Dynamics and level statistics of interacting fermions in the Lowest Landau Level

We consider the unitary dynamics of interacting fermions in the lowest Landau level, on spherical and toroidal geometries. The dynamics are driven by the interaction Hamiltonian which, viewed in the basis of single-particle Landau orbitals, contains correlated pair hopping terms in addition to static repulsion. This setting and this type of Hamiltonian has a significant history in numerical studies of fractional quantum Hall (FQH) physics, but the many-body quantum dynamics generated by such correlated hopping has not been explored in detail. We focus on initial states containing all the fermions in one block of orbitals. We characterize in detail how the fermionic liquid spreads out starting from such a state. We identify and explain differences with regular (single-particle) hopping Hamiltonians. Such differences are seen, e.g. in the entanglement dynamics, in that some initial block states are frozen or near-frozen, and in density gradients persisting in long-time equilibrated states. Examining the level spacing statistics, we show that the most common Hamiltonians used in FQH physics are not integrable, and explain that GOE statistics (level statistics corresponding to the Gaussian Orthogonal Ensemble) can appear in many cases despite the lack of time-reversal symmetry.


Abstract.
We consider the unitary dynamics of interacting fermions in the lowest Landau level, on spherical and toroidal geometries. The dynamics are driven by the interaction Hamiltonian which, viewed in the basis of single-particle Landau orbitals, contains correlated pair hopping terms in addition to static repulsion. This setting and this type of Hamiltonian has a significant history in numerical studies of fractional quantum Hall (FQH) physics, but the many-body quantum dynamics generated by such correlated hopping has not been explored in detail. We focus on initial states containing all the fermions in one block of orbitals. We characterize in detail how the fermionic liquid spreads out starting from such a state. We identify and explain differences with regular (single-particle) hopping Hamiltonians. Such differences are seen, e.g. in the entanglement dynamics, in that some initial block states are frozen or near-frozen, and in density gradients persisting in long-time equilibrated states. Examining the level spacing statistics, we show that the most common Hamiltonians used in FQH physics are not integrable, and explain that GOE statistics (level statistics corresponding to the Gaussian Orthogonal Ensemble) can appear in many cases despite the lack of time-reversal symmetry.

