Melting of generalized Wigner crystals in transition metal dichalcogenide heterobilayer Moiré systems

The promise of a highly tunable lattice system that can allow solid-state-based simulation of strong coupling physics Andrei et al. (2021) has largely driven the explosion of efforts studying Moiré superlattices. The transition metal dichalcogenide (TMD) heterobilayer Moiré systems with zero twist-angle (see Fig. 1(a)) and localized Wannier orbitals form a uniquely simple platform to explore phases driven by strong interactions Tang et al. (2020); Regan et al. (2020). Indeed, upon sweeping the density of electrons per Moiré unit cell, incompressible charge ordered states have been observed at various fractional fillings Regan et al. (2020); Xu et al. (2020); Li et al. (2021a). These charge orders can be viewed as a generalized Wigner crystalline state that reduces the symmetry of the underlying Moiré lattice, as they are driven by the long-range Coulomb interaction. The density controlled melting of Wigner crystals is expected to result in a rich hierarchy of intermediate phases Spivak and Kivelson (2004, 2006); Jamei et al. (2005). While a microscopic theoretical study of Wigner crystal melting is challenging due to the continuous spatial symmetry, the melting of generalized Wigner crystals is more amenable to a microscopic study due to the lattice. The observation of intermediate compressible states with optical anisotropy Jin et al. (2021) (see Fig. 1(b)) and the tunability beyond density Li et al. (2021b); Ghiotto et al. (2021) present a tantalizing possibility to study melting and the possible intermediate phases of the generalized Wigner crystals.

The underlying lattice in the generalized Wigner crystal reduces the continuous rotational symmetry to $D_{3}$ point group symmetry. Ref. Coppersmith et al. (1982) studied the melting of a $1/3$-filled crystalline state on a triangular lattice in the context of Krypton adsorbed on Graphene. Based on the free energy costs of the domain walls and domain wall intersections, they reasoned that the generalized WC would first melt into a hexagonal liquid, and then crystallize into a stripe solid. From the modern perspectives of electronic liquid crystals Kivelson et al. (1998), one anticipates nematic fluid states in the vicinity of crystalline anisotropic states such as the stripe solid. Moreover, the $D_{3}$ point group symmetry of the triangular lattice relevant for hetero-TMD Moiré sytems further enriches the possibilities of the intermediate fluid phases. The triangular lattice admits two types of nematic states due to the nematic order parameter transforming as a 2-dimensional irreducible representation of the lattice point group Serre (2012); Fernandes and Venderbos (2020). The hetero-TMD Moiré systems present an excellent opportunity to study these intermediate liquid phases.

As the quantum melting of charge order is a notoriously difficult problem  Kivelson et al. (2003), in this paper we will take advantage of the small bandwidth in hetero-TMD and use a strong coupling approach. We study Monte Carlo simulations inspired by hetero-TMD Moiré systems at and between commensurate charge ordered states. We analyze our results in terms of the structure factor and a nematic order parameter correlation function. In particular, we distinguish between the two possible types of nematic states illustrated in Fig. 1(c), one associated with the director aligned with a single majority bond orientation (which we dub type-I), and the other perpendicular to a single minority bond orientation (type-II). As shown in Fig. 1(d), we find the type-I nematic and type-II nematic to each robustly appear between 2/5 and 1/2 and 1/3 and 2/5 respectively. We conclude with a discussion.

As the orientation of the nematic director is defined within the angle range $\theta\in[0,\pi)$ (Fig. 1(c)) we define the local nematic field using complex notation $N(\vec{r})=n(\vec{r})e^{i2\theta(\vec{r})}$, where $n(\vec{r})\in\mathbb{R}$. In terms of this nematic order parameter field, the free energy density describing the isotropic-nematic transition in a trigonal system takes the following form Fernandes and Venderbos (2020); Venderbos and Fernandes (2018); Hecker and Schmalian (2018); Little et al. (2020):

The sign of the coefficient of the cubic term, $\gamma$, determines whether the system wants to be in a type-I nematic state ($\langle\cos(6\theta)\rangle=+1$) or type-II nematic state ($\langle\cos(6\theta)\rangle=-1$).

We explore the phase diagram using classical Monte Carlo as a function of temperature $T$ and the number of particles per Moiré site $\nu$. To emulate the experimental setup in refs. Xu et al. (2020); Jin et al. (2021); Li et al. (2021c); Regan et al. (2020), the Hamiltonian that we simulate describes the Coulomb interaction for electrons halfway between two dielectric gates a distance $d$ apart with dielectric constant $\epsilon$:

