Molecular origins of exciton condensation in van der Waals heterostructure bilayers

Abstract

Recent experiments have realized exciton condensation in bilayer materials such as graphene double layers and the van der Waals heterostructure MoSe₂–WSe₂ with the potential for nearly frictionless energy transport. Here we computationally observe the microscopic beginnings of exciton condensation in a molecular-scale fragment of MoSe₂–WSe₂, using advanced electronic structure methods based on reduced density matrices. We establish a connection between the signature of exciton condensation—the presence of a large eigenvalue in the particle-hole reduced density matrix—and experimental evidence of exciton condensation in the material. The presence of a “critical seed” of exciton condensation in a molecular-scale fragment of a heterostructure bilayer provides insight into how local short-range strongly correlated effects may give rise to macroscopic exciton condensation. We find that molecular-scale properties such as layer alignment and interlayer distance can impact the formation of nonclassical long-range order in heterostructure bilayers, demonstrating the importance of geometric considerations for the rational design of exciton condensate materials. Mechanistic insights into the microscopic origins of exciton condensation have potential implications for the design and development of new materials with enhanced energy transport properties.

---

Introduction

Current power transmission technologies are plagued by energy losses due to frictional dissipation of energy, and the search for materials capable of frictionless energy transport at ambient pressure and temperature is ongoing. Excitons, quasi-bosonic bound states of electrons and holes, are capable of quantum condensation. The resulting superfluid is theoretically capable of nondissipative energy transfer,^1,2 which could inspire novel electronic devices and spur tremendous innovation to efficient energy transfer applications. Additionally, exciton condensation is predicted to be possible at higher temperatures than traditional superconductivity.³ While the condensation is difficult to achieve because excitons are susceptible to recombination, especially at room temperature, it has been experimentally realized through the coupling of excitons to polaritons,^3,4 and the spatial separation of electrons and holes in bilayer materials.^5–8 Bilayer systems provide an important platform for exciton condensation due the spatial separation of electrons and holes between layers preventing quick recombination of excitons. Graphene bilayers have proven to be promising candidates for exciton condensation, with the twist-angle dependence of their electronic states affording an avenue for tuneability.^5,7,9–11 Recently, there has been experimental evidence for exciton condensation in bilayers composed of transition metal dichalcogenides such as MoSe₂–WSe₂ heterostructures.^8,12–14 Additionally, recent theoretical work has shown the possibility for the beginnings of exciton condensation in several small-scale molecular systems.^15–19

Here we computationally observe exciton condensation emerging in molecular-scale fragments of a family of transition metal dichalcogenide bilayers. Using variational two-electron reduced density matrix (2-RDM) theory,^20–34 we calculate the populations of the excitonic states—the eigenvalues of the particle-hole reduced density matrix (RDM)—for molecular-scale fragments of several transition metal dichalcogenide bilayers. In the noninteracting limit, all eigenvalues of the modified particle-hole RDM are equal to zero or one, and no condensation is present. In the strongly correlated limit, off-diagonal long-range order emerges in the particle-hole reduced density matrix as a consequence of long-range coherence, leading to eigenvalues greater than one. The presence of a maximal exciton population λ_G greater than one, indicating more than one particle-hole pair in a given excitonic mode, is thus the theoretical signature of exciton condensation.^15,35 Exciton condensation requires both stable, long-lived excitons (such as interlayer excitons) as well as strong correlation for long range order in the material. Factors that disrupt the formation of stable excitons or the coherence of the system will disrupt the exciton condensate and cause a decrease in λ_G, making λ_G ideal for probing the effects of geometric modifications to an exciton condensate.

We find that this signature begins to appear in fragments of MoSe₂–WSe₂, a material which has previously shown experimental evidence for exciton condensation.^8,12,13 The degree of condensation increases with the size of the molecular fragment, suggesting that we are seeing a “critical seed” of the condensation that would be present in a bulk system. Additionally, we show that layer alignment and interlayer distance influence the extent of exciton condensation. Finally, we find that the signature is also present in several transition metal dichalcogenide bilayers similar to MoSe₂–WSe₂. The presence of a “critical seed” suggests that local, short-range strongly correlated effects are significant for the formation of a macroscopic condensate. Furthermore, our results highlight the importance of geometric considerations in the design of exciton condensate materials. The presence of the signature of exciton condensation in several different van der Waals heterostructure bilayers indicates that the identity of the particular transition metal dichalcogenide does not have a large effect on the potential for exciton condensation, suggesting that materials capable of exciton condensation may be more ubiquitous than previously thought.

