Dynamics of quasiholes and quasiparticles at the edges of small lattices

We study quench dynamics of bosonic fractional quantum Hall systems in small lattices with cylindrical boundary conditions and low particle density. The states studied have quasiholes or quasiparticles relative to the bosonic Laughlin state at half filling. Pinning potentials are placed at edge sites (or sites close to the edges) and are then turned off. Because the edges of fractional quantum Hall systems host chiral edge modes, we expect chiral dynamics, with motion in one direction for positive potentials pinning quasiholes, and motion in the other direction for negative potentials pinning quasiparticles. We numerically show that chiral motion of the density distribution is observed and robust for the case with positive potentials (quasiholes), but that there is no noticeable chiral motion for negative potentials (quasiparticles). The comparison of the numerical ground states with model lattice Laughlin wavefunctions suggests that both positive and negative potentials do create and pin anyons that are not necessarily well-separated on small lattices. Initializing the dynamics with the model state also shows the lack of chiral dynamics of quasiparticles. Our results suggest that, in small lattices with low particle density, quasiparticles are strongly adversely affected in dynamical processes, whereas quasiholes are dynamically robust.


I. INTRODUCTION
A fascinating aspect of topologically ordered phases of matter is that they support anyonic quasiparticles with fractional exchange statistics [1].Anyons appear, e.g., as charged excitations of fractional quantum Hall (FQH) systems in two-dimensional electron gases subject to strong magnetic fields [2].The robustness of exchange statistics against local noise makes anyons an interesting platform for topological quantum computation, which has motivated much work towards realizing anyons in different physical systems [3,4].Various models hosting FQH physics and anyons have been proposed for different physical platforms [5][6][7][8][9][10].Both the fractional charge [11][12][13][14] and the fractional statistics [15] of anyons have been observed experimentally.
It is a crucial question whether the characteristic features of FQH systems in two-dimensional electron gases still hold for small lattice systems with low particle densities, or what features should be chosen to characterize FQH states in such systems [31,32].Small lattice systems are ideal testing grounds from both experimental and computational viewpoints: Experimentally, the manipulation of single sites and single atoms [19,35] provides tools to observe the FQH effect in small lattices (if we know what to observe), and computationally, the exponentially growing size of the Hilbert space puts a limit on the system sizes that can be studied.Proposals for probing the fractional charge and fractional statistics in small optical lattices have been made [31,33,34,46].It is, however, an open question whether known results for solid-state systems with macroscopic numbers of electrons are still true for small (or even medium) size lattice systems with dozens of sites and few particles.Finite-size effects can influence the results in various ways, for example by making it impossible to separate the anyons sufficiently, or by altering the energy spectrum.Another issue is that, due to low site occupancies, the process of pinning or localizing quasiparticles could excite an uncontrolled amount of excitations.
One of the characteristic features of FQH states is the presence of chiral modes at the edges.This result was initially derived for large continuum systems (longwavelength limit) [36], but numerical calculations have shown chiral motion of quasiholes in small continuum systems [37,38] and chiral motion of "full" particles in (relatively large) lattice systems [39].In addition, even very small lattice systems retain some degree of similarity to the continuum, which is seen in the counting of the edge states [40].It is therefore relevant to ask whether the chiral motion of charges such as quasiparticles and quasiholes can be observed in small lattice systems, comparable to the ones used in experiments [28].
Here, we investigate whether quasiholes (quasiparticles) can be pinned by positive (negative) potentials in small lattices and whether chiral dynamics along the edge is present in the quench dynamics (i.e. after turning off the pinning potential).We consider hardcore bosonic FQH systems with filling factor ν = 1/2 and place pinning potentials at or near the edges of the lattice.The hardcore constraint is the U/J → ∞ limit, where U is the on-site interaction, and J is the hopping strength.In arXiv:2305.18364v2[cond-mat.mes-hall]9 Feb 2024 practice, when U/J is large enough, the hardcore constraint is a good approximation (for negative pinning potentials V 0 , we also need |V 0 | ≪ U ).This results in the creation of density depletions (increases), which we interpret as quasiholes (quasiparticles) due to the accumulated excess density being close to ±0.5 and due to high overlap with model wavefunctions.For the lattice sizes considered in this work, i.e. dozens of sites, we observe chiral motion of quasiholes.We also find that the chiral motion is robust in the sense that it is observed for a range of lattice sizes and for various strengths and locations of the potentials.On the contrary, we do not observe chiral dynamics for the quasiparticles.
We find that the ground states with potentials (both positive and negative) have high overlaps with model lattice Laughlin states with anyons, suggesting that in both cases we do create anyons.This suggests that the absence of chiral motion for negative potentials stems from complexities arising in the dynamics of quasiparticles, presumably involving a large number of highly excited states.Thus, we show that in these setups (small lattices and low particle densities) there are substantial differences between the dynamics of quasiparticles and the dynamics of quasiholes.
The paper is structured as follows.In Sec.II, we introduce the model.In Sec.III, we consider the case with positive potentials and demonstrate chiral motion of the density distribution in the quench dynamics.In Sec.IV, we study the case with negative pinning potentials and find no chiral motion of the density distribution.In Sec.V we compare the ground state of the Hamiltonian with potentials with model lattice Laughlin states with anyons, showing that the overlap is high.In Sec.VI we expand on the proposed explanation of our main observations and conclude the paper by pointing out implications of our results.Appendix A explains the model lattice wavefunctions used in Sec.V, and Appendix B shows how to relate them to continuum Laughlin wavefunctions.

