Fractional quantum Hall effect in the interacting Hofstadter model via tensor networks

TL;DR

This paper demonstrates, using tensor network methods, that the Harper-Hofstadter model with hard-core bosons on a lattice exhibits a topological fractional quantum Hall phase at filling factor 1/2, with robust properties and measurable topological invariants.

Contribution

It provides the first tensor network evidence of a lattice fractional quantum Hall state in the interacting Hofstadter model, including topological entanglement entropy estimation.

Findings

Supports the $ u=1/2$ fractional quantum Hall phase on the lattice.

Measures the many-body Chern number confirming topological order.

Shows gapless edge modes via Green function analysis.

Abstract

We show via tensor network methods that the Harper-Hofstadter Hamiltonian for hard-core bosons on a square geometry supports a topological phase realizing the $\nu=1/2$ fractional quantum Hall effect on the lattice. We address the robustness of the ground state degeneracy and of the energy gap, measure the many-body Chern number, and characterize the system using Green functions, showing that they decay algebraically at the edges of open geometries, indicating the presence of gapless edge modes. Moreover, we estimate the topological entanglement entropy by taking a combination of lattice bipartitions that reproduces the topological structure of the original proposals by Kitaev and Preskill, and Levin and Wen. The numerical results show that the topological contribution is compatible with the expected value $\gamma = 1/2$. Our results provide extensive evidence that FQH states are within reach of state-of-the-art cold atom experiments.

Figures

Equations

$$H$$

$$P({\theta }_x,{\theta }_y)$$

$${\Omega }_{\Phi \to {\Phi }^{\prime }}=\arg \left(\det ⟨{\Phi }_j^{\prime }|P({\theta }_x,{\theta }_y)|{\Phi }_k⟩\right)$$

$${\mathcal{I}}_x(x,y)=i⟨a_{x+1,y}^{†}{a}_{x,y}-a_{x+1,y}{a}_{x,y}^{†}⟩,$$

$$-\gamma$$

$$\tilde{H}=H+\mu \sum _{(x,y)\in sup.latt.}{n}_{x,y},$$

$${\Gamma }_{\Phi }({\theta }_x,{\theta }_y)=⟨\Phi |{\Psi }_0({\theta }_x,{\theta }_y)⟩⟨{\Psi }_0({\theta }_x,{\theta }_y)|\Phi ⟩,$$

$$E_0=⟨{\Psi }_0|H|{\Psi }_0⟩\stackrel{!}{=}min.$$

$$E_n=⟨{\Psi }_n|H|{\Psi }_n⟩+\sum _{k=0}^{n-1}{ϵ}_k⟨{\Psi }_n|P_k|{\Psi }_n⟩\stackrel{!}{=}min.,$$

Full text

I Introduction

The Harper-Hofstadter model Hofstadter (1976) plays an archetypical role in the current understanding of topological quantum matter on a lattice. It encompasses the basic coupling between particles and a background magnetic field, and supports topological bands with finite Chern number for a broad range of fluxes and tunneling rates Bernevig and Hughes (2013). Those remarkable properties have motivated proposals Dalibard et al. (2011); *Jaksch2003; *Goldman2014 and recent experiments in both solid-state Dean et al. (2013) and cold atom systems Aidelsburger et al. (2013); Miyake et al. (2013); Kennedy et al. (2015); Aidelsburger et al. (2015); Mancini et al. (2015); *Stuhl2015, which have extensively investigated its non-interacting limit, including the observation of its fractal spectrum Dean et al. (2013) and the measurement of a finite Chern number Aidelsburger et al. (2015) in some of its bands. More recently, experiments using bosonic atoms in optical lattices have shown impressive capabilities to approach the strongly interacting regime Kennedy et al. (2015), with the ultimate goal of stabilizing lattice analogues of fractional quantum Hall (FQH) states Bernevig and Hughes (2013); Fradkin (2013). However, despite promising small system results based on exact diagonalization Sørensen et al. (2005); Hafezi et al. (2007a); Möller and Cooper (2009); Sterdyniak et al. (2012), in the large flux regime available to experiments – when the magnetic length is of order of the lattice spacing – theoretical evidence of such states at large scales has been lacking, especially regarding smoking guns of topological order – gapless edge modes, entanglement properties, and many-body Chern numbers.

