# Differentiating Surgeon Expertise Solely by Eye Movement Features

Benedikt Hosp<EMAIL_ADDRESS>https://orcid.org/0000-0001-8259-5463 Human-Computer Interaction, University
of TübingenSand 14TübingenGermany72076 ,  Myat Su Yin  Faculty of ICT, Mahidol
UniversityNakhon PathomThailand<EMAIL_ADDRESS>,  Peter Haddawy
Faculty of ICT, Mahidol UniversityNakon PathomThailand<EMAIL_ADDRESS>,  Paphon Sa-Ngasoongsong<EMAIL_ADDRESS>Faculty of Medicine,
Ramathibodi Hospital, Mahidol UniversityBangkokThailand and  Enkelejda Kasneci
https://orcid.org/0000-0001-8259-5463 Human-Computer Interaction, University
of TübingenSand 14TübingenGermany72076

(2018)

###### Abstract.

Developments in computer science in recent years are moving into hospitals.
Surgeons are faced with ever new technical challenges. Visual perception plays
a key role in most of these. Diagnostic and training models are needed to
optimize the training of young surgeons. In this study, we present a model for
classifying experts, 4th year residents and 3rd year residents, using only eye
movements. We show a model that uses a minimal set of features and still
achieve a robust accuracy of 76.46 % to classify eye movements into the
correct class. Likewise, in this study, we address the evolutionary steps of
visual perception between three expertise classes, forming a first step
towards a diagnostic model for expertise.

surgeon, eye, tracking, diagnostic, model, machine learning

††copyright: acmcopyright††journalyear: 2018††doi:
10.1145/1122445.1122456††conference: ETRA ’21: ACM Symposium on Eye Tracking
Research & Applications ; June 03–05, 2018; Stuttgart, BW††booktitle:
Woodstock ’18: ACM Symposium on Neural Gaze Detection, June 03–05, 2018,
Woodstock, NY††price: 15.00††isbn: 978-1-4503-XXXX-X/18/06††ccs: Social and
professional topics Medical technologies††ccs: Computing methodologies
Supervised learning††ccs: Computing methodologies Biometrics

##  1\. Introduction

Arthroscopy is a popular procedure that promises better patient outcomes on
the one hand and conserve hospital resources on the other. According to Monson
et al. (Monson, 1993), patients experience less pain, have fewer complications
and recover faster. However, a surgeon needs advanced technical skills for
this type of operation (Hermens et al. (Hermens et al., 2013)), for example a
technical understanding and adequate conception of the 3D tissue in a 2D view
and the separation of the field of view (screen) and the instrument plane.
This separation presents a challenge insofar as cognitive projections of
motion from instruments in 3D space must be made into the 2D space of the
screen. Likewise, the instruments determine and constrain the degrees of
freedom of movement. Therefore, in order to reach the appropriate site during
an operation via arthroscopy, navigation plays a central role for the surgeon.
This involves inserting a camera through a portal into the patient’s shoulder.
The camera can rotate in multiple dimensions and casts its image on a screen
placed next to the patient, so that surgeons largely rely on this image during
surgery. The complex anatomy as well as the difficult perception due to
rotation of the camera into 3D tissue, pose a challenge. Due to these
technical challenges, there is a growing interest within the scientific
community to optimize training. In particular,the role of eye movements is
increasingly being investigated. Therefore, diagnostic systems are needed to
improve navigation via perception, especially for non-experts and to
understand the perception of surgeons during surgery and despite that to draw
conclusions for training. In order to find out whether a perceptual-cognitive
training system is useful, first it must be determined whether, and to what
extent, differences in surgeons’ expertise are reflected by their eye
movements. The findings from this diagnostic are significant for the design of
adequate training scenarios for perceptual-cognitive diagnostic and training
systems.

In this work we consider the perception of surgeons using eye movement
patterns from three expertise levels during an arthroscopic surgery of the
shoulder. We use stimulus-independent eye movement patterns to classify the
expertise of the subjects and to investigate differences in the eye movement
behavior of the subjects in detail. These findings are not only useful for
training, but also for assigning surgeons with different skills to different
complex operations.

A fundamental aspect of perception research is the validity of the
presentation as a highly natural presentation allows the participants to apply
their natural gaze behaviour. Many studies were conducted using artificial
forms of presentation like virtual reality (VR) (Law et al., 2004; Zheng et
al., 2011; Wilson et al., 2011; Atkins et al., 2012) or images (Sodergren et
al., 2010; Eivazi et al., 2012). While this kind of presentation allows a
certain level of control for the experimenter, artificially created
presentations are hiding important influences the participant needs to infer.
Zheng et al. (Zheng et al., 2011) for instance report higher levels of
frustration because the physical demands in their study were low, meaning
their simulator might prevent optimal settings for expertise research. To
provide a presentation mode that is as natural as possible, we use so called
soft cadavers that provide natural tactile sensation and therefore increase
the naturality of the scene. Another important aspect in perception research
is the level of control.

In lab studies, remote eye trackers are usually used (Law et al., 2004; Atkins
et al., 2013), since remote eye trackers typically have a screen as stimulus
output and collect data in relation to this input screen, a comparison between
participants is much easier as the control about the stimulus is totally in
the experimenters hands. A week point of remote eye trackers is their
flexibility. As soon as the participant is looking into another direction
(i.e. down at the cadaver), remote eye trackers cannot capture the gaze signal
anymore. To allow the participant to apply a normal gaze behaviour and move
freely without data lost, we use a head-mounted eye tracker in combination
with a 4k-screen. This setup supports natural gaze behaviour as well as high
control on the stimulus. Thus, we can capture highly detailed information of
the tissue on a screen with high resolution and gaze signals on the cadaver,
both with the same field camera.

As eye trackers provide a direct measuring method, they generally allow a high
experimental control. Eye trackers, in fact, provide a sound foundation with
high temporal as well as spatial resolution to research on perceptual
processes. However, studies in surgery also differ in how they evaluated the
gaze signal. The gaze signal on the stimulus was considered, i.e. target gaze
behavior, switching behavior (alternating gaze between target and instrument),
or following behaviour (eye following the instrument) (Law et al., 2004).
Other studies focused on quiet eye periods (Wilson et al., 2011). However,
there are also studies that have gained insights at feature level. For
example, Kocak et al.(Kocak et al., 2005) used stimulus-independent eye
features in their binary classification and found significantly lower saccade
rates as well as significantly higher peak velocities for experts, which was
confirmed by other studies, e.g., by Hermens et al.(Hermens et al., 2013).
Tien et al. (Tien et al., 2011) found a higher fixation rate in experts.
Richstone et al. (Richstone et al., 2010) approved a binary classification
with a mixture of fixation rates, vergence, blink rates and pupil diameter
variability. Eivazi et al. (Eivazi et al., 2012) show differences in time to
first fixation and mean fixation duration. However, theses differences were
not approved by Sondergren et al. (Sodergren et al., 2010), as in both studies
fixation durations are analyzed differently and the choice of regions of
interest plays an important role. These results show that eye movements can be
used to assess the surgical expertise and to define differences between
groups. A wide range of eye tracking features has been shown to reflect
expertise.

In this regard, many studies have focused on the detection of differences in
expertise between experts and novices (Wilson et al., 2011; Zheng et al.,
2011; Tien et al., 2010, 2011; Schulz et al., 2013) and only few studies have
focused on the development of eye movements. However, studies focusing on
development used mostly simulators (Kocak et al., 2005) or images (Sodergren
et al., 2010). By looking at different expertise stages and the
characteristics of their eye movements, we seek to make statements about the
development of gaze behavior through the changes and differences in eye
movements of different classes. We only consider relative features and avoid
absolute features like operation length as such can distort the picture. Thus,
we use only eye movement features to classify perception in isolation and
apply the same methods like Hosp et al. (Hosp et al., 2020) in their study
with athletes.

Machine learning is the state-of-the-art when it comes to big data. As
traditional methods quickly reach their limits when complex data are involved,
machine learning methods are increasingly finding their way into the analysis
of eye tracking data. Another point is that machine learning techniques are
data-driven and therefore can learn patterns from previous data. Usually, the
more data there are available for training, the more robust is the result. So
far, several algorithms have been introduced to eye tracking including
supervised methods like support vector machines (Tong and Chang, 2001; Hosp et
al., 2020; Castner et al., 2018) or neural networks (Castner et al., 2020).
However, hidden markov models (HMM) have also been used for this purpose. In
the study by Sodergren et al. (Sodergren et al., 2010) it was shown that their
HMM can reveal differences in eye movement patterns between high and low
performers. Ahmidi et al. (Ahmidi et al., 2010) mixed instrument movements
with eye movement data and achieved a binary classification accuracy of 82.5%
for skill level classification and confirmed their result with a second study
(Ahmidi et al., 2012). All these studies show that eye movement data can be
used to differentiate between experts and novices and that it is not necessary
to determine where exactly the surgeons were looking at to measure their skill
accurately.

This work is based on the data of Yin et al. (Yin et al., 2020), in the sense
that the same eye tracking data are used for our classification. In contrast,
however, the evolution of eye movement patterns is done over three classes,
namely, novice (residents, 3rd year), intermediates (residents,4th year) and
experts (fellow surgeons). In addition, we look at the characteristics of the
features in order to obtain statements about the perception of the different
groups. Since we want to look at perception in isolation and draw conclusions
about expertise, our evaluation is completely limited to the use of relative
eye movement features. No other features like operation length, number of
saccades/fixations etc. are used.

##  2\. Participants and Methods

###  2.1. Procedure

The data of this study belong to participants from three different groups of
surgeons. The expert group contains n = 5 fellow surgeons with 4-10 years of
experience from the Orthopaedics Faculty of Medicine, Mahidol University. N=5
participants are in the 4th year of the Ortopaedic Surgery Residency Program
and belong to the intermediates group. The novice group contains n=5
participants, who are in the 3rd year of the same program. All residents have
no arthroscopic experience. The study was approved by the Mahidol Unversity
Institutional Review Board. Participants were placed in front of the cadaver
and four feet away from the 4k, 52-inch screen where the output of the
arthroscope is displayed (see Figure 1 ). Each participant was familiarized
with the setup and asked to navigate to 12 anatomical landmarks in the
shoulder, while having a Tobii Pro Glasses 2 eye tacker on. The gaze was
recorded with Tobii software.

Figure 1. Study setup with 4k-screen four feet away from the surgeon and
cadaver directly in front of surgeon.

###  2.2. Data preparation

The Tobii Glasses 2 were set to a frame rate of 100 Hz. This means a gaze
sample should be available every 10 ms and is saved with a timestamp, x- ,and
y-coordinates. These samples are used to calculate metrics like fixations and
saccades. To conduct this process we used the Tobii Fixation Filter, which
behaves mainly like a I-VT filter (velocity threshold based) but uses a
sliding window averaging method. The feature calculation is based on the
classification algorithm of Olson (Olsson, 2007). We used the default values,
which result in a velocity threshold of 0.7 pixels/ms. The raw eye tracking
data, as well as the fixations and saccade metrics, can be exported from the
Tobii Studio software. As our aim is to infer which features contribute to
expertise differences, we first use all the exported features from Tobii
Studio Software and added common metrics, too. Subsequently, we evaluated
their frequency in the model building process and rated the most frequent
features to build a model with this subset of features for expertise
acquisition. For this classification we focus first on the following features:

  * •

Velocity of saccades (average, min, max. std. dev.)

  * •

