Charge density waves and their transitions in anisotropic quantum Hall systems

Abstract

In recent experiments, external anisotropy has been a useful tool to tune different phases and study their competitions. In this paper, we look at the quantum Hall charge density wave states in the N = 2 Landau level. Without anisotropy, there are two first-order phase transitions between the Wigner crystal, the 2-electron bubble phase, and the stripe phase. By adding mass anisotropy, our analytical and numerical studies show that the 2-electron bubble phase disappears and the stripe phase significantly enlarges its domain in the phase diagram. Meanwhile, a regime of stripe crystals that may be observed experimentally is unveiled after the bubble phase gets out. Upon increase of the anisotropy, the energy of the phases at the transitions becomes progressively smooth as a function of the filling. We conclude that all first-order phase transitions are replaced by continuous phase transitions, providing a possible realisation of continuous quantum crystalline phase transitions.

Introduction

Quantum Hall (QH) systems play an important role in understanding correlated phenomena. Because of the Landau level (LL) quantisation, the interaction dominates over the kinetic energy when the ratio ν = N_e/N_ϕ between the electron number N_e and number of states N_ϕ inside an LL is not an integer. Various correlated phases appear depending on the LL index N, such as topological QH liquids^1,2 in the lowest LL (N = 0), non-Abelian QH states^3,4,5 and QH nematics^6,7 in the N = 1 LL. Higher LLs with N > 1 allow for the existence of charge density wave (CDW) states with ordered modulation in electron density^8,9,10. Like the QH liquids, these CDW orders emerge from the inherent strong interaction¹¹.

Recently, interest has focused on QH states perturbed by anisotropy. Anisotropy breaks the rotation symmetry of the system and changes its geometry. It is interesting to investigate how different QH phases, e.g., gapped QH fluid^{12,13,14,15,16,17,18,19} or gapless composite fermion liquid states^20,21,22,23, can be tuned through external anisotropy. These studies greatly enhance the understanding of topological robustness against geometric perturbation. Meanwhile, the reaction of CDW states to external anisotropy has been much less studied.

The CDW instability leads to Wigner crystals (WC), bubble phases, and stripe phases^8,24,25. In experiments, CDW phases have different transport properties as compared to the QH fluid phases. The WC and bubble phases are indeed insulating because they are collectively pinned by disorder and thus do not contribute to the Hall conductivity. This is e.g., at the origin of the re-entrant integer QH effect^25,26, which has also been predicted²⁷ and observed in higher LLs of graphene²⁸. The stripe phase strongly breaks the ‘rotational’ symmetry and exhibits a large anisotropy in the DC diagonal resistance. Meanwhile for the N = 2 LL, no QH liquid phase has been observed experimentally so far²⁹ except for the possible $\tilde{\nu }=1/5$ and $\tilde{\nu }=4/5$ states³⁰ at intermediate temperatures, where $\tilde{\nu }$ is the partial filling factor in the N = 2 LL. (In conventional thin GaAs quantum wells, the filling factor is $\nu ={N}_{{\rm{e}}}/{N}_{\phi }=4+\tilde{\nu }$.) The fewer and clearer candidates for ground states in the N ≥ 2 LLs make the study of their competitions in the presence of anisotropy more feasible and reliable.

Because of the strongly interacting nature, QH systems often suffer from the limited availability of theoretical tools and in many cases, it is necessary to resort to numerical calculations. However, CDW phases, unlike the correlated liquid phases, are easily captured by an analytic Hartree–Fock (HF) approximation^31,32. The validity of the HF approximation improves as N becomes larger^9,33 and it also turns out to be capable of catching the essential physics for intermediate N^10,34 confirmed by experiments^26,35,36. Meanwhile, numerical tools always serve as an important reference in QH problems. In the isotropic case, they have turned out to be feasible to exhibit CDW phases. The exact diagonalization (ED)³⁷ and the density matrix renormalization group (DMRG)^38,39,40 reach a good qualitative agreement with the HF approximation as well as experiments for isotropic systems. Therefore we can use both theoretical and numerical calculations to study how higher-LL QH systems react to anisotropy.

In this paper, we provide a quantitative comparison between analytical HF and numerical DMRG calculations to study the CDW phases of spin-polarised electrons in the N = 2 LL under a mass anisotropy, which can be realised in a 2D electron gas in AlAs quantum wells with a mass anisotropy m_x/m_y ≈ 5⁴¹. A tunable mass anisotropy of a 2D electron gas can also be realised by strain⁴² or moiré pattern⁴³. The HF approximation yields a reliable picture for the appearance of different CDW orders while the DMRG calculation additionally provides quantum fluctuations beyond the mean-field ansatz. The predictions of the two methods reach good agreement. We find that the 2-electron (2e) bubble phase is suppressed by increasing the mass ratio between two orthogonal directions. As a result, the unidirectional stripe phase near half filling and the low-filling WC become adjacent in the phase diagram at intermediate fillings. In the isotropic case, previous studies^{44,45,46,47,48,49} suggest that when going from half to intermediate fillings, the unidirectional stripe phase should behave as a stripe crystal, a highly anisotropic WC with the same transverse period as the unidirectional stripe phase. However, such a stripe crystal is usually covered by the triangular 2e bubble phase. When anisotropy enters the system, our result shows that this stripe crystal regime naturally appears and dominates over other CDW states in intermediate fillings. Its density profile can be directly reflected by our DMRG studies. Our analysis of the CDW periodicity shows that the transition from the WC at low fillings and the stripe crystal should be continuous and at least second order for sufficiently large anisotropy. All first-order transitions, found in the isotropic case, are replaced by continuous ones among the WC regime and the stripe regime.

Results

Model and relevant phases

We consider electrons with anisotropic mass m_y/m_x = α² ≠ 1 in a uniform magnetic field, with isotropic Coulomb repulsive interactions. We restrict the electrons to a single partially filled LL, such that the kinetic energy is quenched. The electron-electron interaction, projected to the N^th LL is (see Methods for derivation)

$$\bar{H}=\frac{1}{2A}\mathop{\sum}\limits_{{\bf{q}}}{V}_{\text{eff}}({\bf{q}})\bar{\rho }({\bf{q}})\bar{\rho }(-{\bf{q}});$$

(1)

$${V}_{{\rm{eff}}}({\bf{q}})={F}_{N}^{2}(\alpha {q}_{x}^{2},{q}_{y}^{2}/\alpha )\frac{2\pi {e}^{2}}{\epsilon \sqrt{{q}_{x}^{2}+{q}_{y}^{2}}},$$

(2)

where A is total area of the system. The projected density operators $\bar{\rho }$ consist only of components in the N^th LL and the form factor F_N keeps track of the associated LL wave functions,

$${F}_{N}({q}_{x}^{2},{q}_{y}^{2})={L}_{N}\left(\frac{{q}_{x}^{2}+{q}_{y}^{2}}{2}{l}^{2}\right){e}^{-({q}_{x}^{2}{l}^{2}+{q}_{y}^{2}{l}^{2})/4},$$

(3)

$$\bar{\rho }({\bf{q}})=\mathop{\sum}\limits_{i}{e}^{-i{\bf{q}}\cdot {{\bf{R}}}_{i}},$$

(4)

where R_i is the guiding-centre coordinate of the electron i and L_N is the Laguerre polynomial. The magnetic length is given as $l=\sqrt{\hslash /eB}$. Notice that the mass anisotropy affects the effective interaction only through the form factor. Our HF and DMRG⁵⁰ calculations are based on the Hamiltonian equation (1).

Let us now briefly review the relevant CDW phases in the N = 2 LL for isotropic interactions before presenting our own work. In the dilute limit $\tilde{\nu }\to 0$ (but away from the $\tilde{\nu }=0$ integer QH plateau), the WC is likely to form^10,38, where the electrons have enough space to stay away from each other due to the Coulomb repulsion. For an isotropic effective interaction, the lattice takes a triangular form, which has a maximal crystalline rotational symmetry⁵¹. As the density is increased, the electrons are squeezed. They tend to cluster around an interaction range where the pure Coulomb repulsion is weakened by the LL projection (see Fig. 1 in ref. ¹⁰). A bubble phase with two electrons at each lattice site is formed^10,37. The 2e bubble phase still lies on a triangular lattice but has a different number of electrons in each unit cell, leading to a discontinuity in the derivative of the energy per particle ${E}_{{\rm{pp}}}(\tilde{\nu })$ which will be elaborated later. Around this first-order transition, a mixed phase that consists of a WC coexisting with the 2e bubble phase can form^10,36. When the filling factor further approaches $\tilde{\nu }\to 1/2$, a particle-hole symmetry (PHS) emerges. In this case, a stripe phase manifesting PHS is the natural candidate. This is confirmed by experiments^52,53 and ED³⁷, while theoretically a 3e bubble phase is in close competition with the stripe^10,34,49. The transition between the 2e bubble phase and the stripe phase is again first-order because of their different periodic structures.