In this work, we study the strongly-interacting bosonic Hofstadter model, and show that it supports a FQH ground state (GS) akin to the $\nu=1/2$ Laughlin state in the continuum Fradkin (2013); Kalmeyer and Laughlin (1987); Regnault and Bernevig (2011). Our analysis is based on a combination of diagnostics, including GS degeneracies on different topologies, the measurement of the many-body Chern number, Green functions’ decay at the edge and in the bulk, and a direct measurement of the topological entanglement entropy. The corresponding results serve as a quantitative guideline to address the stability of such phases against temperature: crucially, we show how the spectral gap of the bulk excitations is of the order of $10\%$ of the tunneling rates, showing how FQH states are within reach for temperatures available to current experiments.

The enabling tool of our analysis are numerical simulations based on the tree-tensor network ansatz Tagliacozzo et al. (2009); Murg et al. (2010); *Nakatani2013; Gerster et al. (2014). This class of tailored variational wave functions extends the matrix-product state (MPS) ansatz – the tensor network class at the heart of the density-matrix renormalization group (DMRG) White (1992); Schollwöck (2011). To obtain the presented results on square lattices with sizes of up to $32\times 32$ sites, we exploit the reduced scaling of computational costs of the loopless geometry of TTNs, a characteristic not shared by other network structures such as PEPS and MERA Verstraete et al. (2006); *Vidal_2007. Additionally, differently from typical DMRG approaches, TTN warrant direct access to the reduced density matrix of a variety of lattice bipartitions. In the following, we show how this feature enables the direct evaluation of the topological entanglement entropy using a specific combination of various partitions, in parallel to the original proposals of Refs. Kitaev and Preskill, 2006; *Levin2006, which are instead not immediately applicable to DMRG studies.

The paper is structured as follows: First, we define the model (Sec. II) and describe the employed numerical method (Sec. III). Then, we proceed to present the numerical evidence for a FQH GS in Sec. IV, based on the low-energy spectrum (Sec. IV.1), the many-body Chern number (Sec. IV.2), the correlation functions (Sec. IV.3), and the topological entanglement entropy (Sec. IV.4). In Sec. V we analyze the robustness of the FQH state upon introducing a superlattice potential. Finally, we conclude our work in Sec. VI.

II Model Hamiltonian

We study spinless bosonic particles hopping on a $L\times L$ square lattice under the influence of an external magnetic field, as illustrated in Fig. 1. In the Landau gauge, the system is described by the following Hamiltonian Sørensen et al. (2005)

[TABLE]