Introduction
In recent years, intense research effort has been directed at the real-time dynamics of many-body quantum systems in the absence of external baths or dissipation mechanisms [1][2][3][4]. Motivated by experiments -most prominently with laser-cooled atoms, but increasingly also in other settings -where experimental measurements are performed at timescales shorter than dissipation timescales, a large and growing body of theoretical work has addressed issues such as relaxation, thermalization, and dissipationless transport in isolated quantum many-body systems.
The bulk of the literature on non-dissipative many-body quantum dynamics focuses on short-range Hamiltonians with density-density interactions and singleparticle hopping. In particular, the dynamics of one-dimensional (and to a lesser extent two-dimensional) lattice Hamiltonians have enjoyed much attention. Models that can be mapped to free-fermion systems, such as the transverse-field Ising model and hardcore bosons on 1D chains, as well as models that are fundamentally interacting, such as the Heisenberg spin- 1 2 chain and the bosonic and fermionic Hubbard models, have served as widely used test beds for understanding non-equilibrium quantum physics in the dissipationless regime [1][2][3][4].
In the present work, we consider the dynamics of a Hamiltonian which also has a venerable history in the condensed matter physics literature, but has some remarkable properties which stand out from the Hamiltonians more usually considered in the non-equilibrium literature. This Hamiltonian stems from the community studying fractional quantum Hall (FQH) physics. We consider interacting fermions on a two dimensional surface subject to a strong magnetic field, such that the dynamics can be considered to be confined to the lowest Landau level (LLL) at finite filling. Following a strong tradition in numerical studies of FQH physics, the fermions are placed on closed geometries, either a sphere [5][6][7][8][9][10] or a torus [10][11][12][13]. The single-particle basis then consists of a finite number of orbitals circling the sphere (along lines of latitude) or circling one direction of the torus. When the interaction potential is projected onto the LLL, we obtain a quartic term of type a † i+k+m a † i−m a i+k a i , where the subscripts are orbital indices. This contains both static density-density repulsion (m = 0 or m = −k) and, more interestingly, correlated hoppings of particle pairs. When focusing on the LLL, the kinetic energy (quadratic) part of the Hamiltonian is quenched. Thus the dynamics are driven by the correlated pair hopping contained in the interaction. Since the orbitals are arranged in a sequence, one can think of the arrangement as a onedimensional 'lattice' system with open or periodic boundary conditions (sphere and torus respectively). However, the quartic terms driving the dynamics are very peculiar from the perspective of lattice models. Viewed as a lattice (chain), our system is described by a Hamiltonian of the form The hopping terms involve two particles hopping symmetrically outward away from, or symmetrically inward toward, each other. This form of the Hamiltonian appears because the system is described in momentum space, and momentum is conserved. Since momentum in one direction is coupled to position in the transverse direction, the process conserves the center of mass position of the pair. The details of the coupling V (i) k,m depend on the exact geometry, but it has the generic feature that it falls off at long distances (as a function of either k or m) in gaussian fashion (faster than exponentially), with a length scale scaling as a square root of the number of orbitals (number of "lattice sites"). The coupling is thus neither completely local nor quite non-local.
To highlight the peculiarities of this Hamiltonian, we may compare/contrast with a paradigmatic single-particle-hopping Hamiltonian that is used in the non-equilibrium literature: This is a tight-binding model for fermions with nearest-neighbor interactions of strength V , and is equivalent to the XXZ spin chain apart from boundary terms. Although this is an interacting quartic Hamiltonian, like the Hamiltonian (1) to be examined in this work, the hopping occurs through a quadratic term (single-particle hopping). This feature is quite typical of models generally considered in the quantum non-equilibrium literature [1][2][3][4]. Clearly, when considered in momentum space, the dynamics in this model will also proceed through symmetric two particle hopping. However, outside of the FQH context, we can not interpret this as hopping in real space. We will discuss this further in section 6. We consider dynamics arising from an inhomogeneous initial state. We pack the N fermions in a block of N contiguous orbitals, so that there is a block of filled orbitals which spreads into the rest of the sphere/torus, the rest of the orbitals being initially empty. With respect to any subdivision of orbitals, this initial state is a zeroentanglement product state. This has an obvious analogy in fermionic tight-binding chains, and is also analogous to a spin chain where some contiguous block start in a spin-up state and the rest start in a spin-down state, e.g. the 'domain wall' state |↑↑↑↑↓↓↓↓ . For tight-binding fermionic chains and spin chains, this type of dynamics by now has been studied in some detail [14][15][16][17][18][19][20][21][22][23][24][25][26][27][28][29][30][31]. Related situations, e.g. when one/both of the regions are not completely full/empty, are also of interest and sometimes involve similar physics [23,[32][33][34][35][36][37][38][39][40][41]. As the dynamics can be regarded as being initiated by abruptly joining two regions in different equilibrium states, this process is sometimes known as an 'inhomogeneous quench'. If both parts are locally ground states of the Hamiltonian, there is a spreading or 'transport' of fermions or magnetization from one region to the other. For local Hamiltonians, an intermediate current-carrying region grows linearly as a 'light cone', within which a non-equilibrium steady state emerges. The form of the non-equilibrium steady state [14,21,23], the nature of the boundary of the light-cone [22,23], spreading of entanglement [23], and the effects of interactions [15,16,21,28] and integrability [25,27] are some of the many aspects that have been investigated. In this literature, the spreading dynamics is driven by single-particle hopping, as opposed to the correlated pair hopping driving our LLL dynamics. We expect such 'transport' dynamics to be well-suited for exposing the differences between the two classes of hopping.
A hallmark of isolated quantum systems is that the dynamics can be governed by any part of the many-body eigenspectrum: the low-energy parts of the spectrum do not necessarily play an important role. The reason is that, if the initial state has overlaps with some set of eigenstates, possibly far from the low-energy part of the spectrum, then (assuming a time-independent Hamiltonian) the system dynamics does not explore any other part of the spectrum, as there is no tendency toward the ground state in the absence of external dissipation mechanisms. Initial states of the inhomogeneous type we use, typically have largest overlap with eigenstates far removed from the lowenergy regime. The interest in the present Hamiltonian and associated geometries has been largely due to the exotic topological and fractionalization properties of the ground state and low-lying excited states [42][43][44][45][46]. In this work, this class of properties is not expected to play any role, since the physics explored is not directly related to ground state physics. Instead of the low-energy part of the spectrum, we present a study of the level spacing statistics of the full many-body spectrum. The level statistics are known to be sensitive to conservation laws and integrability [47][48][49][50][51], which may have nontrivial effects on the non-equilibrium behaviors.
Although the physics explored here is not directly related to FQH physics, we show results at the filling at which the best-known FQH state, namely the Laughlin ν = 1/3 state [42], appears. We place N fermions on a torus with 3N orbitals or on a sphere with 3N − 2 orbitals. The phenomena we describe do not depend strongly on the filling, as long as the filling is not too close to 0 or 1. Our choice of 1/3 filling simply reflects the historical importance of this filling. The extreme case of unit filling represents an integer quantum Hall state, which has no dynamics as long as we are confined to a single Landau level.
The precise form of the interaction we use also reflects a FQH tradition: we use the so-called pseudopotential V 1 . This is a short-range potential in real space. It has the Laughlin ν = 1/3 state as its exact ground state [7,44,52]. The dynamics we will present depend on the correlated-pair-hopping nature of the Hamiltonian, but do not depend strongly on the precise form of the interaction. When appropriate, we comment on differences or similarities between the short-range V 1 potential and the long-range Coulomb potential, with regards to dynamics and level statistics.
In Section 2, we provide a description of the geometries and the form of the Hamiltonian in the orbital basis. In Section 3, we describe our initial states. Considering initial states where the N fermions are placed on a block of N consecutive orbitals, the position of this block can still be varied. On the torus, shifting the block position makes no difference due to translation symmetry, analogous to a periodic tight-binding chain. On the sphere, one could consider having the block starting at one of the poles. In the chain analogy, this corresponds to the block starting at one edge of the open chain, so that there is a single 'domain wall'; this is strongly analogous to the literature on inhomogeneous quenches. We explain however that such a state is inert and has no dynamics whatsoever, due to conservation of angular momentum. Instead, we primarily focus on the initial state having the block of fermions centered at the equator. Section 4 focuses on real-time dynamics. We present details of the 'melting' process of the equatorial block, such as the intermediate and longtime occupancy/density profiles, the evolution of entanglement entropies, relaxation timescales etc. We point out differences and similarities with a 'usual' tight-binding model with single-particle hoppings; for example the evolution of the entanglement entropy profile shows a striking difference. We examine the 'overlap distribution'the overlaps of the initial state with all energy eigenstates. For the torus, the overlaps are found to be dominated by a few eigenstates at the very top of the many-body spectrum. As a result, our initial state on the torus has no substantial dynamics, even less so with increasing sizes. We also examine a shifted initial state on the sphere, which leads to a non-uniform asymptotic occupancy profile whose overall shape (nearly linear) we are able to explain using mean-field-like considerations.
Section 5 discusses the level spacing statistics. By identifying Wigner-Dyson statistics, we demonstrate that the LLL Hamiltonians are not integrable, either with V 1 or with Coulomb interactions. We demonstrate also that the Hamiltonian on (a) The curves are proportional to the real-space density |ψ| 2 for two adjacent orbitals (out of a total of nine) on a torus. The dots represent the y-position of the density peaks of orbitals. The torus has its orbitals ψ k equidistantly spaced at positions y = kLτ 2 /N φ , where k is the momentum quantum number in the x direction. Each orbital wraps around the torus in the x direction. (b) The curves show the real space densities |ψ| 2 for all (nine) orbitals on a sphere, each normalized as |ψ(θ)/ψ(θmax)| 2 for better overall visibility. Dots are θ positions of the density peaks. The orbitals each wrap around the sphere along lines of equal latitude. The orbitals are not spaced equidistantly, but have peaks at θ = arccos( m Q ); where m is the orbital index running from −Q to +Q.
the sphere and in several momentum sectors on the torus follow GOE (Gaussian orthogonal ensemble) statistics, despite the system having broken time-reversal symmetry. We explain this counter-intuitive effect using a result reported in Ref.
[51] for single-particle systems, namely, that GOE statistics may appear in timereversal-broken systems if there is a unitary symmetry which, combined with timereversal, provides an anti-unitary symmetry leaving the Hamiltonian invariant. To our knowledge, this is the first time this effect has been observed in a many-body setting.
Concluding comments and discussion are presented in Section 6.

Geometries and Hamiltonian
In this section we introduce the geometries and the forms that the interaction Hamiltonian takes in the lowest Landau level on the two geometries. Our goal here is to provide enough details for a reader from outside the field of fractional quantum Hall physics to appreciate the geometry and the forms of the Hamiltonian. We will first describe the single-particle problem and single-particle eigenfunctions (orbitals) on both the sphere and torus, followed by discussion of the interacting problem.
2.1. Single-particle orbitals and quantum numbers A charged particle (possibly an electron) with mass m and charge q, confined to move in two dimensions (x and y), in a magnetic field B pointing in the third (z) direction, is subject to the Hamiltonian H = j=x,y Here A is the vector potential. The quadratic Hamiltonian can be diagonalized by introducing ladder operators a, a † and b, b † . The diagonalized hamiltonian is H = ω a † a + 1 2 , with ω = |q|B/m. The absence of the b operator signifies an extensive degeneracy with respect to each eigenvalue of a. Each energy band is called a Landau Level (LL), and the lowest band is the lowest Landau Level (LLL). We are interested in the high-field regime such that the separation between Landau levels is far larger than the interaction energy scale, so that we can focus on the LLL and ignore all other Landau levels. The single-particle eigenstates within the LLL have a characteristic length scale ('magnetic length') B = /|q|B.

Torus
We define the torus as a Bravais lattice with lattice vectors L 1 and L 2 , not necessarily perpendicular to each other. We define L 1 = (L, 0) and L 2 = (τ 1 L, τ 2 L), where τ = τ 1 + iτ 2 is a complex parameter. This is a standard parametrization of the torus geometry [13,[53][54][55][56][57][58][59]. Rectangular tori, for which L 1 and L 2 are orthogonal, are obtained for τ 1 = 0. In this case, for τ 2 → 0 and τ 2 → ∞, one obtains the 'thin torus' limit, which has been used fruitfully to gain insights into FQH physics [60][61][62][63][64][65][66][67][68][69][70]. We will restrict to unit aspect ratio (|L 1 | = |L 2 | and |τ | = 1). In this case, defining τ = e iθ , a square torus is obtained for θ = π/2, while a trapezoidal torus with 'hexagonal' symmetry is obtained when the shape is distorted to θ = π/3. While the torus has translational invariance in both lattice directions, the geometry within a LL is non-commutative, hence the momenta in the two directions cannot be simultaneously specified [12]. Only one of the lattice momenta can be used as a quantum number to distinguish the degenerate single-particle states in the LLL. The torus area is required to be A = 2πN φ 2 B , where N φ is the integer number of magnetic flux quanta that is piercing the torus.
In the LLL, there are N φ linearly independent momentum eigenstates, labeled by the integer k (single-particle momentum in the x direction). The eigenstates are expressed exactly as sums over displaced Gaussian functions [11]; however the structure of the orbitals is more easily visualized from the following approximation: As shown in Fig. 1(a), these orbitals are extended in the x-direction, and localized in the y-direction at y = −Lτ 2 k/N φ with width B . Note that position in the y direction corresponds to momentum in the x direction -each orbital is localized around a yposition determined by the x-momentum quantum number k. This is another feature of the noncommuntative geometry within LLs. One could of course have chosen to work with eigenstates of the momentum in the y direction; in that case the orbitals would wrap around the torus in y direction and be localized in the x direction around a x-position set by the y-momentum quantum number.

Sphere
The so-called Haldane sphere [5][6][7][8][9][10] has a magnetic monopole in its center, such that it is pierced by N φ = 2Q flux quanta, corresponding to a magnetic field B = 2Qφ0 4πR 2r . Here 2Qφ 0 is the flux through the spherical surface. The Dirac quantization condition dictates that N φ is an integer.
The Hamiltonian of a single particle on the sphere may be written as H = We use the notation (φ, θ) for azimuthal and polar angles. The vector potential A = − cQφ0 eR cot θφ has two Dirac strings, one through the north pole and one through the south pole. Unlike the torus, there is no momentum quantum number, however, a magnetic analogue of angular momentum (L) can be defined [5,8,10]. The Hamiltonian can be rewritten as and are localized around θ = 2 arctan( Q+m Q−m ) = arccos(m/Q). A cartoon of the arrangement of orbitals is shown in Fig. 1(b).

Many interacting fermions
We now consider N interacting fermions. In the original FQH literature, the Coulomb interaction was most natural to consider as it is the physical interaction between electrons. Another interaction that also has a long history is the short-range interaction that produces the Laughlin state [42] as its exact ground state -this is the Trugman-Kivelson interaction ∇ 2 δ( r) [52] projected onto the LLL. This interaction penalizes two particles interacting via the p-wave channel (the s-wave channel being anyway forbidden by fermionic antisymmetry). It is thus also an interaction form that could conceivably be realized in cold-atom experiments with spin-polarized fermions in reduced dimensions [71]. When projected onto the LLL, two-body interactions are conveniently described in terms of pseudopotentials V m [5,7,10,44,61]. The Trugman-Kivelson interaction leads exactly to the V 1 pseudopotential; hence this is also known as the V 1 interaction. Given its historical importance and somewhat easier form, we will focus on the V 1 interaction in this work, including occasional comparisons with the Coulomb interaction.
To treat a central isotropic interaction potential V (| r− r |), on the torus one needs to use a periodicized form,Ṽ ( r) =Ṽ ( r + nL 1 + mL 2 ), while on the sphere the arc distance between r and r is used. In the basis of orbitals, such an interaction (either Coulomb or V 1 ) takes the form . This Hamiltonian includes both static repulsion (for m = 0) and long-range correlated pair hopping terms (m = 0).
On the torus, the interaction is translationally invariant, i.e. V The amplitude of correlated hopping falls off rapidly with k and m but nevertheless remains finite at long ranges, which allows the dynamics to get initiated starting from our block initial states. The length scale of the Gaussian in units of B is set by L, the circumference in the transverse direction. For tori of unit aspect ratio, the interaction length scale thus scales as ∼ N φ . Viewed as an interaction on a lattice with ∼ N φ sites, V (i) k,m cannot be viewed either as short-range or as long-range in any conventional meaning of these terms.
Thus, although the Trugman-Kivelson interaction is short range in real space, when the system is viewed as a tight-binding chain of orbitals, it leads to fermions with correlated hoppings that are not very short-range, as shown in Fig. 2(a).
On the sphere, the interaction V (i) k,m depends on i, reflecting the lack of translation invariance in this geometry. The dependence on k and m is more complicated to write than is the case on the torus, but has similar qualitative form, and is shown in Fig. 2 In particular, we expect the length scale of the interaction to scale on average as ∼ N φ .

Symmetry Sectors
On the sphere, there exists a magnetic analogue of angular momentum L, that preserves all the commutation relations of usual angular momentum operators. We may thus choose to use the L z operators to enumerate the basis of single particle orbitals, and diagonalize the Hamiltonian within this basis. In the many-body system, there now exists a total L z and corresponding total L 2 , which both commute with the Hamiltonian. Block diagonalizing the Hamiltonian with respect to L z is trivial, but dividing the Hilbert space into L 2 sectors is tricky and most often not implemented directly in numerical calculations. Our calculations are performed in fixed total L z sectors of the Hilbert space. Our time evolution is performed with initial states that are not eigenstates of total L 2 , thus multiple L 2 sectors are involved. When separation of L 2 sectors is needed (in the level statistics analysis of Section 5), we add a term αL 2 to the hamiltonian and post-process the data by sorting the eigenvectors by L 2 .
On the torus there are two generators t 1 and t 2 which measure the (magnetic) momentum (k 1 and k 2 ) of the single particle orbitals in the two lattice directions.
These two operators do not commute, but satisfy t 1 t 2 = e i 2π N φ t 2 t 1 , which shows that the maximal sets of commuting operators are t 1 , t For the many-body state we may define the operators 2 , which measure the total momentum K 1 and K 2 of the many body state. If the filling fraction is This shows that the maximal sets of commuting operators are T 1 , T q 2 or T q 1 , T 2 . As a consequence, some combinations of K 1 and K 2 are redundant. Choosing e.g. T 1 , T q 2 , means that K 1 = 0, . . . , N φ − 1, whereas K 2 and K 2 +N give the same T q 2 eigenvalue, reducing K 2 to K 2 = 0, . . . , N −1. Finally, since (K 1 , K 2 ) and (K 1 + N, K 2 ) have the same energy spectrum, we may define a reduced many-body Brillouin zone using (K 1 , K 2 ), where K 1 , K 2 = 0, . . . , N − 1, as in previous work [12,59,[72][73][74].

Initial state
The Hamiltonian is number-conserving, and the number of fermions, N , is thus fixed during dynamics. In considering scaling with system size, we keep the filling fixed. Since our Hamiltonian and geometries are motivated by the literature on fractional quantum Hall states, we use the filling corresponding to the most prominent such state, namely the Laughlin state at filling one-third. For the torus, this means that we have N fermions in exactly 3N orbitals. For the sphere geometry, the Laughlin state appears when N fermions are placed on a sphere with 3N − 2 orbitals. Here 2 is the so-called 'shift' and is present due to the curvature of the sphere [5][6][7][8][9][10].
Our choice of filling means that the ground state is a topological (Laughlin) state. However, since we are looking at non-dissipative dynamics, the ground state or lowenergy segment plays no distinguished role; the eigenstates relevant to our dynamics are far from this part of the spectrum. In fact, similar dynamics is expected for slightly altered filling fractions, even though such a sector would have utterly different properties from the point of view of studying ground-state properties.
We consider initial states which are product states in the orbital basis: the N fermions are placed in N orbitals. We will focus on initial filling of N successive orbitals: This is the analog of fermionic dynamics on a tight-binding lattice where all the fermions are initially packed on adjacent sites [22][23][24]27], or spin-1 2 dynamics starting from a configuration in which all the up-spins are on adjacent sites, i.e. starting from a domain-wall state [14-21, 25, 26, 28, 29, 31].
On the torus, the starting point (j in Eq. (5)) of the block of fermions does not matter; translation symmetry guarantees that the dynamics is equivalent for any starting position. For the sphere, however, the positioning of the string affects the dynamics drastically.
For the sphere, we could consider placing the fermions in the N orbitals closest to one of the 'poles'. In the 'tight-binding chain' picture, this would mean unit filling for the first third or last third of the chain, and zero filling for the rest of the chain. This would then be analogous to the dynamics considered in, e.g. Refs. [14][15][16][17][18][19][20][21][22][23][24][25][26][27][28][29][30][31], for spin or tight-binding systems, when the initial state has all fermions or all ↑-spins on one side of the chain. In such cases, the single-particle hopping leads to transport of particles from the filled side to the empty side, with a current-carrying non-equilibrium steady state arising dynamically within the light cone region. In our system, however, this type of initial state leads to an inert state with no dynamics whatsoever.
The lack of dynamics starting from a polar-cap initial state on the sphere can be understood in several ways. The sphere Hamiltonian conserves the L z quantum number (introduced in Sec. 2), which is maximal for this state. That means that this state belongs to a symmetry sector of the Hilbert space composed of a single state, namely the polar-cap state. A 1-dimensional Hilbert space, of course, has no dynamics. Alternatively, one can consider the correlated pair hopping process, which can occur only if two particles can hop symmetrically outwards or inwards, Fig. 2(a). When the particles are packed on the first N orbitals starting from the boundary, only hoppings in one direction are possible, hence there is no configuration space available for pair hopping, resulting in a frozen state. (The situation can be described using terrestrial analogy: 'The polar-cap state is frozen.') Therefore, we consider dynamics where the block of N fermions are near the equator. In fact we will mostly concentrate on the initial state covering a symmetric band around the equator. However, in Section 4.4 we also show some interesting dynamical effects for 'lopsided' initial states.
While the orbital description is very convenient, one might also want to know what our initial states look like in real space. Since the orbitals are loosely localized in space, one expects that, at large N , a region of filled orbitals corresponds to a region with a nearly constant real space density equal to the maximal density possible in the lowest Landau level. This density is the real space density observed in the integer quantum Hall state with a completely filled lowest Landau level, ν = 1. Fig. 3 shows the real space density for various particle numbers. While this correspondence is correct for large N , there are significant deviations for N ∼ 10, i.e. for system sizes accessible numerically. The initial real space density profile then has a bell-shaped form.
In Fig. 3, for large enough N , the particles in the initial state fill up a region near the equator which covers less than one-third of the azimuthal angle. The reason is that orbitals in the equatorial region (which are longer) have smaller width than orbitals in the polar region (which are shorter), because each orbital must cover the same area. If there is a plateau, the real-space density falls off at both edges of the plateau with a fall-off width of the order of a magnetic length B ; the decay is faster than exponential.
Although the large-N behavior has a plateau structure, for the dynamical effects presented in this paper, it is simply more appropriate to think of the initial state having a gaussian-like real-space density profile, corresponding to the numerically accessible N values.

Time evolution
In this section, we present our results on the time evolution of initial states of the type (5). Most of the analysis is for the equatorial initial state on the sphere.
We first describe the evolution of the orbital occupancies (4.1), considering snapshots of evolving orbital occupancy profiles and the deviation from the final uniform state. Subsection 4.2 describes the same evolution in terms of the realspace densities. We then present the evolution of the entanglement entropy Initial density profile Figure 3. The real space density profile ρ(θ) of the equatorial initial state on the sphere, for various system sizes. To make the comparison of different N possible, the density ρ(θ) is re-scaled by the corresponding density ρ ν=1 = 2Q+1 that the filled LLL (ν = 1 ) would have. For larger systems a well defined density plateau emerges near the equator over the angles θ − < θ < θ + , where θ ± ≈ 2 arccos(∓ν) = arccos(∓ 1 3 ). The density decreases from ρ = ρ ν=1 to ρ = 0 over a distance of the order of the magnetic length B . For system sizes numerically treated in this work (N ≈ 10), such a density plateau has not yet formed.
in particular a contrast with the way entanglement entropy grows in local models with single-particle hopping. These three sections focus on the equatorial initial state. Subsection 4.4 considers a 'lopsided' initial state, displaced from the equator. The correlated-pair-hopping nature of the Hamiltonian (or equivalently, L z conservation) now manifests itself in a remarkable manner -the long-time state is also lopsided, and for a range of displacements has a remarkably linear form. We explain the near-linear shape using mean-field arguments for the pair-hopping Hamiltonian.
The section ends with a presentation of the overlaps of the equatorial initial state with all eigenstates of the Hamiltonian, the so-called 'overlap distribution' (4.5). This is compared with a tight-binding model and the torus case. In the torus case, the form of the overlap distribution leads to the conclusion that the block initial state on the torus is frozen for larger sizes; we confirm this by showing orbital occupancy dynamics for a sequence of sizes.
The numerical results presented in this work are all obtained by exact numerical diagonalization. To present the numerical data, we have used the gap ∆ in the ground state sector (the so-called 'neutral gap') as the unit of energy. This is the energy difference between the Laughlin (ground) state and the first excited state in the same L z sector or the same (K 1 , K 2 ) sector. This is a well-studied energy scale and is of the same order of magnitude as the larger interaction terms (static or pair-hopping); hence we regard it as a reasonable analogue of the usual practice in the non-equilibrium literature of using the hopping term to define the units of energy. Time is measured in units of /∆ = 1/∆.

Evolution of orbital occupancies
In Fig. 4(a,b), the dynamics starting from the equatorial initial state on the sphere is shown through successive snapshots of the orbital occupancies. The two panels show the same dynamics at shorter and at longer timescales. The fermions gradually spread out from their initial position and occupy all the orbitals. In the very long-time limit, the occupancies approach a 'flat' distribution in which each orbital is equally occupied; this is shown also through individual orbital occupancy evolutions in Fig. 4(c). Of course, for any finite size, the occupancies cannot relax completely and at long times fluctuate around the average long-time distribution.
It is tempting to compare the melting at either edge to 'domain wall' melting in inhomogeneous quenches driven by a regular tight-binding Hamiltonian with nearest-neighbor single-particle hopping [14,[21][22][23]. Detailed similarity is perhaps not expected, because in our case the melting of the two domain walls are correlated and not independent local processes. In the dynamics here, the two discontinuities in orbital occupancy survive for quite a while, until t ∼ 8, Fig. 4(a), and it is only after this time that a smooth profile is obtained. (With nearest-neighbor singleparticle hopping, the discontinuity generally smoothens out rapidly.) Our interaction is not strictly short-range and hence a clean light cone is not expected. Nevertheless, the spreading out is gradual -the empty orbitals near the band get occupied first, and then the occupancy spreads gradually toward the poles. This presumably reflects the fact that the V km matrix elements generally fall off with m (Fig. 2). As a clear light cone is not expected, it is unclear whether or not at larger sizes a spatially distinguishable region of current-carrying 'nonequilibrium steady state' [14,[21][22][23] emerges. For the sizes accessible in our calculations, this question cannot be meaningfully addressed.
One effect of our peculiar (correlated-hopping) interaction is visible at very short times (t 1) -the orbital occupancy profile is non-monotonic at these times. Since V km falls off with both k (except for very small values of k) and m, hopping of pairs from farther inside the band might be more favorable than hopping from the very edges. As the two effects -dependence on k and dependence on m -compete, one can obtain a non-monotonic effect; this shows up in the short-time profile. At larger times, the effect is washed out. This mechanism for non-monotonicity is described in more detail in Appendix A.
Another contrast with spreading dynamics driven by single-particle hopping is that, in our case, we do not observe any reflection after the fermionic liquid reaches the poles (edges). In propagation dynamics on regular tight-binding chains, there is almost always reflection, even in the presence of strong interactions.
We now turn to the long-term dynamics, Fig. 4(b,c). The occupancies approach a uniform distribution at long times. The occupancy at any orbital approaches the average occupancy value N/(3N − 2) and then oscillates around this value, Fig. 4(c). The relaxation to the average occupancy value is faster in the region near the boundary of the initial state, away from both equator and poles. The orbitals near the pole take a noticeably long time to settle to the final value, Fig. 4(c). This may be regarded as an effect of the correlated-hopping nature of the Hamiltonian. Those V km processes which are symmetric around the equator would need to have a large k or m if they were to populate the polar-region orbitals, and are hence suppressed. For non-symmetric V km processes with smaller k and m, the partner that also hops would have to hop from a smaller density to a larger density region; such processes are thus also not very effective.
The latter effect can be formulated intuitively in 'hydrodynamic' language as follows. For a single-particle-hopping Hamiltonian, density gradients drive particle flow, as hoppings from a high-density to a low-density site is kinematically favored. On the other hand, a correlated-pair-hopping process involves two particles hopping away from (or in toward) a central region; this will be favored if the densities are lower (resp. higher) on both sides of the central region. Thus correlated hopping leads to flow driven by second derivatives of density. This intuitively explains why peaked structures like the short-time density profile lead to rapid spreading out, but the final density inhomogeneity near the ends of the 'chain' are not so easily corrected by the type of hopping we have. Later in Section 4.4, we will explicitly show a situation where a density gradient leads to no flow.
The approach to the final distribution can be visualized by plotting the deviation from the uniform distribution. This deviation is quantified, for example, by the standard deviation σ of the orbital occupancies. The time evolution of this quantity is shown for several sizes in Fig. 5. The initial fast decrease of the variance, Fig. 5(a), is roughly proportional to the system size, i.e. proportional to the number of orbitals the fermions must spread into. This is similar to what one would expect for a short-range hopping, in which case there is a light-cone effect of information/particles spreading out at a linear rate. We show here that the same feature is visible for our longrange correlated hopping; presumably this aspect is not very different from shortrange models because of the rapid decay of matrix elements V km with large k or m. In Fig. 5(b), there appears to be a second slower relaxation period, which lasts until some timescale increasing rapidly (possibly exponentially) with system size; unfortunately the available sizes do not allow a clean investigation of this time scale.
At longer times the orbital occupancies each fluctuate around the average value, which is seen in Fig. 5(b) as the final fluctuating behavior of the standard deviation σ. The strength of the fluctuation decreases with system size. The data is consistent with an exponential decrease with the system size -on the logarithmic plot scale, the average final value of σ (shown using left-pointing arrows) decreases in roughly equal steps as N is increased. The long-time value of the standard deviation is an average measure of the fluctuation strength of the individual orbital occupancies. In the recent literature on the non-equilibrium dynamics of local observables, the longtime fluctuations have been studied in various systems [75][76][77][78][79][80][81][82][83][84][85]. The magnitude of fluctuations around the long-time average value is understood to generically decrease exponentially with system size. This fluctuation strength is determined by the average magnitude of the off-diagonal matrix elements of local observables between eigenstates of the Hamiltonian [75,[86][87][88][89][90][91][92]. At least for non-integrable Hamiltonians, this decreases exponentially with system size (or more precisely, as D −1/2 , where D is the Hilbert space dimension), leading to an exponential decrease of the long-time fluctuation, in our case measured by the long-time average of σ. The size-dependence of the longtime fluctuations in our system, Fig. 5(b), is thus similar to the generically expected behavior for non-integrable systems.

Evolution of real space densities
In Fig. 6 the dynamics are visualized through real-space density profiles. The realspace profiles are obtained straightforwardly from the orbital occupancies n i using the real-space wavefunctions ψ i ( r) of the single-particle orbitals: the real-space density at the location r is ρ( r) = i n i |ψ i ( r)| 2 .
For the system sizes we can reach, the initial profile (real-space profile of the initial state) is a bell-shaped curve without a plateau (Sec. 3). The density spreads out over the full sphere and is eventually near-constant.
As noticed with the orbital density profiles, the density in-between the equator and poles relaxes faster, while the densities at the equator and pole take longer to relax, and also perform persistent long-term oscillations.
Remarkably, the non-monotonicity at short times in the orbital occupancy profile is washed out by the orbital profiles and does not get reflected in the real space density profiles. It is perhaps reassuring that the widely used V 1 potential does not lead to artificial effects when viewed in physical space, despite its singular (derivative of a delta function) form.

Entanglement dynamics
We now examine the entanglement generated in the dynamics.
We calculate the bipartite entanglement entropy between two blocks of orbitals (conventionally referred to as A and B), the first partition containing the first block of l A orbitals starting from one pole of the sphere and the second partition containing the rest of the orbitals. This type of 'orbital partitioning' entanglement has been widely studied since the beginning of entanglement investigations of fractional , starting with an analogous initial state.
quantum Hall systems, including studies of both entanglement entropies [93][94][95][96][97][98][99][100][101][102] and entanglement spectra [96,99,100,[103][104][105][106][107][108][109][110][111][112][113][114]. The orbital entanglement may be regarded as a reasonable proxy for entanglement between real-space partitions of the sphere [115][116][117], because of the sequential spatial arrangement of orbital positions. (However, because the width of orbital wavefunctions is larger than their mutual spacing, this correspondence is nontrivial.) The interest in the orbital entanglement entropy/spectrum of ground state quantum Hall wavefunctions is largely due to topological information contained in these quantities. Topological aspects do not play a role in the present work. On the other hand, the evolution of entanglement is also widely used to characterize non-equilibrium dynamics; this is the spirit of our analysis below. Fig. 7(a) shows snapshots of the spatial profile of the entanglement entropy, obtained by plotting the entanglement entropy against the cut position, l A . As we have done until now, we focus on the equatorial initial state: the N = 8 orbitals nearest to the equator are filled, on a sphere with 3N −2 = 22 orbitals. At t = 0 the entanglement entropy is zero, since the initial state is a product state. The entanglement entropy grows with time; the initial growth for almost all cuts in the equatorial region is approximately linear. This is seen explicitly in Fig. 7(b), where the entanglement entropy for the half-partition (l A = 1 2 (3N − 2)) is plotted as a function of time. After t ∼ 10 the entanglement entropy saturates. The saturation time grows approximately linearly with system size (Fig. 7(b)) and correlates with the time it takes for the quantum fluid to reach the north and south pole of the sphere (Fig. 4).
At long times, the entanglement entropy profile eventually settles to values consistent with the entanglement entropy in high-energy eigenstates, or equivalently, with the entanglement entropy in random or 'thermal' eigenstates [118][119][120][121][122][123][124][125]. This is seen in the right panel, where we have plotted the orbital entanglement entropy profiles (entanglement entropy versus l A ) for the 20 eigenstates that are closest in energy to the expectation energy of the time-evolving (or initial) state. The long-time entanglement profile follows the eigenstate entanglement profile quite closely.
It is interesting to compare the long-time (or eigenstate) entanglement profile with the ground state entanglement entropy of the Laughlin state, studied first in detail in Ref. [93]. The ground state entanglement entropy follows an area law [126][127][128][129] and encodes topological properties through its subleading term, the so-called topological entanglement entropy [130,131]. The long-time entanglement entropy after a quench, in contrast, should follow a volume law, and is more similar to a typical random state than to a low-energy state. Indeed the long-time entanglement profile is seen in Fig. 7(c) to be significantly larger than the Laughlin (ground state) profile. Distinguishing volume law from area law quantitatively in finite size data however requires careful finite-size scaling. For the quantum Hall geometry on the sphere, area law and volume law behavior mean that the top of the entanglement profile scales as ∼ √ N and as ∼ N respectively. Our quench data, Fig. 7(b), shows that the entanglement entropy for the half-system cut indeed saturates at values growing linearly with system size.
We now return to the short-time behavior, which displays a striking effect of the nature of our Hamiltonian. At very short times, the entanglement entropy profile is flat in the equatorial regime where the product state fermions are initially placedthere is no peak at the two edges. This is in sharp contrast to what one finds in the case of regular (uncorrelated) hopping, as we display in Fig. 7(d), where entanglement profiles are shown for a tight-binding chain starting from the corresponding initial state -8 fermions placed in the middle of a 22-site chain. (The chain Hamiltonian is that given in the Introduction, Eq. (2), with open boundary conditions and V = 1.) An uncorrelated-hopping Hamiltonian, with hopping terms of the form c † c rather than c † c † cc, generates entanglement first at the boundaries between filled and unfilled regions. The flat entanglement profile generated in our case is a consequence of the correlated hopping processes governing the dynamics in the lowest Landau level Hamiltonian.
With correlated hopping, a fermion hops from the initial central region to an initially empty site only when there is a partner fermion hopping across the other boundary of the initial filled region. According to intuitive ideas of entanglement, this correlated process clearly will generate entanglement between the two edge regions, and will show up as nonzero entanglement across a cut at the center. Thus correlated hopping automatically generates entanglement across a cut placed between the two edges; hence the flat profile at very short times. In contrast, for uncorrelated hopping, a hop on the left edge and a hop at the right edge are independent; these processes at short times generate entanglement across cuts near either edge but no entanglement across cuts halfway between the edges. Thus the short-time entanglement for a tightbinding chain has separated peaks near either edge, as seen in Fig. 7(d).
This intuitive picture can be made more precise through a toy calculation. Using the notation 0 n 1 l 0 n for our initial state (a block of l = N fermions sandwiched between two blocks of n empty sites), a tight-binding nearest-neighbor uncorrelated hopping Hamiltonian generates the following wavefunction at linear order in time: where α is a small parameter proportional to time, and N 1 = 1 + 2|α| 2 −1/2 normalizes the wavefunction. With this simplified wavefunction, for a cut between the edges we obtain the reduced density matrix for one partition to be ρ A = N 2 1 1 + |α| 2 α α * |α| 2 . This leads to zero entanglement at leading order (order |α| 2 ). A cut at either edge, on the other hand, leads to finite entanglement at leading order. This explains the double-peak structure at small times in Fig. 7(d).
To model the correlated hopping effect tractably, we make the simplification that correlated nearest-neighbor hoppings (m = 1 terms in Fig. 2) dominate. The shortterm wavefunction is then where γ is a small parameter proportional to time, and N 2 = 1 + |γ| 2 −1/2 . Here, hoppings across the two edges do not generate separate terms, but are correlated. Now, for a cut between the edges we obtain the reduced density matrix N 2 2 1 0 0 |γ| 2 , yielding a nonzero entanglement entropy at leading order, O |γ| 2 . A cut at the edge yields an identical entanglement entropy at leading order. This explains the flat nonzero entanglement profile in the region between the edges for the lowest-Landaulevel Hamiltonian.

Lopsided initial state
Until now, we have described dynamics starting from the equatorial initial state, 0 l 1 N 0 l , with l ≈ N , where the filled block of orbitals is placed symmetrically around the equator. In this case the total angular momentum is zero, L z = 0. Since L z corresponds roughly to total 'latitude' (angular position) on the sphere and is conserved, the state at any time is also symmetric around the equator. The long-time state, having uniform filling on all orbitals, is obviously symmetric as well.
We now consider initial states where the block is displaced from the symmetric location. Taken to the extreme, this results in the polar block (maximal L z ) which we have explained is an inert state (Section 3). We now consider initial states intermediate between these two cases, such that L z is nonzero but not maximal. Since the time evolution conserves L z , the state continues to be lopsided and biased toward one pole.
We first examine moderate displacements from the symmetric situation, 0 l1 1 N 0 l2 , with l 1 ≈ 1 2 N and l 2 ≈ 3 2 N . Results for such a case are shown in Fig. 8(a,b). As explained above, the density profile remains lopsided due to L z conservation -the long-time asymptotic profile is not uniform. This is true for both the orbital occupancy and the real-space density.
The asymptotic occupancy profile in Fig. 8(a) is, remarkably, nearly linear. While L z conservation guarantees an inhomogeneous profile, it does not by itself predict this simple shape for the final profile. We can interpret the linear shape using a 'hydrodynamic' or 'rate equation' approach for the correlated hopping term. Let us consider four orbitals such that correlated hopping can occur from the outer two orbitals to the inner two, or vice versa. The situation is depicted in Fig. 8(c) If we parametrize the differences of densities as n 1 = n 2 + x = (n 3 + ) + x = (n 4 + y) + + x, then a linear occupancy profile would correspond to x = y, because the 1-2 and 3-4 distances are equal. We can work in the limit where x, y, are much smaller than n i , i.e. the four orbitals are not very far apart. At first order in this approximation, the rate equation (6) yields n 2 4 (1−3 −2x−4y)+n 4 (x+ +y) = n 2 4 (1−3 −x−5y)+n 4 ( +2y)(8) so that we obtain x = y. Thus the approximate hydrodynamic description predicts a linear occupancy profile. (The difference is not determined by this calculation, which is consistent as the argument is independent of the distance k. Near-linearity of the final distribution is a nontrivial consequence of the symmetric correlated hopping, coupled with the geometry allowing for a hydrodynamic description. Applying the same type of argument to a non-correlated single-particle hopping term (between orbitals 1 and 2) gives n 1 (1 − n 2 ) = n 2 (1 − n 1 ), i.e. n 1 = n 2 , Lopsided N = 11 Figure 9. Snapshots of the (a) orbital occupancy profile and (b) real space density distribution, for N = 11 electrons on a sphere. The initial state is even more lopsided than in Fig. 8.
a flat density profile. Of course, it is difficult to imagine an analog of L z conservation with a single-particle hopping term; thus a non-homogeneous asymptotic state is not expected. By squeezing the initial block of fermions closer to the pole, we can remove the geometric freedom required for a mean-field description. This is shown in Fig. 9. This may be regarded as an intermediate case between the moderately lopsided case where the occupancy profile approaches a linear form, and the polar initial state for which the density profile never changes. In the present case, there is some dynamics (the Hilbert space is smaller than equatorial or moderately lopsided cases, but nevertheless finite). However, the long-time state now has a strong remnant of the initial peak in the occupancy profile.

Overlap distribution and torus dynamics
The dynamics of a quantum state is governed by the eigenstates of Hamiltonian with which the state has significant overlap. The distribution of overlaps -the overlaps with eigenstates as a function of corresponding eigenenergies, thus has been examined in the non-equilibrium literature for various Hamiltonians and initial states [75,86,[132][133][134][135][136][137][138][139][140].
In Fig. 10(a) we show the overlaps of the equatorial initial state with the eigenstates of the Hamiltonian. We consider overlaps with all eigenstates in the L z = 0 sector. Eigenstates with any L 2 can have nonzero overlap with the equatorial initial state. The energy of the equatorial state, indicated with a black vertical line, is far from the edges of the spectrum, but nearer to the upper end of the spectrum than the lower end. The highest overlaps come from eigenstates in this region, but not all eigenstates in this region have a large overlap. The largest overlaps in other energy regions follow an approximately exponential envelope, ∼ e λ(E−Einit) , for E < E init , and an approximately Gaussian envelope, ∼ e −(E−Einit+λ) 2 /σ 2 , for E > E init .
The overlap distribution is similar (in terms of overall shape and exponential decrease) to the overlap distribution for the tight-binding interacting chain, shown for comparison in Fig. 10(b). (This is the same chain of interacting fermions with single-particle uncorrelated hopping, used previously in Fig. 7. The initial state is N fermions in the central N sites among 3N − 2 sites.) This similarity suggests that the overall shape and the exponential envelope is not due to the correlated-hopping part of the Hamiltonian. Rather, it is due to the static interaction part (m = 0 terms of V km ), which prefers to keep the fermions apart from each other. Lower energy eigenstates have smaller contributions from configurations with large 'static' interaction energy.
The Gaussian envelope at larger energies follows a well-distinguished set of eigenstates. Such eigenstate branches appear to be common in overlap distributions, e.g. in Figure 5a of [138] and in Figure 3a of Ref. [140]. In the present case, states in this branch have much smaller overlaps than the dominant eigenstates, and hence do not play any noticeable role in the dynamics, unlike the case of Ref. [140].
In Fig. 10(b) we show the overlap distribution for the analogous initial state on the torus. The dominant eigenstates are now the few (around N ) highest-energy eigenstates. These eigenstates have the structure of linear combinations of product states like our initial state, centered at different positions of the torus. These few eigenstates become increasingly dominant with increasing system size. A result of this extreme overlap distribution is that the dynamics is frozen. This is shown in Fig. 10(d). The effect is more extreme for larger sizes.
Since the static part of the interaction energy is minimized by maximizing the spacing between the fermions, it is expected that our spatially packed initial states should have energy nearer to the top of the spectrum. On the torus, this type of state appear at the very top of the spectrum. On the sphere, one can create states with energies even higher than our equatorial initial states. (For example, since packed configurations near the poles are more costly than packed configurations near the equator, we can pack N/2 fermions near the north pole and the rest near the south pole, obtaining an L z = 0 configuration with energy higher than the equatorial state.) This explains why the energy of the equatorial initial state in Fig. 10(a) is not at the very top of the energy spectrum.

Level statistics -Integrability and universality classes
We now consider the level spacing statistics of the fermionic lowest-Landau-level Hamiltonians.
Considerations of level statistics pervade the quantum dynamics literature, both in the field of single-particle quantum chaos [47][48][49][50] and in the study of many-body quantum dynamics [3]. The main reason is that integrability manifests itself in the level spacing distribution, and quantum dynamics can be strongly affected by the presence or proximity of integrability. The spectrum of an integrable system has Poissonian spacing statistics, whereas the spectrum of a non-integrable system will show Wigner-Dyson (GOE, GUE or GSE) statistics [47][48][49][50]. Of course, it is assumed that the energy spectrum is first sorted according to symmetry sectors. One can thus check for integrability using an analysis of level statistics. As far as we are aware, no study of the level statistics of FQH systems has been published to date. We present such a study now. Our analysis shows explicitly that the LLL-projected Hamiltonians are not integrable in the geometries considered. Since a magnetic field violates time reversal symmetry, GUE (Gaussian Unitary Ensemble) statistics would be the obvious expectation, as GUE is generally associated with the lack of time reversal invariance. However, we show and explain below that GOE statistics appears for the sphere and for several sectors of the torus.
For the LLL-projected Coulomb Hamiltonians, there is no a priori reason to expect integrability. However, it is not easy to rule out the possibility, as the projection to the LLL is a nontrivial operation. (It is well known that projection onto Wannier levels can drastically affect integrability properties, e.g. the continuum Lieb-Liniger model is integrable while the lattice Bose-Hubbard model is not.) There is also the possibility that the pseudopotential (V 1 ) Hamiltonian could be integrable in one of the geometries. Also, it is not clear that (non-)integrability in one of the common geometries (sphere, torus) implies the same for the other. It is worth noting that the similar low-energy properties of these different situations (Coulomb and V 1 , sphere and torus) do not imply that integrability properties are the same. To our knowledge, the integrability of LLL Hamiltonians has not until now been examined explicitly. In our analysis described below and displayed in Fig. 11 and Fig. 12, we show that each symmetry sector of the sphere and the torus displays Wigner-Dyson statistics, for the filling appropriate for the Laughlin-1/3 state. We display these results for the V 1 potential that we have focused on in this work, but we have also checked We characterize the level statistics through calculations of the ratio of consecutive level spacings r [141,142]. This measure is based on s n = E n+1 − E n , the set of level spacings in an ordered list E n of eigenenergies. The ratio r n = min(s n , s n−1 ) max(s n , s n−1 ) is defined for each pair of consecutive level spacings. The statistics of the ratio r n is more convenient than the bare level spacings s n because it bypasses the need to account for varying density of states through unfolding procedures. For Poisson statistics, the probability distribution of r n is P (r) = 2/(1 + r) 2 with mean r = 2 ln 2 − 1 ≈ 0.39. For the Wigner-Dyson ensembles, the probability distributions are well-approximated by the surmise [142] P (r) ∝ (r + r 2 ) β /(1 + r + r 2 ) 1+3β/2 up to normalization, with β = 1 (β = 2) for GOE (GUE). The averages are r GOE ≈ 0.53 and r GUE ≈ 0.60 [142]; these values provide a quick check for the nature of numerically obtained spectra.

Sphere
In Fig. 11(a) we show the distribution of r values for the spectrum of the V 1 Hamiltonian on the sphere. One has to first separate the symmetry sectors. Numerical diagonalization is naturally done in distinct L z sectors. We consider here the L z = 0 sector which contains the equatorial initial state and the Laughlin state. However a single L z sector contains multiple L 2 sectors. The L z = 0 spectrum was post-processed into separate spectra for each L 2 sector (by numerically applying the L 2 operator on the eigenstates), and the r values for each sector were collected separately. Under the assumption that each L 2 sector possesses the same statistics, it is reasonable to combine the groups of r values and present the combined distribution, as we have done. Combining the level spacings (s values) themselves from different sectors would have been far trickier as one has to take into account the different density of states effects in the different sectors. Considering the ratios r thus allows us to obtain reasonable statistics even when the individual sector Hilbert spaces are quite modestsized. The distribution P (r) obtained in this manner very clearly follows the GOE form. The excellent agreement provides a posteriori justification for combining data from different L 2 sectors -if different sectors had different statistics, one would expect the numerical P (r) histogram to be intermediate between two or more of the standard distributions.
We have thus shown that the sphere Hamiltonian has GOE statistics. (We have checked that GOE statistics is also obtained in other L z sectors.) Level statistics following one of the Wigner-Dyson distributions demonstrates the Hamiltonian is non-integrable. To explain why the statistics is GOE rather than GUE, one has to consider symmetries of the system. The conventional expectation is that Hamiltonians with time-reversal symmetry have GOE statistics, while those with broken timereversal symmetry have GUE statistics. In most usual cases, the Hamiltonian matrix elements are real in the first situation and complex in the second situation. This latter distinction however, cannot obviously be rigid as any complex Hamiltonian matrix can be basis-transformed to have real matrix elements, and vice versa. Clearly, since we have GOE statistics in a situation with a magnetic field (which breaks time-reversal symmetry), the current situation is more involved than the simple picture above.
The explanation is that, if the Hamiltonian breaks time-reversal symmetry but is invariant under another anti-unitary transformation (e.g. time-reversal coupled with a reflection), then the statistics is GOE. This was first pointed out by Robnik and Berry [51] in the single-particle context. To our knowledge, this is the first time a many-body example has been discussed. In the case of the Haldane sphere, a combination of time-reversal and a reflection around the equatorial plane keeps the system invariant. Time-reversal reverses the directions of cyclotron motion along each orbital, i.e. it changes the sign of L z on each orbital without changing the position, so that the orbitals in the northern hemisphere now have negative L z . (This is equivalent to changing the sign of the monopole at the center of the sphere.) A reflection around the equatorial plane switches the positions of the orbitals. Thus the combination of the two operations is an anti-unitary operation that keeps the system invariant. In this manner, we obtain GOE statistics in a many-body system violating time-reversal mechanism, through the mechanism of Ref. [51].

Torus
In Fig. 11(b) we summarize some results for the various torus geometries. A different distribution appears in each momentum sector (as defined in 2.   r values for the N = 8 spectrum in a selection of toroidal geometries. Each colored square/hexagon is a momentum sector (K 1 , K 2 ). Dashed black lines show the borders of the many body Brillouin zone. The darkest blue sectors are (K 1 , K 2 ) = (0, 0), (N/2, 0), (0, N/2), (N/2, N/2), the latter three existing only when N is even. Here there is a residual symmetry of orbital inversion, which superficially makes the statistics look Poissionian ( r ≈ 0.386, purple). In the other sectors, the r value is clearly identifiable as either GOE ( r ≈ 0.536, green) or GUE ( r ≈ 0.603, yellow). GOE statistics is found in momentum sectors which are invariant under the product of mirror and timereversal symmetry. These sectors are located along the mirror axes of each geometry. All other momentum sectors have GUE statistics. To illustrate this effect, we have represented tori corresponding to all possible Bravais lattices for the system size N = 8. The red dashed lines are reflection axes.
in the numerical diagonalization.) The r values are close to either the GOE or the GUE values for every momentum sector shown in Fig. 11(b), and for every torus orientation. This demonstrates that the Hamiltonians are non-integrable. Whether r is close to GOE or GUE expectation reflects the presence or absence of an antiunitary combination of operators which preserves the Hamiltonian in that sector.
The torus sectors are presented in more detail in Fig. 12. In the sectors (K 1 , K 2 ) = (0, 0), (N/2, 0), (0, N/2), (N/2, N/2), there is an additional symmetry (orbital inversion) not resolved in the numerical diagonalization. The combination of spectra from two symmetry sectors results in a combined spectrum which has no level repulsion, and the resultant distribution looks closer to Poissonian than Wigner-Dyson. This is reflected in the corresponding squares/hexagons in Fig. 12 having near-Poissonian values of r .
The Hamiltonian in most (K 1 , K 2 ) sectors has GUE statistics, reflecting the broken time-reversal symmetry due to the magnetic field. However, a smaller number of momentum sectors has GOE statistics in spite of the broken time-reversal symmetry: for example the momentum sectors K 1 = ±K 2 for tori of aspect ratio 1. In this case, time-reversal combined with a reflection along the diagonal keeps the Hamiltonian invariant. More generically, we note that the momentum sectors with GOE statistics are determined by the Bravais lattice symmetry of the considered geometry. These sectors indeed coincide with the reflection axes of each torus, such that the product of time-reversal and reflection leaves the Hamiltonian unchanged in these sectors. To verify this statement, we have obtained the level-spacing statistics in geometries representative of all possible types Bravais lattices. The results for the largest size (N = 8) are given in Fig. 12.
A noteworthy feature is that GOE and GUE statistics do not correspond to real and complex Hamiltonians. For the sake of concreteness, let us consider an example in the square geometry. In that case, the sectors (K 1 , K 2 ) = (K, ±K) (K = 0) have GOE statistics even though the Hamiltonian has complex matrix elements.

Discussion and Context
We have presented a detailed exploration of non-equilibrium dynamics and levelspacing statistics of an interacting fermion system in the lowest Landau level, on the commonly used sphere and torus geometries. The dynamics is initiated in a state with a consecutive block of filled orbitals.
Quench dynamics of interacting fermions in the LLL has been recently studied [143], with the quench changing the torus aspect ratio. This type of quench is sensitive to low-energy physics (such as magnetorotons), in contrast to our setup where the part of the many-body spectrum that is relevant is closer to the top of the spectrum than the bottom.
Our setup is specific for the finite geometries (sphere and torus) that are common in numerical diagonalization, so some care is required in interpreting this situation in the thermodynamic limit, and even more care would be needed to relate to experiments. Our initial state roughly mimics a puddle or a stripe of fermions in a magnetic field (a constant density of fermions in the interior of the puddle is achieved for large enough system sizes). The state is designed to be such that the only dynamics is due to the interactions (i.e. the initial state is an eigenstate of the magnetic field and kinetic energy parts of the Hamiltonian). While this last feature of the initial state may seem unnatural in a continuum context, electrons confined to a stripeshaped region will approximately fill the orbitals in that region, and hence our initial state may be a good zeroth approximation. Since we consider the full many-body spectrum of the LLL Hamiltonian, we are focusing on the parameter regime in which the magnetic field is so large that the Landau level splitting is much larger than the many-body bandwidth.
Viewed as a one-dimensional (1D) lattice system, our setup bears resemblance to inhomogeneous quenches starting from domain-wall initial states, which is a topic of extensive interest with more conventional 1D lattice Hamiltonians. The perspective of the FQH interaction being treated as a 1D chain Hamiltonian pair-hopping has been adopted previously, e.g. Refs. [144][145][146][147]. In some cases [145][146][147], the proximity to the thin torus limit allows one to truncate the longer-range terms, leading to a more manageable or even integrable Hamiltonians. Such a truncated Hamiltonian would not give interesting dynamics in our case, as the pair-hopping has to be long-range enough to initiate the melting.
Viewed as a 1D chain, our interactions have the peculiarity of being neither local nor non-local in the conventional sense. The "interaction range" (length scale of the gaussian fall-off) grows as the square root of the number of orbitals. Thus the range grows indefinitely with system size but the fraction of the system covered by the interaction decreases with system size. An interesting aspect of the ∼ N φ growth of the interaction range is that, for large enough system sizes, the equatorial initial state would be frozen (up to exponentially long times) if the filling fraction were kept unchanged. If the width of the initial block grows linearly as ∼ N φ /3, the interactions will not be long-ranged enough to initiate melting at larger sizes. However, as long as the block size is kept N φ , spreading dynamics at reasonable time scales is expected at larger N φ as well.
We have stressed the correlated-pair-hopping nature of our dynamics-driving terms, contrasting the dynamical effects of these quartic terms with dynamics driven by quadratic single-particle-hopping Hamiltonians. Correlated-pair hopping has been rarely addressed explicitly in the quantum dynamics literature. Interacting Sachdev-Ye-Kitaev models (SYK 4 ), for example, contain correlated pair hopping terms of random strengths. The relaxation dynamics and level spacing statistics of SYK 4 models have been studied recently [92,[148][149][150]. However, these models have no spatial structure, so there are no transport/spreading effects. Also, we are not aware of any dynamical features explicitly due to correlated pair hopping.
More generally we may remark, as we did after equation (2), that any model with one or more conserved momentum variables and two-particle interactions can be described in momentum space and then the dynamics necessarily shows symmetric two particle hopping. Since the momentum states usually do not have spatially varying density, this hopping is difficult to interpret in real space. Also, hopping may now occur for pairs of widely separated momenta. Many of the phenomena we describe here are therefore specific to FQH and other chiral systems. For example we would not expect "freezing" of wide bands of momentum states. Nevertheless some phenomena observed here may carry over. For example the asymptotic distribution in momentum space should depend on the conserved total initial momentum (or angular momentum). E.g. for systems on the sphere it would be interesting to see under what conditions a linear profile in angular momentum space emerges, like that observed in section 4.4.
We have reported a (to our knowledge, first) study of the level statistics of the LLL-projected Hamiltonians. Level statistics is considered important for dynamics, primarily because it is sensitive to integrability, and some classes of non-equilibrium phenomena are affected by proximity to integrability. For many-body systems, there is a stark difference between integrable and non-integrable (chaotic or generic) Hamiltonians with regards to the long-time dynamics -integrable systems do not obey the Eigenstate Thermalization Hypothesis, and hence do not thermalize in the same way as chaotic many-body systems [1][2][3]75]. In addition, integrability affects quantum transport, leading in some cases to ballistic behavior even at finite temperatures [151][152][153][154][155][156]. For inhomogeneous quenches, integrability allows one to formulate a generalized hydrodynamics [25,31,156,157]. The type of expansion dynamics we are considering is in spirit reminiscent of transport; hence the presence or absence of integrability is a relevant question in the context of the present work. Through our analysis of level statistics, we have demonstrated that the LLL-projected Hamiltonians are not integrable in either of the two geometries. However, considering the types of Wigner-Dyson statistics in various symmetry sectors, we have found that many sectors show GOE statistics despite the broken time reversal symmetry. These are many-body manifestations of the mechanism proposed in Ref. [51] -the relevant anti-unitary symmetry operation is not pure time reversal but a combination of time reversal and spatial reflection.
While this has been a detailed exploration, we believe we have only scratched the surface of fermionic LLL dynamics, and many questions remain open. The interplay of  Figure A1. Hopping elements for the V 1 potential and the Coulomb potential along the cut V L−2d,d .
the interaction range and system size might give rise to nontrivial size dependencies of dynamical behaviors, which we have not pursued. We expect the dynamical features to depend only loosely on the filling fraction, but this remains to be explicitly explored.
Our block initial states have been chosen in analogy to the domain-wall literaturea range of other initial states could be of interest and would lead to different classes of physics. While studying the level statistics, we have found in some cases GOE statistics in the absence of time reversal invariance. This calls for a more detailed analysis of symmetries in condensed-matter models violating time reversal symmetry.
the hopping term V i+d L−2d,d , as shown pictorially in Fig. A1(a). The assumption here is that hopping predominantly occurs to the first orbital outside of the block, at very short times, which is consistent with the observations of Fig. 4(a).
Inspecting the values of V i+d L−2d,d for a N = 10 particle system (black line) in Fig. A1(b), we see that the maximum is not at d = 1, the most short ranged allowed hopping, but rather at d = 3. This is consistent with the position of the non-monotonic dip in Fig. 4(a). For large N , the position of the maximum of V i+d L−2d,d is found to scale as d ∼ √ N . This means that the distance of the dips from the edges of the initial block should scale as ∼ √ N at large N . Considering instead the Coulomb interaction, Fig. A1(c), we see the maximum of V i+d L−2d,d occurs at d = 1 (d = 0 means no hopping) showing that this potential would not have the same non-monotonic behaviour at early times.