# Exact Results for Interacting Hard Rigid Rotors on a d-Dimensional Lattice

Sushant Saryal<EMAIL_ADDRESS>Department of
Physics,Indian Institute of Science Education and Research, Pune 411008, India
Deepak Dhar<EMAIL_ADDRESS>Department of Physics,Indian Institute of
Science Education and Research, Pune 411008, India

###### Abstract

We study the entropy of a set of identical hard objects, of general shape,
with each object pivoted on the vertices of a d-dimensional regular lattice of
lattice spacing $a$, but can have arbitrary orientations. When the pivoting
point is situated asymmetrically on the object, we show that there is a range
of lattice spacings $a$, where in any orientation, a particle can overlap with
at most one of its neighbors. In this range, the entropy of the system of
particles can be expressed exactly in terms of the grand partition function of
coverings of the base lattice by dimers at a finite negative activity. The
exact entropy in this range is fully determined by the second virial
coefficient. Calculation of the partition function is also shown to be
reducible to that of the same model with discretized orientations. We
determine the exact functional form of the probability distribution function
of orientations at a site. This depends on the density of dimers for the given
activity in the dimer problem, that we determine by summing the corresponding
Mayer series numerically. These results are verified by numerical simulations.

Many organic liquids, on cooling, first freeze into a state where there is
long range crystalline order, but the molecules are partly able to rotate
about their center of mass. In contrast to the liquid crystals that shows
orientational order, but no positional order, these materials show no long
range orientational order, and have been called plastic solids, rotator
crystals or orientationally disordered solids. Timmermans studied it in the
1930’s timmermans , and early work has been summarized in staveley . These are
currently attracting a lot of interest, for applications such as storage
batteries batteries , drug delivery drug-delivery , and refrigeration
refrigeration .

The simplest theoretical description of this state would be to assume that
rigid molecules have their centers of mass fixed at lattice positions, but
have orientational degrees of freedom. Many such models have been studied
since the beginning of these studies, e.g. pople-karasz . Casey and Runnels
Runnels1  studied a system of hard squares with centers fixed on the one
dimensional lattice and also with centers free to move along a 1D line.
Freasier and RunnelsRunnels2  studied a model of rotors with their centers
pivoted on a one dimensional lattice. More recently, Kantor and Kardar Kardar
, and Gurin and Varga Gurin  studied a system of hard rods in isobaric
ensemble with their centers restricted to move only along a 1D line but free
to rotate in a plane.
In an earlier paper, we have studied hard linear rods, with centers fixed on
sites of a 1-d lattice, but free to rotate in a plane saryal1 . We showed that
even in one dimension, this system shows an infinite number of singularities
of entropy per rod as a function of the ratio of length of rods and the
lattice spacing. Plastic solids often show many (three or more) crystalline
phases with varying amount of orientational order guthrie . It is gratifying
that the simple model of rigid hard particles can describe this basic
phenomenology.

Figure 1: An example of hard particles with their pivots fixed at sites of a
square lattice of lattice spacing $a$. The farthest point in the object from
the pivot is at distance $\ell$. The orientation of the rotor at site ${\bf
i}$ is specified by the orientation of the long axis to the x-axis
$\theta_{\bf i}$.

There are few exact results available for this problem in more than one
dimension. In this work, we show that in this model of rigid rotors on a
lattice [an example in two dimensions is shown in Fig. 1], there is a range of
lattice spacings in which the model becomes reducible to the dimer model on
the same lattice with a negative value of the dimer activity. We determine the
exact functional form of the probability distribution function of orientations
at a site. This depends on the density of dimers for the given activity in the
dimer problem, that we determine by summing the corresponding Mayer series
numerically. These results are valid for general shapes of rotors, and for
arbitrary lattices, and dimensions. We have verified these results by
numerical simulations for thin hard linear rods pivoted at one end at the
vertices of a square lattice [Fig. 2].
One important feature of the model is the fact that interaction between two
rotors has a different form in different directions. In some cases, e.g. in
the class of models known as compass models compass , this feature leads to
the model becoming analytically tractable. The best known example of this
class is the Kitaev model kitaev . In the standard version of the Kitaev
model, the degrees of freedom are quantum spins, though the classical case has
also been studied classical-kitaev . We will only study the case of classical
rotors here.

We consider a system of hard bodies, called rotors, of identical shape and
size in $d$ dimensions. Each rotor has a marked point called the pivot, which
is at the same relative position in each rotor. The objects are rigid, of a
given shape. Our treatment is valid for arbitrary shapes. The pivots of the
particles are fixed at the vertices of a regular $d$-dimension lattice ( e.g.
the $d$-dimensional hypercubical lattice), but the rotors can take arbitrary
orientations, subject to the constraint that different rotors cannot overlap.
The orientation of the rotor at the vertex $\bf{r}$ of the lattice may be
specified by a $d\times d$ real orthogonal matrix $M(\bf{r})$. Then, the full
configuration is specified by a lattice model, with a $d\times d$ matrix
belonging to the group $SO(d)$ at each vertex. [Sometimes, the experimentally
prepared crystals are racemic mixtures of the two species of molecules related
by mirror symmetry. These crystals would be described by the group O(d). Our
treatment is easily extended to this case.]

Figure 2: Hard rods pivoted at the one end on a square lattice. Length of each
rod is $l$ , lattice spacing is $a$. The orientation of rotor at site ${\bf
i}$ is given by the angle its long axis makes with the x-axis.

We first consider the case when the pivot is placed asymmetrically on the
rotor, and there is a unique point on the rotor at maximum distance from the
pivot. Let this distance be denoted by $\ell$. Let the distance between
nearest neighbors vertices of the lattice be $a$.The line joining the pivot to
the farthest point will be called the long axis of the rotor. Then, we will
study the entropy per site for this model, as a function of the dimensionless
coupling parameter $\ell/a$, for a specified shape of the object.

We define the indicator function $\eta(\bf{a},\bf{b})$ to be $1$ if the rotors
at sites $\bf{a}$ and $\bf{b}$, with given orientations $M(\bf{a})$ and
$M(\bf{b})$, overlap, and zero otherwise. For notational simplicity, we will
not show the dependence on orientations explicitly. Then the partition
function for this system is defined as

$\mathcal{Z}_{N}=\left[\prod_{\bf{r}}\int
dM(\bf{r})\right]\prod_{<\bf{a},\bf{b}>}\left[1-\eta(\bf{a},\bf{b})\right],$
(1)

where the second product is over all nearest neighbor pairs, and the integral
over the matrix $M(\bf{r})$ is the normalized Haar measure over the group
SO(d), so that $\int dM({\bf r})~{}1=1$. We will only consider the case where
rotors at sites farther than nearest neighbors cannot overlap, and second
product is restricted to only over nearest neighbor pairs. The entropy per
site $s(\ell/a)$ is defined by

$s(\ell/a)=\lim_{N\rightarrow\infty}[\log{Z}_{N}]/N,$ (2)

where $N$ is the number of sites in the lattice. If $\ell/a<1/2$, then
different rotors cannot overlap, and the entropy per rotor takes its maximum
value, zero in our normalization. As the spacing $a$ is decreased, keeping
$\ell$ constant, the entropy will also decrease. When $a$ is decreased to the
minimum allowed value, corresponding to closest packing, entropy per rotor
will tend to $-\infty$.

Now consider the case $\ell/a=1/2+\epsilon$, where $\epsilon$ is a small
positive number. Consider two rotors at nearest neighbor vertices $\bf{i}$ and
$\bf{j}$ along the bond along one of the lattice directions, say
${\bf\hat{e}}_{1}$. Then, these two rotors can overlap only if the long axes
of both are substantially aligned parallel to the direction
${\bf\hat{e}}_{1}$. We define $\omega_{1}$ as set of orientations of the rotor
at $\bf{i}$ such that there is a non-empty set of orientations of the rotor at
$\bf{j}$ that overlap with the rotor at $\bf{i}$. The volume of the set
$\omega_{1}$ decreases to zero, as $\epsilon$ tends to zero.

We similarly construct sets $\omega_{-1},\omega_{\pm 2},...\omega_{\pm d}$,
using other neighbors of $\bf{i}$ (set $\omega_{-\alpha}$ corresponds to the
$-{\bf\hat{e}}_{\alpha}$ direction). The total number of such sets is the
number of neighbors of $\bf{i}$, i.e. $2d$ on the $d$-dimensional hypercubical
lattice. For small $\epsilon$, these sets $\omega_{i}$ are disjoint. As $a$ is
decreased, the sets $\omega_{i}$’s increase in size, and eventually, some will
touch. Let $a^{*}$ be the smallest value of $a$ such that the sets
$\omega_{i}$ are still disjoint. Then for all the range of values
$2l>a>a^{*}$, a rotor at any site $\bf{i}$, in any orientation $M(\bf{i})$ can
have an overlap with at most one of its neighbors, whatever their
orientations. Let ${\bf j}$ and ${\bf j^{\prime}}$ be two neighbors of ${\bf
i}$, with ${\bf j}\neq{\bf j^{\prime}}$. Then we have, for all ${\bf i},{\bf
j},{\bf j}^{\prime}$, and all orientations $M({\bf i}),M({{\bf
j}}),M(\bf{j}^{\prime})$,

$\eta({\bf i},{\bf j})\eta({\bf i},{\bf j}^{\prime})=0.$ (3)

We will call this the at-most-one-overlap ( AOO) condition.

A simple example of systems where AOO condition holds in some region of the
parameter space is hard rigid lines of length $\ell$, where one of the ends is
the pivoting point as shown in Fig. 2. In this case, it is straightforward to
verify that $a^{*}=\sqrt{2}\ell$, and determine explicitly the function
$\eta(\theta,\theta^{\prime})$. The details are given in Appendix A.

