Strongly interacting bosons on a three-leg ladder in the presence of a homogeneous flux

We perform a density-matrix renormalization-group study of strongly interacting bosons on a three-leg ladder in the presence of a homogeneous flux. Focusing on one-third filling, we explore the phase diagram in dependence of the magnetic flux and the inter-leg tunneling strength. We find several phases including a Meissner phase, vortex liquids, a vortex lattice, as well as a staggered-current phase. Moreover, there are regions where the chiral current reverses its direction, both in the Meissner and in the staggered-current phase. While the reversal in the latter case can be ascribed to spontaneous breaking of translational invariance, in the first it stems from an effective flux increase in the rung direction. Interactions are a necessary ingredient to realize either type of chiral-current reversal.


Introduction
Experimental progress with ultracold quantum gases has made feasible engineering the coupling between the different states of the atoms, in order to realize synthetic gauge fields [1,2]. The effective magnetic fields acting on the neutral atoms can be much larger than what is possible in solid-state systems. These advances bring the simulation of a wide range of Hamiltonians into reach that are important in condensed matter physics [3,4,5,6,7,8]. Indeed, some of the most intriguing phenomena in condensed matter physics involve the presence of strong magnetic fields. For instance, topological states of matter are realized in quantum Hall systems [9,10], which are insulating in the bulk, but bear conducting edge states. Remarkably, topological phase transitions were observed in experiments with cold atoms [11,12,13].
Such N-leg ladder systems, similar to their counterparts in quantum magnetism [49], provide an elegant bridge between the physics in one and two dimensions [50]. Interacting bosons on two-leg ladders in homogeneous magnetic fields harbor very rich physics, featuring Meissner and vortex phases [51,17,18,19,52,53,21,23], a biased-ladder phase that breaks the symmetry between the two legs [54,55,22], analogues of Laughlin states [56] as well as an interaction-driven reversal of the direction of the current due to spontaneous breaking of translational invariance in vortex lattices [22].
Here, we use density-matrix renormalization-group (DMRG) [57,58,59] simulations to explore the phase diagram of strongly interacting bosons on a three-leg ladder subjected to a homogeneous flux. We focus on a filling of one-third of a boson per site, which can easily be realized in experiments [26]. At this filling, the system is expected to be in a Mott-insulating state, based on the existence of a magnetization plateau at one-third of the saturation magnetization in the closely related three-leg spin-1/2 ladders [60,61] and on work on the MI-SF transition in bosonic three-ladders at zero flux a zero flux [62]. We present the phase diagram for hard-core bosons (HCBs) in dependence of the inter-leg coupling strength and the magnitude of the magnetic flux. As a main result, we observe vortex phases that cannot be traced back to features in the single-particle dispersion. In the Meissner phase in this system, there is a fascinating reversal of the direction of the current, driven by the magnetic flux. We provide an explanation for this effect and we show that it persists at  (1). (r, ℓ) labels the site on the r-th rung and on the ℓ-th leg. The tunneling strength along the legs(rungs) is J(J ⊥ ). We choose the gauge such that the phase is picked up on the rungs, resulting in a net flux of φ through each plaquette. The on-site interaction strength is given by U .
intermediate interaction strengths. Finally, a staggered-current phase, generalizing the one found on the two-leg ladder [17,18], is observed around φ ∼ π, which also triggers a chiralcurrent reversal.

Model and method
We consider the Bose-Hubbard Hamiltonian on a three-leg ladder of length L (see figure 1): where a † r,ℓ creates a boson on the r-th rung and the ℓ-th leg of the system and n r,ℓ ≡ a † r,ℓ a r,ℓ is the on-site number operator. J and J ⊥ are the tunneling matrix elements along the legs and rungs, respectively, U is the on-site interaction strength, and φ is the magnetic flux through a plaquette of the ladder.
We compute the ground state of (1) numerically with a single-site DMRG algorithm [63]. While exploiting the U(1) symmetry of the Hamiltonian associated to particle number conservation, we keep up to 4000 DMRG states and the data in the figures of the main text are for L = 100, while we have considered L as large as L = 200 in the appendix. For U/J < ∞, we limit the number of bosons per site to six at maximum. For the determination of the Mott insulator (MI) to superfluid (SF) transitions, we use infinite-size DMRG [64] with up to 600 DMRG states [see Appendix B].
The different phases are primarily characterized by their local current configurations. The associated operators are We compute the chiral current from Throughout the paper, we fix the filling to 1/3 bosons per site.