Results

The transition metal dichalcogenide bilayer MoSe₂–WSe₂ is probed for the signature of exciton condensation—λ_G > 1, the number of excitons condensed into a single particle-hole state. Two cluster models are constructed for each system, consisting of bilayers with the primitive cell tiled 2 × 2 and 3 × 3 in the XY plane respectively (Fig. 1). The signature of exciton condensation is obtained for each molecular structure via variational two-electron reduced density matrix theory (v2RDM) as described in the Methods.^20–34 The theoretical density of the excitonic mode corresponding to the largest eigenvalue is plotted to determine if the system is exhibiting the desired interlayer exciton. To gauge the influence of layer alignment and interlayer distance on the extent of exciton condensation, λ_G is calculated for a variety of different layer alignments and interlayer spacings. Finally, a family of related transition metal dichalcogenide bilayers is probed for the signature of exciton condensation.

Fig. 1 Example 2 × 2 and 3 × 3 cluster models for the MoSe₂–WSe₂ heterobilayer, with molybdenum atoms in light blue, tungsten atoms in darker blue, and selenium atoms in orange. Passivating hydrogens omitted for clarity. Black outline indicates one primitive cell.

Exciton condensation in MoSe₂–WSe₂

MoSe₂–WSe₂ is a transition metal dichalcogenide bilayer that has shown strong experimental evidence of exciton condensation.^8,12,13 Calculation of the largest eigenvalue of the particle-hole reduced density matrix for MoSe₂–WSe₂ for the 2 × 2 and 3 × 3 models results in a λ_G of 1.29 and 1.40, respectively. The presence of λ_G > 1 indicates the potential for exciton condensation, aligning with previous experimental studies. While the exciton population in a single quantum state is quite small, the increase in λ_G as system size increases from the 2 × 2 to 3 × 3 model suggests that exciton population will continue to increase as the system approaches the size of a bulk material sample, potentially leading to macroscopic exciton condensation. To determine if the exciton condensate character is corresponding to the desired interlayer exciton, we visualize exciton density by calculating the probabilistic location of an electron corresponding to a specified hole location (see Methods for more information). The exciton density corresponding to λ_G is plotted in Fig. 2. Hole density is constrained to a tungsten 5d_x²−y² orbital (in red) while the probabilistic location of the electron is plotted in purple.

Fig. 2 Visualization of exciton density, where the red atomic orbital indicates the constrained location of the hole and the purple indicates the probabilistic density of the electron corresponding to that particular hole location.

Plotting the exciton density shows both inter- and intra-layer contributions to the exciton mode corresponding to the large eigenvalue. As the interlayer exciton is desirable due to being less prone to quick recombination, the presence of an interlayer exciton state supports the potential for practical realization of exciton condensation in this particular transition metal dichalcogenide bilayer. These results align with results from periodic electronic structure calculations with the Bethe–Salpeter equation, a common method for studying excitons in materials. Calculations of the excitonic spectra with the Bethe–Salpeter equation reveal strong interlayer excitons, with the lowest energy excitonic state composed of holes located on the WSe₂ layer and electrons located on the MoSe₂ layer.^36–40 Excitons with intralayer character are also present. Plots of the excitonic wavefunctions reveal highly delocalized excitons,^36,39 similar to our observation of delocalized exciton density.

As the molecular models are passivated with hydrogens, it is important to determine if the added hydrogens have a significant impact on the excitonic states we are investigating. The probabilistic density corresponding to the excitonic ground state (and thus the condensate) shows essentially no density residing on the hydrogens, demonstrating that the hydrogen atoms contribute very minimally, if at all, to the exciton state containing the condensate character and thus their inclusion does not have significant impact on the results presented here.

Interlayer distance in MoSe₂–WSe₂

Experimental studies of exciton condensation in transition metal dichalcogenide bilayers primarily utilize devices constructed with two transition metal dichalcogenide monolayers separated by one or more layers of hexagonal boron nitride (hBN).^8,13,41 The hBN layers are included as a tunnel barrier to prevent exciton recombination. Previous theoretical work has suggested that interlayer distance can have a large effect on the extent of exciton condensation in bilayer systems,^17,18 suggesting that varying the thickness of the tunnel barrier could serve as a way to tune exciton condensation in transition metal dichalcogenide bilayers.