Now, we expand the product in Eq. (1), and make a graphical representation of
the terms of the expansion, with a bond between vertices $\bf{i}$ and
$\bf{j}$, iff the term contains $\eta(\bf{i},\bf{j})$, see Fig. 3. This is a
lattice version of the well-known Mayer cluster expansion. Then, using the AOO
condition, only graphs with at most one bond coming to a site survive. Thus,
we get

$\mathcal{Z}_{N}=\mathcal{Z}_{dimers}(z),$ (4)

where $\mathcal{Z}_{dimers}(z)$ is the grand partition function of partial
covering of the vertices of the lattice by dimers (see Fig. 3),

$\mathcal{Z}_{dimers}(z)=\sum_{\textrm{dimer-confgs}}z^{\textrm{number of
dimers}},$ (5)

with the fugacity of a dimer is the (negative ) real number, given by

$z=-\int dM({\bf i})\int dM({\bf j})~{}\eta({\bf i},{\bf j}).$ (6)

In fact, $z$ is exactly the second virial coefficient $B_{2}$ of the virial
expansion. It is straight forward to determine it for a a given shape of
molecules. In Fig. 4, we have plotted $z$ as a function of $\ell/a$ for the
case of linear rods, pivoted at one end shown in Fig. 2. Let $g(z)$ be the
logarithm of the grand partition function per site, defined as

$\displaystyle
g(z)=\lim_{N\rightarrow\infty}[\log{\mathcal{Z}}_{dimers}(z)]/N$ (7)

Then, the entropy per site is $s(\ell/a)=g(z)$.

Figure 3: A typical term in the graphical expansion of the partition function
[Eq.(1)]. This can also be thought of as a configuration of a set of dimers on
the lattice.

The evaluation of the partition function of partial covering of dimers is a
well-known hard problem in lattice statistics, and no exact expression for
$g(z)$ is known, except in one dimension and in some special graphs, like the
Bethe lattice. However, the equivalence is not without use, as we see that the
exact entropy $s(\ell/a)$, in the entire range of $a$ satisfying the AOO
condition, depends on shape of the rotors, only through the second virial
coefficient $z$.

Also, we note that the equality of the dimer and rotor model partition
functions holds also for finite lattices of size $L_{1}\times L_{2}\times
L_{3}...$. The dependence of these partition functions on one of the
dimensions is of the form

$\mathcal{Z}\sim\sum_{\alpha}c_{\alpha}\lambda_{\alpha}^{L_{1}},$ (8)

where $\lambda_{\alpha}$ are the eigenvalues of the transfer matrix in the
direction $1$, and $c_{\alpha}$ are some constants. Since this is true for all
$L_{1}$, we conclude that all the nonzero eigenvalues of the transfer matrix
for the rotor model are the same as eigenvalues of the transfer matrix for the
dimer model with activity $z$. In Appendix B and C, we show an explicit
calculation of the model for hard rods on a line and zigzag lattice
respectively, using transfer matrices, and show that the results agree with
the algebraic method used here.

Figure 4: $-z$ vs $l/a$, for a system of hard linear rods of length $l$
pivoted at one end on 2d square lattice.

We will now show that under the AOO condition, our model with continuous
$SO(d)$ degrees of freedom is reducible to a discrete model with $2d$ discrete
degrees of freedom per site. Following the terminology of ice and KDP models
ice , we will call this model the Hydrogen Iodide (HI) model, where the heavy
“iodine”(I) atoms sit at lattice sites, and there is one “hydrogen”(H) atom
near each I atom that sits along one of the nearest neighbor bonds, close to
it. There is a repulsive energy cost if the H atoms of two adjacent I’s sit on
the same bond. Thus, there is a discrete degree of freedom at each site ${\bf
r}$ of the lattice, denoted by $\sigma(\bf r)$, which takes non-zero integer
values between $-d$ and $d$. If the hydrogen atom is along the bond
${\bf\hat{e}}_{\alpha}$, we will call the state $\alpha$, if along
$-{\bf\hat{e}}_{\alpha}$, we call the state $-\alpha$ ($1\leq\alpha\leq d$) .
Therefore, for a bond between the sites ${\bf i}$ and ${\bf
i}+{\bf\hat{e}}_{\alpha}$, the energy cost is $J$ iff the spin at site ${\bf
i}$ is in the state $\alpha$, and that at site $({\bf
i}+{\bf\hat{e}}_{\alpha})$ is in the state $-\alpha$. The interaction along
the bond in the direction ${\bf\hat{e}}_{\alpha}$ only involves the spin
states $\pm\alpha$.

The hamiltonian of the model is

$\displaystyle{\mathcal{H}}_{HI}=J\sum_{{\bf
r}}\sum_{\alpha=1}^{d}\delta(\sigma({\bf r}),\alpha)\delta(\sigma({\bf
r}+{\bf\hat{e}}_{\alpha}),-\alpha),$ (9)

where $\delta(a,b)$ is the Kronecker delta function. We write
$\exp\left[-\beta J\delta(a,b)\right]=1+\tilde{z}\delta(a,b)$, with
$\tilde{z}=e^{-\beta J}-1$, $\beta$ is the inverse temperature, and expand the
partition function in powers of $\tilde{z}$, we again get the partition
functions of dimers with activity $z=\tilde{z}/(2d)^{2}$. Let $Z_{HI}$ denote
the partition function of the HI model, then we have,

$\displaystyle
Z_{HI}=(2d)^{N}{\mathcal{Z}}_{dimers}\Big{(}z=\frac{\tilde{z}}{(2d)^{2}}\Big{)}$
(10)

To complete the correspondence with the rotors models, we break the manifold
of states given by $M$ into $2d$ equal submanifolds. A given matrix $M$ is
associated with the discrete state $\alpha$ of the HI model as follows: we
imagine starting with a large value of the lattice spacing $a$, and decreasing
it slowly, and monitor the overlap integral with neighbour in the direction
$\alpha$, for fixed $M$. At some value of spacings, one of the integrals
becomes non zero for the first time. That is the direction of the HI model
spin associated to this $M$. Values on $M$ which lie at the common boundaries
of two submanifolds, and where this direction is not unique, have zero
measure. For all smaller $a$, this overlap will increase, but by definition,
all other overlaps will remain strictly zero, so long as the AOO condition is
satisfied. As a simple example, in two dimensions, the manifold $M$ is the
full circle with $0\leq\theta\leq 2\pi$, and it is broken into four
submanifolds $\big{(}[-\pi/4,\pi/4),[\pi/4,3\pi/4),[3\pi/4,5\pi/4)$ and
$[5\pi/4,7\pi/4)\big{)}$ which correspond to the four states of the discrete
model.

The graphical expansion of the partition function, for a given dimer
configuration, the sites covered by dimers have a fixed $\sigma$-variable, and
we sum over all orientations of arrows at sites where there is no dimer. The
full partition function has the form of the inclusion-exclusion principle. The
empty dimer configuration corresponds to all possible states of the discrete
spins. The term with one dimer is negative, and subtracts the sum of
contributions where specified one bond is doubly occupied, and others are
aribitrary, summed over all possible positions of the bond. The next term
corrects for the overcounting by adding the contribution where two bonds are
doubly occupied, and others are arbitrary, and so on. Then the integration
over submanifolds for specified discrete state at each site is easily done. In
the integration over the $M$-variables, restricted to specific submanifolds,
sites not covered by same dimer are independent. This implies that all the
$n$-point correlation functions of rotor model with AOO condition can be
expressed as linear functions of the $n$-point correlations of
${\mathcal{H}}_{HI}$. We illustrate this below using the 1-point function.

Figure 5: $f_{1}(\theta)$ for a system of hard linear rods of length 2d square
lattice for $a=1.6l$ . Because of the symmetry of the square lattice other
$f_{\alpha}$ are related to $f_{1}$ such that
$f_{-1}(\theta)=f_{1}({\theta-\pi})$ and $f_{\pm
2}(\theta)=f_{1}({\theta\pm\pi/2})$.

Let $f_{\alpha}(M)dM$ be the conditional probability that variable $M$ of the
rotor at ${\bf O}$ lies in a small volume $dM$ near the value $M$, given that
it overlaps with the neighboring rotor in the direction $\alpha$.
Then we have

$f_{\alpha}(M)=\dfrac{\int dM({\bf\hat{e}}_{\alpha})~{}\ \eta({\bf
O},{\bf\hat{e}}_{\alpha})}{\int dM({\bf O})\int dM({\bf\hat{e}}_{\alpha})~{}\
\eta({\bf O},{\bf\hat{e}}_{\alpha})},$ (11)

In Fig. 5, we have shown the plot of $f_{1}(\theta)$ for the linear rods
pivoted at the end. Clearly $f_{\alpha}$ is zero outside the the submanifold
$\omega_{\alpha}$.
Let ${\rm P}(M)dM$ denote the probability that the matrix at the origin ${\bf
O}$, in the thermodynamical limit, lies in the volume $dM$ centered at $M$. In
the partition sum in Eq.(1), we sum over all sites, except the origin. If
$\bar{n}(z)$ is the number density of dimers per site, then with probability
weight $(1-2\bar{n}(z))$, the site is not covered by a dimer, and with weight
$\dfrac{\bar{n}(z)}{d}$, it is covered by a dimer in direction $\alpha$. Then
${\rm P}(M)$ is given by,

${\rm
P}(M)=\left[(1-2\bar{n})+\frac{\bar{n}}{d}\sum_{\alpha}f_{\alpha}(M)\right].$
(12)

