# Near-field enhancement by waveguide-plasmon polaritons in a nonlocal
metasurface

Xiaorun Zang and Andriy Shevchenko  Department of Applied Physics, Aalto
University, P.O. Box 13500, FI-00076 Aalto, Finland<EMAIL_ADDRESS>
###### Abstract

Localized surface plasmons in metal nanoparticles are widely used in nano-
optics to confine and enhance optical fields. It has been previously shown
that, if the nanoparticles are distributed periodically, an additional
enhancement can be achieved by coupling the localized surface plasmons to the
diffraction orders of the lattice, forming surface lattice resonances. In this
work, we study an even further improvement of the near-field enhancement by
placing a metal-dielectric slab waveguide beneath the lattice of the particles
to excite coupled waveguide-plasmon polaritons. These excitations can extend
over many periods of the lattice, making the metasurface highly nonlocal. We
numerically demonstrate that the approach can provide a significant extra
increase in the near-field intensity – by a factor of 80 over that produced by
a single-particle plasmon resonance and by 7 over the lattice-resonance
enhancement. The described enhancement mechanism can be used to design
extraordinarily efficient nonlocal optical metasurfaces for many applications,
including surface-enhanced Raman spectroscopy, fluorescence spectroscopy,
nonlinear optics, and solar energy harvesting.

Keywords: Metasurfaces, Waveguide-plasmon polaritons, Surface plasmons, Bloch
wave, SERS, Plasmonics, Field enhancement

##  1 Introduction

Local enhancement of optical fields plays an essential role in many
applications of optics and photonics, including optical sensors and detectors
[1, 2], light sources [3, 4, 5, 6], photovoltaic devices [7, 8, 9], and
nonlinear optical components [10, 11, 12]. An efficient near-field enhancement
can be obtained with the help of metal nanostructures that exhibit localized
surface plasmon resonances (LSPRs) [13, 14, 15]. Due to a high density of
metal electrons, light can be strongly squeezed into subwavelength regions
[15]. The squeezing is especially effective near locations of abrupt and
small-volume inhomogeneities of the structures, including interparticle gaps
and sharp edges. At such inhomogeneities, very high field intensities can be
obtained. For example, the optical intensity in the gap of a gold bowtie
antenna consisting of two tip-to-tip facing prisms with sharp corners provides
a higher enhancement than the analogue dimer of nanodiscs with round surfaces
[16, 17]. These "hot spots" are in the core of many applications, for example,
in surface-enhanced Raman spectroscopy (SERS) [18, 19]. The electromagnetic
enhancement factor of the SERS signal is given by
$\eta_{\mathrm{SERS}}\approx|E/E_{0}|^{4}$ [18], where $E$ and $E_{0}$ are,
respectively, the total and incident electric fields in the hot spot.

Nanoparticles of various materials, compositions, shapes, sizes, and
surrounding dielectric environments have been explored to control the local
field intensity [20, 21, 22]. In addition to conventional approaches, high-
index dielectric nanoparticles have been recently used to resonantly enhance
both the electric and magnetic fields while significantly reducing the ohmic
losses [23, 24, 25, 26]. Plasmonic dimers that exhibit a toroidal dipole
resonance have been shown to outperform the conventional disc dimers regarding
the field enhancement and third harmonic generation [27]. It has also been
shown that, reducing the gap size increases the local field intensity up to a
certain limit dictated by quantum effects, such as the effect of electron
tunnelling through the gap and nonlocal screening [28, 29, 30, 31, 32, 33]. On
the other hand, the quantum confinement effect can play a crucial role for
plasmonic particles that are smaller than $10$ nm in size [28]. Recently,
quantum confinement of electrons in a one-dimensional metallic single-walled
carbon nanotube with a diameter of $1$-$2$ nm was reported to form a quantized
Luttinger-liquid plasmon that exhibits a strong spatial confinement and low
attenuation loss. Such quantum plasmons could be useful for plasmonic
integrated devices [34].

Much higher local fields can be obtained by arranging plasmonic nanoparticles
in a lattice [35, 36] or superlattice [37, 38]. In this case, the field
scattered by one particle can arrive at the neighbouring particles in phase
with the other scattered fields, when the periods are appropriate [39, 40,
38]. This results in collective LSPRs, which have also been called plasmonic
surface lattice resonances [41, 38, 42] or collective lattice modes [39, 40,
38, 43]. These excitations lead to a further increase in the local field
around each nanoparticle. It has been estimated that a surface lattice
resonance can increase the local field intensity by an order of magnitude [44,
45, 46]. In the case of a superlattice, the field enhancement can be achieved
simultaneously at multiple modes [37]. Moreover, LSPRs can be coupled to
optical modes guided in the substrate of a planar metasurface [47], leading to
the so-called waveguide-plasmon polaritons (WPPs) [48, 49]. It has been shown
that a hyperbolic metasurface can support such hybrid waves or quasi-TE and
quasi-TM modes [50]. Very recently, an array of lossy plasmonic nanoparticles
on a dielectric slab waveguide was shown to guide a resonant mode of hybrid
polarization states with an ultrahigh quality factor [51]. WPPs hold great
potential for yielding extra near-field enhancement in arrays of plasmonic
nanoparticles, and in our opinion, an in-depth investigation into this
mechanism is highly relevant and important.

In this work, we design and study a nonlocal metasurface [52] that provides a
significant additional near-field enhancement by WPPs. The metasurface
consists of a lattice of silver nanodimers mounted on a metal-dielectric slab
waveguide. Besides a gap enhancement in each nanodimer, an additional
intensity enhancement – by a factor of $80$ – arises from the coupling of
LSPRs to the guided Bloch modes. These modes belong to the second stop band of
the structure, exhibiting a collective resonance at normal incidence.
Metasurfaces of this type can be used for many applications, but for the
purpose of demonstration, we consider a SERS substrate operating at a
wavelength of $780$ nm, one of typical wavelengths in Raman spectroscopy [53].
At this wavelength, the SERS enhancement in each gap of the dimers is more
than $6000$ times higher than that in an isolated dimer on a semi-infinite
dielectric. The metasurface can easily be redesigned to enhance an optical
field at any other wavelength. In addition, nonlocal metasurfaces consisting
of meta-atoms other than nanodisc dimers, such as bowtie antennae or trimers,
are expected to benefit from the extra near-field enhancement by WPPs as well.

The paper is organized as follows. In Section 2, we study the guided modes of
a one-dimensional slab waveguide. This study provides a starting point for the
design of a planar metasurface and a physical insight into the properties of
possible guided Bloch modes. In Section 3, we establish a connection between
the near-field enhancement spectra in a designed planar metasurface and the
empty-lattice Bloch-mode band diagrams. In Section 4, we optimize the designed
metasurface to maximize its near-field enhancement factor. To achieve a
significant additional enhancement of light in the gap hot spots, we match the
frequencies of the lattice resonances and the guided modes by varying the
lattice periods and slab thickness. We draw conclusions in Section 5.

##  2 Guided modes in a slab waveguide

Figure 1: (a) A schematic diagram of a silver-silica-air slab waveguide. From
top to bottom, the air, silica and silver layers have refractive indices
$n_{\mathrm{air}}$, $n_{\mathrm{SiO_{2}}},$ and $n_{\mathrm{Ag}}$,
respectively. The silica slab layer has a thickness of $t_{\mathrm{slab}}$.
Field profiles of the fundamental (b) TM and (c) TE modes when the slab
thickness is $t_{\mathrm{slab}}=200$ nm. The dispersion curves are plotted for
(d) $t_{\mathrm{slab}}=40$ nm and (e) $t_{\mathrm{slab}}=200$ nm. Blue and
orange dots indicate the TM and TE modes, respectively. At the target
frequency with $\omega/\omega_{\mathrm{tg}}=1$, the propagation constants of
the fundamental TM mode $\beta_{\mathrm{TM_{0}}}^{\mathrm{tg}}$ and TE mode
$\beta_{\mathrm{TE_{0}}}^{\mathrm{tg}}$ are indicated, respectively, by
vertical blue and orange dashed lines. The black solid line is the light line
in vacuum, whereas the grey solid curve plots the SPP dispersion curve at the
silver-silica interface.

A metasurface considered in this work is composed of a slab waveguide and a
lattice of nanoparticles (meta-atoms) on its surface. We first study the
properties of the waveguide alone, treating it as a periodic structure with an
empty lattice of meta-atoms [54, 40]. A group of qualitatively similar
dispersion relations can then be found in a metasurface based on the same slab
waveguide [54, 55]. Therefore, it is useful to first investigate the modal
dispersion of a slab waveguide and then apply the obtained knowledge to the
problem of the near-field enhancement in a real 2D metasurface.

The waveguide is made of a silica slab on a silver substrate [see Fig. 1(a)].
The silver substrate can also enhance the field on the surface by simply
reflecting the incident field. The refractive indices of silver, silica, and
air are denoted by $n_{\mathrm{Ag}}$, $n_{\mathrm{SiO2}}$, and
$n_{\mathrm{air}}$, respectively. The dispersive refractive index of silver is
downloaded from the online refractiveindex.info database [56], which uses
Johnson and Christy’s tabulated data [57]. The silica and air regions are
treated as having constant refractive indices of $1.4537$ and $1$,
respectively. The waveguide modes are calculated by solving the following one-
dimensional problem in a weak formulation, using the finite element method
(FEM),

$\int_{C}dl\left[\frac{1}{\alpha}\frac{\partial v_{z}}{\partial
x}\frac{\partial u_{z}}{\partial x}-\left(\frac{\omega}{c}\right)^{2}\gamma
v_{z}u_{z}\right]=-\beta^{2}\int_{C}dl\frac{1}{\alpha}v_{z}u_{z},$ (1)

where $u_{z}$ is the only unknown transverse component of the mode as a
function of $x$, and $v_{z}$ is the corresponding test function. For a TM
mode, $u_{z}$ denotes the transverse magnetic field and the domain-wise
parameters are $\alpha=\epsilon$ and $\gamma=\mu$. For a TE mode, $u_{z}$
represents the transverse electric field and the parameters are $\alpha=\mu$
and $\gamma=\epsilon$. The modal propagation constant $\beta$ is sought by
solving the above eigenvalue problem for $-\beta^{2}$. We have implemented the
weak formulation calculations using both the Equation-based Module in COMSOL
Multiphysics [58] and the open-source software NGSolve [59], which generate
consistent results. In the calculations, the semi-infinite air and silver
regions are truncated by Perfectly Matched Layers (PMLs), which are modelled
as frequency-dependent, anisotropic absorbing media [60, 61]. The outside
boundaries of PMLs are those of a perfect electric conductor for a TE mode and
a perfect magnetic conductor for a TM mode.

The field profiles of the fundamental $\mathrm{TM}$ and $\mathrm{TE}$ modes
are shown in Figs. 1(b) and 1(c), respectively, for a slab thickness of
$t_{\mathrm{slab}}=200$ nm at the wavelength of $780$ nm. The $\mathrm{TM}$
mode has a dominant field component $E_{y}$ due to the presence of silver at
the bottom of the waveguide. It is more concentrated inside the slab and
resembles a surface plasmon polariton (SPP) mode at a silver-silica interface
[14]. However, the $\mathrm{TE}$ mode has only the $E_{z}$ component that is
nearly negligible at the metal surface and has intensity maxima closer to the
silica-air interface. In addition, the $\mathrm{TM}$ mode has a shorter modal
wavelength than the $\mathrm{TE}$ mode, with the former being
$\lambda_{\mathrm{TM_{0}}}=550$ nm and the latter being
$\lambda_{\mathrm{TE_{0}}}=753$ nm. The modal dispersion curves of slab
waveguides with $t_{\mathrm{slab}}=40$ nm and $t_{\mathrm{slab}}=200$ nm are
plotted in Figs. 1(d) and 1(e), respectively. For relatively thin slabs, the
fundamental $\mathrm{TM}$ mode approximately follows the dispersion curve of
the $\mathrm{SPP}$ on a silver-silica interface at high frequencies, while it
approaches the light line in air at low frequencies. If the slab thickness
increases, more modes emerge and the dispersion curve of the fundamental
$\mathrm{TM}$ mode is pushed towards the $\mathrm{SPP}$ curve. The target
wavelength of $780$ nm is shown by the horizontal solid grey line that
intersects with the dispersion curves and determines the modal propagation
constants $\beta_{\mathrm{TM_{0}}}^{\mathrm{tg}}$ and
$\beta_{\mathrm{TE_{0}}}^{\mathrm{tg}}$ at this wavelength.

Figure 2: (a) Mode wavelengths
$\lambda_{\mathrm{wg}}=2\pi/\beta_{\mathrm{wg}}^{r}$ of the slab waveguide;
(b) Normalized decay lengths
$2\pi/\beta_{\mathrm{wg}}^{i}/\lambda_{\mathrm{wg}}$ for a varying slab
thickness. The real and imaginary parts of the propagation constant are
denoted by $\beta_{\mathrm{wg}}^{r}$ and $\beta_{\mathrm{wg}}^{i}$,
respectively.

The wavelengths of the guided modes,
$\lambda_{\mathrm{wg}}=2\pi/\beta_{\mathrm{wg}}^{r}$, with
$\beta_{\mathrm{wg}}^{r}$ denoting the real part of the propagation constant,
are shown in Fig. 2(a) for a varying slab thickness. In addition, the
normalized decay lengths, $2\pi/\beta_{\mathrm{wg}}^{i}/\lambda_{\mathrm{wg}}$
[62], are plotted in Fig. 2(b), where $\beta_{\mathrm{wg}}^{i}$ is the
imaginary part of the propagation constant. The guided modes can therefore
propagate over a distance of $10^{3}-10^{4}$ modal wavelengths, leading to a
long-distance, highly nonlocal interaction of the meta-atoms.

