Direct measurement of discrete valley and orbital quantum numbers in bilayer graphene

The high magnetic field electronic structure of bilayer graphene is enhanced by the spin, valley isospin, and an accidental orbital degeneracy, leading to a complex phase diagram of broken symmetry states. Here, we present a technique for measuring the layer-resolved charge density, from which we directly determine the valley and orbital polarization within the zero energy Landau level. Layer polarization evolves in discrete steps across 32 electric field-tuned phase transitions between states of different valley, spin, and orbital order, including previously unobserved orbitally polarized states stabilized by skew interlayer hopping. We fit our data to a model that captures both single-particle and interaction-induced anisotropies, providing a complete picture of this correlated electron system. The resulting roadmap to symmetry breaking paves the way for deterministic engineering of fractional quantum Hall states, while our layer-resolved technique is readily extendable to other two-dimensional materials where layer polarization maps to the valley or spin quantum numbers.

. Capacitance and dissipation for the data set in Fig. 2 of the main text. Supplementary Fig. 3. CA measured in a different device at B=35T for −4 < ν < 0. The device has nearly identical geometry to that described in the main text, but a larger area, enabling higher resolution measurements. Several fractional quantum Hall features are visible at one third filling of different LLs.

SUPPLEMENTARY NOTE 1: ELECTROSTATIC MODEL OF BILAYER GRAPHENE CAPACITANCE MEASUREMENTS
We model the dual gated graphene bilayer as a four plate capacitor, with the c i corresponding to the geometric capacitances as indicated Fig. 1A of the main text. The n i denote the areal electron densities on the four plates. Equations for the charge stored on each capacitor plate, as well as overall charge neutrality, result in four equations, Note that in the language of the main text, φ 2 − φ 1 = u. Suppl. Eqs. 1-4 are supplemented by the condition of electrochemical equilibrium between the top and bottom layers of the bilayer, where v 0 is the voltage applied to the bilayer, and µ i is the chemical potential on layer i. The µ i depend on both n 1 and n 2 through the constitutive relations that derive from the electronic structure of the bilayer. Capacitance measurements are performed with a small AC signal applied to one of three terminals while the corresponding variation in charge density is read out on another terminal. For small variations, then, the differential versions of Suppl. Eqs. 1-6 are relevant. In particular, Suppl. Eqs. 5,6 can be expressed in terms of the inverse compressibility matrix of the bilayer itself, κ ij = ∂µ i /∂n j , δφ 1 = δv 0 − κ 21 δn 1 − κ 22 δn 2 (7) δφ 2 = δv 0 − κ 11 δn 1 − κ 12 δn 2 (8) Note that κ 12 = κ 21 follows from a Maxwell relation. Experimentally, we measure the elements of the capacitance matrix Where the indices indicate the voltages applied to the top gate (v t ), bottom gate (v b ), or bilayer itself (v 0 ). Three elements of the capacitance matrix are independent, and we choose the most directly experimentally relevant combinations: the penetration field capacitance C P ≡ −C BT = −C TB , and top and bottom gate capacitances C B ≡ C B0 and C T ≡ C T0 . Expressions for these three quantities can be found by varying Suppl. Eqs. 1-4 and using Suppl. Eqs. 7-8 to eliminate φ 1 and φ 2 . All measurable capacitances depend on all three components of the compressibility matrix, As described in the main text, the ultimate quantities of interest are the total density and layer density imbalance of the bilayer, n ≡ n 1 + n 2 and p = n 1 − n 2 , while the most natural control parameters are The partial derivatives of n with respect to n 0 and p 0 follow trivially from the fact that partial derivatives of the n with respect to the gate voltages can be measured directly: Supplementary Fig. 4. Bilayer graphene hopping parameters.
For convenience we have introduced the average geometric capacitance c = c b +ct 2 and the geometric capacitance asymmetry between top and bottom gates δ = c b −ct c b +ct ; we also define here the capacitance observables measured and described in the main text as Derivatives of p can also be computed by varying 1-4 and using Suppl. Eqs. 7-8, to give expressions in terms of the κ ij . After solving Suppl. Eqs. 10-12 for the three κ ij and substituting the results, we arrive at expressions for these derivatives in terms of observable capacitances (C P , C S , and C A ): In our device the top and bottom gate geometric capacitances are nearly symmetric with δ ≈ .029; in addition, we estimate that c 2c0 .0086. Taking the leading order terms in these two small parameters, we finally arrive at simplified expressions, some of which are used in the main text. As described in the main text, calculating a phase diagram of bilayer graphene that correctly captures the experimentally observed effects requires considering a variety of both single particle and correlated electron effects. In this section we describe a minimal model which accounts for these effects, and compute ground state energies using Hartree-Fock and DMRG in order to determine the ν-dependence of the phase boundaries. The presentation is organized as follows: The Hamiltonian is a sum of the single-particle and two-body terms, We first describe the terms we choose to include in H. To make numerical progress, we then project the problem into the ZLL. Finally, in order to compute the energies of the competing phases at arbitrary ν, we perform both Hartree-Fock and DMRG computations.
We start with a tight binding model of bilayer graphene that includes three intersite hopping terms (γ 0 , γ 1 , and γ 4 as well as the dimer on-site energy (∆ ) but neglect the trigonal warping term (γ 3 ). A diagram of the bilayer graphene structure with relevant hopping integrals is shown in Fig. 4). At low energies, the Fermi surface has two disconnected parts near the corners of the Brillouin zone around the two inequivalent K points. The single particle Hamiltonian can then be reduced to a four-band model corresponding to the four sites of the BLG unit cell, resulting (in, e.g., valley K) in where the basis in the K valley consists of the wavefunction weight on the four lattice sites in the bilayer graphene unit cell Here π = p x − ip y and π † = p x + ip y are momentum operators, and u = φ 2 − φ 1 is the potential difference across the bilayer induced by the perpendicular electric field. Velocities are defined in terms of the monolayer graphene lattice constant, a = 2.46Å, as v 0 = To extend this Hamiltonian to the case of large perpendicular magnetic field, we introduce creation and annihilation operators for the scalar Landau level wavefunctions localized on each lattice site, defined asâ ≡ B (q x − iq y ) andâ † ≡ B (q x − iq y ) where q i ≡ k i − e c A i and A is the magnetic vector potential. The operators operate on scalar Landau level wavefunctions such thatâ|n = √ n|n − 1 andâ † |n = √ n + 1|n + 1 . The Hamiltonian in valley K, for example, then becomeŝ where the monolayer graphene cyclotron energy is ω 0 = v0 where ξ labels the valley (henceforth denoted ξ = ±), N labels the orbital quantum number, and the |n are the oscillator states of a, equivalent to the conventional quadratic-band LL-wavefunctions. The coefficients are then determined by the band-structure.
B γ1 γ 0 , while the N = 0 and 1 orbital are nearly degenerate, both having zero energy for u = γ 4 = ∆ = 0. γ 4 , ∆ , and finite u all weakly lift this degeneracy, but still leave an eight fold near-degeneracy between states of different orbital, spin, and valley quantum numbers. The ZLL is well separated from the N = ±2 states at E ±2 ≈ ± √ 2 ω c ≈ ±113meV. While N ≥ 2 LLs have support on all four sublattices in the unit cell, the N = 0, 1 eigenstates vanish on one or more sublattices: Defining the layer polarization as α ≡ |c A | 2 − |c B | 2 + |c B | 2 − |c A | 2 , we note that within the 4 band model the N = 0 orbital is fully layer polarized (α 0 = 1) but the N = 1 is not; using the tight binding parameters above we find α 1 ≈ .63 at B ⊥ = 31T. Wavefunctions in the opposite valley have correspondingly opposite layer polarization. Throughout our experiment u/ ω c 1, and to leading order the single-particle energies of the ZLL are Here, as in the main text, E Z is the Zeeman energy, σ = ± 1 2 denotes the spin projection along the direction of the applied field, ∆ 10 ≈ ω c (2 γ4 γ0 + ∆ γ1 ) is the single particle orbital splitting, u = φ 2 − φ 1 is the E-field induced potential difference across the bilayer, α N is the layer polarization of the orbital, and ξ = ± indexes the valley. The single particle levels are shown in Fig. 3A of the main text. At B = 31T, ∆ 10 = 9.7 meV, E Z = 3.58 meV, and {α 0 , α 1 } = {1, 0.63} as follows from Suppl. Eq. 28. The large splitting to the higher Landau levels ensures that they are not involved in any uor E Z -tuned phase transitions.

