Leonardo Benini

Piero Naldesi

Rudolf A Römer

Tommaso Roscilde

Loschmidt echo singularities as dynamical signatures of strongly localized phases

Keywords: Localization, Disordered systems, Non-equilibrium dynamics, Entanglement

Quantum localization (single-body or many-body) comes with the emergence of local conserved quantities -whose conservation is precisely at the heart of the absence of transport through the system. In the case of fermionic systems and S = 1/2 spin models, such conserved quantities take the form of effective two-level systems, called l-bits. While their existence is the defining feature of localized phases, their direct experimental observation remains elusive. Here we show that strongly localized l-bits bear a dramatic universal signature, accessible to state-of-the-art quantum simulators, in the form of periodic cusp singularities in the Loschmidt echo following a quantum quench from a Néel/charge-density-wave state. Such singularities are perfectly captured by a simple model of Rabi oscillations of an ensemble of independent two-level systems, which also reproduces the short-time behavior of the entanglement entropy and the imbalance dynamics. In the case of interacting localized phases, the dynamics at longer times shows a sharp crossover to a faster decay of the Loschmidt echo singularities, offering an experimentally accessible signature of the interactions between l-bits.

Introduction

Constructive interference of paths bringing a particle back to its initial location in real space is at the heart of single-particle (or Anderson) localization (AL) [START_REF] Akkermans | Mesoscopic Physics of Electrons and Photons[END_REF][START_REF]Anderson Localization and its Ramifications[END_REF]; more recently a similar phenomenon occurring in Hilbert space (MBL) [START_REF] Nandkishore | [END_REF][4][5][6] has been shown to prevent many-body quantum systems from relaxing to thermal equilibrium, undermining the ergodic hypothesis in a large class of models of interacting quantum particles. Localized phases are generally characterized in the negative (absence of transport, of long-range order, of spectral gaps, etc.), while positive characterizations are typically elusive. A crucial aspect of localization is the persistence of initial conditions, which, in the case of AL of non-interacting particles, is related to the conservation of populations in the localized single-particle eigenstates of the Hamiltonian. In the case of MBL, the analog phenomenon would be the appearance of local conserved quantities (called local integrals of motion or l-bits [7][8][9][10]) which are obtained by unitary transformations of local operators; and which, if extensive in number, constrain the dynamics of the system to the point of preventing relaxation.

The existence of l-bits in disordered spin chains can be mathematically proven under the assumption of limited level attraction [9], and approximate l-bits for manybody systems can be constructed with a variety of analytical as well as numerical methods [11][12][13][14][15][16][17][18][19][20]. Much of the phenomenology of MBL dynamics (persistence of traits of the initial state, logarithmic growth of entanglement entropies, etc.) -observed in numerical studies as well as in experiments [4,5] -can be directly explained in terms of the existence of l-bits and interactions between them. Yet observing l-bits directly is an arduous task, given that their expression is highly disorder-dependent (and generally unknown even in theory), and it would require high-precision measurements of local observables in different local bases. An even more ambitious task for experiments is the one of probing directly the existence of interactions among l-bits, which is a defining feature distinguishing MBL from AL. Measuring the consequences of such interactions on correlation and entanglement dynamics is currently the focus of a considerable experimental effort based on state-of-the-art quantum simulation platforms [21][22][START_REF] Chiaro | Direct measurement of non-local interactions in the many-body localized phase[END_REF].

The purpose of this work is to show that, in the case of strongly localized phases, the existence of l-bits can offer striking signatures in the dynamics of the Loschmidt echo, namely in the logarithm of the return probability to the initial state |ψ 0

λ(t) = - 1 L log | ψ 0 |e -iHt |ψ 0 | 2 av . ( 1
)
Here H is the system's Hamiltonian and L the lattice size; [...] av indicates the disorder average. Please note that we employ the "Loschmidt echo" terminology similar to other recent works in which singular behavior has been observed in the quench dynamics of many-body quantum systems [START_REF] Heyl | [END_REF][25][26][27][28][29][30][31]. Nonetheless, In the quantum-chaos literature the Loschmidt echo is more generally defined as the scalar product between the evolution of the same state |ψ 0 with two different Hamiltonians, H 1 and H 2 [32,33]. Starting from L(t) = ψ|e iH 2 t e -iH 1 t |ψ , our definition is retrieved in the case of H 1 = H and H 2 = 0.

When |ψ 0 has a simple factorized form, and in the case of strong disorder, we find that the Loschmidt echo displays periodic singularities, decaying very slowly in amplitude -as illustrated using a model of disordered spinless fermions in 1d (corresponding to the S = 1/2 XXZ model in a fully random or quasi-periodic field) initialized in a charge density-wave (CDW) state. The singularities in the Loschmidt echo are fully explained quantitatively by a simple model of a collection of localized 2-level systems (2LS) undergoing independent Rabi oscillations, and approximating strongly localized lbits. The same minimal model captures quantitatively the dynamics of the entanglement entropy at short times as well as of the number entropy at longer times; and the dynamics of the density imbalance characterizing the initial state. At longer times the deviation of the exact results for the MBL dynamics from the predictions of the 2LS ensemble offers direct evidence of the interactions among the lbits in the form of a faster decay of the Loschmidt-echo singularities and imbalance oscillations. As the Loschmidt echo and the imbalance are generally accessible to quantum simulators, either measuring individual degrees of freedom [26,34] or even global ones [35][36][37], our results show that strong direct signatures of l-bits dynamics and interactions are within the immediate reach of state-of-the-art experiments on disordered quantum systems.

The structure of the paper is as follows: section 2 introduces the XXZ model in a random/quasi-periodic field; section 3 discusses the observation of Loschmidt echo singularities in the exact dynamics, as well as their quantitative understanding via a model of 2LS as well as three-level systems (3LS); section 5 discusses the dynamics of imbalance, and the comparison with the prediction of the 2LS and 3LS model; section 6 discusses the departure of the exact data (for the Loschmidt echo and imbalance dynamics) from the 2LS/3LS predictions as a signature of l-bit interactions; section 7 shows that the 2LS/3LS models capture quantitatively the entanglement dynamics at short time, and of number entropy at longer times; conclusions are drawn in section 8.

Model

Our platform for the investigation of Loschmidt-echo dynamics is given by a paradigmatic model, namely the S = 1/2 XXZ chain in an inhomogeneous magnetic field [38,39], corresponding to a model of spinless fermions with nearest-neighbor interactions in an inhomogeneous local chemical potential [40]

H = L-1 i=1 - J 2 S + i S - i+1 + h.c. + J z S z i S z i+1 - L i=1 h i S z i = L-1 i=1 - J 2 c † i c i+1 + h.c. + J z n i n i+1 - L i=1 h i n i , (2)
where S α i (α = x, y, z) are spin operators and c i , c † i and n i = c † i c i are fermionic operators; the equality between the two Hamiltonians is true up to an additive constant via Jordan-Wigner mapping. In the following the external field/potential h i is taken to be either quasi-periodic (QP) [41,42], namely h i = ∆ cos(2πκi + φ) with κ = 0.721 (inspired by experiments on bichromatic optical lattices [35,43]) and φ a random phase; or to be fully random (FR) and uniformly distributed in the interval [-∆, ∆]. We consider chains of length L (up to L = 22) with open boundaries, and we average our results over ∼ 10 3 realizations of the random phase (QP) or of the full random potential (FR). All the unitary evolutions considered in this study are obtained using exact diagonalization (ED), and they start from the charge-density wave state |ψ 0 = |1010101... , corresponding to a Néel state for the spins. We shall focus on the case of interacting fermions J z = J (corresponding to an SU(2) invariant spin-spin interaction) and contrast it with the limit of free fermions J z = 0. In the latter case, the QP potential leads to a transition to fully localized single-particle eigenstates for ∆ ≥ J, with an energy independent localization length ξ = 1/ log(∆/J); while the FR potential leads to AL of the whole spectrum at any infinitesimal value of disorder. In the interacting case, instead, a QP potential of strength ∆ 4J [44] and a FR potential of strength ∆ 3.5J [45][46][47][48][49][50][51][52] are numerically found to lead to MBL. In the following sections we shall generally start our discussion from the case of the QP potential, which has a simpler spatial structure devoid of rare regions, leading to stronger localization effects; and we shall later discuss how to enrich the picture in the case of the FR potential, in order to account for the existence of rare regions.

