Wigner crystals in two-dimensional transition-metal dichalcogenides: Spin physics and readout

J. Knörzer ,1,2,* M. J. A. Schuetz,3,† G. Giedke ,4,5 D. S. Wild,3 K. De Greve,3,‡ R. Schmidt ,1,2 M. D. Lukin,3 and J. I. Cirac1,2 1Max-Planck-Institut für Quantenoptik, Hans-Kopfermann-Str. 1, D-85748 Garching, Germany 2Munich Center for Quantum Science and Technology (MCQST), Schellingstr. 4, D-80799 München, Germany 3Physics Department, Harvard University, Cambridge, Massachusetts 02318, USA 4Donostia International Physics Center, Paseo Manuel de Lardizabal 4, E-20018 San Sebastián, Spain 5Ikerbasque Foundation for Science, Maria Diaz de Haro 3, E-48013 Bilbao, Spain

Wigner crystals are prime candidates for the realization of regular electron lattices under minimal requirements on external control and electronics. However, several technical challenges have prevented their detailed experimental investigation and applications to date. We propose an implementation of two-dimensional electron lattices for quantum simulation of Ising spin systems based on self-assembled Wigner crystals in transition-metal dichalcogenides. We show that these semiconductors allow for minimally invasive all-optical detection schemes of charge ordering and total spin. For incident light with optimally chosen beam parameters and polarization, we predict a strong dependence of the transmitted and reflected signals on the underlying lattice periodicity, thus revealing the charge order inherent in Wigner crystals. At the same time, the selection rules in transition-metal dichalcogenides provide direct access to the spin degree of freedom via Faraday rotation measurements. DOI: 10.1103/PhysRevB.101.125101

I. INTRODUCTION
Ever since its theoretical inception 85 years ago [1], Wigner crystallization has stimulated both theoretical and experimental research to find unambiguous evidence for this elusive state of matter. Since the earliest indication for quantum Wigner crystals (WCs) obtained from high-magneticfield transport measurements [2,3], it has proven to be a very demanding task to study WCs, especially in a minimally invasive manner without destroying the crystalline order. Recent experimental work demonstrated nondestructive read-out of the charge distribution of one-dimensional WCs in carbon nanotubes [4]. However, it remains an open challenge to find approaches for the noninvasive detection of WCs in twodimensional and broader ranges of one-dimensional quantum systems.
Apart from a fundamental interest in the physics of Wigner crystallization, self-assembled crystals promise a route toward highly ordered and scalable many-body systems under minimal external control. Thus, they meet some of the key requirements posed by quantum computers [5] and simulators [6]. It has therefore been proposed that WCs hosted in semiconductor nanostructures [7,8], trapped above the surface of liquid helium [9,10] or composed of trapped ions [11,12], can be utilized for quantum information processing and simulation. In particular, electrons confined to low-dimensional semiconductors [13] may be brought into the low-temperature regime k B T ε F (Fermi energy ε F ) where quantum phenomena occur and spin-exchange interactions can play an important role. Since solid-state systems also offer a genuine prospect for miniaturization and on-chip integration, the quest for a faithful implementation of solid-state quantum WCs at zero magnetic field remains tantalizing.
In this paper, we demonstrate the potential of scalable quantum simulators based on two-dimensional WCs in TMDs and propose an all-optical detection scheme for charge ordering and partial spin information in these systems (see Fig. 1). In particular, the scheme possesses three key properties: (i) It provides clear evidence for Wigner crystallization in monolayer TMDs. (ii) Under conditions specified below, the detection scheme is noninvasive and leaves charge and spin order intact. (iii) Optical selection rules provide spin-selective addressability, which is a crucial requirement for quantum simulation.
FIG. 1. Schematic illustration of proposed setup and optical detection scheme. Charge ordering of electrons in a lattice (black dots) competes with random disorder-induced dislocations of lattice sites in the presence of impurities and defects (green triangle). The angle-dependent (φ) reflection of a tilted (θ) focused laser beam with wave vector k from a WC probes its lattice geometry. Light polarization provides further information about the spin via optical selection rules of TMDs.