Coulomb Hamiltonian
The Coulomb interactions decompose into a dominant isospin SU(4)-symmetric part and subleading capacitive and valley anisotropies: We now discuss these terms in turn.
Screened Coulomb interaction. The bare Coulomb interaction is screened by the surrounding hBN dielectric, the proximal metallic gates, and filled LLs below the ZLL of the BLG itself. The screening due to the hBN dielectric is incorporated as a dielectric constant in the Coulomb scale assuming BN ≈ 6.6 0 2,3 . The metallic gates, each at distance D ≈ 20nm from the bilayer, exponentially screen the interaction when r D. The 2D Fourier transform of the gate-screened potential is in units of E C and B . Here we neglect the finite width of the bilayer itself, which is an order of magnitude smaller than B ; it will be reincorporated as a capacitive energy below.
Since we wish to work with a model projected into the ZLL, we must account for the residual response of the other LLs, colloquially referred to as "LL-mixing". The dimensionless parameter controlling their response is E C ω0 ∼ 3.14/ B/T ≈ 0.56 at B = 31T, which is comparable to values in GaAs. Motivated by the large number of isospin flavors (four), the standard approach for BLG is the random phase approximation (RPA) , 4,5 in which we replace the bare Coulomb potential V (k) with the effective potential where Π is the polarization response. The RPA result is then further approximated by the static ω = 0 value. However, RPA calculations have only been reported for the PH-symmetric two-band model 5 at ν = 0, which is not quantitatively correct at large electric fields or at the magnetic fields relevant for our experiment 6 . Moreover, Ref. 5 found that the static ω = 0 approximation strongly overestimates screening (IQHE gaps were underestimated by a factor of three). Thus even recalculating the RPA value with a four band model is unlikely to be quantitatively accurate, as it is not possible to incorporate ω dependence into ground state numerical methods like exact diagonalization or DMRG. Thus, at present, there does not appear to be a satisfactory "ab initio" tool for quantitatively predicting the strength of screening.
For these reasons, we use a phenomenological model following the approach of Ref. 7, taking This is motivated by the RPA form when approximating the polarization as Π(k) = a tanh(bk 2 2 B )4 log(4)/2π, with V (k) given in Suppl. Eq.(35). If a, b are chosen to match the low-k and high-k behavior of the ν = 0 two-band RPA calculation, 5 one finds a RPA ≡ E C ω0 and b RPA ≡ 0.62. However, much more generally we must have Π(k) ∝ k 2 at low-k and Π(k) → const at high-k, as captured by the ansatz. Quantum Hall calculations are only sensitive to the form of the interaction in the vicinity of k < ∼ −1 B , so the magnitude of the low-k behavior (k 2 ) forms a one-parameter space of screening behaviors set here by the product ab. Thus we fix b = b RPA , and treat a scr as a phenomenological measure of the screening strength.
ZLL projected Hamiltonian-We then project the screened Coulomb interaction V eff (q) into the eight components of the ZLL. For the moment we neglect the small, lattice-scale valley anisotropies which break the valley-SU(2) symmetry; these effects will be introduced as phenomenological couplings shortly. The SU(4)-isospin symmetric interaction is Here n ZLL (q) = ξσ n ξσ (q) is the total density in the ZLL, which is a sum of the four isospin components, whileN ξ1σ is the electron number in level ξ1σ. As explained in Ref. 8, a shift ∆ Lamb between the N = 0, 1 orbitals arises when projecting into the ZLL, since their Coulomb exchange with the filled LLs below the ZLL differs. Under the approximation of particle-hole symmetry, it was shown that N N is the Coulomb exchange per-electron when fully filling an N level with a Slater-determinant. We evaluate this shift using the screened interaction V eff , and find a near perfect fit to ∆ Lamb − 0.2E C 1+2.73ascr . The density n ξ,σ (q) contains a contribution from both the N = 0 and N = 1 levels. The ratio between the Coulomb scale and the single-particle splitting between the orbitals is Thus, unlike the levels outside the ZLL, it is not well-justified to project the interaction into only one of the two N orbitals; we must keep both. BLG form factors-For completeness, we explain how the BLG "form factors" F N M can be used to compute the Coulomb matrix elements projected into the Landau-level basis, which are required for the actual computations. The density n σξ is not diagonal in the orbital index N , but instead involves "orbital-mixing" contributions: Hereρ is a guiding-center density operator, which in the Landau gauge reads The form factor F is expressed in terms of the conventional quadratic-band form factors F nm as In the four-band model, we refer to Suppl. Eqs.(30)-(31) to find It is thus convenient to parameterize the interaction by At low perpendicular magnetic fields, Θ → 0 and the problem reduces to the two-band model equivalent to the two lowest Landau-levels of the conventional quadratic-band QHE.
Lattice scale effects at order a/ B generate valley-SU(2) breaking perturbations sensitive to the details of the lattice and orbital structure, so must be treated phenomenologically 9,10 . The first anisotropy is a capacitive energy due to the finite thickness of the BLG (d), which modifies the interlayer Coulomb interaction from V (q) → V (q)e −qd . This perturbation is smaller than the long range Coulomb (H SU (4) ) by a factor of d/ B 1, but the q = 0 part is unscreened, so for simplicity we retain the Hartree-type capacitive charging energy The ratio of dielectric constants arises because we included BN in E C . The dielectric constant ⊥ BLG should be the same one used for converting the applied gate voltage p 0 /c to the inter-layer bias u.
Following Refs. 9 and 10, the remaining valley anisotropies are incorporated as short-range (δ-function) interactions which preserve spin SO(3) and valley U (1) × U (1) Z 2 . In contrast to monolayer graphene, symmetry considerations allow a rather large space of possible perturbations since they can take an arbitrary form in the N = 0, 1 orbital index. Here we will assume the perturbations can be expressed using local operators that depend only on the isospin, e.g. O µν (r) = ψ † ξ,σ (r)τ µ ξ,ξ σ ν σ,σ ψ ξ ,σ (r), where ψ ξ,σ = N ψ ξ,N,σ (r) and τ, σ are Pauli operators. This assumption is reasonable when Θ → 0, since the anisotropies arise from lattice structure and the N = 0, 1 orbitals sit on the same site in this limit. While this approximation neglects the Θ = 0 corrections due to the finite polarization of the N=1 levels, we note that the capacitance term H c0 captures at least some of these.
We express the anistropy energy using the total valley number density n ± (r) and the valley spin density S ± (r) = 1 2 M N ab ψ † ±M a (r)σ ab ψ ±N b (r). In addition to the total density (which preserves SU(4)), the most general symmetry preserving perturbation is Note that previous treatments 9 have assumed interactions of the type u ⊥ 2 (τ x τ x + τ y τ y ) and uz 2 τ z τ z . These forms are in fact equivalent: using the anti-symmetry of the fermions, we have E C