Amplitude of saccades (average, min, max. std. dev.)

  * •

Total amplitude of saccades

From the velocity of the saccades, as well as the amplitude of saccades, we
use the person-specific average, minimum, maximum and standard deviation as
features. Next to the total amplitude of saccades Tobii provides metrics about
the first saccade and first fixation. We did not include them since our
participants have been familiarized with the glasses in different length and
therefore we end up with chaotic first saccades, which have no informative
character, when the trial starts.

We added certain typical eye movement features which we calculated by
ourselves.

  * •

Saccade duration (average, min, max. std. dev.)

  * •

Fixation duration (average, min, max. std. dev.)

  * •

Fixation frequency

  * •

Saccade frequency

  * •

Pupil diameter (average, min, max. std. dev.)

  * •

Gyroscope X,Y,Z (average, min, max. std. dev.)

We decided to include the gyroscope values as we consider a lot of information
about head movement between screen and cadaver may be revealed. The
integration of pupil diameter features is based on the assumption, that
experts may have less fluctuating puil diameter since there mental effort is
considered to be smaller. Vice versa, the pupil diameter of intermediates and
novices may reveal expertise differences by such effects. Our first model is
built on these mentioned 35 features.

###  2.3. Machine learning model

We took all 34 features and built a SVM model. We protected the model against
overfitting by partitioning the data set into folds and estimate the accuracy
on each fold by applying a 12-fold cross-validation. On each of 1000 runs we
kept out one participant of each group. These participants have never been
seen by the model before (leave-one out validation). Therefore, we trained the
model with the remaining four participants of each group and evaluated each
run with the three participants of the leave-one out validation data set. The
membership of a participant to validation or training set is decided randomly
in each run.

Firstly, we used all 34 features to check the classifiability of the data set
and afterwards reduced the amount by taking the 4 most frequent features of
1000 runs. The most frequent features are features that have the highest
importance values for a single model prediction. In each run we built a queue
of all 34 features sorted by importance for the current model. Subsequently,
we computed their overall frequency over all models.

##  3\. Results

Table 1. The most important and frequent features on 1000 runs. Feature | derivation
---|---
Peak velocity of saccades | standard deviation
Amplitude of saccades | minimum
Total amplitude of saccades | sum
Gyroscope z | min

Our first classifiability model shows promising results with an average
accuracy of 60%. As a system that would simply guess the class, would only
reach a chance-level of 33.33%, the all feature model can already be
considered as well-performing. But as we want to specify the results to allow
a precise statements about a high performing classification with the least
amount of features, we continued by collecting all features and their
importance values, shown in table 1, on 1000 runs of the all feature model.
With a subset of 4 features, earlier counteracting features may be avoided and
a precise statement about the differences of the groups can be stated.

###  3.1. Performance metrics

Figure 2. Performance values on 1000 runs.

Boxplot in Figure 2 shows in the first boxplot to the left, the distribution
of the accuracy values on 1000 runs. In 430 runs out of 1000 the model reached
an accuracy of 90-100%. The accuracy of 446 out of 1000 runs lays between
60-70 %. Only 112 out of 1000 runs show accuracy values of 30-40% and 12 out
of 1000 runs show 0-12%. This means that 98.8% of the runs were better than
chance level and 87.7% still better than twice the chance-level. The dashed
red line shows the average accuracy of 76.46% of the model, based on the 4
most frequent features. This is a solid value for a three class
classification. The average accuracy of the model is more than twice the
chance-level. As the accuracy of the model is a more general information
provider, we calculated other metrics, too. To know about correct predictions,
we calculated the precision or positive predictive value, shown in Figure 2,
second to left. The precision indicates the proportion of results correctly
classified as positive in the totality of results classified as positive. It
is a metric that tells about how precise the predictions of the model are. In
other words, the probability of a sample being predicted to belong to class x
which really belongs to class x. Out of 1000 runs, 430 runs had a precision of
100%. Only 12 of the runs were between 0-10%. The same behaviour can be seen
for the recall (Figure 2). Recall, or hit rate, is a measure that tells about
the probability that a positive object is correctly classified as positive.
Recall is especially important as a model metric, when there are high costs
associated with false negatives. The fourth boxplot to the left shows the miss
rate of the model 2. Out of 1000 runs, the model had a miss rate between
0.05-0.2% on 480 runs. On 380 runs the miss rate was even smaller at
0.0-0.05%. Accordingly, the false negative rate or miss rate indicates the
proportion of objects falsely classified as negative out of the total number
of positive objects. The fifth boxplot in Figure 2 shows F1-scores. With
f1-scores we have a metric that balances between precision and recall. The
F1-score combines precision and recall by means of the weighted harmonic mean,
where precision and recall are weighted equally. The performance of the model
is balanced well, as can be seen in the boxplot. The average f1-score is at
50%. 480 out of 1000 runs achieved a f1-score of 50%-60%. Another 380 runs
achieved even 90%-100%.

###  3.2. Feature evolution

As we can see, that the model is performing very well and therefore able to
classify the three groups of expertise, we have a look at the feature
characteristics as well to analyze their evolution. We do that with the 4 most
frequent features from Table 1. Table 2 contains the characteristics of the 4
MFF features. The table shows that experts have a smaller standard deviation
of the peak velocity of the saccades (93.26 °/s). This feature is hard to
interpret, but one assumption may be, that experts have a more uniformly
distribution of saccade velocites. This means they do more saccades at the
same speed, compared to intermediates and novices. Interestingly,
intermediates as the middle class between expert and novice show a much more
divers saccade peak velocity behaviour (121.72 °/s). Novices are in the middle
between experts and intermediates. A higher value for the standard deviation
of the saccade peak velocities could a an indicator for a more chaotic gaze
behaviour, but it is hard to draw a conclusion on such feature. When having a
look at the minimum saccade amplitudes, we can see the same differences. The
experts have on average a bigger minimal saccade of the length of 0.86 °,
compared to the intermediates with 0.40°and the novices with 0.64°. Again, we
can see that the novices are in between the experts and intermediates. Only
the total amplitude of all saccades show a uniform evolution. The experts do a
total of 481.32°of saccade length, where intermediates do more than twice the
experts (1120.74°) and novices (1956.21°) even more than five time the experts
and nearly double the intermediates. Another interesting feature evolution can
be seen at the gyroscope minimum measure in z-axis. This corresponds to head
rotations to left and right. With -66.47°the intermediated have the smaller
average value of z-axis rotation, followed by the experts with -72.90°and the
novices with -80.12 °.

Table 2. Average feature evolution between classes. Feature | Expert(Fellow) | Intermediate (R4) | Novice (R3)
---|---|---|---
Sasccade peak velocity (std. dev.) | 93.26 °/s | 121.72 °/ s | 117.45 °/s
Saccade amplitude (min) | 0.86 ° | 0.40 ° | 0.64 °
Total saccade amplitude | 481.32 ° | 1120.74 ° | 1956.21 °
Gyroscope Z (min) | -72.90 ° | -66.47 ° | -80.12 °

##  4\. Discussion

In this work we developed a model with supervised machine learning techniques,
that is able to distinguish three levels of expertise solely on eye movements
during an arthroscopic surgery of the shoulder. With an accuracy of 76.46% the
model can be considered as performing well, since it is two times better than
guessing. Even for a binary classification, this result would be high. We
showed that different performance metrics (accuracy, precision, recall, miss
rate and f1-score) allow to state that the model is robust, precise and
valuable to be used in a training system. To further understand the
differences between the three groups of expertise, we had a look at the 4 most
frequent features of the model and analyzed the evolution of the
characteristics between the groups. Except for the total amount of saccade
amplitudes, the three of the four most frequent features show an uniform
evolution. First, novices tend to have a more chaotic gaze behaviour and
distribute their gaze over a bigger amount of the scene by making lots of
different saccades with different speed. They also tend to look more on the
outside than on the center and move their head more from left to right and
vice versa, which might indicate a problem in the cognitive projection of the
2D view of the arthroscope into the 3D-tissue in front of them. The evolution
to intermediates shows an atypical behaviour, as they tend to still gaze over
a bigger area of the scene than the experts, but do smaller saccades with a
still divers velocity. They move their head less. This might indicate, that
they try to focus on more specific visual clues and start to concentrate on
the center of the scene. In the next evolution step, the saccade velocities
shrink significantly, which signifies a more planned scanning behaviour, with
little longer saccades, concentrated more on specific areas but with more head
movements to the left or right. To summarize our findings, one can state that
the evolution of novices to intermediates first tend to lead to a partly more
chaotic gaze behaviour and then turning to be more precise.

Further steps are to add more participants to the single classes, and refine
the number of classes. This allows a much finer classification and therefore a
better understanding of the differences between the classes. A finer
classification is important to robust assumptions made by the model about gaze
behaviour. In future we plan to add more participants from more classes to
investigate on expertise differences in eye movements.

