Superradiant coupling effects in transition-metal dichalcogenides

Cooperative effects allow for fascinating characteristics in light-matter interacting systems. Here, we study naturally occurring superradiant coupling in a class of quasi-two-dimensional, layered semiconductor systems. We perform optical absorption experiments of the lowest exciton for transition-metal dichalcogenides with different numbers of atomic layers. We examine two representative materials, MoSe$_2$ and WSe$_2$, using incoherent broadband white light. The measured transmission at the A exciton resonance does not saturate for optically thick samples consisting of hundreds of atomic layers, and the transmission varies nonmonotonously with the layer number. A self-consistent microscopic calculation reproduces the experimental observations, clearly identifying superradiant coupling effects as the origin of this unexpected behavior.


INTRODUCTION
With the ability to fabricate them as monolayers, the electronic and optical properties of naturally occurring transition-metal dichalcogenide (TMDC) multilayer structures have attracted renewed interest. While excitons in monolayers are naturally confined within the layer, systematic studies of the band structure as a function of the number of layers show that even though the system becomes energetically indirect, the direct gap at the K -points of the Brillouin zone is preserved even in the bulk limit [1]. The near-K -point bands have a flat dispersion along the K − H axis, indicating that even in a multilayer structure, the corresponding quasi-particles are strongly confined within the individual layers and can be considered as effectively two-dimensional. Indeed, recent theoretical and experimental studies attribute the excitonic series in TMDC multilayer systems as being due to quasi twodimensional intra-and interlayer excitons [2,3].
Up to now, most theoretical and experimental investigations of TMDC structures have focused on the evolution of the bandgap and exciton binding, while changing the dielectric substrate or number of layers. Surprisingly, little or no attention has been paid to changes in the optical spectra of TMDC multilayer structures due to light propagation effects and radiative coupling. Therefore, the studies in this paper will focus on those specific points.
In 1954, Dicke [4] first considered the subtle situation in which two emitters decay in close proximity to each other and affect each other's emission and absorption processes, leading to cooperative behavior. In general, radiative coupling occurs if an ensemble of oscillators interacts with a common light field. The interaction of a single oscillator with the light field (re-)emitted by all oscillators in the ensemble leads to the formation of collective modes, which describe the response of the system as a whole instead of the response of the individual emitters. In fact, the emission dynamics of a single two-level system are altered by the presence of a second one, even if it is in its ground state.
As already shown by Hopfield [5], radiative coupling in a resonantly excited homogeneous semiconductor leads to the formation of exciton-polaritons with infinite radiative lifetimes. The infinite radiative lifetime is an immediate consequence of the infinite system dimensions, inhibiting an emitted photon from escaping from the sample. Instead, it is trapped in an endless cycle of (re-)absorption and (re-)emission processes, forming the hybrid exciton-light state, i.e., the exciton-polariton. Nevertheless, the polariton is damped due to nonradiative dephasing processes in its material part that lead to an irreversible (true) absorption of energy from the optical field. The concomitantly reduced light transmission can be described by the Beer-Lambert law, predicting an exponential decay of the transmitted light as a function of the optical thickness.
In contrast, a photon emitted by an individual localized oscillator escapes from the sample without re-absorption, thus opening a radiative decay-channel. The oscillator strength, f λ , determines both, the reversible light absorption and the emission, and for an oscillator with transition energy E λ , the associated radiative decay rate is Γ λ f λ E λ ∕ℏ [6][7][8][9]. In an ensemble of such localized oscillators, cooperative effects depend strongly on the relative phases of the light (re-)emitted by different oscillators. If N identical excited oscillators are confined to a region small compared to the optical wavelength, the overall emission is in phase, leading to a peak intensity of the emitted light field ∝N 2 . This superlinear intensity increase is accompanied by a corresponding reduction of the radiative lifetime such that the total emitted energy, R I tdt, increases linearly with the number of oscillators, in agreement with the requirement of energy conservation. This phenomenon is denoted in the literature as "superfluorescence" or "superradiance" [4]. For the ensemble of identical localized oscillators, the superradiant mode has the features of a single oscillator mode with N times enhanced coupling strength, i.e., with N times reduced radiative lifetime, and is the only optically active mode, whereas all other N − 1 modes are dark.
In solid state physics, superradiant coupling effects have been studied intensively in semiconductor multiple quantum well (MQW) structures [6][7][8][9][10]. In an MQW, excitons are confined to the quantum-well region, and the lateral spacing between the wells can be controlled within nanometer accuracy. The analogous situation to what was originally considered by Dicke is then realized in a so-called Bragg geometry where the center-to-center distance of neighboring quantum wells (QW) corresponds to an integer multiple of half the optical wavelength, leading to phase coherence between the emitters in the different QWs. Indeed, four-wave mixing experiments on MQW Bragg structures have demonstrated the hallmark of superradiance, a radiative decay rate increasing linearly with the QW number [10].
In the present paper, we analyze the influence of radiative coupling and light propagation effects on the optical absorption of TMDCs as a function of the number of layers in the crystal. Similar to semiconductor QWs, absorption of light in the TMDCs leads to the formation of excitons that are confined reasonably well within the individual layers. However, as compared to typical QW excitons, the excitons in TMDC monolayers display huge oscillator strengths, leading to as much as ∼10% extinction of the incoming light for a typical monoatomic layer of e.g., WSe 2 [11]. As has been demonstrated recently, the transmittance of a single layer can be decreased even further if it is encapsulated, e.g., by hexagonal boron nitride [12,13]. Such an encapsulation protects the TMDC monolayer from surface damage and leads to a considerable reduction of the nonradiative excitonic linewidth, manifesting itself in an observed reflection coefficient as large as 80% for a monolayer MoS 2 . According to the Beer-Lambert law, this strong single-layer extinction should lead to a rapid decay of the transmitted light at the excitonic resonance frequency with an increasing number of atomic layers. In fact, a simple calculation of the optical density indicates that the crystal should be nearly opaque already for a few tens of atomic layers. However, the close proximity of different layers, together with the very large excitonic oscillator strength, provides favorable conditions for the occurrence of superradiant radiative coupling effects.
Whereas superradiant light emission cannot be expected in an optically indirect semiconductor, the absorption characteristics at the direct gaps at the K -points in a multilayer TMDC flake can still experience a superradiant coupling with a correspondingly reduced radiative lifetime. In an absorption measurement, the increased radiative decay rate then leads to a reduced percentage of the incoming energy that is irreversibly absorbed by the sample, i.e., to a reduction of the true absorption. To investigate this behavior, we perform a series of conceptually simple broadband white light transmission and reflection experiments to study the cooperative radiative coupling effects in TMDCs as a function of crystal thickness.

