Discrete Interactions between a few Interlayer Excitons Trapped at a MoSe$_2$-WSe$_2$ Heterointerface

TL;DR

This study investigates the interactions of a few interlayer excitons confined at a heterointerface in TMD bilayers, revealing discrete multiparticle states and quantifying their dipolar and exchange interactions, advancing understanding of quantum correlations in 2D systems.

Contribution

It provides the first detailed experimental measurement of discrete multiparticle interlayer exciton states and their interactions in a nanoscale confined TMD heterobilayer.

Findings

Observation of discrete multiparticle exciton states with N_{IX}~1-5.

Measurement of interlayer biexciton energy splitting and exchange interactions.

Numerical modeling supports the experimental confinement lengthscale of ~3 nm.

Abstract

Interlayer excitons (IXs) in hetero-bilayers of transition metal dichalcogenides (TMDs) represent an exciting emergent class of long-lived dipolar composite bosons in an atomically thin, near-ideal two-dimensional (2D) system. The long-range interactions that arise from the spatial separation of electrons and holes can give rise to novel quantum, as well as classical multi-particle correlation effects. In order to acquire a detailed understanding of the possible many-body effects, the fundamental interactions between individual IXs have to be studied. Here, we trap a tunable number of dipolar within a nanoscale confinement potential induced by placing a MoSe$_2$-WSe$_2$ hetero-bilayer (HBL) onto an array of SiO$_2$ nanopillars. We control the mean occupation of the IX trap via the optical excitation level and observe discrete sharp-line emission from different configurations of interacting IXs. We identify these features as different multiparticle states with $N_{IX}\sim1-5$ via their power dependencies and directly measure the hierarchy of dipolar and exchange interactions as $N_{IX}$ increases. The interlayer biexciton ($N_{IX}=2$) is found to be an emission doublet that is blue-shifted from the single exciton by $\Delta E=(8.4\pm0.6)$ meV and split by $2J=(1.2\pm0.5)$ meV. The blueshift is even more pronounced for triexcitons ($(12.4\pm0.4)$ meV), quadexcitons ($(15.5\pm0.6)$ meV) and quintexcitons ($(18.2\pm0.8)$ meV). These values are shown to be mutually consistent with numerical modelling of dipolar excitons confined to a harmonic trapping potential having a confinement lengthscale in the range $\ell\approx 3$ nm. Our results contribute to the understanding of interactions between IXs in TMD HBLs at the discrete limit of only a few excitations and represent a key step towards exploring quantum correlations between them.

Figures

