This paper was converted on www.awesomepapers.org from LaTeX by an anonymous user.
Want to know more? Visit the Converter page.

\externaldocument

[sm-]sm \externaldocumentsm

Competing Generalized Wigner Crystal States in Moiré Heterostructures

Shu Fay Ung su2254@columbia.edu Department of Chemistry, Columbia University, New York, NY, United States    Joonho Lee joonholee@g.harvard.edu Department of Chemistry, Columbia University, New York, NY, United States Department of Chemistry and Chemical Biology, Harvard University, MA, United States    David R. Reichman drr2103@columbia.edu Department of Chemistry, Columbia University, New York, NY, United States
Abstract

We present a comprehensive study of generalized Wigner crystals across various filling factors for system sizes up to 162 holes employing Hartree-Fock theory and explicitly correlated wave function approaches. While we find broad agreement with the behavior observed in experiments and classical Monte Carlo simulations, we highlight the fact that the Hartree-Fock energy landscape appears to be remarkably complex, exhibiting many competing states, both ordered and disordered, separated by energies of a fraction of \sim 1 meV/hole. We demonstrate which of the located states are metastable by performing a stability analysis at the Hartree-Fock level. Correlated wave function methods furthermore reveal small correlation energies that are nevertheless large enough to tip the balance of state ordering found within Hartree-Fock theory.

I Introduction

The interest in a novel class of physical systems–so-called moiré systems–has recently gained tremendous momentum in the condensed-matter physics and materials science communities [1]. These systems are created from two=dimensional layers stacked with a small misalignment such that a characteristic long wavelength moiré pattern forms. The associated moiré period parametrizes the kinetic and interaction energy scales in the stacked system, providing a control knob to alter the electronic Hamiltonian. In particular, the interaction strength between charges can be tuned by changing the twist angle or the layer constituents, rendering such systems a versatile platform for studying strong-correlation physics. Experiments on moiré superlattices have revealed a rich phase diagram of strongly-correlated electronic states, including Mott insulators, generalized Wigner crystals (GWCs), stripe phases, and liquid crystals by varying the interaction strength together with the electron filling factor [2, 3, 4, 5, 6, 7].

Despite the large moiré unit cell, the low-energy physics can be described by an emergent single-particle Hamiltonian (i.e. the continuum model) sharing the periodicity of the moiré superlattice. We focus on transition metal dichalcogenide (TMD) heterobilayers, where carriers near the Fermi level are localized in one layer only and experience a periodic moiré potential due to the spatial variation in the valence-band maximum of the other layer [8]. Since the valence-band extrema located at the KK and KK^{\prime} momentum points are decoupled and exhibit a large spin-splitting caused by spin-orbit interactions, spin indices are locked with valley degrees of freedom, giving rise to an effective two-fold spin degeneracy in the continuum model [9].

The accurate description of strongly-correlated states, however, requires the treatment of many-body interactions. A natural approach is to incorporate long-range Coulomb interactions into the continuum model description (i.e. the interacting continuum model) [10, 11, 12]. This framework resembles that of the two-dimensional electron gas (2DEG) with a pinning moiré potential, which may harbor exotic physics arising from Coulomb frustration [13] and glassy dynamics [14] in addition to well-known 2DEG phenomena such as standard Wigner crystallization [15, 16, 17].

Previous theoretical work on GWCs have employed Hartree-Fock (HF) [12, 18, 19, 20, 21], density functional theory [22], classical and quantum Monte Carlo (MC) [23, 24, 25, 22], and exact diagonalization (ED) [26, 27] to probe the ground state phase diagram. While classical MC is reasonable deep within the crystalline regime, a comprehensive picture requires a quantum treatment beyond mean-field theory to capture strong correlation effects. Such studies have largely focused on the magnetic behaviour at filling factors ν=1/3,2/3,\nu=1/3,2/3, and 1. Essentially no work has been done to characterize charge order across all relevant filling factors using explicitly quantum mechanical approaches in large systems.

In this paper, we present a comprehensive quantum mechanical study of GWCs in hole-doped TMD heterobilayers at ten distinct filling factors up to ν=1\nu=1. We first characterize the quantum energy landscape of HF solutions for systems up to 162 holes, establishing that multiple metastable states (as defined by a positive-definite spectrum of the HF Hessian matrix) generically exist, which is suggestive of glassy behaviour in the system. The states are nearly degenerate in energy and may exhibit ordered or disordered configurations. Next, we verify that the complexity in the HF energy landscape persists upon incorporating electron correlation effects via correlated wave function approaches. Correlation energies were found to be small yet sufficiently large to alter the energetic ordering of the states revealed by Hartree-Fock theory. Most notably, our results suggest that quantum fluctuations can stabilize partially crystalline ground states. We demonstrate that such approaches–which have reasonable scalings with system size–may be used to accurately estimate the correlation energies of GWC states. Hartree atomic units are used throughout the text unless stated otherwise.

II Model and methods

II.1 Interacting continuum model

The low-energy physics of TMD heterobilayers is described by a continuum model constructed from the topmost valence band [9]. Carriers near the Fermi level are localized in one layer only, with the effect of the other layer appearing through local variations in the band gap that correspond to moiré length-scale variations in the local atomic registry [8]. This manifests in the carriers experiencing a moiré potential Δ(r)\Delta(\vec{r}) that is periodic with respect to the superlattice, giving rise to the single-particle moiré Hamiltonian [9]

h^i\displaystyle\hat{h}_{i} =i22m+Δ(ri),\displaystyle=-\frac{\nabla_{i}^{2}}{2m^{*}}+\Delta(\vec{r}_{i}), (1)

where mm^{*} is the carrier effective mass. To the lowest order in a harmonic expansion, we have

Δ(r)=2Vmj=13cos((gjr+ϕ)),\Delta(\vec{r})=2V_{m}\sum_{j=1}^{3}\cos{\left(\vec{g}_{j}\cdot\vec{r}+\phi\right)}, (2)

where VmV_{m} is the potential strength, ϕ\phi characterizes the shape of Δ(r)\Delta(\vec{r}) in the moiré unit cell, and gj\vec{g}_{j} are the moiré reciprocal lattice vectors in the first shell which satisfy C3C_{3} symmetry. The emergence of Δ(r)\Delta(\vec{r}) leads to zone-folding of the monolayer Brillouin zones, yielding smaller moiré Brillouin zones (mBZ) that host flat bands. VmV_{m} and ϕ\phi are intrinsic material properties that can be obtained from fitting to ab initio band structures [9, 10].

Taking the carriers to be holes, we set m=0.45mem^{*}=-0.45m_{e} (mem_{e} is the bare electron mass), Vm=15V_{m}=-15 meV, and ϕ=45\phi=45^{\circ} for \chWSe2/\chWS2. At zero twist, the 4% lattice mismatch between the \chWSe2 and \chWS2 monolayers produces a moiré superlattice with a lattice constant of Lm=8.3L_{m}=8.3 nm [10]. Spin-valley locking [9] implies an effective two-fold degeneracy in the model, thus we identify the filling factor ν=2\nu=2 (2 carriers per moiré unit cell) as “full-filling”. We only examine ν1\nu\leq 1 since interactions renormalize bands more strongly with higher charge densities, leading to remote band mixing that demands caution with respect to the single-band continuum model for ν>1\nu>1 [27].

The interacting continuum model is obtained by adding Coulomb interactions

H^=i=1Nh^i+12ij=1N1ϵ|rirj|.\hat{H}=\sum_{i=1}^{N}\hat{h}_{i}+\frac{1}{2}\sum_{ij=1}^{N}\frac{1}{\epsilon|\vec{r}_{i}-\vec{r}_{j}|}. (3)

Here, we used a conventional 1/ϵ1/\epsilon screened Coulomb potential with ϵ\epsilon the effective dielectric constant, although the effects of gate-screening [18, 27, 25] and a Rytova-Keldysh-type screening [28, 29] can also be considered (see Appendix F for a discussion). We take ϵ=7\epsilon=7 to reflect the dielectric environment of hexagonal boron nitride (hBN) encapsulating layers often used in experiments [30, 3, 5, 4]

II.2 Generalized Hartree-Fock

Our study employed the generalized Hartree-Fock approximation, where spin-orbitals {χi(r)}\{\chi_{i}(\vec{r})\} comprising the Slater determinant |ΨHF\mathinner{|{\Psi_{\text{HF}}}\rangle} possess both spin-up and spin-down character, i.e. they retain a 2-component spinor structure

χi(r)=[χi(r)χi(r)],\chi_{i}(\vec{r})=\begin{bmatrix}\chi_{i}^{\uparrow}(\vec{r})\\[5.0pt] \chi_{i}^{\downarrow}(\vec{r})\end{bmatrix}, (4)

where χiσ(r)\chi_{i}^{\sigma}(\vec{r}) are spatial orbitals in the spin-σ\sigma sector. By variationally minimizing the HF energy EHF=ΨHF|H^|ΨHFE_{\text{HF}}=\mathinner{\langle{\Psi_{\text{HF}}}|}\hat{H}\mathinner{|{\Psi_{\text{HF}}}\rangle} with respect to the orbitals {χi}\{\chi_{i}\}, we obtain a set of integro-differential equations that can be self-consistently solved for the optimal {χi}\{\chi_{i}\}. We represented the χiσ\chi_{i}^{\sigma} in a plane-wave basis and performed supercell calculations at the Γ\Gamma-point of the mBZ. By imposing periodic boundary conditions only over the supercell, we can naturally obtain solutions that are expected to break the superlattice translational symmetry at fractional filling factors [3, 5, 7, 24, 25, 27]. At each filling factor, we found at least one HF solution starting from solutions of the continuum model with matrix elements in the mixed-spin sectors of the initial density matrix 𝐏\mathbf{P} (i.e. 𝐏\mathbf{P}^{\uparrow\downarrow} and 𝐏\mathbf{P}^{\downarrow\uparrow}) drawn from a standard Gaussian distribution 𝒩(0,1)\mathcal{N}(0,1), normalized by 103×max|𝐏|10^{-3}\times\max|\mathbf{P}|. This “noisy” continuum model guess served to unbias solutions away from collinearity. In addition, we selectively located other solutions via physically motivated initial guesses and determined the stability of each state by diagonalizing the HF Hessian matrix. We considered systems with up to 162 holes and, where possible, performed finite-size analyses as detailed below. Further computational details are available in Appendix C.