## References

  * (1)
  * Ahmidi et al. (2010) Narges Ahmidi, Gregory D Hager, Lisa Ishii, Gabor Fichtinger, Gary L Gallia, and Masaru Ishii. 2010\.  Surgical task and skill classification from eye tracking and tool motion in minimally invasive surgery. In _International Conference on Medical Image Computing and Computer-Assisted Intervention_. Springer, 295–302.
  * Ahmidi et al. (2012) Narges Ahmidi, Masaru Ishii, Gabor Fichtinger, Gary L Gallia, and Gregory D Hager. 2012.  An objective and automated method for assessing surgical skill in endoscopic sinus surgery using eye-tracking and tool-motion data. In _International forum of allergy & rhinology_, Vol. 2. Wiley Online Library, 507–515.
  * Atkins et al. (2012) M Stella Atkins, Xianta Jiang, Geoffrey Tien, and Bin Zheng. 2012\.  Saccadic delays on targets while watching videos. In _Proceedings of the symposium on eye tracking research and applications_. 405–408.
  * Atkins et al. (2013) M Stella Atkins, Geoffrey Tien, Rana SA Khan, Adam Meneghetti, and Bin Zheng. 2013.  What do surgeons see: capturing and synchronizing eye gaze for surgery applications.  _Surgical innovation_ 20, 3 (2013), 241–248.
  * Castner et al. (2018) Nora Castner, Enkelejda Kasneci, Thomas Kübler, Katharina Scheiter, Juliane Richter, Thérése Eder, Fabian Hüttig, and Constanze Keutel. 2018.  Scanpath comparison in medical image reading skills of dental students: distinguishing stages of expertise development. In _Proceedings of the 2018 ACM Symposium on Eye Tracking Research & Applications_. 1–9.
  * Castner et al. (2020) Nora Castner, Thomas C Kuebler, Katharina Scheiter, Juliane Richter, Thérése Eder, Fabian Hüttig, Constanze Keutel, and Enkelejda Kasneci. 2020.  Deep semantic gaze embedding and scanpath comparison for expertise classification during OPT viewing. In _ACM Symposium on Eye Tracking Research and Applications_. 1–10.
  * Eivazi et al. (2012) Shahram Eivazi, Roman Bednarik, Markku Tukiainen, Mikael von und zu Fraunberg, Ville Leinonen, and Juha E Jääskeläinen. 2012\.  Gaze behaviour of expert and novice microneurosurgeons differs during observations of tumor removal recordings. In _Proceedings of the Symposium on Eye Tracking Research and Applications_. 377–380.
  * Hermens et al. (2013) Frouke Hermens, Rhona Flin, and Irfan Ahmed. 2013.  Eye movements in surgery: A literature review.  (2013).
  * Hosp et al. (2020) Benedikt Hosp, Florian Schultz, Enkelejda Kasneci, and Oliver Höner. 2020.  Eye Movement Feature Classification for Soccer Expertise Identification in Virtual Reality.  _arXiv preprint arXiv:2009.11676_ (sep 2020).
  * Kocak et al. (2005) Ergun Kocak, Jan Ober, Necip Berme, and W Scott Melvin. 2005\.  Eye motion parameters correlate with level of experience in video-assisted surgery: objective testing of three tasks.  _Journal of Laparoendoscopic & Advanced Surgical Techniques_ 15, 6 (2005), 575–580.
  * Law et al. (2004) Benjamin Law, M Stella Atkins, Arthur E Kirkpatrick, and Alan J Lomax. 2004.  Eye gaze patterns differentiate novice and experts in a virtual laparoscopic surgery training environment. In _Proceedings of the 2004 symposium on Eye tracking research & applications_. 41–48.
  * Monson (1993) JR Monson. 1993\.  Advanced techniques in abdominal surgery.  _British Medical Journal_ 307, 6915 (1993), 1346–1350.
  * Olsson (2007) Pontus Olsson. 2007\.  Real-time and offline filters for eye tracking.
  * Richstone et al. (2010) Lee Richstone, Michael J Schwartz, Casey Seideman, Jeffrey Cadeddu, Sandra Marshall, and Louis R Kavoussi. 2010.  Eye metrics as an objective assessment of surgical skill.  _Annals of surgery_ 252, 1 (2010), 177–182.
  * Schulz et al. (2013) Christian M Schulz, Mica R Endsley, Eberhard F Kochs, Adrian W Gelb, and Klaus J Wagner. 2013\.  Situation Awareness in AnesthesiaConcept and Research.  _Anesthesiology: The Journal of the American Society of Anesthesiologists_ 118, 3 (2013), 729–742.
  * Sodergren et al. (2010) Mikael Hans Sodergren, Felipe Orihuela-Espina, James Clark, Ara Darzi, and Guang-Zhong Yang. 2010\.  A hidden markov model-based analysis framework using eye-tracking data to characterise re-orientation strategies in minimally invasive surgery.  _Cognitive processing_ 11, 3 (2010), 275–283.
  * Tien et al. (2010) Geoffrey Tien, M Stella Atkins, Bin Zheng, and Colin Swindells. 2010.  Measuring situation awareness of surgeons in laparoscopic training. In _Proceedings of the 2010 symposium on eye-tracking research & applications_. 149–152.
  * Tien et al. (2011) Geoffrey Tien, Bin Zheng, and M Stella Atkins. 2011.  Quantifying surgeons’ vigilance during laparoscopic operations using eyegaze tracking.. In _MMVR_. 658–662.
  * Tong and Chang (2001) Simon Tong and Edward Chang. 2001.  Support vector machine active learning for image retrieval. In _Proceedings of the ninth ACM international conference on Multimedia_. 107–118.
  * Wilson et al. (2011) Mark R Wilson, John S McGrath, Samuel J Vine, James Brewer, David Defriend, and Richard SW Masters. 2011.  Perceptual impairment and psychomotor control in virtual laparoscopic surgery.  _Surgical endoscopy_ 25, 7 (2011), 2268–2274.
  * Yin et al. (2020) Myat Su Yin, Peter Haddawy, Benedikt Hosp, Paphon Sa-ngasoongsong, Thanwarat Tanprathumwong, Madereen Sayo, Supawit Yangyuenpradorn, and Akara Supratak. 2020.  A Study of Expert/Novice Perception in Arthroscopic Shoulder Surgery. In _Proceedings of the 4th International Conference on Medical and Health Informatics_. 71–77.
  * Zheng et al. (2011) Bin Zheng, Geoffrey Tien, Stella M Atkins, Colin Swindells, Homa Tanin, Adam Meneghetti, Karim A Qayumi, and O Neely M Panton. 2011\.  Surgeon’s vigilance in the operating room.  _The American Journal of Surgery_ 201, 5 (2011), 673–677.


# Supergranules as Probes of the Sun’s Meridional Circulation

David H. Hathaway  NASA Marshall Space Flight Center, Huntsville, AL 35812 USA
<EMAIL_ADDRESS>
###### Abstract

Recent analysis revealed that supergranules (convection cells seen at the
Sun’s surface) are advected by the zonal flows at depths equal to the widths
of the cells themselves. Here we probe the structure of the meridional
circulation by cross-correlating maps of the Doppler velocity signal using a
series of successively longer time lags between maps. We find that the
poleward meridional flow decreases in amplitude with time lag and reverses
direction to become an equatorward return flow at time lags $>24$ hours. These
cross-correlation results are dominated by larger and deeper cells at longer
time lags. (The smaller cells have shorter lifetimes and do not contribute to
the correlated signal at longer time lags.) We determine the characteristic
cell size associated with each time lag by comparing the equatorial zonal
flows measured at different time lags with the zonal flows associated with
different cell sizes from a Fourier analysis. This association gives a
characteristic cell size of $\sim 50$ Mm at a 24 hour time lag. This indicates
that the poleward meridional flow returns equatorward at depths $>50$ Mm –
just below the base of the surface shear layer. A substantial and highly
significant equatorward flow (4.6 $\pm$ 0.4 m s-1) is found at a time lag of
28 hours corresponding to a depth of $\sim 70$ Mm. This represents one of the
first positive detections of the Sun’s meridional return flow and illustrates
the power of using supergranules to probe the Sun’s internal dynamics.

convection, Sun: rotation

##  1 INTRODUCTION

Supergranules are cellular flows observed in the Sun’s photosphere as a
surface filling pattern of largely horizontal flows (Rieutord & Rincon, 2010).
Typical supergranules have diameters of $\sim 30$ Mm and lifetimes of $\sim 1$
day – much larger and longer lived than granules (diameters of $\sim$ 1 Mm and
lifetimes of $\sim$ 5 minutes).

Since their discovery by Hart (1954) and their initial characterization by
Leighton et al. (1962), supergranules have been recognized as the principal
organizer of the Sun’s chromospheric/magnetic network and as key to the
transport of magnetic flux in the photosphere (Leighton, 1964; DeVore et al.,
1984).

The rotation rate of the supergranules has been measured by a number of
investigators with interesting results. Duvall (1980) cross-correlated the
Doppler velocity pattern from equatorial spectral scans obtained over several
days and found that the pattern rotates about 3% faster than the photospheric
plasma. Faster rates are found for the 24-hr time lags from day-to-day than
for the 8-hr time lags from the beginning to end of an observing day. Beck &
Schou (2000) measured the rotation of the Doppler velocity pattern using a
Fourier technique and found that the larger cells rotate more rapidly than the
smaller cells but with apparent rotation rates that exceed the peak internal
rotation rate. This “superrotation” was attributed to line-of-sight projection
effects on the Doppler velocity pattern by Hathaway et al. (2006). In fact,
Schou (2003) had earlier found rotation consistent with the internal plasma
flow profile when this line-of-sight projection was removed.

These observations are consistent with the original conclusion by Duvall
(1980) that larger cells dominate the longer time lags and that these larger
cells are more deeply anchored in a surface shear layer in which the rotation
rate increases with depth. While it should be noted that some characteristics
of the motions of the supergranules have been attributed to wave-like
properties (Gizon et al., 2003; Schou, 2003), it is clear that the cellular
structures are advected largely by the zonal flows in the surface shear layer.

The existence of this shear layer was suggested by Foukal & Jokipii (1975) as
a consequence of the conservation of angular momentum by convective elements
moving inward and outward in the near surface layers. Helioseismology
(Thompson et al., 1996; Schou et al., 1998) now indicates that the shear layer
extends to a depth of $\sim 50$ Mm.

Recently Hathaway (2012) found a one-to-one correspondence between the
rotation rate of supergranules with increasing size and the rotation rate in
this surface shear layer with increasing depth. Results, using a Fourier
technique like that of Beck & Schou (2000) with the line-of-sight projection
effects removed using the method of Schou (2003), indicate that supergranules
are advected by the flows in the shear layer at depths equal to their widths.
We expect the advection of the supergranules to be dominated by flows near
their bases since the plasma density increases rapidly with depth. However,
this does imply that the cells must extend even deeper than this anchoring
level.

Note that while cells that are as deep as they are wide are typical in
numerical simulations of solar convection (Stein et al., 2011), results from
local helioseismology have suggested flattened supergranules – but with
conflicting results. Duvall (1998) found that the outflows in the
supergranules reversed at a depth of 5 Mm and then disappeared by 8 Mm depth.
Zhao & Kosovichev (2003) found reversal at a depth of 8 Mm and disappearance
at 15 Mm. Švanda et al. (2009) found disappearance at an even greater depth –
25 Mm.

Here we measure the meridional motion of the pattern of supergranules by
analyzing the same data used in Beck & Schou (2000), Schou (2003), and
Hathaway (2012).

The Sun’s meridional flow is extremely weak when compared to the other
photospheric flows. Both the strength and the direction of the meridional flow
were uncertain until data and data analysis techniques were sufficiently
improved (Komm et al., 1993; Hathaway et al., 1996; Giles et al., 1997;
Ulrich, 2010). The meridional flow is now known to be poleward in each
hemisphere with a peak velocity of about 10-20 m s-1 while the zonal flows
(differential rotation) are an order of magnitude stronger and the convective
flows (granules and supergranules) are yet another order of magnitude
stronger.

The Sun’s meridional flow plays a critical role in the transport of magnetic
flux (DeVore et al., 1984; Wang et al., 2009) and in some dynamo models for
the Sun’s magnetic cycle (Dikpati & Charbonneau, 1999; Nandy & Choudhuri,
2002). Although the nonaxisymmetric convective motions transport magnetic
elements at much higher speeds to the cell boundaries, those motions are in
random directions and give diffusion away from flux concentrations. The net
poleward transport has a larger contribution from the slow, but direct,
meridional flow. This poleward transport is directly responsible for the
reversal of the Sun’s polar magnetic fields every solar cycle and for
producing the polar fields that determine the strength of the following solar
cycle (Schatten et al., 1978; Svalgaard et al., 2005).

Fourier techniques (which do give explicit information on the cell sizes) are,
unfortunately, not available for measuring the meridional motion of the
supergranule Doppler velocity pattern. Here we use cross-correlations with
time lags that vary from 1 hour to 28 hours. The displacement of the peak in
the cross-correlation gives both zonal and meridional velocities which vary
systematically with time lag. We determine the characteristic cell size for
each time lag by comparing the zonal velocities with those obtained with the
Fourier technique (Hathaway, 2012). This then gives the meridional flow at a
series of anchoring depths corresponding to the series of time lags.

##  2 DATA PREPARATION

The data consist of $1024^{2}$ pixel images of the line-of-sight velocity
determined from the Doppler shift of a spectral line due to the trace element
nickel in the solar atmosphere by the Michelson Doppler Imager (Scherrer et
al., 1995, MDI) on the ESA/NASA Solar and Heliospheric Observatory (SOHO).