Evaluation of H (2)
Our model depends on a number of parameters, summarized in Tab. I. Some are known from the literature (e.g., BN and the tight binding parameters shown in Suppl. Eq. 27), some can be derived directly from these assumed values (e.g., ∆ 10 , α 1 , E C ), and some follow from theory (e.g., ζ, ∆ Lamb ). In addition, there are several phenomenological parameters which we constrain from experiment, namely ⊥ BN , g ⊥ , g z , a scr , ⊥ BLG . To get a handle on these parameters, we begin with a Hartree Fock approximation for H (2) , which allows us to conveniently estimate the predicted location of the transitions within our model.

Hartree Fock approximation for H (2)
Within the Hartree-Fock approach, we assume that the ground-state at integer filling is a Slater determinant which successively fills orbitals ξN σ according to the proposed phase. We neglect the possibility of valley off-diagonal coherence (e.g., non-zero ψ † ξ ψ ξ = 0, ξ = ξ ). While such phases have been predicted to occur, for instance in a very narrow range of p 0 at ν = −3, 11 there is no evidence for them in the current data, since there appears to be a single direct transition at p 0 = 0.
The energy of H SU(4) at integer filling can then computed as an integral involving F N M and V eff (see, for example, Ref. 11). In contrast to earlier works, we use both the screened interactions and the form factors appropriate to the B = 31T four-band model, as described above. We interpolate the integer result to fractional ν by generalizing the first term in the interpolation of Fano and Ortolani 12 to the multi-component case: ). Calculations were repeated for a range of a scr , and we find gives an almost perfect interpolation of the result. The Lamb shift, meanwhile, is found to be Note that this expression for the Lamb shift remains true beyond Hartree-Fock, so will also be used in our DMRG calculation. The valley-anisotropies are off-diagonal in the valley index, so they reduce to a Hartree energy, e.g. n + (r)n − (r) → n + (r) n − (r) . Using n ± (r) = 1 2π 2 B ν ± , we arrive at the final expression Since ⊥ BN determines the capacitance between the gates and BLG, we can measure ⊥ BN by fitting the Landau level spacing, as measured in the applied gate voltages n 0 /c, to their known densities n = ν 2π 2