As a concrete example, we calculate ${\rm P}(M)$ for thin hard linear rods of
length $\ell$ pivoted at one end at the vertices of a square lattice [Fig. 2].
Here $M$ is just a simple rotation angle $\theta\in[-\pi,\pi)$. In this case,
it is straight forward to calculate the function $f_{\alpha}(\theta)$, for
$\alpha=\pm 1,\pm 2$, numerically for any given value of $\ell/a$ using
Eq.(11).
In Fig. 6, we have plotted the function ${\rm P}(\theta)$, as determined by
Monte Carlo simulations, against the theoretical form, for different values of
$\ell/a$. We do not have the functional form of $\bar{n}(z)$, as the solution
to the dimer model is not known. However, the experimental data can be fitted
to the functional form in Eq.(12), treating the unknown value $\bar{n}(z)$ as
a single fitting parameter , for a given value of $\ell/a$.This is done in
Fig. 6, and we see that using only a single parameter, we are able to get a
very good fit to the function ${\rm P}(\theta)$ for the entire range of
$\theta$, and for different values of $\ell/a$. In Appendix B and C, we verify
that the distribution of orientations also agrees with that calculated using
the eigenvectors of the transfer matrix for the 1-dimensional chain, and the
zig-zag chain.

In Fig. 7, we have plotted the fitting value of $\bar{n}(z)$ as a function of
the overlap integral $z$ for different values of $\ell/a$. As already
mentioned, exact functional form of $\bar{n}(z)$ is not known. But, for the
low values of activity in the range where AOO holds, it is adequate to use
first few terms in the low density Mayer series for the dimer density as a
function of the activity $z$ :

$\bar{n}(z)=\sum_{n=1}^{\infty}(-1)^{n-1}a_{n}z^{n}$ (13)

Figure 6: The probability distribution of angles ${\rm P}(\theta)$ as a
function of the orientation angle $\theta$, for a system of hard rods pivoted
at the one end on a square lattice (see Fig. 2), for different values of
lattice spacing. The function outside this range can be obtained using the
symmetries of the square lattice.

It was shown by Viennot that the Mayer coefficients $a_{n}$ are integers equal
to the number of heaps one can form on the simple cubic lattice using dimers
vien . For small values of $z$ of interest here, only the first few terms of
this series are enough to determine $\bar{n}(z)$ to three digit accuracy,
which is given by

$\bar{n}(z)=2\left[z-7z^{2}+58z^{3}-521z^{4}....\right]$ (14)

Figure 7: Plot of negative of average density $\bar{n}(z)$ per site obtained
from Monte Carlo simulations and theoretical value given by eq.(14).

At $\ell/a=1/\sqrt{2}$, the value of activity is approximately $z=-0.014$,
which is much less than the radius of convergence of the series in Eq.(14).
Even though the value of $z$ is rather small, the effect on $P(\theta)$ is
quite substantial as is seen in Fig. 6. From the listed few terms above, we
can estimate that the radius of convergence is near $0.1$ ( a much more
precise estimate is clearly possible, but this is adequate for our purpose
here).
This can be understood as follows: We need the AOO condition to establish the
equivalence between the rotors model and the dimer model. The breakdown of the
AOO condition is due to geometrical reasons specific to the rotors model. The
dimer model partition function shows no singularity as a function of activity,
if the lattice spacing is decreased so that the AOO condition is no longer
satisfied.

Of course, the rotors model is well defined, even outside the range of lattice
parameters where the AOO condition is satisfied, the minimum lattice spacing
allowed being set by density of closest packing. The mechanism leading to
geometrical transitions discussed in saryal1  is quite robust, and holds in
all dimensions. If we decrease $a$ below $a^{*}$, we would expect to get a
sequence of geometrical transitions. In addition, there can be other order-
disorder type of transitions.
As $a$ is decreased below $a^{*}=\sqrt{2}l$, the function z(a) is a smooth
function of $a$, but there are configurations there one dimer overlaps with
two others. These give rise to non-zero positive correction to the entropy
$s(\ell/a)$ . If $a=a^{*}(1-\epsilon)$, this correction is of order
$\epsilon^{3/2}$ (see Appendix D).

Figure 8: A typical term in the graphical expansion of the partition function
of rotors, when the objects are reflection symmetric, and the AOO condition is
not satisfied. Here, the terms can be thoughts of as configurations in the
grand partition function of a system of polydisperse straights rods on the
lattice, where the activity of a rod depends on its length.

Now, we discuss the case where the pivot point is a center of inversion of the
rotor. For example, the two-dimensional case, when rotors are elliptical, with
pivot at the center of the ellipse. In this case, AOO condition can not be
satisfied on the square lattice, but there is a range of values of $a$, such
that an ellipse can touch two other ellipses, only if the centers are
collinear. Then, the graphical representation, terms that contribute can be
thought of as configurations of non-overlapping straight needles, whose length
can be any integer. The weight of a needle depends on its length ( see Fig.
8). Thus, in this case, the partition function becomes equivalent to that of a
polydisperse long rods on the same lattice. The latter problem has been
studied numerically, but remains quite intractable analytically.
As a final remark we want to point out that this analysis is also valid for
soft rods. In case of attractive interaction among rods $z$ will be positive.
In summary, we have shown that there is a range of lattice spacings, where the
calculation of partition function of the rotor model becomes exactly
equivalent to the problem of evaluating the dimer model partition function on
the same lattice, analytically continued to negative activities of dimers.
This equivalence allows us to determine the distribution of orientations for
the original model. The results agree with direct Monte Carlo simulation of
the rotor model. The results are valid for arbitrary shapes of molecules, and
all dimensions.
It is hoped that insights from the rather idealized model will help generate a
better understanding of the actual plastic crystals.

Acknowledgments: We thank M. Barma, K. Damle, J. Klamser, S. N. Majumdar, R.
Rajesh, T. Sadhu, and G. J. Sreejith for their comments on the paper. DD’s
work is supported by a Senior Scientist fellowship from the National Academy
of Sciences of India. S.S. acknowledges financial support from Council of
Scientific and Industrial Research (Grant No. 1061651988).

## Appendix A: Calculation of $\eta(\theta,\theta^{{}^{\prime}})$ for thin
hard linear rods pivoted at one end

Consider two thin hard rods of length $l$ placed at adjacent sites, see Fig.
9. $\eta(\theta,\theta^{{}^{\prime}})$ will be 1 when they overlap and 0
otherwise, see Fig. 10. For notational convenience we define $\kappa=l/a$.
Then, if $\kappa<1/2$, each rod can rotate fully freely, and
$\eta(\theta,\theta^{{}^{\prime}})=0$ for all $\theta,\theta^{{}^{\prime}}$.

Figure 9:  Two thin hard rods pivoted at adjacent lattice sites.

While AOO condition in 2d square lattice is valid for $1/2\leq\kappa\leq
1/\sqrt{2}$, for 1d case, which we will discuss in next section, it is valid
for $1/2\leq\kappa\leq 1$. Therefore, we will discuss here the full range
$1/2\leq\kappa\leq 1$.
Because of the symmetry of the problem we have
$\eta(\theta,\theta^{{}^{\prime}})=\eta(-\theta,-\theta^{{}^{\prime}})$, we
will only consider $\theta>0$ to determine $\eta(\theta,\theta^{\prime})$.

For $1/2\leq\kappa\leq 1/\sqrt{2}$ , there is no overlap if
$|\theta|>\cos^{-1}(\frac{1}{2\kappa})$ . If
$0<\theta<\cos^{-1}(\frac{1}{2\kappa})$ rods will overlap only if
$\theta^{{}^{\prime}}\in[\theta^{{}^{\prime}}_{min},\theta^{{}^{\prime}}_{max}]$
where

$\displaystyle\theta^{{}^{\prime}}_{max}=\pi+\theta-\sin^{-1}\Big{[}\frac{\sin{\theta}}{\kappa}\Big{]}$
(A1)
$\displaystyle\theta^{{}^{\prime}}_{min}=\pi-\tan^{-1}\Big{[}\frac{\kappa\sin{\theta}}{1-\kappa\cos{\theta}}\Big{]}$
(A2)

For $1/\sqrt{2}\leq\kappa\leq 1$, there is no overlap if
$|\theta|>\sin^{-1}\Big{(}\kappa\Big{)}$. Here $\theta^{{}^{\prime}}_{max}$
has the same expression as given in eq.(A1) but $\theta^{{}^{\prime}}_{min}$
has different expression for two different ranges of $\theta$. For
$\theta\in[\cos^{-1}\Big{(}\frac{1}{2\kappa}\Big{)},\sin^{-1}\Big{(}\kappa\Big{)}]$

$\displaystyle\theta^{{}^{\prime}}_{min}=\theta+\sin^{-1}\Big{[}\frac{\sin{\theta}}{\kappa}\Big{]}$
(A3)

For $0<\theta<\cos^{-1}(\frac{1}{2\kappa})$ $\theta_{max}$ is given by
eq.(A1). A graphical representation of $\eta(\theta,\theta^{\prime})$ is given
in Fig. 10, for some representative values of $\ell/a$.

Given the functional form of $\eta(\theta,\theta^{\prime})$, the integral in
equation (6) of main text is easily done numerically, to determine $z$ as a
function of $l/a$.

Figure 10: $\eta(\theta,\theta^{{}^{\prime}})$ in
($\dfrac{\theta}{2\pi}$,$\dfrac{\theta^{{}^{\prime}}}{2\pi}$) plane. Shaded
area(red) corresponds to overlap region where $\eta$ is $1$. In the unshaded
area $\eta=0$.

## Appendix B: Transfer matrix calculation for the 1-dimensional chain of hard
thin rods pivoted at one end

We consider a system of thin hard rigid linear rods, each of length $\ell$.
One end of the rods is fixed on regular lattice of spacing $a$, and the rods
can rotate in a plane, see Fig. 11. The orientation of each rod is specified
by the angle it makes with the x-axis. Let the angle the $i$th rod makes be
called $\theta_{i}$, with $-\pi\leq\theta_{i}<\pi$.

Figure 11:  A system of hard rods pivoted at one end on a one dimensional
lattice.

Then we partition function is given by,

