Incommensurate phases of a bosonic two-leg ladder under a flux

A boson two--leg ladder in the presence of a synthetic magnetic flux is investigated by means of bosonization techniques and Density Matrix Renormalization Group (DMRG). We follow the quantum phase transition from the commensurate Meissner to the incommensurate vortex phase with increasing flux at different fillings. When the applied flux is $\rho \pi$ and close to it, where $\rho$ is the filling per rung, we find a second incommensuration in the vortex state that affects physical observables such as the momentum distribution, the rung-rung correlation function and the spin-spin and charge-charge static structure factors.

A remarkable characteristic of charged systems with broken U (1) global gauge symmetry such as superconductors is the Meissner-Ochsenfeld effect 1 . In the Meissner phase, below the critical field H c1 , a superconductor behaves as a perfect diamagnet, i.e. it develops surface currents that fully screen the external magnetic field. In a type-II superconductor, for fields above H > H c1 , an Abrikosov vortex lattice phase is formed in the system, where the magnetic field penetrates into vortex cores. In quasi one-dimensional systems, analogues of the Meissner and Abrikosov vortex lattice have been predicted for the bosonic two-leg ladder 2-5 , the simplest system where orbital magnetic field effects are allowed. It was shown that in this model, the quantum phase transition between the Meissner and the Vortex phase is a commensurate-incommensurate transition [6][7][8] . For ladder systems at commensurate filling, a chiral Mott insulator phase with currents circulating in loops commensurate with the ladder was obtained [9][10][11][12] Initially, Josephson junction arrays [13][14][15][16] were proposed as experimental realizations of bosonic one dimensional systems. 17,18 However, Josephson junctions are dissipative and open systems 19,19-21 that cannot be described using a Hermitian many-body Hamiltonian in a canonical formalism. Moreover, the quantum effects in the vortex phase of the Josephson ladder are weak 22 . Fortunately, with the recent advent of ultracold atomic gases, another route to realize low dimensional strongly interacting bosonic systems has opened [23][24][25] . Atoms being neutral, it is necessary to find a way to realize an artificial magnetic flux acting on the ladder. Alternatively, one can consider the mapping of the two-leg ladder bosonic model to a two-component spinor boson model in which the bosons in the upper leg become spin-up bosons and the bosons in the lower leg spin-down bosons. Under such mapping, the magnetic flux of the ladder becomes a spin-orbit coupling for the spinor bosons. Theoretical proposals to realize either artificial gauge fields or artificial spin orbit coupling have been put forward 26,27 , and an artificial spin-orbit coupling has been achieved in a cold atoms experiment 28,29 Recently, the Meissner effect and the formation of a vortex state have been observed for noninteracting ultracold bosonic atoms bosons on a two leg ladder in artificial gauge fields induced by laser-assisted tunneling 30 . The behavior of the chiral current as a function of the coupling strength along the rungs of the ladder, indicates a diamagnetic phase when it reaches a saturated maximum and a vortex lattice phase when it starts to decrease. This experimental achievement has revived the theoretical interest for bosonic ladders in the presence of magnetic flux and its spinor-boson equivalent in the presence of interaction, where an even richer phase diagram is expected 5,[31][32][33][34][35][36][37][38] .
In the present manuscript, we study the commensurate-incommensurate transitions of the hard-core boson ladder with equal densities in the two legs, for varying interleg coupling and flux 39 and fixed fillings away from half-filling. We confirm that above a threshold in the interleg coupling, the Meissner phase is stable for all fluxes 40 while below that threshold the commensurate-incommensurate phase transition 3 to the vortex phase takes place at large enough flux. However, within the vortex phase, we find that a second incommensuration 39 appears at a flux commensurate with the filling, which we characterize by different observables. The paper is organized as follows. In Sec. I, we present the model and the Hamiltonian and define the observables. In Sec. II we describe the bosonization treatment, the Meissner state and the commensurate-incommensurate (C-IC) transition. In Sec. III, we discuss the second incommensuration as a function of the filling. Finally, in the conclusion we present the phase-diagram emerging for the half-filled case.

