Charge and statistics of lattice quasiholes from density measurements: a Tree Tensor Network study

E. Macaluso, T. Comparin, 2 R. O. Umucalılar, M. Gerster, S. Montangero, M. Rizzi, 7 and I. Carusotto INO-CNR BEC Center and Dipartimento di Fisica, Università di Trento, I-38123 Trento, Italy Univ Lyon, ENS de Lyon, Univ Claude Bernard, CNRS, Laboratoire de Physique, F-69342 Lyon, France Department of Physics, Mimar Sinan Fine Arts University, 34380 Sisli, Istanbul, Turkey Institute for Complex Quantum Systems and Center for Integrated Quantum Science and Technologies, Universität Ulm, D-89069 Ulm, Germany Dipartimento di Fisica e Astronomia “G. Galilei”, Università degli Studi di Padova & INFN, I-35131 Padova, Italy Institute of Quantum Control (PGI-8), Forschungszentrum Jülich, D-52425 Jülich, Germany Institute for Theoretical Physics, University of Cologne, D-50937 Köln, Germany (Dated: October 14, 2019)


I. INTRODUCTION
In three spatial dimensions, quantum particles are typically classified into bosons and fermions, according to the symmetry properties of their many-body wave functions. In particular, bosonic (fermionic) many-body wave functions must be globally symmetric (anti-symmetric) in the particle coordinates, meaning that they take an overall +1 (−1) factor upon particle exchange. This classification is enriched in two dimensions (2D), where exotic particles called anyons have been predicted to exist [1][2][3][4][5][6] so that the effect of particle exchange (resp. braiding) on the many-body wave functions is a generic phase factor e iϕst (resp. e iϕ br = e i2ϕst ), where the statistical phase ϕ st can take any value in [0, 2π). While bosons and fermions are characterized by ϕ st = 0 and ϕ st = π, Abelian anyons have statistical phase ϕ st = απ, with α a non-integer number. In the presence of topologically degenerate ground states, the statistical phase factor is further generalized to non-commuting unitary transformations acting on the ground-state manifold, and anyons are said to be non-Abelian [7][8][9][10][11].
Among the physical systems for which the existence of anyons has been predicted, fractional quantum Hall (FQH) systems are probably the most popular ones [7,10,[12][13][14]. Such strongly correlated quantum fluids can host bulk elementary excitations -called quasiholes (QHs) and quasiparticles (QPs)-which have been theo-rized to carry fractional charge and exhibit anyonic behavior. Although the QH/QP fractional charge was measured in electron experiments [15], the anyonic statistics of these excitations still lack a clear-cut experimental evidence. For this reason, a large ongoing effort is based on controllable analog systems, where magnetic quantum-mechanical effects occur for neutrally charged particles such as atoms and photons prepared in the FQH regime [16,17].
In this context, several theoretical studies investigated the adiabatic preparation of different FCI states and the associated phase diagrams [35][36][37][38][39][40]. Some works focused on the numerical characterization of these states by inspecting key quantities such as the many-body Chern number, the particle entanglement spectrum, the behavior of the correlation functions and the topological entanglement entropy [41][42][43], while others proposed ex-perimentally applicable schemes to identify these elusive strongly correlated phases of matter [44,45]. Finally, growing attention has been given to FCI bulk excitations [46][47][48][49], which (similarly to those characterizing the FQH effect) display fractional charge and anyonic statistics.
In the presence of interactions, finding the exact ground state of a quantum many body problem is often possible only for small system sizes, by means of exact diagonalization (ED) of the Hamiltonian. A more recent class of numerical techniques is based on the optimization of quantum states in the family of Tensor Network (TN) states, which is a generalization of Matrix Product States -see for instance Refs. [50][51][52][53] for general reviews. Tensor Network states, which are distinguished by their different network topologies, notably include Projected Entangled-Pair States (PEPS) [54,55], Entangled-Plaquette States (EPS) [56], the Multi-scale Entanglement Renormalization Ansatz (MERA) [57,58], and the Tree Tensor Network (TTN) states used in this work. States in the Tensor Network family have already proved their usefulness for treating FCIs in several works, see e.g. in Refs. [39,[42][43][44].
In this work we make use of the TTN ansatz used in Ref. [42], to study the properties (the charge and the statistics) of the QH excitations appearing in the HH model with hardcore bosons. The key observable of our analysis are the depletions created by the QHs in the density profile of the system. While previous works inspected the QH density depletions mainly to quantify the charge of the QHs and/or to give an estimate of their size [see e.g. in Refs. [46,49]], only recently it has been discovered that, for FQH states in the lowest Landau level (LLL), these depletions encode information also on the QH anyonic statistics [59]. By considering the interacting HH model as a concrete example, here we provide numerical evidence that the experimental protocol proposed in Ref. [59] to extract the QH anyonic statistics from density measurements can be generalized also to FCIs in lattice geometries. Moreover, in order to provide a description which is as close as possible to experimentally realistic situations, we mainly focused our attention on systems with open boundary conditions (OBC) and mesoscopic numbers of particles. This makes our proposal to measure the anyonic statistics of FCI QHs readily applicable in state-of-the-art experiments with ultracold atoms and superconducting qubits, in which the HH model has already been implemented [60][61][62][63], and the occupation of the different lattice sites can be easily measured [62,[64][65][66][67].
The structure of the article is the following: In Sec. II we introduce the HH Hamiltonian for bosons with on-site interactions, and we briefly review its basic properties together with the main features of the TTN technique. Then, in the same section, we describe the Monte Carlo sampling of discretized Laughlin wave functions, which we use as an auxiliary method to interpret the results obtained through the TTN ansatz. In Sec. III we review the main findings of Refs. [59,68] and, in particular, the relation connecting the QH braiding phase to the depletions that these excitations create in the density profile of continuum systems. This relation is then generalized to the case of lattice systems in Sec. III B. The main results of our work are reported in Sec. IV. First, in Sec. IV A, we discuss how to use localized pinning potentials to selectively stabilize states presenting either a single QH or two QHs on top of each other. Then, in Sec. IV B, we provide numerical evidence that the method discussed in Sec. III is indeed able to reconstruct the QH anyonic statistics even in lattice systems with a relatively small number of particles. Conclusions are finally drawn in Sec. V.

A. Interacting Harper-Hofstadter model and Tree
Tensor Networks The Harper-Hofstadter (HH) model describes noninteracting spinless particles hopping on a twodimensional square lattice and subjected to an external magnetic field [24,25]. In the presence of on-site interactions and local energy offsets, the interacting Hamiltonian reads whereâ i (â † i ) is the bosonic annihilation (creation) operator at the site i of an L × L square lattice, and wherê n i =â † iâ i is the on-site particle density. The hopping between neighboring sites i and j is characterized by an amplitude t and by the Peierls phase θ ij , which is related to the magnetic flux passing through the system [69]. We focus on the case of a uniform magnetic flux per plaquette φ = αΦ 0 , where Φ 0 = 2π /e is the magnetic flux quantum. For a given value of the flux density α, multiple choices of the phases θ ij are possible. Among them, we choose the Landau gauge considered in Ref. [42]. Moreover, we consider the limit of hard-core bosons (that is, with infinite on-site repulsion U ), where each site can host at most one particle. The local energy offsets V i represent additional attractive or repulsive potentials on some lattice sites [see Fig. 1], which can encode the localized impurity potentials to pin the QHs and/or an additional trapping potential for the particles. Even though in the following we focus only on harmonic confinements, more complicated forms could be studied with no additional difficulties.
Our study of this model is based on the Tree Tensor Network technique. This is a variational ansatz in the class of Tensor Network states, with some noticeable advantages and shortcomings compared to other TN ansatz states, all essentially related to their loop-free structure. On the one hand, TTNs do not capture the entanglement area law for arbitrarily large systems (which, instead, PEPS would do by construction). On the other hand, their simple connectivity allows for very efficient computational algorithms [52,53,70]. Thus, the bond dimension D can be pushed to large enough values to compensate for the intrinsic weaknesses of the ansatz, yielding reliable numerical results for system sizes which are beyond reach with exact diagonalization.
As shown in Ref. [42], the TTN technique is capable to fully reproduce the properties of a FCI state without QH excitations, in the interacting Harper-Hofstadter model. These properties include the correct topological degeneracy when the system is placed on a torus, the many-body Chern number, the behavior of correlation functions, and the entanglement-entropy scaling laws. For the case of a FCI state with QH excitations, we benchmark the validity of the TTN ansatz through a comparison to exact results -see Appendix A.
The core results of this work are obtained for systems on a 16 × 16 square lattice, with open boundary conditions (OBC) and in the presence of an additional harmonic trapping. In the absence of trapping potentials, the number of magnetic fluxes would have to be commensurate with the number of particles in order to obtain an appropriate magnetic filling ν = N/N Φ , where the total number of fluxes for systems with OBCs is N Φ = (L − 1) 2 α. However, adding a harmonic confinement relaxes this constraint, and for a given number of particles N it is possible to obtain a FCI state as the ground state for different values of the flux density α [71]. Note that different values of α correspond to different ratios between the lattice constant a and the effective (model dependent) magnetic length l B ≡ A/2π, where A is the area of the magnetic unit cell [46]. In particular, for the model under study (i.e., the HH model on a square lattice), A = a 2 /α and such a ratio reads This expression of the lattice spacing in terms of the effective magnetic length will be of crucial importance both in Sec. II B and for the analysis of the QH density depletions in Sec. IV B. A final remark concerns the competition between the different possible states (FCIs, superfluid states, Mott insulators, charge-density-waves, etc.) which might be the ground state of the HH Hamiltonian in Eq. (1) -see for instance Refs. [72][73][74]. In this respect, superimposing an additional harmonic confinement to the lattice further enriches an already complicated scenario. However, for all the values of α we consider in the following it is possible to identify a finite interval of strengths of the harmonic confinement in which the system ground state is a FCI state. Although a detailed analysis of the different competing phases as a function of the confining strength is an interesting task, we postpone it to a future work and to focus here only on intermediate values of the harmonic potential strength for which the system ground state has a FCI nature.

B. Monte Carlo sampling of discretized
Laughlin-like wave functions As we will see, the TTN approach allows us to capture a very close approximation of the true ground state of the Harper-Hofstadter model for sizes compatible with nearfuture experiments in cold gases or circuit QED systems. However, there is a practical limit on the lattice size and particle number that can be treated.
To address more technical questions related to the effect of the lattice grid, rather than of the system size [see Sec. III B], we employ an auxiliary method which provides more flexibility. This consists in the Monte Carlo sampling of a discretized version of the Laughlin-like wave functions [13], which are evaluated on the sites of a two-dimensional grid [as done for instance in Ref. [20]]. More explicitly, at magnetic filling 1/2, the wave function for N particles and k QHs reads where η µ is the position of the µ-th QH, Ψ L (z 1 , . . . , z N ) is the celebrated Laughlin wave function, i.e.
and, on a lattice, the position of the j-th particle (z j = x j +iy j , with the complex coordinate notation) only takes discrete values (that is, with x j /a and y j /a integer). The use of the Laughlin wave function as a reference state for Harper-Hofstadter systems in the limit of low flux density α is common in literature [18,20]. In this work, we use the wave functions in Eq. (3) to study the properties of QH excitations, and to provide a comparison with the HH model. The usefulness of this comparison relies on the study in Ref. [46], which proves that the density profile around FCI QHs is a discretized version of the one of Laughlin QHs in continuum space, up to a proper rescaling of the length units. Of course, for a given flux density α, a proper comparison with the HH model requires that the lattice spacing of the discretization grid satisfies Eq. (2).
For a discretized Laughlin-like wave function, it is straightforward to generalize the Monte Carlo technique which is typically used to extract observable quantities for the continuum-space Laughlin state [75,76]. This consists in sampling configurations {z 1 , . . . , z N } distributed as in |Ψ(z 1 , . . . , z N )| 2 , which give access to observables like the density profile on the lattice. By changing the discretization of Ψ (that is, by tuning the lattice constant a), we have access to different values of a/l B . This is valuable to have a one-to-one comparison with HH systems [where this ratio is set by α -see Eq. (2)], but also to perform a more systematic study of the convergence towards the continuum limit (a/l B → 0).
Note, to conclude, that other choices are available for the discretized version of Laughlin-like wave functions, see for instance Ref. [77], and we expect that similar results can be obtained for all of them.

III. ANYONIC STATISTICS FROM DENSITY-PROFILE MEASUREMENTS
We now consider FQH and FCI states of N particles, in the presence of some localized QH excitations. We start by reviewing a recent proposal to characterize the braiding phase of these excitations through density profile measurements in the continuum [59,68], and we then conjecture on its generalization to the lattice case.

A. Continuum systems
In Refs. [59,68], some of the current authors developed a scheme to access the anyonic statistics for QH excitations of FQH states. In contrast to other proposals based on interferometric experiments [see for instance Ref. [78]], this proposal only requires static measurements. More precisely, a useful relation was obtained to relate the QH braiding phase ϕ br and the expectation value L z of the angular momentum operator taken on some specific quantum states. Each state hosts two QHs at positions η 1 and η 2 , which are either diametrically opposite with respect to the system center or on top of each other. For FQH states in the LLL, the correspondence between L z and the mean square radius r 2 [79] simplifies our scheme even further, since it implies that one can obtain the braiding phase simply by measuring the density profile in the presence of QHs. More precisely, L z / + N = N r 2 /2l 2 B and the QH braiding phase reads [68] The subscripts of r 2 refer to two aforementioned QH configurations that need to be considered: "opp" indicates diametrically opposite QHs (i.e., η 1 = − η 2 ), while "over" indicates overlapping QHs (i.e., η 1 = η 2 ). We stress that, in Eq. (5), l B is the actual magnetic length l B ≡ /qB, depending on the particle charge q and on the magnetic field B, and not a model dependent quantity.
Several practical issues appear when applying Eq. (5), due to the fact that it involves global properties of the FQH cloud: First, one has to measure r 2 opp and r 2 over on systems with exactly the same number of particles N and with the same value of | η 1 | and | η 2 |. Second, one needs to consider large-enough system sizes, so that η 1 and η 2 (in the "opp" configuration) are far enough from each other and from the cloud boundaries. Third, global properties of the system, like r 2 , are not robust with respect to low-energy perturbations, which in the case of pinned QHs typically consist in the excitation of edge modes [80,81]. These three issues are drastically mitigated once we rewrite Eq. (5) in terms of the depletion d( ρ), which is the change in the FQH density profile at position r induced by a QH at position η, with r = η + ρ. More precisely, we define the depletion profiles d 1QH ( ρ) and d 2QH ( ρ) close to single or double QH as Here n kQH ( r) represents the average density on a state with k QHs at position η, for k ∈ {0, 1, 2}. By computing d 1QH ( ρ) and d 2QH ( ρ), the braiding phase can be expressed as [59]: For finite-size systems the integral is defined up to a cutoff distance | ρ| = R max , which should be significantly larger than the QH size, but also small enough to avoid spurious effects due to the FQH cloud boundaries. This guarantees the mathematical equivalence between Eqs. (5) and (7). The dependence of Eq. (7) on the cutoff R max is characterized by damped oscillations which converge towards the actual value of ϕ br [59], as visible in Figs. 4(b) and 5 (c) and (d).
The new expression for the braiding phase, Eq. (7), has some clear advantages over Eq. (5). First, it only depends on the local density perturbation induced by the QHs, rather than on the global shape of the cloud. This simplifies the measurement, which is now independent on the precise position of the QHs. Second, the constraint on the cutoff distance R max is milder than the one on the distance | η 1 − η 2 | in the "opp" configuration, since now there is no need to consider two spatially separated QHs. Third, local properties like the depletions in Eq. (7) are not modified by perturbations which are localized at the edge of the system, so that this measurement is expected to be more robust against edge excitations and finitetemperature effects.
In previous works, we numerically confirmed the validity of Eqs. (5) and (7) by a Monte Carlo study of two paradigmatic FQH states. We considered the Laughlin [68] and Moore-Read [59] wave functions, and for the latter we focused on the case of two non-Abelian QHs [7,11]. Nevertheless, the procedure is totally general, and it could be applied to any state in the LLL.

B. Lattice generalization
While in the previous works we presented the theory relating the anyonic statistics of QHs to their depletion profiles in a continuum geometry, here we discuss how these relations [Eqs. (5) and (7)] change if one considers lattice systems.
The mean square radius r 2 now reads where r j is the position of the j-th site. With this definition, it is straightforward to generalize Eq. (5) to the lattice case: Similarly, we define the depletions d 1QH ( ρ j ) and d 2QH ( ρ j ) as where n j kQH is the average density on site j for a state with k QHs at position η, and where r j = η + ρ j . Thus Eq. (7) becomes where the sum over j is restricted to sites with | ρ j | < R max , as for Eq. (7). Note that both expressions for the braiding phase ϕ br [see Eqs. (9) and (11)] explicitly depend on the effective magnetic length l B , defined in Eq. (2). Before moving on, we need to stress that the angular momentum operator is not properly defined on a lattice, and therefore the relation between the QH braiding phase and the density profile [Eqs. (9) and (11)] is not mathematically rigorous in this case. However, the idea of generalizing Eqs. (5) and (7) to the case of FCIs on a lattice is motivated (and partially justified) by two observations: First, the study by Liu and co-authors [46] clearly shows that, for systems with periodic boundary conditions (PBCs), the density profile close to a single FCI QH is a discrete sampling of the continuum case, once a suitable (model-dependent) effective magnetic length is introduced. In Sec. IV, we explicitly confirm this result to the case of the HH model with OBC, both for a single and a double QH. Second, we perform a complementary analysis of the lattice case, based on the discretized Laughlin wave function described in Section II B. This flexible ansatz allows us to scan different values of the grid spacing a, and to compute the braiding phase through the discretized versions of Eqs. (5) and (7). In Table I, we report the numerical results obtained via Eq. (9), while we will use Eq. (11) in Section IV B. We find that the braiding phase of the QH excitations of the discretized Laughlin state is in full agreement -up to some small deviations due to finite-size and discretization effects-with the expected value ϕ br /(2π) = 1/2, within the statistical uncertainties of the Monte Carlo method. It is remarkable that this still holds true for a/l B 1.77, which corresponds to the maximum flux density that can be realized in the HH model, i.e. α = 1/2. In the next section we verify the validity of Eq. (11) on the ground states of the interacting HH model, obtained with the TTN technique.

IV. RESULTS
In our study of the HH model through the TTN ansatz, we consider two sets of parameters: N = 12 bosons and α = 0.15 (Case I), or N = 18 particles and α = 0.25 (Case II). These two choices for α correspond to a/l B 0.97 and a/l B 1.25, respectively [see Eq. (2)]. Considering different values of α allows us to modify N in a significant way without changing the size of the lattice, and it also gives us the opportunity to inspect discretization and flux-dependent effects. Note also that for Case II we choose α = 0.25, which is one of the most appealing flux densities for realizing almost flat bands in realistic experiments [60][61][62][63]. We also introduce an additional harmonic confinement in the form V j = Ω| r j − r center | 2 , where the trap center corresponds to the center of the L × L lattice. The value of Ω cannot be chosen arbitrarily, since it is not guaranteed that the ground state is an FCI for any trapping strength. The main FCI signature is the formation of a flat-density region in the center of the trap (with density equal to α/2), which signals its incompressibility. Through the TTN technique, we identify a set of Ω's where the density profile shows an incompressible central region with density α/2. In Case I, this occurs at least for Ω/t between 10 −4 and 3 × 10 −3 , while for Case II this happens for Ω/t between 2 × 10 −3 and 10 −2 . In the following, we choose trapping strengths within these intervals, namely Ω/t = 1 × 10 −3 for Case I and Ω/t = 3×10 −3 for Case II. A systematic study of the stability of the incompressible core as a function of both the magnetic flux density α and the confining strength Ω is certainly interesting for experimental purposes, but it is left for future work.
Concerning the TTN ansatz, we use bond dimensions as large as D = 500 (for the largest N ), and we verify that the systematic error in the relevant observables is negligible [see Appendix B].
The two sets of parameters used in the numerical calculations are summarized in Table II. In the following, we describe how to stabilize QH states for these two cases and we apply the procedure described in Section III B to demonstrate their anyonic nature.  A. Stabilization and charge of the quasiholes As a preliminary step before entering into the discussion of the braiding phase, we show that it is possible to stabilize states with either one or two (overlapping) QHs. To do so, we make use of some localized pinning potentials with different strengths and shapes (e.g. point-like potentials acting on a single-site, plus-shaped potentials, square-shaped potentials, etc.). For systems with PBCs, it is already known that some of these potentials can select states with a single QH as the ground state of the system [46,48,49]. On the other hand, only recently Račiūnas and co-workers have inspected the stabilization of QHs in systems with OBCs [47]. However, the use of exact diagonalization techniques limited their analysis to systems with small number of particles (N = 4), preventing them from clearly observing the expected QH fractional charge.
We consider the Hamiltonian in Eq. (1) with Case I parameters and we add a pinning potential with finite intensity V i on each of the four sites of a plaquette at the center of the trap (squared-shaped potential). Then we compute the charge Q l , defined as: where the sum is restricted to the sites j of a l × l square at the center of the trap. For large V i , the density on the central 2 × 2 plaquette tends to zero, so that (for l = 2) Q l tends to −4 × (α/2) = −0.3, which is minus four times the bulk density in the absence of pinning potentials. This is clearly visible in the lower panel of Fig. 2. Note that for α = 0.25 the plateau in Q 2 would occur exactly at −1/2, independently on the number of QHs pinned by the impurity potential, and this could lead to the misinterpretation of the data. The correct values for the (fractional) charge of the QH excitations are recovered for large enough l. The lower panel in Fig. 2 clearly shows that, by increasing the potential strength, Q l saturates at −1/2 and −1, which are the expected charges for one and two QHs, respectively. Note that the exact location of the transitions between 0, 1 and 2 QHs also depends on the trapping strength Ω, on the flux density α, and on the particle number N . Along this line, there are two important remarks we want to make. The first one concerns the robustness of the bulk against weak perturbations: While for pinning potentials strong enough to create some QH excitations, the density depleted by the QHs is displaced to the system boundary [see middle and right pictures in Fig. 2, top panel], weak repulsive potentials modify the systems density only around the potential position, giving an overall depleted charge Q l 0 for sufficiently large values of l and no accumulated charge on the boundary [see left picture in Fig. 2, top panel]. Note also that the robustness of the bulk is expected to increase with the number of particles, since in the large-N limit the density depleted by the QH must be pushed much farther away to reach the system boundary.
The second point concerns the role played by harmonic confinement, which we found to be very important to efficiently stabilize the QH states. The advantage of introducing an additional harmonic trap is twofold: On the one hand, it shifts the edge modes of the FCI state to higher energies. On the other hand, it allows the system to automatically regulate its spatial extension to minimize the energy. Thus the system can adapt itself to the presence of the pinning potentials by putting QHs in correspondence with the potentials and displacing the extra density to the boundary, without the need of removing particles or modifying the flux density α. In other words, in the presence of a harmonic confinement the system automatically chooses the right amount of flux quanta in order to be an FCI at filling ν = 1/2 with the most suitable number of QH excitations. Note that this behavior is completely different from what usually happens in systems with OBCs (without additional trapping potentials), where the total number N QH of QH excitations is set by the fact that ν = (N + N QH /2)/(L − 1) 2 α. Even though we considered squared-shaped potentials in our analysis of the QH charge, we verified that similar results can be obtained for pinning potentials of different forms [see for instance the single-site and the plus-shaped potentials used in Sec. IV B to stabilize one and two QHs, respectively]. The discussion about the pros and cons of the different pinning potentials, and the effects they have on the QH properties, is postponed to the next section.

B. Quasihole statistics
As discussed in Sec. III A, the method proposed in Ref. [59] has the great advantage of allowing us to determine the QH braiding phase in a setup that would be prohibitive for the more traditional schemes involving real braiding and interference. By lifting the needs of an explicit rotation of the QHs, and of their spatial separation, this method allows us to work with trapped samples with OBCs (rather than with abstract periodic structures) also when the number of particles is relatively small [82]. On this basis, we decided to apply the lattice generalization of that protocol -see Sec. III B-to the FCI QHs under study [see parameter sets in Table II]. We stress that both in Case I and in Case II, the system size is too small to accommodate two QHs far enough from each other and from the system boundary to apply the usual braiding schemes.
The key ingredients to measure the QH braiding phase, through the protocol described by Eq. (11), are the density depletions d 1QH ( ρ j ) and d 2QH ( ρ j ). To take advantage of Eq. (11), the set of vectors { ρ j }, i.e. the distances between the lattice sites and the QH center, must be the same for both d 1QH and d 2QH . As a consequence, both the single QH and the two overlapping QHs must be centered at the same position. To satisfy this constraint we identified two possible configurations: (i) a single-site potential to pin the single QH and a plus-shaped potential (centered on the same site and spread over five sites) to pin the two overlapping QHs; (ii) a 2 × 2 square-shaped potential with different strengths to selectively pin one or two QHs.
In order to extract the depletions d 1QH ( ρ j ) and d 2QH ( ρ j ), one needs the density profiles n j kQH characterizing states without QHs (k = 0), with a single QH (k = 1), and with k = 2 overlapping QHs [see Eq. To simplify the interpretation of the figures, we use the following marker notation: First, blue (orange) markers denote the data obtained through the TTN technique for the depletion due to a single (double) QH. Second, the different marker shapes are associated with different physical parameters and/or pinning potentials, as summarized in Table III Table III] to create the quasiholes in a system with Case II parameters. Panels (c) and (d): Comparison between the quasihole braiding phase, as a function of the cutoff radius Rmax, obtained for the interacting HH model with different pinning potentials (red markers), for the discretized quasihole wave functions in Eq. (3) with a 1.25 lB (black crosses), and for their continuum version (black dashed line). In panels (a) and (c) quasiholes are centered on a lattice site, while in (b) and (d) they are located at the center of a plaquette. The TTN bond dimension is D = 500 for all data sets.
corresponds to the behavior of ϕ br for the Laughlin QHs in the continuum.
The comparison of the QH density depletions obtained through the TTN technique with those characterizing the Laughlin QHs is in general very good [see panel (a) in Fig. 4 and panels (a) and (b) in Fig. 5], up to small deviations. At large distances there are discrepancies due to finite-size effects which affect both d 1QH ( ρ j ) and d 2QH ( ρ j ). To be precise, these are caused by the presence of the density bump close to system boundary, which modifies the tail of the depletions [83]. To mitigate this problem, we restrict the plotted data to distances for which the depletions behave smoothly and the QH charges are close to the expected values (i.e., −0.5 and −1 for the single and double QH, respectively). At small and intermediate distances, instead, the reason why the QH density depletions appear slightly different with respect to their continuum counterparts is twofold. On the one hand, there is the effect of the pinning potentials [see for instance the radial profiles of the 2QH density depletions reported in panel (a) of Fig. 5], which is more evident for α = 0.25 than for α = 0.15. This dependence on the flux density is due to the fact that larger α's correspond to larger values of the ratio a/l B . This implies that, in the case of α = 0.25, the same (finitewidth) pinning potential used for α = 0.15 effectively acts at larger distances from the center of the QH, where the system density is higher and the potential can have stronger effects on the QH depletions d 1QH and d 2QH . On the other hand, there are unexpected deformations which are present also in the density depletion of the single QH pinned by a single-site potential [see blue squares and blue diamonds in panel (a) of Fig. 5 at short and intermediate distances]. Zero-range pinning potentials should not deform the QH density depletions and therefore we consider the observed discrepancies as model-dependent effects, probably due to the specific properties of the HH bands for a given value of α. This would explain why this behavior is more evident for Case II (α = 0.25) than for Case I (α = 0.15).
Finally, we used the method described in Sec. III B to extract the braiding phase ϕ br of the lattice QHs, as a function of the cutoff radius R max [see Figs. 4 (b) and 5 (c) and (d)]. We compare it with the results obtained for both the QHs of the continuum-space Laughlin state and their discretized counterparts. Note that for the sampling of the discretized QH wave functions we chose a different discretization grid spacing for Case I and Case II, since they are characterized by a different value of the ratio a/l B [see Eq. (2)].
For Case I, the behavior of ϕ br obtained for the QHs of the interacting HH model is extremely similar to the one predicted by the Monte Carlo sampling of the discretized Laughlin QHs up to a certain value of R max [see Fig. 4, panel (b)]. After that, at larger cutoff radii, we observe small deviations between the two data sets, reflecting the deformations in the depletion d 2QH ( ρ j ) caused by the plus-shaped pinning potential. However, the results shown in Fig. 4(b) clearly indicate that for Case I the anyonic nature of the QHs in the interacting HH model can be probed through simple density measurements.
The interpretation of the QH braiding phase obtained for Case II is more subtle. As we have already discussed, for α = 0.25 the density depletions of the lattice QHs display more visible discrepancies, with respect to their continuum counterparts, than for α = 0.15. Moreover, due to the | ρ j | 2 factor in Eq. (11), these discrepancies in the depletions translate into even stronger effects affecting the behavior of ϕ br (R max ). This is clearly visible in panel (c) of Fig. 5. Despite discrepancies in the profile of the density depletions, the method proposed in Eq. (11) is still valid and the correct results for ϕ br (R max ) are recovered for large enough integration regions -where the deformations in the density depletions are damped. We attribute this behavior to the topological robustness [84] of the QH braiding properties and we expect further confirmation of this interpretation from future studies of of larger systems. In spite of these additional difficulties, the results obtained for Case II confirm that the anyonic nature of the QHs can be inspected through simple density-profile measurements also in the experimentally promising α = 0.25 case.
Even though our method to measure QH braiding phase seems to be robust against the deformations induced by extended pinning potentials in the depletion profiles, they affect the behavior of ϕ br (R max ) at small and intermediate values of the cutoff radius. To get rid of these effects, one could in principle create multiple QHs at the same position by locally inserting the suitable amount of flux quanta, in the presence of a single-site potential able to pin a single QH [36,47,85]. Despite being experimentally feasible, in the theoretical framework this flux-insertion procedure requires time-dependent simulations of the interacting HH Hamiltonian, which at the moment go beyond the capabilities of our TTN technique. Along this line, a recent application of the TTN ansatz for time-dependent simulations opens interesting perspectives [86].

V. CONCLUSIONS AND OUTLOOK
In this work we used a Tree Tensor Network ansatz to study the properties of the quasihole excitations of the fractional Chern insulator described by the Harper-Hofstadter Hamiltonian with hardcore interactions. The loop-free geometry of the Tree Tensor Network ansatz allowed us to study systems with open boundary conditions, far beyond the typical system sizes manageable by exact diagonalization calculations.
In this way, first we showed that it is possible to use localized pinning potentials to stabilize states hosting either a single or two overlapping quasiholes and that the expected fractional charge of these excitations is already clearly visible in a N = 12 particle system. In this respect, we discovered that superimposing an additional harmonic confinement to the lattice greatly simplifies the stabilization of the quasihole states.
Then, to characterize the statistics of the quasiholes, we applied a lattice version of the equation proposed in Ref. [59] and relating the quasihole braiding phase to the depletions induced by such excitations in the system density. Our results clearly show that these excitations are anyons, namely that they are neither bosons nor fermions (for which ϕ br = 0, 2π), and that their braiding phase is very close to the predicted one, i.e. ϕ br = 2πν. In spite of the obvious limitations in the accuracy of the measurement of ϕ br , mainly due to the size of the state-of-the-art samples, we stress that our results have been obtained for systems which are too small to accommodate two spatially separated quasiholes and adiabatically braid them to inspect their statistical properties, as it would be required by traditional measurement schemes.
As a result, the present study provides numerical evidence that the anyonic statistics can indeed be observed through simple density measurements in state-of-the-art experiments with ultracold atoms and superconducting qubits. First of all, the flux densities we considered in this work -i.e., α = 0.15 and 0.25-are already within the current experimental capabilities. The case of α = 0.25, in particular, is of great interest from the experimental point of view. At such flux density the single-particle spectrum of the HH Hamiltonian is characterized by four energy bands, with a very convenient (low) ratio between the width of the lowest band and its separation from the higher ones. This makes α = 0.25 one of the most appealing flux densities for realizing almost flat bands in realistic experiments [60][61][62][63]. At the same time, the lattice size we looked at -i.e, L = 16-is comparable with the one used in Ref. [63]. On top of that, adding an overall harmonic confinement to the lattice should be possible in both relevant setups: For ultracold atomic sys-tems the harmonic trap is typically already present on top of the lattice potential (and in most cases it is difficult to remove); while platforms based on superconducting qubits should allow one to independently tune the different on-site potentials V i [62,87]. Finally, the hardcore constraint describes well the on-site interactions in superconducting qubits [62]. This condition might be difficult to reach with ultracold atoms; however, as suggested in Refs. [40,43,47], we believe that considering finite-strength interactions instead of the hard-core ones should not modify our results in a considerable way. A comprehensive analysis of the case of soft-core interactions is left for a future work.
From the theoretical point of view, our work, combined with the observations by Liu and co-workers [46], numerically shows that the expression relating the QH braiding phase to the QH density depletions is valid also for more general lattice systems, once a suitable effective magnetic length is introduced. Since on a lattice the angular momentum operator is not properly defined, such an expression can not be rigorously derived as done in Ref. [59] for the continuous system. This opens the interesting question whether it is possible to find a deeper and more general explanation of the link between the statistical properties of the QHs and their density profiles, which seem to be independent of the model under study.
Possible extensions of this work include other FCI states, and in particular those hosting non-Abelian QHs [88][89][90][91][92].  To perform a systematic benchmark of the TTN results against exact diagonalization (ED), we consider N = 3 particles on a 8 × 8 lattice with periodic boundary conditions (PBCs), and with α = 1/8. For these parameters, the magnetic filling is such that the system hosts two QHs. To localize these two QHs at the same position, we also include a five-sites plus-shaped pinning potential centered at site (4,4), with a repulsion V i /t = 1 on each site.
Differently from the OBC cases treated in the main text, the toroidal geometry of a PBC system induces a 1/ν topological degeneracy of the ground state [93], that is, a twofold degeneracy for the case treated in this work. This degeneracy is preserved when the system hosts localized Abelian QHs (up to finite-size effects), in contrast to the general case of non-Abelian or non-localized QHs, in which the number of (quasi-)degenerate states is typically larger [31,48].
For N = 3 and L = 8, the Hilbert space dimension is ≈ 4.2 × 10 4 , and we can obtain the three lowest-energy states Ψ 0 , Ψ 1 and Ψ 2 through the ED technique (via the Lanczos algorithm). The corresponding energies [see Table IV] display a signature of the (quasi-)degeneracy of Ψ 0 and Ψ 1 , while the state Ψ 2 has a larger energy gap -see Ref. [42] for a study of the spectrum of larger systems. For this small system, the TTN ansatz has moderate requirements in terms of the bond dimensions of its tensors. In general, within the hierarchical tree structure of a loop-free TTN state, the bonds of a tensor on the l-th layer may have dimension up to d 2 l [51,53], where d is the dimension of the local Hilbert space. To make this ansatz tractable, one typically introduces a cutoff D on the bond dimension, so that the dimension of each bond becomes min(d 2 l , D). For L = 8, the root tensor (i.e. the one at the top of the tree structure) has a maximum bond dimension of min(d 2 5 , D) = min(d 32 , D).
For the hard-core Bose-Hubbard model treated in this work, the local dimension is d = 2. By considering a fixed number N of bosons, the number of states is reduced, and we can compute the bond-dimension D exact needed to span the entire many-body Hilbert space. For N = 3 and L = 8, a TTN state with D ≥ D exact = 66 can reproduce an arbitrary many-body state. We then set D = 66 in this Appendix, and show that the optimization of the TTN parameters yields the same physical results as ED.
First of all, the TTN energies for Ψ 0 , Ψ 1 and Ψ 2 perfectly match with the ED values in Table IV, up to the precision reported therein. We also compute the density profile close to the pinned QHs, as shown in Fig. 6 for Ψ 0 and Ψ 1 . These two degenerate states have different profiles, but once we average over the two states we obtain a smooth curve, as known from Ref. [46]. Also for the profiles, as for the energies, we observe full agreement between ED and TTN results. This further validates the TTN method, also in the presence of a localized QH excitation.

Appendix B: Bond dimension dependence of observables
In App. A, we showed that for a small system one can reach the saturation of the bond-dimension cutoff (D = D exact ), which makes the TTN ansatz a general parametrization for any state in the Hilbert space. The computational cost of this task, however, becomes prohibitively large for larger systems. For Case I (Case II) in the main text, such saturation would require a bond dimension D 6 × 10 9 (D 2 × 10 13 ). Since this is clearly beyond reach, the TTN ansatz can only represent states in the subset of Hilbert space with low to moderate entanglement content, which introduces a bias in the computed observables.
As a representative example, in this Appendix we consider the parameters of Case II reported in Table II, supplemented with a plus-shaped pinning potential of intensity V i /t = 1. In this setup the system has a localized double QH -see orange diamonds in Fig. 5(a). We first look at the energy, for bond dimension D between 100 and 500 -see data in Table V. Even if the large-D convergence is not reached, the energies for the two largest bond dimensions are very close, with a relative energy difference as small as 1.5 × 10 −4 . This suggests that the TTN ansatz with the currently available bond dimension may be sufficiently accurate to represent the ground state of the model under study. Note that the relative energy differences reached in this work are comparable with the TTN analysis for the case of a homogeneous FCI [42]. Furthermore, we look at how the depletion profile around a QH (the key observable in the current work) depends on the TTN bond dimension. We look again at the same case (a double QH in Case II, with a plusshaped pinning potential), and compute the depletion This quantity is shown in Fig. 7, for D between 100 and 400. We observe that its fluctuations decrease for increasing bond dimension, and that the curve for D = 400 is barely distinguishable from 0, on this scale. We conclude that for a TTN with D = 500 (the bond dimension used for Case II in the main text) the systematic error in the depletion profile is negligible.