Electrically tunable dipolar interactions between layer-hybridized excitons

Transition-metal dichalcogenide bilayers exhibit a rich exciton landscape including layer-hybridized excitons, i.e. excitons which are of partly intra- and interlayer nature. In this work, we study hybrid exciton–exciton interactions in naturally stacked WSe2 homobilayers. In these materials, the exciton landscape is electrically tunable such that the low-energy states can be rendered more or less interlayer-like depending on the strength of the external electric field. Based on a microscopic and material-specific many-particle theory, we reveal two intriguing interaction regimes: a low-dipole regime at small electric fields and a high-dipole regime at larger fields, involving interactions between hybrid excitons with a substantially different intra- and interlayer composition in the two regimes. While the low-dipole regime is characterized by weak inter-excitonic interactions between intralayer-like excitons, the high-dipole regime involves mostly interlayer-like excitons which display a strong dipole–dipole repulsion and give rise to large spectral blue-shifts and a highly anomalous diffusion. Overall, our microscopic study sheds light on the remarkable electrical tunability of hybrid exciton–exciton interactions in atomically thin semiconductors and can guide future experimental studies in this growing field of research.

where φ ξL n,k is the excitonic wave function in state n = 1s, 2s..., valley ξ = (ξ e , ξ h ), and layer L = (l e , l h ) and E ξL n,bind is the exciton binding energy. Here, the reduced exciton mass m ξL red = m ξe le m ξ h l h m ξ h l h +m ξe le , as well as the screened electron-hole Coulomb interaction V c le v l h q enter. The valley-specific electron (hole) masses m ξele (m ξ h l h ) are obtained from DFT calculations [2]. When evaluating the Coulomb matrix elements we explicitly include the finite thickness of the TMD layers as well as the as the dielectric environment via a generalized Keldysh screening [3]. In this work, we explicitly include hybridization of intra (X)-and interlayer (IX) excitons. In particular, the four possible intra-and interlayer exciton states (here expressed as L ≡ IX 1 , IX 2 , X 1 , X 2 focusing on the energetically lowest n = 1s transitions such that the exciton index can be omitted) are generally coupled by electron/hole tunneling. The resulting hybrid exciton states are obtained from diagonalizing the following exciton Hamiltonian [4,5] H x,0 = ξ,L,Q E ξ L,Q X †ξ L,Q X ξ L,Q + ξ,L,L ′ ,Q containing the exciton centre-of-mass dispersion E ξ L,Q = ℏ 2 Q 2 2M ξL + E ξL bind + ∆ ξL , M ξL = m ξ h l h + m ξele being the total exciton mass, X ( †) being excitonic and bosonic ladder operators and ∆ ξL is the valley-specific band gap. The free Hamiltonian also contains a tunneling contribution which takes into account the tunneling of electrons and holes between different layers (l e ̸ = l ′ e or l h ̸ = l ′ h ) via the excitonic tunneling matrix element The excitonic tunneling matrix element crucially depends on electron and hole tunneling strengths, T cξe lel ′ e and T vξ h l h l ′ h respectively, as well as exciton wave function overlaps F ξ LL ′ = k φ * ξL k φ ξL ′ k . The electron and hole tunneling strengths are obtained from ab-initio calculations and are reported in Ref. [4] for common TMD bilayers. For the considered case of 2H-stacked WSe 2 homobilayers we adopt the tunneling strengths T cK = 0, T vK = 66.9 meV and T cΛ = 236.6 meV for the most relevant K/K' and Λ/Λ' valleys in this structure. Note that the electronic tunneling matrix elements are generally stacking-and momentum-dependent, however in this work we focus on naturally stacked (H h h ) homobilayers, and evaluate the matrix elements at the high-symmetry points. The Hamiltonian in Eq. (S2) is now diagonalized via the basis transformation [5] X †ξ where Y ( †) is a new set of hybrid exciton operators and C ξη L (Q) is the mixing coefficient determining the relative intra/interlayer content of the hybrid exciton, enabling us to define a hybrid exciton state as |hX η ⟩ = i=1,2 (C η Xi |X i ⟩+ Electronic Supplementary Material (ESI) for Nanoscale. This journal is © The Royal Society of Chemistry 2023 2 C η IXi |IX i ⟩) with i=1,2 (|C η Xi | 2 + |C η IXi | 2 ) = 1 for a fixed hybrid exciton state η. The mixing coefficients are obtained from solving the following hybrid eigenvalue problem [6] introducing the hybrid exciton eigenenergy E (hX)ξ η,Q . We can now express the exciton Hamiltonian above in the hybrid basis such that H x,0 →H x,0 withH By solving the eigenvalue problem in Eq. (S5) we get microscopic access to the full hybrid exciton landscape in TMD bilayers. Furthermore, we investigate the impact of an electric field on the hybrid exciton landscape. This is done by exploiting the electrostatic Stark shift of interlayer resonances, i.e. by taking E ξ L=IX,Q → E ξ L=IX,Q + ∆E, with ∆E = ±de 0 E z , where d ≈ 0.65 nm is the dipole length (here assumed to be the same as the TMD layer thickness [7]),e 0 is the electric charge and E z is the out-of-plane electric field [8]. In this way, hybrid exciton eigenenergies and mixing coefficients become tunable with respect to electric fields. In Table S1, we report the exciton hybrid energies and intralayer and interlayer mixing coefficients for hybrid excitons composed by electrons in the ξ e = K, K', Λ, Λ' valleys and holes in the ξ h = K, K' valleys in 2H-stacked hBN-encapsulated WSe 2 homobilayers. The energies are given relative to the intralayer A exciton energy and the electric fields E z = 0, ±0.3 V/nm are considered.