Here, $K_{0}$ is the modified Bessel function of the second kind, $a$ is the Moiré lattice constant, and $\rho_{i,j}\in\{0,1\}$ are the occupancies of lattice sites $i,j$. As in refs. Xu et al. (2020); Jin et al. (2021), we take $a=8nm$ and $d=10a$, and we take $e^{2}/(4\pi\epsilon\epsilon_{0}a)$ as our unit of energy for simulation.

Because the interaction is long-ranged, simply simulating a system with periodic boundary conditions would result in ambiguous distance calculations. Thus, we simulate a formally infinite system that is constrained to be periodic in an $\ell\times\ell$ rhombus. Particles interact both within and between copies of the system. The choice of an $\ell\times\ell$ rhombus has the full symmetry of the triangular lattice as, when one considers the infinite system, the action of each element of the point group is a bijective map on the set of unique sites contained within the simulation cell. Thus we do not expect our choice of geometry to artificially promote rotational symmetry breaking. Moreover in each nematic state that we report, our simulations find configurations with each of the three possible director orientations for the relevant nematic type with equal probability. For Monte Carlo updates, we use arbitrary-range, single particle occupancy exchanges with standard Metropolis acceptance rules. However, the prevalence of short-range correlated structures leading to long autocorrelation times complicates our simulations, especially in the incompressible density region. To deal with this, we developed a cluster algorithm in the spirit of the well-known Wolff algorithm Wolff (1989) and the later geometric cluster algorithm Heringa and Blöte (1998). For more detail, see appendix A. We perform a cluster update after every 1000 single particle occupancy exchange updates.

At each point in phase space, we calculate the Monte Carlo average of the structure factor

to assess crystalline order. To assess the degree of rotational symmetry breaking, we also calculate the average of the nematic order parameter correlation function (per site) given by

where $\vec{q}$ denotes Fourier momentum. At high temperatures, when $\langle N(\vec{r})\rangle=0$, $\langle\tilde{C}(\vec{q}=0)\rangle/\ell^{2}$ behaves as $k_{B}T$ times the nematic susceptibility: $\chi(\vec{q}=0)k_{B}T$. Generically we expect this to have some continuous behavior as a function of temperature. However, when the order parameter develops an expectation value in a nematic state, $\langle\tilde{C}(\vec{q}=0)\rangle/\ell^{2}$ should acquire a constant, non-zero value. To determine the type of nematicity exhibited by nematic states, we also calculate $\langle\cos(6\theta)\rangle$, where, as in Fig. 1(c), type-I (type-II) nematic states have $\langle\cos(6\theta)\rangle=+1\ (-1)$. For further details about the calculation of these quantities from our Monte Carlo simulation data, see appendix B. All results that we show are obtained from an $\ell=20$ system, except for exactly at $\nu=1/3$ since $20\times 20/3$ is not an integer. In all cases, we perform $10^{5}$ updates per site for equilibration at each temperature, and then $2\times 10^{5}$ updates per site for data collection.

At $\nu=1/3$, we find the isotropic generalized Wigner crystalline phase, shown in Fig. 2(a) for $\ell=12$. This phase has lattice vectors $\vec{a}^{\text{wc}}_{1}=(0,\sqrt{3})$ and $\vec{a}^{\text{wc}}_{2}=(3/2,\sqrt{3}/2)$ as indicated by the black arrows in Fig. 2(a). The structure factor shows well defined peaks at the reciprocal lattice vectors $\vec{G}^{\text{wc}}_{1}=(-2\pi/3,2\pi/\sqrt{3}))$ and $\vec{G}^{\text{wc}}_{2}=(4\pi/3,0)$ associated with the crystalline state (Fig. 2(b)). Upon increasing the density, this crystalline state starts to melt, but it maintains an isotropic, compressible state to a certain filling. At small fillings away from the $1/3$-state, as shown in Fig. 2(c), the extra particles form domain walls between the three different registries of the generalized WC state. Three domain walls are marked with black lines in Fig. 2(c). The domain walls meet at $2\pi/3$ angles, reminiscent of what was found in ref. Coppersmith et al. (1982).

As in ref. Coppersmith et al. (1982), this domain wall structure is stable while the density of domain walls is dilute (and hence the domain walls are long) because energetics favor $2\pi/3$ angles. Taking into account interactions up to fifth neighbor, the energy contributed to the Hamiltonian by the six particles in the three dimers composing the $2\pi/3$ vertex is