II.3 Estimating the charge gap

The charge gap Δc\Delta_{c} is defined as

Δc\displaystyle\Delta_{c} =limNΔc(N),\displaystyle=\lim_{N\to\infty}\Delta_{c}(N), (5)
Δc(N)\displaystyle\Delta_{c}(N) =E(N+1)+E(N1)2E(N)\displaystyle=E(N+1)+E(N-1)-2E(N) (6)
=[E(N)E(N+1)]+[E(N1)E(N)]\displaystyle=-\left[E(N)-E(N+1)\right]+\left[E(N-1)-E(N)\right]
=EEA+EIP,\displaystyle=-E_{\text{EA}}+E_{\text{IP}}, (7)

where E(N)E(N) is the total energy of a NN-particle system, EEA=E(N)E(N+1)E_{\text{EA}}=E(N)-E(N+1) is the electron affinity and EIP=E(N1)E(N)E_{\text{IP}}=E(N-1)-E(N) is the ionization potential. Assuming that the orbitals in the NN and N±1N\pm 1 ground states are identical within Hartree-Fock theory, Koopman’s theorem [31] allows us to identify

EEA\displaystyle E_{\text{EA}} =εLUMO,EIP=εHOMO,\displaystyle=-\varepsilon_{\text{LUMO}},\quad E_{\text{IP}}=-\varepsilon_{\text{HOMO}}, (8)

where εLUMO\varepsilon_{\text{LUMO}} (εHOMO\varepsilon_{\text{HOMO}}) is the lowest unoccupied (highest occupied) molecular orbital. This gives

Δc(N)εLUMOεHOMO.\displaystyle\Delta_{c}(N)\approx\varepsilon_{\text{LUMO}}-\varepsilon_{\text{HOMO}}. (9)

In the thermodynamic limit, the approximation in (9) should become an exact equality since the addition or removal of a single particle will have minimal effect on the ground state orbitals. Evaluating (9) instead of (6) furthermore has the advantage of ensuring Δc0\Delta_{c}\geq 0 while being computationally cheaper since it only requires one self-consistent loop. However, finite-size errors inevitably arise in practical calculations performed on a finite number NN of particles. While periodic boundary conditions help to reduce finite-size effects, non-negligible errors originating from spurious “self-interactions” between particles and their periodic images still remain–especially for singular long-ranged potentials such as the Coulomb interaction [32].

These contributions can be mitigated by employing corrections involving the triangular lattice Madelung energy, given by [33]

εM=3.9210342ϵπrs1N1/2,\varepsilon_{M}=\frac{-3.921034}{2\epsilon\sqrt{\pi}}r_{s}^{-1}N^{-1/2}, (10)

where ϵ\epsilon and rsr_{s} are the effective dielectric constant and the Wigner-Seitz radius, respectively. For a NN-particle system, the Madelung-corrected charge gap is [34, 32]

Δ~c(N)\displaystyle\tilde{\Delta}_{c}(N) =Δc(N)+2|εM|,\displaystyle=\Delta_{c}(N)+2|\varepsilon_{M}|, (11)

which improves the scaling in the finite-size errors from 𝒪(N1/2)\mathcal{O}(N^{-1/2}) to 𝒪(N3/2)\mathcal{O}(N^{-3/2}) [35]. Given the above considerations, we estimated Δc\Delta_{c} in the thermodynamic limit by first computing the Madelung-corrected charge gap (11) using (9), then extrapolating–where appropriate–Δ~c(N)\tilde{\Delta}_{c}(N) to NN\to\infty assuming the extrapolation formula

Δ~c(N)=γN3/2+Δc(),\displaystyle\tilde{\Delta}_{c}(N)=\frac{\gamma}{N^{3/2}}+\Delta_{c}(\infty), (12)

with γ\gamma being a parameter to be fitted. Illustrative examples are provided in Appendix E.

II.4 Correlated methods

To go beyond mean-field theory and incorporate electron correlation effects, we leveraged several state-of-the-art approaches infrequently used in the study of moiré systems. We investigated systems smaller than our largest HF calculations yet too large to examine by ED–a compromise between system sizes where finite-size errors obscure true correlation effects and computational costs limit the availability of reference solutions for comparison. To find (near) exact energies from our distinct HF solutions, we employed the heat-bath configuration interaction (HCI) [36, 37] and coupled-cluster singles-doubles (CCSD) plus its improvement with perturbative triples (CCSD(T)) approaches [38, 39, 40]. Note that the scaling of CCSD and perhaps CCSD(T) enables their application to the largest systems we examined, making them attractive tools for future investigations of correlation effects in moiré systems [41]. A brief discussion of these approaches can be found in Appendix B.

III Results

III.1 Hartree-Fock survey across ν\nu

We first present a survey of the behavior of HF solutions across ν\nu for the interacting continuum model. The system sizes we considered are such that the GWCs are commensurate with the supercell and are listed in Table 1. All solutions discussed in this section were initiated from the continuum model guess as described previously and found to be stable.

Filling factor ν\nu System size NN
1/7 7, 28
1/4 4, 16, 36, 49, 64, 81
1/3 3, 12, 27, 48, 75, 108
2/5 10, 40
1/2 8, 18, 32, 50, 72, 98, 162
3/5 15, 60
2/3 6, 24, 54, 96, 150
3/4 12, 27, 48, 75, 108, 147
6/7 42, 168
1 9, 16, 36, 64, 81, 144
Table 1: The filling factors and system sizes studied with HF.

Fig. 1a shows charge gaps Δc\Delta_{c} obtained. At ν=1/4, 1/3, 1/2, 2/3, 3/4\nu=1/4,\ 1/3,\ 1/2,\ 2/3,\ 3/4, and 11 where there are sufficient data points, we performed the extrapolation as detailed in Section II.3. At ν=1/7, 2/5, 3/5,\nu=1/7,\ 2/5,\ 3/5, and 6/76/7 where there are only two data points, we estimated Δc()\Delta_{c}(\infty) using Δ~c(N)\tilde{\Delta}_{c}(N) computed with the largest system size. A large enough NbasismaxN_{\text{basis}}^{\text{max}} was used such that all Δ~c(N)\tilde{\Delta}_{c}(N) are converged to <0.5<0.5 meV. At each ν\nu, we also only considered Δc\Delta_{c} obtained from ordered solutions if they were found; the charge gaps otherwise originate from disordered solutions. Larger values of Δc\Delta_{c} can occur in ordered solutions (indicated at ν=1/2\nu=1/2 in Fig. 1a), but require tuning of the initial guess to find.

Refer to caption
Figure 1: Charge gaps Δc\Delta_{c} relative to the moiré scale interaction strength Um=1/ϵLmU_{m}=1/\epsilon L_{m} at various fractional filling factors ν\nu, where LmL_{m} is the moiré superlattice constant. The inset shows experimental measurements of the 2s exciton resonance energy E2sE_{\text{2s}} used to probe the charge gap in untwisted \chWSe2/\chWS2 heterobilayers [3]. All solutions depicted were obtained using the continuum model guess as discussed in the main text, except for the red cross which marks Δc/Um\Delta_{c}/U_{m} calculated from the ordered stripe solution.
Refer to caption
Figure 2: The charge density ρ(r)\rho(\vec{r}) at filling factors ν1\nu\leq 1 obtained from the unpolarized continuum model guess. The densities are plotted over the supercell and the grid lines show the moiré unit cell boundaries. The indices n1,n2+n_{1},n_{2}\in\mathbb{Z}^{+} represent a point r=n1Lm,1+n2Lm,2\vec{r}=n_{1}\vec{L}_{m,1}+n_{2}\vec{L}_{m,2} in the supercell, where Lm,i\vec{L}_{m,i} are the primitive moiré lattice vectors. Disordered configurations were generally observed across ν\nu; ordered phases were seen at ν=1/7,1/3,2/3\nu=1/7,1/3,2/3, and 1. At ν=1/2\nu=1/2, the continuum model guess yielded a disordered phase in contrast to the ordered phase obtained from the stripe guess (Fig. 4a).
Refer to caption
Figure 3: Fast Fourier transform (FFT) of the charge densities in Fig. 2 with the average density subtracted, i.e. |ρ(k)|=|FFT{ρ(r)ρ}||\rho(\vec{k})|=|\mathrm{FFT}\left\{\rho(\vec{r})-\langle\rho\rangle\right\}|. The indices κ1,κ2\kappa_{1},\kappa_{2}\in\mathbb{Z} represent a point g=κ1G1+κ2G2\vec{g}=\kappa_{1}\vec{G}_{1}+\kappa_{2}\vec{G}_{2} in reciprocal space, where Gi\vec{G}_{i} are the primitive reciprocal lattice vectors associated with the supercell.

We plot Δc\Delta_{c} normalized by the moiré scale interaction strength Um=1/ϵLmU_{m}=1/\epsilon L_{m}. As in previous experiments [3, 7] and theoretical work [27], we observed an asymmetry in Δc\Delta_{c} about ν\nu = 1/2 at the level of Hartree-Fock theory. The charge gaps for ν<1/2\nu<1/2 are larger than for ν>1/2\nu>1/2 due to the larger Coulomb-to-kinetic-energy ratio at lower average charge densities, which is reminiscent of the low-density behaviour in the 2DEG.

Compared to the exact diagonalization results of Morales-Durán et. al. [27] obtained with small system sizes (N16N\leq 16 for ν1\nu\leq 1), we found larger charge gaps and more stable GWCs. Our (Lm,m,Vm)(L_{m},m^{*},V_{m}) parameter values give Wm/Vm=0.16W_{m}/V_{m}=0.16, where Wm=1/mLm2W_{m}=1/m^{*}L_{m}^{2} is the moiré scale kinetic energy. While they predicted metallic phases across fractional filling factors at this value of Wm/VmW_{m}/V_{m}, we found insulating states with larger Δc/Um\Delta_{c}/U_{m} values than those computed in their insulating states. The discrepancy with the charge gaps reported in [27] appears not to arise entirely due to the use of Hartree-Fock, but due to differences in the model (i.e. whether a projection onto continuum model orbitals was involved), the system size, the basis set size, and the manner in which the charge gap was calculated (e.g. via (6) or (9)) as well. Replicating the ED calculations of [27], we found a closing of the gap computed using (6) towards metallic behavior with increasing Wm/VmW_{m}/V_{m} (i.e. decreasing VmV_{m}) compared to the values reported in Fig. 1.

