The quadruplon in a monolayer semiconductor

Abstract

The ultimate goal of understanding the structure of matter has spurred a constant search for composite particles, especially high-order correlated entities for nearly all forms of matter, from elementary particles, nuclei, and cold atoms, to condensed matter. So far, composite particles involving two or three constituent particles and their weak-coupling combinations have been experimentally studied, such as the Cooper pairs, excitons, trions, and bi-excitons in condensed matter physics, or diquarks, mesons, and di-mesons in quantum chromodynamics. Although genuine four-particle correlated entities have long been theorized in various materials, alternatively known as quadruplons (Rausch and Potthoff in New J. Phys. 18, 2016), quadrons (Quang et al. in Physica B 602, 2021), or quartets (Jiang et al. in Phys. Rev. B 95, 2017), the only closely related experimental evidence is the tetraquark observation at CERN (LHCb in  Nat. Phys. 18, 751–754, 2022). In this article, we present for the first time the experimental evidence for the existence of a four-body entity in condensed matter, the quadruplon, involving two electrons and two holes in a monolayer of Molybdenum Ditelluride. Using the optical pump–probe technique, we discovered a series of new spectral features in addition to those of excitons and trions. Furthermore, we found that all these spectral features could be reproduced theoretically using transitions between the two-body and four-body complexes based on the Bethe–Salpeter equation. Interestingly, we found that the fourth-order irreducible cluster is necessary and sufficient for the new spectral features by using the corresponding cluster expansion technique. Thus, our experimental results combined with theoretical explanation provide strong evidence for the existence of a genuine four-particle entity, the quadruplon. In contrast to a bi-exciton which consists of two weakly interacting excitons, a quadruplon involves tightly bound four-particle entity without the presence of well-defined excitons. Our results could impact the understanding of the structure of materials in a wide range of physical systems and potentially lead to new photonic applications based on quadruplons.

1 Main

Whereas the search for “elementary” constituent particles of matter has been a never-ending pursuit in high-energy physics, the understanding of the interactions or correlations among these particles dominates studies in lower energy scales that are more relevant to our daily experience and technology. Such interactions lead to the formation of correlated entities or composite particles that determine the basic material properties, underline our fundamental understanding of almost all fields of physical and material sciences, and provide the foundation for all modern technologies. In condensed matter physics, correlated entities, such as excitonic complexes: excitons (X), trions (T), and bi-excitons (BX), are critical to our understanding of the basic material properties, and especially to the ever-richer physics of the celebrated Mott transition beyond the simple exciton-plasma picture. Recently, multi-particle entities have attracted much interest, including the Bose–Einstein condensation (BEC) [5], BEC-BCS crossover via different bosonizations of Fermions in strongly correlated systems, or dropletons [6,7,8] in semiconductors. In superconductivity theories, the charge-4e configuration [3, 9] was proposed as an alternative to the conventional Cooper-pair mechanism. Multiplons were also recently studied theoretically in the 1D Hubbard model [1]. It was shown theoretically that the formation of the quadron or quadruplon was more favorable than the bi-exciton in a strongly confined parabolic quantum dot [2]. The situation is completely analogous in elementary particle physics where the meson molecules [10, 11] (analogous to bi-exciton or exciton molecule) and the genuine tetraquarks [4, 11,12,13] (the quadruplon analogue, see Fig. 1 in Ref. [11] as well as Fig. 6 later on) were both predicted theoretically [10] and observed experimentally [4, 12, 13]. In terms of cluster expansion language, bi-exciton and meson molecules are two weakly interacting irreducible clusters of order 2, or , while quadruplon and tetraquarks correspond to irreducible cluster of order 4, or . To date, experimental evidence for the existence of such a four-body (4B) irreducible cluster (or correlated entity) is still lacking in other fields beyond high-energy physics.

Fig. 1