II. MODEL
We consider an interacting Hofstadter model with hardcore bosons on a two-dimensional square lattice.The sites are at the positions ⃗ r i = a(x i , y i ), where a is the lattice constant, x i ∈ {1, 2, . . ., N x }, and y i ∈ {1, 2, . . ., N y }.The Hamiltonian takes the form where the sum is over all k and j for which the sites at ⃗ r k and ⃗ r j are nearest neighbors.Note that each pair of neighbors appears twice in the sum.The operator c k annihilates a boson at site ⃗ r k , ϕ 0 = h/e is the flux quantum, h is Planck's constant, e is the elementary charge, and ⃗ A is the vector potential corresponding to a uniform magnetic field B perpendicular to the plane of the lattice.The strength of the nearest neighbor hopping is set to J = 1 for simplicity.The interaction between bosons is included implicitly by enforcing the hardcore condition, i.e. allowing at most one boson per site.We consider a square lattice with open boundary conditions in the x direction and periodic boundary conditions in the y direction.This cylinder topology provides two parallel edges without corners, and hence is particularly suitable for studying chiral dynamics along the edges.We use the Landau gauge where x 0 = (N x + 1)/2, which ensures that the vector potential vanishes at the center of the system.Note that a 2 B = αϕ 0 , where α is the number of flux quanta per plaquette.
We consider half filling ν = 1/2 in this paper, where the filling factor ν is the number of particles divided by the number of flux quanta penetrating the lattice.(Note that ν = 1/2 is the filling fraction of the lowest Hofstadter band, and not the lattice filling fraction or average particle density, which is considerably smaller in our case.)If there are quasiholes (quasiparticles) present, they each count for ν particles (−ν particles) when computing ν.The number of flux quanta is α times the number of plaquettes N plaq in the lattice.Altogether, we hence have where M is the number of particles and p k = +1 (p k = −1) if the kth anyon is a quasihole (quaisparticle).At the edge appearing due to the open boundary conditions in the x direction, there are ambiguities in how to count the amount of flux through the lattice, and this ambiguity is significant for small lattices.It was suggested in [31] that the choice N plaq = (N x − 1)N y is best for creating and stabilizing a FQH droplet in small lattice systems, and we hence also use this way of counting the flux here.We compute α from (3) with ν = 1/2, and this gives us the vector potential ⃗ A appearing in the Hamiltonian.Since we are considering the ν = 1/2 Laughlin state, a quasihole corresponds to half a particle missing in a local region, and a quasiparticle corresponds to half a particle extra in a local region.We would like the anyons to be pinned at particular positions to begin with and local pinning potentials are a standard approach to do that.Adding a potential term V j n j to the Hamiltonian, where n j = c † j c j is the particle number operator, results in an energy penalty for a particle to sit at the jth site if V j > 0 and the opposite for V j < 0. A natural way to try to trap two quasiholes is hence to reduce the number of particles M by one and add positive potentials on two sites.Similarly, one could try to trap two quasiparticles by increasing the number of particles M by one and putting negative potentials on two sites.Note that M + ν k p k is unchanged in this process, and hence α is also unchanged.The Hamiltonian with potentials is where only two of the V j are nonzero.Below we shall investigate quench dynamics.We take the initial state at time zero to be the ground state |ψ(0)⟩ ≡ |E 0 ⟩ H of the Hamiltonian (4) with trapping potentials present.At time zero, we turn off the trapping potentials to start the quench dynamics, and after a time t, the state has evolved into |ψ(t)⟩ = exp (−iH 0 t/ℏ)|ψ(0)⟩. ( With the purpose of investigating properties of |ψ(t)⟩, we define the density distribution at time t as The first term ⟨ψ(t)|n j |ψ(t)⟩ on the right hand side is the expectation value of the number of particles on site j in the state |ψ(t)⟩, and the second term is the average number of particles per site when there are M + ν k p k particles in the system.Note that j ρ j (t) = −ν k p k , so j ρ j (t) = −1 for two quasiholes and j ρ j (t) = +1 for two quasiparticles at half filling.FQH systems tend to have uniform density in the bulk, and hence nonzero ρ j (t) can appear due to edge variations or density disturbances due to a potential.To eliminate the edge variations, we shall also sometimes consider the excess particle density ρj (t) = ⟨ψ(t)|n j |ψ(t)⟩ − ⟨n j ⟩ 0 , where ⟨n j ⟩ 0 is the expectation value of n j computed in the ground state of H 0 with M + ν k p k particles.Our computations involve exact diagonalization.We consider 3 particles for the case without anyons to ensure that the dimension of the Hilbert space is not too large.When trapping potentials are present, we have M + ν k p k = 3, i.e., M = 2 bosons for the case of positive potentials (two quasiholes), and M = 4 bosons for the case of negative potentials (two quasiparticles).We present results for the real-time dynamics of the lattice system with N = N x ×N y = 7×7 sites.We have numerically checked that the dynamics is qualitatively the same in other lattices of comparable size, e.g.N = 7 × 11, N = 11 × 9, etc.We concentrate on initial states for which the trapping potentials are located at particular sites.Two cases are considered: case (i), in which the two sites are placed on the middle of the left and right edges, i.e. the sites (1a, 4a) and (7a, 4a) for the 7 × 7 lattice; case (ii), in which the two sites are moved one site into the bulk from the left and right edges, i.e. (2a, 4a) and (6a, 4a) for the 7 × 7 lattice.V 0 between the ground state of the total Hamiltonian H with potentials and the lowest eigenstates of the Hamiltonian H0 without potentials.V0 is in units of hopping strength J.

III. QUASIHOLES
We first consider the case of quasiholes, with positive potentials used to pin (localize) them.We choose the following configuration: the number of sites is N = 7 × 7, the number of particles in the system with positive potentials is two, the number of fluxes per plaquette is α = 1/7, and positive potentials with strength V 0 are introduced on the sites (1a, 4a) and (7a, 4a).In contrast to the case of quasiparticles, which is shown and discussed in the following section, the chiral motion along the y-axis is clearly demonstrated in the case with positive potentials.Moreover, the observed chiral motion is quite robust in the sense that it does not depend on the strength of the potentials or their precise location on the edge.For example, we obtain similar results if we instead put the potentials on the sites (2a, 4a) and (6a, 4a).We also find that the chiral motion persists for a very long time.
We first investigate the effect of the potential strength V 0 on the system by considering the particle density at site i and the overlap between the ground state |ψ(0)⟩ of the total Hamiltonian H with potentials, and the nth eigenstate |E n ⟩ of the Hamiltonian H 0 with the same number of particles.If H 0 has several degenerate eigenstates |E n ⟩ m at a given energy E n , their total overlap is given by The squared overlap O 2 V0 measures the percentage that the eigenstate |E n ⟩ (or, in case of degeneracy, all the states at energy E n ) contributes during the quench dynamics given the initial state is |ψ(0)⟩.In Fig. 1(a) we plot the particle density ⟨n i ⟩ V0 at site (1a, 4a).Due to the symmetry under π rotations (inversion), the density at this site is equal to the density at site (7a, 4a).We observe that the particle density at site (1a, 4a) is quite small even when the potential is not strong, say, V 0 = 0.5.In addition, ⟨n i ⟩ V0 drops quickly as V 0 increases.This is expected, since the particle density on the site should go to zero when there is a large positive potential on the site.
In Fig. 1(b) we show the squared overlap between the initial state |ψ(0)⟩ and the zeroth, first, and second excited state of H 0 .The ground state |E 0 ⟩ of H 0 is nondegenerate, while the first excited state is two-fold degenerate |E 1 ⟩ 1,2 .This degeneracy is not lifted even with a very strong potential, say, V 0 = 10000.We can see that the overlap of the ground state dominates over all other overlaps of the excited states.We find that O 2 V0 > 0.9604 for the ground state for all V 0 > 0. We note that this means that there is a large overlap between states with pinned and unpinned quasiholes, as the ground state with two particles has two unpinned quasiholes.
In Fig. 2 we show the quench dynamics of the system with positive potentials.In Fig. 2(a) we observe that the density distribution moves chirally around the edge (the y-axis) in addition to the dispersion along the edge.Meanwhile, there is little density spreading into the bulk sites.This is different from the dynamics observed in the models studied in Refs.[39,41] where the particles spread into larger regions in the bulk on longer timescales.In our work, however, the density depletion hardly spreads into the bulk anymore after a short time t = 1.The leftmost and rightmost columns (edge columns) do not show much visible variation in density, so it is convenient to use the second column (column next to the edge) to track and visualize the density dynamics.In Fig. 2(b) we take snapshots of the density distribution of the second column of sites counted from the left, and concatenate snapshots at different instants to form a spacetime plot.This visualizes the chiral motion of the density distribution which is seen as a 'staircase' pattern in such a plot.The observed chiral dynamics suggests that we have a quasihole propagating along each of the two edges.
This picture is supported by monitoring the excess density along the edges.We observe that even during long periods of quench dynamics the sum of ( 7) over the sites along the two leftmost (rightmost) columns of sites equals 1/2 with a good accuracy, ρedge = i:xi≤2 ρi (t) ≈ −0.49.The total excess density corresponding to half a particle near each edge is consistent with the expectation that quasiholes are created and trapped.We note however that the system is inversion-symmetric, and there is one particle missing compared to the case without potentials, so a similar density pattern could in principle also arise in non-topological systems.The interpretation of the density depletions as quasiholes is further corroborated by the analysis of model wavefunctions in Sec.V.The nearly-constant value of ρedge over time suggests that the quasiholes stay confined to the edges as they evolve in time.
We find similar results as above when we examine the quench dynamics for various sizes of the lattice, different locations of the trapping potentials (at the edges or one site away from the edges) and various potential strengths V 0 .In fact, we find that ρedge ≈ −0.5 even if no potentials are present (V 0 = 0), although in such a case the excess density is spread uniformly in the y direction and hence chiral motion is not visible.But as soon as the potentials are strong enough to create a discernible "bulge" of the excess charge near their locations, the chiral motion can be observed.We verified that this is the case even for weak potentials, such as V 0 = 0.5.
Moreover, the chiral motion is also observed for values of the flux density, say, α = 6/49, which deviate slightly from the correct value α = 1/7.(In fact, if we replace N plaq with N in Eq. ( 3) we get the filling factor 1/2 for α = 6/49.The ambiguity in defining the magnetic flux through the lattice, when edges are present, is thus not a severe problem for the considered systems.)These results indicate that the chiral dynamics of quasiholes pinned by the positive potentials is quite robust in small lattices.In addition, we observe that the chiral motion of the quasiholes persists for a large timescale t > 700.If the flux density deviates much from the correct value α = 1/7, one would expect that the chiral motion is absent.In Fig. 3, we show the quench dynamics of the system with the same configurations in Fig. 2, but with the flux density α = 3/10.As expected, no chiral motion is observed.
Furthermore, we consider an extreme case in which the magnetic field is turned off, α = 0.The result is shown in Fig. 4. In this case, there is no chiral motion.Instead, we observe that the density distribution at the sites of the left (and right) half of the lattices split into two parts with equal weight, and two regions of depleted density move with opposite chirality (upward and downward directions).

IV. QUASIPARTICLES
In this section, we consider the case with four particles.We use negative potentials for creating the initial state.In this way we attempt to spatially localize quasiparticles.We find that the trapped objects have total excess particle density close to 0.5, as quasiparticles are expected to, but that they do not display clear chiral dynamics.V 0 between the ground state of the total Hamiltonian H with potentials and the lowest eigenstates of the Hamiltonian H0 without potentials.V0 is in units of hopping strength J.
In Fig. 5(a), we show the particle density at site (1a, 4a) as a function of the potential strength −V 0 .As −V 0 increases, ⟨n i ⟩ V0 grows and eventually converges to its maximum possible value, which is 1.This is expected because the stronger the on-site negative potential is, the more the energy is reduced by having a particle on that site.We avoid potentials close to zero as |ψ(0)⟩ has a degeneracy for V 0 = 0.In Fig. 5(b) we show the squared overlap O 2 V0 between the initial state |ψ(0)⟩ and the lowest eigenstates of H 0 with four particles as a function of |V 0 | = −V 0 .The ground states of H 0 are doubly degenerate, and thus the overlap is calculated using (10).As −V 0 increases, more and more higher excited states have appreciable overlap and thus contribute to the quench dynamics.We observe that the overlaps with respect to |E n ⟩ depend substantially on |V 0 |, and many of them have comparable magnitude.The total overlap with respect to the doubly degenerate ground state |E 0 ⟩ 1,2 decreases rapidly with −V 0 ; its magnitude does not dominate over other eigenstates for −V 0 ≳ 1.5.This behavior is very different from the case of quasiholes trapped by positive potentials, shown in Fig. 1(b).In that case the ground state dominates even for large potentials.
In Fig. 6 we demonstrate the quench dynamics for negative potentials located at the edges with strength V 0 = −100 which is huge compared to the hopping strength.In Fig. 6(a) we show the dynamics of the density distribution, and in Fig. 6(b) we plot the density distribution of the first column of sites, counted from the left, as a func-tion of time.The particle density (7) in the initial state has an excess density of ρedge (0) = i:xi≤2 ρi (0) ≈ 0.51 in the two leftmost (two rightmost) columns of sites.That is, at the beginning of the dynamics we have an excess density of approximately a half of a particle (this time with opposite sign) located near each edge, consistent with the expectation that we trapped quasiparticles.Again we note that the system has inversion symmetry, and there is one particle extra compared to the case without potentials, so a similar excess density could also appear without the presence of topology.Further support that quasiparticles are trapped is presented in Sec.V.
A significant difference from the results in Sec.III appears when the density evolves in time.From Figs. 6(a) and 6(b), we observe no chiral motion of the density distribution, but patterns that are different from the case of quasiholes.We find that, during quench dynamics the excess density roughly splits into two parts and each of them move in opposite directions.Meanwhile, in contrast to the quasihole case, the excess density spreads into the bulk with much weight in the beginning, but stops spreading more into the bulk sites after some time.
In addition, we have varied V 0 and also varied the positions of the trapping potentials, e.g., at the edges, or one or two sites away from the edge.We have not observed clear chiral motion in any of these cases.Furthermore, we find that the patterns of quench dynamics for negative potentials vary substantially as we change the strength and/or locations of the pinning potentials, which was not the case for quasiholes trapped by positive potentials.We also tried to use pinning potentials over a block of sites to create quasiparticles as in [31].Also in that case, we have not seen clear chiral motion.

V. MODEL QUASIPARTICLE STATE
Could the difference between the dynamics for positive and negative potentials arise because the former succeeds in pinning anyons, and the latter does not?To further study the question whether the trapped objects are quasiholes and quasiparticles, we here discuss a model wavefunction with anyonic excitations [42][43][44].
We write the position of site j as a complex number ξ j = x j + iy j .Then, we define a mapping from a cylinder to a complex plane: z j = exp(2πξ j /L y ), where L y is the circumference of the cylinder (for a square lattice with unit lattice constant we have L y = N y ).For brevity, we write all z j s as a vector z = [z 1 , z 2 , . . ., z N ].Like any state of a hardcore boson system, the model lattice Laughlin state |ψ model ⟩ can be written in the occupation number basis, where n = [n 1 , n 2 , . . ., n N ] is the vector of the occupation numbers of each site, with n j ∈ {0, 1}, |n⟩ is the corresponding basis state, ψ model (z, n) are the unnormalized wavefunction coefficients and C is the normalization constant.For a model lattice Laughlin state without anyons on a cylinder [43], the coefficients are Here, δ(qM + p + + p − − N η) is the Kronecker delta fixing the charge neutrality, i.e. the number of particles M (with M = j n j ) and γ j = (qn j − η)/ √ q, with η being the flux assigned to each site (in our case, we put η = α).The integer q determines the topological order of the wavefunction and in our case q = 2.The gauge factors χ j are arbitrary complex numbers with |χ j | = 1.We note that in this work we use the word "gauge" in two meanings: the Landau gauge (2) of the hoppings of the Hofstadter model ( 1), and the factors χ j of the model wavefunctions.To avoid confusion, we will always refer to the latter as "the gauge factors χ j ".The real numbers p − , p + , have to fulfill the relation We note that the formulation in [42,43] of the lattice Laughlin states on a cylinder did not contain p − and p + .The origin and meaning of these parameters is described in Appendix A, while the relation of ( 12) to the continuum states is described in Appendix B.
In our calculations, we want to use a state with a given M , N and η = α.Thus, we determine p − from a combination of the charge neutrality condition and (13), and then p + from (13).Some properties of ( 12), such as the particle density and the topological order, are independent from the gauge factors χ j .Its overlap with the ground state of the Hofstadter model, however, strongly depends on χ j .Therefore, we optimize these factors numerically to maximize the overlap.For the 7 × 7 system with M = 3 and α = 1/7, we obtain the squared overlap |⟨ψ(0)|ψ model ⟩| 2 ≈ 0.984 after optimization.
We can add anyonic excitations to the state (12).They are characterized by the integers p k (see Eq. ( 3)), with p k = 1 and p k = −1 corresponding to a basic quasihole and quasiparticle, respectively, and the complex positions ω k ∈ C. Note that the complex positions are not required to coincide with lattice sites [44].After mapping to the complex plane, these locations transform into w k = exp(2πω k /L y ).The wavefunction with anyons is FIG. 7. The locations and strengths of the potentials that we optimize to increase the overlap of the ground state of the Hamiltonian with negative V0 with the model wavefunction (15) with two quasiparticles.During the optimization, V0 = −100 stays constant, while VA and VB are adjusted.V0 is in units of hopping strength J. then The w k have the braiding properties expected for anyons in systems with Laughlin type topology when the w k remain well-separated during the entire braiding process [44].This is true even if p k is negative -while such a wavefunction does not have a consistent continuum limit [45], it is a valid quasiparticle ansatz on lattices with finite lattice constant, yielding correct quasiparticle charge and statistics [44].On small lattices with low densities, the changes in the particle density caused by the anyons may not be fully separated.
Therefore, the ground states of the Hofstadter model states are similar to the model states, but not perfectly equal to them.To check whether the discrepancy is the reason for the lack of chiral dynamics in the quasiparticle case, we also compute the quench dynamics generated by the hopping Hamiltonian H 0 (Eq.( 1)) with the model state (15) with two quasiparticles as the initial state.The results are quite similar to Fig. 6, i.e. we do not observe clear chiral dynamics.
We note that we can increase the overlap between the quasiparticle state of the Hofstadter model and the model lattice Laughlin quasiparticle state by introducing auxiliary potentials V A and V B at the nearest-neighboring sites of the original potentials (see Fig. 7), and optimize them to maximize the overlap.In such a way, we obtained the squared overlap |⟨ψ(0)|ψ model ⟩| 2 = 0.993.However, using the ground state with the optimized potentials as the initial state also does not lead to chiral dynamics.
The results presented in this section suggest that the difference between the cases of positive and negative potentials is not related to the ability to pin anyons, as both quasiholes and quasiparticles seem to be successfully created and pinned by the pinning potentials.

VI. DISCUSSION AND CONCLUSIONS
We have studied quench dynamics in small lattices with cylindrical boundary conditions, using interacting (hardcore) bosons subjected to the Hofstadter Hamiltonian.We considered particle numbers appropriate to both two quasiholes (M = 2) and two quasiparticles (M = 4), relative to a FQH state with filling factor ν = 1/2.In the two cases we used positive and negative pinning potentials so that the initial positions of the quasiholes and quasiparticles are localized near the edges of the lattices.The quench dynamics begins when the trapping potentials are turned off.One expects to observe chiral motion of the excess density due to the presence of chiral edge modes.Our results for the two cases are very different.While it is generally easy to observe chiral motion of quasiholes (negative excess density), and the observed chiral dynamics are quite robust, we have not found clear chiral motion of quasiparticles (positive excess density).
The absence of the chiral motion of quasiparticles was further verified by investigating the quench dynamics with the model lattice Laughlin state as the initial state.Moreover we introduced auxiliary potentials in the Hofstadter model to optimize the overlap between the ground state of the Hofstadter Hamiltonian and the model lattice Laughlin state.The chiral dynamics was still absent.
Why is there such a pronounced difference between the quasihole and quasiparticle cases?Due to the small particle densities, the details of pinning are quite different in the two cases-locally decreasing the density is a very different process compared to locally increasing the density.In the former case (quasiholes, negative excess density) there is a natural limit to how far the local density ⟨n i ⟩ V0 can be reduced, since ⟨n i ⟩ V0 is small already for V 0 = 0.In the latter case (quasiparticles, positive excess density), the site density can increase enormously, all the way up to 1, leading to the excitation of many high-energy modes whose nature might be unrelated to FQH physics.Notice that this is true with the hardcore constraint U/J → ∞.These high-energy contributions can be expected to mask or disturb the topological and chiral physics that FQH quasiparticles should show.This difference between the two cases is mirrored in the overlaps shown in Figs.1(b) and 5(b).This difference raises the possibility that the pinning of quasiparticles might be disrupted, unlike the pinning of quasiholes.However, the comparison with model quasiparticle states (Sec.V) indicates that both positive and negative potentials pin the anyons successfully.We therefore conclude that the difference must come from the contributions of excited states in the dynamics -the quasihole dynamics is dominated by the ground state, while the quasiparticle dynamics involves a number of highly excited states, as indicated in Figs.1(b) and 5(b).The disruption of chiral dynamics is a dynamical effect.
The absence of chiral motion suggests that the dynamics of quasiparticles in small lattices with low particle densities, of the sizes considered in this work, is different from that in large lattices, and from that in twodimensional electron gases.This has important implications for the study of FQH physics in novel platforms such as that of Ref. [28], where the particle density and system size might be naturally small.Our results indicate that, in such situations, quasiparticle dynamics is fragile as there is a natural tendency to create substantial additional (possibly unwanted) excitations when localized quasiparticles are created.Remarkably, the difficulty with observing quasiparticle features shows up primarily in the dynamics, and not necessarily in creating or pinning them.
This work opens up a number of questions: (1) Our conjectured explanation above suggests that quasiparticles and quasiholes would behave more similarly if the average particle density (lattice filling) were close to 1/2.It would be worthwhile to check this explicitly; however in our setup this is computationally too challenging.Systems with larger numbers of sites and particles can be studied using tensor network methods [47].(2) A detailed understanding of how the dynamical behaviors change with system size (e.g., with constant density or with constant particle number) is currently lacking.An intriguing question is whether chiral motion of quasiparticles become readily observable for large lattice sizes even when the particle density remains small.(3) It would also be interesting to investigate how the presence or absence of inversion symmetry affects the excess density and dynamics in small systems.(4) We here considered the case of Abelian anyons.In future work it would be interesting to study the quench dynamics of FQH systems which host non-Abelian anyons.
Appendix A: Inversion-symmetric lattice Laughlin wavefunctions on a cylinder In this Appendix, we derive the expression (12) and condition (13) from the known expressions for lattice Laughlin wavefunctions [42,43] by demanding inversion symmetry.
Let us start by considering a planar system of N sites at locations described by complex numbers z j = x j + iy j , with j = 1, 2, . . .N , populated by M particles: either by hardcore bosons or fermions.To each site, we assign a flux η ∈ R in units of the flux quantum, so the total number of flux quanta in the system is N η.We also place a compensating charge of η/q (in units where a single particle has charge −1) at each site, in analogy to a continuum Landau level where a uniform background charge is used.
According to [42,43], the coefficients of a lattice Laughlin wavefunction on a plane (see Eq. ( 11)) are (see below Eq. ( 12) for an explanation of the notation).
One can add one or more anyonic excitations to (A1), [43,44] which transforms it into where, similarly to (15), p k ∈ Z and p k /q is the charge of the kth anyon.The coordinate w k ∈ C is its position.Both p k and w k are external parameters of the wavefunctions.Note that w k is not required to be equal to one of the z j , but can take any value in the complex plane.Equations (A1) and (A2) can be generalized to a cylinder [42,43].Let us consider a plane of size L x ×L y , whose upper and lower edges are glued to make a cylinder (i.e.y is the periodic direction, and L y is the circumference of the cylinder).If ξ j and ω i are the coordinates of sites and anyons, respectively, within the L x × L y rectangle, then the substitution z j = exp(2πξ j /L y ), w j = exp(2πω j /L y ) turns (A1) and (A2) into cylinder wavefunctions [42].
It was shown that these cylinder wavefunctions have Laughlin-type topological order [42,43].However, while the particle density in the continuum cylinder wavefunction with no anyons seems to be invariant under inversion [48], in the lattice cylinder wavefunction it is not (except from the special case η = q/2).An example for a 5 × 5 square lattice with q = 2 and M = 2 is shown in Fig. 8  The same quantity for the modified wavefunction (12) with p− given by ( 14), and q, M and η the same as in (a).
The lack of inversion symmetry is alarming for two reasons.First, in real quantum Hall systems exchanging the left and right side of the cylinder should not matter.Second, on the plane, the lattice wavefunctions approach the continuum wavefunction in the limit of η → 0 (i.e.infinite number of sites per flux quantum) [42].Similar arguments should apply for a cylinder [43], and if this is indeed the case, we should explain why the inversion invariance does not exist for general η and is restored at η = 0.
In the following, we will show that the inversion invariance of (A1) on a cylinder can be restored at any η by choosing the numbering of the sites appropriately and adding appropriate charges, described by the numbers p − , p + ∈ R, at the coordinates w − = 0 and w + = ∞, respectively.They should be understood as infinitely thin, immobile flux tubes.If we plot the complex numbers z j on the plane, they form rings of different radii.The flux tube at (w − = 0) is at the center of the rings, while the flux tube at (w + = ∞) is located outside the rings, infinitely far away.As a result, the effects of the charge p + appears only in the charge neutrality condition [43].The wavefunction is thus given by ( 12) (but so far, we have not determined the values of p − and p + ).
In our considerations, we will permute and modify the z j and n j .Therefore, we will write the individual arguments explicitly, i.e.
Inversion of the coordinates ξ j around the point ξ c transforms ξ j − ξ c → −(ξ j − ξ c ) and hence where c = e 4πξc/Ly .For now, we will consider a general lattice, which is not necessarily inversion invariant.As √ πl ′ B described in the main text.The many-body wavefunction, with x-dependent particle density shown schematically in (a), is composed out of three single-particle orbitals with m = 0, 1, 2 -the blue ones in (b).The cylinder is drawn schematically in (c), with blue and grey lines denoting the Gaussian centers of the orbitals.The blue region of the cylinder located between the centers of the 0th and 2nd orbitals contains q(M − 1) = 2 flux quanta.The dots ". . ." in (b) and (c) mean that the cylinder extends infinitely in both positive and negative x directions.The blue and green regions contain a total number of 6 flux quanta.In (d) we show the blue and green regions unwound to a plane.The circumference of the cylinder is chosen in such a way that the blue region is a square.In (e) and (f) we show the site positions for the analogous lattice systems with total flux N η = 2 and N η = 6, respectively, for an Nx × Ny = 10 × 10 lattice.
The region between the Gaussian centers of the nearest-neighbor orbitals always contains one flux quantum.Therefore, the x distance between the Gaussian centers depends on L ′ y .As the circumference of the cylinder decreases, the orbitals get more separated, and the nearly-uniform quantum Hall state transforms into a charge density wave [48] (in the example in Fig. 9(a), the inhomogenity of the density is already quite developed, with high density on the m = 0, 2 orbitals and low density on the m = 1 orbital).
We also note that (B1) is in fact defined on an infinite cylinder.However, we can also restrict (B1) to a finite cylinder by truncating the domain of Ξ ′ j to some finite range and normalizing the wavefunction appropriately.One possible choice is to set the cylinder borders to the centers of m = 0 and m = q(M − 1) orbitals -in such a way the system is pierced by q(M − 1) flux quanta, in accordance with (A14).In our M = 2 example, this is the blue region in Fig. 9(c), also shown in Fig. 9(d) after unwinding to a plane.It contains two flux quanta.However, we can also use a bigger system, which would allow us to observe the Gaussian "tails" at the edges.For example, we can study a system thrice as large (six flux quanta), denoted by combined green and blue regions in Fig. 9(c), (d).Since we focus on the systems where the particle density is inversion-invariant, we consider only cases where the lengths added before the center of m = 0 orbital and after the center of m = q(M − 1) orbitals are the same (i.e. that the two green regions in Fig. 9(c), (d) have the same size).We do not restrict the total flux to be an integer, i.e. the edges of the domain do not have to lie at the Gaussian centers as in Fig. 9. Now, knowing the meaning of various parameters of the continuum wavefunction, we propose that to match the particle densities of continuum and lattice wavefunctions, we have to do the following: 1. obtain the magnetic length l B of the lattice system to establish a common length scale, 2. match the circumferences L ′ y /l ′ B = L y /l B (i.e. the distance between the Gaussian centers), 3. match the total flux through the system (i.e. the extent of the system in the x direction), 4. match the center of symmetry (i.e.shift the x coordinates so that the center of both systems is at x = 0), not necessarily in this order (it depends on which parameters of which system we want to vary).
As an example, let us consider finding lattice systems corresponding to the M = 2 continuum wavefunction from Fig. 9 defined on the blue region (two flux quanta).We will consider a rectangular lattice of a fixed size N x × N y in terms of unit cells, with a unit cell of size a × 1, where a will be adjusted in the process.We have N = N x N y and L y = N y .Thus, by doing step 2 (i.e.demanding L ′ y /l ′ B = L y /l B ), and recalling that for the system from Fig. 9 we had L ′ y = 2 √ πl ′ B , we find l B = N y /(2 √ π).From step 3 we get that N η = 2, and hence η = 2/N .Now, let us look again at the magnetic length.Eq. (12) does not refer to it explicitly, but we can infer the relation of l B to the unit cell size by considering a loop and comparing its area to the encircled flux (for simplicity, we treat the flux as uniformly distributed in space).A unit cell of size a × 1 corresponds to η flux quanta, so we must have a = 2πηl 2 B .In this way, we determined the shape of the lattice.
Next, we find p − from (14). Figure 10(a) shows the comparison of the particle densities for continuum and lattice for various N x × N y .In the plot, the x coordinate is plotted in units of magnetic length of the respective system and shifted so that the center of the system is at x = 0.The particle density is normalized so that its integral is equal to M (in the case of lattice systems, we plot ρ(x) = ⟨n(x)⟩l B /a).The densities are calculated exactly for sufficiently small systems and using Monte Carlo for the largest ones.It can be seen that there is a good match between the lattice and continuum wavefunctions, especially in the center.Notable mismatches, especially near the edges, exist only for the cases with small N y .As N y approaches N x , the mismatch vanishes almost completely.
Similarly, we can repeat the procedure for the combined blue and green regions of Fig. 9 (six flux quanta), with the result shown in Fig. 10(b).Again we obtain a good match, with small but visible mismatches for small N x .While in the case of two flux quanta and small N y the density near the edges significantly departed from the continuum values, here we do not observe such an effect, probably because the system is more elongated, or because the density at the edges is very small.If one considers a Hofstadter model (1) with the gauge (2), then one can compare the particle densities in the ground state and the model wavefunction in an analogous way.Matching the wavefunctions to compute the overlap is also possible, but more complicated.The model wavefunction has to be centered at the m = 0 orbital, which can be achieved by a shift of the X ′ j coordinates together with a modification of the phase of (B1) (the phase transformation does not affect the particle density).

FIG. 1 .
FIG. 1. Effects of the pinning potential in the quasihole case with two particles.(a) Particle density ⟨ni⟩V 0 for the site (1a, 4a) (at the edge) versus potential strength V0.(b) Squared overlap O 2V 0 between the ground state of the total Hamiltonian H with potentials and the lowest eigenstates of the Hamiltonian H0 without potentials.V0 is in units of hopping strength J.

FIG. 2 .
FIG. 2. (a) Quench dynamics of the density distribution (6) for the quasihole case.The size of the lattice is 7 × 7, and the number of particles is two.The trapping potentials are located at the sites (1a, 4a) and (7a, 4a), i.e. at the left and right edges, and have strength V0 = 100.The flux density is α = 1/7.(b) The density distribution of the second column of sites, counted from the left, as a function of time for the quench dynamics shown in (a).A staircase pattern is seen, which means that there is chiral motion.Time is in units of inverse of hopping strength, i.e., ℏ/J.
FIG. 3. (a) Snapshots of the time-evolving density distribution (6) for the system with quasiholes and a flux density α = 3/10, significantly different from the 'correct' value α = 1/7.The trapping potentials (V0 = 100) are at the edges, (1a, 4a) and (7a, 4a).(b) The density distribution of the second column of sites, counted from the left, as a function of time for the quench dynamics shown in (a).The lack of a clear staircase pattern signals the absence of chiral motion.Time is in units of inverse of hopping strength, i.e., ℏ/J.

FIG. 4 .FIG. 5 .
FIG. 4. (a) Snapshots of the density distribution (6) for the quasihole case in the absence of the magnetic field, i.e. α = 0.The trapping potentials (V0 = 100) are located at the edge sites (1a, 4a) and (7a, 4a).(b) The density distribution of the second column of sites, counted from the left, as a function of time for the quench dynamics shown in (a).When the magnetic field is turned off, the model has mirror symmetry both between up and down and between left and right, and hence there is no chiral motion.Time is in units of inverse of hopping strength, i.e., ℏ/J.

9 FIG. 6 .
FIG. 6.(a) Snapshots of the density distribution (6) in the quasiparticle case.The trapping potentials (strength V0 = −100) are at the edge sites (1a, 4a) and (7a, 4a).A different color scale has been used in this plot, because the particles spread into the bulk sites on a longer timescale.(b) The density distribution of the leftmost column of sites as a function of time for the quench dynamics shown in (a).No staircase pattern is seen, which signals the absence of chiral motion.Time is in units of inverse of hopping strength, i.e., ℏ/J.
(a).The figure shows the average particle number at a given x coordinate, defined as n(x) = j δ(x − x j )⟨ψ model |n j |ψ model ⟩.

FIG. 8 .
FIG. 8. (a)The average particle number as a function of the (non-periodic) coordinate x for a 5 × 5 square lattice, for the "naive" implementation of a q = 2, M = 2 lattice Laughlin wavefunction on a cylinder, i.e.Eq. (A1) with zj = exp(2πξj/Ly) and η = 4/25.(b) The same quantity for the modified wavefunction(12) with p− given by (14), and q, M and η the same as in (a).

FIG. 9 .
FIG.9.The relation between continuum and lattice systems for the example with q = 2, M = 2, and L ′ y = 2