The images are acquired at a 1 minute cadence and cover the full visible disk
of the Sun. We average the data over 31 minutes with a Gaussian weighting
function which filters out variations on time scales less than about 16
minutes, and sample that data at 15 minute intervals. We then map these
temporally filtered images onto a $1024^{2}$ grid in heliographic latitude
from pole to pole and in longitude $\pm 90\arcdeg$ from the central meridian
(Figure 1). This mapping accounts for the position angle of the Sun’s rotation
axis relative to the imaging CCD and the tilt angle of the Sun’s rotation axis
toward or away from the spacecraft. Both of these angles include modifications
in line with the most recent determinations of the orientation of the Sun’s
rotation axis (Beck & Giles, 2005; Hathaway & Rightmire, 2010). We analyze
data obtained during two 60 day periods of continuous coverage: one in 1996
from May 24 to July 22; the other in 1997 from April 14 to June 17.

Each filtered image is analyzed using the methods described by Hathaway (1992)
to determine and remove the stationary signals due to differential rotation,
convective blue shift, and the meridional flow – as well as the time varying
signal due to the motion of the spacecraft. Following this, an average image
from all of the images in each 60 day period was subtracted from each image to
remove any instrumental artifacts. This leaves behind the non-axisymmetric and
time-varying signal due to the convective flows.

Figure 1:  Heliographic map detail of the line-of-sight (Doppler) velocity
from SOHO/MDI. The map detail extends $90\arcdeg$ in longitude from the
central meridian on the left and about $35\arcdeg$ in latitude from the
equator (the thick horizontal line). The mottled pattern is the Doppler signal
(blue for blue shifts and red for red shifts) due to the supergranules.

##  3 CROSS-CORRELATION ANALYSIS

We determine the zonal and meridional velocity of the supergranule pattern by
cross-correlating strips of data 11 pixels high in latitude by 601 pixels long
in longitude. The strips are centered on each of the 860 latitude positions
from $75\arcdeg$ south to $75\arcdeg$ north of the equator with overlap from
adjacent strips. These data strips are cross-correlated with corresponding
strips from images obtained at time lags of 1, 2, 4, 8, 16, 24, and 28 hours.

This process is done both forward and backward in time. A data strip centered
on the central meridian is correlated with strips offset to the left at later
times and then data centered on the central meridian at later times is
correlated with strips offset to the right in the initial data map. This
minimizes any systematic errors.

The shift of the Doppler pattern is determined to the nearest pixel by finding
the displacement of the strip in the later Doppler map that gives the maximum
correlation. This pixel-wise shift is then refined by fitting parabolas
through the correlation amplitudes at the displacement that gives the maximum
correlation and those at plus and minus one pixel in each direction. The
maxima of these parabolas then give the Doppler pattern shift to a fraction of
a pixel. The shifts in longitude give the apparent zonal velocity while the
shifts in latitude give the meridional velocity. The 60 days of data in each
data set give some 1400 independent measurements of these velocities at each
of the 860 latitude positions for each time lag.

The meridional flow measurements for time lags from 1 to 24 hours are shown in
Figure 2. Here the data are smoothed in latitude with a tapered Gaussian
having a FWHM of $10\arcdeg$. (A tapered Gaussian is a Gaussian with a
quadratic function subtracted so that the smoothing filter and its first
derivative go to zero at the endpoints.) The flow is poleward in each
hemisphere from the equator up to latitudes of at least $60\arcdeg$. The speed
of the flow decreases as the time lag increases and the poleward flow
virtually vanishes when measured with the 24 hour time lags. The slow down for
time lags from 2 to 16 hours was previously noted by Hathaway et al. (2010).
The lack of a poleward meridional flow for 24 hour time lags was previously
noted (and referred to as “anomalous”) by Gizon et al. (2003).

Figure 2:  Meridional flow profiles for time lags from 1 to 24 hours. The
measurements from the 1996 data are shown with the red line while the
measurements from the 1997 data are shown with the blue line. The dotted lines
indicate $2\sigma$ error limits for each year. The meridional flow is poleward
up to latitudes of at least $60\arcdeg$. The amplitude of this poleward
meridional flow decreases as the time lag increases and disappears at 24 hour
time lags.

Measurements with 28 hour time lags are shown in Figure 3. Again, the data are
smoothed in latitude with a tapered Gaussian having a FWHM of $10\arcdeg$.
With this longer time lag the measured meridional flow is equatorward at
virtually all latitudes with a peak of $\sim 4$ m s-1. The $2\sigma$
variations on the measurements at each latitude are $\sim 2$ m s-1 making
nearly all of the individual latitude measurements statistically significant.
Fitting the measurements to a 4th-order polynomial in $\sin\theta$, shown by
the thick black line in Figure 3, gives a peak velocity of 4.6 $\pm$ 0.4 m s-1
with $2\sigma$ errors. As such this represents the first positive measurement
and significant ($>10\sigma$) detection of the Sun’s equatorward return flow.

Figure 3:  Meridional flow profiles for time lags of 28 hours. The
measurements from the 1996 data are shown with the red line while the
measurements from the 1997 data are shown with the blue line. The dotted lines
indicate $2\sigma$ error limits for each year. The 4th-order fit to the
combined data is shown with the thick, black line.

##  4 DEPTH DETERMINATIONS

Duvall (1980) suggested that at different time lags the cells that survive are
advected by flows with different speed (and direction) depending upon the
depth to which those cells extend. The different time lags thus probe the
depth dependence of the meridional flow. Hathaway (2012) found that
supergranules of a given cell size (given by their wavenumber in a Fourier
analysis) are advected by the zonal flows at depths equal to their widths.
This was determined by comparing their zonal velocities from the Fourier
analysis with zonal velocities in the surface shear layer determined with
global helioseismology by Schou et al. (1998).

Small cells have short lifetimes. Granules with typical diameters of $\sim 1$
Mm have lifetimes of $\sim 5$ minutes (Title et al., 1989). Supergranules with
typical diameters of $\sim 30$ Mm have lifetimes of $\sim 1$ day (Simon &
Leighton, 1964; Wang & Zirin, 1989; Hirzberger et al., 2008). Cells of
intermediate (mesogranular) size with diameters 5-10 Mm have intermediate
lifetimes of $\sim 2$ hours (November et al., 1981). As the time lag used in
the cross-correlation increases, the cells with lifetimes shorter than the
time lag present an uncorrelated random pattern that does not influence the
determination of the pattern motion. Longer time lags are associated with
larger cells which extend deeper into the Sun. This assertion is supported by
observations of the monotonic decline in the cross-correlations for granules
(Title et al., 1989), the monotonic decline in the cross-correlations for the
supergranule Doppler pattern (Hathaway et al., 2010),and the increase in
lifetimes for larger supergranule cells determined from local helioseismology
(Hirzberger et al., 2008).

We can determine the characteristic cell size at each time lag by matching the
zonal velocities from the cross-correlations with those from the Fourier
method. This Fourier method is described in detail by Beck & Schou (2000) and
Hathaway (2012). Lines of mapped Doppler data at equatorial latitudes are
Fourier transformed in longitude to give spectral coefficients. Spectral
coefficients from 6 10-day series of Doppler maps obtained at 15-minute
cadence are then Fourier transformed in time. The location of the peak power
gives the zonal (longitudinal) velocity for each spatial wavelength or cell
size.

The line-of-sight projection for the Doppler signal influences the measured
zonal flow (but, as shown by Hathaway et al. (2010), not the meridional flow).
When this projection is divided-out (as can only be done at the equator) the
zonal flow as a function of wavelength from the Fourier analysis matches the
zonal flow as a function of depth from global helioseismology. Figure 4 shows
zonal velocities relative to the Carrington Rotation frame of reference for
supergranules of varying wavelength along with the zonal velocity as a
function of depth from global helioseismology (Schou et al., 1998).

Figure 4:  Equatorial zonal velocities as functions of cell wavelength and
depth within the Sun. Zonal velocities as functions of wavelength from the
1996 MDI data in which the line-of-sight projection is divided out are shown
with the small black circles. Zonal velocities as functions of depth from
global helioseismology (Schou et al., 1998) are shown with the large red
circles. The match between these measurements at virtually all wavelengths and
depths indicates that supergranules are advected by the flows at depths equal
to their wavelengths.

The zonal velocities for the raw data (without the line-of-sight projection
removed) are shown in Figure 5. Comparing these zonal velocities with those
obtained at the equator (also with the raw data) by the cross-correlation
analysis gives a correspondence between time lag and characteristic wavelength
of the cells that dominate at that time lag. This correspondence is given in
Table 1.

Figure 5:  Equatorial zonal velocities as functions of cell wavelength and depth within the Sun. Zonal velocities as functions of wavelength from the 1996 MDI data are shown with the small black circles. Zonal velocities as functions of depth from helioseismology (Schou et al., 1998) are shown with the large red circles. The zonal velocities at the equator from the different time lags are indicated along the vertical axis with lines indicating where they match the zonal velocities from the Fourier analysis. Time Lag | Velocity | Wavelength
---|---|---
1 hour |  $12\pm 1$ m s-1 |  $7\pm 1$ Mm
2 hours |  $24\pm 1$ m s-1 |  $14\pm 1$ Mm
4 hours |  $40\pm 1$ m s-1 |  $22\pm 1$ Mm
8 hours |  $49\pm 2$ m s-1 |  $26\pm 2$ Mm
16 hours |  $59\pm 2$ m s-1 |  $37\pm 2$ Mm
24 hours |  $71\pm 3$ m s-1 |  $50\pm 5$ Mm
28 hours |  $80\pm 3$ m s-1 |  $70\pm 5$ Mm
Table 1: Correspondence between time lag, equatorial zonal velocity, and cell
wavelength from Figure 5.

##  5 CONCLUSIONS

We conclude that the Sun’s poleward meridional flow is confined to it’s
surface shear layer (from the photosphere to a depth of $\sim 50$ Mm). We draw
this conclusion from the fact that the poleward flow weakens and then
disappears as the cross-correlation time lag increases from 1 to 24 hours
while the measured rotation rates indicate that the cross-correlated patterns
are dominated by cells with sizes, and therefor depths, that increase from
$\sim 7$ Mm to $\sim 50$ Mm. A positive detection of the equatorward return
flow is found with time lags of 28 hours corresponding to depths of $\sim 70$
Mm.

The fact that the return flow signal was significantly detected at a longer
time lag in a deeper layer also indicates that the lack of a poleward flow at
the 24 hour time lag is due to an actual lack of poleward flow rather than an
insensitivity to the signal.

Estimating the poleward mass flux at $30\arcdeg$ latitude within the surface
shear layer using the measured flow speeds and the run of plasma density with
depth, indicates that a balancing equatorward mass flux at a flow speed of
$\sim 4$ m s-1 should only extend another 20-30 Mm below the base of this
surface shear layer.

This shallow return flow is somewhat surprising given earlier results from
helioseismology. Braun & Fan (1998) used frequency shifts between northward
and southward propagating waves and concluded that the meridional flow was
poleward to a depth of at least 100 Mm. Giles (2000) used time-distance
helioseismology with a mass flux constraint and concluded that the return flow
was below a depth of about 140 Mm. However, Duvall & Hanasoge (2009) have
found that these measurements are plagued by systematic errors leading Gough &
Hindman (2010) to state that, with helioseismology, the meridional flow below
a depth of 30 Mm remains uncertain.

Recently Zhao et al. (2012A) have characterized these systematic signals,
removed them from the data, and retrieved a meridional flow profile in the
near surface layers that agrees well with measurements obtained from magnetic
feature motions (Hathaway & Rightmire, 2010) and direct Doppler measurements
(Ulrich, 2010). Zhao et al. (2012B) have now reported that they find evidence
for the meridional return flow at depths much shallower than the tachocline –
supporting the conclusions reached here.

