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

Interlayer excitons (IXs) in heterobilayers (HBLs) 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 emergence of quantum correlations in finite 2D systems may open up tantalizing prospects of excitonic Bose-Einstein condensation and Wigner-crystal-like phases at sufficiently low temperatures and in confined geometries. Here, we trap a tunable number of dipolar IXs within a nanoscale confinement potential induced by placing a MoSe$_2$-WSe$_2$-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. The intensities of these features exhibit characteristic near linear, quadratic, cubic and quartic power dependencies, which allows us to identify them as different multiparticle configurations with $N_{IX}=1-4$. We directly measure the hierarchy of dipolar and exchange interactions as $N_{IX}$ increases. The interlayer biexciton 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) and quadexcitons ($(17.0\pm0.5)$ meV). These values are shown to be mutually consistent with numerical modelling of dipolar excitons moving in a harmonic trapping potential having a confinement length scale in the range $l\approx3$ 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 IXs in TMD-HBLs.

Here, we trap a tunable number of dipolar IXs (N IX ∼ 1 − 4) within a nanoscale confinement potential induced by placing a MoSe 2 -WSe 2 hetero-bilayer (HBL) onto an array of SiO 2 nanopillars [18,19]. 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. The intensities of these features exhibit characteristic near linear, quadratic, cubic and quartic power dependencies, [20] which allows us to identify them as different multiparticle configurations with N IX ∼ 1 − 4. We 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 ∆E = (8.4 ± 0.6) meV and split by 2J = (1.2 ± 0.5) meV. The blueshift is even more pronounced for triexcitons ((12.4 ± 0.4) meV) and quadexcitons ((17.0 ± 0.5) meV). These values are shown to be mutually consistent with numerical modelling of dipolar excitons moving in a harmonic trapping potential having a confinement lengthscale in the range ≈ 3 nm. [21,22] Our results contribute to the understanding of interactions between IXs in TMD heterobilayers at the discrete limit of only a few excitations and represent a key step towards exploring quantum correlations between IXs in TMD hetero-bilayers.
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 orienta-tion of the basal plane of the component layers. [7] Such atomically-thin 2D heterointerfaces can host IXs -Coulomb-bound states of electrons (e) and holes (h) located in the different component monolayers. [7] The spatial separation of e and h gives rise to a significantly increased lifetime compared to intralayer excitons [7,9,23] and a static electric out-of-plane dipole moment with a magnitude |p z |/e ∼ 0.7 nm for the MoSe 2 -WSe 2 heterointerface [7].
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 energy [7,24,25] and lifetime [24,25] 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 system [7,9,23]. It has also been demonstrated that applying strain to the HBL can strongly tune the IX energy [26]; a method that is routinely used to define exciton trapping potentials in monolayer TMDs. [18,19,27] Overall, IXs constitute a highly tunable platform for optically exploring and controlling interacting gases of dipolar excitons in the solid-state.
Our sample structure is schematically depicted in fig. 1a. It consists of stacked monolayers of WSe 2 and MoSe 2 that form a HBL which is then transferred onto a SiO 2 substrate patterned with nanopillars having a height of ∼ 90 nm and diameter ∼ 130 nm (see Methods). Fig. 1b shows an optical micrograph of the sample in which the monolayer MoSe 2 , WSe 2 and HBL regions can be seen, as well as the nanopillars in the HBL region. The contours of the WSe 2 and MoSe 2 monolayers are outlined with blue and orange dashed lines, respectively, and pillars that did not pierce the TMD layers appear as black dots. Fig. 1c shows normalized spectra recorded from three different unstrained sample regions, namely the WSe 2 monolayer, the MoSe 2 monolayer and the WSe 2 -MoSe 2 -HBL away from the nanopillars. These spectra were recorded at T = 10 K with cw excitation at 633 nm and a power of P ex =35 µW focused to ∼ 1 µm. Both, WSe 2 and MoSe 2 , have prominent emission features at their characteristic exciton and trion energies at 1.75 eV and 1.72 eV, [28] and 1.66 eV and 1.63 eV, [29] respectively. Furthermore, the WSe 2 component monolayer exhibits a broader band of unresolved emission at lower energies (1.60 − 1.70eV) that has been identified as stemming from dark excitons [30,31] and trions, [32] charged biexcitons [33][34][35][36][37] and defect-bound excitons. [38,39] Each of these emission features are also present in the HBL region, albeit strongly quenched, due to interlayer charge transfer that occurs Besides the power-dependent measurements shown in e,f all spectra were recorded at 10 K subject to strong continuous wave (cw) excitation at 633 nm using a power of 35 µW focused to a spot-size of ∼ 1 µm.
over sub-picosecond timescales, faster than the exciton lifetime. [40][41][42] The most prominent emission feature from the HBL is at ∼1.38 eV and is attributed to IXs formed by an electron in the MoSe 2 and a hole in the WSe 2 . [7] This attribution is supported by the emission energy and the characteristically asymmetric lineshape featuring red-shifted emission from momentum-indirect IXs [23]. Fig. 1e shows the power dependence of the IX emission revealing a significant blueshift (≥ 20 meV) arising from repulsive dipolar interactions in the gas of IX for increasing excitation levels [5,9]. We note that this blueshift does not occur over the whole range of excitation levels, but rather exhibits an onset close to P ex ∼ 30 nW, marking a low-density regime where dipolar interactions are weak (see supplemental - Fig. S2).
Furthermore, the low density IX regime also manifests itself by a linear power dependence of the IX peak intensity that undergoes a transition to a √ P ex dependence at higher power, indicative of non-radiative exciton-exciton annihilation at higher densities [43].
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 eV can still be observed in each emission spectrum since the laser focal volume is much larger than the nanopillar size. However, additional peaks emerge in the vicinity of the nanopillar sites, red-shifted by up to ∼100 meV from the peak of the free IX emission.
The observation of a pronounced red-shift is similar to reports of strain-induced quantum emitters in monolayer WSe 2 [18,19] and MoSe 2 [27]. We attribute these features to local- as the excitation level increases. We note that the center of the spectral weight of the LIX emission shifts in a way that closely mimics the free IX peak (see supplemental information - fig. S3) 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.
Since pillar A exhibited the brightest LIX emission, we recorded detailed excitation-power dependent measurements for much lower excitation levels than presented in fig. 1 on it. The results of these measurements are presented in Fig. 2a for ultra-low excitation powers ranging from 79nW to ≥ 500nW, two orders of magnitude weaker than the power used to record . These observations mirror previous findings for multi-exciton states in semiconductor quantum dots [45] and double quantum wells [14] and are consistent with the Poisson statistics governing the independent capture of excitons from a reservoir into the trapping potential. of the biexciton LIX 2 doublet is significantly blueshifted by (8.4 ± 0.6) meV with respect to LIX 1 and the two components are separated by ∼ (1.2 ± 0.5) meV.
In order to quantitatively understand the strong blueshift of the LIX biexciton from the 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 = /m * Ω 0 . Fig. 3a illustrates how we calculate the equilibrium separation ∆ρ between two IX by minimizing the sum of repulsive (dipolar) and attractive (harmonic) potential energies. For values of ∼ 3 nm-4 nm, reference to fig. 3a shows that ∆ρ ∼ 5 nm-6.5 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. [22] and start with a Hamiltonian describing the biexciton complex of the following form: where describes non-interacting excitons exposed to the weak harmonic confinement potential and represents the Coulomb interactions between the excitons. In equations 2 and 3, the co- 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 , giving rise to two, energetically distinct biexciton states separated by the exchange energy 2J(∆ρ). Here, J(∆ρ) is given by [50] where Ψ(x, y) is the ground state of eqn.
(1) that has been transformed to the coordinate system x = (ρ 1 − ρ 2 − ∆ρ)/ √ 2 and y = (ρ 1 + ρ 2 )/ √ 2. From this expression, we numerically calculate the exchange splitting to obtain the absolute energies of the two biexciton states:  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 K and a HeNe laser (λ = 633 nm) focused to a diameter of ∼ 1 µm. In order to counteract thermal cycling effects (see Supplementary), the measurements presented in fig. 2 were carried out in a setup equipped with a closed-cycle cryostat at T = 4 K. The excitation laser used in this setup was at λ = 532 nm and was defocused due to chromatic aberration in the used objective.
We note that the data presented in figs. 1f and 2 were both recorded from pillar A, but following thermal cycling 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 ∼ ±10 meV. Nevertheless, following each cool-down, the generic form of the group of peaks was retained. The supplemental fig. S4 shows power dependent µP L spectra recorded from pillar A (from the same measurement as the data presented in fig. 1f), illustrating how the single sharp emission feature is observed close to 1.292 eV at the lowest excitation power investigated, with additional peaks emerging in the range ∼ 5 meV-30 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 meV around the respective center energies. The only exception is the lowest-energy peak for which fitting was possible over the whole power range due to its spectral isolation. The data points for LIX 3A and LIX 3B in fig. 2b initially overlapped. Thus, the data for LIX 3A was slightly offset for better visibility. 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 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) and a square (quadexciton), all centered around the minimum of the harmonic potential. As a function of the nearest-neighbor distance ∆ρ, the harmonic and dipole-dipole potential energies are then given by: , where M is the total mass of the exciton and 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Ẽ = d 2 , we find: 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: Due to the universal scaling with both, the strength of the harmonic potential and the dipoledipole interaction, the ratios of blueshifts with respect to the single-exciton transition take fixed values: As discussed in the main text, we observe the localized bi-, tri-and quadexciton emission to be blueshifted by (8.4 ± 0.6) meV, (12.4 ± 0.4) meV and (17.0 ± 0.5) meV from the single exciton, respectively. Evaluating the blueshift ratios introduced above for these values, we obtain: The predictions from the model are slightly below our experimentally determined ratios but 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