$\mathcal{Z}_{L}=\left[\prod_{i}\int\frac{d\theta_{i}}{2\pi}\right]\prod_{i=1}^{L-1}\left[1-\eta(\theta_{i+1},\theta_{i})\right]$
(B1)

We will obtain the entropy per rod using the transfer matrix method, which
provides a check of the algebraic method used in the main text. Also, we can
obtain the probability distribution of the angles ${\rm P}(\theta)$ in terms
of the eigenvector of the transfer matrix, which also agrees with the result
obtained in the main text.

AOO condition is clearly satisfied if $1/2\leq\kappa\leq 1$. Let
$\mathcal{Z}_{L}(\theta)$ be the restricted partition function of the system
of $L$ hard rods with left most rod having angle $\theta$. Then it satisfies
the following recursion relation

$\displaystyle\mathcal{Z}_{L+1}(\theta)=\int_{-\pi}^{\pi}\frac{d\theta^{\prime}}{2\pi}\mathcal{T}(\theta,\theta^{\prime})\mathcal{Z}_{L}(\theta^{\prime})$
(B2)

where ${\mathcal{T}}$ is the transfer matrix. It is $0$ when the adjacent rods
with angles ($\theta,\theta^{{}^{\prime}}$) overlap and 1 when they don’t
overlap see Fig. 10. Then the free energy per site in the thermodynamic limit
is given by $\log(\lambda)$ where $\lambda$ is the largest eigenvalue of the
transfer matrix ${\mathcal{T}}_{\kappa}$ , whose eigenequation is given by

$\displaystyle\int\frac{d\theta^{{}^{\prime}}}{2\pi}{\mathcal{T}}(\theta,\theta^{{}^{\prime}})\psi_{\lambda}^{R}(\theta^{{}^{\prime}})=\lambda\psi_{\lambda}^{R}(\theta)$
(B3)

where $\psi_{\lambda}^{R}(\theta)$ is the right eigenfunction of integral
kernel ${\mathcal{T}}(\theta,\theta^{\prime})$. Clearly

${\mathcal{T}}(\theta,\theta^{\prime})=1-\eta(\theta,\theta^{\prime}).$ (B4)

Now, $\eta(\theta,\theta^{\prime})$ is nonzero only if $\cos\theta>0$, and
$\cos\theta^{\prime}<0$. Thus, we see that the AOO condition is satisfied

$\eta(\theta,\theta^{\prime})\eta(\theta^{\prime},\theta^{\prime\prime})=0,{\rm~{}for~{}all~{}}\theta,\theta^{\prime},\theta^{\prime\prime},$
(B5)

Let us define function $f(\theta)$ given by

$f(\theta)=\int\frac{d\theta^{\prime}}{2\pi}\eta(\theta,\theta^{\prime}).$
(B6)

Clearly, $f(\theta)$ is nonzero only if,
$|\theta|<\cos^{-1}(\frac{1}{2\kappa})$ for $1/2\leq\kappa\leq\sqrt{2}$ and
$|\theta|<\sin^{-1}(\kappa)$ for $1/\sqrt{2}\leq\kappa\leq 1$ , and it equals
to
$[\theta^{{}^{\prime}}_{max}(\theta)-\theta^{\prime}_{min}(\theta)]/(2\pi)$.
Then it is easily seen that the right eigenfunction is given by

$\psi_{\lambda}^{R}(\theta)=1-\frac{1}{\lambda}f(\theta),$ (B7)

where eigenvalue $\lambda$ satisfies the equation

$\displaystyle\lambda^{2}-\lambda+A=0$ (B8)

where $A$ is the fraction area of the forbidden region (red region in Fig. 10)
, given by

$\displaystyle
A=\int\frac{d\theta}{2\pi}\int\frac{d\theta^{{}^{\prime}}}{2\pi}\eta(\theta,\theta^{\prime}).$
(B9)

Thus the larger eigenvalue is given by

$\lambda=\frac{1+\sqrt{1-4A}}{2}$ (B10)

Now it is easily that for dimer model in one dimension with activity $z$,
eigenvalue of the transfer matrix is given by,

$\lambda_{dimers}=\frac{1+\sqrt{1+4z}}{2}$ (B11)

Hence the model of hard rods is equivalent to the dimer model provided $z=-A$.
In order to obtain the probability distribution $P(\theta)$ we need both left
and right eigenfunctions of the transfer matrix. Because of the symmetry of
the problem left eigenfunction is simply,
$\psi_{\lambda}^{L}(\theta)=\psi_{\lambda}^{R}(\theta-\pi)$. Then $P(\theta)$
is given by (see Fig. 12)

$\displaystyle
P(\theta)=\frac{1-\dfrac{1}{\lambda}[f(\theta)+f(\theta-\pi)]}{1+\dfrac{2}{\lambda}z}.$
(B12)

.

Figure 12: Probability distribution of angles for 1-dimensional chain of hard
thin rods pivoted at one end. $a=1.5l$.

## Appendix C: Transfer matrix for the zigzag chain

We consider a slightly more complicated case of a one-dimensional zigzag
lattice, where the angle between adjacent bonds alternates between two values
$\phi$ and $-\phi$ [ Fig. 13]. At each site is a linear rigid thin rod of
length $\ell$ pivoted to lattice site at its midpoint. The orientation of the
rod at site $i$ will be specified with its angle $\theta_{i}$ it makes with
one of the bond directions. We can choose the angle to lie in the range
$[-\pi/2,+\pi/2]$.

Figure 13: System of hard rods pivoted at center on a 1-dimensional zigzag
lattice. Type 1 and Type 2 bonds make ($\pi/2-\phi/2$) and (-$\pi/2+\phi/2$)
angles from the x-axis respectively.

Each rod is free to rotate about its center as long as it do not overlap with
other rods.

In this case, in the transfer matrix method, we get two coupled eigenvalue
equations of the form

$\displaystyle\mathcal{T}_{1}\ket{g}=\lambda\ket{f}$ (C1)
$\displaystyle\mathcal{T}_{2}\ket{f}=\lambda\ket{g}$ (C2)

Here ${\mathcal{T}}_{1}$ and ${\mathcal{T}}_{2}$ are the transfer matrices
along the two alternating bonds. We write $\mathcal{T}_{1}$ and
$\mathcal{T}_{2}$ as

$\displaystyle\mathcal{T}_{1}=\ket{1}\bra{1}-\Delta_{1}$ (C3)
$\displaystyle\mathcal{T}_{2}=\ket{1}\bra{1}-\Delta_{2}$ (C4)

where $\ket{1}$ is a vector whose value in $\theta$ representation is always
one and $\Delta_{1,2}$ is the overlap operator which in $\theta$
representation takes value 1 when there is a overlap between adjacent rods and
takes value zero otherwise. In this case, the AOO condition becomes

$\displaystyle\Delta_{1}\Delta_{2}=\Delta_{2}\Delta_{1}=0.$ (C5)

Then the following vectors are solutions of eq.(C1) and eq.(C2)

$\displaystyle\ket{f}=(1-\frac{1}{\lambda}\Delta_{2})\ket{1}$ (C6)
$\displaystyle\ket{g}=(1-\frac{1}{\lambda}\Delta_{1})\ket{1}$ (C7)

where the eigenvalue $\lambda$ satisfies following quadratic equation

$\displaystyle\lambda^{2}-\lambda+A=0$ (C8)

where $A$ is , as before, given by

$\displaystyle A$ $\displaystyle=$
$\displaystyle\braket{1}{\Delta_{1}}{1}=\braket{1}{\Delta_{2}}{1}$ (C9)
$\displaystyle=$
$\displaystyle\int_{-\pi/2}^{\pi/2}\frac{d\theta}{\pi}\int_{-\pi/2}^{\pi/2}\frac{d\theta^{{}^{\prime}}}{\pi}\bra{\theta^{{}^{\prime}}}\Delta_{1,2}\ket{\theta}$

which is nothing but the area of the overlap region. Again we are able to
obtain the same result obtained by analytic method by using simple algebraic
property of the transfer matrix given by eq.(C5). In Fig. 14 colored area
corresponds to the region where AOO condition is satisfied in the ($a/l,\phi$)
plane.

Figure 14: Colored area corresponds the region where AOO condition is valid.
Figure 15: Probability distribution for system of rods are pivoted at the
center on a one-dimensional zigzag lattice. $a=0.8l$ and $\phi=\pi/2$

Probability distribution $P(\theta)$ can also calculated from eigenvectors,
given by (see Fig. 15),

$\displaystyle
P(\theta)=\frac{1-\dfrac{1}{\lambda}[\braket{\theta}{\Delta_{1}}{1}+\braket{\theta}{\Delta_{2}}{1}]}{1+\dfrac{2}{\lambda}z}.$
(C10)

where $z=-A$ as before.

##  Appendix D: Qualitative behaviour of the entropy for lattice spacing just
below $a^{*}$

In this section we will consider 2-dimensional model of hard rods pivoted at
one end on lattice site (see Fig. 2 of main text), when lattice spacing
$a=a^{*}(1-\epsilon)$ with $\epsilon\ll 1$. We present a heuristic argument
that for this case, we have

$\displaystyle s(l/a)$ $\displaystyle=$ $\displaystyle
g(z),~{}~{}~{}~{}~{}~{}~{}~{}a>a^{*}$ $\displaystyle=$ $\displaystyle
g(z)+k\epsilon^{3/2}+\textrm{higher order terms in}~{}\epsilon,a<a^{*}$

where $z$ is the dimer activity corresponding to lattice spacing $a$ given by
eq.(6) of main text.
If AOO condition is violated, there are configuration where one rotor overlaps
with two neighbours. A typical case is shown in Fig. 16

Figure 16: Rotor at $\bf{O}$ overlaps simultaneously with rotors at
$\bf{O}+\bf{\hat{e}_{y}}$ and $\bf{O}+\bf{\hat{e}_{y}}$.

