Unconventional excitonic states with phonon sidebands in layered silicon diphosphide

Complex correlated states emerging from many-body interactions between quasiparticles (electrons, excitons and phonons) are at the core of condensed matter physics and material science. In low-dimensional materials, quantum confinement affects the electronic, and subsequently, optical properties for these correlated states. Here, by combining photoluminescence, optical reflection measurements and ab initio theoretical calculations, we demonstrate an unconventional excitonic state and its bound phonon sideband in layered silicon diphosphide (SiP2), where the bound electron–hole pair is composed of electrons confined within one-dimensional phosphorus–phosphorus chains and holes extended in two-dimensional SiP2 layers. The excitonic state and emergent phonon sideband show linear dichroism and large energy redshifts with increasing temperature. Our ab initio many-body calculations confirm that the observed phonon sideband results from the correlated interaction between excitons and optical phonons. With these results, we propose layered SiP2 as a platform for the study of excitonic physics and many-particle effects.

n exciton, the electron-hole pair formed via Coulomb interaction, is an ideal platform for understanding many-body effects [1][2][3][4][5][6][7][8] . The properties of excitons strongly depend on the crystal structure and dimensionality of host materials 9,10 . Due to quantum confinement, the electronic properties of quasiparticles (electrons, holes and excitons) in low-dimensional materials can be remarkably different from those in three-dimensional (3D) bulk materials. The Coulomb screening in low-dimensional quantum-confined structures, particularly in one-dimensional (1D) electronic systems, is known to be weaker than that in bulk systems and consequently leads to larger exciton binding energy [11][12][13][14] and other emergent excitonic phenomena 10 . Experimental observations of anisotropic excitons have been demonstrated in two-dimensional (2D) van der Waals (vdWs) materials in which electrons and holes taking part in the formation of 2D excitons are confined in the same monolayer [15][16][17][18][19] . Meanwhile, in 1D materials such as carbon nanotubes (CNTs) 12,20 , 1D excitonic states have also been observed in which constituent electrons and holes are known to be confined within 1D nanostructures. A unique excitonic state with hybrid dimensionality, which is yet elusive, such as a bound electron-hole pair with an electron confined along one dimension (1D-confined electron) and a hole confined along two dimensions (2D-confined hole), or vice versa, would be of great interest and highly desired in terms of its optical properties and interactions with other emergent quasiparticles.
In this work, we demonstrate the observation of an unconventional bright exciton in a layered silicon diphosphide (SiP 2 ) crystal, accompanied by a correlated phonon sideband in the optical spectrum. Based on our ab initio many-body GW and GW plus Bethe-Salpeter equation (GW-BSE) calculations, as well as non-perturbative model calculations, we find that the electrons constituting the excitons are confined within the 1D phosphorus-phosphorus chains of SiP 2, while the correlated holes extend over the 2D SiP 2 atomic plane. Therefore, excitonic states in layered SiP 2 are expected to exhibit hybrid dimensionality properties. Photoluminescence (PL) spectroscopy and reflectance contrast (RC) spectroscopy show that, regardless of the polarization of the excitation laser, the optical response of the excitonic state is always linearly polarized along the x direction of the SiP 2 lattice and is accompanied by a unique sideband feature. Both the excitonic emission and the sideband feature undergo dramatic redshifts as the temperature increases, in contrast to a slight temperature-dependent redshift of the band edge that is mainly influenced by electron-phonon coupling 21,22 . This reveals that in SiP 2 the interaction of the electronic degrees of freedom with the phononic degrees of freedom is strongly enhanced by excitonic effects. The phonon sideband feature can be theoretically modelled using a non-perturbative approach to describe the interaction between the unconventional excitons and optical phonon modes. Note that reduced dimensionality normally leads to excitonic features that are strongly affected by extrinsic Unconventional excitonic states with phonon sidebands in layered silicon diphosphide environmental effects, such as disorder from the substrate and surface additives 10 . Here we provide an investigation on the intrinsic excitonic behaviour in thicker, bulk-like SiP 2 flakes. Such a tightly bound unconventional exciton in SiP 2 not only can be envisioned as a platform for the exploration of exciton-phonon (ex-ph) coupling [23][24][25][26][27][28] and other many-body physics but also may lend itself to potential applications for anisotropic optoelectronic devices.

