{'fraction_non_alphanumeric': 0.04422483791779442, 'fraction_numerical': 0.017083149745074588, 'mean_word_length': 4.278490265133896, 'pattern_counts': {'":': 0, '<': 34, '<?xml version=': 0, '>': 15, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Symmetry-broken states in a system of interacting bosons on a two-leg ladder with a uniform Abelian gauge field

We study the quantum phases of bosons with repulsive contact interactions on a two-leg ladder in the presence of a uniform Abelian gauge field. The model realizes many interesting states, including Meissner phases, vortex-fluids, vortex-lattices, charge-density-waves and the biased-ladder phase. Our work focuses on the subset of these states that break a discrete symmetry. We use density matrix renormalization group simulations to demonstrate the existence of three vortex-lattice states at different vortex densities and we characterize the phase transitions from these phases into neighboring states. Furthermore, we provide an intuitive explanation of the chiral-current reversal effect that is tied to some of these vortex lattices. We also study a charge-density-wave state that exists at 1/4 particle filling at large interaction strengths and flux values close to half a flux quantum. By changing the system parameters, this state can transition into a completely gapped vortex-lattice Mott-insulating state. We elucidate the stability of these phases against nearest-neighbor interactions on the rungs of the ladder relevant for experimental realizations with a synthetic lattice dimension. A charge-density-wave state at 1/3 particle filling can be stabilized for flux values close to half a flux-quantum and for very strong on-site interactions in the presence of strong repulsion on the rungs. Finally, we analytically describe the emergence of these phases in the low-density regime, and, in particular, we obtain the boundaries of the biased-ladder phase, i.e., the phase that features a density imbalance between the legs. We make contact to recent quantum-gas experiments that realized related models and discuss signatures of these quantum states in experimentally accessible observables.


I. INTRODUCTION
An important part of the physics of quantum particles moving in a two-dimensional plane under the action of a magnetic field is related to the quantum Hall effect [1,2]. While the essence of the integer quantum Hall effect can be understood from considering noninteracting electrons, interacting particles in two dimensions and in the presence of an Abelian gauge field provide an ideal playground to explore exotic many-body physics [3][4][5], encompassing, most notably, the fractional quantum Hall effect. Periodic lattice potentials introduce additional intriguing physics. A magnetic field applied perpendicular to the plane of motion of a charged particle in a lattice produces a fascinating structure of energy levels known as the Hofstadter butterfly [6] resolved experimentally in solid state systems only recently [7][8][9].
The experimental progress in the field of ultracold atomic gases with emulating gauge fields or spin-orbit coupling in these systems of neutral particles has opened new prospects for observing many-body physics in the presence of gauge fields in a very clean and highly tunable environment [10][11][12]. While pioneering experiments have demonstrated the successful implementation of spinorbit coupling in Bose gases in the continuum [13,14], the field has seen a tremendous activity in studying the combined effects of optical lattices with artificially engi-neered gauge fields [15][16][17][18][19], accessing the physics of the Hofstadter model [18,19] as well as the famous Haldane model [20,21].
So far, quantum gas experiments have focussed on non or weakly-interacting quantum gases. In this regime, some key hallmark features of topological states of matter [22,23] were measured, such as the Zak phase in one-dimensional systems [24], the Berry curvature in a Floquet system [25], the Chern number [21,26] or the Berry flux in momentum space [27].
Our work is primarily motivated by the experimental realization of ladder systems combined with uniform Abelian gauge fields in ultracold atomic gases. Ladders, which here we have in mind to be arrays of plaquettes as indicated in Fig. 1, are the simplest possible extensive lattices that allow one to study nontrivial orbital effects in the presence of a synthetic magnetic field. Such ladder structures can either be obtained using superlattices or a so-called synthetic lattice dimension [43]. The former approach has been utilized in [44] to study a weaklyinteracting Bose gas loaded into an optical lattice with the two-leg ladder geometry. This experiment heavily relied on the measurement of local currents and could access the chiral edge current predicted to exist in this system in the presence of a uniform flux per unit cell [45,46]. Moreover, the experiment established a similarity to superconductors (also previously discussed in theoretical papers [45,46]), since the chiral current at small flux behaves similar to the screening current in the Meissner phase (with no current in the bulk, which, for a ladder, implies a vanishing of local currents on its rungs) while at large fluxes, finite rung currents emerge, reminiscent of a vortex phase in type-II superconductors.
The synthetic lattice dimension approach combines an actual optical lattice, which is typically one-dimensional, with Raman lasers that drive transitions between a subset of the hyperfine states of bosonic or fermionic atoms [43]. In the first two implementations of this scheme, two-and three-leg ladders have been realized using this (a) Sketch of the two-leg ladder model and interaction-and tunneling terms as defined in Eqs. (1), (3) and (19). Current patterns and onsite density in (b)-(d) the different vortex lattices with vortex density (b) ρv = 1/2, (c) ρv = 1/3, and (d) ρv = 1/4, (e) the biased-ladder phase (BLP) phase, and (f) in the charge-density-wave (CDW 1/4 ) phase at particle density ρ = 1/4. The arrows indicate the direction and by their length, the strength of local currents. The density is represented by the size of the circles and the background shading. method [47,48], and both experiments also succeeded in measuring edge currents (more specifically momentum distributions), for fermions [47] and bosons [48]. More recent experiments realized two-leg ladders using an optical clock transition [49,50] (based on the proposal from [51]) or even using an all-synthetic-lattice approach [52].
These examples demonstrate the experimental study of the effects of synthetic gauge fields in neutral ultracold atomic gases in low-dimensional lattice geometries constitutes an active field of timely research. Despite the enormous recent progress with experimental observations of fascinating single-particle behavior due to synthetic gauge fields, the strongly-interacting regime remains hard to access. One possible reason that is being investigated is heating in driven many-body systems [53][54][55][56] due to the various driving schemes used to mimic the artificial gauge fields [16,18,19,21].
The interest in bosonic ladders dates back to the development of the theory of arrays of Josephson junctions [57][58][59][60], adequately described by high density and weak interactions. Subsequently, a first bosonization study for weakly-coupled legs of a bosonic two-leg ladder explored the strongly-interacting regime, predicting the stability of Meissner-like, vortex fluids and vortex lattices [45], whose existence was initially established for the weaklyinteracting regime only [57][58][59][60]. Then, fueled by the experimental progress with emulating artificial gauge fields, Dhar et al. [61,62] demonstrated the existence of socalled chiral Mott insulators in bosonic two-leg ladders at flux φ = π per plaquette (see the sketch of the model in Fig. 1(a)), which spontaneously break time-reversal symmetry. A bosonization study of the orbital response of a bosonic ladder in the strongly-interacting regime at arbitrary flux was presented in [63], predicting Meissner and vortex phases also in the Mott-insulating case. Understanding the orbital response of interacting fermions on two-leg ladders has been the topic of Refs. [64,65].
Other efforts addressed the fate of the Meissner and vortex-like states inherited from the noninteracting limit of bosonic two-leg ladders. Their origin can be traced back to the existence of either a unique minimum in the single-particle dispersion at quasi-momentum or two equivalent minima at incommensurate momenta, the former corresponding to the Meissner, the latter to the vortex phase [44,46]. Similar states survive even in the limit of strongly-interacting bosons [63,66,67], both on top of Mott-insulators and superfluids. Overall, interactions lead to a suppression of vortex phases.
While most of the theoretical work has investigated the case of pure contact interactions of the type H int = U/2 i n i (n i −1), where n i measures the density in a site of a lattice (U is the interaction strength), in the experimental realizations using a synthetic lattice dimension [47,48], atoms experience long-range interactions. The reason is that the physical interaction is of contact type, meaning that all atoms in a rung interact with each other, regardless of their hyperfine state. The effects of longerrange interactions have been studied in [70,73,[78][79][80][81][82][83][84][85][86].
Returning now to the case of bosonic two-leg ladders with contact interactions in the model, while Meissnerlike and vortex-liquid phases exist in the presence of strong interactions, vortex lattices, on the contrary, have remained elusive in the strongly-interacting regime, with the exception of the chiral Mott insulators [61,62], which in fact represents a vortex lattice state with a maximal possible vortex density ρ v = 1/2, and a unit cell of two plaquettes as shown in [87]. In Ref. [87], we determined the range of stability of the above mentioned vortex lattice at vortex densities of ρ v = 1/2 and of another vortex lattice state ρ v = 1/3 (having a unit cell made of three plaquettes) in the regime of intermediate interactions 1 < U/J < 10 and a low filling ρ ≤ 1. Moreover, we discovered that interactions can lead to a spontaneous reversal of the circulation direction of the boundary chiral current in certain vortex-lattice states: the atoms there behave as if the direction of the external magnetic field (or the sign of the flux) had been inverted. This effect can be understood as resulting from the periodicity of the chiral current with 2π, the increase of the effective flux from φ to qφ in vortex lattices with unit cells of q plaquettes and the quantum nature of atoms. Only those vortex lattices that are stable at values of the flux for which π < qφ < 2π lead to this reversal. Remarkably, this reversal is stable against temperatures that are possible to realize in current experiments and hence it can be used as an experimental probe of the existence of vortex lattices. In bosonic ladders, such a reversal of the chiral current occurs only in the interacting case, while for fermions, in a complicated band structure that results from adding the flux, changing the filling can lead to a chiral-current reversal already in noninteracting systems simply due to the Pauli principle [64].
The purpose of the present work is to provide a comprehensive analysis of quantum phases in bosonic two-leg ladders subject to a uniform flux φ per plaquette that spontaneously break some discrete symmetry of the microscopic model. We primarily focus on the low-density regime ρ ≤ 1 and study such states as a function of interaction strength, flux, and the ratio J ⊥ /J of hopping matrix elements along rungs J ⊥ and legs J (see the sketch of the model shown in Fig. 1(a)).
Most notably, the list of states with broken discrete symmetries includes vortex lattices, which break lattice translation invariance. We study the previously known vortex lattices at vortex densities ρ v = 1/2 [61,62,87] and 1/3 [87] and we report numerical evidence for the existence of an additional vortex lattice at ρ v = 1/4. Typical configurations for the local particle currents and densities in these three vortex lattice states are plotted in Figs. 1(b)-(d). The position of the vortex cores is denoted by the symbol V in the figure. While the screening current in a Meissner phase goes counterclockwise around the boundary of the system, the vortices carry currents of the opposite chirality, thus reducing the overall chiral current. In addition, the vortex lattices with ρ v < 1/2 feature density modulations that are locked to the structure of the vortex lattice. The results shown in Fig. 1(b)-(f) were obtained from density matrix renormalization group (DMRG) simulations [88][89][90], the primary tool in our analysis.
Our main results beyond those of Ref. [87] for the vortex lattices are the discovery of a stable ρ v = 1/4 vortex lattice, the stability analysis of the ρ v = 1/3 against increasing interactions, the analysis of the phase transitions between the vortex lattices and the neighboring phases, and an intuitive discussion of the chiral-current reversal that we develop by comparing systems with spontaneously enlarged unit cells to systems with explicitly larger unit cells. Moreover, we demonstrate that nearestneighbor interactions on the rungs of a two-leg ladder suppress vortex phases in favor of the Meissner state. From a conceptual point of view, it is very important to compare different measures of the vortex density, which does not have a microscopic definition. We also discuss the signatures of vortex lattices in the experimentally accessible quasi-momentum distribution function, comparing different gauges.
Another state that breaks a discrete symmetry is the biased-ladder phase (BLP), first discussed by Wei and Mueller [75] and also studied in [76,87,91]. In this state, the density between the two legs of the ladder is imbalanced, which serves as an order parameter. Thus, this state breaks the Z 2 symmetry associated with inversion of the two legs together with inversion of the sign of flux. A typical configuration of the local currents and density in the BLP state is shown in Fig. 1(e). We determined the phase boundaries of the BLP phase at intermediate interaction strength in [87] from accurate DMRG simulations, thus providing robust evidence for its existence beyond mean-field [75] and bosonization [76,91] predictions. In our present work we also introduce a theory of the emergence of the Meissner state, vortex fluids and the biased-ladder phase based on the limit of a dilute Bose gas that some of us originally developed for the description of frustrated spin chains just below their saturation magnetization [92,93].
Finally, we discuss charge-density-wave (CDW ρ ) states. CDW 1/4 states exist, in particular, in the limit of hard-core bosons at a filling of ρ = 1/4 and for sufficiently large values of the hopping on the rungs J ⊥ > J [67,71]. In Ref. [67], we provided a theoretical explanation for their existence based on a mapping to an effective spin-1/2 Hamiltonian valid in the regime of J ⊥ > J. In the present work, we primarily focus on the stability of the CDW 1/4 state against going to lower values of c ρv q avg|jR| ∆n Quantum phases of bosons with repulsive contact interactions on a two-leg ladder with a uniform Abelian gauge field studied in this work. We list those states that we actually detected in our numerical simulations with no claim of exclusiveness as additional states exist, which we do not explicitly discuss in this work. Meissner, vortex liquid and vortex lattice phases exist either atop superfluid (SF) or Mottinsulating (MI) states. We also list characteristic properties (see the text for details) such as the central charge c, counting the number of gapless modes, the vortex density ρv, the size of the effective unit cell of the groundstate q (plaquettes), the average local rung current in the thermodynamic limit avg|jR| (see Eq. (7)) and the leg-density imbalance ∆n (see Eq. (27)). The "shaded" states break a discrete symmetry and are at the main focus of this study.
U/J < ∞. As a result, we find that the CDW 1/4 state survives down to U/J 30 at quarter filling ρ = 1/4. At smaller values of U/J, the system transitions from the CDW 1/4 state into a Mott insulator that carries a vortex lattice with ρ v = 1/2. The existence of this Mott insulator is very interesting since it results from the combined effects of interactions, flux and filling, unlike other Mott insulators that exist in the bosonic two-leg ladder [61,63,66,67,94] that can be traced back to the limit of φ = 0. A typical configuration of the local currents and density in the CDW 1/4 state is shown in Fig. 1(f). For strong interactions in the synthetic dimension U/J → ∞ and V /J → ∞ we observe the stability of a CDW 1/3 state, i.e., at filling ρ = 1/3. Table I provides an overview over the quantum phases that are realized in the bosonic two-leg ladder with onsite interactions and some of their characteristic properties. These include the Meissner superfluid (M-SF), the Meissner-Mott insulator (M-MI), vortex liquid superfluids (V-SF), and vortex liquids on top of Mott insulators (V-MI), vortex-lattice superfluids at a vortex density ρ v (V ρv -SF) as well as the VL 1/2 -MI that sits on top of a Mott-insulating state, the charge-density-wave state (CDW 1/4 ) and the biased-ladder superfluid phase (BLP). Other states that have also been proposed to exist in this model [70,77] are not at the main focus of the work and are thus not included in the table. The terminology for these phases as well as the acronyms already suggest the existence of transitions in two sectors: the charge sector, in which the Mott-insulator-to-superfluid transition takes place and the antisymmetric sector, in which the Meissner-to-vortex-to-vortex-lattice transitions occur. A classification of these states can be obtained from computing the central charge c or certain order parameters for the vortex lattice (namely the vortex density ρ v or the average value of local currents j R on rungs) or the density imbalance ∆n between the two legs that is nonzero in the BLP phase. We also list the number q of elementary four-site plaquettes that the unit cell in a given phase contains.
The plan of the paper is the following. We first introduce the model in two different gauge conventions as well as some key observables in Sec. II. Details on our main method, the DMRG technique, are provided in Sec. III. Section IV discusses the various vortex lattices that are stable at low particle densities as well as the transitions into neighboring phases. In Sec. V, we provide an intuitive explanation of the chiral-current reversal that is tied to certain vortex lattices and study the effect of temperature. Section VI is devoted to the CDW 1/4 state that exists at filling ρ = 1/4. The BLP state is studied in Sec. VII, both analytically and numerically. We summarize our findings in Sec. VIII. An Appendix contains our results for additional incommensurabilities in the Meissner phase and a discussion of their possible interpretation.

II. MODEL AND CHIRAL CURRENT
The system is described by the following Hamiltonian with the matrix elements corresponding to hopping along the rungs and legs of the ladder J ⊥ and J, respectively. U is the strength of the onsite interaction (we consider repulsive interactions in this work unless stated otherwise). b † ℓ,r creates a particle in the r-th site on the leg ℓ = 1, 2 and n ℓ,r = b † ℓ,r b ℓ,r . The total number of bosons is denoted by N , while the number of sites in each leg is L (i.e., 1 ≤ r ≤ L). We define the particle filling as ρ = N/(2L). In the following, we set J = 1 as the unit of energy ( = 1).
The model exhibits a gauge freedom in choosing different distributions of the Peierls phases as long as the total flux per ladder plaquette remains invariant. For instance, by means of a unitary transformation to new bosonic operators we can make the hopping matrix elements along the rungs real, but instead the hopping matrix elements along the legs become complex. The Hamiltonian (1) is then given by We shall refer to these two gauge choices described by the Hamiltonians Eqs. (1) and (3) as rung gauge and leg gauge, respectively. An important, experimentally accessible [44,47,48] observable in ladder systems is the local current, defined on either bonds in the bulk of the system or on the boundaries. From the continuity equation dnr dt = i [H, n r ] = − s j(r → s) we can define the current j(r → s) from a site r to a neighboring site s (where r = (ℓ, r)). In particular, for the model (1), we obtain the local currents on legs and rungs from Apart from the configuration of local currents, the average current that circulates along the boundary of the system may reveal important properties of the quantum phases. This so-called chiral current (also dubbed edge, screening or Meissner current) is defined as For a two-leg ladder, we may obtain j c from the Hellmann-Feynman theorem as the derivative of the ground-state energy E 0 per particle with respect to the flux φ This, in particular, shows that the expectation value of the chiral current is gauge invariant. The experimental realizations of ladder models with either superlattices [44] or synthetic lattice dimensions [47,48] correspond to the gauge choice with complex hopping matrix elements along the rungs as in Eq. (1). The superlattice experiment [44] used flux values of φ ≃ ±π/2, while the synthetic lattice dimension experiments were operated at φ ≃ ±2π/3 [48] and φ ≃ ±0.37π [47].

III. NUMERICAL METHOD: DENSITY MATRIX RENORMALIZATION GROUP TECHNIQUE
Most results of this work are obtained by means of large scale numerical density matrix renormalization group [88][89][90] simulations (DMRG), which is a standard method for the simulation of one-dimensional chains or ladder-like systems at zero temperature. We simulate the ladder model of Eq. (3) with up to L = 160 rungs, typically using 1000 DMRG states. We control the accuracy by enforcing a sufficiently small discarded weight.
The repulsive interactions allow us to employ a cutoff for the occupation of bosons per site to address large system sizes. We typically use a cutoff of n max = 4 bosons for U J, n max = 3 for U 10J and n max = 2 for U 30J and fillings ρ < 1. By comparison with larger and smaller cutoffs we have ensured the independence of the numerical data on the cutoff, for the quantities shown in this work. A detailed analysis of the dependence on n max is contained in the supplemental material of Ref. [87].
Close to the V-SF to VL ρv -SF boundaries, the DMRG simulations tend to converge to metastable excited states with a varying vortex density. To overcome this problem, we perform several calculations (sometimes, up to sixteen runs) starting from different randomly chosen initial states. Selecting those states with the lowest energy gives the piecewise continuous results for, e.g., the chiral current j c shown in Fig. 10.

IV. VORTEX LATTICES
The existence of vortex lattices in the bosonic twoleg ladder in the presence of a uniform gauge field was initially predicted from the consideration of largecapacitance Josephson-junction arrays in the so-called classical limit [57]. Since vortex lattices break a discrete symmetry of the model, they are robust to finite quantum fluctuations, as confirmed by a bosonization analysis of the J ⊥ ≪ J regime [45]. For the strongly interacting, low-density regime and arbitrary J ⊥ /J, vortex lattices were first seen for filling ρ = 1 and φ = π [61,62]. In those studies [61,62], the emphasis was put on the spontaneous breaking of time-reversal symmetry and thus this state was dubbed a chiral Mott insulator [61,62] (see also [95] for a discussion of chiral Mott insulators in two dimensions). In this state, translation symmetry is also broken spontaneously, which is, however, not independent from time-reversal symmetry breaking, since translation with respect to one ladder plaquette accompanied with time reversal remains intact. Hence, one can talk about breaking of translation symmetry instead of timereversal symmetry in this state, interpreting it as a usual vortex lattice state packed with the maximal number of vortices (siting on every other plaquette) [87].
Our DMRG study of the hard-core boson limit U/J = ∞ did not observe any vortex lattice (at any density ρ) but merely Meissner and vortex-liquid states [67]. Vortex-lattice states are, at low densities ρ ≤ 1, stable at intermediate interactions as we demonstrated in Ref. [87]: there, we reported evidence for the existence of vortex lattices at ρ v = 1/2 and ρ v = 1/3 for inter-action strength 1 U/J ≤ 10. While detailed phase diagrams for such intermediate values of U/J and low densities have been reported and discussed in Ref. [87], we here focus on the properties of vortex lattices at particle fillings ρ = 0.8 and ρ = 0.5 and the transitions between these vortex lattices and other quantum phases (see Secs. IV C and IV D). Moreover, we have detected another stable vortex lattice at ρ v = 1/4 (see Sec. IV F). This vortex lattice is interesting since with such a large unit cell, more information on the location of vortices and their extension can be extracted.
In our discussion, we devote particular attention to various measures of the vortex density ρ v . For instance, one can extract ρ v from the Fourier transform of rungcurrent patterns (in the case of open boundary conditions), the momentum distribution function (at least in the leg gauge as defined above), or from modulations in the particle density (see Secs. IV A and IV B). While these three measures yield consistent results in vortex lattices and (most) vortex liquids, the outcome can differ in, for instance, the Meissner phase for certain parameter regimes. Finally, in Sec. IV G we discuss the stability of vortex lattices against augmenting the Hamiltonian with nearest-neighbor repulsive interactions on the rungs (as appropriate for synthetic-lattice dimension realizations).
A. Overview: quantum phases and vortex lattices at density ρ = 0.8 Figure 2 depicts various observables for the characterization of vortex-lattice and vortex-liquid phases in the model given in Eq. (1) as a function of flux φ. The numerical results were obtained from DMRG simulations for U = 2J and a low filling ρ = 0.8. A sequence of phases starting from the M-SF, a first vortex-liquid, the VL 1/3 -SF, another sliver of the V-SF, and finally, the VL 1/2 -SF is realized for these parameters.
In the Meissner and vortex-lattice phases, the chiral current j c exhibits a characteristic quasi-linear increase with the flux φ as shown in Fig. 2(a). The average rung current (also plotted in Fig. 2 exhibits a stable large plateau in the vortex-lattice phases. The transition from the vortex-lattice phases to the V-SF phases is indicated by a marked drop of avg|j R |. Apart from identifying the vortex-lattice and Meissner phases by their characteristic local current configurations (see the discussion in Ref. [87] and Fig. 1), they may clearly be discriminated from the vortex-liquid phase by calculating the central charge c, which may be extracted from scaling properties of the entanglement entropy S vN (l) of a a subsystem of length l embedded in a where we have omitted non-universal constants and higher-order oscillatory terms due to the finite size of the system. For a more detailed discussion of the behavior of the entanglement entropy in this model and the extrac-tion of c, see the supplemental material of [67,87]. In praxis, we compute S vN for blocks that contain r rungs, i.e., we discard blocks that would cut a rung and thus we plot S vN versus the number of rungs r in the block in the figures.
In Fig. 2 (b), we depict the extracted central charge from fitting Eq. (8) to the numerical data, which is well consistent with c = 1 in the vortex lattices atop the superfluid phase (such as the VL 1/2 -SF and VL 1/3 -SF states) and the M-SF phase and c = 2 in the V-SF phase. Interestingly, as shown in the inset of Fig. 2(b), the entanglement entropy in systems with open boundaries exhibits small oscillations that follow the lattice structure of the vortex-lattice phases.
We next estimate the vortex density ρ v by analyzing the rung-current configurations j ⊥ r . We follow our previous analysis detailed in [67], where we introduced a measure for the vortex density, given by the inverse typical distance between the vortex cores l v We extract this distance from the Fourier transform of the real-space patterns of the rung currents j ⊥ r (which in the vortex fluids are obviously discernible due to finitesize effects for open boundary conditions). The resulting vortex density depicted in Fig. 2(c) shows the typical devil's staircase-like structure predicted in Ref. [45,57], i.e., a sequence of the Meissner phase (ρ v = 0), incommensurate vortex liquids and vortex lattices. We draw the reader's attention to the sharp jump of the vortex density at the boundary of the M-SF to the neighboring vortex-liquid phase located in the vicinity of φ = 0.7π. This jump may be indicative of a first-order transition.
Alternatively, the vortex density may be related to local density fluctuations in systems with open boundaries (as already discussed by Ref. [44] for noninteracting particles). The Fourier transform of such local density fluctuations exhibits a sharp peak at a momentum k n max , which, outside of the VL 1/2 -SF and Meissner phases coincides with ρ v (see Fig. 2(c)). In the VL 1/2 -SF, there are no density fluctuations by symmetry [87] and hence k n max = 0. Interestingly, in the Meissner phase and for a finite flux, we also observe finite, small-amplitude fluctuations in the density as a precursor of the transition to the vortex phases. As we shall see below in Sec. IV B, in the superfluid phase, where single-particle correlation functions decay algebraically, we can also estimate the vortex density from the momentum distribution function. We will continue our discussion of the vortex density in Sec. IV B.

B. Momentum distribution function: Experimental observable and a possible measure of vortex density
In the following, we study the momentum distributions along the legs of the ladder n ℓ (k), with ℓ = 1, 2: which are measurable in time-of-flight experiments [98]. Using Eq. (2) one can see that the momentum distributions in the two different gauges defined in Sec. II are related to each other via Figure 3(a) shows that, in the leg gauge, the central peak of n ℓ (k) perfectly coincides with the vortex density ρ v of the system (compare Fig. 2(c)). For the rung gauge, the peak position exhibits an additional shift and behaves qualitatively similarly to the chiral current j c (compare Fig. 2(a)). We next compute the momentum distributions in Meissner, vortex-lattice and vortex-liquid phases with the help of an effective field-theory approach based on bosonization. We introduce two pairs of bosonic fields (θ ℓ , φ ℓ ), describing the phase and density fluctuations of bosons on leg ℓ, respectively, with . The low-energy properties of the model Eq. (1) are governed by the following Hamiltonian density The expression for θ − depends on the gauge: for the case where the Peierls phases are along the rungs, θ − = (θ 1 − θ 2 −φx/ √ π)/ √ 2, whereas for the gauge where the Peierls phases are along the legs, θ − = (θ 1 − θ 2 )/ √ 2. K ± are Luttinger-liquid parameters corresponding to the total and relative fluctuations on the two-leg ladder and v ± are the corresponding velocities, which in general need to be determined from a comparison with numerics.
Since the flux couples to the topological charge of the sine-Gordon model describing the antisymmetric sector, we explicitly separate out the zero-momentum mode in the field expansion which is related to the vortex density [57] For small values of φ, the most important term in Eq. (12) is the one proportional to λ 1,0 ∼ J ⊥ , which at any filling and U opens a gap in the antisymmetric sector already for an arbitrarily small interchain tunneling and pins θ − (i.e., it locks the relative phase of bosons on the two legs). Thus, the system is in the Meissner phase as long as φ < φ c , where φ c is a soliton gap of the quantum sine-Gordon model describing the antisymmetric sector [57].
With increasing flux, the terms proportional to λ q,0 with q > 1 can become commensurate (for integer values of qρ v ) and relevant (for q 2 ≤ 4K − ) in vortex-lattice states with a q-fold degenerate ground state and a spontaneously enlarged unit cell with q plaquettes [45].
Using the representation (where the +(−) sign corresponds to ℓ = 1(2)) we can calculate the momentum distributions, say along the first leg of the ladder in vortex-lattice phases (the expression also applies to the Meissner phase where ρ v = 0), resulting in On the other hand, in the vortex-liquid phase, where the antisymmetric sector is also described by the Luttinger liquid, we obtain The position of the peak in the momentum distribution is gauge dependent where α = 1 for the rung gauge and α = 0 for the leg gauge. Since the Luttinger-liquid parameters K ± are positive numbers, one can see by comparing Eqs. (15) and (16) that in vortex-lattice states, the momentum distribution has a larger weight at its peak value than in vortex-liquid states, which is confirmed by the numerical data shown in Fig. 3.
In the rung gauge (the case relevant for experiments), the kinematic and canonical momenta of particles coincide and hence the position of the peak in the momentum distribution function is related to the chiral current (all bosons on leg 1(2) can be thought of as being quasicondensed at the momentum k p (−k p )). In the weakcoupling limit we can hence use the approximation where the sin functions appear due to the presence of the lattice and the fact that the currents are related to particle velocities rather than to their quasi-momenta. From Eq. (18), one can also see that the chiral current reverses its circulation direction in certain VL q states and, in particular, in VL 1/2 states. The chiral current can change its sign in other vortex lattices as well. For instance, if a VL 1/3 state is realized for flux values including the point φ = 2π/3 then, from Eq. (18), we infer that a sign change of j c occurs at φ = 2π/3, where j c vanishes. A similar conclusion has been reached previously by Orignac and Giamarchi [45]. The physical mechanism underlying the chiral-current reversal is a spontaneous increase of the effective flux piercing the unit cell of vortex-lattice states, as discussed in our recent work [87].
C. Vortex lattice at ρv = 1/2 and density ρ = 1/2 In the following we study the vortex-lattice phase for ρ v = 1/2 at filling ρ = 1/2. Figure 4 shows the phase diagram as a function of J ⊥ /J and U/J. While in the limit of hard-core bosons U/J → ∞, there is just a direct transition from a V-MI to a M-MI phase [67], for finite values U/J, an intermediate VL 1/2 -MI phase exists.
First indicators of the vortex-liquid (V-SF) to VL 1/2 -SF transition can be detected in the behavior of j c and avg|j R |. Namely, there is a marked increase of avg|j R | and a kink in the chiral current j c (see Fig. 5  the jump in avg|j R |, when transitioning from the vortexliquid to the vortex-lattice state, becomes more and more pronounced. A more precise estimate of the phase transition point is possible via the calculation of the vortex density ρ v . The results are illustrated in Fig. 5(b). For small values of U/J such as U/J = 4 (see Fig. 5(a)), the vortex-lattice phases VL 1/2 -SF and VL 1/2 -MI give rise to an inversion of the sign of the chiral current j c < 0 as discussed in Ref. [87]. The transitions from the SF to MI phases are of the Berezinskii-Kosterlitz-Thouless (BKT) type and their position is the point where the Luttinger-liquid parameter K ρ takes the value K ρ = 1. We extract K ρ from the long-wavelength behavior of the static structure factor, 1 L i,j e i(i−j)k n i n j [99]. For the VL 1/2 -SF to VL 1/2 -MI phases we may verify this estimate by an analysis of single-particle correlation functions, which are predicted to decay as a † l,r a l,r+x ∼ x − 1 4 at the transition. We study the finite-size scaling behavior of peaks in the quasi-momentum distribution function [87] for extracting the BKT-transition point. An example is illustrated in Fig. 6.
For large values of J ⊥ /J, we observe a continuous Ising-type phase transition from the VL 1/2 -MI to the M-MI phase. In order to examine this phase transition, we calculate the ground-state fidelity susceptibility from the overlap of the ground-state wave functions F = Ψ 0 (J ⊥ )|Ψ 0 (J ⊥ + δJ ⊥ ) . Figure 7 depicts the behavior of χ F S (J ⊥ ) in the vicinity of the VL 1/2 -MI to M-MI transition for several system sizes L. The quadratic increase of χ F S (J ⊥ ) ∼ L 2 confirms an Ising-type character [100][101][102].
D. Vortex lattice at ρv = 1/3 and ρ = 0.8 Apart from the vortex lattice at ρ v = 1/2 we have also resolved a VL 1/3 -SF state in the regime of weak interparticle interactions but low density in our previous work [87]. At this point we would like to discuss the fate of the VL 1/3 -SF state as the interaction strength is increased, keeping the value of the flux at φ ≃ 3π/4 and filling ρ = 0.8 fixed. Figure 8 shows the vortex density ρ v (obtained from the Fourier transform of the rung currents, Eq. (9)): it is first constant and pinned at ρ v = 1/3, as expected for a VL 1/3 state, then it increases once the system enters into a vortex-liquid state at U/J ∼ 4.5. Thus, The data shown in Fig. 8 suggest the existence of another state with a commensurate structure in the rung currents. An obvious candidate state would be a vortex lattice at ρ v = 2/5. We show the pattern of local currents and the density modulations in Fig. 9. The finitesize ground-state configuration is seemingly periodic and thus resembles the structure of a vortex lattice. However, a more detailed analysis suggests that this is not the case. First of all, based on arguments from bosonization, the stabilization of a VL 2/5 state at large values of U/J where vortex lattices with smaller periodicities have already melted into vortex liquids is extremely unlikely.
To further elucidate this case we carry out a finite-size analysis of the amplitude of the local density and current modulations for the limit of hard-core bosons (see the Appendix) which indicates that the average rung current vanishes as avg|j R | ∼ 1/L α with α ≈ 0.4 in the thermodynamic limit. Hence, similar to the vortex-liquid and Meissner states, this state would not break a translational symmetry in the thermodynamic limit as local oscillations die out with increasing the system size. A possible interpretation of this region is a Meissner phase with a Luttinger-liquid parameter K ρ < 1/2, which leads to strong correlation effects, a blurred momentum distribution and enhanced density modulations.
There is, however, another possibility that would lead  to a commensurate locking of density, namely the onedimensional analogue of the ν = 1/2 Laughlin state (see [68,70] for a discussion). Here, ν = N/N φ is the ratio of particle number over flux quanta. For this ratio, the results of [77] may hint at the presence of another commensurate phase in the hard-core boson limit. In our case, i.e., for the parameters of Fig. 8, the condition of ν = 1/2 is also approximately fulfilled with φ/π = 3/4 and ρ = 0.8 (the condition of ν = 1/2 would be fulfilled with ρ = 0.8 and φ/(2π) = 0.4, since ρ = 0.8 is equivalent to ρ = 0.2 for hard-core bosons due to particle-hole symmetry). While giving a definite answer to this interesting question is beyond the scope of our work, in the Appendix, we describe a similar situation, which we have encountered studying hard-core bosons close to the boundary between Meissner and vortex-liquid phases for low particle densities [67], which is also the regime discussed by Petrescu et al. [77]. Finally, the case studied here leads us to another conceptual issue, namely, the ambiguity in defining and estimating the vortex density. First, we note that the vortex density has no direct microscopic definition, unlike particle densities and currents. As explained in Secs. IV A and IV B, we extract the vortex density from different quantities, but primarily from the Fourier transform of rung currents from finite-size data with open boundary conditions Eq. (9) or from the position of the maximum in the momentum distribution function computed in the leg gauge of Eq. (3). In most of the cases considered here and, in particular, in vortex-lattice and most vortexliquid states, the vortex densities extracted from these different methods agree with each other. However, this is not necessarily the case in certain other parameter regimes, including the state that we described above at large U/J > 7, ρ = 0.8 and φ = 3/4. Another example is the Meissner phase at small U/J = 2 (compare Fig. 2 and its discussion in Sec. IV A). An analogous behavior in the Meissner phase in the hard-core boson limit will be discussed in the Appendix. For small values of U/J, we resolve another vortexlattice state at vortex density ρ v = 1/4. Figure 10 provides our numerical evidence for the existence of such a vortex lattice. The corresponding configuration of local currents is shown in Fig. 1(d).
We can use the example of the VL 1/4 state to discuss the location and extension of vortices. In general, the current configurations in VL 1/q phases are best interpreted [87] as small portions of Meissner regions (extended over q − 1 consecutive plaquettes and with the screening current circulating around the boundary of the region), separated by vortices, which occupy every q−th plaquette. An analytical estimate of the vortex size l 0 can be obtained in the weak-coupling limit and it gives l 0 ∼ J/(2J ⊥ )a [57], which suggests that the vortices are rather tightly localized objects. Note, however, that in the VL 1/4 -SF state, the rung currents do not fully vanish inside the Meissner portion due to the nonzero screening length.
This observation leads us to the interesting question of how the vortex size compares to the mean distance between the vortex cores. In our previous analysis of VL 1/2 and VL 1/3 states (and vortex fluids in systems with open boundaries) [67,87], we usually observe one length scale. However, VL 1/2 and VL 1/3 states are the most densely packed vortex lattices and hence are not the optimal cases to resolve the size of individual vortices.
Moreover, due to symmetry reasons, in the VL 1/3 states, in the middle rung of each Meissner portion, the rung current vanishes exactly. Such a behavior is expected for any VL q state with odd q, since the Meissner phase of the ladder with an odd number of rungs has a reflection symmetry with respect to the middle rung of the ladder, accompanied with a reversal of the current circulation direction, which implies an exact vanishing of the rung current in the middle rung.
The local particle density shows a strong modulation in the vortex-lattice phases (except for the vortex lattice at the maximal possible vortex density ρ v = 1/2). This effect has not been captured in previous bosonization studies [45]. Moreover, in the VL 1/4 -SF, we also observe a modulation of the absolute values of the local rung currents, also not captured by previous bosonization analyses.
To summarize, our interpretation of certain plaquettes as the position where the vortex cores are localized in vortex lattices is based on two facts: first, the direction of the local particle currents around theses plaquettes is opposite to the direction of the chiral current circulating around the Meissner region (thus vortices reduce the overall chiral current, consistent with the flux and J ⊥ dependence of j c ). The second reason is that local particle densities are reduced in the plaquettes where vortices are localized. In order to numerically study the size of vortices, we would need access to less closely packed vortex lattices, which is left for future studies.

G. Stability of vortex lattices in systems with a synthetic lattice dimensions
For experimental realizations using a synthetic lattice dimension [43,47,48], long-range interactions in the rung direction have to be taken into account. In our case this amounts to The ratio U/V depends on the properties of the atomic species and might be controlled externally by means of Feshbach resonances [98] or lattice modulation techniques [103]. A reasonable first approximation is to set U = V , while here we will allow V to vary between 0 < V /J < 1.
In Fig. 11, we present the phase diagram for finite positive values of V , J ⊥ = 1.6J and φ = 0.9π. The presence of large rung interactions V /J favors the M-MI phases and suppresses vortex phases, in agreement with the observation made for hard-core bosons on a three-leg ladder [78].
For the parameters of Fig. 11 and for the relevant case of U = V , a transition from the VL 1/2 to the Meissner phase would be expected for U = V ≈ J. However, it is important to note that the regime of stability strongly depends on the particle filling since for a low filling, the interaction V becomes less relevant: Anticipating the results of Sec. VI, at quarter filling ρ = 0.25 a stable VL 1/2 phase (as well as other interesting phases discussed there) can be found up to large values U = V ∼ 30J.

V. CHIRAL-CURRENT REVERSAL: SPONTANEOUS VS EXPLICIT SYMMETRY BREAKING
A main result of our previous work [87] is the observation of a sign change of the chiral current (i.e., a reversal of its circulation direction) in certain vortex lattice phases. We explained this via the mechanism of an increase of the effective flux seen by the particles as a result of the spontaneous breaking of lattice translation symmetry in the vortex lattices, which results in a q-fold enlarged unit cell. If φ is the flux per plaquette, then the effective flux is φ eff = qφ and therefore, the chiral current is j c = j c (qφ). Since the chiral current is 2π- periodic, this can correspond to a negative current, if, for instance, π < qφ < 2π since j c (qφ) = j c (qφ − 2π).
Among the examples shown in Fig. 1, the VL 1/2 states exhibit this behavior, since these vortex lattices are typically stable for φ π. Examples for the sign change of the chiral current are shown in Figs. 2(a) and 5(a).
Here, we will provide an intuitive explanation for the effect of the chiral-current reversal in bosonic ladders with a unit cell that is larger than just one plaquette.
As a starting point, consider the uniform two-leg ladder geometry presented in Fig. 12(a) (i.e., the geometry that corresponds to Eq. (1)) and the simplest case of noninteracting bosons. As a modification, in Fig. 12(b) the hopping along the rungs indicated by dashed lines is fully suppressed such that the shortest closed path to pick up a phase is the boundary of q elementary plaquettes that the system in Fig. 12(a) is built up from. Hence, the relevant flux in this case is qφ. We will argue that the chiral current for the case shown in Fig. 12(b) is expected to be related to the chiral current of a two-leg ladder with q times less plaquetts, but with a flux qφ per plaquette as shown in Fig. 12(c). If the ladder sketched in Fig. 12(c) is in the Meissner phase, which is the case for J ⊥ > 2 tan qφ 2 sin qφ 2 , then its chiral current will be given by where the 1/q factor follows from the fact that there are q times less links along the boundary of the ladder (i.e., links contributing to the chiral current) for the ladder with q-times less plaquettes as shown in Fig. 12(c) as compared to the case of Fig. 12(a). As a generalization of the cases shown in Figs along the dashed links in Figs. 12(b). The case δ = 1 corresponds to the uniform ladder shown in Fig. 12(a) and the case δ = 0 applies to the case shown in Fig. 12(b). We will next present concrete examples for the minimal cases of q = 2 and q = 3.
A. Rung-dimerized ladders: The case of q = 2 First, we consider the minimal case of a rung-dimerized ladder, i.e., q = 2. In Fig. 13, we present the chiral current as a function of flux for δ = 1, 0.5 and 0. For comparison, in the same plot we also depict the chiral current corresponding to a uniform ladder with a flux of 2φ per plaquette, but twice less links along the legs corresponding to the situation shown in Fig. 12(c).
Obviously, Fig. 13 shows that the chiral current for δ = 0 changes its circulation direction for 0.5 < φ/π < 1. The agreement between the chiral currents of the cases of Fig. 12(b) and (c) for q = 2 is excellent for values of the flux corresponding to the Meissner phase of the model in Fig. 12(c). Most importantly, the effective uniform ladder with doubled flux reproduces correctly the sign of the chiral current for the case of Fig. 12(b) with q = 2. Clearly, the two curves are not identical because they correspond to two distinct microscopic models.
B. Rung-trimerized ladders: The case of q = 3 Next, we discuss the rung-trimerized ladder, i.e., q = 3. In Fig. 14, we present the chiral current as a function of flux for δ = 1, 0.5, 0. For comparison, in the same plot, we also depict the chiral current corresponding to the uniform ladder with a flux of 3φ per plaquette, but a factor of 1/3 less links along the legs as shown in Fig. 12(c).
The chiral current of noninteracting bosons for δ = 0 and for the case of q = 3 changes its circulation direction for 1/3 < φ/π < 2/3. Similar to the q = 2 case, for q = 3, the agreement between the chiral currents of the cases   Fig. 12(c). Most importantly, the effective uniform ladder with tripled flux reproduces correctly the sign of the chiral current for the case of Fig. 12(b) with q = 3.

C. Finite temperatures: Weak-coupling approach
In order to study the temperature dependence of the chiral current under experimentally realistic conditions [44], i.e., typically large particle numbers per site, one may use a weak-coupling approximation of model Eq. (3), introduced in the context of Josephson-junction ladders [57] and applied in the present context in Ref. [87]. In the regime of a large filling ρ ≫ 1 and weak but finite interactions, a suppression of density fluctuations occurs and one can thus map the system to a frustrated XY-model of classical spins This model has been studied using either the effective potentials method [59,104] or a transfer-matrix approach [60] at finite temperatures. We use the latter approach based on the transfer-matrix method [60] directly in the thermodynamic limit to evaluate the chiral current for finite temperatures T (k B = 1) through the generalization of the Hellman-Feynman-theorem Eq. (6). We compute the derivative of the free energy with respect to the flux Figure 15 illustrates the temperature dependence of the chiral current for different values of φ/π. In particular, in the proximity of the VL 1/2 -SF phase (see the φ = 0.9π curve of Fig. 15), it is possible to observe a chiral-current reversal up to temperatures of the order of T ≈ J/2 (for more details, see Ref. [87]). Interestingly, j c may exhibit a local maximum at finite temperatures, which is related to the frustration of the model and may also be observed for noninteracting particles.
It is important to note that the transfer-matrix technique does not capture the correct high-temperature behavior of the model (21) due to the assumed mapping to a one-dimensional chain. Nonetheless, the leading temperature dependence for T ≫ J still comes our correctly, with a decay of the current as j c ∼ T −3 . In fact, also for noninteracting particles, one finds a similar decay

D. Phenomenological chiral-current curve in vortex-lattice states of weakly interacting bosons
We now use the method discussed in the previous section to study the high-density regime and very low temperatures. We thus consider bosons on a uniform ladder (i.e., the geometry of Fig. 12(a)) but in the weakly-interacting regime ρ ≫ 1 and U ≪ Jρ. In [87], we have shown that, for example, for J ⊥ = 0.5J, there are pronounced vortex-lattice states (for the smallest temperature that we could access in the transfermatrix approach) for the following vortex densities: ρ v = 0, 1/2, 2/5, 1/3, 1/4, 1/5 (for which the unit cell thus consists of q = 1, 2, 5, 3, 4, 5 plaquettes, correspondingly). Other vortex lattices are washed out already by a small temperature, which we cannot avoid in the transfermatrix approach. The vortex-lattice state with q = 1 is the Meissner state. Now, in the weak-coupling regime, we can reconstruct a chiral-current curve for those flux values, for which vortex-lattice states are realized from just knowing the vortex-density curve as a function of flux by the following method: in the vortex-lattice states, we use the expression for the chiral current of free bosons with a correspondingly enlarged unit cell. Hence, in the vortex-lattice states with q times increased unit cell we will use the expression Eq. (20), provided that J ⊥ > 2J sin qφ 2 tan qφ 2 (which happens to be the case for all vortex-lattice states that we observe).
In Fig. 16, we see that in vortex-lattice states with a qtimes enlarged unit cell the behavior of the chiral current is captured well by Eq. (20). Hence, we conclude that the response of the chiral current of bosons to a spontaneous increase of the unit cell in the weak-coupling regime is not only qualitatively but also quantitatively similar to the response of the chiral current of noninteracting bosons to an explicit enlargement of the unit cell.

VI. CDW PHASE
At quarter filling ρ = 1/4 and for sufficiently large J ⊥ , U > J, a fully gapped CDW 1/4 phase with a spontaneously broken translational symmetry and a two-fold enlarged unit cell can be observed. We first reported evidence and a theoretical explanation for this state in [67]. An example for the typical configuration of currents and density with staggered rung-density oscillations is sketched in Fig. 1(f) (the data are compiled from the central part of a system with L = 80 rungs, φ = 0.98π, J ⊥ = 3J). The ground-state currents look Meissner-like, the rung-currents being suppressed.
As initially described in Ref. [67], the emergence of this CDW 1/4 phase is best understood from the limit of strong interchain tunneling J ⊥ /J → ∞. By introducing a pseudo-spin-1/2 degree of freedom on a rung r via one may write down an effective spin-1 2 model. To first order in 1/|J ⊥ |, the Hamiltonian is H 1 2 = J H 0 While for small fluxes and in this effective model, the term H 0 1 2 dominates and describes a usual (c = 1) Luttinger-liquid phase, corresponding to the M-SF phase, for fluxes φ → π, the correlated hopping and nearestneighbor Ising-type interaction terms become relevant. At φ = π, we may simplify the effective model to Since for a large filling and due to the correlated hopping basically all tunneling processes are strongly suppressed, the Ising term S z r S z r+1 induces a transition to a doubly degenerate Néel state at quarter filling (ρ = 1/4) and in the vicinity of φ = π. For the original bosonic particles, this corresponds to a CDW 1/4 phase.
In Fig. 17, we show the phase diagram of the CDW 1/4 phase as a function of the flux φ/π and J ⊥ /J in the limit of hard-core bosons U/J → ∞. We estimate the position of the SF-to-CDW 1/4 transition by calculating the Luttinger-liquid parameter K ρ , which at the transition should be K ρ = 1/2. The CDW 1/4 phase remains stable for fluxes φ 0.8π and J ⊥ 1.5J.
We next analyze the stability of the CDW 1/4 phase at finite interactions U/J < ∞ (see Fig. 18). For small U/J, we expect a VL 1/2 -SF phase for sufficiently large fluxes φ → π. As shown in Fig. 18, remarkably, we observe a large regime of a fully gapped VL 1/2 -MI phase with finite staggered rungs currents and a flat entanglement entropy profile S vN (r) (see the inset of Fig. 18(b)). In this MI, one particle is delocalized in each plaquette, with suppressed charge fluctuations between plaquettes. This state results from applying the flux and is not linked to a trivial band insulator in the absence of the flux.
For J ⊥ = 3.2J, we estimate the BKT-transition point from the VL 1/2 -SF to the VL 1/2 -MI phase to be at U ≈ 10J, again from determining the point at which the Luttinger-liquid parameter becomes K ρ = 1. For U 30J, we finally observe the CDW 1/4 phase, with vanishing avg|j R | and a finite staggered charge-densitywave order, indicated by the peak-value of the static density structure factor S n (k = π) (see Fig. 18(a)). The possibly Gaussian phase transition can be located precisely from the pronounced peak of the fidelity susceptibility χ F S /L, which we find to diverge as max (χ F S /L) ∼ L 3/2 .
We also study the CDW 1/4 phase for systems with a synthetic lattice dimension including the rung interactions Eq. (19) for the case U = V . It is important to note that we observe broad regimes of both stable CDW 1/4 and VL 1/2 -MI phases for similar parameter ranges as in Fig. 18 for the V = 0 case (data for V > 0 not shown here). Our preliminary results, however, suggest the possibility of an intermediate V-MI phase for U = V ∼ 30J for V > 0. These results will be published elsewhere.
In Fig. 19, we present the equation of state ρ = ρ(µ) for the case of strong nearest-neighbor rung and onsite interactions, V /J → ∞ and U/J → ∞. Hence, in this limit, only a maximal occupation of a single particle per rung is allowed. The DMRG simulation shows that, due to the strong interactions, already for small interchain tunneling J ⊥ ∼ J an extended CDW 1/4 phase is stabilized at quarter filling ρ = 1/4. Very interestingly, close to the limit of φ = π (see Fig. 19), a CDW 1/3 also emerges at filling ρ = 1/3. This state has a three-fold enlarged unit cell and a density oscillation corresponding to a Néel-state of type | · · · ↑↑↓↑↑↓ · · · in the rung-singlet basis. A detailed analysis of these CDW-states in the rung-hard-core limit will be published elsewhere.

VII. BIASED-LADDER PHASE
In this section, we turn to the discussion of the biasedladder phase (BLP), first discussed using mean-field the- ory by Wei and Mueller [75]. In the BLP phase, the Z 2 symmetry associated with interchanging the leg index ℓ = 1 → 2 and ℓ = 2 → 1 (or, in other words, reflection symmetry with respect to reflections about the middle of each rung) is spontaneously broken. This results in a density imbalance ∆n between the two legs, which serves as the order parameter for this phase. We define ∆n as ∆n = r | n 1,r − n 2,r |/N .
The BLP phase was studied previously in several works [75,76,87,91]. In [87], we established the existence of the BLP phase at intermediate values of U/J ∼ 2 from DMRG simulations for ρ 1. Here, we will apply a theoretical framework valid in the regime of a dilute Bose gas to describe this state and to obtain the phase diagram in the dilute-gas limit. That theory relies on a mapping of the system to a two-component Lieb-Liniger gas, whose parameters we relate to microscopic parame- ters by studying the scattering problem. This theory is described in Sec. VII A, while we complement the analytical analysis by DMRG results for finite densities ρ 1 presented in Sec. VII B.
A. Biased-ladder phase in the dilute limit In this section, we address the limit of a dilute gas of bosons, which in one dimension is a strong-coupling regime, invalidating a mean-field type approach. Moreover, it is not possible to develop an effective fieldtheory approach based on bosonization because the velocities obtained from linearizing the single-particle dispersion vanish together with the density. We will follow an approach that we developed for frustrated onedimensional spin systems close to their saturation magnetization [92,93], a method that works qualitatively the best in the dilute limit.
For convenience, we use the leg gauge to render the system explicitly translationally invariant. Thus (quasi-)momentum is a good quantum number. In that gauge, the Hamiltonian is given by H leg from Eq. (3).

Single-particle dispersion
First, we study the single-particle dispersion on a ladder with a nonzero flux. The dispersion consists of two branches (or bands), labelled by ±: For small values of fluxes, J ⊥ > 2 tan φ 2 sin φ 2 , the lower single-particle band has a unique minimum at k = 0, corresponding to the Meissner phase [45]. It develops Two branches of the single-particle dispersion with doubly degenerate minima at ±k0 in the lowest band for parameters J ⊥ = J and φ = 0.8π. k0, k1 andk0,k1 are the incoming and outgoing momenta in the two-body scattering problem described in Sec. VII A 3. a double minimum for J ⊥ < 2 tan φ 2 sin φ 2 at two points ±k 0 as depicted in Fig. 20, In the free case U = 0, the many-body ground state is infinitely degenerate for periodic boundary conditions. Interactions can lift this degeneracy, which leads to the many interesting quantum phases detected in this model. The relevant question is whether the interacting system prefers to form a condensate that is a fifty-fifty mixture of condensates at each of the minima or to exclusively populate one of the two equivalent minima in the dispersion.

Two-component Lieb-Liniger model
We next follow an approach that we have developed for the dilute Bose-gas limit [92], using a mapping of the original bosonic system with a doublydegenerate single-particle dispersion to a two-component Lieb-Liniger model where Ψ 1,2 (x) are boson field operators corresponding to particles from left and right dispersion minima and n 1,2 are the corresponding density operators. The coupling constants of the effective two-component Lieb-Liniger model are given by where a andã are intra-and inter-species scattering lengths, respectively, related to the corresponding scattering phase shifts via In Eq. (32), p is the relative momentum of the low-energy two-boson scattering problem modulo k 0 . The effective mass m is the same for both species and it diverges at the Lifshitz transition when the two minima of the single-particle dispersion at ±k 0 merge into a single one at k = 0. From the effective model Eq. (30) it follows that for g < g, the ground state corresponds to a single-component Luttinger-liquid state (thus an immiscible state where only one species of bosons are present) with a spontaneously broken Z 2 symmetry, whereas for g >g, energetically a two-component Luttinger-liquid state (and thus a miscible state with both species of bosons present in the ground state) is preferred. The phase transition line between these two states is given by g =g and the phase transition is first order.

Scattering problem
We will extract the relevant scattering lengths from solving the low-energy two-boson scattering problem onshell, hence neglecting the upper dispersion branch in Fig. 20. The scattering state of two particles (i = 1, 2) with momenta k 1 and k 2 has the following energy (34) The two particle wave-function is represented as The upper indices on the amplitudes C ℓ,ℓ i,j indicate that both particles belong to the same leg ℓ. If there are no such upper indices, then those amplitudes correspond to the case in which one particle is on the first leg and the other one on the second leg.
Introducing the total momentum Λ = k 1 + k 2 , we separate the center-of-mass motion where we introduced a relative coordinate r = j − i. The Schrödinger equation leads to the following system of equations for the amplitudes C ℓ,ℓ r , C r , and C −r for r > 1, and to the following system of equations for r ≤ 1, From the structure of Eqs. (38) and Eqs. (39), it follows that C ℓ,ℓ r are real for all r ≥ 0, whereas C −r = C * r and, in particular, C 0 is real.
To extract the intra-species scattering length we set k 1 = k 0 + p 1 and k 2 = k 0 + p 2 and take the limits p 1 → 0 and p 2 → 0. Next, we define the relative momentum p = (k 1 − k 2 )/2 and construct scattering states for r ≥ 1 as follows where δ is the scattering phase shift. Furthermore, we introduced Bogoliubov coefficients in analogy to the free boson case The scattering states of Eq. (40) should, for large relative distance r ≫ 1 and for δ = 0, reproduce the scattering states of noninteracting bosons U = 0, hence |z i | < 1 for i = 1, 2, 3. The real numbers v 1 , v 2 and v 3 will be fixed later together with the scattering phase shift δ. First, we insert the ansatz given by Eq. (40) into the Eqs. (38) to determine the complex coefficients a, b, c and z, where, for physically acceptable solutions, we require |z i | < 1. This leads to the following set of equations: Taking only physically meaningful solutions of this system of equations (in general, there are exactly three such solutions) for a i , b i , c i and z i for i = 1, 2, 3, we insert the ansatz Eq. (40) into Eq. (39) to solve for the unknown quantities δ, C 1,1 0 , C 2,2 0 , C 0 , v 1 , v 2 , and v 3 . Along similar lines, we extract the scattering phase shift for the inter-species scatteringδ. In that case, we can fix total momentum to zero Λ = 0 by considering one boson with momentum k 1 = k 0 + p, p → 0, and another one with −k 1 . Hence, one can search for solutions satisfying C 1,1 r = C 2,2 r , reducing the number of unknown constants. However, one needs to take into account the fact that low-energy scattering states for interspecies scattering with a given total momentum and energy are doubly degenerate. There is no such degeneracy for the intra-species scattering problem outlined above: for intra-species scattering the scattering state is uniquely characterized by total momentum and energy. From the form of the single-particle dispersion, the origin of the degeneracy of the inter-species scattering problem is obvious. For example, the scattering state of zero total momentum and asymptotic momenta k 1 and −k 1 is clearly degenerate with the scattering state of two bosons with momentak 1 and −k 1 as indicated in Fig. 20. When constructing the scattering state for inter-species scattering one therefore has to admix different asymptotic momenta: the incoming state with momenta k 1 and −k 1 will produce a state with similar outgoing momenta superposed with a scattering state with outgoing momentã k 1 and −k 1 . Denoting k 2 = −k 1 andk 2 = −k 1 , we write the following ansatz for inter-species scattering for r ≥ 1, C l,l r = −2 cos Θ k1 cos Θ k2 cos(k 1 r +δ) + vz r −2v 1 cos Θk 1 cos Θk 2 cos(k 1 r −δ) C r = cos Θ k1 sin Θ k2 e −i(k1r+δ) + sin Θ k1 cos Θ k2 e i(k1r+δ) where v, v 1 , and z are real numbers with |z| < 1 to ensure physically meaningful solutions and v 2 can be any complex number. We note that since k 1 and k 2 (as well ask 1 andk 2 ) have opposite sign one has to take the proper branches of the arctan in the definition of the coefficients Θ ki . The unknowns v 2 and z are fixed by inserting Eq. (43) into the system of Eqs. (38). In particular, for z, we obtain and for v 2 , The remaining five unknownsδ, C 1,1 0 = C 2,2 0 , C 0 , v, and v 1 are determined by inserting the ansatz Eq. (43) into the system of Eqs. (39).
With the help of relations Eq. (31) and Eq. (32) we determine the intra-and inter-species scattering lengths and the corresponding Lieb-Liniger coupling strengths. For g >g, both minima at ±k 0 are equally populated in the ground state and a two-component Luttinger-liquid phase is realized for a finite but small density. This is the vortex-superfluid phase. There is no density imbalance between the legs of the ladder since C 1,1 r = C 2,2 r for the inter-species scattering problem. However, for g <g, only one of the minima in the single-particle dispersion is populated and a one-component phase is selected, where a density imbalance between the two legs exists, which, for the U → 0 limit is given by, This state is the BLP superfluid. 4. Phase diagram with the BLP state in the dilute Bose-gas regime In Figs. 21(a) and (b), we present the ground-state phase diagram as a function of J ⊥ and U obtained in the dilute-gas limit for the case of V = 0 and for two values of the flux φ = 0.5π and φ = 0.8π, respectively. One can see that the region in parameter space, in which the BLP-SF phase exists, grows significantly in the parameter plane U/J versus J ⊥ /J when increasing the flux from φ = π/2 to φ = 0.8π.
Apart from the phase transition lines we also indicate the line (dashed line) above which, in the V-SF phase, the intra-species interaction enters into the socalled super-Tonks regime [105][106][107] with a > 0, meaning that intra-species repulsion (i.e., repulsion between the particles with almost the same momenta) is effectively stronger than the hard-core contact repulsion. There, instead of relating the scattering length a to the Lieb-Liniger coupling constant g through Eq. (31) (which would wrongly imply attractive g), the intra-species scattering length should be interpreted as an excluded volume. Inter-species interactions, on the contrary, never enter the super-Tonks regime. The interesting property of the gauge field is that the super-Tonks regime (for intra-species interactions) is attained for the case of contact interactions for finite values of the repulsive interaction. The critical value of the repulsive interactions for attaining the super-Tonks regime in intra-species scattering goes to zero when approaching the Lifshitz point, where the effective mass of the Lieb-Liniger model given in Eq. (33) diverges.
Including repulsive interactions along the rung with V > 0 increases the region of stability of the BLP phase. For example, for φ = 0.5π and at J ⊥ = 1.2J, the transition from the BLP-SF to the V-SF for V = 0 is at U ≃ 0.37J, whereas for V = U , that transition shifts to U/J ≃ 0.5. For φ = 0.8π and at J ⊥ = 3.5J, the transition from the BLP-SF to the V-SF is, for V = 0, at U/J ≃ 2.65, whereas for V = U , the transition shifts to U/J ≃ 7.
Attractive interactions along the rung V < 0 with |V | ≪ U , on the contrary, shrink the region of stability of the BLP-SF, consistent with expectations on physical grounds. For stronger attractions, bound states can develop and a pair superfluid can get stabilized with no density imbalance between the legs of the ladder. For example, for φ = 0.8π, J ⊥ = 3.5J and U = 3J there is a resonance in inter-species scatteringg = 0 at V = V c ≃ −1.84J and for V < V c , instead of the two-component V-SF, a single-component pair-superfluid phase (P-SF) is stabilized. Further decreasing V , the system can eventually collapse. We estimate the instability to a collapse to occur when the intra-species interaction constants become attractive g < 0. For φ = 0.8π, J ⊥ = 3.5J and U = 3J this happens at V /J ≃ −2.6.

B. BLP at finite densities
We next study the BLP phase using DMRG calculations at finite densities. The current configuration of the BLP phase is very similar to the Meissner phase, as can be seen in Fig. 1(e): the current flows only along the boundary of the ladder while the rung currents are suppressed. The particle density, however, exhibits a marked imbalance between the legs, which we calculate from Eq. (27).
In the thermodynamic limit the ground state is thus two-fold degenerate, spontaneously breaking the Z 2 mirror symmetry between the legs. In order to numerically stabilize the simulation of the BLP phase, we add small potentials at the boundary of the ladder explicitly breaking the symmetry of the system. By comparing to simulations with smaller or larger edge potentials, we verify that their presence does not influence the magnitude of the order parameter.
We identify the phase transition by a sharply increasing particle density imbalance between the legs ∆n, as presented in Fig. 22. The data are consistent with a second-order Ising-type transition between the M-SF and the BLP phase, as one might expect from the spontaneous breaking of a Z 2 symmetry. Still, a weak first-order nature of the transition cannot be excluded.

VIII. SUMMARY
In summary, we presented an extensive study of the ground-state physics of repulsively interacting bosons on a two-leg ladder in the presence of a uniform Abelian gauge field. In particular, we focused on the discussion of quantum phases with (spontaneously) broken discrete symmetries, including various vortex-lattice phases, a charge-density-wave phase at quarter filling and the biased-ladder phase. We analyzed the vortex-lattice phases at vortex densities ρ v = 1/2, 1/3 and 1/4 in detail and studied different properties and observables such as the central charge, the structure of local currents and the momentum distribution function in different gauges. The vortex density ρ v can be extracted numerically in several ways, such as from analyzing the local current structure or the momentum distribution. We furthermore characterized the various phase transitions between these phases and neighboring ones. We investigated the stability of vortex lattices against including nearest-neighbor interactions on the rungs, relevant for synthetic-lattice dimension experiments [47,48]. As we showed in our previous work [87], vortex-lattice phases may feature an exotic chiral-current reversal effect. We here discussed how this phenomenon may be understood intuitively for a simplified model of noninteracting bosons with an explicitly broken translational symmetry. Thus, the effect is clearly related to the effective flux seen by the particles, which either results from spontaneously enlarging the unit cell or by constructing models with intrinsically larger unit cells.
From the limit of strong rung-couplings J ⊥ /J, we may understand the emergence of the CDW phase at quarter filling introduced in Ref. [67]. Here, we presented results for its stability as a function of J ⊥ /J, φ and also U/J. Remarkably, at large values of U/J, we observe a direct transition to a fully gapped VL 1/2 -MI phase. For strong nearest-neighbor interactions on the rungs (as realizable with a synthetic lattice dimension), additional CDW phases are stabilized.
Finally, we discussed the properties of the BLP phase starting from an analytical analysis that is set up for the limit of a dilute Bose gas (i.e., a low filling ρ → 0), which allows for an intuitive understanding of the nature of the BLP phase.
Open questions naturally arising from the discussion of our model (1) concern its connection to the physics studied extensively in two dimensions such as the fractional quantum-Hall effect or more specifically, fractional Chern insulators [108][109][110]. In particular, the proposal that Laughlin-like phases exist in this model as conjectured in Ref. [68][69][70]77] offer exiting further possibilities for studies of this apparently simple, yet ultimately very rich two-leg ladder model.
In this appendix, we address the problem of defining the vortex density in certain parameter regimes. For the sake of simplicity, we focus on the limit of hard-core bosons. We consider the case of a low particle filling ρ = 0.2 and increase the flux starting from the Meissner phase (or alternatively, we increase J ⊥ starting from the vortex-liquid state).
We start by analyzing the momentum distribution function in the symmetric leg gauge of Eq. (3), as a function of either J ⊥ /J (Fig. 23(a)) or φ (Fig. 23(b)). For φic < φ < φc, we observe two broad maxima in the momentum distribution n1(k) shown in Fig. 23(b).
The transition to the vortex-liquid state V-SF is indicated by the dashed-dotted line. From Fig. 23, we see that below a certain value of J ic ⊥ or flux φ ic indicated by the dotted lines, the momentum distribution function becomes blurred in the Meissner phase. Anticipat- ing the results of the following discussion we denote the two regions as commensurate (cM-SF) and incommensurate (icM-SF) regions of the Meissner-superfluid phase. Thus, n 1 (k) ceases to be sharply peaked at zero momentum, which we otherwise would expect for the Meissner phase if we linked the position of the maximum of the momentum distribution to the vortex density. In addition to the shallow maximum at k = 0, another maximum at a k m = 0 appears in the momentum distribution function for J ic ⊥ < J ⊥ < J c ⊥ and φ ic < φ < φ c at a k = k m = 0, respectively, and the weight of the momentum distribution continuously shifts from k = 0 to k = k m when moving towards the vortex-liquid phase. In the same parameter regime, we observe strong modulations of local particle densities and currents, which extend deep into the Meissner phase and diminish smoothly when departing from the phase transition (J c ⊥ , φ c ) (dot-dashed lines in the figure) from the vortex-liquid states into the Meissner regime. These oscillations are the combined effect of open boundaries and finite system sizes and die out with increasing the system size as can be seen in Figs. 24(b), 25 and 26(a). As shown in Fig. 25 the numerical data indicates that the average rung-current avg|j R |/J vanishes as ∼ 1/L α as expected for boundary driven effects. The vortex density extracted from the Fourier transform of the rung currents shows a plateau at k max = 0.4 as depicted in Fig. 26(b). That plateau extends deep into the Meissner phase and its presence can thus not be used as an unambiguous measure of vortex density.
Most importantly, we could not find any trace of an actual phase transition between the icM-SF regime with multiple peaks in the momentum distribution and the conventional cM-SF phase with its single maximum in the momentum distribution at k = 0. In particular, the fidelity susceptibility (see Fig. 26(c)) is featureless and the block entanglement entropy does not indicate the presence of a conventional second-order phase transition (see Fig. 26(d)). The blurring of the momentum distribution may indicate that single-particle excitations become either gapped, which would result in a phase that is thermodynamically distinct from the Meissner phase, or that the single-particle correlations still decay algebraically, but much slower than in the Meissner phase at either small J ⊥ /J or φ (see Fig. 26(f)). Even if we did not observe any drastic change in the ground-state characteristics between the Meissner phase and the regime with pronounced finite-size modulations realized close to the boundary of the vortex-liquid state, we observe a distinct level crossing in excited states with negligible finite-size effects as shown in Fig. 26(e). In the bulk of the Meissner phase, the lowest excitation is a single-particle excitation ∆E c = (E 0 (L, N − 1) − 2E 0 (L, N ) + E 0 (L, N + 1)) /2 with E m (L, N ) being the m-th eigenenergy of a system of L rungs and N particles, while near the bound-ary between the Meissner and vortex-liquid phases, the lowest excitation is a charge-neutral excitation from the subspace that has the same number of particles as the ground state, ∆E n (L, N ) = E 1 (L, N ) − E 0 (L, N ). We mark this level-crossing position, which coincides with K ρ = 0.5 (compare Figs. 26(e) and (f)), with a dotted line in Figs. 23-26. Based on our current data, it remains unclear whether we are dealing with a thermodynamically distinct state from the Meissner state such as, e.g., a ν = 1/2 Laughlin state [70,77], where the vortex density is pinned to two times the particle density.