The contribution of such configurations to the partition function is of the
form,

$\int\frac{d\theta_{o}}{2\pi}\,\int\frac{d\theta_{1}}{2\pi}\,\int\frac{d\theta_{2}}{2\pi}\eta(\theta_{o},\theta_{1})\,\eta(\theta_{o},\theta_{2}).$

It is clear from the geometry of the figure that non-zero contribution to the
integral happens only when, ($\theta_{o}-\pi/4$), ($\theta_{1}-3\pi/4$) and
($\theta_{2}+\pi/4$) must be of order $\epsilon^{1/2}$ and thus integral is of
order $\epsilon^{3/2}$.This integral makes positive contribution to the
partition function. Hence we get the result given in eq.(LABEL:aoov).

## References

  * (1) L. Onsager, The Effects of Shape on the Interaction of Colloidal Particles , Ann. N.Y. Acad. Sci. $\mathbb{51}$, 627 (1949).
  * (2) J. Timmermans, Plastic crystals: a historical review, J. Phys. Chem. Solids, 18 (1961), pp. 1-8.
  * (3) L. A. K. Staveley, Thermodynamic studies of molecular rotation in solids, J. Phys. Chem. Solids, Vol. 18, (1961) pp46-52.
  * (4) Jennifer M. Pringle, Patrick C. Howlett , Douglas R. MacFarlane and Maria Forsyth, Organic ionic plastic crystals: recent advances, J. Mater. Chem., 2010, 20, 2056-2062. DOI: 10.1039/B920406G
  * (5) E. Shalaev, K. Wu, S. Shamblin, J.F. Krzyzaniak, M. Descamps Crystalline mesophases: structure, mobility, and pharmaceutical properties Adv. Drug Deliv. Rev., 100 (2016), pp. 194-211.
  * (6) Li, B., Kawakita, Y., Ohira-Kawamura, S. et al. Colossal barocaloric effects in plastic crystals. Nature 567, 506–510 (2019). https://doi.org/10.1038/s41586-019-1042-5.
  * (7) J. A. Pople and F. E. Karasz, A Theory of fusion of molecular crystals: I. The effect of hindered rotation, J. Phys. Chem. Solids, 18, (1961) pp 28-39.
  * (8) L. M. Casey and L. K. Runnels, Model for Correlated Molecular Rotation ,J. Chem. Phys. $\mathbb{51}$, 5070 (1969).
  * (9) B. C. Freasier and L. K. Runnels , Classical rotators on a linear lattice, J. Chem. Phys. $\mathbb{58}$, 2963 (1973).
  * (10) Y. Kantor and Mehran Kardar , One-dimensional gas of hard needles,Phys. Rev E $\mathbb{79}$, 041109(2009).
  * (11) P. Gurin and S. Varga, Towards understanding the ordering behavior of hard needles: Analytical solutions in one dimension, Phys. Rev. E $\mathbb{83}$, 061710 (2011).
  * (12) S. Saryal, J. U. Klamser, T. Sadhu and Deepak Dhar, Multiple Singularities of the Equilibrium Free Energy in a One-Dimensional Model of Soft Rods, Phys. Rev. Lett. $\mathbb{121}$, 240601(2018).
  * (13) G. B. Guthrie and J. P McCullough, Some observations on phase tranformations in molecular crystals, J. Phys. Chem. Solids, 18 pp53-61(1961).
  * (14) Zohar Nussinov, Jeroen Van Den Brink, Compass models: Theory and physical motivations, Rev. Mod. Phys., Vol. 87 (2015) pp1-59.
  * (15) Alexei Kitaev, Anyons in an exactly solvable model and beyond, Ann. Phys. Vol. 321 (2006) pp2-111.
  * (16) G. Baskaran, Diptiman Sen, and R. Shankar, Spin-S Kitaev model: Classical ground states, order from disorder, and exact correlation functions, Phys. Rev. B 78 (2008) 115116.
  * (17) E. Lieb and Fa Yueh Wu, Two dimensional Ferroelectric models, in Phase transitions and critical phenomena, Vol. 1, Ed. C. Domb and M. S. Green, [Academic Press, London, 1972], pp 332-490.
  * (18) G. X. Viennot, in Proc. of the Colloque de Combinatoire Enumerative, Lecture Notes in Mathematics Vol. 1234, pp 321-350 (Springer-Verlag, Berlin, 1986).
  * (19) Julianne Klamser, Tridib Sadhu, and Deepak Dhar,
Sequence of phase transitions in a model of interacting rods,
arXiv:2202.03065.


# Ultrafast charge-transfer dynamics in twisted MoS2/WSe2 heterostructures

J. E. Zimmermann1, M. Axt1, F. Mooshammer2, P. Nagler2, C. Schüller2, T.
Korn3, U. Höfer1 and G. Mette1 1 Fachbereich Physik und Zentrum für
Materialwissenschaften, Philipps-Universität, 35032 Marburg, Germany  2
Institut für Experimentelle und Angewandte Physik, Universität Regensburg,
93053 Regensburg, Germany  3 Institut für Physik, Universität Rostock, 18059
Rostock, Germany<EMAIL_ADDRESS>
###### Abstract

Two-dimensional transition metal dichalcogenides (TMD) offer a unique platform
for creating van-der-Waals heterojunctions with fascinating physical
properties and promising applications in optoelectronics and valleytronics.
Because of their typical type-II band alignment, photoexcited electrons and
holes can separate via interfacial charge transfer. To understand the nature
and the dynamics of this charge transfer is of utmost importance for the
design and efficiency of potential devices. However, systematic studies
concerning the influence of the stacking angle on the charge transfer remain
sparse. Here, we apply time- and polarization resolved second-harmonic imaging
microscopy to investigate the charge-transfer dynamics for three MoS2/WSe2
heterostructures with different stacking angles at a previously unattainable
time-resolution of $\approx$ 6 fs. For 1.70 eV excitation energy, electron
transfer from WSe2 to MoS2 is found to depend considerably on the stacking
angle with the fastest transfer time observed to be as short as 12 fs. At 1.85
eV excitation energy, ultrafast hole transfer from MoS2 to hybridized states
at the $\Gamma$-point or to the K-points of WSe2 has to be considered.
Surprisingly, the corresponding decay dynamics show only a minor stacking-
angle dependence indicating that radiative recombination of indirect
$\Gamma$-K excitons becomes the dominant decay route for all samples.

Keywords: time-resolved second-harmonic generation, transition metal
dichalcogenides, ultrafast charge transfer, heterostructure, pump-probe
experiment, nonlinear optical spectroscopy, stacking angle

##  1 Introduction

Two-dimensional van-der-Waals materials offer a plethora of opportunities for
the design and investigation of stacked heterostructures (HS) 1, 2. In
particular, heterostructures of transition metal dichalcogenides (TMD) have
revealed fascinating properties stimulating fundamental and applied research
in the field of optoelectronics and valleytronics 3, 4. Many combinations of
different TMD layers form type-II heterojunctions 5 which enables efficient
charge separation and results in spatially separated electron-hole pairs
following an optical excitation (so called charge-transfer or interlayer
excitons) 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17. Due to their van-der-
Waals coupled nature, structures with arbitrary stacking angles can be
designed and manufactured. The relative orientation of the TMDs in turn
influences the coupling between the layers and therefore charge transfer,
recombination and other properties of the interlayer excitons 18, 19, 20, 11,
21, 22, 23, 24. However, recent experimental studies investigating the
influence of the stacking on the ultrafast charge transfer obtained
surprisingly diverse results. For MoS2/WSe2 heterostructures, charge transfer
has been reported to be faster than the experimental time resolution of 40 fs
independent of the examined stacking angles 25. In contrast, much slower
charge transfer times of a few hundred femtoseconds and a significant increase
for larger rotational mismatch were observed for WS2/WSe2 26.

Thus, the underlying mechanism for the charge transfer process for different
van-der-Waals heterostructures still remains elusive and there are several
unresolved questions concerning the influence of the stacking configuration on
the ultrafast interlayer charge transfer. We have recently shown the
capabilities of the experimental method of time- and polarization-resolved
second-harmonic (SH) imaging microscopy giving us direct access to the
directional ultrafast charge transfer in a rotationally mismatched WSe2/MoSe2
heterostructure 8. In particular, our technique allows for pump-probe
experiments in $\mu$m-sized regions of heterostructures with a time-resolution
not accessible in previous studies. Here, we employ this new approach to
examine the charge-carrier dynamics in MoS2/WSe2 heterostructures with
different stacking angles. The three MoS2/WSe2 samples studied in the present
work had been examined beforehand in a systematic study by Kunstmann  _et al._
11. There it has been conclusively shown that for all stacking configurations
photoluminescence from momentum-space indirect $\Gamma$-K interlayer excitons
is visible after optical excitation with 532 nm (2.33 eV) photons. It is
therefore established that for excitation energies significantly exceeding the
respective band gaps interlayer charge transfer occurs and that at least some
of the excited charge carriers recombine radiatively via a $\Gamma$-K
transition. Recently, also momentum-space direct K-K photoluminescence at
$\approx$1.0 eV has been observed for well-aligned MoS2/WSe2 heterostructures
27. In the present study, we want to elucidate the corresponding ultrafast
dynamics in dependence of the stacking angle for three selected stacking
angles: 9° (referred to as quasi-3R), 52° (quasi-2H) and 16° (misaligned).

##  2 Results and Discussion