A. Wigner crystals
At electron densities n below a critical density n cr and in the presence of an external confinement potential, interacting charge carriers (referred to as electrons in the following) arrange themselves in a lattice [25], leading to a periodic modulation of charge density n(r). In this low-density regime, electrostatic interactions dominate over the kinetic energy of electrons. In two dimensions, this regime is characterized by a sufficiently large interaction parameter r s = 1/( √ π na B ), with the Bohr radius a B = 4πεh 2 /(e 2 m), effective electron mass m, and permittivity ε. Monolayer TMDs feature an extraordinarily small Bohr radius a B 0.5 nm and thus render the large-r s regime accessible at experimentally achievable [26,27] densities n n cr . For our calculations, we choose n cr = 10 11 cm −2 [14] and m = 0.5m 0 (representative of MoX 2 monolayers where X = S, Se [28]), where m 0 denotes the bare electron mass. In a square (triangular) lattice, this maximum electron density corresponds to a minimum lattice spacing of a 32 nm(a 34 nm).

B. Model
We consider N electrons trapped at z = 0 in a global harmonic potential such that the total potential reads where r i = (x i , y i , 0) denotes the position of the ith electron. The confinement is characterized by the trapping frequency ω and V int denotes the two-body interaction potential. In TMDs, the former may be induced by strain [29,30] or defined via local gates [31] and the latter is usually modeled by the Keldysh potential [32], with a material-specific length scale r 0 ≈ 5 nm. H 0 and Y 0 are Struve and Bessel functions, respectively. At electron concentrations n < n cr , the interparticle distance |r i − r j | r 0 and hence V int (r i , r j ) ∼ 1/|r i − r j | behaves like a Coulomb potential.
In a WC, the electrons are localized around lattice sites at r 0 i (i = 1, . . . , N) which can be determined from the equilibrium conditions ∇ i V | r i =r 0 i = 0. Numerical calculations show that harmonic confinement potentials, as described in Eq. (1), give rise to triangular lattice geometries while other potentials can give rise to, e.g., square lattices; see Appendix A for details. For any ω, the maximum number of WC electrons can be calculated given a critical density, and vice versa. Small systems containing N ∼ (10 − 100) electrons requirē hω ∼ (1 − 3) meV at n ∼ n cr (see Fig. 2).
The strong interactions in Eq. (2) enable the description of charge excitations in terms of phonons in the WC. These can be expressed as small displacements q i = r i − r 0 i (i = 1, . . . , N) from the lattice sites such that V = (m/2) K αβ i j q α i q β j (α, β ∈ {x, y}) with an elasticity matrix K. All 2N normal modes of the system with eigenfrequencies n (n = 1, . . . , 2N) are readily obtained by diagonalization of K and for the nonzero eigenfrequencies one finds that n ω (cf. Appendix A). Given the relation between ω and N at n ∼ n cr , this indicates that large WCs have low-energy phonon modes. Using anharmonic potentials, there is no limit placed on N by the phonon modes or n cr .

C. Requirements
Wigner crystallization requires low disorder. Disorderinduced potential fluctuations are incorporated based on Eq. (1) by adding further randomly distributed local confinement terms to analyze the impact of impurities (e.g., atomic defects or charges) on the electron lattice. To obtain a regular lattice structure with an approximately equidistant spacing between adjacent electrons (see schematic Fig. 1), the impurity density n imp should be significantly smaller than the electron density, i.e., n imp 0.1n; see Appendix B for details.
To date, atomic and charge defects in TMDs still prevent the realization of systems with sufficiently low disorder [33]. However, both sample quality and deterministic control over defects [34] have been improving rapidly in recent years and defect densities around n cr can already be achieved. Moreover, WCs require sufficiently low temperature. Cooling into the motional ground state requires low temperatures T ∼ 1 K forhω meV, as the thermal occupationn th = 1/[exp(h n /(k B T )) − 1] of the modes increases as ω is decreased (cf. Appendix C).
There are many interesting aspects about the dynamics of strongly correlated electrons that can be studied in the system we describe, including the entanglement properties of the ground state, the nature and dynamics of excitations, and the transitions to neighboring phases. In the following, we focus on the spin physics.