A shallow return flow is also supported by Sivaraman et al. (2010) who
measured the meridional motions of sunspot umbrae. They compared the rotation
rates of sunspot umbrae with rotation rates derived from helioseismology to
find anchoring depths for umbrae of different sizes (similar to what was done
here with supergranules) and found that small sunspots with umbral areas
$<5\mu$Hem (a $\mu$Hem is $10^{-6}$ the area of a solar hemisphere) are
anchored at depths of $\sim 90$ Mm and are moving equatorward at speeds of
$\sim 7$ m s-1. Slightly deeper anchoring depths and slower equatorward
velocities were found for larger sunspots.

##  6 FURTHER IMPLICATIONS

The shallow depth of the meridional return flow has implications for our
understanding of both solar convection zone dynamics and solar dynamo theory.

The axisymmetric meridional flow is maintained by a balance between the
combined effects of three driving forces (latitudinal pressure gradients,
Reynolds stresses resulting from correlations between radial and meridional
motions in the cellular flows, and the Coriolis force acting on the
axisymmetric zonal flow) and the viscous diffusion due to the cellular flows.
Simulations of the effects of the Sun’s rotation on supergranules in plane-
parallel geometry (Hathaway, 1982) indicate that the Reynolds stresses
produced by the supergranules themselves can maintain both the zonal shear
(slower rotation at the surface and faster rotation below) and a meridional
circulation with poleward flow near the surface and equatorward flow below.
The meridional Reynolds stresses (correlated poleward and upward motions in
the cellular flows) are reinforced by the Coriolis force on the axisymmetric
zonal flow in these simulations. (Latitudinal pressure gradients cannot be
maintained or included due to the periodic boundaries in both horizontal
directions.)

The meridional circulations produced by larger convection cells in rotating
spherical shells extending over several density scale heights (but thus far
not including the surface shear layer) are usually highly structured in
latitude and radius and are highly variable in time (Miesch et al., 2000).
More recent spherical shell simulations (Miesch et al., 2008) impose a
latitudinal entropy (temperature) gradient at the base of the convection zone
to give rotation profiles and meridional circulations more in line with
observations but still produce highly time dependent meridional flows. The
source of this structure and variability can be attributed to the small number
($\sim 100$) of convection cells that populate the simulated global volume. A
meridional circulation driven in the surface shear layer by ($\sim 10000$)
supergranules (the convection cells that populate the surface shear layer) is
far more likely to be less structured and variable – as is observed.

The shallow return flow we find here has implications for flux transport
dynamos (Dikpati & Charbonneau, 1999; Nandy & Choudhuri, 2002) – it violates a
key assumption. These dynamo models assume that the return flow begins at a
depth of 100 Mm and reaches a maximum equatorward speed of just 1-2 m s-1 at
the base of the convection zone. This is a critical assumption in that the
equatorward drift of the sunspot zones is directly produced in these dynamo
models by this slow equatorward flow at the base of the convection zone.

The results presented here support the idea of using supergranules to probe
the dynamics of the solar convection zone. This method can complement the
probing done with acoustic waves by the methods of helioseismology and provide
valuable information on the dynamics of solar and stellar convection zones.
The method used here, finding the peak in the cross-correlated Doppler
velocity pattern for data pairs with different time lags, is limited to time
lags less than $\sim 32$ hours. At longer time lags the correlation amplitudes
from individual data pairs are similar to the noise level from the cross-
correlations. However, averaging the cross-correlations over as few as 24
hourly measurements could extend the results to longer time lags and deeper
layers. Local Correlation Tracking using smaller patches might also yield
measurements of the large scale nonaxisymmetric flows – Giant Cells.

Finally, it should be noted that, while these results indicate an equatorward
return flow at a fairly shallow depth, they do not preclude the possible
existence of additional cells of poleward and equatorward flow at even greater
depths.

The author would like to thank NASA for its support of this research through
grants from the Heliophysics Causes and Consequences of the Minimum of Solar
Cycle 23/24 Program and the Living With a Star Program to NASA Marshall Space
Flight Center. He is indebted to Ron Moore and Lisa (Rightmire) Upton and an
anomymous referee whose comments greatly improved the manuscript and to John
Beck who produced the temporally filtered data from the original MDI Doppler
data. He would also like to thank the American taxpayers who support
scientific research in general and this research in particular. SOHO, is a
project of international cooperation between ESA and NASA.

## References

  * Beck & Giles (2005) Beck, J. G., & Giles, P. 2005, ApJ 621, L153
  * Beck & Schou (2000) Beck, J. G., & Schou, J. 2000, Sol. Phys. 193, 333
  * Braun & Fan (1998) Braun, D. C., & Fan, Y. 1998, ApJ 508, L105
  * DeVore et al. (1984) DeVore, C. R., Sheeley, N. R., Jr., & Boris, J. P. 1984, Sol. Phys. 92, 1
  * Dikpati & Charbonneau (1999) Dikpati, M., & Charbonneau, P. 1999, ApJ 518, 508
  * Duvall (1980) Duvall, T. L., Jr. 1980, Sol. Phys. 66, 213
  * Duvall (1998) Duvall, T. L., Jr. 1998, in ESA Special Publication, Vol. 418, Structure and Dynamics of the Interior of the Sun and Sun-like Stars, ed. S. Korzennik, 581
  * Duvall & Hanasoge (2009) Duvall, T. L., Jr., & Hanasoge, S. M. 2009, in ASP Conf. Series 416 Solar-Stellar Dynamos as Revealed by Helio- and Asteroseismology, eds. M. Dikpati, T. Arentoft, I. González Hernández, C. Lindsey and F. Hill, 103
  * Foukal & Jokipii (1975) Foukal, P., & Jokipii, R. 1975, ApJ 199, L71
  * Giles (2000) Giles, P. M. 2000, PhD Thesis, Stanford Univ.
  * Giles et al. (1997) Giles, P. M., Duvall, T. L., Jr., Scherrer, P. H., & Bogart, R. S. 1997, Nature 390, 52
  * Gizon et al. (2003) Gizon, L., Duvall, T. L., Jr., & Schou, J. 2003, Nature 421, 43
  * Gough & Hindman (2010) Gough, D. O., & Hindman, B. W. 2010, ApJ 714, 960
  * Hart (1954) Hart, A. B. 1954, MNRAS 114, 17
  * Hathaway (1982) Hathaway, D. H. 1982, Sol. Phys. 77, 341
  * Hathaway (1992) Hathaway, D. H. 1992, Sol. Phys. 137, 15
  * Hathaway (2012) Hathaway, D. H. 2012, ApJ 749, L13
  * Hathaway et al. (1996) Hathaway, D. H. et al. 1996, Science 272, 1306
  * Hathaway & Rightmire (2010) Hathaway, D. H. & Rightmire, L. 2010, Science 327, 1350
  * Hathaway et al. (2006) Hathaway, D. H., Williams, & Cuntz, M. 2006, ApJ 644, 598
  * Hathaway et al. (2010) Hathaway, D. H., Williams, P. E., Dela Rosa, K., & Cuntz, M. 2010, ApJ 725, 1082
  * Hirzberger et al. (2008) Hirzberger, J., Gizon, L., Solanki, S. K., & Duvall, T. L., Jr. 2008, Sol. Phys. 251, 417
  * Komm et al. (1993) Komm, R. W., Howard, R. F., & Harvey, J. W. 1993, Sol. Phys. 147, 207
  * Leighton (1964) Leighton, R. B. 1964, ApJ 140, 1547
  * Leighton et al. (1962) Leighton, R. B., Noyes, R. W., & Simon, G. W. 1962, ApJ 135, 474
  * Miesch et al. (2000) Miesch, M. S., Elliott, J. R., Toomre, J., Clune, T. L., Glatzmaier, G. A., Gilman, P. A. 2000, ApJ 532, 593
  * Miesch et al. (2008) Miesch, M. S., Brun, A. S., DeRosa, M. L., Toomre, J. 2008, ApJ 673, 557
  * Nandy & Choudhuri (2002) Nandy, D., & Choudhuri, A. R. 2002, Science 296, 1671
  * November et al. (1981) November, L. J., Toomre, J., Gebbie, K. B. & Simon, G. W. 1981, ApJ 245, L123
  * Rieutord & Rincon (2010) Rieutord, M. & Rincon, F. 2010, Living Rev. Solar Phys. 7, 2
  * Schatten et al. (1978) Schatten, K. H., Scherrer, P. H., Svalgaard, L. & Wilcox, J. M. 1978, Geophys. Res. Lett. 5, 411
  * Scherrer et al. (1995) Scherrer, P. H., et al. 1995, Sol. Phys. 162, 129
  * Schou et al. (1998) Schou, J. et al. 1998, ApJ 505, 390
  * Schou (2003) Schou, J. 2003, ApJ 596, L259
  * Simon & Leighton (1964) Simon, G. W. & Leighton, R. B. 1964, ApJ 140, 1120
  * Sivaraman et al. (2010) Sivaraman, K. R., Sivaraman, H., Gupta, S. S. & Howard, R. F. 2010, Sol. Phys. 266, 247
  * Stein et al. (2011) Stein, R. F., LagerFjärd, A., Nordlund, Å., & Georgobiani, D. 2011, Sol. Phys. 268, 271
  * Svalgaard et al. (2005) Svalgaard, L., Cliver, E. W., & Kamide, Y. 2005, Geophys. Res. Lett. 32, L01104
  * Švanda et al. (2009) Švanda, M., Kosovichev, A. G., Klvaňa, K., Sobotka, M, & Duvall, T. L., Jr. 2009 in ASP Conf. Series 416 Solar-Stellar Dynamos as Revealed by Helio- and Asteroseismology, eds. M. Dikpati, T. Arentoft, I. González Hernández, C. Lindsey and F. Hill, 547
  * Thompson et al. (1996) Thompson, M. J., et al. 1996, Science 272, 1300
  * Title et al. (1989) Title, A. M., Tarbell, T. D., Topka, K. P., Ferguson, S. H., Shine, R. A. & The SOUP Team 1989, ApJ 336, 475
  * Ulrich (2010) Ulrich, R.K. 2010, ApJ 725, 658
  * Wang et al. (2009) Wang, Y.-M., Robbrecht, E. & Sheeley, N. R., Jr. 2009, ApJ 707, 1372
  * Wang & Zirin (1989) Wang, H. & Zirin, H. 1989, Sol. Phys. 120, 1
  * Zhao et al. (2012B) Zhao, J., Bogart, R.S., Kosovichev, A.G., & Duvall, T. L., Jr. 2012B, AAS Meeting 220, 109.05
  * Zhao et al. (2012A) Zhao, J., Nagashima, K., Bogart, R.S., Kosovichev, A.G., & Duvall, T. L., Jr. 2012A, ApJ 749, L5
  * Zhao & Kosovichev (2003) Zhao, J. & Kosovichev, A.G. 2003, in ESA Special Publication, Vol. 517, GONG+ 2002. Local and Global Helioseismology: the Present and Future, ed. H. Sawaya-Lacoste, 417


Center for coherence and quantum optics, Department of Physics, University of
Rochester, 480 Intercampus Dr, Rochester, NY 14627, USA Center for coherence
and quantum optics, Department of Physics, University of Rochester, 480
Intercampus Dr, Rochester, NY 14627, USA Materials Science, University of
Rochester, Rochester, NY 14627, USA Department of Physics and Astronomy,
University of Rochester, Rochester, NY 14627, USA IR,NMR,UV