B
. Data were taken at p 0 = 0 and B ⊥ = 2T . Starting from the electrostatic model of Suppl. Eqs. 1-2, we ignore interlayer capacitance (c 0 ∼ ∞) to treat the bilayer as a single 2D electron system, and neglect the finite quantum capacitance, which is reasonable at low fields. Fitting the separation, in n 0 , of two four fold degenerate LLs to , we find ⊥ BN ∼ 3.0 ± .15.
Experimental estimate of g ⊥ By measuring the critical total B-field of the ferromagnet to canted antiferromagnetic transition at ν = p 0 = 0, previous experiments 13 have estimated that B tot * = (2.42 ± 0.21)B ⊥ . This implies u ⊥ = 0.14meV/T, or in our parameterization, 2. While this experiment shows the expected B ⊥ -linear scaling for B ⊥ = 4 − 7 T, we must be careful when extrapolating to higher magnetic fields.
To address this question, we analyze the tilted field measurements of four phase transitions at B ⊥ = 15T, as shown in Fig. 6. From the experimental data, we extract the shift in ∆p 0 /c over the probed magnetic field range. Each transition is expected to  Fig. 6 and Suppl. Eq. 51 ν Transition Canting? shift by ∆u = σ∆E Z α , where α is the effective change in layer polarization across the transition and σ is the effective change in spin, and ∆E Z = 0.116∆B tot meV/T the change in Zeeman energy. Solving for σ, we find where 0.0086 = 0.335nm 39nm is the ratio of geometric capacitances. When calculatingᾱ, we must remember that α 1 = 0.8 at B ⊥ =15T, as extracted from band structure parameters 1 . In the absence of antiferromagnetism, all transitions considered involve reversal of one full electron spin, so that σ = 1. Canting of the spins due to g ⊥ can reduce the effective spin; however, σ > 1 is unphysical as it would imply a larger than unity spin per electron. Taking into account experimental error, each tilted field data point thus imposes a lower limit on ⊥ BLG , the most stringent of which is ⊥ BLG > 2.57. As we will see below, this limit is ultimately consistent and does not influence the final fitted value of ⊥ BLG = 2.76. Using this value for ⊥ BLG , we can then calculate σ for each transition, shown in Table II. Two of the four spin transitions are consistent with absence of canting (σ = 1), and two are consistent with strong canting (σ < 1). Additional data (not shown) shows that the phase transitions at ν = −2 do not depend on in-plane magnetic field, suggesting that antiferromagnetism has not yet set in at this filling.
In a phase characterized by net spin S + and S − in each valley, the threshold for canting in our model arises from the interplay between the Zeeman effect (which favors polarization) and g ⊥ (which favors canting) via where observation of σ < 1 at a phase transition implies that E Z < E c Z for one of the adjacent phases. The canting at T 0< and the absence of canting at both T −2 and T 0> together constrain 0.52 < g ⊥ < 1.04. T −1 and T +1 , meanwhile, seem to give contradictory constraints, 0.69 < g ⊥ and 0.69 > g ⊥ . The lack of consistency likely stems from the dependence of the g ⊥ anisotropy on the orbital filling: T −1 and T +1 involve different partial fillings of an isospin flavor, and g ⊥ need not be identical for different orbital combinations in the different valleys. At T −1 , occupation is transferred from a |0+ ↑ to a |0− ↓ state, while at T −1 occupation is transferred from a |1 + σ to a |1 − σ state. As we discussed, this effect is not captured in Suppl. Eq. 46, which treats orbital components equally. While the precise value and possible orbital substructure of g ⊥ is critical for determining the precise dependence of the antiferromagnetism on filling, within our model we find g ⊥ has negligible impact on the locations of the u-tuned phase transitions, which are the focus of this work. Thus we set g ⊥ = 0.69, and defer a more full analysis of the ν dependence of the antiferromagnetism to future high resolution studies of the tilted field dependence of the phase transitions observed in C A .
ν ξN σ of a filling sequence, the phase boundaries u * (ν) are determined by equating E (HF) The layer bias u * is then related to the experimentally measured gate bias via p0 0.335nm . Experimental data points are provided by the positions, in p 0 /c, of the transitions T −2 , T −1 , T 0> , and T +1 (see Fig. 5 for the labeling scheme). In the absence of canting (we have assumed g ⊥ = 0.69), the predicted phase boundaries are To compare to experimental values of p 0 /c for each transition, we convert V (T i ) = u * (T i ) 2c0 c . An fifth constraint arises from the tilt-field dependence of Suppl. Eq. (54), which implies Given five experimental constraints in three unknowns, we do a least-squares fit for our overconstrained model. The model is validated, simultaneously fitting all data points across a wide range of ν, and providing physically reasonable values for ⊥ BLG = 2.76, a scr = 0.29 = 0.5a RPA , and g z = .082. The resulting phase diagram is plotted in Fig. 3c of the main text, lower panel, and similarly shows good quantitative agreement with most of the features of the experimentally measured phase diagram.