of bosonic particles, $[a_{x,y},a_{x^{\prime},y^{\prime}}^{\dagger}]=\delta_{xx^{\prime}}\delta_{yy^{\prime}}$, $n=a^{\dagger}a$. Here, $\phi$ is the magnetic flux through each plaquette (resulting in a magnetic filling factor $\nu=N/(\phi L^{2})$, with $N$ the number of bosons in the system) and $\theta_{x}$, $\theta_{y}$ implement the twists in the boundary condition Hafezi et al. (2007a). In the dilute limit (small densities and small fluxes) the lattice physics approaches the one of the continuum 111Notice that the continuum limit is also recovered by adding tailored long-range hoppings, which effectively flatten the lowest band Kapit and Mueller (2010).. However, in the large flux limit, available in cold gases experiments, the phase diagram is not set. On small systems, it has been shown by exact diagonalization (ED) that the GS of the model described by Eq. (1) at filling factor $\nu=1/q$ (where $q$ is an even integer) is compatible with a lattice analogue of the (bosonic) Laughlin wave function Hafezi et al. (2007a); Laughlin (1983); Möller and Cooper (2009); Sterdyniak et al. (2012); Hügel et al. (2017), exhibiting topological GS degeneracy and a non-zero Chern number Niu et al. (1985); *Tao1986. However, the overlap with the exact Laughlin wave functions rapidly degrades with system size already for small systems of 6 particles. From a complementary viewpoint, the ladder version of the Hofstadter model has also been shown to share similarities with FQH states Budich et al. (2017); Strinati et al. (2017); Petrescu et al. (2017); Haller et al. (2017). Very recently, iDMRG results on cylinders have shown strong signatures of integer quantum Hall states, and have reported fractional current quantization in regimes different from the one we consider here He et al. (2017). Throughout, we focus on the strongly interacting case $U\rightarrow\infty$ (hard-core bosons) with flux values $\phi=1/8$ and $1/16$ respectively, which correspond to flux setups that are experimentally available. Finally, we fix the energy scale by setting $t=1$.

III Tree-tensor network ansatz

We employ a tree-tensor network (TTN) ansatz Shi et al. (2006) for the GS and the two lowest excited states to verify the properties discussed above. The specific binary TTN used in this work is illustrated in Fig. 2, where the standard graphical notation for tensor networks (TN) is employed Orús (2014); Silvi et al. (2017): tensors are drawn as cubes, with attached lines symbolizing tensor indices (links). Links that are shared by two tensors are contracted, which in TN language means that over their corresponding mutual indices is to be summed. The dimension of each link in the TTN is upper bounded by a constant $m$ (bond dimension), which serves as the refinement parameter of the ansatz: the larger $m$, the more accurate the true many-body state can be approximated. In a binary TTN each tensor has at most three links; therefore, the scaling of the computational resources with the bond dimension is moderate in algorithms using this class of TN states Gerster et al. (2014). Furthermore, we exploit particle number conservation by restricting the ansatz to the $N$ particle symmetry sector Singh et al. (2011).

While it is known that a 2D-TTN is not compatible with the area law for the entanglement entropy Tagliacozzo et al. (2009); Ferris (2013), it possesses several beneficial features which make it a promising tool for the study of intermediate system sizes: (a) The existence of a numerically stable search algorithm for eigenstates Murg et al. (2010); *Nakatani2013; Gerster et al. (2014); Silvi et al. (2017); (b) a low-order polynomial scaling $\mathcal{O}(m^{4}L^{2})$ of the computational cost; (c) easy interchange of various boundary conditions (open, periodic, twisted); (d) access to the entanglement entropy for bipartition shapes that enable the determination of the topological entanglement entropy (TEE) Kitaev and Preskill (2006); *Levin2006. In what follows, we will exploit these properties to gather a number of numerical pieces of evidence supporting a FQH GS of the model Eq. (1) in the case of filling $\nu=1/2$.

IV Numerical evidence for FQH ground state

IV.1 Low-energy spectra

As first evidence we verify the GS degeneracy, intimately connected to the topological order of the system Wen and Niu (1990); Fradkin (2013). On a torus geometry the GS degeneracy at $\nu=1/2$ filling is expected to be twofold (independent of the twist angles $\theta_{x}$, $\theta_{y}$), while the first excited state is to be separated by a bulk gap. We determined the three lowest-energy eigenstates, reported in Fig. 3a: We clearly observe a finite energy gap $E_{2}-E_{0}\approx 0.1t$ which is typically more than two orders of magnitude larger than the energy difference $E_{1}-E_{0}$ between the two states manifold (multiplet) in the GS for the system sizes considered here. Conversely, we verified that for open boundary conditions (OBC) the quasi-degeneracy is removed and we observe that $E_{2}-E_{1}\approx E_{1}-E_{0}$.