##  3 Bloch modes in plasmonic metasurfaces

Before discussing the guided Bloch modes and lattice resonances in our
designed metasurfaces, we extend the empty-lattice dispersion relation of the
free-photon approximation [54, 40] to the free-guided-wave approximation. The
empty-lattice dispersion relation is

$\left|\mathbf{k}_{||}+\mathbf{G}\right|=\beta_{\mathrm{wg}},$ (2)

where $\mathbf{k}_{||}=\hat{y}k_{y}+\hat{z}k_{z}$ is the in-plane Bloch wave
vector, $\mathbf{G}=m\frac{2\pi}{\Lambda_{y}}+n\frac{2\pi}{\Lambda_{z}}$ is
the reciprocal lattice vector (with the integers $m$ and $n$ denoting the
Bloch-mode orders along the $y$ and $z$ directions, respectively), and
$\beta_{\mathrm{wg}}$ is the propagation constant of the slab mode. The slab
waveguide is homogeneous and extends to infinity in the transverse plane.
However, the fictitious periodicity causes the band folding [54, 40, 55, 63].
If we choose the periods to be equal to the guided-mode wavelengths, the
empty-lattice Bloch modes appear to be formed at the $\Gamma$-point with
$k_{y}=k_{z}=0$, which can be used for designing the Bloch modes in a periodic
metasurface to be excited at normal incidence [63].

Figure 3: Schematic diagrams of (a) a metal-dielectric metasurface and (b) the
corresponding empty lattice. The dimer consists of two identical nanodiscs
with a diameter of $w$ and a height of $h$, which are separated by a gap of
$g$. The two in-plane periods are $\Lambda_{y}$ and $\Lambda_{z}$. In the
empty lattice, the dashed outline implies that the nanodiscs are actually
removed. Note that the dimers in our calculations are about two times smaller
than in the pictures. The parameter values used in the calculations are
$\Lambda_{y}=550$ nm, $\Lambda_{z}=600$ nm, and $t_{\mathrm{slab}}=200$ nm.
The near-field enhancement in the metasurface and the empty-lattice dispersion
curves are shown in (c) and (d).

In a metasurface containing silver meta-atoms, the meta-atomic LSPRs are
coupled to the empty-lattice Bloch modes, resulting in modified Bloch modes
[40]. The unit cell of the metasurface is shown in Fig. 3(a). It contains a
meta-atom in the form of a silver-disc dimer with dimensions $w=120$ nm,
$h=20$ nm, and $g=6$ nm. The values of other parameters of the structure are
$t_{\mathrm{slab}}=200$ nm, $\Lambda_{y}=550$ nm, and $\Lambda_{z}=600$ nm. At
$\lambda=780$ nm, the period $\Lambda_{y}$ is equal to the mode wavelength
$\lambda_{\mathrm{TM_{0}}}$, but the period $\Lambda_{z}$ is shorter than
$\lambda_{\mathrm{TE_{0}}}$. Note that plasmonic dimers are among the most
efficient standard nanoparticles used to locally enhance optical fields, and
dimers composed of nanodiscs rather than particles with sharp corners are
relatively easy to nanofabricate and treat numerically.

The near-field enhancement spectra in the nonlocal metasurface are calculated
using the frequency-domain FEM of the Wave Optics Module of COMSOL
Multiphysics [58]. The background excitation field is either a $p$\- or
$s$-polarized plane wave in the scattered-field formulation. A $p$-polarized
incident plane wave (with the electric field lying in the $xy$-plane) couples
effectively to the TM Bloch modes of the structure, yielding Fano-type [64]
asymmetric line shapes in the near-field enhancement spectra. This is revealed
in Fig. 3(c), where high (in black) and low (in white) near-field enhancement
factors appear near the empty-lattice dispersion curves. A band gap showing
low near-field enhancement is formed around the normalized frequency
$\omega/\omega_{\mathrm{tg}}=1$ at $k_{y}=0$, as a result of coupling of the
LSPR with the empty-lattice Bloch modes $\mathrm{TM_{0}^{(+1,0)}}$ and
$\mathrm{TM_{0}^{(-1,0)}}$ that counter-propagate along axis $y$. For higher-
order empty-lattice TM modes with order $n=1$ (the blue dashed curve), another
band gap is formed near the normalized frequency
$\omega/\omega_{\mathrm{tg}}=1.04$ at $k_{y}=\pi/\Lambda_{y}$.

If light is incident in the $zx$-plane and $s$-polarized, it couples
effectively to the Bloch modes propagating in the $z$-direction. The excited
LSPRs in the dimers still provide a high gap enhancement. In Fig. 3(d), one
can observe the near-field enhancement spectra shown in band gaps near
$\omega/\omega_{\mathrm{tg}}=1.0$ and $1.08$ due to the interaction of LSPRs
with the $\mathrm{TM}_{0}^{(1,n)}$ modes at $k_{z}=0$ and
$k_{z}=\pi/\Lambda_{z}$, respectively. The interaction leads to high near-
field enhancement factors (see the dark regions in the figure). Additional
empty-lattice TE and TM modes do not significantly affect the near-field
enhancement and they are not shown in the figure.

##  4 Near-field enhancement in a planar metasurface

Figure 4: Near-field enhancement factor $|E/E_{0}|$ as a function of periods
$\Lambda_{y}$ and $\Lambda_{z}$ for several slab thicknesses
$t_{\mathrm{slab}}$ obtained at the fixed wavelength
$\lambda_{\mathrm{tg}}=780$ nm. Thick black ticks with labels
$\lambda_{\mathrm{TM}_{0}}$, $\lambda_{\mathrm{TM}_{1}}$, and
$\lambda_{\mathrm{TE}_{0}}$ mark the guided-mode wavelengths in the slab
waveguide of the corresponding thickness.

We can now maximize the near-field enhancement by optimizing the periods
$\Lambda_{y}$ and $\Lambda_{z}$, as well as the slab thickness
$t_{\mathrm{slab}}$. As shown in Fig. 4(a-f), both the periods affect the
near-field enhancement, and the fingerprints of band gaps are recognized as
straight stripes of low near-field enhancement near the periods equal to the
guided-mode wavelengths (indicated by thick black ticks with labels
$\lambda_{\mathrm{TM}_{0}}$, $\lambda_{\mathrm{TM}_{1}}$, and
$\lambda_{\mathrm{TE}_{0}}$). The Bloch modes formed above or below the
guided-mode wavelengths contribute the most to the near-field enhancement in
the gaps of the nanoparticles. A maximum field enhancement, $|E/E_{0}|=577$,
is achieved at periods $\Lambda_{z}=595$ nm and $\Lambda_{y}=760$ nm for the
slab thickness $t_{\mathrm{slab}}=200$ nm. Such a configuration allows for a
SERS enhancement factor, $\eta_{\mathrm{SERS}}$, of the order of $10^{11}$ at
a gap size of $6$ nm. At smaller gaps, the enhancement factor is larger, but
we consider 6 nm as a reasonable gap size that can be obtained in
nanofabricated samples.

Figure 5: Near-field enhancement factors of (a) a single dimer as a function
of wavelength $\lambda$ and (b) a lattice of dimers on top of a semi-infinite
silica as a function of periods $\Lambda_{y}$ and $\Lambda_{z}$.

It is interesting to compare the near-field enhancement in the optimized
metasurface to that of a single dimer and a lattice of dimers on a thick
silica substrate. For an isolated single dimer on a thick substrate, the near
field can be resonantly enhanced by a factor of $65$ at $\lambda=780$ nm, as
shown in Fig. 5(a). If such dimers are arranged periodically in a lattice on
the same substrate, an extra factor of $3$ can be obtained for the near-field
enhancement [see Fig. 5(b)]. In our design, with a silver mirror positioned at
a distance of $200$ nm below the surface, an overall near-field enhancement of
$577$ is obtained. Hence, arranging the particles in a lattice on a slab
waveguide resulted in about $80$ times higher local-field intensity in the
gaps compared to the case of isolated dimers on a thick substrate, providing
an extra SERS enhancement of more than $6000$ (calculated as a squared
intensity enhancement). The presence of the waveguide additionally increases
the local-field intensity by more than $7$ times, corresponding to a $50$-fold
increase of $\eta_{\mathrm{SERS}}$. It is well known that naked metal
nanostructures of typical SERS substrates are not suitable for multiple use
because of their fast and irreversible degradation caused, e.g., by humidity
variations, surface oxidation, and contamination. The structures could be
protected with a few-nm-thick dielectric coating, but this usually lowers
$\eta_{\mathrm{SERS}}$ by orders of magnitude. The extra enhancement
considered in this work can compensate for this decrease of SERS signal and
allow making the substrates reusable.

Figure 6: Near-field enhancement spectra for $t_{\mathrm{slab}}=200$ nm and
(a) a fixed period $\Lambda_{y}=760$ nm, (b) varying periods $\Lambda_{y}$ and
$\Lambda_{z}$, and (c) a fixed period $\Lambda_{z}=595$ nm. For three typical
configurations at points D, E, and F, the near-field enhancement spectra as a
function of the incidence angles along the $y$ and $z$ axes are shown in (d),
(e), and (f), respectively. Blue and orange lines indicate the empty-lattice
Bloch modes $\mathrm{TM}_{0}^{(m,n)}$ and $\mathrm{TE}_{0}^{(m,n)}$,
respectively. In both colours, solid lines indicate the zero-order modes,
whereas the dashed lines stand for the $1$st-order modes.

We further explore the near-field enhancement spectra of three designed
metasurfaces with $t_{\mathrm{slab}}=200$ nm, marked by points D, E, and F in
Fig. 6(a)-(c). For the metasurface at point D, the periods are
$\Lambda_{y}=550$ nm and $\Lambda_{z}=595$ nm, and the near-field enhancement
is very low due to the band gap near the normalized frequency
$\omega/\omega_{\mathrm{tg}}=1$ at $k_{y}=0$ [see Fig. 6(d)]. This band gap is
formed because the period $\Lambda_{y}$ matches with the modal wavelength
$\lambda_{\mathrm{TM}_{0}}=550$ nm of the fundamental $\mathrm{TM}$ mode
propagating along the $y$-axis [see Fig. 2(a)].

At point E, the periods are $\Lambda_{y}=760$ nm and $\Lambda_{z}=755$ nm,
with the latter being close to the modal wavelength
$\lambda_{\mathrm{TE}_{0}}=753$ nm of the fundamental $\mathrm{TE}$ mode
propagating along the $z$-axis [see Fig. 2(a)]. This results in a band gap
near the normalized frequency $\omega/\omega_{\mathrm{tg}}=1$ at $k_{z}=0$
[see Fig. 6(e)], which makes the near-field enhancement low.

The maximum near-field enhancement is obtained for the metasurface at point F
with a pair of periods $\Lambda_{y}=760$ nm and $\Lambda_{z}=595$ nm. As shown
in Fig. 6(f), the upper band of the guided Bloch mode is formed at the
normalized frequency $\omega/\omega_{\mathrm{tg}}=1$ at $k_{y}=0$, because
$\Lambda_{y}$ is larger than the modal wavelength $\lambda_{\mathrm{TM}_{0}}$.
On the other hand, the lower band of the guided Bloch mode is formed at the
normalized frequency $\omega/\omega_{\mathrm{tg}}=1$ at $k_{z}=0$, as
$\Lambda_{z}$ is smaller than the modal wavelength
$\lambda_{\mathrm{TE}_{0}}$. The upper and lower bands of the guided Bloch
modes counter-propagating along the $y$ and $z$ axes jointly enhance the near-
field in the gaps of the dimers.

##  5 Conclusions

We have shown that, with the aid of guided Bloch modes in a designed nonlocal
optical metasurface, the local field intensity can be considerably increased
in addition to the enhancement by isolated plasmonic nanoparticles. The
additional enhancement mechanism relies on the coupling of the incident light
to the LSPRs of individual nanoparticles and to the Bloch modes guided in a
metal-dielectric slab waveguide beneath the particles. The coupling takes
place at the second stop band edge for normally incident light. The maximum
enhancement is obtained by optimizing the two periods and the thickness of the
waveguide. The approach clearly outperforms the lattice-resonance-based
enhancement. For plasmonic nanodimers with a gap size of $6$ nm, we obtained a
SERS enhancement of $10^{11}$. By reducing the gap, the enhancement factor can
be increased by further orders of magnitude. In addition, using dimers with
sharp edges, such as bowtie antennae, would further enhance the local field
intensity. We expect that the proposed mechanism of additional near-field
enhancement will find applications not only in SERS, but also in fluorescence-
and scattering-based plasmonic sensing, light-energy harvesting components,
and nonlinear optics.

We acknowledge the Academy of Finland Flagship Programme Photonics Research
and Innovation (PREIN; grant No. 320167). We thank Dr. R. Kolkowski for
fruitful discussions and advice on theoretical simulations of metasurfaces and
slab waveguide modes. We also acknowledge the computational resources provided
by the Aalto Science-IT project.