To explore the influence of varying interlayer distance between the transition metal dichalcogenide layers on the extent of exciton condensation, λ_G is calculated for MoSe₂–WSe₂ (2 × 2) with interlayer distances ranging from 2.745 to 5.745 Å. This distance range was chosen to be in the realm of those achievable with hBN layers, although computational expense prohibits the inclusion of actual hBN layers.

These results, plotted in Fig. 3e, show the signature of exciton condensation maximized around 3.745 Å (the layer spacing of the MoSe₂–WSe₂ model used in the calculations above). Directly above the data point for each interlayer distance is plotted exciton density of the exciton mode corresponding to that particular large eigenvalue (Fig. 3a–d). These plotted exciton densities show interlayer exciton character and electron delocalization maximized around 3.745 Å, the interlayer distance which also maximized λ_G. As interlayer distance increases to 4.745 Å, an interlayer exciton is barely present, and the electron density is highly localized on a selenium p orbital. At 5.745 Å, all interlayer character has vanished and the exciton is localized on adjacent tungsten d orbitals. When interlayer distance is decreased to 2.745 Å, all interlayer exciton character has once again vanished and the exciton is localized on adjacent tungsten d orbitals as for the 5.745 Å case. These calculations appear to demonstrate that interlayer distance impacts the formation of the interlayer exciton desirable for long exciton lifetimes, which is a necessary factor for the formation of an exciton condensate. It should be noted that these calculations do not include the actual hBN spacer that would be used in experimental devices, and so the exact distance trends observed here may not be consistent to those seen experimentally. We believe, however, that the general trend observed in our calculations will hold in the experimental case—that there is a range of ideal layer separations, and distances above and below that range will result in a decrease in observed exciton condensation.

Fig. 3 Visualization of the exciton density of the exciton mode corresponding to the largest eigenvalue of the particle-hole density matrix (λ_G) (a–d) for different interlayer spacings of MoSe₂–WSe₂. The largest eigenvalue corresponding to each density is plotted below each image against the distance between MoSe₂ and WSe₂ layers (e). The atomic orbital corresponding to the constrained location of the hole is plotted in red, while the probabilistic density of the electron is plotted in purple. The density plot and data point for the original interlayer distance are highlighted in grey.

Layer alignment in MoSe₂–WSe₂

The 2D nature of the heterostructure bilayers means that electronic properties can be engineered via different interlayer stacking alignments.^42–46

Due to the broken inversion symmetry in TMD bilayers, two distinct stacking orders emerge: a parallel configuration termed R (rhombohedral, 0° alignment between the layers) and an antiparallel configuration termed H (hexagonal, 180° alignment between the layers). Each configuration gives three high-symmetry stacking orders. To assess the influence of interlayer alignment on the degree of exciton condensation, λ_G is calculated for 2 × 2 models based on all six high-symmetry stacking orders of MoSe₂–WSe₂. These six stacking orders are referred to via the terminology used in ref. 47, wherein the superscript and subscript indicate to the top and bottom layers respectively, and h, X, and M refer to the hexagonal center, the transition metal atom, and the chalcogen atom. Each model is constructed by rotating the MoSe₂ layer of MoSe₂–WSe₂ such that the layer alignment corresponds to one of the six high symmetry stacking orders of MoSe₂–WSe₂. The signature of exciton condensation is plotted for each of the six interlayer alignments in Fig. 4g, with the exciton density of the mode corresponding to the large eigenvalue plotted next to each point. Note that here the location of the electron (purple) is constrained while the probabilistic density of the hole (red) is plotted. The location of the constrained particle differs between stackings, as the atomic-orbital composition of the active molecular orbitals involved in the exciton differs between stackings. A schematic of each alignment is depicted in Fig. 4a–f.

Fig. 4 Schematics of six high-symmetry stacking orders of MoSe₂–WSe₂ (a–f). Passivating hydrogens omitted for clarity. (g) The largest eigenvalue of the particle-hole reduced density matrix (λ_G) is plotted for all high symmetry stackings, with R-type stackings plotted with blue circles and H-type stakings plotted with red triangles for clarity. Visualizations of the exciton density of the exciton mode corresponding to λ_G are provided next to each data point. The atomic orbital corresponding to the constrained location of the particle is plotted in purple, while the probabilistic density of the hole is plotted in red. The schematic and data point corresponding to the original layer alignment (H^h_h) are highlighted in grey.