Summary of results for two-leg ladders
Let us first briefly summarize some results for the corresponding two-leg ladder model in the low density regime (obtained by limiting the number of legs to two in figure 1), which will help to appreciate our results for the three-leg ladder to be presented in the following. The phase diagram for hard-core bosons at a filling of one boson per rung has been reported in [21,23]: at fixed rung coupling J ⊥ /J < 1.6, there is a transition from a phase with Meissner-like currents to a vortex phase when increasing the flux. Above a critical value J ⊥ c ≃ 1.6J, the Meissner phase is stable at any flux. Furthermore, at certain commensurate vortex densities, vortex lattices are expected to form [51], which, however, do not seem to exist in the limit of hard-core bosons [21]. At smaller interaction strength, though, stable vortex lattices at vortex densities of 1/2 and 1/3 have been identified in [22]. In those phases, the translation invariance associated with translations along the legs of the ladder is spontaneously broken and an enlarged unit cell forms. As a consequence, the effective flux sensed by the bosons is modified. This leads to a spontaneous reversal of the direction of the chiral current, under certain conditions on the flux and interaction strength [22].

3.2.
Overview of the phase diagram for the three-leg ladder Figure 2 shows the ground-state phase diagram for HCBs [(a † r,ℓ ) 2 = 0; U/J → ∞] of the three-leg ladder model equation (1) as a function of J ⊥ ∈ [J, 2J] and flux. Results are displayed for φ ∈ [0, π] since all physical quantities are 2π-periodic, odd or even, functions of the flux [22]. We find the following phases: (i) a Meissner phase (M-MI/M-SF), which shows a reversal of the current direction for large values of the flux, (ii) vortex-liquid phases (V), (iii) a vortex-lattice phase (VL), and (iv) staggered-current phases (SC-MI/SC-SF). The corresponding transitions are located by solid lines in Fig. 2, in which the gray '+' symbols indicate the parameter points that were analyzed numerically. For J ⊥ /J 1.6 or φ 0.3π, we observe a finite mass gap [see Appendix B] indicating that the system is in a Mottinsulating state. This is inherited from the J ⊥ /J ≫ 1 limit, in which the ground state is the product of the local ground state on each rung, separated by a gap √ 2J ⊥ from all excited states (similar to the two-leg ladder with one boson per rung [65]). Upon lowering J ⊥ at fixed φ, the mass gap eventually closes. The computation of the Binder cumulant [66] of an appropriate string-order parameter for the Mott insulator yields an upper bound for the transition to a superfluid phase when lowering J ⊥ (see the dotted line in figure 2). This transition is, for φ/π 0.3, compatible with a Berezinskii-Kosterlitz-Thouless transition [see Appendix B] as expected from the theoretical work on magnetization plateaux in N-leg spin ladders [60,61] and numerical work at φ = 0 [62]. The transition at 0.9 < φ/π ≤ 1, on the other hand, is compatible with a second-order phase transition [see Appendix B]. In the following, we refer to states with a Meissner-like current configurations using the label M, and to states with a staggered-current pattern using the label SC (irrespective of whether these are MI or SF phases).

Meissner phases
For J ⊥ /J ≥ 1.6, only a Meissner phase exists, characterized by currents occurring exclusively in the upper and the lower leg [see figures 3(a) and 3(e)], with a finite mass gap. The existence of such a Meissner-Mott insulator (M-MI) on the three-leg ladder has been predicted in reference [56]. For φ/π ≤ 0.3, we find a transition from this Mott-insulator to a superfluid phase with Meissner currents (M-SF) by lowering the inter-chain coupling J ⊥ [dashed line in figure 2]. Intriguingly, the chiral current reverses its chirality from counterclockwise [figure 3(a)] to clockwise [figure 3(e)] for φ/π 0.75, meaning that the atoms flow in the direction opposite to the one favored by the effective magnetic field. As a consequence, at the point at which the reversal occurs, no current flows even though the bosons feel a very strong, non-staggered, magnetic flux. A typical example for the j C = j C (φ) curve in the Meissner

Vortex phases
Upon lowering J ⊥ to J ⊥ /J 1. ]. The transition from the Meissner into the vortex phases is characterized by j R becoming non-zero and by kinks in the j C = j C (φ) curves (see the example of J ⊥ /J = 1.2). The vortex phases can be further divided into vortex-liquid (V) phases, which are incommensurate, gapless phases [see, e.g., figure 3(b)-(c)] and a vortex lattice (VL), which is fully gapped and forms at a commensurate vortex density [51] [see, e.g., figure 3(d)].