I. MODEL AND HAMILTONIAN
The lattice Hamiltonian of the bosonic ladder in a flux 2,3 reads: where the operator b ( †) j,σ destroys (creates) a hard core boson on site j of the σ chain. We have defined σ = ±1/2 as the chain index 3,39 , λ as the flux in each plaquette (corresponding to a Landau gauge with the vector potential parallel to the legs), Ω as the interchain hopping. The te iλσ is the hopping amplitude on the chain σ. A schematic picture of the model and its relevant parameters is shown in Fig.1. This hard-core boson model can be mapped into a spin-ladder model with Dzialoshinskii-Moriya interactions 41,42 , as detailed in Appendix A. As a result of translational invariance and parity, the spectrum of the Hamiltonian (1) is even and 2π-periodic in λ.
The leg-current operator J (j, λ) is defined as: while the rung current is defined as: The average densities of bosons are ρ σ = Nσ L where N σ is the number of particles in chain σ and L is the length of the chain. In the rest of the manuscript, we will be considering a fixed total density ρ = ρ ↑ + ρ ↓ . In the absence of applied flux λ the ground state of the system is a rung-Mott Insulator for ρ = 1 and a superfluid for ρ < 1 43 . This situation is not changed at finite λ so that for ρ = 1 Mott-Meissner and Mott-Vortex phase 11,12,40 are obtained.
For our analysis, we are interested in the following observables: the rung-current correlator C(k) the leg-symmetric density correlator S c (k) the leg-antisymmetric density correlator S s (k) and the leg-resolved momentum distribution n σ (k) The non leg-resolved momentum distribution is n(k) = n ↑ (k)+n ↓ (k). The latter quantity is accessible in time-of-flight spectroscopy 30 . The presence of an artificial magnetic flux λ per plaquette, induces the hopping terms on the chain to acquire a phase that depends on the spin (chain). No double occupancy is allowed due to hard-core interaction.

II. BOSONIZATION OF THE TWO-LEG BOSON LADDER
We apply Haldane's bosonization of interacting bosons 44 to the Hamiltonian (1) assuming that Ω is a perturbation. In the absence of interchain couplings and spin-orbit coupling, the Hamiltonian of the bosons can be written as: , u σ is the velocity of excitations, K σ is the Tomonaga-Luttinger (TL) exponent.
In the case of hard-core bosons, u σ = 2t sin(πρ 0 σ ) and K σ = 1. Introducing the fields θ α = π x Π α , we can represent 44 the boson annihilation operators as: and the density operators 44 as: Here, we have introduced the lattice spacing a, while A m and B m are non-universal coefficients. In the case of hard core bosons at half filling, these coefficients have been found analytically 45 . From Eq.(9), we deduce the bosonized expression of the interchain hopping as: where we have kept only the most relevant term in the renormalization group sense 3 . For a model with equivalent up and down leg as Eq. (1), and in the absence of the spontaneous density imbalance between the chains found for weak repulsion 46,47 , u ↑ = u ↓ and K ↑ = K ↓ , it is convenient to introduce the leg-symmetric and leg-antisymmetric representation: in order to rewrite (8)-(11) as: The Hamiltonian H c describes the gapless leg-symmetric density modes, while H s , which describes the legantisymmetric modes, has the form of a quantum sine-Gordon model [48][49][50] and is gapful for K s > 1/4. In a model of bosons with spin-orbit coupling, H s would describe the spin modes, and H c the total density (i.e. the "charge" in the bosonization litterature) modes. Till now, we haven't considered the effect of the flux λ. We now show that it can be exactly incorporated in the bosonized Hamiltonian. In the absence of interchain hopping Ω, we can perform independent gauge transformations on the upper and the lower leg of the ladder. In particular, the gauge transformation: entirely removes λ from the Hamiltonian. We can then apply the Haldane bosonization (8)- (9) to theb j,σ operators. Combining the resulting expressions with Eq. (17), we see that b j,σ has now a bosonized expression of the form: The boson operators b j,σ can be written in the form (9) with φ σ =φ σ and θ σ (x) =θ σ (x)−σ λx 2a and the Hamiltonian expressed in terms of the fields Π and φ σ now reads: leading to a modified Hamiltonian for the leg antisymmetric modes, As discussed in Ref. 3, when Ω = 0 the λ term is imposing a gradient of θ ↑ − θ ↓ , while the term (11) is imposing a constant value of θ ↑ − θ ↓ . For sufficiently large values of λ it becomes energetically advantageous to populate the ground state with solitons giving rise to an incommensurate phase. In the ladder language, such incommensurate phase is the vortex lattice 3 .