Exciton
Energy Note that the KΛ and K'Λ' exciton states are energetically degenerate at vanishing electric fields. This is a consequence of the H-type stacking, where the individual TMD layers are rotated 180 degrees with respect to each other such that the spin-orbit coupling in one of the layers is inverted. Moreover, these states can be expressed as |KΛ⟩ = |IX 2 ⟩ such that each of the states only mixes contributions from a single intralayer and a single interlayer exciton species. Hence, it follows that the KΛ and K'Λ' hX carry opposite out-of-plane dipole moments via their interlayer components, and therefore the energy of these states shifts in opposite directions under the application of an electric field (cf . Table S1).

II. HYBRID EXCITON-EXCITON INTERACTION HAMILTONIAN
Here, we provide a microscopic derivation of the hybrid exciton-exciton interaction Hamiltonian. The starting-point is the bilayer carrier-carrier Hamiltonian: where λ ( ′ ) = (c, v), ξ, and l ( ′ ) are the band, valley, and layer indices, respectively. Here, the operators λ ( †) annihilate (create) carriers in band λ. Moreover, we note that V λ l λ l ′ q describes an intraband intralayer Coulomb interaction if l = l ′ and an interlayer Coulomb interaction if l ̸ = l ′ . Furthermore, we consider the long-range part of the Coulomb interaction such that V 2ϵ0A|q|ϵinter,q (l ̸ =l). The intra-and interlayer dielectric functions ϵ intra,q and ϵ inter,q can be found in the Supplementary Material of Ref. [9]. Interband Coulomb interactions, which give rise to electron-hole exchange [10] or Auger recombination [11], are not expected to contribute significantly to experimentally accessible density-dependent energy renormalizations (Supplementary Section IV) and are therefore neglected in this work.
Given the carrier-carrier Hamiltonian, we now proceed as follows: i) we find the equation of motion for the transform the equation of motion to the exciton basis [12], iii) make an ansatz for the exciton-exciton interaction Hamiltonian and compute the equation of motion for the polarisation in the exciton picture, iv) read off the exciton-exciton interaction matrix element such that the results from steps ii) and iii) coincide. Starting with the first step i), we obtain the equation of motion for the polarisation directly from the Heisenberg equation of motion [1]. Including only the Coulomb contributions from Eq. (S7) we obtain Next, we transform the entire equation above to the excitonic basis and make use of the pair operator expansions [12] c where the pair operators can be further expressed in the exciton basis as P ξele,ξ h l h k,k ′ = n φ ξL n,β ξL k+α ξL k ′ X ξ n,L,k−k ′ , with φ ξL n,k being the exciton wave function (cf. Supplementary Section I) and the compound indices ξ = (ξ e , ξ h ), L = (l e , l h ) (such that l e = l h corresponds to the intralayer wave function and l e ̸ = l h corresponds to the interlayer wave function). In the following, we will only consider the lowest-lying n = 1s exciton states, so that the index n can be omitted. By doing this, the equation of motion Eq. (S8) separates into two parts, a direct part and an exchange part. The second, fourth, fifth and seventh term in Eq. (S8) gives rise to the direct terms reading Here, we also defined the excitonic form factors F (x ξL q) ≡ k φ * ξL k+x ξL q φ ξL k . We may now construct the corresponding direct exciton-exciton interaction Hamiltonian with with the direct part of the exciton-exciton interaction reading such that a commutation of the excitonic Hamiltonian (S11) with the polarisation gives rise to Eq. (S10). We note that, in the long wavelength limit i.e. we find a vanishing direct interaction between intralayer excitons (L, L ′ = X i , i = 1, 2) and recover the widely used plate capacitor formula when considering interactions between interlayer excitons (L, L ′ = IX i , i = 1, 2) as has been previously confirmed in literature [11,13,14]. Here, the material-specific constants d i,TMD and ϵ ⊥ i,TMD denote individual TMD monolayer thicknesses and out-of-plane components of the TMD dielectric tensors, respectively. In the main manuscript we set d 1,TMD = d 2,TMD ≡ d TMD and ϵ ⊥ 1,TMD = ϵ ⊥ 2,TMD ≡ ϵ ⊥ TMD as we are considering a homobilayer. Note that interactions between different interlayer exciton species IX i and IX j , i ̸ = j are attractive due to the opposite dipole orientations of these excitons. Now, we consider the remaining terms (i.e. the first, third, fifth and eight terms) in Eq. (S8) and find that these give rise to the following exchange terms from which we may construct an exchange matrix element such that where the exchange part of the exciton-exciton interaction reads . (S16) Here, we note that the first term corresponds to hole-hole exchange within the excitons and the second term corresponds to electron-electron exchange. In particular, the Kronecker deltas imply that fermionic exchange of individual charge constituents is only allowed if charges of the same species reside in the same layer and valley. Moreover, the exchange interaction is generally dependent on both centre-of-mass momenta Q 1 , Q 2 as well as the relative momentum q. In the long wavelength limit (q, Q 1 , Q 2 ≪ a −1 B , a B being the exciton Bohr radius) the exchange interaction is nonvanishing for both intra-and interlayer exciton species, and it is the dominating contribution to the exciton-exciton interaction for intralayer excitons [13,15,16]. However, we remark that the resulting density-dependent energy renormalizations due to intralayer exchange interactions are negligible (see Supplementary Section IV). Exchange interactions are therefore not considered in the main manuscript, but included here only for the sake of completeness. Hence, we obtain the exchange part of the exciton-exciton Hamiltonian We now have access to the multilayer exciton-exciton interaction involving both intra-and interlayer excitons. However, generally, excitons are hybridized between the layers due to electron/hole tunneling. To include the effect of hybridization, we transform the excitonic operators to the hybrid basis, cf. Eq. (S4). The exciton-exciton Hamiltonian then transforms into a hybrid Hamiltoniañ where the hybrid exciton-exciton interaction contains of a direct part and an exchange part according tõ with the unhybrised direct (D) and exchange (E) matrix elements defined in Eq. (S12) and (S16), respectively. Having derived the most general form of the hybrid exciton-exciton interaction, we now remark on the considered case of untwisted homobilayers. In this case, it holds that the mixing coefficients are approximately constant in momentum, such that C ξη L,Q ≈ C ξη L [4]. Consequently, the direct hybrid exciton-exciton interaction depends only on the relative momentum q. In the main manuscript, we only consider the lowest-lying hybrid exciton states for each valley configuration and therefore the indices η i , i = 1...4 are omitted therein. Moreover, note that the intra-and interlayer mixing coefficients enter the hybrid exciton-exciton interaction strengths. This provides an intriguing way of tuning the interaction strength with externally applied electric fields.