Crystal structure and electronic property of SiP 2
Layered SiP 2 is chosen as our target material because of its following unique characteristics. Compared with hexagonal layered materials such as graphene and MoS 2 , the cleavable SiP 2 crystal (space group Pnma) possesses an orthorhombic layered structure with a huge in-plane lattice anisotropy, as schematically shown in Fig. 1a and experimentally confirmed by scanning transmission electron microscopy-annular dark-field (STEM-ADF) imaging in Remarkably, based on their atomic surroundings, two types of inequivalent phosphorus atoms P A and P B can be distinguished in the SiP 2 lattice. As shown in Fig. 1a, P A binds to three silicon atoms, while P B binds to one silicon atom and the other two equivalent P B atoms. Note that the P B atoms along the y direction of the crystal lattice can naturally form phosphorus-phosphorus chains (denoted as P B -P B chains) embedded in the bulk SiP 2 (blue shades in Fig. 1a), which play a critical role in realizing the quasi-1D electronic states involved in exciton formation. To identify the variation in the chemical bonding environment around P A and P B atoms and the resulting unique properties of P B -P B chains in layered SiP 2 , we performed arsenic doping experiments (Supplementary Information, section 2.4) and used STEM characterization ( Supplementary Fig. 5). One can see that the doped arsenic atoms only selectively substitute the P B atoms inside the P B -P B chains (more details in Supplementary Fig. 5), indicating that the atomic structure containing P B -P B chains in SiP 2 is distinct from the buckled structure in black phosphorus.
More importantly, the anisotropy induced by quasi-1D P B -P B chains in layered SiP 2 directly results in unique electronic properties. Figure 1e shows the band structure of semiconducting bulk SiP 2 obtained from GW calculations. We found that the conduction band edge states in the X-Γ-Z plane of the first Brillouin zone (BZ) are relatively flat with a large effective mass (Supplementary  Table 1), and the corresponding charge densities are localized on the P B -P B chains (Fig. 1f), behaving like 1D-confined electrons. Importantly, in the direction along the P B -P B chains (the y direction of the crystal lattice), the electron hopping on P B atoms is significantly larger (bandwidth, ∼1.63 eV) than that across the P B -P B chains (bandwidth, ∼0.08 eV) (see details in Supplementary  Fig. 22), confirming the 1D nature of this electronic state on the conduction band edge. On the other hand, the hole states at the valence band edge do not show the same level of anisotropy (Supplementary Information, section 2.2), which, compared with 1D electrons, are relatively extended over the whole atomic plane in a quasi-2D fashion. The hybrid dimensionality of these band edge states in SiP 2 is remarkably different from those of the anisotropic 2D states in black phosphorus 15,29,30 . By analysing the calculated phonon bands given in Supplementary Information, section 13, we identify that the optical phonons localized on P B atoms and neighbouring silicon atoms could have a large coupling with quasi-1D electronic states in layered SiP 2 . exciton with 1D-confined electron and 2D-confined hole Figure 2a,b presents the PL spectrum and the second derivative of the RC (d-RC; see Supplementary Information, section 6) of an SiP 2 flake (228 nm) at 5.5 K, which reflects the light emission and absorption properties, respectively. The PL spectrum shows a main peak A at 2.06 eV (the lowest bright excitonic bound state denoted as the A exciton) and a broadened sideband feature Aʹ at 2.01 eV. The main peak A, obtained from all SiP 2 flakes measured at 5.5 K, is consistently located at an emission energy of 2.06 ± 0.01 eV (here, 0.01 eV is the energy uncertainty obtained from the standard deviation of the emission energies of several measured SiP 2 flakes; see Supplementary Information, section 7). Such a peak A in the PL spectrum matches the peak at 2.05 eV in the d-RC spectrum, as indicated by the red arrow (Fig. 2b). Due to the interference of the RC signals from the different interfaces in the SiP 2 thin films supported by substrates ( Supplementary Information, section 6), the phonon sideband feature is difficult to identify from the d-RC spectrum. Figure 2c shows the absorbance spectra obtained from the GW-BSE calculation and GW calculation with the random phase approximation (GW-RPA). Compared with the calculated absorption spectra based on GW-BSE and GW-RPA, we confirm that the emission peak A at 2.12 eV originates from the recombination of an excitonic state, in which the electronic states for electrons are quasi-1D and related electronic states for holes are quasi-2D (Fig. 1f). As shown in Fig. 2c,d, the calculated binding energy of such an unconventional exciton is approximately 140 meV (for more details about GW-BSE calculations, see Methods and Supplementary Information, section 12). From the modulus squared of the exciton wavefunction in real space shown in Fig. 3c, the observed exciton behaves like a Wannier-type exciton with twofold rotational symmetry, in sharp contrast to 2D excitons in monolayer transition-metal dichalcogenides 13 . More importantly, this exciton is embedded in a bulk layered material with an unusual atomic structure in contrast to those of reported pure 1D excitons in semiconductor nanowires 31 and CNTs 11,12 , leading to strongly anisotropic Coulomb screening for 1D-confined electrons and 2D-confined holes.