Vortex lattice.
The commensurability of the phases can be unveiled by studying the spatial patterns of the rung currents j ⊥ r,ℓ=1 . In general, their Fourier transforms bear two typical wavelengths q A and q B (see the inset of figure 5, we choose |q B | < |q A |). Motivated by this observation, we define two vortex densities ρ A,B v = q A,B a −1 /(2π) (with a the lattice spacing), which are shown in figure 5 for J ⊥ /J = 1.2. In the Meissner and VL phases, both vortex densities are commensurate with ρ α v = 0 and ρ α v = 1/2 (with α = A, B), respectively. These two phases are fully gapped. In the case of the VL, we ascribe this behavior to the incompressibility of the vortex pattern: it costs a finite energy to add a vortex to the system, even in the thermodynamic limit. The transitions to the vortex-liquid phases are continuous commensurate-incommensurate transitions [51]. We therefore expect the vortex-liquid region to surround the VL phase everywhere, even though the proximity of various phase transitions renders it very difficult to resolve numerically.

Vortex liquids.
The incommensurate vortex-liquid region (V) encompasses phases in which both vortex densities are incommensurate [e.g., at φ/π = 0.7 in figure 5] as well as phases where one mode is commensurate (at either ρ α v = 0 or ρ α v = 1/2) and the other is not [e.g., at φ/π = 0.375 and 0.5 in figure 5]. This is fully corroborated by the study of the von Neumann entropy, yielding either a central charge c = 2 or c = 1 [see Appendix A] in the vortex-liquid phases, corresponding to two-and one-component Luttinger liquids, respectively (we do not distinguish between the c = 1 and c = 2 vortex-liquid phase in the figures). In principle, the study of current-current correlations could also permit to distinguish the different vortex phases. However, the superposition of the behaviours of the different components renders the analysis of correlations less conclusive for finite-size systems.
We stress that the emergence of vortex phases for bosons on the three-leg ladder is a many-body effect: the minimum of the single-particle dispersion is always at zero momentum, corresponding to a Meissner phase with counterclockwise chiral-current. The vortex phases are thus not inherited from finite-momentum global minima in the single-particle dispersion relation, unlike in two-leg ladders [11,52,21,22].

Chiral-current reversal in the Meissner phase
Let us further investigate the region φ/π 0.75 where the chirality of the chiral currents is opposed to the one favored by the bare flux. In the strong-rung limit J ⊥ ≫ J, the three-site problem is solved analytically and standard perturbation theory (PT) can be applied. To first order in J/J ⊥ , the chiral current is then given by The flux dependence j C ∝ sin(φ) is directly inherited from that of the two-leg ladder in the Meissner phase [21]. The two terms of (5) can therefore be interpreted as stemming from two contributions, which are sketched in figure 6. The first term (left panel) corresponds to chiral currents flowing on the two sub-ladders (formed by the middle leg and the upper or lower leg) resulting in a zero net-current on the middle leg, whereas the second term (right panel) corresponds to particles propagating only along the top and bottom leg of the threeleg ladder. This can be thought of as a Meissner phase with an effective flux φ eff = 2φ. For + Figure 6. Schematic representation of the two contributions to equation (5) for π/2 < φ < π. π/2 < φ < π, the second contribution is negative, which can lead to a total chiral current with a reversed chirality. This reversal is thus associated with a doubling of the effective flux along the rung direction. It can already be captured in the minimal model of two plaquettes (L = 2), in which a qualitatively similar reversal occurs in the strongly interacting limit. An increase of the effective flux also underlies the chiral-current reversal in the case of two-leg ladders [22], yet there, the increase results from spontaneous breaking of translation symmetry.
The J ⊥ ≫ J limit also permits to understand another numerical observation: the density is imbalanced between the legs, in all phases. PT indeed predicts a density twice as large on the middle leg than on the outer ones. This difference decreases when lowering J ⊥ /J.

Finite interactions U/J < ∞
In the single-particle case (U/J = 0), the Meissner phase has counterclockwise chirality for any choice of parameters. Therefore, there must be a reversal of the chiral current as the interaction strength increases. This is indeed the case as shown in figure 7, displaying j C as a function of U, for φ/π = 0.8 and 0.85 and J ⊥ /J = 1.6. The reversal from a negative to a positive current occurs at a finite interaction strength U whose value depends on φ. For φ/π = 0.85, the system enters into a staggered-current phase at intermediate values of U, as indicated in the figure. Translational invariance is spontaneously broken in the SC phase and the unit cell comprises four plaquettes [see figure 3(f)]. For φ ∈ [3π/4, π], the effective flux is φ eff = 4φ ∈ [−π, 0] modulo 2π, and the current is reversed. This realizes another instance of the chiral-current reversal due to spontaneously enlarged unit cells first presented in reference [22]. The presence of the SC phases stabilizes the current reversal down to much smaller values of U/J compared to parameters for which the SC phase is absent (compare the data for φ/π = 0.8 and 0.85 shown in figure 7). Note that the VL phase leads to the same unit-cell enlargement. However, as it occurs for φ ∈ [π/2, 3π/4] in the HCBs case (see figure 2), one has φ eff = 4φ ∈ [0, π] modulo 2π, and the chiral current is not reversed.