Basics of experiment and key observations. a, Illustration of the e–h “soup” generated by the intense pump pulse, showing the possible 2B and 4B states. b, Schematic of the charge-tunable device by the gate voltage, V_g. c, d, CW RCS (R – R’)/R’ (see Methods S2.1 for definition) at several temperatures under the charge-neutral condition at V_g = –1 V (c) and under the p-doping condition at V_g = –6 V (d) (Device #1). The solid lines are the results of the Gaussian fittings. e, TDAS (see the main text for definition) at different temperatures for the cross- (σ- σ +) circularly polarized pump-probe configuration at the pump-probe delay time t ≈ 0 ps at the charge-neutral voltage, V_g = –1 V. Only spectral region below exciton is shown (compare with c). Each spectrum for a given temperature was fitted with 6 Gaussian peaks marked by P₁ – P₆. The points are experimental values, while the red solid lines are the results of the Gaussian fittings. f, Plots of the fitted central energies of the Gaussian peaks with respect to the temperatures for P₁ – P₆ (e) and T (d). The spectral locations of T marked with the green dashed lines in e were obtained from the CW results in d

Two-dimensional (2D) layered semiconductors, such as monolayer transition metal dichalcogenides (ML-TMDCs), provide a unique platform for the study of multi-particle complexes. The reduced dielectric screening in ML-TMDCs leads to extremely large excitonic binding energies [14, 15] and much more stable multi-particle complexes than those in conventional semiconductors. Indeed, trions [16, 17], bi-excitons [18,19,20,21,22], and even charged bi-excitons [23,24,25,26,27] were experimentally observed in ML-TMDCs with larger binding energies and at higher temperatures than those in bulk semiconductors. In addition, the unique spin-valley locking [28, 29] leads to more varieties of the above entities than in conventional semiconductors. For the same reason, the multi-particle entities with specific spin-valley polarizations are addressable by choosing the helicity of pump or probe light [21, 22]. The combined unique features described above have not been available in other materials for such a study. Therefore, ML-TMDCs with the above combined advantages serve as a unique platform for the further study of the high-order irreducible correlations or the entities.

To take advantages of these unprecedented opportunities, we conducted a combined theoretical and experimental investigation into the possible existence of higher-order correlated complexes in ML-TMDCs beyond the known trions, bi-excitons, and charged bi-excitons (including those that are reducible to excitons, trions, or bi-excitons). Using the helicity-resolved pump–probe technique, we observed a series of spectral features in the transient reflectance in a wide range of probe photon energies in gate-controlled monolayer molybdenum ditelluride (ML-MoTe₂) samples. Up to six spectral peaks were revealed, extending over 40 meV from below T all the way up to X. These spectral features cannot be attributed to defects or phonon-related processes. To understand the new spectral features, the four-body Bethe–Salpeter equation was solved to obtain the transition spectra [30] and found in good agreement with experimental results. By using the cluster expansion technique, interestingly, we found that the clusters corresponding to the trions and bi-excitons cannot produce most of the new spectral features, thus excluding the trions and bi-excitons as their origin. Importantly, such theory–experiment comparison shows that the 4^th-order irreducible 2e2h-cluster, or quadruplon, is necessary and sufficient in producing all the experimental spectral features, thus providing experimentally spectroscopic signatures of the existence of the quadruplons.

1.1 Experimental results

Figure 1a represents schematically the experimental situation where a strong pump produces an e–h “soup”: a combination of various excitonic complexes or correlated entities of various orders and their corresponding excited states. Figure 1b shows our device structure of a ML-MoTe₂ sandwiched between two hexagonal boron nitride (h-BN) layers, with fabrication details presented in Methods S1. A back gate was used to control the background charge of the ML-MoTe₂. The temperature-dependent reflection contrast spectra (RCS) of Device #1 are obtained from the continuous-wave (CW) reflectance spectra with (R) and without (R’) the sample, i.e. defined as (R – R’)/R’ (for the details of data processing see Methods S2.1, Fig. S1a – S1h, and Ref. [31]). The RCS are shown under the charge-neutral condition at V_g = –1 V in Fig. 1c and for the case of p-type doping at V_g = –6 V in Fig. 1d. The grey dots with an additional minus sign (i.e. –(R – R’)/R’) could represent the ground-state absorption (GSA) of the material. According to the fittings for the results at 4 K, the X peak is spectrally positioned at ~ 1.168 eV at the gate-compensated charge-neutral voltage, V_g = –1 V (Fig. 1c), and the T peak appears at ~ 1.149 eV in the doped regime at V_g = –6 V (Fig. 1d).

The GSA spectra show two simple features: 1) The single T peak for ML-MoTe₂ (see also previous papers on MoTe₂ [31,32,33,34,35]) corresponds to the inter-valley spin-singlet trion. This is closely related to the special band-structure of ML-MoTe₂ (or ML-MoSe₂): i.e. intravalley bright exciton has a smaller energy than intravalley dark exciton [27, 33] (in contrast to ML-WS₂ and ML-WSe₂ [23,24,25,26]), and the spin–orbit (SO) splitting of conduction band (CB), ~ 30 – 60 meV [32], is sufficiently large (in contrast to ML-MoS₂, ~ 3 – 15 meV [36, 37]). 2) More importantly, there are no additional absorption features below X in the charge-neutral regime without a pump (see the low-energy side of X in Fig. 1c). The simple GSA spectra can be used as a base or reference for later study of complicated features in the transient differential absorption or reflection spectra (TDAS or TDRS). The fact that there is no more complicated feature below T and X makes the MoTe₂ one of the ideal systems among various TMDCs (for more specific reasons or advantages that we chose to study MoTe₂ see Methods S3).

Typically by measuring the CW RCS and photoluminescence (PL) spectra (at both low and high excitation densities), we pre-screened the samples for the following pump-probe experiments and selected those without visible defect features. Representative CW RCS and PL spectra of ML-MoTe₂ could be seen with or without defect peaks below T, as presented in Methods S2.2, Fig. S2 & S3.