Calculating exciton population for each layer alignment shows that exciton condensate character is maximized for the R^h_h stacking, and minimized for the R^X_h and H^M_h stackings although the eigenvalues of all stackings other than R^h_h are quite similar, falling in the range of 1.23–1.29. These differences in λ_G are likely due to differing interlayer interactions between the different stacking orders. Differing layer alignments lead to differing electrostatic interactions between atomic sites,^42,48 which may impact exciton formation. Differing orbital overlap may also drive the differences in λ_G observed-aligning with previous experimental and theoretical studies which have highlighted the impact of interlayer twist angle, and thus interlayer coupling and hybridization, on exciton formation in TMD van der Waals heterostructures.^49–53 It is notable that the R^h_h stacking, which corresponds to each metal atom aligned directly over a metal atom on the opposite layer and each chalcogen atom directly over a chalcogen atom on the opposite layer, results in greater exciton condensation relative to the other stackings. This suggests that maximizing the alignment between sites on opposing layers may be significant for the emergence of large λ_G.

By comparing the exciton density of the mode corresponding to the large eigenvalue of the different high symmetry stackings, it is apparent that while each stacking has highly delocalized intralayer exciton density, only the R^h_h, H^h_h, R^M_h, and H^X_h stackings have significant interlayer exciton densities. Amongst these, the high-λ_GR^h_h stacking has the most interlayer character, with the density of the hole delocalized among all four Tungsten atoms. The H^M_h and R^M_h stackings, which have no major interlayer exciton density, also have the lowest values of λ_G. These results suggest that the differing interactions between layers in the various stackings make some stackings more favorable for the formation of stable, highly delocalized interlayer excitons, which in turn supports the emergence of exciton condensation.

While it is not clear if this exact trend for the six high symmetry stackings persists in bulk samples of this material, we expect that layer alignment will have a significant impact on exciton formation and thus the ability of the material to form an exciton condensate.

Exciton condensation in different bilayers

It is of interest to see if other TMD bilayers show potential for the exciton condensation observed in MoSe₂–WSe₂. To this end, a family of related homo and heterobilayers composed of different combinations of MoSe₂, WSe₂, MoS₂, and WS₂ were chosen for investigation. To probe the chosen transition metal dichalcogenide bilayers for the beginnings of exciton condensation, λ_G is calculated for both the 2 × 2 and 3 × 3 models of each bilayer. These results are plotted in Fig. 5h. The exciton density of the exciton mode corresponding to λ_G is plotted for each 3 × 3 model.

Fig. 5 (a–g) Visualizations of the exciton density of the exciton mode corresponding to λ_G for the 3 × 3 model of each heterostructure. The atomic orbital corresponding to the constrained location of the hole is plotted in red, while the probabilistic density of the electron is plotted in purple. The largest eigenvalue of the particle-hole density matrix (λ_G) for both the 2 × 2 and 3 × 3 cluster models of several related homo and heterobilayers. Schematics of a top view of the 2 × 2 and 3 × 3 models are provided for reference (passivating hydrogens omitted for clarity).

As can be seen in Fig. 5, all homobilayer and heterobilayer systems exhibit λ_G > 1 and thus demonstrate the potential for exciton condensation. The relative magnitudes of λ_G between the systems are not necessarily indicative of their relative capabilities of exciton condensation, as the above results have shown that structural factors such as layer spacing have considerable impact on the beginnings of exciton condensation. While all models have interlayer distances between 3.0 and 4.0 Å, close to the ideal layer separation for MoSe₂–WSe₂ determined previously, the ideal layer separation for each of these new systems is likely different from that for MoSe₂–WSe₂. However, the plotted exciton densities in Fig. 5a–g give some insight into the trends observed in Fig. 5h. MoSe₂–WS₂ displays a notably higher λ_G than the other bilayers, which, as in the other cases discussed in the work, is accompanied by an increase in interlayer exciton density. All of the bilayers show highly delocalized exciton density with some interlayer component, which is a necessary factor for the emergence of exciton condensation. The fact that λ_G > 1 is observed for each new bilayer system, even in very small fragments of the bulk systems, speaks promisingly to the potential of these materials for exciton condensation.

Conclusions