Figure 1: (a) Optical microscopy image of the studied misaligned (16°)
MoS2/WSe2 sample. MoS2 and WSe2 monolayer as well as the heterostructure
regions are highlighted with green, orange and red outlines, respectively.
(b)-(d) SH images measured with three different probe polarizations sensitive
to both monolayers (b), exclusively to MoS2 (c) or WSe2 (d), respectively. (e)
Sketch of the pump-probe experiment for an example of resonant optical
excitation in WSe2 and selective SHG detection of potential charge transfer in
MoS2. (f) Polarisation-dependent SHG of the misaligned sample shown in (a)
evaluated for MoS2 (green) and WSe2 (orange) monolayers and the
heterostructure (red). The dashed lines mark the polarization angles I and II
applied for the SH images shown in (c) and (d). (g) Schematic band alignment
of MoS2/WSe2 showing the individual K-points of both materials and the
hybridized $\Gamma$-point 11. The dashed ellipses mark the previously reported
K-K and $\Gamma$-K interlayer excitons 11, 27. (h) Corresponding hexagonal
Brillouin zones of the studied heterostructures for the three stacking angles
9° (red), 16° (yellow) and 52° (blue). The dashed rectangles delineate the
regions of interest highlighted in Figs. 2(d) and 3(d) respectively.

Our method of choice is polarization- and time-resolved second-harmonic
imaging microscopy. By careful selection of the polarization angle of the
800-nm (1.55 eV) probe light, we are able to extract changes in the nonlinear
susceptibility of individual monolayers inside the TMD heterostructures. For
TMD monolayers, the pump-induced changes of the SH response upon resonant
optical excitation of intralayer excitons have been shown to be closely
correlated to the observed exciton dynamics in linear optical spectroscopy 28.
In case of TMD heterostructures, our pump-probe experiments schematically
sketched in Fig. 1(e) allow us to access the dynamics of interlayer charge
transfer 8. The three studied MoS2/WSe2 samples were mechanically exfoliated
and transferred onto a SiO2/Si(001) substrate as described in Ref. 11. Figs.
1(a)-(d) compare an optical microscopy image of one of the studied samples
with respective images obtained for different probe polarizations with our SH
imaging microscopy. Due to the inherent structural sensitivity of second-
harmonic generation (SHG), we are able to directly identify the different
monolayer regions and their overlap to determine the position of the
heterostructure.

Fig. 1(f) displays the normalized SH response in dependence of the
polarization of the probe laser evaluated for the monolayer and
heterostructure regions. From these measurements the crystal orientations in
relation to the lab coordinates can be determined. The probe polarization for
our time-resolved measurements is then chosen to maximize the sensitivity for
the MoS2 layer by fully suppressing the SH contribution of WSe2 (polarization
I). For 1.70 eV pump-photon energy intralayer A-excitons of WSe2 are
resonantly excited 29, 30. Fig. 2(a) shows the corresponding SH-transients of
the heterostructures selectively detected by probing the ultrafast nonlinear
response of the MoS2 layer in direct comparison with the respective transients
obtained from the individual monolayers. All SH transients have been
normalized to the signal at negative delays. As expected the MoS2 monolayer
transient (filled black data) shows no pump-induced change since the
excitation energy is below its A-exciton resonance. In contrast, the SH-
intensity of the MoS2 layers of the heterostructures exhibits a clear
modulation as a function of the pump-probe delay time. We observe a rapid
decrease and a subsequent slower recovery. Since any direct excitation of the
MoS2 can be excluded from the monolayer results, we can therefore assign the
observed dynamics in the heterostructures to ultrafast electron transfer to
MoS2 after optical excitation of WSe2 as sketched in Fig. 2(c). This
attribution is also corroborated by the delayed onset of the heterostructure
signals in relation to the WSe2 monolayer response and by the ongoing/delayed
decrease of the SH intensity, which extends even beyond the duration of the
pump pulse (cf. dashed line, CCR).

While there is clear evidence for charge transfer in all heterostructure
samples, the observed ultrafast dynamics differ considerably for the three
different stackings. Apparently, the fastest electron transfer can be found
for the quasi-2H stacked sample (blue data), followed by the quasi-3R (red
data) and the misaligned sample (yellow data). This is confirmed
quantitatively by a rate-equation model consisting of coupled differential
equations in which one energetically higher lying state is filled by a
gaussian pump pulse and then subsequently populates the lower lying state with
a transfer time $\Delta\mathrm{t}_{\mathrm{CT}}$. The extracted transfer times
for the quasi-3R and the misaligned sample are comparable with values of
65$\pm$10 fs and 84$\pm$10 fs, respectively, whereas the transfer time for the
quasi-2H stacked sample is more than five times faster (12$\pm$6 fs). This
disparity might be surprising since the variation of the lattice separation
between the layers, which has a central influence on the tunneling probability
31, is negligible between 2H- and 3R-stacking 11. However, it has also been
reported that the wavefunction overlap between the layers plays an important
role for the charge transfer across the interface as well as for the
subsequent recombination dynamics 19. Accordingly, wavefunction calculations
for the interfacial plane in TMD heterostructures have shown strong
differences for 3R- and 2H-stacking 32, which then in turn can indeed
influence the dynamics 33.

Figure 2: Time-resolved SHG of the MoS2/WSe2 heterostructures with three
different stackings. (a,b) SH transients upon 1.70 eV resonant excitation of
the WSe2 intralayer A-exciton for a sub-picosecond and an extended range of
pump-probe delays, respectively. Colored data points represent the
heterostructure samples, black data points the corresponding monolayer
dynamics as denoted in (b). Solid lines in (a) correspond to rate-equation
fits, those in (b) to exponential fits. Filled (unfilled) data points are
measured at probe polarization I (II), sensitive to MoS2 (WSe2). The dashed
line in (a) represents the cross-correlation (CCR) of the laser pulses. (c)
Sketch of the band alignment for MoS2/WSe2. Ultrafast charge transfer to MoS2
after resonant optical excitation of WSe2 results in the formation of K-K
interlayer excitons (dashed ellipse). SH detection selectively probes the
transient response of MoS2. (d) Cutout from the Brillouin zone illustrating
the electron transfer for the respective samples.

Consequently, we conclude that the faster interlayer electron transfer for the
quasi-2H stacked sample is caused by the stronger interlayer interaction of
the involved monolayer states at the K-point. In accordance, the observed
recombination dynamics shown in Fig. 2(b) also differ significantly for the
three heterostructures as determined from a single-exponential fit of the SH
transients to extract the individual lifetimes. The misaligned sample shows
the longest (467$\pm$16 ps), the quasi-3R stacked sample an intermediate
(309$\pm$12 ps) and the quasi-2H stacked sample the shortest lifetime
(158$\pm$12 ps). The enhanced lifetime for the misaligned sample could be an
indication for the absence of the radiative decay channel as was reported
recently27. The difference in lifetime between quasi-2H and quasi-3R can then
be explained by the same reasoning as the enhanced transfer. Since the charge
transfer via a phonon-assisted tunneling process can be amplified by the
spatial overlap of the involved wavefunctions 6, 31, the same should hold for
the reversed process of interlayer exciton recombination: A stronger spatial
coincidence facilitates the recombination resulting then in a reduced
lifetime. Please note, that all three heterostructure lifetimes are
considerably extended in comparison to the lifetime observed for the WSe2
monolayer. The latter can be described by a bi-exponential recovery behaviour
with much shorter lifetimes of $\tau_{1}$ = 3.5 ps and $\tau_{2}$ = 30 ps.
This enhanced lifetime for the heterostructure signals is in perfect agreement
with our interpretation, since the spatially indirect nature of interlayer
excitons leads to a strong increase of their lifetimes in comparison to
intralayer excitons 34, 15.

In the following, we would like to address open questions concerning the
ultrafast hole transfer and the particular influence that hybridized states at
the $\Gamma$-point have on charge transfer and recombination. In order to
elucidate these processes, the pump-photon energy is tuned to 1.85 eV for
resonant optical excitation of MoS2 intralayer A-excitons 29, 30. While at
this photon energy both materials are excited, the generated exciton density
in MoS2 (5.3$\cdot 10^{12}/$cm2) is about five times larger than in WSe2
(9.3$\cdot 10^{11}/$cm2), as calculated based on energy-dependent absorption
of the two layers (see method section). Thus, the pump-induced effects in WSe2
are not negligible, however, the dynamic response measured in MoS2 is mainly
dominated by its inherently generated charge carriers. Fig. 3 shows the same
time regimes as before using identical colors for the different structures,
but now the systems are pumped at 1.85 eV. By comparing the dynamics of Figs.
2(a) and 3(a) striking differences become obvious: First of all, upon 1.85-eV
photoexcitation the pump-induced effects now already start to occur precisely
at temporal overlap in accordance with a direct excitation in MoS2. With
respect to the delayed filling, the SH transients of the quasi-3R and the
misaligned sample stay similar, whereas the behaviour of the quasi-2H stacked
sample changes drastically. In contrast to the 1.70-eV excitation, where the
SH signal of the quasi-2H sample was observed to recover mono-exponentially
with a slow recovery rate, the SH response after 1.85-eV excitation begins to
recover immediately on a femtosecond timescale. The extracted lifetime
$(\tau_{\mathrm{2H}}\,=\,220\mathrm{\,fs})$ is only slightly slower than the
corresponding monolayer decay times for 1.85-eV excitation
($\tau_{\mathrm{MoS}_{2}}\,=\,75\mathrm{\,fs},\
\tau_{\mathrm{WSe}_{2}}\,=\,146\mathrm{\,fs}$).

Figure 3: Time-resolved SHG of the MoS2/WSe2 heterostructures with three
different stackings. (a,b) SH transients and corresponding fits as in Fig.
2(a,b) but for 1.85 eV pump-photon energy, resonant with the MoS2 intralayer
A-exciton. (c) Sketch of the energy alignment for MoS2/WSe2 following the same
scheme as in Fig 2(c), but for resonant excitation in MoS2 while still probing
MoS2. The higher excitation energy opens up additional relaxation pathways by
ultrafast hole transfer to the hybridized $\Gamma$-point or to the K/K’-points
of WSe2. Three different kinds of conceivable interlayer excitons are marked
by dashed ellipses. For the quasi-2H stacked sample, hole transfer to K’ is
spin-forbidden and only two transfer possibilities remain. (d) Cutout from the
Brillouin zone illustrating the charge recombination via the $\Gamma$-point.