Experimental realizations
The most straightforward experimental approach would be to use a superlattice to split a twodimensional lattice into three-leg ladders [67]. This would create an energy offset between the middle and the outer two legs, which could easily be compensated by another superlattice with a two-site periodicity. For systems with a synthetic lattice dimension [26,25,38], the interaction in the rung direction is not on-site as considered here, but depends on the total density n t r = 3 ℓ=1 n r,ℓ in all three rung sites. Thus, the appropriate form is Our DMRG results for U/J = 10 show that this type of long-range interactions suppresses vortex phases. This is qualitatively consistent with the study of the two-leg ladder, in which this type of interaction also favors the Meissner phase over the vortex phase [19,56]. The chiral-current reversal in the Meissner phase at a filling of one-third survives at small J ⊥ /J 0.3 [see figure 8].

Summary
We presented the phase diagram of strongly-interacting bosons on a three-leg ladder subjected to a homogeneous magnetic flux. We identified several phases, including vortex-liquid phases, a vortex lattice and a Meissner phase. Moreover, there is a state with staggered currents, which leads to a reversal of the chiral current due to the spontaneous increase of the unit cell, similar to the situation discussed in [22]. Fascinatingly, the Meissner phase also shows a chiralcurrent reversal when increasing the strength of the magnetic flux per plaquette, yet in the Meissner phase, translational invariance is not broken. The analysis of the strong-rung and large-interaction limit indicates that a doubled flux is experienced along the rung direction. Therefore, the chiral-current reversal in the Meissner phase is qualitatively different. We argue that interactions are a necessary ingredient to obtain these behaviors.

Acknowledgments
We

Appendix A. Entanglement entropy
The entanglement entropy S vN is defined as the von Neumann entropy corresponding to a bipartition of the wavefunction into two subsystems A and B: where ρ A = Tr B |ψ ψ| is the reduced density matrix of subsystem A. In conformal field theory, for a one-dimensional system with open boundary conditions and total length L, the von Neumann entropy of a subsystem of length r is given by [68,69] S vN (r) = c 6 log 2L π sin πr L + g, where c is the central charge of the conformal field theory and g is a non-universal constant.
The central charge measures the number of gapless modes in the system. In figure A1 we show the entanglement entropy S vN (r) for straight cuts through the three-leg ladder between rungs r and r + 1. By fitting (A.2) to the data we obtain estimates for the central charge. Oscillations of the entropy in the vortex phases make the fit very sensitive to the fitting region used. This applies, in particular, to the vortex-liquid phases. Nevertheless, the behavior of the central charge provides further insights into the nature of the commensurate-incommensurate transition. In the Meissner phase [figures A1(a) and (e)] the values of the central charge are very close to zero in agreement with an entanglement entropy independent of the subsystem size, which corresponds to the area law of a fully gapped phase [70].
The vortex phases can be divided into vortex-liquid and vortex-lattice phases with incommensurate and commensurate vortex densities, respectively. As illustrated in figure 5 of the main text there are two vortex densities that undergo a transition from the Meissner to a vortex-liquid and vortex-lattice phase with increasing flux φ/π and at low enough J ⊥ /J. This picture is consistent with the behavior of the central charge. More precisely, at J ⊥ /J = 1.2 and φ/π = 0.375, we find one incommensurate vortex density in the Fourier transform of the rung currents and a central charge of c ≈ 1 [ figure A1(b)]. At flux φ/π = 0.7, there are two incommensurate vortex densities. The central charge is c ≈ 2 in agreement with the presence of two gapless modes [ figure A1(d)]. At φ/π = 0.6, both vortex densities are commensurate and we find an entanglement entropy independent of subsystem size (up to oscillations due to the increased unit cell), or c = 0 in terms of the central charge [ figure A1(c)].
This suggests that the modes, which are measured by the central charge, can be related to the dominating vortex densities in the Fourier transform of the rung currents. An incommensurate vortex density would be associated with a gapless mode while a commensurate vortex density would correspond to a gapped mode.
In the staggered-current phase the vortex densities are commensurate, while the fit gives a central charge c ≈ 1. But the scaling of the entanglement spectrum indicates that it is two independent c = 0.5 modes, and also the scaling of the correlation function gives the scaling dimension ∆ ≃ 1/8 (not shown), consistent with c = 0.5. This indicates that the mass gap is closed, which is compatible with the analysis of the mass gap and the Binder cumulant presented in the next section.