In addition to the above picture, it is found that away from $\tilde{\nu }=1/2$, the unidirectional stripe phase becomes unstable against fluctuation along the stripe direction^46,48. These fluctuations lead to a modulation where the resulting phase has every stripe broken into a 1D crystal with one electron per unit cell. The Coulomb repulsion requires neighbouring 1D crystals to have a phase difference of π while kinetic and thermal dynamics allow them to slide in the stripe direction. This competition determines whether this array of 1D crystals behaves like a unidirectional stripe or a 2D crystal experimentally. This is reflected by the shear elastic modulus^47,49 or by viewing a sliding process across one period as a soliton and studying the unbinding of soliton/anti-soliton pairs⁴⁵. Both criteria predict that when the filling $\tilde{\nu }$ deviates substantially from $\tilde{\nu }=1/2$, the unidirectional stripe phase eventually behaves as a highly anisotropic crystal, called the stripe crystal^25,44, and their transition should be continuous. However, the filling-factor range where the stripe crystal could be found coincides with the more favoured 2e bubble phase in the isotropic case. This stripe crystal is thus almost entirely covered.

Phase diagram

We compute the CDW phases for a series of mass ratios between 5 ≥ m_x/m_y ≥ 1, i.e., $1\ge \alpha \ge \sqrt{0.2}$. The phase diagrams from HF and DMRG calculations are displayed in Fig. 1a and b. A significant trend under mass anisotropy is the shrinking of the 2e-bubble regime (see Fig. 1). For m_x/m_y = 1, the 2e bubble is dominant between $\tilde{\nu }\simeq 0.22$ and $\tilde{\nu }\simeq 0.36$. By increasing the mass anisotropy, the stripe phase and the WC at small $\tilde{\nu }$ expand their respective regions in the phase diagram and finally become adjacent around m_x/m_y = 5. We also find that the stripe now picks the heavier-mass x-direction, with its periodicity along the lighter-mass y-direction, in accordance with refs. ^40,41. We discuss the origin of this behaviour in section ‘Explanation for the lattice constants’.

Fig. 1: Phase diagrams and an illustration of the lattice.

Phase diagram: (a) Hartree–Fock results as a function of the mass ratio m_x/m_y and the partial filling $\tilde{\nu }$. The stripe phase is assumed to be unidirectional. The regime of stripe crystal [covered by the 2-electron (2e) bubble for m_x = m_y] is taken from the isotropic calculation in refs. ^45,47,49. For m_x/m_y = 2, a small stripe regime sandwiched between the Wigner crystal (WC) and 2e bubble appears around $\tilde{\nu }=0.2$ and it eventually takes over the 2e bubble phase at m_x/m_y = 5. b DMRG results. Both the stripe crystal and the WC take 2D lattice order. The colour represents the modulation anisotropy ${\bar{\rho }}_{\parallel }/{\bar{\rho }}_{\perp }$, where ${\bar{\rho }}_{\parallel }$ is the maximal density wave amplitude $\bar{\rho }({\bf{q}})$ with a non-zero wave vector q in the stripe direction and ${\bar{\rho }}_{\perp }$ is the counterpart in the direction perpendicular to the stripe (see Methods). Lighter colour indicates higher anisotropy, and the highest anisotropy is reached by the unidirectional stripe. c The lattices of the Wigner crystal and the stripe crystal. Left: triangular Wigner crystal in isotropic systems. Right: stretched (S) lattice deformation. When the lattice deformation is sufficient, this becomes a stripe crystal and the stripe phase forms after the modulation along y-direction becomes small, indicated by the orange area.

We notice that the HF and the DMRG calculations are complementary in describing the unidirectional stripe phase and the stripe crystal. The HF method here does not take into account the stripe-direction modulation inside the unidirectional stripe phase. So we indicate the regime of the stripe crystal computed from refs. ^45,47,49, in which it is believed that the unidirectional stripe phase computed here should actually correspond to the stripe crystal²⁵. Meanwhile, the DMRG calculation naturally includes the stripe-direction fluctuations as it captures the exact ground state. In Fig. 1b, we use the Fourier decomposition of the density for non-zero wave vectors along the stripe direction to demonstrate the stripe modulation, as shown in colour plot. It is, however, difficult to distinguish weak modulation from zero modulation, see Supplementary Note 3. The stripe crystal computed here near $\tilde{\nu }=1/2$ can behave as a unidirectional stripe in experiments. This is the case near the isotropic limit, where the modulation has been predicted to be very weak^46,49 and likely to be beyond experimental probes.

In the presence of anisotropy, a feature worth noticing is that the lattice constants have two local minima in energies. One is a high aspect-ratio rhombus, denoted as the stretched (S) lattice, as shown in Fig. 1c. The other is closer to a square lattice, denoted as the compressed (C) lattice, which can be found in Methods. The physics behind this can be illustrated from the deformation of the isotropic triangular lattice. When we add anisotropy with D₂ symmetry, the diagonals of the rhombus are reoriented along the two principal axes of the D₂ anisotropy. There are two choices of reorientation, with the long diagonal along either the x-direction or the y-direction. For the anisotropy considered in our case, the x-direction length should be compressed and y-direction length should be stretched. As a result, we can expect two local optimal configurations due to two different ways of reorientation. These two triangular lattices are degenerate when the interaction is isotropic. As anisotropy is switched on, one rhombus is becoming thinner, while the other becomes closer to a square.

In our high-accuracy DMRG calculations on the cylinder (Supplementary Note 2), we find that the S lattice is slightly favoured, and its dominance becomes more evident with the increase of anisotropy. For m_x/m_y = 5, the two lattices are close in energy for $\tilde{\nu }\;\lesssim \;0.2$. For higher fillings, the S lattice is much more favoured. Its periodicity happens to be closer to that of a stripe phase, as reflected in Fig. 2. As for the C lattice, it tends to form a metastable stripe with a period half of the lowest-energy stripe. The details based on HF computation can be found in Methods. In the following, we assume that the S lattice represents the stable phase and is the relevant one in all phase transitions.

Fig. 2: The energy and the lattice constants of different phases.

a The energy per particle of different phases as a function of filling $\tilde{\nu }$ at mass ratio m_x/m_y = 1 and m_x/m_y = 5, computed in the Hartree–Fock (HF) approximation in the unit of e²/l. b The lattice constants of the Wigner crystal and the stripe from HF calculation. We present the half diagonal length l_b/2 of the rhombus unit cell and the stripe period. The dashed lines are the length scales set by the HF minimum 2π/q^* at m_x/m_y = 1 and $2\pi /{q}_{1}^{* }$ at m_x/m_y = 5. The solid blue line represent the value ${l}_{{\rm{b}}}=\sqrt{4\pi \sqrt{3}/\tilde{\nu }}l$, corresponding to an exact triangular lattice. As the stripe phase takes a unidirectional form in the HF ansatz, there is a discontinuity for m_x/m_y = 5 because of different symmetry orders. But further DMRG calculation reveals that quantum fluctuations lead to stripe crystal at this filling and smooth the discontinuity. c The energy per particle of the ground state at each $\tilde{\nu }$ by DMRG calculations. d The half diagonal length l_b/2 of the Wigner/stripe crystal from DMRG calculations. The colour represents the stripe-direction modulation as in Fig. 1b. The absence of some data points of (d) comparing to (b) is due to that stripe phase cannot be realised as a stable state in this region.

Energy per particle and lattice constants

The results of the energy per particle E_pp from HF calculations are summarised in Fig. 2a. The S lattice is employed here (for the C lattice, see Supplementary Note 1). We can see that upon increasing the mass ratio, the 3e bubble is no longer in close competition with the stripe phase, as the symmetry of the latter fits better the external anisotropy which breaks rotation symmetry down to D₂. To see how the WC evolves into the stripe, we compare the lattice constant l_b [parameterisation in Fig. 1c], the characteristic scale 4π/q^*, which will be discussed in the next section, and the stripe period in Fig. 2b. An important feature is that these quantities approach almost the same value at the transition point for larger anisotropy. We will further elaborate this in the subsequent two sections.