III. DIPOLE-DIPOLE INTERACTION
In here, we show that the real-space representation of the direct interlayer exciton-exciton interaction (Eq. (S12), with L = L ′ = IX and ξ = ξ ′ ) can be interpreted as a classical dipole-dipole interaction at large distances. By considering the two TMD layers forming the homobilayer as two infinitely thin slabs separated by the distance d and approximating the dielectric environment as homogenous, with a single effective dielectric constant ϵ BL , we find an analytic expression for the direct exciton-exciton interaction. Within these approximations, the intra (X)-and interlayer (IX) Coulomb interactions read [3] Now, substituting these simplified expressions into the direct exciton-exciton interaction and setting the excitonic form factors F ≈ 1, we find that the direct interlayer exciton-exciton interaction reduces to ) .
The real-space representation of the interaction is obtained by taking the Fourier-transform: ) .
Finally, we are interested in the asymptotic behavior of the interaction and therefore let r ≫ d. In this limit, we find which is precisely a classical dipole-dipole interaction. Now, as the interlayer components of the mixing coefficients are approximately constant in momentum (cf. Supplementary Section I) for untwisted TMD homobilayers, we note that also the corresponding real-space hybrid exciton-exciton interaction obeys the asymptotic 1/r 3 behavior for large r.
In Fig. S1, we show the direct hybrid exciton-exciton interaction strength between KΛ excitons in hBN-encapsulated WSe 2 homobilayers as a function of distance, r, including the dominant interlayer contributions to the interaction (solid yellow curve). Unlike in the main manuscript, we here plot the logarithm of the interaction strength over large distances. Importantly, we find that the interaction scales as 1/r 3 (dashed black curve) at distances r ≳ 15 nm. This confirms the dipole-dipole-like character of the hybrid exciton-exciton interaction at large distances.
FIG. S1. Real-space representation of the hybrid exciton-exciton interaction for KΛ excitons. Note that we here have taken the logarithm of the interaction strength.