Appendix B. Mott insulator to superfluid transition
The mass gap ∆E M is defined as where E(N) is the total ground-state energy of the system with N particles. We generally observe a finite mass gap that becomes smaller when decreasing the interleg coupling J ⊥ . In the J ⊥ = 0 limit, the mass gap is known to be closed. As already mentioned in the main text, in the φ/π = 0 limit, the transition to the superfluid phase at finite J ⊥ /J can be expected to be of the Berezinskii-Kosterlitz-Thouless type [60,61], and, therefore, is very hard to pinpoint numerically.
In the insets of figure A1 we show the mass gap as a function of the inverse system size 1/L extrapolated with a linear fit to the thermodynamic limit. We find a finite mass gap in the Meissner phase for J/J ⊥ = 2, φ/π = 0.2 and J/J ⊥ = 1. The analysis of the central charge is more difficult in the regions where the mass gap can not be resolved (or is actually zero) with the system sizes studied in this work. This is the case in the Meissner phase for J ⊥ /J 1.6 and φ/π 0.3 and in the staggered-current phase [see the inset of figure A1(f)].
Calculating energy gaps is a quite difficult calculation when the gap is small, so where an order parameter exists it is generally better to use it instead to detect the phase boundary of the Mott-insulating region. A suitable order parameter for a Mott phase is the non-local string order parameter [71], given by where n t r = 3 ℓ=1 n r,ℓ is the total particle density on rung r. This can be expressed in terms of a 'kink' operator,  or, we can construct an extensive order parameter P = i p(i). We cannot evaluate P directly, since the sign is indeterminate, however, we can calculate P 2 , which is simply related to O 2 p , as where L is the length of the lattice. The great advantage of this construction is that the P operator is easy to construct for matrix-product-state (MPS) algorithms using the triangular matrix-product-operator (MPO) formulation [72,64], and the expectation value of such triangular MPO's can be evaluated directly in the thermodynamic limit [73]. This allows us to calculate the Binder cumulant [66] of P , For an infinite matrix-product-state (iMPS), the expectation value of the n-th power of an operator P n is obtained as a degree n polynomial of the lattice size L, which is exact in the asymptotic large-L limit. Hence the quantity O 2 p = P 2 /L 2 is evaluated directly as the coefficient of the degree 2 component of P 2 .
If we evaluate the Binder cumulant directly for an iMPS in the large L limit, then it simply probes where O 2 p is zero or non-zero, so does not give any additional information. Since the coefficient of the n-th degree term in P n is simply O n p , it follows that if O 2 p is non-zero then the Binder cumulant in the large L limit is identically equal to 2/3. Similarly, if O 2 p = 0, then the cumulant expansion of P 4 reduces to 3 P 2 2 and the Binder cumulant is identically zero. Hence we have a step function located at the point where O 2 p becomes non-zero, which is what one would expect if one takes the large-L limit of a finite-size Binder cumulant.
However, for an iMPS the correlation length is always finite, so one might expect that finite-entanglement-scaling with the basis size m can be used instead of finite-size scaling. This is indeed the case, if we evaluate the polynomial P n at L = bξ, where ξ is the correlation length, and b is some scaling factor. This is equivalent to calculating the order parameter over a finite section of size L = bξ of the infinite lattice.
Since ξ depends on the number of states m, we can plot a family of curves of the Binder cumulant B m for different values of m. For a finite system, the value of the Binder cumulant at a second-order critical point is independent of L (up to higher order corrections). For an iMPS, the analogous result is that the critical value of the Binder cumulant is independent of m. Note that for finite size systems the actual value of the Binder cumulant at the critical point is not universal and depends on the boundary conditions [74]. For an iMPS, the value of the m-dependent Binder cumulant depends on the chosen scale factor b. For the calculations here, we use b = 1.
As an example, we show in figure B1(a) the Binder cumulant as a function of J ⊥ for φ/π = 0.2. The data is compatible with a Berezinskii-Kosterlitz-Thouless transition from the Meissner Mott-insulating region into the Meissner superfluid region as J ⊥ /J decreases and permits to give an upper bound for the critical point, which is reported in figure 2 in the main text. For φ/π = 0.95 [see figure B1(a)], a Mott-insulator to superfluid transition also occurs in the SC phase and it is compatible with a second-order critical point at J ⊥ /J ≃ 1.035. The Meissner to SC phase transition, which can be characterized by a rise of the average rung-current j R , occurs at a distinct point J ⊥ /J ≃ 1.135.