The corresponding data for energy and lattice constant from DMRG is showed in Fig. 2c and d. In the isotropic limit, the phase boundaries revealed by Fig. 2c are consistent with the HF result (Fig. 2a) up to a small difference. In Fig. 2c, the interpolated curve of E_pp shows clearly discontinuity in its derivative around $\tilde{\nu }=0.22$ and $\tilde{\nu }=0.36$, corresponding to the first-order transition between the WC and the 2e bubble phase and that between the 2e bubble and the stripe phase. While not plotted in Fig. 2c, we also confirm the competing 3e bubble predicted in ref. ¹⁰ with a energy density difference as small as ~ +10⁻⁴(e²/ϵl).

For m_x/m_y = 5, the stripe crystal and the WC become adjacent. As the stripe crystal arises from the modulation of unidirectional stripes^46,54,55, it has one period (l_b) being locked around the characteristic scale 4π/q^* (see Methods). In Fig. 2d, for filling $\tilde{\nu }\,> \,0.2$, there is a variation in l_b smaller than 10%. We consider this region as the stripe region in our DMRG calculation. For relatively low filling, $0.2\,<\,\tilde{\nu }\,<\,0.3$, our calculation clearly shows that there is a density modulation along the stripe. For $\tilde{\nu } \sim 1/2$, the modulation becomes rather small. Similar to the crystal instability of unidirectional stripes in the isotropic limit^46,49, we see that a modulation develops smoothly within the stripe region, i.e., with no clear signature of discontinuity in the derivatives of E_pp (Fig. 2c) and l_b (Fig. 2d). Furthermore, under the mass anisotropy m_x/m_y = 5, the period data indicates no clear discontinuity of the lattice shape between the WC region and the stripe region. However, the derivative or higher order derivative of period and energy with respect to filling is less smooth near $\tilde{\nu }\approx 0.2$. This indicates a transition or a fast crossover between the WC and the stripe region, as we will explain in more detail below.

Explanation for the lattice constants

Let us first clarify several notions from the HF approximation. Under this approximation, the energy of the system can be expressed as the product of the electron mean-field densities with the HF interaction potential u_HF, which is given by the Coulomb potential minus its Fourier transform (see Methods). The HF potential u_HF(q) has minima as a function q (see Fig. 3). In the isotropic case, the minimum is the same for all directions, denoted as q^*. When anisotropy comes in, we denote the minimum along the y-direction as ${q}_{1}^{* }$ while the minimum along the x-direction as ${q}_{2}^{* }$. We can see clearly that the minimum ${q}_{1}^{* }$ is much lower in energy than ${q}_{2}^{* }$. The former serves as the global minimum and explains why the stripe for m_x/m_y = 5 lies along the x-direction.

Fig. 3: The evolution of the repulsion energy $\bar{E}$ with respect to the partial filling $\tilde{\nu }$.

The curves represent the Hartree–Fock (HF) interaction potential u_HF(q) with E in units of e²/l and q in units of 1/l, where l is the magnetic length. The snowflake represents the average $\bar{E}$ between the origin and its six nearest neighbours in the reciprocal lattice. a Isotropic case. The arrows show how $\bar{E}$ evolves with increasing $\tilde{\nu }$. At small $\tilde{\nu }$, the lattice adopts a triangular shape with the average repulsion energy ${u}_{{\rm{HF}}}(\tilde{{{\Lambda }}})$. When $\tilde{\nu }$ increases, two of the six nearest neighbours have a tendency to stay at a distance q^*, while the other four are repelled away from the origin, staying at a distance q_l. The green line roughly indicates the average repulsion from q_l and q^*. As $\tilde{{{\Lambda }}}$ is approaching q_m, a stretched lattice becomes more and more favoured and eventually causes a transition in the lattice shape from the triangular lattice. b Under large mass ratio m_x/m_y = 10. The dashed/solid curve corresponds to the interaction along the y/x-direction. When $\tilde{\nu }$ increases, two of the six nearest neighbours soon fall to the distance ${q}_{1}^{* }$. The crystal then evolves with this distance fixed while the other four nearest neighbours move away. The state becomes a stripe crystal phase through a continuous phase transition, or even crossover without transition.

Now we analyse the lattice constants for the WC and the stripe crystal. In the isotropic case, previous studies⁴⁹ found that a triangular WC should exhibit a sharp deformation around $\tilde{\nu }\simeq 0.22$ when $\tilde{\nu }$ is increased from $\tilde{\nu }=0$. This is also reflected in our simple HF calculation and DMRG calculation. In Fig. 2b, we find that the stripe period and the Wigner crystal have a discrepancy at their transition point for m_x = m_y. In Fig. 2d, a sharp deformation is also found near $\tilde{\nu }\simeq 0.22$. The point marked with an arrow represents a metastable stripe state which has slightly higher energy than the stable 2e bubble phase. This point is also marked by an arrow in Fig. 1b. The highly anisotropic crystal after this deformation is interpreted as the stripe crystal, whose symmetry is different from a triangular lattice and a phase transition should take place.

For sufficient anisotropy, we find that in our HF calculations the discrepancy in periodicity between the stripe and the Wigner crystal disappears. Our DMRG calculation shows that there is no clear discontinuity in the first derivative of ${E}_{{\rm{pp}}}(\tilde{\nu })$ and no sharp deformation of the lattice structure, but instead a minimum of l_b along the light-mass axis near $\tilde{\nu }\approx 0.2$.

We first provide a simple calculation based on the HF approximation to illustrate why for m_x/m_y = 1 the deformation from a triangular lattice to a stripe crystal is so sharp and how it becomes smooth when mass anisotropy is large enough. The energy E_pp is given as a summation over reciprocal lattice vectors [Eq. (21)]. As a rough estimate of E_pp, we consider only the first shell, the six nearest neighbours of the origin that give the largest contribution. For a triangular lattice under isotropic masses, all nearest neighbours have equal distance $\tilde{{{\Lambda }}}=4\pi /(\sqrt{3}{{\Lambda }})$ to the origin, with ${{\Lambda }}=\sqrt{4\pi {l}^{2}/(\sqrt{3}\tilde{\nu })}$ the triangular lattice constant in coordinate space. The average repulsion energy is given by $\bar{E}={u}_{{\rm{HF}}}(\tilde{{{\Lambda }}})$. As $\tilde{\nu }$ starts to increase from 0, the average repulsion starts to decrease as shown in Fig. 3a. At $\tilde{\nu }={\tilde{\nu }}^{{\rm{c}}}\simeq 0.15$, the distance reaches the HF minimum $\tilde{{{\Lambda }}}={q}^{* }$ and $\bar{E}$ is globally minimal. As $\tilde{\nu }$ increases further, the lattice remains triangular for a finite range of $\tilde{\nu }$. To illustrate this, we compute its energy compared to a deformed one, where we have two nearest neighbours stay at the distance q^* but keep the other four at a larger distance q_l to maintain the area of the unit cell, resulting in a rhombic lattice (see Fig. 3a). The energy of the two configurations is computed up to a quadratic expansion of u_HF around q^*:

$${u}_{{\rm{HF}}}(q)={u}_{{\rm{HF}}}({q}^{* })+c{(q-{q}^{* })}^{2}+\ldots$$

(5)

where c is in general positive as for a minimum. Then the average repulsion for the triangular lattice is given by

$${\bar{E}}_{{\rm{tri}}}={u}_{{\rm{HF}}}(q)\simeq {u}_{{\rm{HF}}}({q}^{* })+c\ \delta {\tilde{{{\Lambda }}}}^{2},$$

(6)

where the deviation $\delta \tilde{{{\Lambda }}}$ is

$$\tilde{{{\Lambda }}}=\sqrt{\frac{4\pi \tilde{\nu }}{\sqrt{3}{l}^{2}}}\Rightarrow \,\delta \tilde{{{\Lambda }}}=\frac{{q}^{* }}{2{\tilde{\nu }}^{{\rm{c}}}}\delta \nu .$$

(7)

For the deformed lattice, the average repulsion is