The ultrafast pump-probe experiment is described in Methods S2.3, with the experimental setup illustrated in Fig. S4. The pump energy of 1.174 eV is ~ 6 meV above X with a fluence of ~ 60 μJ·cm^–2 (corresponding to an e–h pair density n_p estimated to be ~ 3 × 10¹² cm^–2, see Methods S2.4 for the estimate), while the probe energy was tuned from 1.127 to 1.161 eV with a resolution of ~ 0.2 meV. The gate voltage was set at –1 V to maintain charge neutrality. The temperature-dependent transient differential absorption spectra (TDAS) of Device #1 are shown in Fig. 1e. Here, the TDAS (–Δα) is defined as –Δα = –(α_p – α₀) ∝ R_p – R₀, where R_p and R₀ are the RCS from the sample with and without pump (α denotes the corresponding absorption. See Methods S2.3 and Ref. [34] for more details about the relation between TDAS and TDRS). It is worth noting that –Δα < 0 means a pump-induced absorption increase, typically related to those of excited-state absorption (ESA) processes. In Fig. 1e, we observe rich spectral features with absorption increases as marked with P₁ – P₆. Based on their relative positions to T and X, these features and the associated fitted peaks can be divided into three energy intervals: P₁, P₂ (below T), P₃, P₄ (near but below T), and P₅ & P₆ (above T or between T and X).

To obtain more quantitative information about these features, we performed multi-peak Gaussian fittings on the spectra at various temperatures in Fig. 1e (The applicability of such multi-peak fittings for P₁ – P₆ in the TDAS will be discussed then in connection with Fig. 2l & m). The central energies of the obtained Gaussian peaks P₁ – P₆ (Fig. 1e) and the fitted peak T (Fig. 1d) are plotted in Fig. 1f versus temperature. The intervals between the neighboring peaks in P₁ – P₆ are in the range of 4 – 7 meV, while the total spread of these peaks is around 25 meV. The lowest peak (P₁, ~ 1.134 eV at 4 K) is about 35 meV below the original X peak (~ 1.168 eV at 4 K) and ~ 15 meV below the T peak (~ 1.149 eV at 4 K). When temperature increases from 4 to 80 K, P₁ – P₆ in the ultrafast spectroscopy (Fig. 1e) show a redshift of ~ 1 – 3 meV, while peaks T and X extracted from the CW results (Fig. 1c & d) show a redshift of ~ 7 meV. As can be seen in Fig. 1e, P₁–P₆ are well resolved at 4 K, but merged more together with increase in temperature due to increased broadening with temperature. To show the reproducibility of the results (the 4 K case in Fig. 1e), we measured the TDAS with a finer spectral resolution of 0.1 meV. The six similar peaks can be seen quite clearly even without the multi-Gaussian fitting, as presented in Fig. S5, Methods S4.

Fig. 2