Loschmidt echo singularities and imbalance oscillations

Figure 1 shows the dynamics of the Loschmidt echo λ(t), (1), along with that of the imbalance

I(t) = 1 L i (-1) i (2[ n i ] av -1) . (3)
The latter saturates to its maximum value of 1 in the initial state and probes the persistence of the initial density/spin pattern [35][36][37]. We observe that for both the QP and FR potentials, and for disorder strengths compatible with the onset of the MBL regime, the Loschmidt echo displays a sequence of periodic cusp-like peaks at times t m = (2m + 1)π/J (m = 0, 1, 2, ...). These times correspond to minima in the imbalance, as the system reaches instantaneous configurations which are the farthest from the initial spin/density pattern.

A closer inspection shows that, for sufficiently strong disorder, all the peaks become sharp cusps, namely they represent genuine non-analyticities of the Loschmidt echo. They are rather remarkable given that they survive disorder averaging, and they seemingly appear in a finite fraction of disorder realizations (see Appendix C for further details); and in particular they decay very slowly in time, as we shall discuss in detail later on. The rest of this work will be devoted to developing a quantitative understanding of the dynamics of Loschmidt-echo singularities and imbalance oscillations as signatures of the existence of l-bits and of their interacting nature. 4. Quantitative modeling of the Loschmidt echo singularities

Modeling with an ensemble of two-level systems

All the essential details of the short-time evolution of the Loschmidt echo can be captured with a surprisingly simple, yet rather insightful model. This model is best understood (and justified) in the case of the QP potential, as illustrated in figure 2.

In the case of strong disorder, the fastest dynamics in the system starting from a Fock state will be offered by those particles that sit on a site i which is nearly resonant with its unoccupied neighbor (say i+1), because the hopping J/2 is either larger than the energy offset δ i = h i+1 -h i (in the non-interacting case) or larger than the screened offset δ i -J z (in the presence of nearest-neighbor repulsion). These 2-site clusters, representing nearly resonant two-level systems (2LS), have the property of being spatially isolated in the QP potential, because of the strong anticorrelation among two consecutive energy offsets (δ i and δ i±1 -see Appendix B). As a consequence, a nearly resonant 2-site system will be generally surrounded by highly non-resonant pairs of sites, which can be considered as nearly frozen to the initial state. This invites us to write for the evolved state a 2LS Ansatz

|ψ(t) ≈ ⊗ p |ψ (p) 2LS (t) ⊗ (⊗ i |ψ 0,i ) , (4)
where the first tensor product ⊗ p runs over the nearly resonant 2LS, while the second tensor product ⊗ i runs over the leftover sites (we have taken the freedom of reordering the sites arbitrarily in the tensor product; see [53][54][55] for a similar Ansatz to study the long-time dynamics). |ψ (p) 2LS (t) is the evolved state of the p-th (isolated) 2LS system, corresponding to two states split by an energy difference δ p = δ p -J z and connected by a Rabi coupling J; while |ψ 0,i is the (persistent) initial state of the site i belonging to the remainder of the system. The Loschmidt echo for such a system is readily calculated as

λ(t) = - 1 L p log 1 -p(δ p , J, t) , (5)
with p(δ,

Ω, t) = (Ω/Ω ) 2 sin 2 (Ω t/2) (and Ω = √ Ω 2 + δ 2
) the well-known probability of finding the 2LS in the state orthogonal to the initial one while performing Rabi oscillations [START_REF] Scully | Quantum Optics[END_REF] -see also Appendix A. When averaging Eq. ( 5) over disorder, it is immediate to obtain the following simple expression

λ 2LS (t) = -P (δ + J z ) log [1 -p(δ ; J, t)] , (6)
where P (x) is the probability that the energy offset between two neighboring sites takes the value x. Going from Eq. ( 5) to Eq. ( 6) implies that we in fact count all of the L pairs of sites in a chain as nearly resonant 2LS, thereby counting twice every site. The mistake that one makes in doing this is minor, though, because the non-resonant pairs of sites give a very small contribution to the Loschmidt echo; and, if neighboring pairs of sites are not simultaneously resonant, a site will not be counted twice in practice. Eq. ( 6) is an analytical integral formula which depends uniquely on the (known) statistics of the 1) presents a pair of quasi-resonant sites for the particle in orange; in the case of interacting particles, region [START_REF]Anderson Localization and its Ramifications[END_REF] shows two quasi-resonant sites for the orange particle, thanks to the partial screening of disorder offered by the interaction with the red particle.

(1) (2) (a) (b) (1) J z = 0 J/2 J z δ -J z J/2 δ (2) J z ≠ 0
disorder potential via the P distribution. In the case of the QP potential

P (x) = [1 -(x/ ∆) 2 ] -1/2 π ∆ (7)
with ∆ = ∆ sin(πκ) [START_REF] Guarrera | [END_REF]; while for the FR potential P (x) is the normalized triangular distribution defined on the [-2∆, 2∆] interval. The cusp singularities in λ(t) at times t m J = (2m + 1)π, m = 0, 1, 2, . . ., descend from the fact that the integrand function f (δ; t m ) = log[1-p(δ, J, t m )], seen as a function of δ, develops a logarithmic singularity at δ = 0, as shown in figure 3, while it is fully regular at any other time. The singular peak centered at δ = 0 has a support shrinking with t m as t -1/2 m -as seen in figure 3 when plotting the function f (δ; t m ) vs δ √ t m , which leads to a collapse of the peak widths at different times (when m 1). The integral of the f function outside the peak contributes to the regular part of the Loschmidt echo, while the integral of the peak dictates fundamentally the height of the cusps above the regular background (estimated as the long-time average λ), namely the quantity λ P (t m ) = λ(t m ) -λ. The decay of the height of these cusps as t -1/2 m will be verified numerically in section 6. Figure 4(a-b) shows that, for the case of the QP potential, Eq. ( 6) is able to predict with high accuracy the ED results deep in the MBL phase without any adjustable parameter. In particular the cusp singularities of the ED results are easily explained as descending from the divergent singularity of the Loschmidt echo for a fully resonant 2LS with Ω = Ω = J, reaching a state orthogonal to the initial one after odd multiples of half a Rabi oscillation t m = (2m + 1)π/Ω. These divergences are smoothened into cusp singularities due to the fact that such resonant 2LSs are a set of zero measure in the disorder statistics. This result has important consequences. Indeed the nearly resonant 2LSs captured by the model are clearly an ensemble of approximate l-bits with Hamiltonian