$${\bar{E}}_{{\rm{dfm}}}=\frac{{u}_{{\rm{HF}}}({q}^{* })+2{u}_{{\rm{HF}}}({q}_{{\rm{l}}})}{3},$$

(8)

where q_l is related to the filling factor by

$${q}_{{\rm{l}}}=\sqrt{{\left(\frac{{q}^{* }}{2}\right)}^{2}+{\left(\frac{2\pi \tilde{\nu }}{{q}^{* }}\right)}^{2}}\Rightarrow \delta {q}_{{\rm{l}}}=\frac{3{q}^{* }}{4{\tilde{\nu }}^{{\rm{c}}}}\delta \nu .$$

(9)

Inserting these expressions, we find that near $\tilde{\nu }={\tilde{\nu }}^{{\rm{c}}}$

$${\bar{E}}_{{\rm{tri}}}-{\bar{E}}_{{\rm{dfm}}}\simeq -\frac{c{q}^{* 2}}{8{\tilde{\nu }}^{c2}}\delta {\nu }^{2}\,<\,0.$$

(10)

This illustrates that a deformed lattice is energetically unfavourable near ${\tilde{\nu }}_{c}$, because the larger values of q_l lead to a cost in $\bar{E}$, and the crystal is still triangular lattice for $\tilde{\nu }\,\le\, {\tilde{\nu }}^{{\rm{c}}}$.

However, if $\tilde{{{\Lambda }}}$ is approaching q_m, the first local maximum of u_HF, this cost for deformation no longer exists. For q > q_m, the energy curve becomes very flat and the repulsion gains very little when the distance in the reciprocal lattice is further enlarged. In this situation, one can imagine that if four of the six nearest neighbours are pushed further to a larger distance q_l while the other two keep a distance of q^*, the energy is much lowered than for the triangular lattice:

$${u}_{{\rm{HF}}}\left(\tilde{{{\Lambda }}}\right)> \frac{2{u}_{{\rm{HF}}}({q}_{{\rm{l}}})+{u}_{{\rm{HF}}}({q}^{* })}{3},\,{\rm{for}}\ \tilde{{{\Lambda }}} \sim {q}_{{\rm{m}}}.$$

(11)

Such a deformed lattice, if one considers its density profile in coordinate space, is rather similar to an array of 1D crystals equally spaced by 2π/q^*. This is exactly a stripe crystal density profile, which fits well with the deformed lattice found in ref. ⁴⁹ by taking l_a = 2π/q^* and ${l}_{{\rm{b}}}=2\pi {l}^{2}/(\tilde{\nu }{l}_{{\rm{a}}})$. According to this simple calculation, the lattice constants should have a sudden jump between $\tilde{\nu }\simeq 0.15$ and $\tilde{\nu }\simeq 0.44$. The very sharp deformation found numerically⁴⁹ at $\tilde{\nu }\simeq 0.22$ is indeed consistent with this simple picture.

Let us turn to the anisotropic case, which we illustrate for an artificially high mass ratio m_x/m_y = 10, albeit our computation shows that around m_x/m_y = 5, the lattice shape is already continous. We consider the S lattice, whose shape eventually evolves into that of a stripe phase. The HF interaction in the x- and y-directions takes rather different shapes (see Fig. 3b). The lattice is anisotropic from low fillings. We use q_s and q_l to parameterise it as in Fig. 3b. We can again perform the quadratic expansion when ${q}_{1}^{* }$ is reached. In this situation, the constant c is direction dependent, for example c_x for the x-direction and c_y for the y-direction. If c_y ≫ c_x, one can expect that the lattice period should be nearly fixed in the y-direction while letting the x-direction period grow with $\tilde{\nu }$, becoming the shape of a stripe crystal immediately after passing ${q}_{1}^{* }$. Furthermore, a simple calculation shows that when q₁ is reached, q_l is still in the descending regime of u_HF. This even more supports that q_l keeps increasing after the system reaches ${q}_{1}^{* }$ in the y-direction. As a result, the above isotropic sharp deformation does not happen. This is indeed reflected through both our calculations in Fig. 2b and d. We find that as m_x/m_y increases, the nearest-neighbour distance q_s in y-direction increasingly drops to and remains fixed at the minimum q₁. Only the distance q_l now evolves with $\tilde{\nu }$ and the system behaves as the stripe crystal without further large deformation. The scale in the y-direction rapidly goes from the dilute limit, where the lattice constant depends primarily on $\tilde{\nu }$, to the HF minimum dominated limit, where this lattice constant is fixed by q^*. In this situation, the lattice constant is a continuous function of $\tilde{\nu }$.

Analysis for the transition between a dilute-limit WC and a stripe crystal

Let us now investigate how the energy per particle transits between different CDW phases. We focus on the situation where electrons form a 2D lattice structure. The general form of E_pp comes from the inspiration of HF results. From Eqs. (20) and (21), the energy per particle is proportional to the sum of u_HF(q)∣ρ_b(q)∣² over the reciprocal lattice. For a given lattice l_a, l_b, the density profile ρ_b on each lattice site is worked out by minimising E_pp. As ρ_b describes how local electrons at one lattice site interact with neighbouring lattice sites, it should rely on the lattice structure and the electron density, ${\rho }_{{\rm{b}}}={\rho }_{{\rm{b}}}({l}_{{\rm{a}}},{l}_{{\rm{b}}},\tilde{\nu })$. It may be reasonable to regard it as a smooth function on its variables as long as the electrons are centred around each lattice site and the electron number per site is fixed to M. In this case we can further write ρ_b = ρ_b(l_a, l_b, M). Then E_pp can be explicitly parameterised by

$${E}_{{\rm{pp}}}\left({l}_{{\rm{a}}},{l}_{{\rm{b}}},\tilde{\nu }\right)={E}_{{\rm{pp}}}\left[2\pi /(\tilde{\nu }{l}_{{\rm{b}}}),{l}_{{\rm{b}}},\tilde{\nu }\right].$$

(12)

The dependence on l_a, l_b and $\tilde{\nu }$ also enters through the reciprocal lattice summation and an overall factor ${n}_{B}\tilde{\nu }/2$ respectively besides ρ_b. From such a structure, ${E}_{{\rm{pp}}}\left({l}_{{\rm{a}}},{l}_{{\rm{b}}},\tilde{\nu }\right)$ is continuous on its three variables. As ${l}_{{\rm{a}}}=2\pi /(\tilde{\nu }{l}_{{\rm{b}}})$, only l_b is a variational parameter when one is looking for the lowest-energy state at a certain filling. The ground state should satisfy:

$$\frac{\delta }{\delta {l}_{{\rm{b}}}}{E}_{{\rm{pp}}}\left({l}_{{\rm{a}}},{l}_{{\rm{b}}},\tilde{\nu }\right)=0,\,\Rightarrow \frac{\partial {E}_{{\rm{pp}}}}{\partial {l}_{{\rm{b}}}}-\frac{2\pi }{\tilde{\nu }{l}_{{\rm{b}}}^{2}}\frac{\partial {E}_{{\rm{pp}}}}{\partial {l}_{{\rm{a}}}}=0.$$

(13)

The solution to the above equation gives l_b as a function of $\tilde{\nu }$, ${l}_{{\rm{b}}}(\tilde{\nu })$. Thus the first-order derivative $d{E}_{{\rm{pp}}}/d\tilde{\nu }$ is:

$$\begin{array}{lll}\frac{d{E}_{{\rm{pp}}}}{d\tilde{\nu }} =\left(\frac{\partial {E}_{{\rm{pp}}}}{\partial {l}_{{\rm{b}}}}-\frac{2\pi }{\tilde{\nu }{l}_{{\rm{b}}}^{2}}\frac{\partial {E}_{{\rm{pp}}}}{\partial {l}_{{\rm{a}}}}\right)\frac{d{l}_{{\rm{a}}}}{d\tilde{\nu }}+\frac{\partial {E}_{{\rm{pp}}}}{\partial \tilde{\nu }}-\frac{2\pi }{{\tilde{\nu }}^{2}{l}_{{\rm{b}}}}\frac{\partial {E}_{{\rm{pp}}}}{\partial {l}_{{\rm{a}}}}\\ =\frac{\partial {E}_{{\rm{pp}}}\left({l}_{{\rm{a}}},{l}_{{\rm{b}}},\tilde{\nu }\right)}{\partial \tilde{\nu }}-\frac{2\pi }{{\tilde{\nu }}^{2}{l}_{{\rm{b}}}}\frac{\partial {E}_{{\rm{pp}}}}{\partial {l}_{{\rm{a}}}}.\end{array}$$