Across all ν\nu, the HF solutions initialized from the continuum model guess consist of localized charge densities, predominantly of a disordered nature. Figs. 2 and 3 show the largest basis set charge densities across filling factors and their Fast Fourier transforms (FFTs). The extraction of FFT peak locations from Fig. 3 is detailed in Appendix D. The C3C_{3} symmetry of the moiré superlattice is generally preserved in Fourier space, implying that the charge densities also obey it (with ν=2/5\nu=2/5 perhaps the exception). The superlattice translational symmetry, on the contrary, can be broken.

The ordered phases we observed occur at ν=1/7,1/3,2/3,\nu=1/7,1/3,2/3, and 1, which exhibit triangular lattice and honeycomb charge orders as seen in experiment [5] and previous theoretical work [42, 7, 27]. At ν=1/3\nu=1/3 and 2/32/3, FFT peaks at k\vec{k} where |k|/|Gm|=1/3|\vec{k}|/|\vec{G}_{m}|=1/\sqrt{3} are observed, reflecting a unit cell defined by the vectors

L1\displaystyle\vec{L}_{1} =2Lm,1Lm,2,L2=Lm,1+Lm,2,\displaystyle=2\vec{L}_{m,1}-\vec{L}_{m,2},\quad\vec{L}_{2}=\vec{L}_{m,1}+\vec{L}_{m,2}, (13)

where Lm,i\vec{L}_{m,i} are the primitive moiré lattice vectors. Such a cell is constructed from occupied (vacant) moiré sites for ν=1/3\nu=1/3 (2/32/3). At ν=1/7\nu=1/7, FFT peaks occur at k\vec{k} where |k|/|Gm|=1/7|\vec{k}|/|\vec{G}_{m}|=1/\sqrt{7}, indicating a unit cell defined by the vectors

L1\displaystyle\vec{L}_{1} =3Lm,1Lm,2,L2=Lm,1+2Lm,2.\displaystyle=3\vec{L}_{m,1}-\vec{L}_{m,2},\quad\vec{L}_{2}=\vec{L}_{m,1}+2\vec{L}_{m,2}. (14)

In contrast to the case with ν=1/3\nu=1/3 and 2/32/3, the ν=1/7\nu=1/7 and 6/76/7 charge densities we obtained are not dual 111In the sense that swapping the occupied and vacant sites in the ν=1/7\nu=1/7 configuration gives the ν=6/7\nu=6/7 configuration to each other in Fig. 2 and can also be distinguished in Fourier space (Fig. 3). This is likely because our solution at ν=6/7\nu=6/7 is one of several local minima with similar configurations that are close in energy.

Disordered states at other filling factors also resemble MC predictions, e.g. the “columnar dimer crystal” at ν=2/5\nu=2/5 (Fig. 1b) [25]. Notably, our survey revealed a disordered (labyrinthine) stripe configuration at ν=1/2\nu=1/2, in contrast to the ordered stripe phase often discussed in the literature 222A similar labyrinthine solution was also obtained via classical Monte Carlo simulations in the Extended Data of [7]. The fragility of the ordered stripe state has been noted experimentally by Li et. al. [5], who highlight the sensitivity of this phase to lattice strain and local inhomogeneities. In the absence of such experimental considerations, however, we posit that our disordered solution is one of several nearly degenerate local minima (i.e. stable solutions) that differ slightly in their real-space configurations. As discussed further below, we can also find low-lying ordered states at filling fractions where the continuum model guess yields disordered configurations. This is accomplished through selecting physically motivated initial guesses.

III.2 Nearly degenerate solutions at ν=1/2\nu=1/2

We now focus on ν=1/2\nu=1/2 to probe the existence of nearly degenerate yet distinct HF solutions by supplying different initial guesses. In an attempt to recover the ordered stripe state, we explicitly constructed a polarized ordered stripe guess in addition to the continuum model guess employed earlier. We first present results obtained with N=18N=18 and Nbasis=211N_{\mathrm{basis}}=211 before discussing larger system and basis set size studies.

Figs. 4a-c show the charge configurations obtained from each initialization. Unsurprisingly, the solution from the ordered stripe guess is also an ordered stripe configuration (Fig. 4a), whereas the continuum model guess yielded a disordered (labyrinthine) stripe state (Fig. 4b). Both solutions were verified to be stable in the charge sector. Although several other HF states were discovered throughout our investigation, we only present one additional solution (Fig. 4c) that displays an inhomogeneous density across the honeycomb lattice and is unstable (i.e. a saddle point) at the HF level. These solutions can be clearly distinguished in Fourier space (Figs. 5a-c)–markedly, the ordered stripe phase exhibits peaks that obey C2C_{2} symmetry instead of the C3C_{3} symmetry of the moiré superlattice. The finite-size HF energy per particle for these states are εordered=41.233\varepsilon_{\mathrm{ordered}}=-41.233 meV, εdisordered=41.465\varepsilon_{\mathrm{disordered}}=-41.465 meV, and εpartiallypinned=41.114\varepsilon_{\begin{subarray}{c}\mathrm{partially}\\ \mathrm{pinned}\end{subarray}}=-41.114 meV, respectively. Taken together, our results suggest the existence of a rough energy landscape with many closely-lying local minima surrounding the ordered phase, which is reminiscent of a glassy system.

Refer to caption
Figure 4: The charge density ρ(r)\rho(\vec{r}) at ν=1/2\nu=1/2. (a)-(c) HF solutions. (d)-(f) CCSD solutions. Comparing the top and bottom rows, we observed that quantum fluctuations induced by correlation delocalizes ρ(r)\rho(\vec{r}). Note that ρ(r)\rho(\vec{r}) in (a) and (d) are relatively shifted by LmL_{m}, while (e) and (f) are relatively rotated by 6060^{\circ}.
Refer to caption
Figure 5: Fourier-space plots of the charge densities obtained at ν=1/2\nu=1/2. (a)-(f) FFT of the corresponding densities in Fig. 4. (g)-(i) FFT peaks extracted from (a)-(c), plotted directly in reciprocal space k=(k1,k2)\vec{k}=(k_{1},k_{2}). Large dots denote where the FFT signal is maximum and Gm\vec{G}_{m} are the primitive moiré reciprocal vectors.
Refer to caption
Figure 6: Convergence of the energy difference between the ordered and disordered (labyrinthine) stripe phases with system and basis set size. The quantity 2|Gmax/Gm|2|\vec{G}_{\mathrm{max}}/\vec{G}_{m}| is a measure of the basis set size, where |Gmax||\vec{G}_{\mathrm{max}}| is the maximum momentum included in our plane wave basis set. (For reference, Nbasis=211N_{\mathrm{basis}}=211 for N=18N=18 corresponds to 2|Gmax/Gm|2.52|\vec{G}_{\mathrm{max}}/\vec{G}_{m}|\approx 2.5.) As NN increases, the ordered stripe configuration becomes lower in energy but the difference |εorderedεdisordered||\varepsilon_{\mathrm{ordered}}-\varepsilon_{\mathrm{disordered}}| decreases.

The existence of competing states is furthermore revealed by the convergence of energies with system size. Fig. 6 illustrates the convergence of the energy difference per particle ε=E(N)/N\varepsilon=E(N)/N between the ordered and disordered stripe phases with system and basis set size. For the largest NbasisN_{\text{basis}}, while the ordered stripe evolves to be the lower energy state as NN increases, the two phases also become more nearly degenerate.

III.3 Further study with correlated methods

To investigate the effect of quantum fluctuations at ν=1/2\nu=1/2, we performed correlated calculations starting from the mean-field solutions found previously. Again considering the case N=18N=18 and Nbasis=211N_{\mathrm{basis}}=211, we leveraged only spin-collinear variants of HCI, CCSD, and CCSD(T) since the HF solutions we obtained were spin-collinear. We will use “vHCI” and “vHCI + PT2” to distinguish between the variational component of HCI and the full algorithm (variational plus perturbative components).

Our results show that the complexity of the energy landscape persists even after introducing electron correlation. Fig. 7 illustrates the small energy scales associated with the solutions: both the correlation energy and energy differences between the three solutions are 1\lesssim 1 meV. For reference, the relevant energy scales characterizing our system are UmVm𝒪(10)U_{m}\sim V_{m}\sim\mathcal{O}(10) meV and Wm𝒪(1)W_{m}\sim\mathcal{O}(1) meV. Similar energies obtained across vHCI, vHCI+PT2, CCSD, and CCSD(T) moreover provide confidence that these values are nearly exact. The three correlated solutions remain distinct, which suggests that the metastability associated with the stable HF solutions is retained. Each approach also qualitatively predicts charge densities that agree with each other and HF: panels (d)-(f) of Figs. 4 and 5 demonstrate how quantum fluctuations delocalize the charge density but still preserve the essential Fourier-space structure from HF.

Nevertheless, energetic corrections arising from correlation effects are non-trivial. The minuscule correlation energy is sufficient to alter the energetic ordering of states from HF, despite indicating that a single Slater determinant is adequate in providing a quantitative description of each state. For example, while HF predicts the lowest energy state (among our isolated solutions) to be the disordered stripe, all correlated methods predict it to be the partially pinned solution instead–albeit it was unstable at the HF level.

Refer to caption
Figure 7: Energy per particle ε\varepsilon relative to the HF ordered stripe solution, denoted with εref\varepsilon_{\mathrm{ref}}. We use vHCI to denote only the variational component of HCI, while vHCI + PT2 denotes the full HCI algorithm. The energy differences are <1<1 meV across all methods. Whereas the energetic ordering found by HF disagrees with correlated approaches, all correlated approaches yield orderings and energies in remarkable agreement with each other.

IV Conclusions

We have numerically studied GWCs in hole-doped TMD heterobilayers at various fractional filling factors using Hartree-Fock and correlated theories. As a concrete example, we focused on the untwisted \chWSe2/WS2 heterobilayer embedded in a hBN dielectric environment, which is an experimentally popular system for probing GWCs. Our results show qualitative agreement with experiment in the asymmetry of charge gaps about ν=1/2\nu=1/2 and in the charge orderings observed, particularly at ν=1/3,2/3,\nu=1/3,2/3, and 1. We moreover discovered multiple potentially disordered solutions that suggest the existence of a rugged energy manifold of nearly degenerate states. To probe this further, we focused on the ν=1/2\nu=1/2 case and found that the energy differences between distinct solutions are <1<1 meV at the mean-field level and beyond. Correlation energies are nevertheless non-trivial as demonstrated by a reversal in the energetic ordering of states between HF and correlated methods. Due to quantum fluctuations introduced by such approaches, HF charge densities become more delocalized while still retaining the essential HF features. In addition, we show for relatively small system sizes that fluctuations can stabilize partially pinned states relative to disordered stripe configurations [45, 46, 47, 48, 49, 50]. Future work will be devoted to characterizing these states. Our observations reveal the complexity of the moiré energy landscape–defined by several closely-lying metastable states–in both mean-field and correlated theories.