H ≈ p K p τ p , (8)
-1 -0.5 0 0.5 1 where τ p = δ p /K p σ z p -J/K p σ x p is a Pauli matrix expressed as a rotation of the Pauli operators σ z p = S z i+1 -S z i (when projected onto the subspace with S z i + S z i+1 = 0) and σ x = S + i S - i+1 + h.c., built from the original spin operators for the pair p = (i, i + 1); and K p = δ 2 p + J 2 is the l-bit splitting. Hence the Loschmidt-echo singularities are a striking manifestation of the existence of such (nearly free) l-bits, to be found in the short-time dynamics of the system.

δ/J 0 2 4 6 8 f (δ; t) (a) -4 -2 0 2 4 √ Jtδ/J (b) Jt = π Jt = 5π Jt = 15π Jt = 25π

Relationship to dynamical quantum phase transitions

It is worthwhile to mention at this point that the existence of singularities in the quench dynamics of the Loschmidt echo is currently the subject of several theoretical and experimental investigations, as they represent the main signature of so-called dynamical quantum phase transitions, studied both in non-random systems [START_REF] Heyl | [END_REF][25][26]30,31] as well as in disordered quantum systems [27,28,58]. Nonetheless our observation of Loschmidtecho singularities is fully explained by a model of individual 2LS, without the need of any many-body effect. Therefore we shall refrain from associating them to any form of time-dependent transition.

From two-level systems to three-level ones

Figure 4(d-f) shows that, in the case of the FR potential, the 2LS model of ( 6) still predicts the correct frequency of the Loschmidt echo singularities, but not the correct height; and it also misses a global offset. This is not surprising, as in the case of the FR potential the assumption of anti-correlation between the energy offset of contiguous pairs is no longer valid, namely the potential can host "rare" regions in which contiguous pairs of sites -(i-1, i) and (i, i+1) -are nearly resonant at the same time. To take those regions into account (at least partially) one can easily promote the 2LS model to a model of 3-site systems (amounting to effective three-level systems -3LS), and approximate the evolved state as that of a collection of independent 3LS. The Hamiltonian of a threesite system (i, i + 1, i + 2) containing two interacting fermions in an initial |101 state is explicitly given by

H 101 = - J 2 c † i c i+1 + c † i+1 c i+2 + h.c. + δ i n i+1 + (δ i+1 + δ i )n i+2 + J z (n i n i+1 + n i+1 n i+2 ) + const. ( 9
)
Here all single-site energies are referred to the energy of site i, and δ i = h i+1h i . The above Hamiltonian assumes that the sites i -1 and i + 3 remain empty during the time evolution. The Hilbert space of the 3-site system is restricted to the three states |101 , |110 and |011 , making of it a three-level system (3LS), with a generic time-dependent wavefunction

|ψ 101 (t) = α(t) |011 + β(t) |101 + γ(t) |110 . (10)
Its explicit form can be easily calculated numerically for any specific choice of the energy differences δ i .

A similar calculation can be done for a three-site system hosting a single particle, and starting from the |010 configuration, with Hamiltonian

H 010 = - J 2 c † i c i+1 + c † i+1 c i+2 + h.c. + J z n i + δ i n i+1 + (δ i+1 + δ i + J z )n i+2 , (11)
which assumes that the sites i -1 and i + 3 host two pinned particles. The Hilbert space |100 , |010 , |001 defines a 3LS, whose instantaneous state takes the generic form

|ψ 010 (t) = α(t) |100 + β(t) |010 + γ(t) |001 . (12)
For the two types of clusters the Loschmidt echo can be readily evaluated as

λ 101 (t; δ i , δ i+1 ) = -log |β(t)| 2 and λ 010 (t; δ i , δ i+1 ) = -log | β(t)| 2 .
We can then model a chain in a QP or FR potential as an ensemble of independent 3LSs by generating sequences of energy offsets δ i , δ i+1 between adjacent site pairs according to the distribution P (δ i , δ i+1 ). The 3LS prediction for the Loschmidt echo of the ensemble is

λ 3LS = 1 2 dδ 1 dδ 2 P (δ 1 , δ 2 ) [λ 101 (t; δ 1 , δ 2 ) + λ 010 (t; δ 1 , δ 2 )] . (13)
In practice, the above integral can be sampled numerically by simply averaging over a large number of different realizations of the potential on 3-site systems, such as those offered by a very long chain, namely

λ 3LS (t) ≈ 1 L L i=1 λ α i (t; δ i , δ i+1 ), (14)
where α i = 101 if i is odd and 010 if i is even, and L 1. Eq. ( 14) for the Loschmidt echo has the apparent drawback of triple-counting each site. Nonetheless, similarly to what was argued for the 2LS case, it is fair to assume (and it can be numerically tested) that, out of the three clusters containing each site, only one at most will contribute significantly to the Loschmidt echo. As a consequence the triple counting has only a mild effect on the final result. One could avoid triple counting by thoughtfully decomposing a chain into non-overlapping clusters of up to 3 sites, in such a way as to maximize the Loschmidt echo; yet this procedure introduces significant complications which are not justified a posteriori, given the quality of the results offered already by the naive ensemble average. As shown by figure 4(c-d), the improvement offered by the 3LS model for the FR potential is substantial; these results can further be improved by moving to 4-site clusters etc., albeit at an exponential cost.

Imbalance dynamics

The 2LS model prediction for the imbalance is very similar to that of the Loschmidt echo, as the imbalance is simply related to the persistence probability of the initial state (|10 or |01 ) on the 2-site cluster -given that the orthogonal state contributes zero to the imbalance. Therefore the 2LS expression for the imbalance simply reads

I 2LS (t) = dδ P (δ) [1 -p(δ, t)] . (15)
The times t m giving cusp singularities in the Loschmidt echo correspond to dips in the imbalance, and these dips come from local dips in the g(δ; t) = 1p(δ, t) function centered around δ = 0 and touching zero for t = t m . The width of these dips is also shrinking in time as t

-1/2 m
. Therefore one expects the depth of the minima in the fluctuations of the imbalance to decay to the long-time average as t -1/2 m as well -this prediction will be verified in section 6.

Extending the 2LS model to 3LS, the imbalance can be calculated as

I 3LS (t) ≈ 1 3L i -|α i | 2 + 3|β i | 2 -|γ i | 2 (16)
with α i , β i , γ i = α(t), β(t), γ(t) or α(t), β(t), γ(t) depending on whether i is odd or even.

Figure 5 shows the comparison between the ED results for the imbalance dynamics of interacting fermions immersed in a QP and fully random potentials of variable strength, compared with the predictions of the 2LS and the 3LS models. For sufficiently strong disorder (∆ 6J) the 2LS predictions are already rather accurate in the case of the QP potential, and the 3LS model offers further improvement. On the other hand in the case of the FR potential the 3LS model offers a more substantial improvement, fixing an overall offset (for sufficiently strong disorder) which is seen in the 2LS predictions.

Dephasing in the Loschmidt echo and imbalance oscillations: evidence of l-bits interactions