## References

  * [1] Bo Liedberg, Claes Nylander, and Ingemar Lunström.  Surface plasmon resonance for gas detection and biosensing.  Sensors and Actuators, 4:299–304, January 1983.  doi:10.1016/0250-6874(83)85036-7.
  * [2] Jiří Homola, Sinclair S. Yee, and Günter Gauglitz.  Surface plasmon resonance sensors: review.  Sensors and Actuators B: Chemical, 54(1):3–15, January 1999.  doi:10.1016/S0925-4005(98)00321-9.
  * [3] Rupert F. Oulton, Volker J. Sorger, Thomas Zentgraf, Ren-Min Ma, Christopher Gladden, Lun Dai, Guy Bartal, and Xiang Zhang.  Plasmon lasers at deep subwavelength scale.  Nature, 461(7264):629–632, October 2009.  doi:10.1038/nature08364.
  * [4] Rupert F. Oulton.  Surface plasmon lasers: sources of nanoscopic light.  Materials Today, 15(1-2):26–34, January 2012.  doi:10.1016/S1369-7021(12)70018-4.
  * [5] Ke Liu, Ning Li, Devendra K. Sadana, and Volker J. Sorger.  Integrated nanocavity plasmon light sources for on-chip optical interconnects.  ACS Photonics, 3(2):233–242, February 2016.  doi:10.1021/acsphotonics.5b00476.
  * [6] Yin Liang, Chun Li, Yong-Zhen Huang, and Qing Zhang.  Plasmonic nanolasers in on-chip light sources: prospects and challenges.  ACS Nano, 14(11):14375–14390, November 2020.  doi:10.1021/acsnano.0c07011.
  * [7] D. M. Schaadt, B. Feng, and E. T. Yu.  Enhanced semiconductor optical absorption via surface plasmon excitation in metal nanoparticles.  Applied Physics Letters, 86(6):063106, February 2005.  doi:10.1063/1.1855423.
  * [8] S. Pillai and M. A. Green.  Plasmonics for photovoltaic applications.  Solar Energy Materials and Solar Cells, 94(9):1481–1486, September 2010.  doi:10.1016/j.solmat.2010.02.046.
  * [9] Harry A. Atwater and Albert Polman.  Plasmonics for improved photovoltaic devices.  Nature Materials, 9(3):205–213, March 2010.  doi:10.1038/nmat2629.
  * [10] A. Bouhelier, M. Beversluis, A. Hartschuh, and L. Novotny.  Near-field second-harmonic generation induced by local field enhancement.  Physical Review Letters, 90(1):013903, January 2003.  doi:10.1103/PhysRevLett.90.013903.
  * [11] Johann Berthelot, Guillaume Bachelier, Mingxia Song, Padmnabh Rai, Gérard Colas des Francs, Alain Dereux, and Alexandre Bouhelier.  Silencing and enhancement of second-harmonic generation in optical gap antennas.  Optics Express, 20(10):10498–10508, May 2012.  doi:10.1364/OE.20.010498.
  * [12] Robert Czaplicki, Antti Kiviniemi, Mikko J. Huttunen, Xiaorun Zang, Timo Stolt, Ismo Vartiainen, Jérémy Butet, Markku Kuittinen, Olivier J. F. Martin, and Martti Kauranen.  Less is more: Enhancement of second-harmonic generation from metasurfaces by reduced nanoparticle density.  Nano Letters, 18(12):7709–7714, December 2018.  doi:10.1021/acs.nanolett.8b03378.
  * [13] T. Klar, M. Perner, S. Grosse, G. von Plessen, W. Spirkl, and J. Feldmann.  Surface-plasmon resonances in single metallic nanoparticles.  Physical Review Letters, 80(19):4249–4252, May 1998.  doi:10.1103/PhysRevLett.80.4249.
  * [14] Stefan A. Maier.  Plasmonics: fundamentals and applications.  Springer, New York, 2007.
  * [15] C. Ciracì, R. T. Hill, J. J. Mock, Y. Urzhumov, A. I. Fernández-Domínguez, S. A. Maier, J. B. Pendry, A. Chilkoti, and D. R. Smith.  Probing the ultimate limits of plasmonic enhancement.  Science, 337(6098):1072–1074, August 2012.  doi:10.1126/science.1224823.
  * [16] Encai Hao and George C. Schatz.  Electromagnetic fields around silver nanoparticles and dimers.  The Journal of Chemical Physics, 120(1):357, 2004.  doi:10.1063/1.1629280.
  * [17] Stephanie Dodson, Mohamed Haggui, Renaud Bachelot, Jérôme Plain, Shuzhou Li, and Qihua Xiong.  Optimizing electromagnetic hotspots in plasmonic bowtie nanoantennae.  The Journal of Physical Chemistry Letters, 4(3):496–501, February 2013.  doi:10.1021/jz302018x.
  * [18] Katrin Kneipp, Yang Wang, Harald Kneipp, Lev T. Perelman, Irving Itzkan, Ramachandra R. Dasari, and Michael S. Feld.  Single molecule detection using surface-enhanced Raman scattering (SERS).  Physical Review Letters, 78(9):1667–1670, March 1997.  doi:10.1103/PhysRevLett.78.1667.
  * [19] Shuming Nie and Steven R. Emory.  Probing single molecules and single nanoparticles by surface-enhanced Raman scattering.  Science, 275(5303):1102–1106, February 1997.  doi:10.1126/science.275.5303.1102.
  * [20] Amanda J. Haes, Christy L. Haynes, Adam D. McFarland, George C. Schatz, Richard P. Van Duyne, and Shengli Zou.  Plasmonic materials for surface-enhanced sensing and spectroscopy.  MRS Bulletin, 30(5):368–375, May 2005.  doi:10.1557/mrs2005.100.
  * [21] Younan Xia and Naomi J. Halas.  Shape-controlled synthesis and surface plasmonic properties of metallic nanostructures.  MRS Bulletin, 30(5):338–348, May 2005.  doi:10.1557/mrs2005.96.
  * [22] Arash Ahmadivand, Burak Gerislioglu, Rajeev Ahuja, and Yogendra Kumar Mishra.  Toroidal metaphotonics and metadevices.  Laser & Photonics Reviews, 14(11):1900326, November 2020.  doi:10.1002/lpor.201900326.
  * [23] Martín Caldarola, Pablo Albella, Emiliano Cortés, Mohsen Rahmani, Tyler Roschuk, Gustavo Grinblat, Rupert F. Oulton, Andrea V. Bragas, and Stefan A. Maier.  Non-plasmonic nanoantennas for surface enhanced spectroscopies with ultra-low heat conversion.  Nature Communications, 6(1):7915, August 2015.  doi:10.1038/ncomms8915.
  * [24] Arseniy I. Kuznetsov, Andrey E. Miroshnichenko, Mark L. Brongersma, Yuri S. Kivshar, and Boris Luk’yanchuk.  Optically resonant dielectric nanostructures.  Science, 354(6314):aag2472, November 2016.  doi:10.1126/science.aag2472.
  * [25] Ivano Alessandri and John R. Lombardi.  Enhanced Raman scattering with dielectrics.  Chemical Reviews, 116(24):14921–14981, December 2016.  doi:10.1021/acs.chemrev.6b00365.
  * [26] D. A. Weitz, S. Garoff, J. I. Gersten, and A. Nitzan.  A comparison of Raman scattering, resonance Raman scattering, and fluorescence from molecules adsorbed on silver island films.  Journal of Electron Spectroscopy and Related Phenomena, 29(1):363–370, January 1983.  doi:10.1016/0368-2048(83)80089-9.
  * [27] Arash Ahmadivand, Michael Semmlinger, Liangliang Dong, Burak Gerislioglu, Peter Nordlander, and Naomi J. Halas.  Toroidal dipole-enhanced third harmonic generation of deep ultraviolet light using plasmonic meta-atoms.  Nano Letters, 19(1):605–611, January 2019.  doi:10.1021/acs.nanolett.8b04798.
  * [28] F. Javier García de Abajo.  Nonlocal effects in the plasmons of strongly interacting nanoparticles, dimers, and waveguides.  The Journal of Physical Chemistry C, 112(46):17983–17987, November 2008.  doi:10.1021/jp807345h.
  * [29] Ruben Esteban, Andrei G. Borisov, Peter Nordlander, and Javier Aizpurua.  Bridging quantum and classical plasmonics with a quantum-corrected model.  Nature Communications, 3(1):825, May 2012.  doi:10.1038/ncomms1806.
  * [30] Yu Luo, A. I. Fernandez-Dominguez, Aeneas Wiener, Stefan A. Maier, and J. B. Pendry.  Surface plasmons and nonlocality: a simple model.  Physical Review Letters, 111(9), August 2013.  doi:10.1103/PhysRevLett.111.093901.
  * [31] Wenqi Zhu and Kenneth B. Crozier.  Quantum mechanical limit to plasmonic enhancement as observed by surface-enhanced Raman scattering.  Nature Communications, 5:5228, October 2014.  doi:10.1038/ncomms6228.
  * [32] Christian Girard, Aurélien Cuche, Erik Dujardin, Arnaud Arbouet, and Adnen Mlayah.  Molecular decay rate near nonlocal plasmonic particles.  Optics Letters, 40(9):2116–2119, May 2015.  doi:10.1364/OL.40.002116.
  * [33] Wenqi Zhu, Ruben Esteban, Andrei G. Borisov, Jeremy J. Baumberg, Peter Nordlander, Henri J. Lezec, Javier Aizpurua, and Kenneth B. Crozier.  Quantum mechanical effects in plasmonic structures with subnanometre gaps.  Nature Communications, 7(1):11495, June 2016.  doi:10.1038/ncomms11495.
  * [34] Zhiwen Shi, Xiaoping Hong, Hans A. Bechtel, Bo Zeng, Michael C. Martin, Kenji Watanabe, Takashi Taniguchi, Yuen-Ron Shen, and Feng Wang.  Observation of a Luttinger-liquid plasmon in metallic single-walled carbon nanotubes.  Nature Photonics, 9(8):515–519, August 2015.  doi:10.1038/nphoton.2015.123.
  * [35] G. Vecchi, V. Giannini, and J. Gómez Rivas.  Shaping the fluorescent emission by lattice resonances in plasmonic crystals of nanoantennas.  Physical Review Letters, 102(14), April 2009.  doi:10.1103/PhysRevLett.102.146807.
  * [36] S. R. K. Rodriguez, A. Abass, B. Maes, O. T. A. Janssen, G. Vecchi, and J. Gómez Rivas.  Coupling bright and dark plasmonic lattice resonances.  Physical Review X, 1(2), December 2011.  doi:10.1103/PhysRevX.1.021019.
  * [37] Danqing Wang, Ankun Yang, Weijia Wang, Yi Hua, Richard D. Schaller, George C. Schatz, and Teri W. Odom.  Band-edge engineering for controlled multi-modal nanolasing in plasmonic superlattices.  Nature Nanotechnology, 12(9):889–894, September 2017.  doi:10.1038/nnano.2017.126.
  * [38] V. G. Kravets, A. V. Kabashin, W. L. Barnes, and A. N. Grigorenko.  Plasmonic surface lattice resonances: a review of properties and applications.  Chemical Reviews, 118(12):5912–5951, June 2018.  doi:10.1021/acs.chemrev.8b00243.
  * [39] Weijia Wang, Mohammad Ramezani, Aaro I. Väkeväinen, Päivi Törmä, Jaime Gómez Rivas, and Teri W. Odom.  The rich photonic world of plasmonic nanoparticle arrays.  Materials Today, 21(3):303–314, April 2018.  doi:10.1016/j.mattod.2017.09.002.
  * [40] Charles Cherqui, Marc R. Bourgeois, Danqing Wang, and George C. Schatz.  Plasmonic surface lattice resonances: Theory and computation.  Accounts of Chemical Research, 52(9):2548–2558, September 2019\.  doi:10.1021/acs.accounts.9b00312.
  * [41] F. J. García de Abajo.  Colloquium: Light scattering by particle and hole arrays.  Reviews of Modern Physics, 79(4):1267–1290, October 2007.  doi:10.1103/RevModPhys.79.1267.
  * [42] Per Lunnemann and A. Femius Koenderink.  Dispersion of guided modes in two-dimensional split ring lattices.  Physical Review B, 90(24):245416, December 2014.  doi:10.1103/PhysRevB.90.245416.
  * [43] Anton D. Utyushev, Vadim I. Zakomirnyi, and Ilia L. Rasskazov.  Collective lattice resonances: Plasmonics and beyond.  Reviews in Physics, 6:100051, June 2021.  doi:10.1016/j.revip.2021.100051.
  * [44] Shengli Zou, Nicolas Janel, and George C. Schatz.  Silver nanoparticle array structures that produce remarkably narrow plasmon lineshapes.  The Journal of Chemical Physics, 120(23):10871, 2004.  doi:10.1063/1.1760740.
  * [45] Michael B. Ross, Chad A. Mirkin, and George C. Schatz.  Optical properties of one-, two-, and three-dimensional arrays of plasmonic nanostructures.  The Journal of Physical Chemistry C, 120(2):816–830, January 2016\.  doi:10.1021/acs.jpcc.5b10800.
  * [46] Frédéric Laux, Nicolas Bonod, and Davy Gérard.  Single emitter fluorescence enhancement with surface lattice resonances.  The Journal of Physical Chemistry C, 121(24):13280–13289, June 2017\.  doi:10.1021/acs.jpcc.7b04207.
  * [47] Fei Ding, Yuanqing Yang, Rucha A. Deshpande, and Sergey I. Bozhevolnyi.  A review of gap-surface plasmon metasurfaces: fundamentals and applications.  Nanophotonics, 7(6):1129–1156, June 2018.  doi:10.1515/nanoph-2017-0125.
  * [48] A. Christ, S. G. Tikhodeev, N. A. Gippius, J. Kuhl, and H. Giessen.  Waveguide-plasmon polaritons: strong coupling of photonic and electronic resonances in a metallic photonic crystal slab.  Physical Review Letters, 91(18):183901, October 2003.  doi:10.1103/PhysRevLett.91.183901.
  * [49] S. R. K. Rodriguez, S. Murai, M. A. Verschuuren, and J. Gómez Rivas.  Light-emitting waveguide-plasmon polaritons.  Physical Review Letters, 109(16):166803, October 2012.  doi:10.1103/PhysRevLett.109.166803.
  * [50] O. Y. Yermakov, A. I. Ovcharenko, M. Song, A. A. Bogdanov, I. V. Iorsh, and Yu. S. Kivshar.  Hybrid waves localized at hyperbolic metasurfaces.  Physical Review B, 91(23):235423, June 2015.  doi:10.1103/PhysRevB.91.235423.
  * [51] Radoslaw Kolkowski and Andriy Shevchenko.  Enabling infinite _Q_ factors in absorbing optical systems.  Nanophotonics, 12(17):3443–3454, August 2023.  doi:10.1515/nanoph-2023-0281.
  * [52] Adam Overvig and Andrea Alù.  Diffractive nonlocal metasurfaces.  Laser & Photonics Reviews, 16(8):2100633, 2022.  doi:10.1002/lpor.202100633.
  * [53] Andriy Shevchenko, Victor Ovchinnikov, and Anna Shevchenko.  Large-area nanostructured substrates for surface enhanced Raman spectroscopy.  Applied Physics Letters, 100(17):171913, April 2012.  doi:10.1063/1.4707158.
  * [54] K. T. Carron, W. Fluhr, M. Meier, A. Wokaun, and H. W. Lehmann.  Resonances of two-dimensional particle gratings in surface-enhanced Raman scattering.  Journal of the Optical Society of America B, 3(3):430–440, March 1986.  doi:10.1364/JOSAB.3.000430.
  * [55] S. G. Tikhodeev, A. L. Yablonskii, E. A. Muljarov, N. A. Gippius, and Teruya Ishihara.  Quasiguided modes and optical properties of photonic crystal slabs.  Physical Review B, 66(4):045102, July 2002.  doi:10.1103/PhysRevB.66.045102.
  * [56] Mikhail Polyanskiy.  refractiveindex.info database, 2008.
  * [57] P. B. Johnson and R. W. Christy.  Optical constants of the noble metals.  Physical Review B, 6(12):4370–4379, 1972.  doi:10.1103/PhysRevB.6.4370.
  * [58] Introduction to COMSOL Multiphysics® v. 6.0.  COMSOL AB, Stockholm, Sweden, 2022.  URL: www.comsol.com.
  * [59] Netgen/ngsolve is a high performance multiphysics finite element software. it is widely used to analyze models from solid mechanics, fluid dynamics and electromagnetics. due to its flexible python interface new physical equations and solution algorithms can be implemented easily.  URL: ngsolve.org.
  * [60] Z.S. Sacks, D.M. Kingsland, R. Lee, and Jin-Fa Lee.  A perfectly matched anisotropic absorber for use as an absorbing boundary condition.  IEEE Transactions on Antennas and Propagation, 43(12):1460–1463, Dec./1995.  doi:10.1109/8.477075.
  * [61] Stephen D. Gedney.  Introduction to the finite-difference time-domain (FDTD) method for electromagnetics.  Synthesis Lectures on Computational Electromagnetics. Springer International Publishing, Cham, 2011.  doi:10.1007/978-3-031-01712-4.
  * [62] L. Solymar and E. Shamonina.  Waves in metamaterials.  Oxford University Press, Oxford ; New York, 2009.
  * [63] V. Yannopapas and N. Stefanou.  Optical excitation of coupled waveguide-particle plasmon modes: A theoretical analysis.  Physical Review B, 69(1):012408, January 2004.  doi:10.1103/PhysRevB.69.012408.
  * [64] Dmitry Khlopin, Frédéric Laux, William P. Wardley, Jérôme Martin, Gregory A. Wurtz, Jérôme Plain, Nicolas Bonod, Anatoly V. Zayats, Wayne Dickson, and Davy Gérard.  Lattice modes and plasmonic linewidth engineering in gold and aluminum nanoparticle arrays.  Journal of the Optical Society of America B, 34(3):691–700, March 2017.  doi:10.1364/JOSAB.34.000691.