Phase boundaries from iDMRG
In order to reproduce the kink in u * (ν) observed in experiment around −3 < ν < −1, we replace the Fano-Ortolani interpolation of the HF result (Suppl. Eq. (47)) with the infinite-DMRG calculation, which takes full account of correlations. As before, we must compute the energies of the competing phases. Since the tilt-field dependence shows that the spins are polarized across the transition, we label the isospin only by its valley ξ = ±. At the transition, density is transferred between isospin components. The cartoon picture is that in phase "001" (low p 0 ), the orbitals fill in order +0, −0, +1, while in phase "010" (higher p 0 ), the fill in order +0, +1, −0. More generally, valley U(1)×U(1) symmetry allows us to assign separately conserved fillings ν + , ν − to the two valleys. The two competing phases are defined by their fillings ν + , ν − (we writeν = ν + 2 in this regime): (ν + , ν − ) 010 = (2 +ν, 0)ν < 0 (2,ν)ν > 0 (58) Supplementary Fig. 7. The densities ν±,N = 2π 2 B n±,N (we assume spin-polarization for simplicity) as the total filling ν is increased. p0 < 0 is small, so the orbitals fill in the order −0, +0, −1, +1. We see that while charge does fluctuate between N = 0, 1 orbitals of the same isospin, this deviation is at most about 10% of the filling. This explains why blue/red (N = 0) and cyan/orange (N = 1) appear as a sharp contrast in our CA data.
Note that for fractional ν very close to the transition, there are presumably a multitude of phases in which fractional filling has transferred between the valleys, which we do not consider above. However, the current experiment isn't sensitive to these more delicate states, which are characterized by a lower energy scale, so they will appear as a small rounding of the phase boundary u * .
The iDMRG conserves U(1)×U(1), so can be used to find the lowest energy state at filling (ν + , ν − ) under the Hamiltonian H = H (1) + H (2) . While in principle the full Hamiltonian could be simulated in DMRG, for technical reasons it greatly simplifies matters to decompose the Hamiltonian as The dominant part, [·] DMRG , will be evaluated in DMRG. The u-dependence [·] P.T. is evaluated in first order perturbation theory, by evaluating N ξ,N using the ground states found in DMRG. When α 1 = 1 first order perturbation theory is exact, since the layer polarization commutes with the Hamiltonian, and our tests indicate that more generally this introduces negligible error. We make this approximation so that u * (ν) can be determined from a single DMRG run at u = 0, rather having to re-run DMRG for each value of u. Finally, [·] H.F. is evaluated by neglecting inter-valley correlations, e.g. by taking n + (r)n − (r) → n + (r) n − (r) : Since one of the two isospins is always at integer filling, and hence is largely inert, the inter-valley correlations are expected to be small (on top of the already small scale d B ). For example, at ν + = ν − = 1, DMRG shows the pair correlation is around n + (r)n − (r) ∼ 0.8 n + (r) n − (r) , a slight suppression from the Hartree value. Thus while largely negligible for determining the phase boundary u * , effects of this form could be relevant for understanding the smaller energy scale governing spin physics, an interesting subject for future work.
All parameters except a scr are those used or determined by the Hartree-Fock analysis. As will be discussed, we find that a scr = 0.22 must be adjusted slightly to preserve the location of the ν = −2 transition. iDMRG numerics. In the infinite-DMRG method we place the quantum Hall problem [·] DMRG on an infinitely long cylinder of circumference L, for which we compute the ground state energy density including the full effect of correlations. When both 2 > ν + , ν − > 0 (for instance, in phase 001), it is necessary to keep four ZLL components in the iDMRG, +0, +1, −0, −1. This is because even when ν ξ = 1, isospin ξ acts as a polarizable medium due to the very small splitting between the N = 0, 1 orbitals. While computationally very expensive, 14 the ZLL orbital mixing is thus fully accounted for.
In principle H has a delicate ν-dependence sensitive to all the fractional competing phases (which may be distinguished at the level of 10 −3 or 10 −4 E C ), which would require finite-scaling analysis to fully resolve. However, given the resolution of the present experiment, we focus on the much larger and slowly varying background (at the level of 10 −1 E C ). For this purpose, we work on cylinders of circumference L = 16 B and use a DMRG-bond dimension of χ = 1600, which results in an error in the energy per particle of around 10 −4 E C , much smaller than the experimental features to be modeled. Specifying the valley fillings ν + , ν − according to Suppl. Eq. (59), iDMRG was used to compute the energy under [·] DMRG for the two competing phases at 1 +ν ∈ {0, 1/5, 1/3, 2/5, 1/2, 3/5, 2/3, 4/5, 1, · · · , 2}. Before discussing the resulting DMRG energies, we first address the issue of why blue/red and cyan/orange appear as distinct scales in the experimental C A data, which requires that the N = 0, 1 levels to largely fill sequentially rather than as a mixture. Consider, for instance, the low p 0 < 0 phase in which (naively) the orbitals fill in the order −0, +0, −1, +1. With interactions, charge can fluctuate between the N = 0, 1 orbitals of the same valley. In Fig. 7, we show the iDMRG result for the densities ν ±,N = 2π 2 B n ±,N as the total filling −4 < ν < 0 is varied. The deviation from the non-interacting expectation is at most around 10%, explaining the sharp contrast.
In Fig. 8a, we show the interacting part of the energy E SU (4) for both phases. The difference ζ = E SU (4) 010 (ν = −2) − E SU (4) 001 (ν = −2) is slightly smaller than its Hartree-Fock value, so to preserve the location of the transition we adjust the screening to a scr = 0.22 (this adjusts ∆ Lamb accordingly). In Fig. 8b, we show their difference after extracting out the linear part, f (ν) = E SU (4) 010 (ν) − E SU (4) 001 (ν) − ζ(1 − |ν + 2|). The significantly different curvature on the two sides of ν = −2 leads to the kink in the phase boundary discussed in the main text.
After calculating the DMRG and HF part of the energy at u = 0, the remaining polarization energy is determined from the DMRG expectation values u 2 ξα N n ξN shown in Fig. 7. From this we compute the boundary u * (ν) shown in Fig. 3e-g.