A significant feature of the Loschmidt echo singularities is their slow decay in timewhich is remarkable given that they result from the Rabi oscillations of a collection of 2LS with a distribution of frequencies that can be a priori expected to lead to fast dephasing. The reason behind the slow decay is also captured by the 2LS model, Eq. ( 6) -namely by the fact that the integral expressing the Loschmidt echo takes contributions from a small window of detunings δ around zero, the smaller the longer the time, as mentioned in section 4.1. When looking at the singularity times t = t m , a direct inspection of the function log(1p(δ , Ω, t m )) seen as a function of δ shows that it has a large peak centered on δ = 0 with a width depending on time as t -1/2 m . The singularity in the average Loschmidt echo comes from the integral of this peak, while the rest of the integral contributes essentially to the regular part of the Loschmidt echo; hence it is immediate to predict that the height of the cusp singularity should decay as the peak width, (namely as t -1/2 m ). Figure 6(a) shows the time evolution of singularity peaks in the Loschmidt echo for free as well as interacting fermions in the QP potential, compared to the prediction of the 2LS model (for the interacting case): we observe that the t -1/2 decay is indeed confirmed by the ED data for free fermions, as well as by the ED data for interacting fermions at sufficiently short times (tJ t * ≈ 100 for ∆ = 8J). On the other hand, at longer times the interacting data are found to display a strong deviation from the 2LS model prediction, exhibiting a much faster decay. This crossover to an interactioninduced dephasing (IID) regime clearly shows the limits of the 2LS model as a model Comparison between the imbalance I(t) for a L = 22 chain and the predictions of the 2LS and 3LS models given by ( 15), ( 16): (a-c) QP potential; (d-f) FR potential. (c-d) Decay of the depth of the imbalance minima I(t m ) with respect to the average value I. ED data in both absence (J z = 0) and presence (J z = 1) of interactions are compared with 2LS and 3LS predictions. The data are obtained for ∆ = 8J; the 2LS and 3LS predictions are for J z = J. The grey-shaded area marks the interaction-induced dephasing (IID) regime exhibited by the exact data for J z = J.

of free l-bits, and it marks a fundamental difference between AL and MBL in the QP system. Indeed the faster decay of the Loschmidt echo must be related to the effect of l-bit interactions, which are a defining feature of MBL, and which add terms of the kind pq U pq τ p τ q + pql V pql τ p τ q τ l + ...

to the effective l-bit Hamiltonian (8). Such terms are responsible for the persistent growth of entanglement entropy in the system [7] as the logarithm of time, and indeed the onset of the log t growth of entanglement occurs at a time compatible with t * (see section 7 -figure 7). A similar crossover from a slow power-law decay of the Loschmidtecho peak height to a faster decay, dictated by the presence of interactions, is also exhibited by the comparison between the ED data for interacting fermions in the FR potential with the same data for non-interacting fermions and for the 3LS model -as shown in figure 6(b).

Remarkably, the same crossover between the dynamics of effectively independent l-bits to a regime of interacting ones is observed in the decay of imbalance oscillations. Figure 6(c) shows the evolution of the depth of the minima of the imbalance at times t = t m , taken with respect to the long-time average, namely the quantity I M (t m ) = Ī -I(t m ). We observe that the predictions of the 2LS system for the fermionic chain immersed in the QP potential shows a clear, slow power-law decay at long times, compatible with t -1/2 , which is indeed reproduced in the case of noninteracting fermions. In the case of interacting fermions, on the other hand, a crossover is observed at long times (t t * ≈ 100) to a faster decay, marking the IID regime. A similar picture is offered by the case of the FR potential. There the ED results are compared with the predictions from the 3LS model; the latter model predicts correctly the decay of the minima depth in the non-interacting case at all times, while the exact results for the interacting system show a clear crossover towards a faster decay for times t t * ≈ 50. For both kinds of disorder, the crossover time t * is compatible both with what is observed in the decay of the Loschmidt echo as well as with the evolution of the entanglement entropy (see again figure 7 in the next section). Therefore we conclude that the crossover to the IID regime is a robust feature of MBL dynamics, clearly exposing the interactions among l-bits.

Entanglement dynamics

Entanglement entropy from the 2LS and 3LS model

The 2LS and 3LS models allow for a simple calculation of the entanglement entropy of a A/B bipartition of the system into two adjacent chains, defined as the von Neumann entropy of the reduced density matrix

S A (t) = -Tr [ρ A (t) log ρ A (t)] , (18)
where ρ A (t) = Tr B |ψ(t) ψ(t)| is the partial trace (over the degrees of freedom in B) of the instantaneous pure-state density matrix associated with the evolved state |ψ(t) . The entanglement associated with such a bipartition simply comes from the entanglement inside the 2-site or 3-site cluster which contains the cut defining the bipartition. In the case of the 2LS model the disorder-averaged entanglement entropy of a bipartition is simply predicted as the entropy of the reduced state of one site in the 2-site cluster, namely

S A (t) = dδ P (δ) h[p(δ; t)], (19)
and

h[x] = -x log x -(1 -x) log(1 -x).
Notice that, unlike for the formulas of the Loschmidt echo and of the imbalance, no-double counting is implied in the above formula, since the entanglement is referred to a cut of the chain, and there is one unique cut per 2-site cluster.

On the other hand a 3-site cluster can be cut in two different ways, that we will indicate as • | • • and • • | • in the following (where • stands for a site and | stands for the cut). The reduced density matrices for the two cuts are readily obtained from the cluster wavefunctions described in section 4.3; e.g. for a 101 cluster the reduced density matrix associated with the

• | • • cut reads ρ •|•• 101 (t) = |α(t)| 2 0 0 |β(t)| 2 + |γ(t)| 2 , (20)
with associated entanglement entropy

S •|•• the one associated with the • • | • cut reads ρ ••|• 101 (t) = |γ(t)| 2 0 0 |α(t)| 2 + |β(t)| 2 . ( 22
)
The density matrices ρ

•|•• 010 and ρ ••|•
010 and related entropies associated with a 010 cluster can be calculated similarly. The disorder-averaged entanglement entropy of the whole system within the 3LS model is then given by

S(t) = 1 2L i S i (t), (23)
where

S i = S •|•• 101 (t) + S ••|• 101 (t) if i odd, S •|•• 010 (t) + S ••|• 010 (t) if i even. ( 24
)
The factor 1/2 in ( 23) comes from the double counting of each cut (which is contained both in a 101 cluster as well as in a 010 cluster).

Figure 7 shows a comparison between the entanglement entropy of interacting fermions in a QP potential and the 2LS prediction. We observe that at moderate disorder in the MBL phase (∆ = 8J) the 2LS and 3LS models only capture the initial rise of the entanglement entropy and (partly) the first maximum; in particular the very existence of a maximum is explained by the models as the result of nearly resonant small clusters returning close to the initially factorized state -albeit at different times due to the inhomogeneously broadened local frequencies, which explains why the entanglement entropy does not come back to (nearly) zero. The 2LS and 3LS models on the other hand completely miss the long-time logarithmic growth of the entanglement entropy -something which is fully expected, given that such a growth is the consequence of interactions between l-bits, not included in the 2LS and 3LS models by construction. On the other hand, at stronger disorder (∆ = 15J) the interactions between l-bits are parametrically suppressed, and the 2LS and 3LS description of entanglement becomes accurate up to very long times.

Entanglement entropy vs. number entropy

The 2LS and 3LS models picture the entanglement between two adjacent subsystems as arising uniquely from the coherent motion of particles within the restricted size of the clusters they describe. When starting from a factorized state, this picture is certainly valid at short times. At long times it remains valid only if particles remain localized within the size of the clusters (namely if the localization length is smaller than the cluster size), and if this is a sufficient condition for entanglement not to spread any further. The latter aspect is true in the case of non-interacting fermions, for which the only mechanism behind entanglement of different spatial partitions is particle motion between them. On the other hand, in the case of interacting fermions in the MBL regime, entanglement keeps growing due to the interactions between l-bits, and distant degrees of freedom can become entangled even without any net particle exchange. In this context it is useful to decompose the entanglement entropy of a subsystem A into a number entropy contribution, and a remainder part (called the configurational entropy), S A = S A,N + S A,c [22,59]. The number entropy is given by