# Diffusion of brightened dark excitons in a high-angle incommensurate Moiré
homobilayer

Arnab Barman Ray  The Institute of Optics, University of Rochester, 480
Intercampus Dr, Rochester, NY 14627, USA  Trevor Ollis  Department of Physics
and Astronomy, University of Rochester, Rochester, NY 14627, USA  Sethuraj K.
R  The Institute of Optics, University of Rochester, 480 Intercampus Dr,
Rochester, NY 14627, USA  Anthony Nickolas Vamivakas  The Institute of Optics,
University of Rochester, 480 Intercampus Dr, Rochester, NY 14627, USA
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

The last few years have witnessed a surge in interest and research efforts in
the field of twistronics, especially in low-angle twisted bilayers of
transition metal dichalocogenides. These novel material platforms have been
demonstrated to host periodic arrays of excitonic quantum emitters, interlayer
excitons with long lifetimes, and exotic many-body states. While much remains
to be known and understood about these heterostructures, the field of high-
angle, incommensurate bilayers is even less explored. At twist angles larger
than a few degrees, the presence of periodicity in these bilayers becomes
chaotic, making the systems essentially aperiodic and incommensurate in nature
due to the limitations of fabrication techniques. In this work, we demonstrate
the emergence of a brightened dark intralayer exciton in twisted molybdenum
diselenide homobilayer. We show that this dark exciton diffuses across the
excitation spot more efficiently as compared to trions or excitons, reaching
diffusion lengths greater than 4 microns. Temperature-dependent spectra
provide corroborative evidence and reveal a brightened dark trion. Our results
reveal some of the richness of the physics of these high-angle systems.

###### keywords:

2D materials, Moiré superlattice, Excitons, Optoelectronics

Keywords: 2D materials, Moiré superlattice, Intralayer excitons,
Optoelectronics

##  1 Introduction

Moiré hetero- and homo-bilayers of transition metal dichalcogenides(TMDCs)
have been shown to host correlated electronic phenomena 1, 2, 3, 4, 5, 6, 7, 8
and arrays of programmable quantum emitters 9, 10, 11, 12, 13, 14, 15. Less
explored has been the physics of high-angle twisted homobilayers
($10^{\circ}<\theta<50^{\circ})$. While the Moiré superlattice is active only
at small angles 16 and allows for the trapping of excitons, at high angles,
this periodicity is broken. The condition for commensurability or periodicity
in a twisted bilayer system of two honeycomb lattices is provided by the
equation17, $\cos{\theta}=\frac{3m^{2}+3mn+n^{2}/2}{3m^{2}+3mn+n^{2}}$, where
$\theta$ is the twist angle and $m$ and $n$ are a pair of coprime positive
integers. Immediately, it can be discerned that angles where the value of the
cosine is irrational do not exhibit periodicity. However, given that there are
an infinite number of coprime pairs of positive integers, it is possible to
have a rational number expressed in the form of the previous equation
arbitrarily close to the irrational number in question. This leads to a highly
sensitive, if not wholly chaotic period at large angles18. With the
limitations of current fabrication techniques19, this leads to an essentially
aperiodic structure without the band modulating and flattening effects that a
Moiré superlattice usually endows.

Without an active and periodic Moiré superlattice, these high-angle bilayers
may seem to hold little promise in terms of the originally intended
applications in quantum computation for these systems. A high twist angle
pushes diffusive interlayer excitons20 out of the light cone while limiting
the possibility of simulating correlated states in periodic lattices21 due to
an effective uncoupling of electronic communication between the two layers, as
is seen in graphene. However, as we show in this work, this class of aperiodic
bilayers can be interesting in its own right. We focus on a high-angle twisted
homobilayer of molybdenum diselenide. These systems have been explored in a
recent work22 over a range of twist angles. We show, using standard micro-
photoluminescence (PL) experiments, that at large twist angles, the proximity
of a highly polarizable monolayer to the other alters its optical properties.
As reported previously 22, we observe a large redshift of the trionic and
excitonic resonances in the bilayer. This redshift can be partly attributed to
the proximal monolayer’s large polarizability, similar to the much smaller
redshift experienced by a free monolayer when encapsulated by a high
dielectric constant insulator hBN 23. Thus, the large-angle twist of the
bilayer serves to prevent the heterostructure from becoming an indirect
bandgap semiconductor which is the case for a Bernal-stacked ($0^{\circ}$)
bilayer 24, while at the same time having the Moiré superlattice potential
inactive.

Upon investigating the spatial diffusion of bound complexes under steady-state
CW excitation at low intensities across the PL spectrum, we uncover the
quenching of diffusion lengths at energies of large quantum yield or maximum
PL intensity. Investigating diffusion lengths points us toward the presence of
a more diffusive excitonic species at a slightly higher emission energy than
the bright spin-singlet exciton. By carefully deconstructing the spectra, we
uncover that this new species of exciton diffuses more efficiently than the
bright excitons. Investigating the temperature dependence of the PL reveals
evidence for a population transfer to this brightened dark exciton. Population
transfer to higher-energy dark excitonic states is what is responsible for the
decrease in quantum yield with increasing temperature in single monolayers of
$\text{MoSe}_{2}$ and $\text{MoS}_{2}$ with the point group symmetry $D_{3h}$,
and is well documented and understood25, 26, 27. However, we find that for a
high-angle twisted bilayer this population transfer still allows us to capture
some of the PL emitted from the spin-forbidden dark exciton, as a result of
the removal of the out-of-plane reflection symmetry - causing the point group
of both the two monolayers to reduce to $D_{3v}$. The emission due to these
dark excitons is usually suppressed due to the vastly smaller radiative rate
that accompanies spin-flip electronic transitions in monolayers28. This
brightening of this triplet exciton is similar to what has been predicted and
observed for the case of interlayer excitons in
$\text{WSe}_{2}-\text{MoSe}_{2}$ heterobilayers29, 30. Furthermore, we show
that these dark triplet excitons are more diffusive than their bright
counterparts, with diffusion lengths exceeding 4 microns. Dark excitons are
usually long-lived, optically decoupled from the environment, and serve as a
reservoir for their bright counterparts, playing a crucial role in the
condensation of excitons in other well-studied semiconductor platforms31, 32.
Hence it is important to understand their properties.

##  2 Results and Discussion

Figure 1:  (a) Optical micrograph of the homobilayer,(b) co-polarized second
harmonic signal from the bottom monolayer and the bilayer used to estimate the
twist angle, (d) PL spectrum of the monolayer and bilayer with the shifted
trionic (T) and excitonic (X) peaks at 4K, 532 nm CW excitation.

Figure 1 (a) provides an optical micrograph of the homobilayer investigated in
this work. We mechanically exfoliate and assemble monolayers of n-type
$\text{MoSe}_{2}$ (2D semiconductors) and thin, flat flakes of hBN (2D
semiconductors) on a distributed Bragg reflector (see Methods) with its
reflection band centered at 770 nm to optimize the PL signal collection. Using
a pulsed ti-sapphire laser at 790.5 nm, we collect the co-polarized second
harmonic signal generated from the more accessible bottom monolayer and the
homobilayer as a function of the laser polarization angle. After accounting
for the effects of the beamsplitters in the signal and collection path, the
corrected SHG signal is presented in Fig. 1 (b) with their respective fits,
revealing the twist angle to be: $41.27^{\circ}\pm 1.55^{\circ}$. We note that
the SHG signal thus collected also validates the high quality of the
fabricated sample and that there is minimal strain present33, 34 away from the
visible bubbles in the micrograph. We probe the photoluminescence spectra in a
confocal microscopy setup where the sample is cooled to cryogenic temperatures
(all measurements are at 6 K unless specified). Figure 1 (c) shows the PL
signal from the monolayer and bilayer. We observe a large redshift of the
trionic and excitonic resonances of $\Delta_{T}=30.4\text{ meV}$ and
$\Delta_{X}=32.8\text{ meV}$.

Next, we shift our attention to the diffusion of the excitonic and trionic
complexes across the emission spot on both the monolayer and the bilayer under
steady-state excitation. While it has been demonstrated that at small angles
the localizing effects of the Moiré potential impedes the diffusion for
interlayer and intralayer excitons35, 36, 37, the case for high-angle bilayers
is less explored either theoretically or experimentally. We focus on the
diffusion lengths obtained for different species in the monolayer and bilayer.
We note that further experimental work involving measurements of PL lifetime
would help calculate the diffusion coefficients. However, the focus of this
work is to chronicle how diffusion lengths were used to identify the
brightened spin-forbidden dark exciton.

Figure 2:  (a) Spatio-spectrum of the monolayer, PL spectrum and the spectral
variation of diffusion lengths of (b) the trion (T), and (c) the exciton (X).
(d) Spatio-spectrum of the homobilayer and the spectral variation of diffusion
lengths of (b) the trion ($\text{M}_{\text{T}}$), and (c) the exciton
($\text{M}_{\text{X}}$). Deep blue lines trace the obtained diffusion lengths
while light blue shaded regions demarcate the 95 % confidence intervals from
the fits.

Upon excitation with a power of $10\,\mu W$ in a spot of radius $\sim 1\,\mu
m$, we image the emission spot on the CCD camera of our spectrometer setup and
record the spatio-spectrum of the monolayer and the bilayer in Fig. 2(a) and
(d). Under steady-state excitation, neglecting the effect of exciton-exciton
interactions, we fit the spatial extent of the PL intensity as a function of
emission energy using the diffusion equation,
$n(x)\propto\int_{-\infty}^{\infty}dx^{\prime}K_{0}\left(\frac{x^{\prime}}{L_{D}}\right)e^{-\frac{(x-x^{\prime})^{2}}{L^{2}}}$,
where $K_{0}$ is the modified Bessel function of the second kind35, 38, 39. We
determine the laser spot linewidth $L$ by fitting it with a Gaussian function.
This fit provides us with the diffusion lengths $L_{D}$ as a function of
spectra. We verify that the diffusion lengths thus obtained do not change
considerably over three orders of magnitudes of the excitation intensity in
the Supporting Information (see Fig. S1). Note that the diffusion equation is
less appropriate for modeling trion diffusion due to local electric field
effects from donor atoms which can modify their dynamics 40. However, the
quenching of the diffusion lengths for trions at energies corresponding to
high quantum yield (high PL intensity) denoted by green arrows in Fig. 2 (b),
(c), (e) and (f) highlights the accuracy of our measurements. Moreover, our
data captures the fact that trions diffuse less as compared to excitons due to
their larger effective mass41 and the aforementioned effects which is well
documented in the literature40, 42. Figures 2 (b), (c), (e) and (f) show that
the spectral variation of the diffusion lengths qualitatively resembles a
reflection of the PL spectra about the horizontal axis.

For any diffusive bound complex, the rate of population decay is given as a
sum of the radiative and non-radiative rates,
$\Gamma_{tot}=\Gamma_{r}+\Gamma_{nr}$. The PL lifetime is given as
$\tau_{tot}=\frac{\tau_{r}\tau_{nr}}{\tau_{r}+\tau_{nr}}$, where
$\tau_{r}=\Gamma_{r}^{-1}$ and $\tau_{nr}=\Gamma_{nr}^{-1}$. As the diffusion
lengths theoretically are given by $L_{D}=\sqrt{D\,\tau_{tot}}$ (where $D$ is
the diffusion coefficient), substituting this gives,
$L_{D}=\sqrt{D\frac{\tau_{r}\tau_{nr}}{\tau_{r}+\tau_{nr}}}$. Using the
relation satisfied by the intrinsic PL quantum yield given as
$\text{QY}=\frac{\Gamma_{r}}{\Gamma_{r}+\Gamma_{nr}}$, it is straightforward
to arrive at the equation, $L_{D}=\sqrt{\frac{D\,\tau_{r}}{QY}}$. This
relation explains the quenching of the diffusion lengths at energies of high
PL intensity or quantum yield.