A. Gapful excitations in the Meissner state
The quantum sine-Gordon model (16) is integrable 51,52 and its spectrum is fully determined. The Hamiltonian (16) for λ = 0 has a gap ∆ s ∼ us a |aΩ/u s | 2Ks 4Ks −1 , where a is the lattice spacing, for K s > 1/4. In its ground state θ s ≡ 0[π √ 2]. For 1/4 < K s < 1/2, the excitations above the ground state are solitons and antisolitons with the relativistic dispersion E s (k) = (u s k) 2 + ∆ 2 s . The soliton and the antisoliton are topological excitations of the field θ s that carry a leg current j z s = ±u s K s . In the case where one is considering the gap between the ground state and an excited state of total spin current zero (i. e. containing at least one soliton and one antisoliton), the measured gap will be 2∆ s . When K s > 1/2, the solitons 49 and the antisolitons attract each other and can form bound states called breathers that do not carry any spin current. The measured gap between the ground state and the lowest zero current state will be the mass of the lightest breather 53 In the case of hard core bosons 54 , which is the one considered in the numerical analysis here, we have K c = K s = 1, so ∆ s ∼ Ω 2/3 . In that limit, the Hamiltonian (16) has been studied in relation with spin-1/2 chain materials with staggered Dzialoshinskii-Moriya in a magnetic field [55][56][57][58][59][60] . With a weak spin-spin repulsion logarithmic corrections 55,56 are actually obtained as a result a marginal flow, and ∆ s ∼ Ω 2/3 | ln Ω| 1/6 . Besides the solitons and antisolitons, there are two breathers 61-63 , a light breather of mass ∆ s and a heavy breather of mass √ 3∆ s . The amplitude A 0 in Eq. (21) can be estimated for hard core bosons in the case of low density, using the continuum limit 64,65 or in the case of half-filling 45 . In the first case, where G is the Barnes G function and n 0 is the number of particles per site, while in the second case, for low density, and: for half-filling.

B. Correlation functions in the Meissner state
As for K s > 1/4 the ground state of H s has θ s long-range ordered and the excitations above the ground state are gapped, the system described by (14) is a Luther-Emery liquid 66 . In such a phase, . This behavior is a remnant of the single condensate obtained in the non-interacting case 5,30 . Since we have J ⊥ (x) = 0 and as |x − x ′ | → ∞. Thus, the average rung-current vanishes and its fluctuations are short ranged and commensurate, so that C(k) takes a Lorentzian shape in the vicinity of k = 0.
In the case of density-density and spin-spin correlation functions, we have: Since the field θ s is long-range ordered, exponentials e iβφs and derivatives ∂ n x φ s of its dual field φ s are short-range ordered. As a result, the density correlations decay as (x − x ′ ) −2 at long distance leading to S c (k) = K c |k|/(2π), while the spin-spin correlations are decaying exponentially giving a Lorentzian shape for S s (k). Finally, if we consider the longitudinal spin current, the obtained bosonized expression is: In the Meissner phase, the linear behavior is obtained, with J (λ) = λu s K s (2πa) −1 .