EXPERIMENT
Using the experimental setup shown schematically in Fig. 1, we carefully measured the optical transmission and reflection of MoSe 2 and WSe 2 samples of different thicknesses at 5 K . We have collected atomic force microscopy (AFM) images of all the samples studied. The thicknesses of the samples and the resulting numbers of atomic layers were obtained using AFM for each sample. In addition, for thinner crystals, Raman scattering and photoluminescence spectroscopy were used. Transmission and reflection spectra were measured quantitatively with respect to the total incoming light, and all the possible losses were carefully taken into account. The true absorption was obtained as where I 0 is the intensity of the incoming signal, I T is the intensity of the transmitted signal, and I R is the intensity of the reflected signal.
The experimental transmitted, reflected, and absorption spectra are shown in Supplement 1, together with the corresponding theoretically obtained spectra. The deduced optical density data at the peak of the A exciton resonance are summarized in Fig. 2 as a function of the number of layers. We notice remarkable deviations of the measured absorption from the Beer-Lambert law, which predicts that at the A exciton resonance frequency, the sample should be essentially opaque for thicker crystals. Instead, we observe an optical density well below saturation. This phenomenon occurs in both studied systems, MoSe 2 and WSe 2 , indicating that it is not unique to a specific material, but likely characteristic for the entire class of TMDC materials. In order to measure the lateral extent of the excitonic wave function, we measured the diamagnetic frequency shift of the A exciton under extremely high magnetic fields. For this purpose, we utilized the 65 T pulsed magnet at the National High Magnetic Field Laboratory at Los Alamos National Laboratory. Magneto-reflectance studies were performed with the samples at cryogenic temperatures down to 4 K using a home-built fiber-coupled optical probe. Broadband white light from a xenon lamp was coupled to the samples using a 100 μm diameter multimode optical fiber. The light was focused onto the sample at nearnormal incidence using a single aspheric lens, and the reflected light was refocused by the same lens into a 600 μm diameter collection fiber. The collected light was dispersed in a 300 mm spectrometer and was detected with a liquid nitrogen cooled charge-coupled device (CCD) detector [14]. The use of very large magnetic fields allows for the observation of the small quadratic diamagnetic shift of the excitons in the TMD materials. From the diamagnetic shift, we deduce the lateral extent of the excitonic wave function to be 3 to 5 times larger than the interlayer distance. We further estimate the excitonic wave function to be well contained within the individual layers, leading to negligible surface effects for thicker samples, supported by measurements and calculations on MoS 2 , WS 2 , and MoTe 2 [2,15]. Theoretical calculations on bulk MoTe 2 estimate >90% of the excitonic wave function to be contained within the individual layers [2].