S A,N = - N A p N A log p N A , (25)
where p N A is the probability of having N A particles in subsystem A. Eq. ( 25) accounts for the particle number uncertainty appearing in subsystem A because of the coherent exchange of particles with its complement B. On the other hand the configurational entropy accounts for correlations establishing between the particle arrangements in A and B once the partitioning of the particles between A and B has been fixed. The 2LS and 3LS models, completely lacking any form of correlations among the clusters, can only capture the number entropy contribution in systems with a localization length smaller than the cluster size. Nonetheless, this limited picture still offers a faithful description of entanglement in the MBL regime for short times (the longer the stronger disorder is), while it can describe entanglement at all times for strongly localized noninteracting particles. Thus, as suggested above, a more appropriate comparison with the entanglement entropies of the 2LS and 3LS models would involve the number entropy from the ED data -shown in figure 8. For non-interacting fermions in a QP potential of strength ∆ = 8J, S A,N is found to nearly coincide with the full entanglement entropy, and to be very well described by the 2LS prediction -see figure 8(a). When adding the interactions, the agreement between the number entropy and the 2LS entropy deteriorates mostly at long times, seemingly due to the ∼ log log t growth of the number entropy observed in the MBL phase [59]. Similar considerations can be made in the case of the FR potential, i.e. the 3LS models describe well the entropies in the non- interacting case, and they miss the slow long-time growth of the number entropy in the interacting case.

Conclusions

In this work we have shown that sharp cusp-like singularities in the Loschmidt echo are a generic feature of the localized dynamics of an extended quantum system initialized in a factorized state. These features can be fully explained by the dynamics of a simple model, describing an ensemble of effective independent two-level (or even three-level) systems, offering an explicit approximation to the conserved l-bits in the AL and MBL regime. Such a model predicts very accurately the Loschmidt echo singularities for strongly disordered systems as well as their decay, along with the imbalance oscillations. A faster decay in the Loschmidt echo and imbalance dynamics compared to that predicted by the model is a direct manifestation of the dephasing effect of interactions between the l-bits, and it intervenes at a time consistent with the onset of the logarithmic growth of entanglement entropy: hence it represents a defining feature of many-body localization (MBL) with respect to Anderson localization (AL). Based on our results, we can conclude that experimental evidence of l-bit dynamics and of their interactions is readily accessible to state-of-the-art quantum simulators which have direct access to the Loschmidt echo and imbalance dynamics, such as e.g. trapped ions [26,34], cold-atom simulators [22,35,37] or superconducting circuits [START_REF] Chiaro | Direct measurement of non-local interactions in the many-body localized phase[END_REF]60].

Technology Platform and HPC Midlands+ (Athena) for computing time and support (EPSRC on grant EP/P020232/1). Part of the exact diagonalization simulations were performed using routines contained in the QuSpin [61,62] and quimb [63] open-source libraries.

Figure 1 .

 1 Figure 1. Loschmidt echo and imbalance dynamics for an L = 22 chain with QP potential (a-c) and FR potential (b-d), for various disorder strengths ∆ = 2, . . . , 8 as indicated by the colors.

  δ

Figure 2 .

 2 Figure 2. (a) Example of a L = 22 chain in a QP potential (lines) in the initial CDW state |1010... . Particles are denoted as coloured balls. (b) Zoom on two quasiresonant regions (shaded areas): in the case of non-interacting particles (J z = 0) the region (1) presents a pair of quasi-resonant sites for the particle in orange; in the case of interacting particles, region[START_REF]Anderson Localization and its Ramifications[END_REF] shows two quasi-resonant sites for the orange particle, thanks to the partial screening of disorder offered by the interaction with the red particle.

Figure 3 .Figure 4 .

 34 Figure 3. Function f (δ; t) vs. δ at different singularity times t m J = (2m + 1)π; as shown in the right panel, for large t m the width of the central peak becomes timeindependent when f is plotted as a function of √ t δ.

Figure 5 .

 5 Figure 5. Comparison between the imbalance I(t) for a L = 22 chain and the predictions of the 2LS and 3LS models given by (15), (16): (a-c) QP potential; (d-f) FR potential.

Figure 6 .

 6 Figure 6. (a-b) Decay of peak heights of the Loschmidt echo, λ(t m ) -λ ( λ stands for the time-averaged Loschmidt echo).(c-d) Decay of the depth of the imbalance minima I(t m ) with respect to the average value I. ED data in both absence (J z = 0) and presence (J z = 1) of interactions are compared with 2LS and 3LS predictions. The data are obtained for ∆ = 8J; the 2LS and 3LS predictions are for J z = J. The grey-shaded area marks the interaction-induced dephasing (IID) regime exhibited by the exact data for J z = J.

10 °1 10 0 10 1 10 2 10 3 Jt10 °1 10 0 10 1 10 2 10 3 JtFigure 7 .

 337 Figure 7. Half-chain entanglement entropy of interacting fermions in a QP potential (a-b) and FR potential (c-d) for a chain of size L = 16, compared with the prediction for the 2LS/3LS models for two different disorder strengths (∆/J = 8 and 15).

10 -1 10 0 10 1 10 2 10 3 Jt10 -1 10 0 10 1 10 2 10 3 JtQP, J z = 1 10 -1 10 0 10 1 10 2 10 3 JtFR, J z = 0 10 -1 10 0 10 1 10 2 10 3 Jt

 331303

Figure 8 .

 8 Figure 8. Half-chain entanglement entropy and number entropy of free and interacting fermions on a L = 16 chain, compared with the 2LS and 3LS predictions: (a) noninteracting fermions in a QP potential; (b) interacting fermions (J z = 1) in a QP potential; (c) non-interacting fermions in a FR potential; (d) interacting fermions (J z = 1) in a FR potential. For all the panels the disorder strength is ∆ = 8J.

Figure C1 .

 C1 Figure C1. Averaged λ(t) (black line) at different disorder strengths ∆ = 1, . . . , 10 for a chain of L = 22 sites in a QP, plotted along with all the realizations used for the averaging procedure (grey lines). The dotted line represents a typical individual realization exhibiting singular behavior.

Figure C2 .

 C2 Figure C2. Averaged λ(t) (black line) at different disorder strengths ∆ = 1, . . . , 10 for a chain of L = 22 sites in a FR potential, plotted along with all the realizations used for the averaging procedure (grey lines). The dotted line represents a typical individual realization exhibiting singular behavior.

Acknowledgments L.B. gratefully acknowledges hospitality and financial support from the Laboratoire de Physique of the ENS Lyon. R.A.R. and L.B. acknowledge funding from the CY Initiative of Excellence (grant"Investissements d'Avenir" ANR-16-IDEX-0008) where this work developed during R.A.R.'s stay at the CY Advanced Studies. L.B. thanks the EPSRC for DTP funding. We thank Warwick's Scientific Computing Research

Appendix A. Two-site cluster as a two-level system and its Rabi oscillations Let us isolate a two-site system (i, i + 1) hosting one particle in the fermionic chain, with Hamiltonian

where we assume that the site i + 2 is occupied by a (pinned) particle, while size i -1 is empty (or occupied by a pinned hole). Introducing the spin operators

the Hamiltonian becomes simply

namely, a two-level system (2LS) with splitting δ and Rabi frequency J. If the system starts from the |10 state, the return probability is given by the well-known formula for the probability of persistence in the initial state during Rabi oscillations [START_REF] Scully | Quantum Optics[END_REF], namely 1p(δ; t), where p(δ; t) = 1 1 + (δ/J) 2 sin 2 1 + (δ/J) 2 2 tJ .

(A.4)

Appendix B. Spatial correlations in the quasi-periodic vs. fully random potential