Acknowledgements.
We acknowledge funding support from NSF CHE–2245592 and thank Y. Zhang for help concerning the implementation of the model and A. Millis for useful discussions. Computations in this paper were primarily executed on the Ginsburg High Performance Computing Cluster at Columbia University. Some calculations were run on the FASRC Cannon cluster supported by the FAS Division of Science Research Computing Group at Harvard University. This work also used Bridges-2 Regular Memory (RM) nodes at the Pittsburgh Supercomputing Center and the Expanse cluster at the San Diego Supercomputer Center through allocation CHE220021 from the Advanced Cyberinfrastructure Coordination Ecosystem: Services & Support (ACCESS) program, which is supported by National Science Foundation grants #2138259, #2138286, #2138307, #2137603, and #2138296.

References

Appendix A Model details

The interacting continuum model is described by the Hamiltonian

H^=i=1N{i22m+Δ(ri)}+12ij=1N1ϵ|rirj|,\displaystyle\hat{H}=\sum_{i=1}^{N}\left\{-\frac{\nabla_{i}^{2}}{2m^{*}}+\Delta(\vec{r}_{i})\right\}+\frac{1}{2}\sum_{ij=1}^{N}\frac{1}{\epsilon|\vec{r}_{i}-\vec{r}_{j}|}, (15)

where mm^{*} is the carrier effective mass, Δ\Delta is the moiré potential, and ϵ\epsilon is the effective dielectric constant. In our calculations, we instead work with a scaled Hamiltonian H^\hat{H}^{\prime} defined by

H^\displaystyle\hat{H}^{\prime} =i=1N{i22+Δ(ri)}+12ij=1N1|rirj|,\displaystyle=\sum_{i=1}^{N}\left\{-\frac{\nabla_{i}^{\prime 2}}{2}+\Delta^{\prime}(\vec{r}\mkern 2.0mu\vphantom{r}^{\prime}_{i})\right\}+\frac{1}{2}\sum_{ij=1}^{N}\frac{1}{|\vec{r}\mkern 2.0mu\vphantom{r}^{\prime}_{i}-\vec{r}\mkern 2.0mu\vphantom{r}^{\prime}_{j}|}, (16)
H^\displaystyle\hat{H}^{\prime} =αH^,Δ=αΔ,r=βr,\displaystyle=\alpha\hat{H},\qquad\Delta^{\prime}=\alpha\Delta,\qquad\vec{r}\mkern 2.0mu\vphantom{r}^{\prime}=\beta\vec{r}, (17)
α\displaystyle\alpha =ϵ2m,β=mϵ,\displaystyle=\frac{\epsilon^{2}}{m^{*}},\qquad\beta=\frac{m^{*}}{\epsilon}, (18)

which effectively maps the interacting continuum model to a 2DEG in the scaled moiré potential Δ\Delta^{\prime} [11]. This form allows the interacting continuum model to be treated with DFT employing exchange-correlation functionals derived from the 2DEG. Note that if |Ψ\mathinner{|{\Psi^{\prime}}\rangle} is an eigenstate of H^\hat{H}^{\prime} with energy EE^{\prime}, then

H^|Ψ=αH^|Ψ\displaystyle\hat{H}^{\prime}\mathinner{|{\Psi^{\prime}}\rangle}=\alpha\hat{H}\mathinner{|{\Psi^{\prime}}\rangle} =E|Ψ,\displaystyle=E^{\prime}\mathinner{|{\Psi^{\prime}}\rangle},
H^|Ψ\displaystyle\implies\hat{H}\mathinner{|{\Psi^{\prime}}\rangle} =Eα|Ψ,\displaystyle=\frac{E^{\prime}}{\alpha}\mathinner{|{\Psi^{\prime}}\rangle}, (19)

i.e. |Ψ\mathinner{|{\Psi^{\prime}}\rangle} is also an eigenstate of the original Hamiltonian H^\hat{H} with energy E=E/αE=E^{\prime}/\alpha. Thus, we divided energies obtained from the scaled Hamiltonian by α\alpha to recover energies associated with the physical moiré system.

Appendix B Correlated methods

Any antisymmetric NN-particle state |Ψ\mathinner{|{\Psi}\rangle} can be represented in the basis consisting of a reference Slater determinant |Φ\mathinner{|{\Phi}\rangle} and all of its excitations, i.e.

|Ψ\displaystyle\mathinner{|{\Psi}\rangle} =(c0+iaciaa^aa^i+i<ja<bcijaba^aa^ba^ia^j+)|Φ,\displaystyle=\left(c_{0}+\sum_{ia}c_{i}^{a}\hat{a}_{a}^{\dagger}\hat{a}_{i}+\sum_{\begin{subarray}{c}i<j\\ a<b\end{subarray}}c_{ij}^{ab}\hat{a}_{a}^{\dagger}\hat{a}_{b}^{\dagger}\hat{a}_{i}\hat{a}_{j}+\cdots\right)\mathinner{|{\Phi}\rangle}, (20)

where cia,cijab,c_{i}^{a},c_{ij}^{ab},\dots are coefficients and a^i(a^i)\hat{a}_{i}^{\dagger}(\hat{a}_{i}) are creation (annihilation) operators defined over the orbital basis set {ϕi}\{\phi_{i}\} from which |Φ\mathinner{|{\Phi}\rangle} is constructed [31]. This representation becomes exact as the orbital basis approaches completeness. In principle, the coefficients can be determined through exact diagonalization, but the problem quickly becomes intractable as the Hilbert space scales exponentially with system size. Correlated wave function methods thus seek to approximate (20)\eqref{ci_wfn} starting from some known reference |Φ\mathinner{|{\Phi}\rangle}–often obtained from Hartree-Fock theory (hence by correlation, we refer to contributions beyond Hartree-Fock). HCI [36, 37] and coupled-cluster (CC) theories [38, 39, 40] are two such methods that we employ in our paper.

B.1 Heat-bath configuration interaction

Given a Hamiltonian H^\hat{H}, HCI finds an approximation to the ground state |Ψ\mathinner{|{\Psi}\rangle} consisting of a selected set of the most important determinants only. This is achieved in two stages: first, a variational wave function is constructed from a set of chosen determinants; second, perturbative corrections to the variational energy determine additional determinants to be added. The selection of determinants in each stage is tuned by two parameters, ϵ1\epsilon_{1} and ϵ2\epsilon_{2}, which control the tradeoff between speed and accuracy. We use vHCI and vHCI + PT2 to denote the variational stage and its correction with perturbation theory (i.e. the full HCI algorithm), respectively.

vHCI consists of iteratively expanding the selected space (which usually only includes the Hartree-Fock determinant initially) while minimizing the selected-space energy. At each iteration, new determinants |Φi\mathinner{|{\Phi^{\prime}_{i}}\rangle} connected to the current space by nonzero Hamiltonian matrix elements Hij=Φi|H^|ΦjH_{ij}=\mathinner{\langle{\Phi^{\prime}_{i}}|}\hat{H}\mathinner{|{\Phi_{j}}\rangle} are added according to the importance measure

f(|Φi)=maxj(|Hijcj|),f(\mathinner{|{\Phi^{\prime}_{i}}\rangle})=\max_{j}(|H_{ij}c_{j}|), (21)

where {cj}\{c_{j}\} are determinant coefficients. The variational stage is summarized as follows [36]:

  1. 1.

    Diagonalize the Hamiltonian in the selected space. Let {cj}\{c_{j}\} denote the determinant coefficients of the lowest eigenvector.

  2. 2.

    Generate all exterior determinants {|Φi}\{\mathinner{|{\Phi^{\prime}_{i}}\rangle}\} that satisfy |Hijcj|>ϵ1|H_{ij}c_{j}|>\epsilon_{1} for at least one determinant |Φj\mathinner{|{\Phi_{j}}\rangle} in the current space. Add these determinants to the selected space.

  3. 3.

    Repeat steps 1 and 2 until the number of new determinants generated is <1%<1\% of the number of determinants in the selected space.

The generation of {|Φi}\{\mathinner{|{\Phi^{\prime}_{i}}\rangle}\} in step 2 is achieved via deterministic heat-bath sampling. Starting from a reference determinant |Φj\mathinner{|{\Phi_{j}}\rangle}, only determinants |Φi\mathinner{|{\Phi^{\prime}_{i}}\rangle} connected to |Φj\mathinner{|{\Phi_{j}}\rangle} through matrix elements Φi|H^|Φj\mathinner{\langle{\Phi^{\prime}_{i}}|}\hat{H}\mathinner{|{\Phi_{j}}\rangle} that satisfy |Φi|H^|Φj|>ϵ|\mathinner{\langle{\Phi^{\prime}_{i}}|}\hat{H}\mathinner{|{\Phi_{j}}\rangle}|>\epsilon are sampled. The threshold ϵ\epsilon here is chosen to be ϵ=ϵ1/|cj|\epsilon=\epsilon_{1}/|c_{j}|. Duplicates are removed as the sampling algorithm is applied to each determinant in the current space.

Once the variational wave function |Ψ0=ici|Φi\mathinner{|{\Psi_{0}}\rangle}=\sum_{i}c_{i}\mathinner{|{\Phi_{i}}\rangle} is obtained, the partitioning of the Hamiltonian used in the perturbative stage is defined as

H^=H^0+V^,\displaystyle\hat{H}=\hat{H}_{0}+\hat{V},
H^0=ijHij|ΦiΦj|,H^0|Ψ0=E0|Ψ0.\displaystyle\hat{H}_{0}=\sum_{ij}H_{ij}\mathinner{|{\Phi_{i}}\rangle}\mathinner{\langle{\Phi_{j}}|},\quad\hat{H}_{0}\mathinner{|{\Psi_{0}}\rangle}=E_{0}\mathinner{|{\Psi_{0}}\rangle}. (22)