f and 1e, respectively. Supplementary Figure [5](#Sx1.F5)a displays example spectra at pillar A for different powers spanning more than three orders of magnitude. As for the dataset presented in Fig. 2, we observe a series of emission features in the window $1.28\text{\,}\mathrm{eV}1.35\text{\,}\mathrm{eV}$. As the excitation power increases from the lowest levels, the spectrum clearly evolves from a single sharp emission line to a multiplet of excitonic features. Since the overlap of the peaks in these measurements are enhanced due to the broader linewidths compared to the data presented in F

f and 1e, respectively. Supplementary Figure [5](#Sx1.F5)a displays example spectra at pillar A for different powers spanning more than three orders of magnitude. As for the dataset presented in Fig. 2, we observe a series of emission features in the window $1.28\text{\,}\mathrm{eV}1.35\text{\,}\mathrm{eV}$. As the excitation power increases from the lowest levels, the spectrum clearly evolves from a single sharp emission line to a multiplet of excitonic features. Since the overlap of the peaks in these measurements are enhanced due to the broader linewidths compared to the data presented in F

b). We thus plot, as a function of excitation power, the ratios of the intensities of the multi-IX states with one state whose particle number we assume to be one below the examined state. This is shown in Suppl. Fig. [9](#Sx1.F9)b. In all cases, we obtain an approximately linear excitation-power dependence which is not influenced by the onset of saturation effects (marked by the gray dashed line), confirming our initial exciton number assignments.

Tables

Table 1. Supplementary Table 1: Parameters extracted from biexponential fits to the time-resolved PL data presented in fig.  4 . The fitted model I ​ ( t ) = A s ​ l ​ o ​ w ​ e − t / τ s ​ l ​ o ​ w + A f ​ a ​ s ​ t ​ e − t / τ f ​ a ​ s ​ t 𝐼 𝑡 subscript 𝐴 𝑠 𝑙 𝑜 𝑤 superscript 𝑒 𝑡 subscript 𝜏 𝑠 𝑙 𝑜 𝑤 subscript 𝐴 𝑓 𝑎 𝑠 𝑡 superscript 𝑒 𝑡 subscript 𝜏 𝑓 𝑎 𝑠 𝑡 I(t)=A_{slow}e^{-t/\tau_{slow}}+A_{fast}e^{-t/\tau_{fast}} directly includes the slow/fast timescale τ s ​ l ​ o ​ w / f ​ a ​ s ​ t subscript 𝜏 𝑠 𝑙 𝑜 𝑤 𝑓 𝑎 𝑠 𝑡 \tau_{slow/fast} . The slow/fast proportion is calculated as A s ​ l ​ o ​ w / f ​ a ​ s ​ t ⋅ τ s ​ l ​ o ​ w / f ​ a ​ s ​ t A s ​ l ​ o ​ w ⋅ τ s ​ l ​ o ​ w + A f ​ a ​ s ​ t ⋅ τ f ​ a ​ s ​ t ⋅ subscript 𝐴 𝑠 𝑙 𝑜 𝑤 𝑓 𝑎 𝑠 𝑡 subscript 𝜏 𝑠 𝑙 𝑜 𝑤 𝑓 𝑎 𝑠 𝑡 ⋅ subscript 𝐴 𝑠 𝑙 𝑜 𝑤 subscript 𝜏 𝑠 𝑙 𝑜 𝑤 ⋅ subscript 𝐴 𝑓 𝑎 𝑠 𝑡 subscript 𝜏 𝑓 𝑎 𝑠 𝑡 \frac{A_{slow/fast}\cdot\tau_{slow/fast}}{A_{slow}\cdot\tau_{slow}+A_{fast}\cdot\tau_{fast}} .

Position	Spectral range	Slow timescale	Slow proportion	Fast timescale	Fast proportion
Pillar	Low energy	$118.8\pm 2.1\mathrm{ns}$	$\left(95\pm 3\right)\%$	$8.1\pm 0.5\mathrm{ns}$	$\left(5\pm 3\right)\%$
Flat HBL	Low energy	$111\pm 4\mathrm{ns}$	$\left(76\pm 5\right)\%$	$6.16\pm 0.12\mathrm{ns}$	$\left(24\pm 5\right)\%$
Pillar	High energy	$74.5\pm 1.5\mathrm{ns}$	$\left(80\pm 3\right)\%$	$9.0\pm 0.3\mathrm{ns}$	$\left(20\pm 3\right)\%$
Flat HBL	High energy	$76.5\pm 1.4\mathrm{ns}$	$\left(78\pm 3\right)\%$	$9.48\pm 0.16\mathrm{ns}$	$\left(22\pm 3\right)\%$

Equations

$$H=H_0+H_{\mathrm{int}},$$

$$H_0=\sum _{i=1,2}\left[-\frac{{\hslash }^2}{2m^{\ast }}\frac{1}{{\rho }_i}\frac{\partial }{\partial {\rho }_i}\left({\rho }_i\frac{\partial }{\partial {\rho }_i}\right)-\frac{1}{4\pi \epsilon }\left(\frac{1}{\sqrt{{\rho }_i^2+d^2}}+\frac{1}{\sqrt{({\rho }_i+(-{)}^i\Delta \rho {)}^2+d^2}}\right)+\frac{1}{2}{m}^{\ast }{\Omega }_0^2{\rho }_i^2\right]$$

$$H_{\mathrm{int}}$$

$$J(\Delta\rho)=\frac{1}{3}\int_{-\Delta\rho/\sqrt{2}}^{\Delta\rho/\sqrt{2}}\mathrm{d}y\big{|}\Psi(x,y)\frac{\partial\Psi(x,y)}{\partial x}\big{|}_{x=0},$$

$$E_{LIX2B/A}=2E_{LIX1}+\frac{2d^2}{4\pi \epsilon (\Delta \rho {)}^3}+\frac{{\hslash }^2{\left(\Delta \rho \right)}^2}{4m^{\ast }{l}^4}\pm J\left(\Delta \rho \right)=2E_{LIX1}+E_{\mathrm{binding}}\pm J\left(\Delta \rho \right)$$

$$\frac{d{n}_N}{dt}=\frac{(N+1)n_{N+1}}{{\tau }_{N+1}}-\frac{N{n}_N}{{\tau }_N}+G_{N-1}{n}_{N-1}-G_N{n}_N$$

$$n_N=\frac{{\left(G\tau \right)}^N}{N!}{e}^{-G\tau }$$

$$\frac{n_{N+1}}{n_N}=\frac{G_N{\tau }_{n+1}}{N+1}\propto P_{exc}$$

$$U^{2X}$$

$$U^{3X}$$

$$U^{4X}$$

$$U^{5X}$$

$$U_0^{2X}$$

$$U_0^{3X}$$

$$U_0^{4X}$$

$$U_0^{5X}$$

$$E^{2X\to 1X}$$

$$E^{3X\to 2X}$$

$$E^{4X\to 3X}$$

$$E^{5X\to 4X}$$

$$\frac{E^{3X\to 2X}-E^X}{E^{2X\to 1X}-E^X}$$

$$\frac{E^{4X\to 3X}-E^X}{E^{2X\to 1X}-E^X}$$

$$\frac{E^{5X\to 4X}-E^X}{E^{2X\to 1X}-E^X}$$

$${\left(\frac{E^{3X\to 2X}-E^X}{E^{2X\to 1X}-E^X}\right)}_{exp}$$

$${\left(\frac{E^{4X\to 3X}-E^X}{E^{2X\to 1X}-E^X}\right)}_{exp}$$

$${\left(\frac{E^{5X\to 4X}-E^X}{E^{2X\to 1X}-E^X}\right)}_{exp}$$

Full text

1 Introduction

HBLs of van-der-Waals bonded 2D TMDs feature an atomically sharp interface with a type-II band alignment and frontier orbital couplings defined by the mutual angular orientation of the basal plane of the component layers.1 Such atomically thin 2D heterointerfaces can host IXs – Coulomb-bound states of electrons ($e$) and holes ($h$) located in the different component monolayers.1 The spatial separation of $e$ and $h$ gives rise to a significantly increased lifetime compared to intralayer excitons 1, 3, 17 and a static electric out-of-plane dipole moment with a magnitude $|p_{z}|/e\sim$0.7\text{,}\mathrm{nm}$$ for the MoSe2-WSe2 heterointerface1. The enhanced radiative lifetime of IXs allows them to cool toward the lattice temperature before recombination takes place, whereas the large static dipole moment facilitates tuning of the IX energy1, 18, 19 and lifetime 18, 19 using static electric fields. In optical experiments, dipole-dipole repulsion between IXs gives rise to strongly blueshifting emission with increasing IX density and the potential to explore the many-body physics of dipolar composite bosons in a solid-state system1, 3, 17. It has also been demonstrated that applying strain to the HBL can strongly tune the IX energy20; a method that is routinely used to define exciton trapping potentials in monolayer TMDs. 12, 13, 21 Overall, IXs constitute a highly tunable platform for optically exploring and controlling interacting gases of dipolar excitons in the solid state.

2 Results and Discussion

Trapping interlayer excitons

Our sample structure is schematically depicted in Fig. 1a. It consists of stacked monolayers of WSe2 and MoSe2 that form a HBL which is then transferred onto a SiO2 substrate patterned with nanopillars having a height of $\sim$90\text{,}\mathrm{nm} and diameter $\sim$130\text{\,}\mathrm{nm} (see Methods). Figure 1b shows an optical micrograph of the sample in which the monolayer MoSe2, WSe2 and HBL regions can be seen, as well as the nanopillars in the HBL region. The contours of the WSe2 and MoSe2 monolayers are outlined with blue and orange dashed lines, respectively, and pillars that did not pierce the TMD layers appear as black dots. Figure 1c shows normalized spectra recorded from three different unstrained sample regions, namely the WSe2 monolayer, the MoSe2 monolayer and the WSe2-MoSe2-HBL away from the nanopillars. These spectra were recorded at $T=$10\text{,}\mathrm{K} with cw excitation at $633\text{\,}\mathrm{nm}$ and a power of $P_{ex}$=$35\text{\,}\mathrm{\SIUnitSymbolMicro W}$ focused to $\sim$1\text{\,}\mathrm{\SIUnitSymbolMicro m}. Both, WSe2 and MoSe2, have prominent emission features at their characteristic exciton and trion energies at $1.75\text{\,}\mathrm{eV}$ and $1.72\text{\,}\mathrm{eV}$, 22 and $1.66\text{\,}\mathrm{eV}$ and $1.63\text{\,}\mathrm{eV}$, 23 respectively. Furthermore, the WSe2 component monolayer exhibits a broader band of unresolved emission at lower energies ($1.60\text{\,}\mathrm{eV}1.70\text{\,}\mathrm{eV}$) that has been identified as stemming from dark excitons 24, 25 and trions, 26 charged biexcitons 27, 28, 29, 30, 31 and defect-bound excitons.32, 33 Each of these emission features are also present in the HBL region, albeit strongly quenched, due to interlayer charge transfer that occurs over sub-picosecond timescales, faster than the exciton lifetime.34, 35, 36 The most prominent emission feature from the HBL is at \sim$$1.38\text{\,}\mathrm{eV} and is attributed to IXs formed by an electron in the MoSe2 and a hole in the WSe2. 1 This attribution is supported by the emission energy and the characteristically asymmetric lineshape featuring red-shifted emission from momentum-indirect IXs17. Figure 1e shows the power dependence of the IX emission revealing a significant blueshift ($\geq$20\text{,}\mathrm{meV}) arising from repulsive dipolar interactions in the gas of IX for increasing excitation levels[9](#bib.bib9), [3](#bib.bib3). We note that this blueshift does *not* occur over the whole range of excitation levels, but rather exhibits an onset close to $P_{ex}\sim$30\text{\,}\mathrm{nW}, marking a low-density regime where dipolar interactions are weak (see Suppl. Fig. 2). Furthermore, the low-density IX regime also manifests itself by a deviation from the $\sqrt{P_{ex}}$ dependence of the IX emission intensity, which might be indicative of non-radiative exciton-exciton annihilation37 only above a threshold density.

The nanopillars have a strong impact on the IX emission spectra. Typical data is presented in Fig. 1d that compares emission of unstrained regions of the HBL with that recorded from three pillar sites (labelled A, B and C as in Fig. 1a). The characteristic IX emission at $1.38\text{\,}\mathrm{eV}$ can still be observed in each emission spectrum since the probed sample volume is much larger than the nanopillar size. However, additional peaks emerge in the vicinity of the nanopillar sites, red-shifted by up to \sim$$100\text{\,}\mathrm{meV} from the peak of the free IX emission. The observation of a pronounced red-shift is similar to reports of strain-induced bandgap modulation in monolayer MoSe2 21, MoS2 38 and particularly WSe2 12, 13, where this effect is routinely used to generate quantum emitters. We attribute these features to localized interlayer excitons (LIXs) at the pillar sites, most likely due to strain-related trapping potentials around the nanopillar. The interlayer exciton character is confirmed by PL-excitation spectroscopy detected on the LIX features, that clearly reveals free exciton fingerprints of both MoSe2 and WSe2 component monolayers, and time-resolved PL measurements, that show similar radiative decay times for localized and free IXs (see Supplementary Notes 3 and 4). To illustrate the confinement of the LIX features to the pillar sites, the inset of Fig. 1d shows the spatially resolved emission intensity of an example LIX peak at pillar site A. The spatial extent of this feature is clearly limited by the spatial resolution of our $\mu$PL setup ($1\text{\,}\mathrm{\SIUnitSymbolMicro m}$), indicative of a point-like source. Figures 1e and 1f compare normalized IX $\mu$PL spectra plotted as a function of excitation power from the flat region of the HBL and on pillar A, respectively. In both cases, free IX emission is observed that blueshifts by up to $\sim$20\text{,}\mathrm{meV}$$ over the examined power range. The spectra recorded from pillar A additionally exhibit discrete LIX emission lines. Individually, these features do not blueshift with excitation power, in strong contrast to the broad IX peak. Instead, changes of excitation power result in a redistribution of the PL intensity among the different LIX peaks. However, if one analyzes the entirety of all LIX emission peaks and determines its center of the spectral weight as a function of excitation power (see Suppl. Fig. 5), one finds a collective LIX blueshift which is slightly stronger than that of the free IX, consistent with stronger interactions for IXs spatially confined to a trapping potential. This effect of the LIX emission shifting among several discrete peaks with varying excitation power was observed for all pillars exhibiting LIX emission and was reproducible following thermal cycling.

Spectroscopy of trapped few-exciton states

The LIX emission spectra of pillar A for excitation powers ranging from $79\text{\,}\mathrm{nW}$ to $\geq$500\text{,}\mathrm{nW}, two orders of magnitude weaker than the power used to record the results depicted in Fig. [1](#S2.F1), are presented in Fig. [2](#S2.F2)a. At the lowest excitation power, the emission spectrum is found to be dominated by a single sharp line at $1.302\text{\,}\mathrm{eV}$, labelled LIX1 in the figure. The spectral linewidth of this peak is $550\text{\,}\mathrm{\SIUnitSymbolMicro eV}$ and it is broadened on the low-energy flank, most likely due to coupling to acoustic phonons[39](#bib.bib39). The integrated intensity of LIX1 is plotted in Fig. [2](#S2.F2)b as a function of the excitation power. From the lowest excitation levels investigated it increases linearly with power, as can be seen from the solid line on the figure that shows a fit to the data with $I=I_{0}\cdot P^{\alpha}$ yielding $\alpha=0.9\pm 0.1$. The dotted line on the figure shows our expectations for a linear power dependence. This observation, combined with the saturation and eventual decrease of the peak intensity at the highest excitation levels investigated, clearly identifies this peak as stemming from the recombination of a single IX in the local trapping potential, as is commonly observed in the emission characteristics of III-V quantum dots and other discrete quantum emitters[40](#bib.bib40). As the excitation power is increased from the minimum of $79\text{\,}\mathrm{nW}$, additional emission features emerge at higher energies, the most prominent of which is a pair of lines labelled LIX${}_{2\text{A}}$ and LIX${}_{2\text{B}}$ in Fig. [2](#S2.F2). In contrast to LIX1, LIX${}_{2\text{A/B}}$ *both* exhibit a near quadratic power dependence ($\alpha_{2\text{A}}=1.96\pm 0.15$, $\alpha_{2\text{B}}=1.57\pm 0.12$) and each saturates and quenches at the highest pump powers investigated (see Fig. [2](#S2.F2)b). For excitation power in excess of $\sim$300\text{\,}\mathrm{nW}, additional manifolds of emission lines emerge, labelled LIX${}_{3\text{A/B}}$, LIX4 and LIX5 on Fig. 2a. The power dependence of each emergent set of lines clearly becomes successively more superlinear in behaviour - LIX${}_{3\text{A/B}}$ have a power dependence that is very close to cubic ($\alpha_{3\text{A}}=3.48\pm 0.27$, $\alpha_{3\text{B}}=3.6\pm 0.7$), LIX4 close to quartic (see dotted lines on figure). These observations mirror previous findings for multi-exciton states in semiconductor quantum dots40 and double quantum wells.8 Since there were only a few data points before saturation for some of the peaks, we substantiate our particle number assignment by comparing ratios of intensities of two peaks to an independent capture model (see Suppl. Fig. 6). This method uses data from the entire available power range.

Figure 2c shows the emission energies of the various LIX peaks as extracted from multiple peak fitting. The energy of the individual LIXi peaks is not significantly power dependent, beyond spectral wandering. This can be prominently seen in the LIX1 peak, which jitters by up to $\sim 1\text{\,}\mathrm{meV}$ between measurements. Similar spectral wandering has also been reported from discrete emitters in monolayer WSe2, where it is commonly attributed to time-dependent fluctuations in the local electrostatic environment. 41, 42 The mid-point of the biexciton LIX2 doublet is significantly blueshifted by $8.4\pm 0.6\text{\,}\mathrm{meV}$ with respect to LIX1 and the two components are separated by $\sim$1.2\pm 0.5\text{,}\mathrm{meV}$$.

Modelling interactions in multi-exciton complexes

In order to quantitatively understand the strong blueshift of the localized interlayer biexciton from the single exciton and the splitting of the biexciton into a doublet, we modelled the interaction of two IX confined to a harmonic trapping potential, as characterised by a confinement lengthscale $\ell=\sqrt{\hbar/m^{*}\Omega_{0}}$. Here, $m^{*}$ denotes the total effective mass of the exciton and $\Omega_{0}$ is the angular frequency of an exciton in the harmonic potential. Figure 3a illustrates how we calculate the equilibrium separation $\Delta\rho$ between two IX by minimizing the sum of repulsive (dipolar) and attractive (harmonic) potential energies. For values of $\ell\sim$3\text{,}\mathrm{nm}4\text{,}\mathrm{nm}, reference to Fig. [3](#S2.F3)a shows that $\Delta\rho\sim$5\text{\,}\mathrm{nm}6.5\text{\,}\mathrm{nm}. We now continue to discuss the calculation of the direct and exchange Coulomb interaction energies for the $LIX_{2}$ and $LIX_{1}$ states. For this purpose, we modified the method outlined in Ref. 16 and start with a Hamiltonian describing the biexciton complex of the following form:

[TABLE]

where

[TABLE]

describes non-interacting excitons exposed to a weak harmonic confinement potential and

[TABLE]

represents the Coulomb interactions between the excitons. In equations 2 and 2, the coordinates $\rho_{1}$ and $\rho_{2}$ are in-plane projections of the relative electron-hole coordinates of the interlayer excitons, as depicted in Fig. 3b, and $\Delta\rho$ was introduced before. Moreover, d\approx 0.7\leavevmode\nobreak\nm is the interlayer separation of the HBL that determines the strength of the dipole-dipole repulsion and $\epsilon=\epsilon_{0}\epsilon_{r}$ describes the effective dielectric screening by the surrounding vacuum (upper) and SiO2 (lower) media. We note that, in the simplest picture in which we average the values for vacuum (upper cladding) and SiO2 (lower substrate), the effective dielectric constant would be $\epsilon_{r}=(\epsilon_{vac}+\epsilon_{SiO_{2}})/2\sim 2.5$. For the effective masses $m^{*}$, we use $m^{*}=0.8m_{0}$ ($0.45m_{0}$) for the MoSe2 electrons43 (WSe2 holes44).

The indistinguishability of the two LIX composite bosons forming the biexciton implies that the total wavefunction must be symmetric with respect to IX exchange. This condition can only be satisfied when both spin and spatial parts of the two-exciton wavefunction are simultaneously symmetric or antisymmetric, respectively, as illustrated schematically in Fig. 3c. Thus, the spatial part of the two-IX (biexciton) wavefunction can take the form $\Psi\sim\frac{1}{\sqrt{2}}[\psi_{IX}(\rho_{1})\psi_{IX}(\rho_{2})\pm\psi_{IX}(\rho_{2})\psi_{IX}(\rho_{1})]$, giving rise to two, energetically distinct biexciton states separated by the exchange energy $2J(\Delta\rho)$. Here, $J(\Delta\rho)$ is given by 45

[TABLE]

where $\Psi(x,y)$ is the ground state of eqn. (1) that has been transformed to the coordinate system $x=(\rho_{1}-\rho_{2}-\Delta\rho)/\sqrt{2}$ and $y=(\rho_{1}+\rho_{2})/\sqrt{2}$. From this expression, we numerically calculate the exchange splitting to obtain the absolute energies of the two biexciton states:

[TABLE]

Here, a positive binding energy refers to an increase in energy of the biexciton state compared to two isolated excitons. Representative results of our calculations are presented in Fig. 3d that shows the numerically evaluated biexciton binding energy and exchange splitting as a function of the effective confinement lengthscale $\ell$. The three curves presented on the figure correspond to effective dielectric constants $\epsilon=2-4$. The experimentally observed blueshift of LIX2 of $8.4\pm 0.6\text{\,}\mathrm{meV}$ is consistent with a confinement lengthscale in the range $\ell=$2.8\text{,}\mathrm{nm}3.3\text{,}\mathrm{nm}, with relatively weak sensitivity to the effective dielectric constant. We interpret the existence of two biexciton states, split by $1.2\pm 0.5\text{\,}\mathrm{meV}$, as reflecting the existence of two distinct spatial wavefunctions with even and odd symmetry with respect to IX exchange, as discussed above. The red curves presented in Fig. [3](#S2.F3)d show the calculated exchange splittings $J(\ell)$ for the same range of effective dielectric constants. The red-shaded region indicates the measured splitting of the biexciton peak. Remarkably, our model calculations simultaneously reproduce *both* the observed biexciton binding energy and the exchange splitting for a dielectric constant in the range $\epsilon\gtrapprox 3$, very close to the average dielectric constant of the environment $\epsilon_{vac-SiO_{2}}=2.5$, and an effective confinement lengthscale of $\ell\approx$3\text{\,}\mathrm{nm}. Thus, we conclude that the strain-induced trapping potentials in the vicinity of the nanopillars give rise to nanometer-scale confinement traps in which the LIX interact via direct and exchange Coulomb interactions.

Figure 2c also shows that the triexciton (LIX3), quadexciton (LIX4) and quintexciton (LIX5) transition lines are even more strongly blueshifted from LIX1, as compared to LIX2. We measure blueshifts of $12.4\pm 0.4\text{\,}\mathrm{meV}$ for LIX3, $15.5\pm 0.6\text{\,}\mathrm{meV}$ for LIX4 and $18.2\pm 0.8\text{\,}\mathrm{meV}$ for LIX5. We extended our calculations for binding energies based on dipole-dipole repulsion in the harmonic trap to three to five excitons (see Supplementary Note 10). The total energy is minimized by maximizing the equilibrium spacing between the particles. Thus, the IXs in LIX3 (LIX4/LIX5) would be expected to be located at the corners of an equilateral triangle (square/regular pentagon) centered around the trap minimum. We find that independent of the system parameters (e.g. dipole moment $d$ or confinement length $\ell$) the blueshift of the LIX3 (LIX4/LIX5) is expected to be about $1.36\times\Delta E_{\mathrm{LIX2}}$ ($1.92\times\Delta E_{\mathrm{LIX2}}$/$2.53\times\Delta E_{\mathrm{LIX2}}$), where $\Delta E_{\mathrm{LIX2}}$ is the localized biexciton blueshift. Remarkably, the experimentally observed shifts amount to $\left(1.48\pm 0.12\right)\times\Delta E_{\mathrm{LIX2}}$ ($\left(1.85\pm 0.15\right)\times\Delta E_{\mathrm{LIX2}}$, $\left(2.17\pm 0.18\right)\times\Delta E_{\mathrm{LIX2}}$), respectively, in excellent agreement with this simple prediction.

In summary, we have directly probed interactions between a discrete number $N_{IX}=$15 of interlayer excitons moving within a localized trapping potential. Power-dependent measurements allow us to identify different few-body states and, moreover, the observed interactions are found to be all mutually in agreement with a localization lengthscale $\ell\approx$3\text{\,}\mathrm{nm}. We reproduce the observed blueshift of the biexciton by $8.4\pm 0.6\text{\,}\mathrm{meV}$ as well as the observed biexciton fine structure with a splitting of $1.2\pm 0.5\text{\,}\mathrm{meV}$ in a model including the dipolar and exchange interaction of two IX in a harmonic trapping potential. The relative energies of states with one to five trapped IXs are well captured by our model, which predicts maximal spacing of the IXs within the limitations of the trapping potential. The significant interaction energies observed in this paper suggest that strongly correlated multiexciton phenomena, such as Wigner crystallization, may be realized in layered van der Waals heterostructures, opening up new avenues for nonlinear coherent optical control and spin-optronics with interlayer excitons.

After preparation and submission of the manuscript, we became aware of similar results by another group. 46

3 Methods

Sample fabrication: We fabricated the sample on a Si-[100] wafer with a $285\text{\,}\mathrm{nm}$ thick SiO2 layer on top. We defined the pillar grid via electron beam lithography on a negative resist (AR-N 7520.07) and etched it into the substrate using F-based reactive ion etching. The pillars have a diameter and height of \sim$$130\text{\,}\mathrm{nm} and \sim$$90\text{\,}\mathrm{nm}, respectively. We exfoliated crystals of WSe2 and MoSe2 (hqgraphene.com) onto polydimethylsiloxane (PDMS, Gel-Pak) stamps, identified monolayers in an optical microscope via reflectivity contrast and then subsequently transferred the WSe2 and MoSe2 using a mask aligner in a cleanroom environment. Each transfer was followed by an annealing step in vacuum ($200\text{\,}\mathrm{\SIUnitSymbolCelsius}$, $30\text{\,}\mathrm{min}$ after the WSe2 transfer; $150\text{\,}\mathrm{\SIUnitSymbolCelsius}$, $20\text{\,}\mathrm{min}$ after the MoSe2 transfer which completed the heterostructure). During assembly of the HBL, we took efforts to orientate the crystal axes of the two materials by aligning long crystalline edges of the monolayers.

Optical spectroscopy measurements: Photoluminescence measurements were carried out in two different setups designed for low-temperature confocal microscopy. The data from Fig. 1 was acquired at a setup using a liquid-He flow cryostat at $T=$10\text{,}\mathrm{K} and a HeNe laser ($\lambda=$633\text{\,}\mathrm{nm}) focused to a diameter of $\sim$1\text{,}\mathrm{\SIUnitSymbolMicro m}. In order to counteract thermal cycling effects (see Supplementary Note 11), the measurements presented in Fig. [2](#S2.F2) were carried out in a setup equipped with a closed-cycle cryostat at $T=$4\text{\,}\mathrm{K}. The excitation laser used in this setup was at $\lambda=$532\text{,}\mathrm{nm}$$ and was defocused due to chromatic aberration in the used objective.

We note that the data presented in Fig. 1f and 2 were both recorded from pillar A, but following thermal cycling of the sample to room temperature. This procedure was typically found to result in a shift of the absolute energy of different groups of sharp-line emission features by $\sim$\pm 10\text{,}\mathrm{meV}. Nevertheless, following each cool-down, the generic form of the group of peaks was retained. Supplementary Figure 5 shows power-dependent $\mu PL$ spectra recorded from pillar A (from the same measurement as the data presented in Fig. [1](#S2.F1)f), illustrating how the single sharp emission feature is observed close to $1.292\text{\,}\mathrm{eV}$ at the lowest excitation power investigated, with additional peaks emerging in the range $\sim$5\text{\,}\mathrm{meV}30\text{\,}\mathrm{meV} to higher energy as $P_{ex}$ increases.

Data analysis: The emission energies shown in Fig. 2c were extracted from Gaussian fits to the PL spectra. In some power regimes, the large peak overlap combined with the low signal-to-noise ratio hindered proper fitting of the data. For this reason, the PL intensities presented in Fig. 2b were determined by integrating the PL signal in an interval of $0.12\text{\,}\mathrm{meV}$ around the respective center energies. The only exception is the LIX1 peak for which fitting was possible over the whole power range due to its spectral isolation. The data points for LIX3A and LIX3B in Fig. 2b initially overlapped. Thus, the data for LIX3A was slightly offset for better visibility.

4 Data availability

The data that support the findings of this study are available from the corresponding author upon reasonable request.

5 Acknowledgements

We gratefully acknowledge funding from the International Max Planck Research School for Quantum Science and Technology (IMPRS-QST), the DFG Clusters of Excellence MCQST, e-Conversion and NIM, DFG Project FI 947/5-1 SQAM, the European Union’s Horizon 2020 research and innovation programme under grant agreement No. 820423 (S2QUIP), the BMBF Verbundprojekt QLink-X and the Alexander von Humboldt Foundation. K.M. acknowledges support from the Bavarian Academy of Sciences and Humanities. In addition, we thank Richard Schmidt and his team for very helpful discussions.

6 Author contributions

M.K., B.D.G., K.M. and J.J.F. conceived and designed the experiments. M.K. and J.G. prepared the sample. M.K., M.B.-G., J.G. and M.M. performed the experiments. M.K. analyzed the data. J.K. modelled the exchange interaction. M.K., K.M. and J.J.F. wrote the manuscript with input from all authors.

7 Competing interests

The authors declare that there are no competing interests.

Supplementary Information to: Discrete Interactions between a few Interlayer Excitons Trapped at a MoSe2-WSe2 Heterointerface

Supplementary Note 1: Enhanced micrograph of the sample

Supplementary Note 2: Excitation-power dependence of free IX emission

Supplementary Figure 2a shows exemplary PL spectra for different excitation powers recorded from a HBL region of the sample away from SiO2 nanopillars, together with comprehensive power-dependent data shown as a false-color plot in Suppl. Fig. 2b. The power-dependence of the integrated IX intensity is plotted in Suppl. Fig. 2c on a double-logarithmic representation. We observe a typical $\sqrt{P_{ex}}$ dependence with a power factor $\alpha=0.43\pm 0.05$ (red line) consistent with exciton-exciton annihilation at elevated excitation levels1. Careful examination of fig. 2c for $P_{ex}\leq P_{ex}^{0}\approx$20\text{,}\mathrm{nW} shows a higher power factor. A power factor of $\alpha=0.66$ would be consistent with Auger recombination. The increased power factor could also hint at a transition towards a regime dominate by radiative recombination ($\alpha=1$). An unambiguous trend cannot be extracted from the data. Supplementary Figure [2](#Sx1.F2)d shows the peak energy of the IX emission as determined from Gaussian fits, illustrating a strong $P_{ex}$-dependent blueshift beyond the threshold power $P_{ex}^{0}$, most likely reflecting the dipolar exciton-exciton interactions. For $P_{ex}\leq P_{ex}^{0}$ a much weaker blueshift is observed in our experiments. The simplest analysis of the blueshift based on a ”plate capacitor formula” would indicate that $\Delta E=(4\pi ne^{2}d/\kappa)$[2](#biba.bib2), where $n$ is the areal concentration of IX, $d$ is the effective separation of charge along the z-axis and $\kappa$ the effective dielectric constant. The righthand axis of Suppl. Fig. [2](#Sx1.F2)d shows the mapping of the observed blueshift onto the areal IX density using $d=$0.7\text{\,}\mathrm{nm} and $\kappa=4$. In more realistic models the effects of spatial inhomogenities in the IX density, finite temperature effects and multi-particle correlations must be properly taken into account3.

Supplementary Note 3: Photoluminescence excitation spectroscopy

We recorded the PL intensity of an example LIX peak (marked in Suppl. Fig. 3a) as a function of excitation energy. This dependence is shown in Suppl. Fig. 3b and exhibits two clear peaks in the spectral range of MoSe2 ($1.65\text{\,}\mathrm{eV}$) and WSe2 ($1.75\text{\,}\mathrm{eV}$) excitons. This behaviour was also observed for other LIX peaks, including peaks identified as multi-IX at higher excitation powers, and is consistent with free IXs4, 5, 6.

Supplementary Note 4: Time-resolved photoluminescence

The single- to few-exciton states discussed in detail in the main manuscript yielded only small PL signals on top of a non-neglible background from the the broad IX emission. It was thus not possible to conduct conclusive photon auto- or cross-correlation or time-resolved PL measurements on these states. In order to obtain a decent signal-to-background ratio of localized IX emission for time-resolved measurements, we instead performed experiments at higher excitation powers: We illuminated the sample with a PicoQuant LDH-P-C-640B laser diode driven by a PicoQuant PDL 800-B driver (wavelength: $640\text{\,}\mathrm{nm}$; Pulse width $<$90\text{,}\mathrm{ps}$$; Repitition rate: $2.5\text{\,}\mathrm{MHz}$) with a cw power equivalent of $43\text{\,}\mathrm{\SIUnitSymbolMicro W}$ (power per pulse: $17.2\text{\,}\mathrm{pJ}$). The resulting PL spectra on pillar A and a reference position in the flat HBL region are given in Suppl. Fig. 4a and b, respectively. At the pillar position, broad LIX emission, composed of several overlapping peaks, can clearly be seen on top of the free IX background. We cannot identify distinct states (and in particular exciton numbers) at these excitation powers anymore, but we still attribute this emission to IX trapped in the pillar-defined potential due to its spatially confined nature and the reduced energy.

We use bandpassfilters (Thorlabs FB930-10 / FB890-10; Center wavelength: $930\text{\,}\mathrm{nm}$ / $890\text{\,}\mathrm{nm}$; Transmission FWHM: $10\text{\,}\mathrm{nm}$) to select emission in the spectral range of the LIX or the free IX main peak. The transmission windows are indicated by the colored bars in Suppl. Fig. 4a and b. We record the time-resolved PL by measuring the time difference between signal events as recorded by a SPAD and the laser trigger with a PicoHarp 300. Supplementary Figures 4c-f show the time-resolved PL measured at the two different positions and spectral ranges. Each of the panels includes the instrument response as a gray dashed line and a fit to the data with a biexponential decay ($I(t)=A_{slow}e^{-t/\tau_{slow}}+A_{fast}e^{-t/\tau_{fast}}$) as a solid black line. The color codes of the data match the indicators in panel a and b. We give a summary of the derived parameters in Supplementary Table 1. The free IX PL recorded at the pillar position (with a detection volume far greater than the strain-defined potential) and at the reference position (panel e and f, respectively) show the same time dependence within the error intervals obtained from the fit.

The time-resolved PL response in the low-energy (LIX) regime looks very different on the pillar and on the reference position (panels c and d) at first glance; the fits to the data reveal, however, that the timescales of the biexponential decay are very similar. The reason for the different appearance of the curves in panels c and d is that the fast decay accounts for 24% of the emission on the flat HBL, whereas this value is only 5% at the pillar site. The general presence of a fast and a slowly decaying component in all four cases may point towards exciton-density-dependent, nonradiative decay channels, as also indicated by the sublinear power-dependence of the free IX intensity (see Suppl. Fig. 2). A definite conclusion about decay mechanisms would require additional measurements, such as time-resolved PL at different excitation powers. However, the very similar timescales support the attribution of the LIX peaks as stemming from interlayer excitons.

Supplementary Note 5: Comparison of interaction-induced blueshifts of free and localized IX

Supplementary Figures 5a and 5b show example spectra from the datasets that are presented in Fig. 1f and 1e, respectively. Supplementary Figure 5a displays example spectra at pillar A for different powers spanning more than three orders of magnitude. As for the dataset presented in Fig. 2, we observe a series of emission features in the window $1.28\text{\,}\mathrm{eV}1.35\text{\,}\mathrm{eV}$. As the excitation power increases from the lowest levels, the spectrum clearly evolves from a single sharp emission line to a multiplet of excitonic features. Since the overlap of the peaks in these measurements are enhanced due to the broader linewidths compared to the data presented in Fig. 2, a thorough quantitative analysis of the individual peaks was not performed for this measurement. Supplementary Figure 5b shows PL spectra recorded at a reference position inside the unstrained HBL region at comparable excitation levels. Clearly, the free IX emission shifts in a similar way to the LIX emission band in the limit of strong excitation. The doublet structure of the IX emission, clearly identified at the highest excitation levels consists of two emission features separated by $\sim$25\text{,}\mathrm{meV}$$, most likely arises from the spin-orbit split conduction band minima in the MoSe2.

In order to quantitatively compare the blueshifts of LIXs and free IXs, we fit the free IX emission at a pillar site outside the LIX spectral range using multiple Gaussians. This fit is used to estimate the free IX contribution inside the LIX spectral range, as illustrated in Suppl. Fig. 5c. We calculated the spectral ”center of mass” of the fitted IX and of the LIX (everything that was not accounted for by the fit). The calculated centers of mass over almost two orders of magnitude in power are superimposed onto a false-color plot of the excitation-power-dependent PL spectra in Suppl. Fig. 5d. It can be seen that both show a similarly increasing trend but the blueshift is more pronounced for the LIXs. This is expected since the dipolar interactions between IXs should be enhanced by the confinement inside a trapping potential.

Supplementary Note 6: Excitation-power dependence of LIX emission at additional pillar sites

During our studies, pillar A consistently showed the clearest signal with the smallest emission linewidths and we therefore focused on it for the quantitative analysis within the main manuscript. We present emission spectra of pillar B and C for different excitation powers in Suppl. Fig. 6. They confirm the trend of the LIX emission being distributed towards higher-energy peaks for increasing excitation powers, which we saw in pillar A. For pillar C, one can see that the main LIX peak at low powers is at $1.28\text{\,}\mathrm{eV}$, whereas the peak at $1.31\text{\,}\mathrm{eV}$ becomes more relevant at higher powers. Simultaneously, the two highest-power LIX peaks at $1.337\text{\,}\mathrm{eV}$ and $1.346\text{\,}\mathrm{eV}$ (indicated by arrows) only emerge as the excitation power increases. The spectra of pillar B show the least clear trend. Nonetheless, a peak at $1.33\text{\,}\mathrm{eV}$ (indicated by the arrow), blueshifted by $9\text{\,}\mathrm{meV}$ from the previous main LIX peak, appears towards the highest excitation powers.

Supplementary Note 7: Polarization-resolved photoluminescence measurements

We analyzed the linear polarization of the emission of several LIX peaks. Two example spectra of the same pillar for perpendicular polarizations are shown in Suppl. Fig. 7a. While the free IX background does not change notably, the LIX peaks show strong intensity variations. The intensities of the marked peaks as a function of detected polarization angle are shown in Suppl. Fig. 7b-d. All peaks show preferential emission towards some angle. These angles do not coincide between the different peaks. The polarization of the excitation laser does not have any influence, ruling out valley coherence as the underlying effect. The highest degree of polarization (68%) was observed for the lowest-energy peak. The reason for this polarization is not clear. We note that due to the sample geometry (heterobilayer forming a tent over a pillar) we may not probe the HBL in normal incidence.

Supplementary Note 8: Enhanced excitation-power-dependent data for Fig. 2

Supplementary Figure 8 presents additional LIX PL spectra at varying excitation powers from the dataset presented in Fig. 2a. These data plus two additional spectra acquired at each of the presented powers were used to extract PL peak intensities and positions presented in Fig. 2b and c.

Supplementary Note 9: Detailed power dependencies of multi-IX complexes

In order to verify our assignment of the single exciton peak (LIX1 in Fig. 2), we performed PL measurements at even lower excitation powers. The integrated peak intensities in these measurements are presented in Suppl. Fig. 9a. Measurements at lower powers than the ones presented did not yield any usable signal. It can be seen that for almost two orders of magnitude in excitation power, the intensity increases linearly with the excitation before finally saturation effects occur. This is a clear indication that the LIX1 stems from the recombination of a single exciton.

For some of the other peaks, we could only extract an intensity for few (or in the case of LIX5 no) excitation powers before the onset of saturation. We used intensity ratios to substantiate our peak assignments with data from the saturated power regime. In order to justify this procedure, we consider a case in which excitons are randomly trapped in the potential at a rate $G$, which (like the exciton generation rate) is proportional to the excitation power, and decay with a radiative lifetime $\tau$. Both parameters could in principle depend on the occupation number $N$ of the trap. The occupation probabilities $n_{N}$ of the N-exciton states are then described by the following set of rate equations:

[TABLE]

If trapping and recombination rate are independent of N, the steady-state solution for these equations is known to be the Poissonian distribution:

[TABLE]

For a more general case, the exact distribution can take a different form. Focusing only on the ratio of populations of two states with occupation numbers $N+1$ and $N$, a quick, induction-based proof can show that:

[TABLE]

This linear power dependence is robust against occupation-dependent capture and recombination rates and persists into the saturated regime (in contrast to the power law fits presented in Fig. 2b). We thus plot, as a function of excitation power, the ratios of the intensities of the multi-IX states with one state whose particle number we assume to be one below the examined state. This is shown in Suppl. Fig. 9b. In all cases, we obtain an approximately linear excitation-power dependence which is not influenced by the onset of saturation effects (marked by the gray dashed line), confirming our initial exciton number assignments.

Supplementary Note 10: Calculations of the exchange energy of the localized interlayer biexciton

We compare our experimental results for the energetic splitting of the observed biexciton states to a theoretical model of exchange interactions between bilayer excitons as outlined in detail in Suppl. Ref. 7. In this model, the tunnel exchange energy for interlayer biexciton complexes is computed in the absence of any external potential. The derivation is based on an Ansatz for the ground-state excitonic wavefunction of the (in-plane) relative electron-hole motion8. In our calculations, we adapted this model to take into account the presence of a weak symmetric harmonic confinement potential in order to model the strain-related trapping potential. For this, we added a harmonic trapping potential to the non-interacting part of the model Hamiltonian in Eq. (2).

Two limiting cases can be considered for the trial wavefunction: (i) For small interlayer distances $d$, the Hamiltonian resembles that of a two-dimensional hydrogen atom. The additional harmonic term then acts as an effective magnetic field that due to its low strength leaves the wavefunction unperturbed. (ii) For large interlayer distances, the Coulomb interaction of electron and hole can be expanded in a Taylor series in terms of the in-plane distance $\rho$ of the two carriers. Upon truncating this expansion at second order, this gives rise to a term proportional to the square of the in-plane distance $\propto\rho^{2}$, which can be absorbed into the harmonic confinement potential introduced to the problem before. In doing so, we can closely follow the remaining steps of the derivation detailed in Suppl. Ref. 7 to evaluate the exchange energy of the interlayer biexciton.

Note that we implicitly assumed the validity of a couple of approximations for the derivation of binding and exchange energies presented in this work. Most notably, we assumed that (i) the effective electron and hole masses are equal, (ii) in-plane electrostatic interactions are well-described by Coulomb interaction potentials and (iii) image-charge effects are negligibly small. Further refining our model in order to take into account effects beyond the approximations (i)-(iii) is expected to result in more accurate results, of course. Despite these approximations, however, we expect the derived values in the main text to provide us with the correct orders of magnitude for the energy scales of interest. This is supported by the good agreement between theoretical predictions and performed measurements. Further note that, by construction, the employed variational Ansatz yields an upper limit for the exchange energy.

Supplementary Note 11: Calculations of the binding energy of the localized multi-IXs

The binding energies of trapped multi-excitonic states were calculated classically in terms of point-like dipoles in an isotropic 2D harmonic potential. Note that we define positive binding energies as an increase of the multi-particle-state energy compared to the energy of the well-separated constituent particles. In order to minimize dipole-dipole repulsion, the IX were maximally spaced around an equipotential of the trapping potential. The constellations that arise from this condition are a line (biexciton), an equilateral triangle (triexciton), a square (quadexciton) and a regular pentagon (quintexciton), all centered around the minimum of the harmonic potential, as illustrated in Suppl. Fig. 10. As a function of the nearest-neighbor distance $\Delta\rho$, the harmonic and dipole-dipole potential energies are then given by:

[TABLE]

, where $M$ is the total mass of the exciton and $\ell$ is the same effective confinement lengthscale defined in the main text. We minimize each of these expressions to find the equilbrium interparticle separation and equilibrium potential energy. Introducing $\tilde{E}=\left(\frac{d^{2}}{4\pi\epsilon}\right)^{2/5}\left(\frac{\hbar^{2}}{M\ell^{4}}\right)^{3/5}$, we find:

[TABLE]

Combined with the energy of a single exciton inside the trap $E^{X}$, the transition energies of the multi-exciton states can then be expressed as:

[TABLE]

Due to the universal scaling with both, the strength of the harmonic potential and the dipole-dipole interaction, the ratios of blueshifts with respect to the single-exciton transition take fixed values:

[TABLE]

As discussed in the main text, we observe the localized bi-, tri-, quad- and quintexciton emission to be blueshifted by $8.4\pm 0.6\text{\,}\mathrm{meV}$, $12.4\pm 0.4\text{\,}\mathrm{meV}$, $15.5\pm 0.6\text{\,}\mathrm{meV}$ and $18.2\pm 0.8\text{\,}\mathrm{meV}$ from the single exciton, respectively. Evaluating the blueshift ratios introduced above for these values, we obtain:

[TABLE]

The predictions from the model show good quantitative agreement, especially considering that the calculations are strictly classical and the exact trapping potential shape of our sample is not known. We emphasize that other spatial configurations of the excitons can easily be shown to be less favorable than the ones discussed above. For example, arranging three excitons in a line would result in $U^{3X}_{0,\textrm{line}}\approx 3.496\tilde{E}>U^{3X}_{0}$, a quadexciton with three excitons on a triangle and the fourth one in the center would give $U^{4X}_{0,\textrm{triangle}}\approx 5.493\tilde{E}>U^{4X}_{0}$ and a quintexciton consisting of four excitons on a square surrounding the fifth in the trap center would have $U^{5X}_{0,\textrm{square}}\approx 7.845\tilde{E}>U^{5X}_{0}$.

In addition to the energies of multi-excitonic states, the model also predicts equilibrium distances of the constituent excitons. These distances are not directly accessible in our experiments but we use the biexciton equilibrium distance $\Delta\rho_{0}^{2X}=\left(\frac{12d^{2}M\ell^{4}}{4\pi\epsilon\hbar}\right)^{1/5}$ as $\Delta\rho$ in the calculations of the biexciton exchange energy.

Supplementary Note 12: Thermal cycling

We investigated the effect of thermal cycling on the localized IX by repeatedly cycling the sample temperature between $10\text{\,}\mathrm{K}$ and $300\text{\,}\mathrm{K}$. Each time, we performed PL measurements at $10\text{\,}\mathrm{K}$ under the same excitation conditions (cw excitation with a HeNe laser at $633\text{\,}\mathrm{nm}$, $30\text{\,}\mathrm{\SIUnitSymbolMicro W}$ focused to a spot of $1\text{\,}\mathrm{\SIUnitSymbolMicro m}$). Supplementary Figure 11 compares typical emission spectra from the same pillar following successive cooldowns. We observe a clear change in the energy at which the LIX features appear. This effect was found to be even more pronounced than with typical strained WSe2-monolayer samples tested under similar conditions. In order to mitigate the effects of thermal cycling, we performed the comprehensive power-dependent studies presented in the main manuscript (data in Fig. 2) with the sample placed inside a closed-cycle cryostat that facilitated all measurements during a single cool down. Under these conditions, we found the absolute energy of the LIX states to be stable over several weeks.