(14)

The second line is obtained by inserting Eq. (13).

First, we verify the above expression by analysing the transition between different bubbles phases. At the transition point, the number of electrons at each lattice site changes abruptly. Therefore ρ_b is discontinuous on the two sides of the transition. Meanwhile Eq. (14) is still valid for either side. Since ρ_b and l_a, l_b are different for different bubble phases, $\partial {E}_{{\rm{pp}}}/\partial \tilde{\nu }$ is discontinuous at the transition point, signifying a first-order transition.

Then we turn to the WC and the stripe crystal. At low fillings, ${l}_{{\rm{b}}}(\tilde{\nu })$ is controlled by the long distance Coulomb tail, and the system forms the dilute-limit WC. As the electrons become denser, our picture in the last section tells us that at intermediate fillings the lattice structure is dominated by the HF minimum q^* as a stripe crystal. We denote the two kinds of dependence as ${l}_{{\rm{b}}}^{{\rm{w}}}(\tilde{\nu })$ and ${l}_{{\rm{b}}}^{{\rm{s}}}(\tilde{\nu })$. The former is controlled by the electron density and the ~1/r repulsion. It slightly deviates from the triangular isotropic result ${l}_{{\rm{b}}}=\sqrt{4\pi \sqrt{3}/\tilde{\nu }}l$, while the later is fixed around 4π/q^*. Analysing Eq. (14) in the isotropic situation, around the transition point, ${l}_{{\rm{b}}}^{{\rm{w}}}(\tilde{\nu })$ deforms sharply to ${l}_{{\rm{b}}}^{{\rm{s}}}(\tilde{\nu })$. As the anisotropy increases, such a sharp deformation should disappear. The lattice structure is continuous at the transition point ${\tilde{\nu }}^{* }$ between the WC and the stripe crystal, ${l}_{{\rm{b}}}^{{\rm{w}}}({\tilde{\nu }}^{* })={l}_{{\rm{b}}}^{{\rm{s}}}({\tilde{\nu }}^{* })$. The first-order derivative $d{E}_{{\rm{pp}}}/d\tilde{\nu }$ is continuous, but in the second-order derivative, the expression $d{l}_{{\rm{b}}}/d\tilde{\nu }$ appears:

$$\begin{array}{lll}\frac{{d}^{2}{E}_{{\rm{pp}}}}{d{\tilde{\nu }}^{2}}= \frac{d{l}_{{\rm{b}}}}{d\tilde{\nu }}\bigg(\frac{\partial }{\partial {l}_{{\rm{a}}}}-\frac{2\pi }{\tilde{\nu }{l}_{{\rm{b}}}^{2}}\frac{\partial }{\partial {l}_{{\rm{a}}}}\bigg)\bigg(\frac{\partial {E}_{{\rm{pp}}}}{\partial \tilde{\nu }}-\frac{2\pi }{{\tilde{\nu }}^{2}{l}_{{\rm{b}}}}\frac{\partial {E}_{{\rm{pp}}}}{\partial {l}_{{\rm{a}}}}\bigg)\\ +\,\frac{4\pi }{{\tilde{\nu }}^{3}{l}_{{\rm{b}}}}\frac{\partial {E}_{{\rm{pp}}}}{\partial {l}_{{\rm{a}}}}-\frac{2\pi }{{\tilde{\nu }}^{2}{l}_{{\rm{b}}}}\frac{{\partial }^{2}{E}_{{\rm{pp}}}}{\partial \tilde{\nu }\partial {l}_{{\rm{a}}}}+\frac{{\partial }^{2}{E}_{{\rm{pp}}}}{\partial {\tilde{\nu }}^{2}}.\end{array}$$

(15)

Because ${l}_{{\rm{b}}}^{{\rm{w}}}({\tilde{\nu }}^{* })$ and ${l}_{{\rm{b}}}^{{\rm{s}}}({\tilde{\nu }}^{* })$ are controlled by different ranges of the interaction, their first-order derivatives could be different (see Fig. 2d). In that case the transition is second order. This does not rule out the possibility that the transition can be higher orders, or no phase transition in the strict sense separating the WC and stripe region, in contrast to the isotropic limit. There, a phase transition must take place when the WC is adjacent to the stripe crystal because of the symmetry difference of the two crystalline CDWs.

Experimental indications

As for experimental implications of our results, notice first that the mass anisotropy of AlAs quantum wells⁴¹ are suggested to be m_x/m_y ≈ 5.

As we have shown, the mass anisotropy leads to the disappearance of the 2e bubble phase. The 3e bubble phase ceases to be in close competition with the stripe phase. As a result, the region of the stripe phase is greatly enhanced and stabilised. In transport measurements, the stripe direction yields the easy direction while the periodic direction is the hard direction. If one measures the longitudinal resistivities along two directions, ρ_xx and ρ_yy, the result would manifest incommensurate behaviours in the stripe phase. In the isotropic case, such a large anisotropy in the longitudinal resistivity is observed near $\tilde{\nu } \sim 1/2$^52,53. As mass anisotropy is increased, we expect that this behaviour will extend down to lower fillings.

We also find that first-order phase transitions between different CDW phases are replaced by continuous phase transitions (or even no transitions) with the increase of mass ratio. As for the first-order phase transition between the WC and the 2e bubble phase, it can be detected through transport measurement under microwave irradiation^36,56,57,58. The pinning mode due to disorder manifests itself as the resonance peak of the real part of the longitudinal conductivity ${\rm{Re}}{\sigma }_{xx}(\omega )$ at finite frequency. It detects the periodic structure of CDW phases. For the isotropic situation ${\rm{Re}}{\sigma }_{xx}(\omega )$ was found to exhibit two peaks corresponding to coexistence between a WC and a 2e bubble phase around the first-order transition point³⁶. As the filling is increased, the weight of the WC is lowered and finally disappears. When mass anisotropy comes into play, such a first-order phase competition is replaced by a continuous phase transition between the WC and the stripe crystal as we showed in the last section. The periodic structure deforms smoothly through the two phases. Therefore we expect that only one peak appears throughout the transition region for intermediate $\tilde{\nu }$, corresponding to that of the WC or stripe crystal. The position of the microwave resonance smoothly changes throughout this phase transition point. The pinning mode is also feasible for the transition between the crystal phase and the unidirectional stripe⁵⁸. For the unidirectional stripe, there is a resonance peak for the longitudinal conductance along the easy direction while no resonance occurs along the hard direction. As now the 2e and 3e phases are completely removed, the system has fewer candidates and it may be interesting to see how the pinning modes in the stripe crystal should eventually evolve into that of a unidirectional stripe. This can help us understand better the nature of the QH stripe phase.

Ref. ⁵⁹ suggests that the magnetic susceptibility may be used as a tool to detect CDW phase transitions. This quantity is related to the second derivative of ${E}_{{\rm{pp}}}(\tilde{\nu })$. We have shown that at the transition between the dilute-limit WC and the stripe crystal, the energy ${E}_{{\rm{pp}}}(\tilde{\nu })$ can be at most discontinuous in its second order derivative. Therefore such an experiment on susceptibility should be able to uncover the WC-stripe crystal transition.

Discussion

Our analytic HF and numerical DMRG calculations reach a remarkable agreement in studying the CDW phases under mass anisotropy in the N = 2 LL. The mass anisotropy suppresses the 2e bubble phase and enlarges the region of the stripe phase in the phase diagram. In particular, the previously predicted stripe crystal now dominates at intermediate fillings. The first-order phase transitions between the WC and the 2e bubble phase are replaced by that between the WC and the stripe crystal. While they are separated by a sharp deformation in the isotropic limit, in the anisotropic case, no sharp deformation between them is observed, but a second-order phase transition is likely to take place. We nevertheless do not rule out the possibility of a crossover due to the discreteness of our numerical data. Our results can lead to many interesting experimental phenomena that enhance the understanding of strong correlation in QH systems.