The second-order correction to the energy is thus

ΔE2\displaystyle\Delta E_{2} =a(iHaici)2E0Ea,\displaystyle=\sum_{a}\frac{\left(\sum_{i}H_{ai}c_{i}\right)^{2}}{E_{0}-E_{a}}, (23)

where the sum is over new determinants {|Φa}\{\mathinner{|{\Phi^{\prime}_{a}}\rangle}\} that give the matrix elements Hai=Φa|H^|ΦiH_{ai}=\mathinner{\langle{\Phi^{\prime}_{a}}|}\hat{H}\mathinner{|{\Phi_{i}}\rangle} and Ea=HaaE_{a}=H_{aa}. Since (23) is expensive to evaluate, vHCI + PT2 approximates ΔE2\Delta E_{2} using only terms in the sum that satisfy |Haici|>ϵ2|H_{ai}c_{i}|>\epsilon_{2}, i.e.

ΔE2a(i,|Haici|>ϵ2Haici)2E0Ea.\displaystyle\Delta E_{2}\approx\sum_{a}\frac{\left(\sum_{i,|H_{ai}c_{i}|>\epsilon_{2}}H_{ai}c_{i}\right)^{2}}{E_{0}-E_{a}}. (24)

Similar to vHCI, the set {|Φa}\{\mathinner{|{\Phi^{\prime}_{a}}\rangle}\} is generated by applying the heat-bath sampling algorithm to each determinant in |Ψ0\mathinner{|{\Psi_{0}}\rangle} using a threshold of ϵ=ϵ2/|ci|\epsilon=\epsilon_{2}/|c_{i}|.

B.2 Coupled-cluster theory

CC theory [51] assumes an ansatz for the NN-particle ground state of the form

|ΨCC=eT^|Ψ0,\mathinner{|{\Psi_{\text{CC}}}\rangle}=e^{\hat{T}}\mathinner{|{\Psi_{0}}\rangle}, (25)

where |Ψ0\mathinner{|{\Psi_{0}}\rangle} is a reference Slater determinant (often taken to be the Hartree-Fock solution). The cluster operator T^\hat{T} is given by

T^\displaystyle\hat{T} =T^1+T^2++T^N,\displaystyle=\hat{T}_{1}+\hat{T}_{2}+\cdots+\hat{T}_{N}, (26)
T^1\displaystyle\hat{T}_{1} =iatiaa^aa^i,\displaystyle=\sum_{ia}t_{i}^{a}\hat{a}_{a}^{\dagger}\hat{a}_{i}, (27)
T^2\displaystyle\hat{T}_{2} =i<ja<btijaba^aa^ba^ia^j,\displaystyle=\sum_{\begin{subarray}{c}i<j\\ a<b\end{subarray}}t_{ij}^{ab}\hat{a}_{a}^{\dagger}\hat{a}_{b}^{\dagger}\hat{a}_{i}\hat{a}_{j}, (28)
\displaystyle\qquad\vdots

where T^i\hat{T}_{i} are cluster operators associated with ii-particle excitations and tia,tijab,t_{i}^{a},t_{ij}^{ab},\dots are cluster amplitudes to be determined. These amplitudes together with the CC energy can be computed through a projection of the Schrödinger equation onto |Ψ0\mathinner{|{\Psi_{0}}\rangle} and its excitations:

Ψ0|eT^(H^E)eT^|Ψ0\displaystyle\mathinner{\langle{\Psi_{0}}|}e^{-\hat{T}}\left(\hat{H}-E\right)e^{\hat{T}}\mathinner{|{\Psi_{0}}\rangle} =0,\displaystyle=0, (29)
Ψ0|a^ia^aeT^(H^E)eT^|Ψ0\displaystyle\mathinner{\langle{\Psi_{0}}|}\hat{a}_{i}^{\dagger}\hat{a}_{a}e^{-\hat{T}}\left(\hat{H}-E\right)e^{\hat{T}}\mathinner{|{\Psi_{0}}\rangle} =0,\displaystyle=0, (30)
Ψ0|a^ja^ia^ba^aeT^(H^E)eT^|Ψ0\displaystyle\mathinner{\langle{\Psi_{0}}|}\hat{a}_{j}^{\dagger}\hat{a}_{i}^{\dagger}\hat{a}_{b}\hat{a}_{a}e^{-\hat{T}}\left(\hat{H}-E\right)e^{\hat{T}}\mathinner{|{\Psi_{0}}\rangle} =0,\displaystyle=0, (31)
\displaystyle\ \ \vdots

The resulting set of coupled polynomial equations (i.e. the coupled-cluster equations) can be solved (often iteratively) for the energy and amplitudes.

Note that |ΨCC\mathinner{|{\Psi_{\text{CC}}}\rangle} with T^\hat{T} untruncated is simply a non-linear reparametrization of (20), i.e. it is equivalent to exact diagonalization. The power of the CC approach, however, lies in the truncated form of the ansatz. The exponential form allows higher order excitations of |Ψ0\mathinner{|{\Psi_{0}}\rangle} to be included in the ansatz despite truncating T^\hat{T}. One of the most widely used variants is the truncation T^=T^1+T^2\hat{T}=\hat{T}_{1}+\hat{T}_{2}, known as CCSD. The accuracy of CCSD can be naïvely improved by including T^3\hat{T}_{3}, but this approach is often too costly. A more feasible way to improve CCSD solutions is to include contributions from T^3\hat{T}_{3} perturbatively, leading to the CCSD(T) method which is often called the “gold standard” of quantum chemistry.

Appendix C Computational details

C.1 Generalized Hartree-Fock

The spinor structure of generalized Hartree-Fock (GHF) orbitals breaks S^2\hat{S}^{2} and S^z\hat{S}_{z} symmetry, hence |ΨGHF\mathinner{|{\Psi_{\text{GHF}}}\rangle} is not constrained to be an eigenstate of S^2\hat{S}^{2} and S^z\hat{S}_{z}. Sacrificing symmetries of the exact solution provides the additional flexibility necessary to obtain the true variational minimum, which can be advantageous in studying magnetically frustrated systems that require a non-spin-collinear description [52]. As previous work [27] on TMD heterobilayers has suggested GWC ground states with non-collinear spins at some fractional filling factors, GHF is an appropriate mean-field method for our problem.

We performed all GHF calculations using a development version of Q-Chem [53]. Solutions were optimized through a hybrid scheme of first using the Direct Inversion in the Iterative Subspace (DIIS) method [54, 55] for 𝒪(1)\mathcal{O}(1) iterations before switching to the Geometric Direct Minimization (GDM) algorithm [56]. This approach combines DIIS’s ability to efficiently converge to the global minimum during the early iterations together with GDM’s ability to traverse challenging energy landscape features and robustly arrive at a local minimum. We executed supercell calculations composed of Ncell×Ncell\sqrt{N_{\text{cell}}}\times\sqrt{N_{\text{cell}}} unit cells at the Γ\Gamma point of the mBZ, which is analogous to performing unit cell calculations at Nk=NcellN_{k}=N_{\text{cell}} k\vec{k}-points adequately sampled from the mBZ. Periodic boundary conditions are only imposed across the supercell, which allow us to naturally obtain solutions where discrete translational symmetry on the moiré superlattice is broken. To break C3C_{3} symmetry, however, requires explicit encoding of the broken symmetry in the initial guess.

For our survey of GHF solutions across filling factors in Fig. 1 in the main text and Figs. 2, 8 and 3, we started with a “noisy” spin-unpolarized initial guess constructed from continuum model orbitals, i.e. eigenstates of (1) in the main text (henceforth referred to as the “continuum model guess”). By “noisy”, we mean that matrix elements of the spin \uparrow\downarrow and \downarrow\uparrow blocks of the density matrix 𝐏\mathbf{P} are randomly sampled from the standard Gaussian distribution 𝒩(0,1)\mathcal{N}(0,1), normalized by 103×max|𝐏|10^{-3}\times\max|\mathbf{P}|. At ν=1/2\nu=1/2, we additionally considered a spin-polarized stripe guess constructed from orthogonal Bloch orbitals of the form

ψk(r)i=1NeikRiϕi(r),\psi_{\vec{k}}(\vec{r})\propto\sum_{i=1}^{N}e^{i\vec{k}\cdot\vec{R}_{i}}\phi_{i}(\vec{r}), (32)

where NN is the particle number, ϕi(r)\phi_{i}(\vec{r}) are Gaussians centered at stripe lattice sites Ri\vec{R}_{i}, and k\vec{k} are crystal momenta that:

  1. (i)

    are restricted to the first Brillouin zone defined by the reciprocal space of the stripe lattice;

  2. (ii)

    satisfy the periodic boundary conditions imposed on the supercell.

To evaluate the basis set convergence of solutions, we performed a series of calculations with increasing basis set size NbasisN_{\text{basis}} starting from a small NbasisminN_{\text{basis}}^{\text{min}}. The converged solution at each NbasisN_{\text{basis}} was input as the initial guess for the next calculation using a larger NbasisN_{\text{basis}}, where the NbasisminN_{\text{basis}}^{\text{min}} solutions were first ensured to be stable.

C.2 Stability analysis

The self-consistent procedure used in GHF only guarantees solutions to be stationary points. To determine whether they are local minima or saddle points, we inspect the lowest eigenvalues of the electronic Hessian following the approach of [57]. Diagonalization is executed via a finite-difference implementation of the Davidson method, where for a vector v\vec{v} in the subspace, the Hessian-vector product is computed as

𝐇v[E(𝐂+)E(𝐂)2ξ]𝖳.\mathbf{H}\vec{v}\approx\left[\frac{\nabla E(\mathbf{C}_{+})-\nabla E(\mathbf{C}_{-})}{2\xi}\right]^{\mathsf{T}}. (33)

E\nabla E is the analytic energy gradient, while 𝐂±\mathbf{C}_{\pm} are the GHF solution coefficient matrices perturbed in opposite directions determined by the finite step size ξ\xi. If the solution is deemed unstable, we can displace the solution along the lowest eigenvector and use this as an initial guess in a new GHF calculation.

C.3 Spin-collinearity test