S2. 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 Ref. [22]. 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 motion [51]. 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 ρ 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 ∝ ρ 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 Ref. [22] 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.

S3. Thermal cycling
We investigated the effect of thermal cycling on the localized IX by repeatedly cycling the sample temperature between 10 K and 300 K. Each time, we performed PL measurements at 10 K under the same excitation conditions (cw excitation with a HeNe laser at 633 nm, 30 µW focused to a spot of 1 µm). Fig. S1 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 WSe 2 -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. with exciton-exciton annihilation at elevated excitation levels [43]. Careful examination of Fig.S2c for P ex ≤ P 0 ex ≈ 20 nW shows a higher power factor, closer to unity, consistent with predominant radiative recombination. Figure S2d shows the peak energy of the IX emission, illustrating a strong P ex dependent blueshift beyond the threshold power P 0 ex , most likely reflecting the dipolar exciton-exciton interactions. For P ex ≤ P 0 ex a much weaker blueshift is observed in our experiments. Whilst the simplest analysis of the blueshift based on a "plate capacitor formula" would indicate that ∆E = (4πne 2 d/κ) [52], where n is the areal concentration of IX, d is the effective separation of charge along the z-axis and κ the effective dielectric constant. The righthand axis of Fig. S2d shows the mapping of the observed blueshift onto the areal IX density using d = 0.7 nm and κ = 4. In more realistic models the effects of spatial inhomogenities in the IX density, finite temperature effects and a b c FIG. S3. Excitation-power-dependet IX and LIX PL spectra. All spectra in this figure were recorded at T = 10 K under excitation with a HeNe laser at 633 nm focused to a spot of ∼ 1 µm diameter. a, Example PL spectra at varying excitation powers recorded on pillar A. b, Example PL spectra at similar powers as presented in S3a recorded on a reference position on the unstrained HBL. c, False color plots of the normalized PL spectra on pillar A and the same reference position as in S3b as a function of excitation power. multi-particle correlations must be properly taken into account [5].  fig. 1f and 1e, respectively. Fig. S3a 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 emission features in the window 1.29 eV-1.35 eV. As the excitation power increases from the lowest levels, the spectrum clearly evolves from a single sharp emission line to a multiplet of excitionic 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. However, we note that the global shift of the spectral center of weight of the LIX and IX peaks becomes similar in the limit of large IX density (P ex ≥ 10 µW). Fig. S3b shows PL spectra recorded from reference position inside the unstrained HBL region at comparable excitation levels. The rightmost panels of Fig. S3 show the same date plotted as a colormap, with the maximum intensity normalized to unity. 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 ∼ 25 meV, most likely arises from the spinorbit split conduction band minima in the MoSe 2 .