C. Commensurate-Incommensurate transition
Adding the spin-orbit coupling λ in (16) gives a Hamiltonian for the spin modes: Expanding (πΠ s + λ/ √ 2a) 2 and using πΠ s = ∂ x θ s , up to a constant shift, the spin orbit coupling adds a term: to the Hamiltonian (16). Now, if we call N s is the number of sine-Gordon solitons and Ns the number of antisolitons, we have: and the contribution of the spin-orbit coupling is rewritten showing that λ acts as a chemical potential for solitons or antisolitons. On the other hand, the energy cost of forming N s solitons and Ns antisolitons is ∆ s (N s + Ns). When |λ| > λ c = a∆s usKs , there is an energy gain to create solitons (or antisolitons depending on the sign of λ) in the ground state. Because of the fermionic character of solitons 67 , their density remains finite, and we obtain another Luttinger liquid. This is the commensurate-incommensurate transition [6][7][8]68 . A detailed picture can be obtained for K s = 1/2, where solitons can be treated as non-interacting fermions as discussed in App.B.
In the incommensurate (IC) phase, the Hamiltonian describing the Luttinger liquid of solitons is: with θ s =θ s − sign(λ)q(λ)x/ √ 2. The density of solitons is proportional to q(λ), while u * s (λ) is the renormalized velocity of excitations and K * s (λ) is the renormalized Luttinger exponent. We now address the behavior of the observables in the IC phase. Near the transition 8,69 The expression of the spin current in the IC phase now becomes: namely the existence of a finite soliton density reduces the average spin current. This justifies the interpretation of these solitons as vortices letting the current to flow along the legs. For large λ, we have q(λ) ∼ |λ|/(πa), so that the expectation value of the spin current eventually vanishes for large flux values. Let us turn to the momentum distribution. In the IC phase and for finite size L with periodic boundary conditions one has: As a result, for 1/(4K c ) + 1/(4K * s ) < 1 one has: so that now n σ (k) has a peak for k = σq(λ), whose height scales as L 1−1/(4Kc)−1/(4K * s ) . That peak becomes a powerlaw divergence in the limit of L → ∞. Comparing with the non-interacting case 30 , these power-law divergences are the remnant of the Bose condensate 5 formed at k = 0 in the Meissner phase or at k = ±q(λ)/2 in the vortex phase.
Turning to the spin-current correlation function, in the IC phase we have 3,4 : Since 1/2 ≤ K * s ≤ 1, the correlation function C(k) presents in the IC phase two cusps at k = ±q(λ). Turning now to the density correlation function, we have: Since 1 ≤ K c +K * s ≤ 2, we find,after taking the Fourier transform, that both S s (k) and S c (k) possess cusp singularities in the vortex phase, with evident notation for the subscript c/s. In the hard core boson system, with K c = K * s = 1, the cusp singularities become slope discontinuities.
Moreover, the behaviors S c (k) ∼ 2Kc|k| π and S s (k) ∼ Ks|k| 2π as k → 0 signal that both charge and spin excitations are gapless in the vortex phase.
We performed numerical simulations for the hard-core spinless bosons on a two-leg ladder as a function of flux and interchain hopping and for different fillings by means of DMRG simulations 70,71 with Periodic Boundary Conditions (PBC). Simulations are performed for sizes up to L = 64, keeping up to M = 1256 states during the renormalization procedure. The truncation error, that is the weight of the discarded states, is at most of order 10 −6 , while the error on the ground-state energy is of order 5 × 10 −5 at most.
A summary for the behavior of observables and correlation functions across the commensurate-incommensurate transition at two different fillings is shown in Fig. 2 for ρ = 0.75 and in Fig. 3 for ρ = 0.5. In both cases, no spontaneous density imbalance 46,47 between the chains is present. In each panel a) of the two figures we compare the behavior of the Fourier Transform (FT) of the rung-current correlation function C(k) in the Meissner phase and in the Vortex phase. The numerical data confirm the prediction of a structureless shape in the Meissner phase and the appearance of two cusp-like peaks in the Vortex phase, respectively at k = q(λ) and k = 2π − q(λ). Since we show data in the vortex phase far from the transition, q(λ) = λ, as expected. The spin gap closure in the Vortex phase is visible also in the low-momentum behavior of the spin static structure factor S s (k) displayed in each panel b) of the two figures 2 and 3: in the Vortex phase S s (k) = K s |k|/2π while in the Meissner phase S s (k) = S s (0) + ak 2 with S s (0) > 0. In these cases K s = 1 as expected for a hard-core boson system. At large momenta the Lorentzian profile centered at k = π, characteristic of the Meissner phase, is replaced by two slope discontinuities at k = πρ and k = 2π − πρ as expected in the Vortex phase for K s = K c = 1. The same evolution can be seen in the charge static structure factors shown in the c) panels of Figs.2 and 3. The commensurate-incommensurate transition is clearly visible in the momentum distribution shown in panels d) of Figs.2 and 3: in the Meissner phase it presents only one cusp-like peak at k = 0 as expected in a bosonic Tomonaga-Luttinger liquid, while the Vortex phase it is characterized by two peaks with same shape, centered at k = ±q(λ)/2.
For K s = 1/2, the sine-Gordon Hamiltonian (30) can be rewritten as a free fermion Hamiltonian 48,66 allowing a more detailed treatment of the commensurate-incommensurate transition 3,6 . Such a treatment sheds additional light on the physics of this commensurate-incommensurate transition, providing an overall alternative description considering that no differences are expected at a qualitative level away from the K s = 1/2 case. The details of such derivation are accounted for in Appendix B.