We evaluated the spin collinearity of our GHF solutions following the approach of [58]. A wave function |Ψ\mathinner{|{\Psi}\rangle} is considered to be collinear if it is an eigenstate of an axial spin operator S^c=cS^\hat{S}_{c}=\vec{c}\cdot\hat{\vec{S}} associated with some Cartesian coordinate, where c3\vec{c}\in\mathbb{R}^{3} is a unit vector and S^=(S^x,S^y,S^z)\hat{\vec{S}}=(\hat{S}_{x},\hat{S}_{y},\hat{S}_{z}) is the usual many-body spin vector operator. In other words, |Ψ\mathinner{|{\Psi}\rangle} is collinear if it satisfies the relation

S^c|Ψ=sc|Ψ,\displaystyle\hat{S}_{c}\mathinner{|{\Psi}\rangle}=s_{c}\mathinner{|{\Psi}\rangle}, (34)
sc{N/2,N/2+1,,N/21,N/2}.\displaystyle s_{c}\in\{-N/2,-N/2+1,\dots,N/2-1,N/2\}.

The collinearity test for a state consists of computing the quantities μ\mu and ε0\varepsilon_{0}, defined as

μ\displaystyle\mu =(S^c)2S^c2,\displaystyle=\langle(\hat{S}_{c})^{2}\rangle-\langle\hat{S}_{c}\rangle^{2}, (35)
ε0\displaystyle\varepsilon_{0} =S^=(S^x,S^y,S^z),\displaystyle=\norm{\langle\hat{\vec{S}}\rangle}=\norm{\left(\langle\hat{S}_{x}\rangle,\langle\hat{S}_{y}\rangle,\langle\hat{S}_{z}\rangle\right)}, (36)

where \norm{\cdot} denotes the Euclidean norm. The conditions for collinearity are thus:

  1. (i)

    μ=0\mu=0.

  2. (ii)

    ε0{N/2,N/2+1,,N/21,N/2}.\varepsilon_{0}\in\{-N/2,-N/2+1,\dots,N/2-1,N/2\}.

Condition (i) is obvious since it can easily be proved that μ=0\mu=0 if and only if |Ψ\mathinner{|{\Psi}\rangle} is an eigenstate of S^c\hat{S}_{c}. Condition (ii) is less straightforward and we refer the reader to [58] for its justification. In practice, we deem μ104\mu\lesssim 10^{-4} as describing a collinear state.

C.4 Correlated methods

We performed correlated calculations using only spin-collinear variants of HCI, CCSD, and CCSD (T). CCSD and CCSD(T) calculations were carried out using PySCF [59, 60, 61], whereas HCI calculations with 𝒪(106)\sim\mathcal{O}(10^{6}) determinants and parameters ϵ1,ϵ2107\epsilon_{1},\epsilon_{2}\approx 10^{-7} were executed using a locally modified version of Dice [36, 37, 62].

Recall that our HF calculations utilized a complex plane wave basis {eiGpr}\left\{e^{i\vec{G}_{p}\cdot\vec{r}}\right\} which yields two-electron integral tensors pq|rs\innerproduct{pq}{rs} possessing a fourfold symmetry, i.e.

pq|rs\displaystyle\innerproduct{pq}{rs} =rs|pq=sr|qp=qp|sr.\displaystyle=\innerproduct{rs}{pq}=\innerproduct{sr}{qp}=\innerproduct{qp}{sr}. (37)

Quantum chemistry codes including PySCF, however, often assume the basis sets to be real and hence pq|rs\innerproduct{pq}{rs} to be eightfold symmetric, i.e.

pq|rs\displaystyle\innerproduct{pq}{rs} =rq|ps=ps|rq=rs|pq\displaystyle=\innerproduct{rq}{ps}=\innerproduct{ps}{rq}=\innerproduct{rs}{pq}
=qp|sr\displaystyle=\innerproduct{qp}{sr} =sp|qr=qr|sp=sr|qp.\displaystyle=\innerproduct{sp}{qr}=\innerproduct{qr}{sp}=\innerproduct{sr}{qp}. (38)

Thus, to interface with such codes, we applied a unitary transformation to rotate the plane wave basis into a real-valued basis composed of sines and cosines (up to a normalization factor):

[eiGpreiGpr]\displaystyle\begin{bmatrix}e^{i\vec{G}_{p}\cdot\vec{r}}\\ e^{-i\vec{G}_{p}\cdot\vec{r}}\end{bmatrix} \xlongrightarrow𝐔2[cos((Gpr))sin((Gpr))],\displaystyle\xlongrightarrow{\quad\mathbf{U}\quad}\sqrt{2}\begin{bmatrix}\cos{(\vec{G}_{p}\cdot\vec{r})}\\ \sin{(\vec{G}_{p}\cdot\vec{r})}\end{bmatrix}, (39)
𝐔\displaystyle\mathbf{U} =12[1i1i].\displaystyle=\frac{1}{\sqrt{2}}\begin{bmatrix}1&-i\\ 1&i\end{bmatrix}. (40)

The resulting two-electron integral tensor has eightfold symmetry recovered.

Appendix D FFT of the charge density

Refer to caption
Figure 8: Peaks extracted from Fig. 3. We use the term “primary” to describe peaks with the largest magnitudes, while “secondary” describes any feature of secondary importance. Gm\vec{G}_{m} denotes the primitive moiré reciprocal vector.

Prior to computing the FFT of the charge densities ρ(r)\rho(\vec{r}) in Fig. 2, we subtracted the average density ρ\langle\rho\rangle to remove the k=0\vec{k}=0 Fourier component. To yield better insight into the periodic structure of the charge densities in Fig. 2, we extracted the locations where the FFT signal was maximum (i.e. primary) and where the signal was sub-maximum (i.e. secondary). The results are shown in Fig. 8. At ν=1/7, 1/3, 2/3, 6/7,\nu=1/7,\ 1/3,\ 2/3,\ 6/7, and 11, the FFT exhibits clean peaks that can be easily extracted. At ν=1/4, 2/5, 1/2, 3/5,\nu=1/4,\ 2/5,\ 1/2,\ 3/5, and 3/43/4, however, in addition to the clear peaks located at the moiré reciprocal lattice sites, the FFT displays secondary disperse features in reciprocal space associated with charge delocalization in real space. We estimated the locations of possible secondary peaks suppressed by the disperse features but nonetheless can still be gauged by eye. An example is shown for ν=1/4\nu=1/4 in Fig. 9(a), where we first removed the primary peaks at the moiré reciprocal lattice sites before applying a Gaussian blur filter to smooth out |ρ(k)||\rho(\vec{k})| near the secondary peaks. The locations (marked by the red dots) where |ρ(k)||\rho(\vec{k})| is above some threshold value are then extracted. In cases where there are clusters of similar |ρ(k)||\rho(\vec{k})| above the threshold value, as for ν=3/4\nu=3/4 in Fig. 9(b), we identified the cluster center by kk-means clustering.

Refer to caption
(a) ν=1/4\nu=1/4
Refer to caption
(b) ν=3/4\nu=3/4
Figure 9: Illustration of our FFT secondary peak extraction by: (i) removing the primary peaks at the moiré reciprocal lattice sites; (ii) applying a Gaussian blur filter to smooth out |ρ(k)||\rho(\vec{k})| near the secondary peaks; and (iii) extracting the locations where |ρ(k)||\rho(\vec{k})| is above some threshold value directly or via kk-means clustering.

Appendix E Finite-size effects

For a NN-particle system, the Madelung-corrected total energy per particle ε(N)=E(N)/N\varepsilon(N)=E(N)/N and charge gap Δc(N)=E(N+1)+E(N1)2E(N)\Delta_{c}(N)=E(N+1)+E(N-1)-2E(N) are computed as [34, 32]

ε~(N)\displaystyle\tilde{\varepsilon}(N) =ε(N)+εM,\displaystyle=\varepsilon(N)+\varepsilon_{M}, (41)
Δ~c(N)\displaystyle\tilde{\Delta}_{c}(N) =Δc(N)+2|εM|,\displaystyle=\Delta_{c}(N)+2|\varepsilon_{M}|, (42)

which improves the scaling in the finite-size errors from 𝒪(N1/2)\mathcal{O}(N^{-1/2}) to 𝒪(N1)\mathcal{O}(N^{-1}) for ε\varepsilon [32] and to 𝒪(N3/2)\mathcal{O}(N^{-3/2}) for Δc\Delta_{c} [35]. Our extrapolations of ε~\tilde{\varepsilon} and Δ~c\tilde{\Delta}_{c} to NN\to\infty thus assume 𝒪(N1)\mathcal{O}(N^{-1}) and 𝒪(N3/2)\mathcal{O}(N^{-3/2}) scalings, respectively. As illustrative examples, Fig. 10 shows the convergence of ε~\tilde{\varepsilon} and Δ~c\tilde{\Delta}_{c} with NN at ν=1/3,1/2,3/4,\nu=1/3,1/2,3/4, and 11. It should be noted that the values at ν=1/3\nu=1/3 and 1 were derived from ordered phases, while those at ν=1/2\nu=1/2 and 3/4 were derived from disordered phases, as depicted in Figs. 10a-d. This likely contributes to the more monotonic convergence in ε~\tilde{\varepsilon} and Δ~c\tilde{\Delta}_{c} at ν=1/3\nu=1/3 and 1 compared to ν=1/2\nu=1/2 and 3/4. For the latter case, the greater variation may be associated with disordered configurations that differ slightly from one another. The variation in ε~\tilde{\varepsilon} as a function of NN is also generally smaller than in Δ~c\tilde{\Delta}_{c}.

Furthermore, a manifestation of finite-size effects was observed in the spin-collinearity of solutions initialized from the continuum model guess, as shown in Fig. 11. While only collinear polarized solutions were realized for the largest system size studied at each ν\nu, the smallest system may be non-collinear. This suggests that–at least at the Hartree-Fock level–the spin properties of low-lying states may be sensitive to finite-size effects. However, we note that our initial guess was occasionally found to be insufficient in biasing solutions away from ferromagnetic behavior; lower energy solutions may indeed be non-collinear. As the primary goal of this paper is to examine charge orders, we do not believe any of our main conclusions are affected.