A fundamental assumption of the 2LS model described in the main text is that quasiresonant two-site systems are spatially isolated in a (quasi)-disordered chain -namely, if a pair of sites (i, i + 1) is quasi-resonant for the motion of a particle, the two adjacent pairs of sites (i -1, i) and (i + 1, i + 2) are not resonant. Defining δ 1 = h i+1h i as the energy difference of the two sites in question, and δ 2 = h i+2h i+1 as that of the following pairs of sites, in the case of non-interacting fermions, the above condition requires that the two energy differences do not vanish simultaneously. Such a form of correlation is indeed observed in the case of the quasiperiodic (QP) potential: figure B1(a) shows the joint probability P (δ 1 , δ 2 ) for two adjacent energy differences, displaying a dip for δ 1 = δ 2 = 0 -an aspect which prevents two successive pairs of sites from being resonant simultaneously. In the case of interacting fermions, on the other hand, the above condition requires that if, e.g. J z ± δ 1 ≈ 0, then J z ∓ δ 2 is non-zero, or vice versa -this prevents a state of the type |101 on the sites (i, i + 1, i + 2) from being simultaneously (quasi-)resonant with |110 and |011 , or, similarly, the state |010 from being quasi-resonant with |100 and |001 . This is indeed guaranteed by the fact that P (δ, -δ) is nearly vanishing for any finite δ, except for δ ≈ 1.25∆ -but the latter situation does not lead to consecutive resonances when ∆ > 1.25J, which is always the case in our study.

On the other hand the uniform potential has no correlations between two consecutive energy differences, and the P (δ 1 , δ 2 ) distribution is the product of two triangular distributions for δ 1 and δ 2 -shown in figure B1(b). This implies that a fundamental assumption behind the 2LS model description is not guaranteed to be satisfied -while it is more likely to have two adjacent pairs of sites with different energy offsets than with similar ones, one cannot exclude the existence of "rare" regions with consecutive nearly resonant pairs. This requires to improve the 2LS model to a three-site (three-level) one, as detailed in section 4.3.

Appendix C. Loschmidt-echo dynamics for different disorder realizations

Figs. C1 andC2 show the disorder average of the Loschmidt echo for a chain of L = 22 sites, along with all the disorder realizations (> 10 3 ) contributing the average, for various strengths (∆/J = 1, 2, ..., 10) of the QP and FR potential, respectively. We observe that sharp cusp singularities are exhibited by a signification portion of the disorder realizations, and that for sufficiently strong disorder these realizations are a finite fraction of the disorder statistics (in the asymptotic limit), so that cusp singularities persist in the disorder-averaged results as well. These plots also suggest the fact that cusp singularities can be observed with a limited disorder statistics, under realistic experimental conditions.


Josef Oswald
email: josef.oswald@unileoben.ac.at

Rudolf A Römer
email: r.roemer@warwick.ac.uk

Microscopic details of stripes and bubbles in the quantum Hall regime

Keywords: numbers: 73.43.-f, 73.43.Nq, 73.23.-b

We use a fully self-consistent laterally resolved Hartree-Fock approximation for numerically addressing the electron configurations at higher Landau levels in the quantum Hall regime for nearmacroscopic sample sizes. At low disorder we find, spatially-resolved, stripe-and bubble-like charge density modulations and show how these emerge depending on the filling factor. The microscopic details of these boundary regions determine the geometrical boundary conditions for aligning the charge density modulation either as stripes or bubbles. Transport is modelled using a non-equilibrium network model giving a pronounced anisotropy in direction of the injected current in the stripe regime close to half filling. We obtain a stripe period of 2.9 cyclotron radii. Our results provide an intuitive understanding of its consequences in strong magnetic fields and indicate the dominance of many particle physics in the integer quantum Hall regime when studied at legnth scales .

In May 2019, the quantum Hall effect [1,2] was formally included among the select group of high-precision experiments to form the basis of a new SI system of units based on the "Planck constant h, the elementary charge e, the Boltzmann constant k, and the Avogadro constant N A " [3]. This had long been awaited and certainly represents a great achievement and fitting fulfillment of the vision for "natürliche Masseinheiten" [4] proposed by Max Planck [5]. In his essay to celebrate this achievement [3], von Klitzing also points out "that a microscopic picture of the quantum Hall effect for real devices with electrical contacts and finite current flow is still missing." Prominent examples of such microscopic details are so-called "bubble" and "stripe" phases [START_REF] Fogler | High Magnetic Fields[END_REF]. They have been identified, e.g., by transport experiments in higher Landau levels (LLs) of ultra-high mobility samples [START_REF] Lilly | [END_REF][8][9] and are characterized by strong transport anisotropies (stripes) or reentrance effects (bubbles). It is believed that the phases correspond to density modulations with characteristic geometric non-uniformities due to the interplay of Coulomb interaction and the wave functions in higher Landau levels.

Early work in modelling density modulations in the quantum Hall regime, starting from the celebrated Chklovskii, Shklovskii and Glazman picture [10], assumed uni-directional charge-density waves (CDWs) [11,12] while mean field treatments established the possibility of anisotropic phases in a Fermi liquid [13,14]. However, spatially resolved information does not yet exist of these phases. Experimentally this is due to the intrinsic challenge of using local scanning probes in low temperatures for such remotely doped systems [15][16][17]. Nev-ertheless, much indirect experimental evidence for the existence of bubble and stripe phases has now been accumulated [17][18][19][20][21][22][23][24][25][26][27][28]. Theoretical modelling has likewise concentrated on transport signatures of these phases [29][30][31][32] while spatially resolved models of bubbles and stripes are only available in clean systems [33]. The period of the stripe patterns has been predicted to follow d ∝ R c with R c = l B √ 2n + 1 the cyclotron radius in LL n at the Fermi energy [11,34] and l B = /eB the magnetic length. Experimentally, 1.5R c [17] and 3.6R c [18] have been reported while ∼ 2.7R c is predicted theoretically [34].

In the present work, we show how stripes and bubbles emerge at weak disorder as self-consistent solutions of the Hartree-Fock (HF) equations, i.e. in the experimentally relevant regime and without any ad hoc assumptions beyond a smooth disorder. We provide the full spatial resolution of both phases from the length scale of l B to near macroscopic sample sizes. This high resolution allows quantitative comparison with current experimental efforts [17,18,24,26,27]. A central insight provided by our work is the importance of many-body aspects. It should be clear that a fully self-consistent HF approach in a disordered environment goes well beyond earlier Thomas-Fermi-based (non-)linear screening models. The inclusion of a converged exchange interaction term essentially alters the physics. This not only changes the spatial distribution of stripes and bubbles, but rather is the main reason of their emergence: neither pure Hartree nor a non-interacting model leads to emerging stripes/bubbles unless coupled with additional assumptions. The key-mechanism is a Hund's rule behaviour for the occupation of the spin-split LLs. The resulting g-factor enhancement is then a local quantity depending on the local filling factor ν( r) = 2πl 2 B ρ( r) [35,36] with ρ( r) the local carrier density. This exchange- enhanced g-factor is a concept that allows to discuss the exchange interaction within the single electron picture [37,38]. By doing so, we find that the local variation of the enhancement of the Zeeman energy due to ν( r) has to be considered in addition to the laterally varying Hartree potential, leading to a modified effective potential for the electrons that also strongly modifies the screening behaviour. In addition to the largely repulsive Hartree part of the self-consistent Thomas-Fermi screening, the ν( r) dependence of enhanced Zeeman energy leads to a positive feedback loop in the self-consistent carrier redistribution and produces an instability of the electron density ρ, which may lead to jumps either to a locally full or locally empty LL [35,36], resulting in a clustering of the filling factor that is triggered by the disorder or edge potential. The boundaries of those clusters finally create narrow channels that align mainly along the edge or random potential fluctuations. In case of a very clean high mobility electron system such a trigger effect for the cluster formation by the random potential is missing and the electron system has to find such a cluster structure by self organization, resulting in the formation of stripes or bubbles.