MICROSCOPIC THEORY
To analyze the observed phenomena, we employ a semiclassical theory of radiatively coupled two-dimensional excitons based on the semiconductor Maxwell-Bloch equations for a classical optical field interacting with equidistantly spaced two-dimensional layers. Treating the Hamiltonian of the individual layers within an effective four-band model, screening and radiative coupling of the bands under consideration is included dynamically, whereas contributions of all the other bands and the dielectric environment are contained in a background dielectric tensor [16]. For the dielectric displacement field of the layered material, we make the ansatz where ϵ k and ϵ ⊥ represent the in-plane and out-of-plane background dielectric constants, respectively, and P describes the nonlocal, time-and frequency-dependent resonant contributions of the bands under consideration. (We use c. g. s. units throughout the theory.) Using the generalized Coulomb gauge, ϵ k ∇ k · A k ϵ ⊥ ∂ z A z 0, and restricting the analysis to light propagation along the c axis, a division into transverse and longitudinal contributions yields the wave equation for the transverse optical field and Poisson's equation for the scalar potential: where E T − 1 c _ A; E L −∇ϕ; A and ϕ are the vector and scalar potentials; and the superscripts T and L denote transverse and longitudinal parts, respectively.
On the length scale of the optical wavelength, the resonant quasi-two-dimensional contributions of the individual layers can be approximated by where z n n − 1∕2d is the central position of the nth layer, d is the natural layer-to-layer distance, and P nk is the optically induced polarization in the nth layer. The radiative coupling effects result from the interaction of the individual layers with the light emitted by all layers in the sample [9]. The total optical field at the position of the nth layer can then be expressed as where E in k z, ω represents the incoming field with a dispersion determined by the dielectric environment of the nonresonant background and substrate, and G 0 z, z 0 , ω describes the radiation of a localized oscillator embedded within this environment. As the material polarization is induced by the total transverse Research Article field, the second part of Eq. (4) results in a dynamical dipoledipole interaction mediated by the optical field.
In the linear regime, the transverse part of the optically induced layer polarization is related to the transverse optical field via P T nk ω χ T n ωE T k z n , ω, where χ T n ω is the transverse susceptibility of the nth layer. To compute this susceptibility, we employ the widely used massive Dirac-fermion (MDF) model [17], where the Hamiltonian describes an effective four-band system adapted to the symmetry properties of the hexagonal lattice. In reciprocal space, the first Brillouin zone has the same honeycomb geometry as in real space, and a direct gap occurs at the corners of the Brillouin zone, i.e., the K Dirac points. The nonequivalent Dirac points are identified by the so-called valley index, τ 1, and they are related by the parity transformation. The near-K -point quasi-particles are described as relativistic Dirac fermions with pseudospins that couple to the light field. In a naturally A-A 0 -stacked multilayer structure, adjacent layers are related by the parity transformation. Thus, the joint Brillouin zone of the multilayer structure can be considered to be composed of mirror-identical copies of the 2D monolayer Brillouin zone, leading to a layer-dependent valley index, τ n −1 n .
For the individual layers, we compute the layer susceptibility microscopically from the coupled gap equations and the Dirac-Bloch equations (DBE), as described in detail in Ref. [16]. Here, we have a single free parameter, which can be related to the physical thickness of a single monolayer. All other parameters are obtained from density functional theory (DFT) calculations. Since the coupled equations contain the screened Coulomb interaction, which is obtained as solution of Poisson's Eq. (3) for the electron and hole charge density, ρ, they contain a Coulomb-mediated interlayer interaction, allowing for a seamless description of the transition from monolayer to bulk. For the numerical evaluations, we use the material parameters listed in Ref. [16].
As screening by the surrounding layers influences the spectral position of the exciton and its oscillator strength, the different dielectric environments of the individual layers lead to slightly different resonance frequencies. Therefore, the response of a multilayer structure is generally inhomogeneously broadened due to the individual contributions from the different layers. However, this effect becomes negligibly small for sufficiently thick structures where the distance of most layers from the sample surfaces well exceeds the in-plane exciton Bohr radius. In such cases, the linear susceptibilities of the individual layers within the sample become identical, influenced only by intrinsic screening effects.
Alltogether, the transmission and absorption coefficients of a multilayer sample reflect the interplay between the individual layer susceptibilities, the dielectric environment, and the propagation effects. For a monolayer embedded between two dielectrics with refractive indices n T The transmitted and reflected intensities are related to the amplitude coefficients via I T n B jtE 0 j 2 , I R n T jrE 0 j 2 , and the incoming intensity is given by I 0 n T jE 0 j 2 . One thus finds for the true absorption coefficient of a monolayer, αω 4πn T ω c Imχ T ω jn T n B ∕2 − 2πi ω c χ T ωj 2 : For a single 2D oscillator with resonance frequency ω 0 and nonradiative decay rate γ, one finds a true absorption coefficient: αω n T n γΓ∕n∕ω − ω 0 2 γ Γ∕n 2 , where Γ is the radiative decay rate for the monolayer in vacuum, andn n T n B ∕2 is the average refractive index of the top and bottom materials. As a true absorption process not only requires an initial excitation but also a subsequent scattering into an optically dark state, the peak absorption is determined by the interplay between radiative and nonradiative decay rates and reaches its maximum value if both are equal. For a typical nonencapsulated monolayer on top of a fused silica substrate, this gives an extinction coefficient of 10% to 20% of the incoming intensity, depending on the specific material and the nonradiative homogeneous linewidth. Neglecting the changes in the dielectric environment and radiative coupling, the Beer-Lambert law would predict an exponential decay of the transmitted intensity, T N n B n T jtj 2N , leading to an extinction coefficient of roughly 90% for a sample consisting of only five layers, in strong disagreement with the experimentally observed opacity.
An increase in the number of layers modifies the transmission and absorption characteristics in two ways: first, it changes Research Article the effective background refractive index experienced by the individual layers, and, second, it leads to radiative coupling. The effect of changing the dielectric environment alone is best illustrated by comparing the transmitted intensity of a monolayer on a substrate with that of a monolayer with the same linear susceptibility embedded in bulk. Such a comparison is shown in Fig. 3(a) for the example of MoSe 2 , using the bulk layer susceptibility. For the background in-plane dielectric constant of MoSe 2 we use ϵ k 5.01, which has been extracted by combining bulk DFT calculations, with the analytically computed contributions, with the longitudinal in-plane dielectric constant [16]. As can be recognized, the increased background refractive index in the bulk environment reduces the layer extinction by roughly a factor of 2.5. The predictions for the transmitted intensity based on Beer's law, once using the monolayer-on-substrate transmission coefficient and once using the monolayer-in-bulk transmission coefficient, are shown as dark yellow dashed and brown dashed-dotted curves in Fig. 2, respectively. As can be recognized, the increased background refractive index of the bulk environment decreases the opacity of a multilayer sample predicted by Beer's law considerably, yielding saturation only after roughly 100 layers.
Apart from simple dielectric effects, the optical spectra are modified by radiative coupling features, which depend sensitively on the phase coherence between the emitters in different layers. Moderately increasing the number of layers and ensuring that the total sample thickness, N × d , remains much smaller than the optical wavelength in the material, all emitters/absorbers are in phase and the emitted fields from the different layers interfere constructively. Here, the superradiant mode is the only bright one, and the system effectively behaves as a single emitter with an N times enhanced coupling strength.
Since the radiative lifetime is inversely proportional to the coupling strength, the enhanced coupling in the superradiant geometry reduces the radiative lifetime and increases the radiative broadening. While the effect of the modified dielectric environment is clearly observable in the transmission characteristics, the radiative broadening is best observed in the true absorption spectra. In Fig. 3(b), the dashed lines show the absorption per layer for a sample containing 10, 20, and 30 layers in a superradiant geometry on a BK7 glass substrate. As we can see, the superradiant coupling reduces the absorption per layer with an increasing number of layers, with a correspondingly increased linewidth (FWHM). The full lines in Fig. 3(b) show the theoretically predicted absorption/layer computed for the corresponding samples on a BK7 substrate using the full geometry. The strong similarity of the spectra shows that, within this thickness range, the radiative coupling is essentially superradiant.
If the total sample thickness is further increased such that it exceeds the optical wavelength, the system essentially develops into an effectively homogeneous medium with a polariton dispersion