IV.2 Many-body Chern number

As second evidence we determine the many-body Chern number (MBCN) Niu et al. (1985); *Tao1986, which is a direct signature of topological order. For the numerical calculation of the MBCN we follow the prescription put forward by Hatsugai Hatsugai (2004); *Hatsugai_2005, which requires the knowledge of the GS manifold on a 2D grid of twist angles $(\theta_{x},\theta_{y})\in[0,1]\times[0,1]$ (note that this grid also has a torus geometry). From a practical point of view, this method relies on the calculation of overlaps, a task which can be easily accomplished with TN states: More specifically, we need to choose two reference multiplets $\Phi$, $\Phi^{\prime}$ which have to be non-parallel but otherwise can be arbitrary. (Commonly, one picks two GS multiplets at twist angles far from each other Hafezi et al. (2007b).) This corresponds to two different gauge references, which can be used to define two scalar fields: $\Lambda_{\Phi^{(\prime)}}=\det\langle\Phi^{(\prime)}_{j}|P(\theta_{x},\theta_{y})|\Phi^{(\prime)}_{k}\rangle$, where the determinant runs over the indexes $j,k\in\{0,1\}$ of the GS multiplet and

[TABLE]

is the projector on the GS manifold at $(\theta_{x},\theta_{y})$. If the gauge has been fixed appropriately, $\Lambda_{\Phi}$ is well-defined (i.e. non-vanishing) where $\Lambda_{\Phi^{\prime}}$ is not and vice versa, thus forming two complementary regions on the boundary condition torus (an example is shown in Fig. 3b-c). Finally, the MBCN is given by the number of branch vortices in the argument field

[TABLE]

to be counted (with sign) in any of these two regions Hatsugai (2004); *Hatsugai_2005; Hafezi et al. (2007b). We performed this procedure using TTN states; the result is shown in Fig. 3d-e. We obtain a MBCN of 1, clearly demonstrating the topological order in the GS of the model. Combined with the two-fold degeneracy discussed above, this shows that the system has a Chern number per state $C=1/2$, as expected for the $\nu=1/2$ Laughlin state.

IV.3 Correlation functions and edge modes

Indirect evidence of a fractional quantum Hall state can also be gathered by monitoring the behavior of the Green function in the system $\mathcal{G}(x,x^{\prime};y,y^{\prime})=\langle a^{\dagger}_{x,y}a_{x^{\prime},y^{\prime}}\rangle$. The Green function in periodic systems contains information about bulk excitations: in our case, it is expected to decay exponentially as a function of distance. In contrast, for open boundaries the Green function is expected to reveal the existence of gapless edge modes. Indeed, in the case of $\nu=1/2$, the phase field operator of the corresponding chiral Luttinger liquid edge mode is expected to have considerable overlap with the creation and annihilation operators on the lattice. In Fig. 4a, we plot the normalized $\mathcal{G}(x,x;1,y)$ as a function of $y$ both for PBC and OBC. In the case of PBC, the results clearly show that $\mathcal{G}$ decays exponentially as a function of distance, with a correlation length of approximately 2 sites - and thus, considerably smaller than our system sizes. Since for OBC translational invariance is broken, we consider different values of $x$ (marked by different colors as also illustrated in Fig. 4c). Here, two distinct regimes are visible: along the edge (here for $x\lesssim 4$) the correlation decays as a power law, indicating a gapless mode localized close to the boundary. In sharp contrast, far from the edges ($x\simeq L/2$) the decay becomes exponential, consistent with the PBC results. Our results on the correlation functions thus strongly confirm a finite bulk gap, coexisting with gapless edge modes.

An alternative route to detect chiral edge modes, likely more accessible in experiments, is given by the particle current density. Specifically, we consider here the current in the $x$-direction

[TABLE]