Refer to caption
Figure 10: Convergence of the Madelung-corrected total energy per particle ε~\tilde{\varepsilon} and charge gap Δ~c\tilde{\Delta}_{c} with system size NN. The values at ν=1/3\nu=1/3 and 1 were derived from ordered phases, while those at ν=1/2\nu=1/2 and 3/4 were derived from disordered phases, as depicted in panels (a)-(d).
Refer to caption
Figure 11: The quantities μ\mu and ε0\varepsilon_{0} of GHF solutions initialized from the unpolarized continuum model guess. While the smallest NN-particle system we studied at each ν\nu may not be collinear, the largest NN-particle system is always collinear. This suggests that the emergence of non-collinear solutions could be an artifact of finite-size effects. In the order of increasing ν\nu, the smallest (largest) NN we considered were 7 (28), 4 (81), 3 (108), 10 (40), 8 (162), 15 (60), 6 (150), 12 (147), 42 (168), and 9 (144).

Appendix F Screened Coulomb potential

Refer to caption
Figure 12: The setup for deriving the screened Coulomb potential.

In experiments on TMD moiré systems that we consider, the setup often involves encapsulating the moiré system between dielectric materials and metallic gates (Fig. 12). The encapsulating material and gates alter the electrostatic environment for a charge in the TMD layer, resulting in a screened Coulomb interaction between in-plane charges. Here, we derive a screened Coulomb potential that is a generalization of the double gate-screened [27, 41, 63] and Rytova-Keldysh [28, 29] potentials considered in moiré and TMD materials, respectively.

Regarding the setup in Fig. 12, our task is to find the electric potential φ2(r;r)\varphi_{2}(\vec{r};\vec{r}\mkern 2.0mu\vphantom{r}^{\prime}) at a point r=(x,y,z)\vec{r}=(x,y,z) due to a charge qq located at r=(0,0,z)\vec{r}\mkern 2.0mu\vphantom{r}^{\prime}=(0,0,z^{\prime}), where z,z[l/2,l/2]z,z^{\prime}\in[-l/2,l/2]. Let φi\varphi_{i} and ϵi\epsilon_{i} denote the potential and dielectric constant in layer ii, respectively. Specializing to the case where ϵ1=ϵ3\epsilon_{1}=\epsilon_{3}, the electrostatic equations to solve are

r2φ1(r;r)\displaystyle\nabla_{\vec{r}}^{2}\ \varphi_{1}(\vec{r};\vec{r}\mkern 2.0mu\vphantom{r}^{\prime}) =0,\displaystyle=0, (43)
r2φ2(r;r)\displaystyle\nabla_{\vec{r}}^{2}\ \varphi_{2}(\vec{r};\vec{r}\mkern 2.0mu\vphantom{r}^{\prime}) =4πqϵ2δ(rr),\displaystyle=-\frac{4\pi q}{\epsilon_{2}}\delta(\vec{r}-\vec{r}\mkern 2.0mu\vphantom{r}^{\prime}), (44)
r2φ3(r;r)\displaystyle\nabla_{\vec{r}}^{2}\ \varphi_{3}(\vec{r};\vec{r}\mkern 2.0mu\vphantom{r}^{\prime}) =0,\displaystyle=0, (45)

with BCs (dropping the r\vec{r}\mkern 2.0mu\vphantom{r}^{\prime} parameter for notational clarity)

φ1(ρ,d2)\displaystyle\varphi_{1}\left(\vec{\rho},\frac{d}{2}\right) =0,\displaystyle=0, (46)
φ1(ρ,l2)\displaystyle\varphi_{1}\left(\vec{\rho},\frac{l}{2}\right) =φ2(ρ,l2),\displaystyle=\varphi_{2}\left(\vec{\rho},\frac{l}{2}\right), (47)
φ2(ρ,l2)\displaystyle\varphi_{2}\left(\vec{\rho},-\frac{l}{2}\right) =φ3(ρ,l2),\displaystyle=\varphi_{3}\left(\vec{\rho},-\frac{l}{2}\right), (48)
φ3(ρ,d2)\displaystyle\varphi_{3}\left(\vec{\rho},-\frac{d}{2}\right) =0,\displaystyle=0, (49)
ϵ1φ1z(ρ,l2)\displaystyle\epsilon_{1}\frac{\partial\varphi_{1}}{\partial z}\left(\vec{\rho},\frac{l}{2}\right) =ϵ2φ2z(ρ,l2),\displaystyle=\epsilon_{2}\frac{\partial\varphi_{2}}{\partial z}\left(\vec{\rho},\frac{l}{2}\right), (50)
ϵ2φ2z(ρ,l2)\displaystyle\epsilon_{2}\frac{\partial\varphi_{2}}{\partial z}\left(\vec{\rho},-\frac{l}{2}\right) =ϵ1φ3z(ρ,l2).\displaystyle=\epsilon_{1}\frac{\partial\varphi_{3}}{\partial z}\left(\vec{\rho},-\frac{l}{2}\right). (51)

ρ2\vec{\rho}\in\mathbb{R}^{2} is a vector in the xyxy-plane. Furthermore, we assume dld\ll l so we can treat the layers as possessing infinite extent in the xyxy-plane. Fourier transforming (43)-(45) in ρ\vec{\rho}, we obtain a system of second-order partial differential equations (PDEs)

2z2φ1(k,z)k2φ1(k,z)\displaystyle\frac{\partial^{2}}{\partial z^{2}}\varphi_{1}(\vec{k},z)-k^{2}\varphi_{1}(\vec{k},z) =0,\displaystyle=0, (52)
2z2φ2(k,z)k2φ2(k,z)\displaystyle\frac{\partial^{2}}{\partial z^{2}}\varphi_{2}(\vec{k},z)-k^{2}\varphi_{2}(\vec{k},z) =4πqϵ2δ(zz),\displaystyle=-\frac{4\pi q}{\epsilon_{2}}\delta(z-z^{\prime}), (53)
2z2φ3(k,z)k2φ3(k,z)\displaystyle\frac{\partial^{2}}{\partial z^{2}}\varphi_{3}(\vec{k},z)-k^{2}\varphi_{3}(\vec{k},z) =0.\displaystyle=0. (54)

Let us define the operator =2z2k2\mathcal{L}=\frac{\partial^{2}}{\partial z^{2}}-k^{2}. Since φ2(k,z)\varphi_{2}(\vec{k},z) is the solution to an inhomogeneous PDE with inhomogeneous BCs, we can write it as

φ2(k,z)\displaystyle\varphi_{2}(\vec{k},z) =φ~2(k,z)+G(k,z),\displaystyle=\tilde{\varphi}_{2}(\vec{k},z)+G(\vec{k},z), (55)

where φ~2(k,z)\tilde{\varphi}_{2}(\vec{k},z) is the solution to the homogeneous PDE φ(k,z)=0\mathcal{L}\varphi(\vec{k},z)=0 with inhomogeneous BCs, while G(k,z)G(\vec{k},z) is the solution to the inhomogeneous PDE φ(k,z)=(4πq/ϵ2)δ(zz)\mathcal{L}\varphi(\vec{k},z)=-(4\pi q/\epsilon_{2})\delta(z-z^{\prime}) with homogeneous BCs. The screened potential we obtain from solving the system of PDEs is given by

φ~2(k,z;r)\displaystyle\tilde{\varphi}_{2}(\vec{k},z;\vec{r}\mkern 2.0mu\vphantom{r}^{\prime})
=Ξ{[ek(z+l/2)sinh[k(l2+z)]+ek(zl/2)sinh[k(l2z)]][ϵ2ϵ1(1cosh([k(dl)]))sinh[k(dl)]]\displaystyle=\Xi\Bigg{\{}\left[e^{k(z+l/2)}\sinh{\left[k\left(\frac{l}{2}+z^{\prime}\right)\right]}+e^{-k(z-l/2)}\sinh{\left[k\left(\frac{l}{2}-z^{\prime}\right)\right]}\right]\left[\frac{\epsilon_{2}}{\epsilon_{1}}\left(1-\cosh{[k(d-l)]}\right)-\sinh{[k(d-l)]}\right]
+[ek(zl/2)sinh[k(l2z)]+ek(z+l/2)sinh[k(l2+z)]][ϵ2ϵ1(1cosh([k(dl)]))+sinh[k(dl)]]},\displaystyle\qquad+\left[e^{k(z-l/2)}\sinh{\left[k\left(\frac{l}{2}-z^{\prime}\right)\right]}+e^{-k(z+l/2)}\sinh{\left[k\left(\frac{l}{2}+z^{\prime}\right)\right]}\right]\left[\frac{\epsilon_{2}}{\epsilon_{1}}\left(1-\cosh{[k(d-l)]}\right)+\sinh{[k(d-l)]}\right]\Bigg{\}}, (56)
Ξ\displaystyle\Xi =4πqϵ1ksinh(kl){(ϵ2ϵ1+1)2sinh(kd)(ϵ2ϵ11)2sinh[k(d2l)]2[(ϵ2ϵ1)21]sinh(kl)}1,\displaystyle=-\frac{4\pi q}{\epsilon_{1}k\sinh{(kl)}}\left\{\left(\frac{\epsilon_{2}}{\epsilon_{1}}+1\right)^{2}\sinh{(kd)}-\left(\frac{\epsilon_{2}}{\epsilon_{1}}-1\right)^{2}\sinh{[k(d-2l)]}-2\left[\left(\frac{\epsilon_{2}}{\epsilon_{1}}\right)^{2}-1\right]\sinh{(kl)}\right\}^{-1}, (57)
G(k,z;r)\displaystyle G(\vec{k},z;\vec{r}\mkern 2.0mu\vphantom{r}^{\prime}) =2πqϵ2ksinh(kl){cosh([k(z+z)])cosh([k(|zz|l)])}.\displaystyle=-\frac{2\pi q}{\epsilon_{2}k\sinh{(kl)}}\left\{\cosh{[k(z+z^{\prime})]}-\cosh{[k\left(|z-z^{\prime}|-l\right)]}\right\}. (58)

F.1 Double gate-screened potential

In the limit ϵ1=ϵ2\epsilon_{1}=\epsilon_{2}, the BCs match those of the double gate-screened potential and so we recover this solution. Setting ϵ1=ϵ2\epsilon_{1}=\epsilon_{2} in (56)-(58),