Figure 2 (d) and (f) hints at the presence of a less bright species of exciton
(which we label $\text{M}_{\text{D}}$) at a slightly higher energy than the
bright exciton $\text{M}_{\text{X}}$. At this point, it is impossible to
ascertain whether this species is fundamentally different from the bright
exciton or whether it is a result of inhomogeneity or dielectric disorder42 in
the sample. Moreover, while Fig. 2 (e) and (f) seem to indicate that both
$\text{M}_{\text{X}}$ and $\text{M}_{\text{D}}$ diffuse much more efficiently
as compared to $\text{M}_{\text{T}}$ or even the monolayer exciton X, the
spectral proximity of these two species may cause leakage of the tails of
their respective spectrums at the peak energies of each other, and thus
artificially inflating their actual diffusion lengths. To circumvent this
problem, we try to disentangle the contributions of each of these two species.
We find that the spectral slices that make up the spatio-spectrums lend
themselves well to double-lorentzian fits (see supporting Figure S2). We are
thus able to investigate the diffusion lengths of both of these species
separately in Fig. 3.

Figure 3:  (a) Corrected spatio-spectrum of the bilayer trion and bright
exciton and (b) the spectral variation of diffusion length for the bright
exciton. (c) Corrected spatio-spectrum of the bilayer dark exciton and (d) the
spectral variation of diffusion length for the dark exciton. Deep blue lines
trace the obtained diffusion lengths while light blue shaded regions demarcate
the 95 % confidence intervals from the fits.

We next compare and contrast the diffusion lengths of different species in
Fig. 4. For the monolayer, the exciton diffuses more efficiently than the
trions, which face an inward electrostatic force from donor atoms, altering
their diffusion significantly. For the bilayer, our results indicate that
despite the absence of a periodic Moiré superlattice, the diffusion of bright
excitons is suppressed. We suggest that the suppression of excitonic diffusion
in these bilayers may arise from induced dipole interactions between the
bright excitons and the donor atoms across both monolayers, leading to a
qualitatively different behavior than the bright excitons in monolayers. We
note that the less bright excitonic species $\text{M}_{\text{D}}$ diffuses
more efficiently as bright excitons in the monolayer, which may arise from a
relatively longer lifetime43, 44 or perhaps due to an unusual interaction with
chiral phonons45.

Figure 4:  Diffusion lengths of different species in the monolayer and bilayer
in a spectral width of 1 meV across their respective peak PL intensities. Deep
blue lines trace the obtained diffusion lengths while light blue shaded
regions demarcate the 95 % confidence intervals from the fits.

To determine the nature of the more diffusive species, we trace the PL signal
as a function of the sample temperature. The evolution of diffusion lengths
with temperature is provided in the Supporting Information. Figure 5 (a)
exhibits the temperature dependence of the PL from the monolayer. We notice
the monotonic redshift with increasing temperature46 and the decrease in PL
yield. This is due to the presence of higher energy dark states in
$\text{MoSe}_{2}$. This decrease in quantum yield is opposite to that of
tungsten-based TMDC monolayers, where the presence of low-lying dark states
leads to an increase in PL yield with increasing temperature. We trace the PL
from the bilayer in Fig. 5 (b). We note evidence for a visible population
transfer to the now-brightened dark states at around $30\,\text{K}$, which
corresponds to a thermal energy of $\sim 2.5\,\text{meV}$, about half of the
difference in the peak energies of $\text{M}_{\text{X}}$ and
$\text{M}_{\text{D}}$. Around that temperature, we detect evidence of a
brightened dark trionic state $\text{M}_{\text{D}}^{-}$47 (see Supporting
Figure 5). The extra binding energy of the dark trion at $30\,\text{K}$ is
$24\,\text{meV}$ and is close to the binding energy of the bright trion
($27\,\text{meV}$) at the same temperature. The four excitonic and trionic
species in question are clearly identifiable in the spectra as four separate
peaks at $16.5\,\text{K}$ in Fig. 5(c). The dark species investigated in this
work are intravalley direct dark excitons and trions, and the PL emission is
not phonon-assisted, which can be surmised from the relative positions of
their peak energies from that of their bright counterparts. Finally, we report
the unusual, non-monotonic behavior of the energy of the dark exciton in Fig.
5(d). In contrast to the continuous redshift of the bright exciton with
temperature, the dark exciton initially undergoes a considerable blue shift in
its energy before it starts to redshift. This, too, points at a different band
origin of the electron in the dark exciton. The cause for this behavior may be
analogous to the anti-funneling effects observed for momentum-forbidden dark
excitons48 in single tungsten-based monolayers, which also arises from a
difference in how the electron bands evolve under strain, or in this case,
temperature.

Figure 5: Evolution of PL with temperature for the (a) monolayer and (b)
bilayer. (c) PL spectra at different temperatures exhibitng four separate
peaks. (d) Dependence of peak energies of bright exciton (yellow arrow) and
dark exciton (dark arrow) with temperature.

##  3 Conclusion

To summarize, we uncover the brightening of the spin-forbidden dark exciton
and dark trion in a high-angle incommensurate Moiré homobilayer. We identify a
more diffusive species by analyzing the spectral variation of diffusion
lengths in the PL spectrum, which we assign to the dark exciton. Investigating
the temperature dependence of the PL spectrum leads us to discover the
population transfer effects that are otherwise undetectable in the monolayer.
Hence, we show that despite the absence of a trapping Moiré superlattice,
these high-angle systems may host interesting and potentially useful physical
phenomena.

##  4 Methods

The first author prepared the manuscript, and all the co-authors participated
in discussing the draft. After the draft was near its final stage, ChatGPT 4o
was used to identify grammatical errors and consistency issues and to improve
the general readability of the text. The prompt and the response from ChatGPT
4o are provided in the supplementary information.

###  4.1 Fabrication

The monolayers and hBN (high-pressure anvil growth) were mechanically
exfoliated from high-quality bulk samples obtained from 2D Semiconductors. The
indivdual flakes were then assembled step-by-step under an optical microscope
using dome-shaped windows constructed from cured PDMS with a thin pane of PPC
(poly-propylene carbonate). After the device was constructed the assembly was
heated to release it when in contact with the DBR chip (with an additional 98
nm of $\text{SiO}_{2}$ on top). The SiN-terminated distributed Bragg reflector
was fabricated using a PECVD method, where 10.5 pairs of SiN/$\text{SiO}_{2}$.

###  4.2 Measurements

The measurements were carried out using a custom-made confocal microscope. A
532 nm DPSS laser is focused into a submicrometer diameter spot using a 0.70
NA objective lens in a closed-cycle cryostat (Montana Instruments) at 6 K. The
emission spot is relayed through the objective and imaged onto the CCD camera
using an achromatic lens system. One of the achromatic lenses is stepped
longitudinally to minimize the effects of longitudinal color in the system
between the measurements of the laser spot and the photoluminescence as they
differ in their wavelengths considerably. The collected PL is then analyzed
using a Princeton Instruments spectrometer (Acton SP-2750i) and an LN2 cooled
Pylon CCD camera. An Msquared Sprite XT femtosecond pulsed ti-saph (80 MHz
repetition rate) is used at 790.5 nm for the second harmonic generation
measurements.

##  5 Supporting Information

More data and analysis complementing the results presented in the main paper
are presented in the supporting information file.

##  6 Acknowledgments

This work was supported by AFOSR FA9550-19-1-0074 from the Cornell Center for
Materials Research with funding from the NSF MRSEC program (DMR-1719875). S.T
acknowledges NSF CMMI 2129412 and NSF DMR 2111812.