Here we observe the microscopic origins of exciton condensation in a family of transition metal dichalcogenide bilayers. We find that the signature of exciton condensation is present in a fragment of MoSe₂–WSe₂, aligning with the experimental evidence for exciton condensation in this material. Visualizing the exciton density of the exciton mode corresponding to the largest eigenvalue shows that the probabilistic electron density (with respect to a hole constrained to a tungsten 5d_x²−y² orbital) is highly delocalized, with both intra and interlayer components. The presence of interlayer exciton density aligns with the experimental results, which specifically target the formation of an interlayer exciton condensate. The degree of condensation increases for a larger-sized molecular fragment, indicating the presence of a “critical seed” of exciton condensation that may become a macroscopic condensate in a bulk system. We also find that layer alignment and interlayer distance influence the extent of exciton condensation. Out of six high-symmetry alignments of MoSe₂–WSe₂, the R^h_h alignment shows the largest signature of exciton condensation, likely due to favorable electrostatic interactions and orbital overlap that aid the formation of robust interlayer excitons. Varying interlayer distance between 2.745 and 5.745 Å shows that maximal exciton condensation occurs at an intermediate distance of 3.745 Å, while bringing the layers too close together or too far apart diminishes the extent of exciton condensation. Finally, we find that the signature of exciton condensation is also present in several transition metal dichalcogenide bilayers similar to MoSe₂–WSe₂, suggesting the potential for exciton condensation. While the magnitudes of λ_G differ between the systems, this is not necessarily indicative of their relative capabilities of exciton condensation, as ideal layer alignment and interlayer separation was not taken into consideration. However, as calculated the MoSe₂–WS₂ structure produces a robust interlayer exciton with a notably large λ_G.

While it has been understood that the large eigenvalue signature is a necessary condition for long-range order and thus the beginnings of exciton condensation,^15,35 previous work has primarily focused on novel molecular systems that are not direct analogues to those explored experimentally.^16–18 In this work, by theoretically probing a molecular fragment of a material, MoSe₂–WSe₂, that has already shown experimental evidence for exciton condensation, we are able to draw a direct connection between the necessary condition for the beginning of exciton condensation and actual experimental observation of condensation. While the large eigenvalues we see in our molecular-scale fragments do not necessarily suggest condensation in the macroscopic sense, we believe they show us the beginnings of exciton condensation which may grow into macroscopic exciton condensation as the system grows to approach the thermodynamic limit. This idea is supported by the fact that λ_G increases as we increase system size from the 2 × 2 to the 3 × 3 model. Our results also give us unique insight into the molecular origins of the phenomena that is observed in material-scale samples, as the presence of a “critical seed” of condensation in a very small fragment of the material, which can be influenced by geometric modifications of the molecular models, suggests that local short-range strongly correlated effects may be significant for the formation of a macroscopic condensate. The fact that layer alignment and interlayer distance impact the magnitude of the signature of exciton condensation also highlights the importance of geometric considerations in the design of exciton condensate materials. We also observe the signature of exciton condensation in several other transition metal dichalcogenide bilayers in addition to MoSe₂–WSe₂, suggesting that the selection of transition metal dichalcogenide involved may not have a large impact on the potential for exciton condensation. This observation may imply that materials capable of exciton condensation are more ubiquitous than previously thought, and motivates further investigation into new bilayer systems for realizing exciton condensation.

Methods

Theory

A definitive theoretical signature for exciton condensation has recently been developed from reduced density matrix (RDM) theory: a large eigenvalue in the particle-hole reduced density matrix

²G_k,l^i,j = 〈Ψ|â^†_iâ_jâ^†_lâ_k|Ψ〉 (1)

where |Ψ〉 is the N-electron wavefunction and â^†_i and â_i are the fermionic creation and annihilation operators.^15,35,54 The particle-hole RDM can be obtained from the two-electron reduced density matrix (2-RDM) by a linear mapping given by:

²G_k,l^i,j = ¹I_l^j1D_kⁱ − ²D_k,j^i,l (2)

where ²D_k,j^i,l is the 2-RDM, ¹I_l^j is the one-electron identity matrix, and ¹D_kⁱ is the one-electron reduced density matrix (1-RDM). Additionally, one must subtract from the particle-hole RDM the components corresponding to a ground-state-to-ground-state transition, which otherwise creates one extraneous large eigenvalue unrelated to exciton condensation:

²G̃_k,l^i,j = ²G_k,l^i,j − ¹D_jⁱ¹D_k^l. (3)