There are a few perspectives from the work in this paper. Notice that in the isotropic case, there can be a Kosterlitz-Thouless transition from the stripe crystal to the unidirectional stripe phase⁴⁵ due to the proliferation of soliton-anti-soliton pairs. How such a transition point evolves under anisotropy remains interesting. This however requires including stripe-direction quantum fluctuations that are beyond the simple HF method used in this work, where we limit our calculations to the classical density profiles. One needs to resort to quantum density profiles. For example in ref. ⁴⁶, the density profile is assumed to take ρ_b as that of a M-particle $\tilde{\nu }=1$ wave function⁶⁰. But there is no such a simple quantum density profile for the stripe phase. A systematic approach is the more complicated self-consistent HF approximation⁵⁴. On the other hand, our DMRG calculation with the implementation of symmetry breaking and quasi-momentum conservation allows us accurately to compute the lattice shapes and energy. The data also reveal the instabilities of a unidirectional stripe to form stripe crystals at intermediate filling. But closer to the half filling, the density modulation in the stripe direction becomes very weak and it is difficult to conclusively distinguish the unidirectional stripe from the stripe crystal. The precise determination of weak stripe crystals or unidirectional stripes calls for more sophisticated data extrapolation.

Methods

LL projection for anisotropic masses

The total Hamiltonian is given by the kinetic part H_k plus the interaction part (1/2A)∑_qV(q)ρ(q)ρ( − q), where V(q) is the Coulomb interaction and $\rho ({\bf{q}})=\mathop{\sum }\nolimits_{i}^{{N}_{{\rm{e}}}}\exp (-i{\bf{q}}\cdot {{\bf{r}}}_{i})$ with r_i being the coordinate of the electron. The kinetic energy can be expressed in terms of the following one-particle Hamiltonian with anisotropic masses:

$${H}_{k}=\frac{{{{\Pi }}}_{x}^{2}}{2m/\alpha }+\frac{{{{\Pi }}}_{y}^{2}}{2m\alpha }.$$

(16)

where the mass ratio is characterised by m_y/m_x = α². The covariant momentum is given as Π = p + ∣e∣A in a magnetic field ${\bf{B}}=\nabla \times {\bf{A}}=-B\hat{{\bf{z}}}$. The LL projection can be achieved by decomposing the electron coordinates into the cyclotron and guiding centre coordinates as η_μ = −∑_νl²ϵ_μνΠ_ν, R_μ = r_μ − η_μ, where $l=\sqrt{\hslash /eB}$ and ϵ_μν is the anti-symmetric tensor. The cyclotron coordinates completely determine the LL structure. In the presence of mass anisotropy the ladder operators for LLs are related to η_x and η_y through

$$a=\frac{{\eta }_{x}}{\sqrt{2\alpha }l}+i\frac{{\eta }_{y}\sqrt{\alpha }}{\sqrt{2}l}\quad {a}^{\dagger }=\frac{{\eta }_{x}}{\sqrt{2\alpha }l}-i\frac{{\eta }_{y}\sqrt{\alpha }}{\sqrt{2}l}.$$

(17)

The LLs are defined by the ladder operators through $\left|N\right\rangle ={a}^{\dagger N}/\sqrt{N!}\left|0\right\rangle$.

The projection to the LL N is done by averaging the density operator $\langle N| \rho ({\bf{q}})| N\rangle ={\sum }_{i}\langle N| \exp [-i{\bf{q}}\cdot ({{\boldsymbol{\eta }}}_{{\boldsymbol{i}}}+{{\bf{R}}}_{i})]| N\rangle$ in the interaction Hamiltonian. After this procedure, the kinetic energy can be dropped.

Hartree–Fock approximation

The analytic HF method provides a physical picture in understanding the reaction of CDWs to anisotropy at a mean-field level, whereas the DMRG numerical calculations incorporate quantum fluctuations absent in the mean-field theory. Different CDW orders, typically the unidirectional stripe order and the 2D crystalline order, are essentially captured by the analytic HF method. On the other hand, the modulation induced by quantum fluctuations in the stripe direction, which leads to the formation of the stripe crystal, is beyond our simple HF approximation (in contrast to the self-consistent HF approximation⁵⁴). Such a stripe crystal can be manifested in the DMRG result. It probes into the regime when the stripe crystal is adjacent to the low-filling WC.

With the CDW orders, the HF approximation enables us to extract the energy of different states. The WC and the M-electron bubble phase possess the 2D lattice order, while the stripe phase takes the unidirectional order. Within the HF approximation, the energy of the system is given by:

$${E}_{{\rm{HF}}}=\frac{1}{2A}\mathop{\sum}\limits_{{\bf{q}}}{u}_{{\rm{HF}}}({\bf{q}})\langle \bar{\rho }({\bf{q}})\rangle \langle \bar{\rho }(-{\bf{q}})\rangle .$$

(18)

In our calculation, we will jump between the discrete sum ∑_q and the continuous integral ∫dq whenever convenient, related through ∑_q = A∫dq/(2π)². The HF potential u_HF is given by the effective potential minus its Fourier transformation⁹

$${u}_{{\rm{HF}}}({\bf{q}})={V}_{{\rm{eff}}}({\bf{q}})-\mathop{\sum}\limits_{{\bf{p}}}\frac{{V}_{{\rm{eff}}}({\bf{p}})}{{N}_{\phi }}\exp (-i{\bf{p}}\times {\bf{q}}{l}^{2}),$$

(19)

The first and second terms are usually called the direct and the exchange interaction, respectively. As common in the HF approximation, the direct interaction is repulsive while the exchange interaction is attractive. Looking back at Eq. (2), because of the Laguerre polynomial in the form factor F_N, the direct interaction has a series of zeros. At its first zero, the repulsive interaction vanishes and a large attractive exchange interaction dominates. This leads to a global minimum q^* in the HF potential²⁵ (also Fig. 3). For CDWs, they will greatly benefit from this minimum if the periodicity is consistent with q^*. In general the CDW period scales indeed as $1/{q}^{* } \sim l\sqrt{2N+1}$ for the LL index N²⁵. This is the reason why the HF approximation becomes increasingly accurate for larger N: the quantum fluctuations arise at the edges of bubbles and stripes, taking place at a length scale l. For large N these edge fluctuations are small compared to the CDW periods.

Now we work out the energy per particle E_pp = E_HF/N_e. The WC and bubble phases can be considered together as the former is equivalent to an 1e-bubble phase. For a system on a lattice, the expectation value of the density takes a decomposition:

$$\langle \bar{\rho }({\bf{q}})\rangle ={\rho }_{{\rm{b}}}({\bf{q}})\mathop{\sum}\limits_{j}{e}^{-i{\bf{q}}\cdot {{\bf{R}}}_{j}},$$

(20)

in which R_j are lattice vectors labelled by j. The function ρ_b represents the density profile at each lattice site. In QH systems, the non-commutativity of the x- and y-guiding centre coordinates requires that each electron must occupy a minimal surface smeared over an area ~ l². So at each lattice site the density distribution cannot be point-like (also see Fig. 4a). For an M-electron bubble phase, the normalisation requires ρ_b(0) = M. The area A_u of the unit cell and the highest partially filled LL filling $\tilde{\nu }$ satisfy $2\pi {l}^{2}M=\tilde{\nu }{A}_{u}$. Using the identity ${\sum }_{j}\exp (-i{\bf{q}}\cdot {{\bf{R}}}_{j})=(A/{A}_{u}){\sum }_{{\bf{Q}}}{\delta }_{{\bf{q}},{\bf{Q}}}$, where Q is a reciprocal lattice vector, one can find that the energy per particle is given by:

$${E}_{{\rm{pp}}}\equiv \frac{{E}_{{\rm{HF}}}}{{N}_{{\rm{e}}}}=\mathop{\sum}\limits_{{\bf{Q}}}\frac{{n}_{B}\tilde{\nu }}{2{M}^{2}}{u}_{{\rm{HF}}}({\bf{Q}}){\left|{\rho }_{{\rm{b}}}({\bf{Q}})\right|}^{2},$$

(21)

where n_B = 1/(2πl²) is the flux density. In the above summation we need to put V_eff(0) = 0 in u_HF(Q) as a result of charge neutrality. This is different from the cohesive energy used in refs. ⁸ and¹⁰, where the exchange energy at Q = 0 is also subtracted. For practical calculation, the summation in Eq. (21) is usually taken over a dozen shells in the reciprocal space for sufficient convergence of the energy.

Fig. 4: Parameterisation of the lattice and the numerical density profile.