III. THE SECOND INCOMMENSURATION APPEARING AT λ ≃ πρ
As λ gets close to πρ, with ρ = N/L is the density per rung and for N/L not small compared to unity, λu s /a becomes of the order of the energy cutoff u s /a and the form (30) for the Hamiltonian cannot be used. In order to describe the low-energy physics at λ = ρπ, it is necessary to choose a gauge with the vector potential along the rungs of the ladder, so that the interchain hopping reads:  Applying bosonization to (41), we obtain from (9) the following representation for the interchain hopping: The latter can be rewritten in terms of SU(2) 1 Wess-Zumino-Novikov-Witten (WZNW) currents 72 : In the case of N/L = 1, the complex exponential of Eq.(43) is replaced 39 by a cosine cos √ 2φ c . At commensurate fillings, N/L is a rational number p/m with p, m mutually prime and a term cos m √ 2φ c is also present in the Hamiltonian. In the presence of such term, the symmetry of the U(1) charge Hamiltonian is lowered to Z m and a spontaneous symmetry breaking giving rise to a charge gap and a long-range ordered e i √ 2φc becomes possible in the presence of long ranged interactions 73 . In such case, an insulating phase with a second incommensuration is obtained. 39 At generic filling, or when the term cos m √ 2φ c is irrelevant , we have an unbroken U (1) symmetry φ c → φ c + γ and θ s → θ s + γ. In such case, the Mermin-Wagner theorem 74,75 precludes long range ordering for φ c + θ s . However, since the perturbation in (43) is relevant in the renormalization group sense and has non-zero conformal spin, it is still expected to give rise to incommensurate correlations at the strong coupling fixed point. To give a qualitative picture of such incommensuration, we turn to a mean-field treatment. Compared with the half-filled case, the assumption φ c = γ would correspond to a spontaneously broken U (1) symmetry, not permitted by the Mermin-Wagner theorem. The Gaussian fluctuations of the φ c modes around the saddle point would in fact restore the U (1) symmetry that one has to assume broken to use a mean-field theory. To partially take into account the effect of these fluctuations, we will first solve the mean-field theory for an arbitrary value of γ, and we will then average the obtained correlation functions over γ. Such averaging procedure ensures that e i √ 2φc = 0, and more generally that the obtained correlation functions respect the U (1) symmetry of the Hamiltonian. Of course, that procedure is not expected to produce quantitative estimates, since the fluctations of φ c are underestimated. In particular, the amplitude of the incommensuration can be less than the one expected from the mean field theory, and the decay exponents of the correlations can be larger. But the mean field treatment is providing some insight on the correlation functions that can reveal the presence of a second incommensuration at the fixed point. Assuming φ c = γ, after the transformation θ s → θ s +γ and φ c → φ c +γ the Hamiltonian H c + H s + H hop can be treated in mean-field theory 76-79 , After defining : using a π/2 rotation around the x axis, J y ν =J z ν , J z ν = −J y ν , and applying abelian bosonization 80 , we rewrite: which allows us to write: In turn, this allows us to solve (44) with h s ∼ Ω 2 and g c ∼ Ω 3 . We obtain a gap in the total density excitations, ∆ c ∼ Ω 2 , while the antisymmetric modes remain gapless and develop an incommensuration. To characterize the incommensuration, we first make a shift of the fieldφ s →φ s + hsx us √ 2 whileθ s →θ s , thus Since for the rung current in the mean-field approximation we have: after the rotation we find: we will have to take the average with respect to φ s and θ s and with respect to γ. The latter averaging partially takes into account the restablishment of the full U (1) symmetry by fluctuations around the mean-field. Averaging over γ gives expressions that are translationally invariant. We find: We therefore see that an incommensuration develops in the k ≃ 0 and k ∼ 2λ = 2πρ component of the rung-current and density-wave correlations. In the Fourier transform peaks are located at πρ , πρ ± h s /u s while the singularities at 2πρ and 2πρ ± h s /u s are discontinuities of slope. Since h s ∼ Ω 2 , the incommensuration increases with interchain hopping. One can repeat the calculation also for the S z operator and its correlation function S s (j − j ′ ) ∼ cos(λj) cos(λj ′ ) 1 |j−j ′ | , giving rise to a peak at k ∼ ±πρ. We note that since we have made very crude approximations to treat the φ c fluctuations, we cannot make accurate predictions on the correct value of the exponents. Regarding the calculation of the momentum distribution, since the boson annihilation operators do not correspond to primary fields of the SU(2) 1 WZNW model, we cannot derive their expression in terms ofθ s andφ s using the SU(2) rotation. However, since e iθs has conformal dimensions (1/16, 1/16) its expression in terms of the fieldsφ s andθ s can be expressed as a sum of operators of conformal dimensions (1/16, 1/16). A general expression for the case λ = π has been derived previously 81 . The general form on n k consisted of three peaks centered at πσ ± h s /u s and πσ or a single broad peak πσ, depending on the value of Ω. In the present case, a broad peak centered at 2λσ or a narrow peak at 2λσ plus satellites centered at 2λσ ± h s /u s are expected. As we noted above, these results can be derived rigorously provided we are at a commensurate filling and a charge gap is formed. At generic filling, or when interactions in the charge sector are insufficiently repulsive, the U (1) symmetry of the term (43) is reestablished by quantum fluctuations. In such case, the mean field treatment is only a suggestion that incommensurate fluctuations will be present in a fully gapless state.
In Fig. 4 we follow the appearance of the second incommensuration in the momentum distribution for the system at ρ = 0.75, spanning from below the critical λ c = πρ = 0.75π up to λ = π, as from panels a) to f). At λ = 0.75π the appearance of the secondary peaks are clearly detectable. The positions of these peaks move towards zero with increasing the flux, and disappear completely at λ = π, where the system is back in the standard Vortex phase. In the presence of the second incommensuration, the position of the peaks is no longer proportional to the applied flux: this is apparent in Fig. 5, where the position k max of the peaks in the momentum distribution is displayed as a function of λ for the case with ρ = 0.75 and Ω/t = 1.25 (red solid dots).
For lower values of Ω/t, the relation q(λ) = λ is valid on a large range, and the possible deviation at critical λ is not appreciable, as it is seen in the Fig. 5 (open black dots). We can also follow the evolution of the momentum distribution at the critical λ c = πρ while varying the interchain hopping Ω. In Fig. 6 we show n(k) for the system at ρ = 0.5 and fixed applied flux λ = π/2 on increasing the interchain hopping. For the case with Ω/t = 0.0625, represented by the dashed black line, the gap in total density is too small to be detected in the present numerical simulation at the L = 64 system size. Yet, at Ω = 0.5 and 0.75 the second incommensuration becomes clearly visible with the predicted appearance of the secondary peaks. As mentioned above, the second incommensuration also shows up in the correlation function for the rung-current. In Fig. 7 we show the FT of this quantity at λ = 0.75π for different fillings. The left panel of Fig. 7 displays the data at a small value Ω/t = 0.0625: here, the system is in the standard Vortex phase (first incommensuration) characterized by peaks located at q(λ) = 0.75π. For the ρ = 0.75 the small interchain hopping leads to second incommensuration too small to be detected for the system size of the present simulation. In all filling cases the peaks are located at k = q(λ) = 0.75π and k = 2π −q(λ) = 5/4π. The right panel of Fig. 7 displays the data at Ω/t = 0.75, which is instead a sufficiently large value so that the second incommensuration becomes sizable: indeed, C(k) gets the expected second incommensuration at the predicted filling ρ = 0.75, while at smaller values of the filling no qualitative differences are seen with respect to the behavior shown in the left panel. At ρ = 1.0 the second incommensuration appears at λ = π and the peak at k = π, and it is still detectable for this applied flux 39 . In our previous study, we found 39 a large region of stability of the second incommensuration near the critical value of λ. In order to see how this region evolves with the filling, we summarize in Fig. 8 the phase diagram obtained from DMRG simulations in PBC for a system size L = 64 at filling ρ = 0.5, in which the boundary in the transition from Meissner to Vortex phase and the extension of the region near λ = π/2 with the second incommensuration are visible. In Fig. 9 we show the behavior of S s (k) for different fillings in two different situations in which only the first or also the second incommensuration appear. In the left panel, the behavior in the standard vortex phase is displayed, after picking small values of λ and Ω/t. In the right panel, we show the behavior at λ = 0.25π: here, the appearance of the second incommensuration is expected at ρ = 0.25, characterized by two peaks develop at k = 0.25π and k = 2π − 0.25π, and with S s (k → 0) getting a sizable finite value and a linear momentum dependence. At the other fillings, the spin correlation function gets small finite values at k = 0 and very low peaks at k = ρπ and k = 2π − ρπ, apart from the case ρ = 0.125 already in the Meissner phase. We conclude this section summarizing in Fig. 10 the effects that the appearance of the second incommensuration produces in the different observables and quantities analyzed in the text. We see that there is almost no effect on the charge static structure factor S c (k): as shown in Fig. 2 or Fig. 3, in fact no sharp slope discontinuity at k = ρπ and k = 2π − ρπ emerges in the second incommensuration with respect to the first, that is the standard vortex case. We remark the difference between the ρ < 1 cases analyzed in the present paper and the ρ = 1.0 case 39 , which instead corresponds to a Mott-insulator, i.e. quadratic behavior at low momenta. The DMRG data show gapless leg-symmetric modes for ρ < 1 in agreement with the Mermin-Wagner theorem 74,75 that implies no breaking of the U (1) symmetry φ c → φ c + γ, θ s → θ s + γ.