The numerical simulations are performed as described in Refs. [35,36,39] via a variational minimization of the self-consistent Hartree-Fock-Roothaan equation [40][41][42][START_REF] Macdonald | [END_REF]. To simulate ultra-high mobility samples, the random potential strength is kept low. The filling factor and the lateral system size are made as large as possible with respect to the available computing power. Configurations up to 1 µm 2 are achievable as shown in Fig. 1 with spatial resolution of ∼ 4.4 nm well below l B for a magnetic field B varying from, e.g., 1 (l B ∼ 26 nm) to 6.5 Tesla (l B ∼ 10 nm). The random potential is generated by Gaussian impurity potentials of radius 40 nm, the number of impurities is N = 2000, their random placement results in a fluctuating potential of V max = 0.43 mV and V min = -0.50 mV [START_REF]ν( r) for pure Hartree interaction and without any interaction, a plot[END_REF]. At a total filling factor of, e.g., ν = 6.54 (= ν ↓ + ν ↑ = 3.54 + 3.0) and B = 1.5 Tesla, this corresponds to more than 2000 electrons. In order to generate transport data, we employ a non-equilibrium network model (NNM) introduced previously [45]. A very large number of step by step calculations are required [START_REF] Sohrmann | Quantum and Semi-classical Percolation and Breakdown in Disordered Solids[END_REF] in the NNM. Hence, for keeping within the available computing time, the transport simulations have been performed for smaller sample size such as 500 × 500 nm 2 .

In Fig. 1 we depict the variation in ν ↓ ( r) for three different densities at fixed magnetic field. Fig. 1(a) shows the situation far from any half-odd average ν ↓ . The required area of the filled (spin-down) clusters is only a minor part of the total area. As seen in the figure this can be achieved by a nearly evenly-spaced distribution of "bubble"-like shaped clusters. Since we are considering already the 4th partly-filled spin-split LL, the boundaries of the bubbles consist of three sub-stripes because of the 3 nodes of the Landau basis function for the 4th LL [16].

A boundary with such an internal structure takes up a substantial area as well and even tends to dominate the region of a single bubble as a whole. In other words, the area that one bubble needs is dominated by the width of the boundaries and not by the region of an assumed idealized full LL. When increasing the filling factor towards half filling, as shown in Fig. 1(b), it becomes clear that such a bubble-like geometry of the clusters is impossible to achieve with the given total ν, because round bubbles would leave too much unfilled space between the bubbles even if touching each other. The only way to remove the unfilled space between bubbles is a change of the geometry so that the boundaries of different clusters arrange almost in parallel, which means a transfer to a "stripe"like geometry. In order to demonstrate that the width of the boundaries of the stripes and bubbles depend on the LL index, Fig. 1(c) shows the half-filled LL at filling factor ν ≈ 4.5 instead of ν = 6.5 and one can see that there is one sub-stripe less at the boundary, allowing a more dense arrangement of stripes or bubbles than in the higher LLs. These results suggest that the higher the filling ν the higher the tendency for creating the stripe pattern near a half-filled LL. Overall, the figures show that our HF calculations can provide a fascinating insight in the spatial behaviour of bubble and stripe configurations.

Fig. 2 shows the transport data obtained by applying the NNM at different carrier densities. Since the stripe alignment tends to be more horizontal, i.e. along the x- direction as shown in Fig. 1, the longitudinal resistance R xx appears higher for vertical sample current and lower for horizontal current flow [START_REF]ν( r) for pure Hartree interaction and without any interaction, a plot[END_REF]. This is consistent with experimental observations resulting in large R xx peaks for vertical current, while for a horizontal current the R xx peaks are hardly visible. When the difference between horizontal and vertical R xx is no longer prominent, we find that we have reached a bubble phase.

The characteristics of the microscopic structure of the stripes consist on one hand in the periodicity of the stripe pattern and on the other hand in the microscopic details of the boundaries of the stripes. Together with the geometric shape of the stripes (and, indeed, the bubles), these determine the relation between the areas of full and empty LLs. Focusing first on the periodicity, Fig. 3 presents the 2D-Fourier transformation of the lateral carrier density at various magnetic fields from B = 1 Tesla to B = 6.5 Tesla at a fixed filling factor of ν = 4.5. While Fig. 3 (a) shows only a selection of a few spectra, Fig. 3 (b) displays the trend of all evaluated spectra. As can be seen in x-direction there is a smooth and somewhat smeared out distribution that starts from zero, indicating that there is no clear periodicity in x-direction. In contrast, in y-direction, there is a well pronounced maximum which matches the reciprocal stripe period in y-direction as also seen in Figs. retical results on the period d of the stripes have led to the expectation d = αR c = α (2n + 1) B/e [11,34]. Experimentally, α = 1.5 [17] and 3.6 [18] have been estimated. For n = 3 and B = 1.5 Tesla, these give 88 nm and 210 nm, respectively. Obviously, only the latter one is compatible with our result at 1.5 Tesla. Conversely, from Fig. 3 (b), we extract a value α = 2.9 ± 0.1. This agrees very well with previous straight-line CDW-based predictions of 2.7 [34] and 2.8 [12].

Friess et al. [17] investigated the Knight shift in nuclear magnetic resonance (NMR) spectra to demonstrate the co-existence of regions with different spin polariza- tion due to the periodic variation of the filling factor in stripe and bubble phases. Their method provides direct information about the area fractions, although there is no direct information about geometry and periodicity. In order to extract also microscopic details of the structural information, they use a semiclassical model [17,[START_REF] Tiemann | [END_REF] based on superpositions of the single electron densities obtained from the Landau basis functions. In our case, we can similarly model the NMR intensity as

I ν (f ) = G{f -[f 0,ν-1 -ν( r) • K max,ν ]
}) dr 2 with Gaussian G describing the absorption spectrum of individual nuclei [START_REF] Tiemann | [END_REF], f the NMR frequency, K max the maximal Knight shift for the fully spin-polarized LL at odd ν and f 0 the frequency of the non-shifted NMR line for the non spin-polarized situation at even filling factor ν -1.

Numerically, I ν (f ) is calculated by evaluating our interacting ν( r) at each r and summing over all points of a typically 229 × 229 grid. In Fig. 4, we show that I ν (f ) exhibits features that red-shift to lower frequencies when increasing ν, e.g., from the non-spin polarized ν = 4 to a fully spin polarized ν = 5 (see supplement [START_REF]ν( r) for pure Hartree interaction and without any interaction, a plot[END_REF] for ν = 2 → 3). At intermediate filling factors the spectral response splits, mainly into double/triple peak structures. This indicates the co-existence of regions with different filling factor as expected from a stripe or bubble like electron distribution (see supplement [START_REF]ν( r) for pure Hartree interaction and without any interaction, a plot[END_REF]). Around ν ∼ 4.4 one can recognise three peaks with changing weights in the superposition while varying the total filling factor. This is in good agreement with the experimental results [17], but naturally less so when compared to a single-particle modelling.

