Capacitively-coupled and inductively-coupled excitons in bilayer MoS$_2$

The interaction of intralayer and interlayer excitons is studied in a two-dimensional semiconductor, homobilayer MoS$_2$. It is shown that the measured optical susceptibility reveals both the magnitude and the sign of the coupling constants. The interlayer exciton interacts capacitively with the intralayer B-exciton (positive coupling constant) consistent with hole tunnelling from one monolayer to the other. Conversely, the interlayer exciton interacts inductively with the intralayer A-exciton (negative coupling constant). First-principles many-body calculations show that this coupling arises via an intravalley exchange-interaction of A- and B-excitons.

The elementary excitation at energies close to the band gap of a semiconductor is the exciton, a bound electronhole pair. The exciton is a prominent feature of the linear optical response of a two-dimensional (2D) semiconductor, for instance monolayer MoS 2 , even at room temperature [1,2]. This is a consequence of the giant exciton binding energies, hundreds of meV, in these materials [3]. Excitons interact with each other. These interactions lead to nonlinear optical effects [4,5], potentially useful in opto-electronics, and to condensation phenomena, for instance the creation of states with macroscopic quantumcorrelations [6]. At the few-exciton level, exciton-exciton repulsion in an appropriate trap is a potential way to engineer a single-photon emitter [7,8].
Here, we probe exciton-exciton interactions in gatedhomobilayer MoS 2 . This is a rich system for probing exciton-exciton interactions: an interlayer exciton (IE) interacts with both the B-and A-excitons of the constituent layers. We show that the optical susceptibility reveals not just the magnitude but also the sign of the exciton-exciton coupling. Remarkably, the IE-A and IE-B couplings have opposite signs. The IE-B coupling has a positive sign due to hole tunnelling from one monolayer to the other. Conversely, the IE-A coupling has a negative sign, showing that the coupling mechanism has a different microscopic origin. Beyond density functional theory (post-DFT) calculations enable us to ascribe the IE-A coupling to spin, specifically an exchange-based hy- * lukas.sponfeldner@unibas.ch bridisation of A-and B-excitons.
The device is constructed using bilayer MoS 2 and sketched in Fig 1a. The naturally 2H-stacked bilayer MoS 2 is embedded in h-BN; few-layer graphene layers provide a back-contact and a top-gate. In addition, the MoS 2 layer is contacted. Results are presented as a function of the vertical electric field (F z ) for the smallest possible electron density (n). The device is illuminated locally with a very weak broadband source. The reflectivity is measured, using the response at large n as reference, and converted into the optical susceptibility using a Kramers-Kronig relationship (see supplement of Ref. [24] for details). At F z = 0, the imaginary part of χ (which determines the absorption) shows three peaks: at low energy the intralayer A-excitons, at high energy the intralayer B-excitons, and in between, the interlayer exciton (IE), see Fig. 2b. These momentum-direct excitonic transitions are shown schematically in Fig. 1a; the relevant band structure at the K-point is shown in Fig. 1c and Fig. 1d. The current understanding is that the IE consists of a hole state delocalised across the two layers bound to an electron localised in one of the two layers [25,26]. For F z = 0, the IE splits into two lines [27][28][29]. At high F z , there is a clear avoided crossing between the IE upper branch and the B-exciton; and a weak avoided crossing between the IE lower branch and the A-exciton, as reported in Ref. [27].
We focus here on the avoided crossings. At the IE-B avoided crossing, one of the B-excitons couples to the IE, the other does not. In Fig. 2a and Fig. 2b, we subtract the peak arising from the uncoupled B-exciton. A simple two-level model describes the peak energies convincingly but not the relative intensities. Strikingly, at the F z for which the two transitions are closest together in energy (F z,0B ), the intensities are quite different: the lower-energy transition is considerably stronger than the higher-energy transition. The IE-A coupling shown in Fig. 2c is weaker than the IE-B coupling. One of the A-excitons couples to the IE, the other does not. The absorption from the uncoupled A-exciton is substracted in Fig. 2c. The minimum energy separation of the peaks is comparable to the peak broadening. Nevertheless, the avoided crossing has a strong effect on the intensities: as F z increases, the IE-like branch enters the avoided crossing with a relatively large intensity, but emerges surprisingly with a much lower intensity. Building on recent studies [27][28][29], we show here that the IE can be tuned energetically below the A-exciton. The IE-A intensity behaviour mimics the behaviour at the IE-B avoided crossing but with one crucial difference. The upper-IE is strong on the low-energy side of the B-exciton; the lower-IE is strong on the high-energy side of the A-exciton. Our target is to understand the origin of the different coupling behaviours.
The optical susceptibility of a quantum well can be calculated from the semiconductor Bloch equations [30][31][32][33]. In this approach, the quantum well is treated quantum mechanically. The final result is identical to a completely classical approach in which the quantum well is treated as a 2D array of optical dipoles [34]. Inspired by the success of the purely classical approach, we set up a heuristic description of the IE-B avoided crossing as sketched in Fig. 1b. The IE and B-excitons are treated as optical dipoles, each driven by the same driving field E, but with different oscillator strengths F 1,2 [25,35]. A linear coupling term κ is introduced: an IE-dipole induces a B-dipole, and vice versa. Solving the equations of motion of this system yields an analytic equation for both the eigenenergies Ω ± and the absorption strengths of the two eigenmodes (see Supplement). The calculated absorption strength I 1,2 of each eigenmode depends on the energy detuning, the coupling κ, and the oscillator strengths F 1,2 .
The classical model reproduces the experimental results extraordinarily well provided the ratio of the oscillator strengths and the coupling constant are well chosen. The measured peak energies (Fig. 2a) and integrated absorption (Fig. 2d) of the IE-B avoided crossing are fitted by the coupled dipole model (light-and dark-blue dashed lines). The fits yield an oscillator strength ratio F 1 /F 2 ≈ 5 and an IE-B coupling strength of κ = +35.8 ± 3.6 meV. The IE-A avoided crossing is described with the same model but with the energies appropriate to the lower IE-branch and the A-exciton along with a different choice of coupling constant. For the Aexciton energy, a quadratic Stark shift is included [36]. Our model reproduces also the IE-A avoided crossing very convincingly, capturing the F z -dependence of the intensities (see Fig. 2c and Fig. 2e). The fits yield an oscillator strength ratio F 1 /F 2 ≈ 5 and an IE-A coupling strength of κ = −3.5 ± 0.4 meV. The absorption spectra for the IE-A and IE-B couplings are calculated separately and added together to describe the full F z -dependence, as shown in Fig. 3. The calculated absorption reproduces very closely the experiments (Fig. 2b).
This model unearths a crucial result: the IE-A coupling constant is of opposite sign to the IE-B coupling constant. The sign of the coupling constant is particularly important in determining the relative strengths of the absorption peaks (see Supplement). We interpret the sign of the coupling by making an analogy to driven, coupled RLC-circuits (see Supplement). Two such circuits can be coupled via an impedance. The equations of motion are analogues of those describing the driven optical dipole. The nature of the coupling impedance determines the sign of the coupling constant: coupling via a capacitance results in a positive coupling constant; coupling via an inductance results in a negative coupling constant. In this analogy, the IE-B coupling corresponds to a capacitive coupling. This suggests that, at a microscopic level, FIG. 2. Absorption (σ + polarisation) of homobilayer MoS2 as a function of applied electric field Fz. (b) Absorption over the whole studied energy range. The green and purple rectangles correspond to the measurement regions centred around the Bexciton (a) and the A-exciton (c), respectively. The measurements were carried out at T = 4.2 K and magnetic field Bz = 9 T. For clarity, the non-interacting BL1 exciton is removed from (a) and (b) and the non-interacting AL2 exciton is removed from (b) and (c) in a data processing step (see Supplement). In (a), Fz,0B marks the energy crossing for zero coupling. The coloured dashed lines in (a) and (c) are fits to the extracted peak energies according to equations S20 and S47 for the eigenenergies Ω±, respectively (see Supplement). (d) Integrated absorption extracted from the spectra shown in (a). The coloured dashed lines are fits to the absorption strength of the eigenmodes I1,2 according to equation S26. (e) Integrated absorption extracted from the spectra shown in (c). As the coupling strength is much smaller than the linewdith of the A-exciton, only the spectra far from the crossing point can be fitted unambiguously. The coloured dashed lines are model fits with the same function as in (d) including the quadratic Stark shift of the A-exciton in the initial equations of motion (see supplement section IV). All parameters extracted from the measured data are summarised in SI Table I is consistent with hole tunnelling from one layer to the other [25]. Conversely, the IE-A coupling corresponds to an inductive coupling. This points to a completely different coupling mechanism, as discussed below.
The model provides an explanation for the F zdependent absorption strengths. We take the IE-A coupling as an example. Without a coupling, both IE and A respond directly to the driving field. A has the stronger response: the induced dipole moment is in-phase with the drive for energies well below the bare A-energy and out-of-phase for energies well above the bare A-energy (the standard behaviour for a driven harmonic oscillator). With a coupling, each eigenmode is a dressed state of IE and A. At detunings far from the avoided crossing, there is an A-like and an IE-like eigenmode. The IE-like mode is driven by two sources: the field acting directly on the IE, and via its coupling to A. Now the sign of the coupling plays a crucial role. If the coupling has a negative sign, then for energies above (below) the bare A-energy, these two terms have the same (opposite) sign and interfere constructively (destructively). At a particular detuning, the destructive interference is complete and the absorption of one mode disappears. Even at energies far from the avoided crossing, this interference has a strong effect on the IE-absorption. The picture inverts for a positive coupling, the IE-B interaction. In this case, the IE-like mode is boosted (suppressed) when it lies below (above) the bare B-energy. At F z = 0, the IE is far from the avoided crossing with both A-and B-excitons but its absorption strength is boosted by its coupling to both A and B. In simple terms, the dielectric constant at the IE-resonance is strongly influenced by the strong Aand B-resonances.
We look for a microscopic explanation for the different IE-A and IE-B couplings. To do this, we describe the band structure of bilayer MoS 2 at the GW-level (one-particle Green's function G, dynamically screened Coulomb interaction W), and use these states to construct excitons by solving the Bethe-Salpeter equation BSE (see Fig. 4 and Supplement). The results describe the general behaviour of the experiments very well as revealed by a comparison of the calculated relative absorption strengths in Fig. 4b and Fig. 4c with the measured integrated absorption in Fig. 2e. (Exact quantitative agreement is not expected as the model assumes that the MoS 2 bilayers are located in vacuum -it does not take into account the full dielectric environment of the sample.) As in the experiment, the IE are relatively strong when they lie energetically between the bare Aand the bare B-resonances but weaker when they lie out of this energy window ( Fig. 4b and Fig. 4c). The post-DFT results can be parameterised by the coupled-dipole model (see Supplement).
An analysis of the spin and orbital composition of the excitons provides an explanation for both the IE-B and IE-A couplings. The IE-B coupling arises via hole tunnelling (see Fig. 1c). In the bare B-state, the electron and hole are localised in one monolayer; in the bare IEstate, the electron is localised in the same monolayer but the hole is localised in the other layer. Hole tunnelling couples these two states (Fig. 1c). Consequently, the excitons observed in the experiment are mixed excitonic states rather than states consisting of a single excitonic transition. The large IE-B coupling constant, +35.8 meV, reflects the efficient hole tunnelling. Conversely, the IE-A coupling arises via a weak admixture between the A-and B-excitons in the same valley and layer, as shown schematically in Fig. 1d. This exciton admixture was proposed in Ref. [37] for MoS 2 monolayers and is found also in our bilayer calculations (Fig. 4). This means that both IE and A couple to B. IE and A then acquire an effective coupling, a second-order effect.
In this analysis, depicted in Fig. 1d, the IE-A coupling is determined by (with the sign convention employed here) where κ IE-B is the IE-B coupling, κ A-B the A-B coupling and ∆E A-B is the energy splitting between A and B [38]. The experiment measures κ eff = −3.5 meV, κ IE-B = 35.8 meV and ∆E A-B = −170 meV. In turn, this determines κ A-B . We find κ A-B = −16.6 ± 2.5 meV. In other words, by using the IE as a tunable probe, we are able to determine the A-B coupling energy. This is an important quantity -it arises even in a MoS 2 monolayer and determines to what extent spin is a good quantum number in the fundamental exciton [37,39]. We state two main conclusions. First, in homobilayer MoS 2 , the IE and B-excitons couple via hole tunnelling; the IE and A-excitons couple via an exchange-induced A-B admixture. The experiment enables the A-B coupling to be determined quantitatively. We find −16.6 meV. Its existence implies that spin is an imperfect quantum number for the A-and B-excitons in the same valley. Second, a measurement of the optical susceptibility enables not just the magnitude but also the sign of the excitonexciton couplings to be determined.
As a final comment, we point out that the driven coupled-oscillator model shows that a weak resonance can be made visible by bringing it into near-degeneracy with a strong resonance. All that is required is a coupling between the two resonances. This point is not restricted to driven dipoles at optical frequencies: it applies quite generally.
We thank Christoph Bruder for fruitful discussions and Jonas Roch for expert experimental assistance. Basel acknowledges funding from the Swiss Nanoscience Institute (SNI), PhD School Quantum Computing and Quantum Technology (QCQT), SNF (project 200020 175748), and NCCR QSIT. Toulouse acknowledges funding from ANR MagicValley, ANR IXTASE, ANR HiLight, ITN 4PHOTON Marie Sklodowska Curie Grant Agreement No. 721394 and the Institut Universitaire de France.
K.W. and T.T. acknowledge support from the Elemental Strategy Initiative conducted by the MEXT, Japan (Grant Number JPMXP0112101001) and JSPS KAK-ENHI (Grant Numbers 19H05790 and JP20H00354). I.C.G. thanks the CALMIP initiative for the generous allocation of computational time, through Project No. p0812, as well as GENCI-CINES and GENCI-IDRIS, Grant No. A010096649.

I. SAMPLE FABRICATION
The studied van der Waals heterostructure was assembled with a dry-transfer technique [1].
The individual layers were mechanically exfoliated from bulk crystals (natural MoS 2 crystal from SPI Supplies, synthetic h-BN [2], and natural graphite from Graphene Supermarket) onto SiO 2 (300 nm)/Si substrates. The flakes were picked up and stacked with a stamp made of polydimethylsiloxane (PDMS) covered with a thin layer of polycarbonate (PC). The naturally 2H-stacked bilayer MoS 2 and few-layer graphene (FLG) layers were contacted by Ti/Au (5 nm/45 nm) electrodes defined by electron-beam lithography. Using atomic force microscopy (AFM), the top and bottom h-BN thicknesses were determined to be d T = 16.2 nm and d B = 21.5 nm, respectively. The studied sample is the same as device 1 in Ref. [3] albeit in a second cooldown cycle. An exception is the measurement shown in Main Fig. 2b and SI Fig. 10b which was recorded in the first cooldown cycle.
The dual-gating scheme allows for independent control of the electric field perpendicular to the MoS 2 plane F z and the total charge carrier density n in the MoS 2 bilayer. F z and n can be changed by applying voltages V TG and V BG to the top and bottom FLG gates while the MoS 2 bilayer is electrically grounded. Assuming no free charge carriers in the bilayer, the electric field across the MoS 2 bilayer can be written as with the dielectric constant of h-BN hBN ≈ 3.76 and of bilayer MoS 2 BL ≈ 6.8 [4]. For a more detailed description of the electrostatic model see Supplementary Information of Ref. [3]. In the experimental electric field sweeps shown in the main paper, the applied voltages only change F z . The carrier density in the MoS 2 bilayer stays constant and at a low value. This can be deduced from the stable absorption strength at the A:1s resonance.

II. EXPERIMENTAL SETUP
Absorption measurements were conducted with the setup sketched in SI Fig. 1. A home-built confocal microscope facilitates a local measurement of the sample's optical susceptibility at a temperature of 4.2 K. White (Osram warm white) light is coupled into a multi-mode fiber and routed to the excitation arm of the microscope. A combination of a linear polariser and a quarter wave plate (λ/4) is used to create a circularly polarised excitation beam. With the help of a computer controlled liquid-crystal (LC) variable retarder, the excitation light can be either set to the σ + or the σ − circular polarisation state. The incident beam is then focused on the sample surface using a microscope objective (NA = 0.45). Piezoelectric nano-positioners are used to control the position of the spot on the sample. The reflected signal from the sample is collected through the same microscope objective and coupled into a single-mode fibre at the end of the collection arm of the confocal arXiv:2108.04248v1 [cond-mat.mes-hall] 9 Aug 2021 Sketch of two masses m1 and m2 coupled by a spring with spring constant κ. m1 has a spring constant of k, m2 has a tunable spring constant of k + ∆k. Each mass is driven by an external force F1,2 with different amplitudes. x1,2 describe the displacements of each mass.
microscope. The collection fiber is connected to a spectrometer where the reflected signal is dispersed by a 600 g/mm (IE-B interaction measurement), a 1500 g/mm grating (IE-A interaction measurement), or a 300 g/mm grating (Main Fig. 2b and SI Fig. 10b) and detected by a liquid nitrogen-cooled charge coupled device (CCD) camera.
The absorption of the bilayer MoS 2 can be described by the imaginary part of the optical susceptibility Im(χ). Im(χ) is determined from the differential reflectivity (R − R 0 )/R 0 = ∆R/R 0 using a Kramers-Kronig relation (see Supplementary Information of Ref. [5] for a detailed calculation). R is the reflected signal from the bilayer MoS 2 and R 0 is a reference reflectivity spectrum. For R 0 we interpolate the raw reflectivity spectrum at high electron densities where the total absorption strength is spread over a large spectral range [6].

III. EXCITON INTERACTION MODEL
The excitonic interactions are modelled as the coupling of two driven oscillating optical dipoles. The dipole oscillation can be described by the oscillation of the dipole extension r or by the oscillation of the charge q. One can think of two frameworks which provide a classical analogue. In a mechanical system, two masses represent the dipoles, a spring represents the coupling. In an electrical system, two RLC circuits represent the dipoles, an impedance (either capacitor or inductor) represents the coupling. Both, the mechanical and the electrical system, can be described by equivalent equations of motion as shown below.

A. Equations of motion
First, we look at a mechanical system. We consider two driven and damped harmonic oscillators with a mass m 1 and m 2 and a spring constant of k 1 = k and k 2 = k + ∆k, respectively. The term ∆k allows for the detuning of one oscillator with respect to the other. The two oscillators are coupled by a spring with spring constant κ. The damping of each oscillator is described by γ 1 and γ 2 . The oscillators are externally driven by the forces F 1 and F 2 , respectively. A sketch of the mechanical system is shown in SI Fig. 2.
In the mechanical system, an oscillating force F 1,2 induces a displacement x 1,2 of the two masses m 1,2 . For the following derivation, we assume m 1 = m 2 = m and take F 1,2 = F 1,2 E 0 Re(e iωt ). F 1,2 is the mechanical oscillator strength of each mass. The matrix form of the equations of motion for this system is

B. Decoupling the equations of motion
For an analytical solution to equation S2, we need to transform the equations of motion into the eigenvector basis. The two eigenmodes r 1,2 are given by [7]  where the rows of the transformation matrix U are the normalised eigenvectors of the matrix in equation S2. For simplicity, only the non-normalised eigenvectors will be shown in the following derivation. All calculations and plots shown in the main text and the Supplement use the normalised eigenvectors. The eigenvectors take the form 3 Introducing the ratio ρ between the detuning and coupling constant ρ = ∆k κ , we can write We can write the transformation matrix U as The angular eigenfrequencies of the system are given by (S8) It is important to note that Ω ± do not depend on F 1,2 .
After transforming equation S2 with the matrix U, we find two uncoupled differential equations for each of the eigenmode displacements r 1 and r 2 : (S9)

C. Frequency response of the coupled masses
Solving equation S9 for the eigenmodes r 1 and r 2 gives The mechanical response f 1,2 of each eigenmode is modified as a function of the ratio ρ through U ij : For large detuning compared to the coupling strength (ρ 1), the mechanical oscillator strength F 1,2 is recovered for each mass: For a large positive detuning, the oscillator strengths shown in equation S12 are correct. For a large negative detuning, the oscillator strengths of the eigenmodes are swapped. To accomodate this problem, we need to change the transformation matrix U upon changing the sign of the detuning ∆k. Switching the rows in the transformation matrix U yields the correct oscillator strengths at large negative detuning. The responses at large detuning are then: The full description of the transformation matrix U is 4 then: (S14)

D. Calculation of the dissipated power
The dissipated power P in the system is given by the product of the velocity of the eigenmodes and the right hand side of the equations of motion, L, in the eigenvector basis: (S15) The driving force oscillates as Re(e iωt ). Therefore, only the real part of the solution for r 1,2 is relevant. After taking the real part and calculating the derivative of r in time, we derive an analytical solution for the timedependent dissipated power, P . We are interested in the steady-state solution of P . After averaging the total power P over time, the total average power P is given by The total average power can be seen as the sum of the average power of each of the two eigenmodes: Assuming that κ k then k+κ m k m . Using this we can rewrite the term Ω 2 ± − ω 2 = (Ω ± − ω) (Ω ± + ω) (Ω ± − ω) 2ω. Inserting this approximation into equation S16 yields: (S18) The average power is the sum of two Lorentzians centred around each of the eigenvalues Ω ± . This is a significant result: the model predicts straightforward lineshapes.
E. Linking the mechanical system to the optical response of coupled excitons The experimental results shown in the main paper study excitonic interactions. Therefore, we forge an analogy between the mechanical and optical systems. We show that excitonic coupling can be readily described by the mechanical model assuming a classical description of the driven excitons.
We model the interaction of inter-and intralayer excitons as two coupled oscillating optical dipoles driven by an external light field. The optical dipole representing the intralayer exciton (index 1) has a constant energy of ω 0 . The IE (index 2) is modelled by a dipole that has a tunable energy through the term ∆ω with a total energy of ω 0 + ∆ω. The two dipoles are coupled to each other through a coupling constant κ. The detuning ∆ω is a linear function of the applied electric field F z on account of the Stark effect. From experiments we know that the energy crossing point occurs at F z = 0. Therefore, where µ is the static dipole moment and F z,0 is the field at which the two bare energies match.
From the experiment we know that the sign of µ is positive for IE 1 and negative for IE 2 . The two dipoles are driven by the external light field E = E 0 Re(e iωt ). The force of the field on the charge e is e E. Each dipole has a different coupling strength F 1,2 to the light field and a dephasing constant γ 1,2 .
The resulting equations of motion have the same form as the equations of motion of the mechanical system shown in equation S2. Therefore, the calculations and solutions of the mechanical system also hold for the excitonic system. The eigenenergies of the system are given by assuming ω 0 + κ ω 0 , as ω 0 κ. For the eigenvectors, the ratio ρ is in this case defined as ρ = ∆ω κ . The excitonic eigenmodes can be written as The result for the eigenmodes r 1,2 from equation S21 is related to the dipole moment via So far we have described the response of two individual dipoles interacting with each other. In practice, an array of dipoles exists in a semiconductor. In a completely classical approach, the quantum well is treated as a twodimensional array of optical dipoles [8]. Following this approach, the linear absorption α(ω) of a quantum well is given by [8] The definition of the radiative coupling constant Γ rad where f is the oscillator strength, n is the refractive index of the quantum well, c is the speed of light, and 1/A is the dipole area density. From equation S23 we see that the solution of the quantum well absorption is additionally broadened by Γ rad 0 due to dipole-dipole interactions in the quantum well [8,9].
The absorption of an array of dipoles shown in equation S23 has a similar form as the equation for the average dissipated power P of two coupled masses in the mechanical system given in equation S18. Both equations are equivalent if we take the experimental dephasing γ 1,2 to be equal to the dephasing of the dipole array γ + Γ rad 0 . By comparing the two equations we find a linear relation between the absorption Im(χ) of two coupled dipoles and the dissipated power P : The absorption strength I 1,2 of each coupled optical dipole is then dependent on the energy detuning ∆ω, the coupling κ, and the oscillator strengths F 1,2 . For all calculations in the main paper and in this supplement, the prefactor 4 e 2 A 0 cmn is set to 1 eV. I 1,2 is given by The absorption is then modelled as F. Semi-classical frequency response and absorption of excitons in a semiconductor In a semi-classical description, the coupling of a dipole to a light field is described by a quantity called the oscillator strength f . The oscillator strength is proportional to the expectation value of the dipole operator squared | φ n | p| φ 0 | 2 [10,11]. Usually f is a sum over all n resonances of the oscillator. In the following we will only consider the ground state transition from state |φ 0 to state |φ 1 with a resonance angular frequency ω 0 . The oscillator strength f is defined as [10,11] Then the semi-classical solution for the response of an oscillating dipole with a dephasing rate Γ can be written as [8,11] In the linear optical regime, the optical susceptibility χ(ω) relates the probe field E to the probe-induced macroscopic polarisation P (ω) [9]: The Coulomb interaction between charges (exciton effects and screening) needs to be considered for a quantum mechanical description of the optical response of semicondcutors. These effects are described by the socalled semiconducting Bloch equations (SBEs) [11,12]. The SBEs are a set of coupled equations that include the microscopic polarisation P k as well as electron and hole distributions [13]. The index k describes the momentum dependence of the charge carriers. The microscopic polarisation is related to the macroscopic polarisation as follows [14]: where the indices c and v stand for valence and conduction band, respectively.
For a classical light field, the optical response of the semiconductor can be calculated by linking the SBEs to Maxwell's wave equation (MSBE) [9]. One can calculate the macroscopic polarisation as [14] where Ψ α is the excitonic wavefunction, ω α is the excitonic resonance angular frequency, and γ α is the polarisation dephasing rate. We now consider the excitonic ground state Ψ 0 and rewrite equation S32 with the definition of the oscillator strength f in equation S28 to:  Remarkably, the result for the macroscopic polarisation following the SBEs is the same as the result for the polarisation of an optical dipole driven by an oscillating field, compare equation S33 with equation S21 and equation S29.
M. Kira and S. W. Koch calculate the linear absorption α(ω) of a quantum well following the MSBEs [9]. Their result for the absorption is the same as the absorption of an array of oscillating dipoles in a two-dimensional plane derived in a classical way shown in equation S23 [8].

G. Influence of the sign of the coupling strength on the excitonic absorption
As mentioned in the main text, comparing the calculated excitonic absorption according to equation S27 with the measured experimental absorption reveals the sign of the coupling strength. SI Fig. 3a shows the calculated average dissipated power of the IE-A and IE-B interactions with the same parameters as the calculations shown in Main Fig. 2 and SI Fig. 11. SI Fig. 3b shows the calculated dissipated power if the two coupling strengths have opposite sign. The two individual interactions are added together to highlight the IE "cage" (SI Fig. 3a) and IE "anti-cage" (SI Fig. 3b) behaviour. For the IE cage, the IE is weak "inside", i.e. at energies between the A-like and B-like excitons, and strong on the "outside"; the IE anti-cage describes the opposite behaviour. By comparing the two calculations with the measured absorption shown in Main Fig. 2 it is clear that a negative coupling models the IE-A interaction, and a positive coupling models the IE-B interaction. SI Fig. 4 shows the calculated absorption spectra for the IE-B and IE-A interaction corresponding to the measurement regions shown in Main Fig. 2a and Main Fig. 2c, respectively. The calculated maps reproduce the experimental data quite well.

H. Coupled RLC circuits
Applying the interaction model to the two experimentally observed excitonic interactions yields a positive coupling constant for the IE-B interaction and a negative coupling constant for the IE-A interaction. A negative coupling constant in a mechanical framework is of course counter-intuitive. Therefore, we will look at an alternative, an electrical system consisting of two coupled RLC circuits (see SI Fig. 5).
First, we derive the equation of motion of a single RLC circuit. The complex impedance of a RLC circuit is given

as
where C is the capacitance, R is the resistance, and L is the inductance.
The voltage can be related to the current I as The current is equal to the change of the charge q on the capacitor in time: We take the voltage driving the circuit as and the corresponding charge as We then calculate the equation of motion of the RLC circuit: can be directly compared to the equation of motion of a mechanical oscillator system: The inductance L corresponds to the mass m, the resistance R corresponds to the damping γ, and the inverse of the capacitance 1 C corresponds to the spring constant k. Alternatively we can write equation S39 as Next, we look at two RLC circuits coupled together by a shared impedance Z c (see SI Fig. 5). Both circuits are driven by an AC voltage with different magnitude 8 but with the same angular frequency and phase. We consider two distinct cases: Z c is either a capacitor or an inductor. Let us take Z c = 1 iωC c first. One can then write the equations of motion as: with C 1,tot = C 1 C c C 1 + C c and C 2,tot = C 2 C c C 2 + C c being the series capacitance of C 1,2 and C c . The matrix form of this coupling is: We define the sign of the electrical coupling in the same way as for the mechanical system shown in equation S2. For a capacitive coupling of the two RLC circuits, the coupling constant is positive. This means the IE-B interaction can be seen as a capacitive coupling.
Next, Z c is set to iωL c . The matrix form of the equations of motion can be written as: For an inductive coupling of the two circuits, the coupling constant is negative. The IE-A interaction can be described as an inductive coupling.
I. Constructive and destructive interference SI Fig. 6 shows the calculated absorption strengths I 1,2 according to equation S26 using the same parameters as extracted for the IE-A interaction in the main text. Here, the absorption strengths are plotted as a function of energetic detuning ∆ω as defined in equation S19. For large detuning, both absorption strengths approach the value of the respective oscillator strengths. This model also predicts destructive interference in eigenmode 1 (SI Fig. 6a), constructive interference in eigenmode 2 (SI Fig. 6b), for negative detuning and a negative coupling constant. The sign of κ determines which eigenmode experiences destructive interference. For κ > 0 (κ < 0) the higher (lower) energy mode shows destructive interference.
The maximum absorption strength in the eigenmode with constructive interference corresponds to the sum of the squared oscillator strengths F 1,2 . The minimum absorption strength in the eigenmode with destructive interference is zero (see SI Fig. 6).
Next, we look at the influence of the coupling strength and the oscillator strength ratio on the constructive/destructive interference. SI Fig. 7a shows the energy detuning ∆ω of the maximum/minimum intensity due to constructive/destructive interference in each eigenmode for various coupling strength values κ. The sign of κ is the same as the sign of the detuning ∆ω needed to reach the extrema. A larger absolute value for κ shifts the absorption strength extrema towards larger detuning. Additionally, the constructively/destructively interfered absorption strengths are spread over a larger detuning range for larger coupling strengths. Let's turn to the large positive coupling strength of the IE-B interaction with κ = 35.8 meV. The effect of the constructive/destructive interference would only be visible at a high detuning of roughly −175 meV = − 5 MV/cm (see SI Fig. 7a and SI Fig. 9 top row). Reaching such high electric fields is experimentally not feasible.
The dependence of the detuning position of the extreme points on κ is fitted with a linear function yielding a slope of ρ = 5.0 for F 1 /F 2 = 5.159. SI Fig. 7b shows the slope values as a function of the oscillator strength ratio F 1 /F 2 . For a large imbalance between F 1 and F 2 , the slopes increase. This means that even for a small κ the detuning needed to reach the extrema is quite large (compare vertical dashed green lines in SI Fig. 9 top row). Towards F 1 /F 2 = 1 the slopes approach zero. A smaller imbalance of F 1 and F 2 pushes the extrema towards zero detuning (compare vertical dashed green lines in SI Fig. 9 top and bottom row). From an experimental point of view, the interplay between coupling strength and oscillator strengths is important for the constructive/destructive interferences to be visible for reasonable experimental parameters (compare SI Fig. 9).

IV. QUANTUM-CONFINED STARK EFFECT OF THE A-EXCITON
In an external electric field perpendicular to the MoS 2 layers F z , excitons experience an energy shift ∆E given by where µ z is the excitonic dipole moment and β z is the excitonic polarisability. This effect is called the quantum-confined Stark effect (QCSE) [15]. While the interlayer excitons have a significant out-of-plane dipole moment, the intralayer excitons have a near-zero out-of-plane dipole moment (µ z 0) leading to a quadratic QSCE [16].
The quadratic shift of the A-exciton is included in the model when fitting the IE-A interaction in the main text. This is carried out by adding the term −β z F 2 z to the uncoupled energy of the optical dipole 1 representing the A-exciton. Then, the eigenenergies have an additional term depending on the quadratic energy shift of the optical dipole 1 (A-exciton). The eigenenergies are The matrix elements of U that are not equal to one become The fits of the IE-A interaction shown in Main Fig. 2  and in SI Fig. 11 use the equations modified by the quadradic Stark shift of the A-exciton.
The fit yields a polarisability of the A-exciton of β z = 6.4 × 10 −9 D m V −1 . The A-exciton polarisability in the bilayer is roughly a factor of 10 larger than the polarisability of A-excitons in a monolayer MoS 2 [16]. In homobilayer MoS 2 the valence band of both layers are mixed [17]. This leads to an increased effective quantum well width as compared to the monolayer which in turn also leads to an increased polarisability [18].
The B-exciton can't be described nicely with a simple peaks. The raw data of those three measurements are shown in SI Fig. 10. Each spectrum is fitted by a sum of three Gaussians. First, the removal procedure of the B L1 in SI Fig. 10a and SI Fig. 10b is described. The linewidth and amplitude of the substracted B L1 exciton at all electric fields is set to the fitted values at the highest electric field. The energy of the B L1 exciton for the scan in SI Fig. 10a (SI Fig. 10b)    removed with these fitting parameters.
Second, the substraction of the A L2 exciton in SI Fig. 10b and SI Fig. 10c   field. The energy is taken to be the fitted zero-field energy modified by the quadratic Stark shift for finite fields. The polarisability is set to the value extracted from the coupled dipole model shown in SI Table I.
The exciton energies and peak areas shown in Main Fig. 2 are extracted from the unmodified measured absorption spectra shown in SI Fig. 10a and SI Fig. 10c. The fitting routines yielding the energies and peak areas are described. The IE-A interaction shown in SI Fig. 10c is fitted as the sum of three Gaussians where the two Gaussians describing the A-excitons in each layer are set to the same energy and intensity (see SI Fig. 10e). This is a reasonable assumption if the spectra are energtically far away from the energy crossing point compared to the coupling strength. For the IE-B measurement shown in SI Fig. 10a, spectra at electric fields higher than −0.81 MV/cm are fitted by a sum of three independent Gaussians while for lower electric fields the two B-excitons in each layer are set to the same energy and intensity for reasonable fitting results (see SI Fig. 10d).
The absorption maps shown in Main Fig. 2a and Main Fig. 2c are overlaid by the eigenenergy evolution as determined by the coupled dipole model. SI Fig. 11 shows the extracted peak energies from the absorption spectra shown in SI Fig. 10a and SI Fig. 10c. The coloured dashed lines in SI Fig. 11a are fits to the IE-and B-energies according to equation S20. The IE-and Aenergies shown in SI Fig. 11b are fitted by equation S47.

A. Computational details
The atomic structures, the quasi-particle band structures, and the optical spectra were obtained using the VASP package [19,20]. Core electrons were treated by the plane-augmented wave scheme [21,22]. A lattice parameter value of 3.22Å was set for all calculation runs. For all the calculation cells, a grid of 15×15×1 k-points was used together with a vacuum height of 21.9Å. Van der Waals interactions between the layers were included by perfoming the geometry's optimisation process at the PBE-D3 level [23]. All atoms were allowed to relax with a force convergence criterion below 0.005 eV/Å. A Heyd-Scuseria-Ernzerhof (HSE) hybrid functional [24][25][26] was used as an approximation of the exchange-correlation electronic term, including spin-orbit coupling (SOC). It was used to determine the eigenvalues and wave functions as an input for the full-frequencydependent GW calculations [27] performed at the G 0 W 0 level. The electric field was applied at this step, just before GW calculation process. Its application was only a small perturbation to the band structures, considering small/moderate electric field values. For partial occupancies, an energy cutoff of 400 eV and a gaussian smearing of 0.05 eV width were chosen. A tight electronic minimisation tolerance of 10 −8 eV was set to determine the corresponding derivative of the orbitals with respect to k needed in quasi-particle band structure calculations with a good precision. The direct band gap convergence was carefully checked (smaller than 0.1 eV as a function of k-points sampling). Afterwards, the total number of states included in the GW procedure was set to 1280 together with an energy cutoff of 100 eV for the response function. A Wannier interpolation proce- dure performed by the WANNIER90 program [28] was used to obtain the band structures. All optical excitonic transitions were calculated by solving the Bethe-Salpeter Equation [29,30]. The twelve highest valence bands and the sixteen lowest conduction bands were used to obtain eigenvalues and oscillator strengths on all systems. From these calculations, we use the imaginary part of the complex dielectric function to construct the absorbance with an artificial broadening of 12 meV.

B. Oscillator strength evolution and interband transition composition of the true excitonic eigenstates
Excitons in TMD monolayers do not consist of a single excitonic transition [31]. In fact, the true excitonic eigenstates are mixed states of multiple excitonic transitions. Let's start by looking at the composition  [31]. This mixing is due to a strong intravalley Coulomb exchange interaction and is also present in the bilayer case, as shown in Main Fig. 4. In the following, we will further discuss the GW+BSE results on the MoS 2 bilayer.
1. IE2-BL2 interaction SI Fig. 12 summarises the results from the GW+BSE calculations for eigenmode 1 and eigenmode 2 of the IE 2 -B L2 interaction. As already discussed in the main text, these two excitons interact strongly with each other. The evolution of the calculated relative absorption strength I 1 and I 2 as a function of electric field is shown as the cyan points in SI Fig. 12b and SI Fig. 12c, respectively. Upon tuning through the interaction region, a clear transfer of absorption strength between the two eigenmodes is visible. These beyond standard DFT results agree quite well with the measured absorption shown in Main Fig. 2a  SI Fig. 12b and SI Fig. 12c also show the interband transition composition of the true excitonic eigenstates at each electric field step. The composition percentage indicates how much each interband transition contributes to the overall excitonic absorption strength. The involved interband transitions are colour-coded and shown in SI Fig. 12a. The calculations reproduce nicely the layer character of each eigenmode upon tuning through the interaction region. Initially, eigenmode 1 (eigenmode 2) has a large interlayer IE 2 , yellow (intralayer B L2 , blue) component. Upon increasing the electric field the IE 2 and B L2 percentages in each eigenmode change: one increases at the cost of the other one. Finally, at high electric fields eigenmode 1 (eigenmode 2) has changed its layer character to an intralayer (interlayer) character. The swapping of the layer character of each eigenmode confirms the interaction of IE 2 and B L2 [3].

BL1 exciton
To confirm that IE 2 only interacts with B L2 , we will now look at the calculation results for the B L1 exciton. SI Fig. 13 summarises the calculated relative absorption strength of B L1 and its interband transition composition. The B L1 absorption strength stays constant over the whole electric field range. The interband transition composition shows several interesting features. First, the intralayer B L1 (blue) component dominates the oscillator strength and stays rather constant for all studied electric fields. These results are a confirmation that the B L1 exciton does not interact with the IE 2 exciton. Second, at 4 MV/cm the B-interlayer exciton (BIE, pink) component is increased as compared to lower and higher electric fields. An explanation could be that B L1 and BIE with opposite spin interact in a similar way as IE 1 and A L1 . From measurements we know that the BIE is energetically located just below the B-excitons at around 2.07 meV [3]. An energetic crossing of B L1 and BIE could therefore be expected around this electric field range. Experimentally, we can't resolve this interaction due to the very small oscillator strength of the BIE and the probably small coupling strength to B L1 . Last, the true B L1 exciton retains a very large and rather constant IE 1 (yellow) component over the whole electric field range even though IE 1 is energetically detuned far away from the B-exciton. This indicates that the orbital hybridisation and consequently the large oscillator strength of the IE is not affected much in the studied electric field range. The IE oscillator strength is not diminished even after energetically tuning it below the A-excitons.

AL2 exciton
Next, we will discuss the GW+BSE results for the true A L2 exciton (see SI Fig. 14). Similarly to the B L1 exciton discussed before, the absorption strength of A L2 and its composition stay constant over the studied electric field range. This confirms that the A L2 exciton does not interact with IE 1 .

Applying coupled dipole model to DFT results
The calculated GW+BSE excitonic peak energies and relative absorption strengths can also be fitted with the model presented in Section III, similarly as the measured data shown in Main Fig. 2. The energies are fitted with equation S20 and the relative absorption strengths are fitted with equation S26. The energy fit of the IE-B interaction shown in SI Fig. 15a yields a coupling strength of κ = +14.2 meV, a B-exciton energy ω 0 of 2.11 eV, a crossing point F z,0 of 4.74 MV/cm, and a linear Stark shift µ of 0.21 e nm. The oscillator strength fit (SI Fig 15b) yields an oscillator strength ratio of F 1 /F 2 = 6.26 and confirms the positive sign of the IE-B coupling.
The energy fit of the IE-A interaction shown in SI Fig. 15c yields an absolute coupling strength of  In addition to the measurement results, the coupled dipole model is also able to parameterise the results from the post-DFT calculations quite well.