Anisotropic exciton and exciton-phonon coupling
Since the unconventional A exciton is mainly contributed by electrons and holes localized along the X-Γ-Z direction in the first BZ (Supplementary Information, section 12.2), we use the band edge states at the X point as the representative k-point to explore the influence of electron-phonon interactions on its electronic structures and optical response. Figure 2f shows the zero-point energy   individually induced by each optical phonon mode with momentum q = 0 at a temperature of zero. L 0 represents the lattice structure without displacement of phonon modes. ±L shift stands for atomic displacements of phonon modes. The change in the band gap is estimated by averaging the energy shifts of the band gap E X g with positive (red bar) and negative (dark blue bar) atomic displacements. Black circles denote energies of the corresponding phonon modes. g, The phonon density of states (DOS) for optical phonon modes, which is projected to the P B atoms in the embedded P B -P B chains (top) and their neighbouring silicon atoms (bottom). Insets: the P B atoms and their neighbouring silicon atoms.
shifts of the band gap at the X point induced by all optical phonon modes with momentum q = 0 at zero temperature. Here, we use the frozen-phonon approximation 32 to estimate the influence of optical phonon vibrations on the electronic states at the X point (see Supplementary Information, section 13.2 for more details). Since the electron wavefunctions of the A exciton are localized on the P B -P B chains, this unconventional exciton couples most strongly to optical phonons, whose vibrational modes are in the X-Γ-Z plane and involve P B atoms and neighbouring silicon atoms (Fig. 2g). These optical phonon modes dramatically modify the electronic structures of the quasi-1D states (Supplementary Information, section 13.4), indicating significant electron-phonon coupling within the P B -P B chains. Comparing the results in Fig. 2f,g, one can see that the prominent energy shifts are from the optical phonon modes with eigenenergies of ∼50-60 meV. More details are given in Supplementary Information, sections 13.2 and 13.4.
The experimental observation of the sideband feature Aʹ also indicates that ex-ph interaction on the quasi-1D P B -P B chains is at least moderately strong ( Supplementary Information, section 8). Therefore, we use a non-perturbative model to simulate the emergence of the sideband feature Aʹ, where a 'generalized Holstein Hamiltonian' is used with inputs from first-principles calculations, and the self-energy effects are included beyond the first-order Fan-Migdal diagram (Methods). In this model, we found that the fitted ex-ph coupling constant M of 30 meV is comparable to the relatively small bandwidth (or hopping, t ex = −20 meV) of the unconventional A exciton (for the estimate of t ex , see Methods and Supplementary  Information, section 13.3). Our approach is similar to the cumulant method considering the ex-ph coupling within the perturbative limit and makes use of the exponential assumption to include the self-energy effects from higher-order diagrammatic terms [33][34][35] . As shown in Fig. 2e, the appearance of the phonon sideband peak in the simulated spectrum agrees with the experimental results, indicating that sideband Aʹ originates from the ex-ph coupling between the unconventional exciton and the abovementioned optical phonon modes. Figure 3a,b shows the contour plots of the PL and d-RC intensity of bulk SiP 2 as a function of emission energy at different detection polarization angles θ (θ = 0° is set along the x direction), suggesting that the linear dichroic absorption and PL emission have similar twofold symmetry characteristics (see Supplementary Information, sections 4 and 12 for more details). Note that the observed linearly polarized PL emission remains along the x direction regardless of the incident laser polarization direction or the sample temperature, as shown in Supplementary Figs. 10 and 11. Our GW-BSE calculations ( Fig. 3d and Supplementary Fig. 21b) show that the absorption peak of the quasi-1D A exciton appears only when the polarization is along the x direction. The absorption signal inside the band gap along the y direction is forbidden by the SiP 2 crystal symmetry, which results in relevant optical excitonic matrix elements being zero ( Supplementary Information, section 12.3). We also performed pump-probe transient optical measurements to characterize the dynamics of the observed bright exciton in bulk SiP 2 (see Supplementary Information, section 9 for more details). The lifetime for the exciton in SiP 2 is as short as 250 fs, which is probably related to an ultrafast process that dissociates these linearly polarized bound excitonic states into unbound and unpolarized states. We further compare the energy shift of the A exciton peak with the energy shift of the quasiparticle band edge as the temperature changes. Figure 4a,b shows the temperature-dependent PL and d-RC spectra for bulk SiP 2 . As shown in Fig. 4c, the optical absorption of the band edges and the exciton peak A, and the sideband feature Aʹ, all exhibit clear redshifts with increasing temperature. The redshifts of the band edge can be fitted with the Bose-Einstein model (see Supplementary Information, section 7 for more details), suggesting that the interaction between electrons and phonons plays an important role in the energy shifts. The redshift of the band edge absorption resulting from the electron-phonon coupling 21,22,36 is approximately 20 meV at 300 K. On the other hand, the redshift of both peaks A and Aʹ is approximately 90 meV at 300 K, much larger than the energy shift of the direct band edge, indicating an additional contribution from the large coupling between the bound exciton and optical phonons. Such a result is consistent with the analysis of temperature-dependent linewidth broadening of the peak for unconventional A exciton (see Supplementary Information, section 5 for more details).