again both for PBC and OBC. In Fig. 4b we show the normalized $\mathcal{I}_{x}(x,y)$ as a function of $y$ and for different values of $x$ (again with the same color code). The results are in strong agreement with the ones obtained from the Green function analysis: For PBC, where there is only bulk, the current vanishes throughout the system, while for OBC we observe strongly enhanced currents along the edges, i.e. at $y\approx 1$ or $y\approx L$. In Fig. 4c we provide an illustration of the edge current present in an OBC lattice, giving an intuitive (albeit over-simplified) explanation for the behavior of the currents plotted in Fig. 4b.

IV.4 Topological entanglement entropy

An unambiguous indicator for topological order is the topological entanglement entropy (TEE) $\gamma$, as introduced in Refs. Kitaev and Preskill, 2006; *Levin2006. A finite value of $\gamma$ signals the presence of anyonic excitations above the GS manifold - remarkably, without having to directly access excited states. In particular, $\gamma$ is directly related to the quantum dimensions of the excitations. Extracting the TEE in numerical simulations is however challenging. System sizes available to exact diagonalization are typically too small to measure $\gamma$ in a meaningful fashion. DMRG simulations can reach considerable system sizes: in that context, one tries to extract the TEE $\gamma$ as a correction to the perimeter-law scaling of the entanglement entropy: $S(\ell)=c\,\ell-\gamma$, where $S(\ell)$ is the von Neumann entropy $S=-\mathrm{tr}[\rho\log_{2}\rho]$ for a spatial partition of perimeter $\ell$. In Fig. 5a, we show our results for $S(\ell)$ for PBC with $L=16$. In order to avoid strong finite size effects, we excluded from our fits the data with $\ell<10$. After extrapolation to $m\rightarrow\infty$, we get $\gamma=0.47\pm 0.19$, which is close to the value obtained with our maximum bond dimension. Despite the large error bar, mostly due to the fact that we have only access to few points in $\ell$, our result show that the system has $\gamma>0$, and the result is compatible with the expected TEE for a Laughlin state with $\nu=1/q$, $\gamma=\log_{2}\sqrt{q}$. While this procedure has been well tested in several models Depenbrock et al. (2012); Jiang et al. (2012), it is desirable to apply a method for obtaining the TEE which follows directly the original prescription: Such technique operates on regions with different shapes in order to exclude spurious effects. Within our TN ansatz this is possible considering the geometry in Fig. 5b, where it can be shown that the following relation holds:

[TABLE]

The basic building block of this procedure are square partitions of size $a\times a$. The results of the corresponding simulations, extrapolated in $m\rightarrow\infty$, are shown in Fig. 5c: for $a=4$ the results indicate a finite TEE of $\gamma\simeq 0.5$, again with rather large error bars (owing to the uncertainty in the $m$ extrapolation), but otherwise in good agreement with the expected behavior for $\nu=1/2$ FQH states. We note that for $a=1,2$ the results are far from this prediction, as expected, since for those values the partition size is smaller than the correlation length $\xi$ as computed from the Green function decay.

V Phase transition from FQH state to trivial insulator

Our diagnostics also allows for the detection of possible phase transitions. Here we demonstrate this by investigating the robustness of the FQH state upon introducing a superlattice potential into the model from Eq. (1), which, in turn, can be engineered in cold atoms experiments Sebby-Strabley et al. (2006); *Foelling2007. We consider the Hamiltonian

[TABLE]

where $\mu\leq 0$ denotes the depth of the superlattice potential. We focus on the case where the number of sites of the superlattice is equal to the number of bosons $N$ in the system. For sufficiently large $|\mu|$, the GS of $\tilde{H}$ is a product state with the particles pinned at the superlattice minima (see Fig. 6) which is, by construction, topologically trivial. The transition between the two qualitatively different GSs occurs at some critical potential depth $\mu_{c}<0$. In order to detect this transition we use the energy spectra on the boundary condition torus, and the MBCN, as shown in Figs. 7 and 8. With increasing depth of the potential, we observe how the GS quasi-degeneracy is removed (see Fig. 7).