while that of the particles composing three straight domain wall dimers is

where $V_{i}$ denotes the energy of two $i$’th neighbor particles. Using eq.( 2) it is easy to check that $E_{DW}-E_{v}>0$, and thus (at least to this order of interaction), it is energetically favorable to have $2\pi/3$-vertices, even at low temperatures. However, the densest possible hexagonal domain state consists of a close packing of $2\pi/3$-vertices, which has density $\nu=3/8$. Thus, this state certainly cannot exist at densities $\nu>3/8$. In Fig. 2(d), the structure factor of this compressible state is still dominated by the generalized WC, but with the domain size becoming finite, the superstructure peaks are broadened. Looking at the ensemble averaged nematic correlation function $\langle\tilde{C}(\vec{q}=0)\rangle$ in Fig. 2(e) confirms that this state is isotropic as it drops to zero below the transition.

The state is dramatically different at $\nu=1/2$. We have the charge stripe state shown in Fig. 3(a) with lattice vectors $\vec{a}^{\text{cs}}_{1}=(1,0)$ and $\vec{a}^{\text{cs}}_{2}=(0,\sqrt{3})$. There are two degenerate charge stripe states whose lattice vectors are related by $\pi/3$ and $2\pi/3$ rotations of $\vec{a}^{\text{cs}}_{1,2}$. The structure factor in Fig. 3(b), averaged over configurations with the same orientation as the one shown in Fig. 3(a), contains peaks at the reciprocal lattice vectors $\vec{G}^{\text{cs}}_{1}=(2\pi,0)$ and $\vec{G}^{\text{cs}}_{2}=(0,2\pi/\sqrt{3})$. As expected, $\vec{G}^{\text{cs}}_{i}\cdot\vec{a}^{\text{cs}}_{j}=2\pi\delta_{ij}$. Diluting the $1/2$-filled state, the stripes become shorter via the introduction of dislocations, as shown in Fig 3(c). The structure factor reflects the finite length of these stripe domains in the splitting of the stripe peak over the span of the stripe domain size scale. The peak at $\vec{a}^{\text{cs}}_{2}$ is split into two peaks separated by $2\pi/L_{N}$ where $L_{N}\approx 4.63$ is the average stripe domain size. The nematic correlation function reveals that, unlike the isotropic phases in Fig. 2(e), $\langle\tilde{C}(\vec{q}=0)\rangle$ shows a sharp jump at $T_{c}$ to a finite value in Fig. 3(e). This indicates that these are nematic states. By examining $\langle\cos(6\theta)\rangle$ in Fig. 3(f), we can see that both of these phases are type-I nematic.

Upon further diluting, the system maintains the same type of anisotropy and forms the columnar dimer crystal state at $\nu=2/5$, shown in Fig. 4(a) with director orientation $\theta=0$. This is the limit of the shortest stripe length, evolving from $\nu=1/2$. The $\nu=2/5$ state is a crystalline state with lattice vectors $\vec{a}^{\text{cdc}}_{1}=(0,\sqrt{3})$ and $\vec{a}^{\text{cdc}}_{2}=(5/2,\sqrt{3}/2)$. We mark the reciprocal lattice vectors $\vec{G}^{\text{cdc}}_{1}=(-2\pi/5,2\pi/\sqrt{3})$ and $\vec{G}^{\text{cdc}}_{2}=(4\pi/5,0)$ in the structure factor shown in Fig. 4(b). Note that the peak at $2\vec{G}^{\text{cdc}}_{2}$ is more intense than the one at $\vec{G}^{\text{cdc}}_{2}$. This is due to the form factor from the lattice basis. As we dilute further, the length of the columns get shorter as dimers get broken up. At lower densities, the columns do not extend over the entire system, so there are finite length segments of columns that can have different orientations as illustrated in Fig. 4(c). The broken pieces of dimers form short-range correlated domains of generalized WC. This is shown by the broad peaks in the structure factor in Fig. 4(d). This compressible state no longer has the mirror symmetries of the columnar state. It is still anisotropic as we can see from the nematic correlation function in Fig. 4(e). Interestingly, the columnar fragments intersect at $2\pi/3$ angles, as well as $\pi/3$ angles, one of which is circled in red in fig. 4(c). While the $2\pi/3$ intersections are isotropic, the $\pi/3$ intersections consist primarily of only two of the three possible nearest-neighbor bond orientations, and hence this state is a type-II nematic phase. We confirm this by observing that $\langle\cos(6\theta)\rangle=-1$ at low temperatures in Fig. 4(f). Thus we predict microscopic mechanism for the type-II nematic phase.