Outlook
Using optical spectroscopic measurements with the support of ab initio many-body calculations, we demonstrated the observation of an unconventional bright exciton in layered SiP 2 . In contrast to those reported 1D and 2D excitons truly confined in CNTs and monolayer transition-metal dichalcogenides, the bound excitonic states in layered SiP 2 exhibit hybrid low dimensionality due to the intrinsic 1D and 2D nature of the constituent electrons and holes, respectively. Interestingly, we envision that SiP 2 can host peculiar trion states, including a negatively charged trion (composed of two 1D-confined electrons and one 2D-confined hole) and a positively charged trion (composed of one 1D-confined electron and two 2D-confined holes). Once we couple layered SiP 2 to other vdWs semiconductors, such as monolayer MoS 2 , to form heterostructures, the interfacial layer coupling can change the rotational symmetry  Fig. 4 | Temperature-dependent spectra and energy evolution for exciton peak A and side peak A′ in SiP 2 . a,b, Temperature-dependent PL (a) and d-RC (b) spectra (second derivative of RC) from 5.5 K up to 300 K. The thick solid lines represent the spectra results, while the thin dashed and dotted lines represent fitting results of the PL A and Aʹ peaks. The solid triangle arrows highlight the redshifts of peak A (red), Aʹ (black) and the band edge (yellow) with increasing temperature. c, Temperature-dependent band edge (yellow circles) and A exciton (red circles) energies extracted from d-RC spectra, and A exciton (red spheres) and Aʹ (grey spheres) energies extracted by multipeak fitting of PL spectra. The error bars indicates the full width at half maximum.
of the semiconducting layers and could bring optical and optoelectronic functionalities via symmetry engineering at the heterointerfaces. Through the doping modulation of carrier polarity in SiP 2 or its heterostructures, rich excitonic physics with exotic dynamic behaviour can be realized in this material platform, such as interlayer excitons and Moiré excitons with tunable dimensionality. Furthermore, the interaction between this unconventional bound exciton and the optical phonon leads to an accompanying phonon sideband. Since a phonon and an exciton fall within the same energy range from zero to several hundred meV, we speculate that such many-body interactions may even lead to the emergence of elementary excitations beyond the Born-Oppenheimer limit in atomic 2D thin films or nanostructures of SiP 2 . Our work will provide a platform to further understand ex-ph coupling and other essential many-body physics and inspire follow-up studies and calculation method developments therein.

Online content
Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/ s41563-022-01285-3.

Methods
SiP 2 crystal growth with flux method. Single-crystalline samples were synthesized by using the tin flux method 37 . Silicon, phosphorus, gadolinium and tin were mixed at a Si:P:Gd:Sn ratio of 1:6:0.03:5 and sealed into evacuated quartz tubes. The mixture was slowly heated to 1,100 °C to avoid bumping phosphorous and kept for 48 h. Subsequently, the sample was cooled to 400 °C in 140 h and then cooled to room temperature by switching off the electric furnace. The tin flux was removed by using diluted HCl (aq). The obtained black crystals were then ultrasonicated in distilled water and ethanol to remove the residuals (such as phosphorous, adhered) on the crystal surface. This procedure was repeated until the water (or ethanol) became transparent enough after ultrasonication.