Gate-voltage and polarization dependence of the TDAS. a, CW RC contour in the plane of photon energy and gate voltage (Device #2 at 4 K). b–k, TDAS contours in the plane of probe energy and delay time for the cross- (σ− σ +) (b – f) and co- (σ + σ +) (g–k) circularly polarized pump-probe configurations. The voltages (5 V, 1 V, –3 V (charge-neutral), –7 V, –11 V) applied for observing the TDAS in b–k are marked with the five cyan dashed lines in a. In d & i, the similar new spectral features (P₁–P₆) to those in Fig. 1e are marked accordingly. l–m, Schematics of the relationship between the absorption spectra (upper panels) and the TDAS (lower panels) under the charge-neutral condition. In the upper panels, the curves filled with red and blue denote the absorption spectra with (α_p) and without (α₀) pump, respectively. The zoomed-in inset in l shows the X with (X_p) and without (X₀) pump to illustrate GSB and BGR. l & m show the GSA (l) and ESA (m) cases. ESA includes the potential many-body effect of high-order correlation, thus corresponding to the presence and absence of P₁ – P₆ in m and l, respectively

We notice that previous studies have also observed a peak below T attributed to charged bi-excitons [23,24,25,26,27]. To examine the charge dependence of our new spectral peaks, we performed a gate-dependent experiment on a device of the same design (Device #2) at 4 K. The gate-dependent CW RC (or absorption) map is shown in Fig. 2a. Similar to the case of Device #1 (Fig. 1c), there are no observable features below X in the charge-neutral regime (V_g = –3 V) without pump. A few selected TDAS are shown in Fig. 2b – k. From top to bottom, the system was gated into the charge-negative, neutral, and positive regimes corresponding to the five voltages as marked by the cyan dashed lines in Fig. 2a. In Fig. 2b – k, the spectral features of P₁ – P₆ similar to those in Fig. 1e are visible and are the strongest in the charge-neutral situations of V_g = –3 V (Fig. 2d & i). The features fade away as the system deviates from charge neutrality. Such gate-dependent behavior was also observed for another device of the same design (Device #5) (see Methods S6, Fig. S9 for details). Contrary to the previous observations of the PL of charged bi-excitons [23,24,25,26] that are the weakest in the charge-neutral regimes, our new peaks P₁ – P₆ are the strongest in these regimes, pointing to the existence of charge-neutral entities.

As can be seen from the polarization-dependent results, the spectral features are stronger for the case of (σ– σ +) (Fig. 2d) than for the case of (σ + σ +) (Fig. 2i) and this is also true for two other devices we measured (see Methods S4, Fig. S5 for Device #1, and Methods S5, Fig. S6 for Device #4). Compared to P₁ – P₆, the signals around X show less such a polarization contrast. The polarization contrast for P₁ – P₆ was also observed for bi-exciton signals in previous experiments [21, 22, 27]. It means that the inter-valley (σ– σ +) configuration is always more favorable than the intra-valley (σ + σ +) one [21, 22, 27], which is also similar to the case of bounding and anti-bounding states of a Hydrogen molecule. Such a polarization contrast is also consistent with the recent theoretical results [30]. In contrast to the six peaks observed here, previous studies have shown one peak for BX [18,19,20,21, 23,24,25,26,27, 38,39,40] and no more than three peaks for fine structure of BX (BXFS) [22, 41]. The spectral positions (or the corresponding binding energies) of BX and BXFS have been calculated for ML-TMDCs previously [22, 38,39,40,41], and widely accepted to be between T and X (14.4 meV below X for ML-MoTe₂ [38]). Therefore, the sequence of our experimental spectral features, i.e. P₁ – P₆, which extend ~ 40 meV from below T all the way up to X, could not be explained by BXFS.

To better understand the new spectral features P₁ – P₆ & X in Fig. 2d & i (also Fig. S6), we schematically illustrate the relationship between the absorption spectra (α_p & α₀) and the TDAS (–Δα), as shown in Fig. 2l & m. Under the charge-neutral condition without pump, α₀ has only a single X peak, X₀ (the corresponding RCS can be seen in Fig. 1c). Under optical pumping, the X peak in α_p (X_p) shows a reduced oscillation strength due to ground-state bleaching (GSB) and a redshifted resonance energy due to bandgap renormalization (BGR), as illustrated in the inset in Fig. 2l. To obtained TDAS, a subtraction of α₀ from α_p with a redshifted X peak leads to the typical anti-symmetric feature shown in the lower panel of Fig. 2l (also Ref. [42]) schematically or in actual experimental data (see Methods S5, Fig. S6). In addition to GSB and BGR, a strong pump could lead to the excited-state absorption (ESA) corresponding to transitions from existing 2B states to 4B entities such as bi-exciton [21, 22, 39, 40], etc. Figure 2l & m compare GSA and ESA cases. For the ESA case, new peaks in α_p emerge with pump in addition to X. A subtraction of α₀ from α_p produces the corresponding pump-induced peaks in the TDAS, i.e. P₁ – P₆ as shown in the lower panel of Fig. 2m (see also Ref. [21, 22, 39, 40] for the BX peak in TDAS). In other words, peaks (P₁ – P₆) resulting from ESA remain in TDAS as in the absorption spectrum, in contrast to the anti-symmetric feature that corresponds to GSA processes.

To fit the TDAS, we assign six negative Gaussian peaks for P₁ – P₆ and a combination of negative (for X_p) and positive (for X₀) peaks for the antisymmetric line-shape around X (see Methods S5 for the fitting details, the fitting method can be seen also in Ref. [22]). We notice that X_p is red-shifted from X₀ by ~ 5 – 7 meV, while peaks P₁ – P₆ are red-shifted by from X_p ~ 6 – 40 meV. Since peaks P₁ – P₆ are sufficiently narrow and well-separated from X_p, they can be independently fitted.

To study the pump-density dependence of P₁ – P₆, we performed a series of pump-probe experiments with a varying pump fluence on a device of the same design (Device #3) at 4 K. Based on the results of Device #2 presented in Fig. 2, the system was gated into the charge-neutral regime to have the maximum effects of the spectral features. The pump-fluence-dependent TDAS are shown in Fig. 3a. From top to bottom, the pump fluence was varied from ~ 4 μJ·cm^–2 to ~ 320 μJ·cm^–2 (n_p estimated to be 2.0 × 10¹¹ cm^–2 to 1.6 × 10¹³ cm^–2, see Methods S2.4 for the estimates). Figure 3b shows the spectra in x–y plot (corresponding to those of Fig. 3a at the delay time of ~ 0.9 ps, where the signals are the strongest) with the multi-Gaussian fittings similar to those shown in Fig. 1e. The integrated intensities and central energies of the fitted Gaussians are plotted in Fig. 3c & d, respectively (We re-plotted Fig. 3c in a linear scale, as can be seen in Methods S7). As can be seen in Fig. 3a & b, the spectrum is nearly featureless when the pump density is below 4.0 × 10¹¹ cm^–2. As the pump density is above 4.0 × 10¹¹ cm^–2, P₆ (between T and X) is the first to appear in the spectrum, followed by P₂ – P₄ (below or near T) & P₅ (between T and X) with blurred features. With the further increase of the pump density beyond 1.6 × 10¹² cm^–2, P₂ – P₄ & P₅ become increasingly visible and better distinguishable, and P₁ (below T) starts to appear (better visible in Fig. 3b than in Fig. 3a). When the pump density reaches 8.0 × 10¹² cm^–2, all the peaks of P₁ – P₆ are visible and can be fitted within relatively small errors (Fig. 3c & d). At any pump level as shown in Fig. 3c, P₁ is always the weakest among the 6 peaks while P₆ is the strongest. The pump dependence described here will be further explained in connection with the discussions of Fig. 5. With the increase of the pump density, P₁, P₂, & P₆ exhibit larger redshifts of ~ 2.5, 2.2, and 2.8 meV than P₃ – P₅ of ~ 1.5 meV or smaller (Fig. 3d).

Fig. 3

Pump-density dependence of the TDAS. a, TDA contours in the plane of probe energy and delay time for the cross- (σ− σ +) circularly polarized pump-probe configuration (Device #3, measured at 4 K in the charge-neutral regime). b, Gaussian fittings (red solid lines) of the TDA extracted from a (dots) at a delay time of ~ 0.9 ps. Each spectrum is fitted with 6 Gaussian peaks marked by P₁ – P₆ as in Fig. 1e. The spectral location of T is marked with the green dashed lines in a & b. The total e–h pair density, n_p, generated by the pump was marked in each panel of a & b. c, d, The total areas (corresponding to absorption increment) (c) and central energies (d) of P₁ – P₆ with respect to the above pump densities

Spectral features represented by P₁ – P₆ as we characterized so far have not been seen before in other semiconductors. In general, many effects, intrinsic or extrinsic, could cause additional spectral peaks such as defects [18,19,20, 23, 26], phonons [43,44,45,46] (e.g. exciton phonon replicas, etc.), any emission signals, second-harmonic generation, and other non-linear mixing effects typically related to e.g. χ⁽³⁾, χ⁽⁵⁾ or higher order contributions, or charge-neutral entities involving more than two excitons or dropletons. We have analyzed all the above effects in detail (see in Methods S8) and concluded that it is unlikely that these effects could cause the new spectral peaks. The spectral peaks have been observed in 5 of the best samples with the least amounts of defects and the spectral features are robust against change of spectral resolutions from 0.1 meV to 2 meV. The significant difference between the cross- and co-polarized pump-probe configurations is also a strong indication that many of the above origins can be excluded. In addition, the good agreement between the experiment and the theory [30] as presented in the following also favors the intrinsic origins of 4B states, since the theory does not include any of the above effects. Of course, other hitherto unforeseen origins or effects should not be ruled out. Continued searches for new unique theoretical and experimental evidence for the quadruplon are of great importance.

To understand the origin of P₁ – P₆, we resort to modern theories of many-body interactions for the calculation of the absorption spectra. We notice that the absorption changes in the presence of a strong pump necessarily involve transitions between two-body states and four-body states. The Bethe–Salpeter equation (BSE) is typically used to calculate the wavefunction of a many-body system. In the following, the two-body and four-body wavefunctions will be calculated by using the BSEs and then used to calculate transition spectra [30] corresponding to ESA.

1.2 Four-body interactions and the cluster expansion approach

Another alternative to BSE and more intuitive approach is the cluster expansion technique, the most natural way of visualizing correlated entities of various orders and associated many-body interactions [47]. This theoretical method has more recently been proven successful in describing the many-body correlations in e–h dropletons [6, 7]. Generalizing such cluster expansion to the case of individual electrons and holes with the spin and valley degrees of freedom, we can write down the complete sequence of irreducible clusters for the 2e2h system as in Fig. 4a–c for ML-MoTe₂. As also described in Ref. [30], represents an irreducible cluster of the n^th-order with n Fermions. Therefore, is a quasi-free electron or hole or a singlon. represents an e–h pair, or a 2-body (2B) state or a doublon, including all the excitonic Rydberg series: 1s-X, 2p-X, 2s-X, …. represents an e-e–h or e–h-h 3-body (3B) state or a triplon. Note that we used the language of multiplons proposed by Rausch et al [1] to name , and , doublon, triplon, and quadruplon, respectively. Each time we include one more cluster of higher order into a truncated expansion, it re-introduces weak interactions (indicated by the wavy lines in, e.g. Fig. 1a & Fig. S11) among irreducible clusters of all the lower orders, as explained in more details in Methods S9. Here, the cluster expansion is truncated up to the 4th order. In the presence of higher-order clusters, the original non-interacting cluster shown in Fig. 4b, c (also Fig. S11d) becomes weakly interacting, i.e. ~ (Fig. S11f), representing a bi-exciton configuration, not only (1s-X) ~ (1s-X) but also all its excited states such as (1s-X) ~ (2p-X), (1s-X) ~ (2s-X) …. Obviously, cluster or quadruplon (or quadron [2]) represents generally a distinct physical entity than the bi-exciton ~. The most fundamental difference between a bi-exciton and a quadruplon is the absence of an exciton in the latter and the lack of a clear association of any one of the two electrons to a given hole. But it is easy to imagine that a certain disassociation event (e.g. an excitation) of a quadruplon could possibly lead to the formation of a bi-exciton. In this sense, a bi-exciton could be an excited state of a quadruplon.

Fig. 4

A cluster expansion picture of a four-body system. a, Decomposition of the 2e2h 4B system into three parts according to the time reversal (K ↔ K’) symmetry, as indicated by the blue-red symmetry of the square boxes. The first two terms/boxes are mutually exchanged under the time reversal and only term is presented in detail in b. The third term/box is time reversal invariant, as presented in c. b, c, Cluster expansions of the 4B system into irreducible clusters (represented by the filled or unfilled circles connected by the straight lines) of various orders (sizes). The spin-up and spin-down electrons are colored in blue and red, respectively. For brevity, we exclude those 2B clusters with the same charges, such as e-e and h–h, in the figures

1.3 Theory–experiment comparison

Figure 5a presents the absorption spectrum calculated using the three-body Bethe–Salpeter Equation (3B-BSE) [48] for ML-MoTe₂, where we see clearly only the T peak in its neighborhood. Obviously, such 3B-BSE does not explain the features related to P₁ – P₆. The BX-related spectral features have been calculated theoretically [22, 38,39,40,41] to be between T and X. Specifically, the BX resonance of 14.4 meV below X for ML-MoTe₂ was calculated to be between T and X by combining the density functional theory with the path-integral Monte Carlo method [38]. But so far, there has been no experimental observation of the BXs in MoTe₂. In addition, our spectral features P₁ – P₆ contain more spectral peaks in a much larger spectral range than those BX-related peaks [22, 39,40,41]. Following Ref. [30], we calculate the two-body wavefunction, $\left|{e}_{3}{h}_{3}\right.\rangle$ and the four-body wavefunction $\left|{e}_{1}{h}_{1}{e}_{2}{h}_{2}\right.\rangle$ by solving the corresponding BSE. The dielectric function for the 2B-4B-transitions can be then calculated by summing over all the transitions as follows,

$$\varepsilon \propto \frac{2\pi }{\hbar }\mathop \sum \limits_{\alpha } \left| {\mathop \sum \limits_{{\left( {v_{i} ,c_{i} ,i = 1,2,3} \right)}} \left. \langle {e_{3} h_{3} } \right|{\mathbf{e}} \cdot \overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{p} \left| {e_{1} h_{1} e_{2} h_{2} } \rangle \right._{\alpha } } \right|^{2} \Gamma$$

(1)

Fig. 5

Theory–experiment comparison of absorption spectra. a–e, Absorption spectra calculated based on the 3B-BSE (a) and the 4B-BSE (b – e) for ML-MoTe₂. d shows the spectra of transition between the 2B state (1 s-X) and the 4B states solved from the full 4B-BSE. b, c, & e show the cases with the same 2B state (plasma), but the 4B states solved from the 4B-BSE truncated up to (b), (c), and (e). The vertical dashed line is the calculated trion (T) energy. The vertical black lines underneath the spectral profiles mark the calculated spectral positions with the height representing the strength of the dipole transitions. The broadening parameter, γ, chosen to be 0.5 meV (unfilled profiles) or 2 meV (shaded profiles). We label the Gaussian-broadened peaks in each spectrum with italic numerals from low to high energy. The solid black profile in f represents Δα at the zero delay for Device #2 at 4 K in the charge-neutral regime, where we also overlay the contour of –Δα in the plane of probe energy and delay time (see also Fig. 2d). The features are labelled with P₁ – P₆ (and X), consistent with those in Fig. 2d. The relationship between the total-energy spectra and the optical spectra is shown in g – j (for the explanation in more details see the main text). The calculated total energies are shown beside the upward axis in g for those typical 2B and 4B states (not to the scale). The energy range of 1.432–1.770 represents the continuous absorption band. Transition ① and ⑥ correspond to d & e, respectively. The flag notation in g marks the 4B ground state, whose spectral feature corresponds to peak 1 in d. b, c, d, e, & g are reproduced from Fig 2, 4 in Ref. [30]

Where $\Gamma$ is a lineshape function taking into account of the broadening of the transitions. For the 2B states, there are a few discrete states from 1s-X, 2p-X, 2s-X, … plasma for ML-MoTe₂. For the 4B states, the following three cases of truncation (in connection with Methods S9, Fig. S11) were considered there, enabling to extract and identify the effects of each: Case 1, up to ; Case 2, up to ; Case 3, up to . For each case, a series of 4B states are solved from the corresponding truncated 4B-BSE and then used to calculate the dipole matrix elements of all the possible 2B-4B transitions.

To obtain a correct and complete picture of the possible 2B-4B transitions, especially the relative positions of spectral lines and their origins, we emphasize here the difference between the typical optical spectrum and the total-energy spectrum. Such total-energy spectrum of the 2e2h 4B system is shown in Fig. 5g for the 4B system truncated to . It is important to realize that the optical spectrum measured in an absorption or emission experiment does not reflect the absolute energies on the total energy scale. As shown in Fig. 5g, a series of transitions occur between different 2B states and the corresponding 4B states on the total energy scale, as marked by the vertical double-arrowed lines with ① – ⑥ (6 examples, selected out of all the 2B states). The actual spectrum of the transitions corresponding to the superposition (Fig. 5j) of all possible individual transitions at different total energies are shown schematically in Fig. 5g. Unfortunately, such total-energy spectra are not easy to obtain in an optical experiment, where only the energy difference (photon energy) between the initial and final states of these transitions is measured. Or equivalently, all these total energy spectra are shifted relative to a common reference (e.g. the 1s-X energy), as shown in Fig. 5i. The actual optical spectrum is the superposition (Fig. 5j) of 15 of these spectra shown in Fig. 5i, or 15 vertically “collapsed” spectra shown in Fig. 5h.

Figure 5b & c show the optical spectra with the 4B states calculated from the 4B-BSE truncated up to (Fig. 5b) and (Fig. 5c) and the 2B state given the plasma state. We notice only one main peak between T and X in the case of (peak 1 in Fig. 5c) versus the featureless background in (Fig. 5b). Clearly, we do not see any spectral features corresponding to P₁ – P₄, especially the features below T. This means that P₁ – P₄ do not originate from cluster or .

Our next approximation is to include all the 4B states calculated from the full 4B-BSE (truncated up to ). The relationship between the total-energy spectra and the optical spectra in this case is shown in Fig. 5g–j. The optical spectra calculated with the 2B states given the 1 s-X (transition ①) and plasma (transition ⑥) are reproduced from Fig. 2 in Ref. [30] and shown here in Fig. 5d & e, respectively. A direct comparison is shown in Fig. 5f (It is very important to point out that the experimental result should not be compared with only one of these calculated spectra, because the actual spectrum is a superposition of all possible 2B-4B transition spectra as we illustrated in Fig. 5g–j). We see that the spectral features (peak 1 – 6 in Fig. 5e) in the entire spectral region resemble closely to those of P₁ – P₆. In other words, the spectral features of P₁ – P₄ originate from cluster . By comparing various 2B-4B transitions in Fig. 5g–j for cluster , clusters , and both on the total-energy scale and optical spectrum, we could make the following conclusion: The states from cluster for each of the same 2B states have lower energies and are thus more stable than those of ~ (BX) and ~ (T⁺ ~ e and T^– ~ h), indicating the most stable existence of cluster or quadruplon. Importantly, we notice from the transitions in the total-energy scale that the lower-frequency features (below T) in the optical spectrum such as P₁ – P₃ correspond to the transitions (such as ④, ⑤, and ⑥ in Fig. 5g) between the more-excited states in the 4B manifold. This explains why P₁ – P₃ decays faster as we mentioned in connection with the discussions of Fig. S6a & S6d for Device #4. Nearly all the 2B-4B transitions shown in the total-energy spectrum in Fig. 5g, e.g. from ① (low energy) to ⑥ (high energy), can contribute to the spectral peaks between T and X. Especially, the spectral peak related to the 4B ground state is calculated to be between T and X (close to P₅ or P₆, see the flag notation in Fig. 5d). This explains why P₆ has the largest intensity and appears the earliest with increasing pump (Fig. 3). However, only the transitions between the highly excited 2B and 4B states such as those of ③ – ⑥ in Fig. 5g, can contribute to the peaks well below T. The fact that P₁ only emerges at high pump levels signifies the existence of highly excited states of both the 2B & 4B entities. Through such 2B-4B spectroscopy, we observe a much more complex many-body system with more refined interplays of these 2B, 3B, and 4B complexes than the simple picture of the Mott transition. The similar spectra (to e.g. Figure 5d & e) were also calculated for the inter− (σ–σ +) and intra−(σ + σ +) valley 4B states in Ref. [30].

One more experimental fact supporting the picture of 2B-4B transitions is that the spectral positions of P₁ – P₆ are much less sensitive to temperature than the T and X peaks (Fig. 1e & f). We notice that the transition energies of T and X follow the same sensitive temperature dependence as bandgap, given by the Varshni formula. The total energy of a 4B system also follows a similar temperature dependence. But since the transition energies of P₁ – P₆ are determined by the differences between the total energies of the 4B system and the 2B system, the difference becomes much less sensitive to temperature change. This explains the much weaker temperature dependence of the P₁ – P₆.

It is important to note that represents a 2e2h 4B irreducible cluster, where none of the 2 electrons or 2 holes belongs to a specific exciton, or there is not a well-defined exciton in such a 4B system. The entity is therefore not a bi-exciton or an excited-state bi-exciton anymore, rather a quadruplon, as schematically shown as the last cluster in Fig. 4b – c. We notice that such a system was recently calculated for a quantum dot system and was called a “quadron” [2]. Therefore, the consistent theoretical and experimental results show that the spectral features corresponding to P₁ – P₄ indeed originate from the quadruplons. Similar quadruplon consisting of 4 identical fermions was recently theoretically studied in a 1D Hubbard model [1]. Recent study showed that the e–h/e-e exchange interaction [22, 41] can result in fine structure of bi-exciton. Therefore, quadruplon features from irreducible cluster should not have been called BXFS [49] (see [49] for a more detailed explanation).

2 Conclusion and outlook

The focus of this paper is the experimental observation of several new spectral peaks that appeared below the trion features and the theoretical explanations using a microscopic theory based on the 4B-BSE. The observation of these peaks covering a large energy range of ~ 40 meV was made experimentally possible through a helicity-resolved pump-probe spectroscopy in gate-controlled ML-TMDC samples. It was at the charge-neutral point that we observed the new rich spectral features. This might be one of possible reasons that they were not observed in previous ultrafast-spectroscopic experiments on ML-TMDCs (for more details see Methods S11).

By using the correspondence between the cluster expansion technique and the BSE approach, we were able to pinpoint the roles of the fourth-order irreducible cluster.

These new spectral features could not be explained theoretically by clusters associated with trions or bi-excitons, while the irreducible cluster of order 4 associated with 2e2h was necessary and sufficient in producing these new spectral features. The agreement between theory and experiment and the roles played by the 4th order irreducible cluster provide strong evidence for the existence of a new correlated 4B entity, the quadruplon. The existence of the 4B entity is further corroborated by the temperature and pump dependent experiment where systematic changes of the absorption spectra are consistent with the possible occupation of higher 2B states and the existence of the final 4B states upon the absorption of the probe photons. Interestingly, these 4B states can stably exist even at a relative high temperature of ~ 60 K.

We would like to make a few comments about other multi-particle states that are products of lower order clusters. Examples include the Suris tetron [50,51,52] that contains a conventional exciton and another e–h pair formed near the Fermi surface, or ~ which are bi-exciton-like, or trion-hole-like [53] entity in the language of this paper. Other examples include the quarternion [54], charged bi-exciton or trion-exciton [23,24,25,26,27], and even 6B and 8B ones: the hexciton and oxciton [55], or tri-exciton [56] ~~ and quad-exciton ~~~ and dropleton ~~~ … ~ as exciton ladders, etc. All these many-body complexes are reducible to product states of excitons, trions, and charges. No association with higher irreducible clusters than 3 has been made. Thus the possibility of 4th order irreducible clusters has not been considered so far [23,24,25,26,27, 50,51,52, 54,55,56]. The difference between product of lower order irreducible clusters and the corresponding higher order irreducible clusters is also important in the field of elementary particles, where meson molecules and genuine tetra-quarks represent separate milestone discoveries [4, 10,11,12,13], as shown in Fig. 6. Since quarks have more quantum numbers and the quark system is much more complex than the electron–hole system, there are many more possibilities of tetraquarks. But one particular type of tetraquark, $(q\,q\,\overline{q}\,\overline{q})$: (such as $c\,c\,\overline{c}\,\overline{c}$) [57] resembles the quadruplon the most. The $2q2\overline{q}$ four-quark system can form either a meson molecule $(q\,\overline{q}) \sim (q\,\overline{q})$ or a genuine tetraquark $(q\,q\,\overline{q}\,\overline{q})$. In the case of quadruplon, both spins and charges are balanced, while in the case of tetraquarks of the type $(q\,q\,\overline{q}\,\overline{q})$, both charges and spins are balanced as follows: (e,e,-e,-e), (1/2,1/2,−1/2,−1/2) for quadruplon and (2e/3,2e/3,-2e/3,-2e/3), (1/2,1/2,-1/2,-1/2) for tetraquark. In the context of ultracold gas, the 4B states (related to the Efimov effects) were theoretically discussed [57] in terms of all irreducible clusters such as B₄, B₃ + B, B₂ + B + B, B₂ + B₂, or B + B + B + B. Attempts were also made to explain an earlier experiment [58] in terms of the above theory [57]. Also, the quadruplons consisting of 4 holes were theoretically studied in the sequence of Q (quadruplon), T (triplon), D (Doublon), and B (Band-like part) in the calculated two-hole excitation spectrum at different filling factors in a 1D Hubbard model [1]. Finally, we mention that a more recent study during the review of this paper provides interesting evidence of 4e and 6e superconductivity [59], although the issue of irreducibility remains to be addressed.

Fig. 6

| A schematic similarity between two configurations of $2\,q\,2\,\overline{q}$ four-quark system (left) and 2e2h electron–hole system (right). For the quark system, the filled circles $(q)$ and the unfilled circles $(\overline{q})$ represent quarks and anti-quarks, respectively. For the electron–hole system, the filled and unfilled circles represent electrons and holes, respectively. The schematic representation of the tetraquark and the meson-molecule were reproduced from Liu’s paper [11]. We illustrated the 2e2h entities by using the Jacobi-like coordinate. The similarities between tetraquarks and quadruplon and between Bi-meson (or meson molecule) and bi-exciton (also called exciton molecule) becomes apparent

To conclude, the observation of quadruplons would also contribute in a very important way to the fundamental understanding of the Mott transition, one of the most celebrated aspects of condensed matter physics. The existence of a new 4B entity adds to the complexity of the Mott physics and would definitely lead to more new physical phenomena. For example, the co-existence and mutual conversion of various known species have been the subjects of extensive studies for their interesting optical transitions, leading to bi-exciton gain [60] or trion gain [31]. Bi-exciton emission processes are also related to the generation of entangled photons. It would be of great interests to study the co-existence and mutual conversion of quadruplons with other known species. The existence of quadruplons opens many new exciting opportunities to study its consequences on optical gain, generation of quantum states, and nonlinear optics, beyond the Mott transitions.

Finally, our study may stimulate more similar studies to search for clusters or many-body complexes of even higher orders beyond the quadruplon and their excited states. Another interesting issue that arises from this study is the experimental study of the total energy spectrum of a 2e2h 4B system. As shown in Fig. 5g, the total energy of the 2e2h system extending over ~ 0.3 – 0.6 eV shows very rich spectral features that cannot be determined through optical spectroscopy. It would be highly interesting to see if techniques such as angle-resolved photoemission spectroscopy (ARPES) could provide more direct evidence of quadruplons in the total energy spectrum. A comparison of our 4B-BSE theory with the total energy spectrum would allow us to study many different states of the 4B irreducible clusters in a more unique fashion. More importantly, it would allow the determination of the relative energetic stabilities of various states. Among all of these states, the experimental verification of the relative stability of bi-excitons versus quadruplons would be of great special interest. Importantly, quadruplon involves very different roles played by individual electrons and holes from bi-exciton. It would be of great interest to see how these roles are reflected in emission processes and properties, especially related to the quantum statistical properties or entanglement of emitted photons. Finally, similar to bi-exciton gain, possible existence of optical gain involving quadruplon as the excited state would be another important subject for future studies.