φ~2(k,z;r)\displaystyle\tilde{\varphi}_{2}(\vec{k},z;\vec{r}\mkern 2.0mu\vphantom{r}^{\prime})
=Ξ{[ek(z+l/2)sinh[k(l2+z)]+ek(zl/2)sinh[k(l2z)]][1cosh([k(dl)])sinh[k(dl)]]\displaystyle=\Xi\Bigg{\{}\left[e^{k(z+l/2)}\sinh{\left[k\left(\frac{l}{2}+z^{\prime}\right)\right]}+e^{-k(z-l/2)}\sinh{\left[k\left(\frac{l}{2}-z^{\prime}\right)\right]}\right]\left[1-\cosh{[k(d-l)]}-\sinh{[k(d-l)]}\right]
+[ek(zl/2)sinh[k(l2z)]+ek(z+l/2)sinh[k(l2+z)]][1cosh([k(dl)])+sinh[k(dl)]]},\displaystyle\qquad+\left[e^{k(z-l/2)}\sinh{\left[k\left(\frac{l}{2}-z^{\prime}\right)\right]}+e^{-k(z+l/2)}\sinh{\left[k\left(\frac{l}{2}+z^{\prime}\right)\right]}\right]\left[1-\cosh{[k(d-l)]}+\sinh{[k(d-l)]}\right]\Bigg{\}}, (59)
Ξ\displaystyle\Xi =4πqϵksinh(kl){4sinh(kd)}1,\displaystyle=-\frac{4\pi q}{\epsilon k\sinh{(kl)}}\left\{4\sinh{(kd)}\right\}^{-1}, (60)
G(k,z;r)\displaystyle G(\vec{k},z;\vec{r}\mkern 2.0mu\vphantom{r}^{\prime}) =2πqϵksinh(kl){cosh([k(z+z)])cosh([k(|zz|l)])}.\displaystyle=-\frac{2\pi q}{\epsilon k\sinh{(kl)}}\left\{\cosh{[k(z+z^{\prime})]}-\cosh{[k\left(|z-z^{\prime}|-l\right)]}\right\}. (61)

Since Coulomb interactions in the GWC states occur over the length scale of the moiré lattice constant LmlL_{m}\gg l, we are interested only in the asymptotic behaviour of φ2(r;r)\varphi_{2}(\vec{r};\vec{r}\mkern 2.0mu\vphantom{r}^{\prime}) for distances ρl\rho\gg l. Using the fact that k1/ρ1/lk\propto 1/\rho\ll 1/l in this case, it is sufficient for us to know the Fourier component φ2(k,z;r)\varphi_{2}(\vec{k},z;\vec{r}\mkern 2.0mu\vphantom{r}^{\prime}) for kl1kl\ll 1. It then follows that k|z|,k|z|1k|z|,k|z^{\prime}|\ll 1 and k|zz|0k|z-z^{\prime}|\approx 0 since l/2z,zl/2-l/2\leq z,z^{\prime}\leq l/2. Taylor expanding φ~2\tilde{\varphi}_{2} and GG to first order in klkl and making the approximation zzz\approx z^{\prime} (note that we still have k|z|𝒪(kl)k|z|\sim\mathcal{O}(kl)), we obtain

φ~2(k,z)\displaystyle\tilde{\varphi}_{2}(\vec{k},z) Ξ{2kl(1cosh(kd))+(kl)2sinhkd4(kz)2sinhkd}\displaystyle\approx\Xi\left\{2kl(1-\cosh{kd})+(kl)^{2}\sinh{kd}-4(kz)^{2}\sinh{kd}\right\}
=πqϵksinhkd{2(1cosh(kd))+klsinhkd4(kz)2klsinhkd},\displaystyle=-\frac{\pi q}{\epsilon k\sinh{kd}}\left\{2(1-\cosh{kd})+kl\sinh{kd}-4\frac{(kz)^{2}}{kl}\sinh{kd}\right\}, (62)
G(k,z)\displaystyle G(\vec{k},z) 2πqϵk[2(kz)2kl12kl],\displaystyle\approx-\frac{2\pi q}{\epsilon k}\left[2\frac{(kz)^{2}}{kl}-\frac{1}{2}kl\right], (63)
φ2(k,z)\displaystyle\implies\varphi_{2}(\vec{k},z) πqϵksinhkd{2(1cosh(kd))+klsinhkd4(kz)2klsinhkd+4(kz)2klsinhkdklsinhkd}\displaystyle\approx-\frac{\pi q}{\epsilon k\sinh{kd}}\left\{2(1-\cosh{kd})+\cancel{kl\sinh{kd}}-\cancel{4\frac{(kz)^{2}}{kl}\sinh{kd}}+\cancel{4\frac{(kz)^{2}}{kl}\sinh{kd}}-\cancel{kl\sinh{kd}}\right\}
=2πqϵktanh((kd2)).\displaystyle=\frac{2\pi q}{\epsilon k}\cdot\tanh{\left(\frac{kd}{2}\right)}. (64)

F.2 Rytova-Keldysh potential

In the limit dd\to\infty, the BCs are similar to the setup in [28] and we recover Eq. (5) therein. We first shift our coordinates zzl/2z\to z-l/2 and zzl/2z^{\prime}\to z^{\prime}-l/2 in (56)-(58) to match the coordinate system of [28]:

φ~2(k,z;r)\displaystyle\tilde{\varphi}_{2}(\vec{k},z;\vec{r}\mkern 2.0mu\vphantom{r}^{\prime}) =Ξ{[ekzsinh(kz)ek(zl)sinh[k(zl)]][ϵ2ϵ1(1cosh([k(dl)]))sinh[k(dl)]]\displaystyle=\Xi\Bigg{\{}\left[e^{kz}\sinh{(kz^{\prime})}-e^{-k(z-l)}\sinh{\left[k(z^{\prime}-l)\right]}\right]\left[\frac{\epsilon_{2}}{\epsilon_{1}}\left(1-\cosh{[k(d-l)]}\right)-\sinh{[k(d-l)]}\right]
+[ekzsinh(kz)ek(zl)sinh[k(zl)]][ϵ2ϵ1(1cosh([k(dl)]))+sinh[k(dl)]]},\displaystyle\qquad+\left[e^{-kz}\sinh{(kz^{\prime})}-e^{k(z-l)}\sinh{\left[k(z^{\prime}-l)\right]}\right]\left[\frac{\epsilon_{2}}{\epsilon_{1}}\left(1-\cosh{[k(d-l)]}\right)+\sinh{[k(d-l)]}\right]\Bigg{\}}, (65)
Ξ\displaystyle\Xi =4πqϵ1ksinh(kl){(ϵ2ϵ1+1)2sinh(kd)(ϵ2ϵ11)2sinh[k(d2l)]2[(ϵ2ϵ1)21]sinh(kl)}1,\displaystyle=-\frac{4\pi q}{\epsilon_{1}k\sinh{(kl)}}\left\{\left(\frac{\epsilon_{2}}{\epsilon_{1}}+1\right)^{2}\sinh{(kd)}-\left(\frac{\epsilon_{2}}{\epsilon_{1}}-1\right)^{2}\sinh{[k(d-2l)]}-2\left[\left(\frac{\epsilon_{2}}{\epsilon_{1}}\right)^{2}-1\right]\sinh{(kl)}\right\}^{-1}, (66)
G(k,z;r)\displaystyle G(\vec{k},z;\vec{r}\mkern 2.0mu\vphantom{r}^{\prime}) =2πqϵ2ksinh(kl){cosh([k(z+zl)])cosh([k(|zz|l)])}.\displaystyle=-\frac{2\pi q}{\epsilon_{2}k\sinh{(kl)}}\left\{\cosh{[k(z+z^{\prime}-l)]}-\cosh{[k\left(|z-z^{\prime}|-l\right)]}\right\}. (67)

Applying the hyperbolic relations

cosh([k(dl)])\displaystyle\cosh{[k(d-l)]} =cosh((kd))cosh((kl))sinh(kd)sinh(kl),\displaystyle=\cosh{(kd)}\cosh{(kl)}-\sinh{(kd)}\sinh{(kl)}, (68)
sinh[k(dl)]\displaystyle\sinh{[k(d-l)]} =sinh(kd)cosh((kl))cosh((kd))sinh(kl),\displaystyle=\sinh{(kd)}\cosh{(kl)}-\cosh{(kd)}\sinh{(kl)}, (69)
sinh[k(d2l)]\displaystyle\sinh{[k(d-2l)]} =sinh(kd)cosh((2kl))cosh((kd))sinh(2kl),\displaystyle=\sinh{(kd)}\cosh{(2kl)}-\cosh{(kd)}\sinh{(2kl)}, (70)

and rewriting in terms of

δ\displaystyle\delta =ϵ2/ϵ11ϵ2/ϵ1+1,\displaystyle=\frac{\epsilon_{2}/\epsilon_{1}-1}{\epsilon_{2}/\epsilon_{1}+1}, (71)

we have upon taking dd\to\infty:

φ~2(k,z;r)\displaystyle\tilde{\varphi}_{2}(\vec{k},z;\vec{r}\mkern 2.0mu\vphantom{r}^{\prime}) 2πqϵ2ksinh(kl)1+δe2klδ2ekl{eklcosh([k(z+zl)])cosh([k(zz)])\displaystyle\approx\frac{2\pi q}{\epsilon_{2}k\sinh{(kl)}}\cdot\frac{1+\delta}{e^{2kl}-\delta^{2}}\cdot e^{kl}\Big{\{}e^{kl}\cosh{[k(z+z^{\prime}-l)]}-\cosh{[k(z-z^{\prime})]}
+δ[cosh(k(zz))]eklcosh([k(z+zl)])},\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad+\delta[\cosh{k(z-z^{\prime})}]-e^{-kl}\cosh{[k(z+z^{\prime}-l)]}\Big{\}}, (72)
φ2(k,z;r)\displaystyle\implies\varphi_{2}(\vec{k},z;\vec{r}\mkern 2.0mu\vphantom{r}^{\prime}) =φ~2(k,z;r)+G(k,z;r)\displaystyle=\tilde{\varphi}_{2}(\vec{k},z;\vec{r}\mkern 2.0mu\vphantom{r}^{\prime})+G(\vec{k},z;\vec{r}\mkern 2.0mu\vphantom{r}^{\prime})
2πqϵ2k{ek|zz|+2δe2klδ2[δcosh([k(zz)])+eklcosh[k(z+zl)]]}.\displaystyle\approx\frac{2\pi q}{\epsilon_{2}k}\left\{e^{-k|z-z^{\prime}|}+\frac{2\delta}{e^{2kl}-\delta^{2}}\left[\delta\cosh{[k(z-z^{\prime})]}+e^{kl}\cosh[k(z+z^{\prime}-l)]\right]\right\}. (73)