a The unit cell of bubble phases. The bubble at each lattice site has a smearing region 2πl² = h/eB per electron, where e is the electron charge and B is the magnetic field. b Charge density waves on the infinite cylinder geometry: visualisation of the DMRG data for the density profile of 2 electron bubble crystal at 2/7 filling and the isotropic limit.

In our HF approximation, there are several assumptions for the density profile ρ_b(q) at each lattice site. Here we take a classical assumption of a uniform distribution of $\tilde{\nu }=1$ inside a disc^10,24. In coordinate space this can be represented as:

$${\rho }_{{\rm{b}}}(r)=\frac{1}{2\pi {l}^{2}}{{\Theta }}({r}_{{\rm{b}}}-r),$$

(22)

with r_b the radius of the disc (see Fig. 4a), which will be determined later. Transforming it to momentum space, we find

$${\rho }_{{\rm{b}}}(q)=\frac{{r}_{{\rm{b}}}}{{l}^{2}q}{J}_{1}(q{r}_{{\rm{b}}}),$$

(23)

where J₁ is a Bessel function. The normalisation requirement ρ_b(0) = M leads to ${r}_{{\rm{b}}}=\sqrt{2M}l$. So for the disc density profile, the energy per particle is

$${E}_{{\rm{pp}}}=\mathop{\sum}\limits_{{\bf{Q}}}\frac{{n}_{B}\tilde{\nu }}{M}{u}_{{\rm{HF}}}({\bf{Q}})\frac{{J}_{1}^{2}(\sqrt{2M}| {\bf{Q}}| l)}{| {\bf{Q}}l{| }^{2}}.$$

(24)

For the stripe phase, it is uniform in one direction but periodic in the perpendicular one. A classical density profile is to have alternating uniform $\tilde{\nu }=1$ and $\tilde{\nu }=0$ unidirectional stripes along the periodic direction¹⁰. The distribution is written as

$$\langle \bar{\rho }({\bf{r}})\rangle =\mathop{\sum}\limits_{j}\frac{1}{2\pi {l}^{2}}{{\Theta }}\left(\frac{a}{2}-| {\bf{r}}\cdot \hat{{\bf{e}}}-j{{{\Lambda }}}_{{\rm{s}}}| \right),$$

(25)

where $\hat{{\bf{e}}}$ is the direction of the CDW order, perpendicular to the stripe direction, and Λ_s is the stripe period. The width of the stripe a has to satisfy the filling factor $\tilde{\nu }=a/{{{\Lambda }}}_{{\rm{s}}}$. Its energy per particle is given by:

$${E}_{{\rm{pp}}}=\frac{{n}_{B}}{2{\pi }^{2}\tilde{\nu }}\mathop{\sum}\limits_{j}{u}_{{\rm{HF}}}\left(\frac{2\pi \hat{{\bf{e}}}}{{{{\Lambda }}}_{{\rm{s}}}}j\right)\frac{{\sin }^{2}(\pi \tilde{\nu }j)}{{j}^{2}}.$$

(26)

The stripe period is a variational parameter fixed by minimising the above E_pp. As the stripe phase only has this periodic structure, we can immediately link it to the minimum of the HF potential. It is found that this period almost coincides with the inverse of the HF minimum^8,10,24:

$${{{\Lambda }}}_{{\rm{s}}} \sim 2\pi /{q}^{* }.$$

(27)

Looking at Eq. (2), we deduce that the zeros of the direct potential are scaling as $1/\sqrt{\alpha }$ in the q_x-direction and $\sqrt{\alpha }$ in the q_y-direction. The HF minimum q^* scales in the same way. This suggests that for stripes arrayed in the x-direction, the period Λ_s behaves as $\sqrt{\alpha }$ and that for stripes arrayed in the y-direction, the period Λ_s behaves inversely. Such behaviours of the minima are also reflected in u_HF for m_x/m_y = 5 presented in Fig. 3b. Therefore for a solid phase, in coordinate space its periodicity in the y-direction is stretched while that in the x-direction is compressed.

Now we briefly mention the ansatz states in the case of an anisotropic mass. The triangular lattice is no longer the optimal one in this case. We restrict ourselves to rhombus lattices whose diagonals are along the x- and y-directions, the principal axes of the anisotropy, as parametrised in Fig. 4a. The lengths l_a and l_b satisfy

$${l}_{{\rm{a}}}{l}_{{\rm{b}}}=\frac{2\pi {l}^{2}M}{\tilde{\nu }},$$

(28)

such that a unit cell contains M electrons. Therefore l_a, l_b are not independent and we use l_a as the variational parameter. At each $\tilde{\nu }$, we work out the l_a with the lowest energy.

In addition to the deformation of the lattice, in principal we also need to consider the deformation of ρ_b(q), the smearing region of the guiding centres around each lattice site. In the isotropic case, this smearing region is assumed to be a disc. In the presence of anisotropy, this region should also be deformed like the lattice shape. Since the HF computation is to search for the lowest-energy state variationally, it is not very efficient when both deformations are included so that the searching dimensions become 2. The optimal smearing region is more accessible through numerical methods, such as the above-mentioned self-consistent HF approximation calculation^48,54 and our DMRG calculation. In order to get a good estimate of this deformation, we adapt two kinds of trial density profiles. The first is still assuming the profile ρ_b(q) to be a round shape. The second is to deform ρ_b(q) like the cyclotron orbitals, namely rescaling it by $\sqrt{\alpha }$ and $1/\sqrt{\alpha }$ in two perpendicular directions. We expect that the energy of the true deformation of ρ_b(q) can be approximated by the lower one of these two configurations. By varying $\tilde{\nu }$ and M, we find that the energy of these two configurations crosses several times. The round disc is usually more favourable when the bubbles have the tendency to cluster into a unidirectional CDW state. Notice that in this HF calculation, the stripe phase and the crystal phases always possess different symmetry orders. This is why in Fig. 2b there is a discontinuity at m_x/m_y = 5. Our DMRG calculation supplements this drawback and reveals that the actual phase near the phase boundary takes an intermediate density profile between unidirectional stripe and a well separated crystal phase. So the discontinuity in Fig. 2b is an artificial result of the simple HF ansatzes.

Details of DMRG calculations

We use the density matrix renormalization group (DMRG) to determine the ground states of the Hamiltonian on an infinitely-long cylinder geometry with circumference L^38,39,50,61. In our calculation, there are at most 4-8 periods (depends on lengths of lattice vectors with L ~ 30l) wrapped around the cylinders. We do not attempt to do extrapolation to infinite L. Nevertheless, the finite-size errors for estimating energies and lattice shapes are small, demonstrated in Fig. 5.

Fig. 5: Estimate energy density by DMRG and determine lattice shape from energy data.

a Estimate energy by extrapolating DMRG data with different bond dimension χ. The unit of energy is e²/ϵl. Our scheme is fitting the ansatz E(χ) = E − b/χ^β by fitting E(χ) as a linear function of 1/χ^β; −b is the fitted slope and an estimation of E is obtained from the intercept. We adjust the other fitting parameter β to do multiple linear fittings and determine β as that maximises the r value. The plot with β = 1.44, is an example of such fit for the ground state at cylinder circumference L = 26, fold of ‘rotational’ symmetry N_R = 4 and 1/6 filling. b Determining lattice shape and energy density by searching optimal L_op defined as the location of energy minimum for nearby L. The energies are extracted using the method illustrated in (a). We use the lattice shapes of systems with L_op to estimate that of the infinite 2D system. The fold of ‘rotational’ symmetry N_R is also the number of unit cells along the tangential direction. For a simple rhombus lattice, the length of one diagonal of the rhombus unit cell l_R can be determined as N_Rl_R = L_op. The plot shows the calculation for systems at filling 1/6 and the isotropic limit. We calculate systems with three different N_R. Two of them have the short rhombus axis along the tangential direction and one has the long rhombus axis along the tangential direction. In the thermodynamic limit, a triangular lattice Wigner crystal is expected. For exact triangular lattices, the optimal should be at where denoted by star symbols (L^*). The data shows that L_ops deviate from L^*s by small values. Accordingly, the estimated l_R deviates from ideal values ${l}_{{\rm{R}}}^{* }$. The deviation is around 1.5% for N_R = 4, 0.7% for N_R = 5 (short); around 0.02% for N_R = 3 (long). Also, we see the estimated energy densities differ by ~10⁻⁵.