III. SPIN PHYSICS
TMD monolayers exhibit strong spin-orbit coupling and an intricate interplay between spin and valley degrees of freedom. Here we focus on the case where, by energetic isolation of the lower spin states of the conduction band, spin and valley become locked [35]. For this reason, we require that the electron density be sufficiently low such that the Coulomb interaction energy E int ∼ r s · ε F = r s πh 2 n/m is small compared to the spin-orbit splitting in the conduction band, c SO . At n n cr , one typically finds E int 10 meV, such that the above condition is readily satisfied in MoSe 2 ( c SO ≈ 23 meV), though not necessarily in MoS 2 ( c SO ≈ 3 meV) [36]. Nevertheless, the requirement can be met in all TMDs by considering holes instead of electrons, since the spin-orbit splitting in the valence band v SO is on the order of a hundred meV [21].
At low temperature and small displacements q i , we assume that the electron spins are localized around the lattice sites at r 0 i . Adjacent spins are coupled via exchange interactions that can be either ferromagnetic or antiferromagnetic, depending on the density n [37,38]. Here, we provide an estimate for the magnitude of the spin-spin coupling, demonstrating the potential of TMD-based electron lattices as a platform for quantum simulation of prototypical spin systems. As exchange couplings decay exponentially with a 2 , where a denotes the interparticle distance, the low-density regime necessary for WCs stands in contrast with the strong couplings of interest for spin physics. However, at intermediate densities n n cr , we still find significant exchange couplings which exceed predicted spin relaxation rates [39,40]. Due to the spin polarization in each of the K and K valleys, we find that the effective spin model in the spin-valley locked, low-temperature regime reduces to an Ising Hamiltonian (cf. Appendix D for details) of the form Here σ z i is a Pauli operator and J i j denotes the coupling strength between spins at sites i and j. In a tight-binding approximation, we calculate J i j (1 i < j N) using Gaussian ansatz wave functions centered around the sites r 0 i . The width of these wave functions is expressed in terms of the normal mode frequencies n and, upon inserting typical material parameters, we find for the magnitude J of the spin-spin interaction between nearest neighbors typical values in the range J ∼ (5 − 30) μeV for n n cr . Due to the exponential decay of J i j with distance, nearest-neighbor interactions are dominant and typically roughly one order of magnitude larger than next-nearest-neighbor interactions. In Fig. 2(b), we show the resulting spin-coupling constant J as a function of ω for different particle numbers 10 N 100. At the intermediate densities n n cr considered here, we find antiferromagnetic exchange couplings which can result in geometrical frustration [41], depending on the lattice structure.

IV. OPTICAL READOUT
We now address the optical detection of charge ordering in TMD-based WCs and consider an incoming (z < 0) Gaussian laser beam E in (r) with wavelength λ focused to a spot on the electron lattice (z = 0) at a tilt angle θ (see Fig. 1). Our approach is similar in nature to the one taken in Refs. [42,43], where the reflection and transmission of arrays of discrete atomic emitters in a lattice configuration were analyzed. Such an approach is valid for highly localized charges [44], in contrast to the study of mobile polarons [45]. Due to optical transition selection rules in monolayer TMDs, specific electron spin states can be addressed using circularly polarized σ + and σ − light. For example, σ − (σ + ) light may couple a WC electron in a |↑ K (|↓ K ) spin state to a trionic state |↑ K , ↓ K ⇓ K (|↓ K , ↑ K ⇑ K ) with a hole spin ⇑ (⇓) in the K (K) valley. For our calculations, we assume a low-amplitude light beam with sufficiently small detuningh 0 E b , E g from the trion resonance such that other quasiparticle excitations and transitions can be neglected. Prototypical values for the trion binding energy E b ∼ 20 meV and quasiparticle band gap E g ∼ 500 meV are given in Ref. [21]. When the incoming beam is sufficiently close to resonance with a dipole transition at lattice points r 0 n , the scattered light field E(r) at position r is obtained by solving a set of coupled linear equations, with the detuning from resonance 0 , the dyadic Green's function G evaluated at k = 2π/λ, and the polarizability tensor α n . The magnitude of the polarizability tensor is given by the scalar polarizability α( 0 ), while the orientation depends on the electron spin at site n; see Appendix E for more details.
To probe charge ordering, it is advantageous to address all WC electrons equally. To this end, we assume for the following discussion that the WC is fully spin polarized, which could be achieved by applying a large magnetic field or via optical pumping [46]. Alternatively, one could consider a TMD heterobilayer system where an electron-hole pair excited in one layer forms a trion state with a WC electron in the other layer, such that both valleys can be addressed independent of the spin of the resident electron [21]. The total power P transmitted by the WC to z > 0 is obtained by integrating the transmitted signal (μ 0 = 1), with the electric and magnetic fields E and B, respectively, and B * denotes the complex conjugate of B. The transmission T = P wc /P 0 is calculated as a function of density n, incidence angle θ , and rotation angle φ (see Fig. 1) by comparing the transmitted power P wc in the presence of a WC with a reference signal P 0 obtained in the absence of localized dipoles [42], e.g. in a system with no doping at n = 0. In Fig. 3, T is shown as a function of the electron density n ∼ 1/a 2 for square [ Fig. 3(a)] and triangular [ Fig. 3(c)] lattices with a lattice constant a. Here we consider 0 = 0, which corresponds to a wavelength λ ∼ (700-800) nm in state-of-the-art TMD setups [47,48]. We choose θ such that the cross section of the Gaussian beam is small enough and does not exceed the size of the WC. Varying the twist angle φ of the laser beam leads to smooth variations in T (φ). The periodic modulation of T (φ) reflects the rotational symmetry of the WC.  (cf. Appendix E). Momentum transfer onto the WC can be safely neglected since the recoil energy E R =h 2 k 2 /(2m) ∼ (5 − 10)μeV is much smaller than interaction energy and trapping potential. This approach already incorporates spin information, as it can be used to detect ferromagnetic ground states and may pick up signatures of the lattice constant 2a prevailing in an antiferromagnetic ground state.