The largest eigenvalue of this modified particle-hole RDM is the signature of exciton condensation, and is referred to as λ_G. The eigenvalues of the modified particle-hole reduced density matrix correspond to the occupation of each particle-hole state, and so the magnitude of λ_G describes the extent of condensation.

Visualization of exciton density

Visualization of the excitonic state corresponding to λ_G is obtained by constraining the location of the excitonic hole to a particular atomic orbital and plotting the probabilistic density of the electron. A matrix of molecular orbitals in terms of atomic orbitals (M_MO,AO) is obtained from the output of v2RDM calculation performed in the Quantum Chemistry Toolbox for Maple.^55,56 This matrix is then used to calculate a matrix of molecular orbitals in terms of atomic orbitals, M_AO,MO

M_AO,MO = (M^T_MO,AO)⁻¹. (4)

A submatrix M^active_AO,MO corresponding to the active orbitals is isolated from M_AO,MO, and the eigenvector or the modified particle-hole RDM corresponding to λ_G is reshaped as a matrix in the basis of the active molecular orbitals. This matrix, called V_max, is multiplied along with the submatrix

(M^active_AO,MO)(V_max)(M^active_AO,MO)^T (5)

to create a matrix representing hole atomic orbitals in terms of the coefficients of the contributions of the electrons to other molecular orbitals.

Molecular models

Several different transition metal dichalcogenide bilayers are chosen for investigation. These include homobilayers composed of two layers of MoSe₂, WSe₂, MoS₂, and WS₂ respectively, as well as the heterobilayers MoSe₂–WSe₂, MoS₂–WS₂, and MoSe₂–WS₂. In order to make models that can be explored via current electronic structure techniques for strong correlation at reasonable computational cost, cluster models are built from bulk experimental crystal structures retrieved from the Inorganic Crystal Structure Database,^57–61 with dangling bonds saturated by hydrogens. Two cluster models are constructed for each system, consisting of bilayers with the primitive cell tiled 2 × 2 and 3 × 3 in the XY plane respectively. As the primitive cells of the models do not contain a C3 rotational axis, the molecular models do not exhibit the C3 rotational symmetry that may be found in the bulk material. The layers of all heterobilayer structures were spaced arbitrarily at 3.474 Å apart. Molecular coordinates for all structures are available in the ESI.†

The 2-RDM is obtained directly from the molecular structure via variational two-electron reduced density matrix theory (v2RDM).^20–34 The v2RDM calculations are performed with the Quantum Chemistry Toolbox for Maple.^55,56 We specifically utilize v2RDM-CASSCF, a combination of v2RDM with the complete active space self-consistent-field (CASSCF) method which allows for optimization of orbitals in the chosen active space.²⁴ An [8,8] active space is used for all 2 × 2 structures, while an [18,18] active space is used for all 3 × 3 structures, where [X,Y] indicates an active space of X electrons in Y orbitals. These active spaces correspond to an active electron in one orbital on each of the metal atoms for both structures. The def2-SVP basis set is used for all calculations.^62,63 The particle-hole RDM is then constructed from the 2-RDM via the method described above. The eigenvalues of the particle-hole RDM are then calculated, with the largest of which referred to as λ_G.

Data availability

Data will be provided by the corresponding author upon reasonable request.

Author contributions

DM conceived of the project. LP performed theoretical and computational calculations. LP, AS, and DM interpreted the results and data and wrote the manuscript.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

DAM gratefully acknowledges support from the U. S. National Science Foundation Grant No. CHE-2155082, the Department of Energy, Office of Basic Energy Sciences, Grant DE-SC0019215, and the ACS Petroleum Research Fund Grant No. PRF No. 61644-ND6.

Notes and references