As long as the GS multiplet is still quasi-degenerate, i.e. separated from the first excited state ($E_{2}-E_{1}\gtrsim E_{1}-E_{0}$) on the whole torus of twist angles, we measure a MBCN of one (see Fig. 8).

At the critical point $-1.3\lesssim\mu_{c}\lesssim-1.0$ this condition is no longer met and the MBCN vanishes, signaling a disappearance of the topological order. The topologically trivial character of such a GS can be evidenced by considering the quantity

[TABLE]

where $\Phi$ represents a gauge reference, formed by a single GS at an (arbitrary) reference point $(\theta^{[r]}_{x},\theta^{[r]}_{y})$. If the GS of the system is in fact topologically trivial, $\Gamma_{\Phi}$ is well-defined (i.e. non-vanishing) on the whole boundary condition torus. On the contrary, this is not the case if the GS has non-trivial topological character: in that case there is a finite region on the torus where $\Gamma_{\Phi}$ vanishes. Whenever this happens, the GS is quasi-degenerate and Hatsugai’s method for determining the MBCN (as described in Sec. IV.2) is applicable. In Fig. 6 we show $\Gamma_{\Phi}$ for different values of $\mu$, demonstrating how it gradually becomes well-defined on the whole boundary condition torus as the depth of the superlattice potential is increased.

VI Conclusions and outlook

We have presented extensive numerical evidence supporting the existence of a fractional quantum Hall phase in the Harper-Hofstadter Hamiltonian for hard-core bosons on a square lattice. Our analysis considered a wide range of independent diagnostics, including spectral properties, the many-body Chern number, correlation functions, currents, and entanglement entropies, and shows how TTN algorithms provide a flexible tool to address the interplay of gauge fields and interactions in two-dimensional systems. The results indicate that the correlation length is typically of the order of a few lattice sites, with corresponding gaps of order $0.1t$: these signatures point toward the fact that lattice fractional quantum Hall states can be realized in present cold atom experiments, as well as potentially in cavity array experiments Kapit et al. (2014); Owens et al. (2017). The temperatures and sizes required match those already available in experiments, and for which adiabatic state preparation protocols have very recently been shown to be applicable He et al. (2017); Motruk and Pollmann (2017).

Acknowledgements.

We thank M. Burrello and A. Sterdyniak for a careful reading of the manuscript, J. Budich and H.-H. Tu for discussions, and F. Tschirsich for contributing numerical libraries. Numerical calculations have been performed with the computational resources provided by the bwUniCluster project 222bwUniCluster: funded by the Ministry of Science, Research and Arts and the universities of the state of Baden-Württemberg, Germany, within the framework program bwHPC., the JUSTUS project, CINECA via the TEDDI project, and the MOGON cluster at the JGU Mainz. We acknowledge financial support from EU projects RYSQ and UQUAM, the German Research Foundation (DFG) through the SFB/TRR21, OSCAR and TWITTER, and the Baden-Württemberg Stiftung via Eliteprogramm for postdocs. S.M. gratefully acknowledges the support of the DFG via a Heisenberg fellowship.

Appendix A Algorithm for ground and low excited states with tree-tensor networks

The ground state search is performed according to a variational scheme, first developed in the context of DMRG/MPS ground state algorithms Schollwöck (2011), and which in the meantime has been generalized to arbitrary loopless TN geometries Murg et al. (2010); *Nakatani2013; Gerster et al. (2014); Silvi et al. (2017). The basic motivation is the following: Although the TN representation of the many-body state already reduces dramatically the number of coefficients in the state vector (as compared to the full Hilbert space dimension), solving the minimization problem

[TABLE]