## References

  * Chen et al. 2022 Chen, D.; Lian, Z.; Huang, X.; Su, Y.; Rashetnia, M.; Ma, L.; Yan, L.; Blei, M.; Xiang, L.; Taniguchi, T.; et. al. Excitonic insulator in a heterojunction moiré superlattice. _Nature Physics_ 2022, _18_ , 1171–1176
  * Tang et al. 2020 Tang, Y.; Li, L.; Li, T.; Xu, Y.; Liu, S.; Barmak, K.; Watanabe, K.; Taniguchi, T.; MacDonald, A. H.; Shan, J.; Mak, K. F. Simulation of Hubbard model physics in WSe2/WS2 moiré superlattices. _Nature_ 2020, _579_ , 353–358
  * Gu et al. 2022 Gu, J.; Ma, L.; Liu, S.; Watanabe, K.; Taniguchi, T.; Hone, J. C.; Shan, J.; Mak, K. F. Dipolar excitonic insulator in a moiré lattice. _Nature Physics_ 2022, _18_ , 395–400
  * Zhao et al. 2023 Zhao, W.; Shen, B.; Tao, Z.; Han, Z.; Kang, K.; Watanabe, K.; Taniguchi, T.; Mak, K. F.; Shan, J. Gate-tunable heavy fermions in a moiré Kondo lattice. _Nature_ 2023, _616_ , 61–65
  * Miao et al. 2021 Miao, S.; Wang, T.; Huang, X.; Chen, D.; Lian, Z.; Wang, C.; Blei, M.; Taniguchi, T.; Watanabe, K.; Tongay, S.; Wang, Z.; Xiao, D.; Cui, Y.-T.; Shi, S.-F. Strong interaction between interlayer excitons and correlated electrons in WSe2/WS2 moiré superlattice. _Nature Communications_ 2021, _12_ , 3608
  * Zhang et al. 2022 Zhang, Z.; Regan, E. C.; Wang, D.; Zhao, W.; Wang, S.; Sayyad, M.; Yumigeta, K.; Watanabe, K.; Taniguchi, T.; Tongay, S.; Crommie, M.; Zettl, A.; Zaletel, M. P.; Wang, F. Correlated interlayer exciton insulator in heterostructures of monolayer WSe2 and moiré WS2/WSe2. _Nature Physics_ 2022, _18_ , 1214–1220
  * Xu et al. 2020 Xu, Y.; Liu, S.; Rhodes, D. A.; Watanabe, K.; Taniguchi, T.; Hone, J.; Elser, V.; Mak, K. F.; Shan, J. Correlated insulating states at fractional fillings of moiré superlattices. _Nature_ 2020, _587_ , 214–218
  * Li et al. 2021 Li, T.; Jiang, S.; Li, L.; Zhang, Y.; Kang, K.; Zhu, J.; Watanabe, K.; Taniguchi, T.; Chowdhury, D.; Fu, L.; Shan, J.; Mak, K. F. Continuous Mott transition in semiconductor moiré superlattices. _Nature_ 2021, _597_ , 350–354
  * Yu et al. 2017 Yu, H.; Liu, G.-B.; Tang, J.; Xu, X.; Yao, W. Moiré excitons: From programmable quantum emitter arrays to spin-orbit–coupled artificial lattices. _Science Advances_ 2017, _3_ , e1701696
  * Seyler et al. 2019 Seyler, K. L.; Rivera, P.; Yu, H.; Wilson, N. P.; Ray, E. L.; Mandrus, D. G.; Yan, J.; Yao, W.; Xu, X. Signatures of moiré-trapped valley excitons in MoSe2/WSe2 heterobilayers. _Nature_ 2019, _567_ , 66–70
  * Baek et al. 2020 Baek, H.; Brotons-Gisbert, M.; Koong, Z. X.; Campbell, A.; Rambach, M.; Watanabe, K.; Taniguchi, T.; Gerardot, B. D. Highly energy-tunable quantum light from moire-trapped excitons. _Science Advances_ 2020, _6_ , eaba8526
  * Ray et al. 2023 Ray, A. B.; Mukherjee, A.; Qiu, L.; Sailus, R.; Tongay, S.; Vamivakas, A. N. Interplay of Trapped Species and Absence of Electron Capture in Moiré Heterobilayers. _Nano Letters_ 2023, _23_ , 5989–5994
  * Soltero et al. 2024 Soltero, I.; Kaliteevski, M. A.; McHugh, J. G.; Enaldiev, V.; Fal’ko, V. I. Competition of Moiré Network Sites to Form Electronic Quantum Dots in Reconstructed MoX2/WX2 Heterostructures. _Nano Letters_ 2024, _24_ , 1996–2002
  * Liu et al. 2021 Liu, E.; Barré, E.; van Baren, J.; Wilson, M.; Taniguchi, T.; Watanabe, K.; Cui, Y.-T.; Gabor, N. M.; Heinz, T. F.; Chang, Y.-C.; Lui, C. H. Signatures of moiré trions in WSe2/MoSe2 heterobilayers. _Nature_ 2021, _594_ , 46–50
  * Wang et al. 2021 Wang, X.; Zhu, J.; Seyler, K. L.; Rivera, P.; Zheng, H.; Wang, Y.; He, M.; Taniguchi, T.; Watanabe, K.; Yan, J.; Mandrus, D. G.; Gamelin, D. R.; Yao, W.; Xu, X. Moiré trions in MoSe2/WSe2 heterobilayers. _Nature Nanotechnology_ 2021, _16_ , 1208–1213
  * Li et al. 2024 Li, Y. et al. Tuning commensurability in twisted van der Waals bilayers. _Nature_ 2024, _625_ , 494–499
  * Lopes dos Santos et al. 2012 Lopes dos Santos, J. M. B.; Peres, N. M. R.; Castro Neto, A. H. Continuum model of the twisted graphene bilayer. _Physical Review B_ 2012, _86_ , 155449
  * Shallcross et al. 2010 Shallcross, S.; Sharma, S.; Kandelaki, E.; Pankratov, O. A. Electronic structure of turbostratic graphene. _Phys. Rev. B_ 2010, _81_ , 165105
  * Lau et al. 2022 Lau, C. N.; Bockrath, M. W.; Mak, K. F.; Zhang, F. Reproducibility in the fabrication and physics of moiré materials. _Nature_ 2022, _602_ , 41–50
  * Unuchek et al. 2018 Unuchek, D.; Ciarrocchi, A.; Avsar, A.; Watanabe, K.; Taniguchi, T.; Kis, A. Room-temperature electrical control of exciton flux in a van der Waals heterostructure. _Nature_ 2018, _560_ , 340–344
  * Mak and Shan 2022 Mak, K. F.; Shan, J. Semiconductor moiré materials. _Nature Nanotechnology_ 2022, _17_ , 686–695
  * Villafañe et al. 2023 Villafañe, V.; Kremser, M.; Hübner, R.; Petrić, M. M.; Wilson, N. P.; Stier, A. V.; Müller, K.; Florian, M.; Steinhoff, A.; Finley, J. J. Twist-Dependent Intra- and Interlayer Excitons in Moiré ${\mathrm{MoSe}}_{2}$ Homobilayers. _Phys. Rev. Lett._ 2023, _130_ , 026901
  * Cadiz et al. 2017 Cadiz, F. et al. Excitonic Linewidth Approaching the Homogeneous Limit in ${\mathrm{MoS}}_{2}$-Based van der Waals Heterostructures. _Phys. Rev. X_ 2017, _7_ , 021026
  * Splendiani et al. 2010 Splendiani, A.; Sun, L.; Zhang, Y.; Li, T.; Kim, J.; Chim, C.-Y.; Galli, G.; Wang, F. Emerging Photoluminescence in Monolayer MoS2. _Nano Letters_ 2010, _10_ , 1271–1275
  * Malic et al. 2018 Malic, E.; Selig, M.; Feierabend, M.; Brem, S.; Christiansen, D.; Wendler, F.; Knorr, A.; Berghäuser, G. Dark excitons in transition metal dichalcogenides. _Physical Review Materials_ 2018, _2_ , 014002
  * Brem et al. 2020 Brem, S.; Ekman, A.; Christiansen, D.; Katsch, F.; Selig, M.; Robert, C.; Marie, X.; Urbaszek, B.; Knorr, A.; Malic, E. Phonon-Assisted Photoluminescence from Indirect Excitons in Monolayers of Transition-Metal Dichalcogenides. _Nano Letters_ 2020, _20_ , 2849–2856
  * Zhang et al. 2015 Zhang, X.-X.; You, Y.; Zhao, S. Y. F.; Heinz, T. F. Experimental Evidence for Dark Excitons in Monolayer ${\mathrm{WSe}}_{2}$. _Phys. Rev. Lett._ 2015, _115_ , 257403
  * Slobodeniuk and Basko 2016 Slobodeniuk, A. O.; Basko, D. M. Spin–flip processes and radiative decay of dark intravalley excitons in transition metal dichalcogenide monolayers. _2D Materials_ 2016, _3_ , 035009
  * Yu et al. 2018 Yu, H.; Liu, G.-B.; Yao, W. Brightened spin-triplet interlayer excitons and optical selection rules in van der Waals heterobilayers. _2D Materials_ 2018, _5_ , 035021
  * Zhang et al. 2019 Zhang, L.; Gogna, R.; Burg, G. W.; Horng, J.; Paik, E.; Chou, Y.-H.; Kim, K.; Tutuc, E.; Deng, H. Highly valley-polarized singlet and triplet interlayer excitons in van der Waals heterostructure. _Phys. Rev. B_ 2019, _100_ , 041402
  * Anankine et al. 2017 Anankine, R.; Beian, M.; Dang, S.; Alloing, M.; Cambril, E.; Merghem, K.; Carbonell, C. G.; Lemaître, A.; Dubin, F. m. c. Quantized Vortices and Four-Component Superfluidity of Semiconductor Excitons. _Phys. Rev. Lett._ 2017, _118_ , 127402
  * Combescot et al. 2007 Combescot, M.; Betbeder-Matibet, O.; Combescot, R. Bose-Einstein Condensation in Semiconductors: The Key Role of Dark Excitons. _Phys. Rev. Lett._ 2007, _99_ , 176403
  * Mennel et al. 2018 Mennel, L.; Paur, M.; Mueller, T. Second harmonic generation in strained transition metal dichalcogenide monolayers: MoS2, MoSe2, WS2, and WSe2. _APL Photonics_ 2018, _4_ , 034404
  * Mennel et al. 2018 Mennel, L.; Furchi, M. M.; Wachter, S.; Paur, M.; Polyushkin, D. K.; Mueller, T. Optical imaging of strain in two-dimensional crystals. _Nature Communications_ 2018, _9_ , 516
  * Choi et al. 2020 Choi, J. et al. Moiré potential impedes interlayer exciton diffusion in van der Waals heterostructures. _Science Advances_ 2020, _6_ , eaba8866
  * Brem et al. 2020 Brem, S.; Linderälv, C.; Erhart, P.; Malic, E. Tunable Phases of Moiré Excitons in van der Waals Heterostructures. _Nano Letters_ 2020, _20_ , 8534–8540
  * Li et al. 2021 Li, Z.; Lu, X.; Cordovilla Leon, D. F.; Lyu, Z.; Xie, H.; Hou, J.; Lu, Y.; Guo, X.; Kaczmarek, A.; Taniguchi, T.; Watanabe, K.; Zhao, L.; Yang, L.; Deotare, P. B. Interlayer Exciton Transport in MoSe2/WSe2 Heterostructures. _ACS Nano_ 2021, _15_ , 1539–1547
  * Rapaport et al. 2006 Rapaport, R.; Chen, G.; Simon, S. H. Nonlinear dynamics of a dense two-dimensional dipolar exciton gas. _Phys. Rev. B_ 2006, _73_ , 033319
  * Ivanov 2002 Ivanov, A. L. Quantum diffusion of dipole-oriented indirect excitons in coupled quantum wells. _Europhysics Letters_ 2002, _59_ , 586
  * Cheng et al. 2021 Cheng, G.; Li, B.; Jin, Z.; Zhang, M.; Wang, J. Observation of Diffusion and Drift of the Negative Trions in Monolayer WS2. _Nano Letters_ 2021, _21_ , 6314–6320
  * Perea-Causin et al. 2022 Perea-Causin, R.; Brem, S.; Malic, E. Trion-phonon interaction in atomically thin semiconductors. _Phys. Rev. B_ 2022, _106_ , 115407
  * Uddin et al. 2020 Uddin, S. Z.; Kim, H.; Lorenzon, M.; Yeh, M.; Lien, D.-H.; Barnard, E. S.; Htoon, H.; Weber-Bargioni, A.; Javey, A. Neutral Exciton Diffusion in Monolayer MoS2. _ACS Nano_ 2020, _14_ , 13433–13440
  * Robert et al. 2017 Robert, C.; Amand, T.; Cadiz, F.; Lagarde, D.; Courtade, E.; Manca, M.; Taniguchi, T.; Watanabe, K.; Urbaszek, B.; Marie, X. Fine structure and lifetime of dark excitons in transition metal dichalcogenide monolayers. _Phys. Rev. B_ 2017, _96_ , 155423
  * Slobodeniuk and Basko 2016 Slobodeniuk, A. O.; Basko, D. M. Spin–flip processes and radiative decay of dark intravalley excitons in transition metal dichalcogenide monolayers. _2D Materials_ 2016, _3_ , 035009
  * Zhu et al. 2018 Zhu, H.; Yi, J.; Li, M.-Y.; Xiao, J.; Zhang, L.; Yang, C.-W.; Kaindl, R. A.; Li, L.-J.; Wang, Y.; Zhang, X. Observation of chiral phonons. _Science_ 2018, _359_ , 579–582
  * Arora et al. 2015 Arora, A.; Nogajewski, K.; Molas, M.; Koperski, M.; Potemski, M. Exciton band structure in layered MoSe2: from a monolayer to the bulk limit. _Nanoscale_ 2015, _7_ , 20769–20775
  * Arora et al. 2020 Arora, A.; Wessling, N. K.; Deilmann, T.; Reichenauer, T.; Steeger, P.; Kossacki, P.; Potemski, M.; Michaelis de Vasconcellos, S.; Rohlfing, M.; Bratschitsch, R. Dark trions govern the temperature-dependent optical absorption and emission of doped atomically thin semiconductors. _Phys. Rev. B_ 2020, _101_ , 241413
  * Rosati et al. 2021 Rosati, R.; Schmidt, R.; Brem, S.; Perea-Causín, R.; Niehues, I.; Kern, J.; Preuß, J. A.; Schneider, R.; Michaelis de Vasconcellos, S.; Bratschitsch, R.; Malic, E. Dark exciton anti-funneling in atomically thin semiconductors. _Nature Communications_ 2021, _12_ , 7221