Faraday rotation
While we have focused on the detection of charge ordering in a spin-polarized WC before, we now further examine the spin degree of freedom by analyzing the polarization of the scattered field. With the probe beam E in detuned far enough from the trionic resonance, the presence of the optical transition merely imprints a state-dependent phase shift on the incoming field. According to selection rules of monolayer TMDs [49,50], σ + (σ − ) polarized light couples to the resident electron density n ↑ (n ↓ ) in the K (K ) valley [see Fig. 4(b)]. In optical Faraday (Kerr) rotation using linearly polarized light, the polarization of the transmitted (reflected) part of the light is rotated by an angle θ F which depends on the spin imbalance n ↑ − n ↓ [46,51]. Here we inspect the Faraday rotation of an incident s-or p-polarized beam, which is given by [52] where χ F = t ps /t ss for s-polarized light (χ F = −t sp /t pp for p-polarized light) depends on the Jones matrix elements t ss , t ps (t pp , t sp ) encoding the polarization state of the scattered light [53]. We consider N ↑ (N ↓ ) electrons in the |↑ K (|↓ K ) conduction band and numerically calculate θ F as a function of spin imbalance N ↑ − N ↓ and detuning 0 . Here we assume that the electron sites r 0 i are distributed in a square-lattice configuration in the spot of the beam with N ↑ (N ↓ ) randomly assigned |↑ (|↓ ) states. We average over many such configurations. In Fig. 4(a), the resulting Faraday rotation is depicted for a p-polarized input field, yielding the strongest signal at | 0 | = γ r /2 with the radiative linewidth γ r . For the strongly localized quantum emitters considered here, we estimatehγ r ∼ 10 −2 μeV. Nonradiative decay processes can also be taken into account in our framework, yielding weaker Faraday signals for larger nonradiative decay rates γ nr (cf. Appendix E). Since the Faraday rotation is proportional to N ↑ − N ↓ , it provides a measure for the spin imbalance in the system. With this tool, one may distinguish between ferromagnetic and antiferromagnetic configurations or even locally probe domain walls in the spin system, where the spatial resolution would be limited by the spot size ∼λ 2 .

V. SUMMARY AND OUTLOOK
In conclusion, we have proposed an all-optical detection scheme for TMD-based WCs, highlighting their potential as a platform for the quantum simulation of geometrically frustrated magnetism with adjustable and self-assembled lattice structures. Beyond the Ising model considered here, richer spin physics with multi-spin-exchange interactions has been predicted for these systems, potentially offering a platform to study three-and four-body interactions [54,55]. Moreover, recent results show that multielectron quantum dots hold promise as exchange-based mediators of quantum information [56]. In this context, intermediate-scale WCs in 2D semiconductors could be interesting for achieving long-range spin coupling with minimal external control requirements [8]. Control over the spin degree of freedom may be provided via magnetic fields or optical pumping into a specific valley, e.g., in parts of the system to study the formation of domain walls. Inversion symmetric TMD bilayers, whose bands are spin degenerate, may further give rise to a wider range of spin Hamiltonians and allow for coherent optical control of the electron spin as no momentum is required to flip the spin. High-quality samples of monolayer TMDs should provide access to first proof-of-principle experiments with small system sizes. Local spin probes may be enabled by illuminating only parts of the WC. Besides the optical techniques we propose, which we believe can be readily implemented given sufficiently clean samples, we envisage that it might become possible in the future to extend existing and developing work on high-resolution electron beam imaging with (close-to) single site resolution [57,58] to the point that a single electron charge can be directly spatially probed. Furthermore, other detection schemes could be considered like magnetic noise spectroscopy [59], microwave spectroscopy [60], or using surface acoustic waves in piezoelectric TMD monolayers [61].