The nematic-II state in the region of $3/8<\nu<2/5$ is supported energetically. Upon increasing density beyond $\nu>3/8$, columnar fragments have to either intersect also at $\pi/3$ or be parallel to each other. Since the distance between columnar fragments increases away from the $\pi/3$ intersection, we expect $\pi/3$ intersections to be favored. See appendix C for a schematic calculation demonstrating this. Such $\pi/3$ intersections involve two nearest-neighbor bond orientations, promoting a nematic-II state.

One could experimentally probe our predicted nematic states by performing optical measurements similar to those done at $\nu=1/2$ in ref. Jin et al. (2021). As one lowers the density from $\nu=1/2$ to $\nu=1/3$ we would anticipate a rotation of the nematic director and consequently a shift in the peaks of the measured optical anisotropy axis. In particular, as one decreases the density from between $\nu=1/2$ and $\nu=2/5$, we predict that the measured anisotropy axis should have peaks along the nematic-I orientations $0,\pi/3,2\pi/3$. Below $\nu=2/5$, when the director rotates into the nematic-II state, the peaks should be at $\pi/6,\pi/2,5\pi/6$. Finally, below $\nu=3/8$ when the system becomes isotropic, we expect that there should be no preferred anisotropy axis at all. For $1/3<\nu<3/8$, one could also look for signatures of the hexagonal WC domain state using Umklapp spectroscopy experiments like those done in ref. Shimazaki et al. (2021). The short-range correlated nature of this state should show up as broadened Umklapp resonances around the $\nu=1/3$ generalized Wigner crystal lattice vectors.

In summary, we studied the electronic states of a system of strongly correlated electrons on a triangular lattice in the region $1/3\leq\nu\leq 1/2$ particles per Moiré site. At $\nu=1/2$, we find the charge stripe state. Upon dilution, the charge stripe state melts into a nematic-I short-ranged charge stripe state via the introduction of dislocations. Once the stripes become short enough, the columnar dimer crystal state emerges at $\nu=2/5$. At even lower densities, the remaining columnar fragments space themselves out to lower their energy by intersecting at $\pi/3$ and $2\pi/3$ angles, resulting in a nematic-II state. Below $\nu=3/8$, the system can again lower its energy by using only $2\pi/3$ columnar fragment intersections to form an isotropic, hexagonal network of domain walls between regions of the $\nu=1/3$ generalized Wigner crystal. Finally, at $\nu=1/3$, the pure, isotropic generalized Wigner crystal state emerges. Our intermediate states were not only promoted by entropy, but we also found them to have lower energy compared to macroscopically phase separated states. Accordingly, we suspect that our proposed states are relevant for finite experimental temperatures where fluctuations due to entropy also play a role. We leave the determination of the classical ground state at $T=0$ as a subject of future work.

Studying the intermediate phases of melted density waves has been of interest since considerations of Krypton adsorbed on graphene Coppersmith et al. (1982). However, limitations in computational resources and experimental methods caused difficulties in probing the intermediate states. With advances in computing power and the advent of the TMD Moiré platform, however, detailed phase diagrams can now be predicted computationally and probed experimentally. Our work demonstrates this capacity to explore intermediate phases and the richness of the phase diagram one can obtain with a classical model, even without considering quantum effects. We found the striped phase predicted upon increasing density in ref. Coppersmith et al. (1982) refines into two distinct stripe crystal states neighboring two distinct types of nematics. In particular, we presented a microscopic mechanism for the formation of the nematic-II state via $\pi/3$ intersections between columnar fragments. As a subject of future work, it would be interesting to study the implications of our findings for the melting of WCs without a lattice potential, such as those recently observed in Zhou et al. (2021); Smoleński et al. (2021).

Acknowledgements: The authors acknowledge support by the NSF [Platform for the Accelerated Realization, Analysis, and Discovery of Interface Materials (PARADIM)] under cooperative agreement no. DMR-U638986. We would also like to thank Steven Kivelson, Kin-Fai Mak, Jie Shan, Sue Coppersmith, Ataç Imamoğlu, and Samuel Lederer for helpful discussions.