DMRG optimises the matrix product state (MPS) variational wave functions^62,63 to approximate the ground states. The ‘quality’ of the MPS ansatz is controlled by the bond dimension χ, where χ is the size of the matrices used in the MPS. The optimised MPS is expected to approach the exact ground state in the χ → ∞ limit. We extrapolate the data of different χ to estimate energy densities (see e.g., Fig. 5a).

Using Landau gauge, the orbitals (single-particle eigenstates) of the kinetic Hamiltonian [Eq. (16)] are localised and aligned along the axis of the cylinder and labelled by the orbital centre/ canonical momentum ($2\pi n/L,n\in {\mathbb{Z}}$). Working with the infinite cylinders (iDMRG), we need to specify the unit cell size N_u, namely to set the MPS to be translational invariant by N_u orbitals.

Unlike working with liquid phases, the circumference L has to be chosen discretely to ensure that the crystal lattices are neither stretched nor compressed by the periodic geometry ^38,39. In other words, for those L (denoted as L_op), the ground state energy density is a minimum for nearby L.

When implementing DMRG, we need to find the pairs (N_u, L_op) as we want to obtain states with broken spatial symmetry instead of cat states. An infinite quasi-one dimensional crystal can be considered as an infinite repetition of a unit cylinder of crystal. By magnetic flux quantisation, this unit cylinder contains an integer number of orbital centres/ fluxes, denoted by N_nu. By definition, N_nu = qN_e/p, where N_e is the number of unit cell in the unit cylinder and p/q is the partial filling fraction of the projected LL ($=\tilde{\nu }$). As N_nu is an integer, N_e is an integer multiple of p. Therefore, N_nu is an integer multiple of q. The symmetry broken quantum state is thus transnational invariant by N_nu and we set N_u = N_nu.

In our situation, we find that the system always takes a rhombus lattice embedded on the cylinder with one rhombus diagonal parallel to the axis (Figs. 4 and 6). Define N_R as the number of unit cells wrapped around the cylinder circumference. We have N_e = 2MN_R and thus N_u = 2qMN_R/p, where M is the number of electrons each bubble contains. As N_u and N_R are both integers, N_u should also be an integer multiple of 2Mq; N_R is an integer multiple of p. Notice that the above condition is consistent with that to implement charge conservation in DMRG, which requires N_u to be an integer multiple of q.

Given a permissible N_u (N_R), by finding L_op, the lengths of the diagonals of the rhombus unit cell can be estimated: l_R = L_op/N_R for the diagonal along the tangential direction. In the following, we show that the estimated l_R, as well as the energy densities, shifts very slightly for the values N_R we study in DMRG.

We compare the energy densities and lattice shapes of the data of different N_R to estimate the finite size error for evaluating those quantities. As the computation cost grows exponentially with N_R, we limit our calculation to a small N_R. Furthermore, since the Hamiltonian is anisotropic, we can choose the axial direction along either the heavy-mass or the light-mass direction. Comparing the results from the two choices gives a further evaluation of the finite-size effect. Here, we show our results (Fig. 5) for the system at filling $\tilde{\nu }=1/6$ in the isotropic limit. The expected triangular lattice has two natural ways to be embedded on the cylinder surface, which serves as two starting points of deformation under the two choices of introducing anisotropy (config1, config2) respectively, see Fig. 6. By computing l_R for different N_R, we find that l_R converges to the length of the short/long diagonal of the ideal triangular lattice for the two types of embedding. Similarly, the data of energy densities shows that the estimation based on finite-size data may be accurate.

Fig. 6: Mass anisotropy setting and anisotropy induced lattice deformation on cylinder geometry.

One can either choose the axial or the tangential direction as the heavy axis, denoted as config1 and config2. In the isotropic limit, the difference between heavy and light axes vanishes; the ground state is a triangular lattice which has two natural ways to be embedded on the cylinder surface. Those serve as two starting points of anisotropy induced deformation for config1 and for config2. The lattice is expected to get compressed along the heavy axis; the two starting points at the isotropic limit lead to two branches of lattice, i.e., the stretched (S) lattice (a) and the metastable compressed (C) lattice (b).

Spatial symmetry breaking in our DMRG calculation

In this part, we discuss the issue of obtaining symmetry broken states to directly detect CDW orders. To do this, we need to figure out if a spatial-symmetry broken state can be an exact ground state on infinite cylinders, and that if CDW orders can be overestimated or induced due to the finite bond dimension of MPS approximation.

Strictly speaking, spatial symmetry can only break in the thermodynamic limit. For finite systems, the exact ground state is the unique cat state. As we work with infinite cylinder geometry, the system is infinite in one direction but finite in the other. There are two kinds of spatial symmetries on the cylinder: the translational and ‘rotational’ symmetries. Here, we argue that the translational symmetry along the cylinder axis can be broken in the strict sense; the ‘rotational’ symmetry cannot be broken for the exact ground states but can be broken for the DMRG finite-χ approximated ground states.

As pointed out in ref. ³⁹, the translational symmetry can be broken in the strict sense for a system on an infinite cylinder with a finite circumference. The observation is that albeit a continuous symmetry is forbidden to be broken in one dimension, the translational symmetry of the Hamiltonian is discrete. On the other hand, the ‘rotational’ symmetry is a continuous symmetry, which is expected to be preserved.

The DMRG calculation is known to usually overestimate the symmetry broken orders once they are not enforced to vanish. The overestimation gets corrected by increasing the bond dimension. One example is that a one-dimensional quasi-order leads to a finite-bond-dimension long-range order which decays with bond dimension^64,65,66. In the current problem of QH CDW on an infinite cylinder, we find that translational symmetry breaking orders are overestimated; the ‘rotational’ symmetry breaking order is induced because of the finite bond dimension MPS.

There are two motivations to look for symmetry breaking states instead of cat states. The first motivation is numerical efficiency. We expect that there is extra entanglement entropy for a cat state comparing to corresponding symmetry broken states. For the breaking of continuous internal symmetry, it has been shown^67,68 that the extra bipartite entanglement entropy is logarithmic in the length of partition boundary. Similar result is expected to apply to spatial symmetry breaking. For the orbital basis bipartition, it is clear that translational symmetry breaking reduces the entanglement entropy by ${\mathrm{log}}\,({N}_{\text{u}})$. Notice that N_u is proportional to the cylinder circumference L. The efficiency of DMRG benefits from a less entangled target state. Working with fully symmetry-broken states is the fastest, even if we can only utilise a quasi-momentum conservation in the DMRG without continuous ‘rotational’ symmetry. The second motivation is that obtaining symmetry broken states allows for the direct confirmation of CDW phases and also the calculation of physical quantities such as densities and certain two-point correlations. In the following, we discuss details of determining QH CDW states by DMRG.

As the translational symmetry breaking along the cylinder axis is exact, it is most straightforward to measure the order along this direction, for example, the (weak) density modulation along the stripe. However, the approximated ground state by MPS with a fixed bond dimension χ may not tell us whether there is a translational symmetry breaking. Even if the exact ground state is uniform in the axial direction, the finite-χ approximation may still have some density modulation. One needs to extrapolate data of different χ to see if the density modulation vanishes in the χ → ∞ limit.

At the same time, as there is no ‘rotational’ symmetry breaking for the exact ground state, we need to interpret the ‘rotational’ symmetry breaking of DMRG data. Each data point comes from DMRG optimisation on a fixed bond dimension MPS, with fixed unit cell size, and the matrix elements are restricted to be real. If the bond dimension approaches infinity, a well-optimised state, of course, should not break the ‘rotational’ symmetry. However, restricting the variational space within the MPS ansatz with a fixed bond dimension, there exists effective pinning which decays with increasing bond dimension. This could be similar to the role of the pinning field in an exact calculation of a finite system. For a range of small pinning strength, the symmetry-breaking order can serve as an estimation of the exact result. If the pinning strength is too small, the symmetry can restore; the threshold strength should be related to the spacing of the low-lying states of the exact spectra. In our calculation of isotropic WC, we observe that for some relatively small systems ~20l, the symmetry restores, for large enough bond dimensions with estimated energy accuracy 10⁻⁶. On the other hand, if the density modulation of the exact state is too weak, the corresponding ‘rotational’ symmetry cannot break even for moderate bond dimension.