APPENDIX A: CALCULATION OF LATTICE STRUCTURE AND NORMAL MODES
We consider a general potential of the form where μ p is the strength of the potential and the interaction potential V int is modeled by the Keldysh interaction potential given in Eq. (2). The results presented in the main text are derived for the special cases p = 2 and μ 2 = mω 2 /2.

Lattice structure
The lattice sites r 0 i are calculated by solving the equations for each electron i ∈ {1, . . . , N}. This leads to a set of 2N coupled equations which are of the form with α ∈ {x, y}, ξ = π e 2 /(2r 3 0 ) and the function which is obtained by making use of recurrence relations for the Struve and Bessel functions of the second kind H ν and Y ν (ν ∈ N), respectively. To solve Eqs. (A3), it is instructive to introduce dimensionless variables scaled by a length scale = [e 2 /(4πεpμ p )] 1/(p+1) . For r 0 , we find that the obtained lattice configurations agree very well with the corresponding results obtained with a Coulomb interaction potential, V int (r i , r j ) ∼ 1/|r i − r j |. Since ≈ 30 nm athω = 1 meV, this condition is typically well satisfied in the situations considered in the main text. The resulting lattice structure

Normal modes
A two-dimensional lattice with N electrons has 2N elementary excitations, the so-called normal modes of the crystal. The normal-mode excitation spectrum of WCs can be calculated from the the system's elasticity matrix K.
Starting with Eq. (A1), the elasticity matrix is obtained from the second-order derivatives of V p with respect to the spatial coordinates. In the general case of arbitrary p 2 and the interaction potential in Eq. (2), we find that and where α, β ∈ {x, y}, α = β and the function g is given by The eigenmodes of the system are then calculated from the eigenvalues of the elasticity matrix K αβ mn = ∂ 2 V p /(∂α m ∂β n ).

APPENDIX B: IMPURITY-INDUCED POSITIONAL DISORDER: EQUIDISTANCE MEASURE
Random dislocations of single electrons from their lattice sites r 0 i may not only affect the lattice structure of a WC, but also the detection scheme and spin couplings discussed in the main text. For a simple estimate of how severe the impact of impurities on the lattice is, we consider N imp randomly distributed Gaussian confinement potentials in addition to the potential in Eq. (A1) and draw both size and depth of these local confinement potentials from normal distributions. For our calculations, we assume that they are localized on a nanometer length scale and have a depth of the order of ∼meV. In a monolayer TMD, such defects could be, e.g., atomistic defects [62]. Starting from Eq. (1), we take these into account by adding a disorder term, with V rand (r 1 , . . . , r N ; with random variables D j ∼ meV and σ j ∼ nm (where both means and standard deviations are of these orders), where {s j } 1 j N denote the positions of the impurities. To illustrate how this impurity model affects the lattice site distribution r 0 i (i = 1, . . . , N) of a small system, an exemplary numer-ical result obtained with N = 8 is shown in Fig. 6(a). The same result, but obtained in the presence of two randomly located (in the lattice) local harmonic potentials, is shown in Fig. 6(b). Averaging over many such instances and calculating the density-density correlations in the WC yields a measure of how much the crystal structure is affected by the presence of disorder. Similarly, here we look at another measure, χ , which quantifies how equidistantly the lattice sites r 0 i are distributed in the x-y plane by summing up the distances between nearest neighbors: Below we show that χ = 2 2/ √ 3 (χ = 1) for an equidistantly (completely randomly distributed) set of points r 0 1, . . . , N). By increasing the number of impurities for a given system size, i.e., increasing the impurity density n imp as compared to the electron density n, χ drops from its maximum value very fast, see Fig. 6. As would be intuitively expected, this underlines that n imp n should be fulfilled in any experiment to maximize the chances to observe charge ordering in regular electron lattices.
We briefly show that χ is upper bounded by χ max = r m /r ∞ = 2 √ 2/3 1/4 ≈ 2.15 [63]. This can be achieved by (i) calculating an upper bound for r m = i min j =i |r i − r j |/N and (ii) estimating r ∞ = 1/(2 √ n) as a function of the average electron density n: (a) In a close-packed lattice with an average nearest-neighbor distance r m , the unit cell occupies an area in A uc = √ 3r 2 m /2. The electron density is then given by n = 2/( √ 3r 2 m ). (b) The mean number of lattice sites in a sector of area A k = π r 2 /k is m = nA k . The probability of finding N sites in A k is given by a Poisson distribution P(N sites in A k ) = m N e −m /N!. Hence, we obtain the probability that two lattice sites are separated by a distance |r 0 i − r 0 j | smaller than a given r, P < (r) := P(|r 0 i − r 0 j | < r) = 1 − exp(−nπ r 2 /k). Therefore, we obtain for the mean of the distance distribution (k = 1): Combining the findings from (i) and (ii), we obtain an upper bound for χ , χ max = 2 √ 2/3 1/4 ≈ 2.15. Similarly, it can be shown that χ = 1 for a random distribution of lattice sites.
In our numerical calculations, we have seen that the detection scheme is only weakly affected by disorder if the impurity density n imp /n 0.1. The influence of disorder on cooperative resonances such as the ones discussed in the main text has also been investigated in Ref. [43].