* Peng et al. (2010) Peng, Y.-j., Lilly, S. J., Kovač, K., et al. 2010, ApJ, 721, 193, doi: 10.1088/0004-637X/721/1/193
  * Planck Collaboration et al. (2011) Planck Collaboration, Abergel, A., Ade, P. A. R., et al. 2011, A&A, 536, A25, doi: 10.1051/0004-6361/201116483
  * Popescu et al. (2000) Popescu, C. C., Misiriotis, A., Kylafis, N. D., Tuffs, R. J., & Fischera, J. 2000, A&A, 362, 138.  https://arxiv.org/abs/astro-ph/0008098
  * Popescu et al. (2011) Popescu, C. C., Tuffs, R. J., Dopita, M. A., et al. 2011, A&A, 527, A109, doi: 10.1051/0004-6361/201015217
  * Schlegel et al. (1998) Schlegel, D. J., Finkbeiner, D. P., & Davis, M. 1998, ApJ, 500, 525, doi: 10.1086/305772
  * Shao et al. (2007) Shao, Z., Xiao, Q., Shen, S., et al. 2007, ApJ, 659, 1159, doi: 10.1086/511131
  * Shen et al. (2003) Shen, S., Mo, H. J., White, S. D. M., et al. 2003, MNRAS, 343, 978, doi: 10.1046/j.1365-8711.2003.06740.x
  * Simard et al. (2011) Simard, L., Mendel, J. T., Patton, D. R., Ellison, S. L., & McConnachie, A. W. 2011, ApJS, 196, 11, doi: 10.1088/0067-0049/196/1/11
  * Stevens et al. (2005) Stevens, J. A., Amure, M., & Gear, W. K. 2005, MNRAS, 357, 361, doi: 10.1111/j.1365-2966.2005.08668.x
  * Sun et al. (2021) Sun, M., Jiang, B., Zhao, H., & Ren, Y. 2021, ApJS, 256, 46, doi: 10.3847/1538-4365/ac1601
  * Tremonti et al. (2004) Tremonti, C. A., Heckman, T. M., Kauffmann, G., et al. 2004, ApJ, 613, 898, doi: 10.1086/423264
  * Tuffs et al. (2004) Tuffs, R. J., Popescu, C. C., Völk, H. J., Kylafis, N. D., & Dopita, M. A. 2004, A&A, 419, 821, doi: 10.1051/0004-6361:20035689
  * van der Giessen et al. (2022) van der Giessen, S. A., Leslie, S. K., Groves, B., et al. 2022, arXiv e-prints, arXiv:2201.10568.  https://arxiv.org/abs/2201.10568
  * Városi & Dwek (1999) Városi, F., & Dwek, E. 1999, ApJ, 523, 265, doi: 10.1086/307729
  * Vogler et al. (2005) Vogler, A., Madden, S. C., Beck, R., et al. 2005, A&A, 441, 491, doi: 10.1051/0004-6361:20042342
  * Weingartner & Draine (2001) Weingartner, J. C., & Draine, B. T. 2001, ApJ, 548, 296, doi: 10.1086/318651
  * Wild et al. (2011) Wild, V., Charlot, S., Brinchmann, J., et al. 2011, MNRAS, 417, 1760, doi: 10.1111/j.1365-2966.2011.19367.x
  * Witt & Gordon (1996) Witt, A. N., & Gordon, K. D. 1996, ApJ, 463, 681, doi: 10.1086/177282
  * Witt & Gordon (2000) —. 2000, ApJ, 528, 799, doi: 10.1086/308197
  * Witt et al. (1992) Witt, A. N., Thronson, Harley A., J., & Capuano, John M., J. 1992, ApJ, 393, 611, doi: 10.1086/171530
  * Xilouris et al. (1999) Xilouris, E. M., Byun, Y. I., Kylafis, N. D., Paleologou, E. V., & Papamastorakis, J. 1999, A&A, 344, 868.  https://arxiv.org/abs/astro-ph/9901158
  * Yip et al. (2010) Yip, C.-W., Szalay, A. S., Wyse, R. F. G., et al. 2010, ApJ, 709, 780, doi: 10.1088/0004-637X/709/2/780
  * Yuan et al. (2021) Yuan, F.-T., Lu, J., Shen, S., & Boquien, M. 2021, ApJ, 911, 145, doi: 10.3847/1538-4357/abec76
  * Zahid et al. (2017) Zahid, H. J., Kudritzki, R.-P., Conroy, C., Andrews, B., & Ho, I. T. 2017, ApJ, 847, 18, doi: 10.3847/1538-4357/aa88ae


# Most Probable Transitions
from Metastable to Oscillatory Regimes in a Carbon Cycle System

Wei Wei1,222Email<EMAIL_ADDRESS>Jianyu Hu1,111Email:
<EMAIL_ADDRESS>the corresponding author, Jianyu Chen1,333Email:
<EMAIL_ADDRESS>and Jinqiao Duan2,444Email<EMAIL_ADDRESS>
1Center for Mathematical Sciences, Huazhong University of Science and
Technology,
Wuhan, Hubei 430074, China.
2Departments of Applied Mathematics & Physics, Illinois Institute of
Technology, Chicago, IL 60616, USA.

###### Abstract

Global climate changes are related to the ocean’s store of carbon. We study a
carbonate system of the upper ocean, which has metastable and oscillatory
regimes, under small random fluctuations. We calculate the most probable
transition path via a geometric minimum action method in the context of the
large deviations theory. By examining the most probable transition paths from
metastable to oscillatory regimes for various external carbon input rates, we
find two different transition patterns, which gives us an early warning sign
for the dramatic change in the carbonate state of the ocean.

Keywords: Large deviation principle, geometric minimum action method, the most
probable transition path, carbon cycle, early warning.

Lead paragraph:

Human activities have been producing more and more carbon dioxide into the
carbon cycle, which in turn, significantly influence the climate nowadays. We
investigate an oceanic carbonate system that plays a major role in the global
carbon cycle. This system has metastable and oscillatory regimes. Due to the
small random fluctuations to the external carbon input rate, this system
undergoes a transition between these two regimes. We use a geometric minimum
action method to capture this transition phenomenon and compute the most
probable transition path. We uncover that as the external CO_2 input rate
$\nu$ changes, two different transition patterns occur with $\nu\approx 0.2$
as a critical value: (i) At a lower level external CO_2 input rate (lower than
$20\%$), the concentration of CO^2-_3 undergoes a larger excursion to shift to
the oscillatory state. (ii) But at a higher level external CO_2 input rate
(larger than $20\%$), a much smaller excursion of the concentration of CO^2-_3
leads to a transition to the oscillatory state. Moreover, as the external CO_2
input rate $\nu$ increases, the arrival value of the concentration of CO^2-_3
decreases from over 150 $\mathrm{\mu mol\cdot kg^{-1}}$ to around 50
$\mathrm{\mu mol\cdot kg^{-1}}$.

##  1 Introduction

The carbon cycle is the biogeochemical cycle in which photosynthesis converts
carbon dioxide(CO_2) to organic carbon and respiration converts organic carbon
back to CO_2. Human activities are adding more and more carbon into the
atmosphere and become an important part of the carbon cycle. Recently, due to
global warming, extreme weather has became more and more common. As a main
part of Earth’s carbon cycle, the oceanic carbonate system worths our
attention.

Rothman [17] introduced a model to describe the dynamical behaviors of the
carbonate system of the upper ocean. He took human activities and volcanic
emissions as an external source of carbon dioxide (CO_2). Using the model, he
demonstrated that as the external CO_2 inputs increases, both the dissolved
inorganic carbon $w$ and carbonate ions CO^2-_3 near a stable state may be
disrupted or excited to be near a higher oscillatory level, and eventually
come back to the stable state. These large amplitude oscillations are often
associated to climate change and mass extinctions. However, Rothman did not
take random fluctuations in the strength of the external source of CO_2. But
due to unexpected events and measurement error, it is reasonable to take such
fluctuations into consideration.

We look into the same model as Rothman, but include inevitable random
fluctuations to the external carbon input. We are interested in the bistable
case of this model. Because in original Rothman’s model, if the system starts
near the stable state, it will never reach the oscillatory state. But in the
presence of fluctuations, no matter how small it may be, in a long run, this
system will stop staying at the original stable state and travel to an
oscillatory state [4, 16]. These fluctuations make the Rothman’s deterministic
model a stochastic dynamical system.

Stochastic dynamical systems are mathematical models for complex phenomena in
physical, chemical and biological sciences [1, 5, 6, 7]. A stochastic system
may possess multistable regimes, such as one stable state and one stable limit
cycle (the oscillatory state), with a sandwiched unstable limit cycle. There
are many literature concerning the first passage time of transition from the
stable fixed point to the limited cycle [8, 9, 10, 11, 18]. The transition
behaviors are also related to early warning signs [13, 14]. We study the
transitions between stable regimes under random fluctuations [2, 3, 12]. Large
deviations theory is a useful tool that can be used to capture transition
behaviors of stochastic dynamical system. Roughly speaking, large deviations
theory measures how small the probability of rare events is. It generalized
the concept of potential of the gradient system to quasi-potential of the non-
gradient system, which is a minimizer of the action of paths connecting two
ends. We regard such a minimizer as our most probable transition path that
connects the metastable state and oscillatory state. We use the geometric
minimum action method, which was proposed by M. Heymann and E. Vanden-Eijnden
[15], to find the minimizer. It used a reparametrization method to transfer
infinity time scale to finite time scale and then numerically solved an Euler-
Lagrange equation to obtain the minimizer.