Intuitively, the driving force for the formation of bubbles and stripes can be understood as follows: the Hund's rule behaviour causes a g-factor enhancement that maintains the tendency to fill up just one spin-level as much as possible before starting to fill up the next spin-level. When entering a stripe/bubble from outside (e.g. ν = 4.0, not spin polarized) to inside (e.g. ν = 5.0, fully spin polarized) the g-factor enhancement also drastically increases. Inside the stripe/bubble the spin-levels are strongly pushed apart, corresponding to a significantly lower energy of the occupied lower spin state, while outside there is minimal g-factor enhancement at even filling factor ν = 4.0. In terms of an effective single particle picture the electrons in the stripes/bubbles encounter a potential well established by the local variation of the g-factor enhancement and which dominates over the repulsive Hartree interaction. From a qualitative point of view, the stripes can be understood as leaky non-coupled electron wave guides that consist of self assembled onedimensional potential wells and the leakage shows up as the boundary region as discussed before.

In conclusion, we have computed the microscopic picture of stripe and bubble phase for the IQH regime at high LLs and in weak disorder. Our results rely on spatially-resolved, self-consistent HF calculations of nearly macroscopic sizes of O(1 µm 2 ) and, for transport calculations, are coupled to device contacts with finite currents. The existence of microscopic stripes and bubbles at weak disorder is thus confirmed. Their spatial features show intriguing extended oscillations along the stripes and surrounding the bubbles. These are clearly due to the structure of the underlying Landau states. We find that the stripe period scales with B -1/2 as expected and agrees in detail with previous experimental measurements. Overall, together with results from HF calculations in the strong disorder regime [35,36,39], this shows that the IQH regime can now be described spatially resolved with high accuracy from microscopic to near-macroscopic length scales. Our results shine new light on the understanding of the microscopic picture of the IQH effect and demonstrate the permanent dominance of many particle physics for quantum Hall physics. The demonstrated Hund's rule behavior in context with the g-factor enhancement allows to incorporate the exchange interaction into an intuitive understanding of the major effects driving weakly and strongly disordered quantum Hall systems.

Supplemental Material

In Fig. S1 we show the variation in ν( r) for three different densities at fixed magnetic field B. These results are identical to Fig. 1 of the main text for ν ↓ ( r) but also show ν ↑ ( r) in the left column. The choice of colors is as in Refs. [1,2]. In Fig. S2(a) and (c) one can see that for pure Hartree interaction there is no stripe formation. The density modulation in ν( r) is much less than in Figs. 1 andS1 and roughly follows the random potential shown in Fig. S3. Furthermore, the charge density modulation in the spin-up and spin-down levels "repel" each other due to the Hartree interaction. In Fig. S2(b) and (c) we show the situation without interaction. Clearly, there is also no stripe formation. The ν( r) modulation rather closely follows the random potential of Fig. S3. The charge densities in the spin-up and spin-down levels do not influence each other and follow nearly identically the disorder potential because of any missing interaction.

Fig. S4 complements Fig. 2 by showing in addition the Hall resistance R xy .

In order to show that the presence of the remnants of the LL wave functions around each stripe is significant, we perform the calculations of I ν (f ) for three test patterns for ν( r). The results are given in Fig. S5. We find that only the variation given by ν( r) as calculated in HF can reproduce essential global features of the experimental NMR results presented in Ref. [3].

In Fig. S6 we show the behavior of I ν (f ) for stripes and bubble-like charge density waves. As in Fig. S5, the HF results for ν( r) lead to a reasonable qualitative agreement with the non-interacting model used in Ref. [3]. Nevertheless, the details around, e.g., ν = 2.5 are rather different, highlighting the importance of interactions. We note that normally stripes appear only starting with filling factor ν = 4.5. This is known also experimentally, but for experimental reasons the authors of Ref. 3 could not go to that filling factor. Instead, they used filling factor ν = 2.5 and forced, by using an in-plane component of the magnetic field, the electron system to form a stripe pattern. Clearly, there is no need for our simulations to also model this experimental "trick". In order to compare the effect of stripe patterns on the NMR Knight shift, we therefore use the stripe pattern in the "correct" range ν = 4-5. The result in Fig. S6 (b) has striking similarities but seems indeed a bit richer in features than the experimental curve for ν = 2-3. However, since the Knight shift spectrum looses its local information due to the spatial integration, we can also evaluate the range ν = 2-3 as shown in Fig. S6. Indeed the agreement with the experiments of Ref. 3 becomes even better in this filling factor range. We can still see in total 3 peaks, two of them clearly separated and a third one as a shoulder on the high frequency flank, just as shown for the experiments in Fig. 2b of Ref. 3. This makes the agreement almost perfect for the experiments as shown in Fig. S7 albeit not with the non-interacting modelling of Fig.

2c [3]. Filling factor 6-7 requires much more computing time and, although a most interesting question, there are currently no experiments available for comparison.

FIG. 1 .

 1 FIG. 1. Spatially resolved filling factor distribution ν ↓ ( r) for different total filling factors (a) ν = 6.20, (b) ν = 6.54 and (c) ν = 4.54 with B = 1.5 Tesla. The colours denote different ν ↓ ( r) values as indicated in the legends. The thin black lines are contours.

FIG. 2 .

 2 FIG. 2. Longitudinal resistance Rxx for easy (horizontal) and hard (vertical) direction (cp. Fig. 1) as a function of carrier density ρ. The sample size is 500 × 500 nm 2 , B = 1.5 Tesla, the impurity radius 40 nm, and the number of impurities N = 100. The random impurity placement results in a fluctuating potential from Vmin = -1.7 mV to Vmax = 1.7 mV.

FIG. 3 .

 3 FIG. 3. (a) Two-dimensional Fourier spectra at ν = 4.5 for different B fields. The different spectra are shifted to the right for better visibility. For a sample size of 1000 × 1000 nm 2 , the frequency in units of 0.001 nm -1 (left axis) equals the number of stripes that can be accommodated within the boundaries. The corresponding period length is shown on the right axis and the color bar indicates the Fourier intensity in arb. units. The horizontal red lines incidate selected values. (b) Variation of the reciprocal stripe period 1/d ( ) as function of √ B. The solid lines corresponds to d = α √ 2n + 1Rc for LL index n = 3 with α = 2.7 (blue line, theory) [34], 3.6 (red line, experiment) [18] and 2.9 (black line, our fit), respectively. The grey dotted grid lines highlight selected B and d values.

  FIG. 4. (a) Simulation results of Iν (f ) using our HF results for ν ↓ ( r) as input. The traces are shifted vertically for better visibility and varied in steps of ∆ν = 0.1. (b) Color plot of Iν (f ) for very fine resolution ∆ν = 0.01. The color legend gives the numerical values of Iν (f ).

  FIG. S1. Lateral carrier density distribution mapped on the filling factor scale ν( r) for different total filling factors (a) ν = 6.20, (b) ν = 6.54 and (c) ν = 4.54. The left and right columns contains results for ν ↑ and ν ↓ , respectively. The color shades represent the filling factor range, where green indicates the second LL for ν ↓,↑ = 1 → 2, red the third LL, light blue the fourth and yellow denoting LL 5. These colours are as in Refs. 1 and 2 and the value indicated for ρ, B denote the electron density in units of 10 11 cm -2 and 1 Tesla, respectively..

FIG. S4. Hall resistance Rxy for the longitudinal resistance Rxx shown in Fig. 2 with the same set of parameters.

This work received funding by the CY Initiative of Excellence (grant "Investissements d'Avenir" ANR-16-IDEX-0008) and developed during RAR's stay at the CY Advanced Studies, whose support is gratefully acknowledged. We thank Warwick's Scientific Computing Research Technology Platform for computing time and support. UK research data statement: Data accompanying this publication are available from the corresponding authors.