Charge density waves and their transitions in anisotropic quantum Hall systems

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.


I. 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 [11].
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 * yuchi.he@rwth-aachen.de † kang.yang@fysik.su.se 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 [27] and observed in higher LLs of graphene [28]. 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 [29] except for the possibleν = 1/5 andν = 4/5 states [30] at intermediate temperatures, whereν is the partial filling factor in the N = 2 LL. (In conventional thin GaAs quantum wells, the filling factor is ν = N e /N φ = 4 +ν.) 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) [37] 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-polarized 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 [41]. A tunable mass anisotropy of a 2D electron gas can also be realised by strain [42] or moiré pattern [43]. 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.

A. Model and relevant phases
We consider electrons with anisotropic mass m y /m x = α 2 = 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) where A is total area of the system. The projected density operatorsρ consist only of components in the N th LL and the form factor F N keeps track of the associated LL wave functions, 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 = /eB. Notice that the mass anisotropy affects the effective interaction only through the form factor. Our HF and DMRG [50] 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ν → 0 (but away from theν = 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 [51]. 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. 10). 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 pp (ν) 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ν → 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 [37], 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ν = 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 [45]. Both criteria predict that when the fillingν deviates substantially fromν = 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.

B. Phase diagram
We compute the CDW phases for a series of mass ratios between 5 ≥ m x /m y ≥ 1, i.e., 1 ≥ α ≥ √ 0.2. The phase diagrams from HF and DMRG calculations are displayed Fig. 1(a) 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ν 0.22 andν 0.36. By increasing the mass anisotropy, the stripe phase and the WC at smallν 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 Ref. 40 and 41. We discuss the origin of this behaviour in section "Explanation for the lattice constants".
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 Ref. 45, 47, and 49, in which it is believed that the unidirectional stripe phase computed here should actually correspond to the stripe crystal [25]. Meanwhile, the DMRG calculation naturally includes the stripe-direction fluctuations as it captures the exact ground state. In Fig. 1(b), we use the Fourier decomposition of the density for nonzero 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ν = 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. 1(c). 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 2 symmetry, the diagonals of the rhombus are reoriented along the two principal axes of the D 2 anisotropy. There are two choices of reorientation, with the long diagonal along either the x-direction or the ydirection. 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ν 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.