IV. CONCLUSION
In conclusion, we have studied the commensurate-incommensurate transition between the Meissner and the vortex state of a two-leg bosonic ladder in an external flux, and the formation of a second incommensuration in the vortex state when the flux is matching the particle density. The predictions of the bosonization treatment and the results of DMRG simulations on the commensurate-incommensurate transition from a commensurate Meissner to a standard incommensurate Vortex phase, and the second incommensuration, have been discussed. As expected from previous results at half-filling 39 , the occurrence of a second incommensuration has been found by the DMRG simulations whenever the ratio between the flux and the filling is equal to π. The developing of the second incommensuration can be followed in the momentum distribution, that can be readily measured in experiments 30 . A qualitative picture of the second incommensuration, based on a phase averaging of mean-field approximation of the bosonized theory has been presented. Our DMRG results have been summarized in the interchain hopping-flux phase diagram Fig. 8 at quarter-filling, displaying the Meissner phase, as well as Vortex phase and second incommensuration. The signatures of the second incommensuration on observables and correlation functions have been summed up in Fig. 10. Our predictions can be tested in current experiments, where observables that we have analyzed and discussed can be accessed. A few questions remain open for future investigations. For example, one could investigate whether a second incommensuration would also be observed in multi-chain systems, such as a ladder with a few legs or a twodimensional array of bosonic chains. From the point of view of bosonization, a more rigorous derivation of the second incommensuration in the case of incommensurate filling would be valuable.
Appendix A: Mapping to a spin ladder In the hard core boson case, a representation 82 in terms of (pseudo) spins 1/2 can be introduced: With such mapping, we can rewrite the Hamiltonian (1) as a two-leg ladder Hamiltonian: where J = t cos λ, D = t sin λ, J ⊥ = Ω 2 , J z ⊥ = U ↑↓ and h = δ/2. The term D is a uniform Dzyaloshinskii-Moriya (DM) 41,42,83,84 interaction, with the DM vector parallel to z. For δ = 0, the two legs of the ladder are exposed to a different magnetic field.