This ultrafast decay process is associated with the coherent radiative
recombination of intralayer excitons 35, 36. On the one hand, the slightly
slower recovery for the heterostructure signal compared to the monolayer can
be explained by the change of the dielectric environment: The enhanced
screening for the heterostructure decreases the binding energy of the excitons
and therefore increases their Bohr radius and the recombination time. On the
other hand, the coherent radiative recombination competes with additional
relaxation mechanism in the heterostructure such as hole transfer to the
hybridized $\Gamma$-point or to the WSe2 layer as illustrated in Fig. 3(c).
Since at 1.85-eV excitation electron transfer from the WSe2 into the MoS2 is
feasible as well, the detected SH signal in the MoS2 is a mixture of various
processes and a clear interpretation of our SH transients becomes more
challenging. The competing processes might also explain why we do not observe
any clear ultrafast component for the quasi-3R and the misaligned
heterostructure sample in our data. Another plausible explanation for the
absence of the ultrafast decay in these two heterostructures could also be
related to the energy alignment and the spin structure of the involved states.
After photoexcitation inside the MoS2 the generated hole has two obvious
relaxation pathways alternative to intralayer recombination: it can either be
scattered to the hybridized states around the $\Gamma$-point or to the
K/K’-point in WSe2. In the case of 2H-stacking, transfer to the valence band
maximum of WSe2 at K’ is spin forbidden as sketched in Fig. 3(c). For the
misaligned sample, scattering to K/K’ should also be reduced due to the large
momentum mismatch. In contrast, hole transfer to the valence band maximum at
the K-point should be most efficient in case of 3R-stacking. Assuming the hole
transfer to K/K’ to be a competing process to the coherent recombination and
taking place on a similar or faster timescale, the decreasing efficiency of
the hole transfer to K/K’ from quasi-3R, misaligned to quasi-2H could then
explain the increasing manifestation of the coherent recombination.

The comparison of the SH transients at large pump-probe delays for both
excitation energies in Figs. 2(b) and 3(b) reveals a surprising change in the
dynamics caused by the higher excitation energy. For 1.70-eV excitation the
recovery rates of the three heterostructure samples differ considerably as
discussed above. In the case of 1.85-eV excitation, however, the lifetimes
become very similar for all samples (quasi-2H: 178$\pm$9 ps, quasi-3R:
217$\pm$8 ps, misaligned: 219$\pm$7 ps). For 1.70-eV excitation, the quasi-2H
configuration showed the fastest recovery time due to the larger wavefunction
overlap between the layers. For 1.85 eV, however, the lack of interaction
between the layers for the other two samples is compensated by the additional
decay routes of the charge carriers generated in the MoS2. The central
difference between the two excitation energies is that for 1.85-eV excitation
the hole states around the $\Gamma$-point become energetically available.
Photoluminescence from momentum-space indirect $\Gamma$-K interlayer excitons
has been observed for all examined stacking configurations of MoS2/WSe2 11.
Therefore, we attribute our observed decay dynamics to the radiative
recombination of $\Gamma$-K interlayer excitons. Due to the fact that the
states around the $\Gamma$-point consist inherently of orbitals which are
delocalized and therefore spread out over both layers the overlap of the
wavefunctions of electron and hole is enhanced and, thus, radiative
recombination is facilitated. This overall enhancement of the recombination
compensates the reduced overlap for the quasi-3R and misaligned sample
equalizing the lifetimes of the excitation independent of the stacking
configuration.

##  3 Conclusion

In conclusion, we have employed time- and polarization-resolved second-
harmonic imaging microscopy to study the ultrafast charge-carrier dynamics
across the MoS2/WSe2 heterostructure interface for different stacking
configurations. For lower excitation energies of 1.70 eV, electron transfer
from WSe2 to MoS2 is found to depend considerably on the stacking angle and
the transfer time is reduced by a factor of seven when going from a larger
rotational mismatch towards 2H-stacking. At higher excitation energies,
ultrafast hole transfer from MoS2 to hybridized states at the $\Gamma$-point
and to the K-points of WSe2 has to be considered. The respective decay
dynamics, however, does not show a significant dependence on the stacking
angle indicating that radiative recombination of indirect $\Gamma$-K excitons
becomes the dominant decay route for all samples.

Funding was provided by the Deutsche Forschungsgemeinschaft (DFG, German
Research Foundation), Project-ID 223848855-SFB 1083. T. Korn gratefully
acknowledges funding by the DFG via KO3612/4-1. C. Schüller gratefully
acknowledges funding by the DFG via Project-ID 314695032-SFB 1277 as well as
SPP2244 (SCHU1171/10-1).

##  4 Methods

The experiments were performed under ambient conditions using a Yb:KGW-based
ultrashort pulse laser system (Light Conversion Carbide) feeding two optical
parametric amplifiers (Orpheus-F Twin, Orpheus-N) providing photon energies in
the range from 1.3-1.9 eV at pulse lengths routinely shorter than 35 fs at a
repetition rate of 200 kHz. The photon energy of the probe pulse was chosen to
be 1.55 eV well below the bandgap of both materials in order to exclude
optical excitation by the probe beam. Our temporal resolution is only limited
by the laser system and can be estimated by measurement of the cross-
correlation on the sample. The CCR is measured to be 30 fs. The minimal time-
resolution is estimated to be 1/5th of the CCR and corresponds to 6 fs. Both
beams are focussed collinearly onto the sample under an angle of 18°. After
passing a 400-nm dielectric filter, the specular reflected SH response of the
probe beam is imaged optically magnified by a camera lens on an electron-
multiplied CCD chip (Princeton Instruments ProEM-HS). Typical exposure times
are 10 seconds. The applied magnification was M$\approx$35-40. The overall
resolution of our imaging microscopy setup is better than 4 $\mu$m. The time-
delay between pump and probe beam is varied by a motorized delay stage. The
polarization of pump and probe beam can be varied by means of
$\lambda$/2-plates. The typically detected p-polarization of the second-
harmonic light is separated by an analyzer. The applied pump fluence was fixed
to 50 $\mu$J/cm2 for 1.85 eV and 38 $\mu$J/cm2 for 1.70 eV on the sample (spot
diameter 1/e). The applied probe fluence was fixed to 250 $\mu$J/cm2. These
fluences lead to an initial exciton density of 3.65$\cdot$1012/cm2 (WSe2) for
1.70 eV and 9.3$\cdot$1011/cm2 (WSe2) as well as 5.3$\cdot$1012/cm2 (MoS2) for
1.85 eV pump energy (calculated with data from 30 with the tmm package 37).
These excitation densities are below the calculated Mott-density 38 ensuring
that we are measuring exciton and not plasma dynamics. Long term measurements
with these fluences applied did not exhibit any multishot damage. A more
detailed description of the setup can be found elsewhere 28. Sample
preparation and characterization is described in detail in Ref. 11.