for the entire TN state $|\Psi_{0}\rangle$ directly is in practice not feasible, because the number of coefficients in $|\Psi_{0}\rangle$ is typically still far too large to be handled by numerical eigensolvers. Instead, one applies an iterative strategy: out of all the tensors, which together constitute the TN state, only the coefficients of one (or two, dependent on the update scheme Hubig et al. (2015); Silvi et al. (2017)) tensor(s) are considered variational, while the others are taken to be fixed. This allows one to contract the physical Hamiltonian $H$ to a reduced, effective Hamiltonian $H_{\mathrm{eff}}$, only acting on the degrees of freedom of the variational tensor(s). For loopless TNs the resulting reduced optimization problem can again be formulated as a standard eigenvalue problem, whose size is now manageable by an eigensolver. One then successively targets all the tensors in the TN, thereby gradually decreasing the energy expectation value $\langle\Psi_{0}|H|\Psi_{0}\rangle$. This procedure is called sweeping. A sweep is completed after all the tensors in the TN have been updated once. For the TTN architecture employed here (sketched in Fig. 2), the ansatz contains $N_{s}-2$ tensors ($N_{s}=L^{2}$: number of lattice sites) and all involved contractions (both computation of $H_{\mathrm{eff}}$ and solving the reduced eigenvalue problem) have computational complexity $\mathcal{O}(m^{4})$ or less. Therefore, the GS algorithm runtime scales as $\mathcal{O}(N_{s}\,m^{4})$. Convergence of the GS energy is typically reached after less than ten sweeps.

To determine excited states, we employ a very similar algorithm to the one just described, with the small modification that orthogonality to all previously determined eigenstates is enforced. This can be achieved by penalizing overlap with these eigenstates; more in detail, in order to obtain the $n$-th excited state of $H$ (orthogonal to all lower-lying eigenstates $|\Psi_{k}\rangle$, $k\in[0,\,n-1]$) we solve the optimization problem

[TABLE]

where $P_{k}=|\Psi_{k}\rangle\langle\Psi_{k}|$ is the projector on the $k$-th eigenstate and $\epsilon_{k}$ is as an energy penalty which has to be chosen large enough, which means larger than the energy difference $|E_{k}-E_{n}|$ to the target state. Of course this energy difference is not known at the start of the algorithm, but in practice one can simply estimate a value which is guaranteed to be large enough, e.g. one can set $\epsilon_{k}$ to be one order of magnitude larger than a typical energy scale in the system. The scaling of the computational complexity of the algorithm is not changed by the additional projective terms in Eq. (9): In complete analogy to the effective Hamiltonian these terms lead to effective projectors, which only contribute a (typically small) overhead to the algorithm. In particular, the runtime scaling for determining excited states remains at $\mathcal{O}(N_{s}\,m^{4})$.

Appendix B Procedure for determining the location of branch vortices

As described in Sec. IV.2, we obtain the branch vortex count of a point $(\theta_{x},\theta_{y})$ on the boundary condition grid by summing up the angle differences $\Delta\Omega$ of the argument field $\Omega_{\Phi\rightarrow\Phi^{\prime}}$ along the four nearest neighbors of $(\theta_{x},\theta_{y})$. This procedure, following Refs. Hatsugai, 2004; *Hatsugai_2005, is exemplified in Fig. 9 for four different illustrative configurations of $\Omega_{\Phi\rightarrow\Phi^{\prime}}$, displaying different branch vorticities. Moreover, in Fig. 9e we provide an illustration for the color code used in Figs. 3d and 8b.

Appendix C Bond dimension convergence of observables

We exemplify the convergence of the energy expectation value as a function of the bond dimension $m$ in Fig. 10. Estimates for finite bond dimension errors can be obtained by extrapolating to the limit $m\rightarrow\infty$. For the system sizes considered here, we reached bond dimensions of up to $m\approx 500$, typically displaying truncation errors of order of $10^{-6}$. In Fig. 11 we demonstrate bond dimension convergence for the Green functions shown in Fig. 4a.