Appendix B: Fermionization approach
We have: where m = πΩA 2 0 a,h = λus 2a and the fermion annihilation operators ψ † R,L are destroying the solitons. The detailed correspondence between the fermionic and bosonic expression of the lattice operators is derived below.
The fermionized Hamiltonian (B1) is obtained by the following correspondence with the boson operators: Within fermionization approach the single-particle correlation function ψ † R ψ L can be evaluated by the singleparticle Green's function of the operators that diagonalize the Hamiltonian (B33) there is no average current between the legs of the ladder. However, it is possible to find fluctuations of the current as shown in Eq.B38. In the commensurate phase the correlators can be evaluated using (B33) and one obtains so that: where K 0 and K 1 are modified Bessel functions. This leads to the result (B40) and one expects an exponential decay of the rung current correlation with correlation length u/(2m).
In the incommensurate phase, the fermion correlation functions are expressible instead in terms of incomplete Bessel functions 85 : Indeed, we have (for T = 0): where (uk F ) 2 + m 2 = h and we have noted that ψ † In Eq.(B16), we have two contributions, one coming from the partially filled upper band, and the other from the filled lower band which was the only contribution in the commensurate case. We see that when h → +∞, k F → ∞ and ψ † R (x)ψ L (x ′ ) → 0 uniformly. We have: and thus: For large distances, |x−x ′ | ≫ u/m, we can neglect the contribution from the lower band. The contribution from the upper band can be obtained from the asymptotic expansions given in Ref. 85 on p. 146, while the simpler derivation can be obtained from physical arguments and is presented in the main text. Indeed, in the case uk F ≪ m, we can make the approximations: giving: Second, in the case of uk F ≫ m, we can linearize the dispersion in the upper band around the points ±k F . We can then make the approximations: This time, we find: We see that the correlator ψ † R (x)ψ L (x ′ ) is smaller by a factor m/(uk F ) ∼ m/h ≪ 1 in that limit. If we had instead written a bosonized Hamiltonian, we would have found that ψ † R (x)ψ L (x ′ ) = 0. With Eqs.(B29), we obtain the expression for the rung-current correlator (B43).
In the fermionic representation (B1), the Hamiltonian is readily diagonalized in the form H = k,r=± (r (u s k) 2 + m 2 − h)c † k,r c k,r , by writing: with: e 2iϕ k = uk+im √ (uk) 2 +m 2 . The commensurate phase 6 is obtained for |h| < |m| and the incommensurate phase for |h| > |m|.
We can express the currents as: and the q ∼ 0 component of n j↑ − n j↓ as: Using (B34), one has 3 J = usλ 4πa in the commensurate phase, and in the incommensurate phase. The finite-size scaling of the leg current has been derived in 86 . As ψ † R ψ L is real, the average rung current vanishes.
However, rung-current fluctuations are non-vanishing. Indeed, with the help of Wick's theorem we obtain: In the commensurate phase, the correlators in (B38) can be evaluated using (B33). One obtains: where K 0 and K 1 are the modified Bessel functions. The exponential decay is thus recovered for |x − x ′ | ≫ u/m. Taking the Fourier transform, we find where E and K are complete elliptic integrals 87 . Using the fermion representation, we can also show that: In the incommensurate phase, the fermion correlation functions are expressible instead in terms of incomplete Bessel functions 85 . The detailed expressions are reported in the Appendix B). For large distances, |x − x ′ | ≫ u/m, we can neglect the contribution from the lower band. The contribution from the upper band can be obtained from the asymptotic expansions given in Ref. 85. In the limit uk F ≫ m, simple physical arguments give: so the Fermi wavevector k F = √ h 2 − m 2 /u s = q(λ)/2. Taking the Fourier transform (B43), we deduce that C(k) has slope discontinuities at k = ±2k F . By contrast, in that limit, we find that (n j↑ − n j↓ )(n j ′ ↑ − n j ′ ↓ ) ∼ (j − j ′ ) 2 as expected from the bosonization arguments.