We will look into the bistable case of the stochastic carbonate system. We
will use the geometric minimum action method to find the most probable
transition path for different external CO_2 input rate. Based on the numerical
simulations, we uncover two different transition patterns.

This paper is arranged as follows. In section 2, we will briefly introduce the
carbon cycle model. The theoretical background of the geometric minimum action
method will be provided in section 3 and in section 4, we will show our
numerical experiments and explanation.

##  2 An oceanic carbonate system

In the marine carbon cycle, we look into the evolution of the carbonate system
in the upper ocean, which is illustrated in Figure 1. We consider a well mixed
open system and investigate the carbonate ions in the form of CO^2-_3, HCO_3^-
and CO_2. A main carbonate input is the dissolved CaCO_3 carried in by rivers.
Respiration is a main source of CO_2. Volcanic emissions, human activities are
regarded as an external source of CO_2. The precipitation and transportation
to lower ocean of CaCO_3 is a main carbonate output of the system.
Photosynthesis is a main export of CO_2 of the system.

Figure 1: An oceanic carbonate system

Rothman [17] introduced a model to describe dynamical behaviors of the
carbonate system and it is formulated as follows.

$\displaystyle dc$
$\displaystyle=\left[\mu\left[1-bs\left(c,c_{p}\right)-\theta\bar{s}\left(c,c_{x}\right)-\nu\right]+w-w_{0}\right]f(c)dt,$
(2.1) $\displaystyle dw$
$\displaystyle=\left[\mu\left[1-bs\left(c,c_{p}\right)+\theta\bar{s}\left(c,c_{x}\right)+\nu\right]-w+w_{0}\right]dt.$
(2.2)

Here, $w$ denotes the total dissolved inorganic carbon, that is

$w=\left[\mathrm{CO}_{2}\right]+\left[\mathrm{HCO}_{3}^{-}\right]+\left[\mathrm{CO}_{3}^{2-}\right],$
(2.3)

and $c$ is the concentration of CO^2-_3, which are expressed in units of μmol
. kg^-1. The strength of external CO_2 injection is $\nu$ and is of our main
interest in this study. Time $t$ is nondimensionalized by dividing the
homeostat’s dominant characteristic timescale $\tau_{w}$ which is about
$10^{5}$ years. Function $\bar{s}=1-s$. $s$ is a sigmoidal function, and $f$
is the “buffer function” as follow

$s\left(c,c_{p}\right)=\frac{c^{\gamma}}{c^{\gamma}+c_{p}^{\gamma}},\ \
f\left(c\right)=f_{0}\frac{c^{\beta}}{c^{\beta}+c_{f}^{\beta}}.$ (2.4)

The parameter $\mu$ is a characteristic concentration, $b$ is the maximum
CaCO_3 burial rate, $\theta$ is the maximum respiration feedback rate, and
$\nu$ is the CO_2 injection rate. The parameter $c_{x}$ is the crossover
CO^2-_3(respiration), $c_{p}$ is the crossover CO^2-_3 (burial), $c_{f}$ is
the crossover CO^2-_3 (buffering), $f_{0}$ is the maximum buffer factor, and
$\gamma,\beta$ are the sigmoid sharpness indexes. The parameters $\mu$, $b$,
$\theta$, $c_{p}$, $c_{x}$, $c_{f}$, $w_{0}$ $\gamma$ and $\beta$ are
constants and are set to fit in properties of the modern ocean. The value of
these parameters are listed in [17, SI Appendix, Table S1]. We only change the
value of $c_{x}$ and $\nu$ in this article.

System (2.1-2.2) behaves as a Hopf bifurcation when parameter $c_{x}$ varies.
For $c_{x}<55.89$ $\mathrm{\mu mol\cdot kg^{-1}}$, there is only one stable
fixed point $(c^{*},w^{*})$ of the system. For $55,89$ $\mathrm{\mu mol\cdot
kg^{-1}}<c_{x}<62.61$ $\mathrm{\mu mol\cdot kg^{-1}}$, a stable fixed point
and a stable limit cycle appear and an unstable limit cycle lies in between.
For $c_{x}>62.61$ $\mathrm{\mu mol\cdot kg^{-1}}$, there is a stable limit
cycle and an unstable fixed point. Rothman draw a bifurcation diagram for the
parameter $c_{x}$ in [17, Fig 7]. For detailed analysis and explanation of the
model, we refer readers to [17].

Figure 2: Phase-space trajectories in the bistable regime. Here
$c_{x}=57\mathrm{\mu mol\cdot kg^{-1}}$. The yellow dashed limit cycle is
unstable. Trajectories (purple arrow lines) initialized inside the unstable
limit cycle return to the stable fixed point (the reddish orange point).
Trajectories (purple arrow lines) initialized outside the unstable limit cycle
evolve to the stable limit cycle (the blue arrow loop).

In the real world, there are unexpected events that would slightly change the
carbon input rate $\nu$. For example, the Anak Krakatau volcano erupted
suddenly in April, 2020 and released massive carbon dioxide compared to a
common volcanic eruption. But this amount of carbon dioxide is still very
small compared to the total released carbon dioxide characterized by $\nu$.
The instrument error may also lead to disruption to the carbon input rate
$\nu$. So, we include a random perturbation to the parameter $\nu$ with
sufficient small amplitude $\epsilon$. This gives us the following stochastic
system:

$\displaystyle dc$
$\displaystyle=\left[\mu\left[1-bs\left(c,c_{p}\right)-\theta\bar{s}\left(c,c_{x}\right)-\nu\right]+w-w_{0}\right]f(c)dt-\epsilon\mu
f(c)dB_{t}^{1},$ (2.5) $\displaystyle dw$
$\displaystyle=\left[\mu\left[1-bs\left(c,c_{p}\right)+\theta\bar{s}\left(c,c_{x}\right)+\nu\right]-w+w_{0}\right]dt+\epsilon\mu
dB_{t}^{2},$ (2.6)

where $B_{t}^{1}$ and $B_{t}^{2}$ are two standard Brownian motions.

No matter how small the random perturbation is, the original stable state
would become a metastable state. And after a sufficient long time, the
solution to (2.5-2.6) would transport from one metastable state to another
metastable state. That means, for the parameter $c_{x}$ lying in the bistable
regime $(55,89$ $\mathrm{\mu mol\cdot kg^{-1}},62.61$ $\mathrm{\mu mol\cdot
kg^{-1}})$, the solution to system (2.5-2.6) will go from the stable state
$(c^{*},w^{*})$ to the oscillatory state, that is the stable limit cycle of
the deterministic system (2.1-2.2).

We will focus on the bistable case, that is $55,89$ $\mathrm{\mu mol\cdot
kg^{-1}}<c_{x}<62.61$ $\mathrm{\mu mol\cdot kg^{-1}}$. In the next section, we
will introduce a large deviation approach that enables us to capture this
transition phenomena.

##  3 A large deviation approach

We will briefly introduce the large deviations theory and the geometric
minimum action method based on it.

We consider a generalized stochastic differential equation in $\mathbb{R}^{2}$
of equation (2.5-2.6) as follows,

$dX^{\epsilon}(t)=\kappa(X^{\epsilon}(t))dt+\epsilon\eta(X^{\epsilon}(t))dB(t),\quad
X^{\epsilon}(0)=x_{0},$ (3.1)

where $\kappa:\mathbb{R}^{2}\rightarrow\mathbb{R}^{2}$ is a regular function,
$\eta:\mathbb{R}^{2}\rightarrow\mathbb{R}^{2\times 2}$ is a $2\times 2$
matrix-valued function and $B$ is a standard Brownian motion in
$\mathbb{R}^{2}$.

Under some assumptions on the coefficients $\kappa$ and $\eta$, for instance,
the Lipschitz continuous conditions, the solution $X^{\epsilon}$ to equation
(3.1) satisfies a large deviation principle with an action functional $S_{T}$.
It enables us to estimate the probability that $X^{\epsilon}$ stays in a
$\delta$-neighborhood of a path $\psi$ though the following asymptotic
relation.

$\mathbb{P}(d(X^{\epsilon},\psi)<\delta)\sim\exp\left(-\frac{1}{\epsilon}S_{T}(\psi)\right),\quad\textrm{as
}\epsilon\to\infty,$ (3.2)

where $d$ is the distance in $\mathbb{R}^{2}$. The action functional $S_{T}$
is

$S_{T}(\varphi)=\int_{0}^{T}L(\varphi,\dot{\varphi})dt,$ (3.3)

with the Lagrange functional

$L(\varphi,\dot{\varphi})=\frac{1}{2}\left(\dot{\varphi}-\kappa(\varphi)\right)(\eta\eta^{T}(\varphi))^{-1}[\dot{\varphi}-\kappa(\varphi)].$

The large deviation theory also equips us to capture the long time behaviors
of $X^{\epsilon}$ with quasi-potential $V$. That is

$V(x_{1},x_{2})=\inf_{T>0}\inf_{\varphi\in\bar{C}_{x_{1}}^{x_{2}}(0,T)}S_{T}(\varphi),$
(3.4)

where $\bar{C}_{x_{1}}^{x_{2}}(0,T)$ is the space of all absolutely continuous
functions that start at $x_{2}$ and end at $x_{2}$. Assume that there exists a
path $\tilde{\varphi}$ and time $\tilde{T}\in[0,\infty]$ satisfying,

$V(x_{1},x_{2})=S_{\tilde{T}}(\tilde{\varphi}).$

Then, in probability $1$, the sample paths of $X^{\epsilon}$ from $x_{1}$ to a
$\delta$-neighborhood of $x_{2}$, converge to $\tilde{\varphi}$ when the noise
intensity $\epsilon$ and the neighborhood size $\delta$ tend to $0$ as well as
the time $T$ tends to $\tilde{T}$. This result means that the sample path
$\tilde{\varphi}$ is the most probable transition path we are looking for.
Readers may refer to [15, Proposition 2.3] for a more precise statement.

In summary, we only need to find the minimizer of equation (3.4) to obtain a
most probable transition path between two metastable states. For a most
probable transition path between a metastable state and an oscillatory state
$D$, for example, a limit cycle, we need to find the minimizer of the
following equation.

$V(x,D)=\inf_{y\in D}V(x,y)=\inf_{y\in
D}\inf_{T>0}\inf_{\varphi\in\bar{C}_{x}^{y}(0,T)}S_{T}(\varphi).$ (3.5)

We will conduct a geometric minimum action method introduced in [15] to find
the most probable transition path between the stable fixed point and points
near the limit cycle. Then, we pick up the path which has the minimum quasi-
potential as the most probable transition path between the stable state and
the limit cycle. That is the minimizer to equation (3.5).

##  4 Numerical simulations and discussions

In this section, we will present the most probable transition paths that are
computed through the geometric minimum action method. We are interested in how
the most probable transition path changes as the external CO_2 input rate
$\nu$ varies, in order to better understand this carbonate system.

We choose $c_{x}=62$ and use the Euler scheme to simulate the limit cycle
value, that is the oscillatory state of the carbonate system. We calculate the
most probable transition path from the stable state $(c^{*},w^{*})$ to a point
near the oscillatory state, for $\nu$ from $0$ to $0.9$. $3000$ points on the
transition path are made to be equidistant. The quasi-potential between two
points on the limit cycle is zero, because there is a deterministic orbit
connecting each other though the limit cycle and this makes the action achieve
its minimum $0$. In practice, as we do not actually choose points on the limit
cycle, we can obtain a most probable transition path though the geometric
minimum action method. Then, we only need to find the path with the least
quasi-potential as the most probable transition path connecting the stable
state $(c^{*},w^{*})$ and the oscillatory state. The most probable transition
paths for $\nu=0,0.19,0.2$ and $0.4$ are shown in Figure 3 in the $c$-$w$
plane.

The most probable transition path has to travel a longer way to achieve the
oscillatory state for $\nu<0.2$ as is shown in Figure 4. The dramatic change
happens near $\nu=0.2$. The most probable transition path stops ending at the
“southeast” part of the limit cycle and starts to end at the “southwest” part
of the limit cycle as is shown in Figure 3. Consequently, the oscillatory
ranges for $c$ and $w$ on the most probable transition path shrink. And the
arrival value of concentration of CO^2-_3 decreases from above 150
$\mathrm{\mu mol\cdot kg^{-1}}$ to around 50 $\mathrm{\mu mol\cdot kg^{-1}}$.

(a) $\nu=0$

(b) $\nu=0.19$

(c)  $\nu=0.2$

(d) $\nu=0.4$

Figure 3: Most probable transition paths ( reddish orange arrow lines) from
the metastable state $(c^{*},w^{*})$ (orange point) to the oscillatory state
(blue limit cycle) for external CO_2 input rate $\nu=0,0.19,0.20.4$ and
parameter $c_{x}=62$.

(a) The length for the most probable transition path

(b) The arrival value of $c$

Figure 4: The length for the most probable transition path and the arrival
value of concentration of CO^2-_3, $c$, with respect to different external
CO_2 input rate $\nu$. Parameter $c_{x}=62.$ Critical changes happen around
$\nu=0.2$. The length of most probable transition path decreases to lower than
400. The arrival value of $c$ decreases from 150 $\mathrm{\mu mol\cdot
kg^{-1}}$ to around 50 $\mathrm{\mu mol\cdot kg^{-1}}$.