THEORY-EXPERIMENT COMPARISON
In order to test our model assumptions, we probe the lateral extent of the excitonic wave function by measuring the diamagnetic frequency shift of the A exciton under very high magnetic fields. The energy position of the A exciton resonance as a function of energy for magnetic fields up to 60 T perpendicular to the layers or parallel to the c axis is shown in Figs. 4(a) and 4(c) for WSe 2 and MoSe 2 , respectively. From the energy position of the A exciton at 60 T and −60 T, we obtain the diamagnetic shift shown in Figs. 4(b) and 4(d) for WSe 2 and MoSe 2 , respectively. The diamagnetic shift is fitted using a simple quadratic formula, which yields shifts of the average position of ∼0.77 μeV∕T 2 and ∼0.12 μeV∕T 2 for WSe 2 and MoSe 2 , respectively. The lateral exciton size is obtained from the diamagnetic shift via r X ffiffiffiffiffiffiffi ffi hri 2 p ffiffiffiffiffiffiffiffiffiffi 8m r σ p ∕e, where m r is the exciton reduced mass. Using m r ∕m 0 0.17 for WSe 2 and m r ∕m 0 0.27 for MoSe 2 , we obtain r X 2.44 nm and r X 1.21 nm for WSe 2 and MoSe 2 , respectively. Using the natural layer-tolayer distance of d 6.5 Å for both materials, the exciton radius is roughly 2 to 4 times larger than the interlayer distance, resulting in negligible surface effects for samples with N ≳ 7.
To analyze the experimental results, we use the values ϵ k 5.01 and ϵ k 13.6 for the bulk in-plane dielectric constants for MoSe 2 and WSe 2 , respectively. These are determined by subtracting the resonant contributions to the dielectric constant from those obtained from DFT calculations for the respective bulk parent materials [18] and from the measured values at 1024 nm [19], where the resonant contributions have been computed analytically within the MDF model [16]. To avoid rapid oscillations from the substrate/vacuum interface that cannot be resolved experimentally, we perform simulations for a half-space geometry, i.e., air/sample/substrate. Furthermore, we use a phenomenological nonradiative linewidth of 15 meV for the A-and 25 meV for the B-excitonic series for all samples. Though the assumption of equal nonradiative decay rates in different samples is questionable in principle, this is a reasonable assumption for the larger samples where surface effects become negligible, and it allows us to identify changes in the absorption characteristics due to radiative coupling effects only. The resulting calculated peak absorption and transmission as a function of the number of layers (sample thickness) is plotted in Fig. 2 for MoSe 2 (Top) and WSe 2 (Bottom) as green and blue-dashed lines, respectively. In particular, the true peak absorption (green lines) exhibits a behavior that differs significantly from the calculated absorption following the Beer-Lambert law, be it based on the monolayer-on-substrate transmission (dark yellow dashed line) or based on the monolayer-in-bulk transmission (brown dashed-dotted line). As we can see, the experimental data are reasonably well reproduced by the calculated true absorption. Furthermore, we experimentally observe a weak, nearly constant absorption of the multilayer structure over the range N 15 to N ∼ 50, in stark contradiction to the exponential decay of the transmitted intensity prescribed by the Beer-Lambert law.
Under superradiant coupling conditions, the peak absorption at resonance is proportional to γN Γ∕γ N Γ 2 , where γ and Γ are the nonradiative and radiative decay rates, respectively. In Fig. 2, the superradiant peak absorption is plotted as a red dashed line. For MoSe 2 , the numerical evaluations give a radiative decay rate of ℏΓ 0.6 meV. Estimating the nonradiative linewidth Research Article to be 15 meV, the maximum peak absorption is expected for N ≈ 20, which decays slowly in the region between N 20 and N 50. These features are characteristic for superradiant coupling and in good agreement with the experimental observations. For samples containing more than roughly 60 layers, we observe deviations from the superradiant coupling. Here, i.e., at a total sample thickness of roughly a quarter-wavelength, the system transitions from the superradiant to the bulk regime.
Due to the larger background in-plane dielectric constant, the optical wavelength in the WSe 2 sample at the resonance frequency is only roughly 120 nm. Correspondingly, the transition from the superradiant regime to the bulk limit occurs for smaller sample sizes with N ≈ 45, which is again in good agreement with the experimental observations.
In Fig. 5, we plot the full absorption spectra measured experimentally and calculated theoretically for both MoSe 2 and WSe 2  for different numbers of layers. We obtain quantitative agreement between experiment and theory in terms of the absolute absorption strength. The experimental line shapes are well reproduced by the theoretical calculations that include superradiant coupling, which further reaffirms the validity of the theoretical approach and the superradiant nature of the light propagation through layered TMDCs.