APPENDIX C: FINITE TEMPERATURE EFFECTS
We first provide a simple estimate of the melting temperature T m of a WC by employing the Lindemann criterion, which has been used extensively in the literature [64,65]. It states that, in a lattice with charge-carrier density n, melting occurs if the root-mean square (RMS) displacement of a charge carrier from its lattice site r 0 i exceeds a certain fraction of the interparticle distance a. The RMS displacement can be obtained from the thermally occupied vibrational (normal) modes of the system at thermal equilibrium. Accordingly, the melting temperature T m and electron density n can be related. Although it is only a phenomenological criterion, it provides an efficient tool for estimating the melting temperature of a lattice. The thereby numerically calculated melting curves, obtained using typical material parameters of GaAs and monolayer TMD systems, respectively, are shown in Fig. 7. For the latter, we estimate melting temperatures of the order of T m ∼ 5 K, which is in agreement with previous estimates [14].
Cooling the system into its motional ground state puts more demanding constraints on temperature than considering melting only. We compare the thermal energy set by k B T to the mode frequencies n and calculate the thermal mode occupationn th = [exp(h n /(k B T ) − 1] −1 . Figure 8 shows that for the center-of-mass (COM) mode, it isn th 1 at T (1 − 5) K andhω 0.5 meV.

APPENDIX D: SPIN-SPIN INTERACTIONS: DERIVATION OF COUPLING CONSTANT
We estimate the spin-coupling strength J as given by Eq. (4) in the main text. For this, we model the interaction potential V int (r i , r j ) between two electrons at r i and r j with a Coulomb potential ∼1/|r i − r j |. In the parameter regime considered here, this (i) simplifies the calculation and (ii) yields the same results as obtained with the Keldysh interaction potential from Eq. (2) to a very good approximation, as confirmed by our numerical calculations.

Estimate of spin-coupling constant
We calculate the spin-exchange interaction between two electrons from the energy difference between the spin-singlet 1 |r 1 − r 2 | b (r 1 ) a (r 2 ), where a/b (r) = φ a/b (r) · χ a/b (r) denotes the electronic wave function and the labels a and b refer to the two electrons located at around r 0 a/b = (x 0 a/b , y 0 a/b ). We model the wave functions with a Gaussian wave packet of width wr ZPF (see Sec. D.2), where i ∈ {a, b}, and take into account spin-valley locking by setting the Bloch wave if the i electron, i ∈ {a, b}, is in a spin-|↑ (spin-|↓ ) state. Next, we evaluate the exchange integral in the spin basis spanned by |↑↑ , |↑↓ , |↓↑ , |↓↓ . With the electrons in different valleys (i.e., opposite spins), by performing some of the integrations analytically, we find for J ab in Eq. (D2) that with a TMD lattice constant a TMD ≈ 0.3 nm for MoX 2 (X = S, Se) [67] and |K − K | = 8π/(3a TMD ). K 0 denotes the modified Bessel function of the second kind. Inserting our numerical results for w, and in particular with wr ZPF a TMD , we find numerically that J KK ab evaluates to negligibly small values as compared with J KK ab , with which we denote the case where the two electrons are in the same valley. We find that J KK ab /J KK ab ∼ a TMD /(wr ZPF ) and that typically J KK ab is several orders of magnitude smaller than J KK ab . For J KK ab , we find an analytical expression and insert TMD parameters such that where we have expressed the electron density as n = 2/( √ 3a 2 ) for a triangular lattice. Similarly, we find for the overlap integral S in Eq. (D2) that In the low-density regime considered here, we find S 1 such that J ≈ J ab in Eq. (D1) to a very good approximation.
Also the Coulomb integral C in Eq. (D2) can be calculated analytically by employing the convolution and Parseval's theorems. Defining f a/b (r) := |φ a/b (r)| 2 and g(x) = 1/|x|, we insert TMD parameters and find that n[10 10 cm −2 ]w 2 where I 0 is the modified Bessel function of the first kind.
Putting our results together, we find that J KK is several orders of magnitude smaller than J KK for realistic parameters. Evaluating the Coulomb interaction Hamiltonian in the spin basis, with these results we obtain the spin model from Eq. (3) in the main text. Finally, putting the results from Eqs. (D3)-(D4) and Eq. (D1) together, we obtain coupling strengths in the range ∼(5 − 30) μeV at densities n n cr , as presented in Fig. 2(a) of the main text.

Width of ansatz wave function
We have considered two approaches to calculate w, for which we have found good agreement. (i) Mean-field approximation: First, we (iteratively, until the result is found to be converged) calculate the effective potential seen by a single electron due to the neighboring electrons by summing up the Coulomb interaction terms. From this potential, we calculate the wave function with a Gaussian ansatz, which yields the width of the wave function ∼w. (ii) Harmonic model: Second, we consider an expansion of the individual electron displacements in the set of collective displacement modes. In this way, we relate w to the normal modes which we have calculated before, where the mode frequencies n are expressed in units of the external confinement ω. For a confinementhω = 3 meV, we obtain r ZPF ≈ 5 nm.

APPENDIX E: OPTICAL READOUT: NUMERICAL AND ANALYTICAL TREATMENT
Here we first briefly summarize how we solve the scattering problem of light incident on a finite WC and then continue with an analytical treatment of the scattering problem for an infinite lattice. The latter provides us with more physical insight into the problem and is useful for optimizing the beam parameters to maximize the transmission or reflection contrast of the readout scheme.

Finite arrays
The principle behind the optical readout scheme discussed in the main text is based on a cooperative resonance effect as described in detail in Refs. [42,43]. As depicted in Fig. 1, we consider a Gaussian beam E in (x , y , z ) incident on the xy 125101-8 plane with a tilt angle θ and azimuthal angle φ, where which is scattered from a lattice of dipoles. Here we have introduced the coordinates ⎛ In Eq. (E1), E 0 denotes the the beam amplitude, w 0 and w(z) = w 0 1 + (z/z R ) 2 are beam waist and radius at z, respectively, z R = π w 2 0 /λ is the Rayleigh length and ϕ = arctan z/z R refers to the Gouy phase of the laser beam [68]. e pol encodes the polarization of the beam. For the results presented in Fig. 3 in the main text, we consider elliptically polarized light with At small detunings 0 from the transition frequency ω 0 , | 0 | ω 0 , each lattice site is modeled as a dipole with polarizability with the radiative (nonradiative) linewidth γ r (γ nr ). In general, the radiative linewidth γ r can be enhanced by the presence of a medium [69], especially for high refractive-index materials like TMDs [70]. At low temperatures as considered here, hexagonal boron nitride (hBN) encapsulated TMD monolayers feature optical transitions with a radiative linewidthhγ 0 ∼ meV [71]. In our calculations, we assume that the excitons are localized on a length scale much smaller than the wavelength, i.e., a B λ. Those spatially localized quantum emitters show much narrower linewidths ∼100 μeV [72][73][74][75]. Using Fermi's golden rule, the increased radiative lifetime of such localized excitons can be calculated, yielding a significantly enhanced emission time as compared to free excitons [76]. We estimate the radiative linewidth of a localized exciton to be of the order ofhγ r ≈ 4π/3(a B /λ) 2 γ 0 ≈ 10 −5 γ 0 ≈ 10 −2 μeV. In the results presented in the main text, we have considered γ nr = 0.
Given the Gaussian input field from Eq. (E1), we solve the Lippmann-Schwinger Eq. (4), with the Green's function [77] G αβ k, r, r 0 with α, β ∈ {x, y, z}. We solve Eq. (4) self-consistently for various angles of incidence θ and φ, beam profiles, detunings, FIG. 9. Transmission at normal incidence (θ = 0) for a square lattice. Other numerical parameters as in Fig. 3 in the main text. and electron lattices. At normal incidence, i.e., θ = 0, the resulting transmission and reflection signals depend on the lattice constant (see Fig. 9) but clearly not on φ. For 0 < θ < π/2, the transmission and reflection contrasts can be of the order of a few percent. An analytical derivation of the maximum contrast for an infinite lattice, depending on the angle of incidence θ and detuning 0 , is presented in Sec. E 3.

Faraday rotation
In the main text, we investigate the Faraday rotation angle according to Eq. (6). For the results in Fig. 4, we consider an incoming beam at normal incidence (θ = 0) with We consider N ≡ N ↑ + N ↓ dipoles which are located at lattice sites r 0 i with the spins assigned randomly to these lattice lattices for fixed N ↑ and N ↓ . Next we average over sufficiently many (∼10 4 ) instances of such configurations to calculate the Faraday rotation.
In Fig. 4, we show results for γ nr = 0. For γ nr > 0, the maximum Faraday rotation decreases and shifts toward more highly detuned frequencies, cf. Fig. 10.

Infinite arrays
Here we consider light scattering off an (infinite) twodimensional lattice of dipoles. If the transition dipole is parallel to the unit vectorê, the electric field at position r satisfies the equation where k 0 denotes the wave number of the transition. This equation can be readily solved using a Fourier transform, 125101-9 FIG. 10. Faraday rotation for different γ nr and the same numerical parameters as in Fig. 3 at N ↑ = 15, N ↓ = 5. Also shown is the maximum Faraday signal as a function of γ nr /γ r . assuming that the medium surrounding the lattice is translationally invariant in the plane of the lattice. One obtains where A is the area of the unit cell and where the sum runs over all reciprocal lattice vectors, denoted by B.
For an incident plane wave with momentum k, only a single term contributes to the sum in Eq. (E7). The plane wave will be Bragg scattered to momenta k + B. However, for sufficiently small lattice constants, |k + B| > k 0 for any B = 0, such that all nonzero scattering orders are evanescent.
In this case, the far field is completely described by We were able to turn the matrix inversion into a simple division by using the fact thatêê † is a projector. It is straightforward to show that the condition |k + B| > k 0 is equivalent to |B| > 4π/λ. For a square lattice, one obtains a < λ/2, while for a triangular lattice a < λ/ √ 3. To simplify Eq. (E9) further, we consider the special case that the array is placed in free space. The free-space Green's function is given by where and P ± (k) denotes the projector onto transverse polarizations for waves propagating up (+, z > 0) or down (−, z < 0). Explicitly, the P ± (k) projects onto the two-dimensional space spanned bŷ where we defined the angles θ and φ according to k x = k 0 sin θ cos φ, k y = k 0 sin θ sin φ, k z = k 0 cos θ.
(E13) We note that k z is always taken to have a positive real (|k| < k 0 ) or imaginary (|k| > k 0 ) part. When |k| < k 0 , all angles are real, and the vector (k x , k y , −k z ) is simply the wave vector of the incident wave. We also point out that the Green's function is discontinuous at z = 0. Right at z = 0, one should take 1 G(k, 0) = i 4k z e ik z |z| [P + (k) + P − (k)].
We focus on a circularly polarized transition, that is, When there is no Bragg scattering, it is easy to see that ImG(k) = Im G(k, 0) such that e † ImG(k)ê = i 4k 0 1 + cos 2 θ cos θ .
Equation (E19) has a simple physical interpretation. The light probes a collective resonance with energy˜ and radiative linewidth . The vectorsv ± correspond to projections of the transverse polarizations onto the transition dipole. The response of the lattice is maximized when E in ∝v − , which corresponds to an elliptic polarization whose projection onto the xy plane is circular. The expression allows us to immediately read off the reflection and transmission coefficients: Both r and t should be thought of as 2 × 2 matrices acting on the subspaces of transverse polarizations. For a fixed incident polarizationê in , we may further compute the intensity reflection and transmission coefficients. They are given by The intensity coefficients satisfy R + T = 1 when γ nr = 0 as required.
In practice, we would like to infer the rotational symmetry of the lattice via the dependence of˜ on φ. Choosing the optimal polarizationê in =v − (θ, φ), the maximum contrast in reflection for a fixed value of θ is given by where˜ min = min φ˜ (θ, φ) and similarly for˜ max . For simplicity, we set γ nr = 0, which implies that the contrast in transmission is equal to the contrast in reflection. We are free to choose 0 to maximize the contrast. Writing 0 = −(˜ min +˜ max )/2 + δ, the contrast can be expressed 125101-11 as R = δ¯ / 2 (δ 2 / 2 −¯ 2 / 2 + 1/4) 2 +¯ 2 / 2 , (E27) where¯ = (˜ max −˜ min )/2. In the limit¯ , the expression simplifies to R ≈ δ¯ / 2 (δ 2 / 2 + 1/4) 2 . (E28) It is easy to show that the contrast is maximized by choosing The value of¯ can be computed numerically. The results for a square and triangular lattice are shown in Figs. 11 and 12. As a final remark, we mention that by measuring the transmission coefficient for a component of the electric field that is neither parallel nor perpendicular to the incident field, it is possible to observe dispersive (asymmetric) line shapes. Such features could potentially enhance the sensitivity.