(a)  Transition in view of $c$ for $\nu=0.19$

(b)  Transition in view of $c$ for $\nu=0.2$

Figure 5: Illustrations for transitions in view of $c$, the concentration of
CO^2-_3, for CO_2 input rate $\nu=0.19,0.2$ and parameter $c_{x}=62$. The
orange line is the value of $c$ on the most probable transition path. The
green line represents the value of $c$ near the metastable state and the blue
line represents the value of $c$ on the oscillatory state.

If we look into how $c$, the concentration of CO^2-_3, changes with respect to
time, we find two different patterns for transitions. Again, we use Euler
scheme to generate the data of $c$ before and after transition for noise
strength $\epsilon=0.01$. Then we connect them with the most probable
transition path and illustrate the results in in Figure 4. When the the
external CO_2 rate $\nu$ is smaller than a threshold near $0.2$, the
concentration of CO^2-_3 will undergo an excursion, whose size is slightly
smaller than the oscillatory state (the limit cycle shown in Figure 3, to
shift to the oscillatory state. When $\nu$ exceeds the threshold near $0.2$,
the concentration of CO^2-_3 will undergo a much smaller excursion to shift to
the oscillatory state. This result implies that at more that $20\%$ external
carbon input rate, much less changes to the concentration of CO^2-_3 may lead
to a transition.

Our calculation shows the early warning sign of the carbonate system. If the
solution $(c,w)$ stays near the stable state $(c^{*},w^{*})$, the carbonate
system of the ocean and the carbon cycle are in dynamic equilibrium. But when
the solution $(c,w)$ enters the oscillatory state (the stable limit cycle),
the carbonate system of the ocean will change dramatically and may cause
massive extinction. It is important to identify the transition between these
two states. Our calculation shows two transition patterns of $c$, the
concentration of CO^2-_3(See Figure 5): (i) At a lower level external CO_2
input rate (lower than $20\%$), $c$ will first drop to lower than 50
$\mathrm{\mu mol\cdot kg^{-1}}$ and then increase to higher than 150
$\mathrm{\mu mol\cdot kg^{-1}}$, before transferring to the oscillatory state.
(ii) But at a higher level external CO_2 input rate (larger than $20\%$), $c$
will only increase a little to around to 100 $\mathrm{\mu mol\cdot kg^{-1}}$
and then drop to around 50 $\mathrm{\mu mol\cdot kg^{-1}}$ to enter the
oscillatory state. Given the external CO_2 input rate, we can identify the
transition before the system enters the oscillatory state.

In conclusion, we have computed the most probable transition path connecting a
metastable state and an oscillatory state for a carbonate system under random
fluctuations in the external CO_2 input rate. As the external CO_2 input rate
$\nu$ changes, two different transition patterns occur with $\nu\approx 0.2$
as a critical value (See Figure 5): (i) At a lower level external CO_2 input
rate (lower than $20\%$), the concentration of CO^2-_3 undergoes a larger
excursion to shift to the oscillatory state. (ii) But at a higher level
external CO_2 input rate (larger than $20\%$), a much smaller excursion of the
concentration of CO^2-_3 leads to a transition to the oscillatory state.
Moreover, as the external CO_2 input rate $\nu$ increases, the arrival value
of the concentration of CO^2-_3 decreases from over 150 $\mathrm{\mu mol\cdot
kg^{-1}}$ to around 50 $\mathrm{\mu mol\cdot kg^{-1}}$.(Figure 4b).

## Acknowledgements

This work was partly supported by the NSFC grants 11771449.

## Data Availability

The data that support the findings of this study are openly available in
GitHub.
https://github.com/JayWeiess/CarbonCycle-gMAM

## References

  * [1] L. Arnold.  Random Dynamical Systems.  Springer, Berlin, 2003.
  * [2] J. L. Creaser, K. Tsaneva-Atanasova, and P. Ashwin.  Sequential noise-induced escapes for oscillatory network dynamics.  SIAM Journal on Applied Dynamical Systems, 17:500–525, 2018.
  * [3] R. de la Cruz, R. Perez-Carrasco, P. Guerrero, T. Alarcon, and K. M. Page.  Minimum action path theory reveals the details of stochastic transitions out of oscillatory states.  Physical Review Letters, 120:128102, Mar 2018.
  * [4] P. D. Ditlevsen.  Observation of $\alpha$-stable noise induced millennial climate changes from an ice-core record.  Geophysical Research Letters, 26(10):1441–1444, 1999.
  * [5] J. Duan.  An Introduction to Stochastic Dynamics.  Cambridge University Press, 2015.
  * [6] J. Duan and W. Wei.  Effective Dynamics of Stochastic Partial Differential Equations.  Elsevier, 2014.
  * [7] P. Imkeller and J.-S. Von Storch.  Stochastic Climate Models, volume 49.  Birkhäuser, 2012.
  * [8] Y. Li, R. Mei, Y. Xu, J. Kurths, J. Duan, and R. Metzler.  Particle dynamics and transport enhancement in a confined channel with position-dependent diffusivity.  New Journal of Physics, 22(5):053016, 2020.
  * [9] Y. Li, Y. Xu, and J. Kurths.  First-passage-time distribution in a moving parabolic potential with spatial roughness.  Physical Review E, 99(5):052203, 2019.
  * [10] Y. Li, Y. Xu, J. Kurths, and J. Duan.  The influences of correlated spatially random perturbations on first passage time in a linear-cubic potential.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 29(10):101102, 2019.
  * [11] Y. Li, Y. Xu, J. Kurths, and X. Yue.  Lévy-noise-induced transport in a rough triple-well potential.  Physical Review E, 94(4):042222, 2016.
  * [12] L. Lin, H. Yu, and X. Zhou.  Quasi-potential calculation and minimum action method for limit cycle.  Journal of Nonlinear Science, 29:961–991, 2019.
  * [13] J. Ma, Y. Xu, Y. Li, R. Tian, G. Chen, and J. Kurths.  Precursor criteria for noise-induced critical transitions in multi-stable systems.  Nonlinear Dynamics, 101(1):21–35, 2020.
  * [14] J. Ma, Y. Xu, Y. Li, R. Tian, and J. Kurths.  Predicting noise-induced critical transitions in bistable systems.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 29(8):081102, 2019.
  * [15] M.Heymann and E. Vanden-Eijnden.  The geometric minimum action method: A least action principle on the space of curves.  Communications on Pure and Applied Mathematics, LXI:1052–1117, 2008\.
  * [16] B. Qiu and W. Miao.  Kuroshio path variations south of Japan: Bimodality as a self-sustained internal oscillation.  Journal of Physical Oceanography, 30(8):2124–2137, 2000.
  * [17] D. H. Rothman.  Characteristic disruptions of an excitable carbon cycle.  Proceedings of the National Academy of Sciences, 116(30):14813–14822, 2019.
  * [18] Z. Wang, Y. Xu, and H. Yang.  Lévy noise induced stochastic resonance in an fhn model.  Science China Technological Sciences, 59(3):371–375, 2016.


††thanks: Contributed equally to this work.††thanks: Contributed equally to
this work.

# Twist-dependent intra- and interlayer excitons in moiré MoSe2 homobilayers

Viviana Villafañe  Malte Kremser  Walter Schottky Institut and Physik
Department, Technische Universität München, Am Coulombwall 4, 85748 Garching,
Germany  Ruven Hübner  Institut für Theoretische Physik, Universität Bremen,
P.O. Box 330 440, 28334 Bremen, Germany  Marko M. Petrić  Walter Schottky
Institut and Department of Electrical and Computer Engineering, Technische
Universität München, Am Coulombwall 4, 85748 Garching, Germany  Nathan P.
Wilson  Andreas V. Stier  Walter Schottky Institut and Physik Department,
Technische Universität München, Am Coulombwall 4, 85748 Garching, Germany  Kai
Müller  Walter Schottky Institut and Department of Electrical and Computer
Engineering, Technische Universität München, Am Coulombwall 4, 85748 Garching,
Germany  Matthias Florian  University of Michigan, Dept. of Electrical
Engineering and Computer Science, 48109 Ann Arbor, MI, USA  Alexander
Steinhoff  Institut für Theoretische Physik, Universität Bremen, P.O. Box 330
440, 28334 Bremen, Germany  Jonathan J. Finley<EMAIL_ADDRESS>Walter
Schottky Institut and Physik Department, Technische Universität München, Am
Coulombwall 4, 85748 Garching, Germany

###### Abstract

Optoelectronic properties of van der Waals homostructures can be selectively
engineered by the relative twist angle between layers. Here, we study the
twist-dependent moiré coupling in MoSe2 homobilayers. For small angles, we
find a pronounced redshift of the K-K and $\Gamma$-K excitons accompanied by a
transition from K-K to $\Gamma$-K emission. Both effects can be traced back to
the underlying moiré pattern in the MoSe2 homobilayers, as confirmed by our
low-energy continuum model for different moiré excitons. We identify two
distinct intralayer moiré excitons for R-stacking, while H-stacking yields two
degenerate intralayer excitons due to inversion symmetry. In both cases,
bright interlayer excitons are found at higher energies. The performed
calculations are in excellent agreement with experiment and allow us to
characterize the observed exciton resonances, providing insight about the
layer composition and relevant stacking configuration of different moiré
exciton species.

Van der Waals (vdW) homo- and heterostructures formed from monolayer
transition metal dichalcogenides (TMDs) are unique semiconductor systems in
which light couples to electronic and spin excitations with potential for
novel optoelectronic and valleytronic applications Geim and Grigorieva (2013);
Liu et al. (2016). In angle-aligned TMD heterostructures, correlated
insulating states have been shown to emerge for different fractional charge
fillings of moiré superlattice sites, enabling investigation of quantum many-
body states Tang et al. (2020); Regan et al. (2020); Xu et al. (2020); Huang
et al. (2021); Jin et al. (2021); Li et al. (2021) with the potential for
optical measurement and coherent control of strongly-correlated phases Jin et
al. (2019); Seyler et al. (2019); Tran et al. (2019); Alexeev et al. (2019);
Baek et al. (2020); Wang et al. (2020). Twisted homostructures not only allow
for controlled tuning of the underlying moiré superlattice potential but also
enhanced formation of hybridized minibands due to the absence of lattice and
energy mismatch in the constitutive monolayers Plochocka (2020). Until now,
few studies have appeared on the angle-dependent optical and electronic
properties of twisted $\mathrm{MoSe}_{2}$ homobilayers including
photoluminescence (PL) experiments exhibiting moiré-trapped trions for twist
angles close to $None$ and $None$ Marcellina et al. (2021); and a static
electric dipole moment characterization of different excitonic species on a
single $\mathrm{MoSe}_{2}$ bilayer with a fixed twist angle of $None$ Sung et
al. (2020). Recent reports on $\mathrm{WS}_{2}$ Yan et al. (2019) and
$\mathrm{MoS}_{2}$ Van Der Zande et al. (2014); Liu et al. (2014) showed that
the K-K exciton transition Yan et al. (2019); Van Der Zande et al. (2014); Liu
et al. (2014) is insensitive to twist angle.

We combine optical spectroscopy on hBN-encapsulated $\mathrm{MoSe}_{2}$
homobilayers with theory to obtain new information about the twist angle
dependent optoelectronic response. Both, direct K-K and indirect $\Gamma$-K
excitons, exhibit an abrupt decrease of the emission energy in the vicinity of
$None$ or $None$. This rapid change of exciton energy is accompanied by the
appearance of indirect exciton PL below $None$ and above $None$ and vanishing
direct exciton PL. Our theoretical predictions based on an ab-initio-based
continuum model are in excellent agreement with experiment, showing that the
transition between different regimes of emission energy can be understood in
terms of increasing exciton localization in moiré sites at small twist angles.
Furthermore, our model allows us to characterize the observed exciton
resonances, providing insight into their layer composition and resulting
binding energies.

Figure 1:  Photoluminescence (a) and reflectivity (b) spectra for
$\mathrm{MoSe_{2}}$ homobilayers with different twist angles. The energies are
given as the detuning with respect to the neutral exciton in the monolayer
region of each sample, labeled by dotted lines. Black spectra in both panels
show an exemplary monolayer signal. (c) Energy shift between monolayer and
bilayer regions for excitons and trions extracted from the data presented in
panels (a) and (b). Green (blue) data points represent the neutral exciton
shift extracted from reflectivity (PL) measurements. Red data points represent
trion shifts extracted from PL measurements. Error bands for the extracted
spectral positions are respectively color coded with the reflectivity
measurements error band being extended across all panels for comparison.
Negative shifts correspond to an energy reduction in the bilayer region with
respect to the monolayer region. Indirect exciton (Xind) PL spectra (d) and
extracted energies (e) as a function of twist angle of the MoSe2 homobilayers.

Samples were prepared using a tear-and-stack technique Kim et al. (2016)
combined with a modified version of the hot-pick-up method Wang et al. (2013);
Pizzocchero et al. (2016) to assemble twisted MoSe2 homobilayers with relative
stacking angle $\Delta\theta$ (see the Supplemental Material (SM)). Fig. 1(a)
presents optical PL emission from $\mathrm{MoSe_{2}}$ twisted bilayers in the
range
$$0\text{\,}\mathrm{\SIUnitSymbolDegree}$\leq\Delta\theta\leq$60\text{\,}\mathrm{\SIUnitSymbolDegree}$$,
probed at $10\text{\,}\mathrm{K}$ using $500\text{\,}\mathrm{n}\mathrm{W}$ of
continuous-wave $532\text{\,}\mathrm{n}\mathrm{m}$ excitation laser focused
onto a diffraction-limited spot (100$\times$ objective, NA=0.7). PL spectra in
Fig. 1(a) are plotted as a function of neutral K-K exciton ($\mathrm{X}^{0}$)
detuning between homo- and monolayer regions and normalized with respect to
$\mathrm{X}^{0}$ intensity measured in the respective monolayers. We observe
that $\mathrm{X}^{0}$ and trion ($\mathrm{X}^{\mathrm{T}}$) energies are
redshifted relative to the monolayer $\mathrm{X}^{0}$. One striking feature of
the data presented in Fig. 1(a) is the complete absence of $\mathrm{X}^{0}$
and $\mathrm{X}^{\mathrm{T}}$ emission for bilayers having twist angles in the
ranges $None$ to $None$ and $None$ to $None$.

Figure 2:  (a) Top panel: Moiré pattern formed by superposition of two
hexagonal lattices twisted by $\Delta\theta=$$$ showing alternating regions of
AABB, AB and BA stacking order. Bottom panel: Side view of the local stacking
configurations marked in the left panel. (b) DFT band structure for AB
stacking including spin-orbit coupling. The green (orange) arrow indicates the
studied direct (indirect) transitions. (c) Variation of the conduction band at
the K-point as well as the valence bands at the K\- and $\Gamma$-point as a
function of interlayer translation. Up (down) triangles correspond to spin up
(down) bands obtained from individual DFT calculations for stacking
configurations defined in (a). The color encodes the layer contribution to
Bloch states. The solid (dashed) lines show the result of our fit model for
the variation of spin up (down) bands. (d) Moiré potentials $\Delta
U(\textbf{r})$ for the indirect exciton $\mathrm{X_{ind}}$ and the direct
exciton $\mathrm{X_{1}}$ shown relative to the potential minimum within the
moiré unit cell.

Fig. 1(b) presents differential reflectivity measurements for the same series
of samples discussed in Fig. 1(a). We note that the reflectivity data reveals
clear signatures of the $\mathrm{X}^{0}$ transition from samples having small
twist angles ($None$ to $None$ and $None$ to $None$) which could not be
observed in PL experiments due to the low quantum yield of these transitions
Frisenda et al. (2017); Evans and Hazelwood (1971); Anedda and Fortin (1980).
To analyze the reflectivity spectra, we modeled the refractive index of the
MoSe2 monolayer as a single Lorentz oscillator, except in the case of $None$
where we used two oscillators, and calculated the differential reflectivity
using the transfer-matrix method (see the SM). The simulated reflectivity
curves are shown as diffuse color lines in Fig. 1(b) with good agreement
between experiment and theory. The reflectivity of $None$, $None$ and $None$
show a double-peak structure predicted by our theoretical model steaming from
different hybridized K-K excitons, as we will see below. We evaluated the
average energies of $\mathrm{X}^{0}$ and $\mathrm{X}^{\mathrm{T}}$, as
obtained from PL and reflectivity measurements. In Fig. 1(c), we present the
energy shift as a function of the stacking angle, relative to the energy of
the corresponding excitonic transition ($\mathrm{X}^{0}$ or
$\mathrm{X}^{\mathrm{T}}$) in the monolayer. Green data points within the
green error band in the uppermost panel in Fig. 1(c) denote relative
$\mathrm{X}^{\mathrm{0}}$ energies measured via differential reflectivity.
Blue and red data points in the middle and bottommost panel show the energy
shift of $\mathrm{X}^{0}$ and $\mathrm{X}^{\mathrm{T}}$ respectively,
determined from PL measurements. At large twist angles (away from $None$ and
$None$), $\mathrm{X}^{0}$ and $\mathrm{X}^{\mathrm{T}}$ exhibit similar
behavior: Both are redshifted by
$\sim$$20\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$ with respect to the
monolayer emission with only a minor twist angle dependence. The observed
energetic reduction for $\mathrm{X}^{0}$ and $\mathrm{X}^{\mathrm{T}}$
transitions in the homobilayer region can be explained considering two
fundamentally different effects: (i) hybridization of electronic states
between layers Li and Galli (2007); Splendiani et al. (2010) and (ii) static
dielectric screening induced by hBN encapsulation and the mutual proximity of
the second MoSe2 monolayer Ugeda et al. (2014); Raja et al. (2017); Florian et
al. (2018); Cho and Berkelbach (2018). In fact, $\mathrm{X}^{0}$ energy has a
strong sensitivity to its dielectric environment enhanced by the bilayers
reduced dimensionality and dielectric constant mismatch between hBN and MoSe2
Segura et al. (2018). Consequently, for stacking angles around
$\Delta\theta\simeq$ $None$, static screening is the dominant effect inducing
a constant redshift in the $\mathrm{X}^{0}$ and $\mathrm{X}^{\mathrm{T}}$
energies Cho and Berkelbach (2018); Lin et al. (2014); Stier et al. (2016);
Florian et al. (2018). Remarkably, as twist angles of $None$ or $None$ are
approached the redshift increases to $\sim$40\text{\,}\mathrm{meV}$$. In this
case, hybridization effects between the layers are dominant and twist-angle
dependent. Figs. 1(d) and (e) present measured PL spectra and energies for
indirect $\Gamma$-K excitons (Xind), appearing $200\text{\,}\mathrm{meV}$
below the observed $\mathrm{X}^{0}$ transition. We observe a strong energy
dependence as a function of stacking angle, similarly as the one found for
$\mathrm{X}^{0}$ and $\mathrm{X}^{\mathrm{T}}$. For twist angles in the ranges
$None$ to $None$ and $None$ to $None$, the decrease in the intensity of the
$\mathrm{X}^{0}$ PL emission is accompanied by an increase on PL emission into
these Xind momentum indirect states at lower energies Sung et al. (2020). As
we now continue to show, the experimentally-observed $\mathrm{X}^{0}$ and
$\mathrm{X}^{\mathrm{T}}$ energetic shifts and Xind PL emission increase for
low stacking angles originate from moiré physics, a behaviour predicted by our
low-energy continuum model.

For small stacking angles, the moiré pattern can be seen as a smooth variation
of the stacking configurations known from untwisted bilayers (see Fig. 2(a)).
In this case, local stacking configurations found on a path across the
supercell can be approximated by translations of the top layer with respect to
the bottom layer, starting from the AABB configuration and passing through
stacking configurations AB and BA. We use the generic term R-stacking for all
configurations of this type. Likewise, for the near-$None$ case (H-stacking),
we introduce a translation through AA, ABBA and BB stacking. Contrary to
rotations, translations do not change the periodicity of the bilayer so that
density functional theory (DFT) calculations are readily performed for the
different stacking configurations. For both, R- and H-stacking, we choose six
translations corresponding to six equidistant points within one supercell. We
use a fixed lattice constant of
$a=$3.336\text{\,}\mathrm{\SIUnitSymbolAngstrom}$$, which accounts for
$1.4\text{\,}\mathrm{\char 37\relax}$ of strain compared to a naturally-grown
homobilayer. Similar values have already been proposed for MoSe2 homobilayers
Sung et al. (2020) and are consistent with strain maxima of up to
$2.1\text{\,}\mathrm{\char 37\relax}$ recently calculated for twisted
homobilayers Linderälv et al. (2022) (see the SM). Fig. 2(b) introduces the
lowest bright K-K transitions $\mathrm{X_{1}}$, $\mathrm{X_{2}}$,
$\mathrm{X_{3}}$, and $\mathrm{X_{4}}$ and the lowest indirect $\Gamma$-K
exciton ($\mathrm{X_{ind}}$) calculated for AB stacking including spin-orbit
coupling. The results for the twist-dependent $\Gamma$\- and K-point energies
are presented in Fig. 2(c). We find that electron localization at the K-point
depends strongly on the stacking configuration. While for high-symmetry
configurations AABB we find perfect electronic delocalization over both
layers, electrons are well localized in one of the two layers for the most
stable configurations AB and BA. However, we note that this phenomenon is
exclusive for R-stacking and limited by spin-layer polarization in case of
H-stacking. For the latter, well-localized states are found for all
configurations (see the SM). Compared to the K-point, the $\Gamma$-point
energy shows a much stronger stacking dependence. This modulation is directly
related to the varying interlayer distance for different stacking
configurations, which has particularly high influence on the $\Gamma$-point
due to strong interlayer hopping. Moreover, interlayer hopping leads to a
pronounced spread of the electronic states over both layers as indicated by
the constant light gray color of the corresponding data points in Fig. 2(c).
Even though there are several indirect excitons involving higher conduction
bands at the K-point, the variation of the $\Gamma$-point dominates the
corresponding moiré potential, which in turn leads to negligible variations
between the different indirect excitons.

To calculate the individual moiré potentials, we model the band variations of
Fig. 2(c) by symmetry considerations using a similar approach to Refs. Wang et
al. (2017); Hagel et al. (2021) (see the SM). Thereby, we have direct access
to a set of periodic band gap variations $U(\boldsymbol{r})$, each defined by
two high-symmetry points and corresponding bands. Excitons inside the moiré
superlattice experience this energetic modulation effectively as the moiré
potential, which is visualized for the exciton species $\mathrm{X_{ind}}$ and
$\mathrm{X_{1}}$ in Fig. 2(d). In this picture, the Hamiltonian for the
center-of-mass (COM) motion of the exciton is given by Wu et al. (2018)

$\hat{H}_{0}=E^{0}_{\mathrm{gap}}+\frac{\hbar^{2}\boldsymbol{\hat{Q}}^{2}}{2M}+U(\boldsymbol{\hat{r}})\,,$
(1)

where $E^{0}_{\mathrm{gap}}$ is the local band gap corresponding to
$\boldsymbol{r}=0$ and $M$ is the total effective mass of electron and hole
combined. $\boldsymbol{Q}$ and $\boldsymbol{r}$ are the COM wavevector and
position vector, respectively. In our model, the twist angle only enters via a
scaling factor inside the moiré potential $U(\boldsymbol{r})$ acting on the
argument $\boldsymbol{r}$. This dependence results from describing every point
inside the supercell by its local and untwisted band structure configuration.
Thus, for the near-$None$ case, we assume the moiré supercell to scale with
$1/\sin(\theta/2)$ and for the near-$None$ case with
$1/\sin((\frac{\pi}{3}-\theta)/2)$ Moon and Koshino (2013). Outside the moiré
regime, approximately between $None$ and $None$, our calculations can be
interpreted as an extrapolation of the moiré case. Even though this
approximation is not strictly valid for large twist angles in the vicinity of
$30\text{\,}\mathrm{\SIUnitSymbolDegree}$, our calculations are useful to
identify differences between R- and H-stacking. The exciton COM motion is
accompanied by a relative motion of electron and hole described by the Wannier
equation, which accounts for the screened electron-hole Coulomb potential. The
latter depends on the actual charge distribution in stacking direction and is
therefore calculated separately for the case of intra- and interlayer
excitons.Wang et al. (2019) Then, the admixture of intra- and interlayer
character is taken into account according to the layer contributions of
valence and conduction band given by the colormap in Fig. 2(c) (see the SM).

Figure 3:  (a) Angle-dependent exciton energies. Solid lines show the result
of our continuum model based on DFT calculations. Each line corresponds to an
individual exciton species, while the color encodes the layer character. The
intermediate case (gray color) indicates delocalization of electron and/or
hole over both layers. The green and orange data points are taken from Fig.
1(c) (DR) and Fig. 1(e) (PL), respectively. (b) Integrated PL intensity for
$\mathrm{X}^{0}$ and $\mathrm{X_{ind}}$ relative to the corresponding
monolayer intensity as a function of twist angle.

The main findings of our work are presented in Fig. 3. Panel (a) compares the
measured ($\mathrm{X}^{0}$, $\mathrm{X_{ind}}$, data points) and calculated
($\mathrm{X_{1,2,3,4}}$, $\mathrm{X_{ind}}$, solid lines) transition energies
as a function of twist angle. The theoretical curves are uniformly shifted by
$554\text{\,}\mathrm{meV}$ to match the experimental data point of the
indirect exciton for small twist angles, which is in good agreement with
expected GW corrections for $\mathrm{MoSe_{2}}$ homobilayers He et al. (2014)
(see the SM). The energy discontinuity that can be seen in Fig. 3(a) around
$None$ arises from the delicate distinction between intra- and interlayer
excitons, which is clearer for H-stacking due to the aforementioned spin-layer
polarization. This effect leads to a strong intralayer exciton binding energy
and thereby results in a slight asymmetry when comparing the energies
approaching from R- and H-stacking configurations at $None$. Since the layer
contributions to each exciton directly influence the binding energies, we
visualized them separately via the line colors. We find that for both, R- and
H-stacking, the two lower direct excitons $\mathrm{X}_{2}$ and
$\mathrm{X}_{3}$ have mostly intralayer character, while interlayer excitons
$\mathrm{X}_{1}$ and $\mathrm{X}_{4}$ have higher energies due to a reduced
Coulomb interaction. In case of R-stacking, intralayer excitons
$\mathrm{X}_{2}$ and $\mathrm{X}_{3}$ are split by $12\text{\,}\mathrm{meV}$,
closely followed by the lowest interlayer exciton $\mathrm{X}_{1}$. Thereby,
$\mathrm{X}_{1}$, $\mathrm{X}_{2}$ and $\mathrm{X}_{3}$ are in agreement with
the experimental findings of Ref. Sung et al. (2020) for angles near $None$.
For H-stacking, inversion symmetry leads to double degeneracy of intra- and
interlayer excitons. The model explains our experimental findings, and we
could assign the measured $\mathrm{X}^{0}$ to $\mathrm{X}_{2}$. However, the
reflectivity data presented in Fig. 1(b) for $None$, $None$ and $None$ samples
hints at a double peak structure of $\mathrm{X}^{0}$ that might originate from
the R-type $\mathrm{X}_{2}$ and $\mathrm{X}_{3}$ intralayer excitons. In that
case, the measured $\mathrm{X}^{0}$ could arise from an hybridized mixture of
these two species Ruiz-Tijerina and Fal’ko (2019). For small stacking angles,
our model predicts a shift of $17\text{\,}\mathrm{meV}$ spanned by both
exciton species, a value that matches very well the experimental shift
obtained for the $\mathrm{X}^{0}$ resonance of $18\text{\,}\mathrm{meV}$. The
small discrepancy between theory and experiment may be due to lattice
reconstruction effects Weston et al. (2020) that are not taken into account in
our models.

As depicted in Fig. 3(a), the indirect $\Gamma$-K transition
$\mathrm{X_{ind}}$ possesses a very strong angle dependence for small stacking
angles, resulting from the interaction with the deep moiré potential shown in
the left panel in Fig. 2(d). Towards $None$ the moiré supercell becomes
incommensurately small and the interlayer distance in the homobilayer
increases. In this regime, the kinetic energy dominates the COM Hamiltonian
(1), resulting in a quasi-homogeneous spreading of the wave function across
the whole superlattice and therefore a weak angular dependence. The critical
angle for delocalization is smaller for direct excitons
($\approx$3\text{\,}\mathrm{\SIUnitSymbolDegree}$$) and comparable to reported
findings in twisted heterobilayers Brem et al. (2020).

Fig. 3(b) presents the integrated experimental PL intensity for
$\mathrm{X}^{0}$ and $\mathrm{X_{ind}}$ as a function of twist angle,
normalized to the intensity of monolayer $\mathrm{X}^{0}$ for each sample.
Interestingly, for angles in the range from $None$ to $None$ and $None$ to
$None$, we observe a PL intensity downconversion from direct K-K excitons to
indirect $\Gamma$-K. For these small stacking angles, we calculate an increase
of the splitting between $\mathrm{X}^{0}$ and $\mathrm{X_{ind}}$ from
approximately $150\text{\,}\mathrm{meV}$ to $200\text{\,}\mathrm{meV}$,
favouring a nonradiative decay of $\mathrm{X}^{0}$ to the lowest
$\mathrm{X_{ind}}$ excitons and thus inducing the shift of the PL emission
intensity.

In conclusion, we presented a systematic investigation of optical properties
of hBN-encapsulated $\mathrm{MoSe}_{2}$ homobilayers as a function of the
interlayer twist angle. To explain our experimental observations, we performed
moiré superlattice exciton calculations using an ab-initio-based continuum
model. Entirely new information was obtained on the optical and electronic
properties of different excitonic species in MoSe2 homobilayers, providing
insight into layer composition and relevant stacking configuration of
different moiré excitonic species. Controlled engineering of these properties
could lead towards the realization of on-chip quantum simulators based on
periodic excitonic arrays. Novel strongly-correlated excitonic states could
become readily available for exploration, such as exciton liquids and two-
dimensional condensates.

###### Acknowledgements.

We gratefully acknowledge financial support from the German Science Foundation
(DFG) via projects FI 947/8-1, DI 2013/5-1 and SPP-2244. Moreover, research
was partly funded via the Excellence Initiative in the framework of the
projects e-conversion (EXS 2089) and the Munich Center for Quantum Science and
Technology (MCQST EXS 2111). J.J.F. thanks the DFG for financial support via
the projects INST95-1496-1 and INST95-1642. M.K. gratefully acknowledges
funding from the International Max Planck Research School for Quantum Science
and Technology (IMPRS-QST). M.M.P. acknowledges TUM International Graduate
School of Science and Engineering (IGSSE). V.V. and M.F. acknowledge support
by the Alexander von Humboldt foundation. The authors would also like to thank
Prof. M. Kira for valuable discussions.

## References

  * Geim and Grigorieva (2013) A. K. Geim and I. V. Grigorieva, Nature 499, 419 (2013).
  * Liu et al. (2016) Y. Liu, N. O. Weiss, X. Duan, H.-C. Cheng, Y. Huang, and X. Duan, Nature Reviews Materials 1, 1 (2016).
  * Tang et al. (2020) Y. Tang, L. Li, T. Li, Y. Xu, S. Liu, K. Barmak, K. Watanabe, T. Taniguchi, A. H. MacDonald, J. Shan, et al., Nature 579, 353 (2020).
  * Regan et al. (2020) E. C. Regan, D. Wang, C. Jin, M. I. B. Utama, B. Gao, X. Wei, S. Zhao, W. Zhao, Z. Zhang, K. Yumigeta, et al., Nature 579, 359 (2020).
  * Xu et al. (2020) Y. Xu, S. Liu, D. A. Rhodes, K. Watanabe, T. Taniguchi, J. Hone, V. Elser, K. F. Mak, and J. Shan, Nature 587, 214 (2020).
  * Huang et al. (2021) X. Huang, T. Wang, S. Miao, C. Wang, Z. Li, Z. Lian, T. Taniguchi, K. Watanabe, S. Okamoto, D. Xiao, et al., Nature Physics 17, 715 (2021).
  * Jin et al. (2021) C. Jin, Z. Tao, T. Li, Y. Xu, Y. Tang, J. Zhu, S. Liu, K. Watanabe, T. Taniguchi, J. C. Hone, et al., Nature Materials 20, 940 (2021).
  * Li et al. (2021) H. Li, S. Li, E. C. Regan, D. Wang, W. Zhao, S. Kahn, K. Yumigeta, M. Blei, T. Taniguchi, K. Watanabe, et al., Nature 597, 650 (2021).
  * Jin et al. (2019) C. Jin, E. C. Regan, A. Yan, M. I. B. Utama, D. Wang, S. Zhao, Y. Qin, S. Yang, Z. Zheng, S. Shi, et al., Nature 567, 76 (2019).
  * Seyler et al. (2019) K. L. Seyler, P. Rivera, H. Yu, N. P. Wilson, E. L. Ray, D. G. Mandrus, J. Yan, W. Yao, and X. Xu, Nature 567, 66 (2019).
  * Tran et al. (2019) K. Tran, G. Moody, F. Wu, X. Lu, J. Choi, K. Kim, A. Rai, D. A. Sanchez, J. Quan, A. Singh, et al., Nature 567, 71 (2019).
  * Alexeev et al. (2019) E. M. Alexeev, D. A. Ruiz-Tijerina, M. Danovich, M. J. Hamer, D. J. Terry, P. K. Nayak, S. Ahn, S. Pak, J. Lee, J. I. Sohn, et al., Nature 567, 81 (2019).
  * Baek et al. (2020) H. Baek, M. Brotons-Gisbert, Z. Koong, A. Campbell, M. Rambach, K. Watanabe, T. Taniguchi, and B. D. Gerardot, Science advances 6, eaba8526 (2020).
  * Wang et al. (2020) L. Wang, E.-M. Shih, A. Ghiotto, L. Xian, D. A. Rhodes, C. Tan, M. Claassen, D. M. Kennes, Y. Bai, B. Kim, et al., Nature materials 19, 861 (2020).
  * Plochocka (2020) P. Plochocka, Nature Nanotechnology 15, 727 (2020).
  * Marcellina et al. (2021) E. Marcellina, X. Liu, Z. Hu, A. Fieramosca, Y. Huang, W. Du, S. Liu, J. Zhao, K. Watanabe, T. Taniguchi, et al., Nano Letters 21, 4461 (2021).
  * Sung et al. (2020) J. Sung, Y. Zhou, G. Scuri, V. Zólyomi, T. I. Andersen, H. Yoo, D. S. Wild, A. Y. Joe, R. J. Gelly, H. Heo, et al., Nature Nanotechnology 15, 750 (2020).
  * Yan et al. (2019) W. Yan, L. Meng, Z. Meng, Y. Weng, L. Kang, and X.-a. Li, The Journal of Physical Chemistry C 123, 30684 (2019).
  * Van Der Zande et al. (2014) A. M. Van Der Zande, J. Kunstmann, A. Chernikov, D. A. Chenet, Y. You, X. Zhang, P. Y. Huang, T. C. Berkelbach, L. Wang, F. Zhang, et al., Nano letters 14, 3869 (2014).
  * Liu et al. (2014) K. Liu, L. Zhang, T. Cao, C. Jin, D. Qiu, Q. Zhou, A. Zettl, P. Yang, S. G. Louie, and F. Wang, Nature communications 5, 1 (2014).
  * Kim et al. (2016) K. Kim, M. Yankowitz, B. Fallahazad, S. Kang, H. C. P. Movva, S. Huang, S. Larentis, C. M. Corbet, T. Taniguchi, K. Watanabe, et al., Nano Letters 16, 1989 (2016).
  * Wang et al. (2013) L. Wang, I. Meric, P. Y. Huang, Q. Gao, Y. Gao, H. Tran, T. Taniguchi, K. Watanabe, L. M. Campos, D. A. Muller, et al., Science 342, 614 (2013).
  * Pizzocchero et al. (2016) F. Pizzocchero, L. Gammelgaard, B. S. Jessen, J. M. Caridad, L. Wang, J. Hone, P. Bøggild, and T. J. Booth, Nature communications 7, 1 (2016).
  * Frisenda et al. (2017) R. Frisenda, Y. Niu, P. Gant, A. J. Molina-Mendoza, R. Schmidt, R. Bratschitsch, J. Liu, L. Fu, D. Dumcenco, A. Kis, et al., Journal of Physics D: Applied Physics 50, 074002 (2017).
  * Evans and Hazelwood (1971) B. Evans and R. Hazelwood, physica status solidi (a) 4, 181 (1971).
  * Anedda and Fortin (1980) A. Anedda and E. Fortin, Journal of Physics and Chemistry of Solids 41, 865 (1980).
  * Li and Galli (2007) T. Li and G. Galli, The Journal of Physical Chemistry C 111, 16192 (2007).
  * Splendiani et al. (2010) A. Splendiani, L. Sun, Y. Zhang, T. Li, J. Kim, C.-Y. Chim, G. Galli, and F. Wang, Nano letters 10, 1271 (2010).
  * Ugeda et al. (2014) M. M. Ugeda, A. J. Bradley, S.-F. Shi, H. Felipe, Y. Zhang, D. Y. Qiu, W. Ruan, S.-K. Mo, Z. Hussain, Z.-X. Shen, et al., Nature materials 13, 1091 (2014).
  * Raja et al. (2017) A. Raja, A. Chaves, J. Yu, G. Arefe, H. M. Hill, A. F. Rigosi, T. C. Berkelbach, P. Nagler, C. Schüller, T. Korn, et al., Nature communications 8, 1 (2017).
  * Florian et al. (2018) M. Florian, M. Hartmann, A. Steinhoff, J. Klein, A. W. Holleitner, J. J. Finley, T. O. Wehling, M. Kaniber, and C. Gies, Nano letters 18, 2725 (2018).
  * Cho and Berkelbach (2018) Y. Cho and T. C. Berkelbach, Physical Review B 97, 041409 (2018).
  * Segura et al. (2018) A. Segura, L. Artús, R. Cuscó, T. Taniguchi, G. Cassabois, and B. Gil, Physical Review Materials 2, 024001 (2018).
  * Lin et al. (2014) Y. Lin, X. Ling, L. Yu, S. Huang, A. L. Hsu, Y.-H. Lee, J. Kong, M. S. Dresselhaus, and T. Palacios, Nano letters 14, 5569 (2014).
  * Stier et al. (2016) A. V. Stier, N. P. Wilson, G. Clark, X. Xu, and S. A. Crooker, Nano letters 16, 7054 (2016).
  * Linderälv et al. (2022) C. Linderälv, J. Hagel, S. Brem, E. Malic, and P. Erhart, arXiv preprint arXiv:2205.15616 (2022).
  * Wang et al. (2017) Y. Wang, Z. Wang, W. Yao, G.-B. Liu, and H. Yu, Physical Review B 95, 115429 (2017).
  * Hagel et al. (2021) J. Hagel, S. Brem, C. Linderälv, P. Erhart, and E. Malic, Physical Review Research 3, 043217 (2021).
  * Wu et al. (2018) F. Wu, T. Lovorn, and A. MacDonald, Physical Review B 97, 035306 (2018).
  * Moon and Koshino (2013) P. Moon and M. Koshino, Physical Review B 87, 205404 (2013).
  * Wang et al. (2019) J. Wang, J. Ardelean, Y. Bai, A. Steinhoff, M. Florian, F. Jahnke, X. Xu, M. Kira, J. Hone, and X.-Y. Zhu, Science advances 5, eaax0145 (2019).
  * He et al. (2014) J. He, K. Hummer, and C. Franchini, Physical Review B 89, 075409 (2014).
  * Ruiz-Tijerina and Fal’ko (2019) D. A. Ruiz-Tijerina and V. I. Fal’ko, Physical Review B 99, 125424 (2019).
  * Weston et al. (2020) A. Weston, Y. Zou, V. Enaldiev, A. Summerfield, N. Clark, V. Zólyomi, A. Graham, C. Yelgel, S. Magorrian, M. Zhou, et al., Nature Nanotechnology 15, 592 (2020).
  * Brem et al. (2020) S. Brem, C. Linderälv, P. Erhart, and E. Malic, Nano letters 20, 8534 (2020).