IV. DENSITY-DEPENDENT ENERGY RENORMALIZATIONS FOR HYBRID EXCITONS
Having access to the microscopic hybrid exciton-exciton Hamiltonian (Eq. (S18) and Supplementary Section II) implies that we are now also able to investigate density-dependent energy renormalizations of hybrid excitons. These energy renormalizations can be derived from the corresponding Heisenberg equation of motion for the (hybrid) polarisation ⟨Y † ζ,Q ⟩ reading where we introduced the compound index ζ = (ξ, η) including the valley index ξ and the hybrid exciton index η. Now, we consider the equation above on a Hartree-Fock level, i.e. we expand the appearing bosonic three-operator expectation value into single-particle expectation values and neglect two-particle correlations. We also make use of the random phase approximation (RPA) [1] such that the equation above becomes where we defined the hybrid exciton occupation N ζ Q ≡ ⟨Y † ζ,Q Y ζ,Q ⟩. Here we also approximated ⟨Y † ζ,Q Y ζ ′ ,Q ′ ⟩ ≈ δ ζ,ζ ′ Q,Q ′ N ζ Q making use of the RPA. Note that the energy renormalization consists of two terms, the first term being a direct term and the second being an exchange term, with the latter reflecting exciton exchange [13]. In contrast to the exchange interactionẼ which includes exchange of individual carriers, the exciton exchange takes into account the exchange of individual excitons. Furthermore, we consider low temperatures in this work such that the exciton distribution N q is strongly peaked around q = 0 and assume the centre-of-mass momentum |Q| ≪ a −1 B , where a B is the exciton Bohr radius. This reduces the equation above to d dt ⟨Y † ζ,Q ⟩ ≈ i ℏ δE ζ ⟨Y † ζ,Q ⟩ introducing the energy renormalization where n ζ1 x ≡ 1 A Q N ζ1 Q with A being the crystal area (cancelling out with the area A in the electronic Coulomb matrix elements). Finally, we restrict ourselves to the energetically lowest hybrid exciton states in this work such that the compound index ζ reduces to the valley index ξ. Generally, the energy renormalization of a hybrid exciton ξ is obtained by taking into account the interactions between all the different exciton species. By assuming that n ξ x ∝ n x , i.e. assuming a thermalized Boltzmann distribution for the hybrid excitons, where n x = ξ n ξ x is the total exciton density, the energy renormalization of a single exciton species can be quantified by a valley-specific effective dipole length d ξ obtained from In the evaluation of valley-specific dipole lengths we here only include the direct (dipole-dipole) contributions to the hybrid exciton-exciton interaction, cf. Eq. (S18). This is done as interlayer exchange interactions are seen to provide a small quantitative correction to the dipole-dipole interaction [9]. Moreover, although intralayer exchange interactions taking into account individual exchange of carriers are dominant in TMD monolayers [11,14], they have negligible impact on the energy renormalizations, as their contributions are largely cancelled out against contributions due to higher-order correlations such as biexcitons [10]. This goes beyond the scope of the Hartree-Fock theory presented in this work. In Fig. S2 we illustrate the valley-specific effective dipole length as obtained from Eq. (S28) for the case of naturally stacked hBN-encapsulated WSe 2 homobilayers as a function of electric field, E z at low temperatures, T = 10 K. Intriguingly, we find a drastic increase in the dipole length at around E z ± 0.15 V/nm. This reflects the transition from a mostly intralayer KΛ (K'Λ') state to a mostly interlayer KΛ' (K'Λ) state under the application of a positive (negative) electric field. At vanishing electric fields, KΛ and K'Λ' excitons coexist. These excitons independently interact via weak repulsive dipole-dipole interactions and mutually interact with each other via attractive dipolar interactions as they exhibit opposite dipole orientations, giving rise to negligible effective dipole lengths. At the largest considered electric fields, E z = 0.3 V/nm, only KΛ' excitons are relevant and the impact of other excitons is negligible due to the large energy separations between exciton states (cf. Table S1). These excitons are mostly interlayer-like in nature (|C IX | 2 = 0.8, cf. Table S1) and interact strongly via dipole-dipole repulsion. Note that the large effective dipole moment of KΛ' excitons directly translates into large energy renormalizations (Eq. (S28)) and give rise to sizable blue-shifts of phonon sidebands with exciton density, as schematically illustrated by the inset in Fig. S2.