C. Energy per particle and lattice constants
The results of the energy per particle E pp from HF calculations are summarised in Fig. 2(a). 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 2 . To see how the WC evolves into the stripe, we compare the lattice constant l b [parameterisation in Fig. 1(c)], the characteristic scale 4π/q * , which will be discussed in the next section, and the stripe period in Fig. 2(b). 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. 2(c)(d). In the isotropic limit, the phase boundaries revealed by Fig. 2(c) are consistent with the HF result ( Fig. 2(a)) up to a small difference. In Fig. 2(c), the interpolated curve of E pp shows clearly discontinuity in its derivative aroundν = 0.22 and ν = 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 the Fig. 2(c), we also confirm the competing 3e bubble predicted in Ref. 10 with a energy density difference as small as ∼ +10 −4 (e 2 / 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. 2(d), for fillingν > 0.2, there is a variation in l b smaller than 10%. We consider this region as the stripe region in our DMRG calculation. 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.
For relatively low filling, 0.2 <ν < 0.3, our calculation clearly shows that there is a density modulation along the stripe. Forν ∼ 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 . 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ν ≈ 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.

D. 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. Now we analyse the lattice constants for the WC and the stripe crystal. In the isotropic case, previous studies [49] found that a triangular WC should exhibit a sharp deformation aroundν 0.22 whenν is increased from ν = 0. This is also reflected in our simple HF calculation and DMRG calculation. In Fig. 2(b), we find that the stripe period and the Wigner crystal have a discrepancy at their transition point for m x = m y . In Fig. 2(d), a sharp deformation is also found nearν 0.22. The point marked with an arrow represents a metastable stripe √ 3/νl, corresponding to an exact triangular lattice. As the stripe phase takes a unidirectional form in the HF ansatz, there is a discontinuity for mx/my = 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ν 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. 1(b). 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. state which has slightly higher energy than the stable 2e bubble phase. This point is also marked by an arrow in Fig. 1(b). 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 pp (ν) and no sharp deformation of the lattice structure, but instead a minimum of l b along the light-mass axis nearν ≈ 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 [equation (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Λ = 4π/ (  √  3Λ) to the origin, with Λ = 4πl 2 /( √ 3ν) the triangular lattice constant in coordinate space. The average repulsion energy is given byĒ = u HF (Λ). Asν starts to increase from 0, the average repulsion starts to decrease as shown in Fig. 3(a). Atν =ν c 0.15, the distance reaches the HF minimumΛ = q * andĒ is globally minimal. Asν increases further, the lattice remains triangular for a finite range ofν. 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. 3(a)). The energy of the two configurations is computed up to a quadratic expansion of u HF around q * : where c is in general positive as for a minimum. Then the average repulsion for the triangular lattice is given by: where the deviation δΛ is For the deformed lattice, the average repulsion is where q l is related to the filling factor by Inserting these expressions, we find that nearν =ν c This illustrates that a deformed lattice is energetically unfavourable nearν c , because the larger values of q l lead to a cost inĒ, and the crystal is still triangular lattice forν ≤ν c . However, ifΛ 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: 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. 49 by taking l a = 2π/q * and l b = 2πl 2 /(νl a ). According to this simple calculation, the lattice constants should have a sudden jump betweeñ ν 0.15 andν 0.44. The very sharp deformation found numerically [49] atν 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 continuous. 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. 3(b)). The lattice is anisotropic from low fillings. We use q s and q l to parameterize it as in Fig. 3(b). 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ν, becoming the shape of a stripe crystal immediately after passing q * 1 . Furthermore, a simple calculation shows that when q 1 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 ydirection. As a result, the above isotropic sharp deformation does not happen. This is indeed reflected through both our calculations Fig. 2(b) and (d). We find that as m x /m y increases, the nearest-neighbour distance q s in ydirection increasingly drops to and remains fixed at the minimum q 1 . Only the distance q l now evolves withν 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ν, 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 ν.
E. 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 equation (20) and equation (21), the energy per particle is proportional to the sum of u HF (q)|ρ b (q)| 2 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, ρ b = ρ b (l a , l b ,ν). 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 The dependence on l a , l b andν also enters through the reciprocal lattice summation and an overall factor n Bν /2 respectively besides ρ b . From such a structure, E pp (l a , l b ,ν) is continuous on its three variables. As l a = 2π/(νl 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: The solution to the above equation gives l b as a function ofν, l b (ν). Thus the first-order derivative dE pp /dν is: The second line is obtained by inserting equation (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 equation (14) is still valid for either side. Since ρ b and l a , l b are different for different bubble phases, ∂E pp /∂ν 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 b (ν) 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 w b (ν) and l s b (ν). The former is controlled by the electron density and the ∼ 1/r repulsion. It slightly deviates from the triangular isotropic result l b = 4π √ 3/νl, while the later is fixed around 4π/q * . Analysing equation (14) in the isotropic situation, around the transition point, l w b (ν) deforms sharply to l s b (ν). As the anisotropy increases, such a sharp deformation should disappear. The lattice structure is continuous at the transition pointν * between the WC and the stripe crystal, l w b (ν * ) = l s b (ν * ). The first-order derivative dE pp /dν is continuous, but in the second-order derivative, the expression dl b /dν appears: Because l w b (ν * ) and l s b (ν * ) are controlled by different ranges of the interaction, their first-order derivatives could be different [see Fig. 2(d)]. 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.

F. Experimental indications
As for experimental implications of our results, notice first that the mass anisotropy of AlAs quantum wells [41] 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 stabilized. 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 behaviors in the stripe phase. In the isotropic case, such a large anisotropy in the longitudinal resistivity is observed nearν ∼ 1/2 [52,53]. As mass anisotropy is increased, we expect that this behavior 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 Re σ xx (ω) at finite frequency. It detects the periodic structure of CDW phases. For the isotropic situation Re σ xx (ω) was found to exhibit two peaks corresponding to coexistence between a WC and a 2e bubble phase around the first-order transition point [36]. 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ν, 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 [58]. 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. 59 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 pp (ν). We have shown that at the transition between the dilutelimit WC and the stripe crystal, the energy E pp (ν) 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.

III. 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 [45] due to the proliferation of soliton-antisoliton 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. 46, the density profile is assumed to take ρ b as that of a M -particleν = 1 wave function [60]. 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 [54]. 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.

A. LL projection for anisotropic masses
The total Hamiltonian is given by the kinetic part H k plus the interaction part (1/2A) q V (q)ρ(q)ρ(−q), where V (q) is the Coulomb interaction and ρ(q) = Ne i exp(−iq · 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: where the mass ratio is characterized by m y /m x = α 2 . The covariant momentum is given as Π = p + |e|A in a magnetic field B = ∇ × A = −Bẑ. The LL projection can be achieved by decomposing the electron coordinates into the cyclotron and guiding centre coordinates as η µ = − ν l 2 µν Π ν , R µ = r µ − η µ , where l = /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 The LLs are defined by the ladder operators through |N = a †N / √ N !|0 . The projection to the LL N is done by averaging the density operator N |ρ(q)|N = i N | exp[−iq · (η i + R i )]|N in the interaction Hamiltonian. After this procedure, the kinetic energy can be dropped.

B. 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 [54]). 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: In our calculation, we will jump between the discrete sum q and the continuous integral dq whenever convenient, related through q = A dq/(2π) 2 . The HF potential u HF is given by the effective potential minus its Fourier transformation [9] u HF (q) (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 equation (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 [25] (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 * ∼ l √ 2N + 1 for the LL index N [25]. 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: in which R j are lattice vectors labeled by j. The function ρ b represents the density profile at each lattice site. In QH systems, the non-commutativity of the x-and yguiding centre coordinates requires that each electron must occupy a minimal surface smeared over an area ∼ l 2 . So at each lattice site the density distribution cannot be point-like [also see Fig. 4(a)]. 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ν satisfy 2πl 2 M =νA u . Using the identity j exp(−iq · R j ) = (A/A u ) Q δ q,Q , where Q is a reciprocal lattice vector, one can find that the energy per particle is given by: where n B = 1/(2πl 2 ) 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. 8 and 10, where the exchange energy at Q = 0 is also subtracted. For practical calculation, the summation in equation (21) is usually taken over a dozen shells in the reciprocal space for sufficient convergence of the energy. 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ν = 1 inside a disc [10,24]. In coordinate space this can be represented as: with r b the radius of the disc [see Fig. 4(a)], which will be determined later. Transforming it to momentum space, we find where J 1 is a Bessel function. The normalization requirement ρ b (0) = M leads to r b = √ 2M l. So for the disc density profile, the energy per particle is 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ν = 1 andν = 0 unidirectional stripes along the periodic direction [10]. The distribution is written as whereê 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ν = a/Λ s . Its energy per particle is given by: The stripe period is a variational parameter fixed by minimizing 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]: Looking at equation (2), we deduce that the zeros of the direct potential are scaling as 1/ √ α in the q xdirection and √ α 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 √ α 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. 3(b). 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 parametrized in Fig. 4(a). The lengths l a and l b satisfy 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ν, 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 lowestenergy 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 abovementioned 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 √ α and 1/ √ α 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ν 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. 2(b) 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. 2(b) is an artificial result of the simple HF ansatzes.

C. 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.
DMRG optimises the matrix product state (MPS) variational wavefunctions [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 optimized 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. 5(a)).
Using Landau gauge, the orbitals (single-particle eigenstates) of the kinetic Hamiltonian [equation (16)] are localised and aligned along the axis of the cylinder and labeled by the orbital centre/ canonical momentum (2πn/L, n ∈ 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 (=ν). 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 (Fig. 4 and 6). Define N R as the number of unit cells wrapped around the cylinder circumference. We have N e = 2M N R and thus N u = 2qM N 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 2M q; 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ν = 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.

D. 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. 39, the translational symmetry The unit of energy is e 2 / 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 maximizes 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 NR = 4 and 1/6 filling. (b) Determining lattice shape and energy density by searching optimal Lop 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 Lop to estimate that of the infinite 2D system. The fold of "rotational" symmetry NR 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 lR can be determined as NRlR = Lop. The plot shows the calculation for systems at filling 1/6 and the isotropic limit. We calculate systems with three different NR. 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 Lops deviate from L * s by small values. Accordingly, the estimated lR deviates from ideal values l * R . The deviation is around 1.5% for NR = 4, 0.7% for NR = 5 (short); around 0.02% for NR = 3 (long). Also, we see the estimated energy densities differ by ∼ 10 −5 .
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 onedimensional 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 log(N 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 utilize 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 y (light) x (heavy) x (heavy) Increase from 1 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). 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 −6 . 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.

DATA AVAILABILITY
Additional results supporting the findings of this study are included in Supplementary Information. The data that supports the plots within this paper is available from the corresponding author on request.

CODE AVAILABILITY
Code is available upon reasonable request. Two-dimensional magnetotransport in the extreme quantum limit.  S2. The energy density of the metastable C lattice with reference to the energy density of the stable S Wigner/stripe lattice for anisotropy mx/my = 5, extracted from the DMRG calculation. The unit of energy is e 2 / l. The error bars are smaller than size of points. On the right side of the filling 1/6, the slope turns to be larger. The turning point coincides with the "boundary" between the Wigner and stripe region of the S lattice. In the small filling limit, S lattice and C lattice all turn into triangular lattice and the energy difference should vanish.
especially the situation of the stripe crystal. For the S lattice, one already sees in the main text that it eventually transforms to the stripe phase. For the C lattice, in isotropic situation, it should be linked to the stripe along the y-direction. But for highly anisotropic case, it tends to form a stripe along the x-direction with π/q * 1 as the periodicity. In Fig. S1, one can clearly observe that the stripe period converges to q * 1 when the electrons are dense enough for all mass ratios considered.

Supplementary Note 2: DMRG study of metastable phases
With the control of N u and L, we can obtain not only states representing the stable phase (e.g., S lattice) but also the metastable phase (e.g., C lattice, 3e bubble phase, and stripe crystal (ν = 2/9) in the isotropic limit) if the energy difference is small enough. This is possible because, for an infinite cylinder geometry, a state representing the metastable phase fits in a different optimal geometry (L and N u ) from those of the states representing the stable phase. In this case, the state representing the metastable phase can be indeed the ground state of that geometry, once it has a lower energy than the stable-phase states which are artificially deformed by that geometry. Such an artificial result disappears for the large L limit, because the deformation brought by the embedding geometry vanishes for any large enough L, even if L is not optimal.
We estimate the energy difference between the stable S lattice and the metastable C lattice in Fig. S2.
Supplementary Note3: Difficulty to determine weak density modulations The goal is to explain why our current data analysis cannot conclude whether a strict unidirectional stripe phase exists in the phase diagram.
Recall that a stripe crystal has weak density modulation along the stripe; the modulation vanishes for unidirectional stripe. We defineρ ⊥ (ρ ) as the largest ρ(q) with a non-zero wave vector along the direction perpendicular (parallel) to the stripe direction. The ratioρ /ρ ⊥ is zero in the unidirectional limit. It is not easy to estimate whether some state is exactly unidirectional. Firstly, DMRG only works with half-infinite systems, and the density modulation ratio of infinite systems in principle needs extrapolation. Secondly, DMRG works within finite bond dimension approximation, with the overestimation of order decays over increasing bond dimension. As discussed in the Methods, we need to align the direction with weak density modulation along the axial direction. It is possible that the modulation is a pure finite bond dimension error. For data with limited bond dimension, this possibility cannot be distinguished from the possibility of a very weak modulation. This is the case near the half filling, see the data ( Fig. S3) with bond dimension χ = 1200. We see the modulation along the stripe direction is very close to zero at half filling. We also see the ratio changes fast overν nearν ∼ 0.2.