Sample preparation for optical measurements and STEM-ADF measurements.
SiP 2 flakes with thicknesses of 5-200 nm were prepared by mechanical exfoliation onto SiO 2 /Si wafers (300-nm-thick SiO 2 layer) or fused silica substrate. The thickness was identified by atomic force microscopy (integrated with a WITec Alpha 300 Raman system) after all optical measurements were finished. SiP 2 is stable in a nitrogen atmosphere and in vacuum and can gradually degrade when exposed to air within several hours. To avoid sample degradation, the whole sample preparation was processed in a glovebox. Atomic-resolution STEM-ADF imaging was performed on an aberration-corrected ARM200F equipped with a cold field-emission gun operating at 80 kV. The STEM-ADF images were collected using a half-angle range from ∼81 to 280 mrad. The convergence semi-angle of the probe was ∼30 mrad.
Optical measurements. Optical measurements, including the PL spectra and RC and Raman spectra, were performed using a confocal Raman system (WITec Alpha 300). Thickness-dependent PL measurements (Supplementary Information, section 3) were carried out at room temperature using a ×100 objective lens with an incident laser (laser power, 0.2 mW) focused to an ∼1 μm spot. Nitrogen conditions were accomplished by protecting samples using continuous nitrogen gas flow. Low-temperature PL and Raman measurements were performed under vacuum conditions with samples installed in a cryostat (Cryo Instrument of America RC102-CFM microscopy cryostat) using a long working distance ×50 objective lens (laser power, 3 mW). For RC measurement, we recorded the subtracted reflectance of the sample normalized by the reflectance of the substrate, that is, RC = 1 − Rsample Rsub , where R sample represents the reflectance of the SiP 2 sample on the silicon dioxide or quartz substrate and R sub represents the reflectance of the bare substrate.
First-principles calculations. First-principles density functional theory (DFT) calculations were performed by using the projector-augmented wave (PAW) 38,39 method implemented in the Vienna Ab initio Simulation Package (VASP) 40 . The energy cut-off for the plane wave basis is set to 500 eV. To test the lattice constants to compare with the experimental value, we used the exchange-correlation functionals of generalized gradient approximation (GGA) with Perdew-Burke-Ernzerhof (PBE) type, local density approximation (LDA), and the PBE functional with vdWs corrections to fully relax the lattice structures. The vdWs interactions were included by using the methods proposed by Dion et al. 41 with the optB88-vdW functional. We found that lattice constants obtained from the method including the vdWs corrections are closest to the experimental values (Supplementary Information, sections 10 and 11), which are used in the following phonon bands, GW and GW-BSE calculations. During the lattice relaxations, the force convergence criterion was 10 -3 eV Å −1 , and a 9 × 21 × 7 k-point mesh was sampled over the BZ. For the self-consistent electronic structure calculations, we set the energy convergence criterion to 10 -6 eV and the k-point mesh to 11 × 25 × 9 over the whole BZ. The phonon spectrum was calculated by the PHONOPY package 42 in the framework of density functional perturbation theory with the finite-displacement approach, in which a 2 × 4 × 1 supercell was employed.
Using VASP 43 , GW calculations 44 were performed using Kohn-Sham DFT wavefunctions (GGA-PBE) calculated on a 4 × 16 × 4 k-point mesh as the initial mean field. The dielectric response function used for the fully frequency-dependent eigenvalue-self-consistent GW calculation is summed over 1,240 Kohn-Sham states (corresponding to a 100 eV cut-off). The frequency grid is divided into a dense part ranging from 0 to 13.75 eV and a coarse grid tail ranging from 13.75 to 178.78 eV. The grid sampling is non-uniform with 80 frequency grid points, resulting in step sizes ranging from 0.31 eV in the dense grid up to 46.35 eV in the tail 45 . We use multiple iterations in the GW calculation to update the eigenvalues of the Kohn-Sham states when calculating both Green's function G and the screened interaction W while keeping the initial Kohn-Sham wavefunctions unchanged. Full convergence was reached after five iterations. This procedure results in better agreement with the experimental results because the standard G 0 W 0 approach underestimates the band gap by approximately 220 meV. The maximally localized Wannier functions obtained from the Wannier90 packages 46,47 were used to plot the GW quasiparticle band structure (Fig. 1e). In the construction of the Wannier functions, the s and p orbitals of both the silicon and phosphorus atoms were used as initial trial wavefunctions. The GW quasiparticle energies and Kohn-Sham wavefunctions are used to construct the kernel of the BSE 48,49 . We employed the standard Tamm-Dancoff approximation and included ten conduction bands and ten valence bands during the calculation of the GW-BSE Hamiltonian.
Calculation of the spectral function of the phonon sidebands. To model the spectral function of the phonon sidebands, we solved for the dressed polaron Green's function of the generalized Holstein Hamiltonian 50,51 , H = H 0 + V, where H0 = ∑ q ωqb † q bq + ∑ k ϵ k c † k c k is the unperturbed single-particle Hamiltonian is the interaction Hamiltonian. For the first term constituting H 0 , q is the crystal momentum of the phonon, b † q and b q are the phonon creation and annihilation operators, and ω q (which shall be taken as a constant independent of q) is the average phonon energy of the dominant optical branch responsible for ex-ph coupling. For the second term constituting H 0 , k is the centre-of-mass momentum of the exciton, c † k and c k are the exciton creation and annihilation operators, and ∈ k is its energy dispersion. We obtain these values from our GW calculations. The interaction Hamiltonian, V, represents the ex-ph interaction, with M k,q being the ex-ph coupling matrix element. Since there is only one exciton in the model Hamiltonian, its solution is independent of the statistics of particle 1 . The same solution would be obtained for any fermion or boson, such as electrons (which is more common), as long as the particles are free to move. To construct the generalized Holstein Hamiltonian, we first calculate most of its parameters from first-principles calculations and finally fit the ex-ph coupling matrix elements, M k,q (also to be taken as a constant), to the experimental results.
First, we note that although more than one phonon mode contributes to the renormalization of the band gap, the dominant contributing phonon modes fall within the energy range of 50-60 meV ( Fig. 2f and Supplementary Fig. 24). Since the phonon bands that project strongly onto the P B -P B chain are relatively flat within the X-Γ-Z plane in reciprocal space ( Supplementary Fig. 23b), we assume the representative phonon band to have negligible phonon bandwidth, as in the Einstein model. Hence, we use the energy of a representative longitudinal optical (LO) phonon mode, ω q ≡ ω LO = 55 meV, to model the ex-ph interaction, as obtained from the ab initio phonon calculations. We also assumed that the exciton has a free-particle dispersion of a periodic 1D chain, namely, ϵ k = −2tex cos(ka), where 4t ex = −80 meV is the exciton bandwidth. The exciton hopping term, t ex , was estimated using the hole bandwidth (4t h ≈ −80 meV) and the electron bandwidth (4t h ≈ 640 meV), which are calculated along the X-Γ-Z direction of the Brillouin zone from our GW calculations ( Supplementary Fig. 22), with the formula t −1 ex = t −1 h + t −1 e . Finally, using the above parameters obtained from the first-principles calculations, we only fit the ex-ph parameter, M k,q ≡ M = 30 meV, so that the spectral function of the calculated dressed Green's function reproduces the PL spectrum shown in the optical experiments (Fig. 2e).
In the calculation of the dressed interacting polaron Green's function, G (k,ω), Dyson's identity, [G (k, ω)] −1 = [G0 (k, ω)] −1 − Σ(k, ω) is used, where G 0 (k,ω) is the free-particle Green's function, given by G0 (k, ω) = (ω − ϵ k + iη) −1 and Σ(k,ω) is the ex-ph self-energy, which consists of an infinite sum of all proper self-energy diagrams. Written more explicitly, G(k,ω) can be written as a continued fraction, , such that Σ(k,ω) is the second term in the denominator given by , that when expanded in powers of M 2 reproduces the Feynman diagrams of each order [52][53][54] . In the calculation of the self-energy, we used the momentum-averaged non-interacting Green's function, as introduced by Berciu 53 and extended by Goodvin, Berciu and Sawatzky 54 . In this approximation, the momentum-dependent non-interacting Green's function, G 0 (k,ω), in the expression of the self-energy, was replaced by its momentum average, Ḡ 0 (ω), given by where N k is the number of k-points and ρ 0 (ϵ) is the density of states. The momentum-averaged self-energy, Σ MA (ω), is now momentum independent, and the interacting Green's function is now G (k, ω) = 1 G −1 0 (k,ω)−ΣMA (ω) . Finally, the spectral function was given by the imaginary part of the interacting Green's function, the main peak of which is fitted to the excitation energy of excitonic state A as obtained from GW-BSE calculations.

Data availability
Source data are provided with this paper. The authors declare that data generated or analysed during this study are provided as source data or included in the