## References

  * 1 A. K. Geim and I. V. Grigorieva.  Van Der Waals Heterostructures.  Nature, 499:419–25, 2013.
  * 2 H. Lim, S. I. Yoon, G. Kim, A. R. Jang, and H. S. Shin.  Stacking of Two-Dimensional Materials in Lateral and Vertical Directions.  Chem. Mater., 26:4891–903, 2014.
  * 3 K. F. Mak and J. Shan.  Photonics and Optoelectronics of 2d Semiconductor Transition Metal Dichalcogenides.  Nat. Photonics, 10:216–26, 2016.
  * 4 X. D. Xu, W. Yao, D. Xiao, and T. F. Heinz.  Spin and Pseudospins in Layered Transition Metal Dichalcogenides.  Nat. Phys., 10:343–50, 2014.
  * 5 V. O. Özcelik, J. G. Azadani, C. Yang, S. J. Koester, and T. Low.  Band alignment of two-dimensional semiconductors for designing heterostructures with momentum space matching.  Phys. Rev. B, 94:035125, 2016.
  * 6 Z. Wang, P. Altmann, C. Gadermaier, Y. Yang, W. Li, L. Ghirardini, C. Trovatello, M. Finazzi, L. Duò, M. Celebrano, R. Long, D. Akinwande, O. V. Prezhdo, G. Cerullo, and C. S. Dal.  Phonon-Mediated Interlayer Charge Separation and Recombination in a MoSe2/WSe2 Heterostructure.  Nano Letters, 21:2165–2173, 2021.
  * 7 S.-W. Zheng, H.-Y. Wang, L. Wang, H. Wang, and H.-B. Sun.  Layer-Dependent Electron Transfer and Recombination Processes in MoS2/WSe2 Multilayer Heterostructures.  The Journal of Physical Chemistry Letters, 11(22):9649–9655, 2020\.
  * 8 J. E. Zimmermann, Y. D. Kim, J. C. Hone., U. Höfer, and G. Mette.  Directional ultrafast charge transfer in a WSe2/MoSe2 heterostructure selectively probed by time-resolved SHG imaging microscopy.  Nanoscale Horiz., 5:1603–1609, 2020.
  * 9 E. V. Calman, L. H. Fowler-Gerace, D. J. Choksy, L. V. Butov, D. E. Nikonov, I. A. Young, S. Hu, A. Mishchenko, and A. K. Geim.  Indirect Excitons and Trions in MoSe2/WSe2 Van Der Waals Heterostructures.  Nano Lett., 20:1869–75, 2020.
  * 10 S. Ovesen, S. Brem, C. Linderalv, M. Kuisma, T. Korn, P. Erhart, M. Selig, and E. Malic.  Interlayer Exciton Dynamics in Van Der Waals Heterostructures.  Commun. Phys., 2:23, 2019.
  * 11 J. Kunstmann, F. Mooshammer, P. Nagler, A. Chaves, F. Stein, N. Paradiso, G. Plechinger, C. Strunk, C. Schüller, G. Seifert, D. R. Reichman, and T. Korn.  Momentum-Space Indirect Interlayer Excitons in Transition-Metal Dichalcogenide Van Der Waals Heterostructures.  Nat. Phys., 14:801–805, 2018.
  * 12 P. Rivera, H. Yu, K. L. Seyler, N. P. Wilson, W. Yao, and X. Xu.  Interlayer Valley Excitons in Heterobilayers of Transition Metal Dichalcogenides.  Nat. Nanotechnol., 13:1004, 2018.
  * 13 B. Miller, A. Steinhoff, B. Pano, J. Klein, F. Jahnke, A. Holleitner, and U. Wurstbauer.  Long-Lived Direct and Indirect Interlayer Excitons in Van Der Waals Heterostructures.  Nano Lett., 17:5229–37, 2017.
  * 14 J. R. Schaibley, P. Rivera, H. Y. Yu, K. L. Seyler, J. Q. Yan, D. G. Mandrus, T. Taniguchi, K. Watanabe, W. Yao, and X. D. Xu.  Directional Interlayer Spin-Valley Transfer in Two-Dimensional Heterostructures.  Nat. Commun., 7:13747, 2016.
  * 15 P. Rivera, J. R. Schaibley, A. M. Jones, J. S. Ross, S. F. Wu, G. Aivazian, P. Klement, K. Seyler, G. Clark, N. J. Ghimire, J. Q. Yan, D. G. Mandrus, W. Yao, and X. D. Xu.  Observation of Long-Lived Interlayer Excitons in Monolayer MoSe2-WSe2 Heterostructures.  Nat. Commun., 6:6242, 2015.
  * 16 H. Fang, C. Battaglia, C. Carraro, S. Nemsak, B. Ozdol, J. S. Kang, H. A. Bechtel, S. B. Desai, F. Kronast, A. A. Unal, G. Conti, C. Conlon, G. K. Palsson, M. C. Martin, A. M. Minor, C. S. Fadley, E. Yablonovitch, R. Maboudian, and A. Javey.  Strong Interlayer Coupling in Van Der Waals Heterostructures Built from Single-Layer Chalcogenides.  P. Natl. Acad. Sci. USA, 111:6198–202, 2014.
  * 17 X. P. Hong, J. Kim, S. F. Shi, Y. Zhang, C. H. Jin, Y. H. Sun, S. Tongay, J. Q. Wu, Y. F. Zhang, and F. Wang.  Ultrafast Charge Transfer in Atomically Thin MoS2/WS2 Heterostructures.  Nat. Nanotechnol., 9:682–86, 2014.
  * 18 M. Brotons-Gisbert, H. Baek, A. Molina-Sánchez, A. Campbell, E. Scerri, D. White, K. Watanabe, T. Taniguchi, C. Bonato, and B. D. Gerardot.  Spin–layer locking of interlayer excitons trapped in moiré potentials.  Nat. Materials, 19:630–636, 2020.
  * 19 P. Merkl, F. Mooshammer, S. Brem, A. Girnghuber, K.-Q. Lin, L. Weigl, M. Liebich, C.-K. Yong, R. Gillen, J. Maultzsch, J. M. Lupton, E. Malic, and R. Huber.  Twist-Tailoring Coulomb Correlations in Van Der Waals Homobilayers.  Nat. Commun., 11:2167, 2020.
  * 20 P. K. Gogoi, Y. C. Lin, R. Senga, H. P. Komsa, S. L. Wong, D. Z. Chi, A. V. Krasheninnikov, L. J. Li, M. B. H. Breese, S. J. Pennycook, A. T. S. Wee, and K. Suenaga.  Layer Rotation-Angle-Dependent Excitonic Absorption in Van Der Waals Heterostructures Revealed by Electron Energy Loss Spectroscopy.  ACS Nano, 13:9541–50, 2019.
  * 21 W.-T. Hsu, L.-S. Lu, P.-H. Wu, M.-H. Lee, P.-J. Chen, P.-Y. Wu, Y.-C. Chou, H.-T. Jeng, L.-J. Li, M.-W. Chu, and W.-H. Chang.  Negative Circular Polarization Emissions from WSe2/MoSe2 Commensurate Heterobilayers.  Nat. Commun., 9:1356, 2018.
  * 22 P. C. Yeh, W. Jin, N. Zaki, J. Kunstmann, D. Chenet, G. Arefe, J. T. Sadowski, J. I. Dadap, P. Sutter, J. Hone, and R. M. Osgood.  Direct Measurement of the Tunable Electronic Structure of Bilayer MoS2 by Interlayer Twist.  Nano Lett., 16:953–59, 2016.
  * 23 S. J. Zheng, L. F. Sun, X. H. Zhou, F. C. Liu, Z. Liu, Z. X. Shen, and H. J. Fan.  Coupling and Interlayer Exciton in Twist-Stacked WS2 Bilayers.  Adv. Opt. Mater., 3:1600–05, 2015.
  * 24 K. H. Liu, L. M. Zhang, T. Cao, C. H. Jin, D. A. Qiu, Q. Zhou, A. Zettl, P. D. Yang, S. G. Louie, and F. Wang.  Evolution of Interlayer Coupling in Twisted Molybdenum Disulfide Bilayers.  Nat. Commun., 5:4966, 2014.
  * 25 H. M. Zhu, J. Wang, Z. Z. Gong, Y. D. Kim, J. Hone, and X. Y. Zhu.  Interfacial Charge Transfer Circumventing Momentum Mismatch at Two-Dimensional Van Der Waals Heterojunctions.  Nano Lett., 17:3591–98, 2017.
  * 26 P. Merkl, F. Mooshammer, P. Steinleitner, A. Girnghuber, K. Q. Lin, P. Nagler, J. Holler, C. Schüller, J. M. Lupton, T. Korn, S. Ovesen, S. Brem, E. Malic, and R. Huber.  Ultrafast Transition between Exciton Phases in Van Der Waals Heterostructures.  Nat. Mater., 18:691–696, 2019.
  * 27 O. Karni, E. Barré, S. C Lau, R. Gillen, E. Y. Ma, B. Kim, K. Watanabe, T. Taniguchi, J. Maultzsch, K. Barmak, R. H. Page, and T. F. Heinz.  Infrared Interlayer Exciton Emission in MoS2/WSe2 Heterostructures.  Phys. Rev. Lett., 123:247402, 2019.
  * 28 J. E. Zimmermann, B. Li, J. C. Hone, U. Höfer, and G. Mette.  Second-harmonic imaging microscopy for time-resolved investigations of transition metal dichalcogenides.  Journal of Physics: Condensed Matter, 32(48):485901, 2020.
  * 29 D. Kozawa, R. Kumar, A. Carvalho, K. K. Amara, W. J. Zhao, S. F. Wang, M. L. Toh, R. M. Ribeiro, A. H. C. Neto, K. Matsuda, and G. Eda.  Photocarrier Relaxation Pathway in Two-Dimensional Semiconducting Transition Metal Dichalcogenides.  Nat. Commun., 5:4543, 2014.
  * 30 Y. L. Li, A. Chernikov, X. Zhang, A. Rigosi, H. M. Hill, A. M. van der Zande, D. A. Chenet, E. M. Shih, J. Hone, and T. F. Heinz.  Measurement of the Optical Dielectric Function of Monolayer Transition-Metal Dichalcogenides: MoS2, MoSe2, WS2, and WSe2.  Phys. Rev. B, 90:205422, 2014.
  * 31 H. Z. Zhou, Y. D. Zhao, W. J. Tao, Y. J. Li, Q. H. Zhou, and H. M. Zhu.  Controlling Exciton and Valley Dynamics in Two-Dimensional Heterostructures with Atomically Precise Interlayer Proximity.  ACS Nano, 14:4618–25, 2020.
  * 32 H. Wang, J. Bang, Y. Sun, L. Liang, D. West, V. Meunier, and S. Zhang.  The Role of Collective Motion in the Ultrafast Charge Transfer in Van Der Waals Heterostructures.  Nat. Commun., 7:11504, 2016.
  * 33 H. Zeng, X. Liu, H. Zhang, and X. Cheng.  New theoretical insights into the photoinduced carrier transfer dynamics in WS2/WSe2 van der Waals heterostructures.  Phys. Chem. Chem. Phys., 23:694–701, 2021.
  * 34 F. Ceballos, M. Z. Bellus, H. Y. Chiu, and H. Zhao.  Ultrafast Charge Separation and Indirect Exciton Formation in a MoS2-MoSe2 Van Der Waals Heterostructure.  ACS Nano, 8:12717–24, 2014.
  * 35 M. Selig, G. Berghäuser, M. Richter, R. Bratschitsch, A. Knorr, and E. Malic.  Dark and Bright Exciton Formation, Thermalization, and Photoluminescence in Monolayer Transition Metal Dichalcogenides.  2D Mater., 5:035017, 2018.
  * 36 C. Poellmann, P. Steinleitner, U. Leierseder, P. Nagler, G. Plechinger, M. Porer, R. Bratschitsch, C. Schüller, T. Korn, and R. Huber.  Resonant Internal Quantum Transitions and Femtosecond Radiative Decay of Excitons in Monolayer WSe2.  Nat. Mater., 14:889–893, 2015.
  * 37 Steven J. Byrnes.  Multilayer optical calculations.  arXiv:1603.02720, 2020.
  * 38 L. Meckbach, T. Stroucken, and S. W. Koch.  Giant Excitation Induced Bandgap Renormalization in TMDC Monolayers.  Appl. Phys. Lett., 112:061104, 2018.