SUMMARY AND CONCLUSIONS
In conclusion, we analyze the optical transmission/absorption properties of exfoliated TMDC samples for different numbers of layers, ranging all the way from monolayers to bulk-like configurations with more than a hundred layers. We carefully characterize the respective sample thicknesses using AFM; for thinner samples, photoluminescence and Raman spectroscopy were additionally used. The absorption was quantitatively measured for the A exciton resonance by carefully measuring transmission and reflection simultaneously, and by taking into account all the possible losses. The resulting optical density as a function of sample thickness deviates significantly from the Beer-Lambert law. We model the observed effects employing the semiconductor Maxwell-Bloch equations for a classical optical field interacting with equidistantly spaced two-dimensional layers, where the A excitons are assumed to be well localized within the individual layers. To validate the model assumptions, we carefully measure the excitonic Bohr radii in the two representative bulk TMDs using the diamagnetic energy shift of the excitonic resonance at magnetic fields up to 65 T. The exciton wave functions obtained from these measurements are well localized within the layers and are also in agreement with recent angle-resolved photoemission studies performed in bulk WSe 2 , which reveal a two-dimensional character of the bands at the K -point [20]. The theoretical calculations reveal that the experimentally observed variation of the resonant optical absorption can be uniquely attributed to the superradiant coupling between the excitons in the different TMDC layers.