1. D. V. Fil and S. I. Shevchenko, Low Temp. Phys., 2018, 44, 867–909.
2. L. V. Keldysh, Phys.-Usp., 2017, 60, 1180–1186.
3. M. S. Fuhrer and A. R. Hamilton, Physics, 2016, 9, 80.
4. J. Kasprzak, M. Richard, S. Kundermann, A. Baas, P. Jeambrun, J. M. J. Keeling, F. M. Marchetti, M. H. Szymańska, R. André, J. L. Staehli, V. Savona, P. B. Littlewood, B. Deveaud and L. S. Dang, Nature, 2006, 443, 409–414.
5. H. Min, R. Bistritzer, J.-J. Su and A. H. MacDonald, Phys. Rev. B, 2008, 78, 121401.
6. A. Kogar, M. S. Rak, S. Vig, A. A. Husain, F. Flicker, Y. I. Joe, L. Venema, G. J. MacDougall, T. C. Chiang, E. Fradkin, J. Van Wezel and P. Abbamonte, Science, 2017, 358, 1314–1317.
7. J. I. A. Li, T. Taniguchi, K. Watanabe, J. Hone and C. R. Dean, Nat. Phys., 2017, 13, 751–755.
8. Z. Wang, D. A. Rhodes, K. Watanabe, T. Taniguchi, J. C. Hone, J. Shan and K. F. Mak, Nature, 2019, 574, 76–80.
9. X. Hu, T. Hyart, D. I. Pikulin and E. Rossi, Phys. Rev. B: Condens. Matter Mater. Phys., 2022, 105, L140506.
10. X. Liu, K. Watanabe, T. Taniguchi, B. I. Halperin and P. Kim, Nat. Phys., 2017, 13, 746–750.
11. J.-J. Su and A. H. MacDonald, Phys. Rev. B: Condens. Matter Mater. Phys., 2017, 95, 045416.
12. L. Sigl, F. Sigger, F. Kronowetter, J. Kiemle, J. Klein, K. Watanabe, T. Taniguchi, J. J. Finley, U. Wurstbauer and A. W. Holleitner, Phys. Rev. Res., 2020, 2, 042044.
13. L. Ma, P. X. Nguyen, Z. Wang, Y. Zeng, K. Watanabe, T. Taniguchi, A. H. MacDonald, K. F. Mak and J. Shan, Nature, 2021, 598, 585–589.
14. M. M. Fogler, L. V. Butov and K. S. Novoselov, Nat. Commun., 2014, 5, 4555.
15. S. Safaei and D. A. Mazziotti, Phys. Rev. B, 2018, 98, 045122.
16. A. O. Schouten, L. M. Sager and D. A. Mazziotti, J. Phys. Chem. Lett., 2021, 12, 9906–9911.
17. L. M. Sager, A. O. Schouten and D. A. Mazziotti, J. Chem. Phys., 2022, 156, 154702.
18. A. O. Schouten, J. E. Klevens, L. M. Sager-Smith, J. Xie, J. S. Anderson and D. A. Mazziotti, Phys. Rev. Mater., 2023, 7, 045001.
19. L. M. Sager, S. E. Smart and D. A. Mazziotti, Phys. Rev. Res., 2020, 2, 043205.
20. Reduced-density-matrix mechanics: with applications to many-electron atoms and molecules, ed. D. A. Mazziotti, Wiley-Interscience, Hoboken, NJ, 2007.
21. D. A. Mazziotti, Phys. Rev. A, 1998, 57, 4219–4234.
22. M. Nakata, H. Nakatsuji, M. Ehara, M. Fukuda, K. Nakata and K. Fujisawa, J. Chem. Phys., 2001, 114, 8282–8292.
23. D. A. Mazziotti, Phys. Rev. Lett., 2004, 93, 213001.
24. G. Gidofalvi and D. A. Mazziotti, J. Chem. Phys., 2008, 129, 134108.
25. D. A. Mazziotti, Chem. Rev., 2012, 112, 244–262.
26. D. A. Mazziotti, Phys. Rev. A, 2016, 94, 032516.
27. D. A. Mazziotti, Phys. Rev. Lett., 2011, 106, 083001.
28. Many-Electron Densities and Reduced Density Matrices, ed. J. Cioslowski and P. G. Mezey, Springer US, Boston, MA, 2000.
29. Z. Zhao, B. J. Braams, M. Fukuda, M. L. Overton and J. K. Percus, J. Chem. Phys., 2004, 120, 2095–2104.
30. E. Cancès, G. Stoltz and M. Lewin, J. Chem. Phys., 2006, 125, 064101.
31. N. Shenvi and A. F. Izmaylov, Phys. Rev. Lett., 2010, 105, 213003.
32. D. A. Mazziotti, Phys. Rev. Lett., 2012, 108, 263002.
33. B. Verstichel, H. Van Aggelen, W. Poelmans and D. Van Neck, Phys. Rev. Lett., 2012, 108, 213001.
34. M. Piris, Phys. Rev. Lett., 2017, 119, 063002.
35. C. Garrod and M. Rosina, J. Math. Phys., 1969, 10, 1855–1861.
36. E. Torun, H. P. C. Miranda, A. Molina-Sánchez and L. Wirtz, Phys. Rev. B, 2018, 97, 245427.
37. S. Gao, L. Yang and C. D. Spataru, Nano Lett., 2017, 17, 7809–7813.
38. K. Sadecka, Acta Phys. Pol., A, 2022, 141, 110–115.
39. R. Gillen and J. Maultzsch, Phys. Rev. B: Condens. Matter Mater. Phys., 2018, 97, 165306.
40. L. L. Li, R. Gillen, M. Palummo, M. V. Milošević and F. M. Peeters, Appl. Phys. Lett., 2023, 123, 033102.
41. H. Kalt and C. F. Klingshirn, Semiconductor Optics 2, Springer International Publishing, Cham, 2024, pp. 421–465.
42. J. He, K. Hummer and C. Franchini, Phys. Rev. B, 2014, 89, 075409.
43. T. Jiang, H. Liu, D. Huang, S. Zhang, Y. Li, X. Gong, Y.-R. Shen, W.-T. Liu and S. Wu, Nat. Nanotechnol., 2014, 9, 825–829.
44. S. Wang, X. Cui, C. Jian, H. Cheng, M. Niu, J. Yu, J. Yan and W. Huang, Adv. Mater., 2021, 33, 2005735.
45. X. Hu, L. Kou and L. Sun, Sci. Rep., 2016, 6, 31122.
46. Y. Wang, W. Wei, B. Huang and Y. Dai, J. Phys.: Condens. Matter, 2019, 31, 125003.
47. S. Zhao, Z. Li, X. Huang, A. Rupp, J. Göser, I. A. Vovk, S. Y. Kruchinin, K. Watanabe, T. Taniguchi, I. Bilgin, A. S. Baimuratov and A. Högele, Nat. Nanotechnol., 2023, 18, 572–579.
48. N. Marom, J. Bernstein, J. Garel, A. Tkatchenko, E. Joselevich, L. Kronik and O. Hod, Phys. Rev. Lett., 2010, 105, 046801.
49. Y. Shimazaki, I. Schwartz, K. Watanabe, T. Taniguchi, M. Kroner and A. Imamoğlu, Nature, 2020, 580, 472–477.
50. G.-B. Liu, D. Xiao, Y. Yao, X. Xu and W. Yao, Chem. Soc. Rev., 2015, 44, 2643–2663.
51. M. Förg, A. S. Baimuratov, S. Y. Kruchinin, I. A. Vovk, J. Scherzer, J. Förste, V. Funk, K. Watanabe, T. Taniguchi and A. Högele, Nat. Commun., 2021, 12, 1656.
52. Y. Jiang, S. Chen, W. Zheng, B. Zheng and A. Pan, Light: Sci. Appl., 2021, 10, 72.
53. L. Zhang, Z. Zhang, F. Wu, D. Wang, R. Gogna, S. Hou, K. Watanabe, T. Taniguchi, K. Kulkarni, T. Kuo, S. R. Forrest and H. Deng, Nat. Commun., 2020, 11, 5888.
54. W. Kohn and D. Sherrington, Rev. Mod. Phys., 1970, 42, 1–11.
55. Maplesoft, Maple, Maplesoft, Waterloo, Ontario, 2024.
56. RDMChem, Quantum Chemistry Toolbox for Maple, Maplesoft, Waterloo, Ontario, 2024.
57. G. Bergerhoff, R. Hundt, R. Sievers and I. D. Brown, J. Chem. Inf. Comput. Sci., 1983, 23, 66–69.
58. J. Champion, Br. J. Appl. Phys., 1965, 16, 1035–1037.
59. V. Kalikhman, E. Gladchenko and L. Pravoverova, Inorg. Mater., 1972, 8, 1020–1022.
60. B. Schönfeld, J. J. Huang and S. C. Moss, Acta Crystallogr., Sect. B: Struct. Sci., 1983, 39, 404–407.
61. A. E. Van Arkel, Recl. Trav. Chim. Pays-Bas, 1926, 45, 437–444.
62. F. Weigend and R. Ahlrichs, Phys. Chem. Chem. Phys., 2005, 7, 3297.
63. K. Eichkorn, F. Weigend, O. Treutler and R. Ahlrichs, Theor. Chim. Acta, 1997, 97, 119–124.

---

Footnote

† Electronic supplementary information (ESI) available: XYZ coordinates of the 7 bilayers. See DOI: https://doi.org/10.1039/d4sc04149f

---