R Citro

S De Palo

M Di Dio

E Orignac

Quantum phase transitions of a two-leg bosonic ladder in an artificial gauge field

We consider a two leg bosonic ladder in a U (1) gauge field with both interleg hopping and interleg repulsion. As a function of the flux, the interleg interaction converts the commensurateincommensurate transition from the Meissner to a Vortex phase, into an Ising-type of transition towards a density wave phase. A disorder point is also found after which the correlation functions develop a damped sinusoid behavior signaling a melting of the vortex phase. We discuss the differences on the phase diagram for attractive and repulsive interleg interaction. In particular, we show how repulsion favors the Meissner phase at low-flux and a phase with a second incommensuration in the correlation functions for intermediate flux, leading to a richer phase diagram than in the case of interleg attraction. The effect of the temperature on the chiral current is also discussed.

I. INTRODUCTION

Trapped ultracold atoms have provided experimentalists with a unique ability to realize highly tunable quantum simulators of many-body model Hamiltonians, [1][2][3] including quasi-one dimensional systems. 4 Moreover, it has recently become possible to simulate the effect of an applied magnetic field using two-photon Raman transitions [5][6][7] , spin-orbit coupling 8 or optical clock transitions 9 . Such situation gives access to a regime where the interplay of low dimensionality, interaction and magnetic field generates exotic phases such as bosonic analogues of the Fractional Quantum Hall Effect. 10 The simplest system to observe nontrivial effects of an artificial gauge field is the bosonic two-leg ladder. 11 Originally, such systems were considered in the context of Josephson junction arrays in magnetic field [12][13][14] and a commensurate-incommensurate (C-IC) phase transition between a Meissner-like phase with currents along the legs and a Vortex-like phase with quasi long range ordered current loops was predicted. However, in Josephson junction systems, ohmic dissipation 15,16 spoiled the quantum coherence required to observe such transition. In cold atom systems, the Meissner and Vortex states have been observed in a non-interacting case. 11 Moreover, recent progress in superconducting qubits 17 engineering offer another promising path 18,[START_REF] Romero | Quantum Simulations with Photons and Polaritons: Merging Quantum Optics with Condensed Matter Physics[END_REF] for realization of low dimensional bosons in artificial flux. Finally, another relevant quantum platform that can implement artificial classical magnetic fields is that of trapped ions 20 .

The availability of experimental systems has thus renewed theoretical interest in the two leg bosonic ladder in a flux . These works have revealed in that deceptively simple model a zoo of ground state phases besides Meissner-like and Vortex-like ones. At commensurate filling, Mott-Meissner and Mott-Vortex phases 30 as well as chiral Mott insulating phases [21][22][23][START_REF] Romen | Chiral Mott insulators in frustrated Bose-Hubbard models on ladders and two-dimensional lattices: a combined perturbative and density matrix renormalization group study[END_REF] have been predicted. Meanwhile, with strong repulsion and a flux Φ = 2πn with n the number of particles per rung, bosonic analogs of the Laughlin states 49 are expected 34,39,41 . Interactions also affect the C-IC transition between the Meissner-like and the Vortex-like phase 50 . In a previous work 51 , we have considered the effect of attractive interchain interactions on the C-IC transition. Using an analogy with statistical mechanics of classical elastic systems on periodic substrates [52][53][54][55][56] , we have shown that interchain attraction split the single commensurate-incommensurate (C-IC) transition point into (a) an Ising transition point between the Meissner-like phase and a density-wave phase, (b) a disorder point 57,58 where incommensuration develops inside the density-wave phase, and (c) a Berezinskii-Kosterlitz-Thouless (BKT) transition 59,60 where the density wave with incommensuration turns into the Vortexlike phase. The density wave phase with incommensuration can be identified as a melted vortex state while the transition (c) can be seen as a melting of the vortex phase. The density wave competing with the Meissner phase at Ising point (a) is induced by interchain interaction [61][62][63] even in the absence of flux. We have verified the existence of those phases in DMRG simulations of hard core bosons. 51 Since the analogy with classical elastic systems holds irrespective of the sign of the interchain interaction, a similar splitting of the C-IC point should also be present in the repulsive case. Differently from the attractive case, we will show that repulsive interactions stabilize the Meissner phase and make the splitting of the C-IC transition occur in a much narrower region. Moreover, they also favor the appearance of an extra periodic oscillation of the correlation function at a wavevector depending on both flux and interchain hopping, even for fluxes non commensurate with density 38 . Such prediction could be easily traced in current experiments by Bragg spectroscopy, that measures the static structure factor, or by time of flight spectroscopy that measures momentum distribution.

The paper is organized as follows: In Section II we introduce the model and its bosonized version, here we also introduce the observables and their correlation functions. In Section III we discuss the Ising transition and the disorder point by using a fermionization approach based on the Majorana fermion representation. Here we also briefly discuss the effect of the temperature on the spin current and momentum distribution. In Section IV we discuss the emergence of the second incommensuration by using a unitary transformation approach and nonabelian bosonization. Section V presents the numerical results for the hard-core limit in the legs. In Section VI we discuss the major results and give some conclusions.

II. MODEL

We consider a model of bosons on a two-leg ladder in the presence of an artificial U(1) gauge field 40,41 :

H = -t j,σ (b † j,σ e iλσ b j+1,σ + b † j+1,σ e -iλσ b j,σ ) + Ω 2 j,α,β b † j,α (σ x ) αβ b j,β + j,α,β U αβ n jα n jβ , (1)
where σ =↑, ↓ represents the leg index or the internal mode of the atom 8,9,64 , b j,σ annihilates a boson on leg σ on the j-th site, n jα = b † jα b jα , t is the hopping amplitude along the chain, Ω is the tunneling between the legs or the laser induced tunneling between internal modes, λ is the Peierls phase of the effective magnetic field associated to the gauge field, U ↑↑ = U ↓↓ is the repulsion between bosons on the same leg, U ↓↑ = U ⊥ the interaction between bosons on opposite legs. In synthetic dimensions, the contact interaction between atoms of different internal quantum numbers becomes the interchain interaction. This model can be mapped to a spin-1/2 bosons with spin-orbit interaction model 38 , where Ω is the transverse magnetic field, λ measures the spin-orbit coupling, U ↑↑ = U ↓↓ is the repulsion between bosons of identical spins, U ↓↑ = U ⊥ the interaction between bosons of opposite spins.

A. Bosonized description

Let us derive the low-energy effective theory for the Hamiltonian (1), treating Ω and U ⊥ as perturbations, and using Haldane's bosonization of interacting bosons. 65 Introducing 65 the fields φ α (x) and Π α (x) satisfying canonical commutation relations [φ α (x), Π β (y)] = iδ(xy) as well as the dual θ α (x) = π

x dyΠ α (y) of φ α (x), we can represent the boson annihilation operators as:

b jσ √ a = ψ σ (x) = e iθσ(x) +∞ m=0 A (σ) m cos(2mφ σ (x)-2mπρ (0) σ x), (2)
and the density operators 65 as:

n jσ a = ρ σ (x) = ρ (0) σ - 1 π ∂ x φ σ + ∞ m=1 B (σ) m cos(2mφ σ (x) -2mπρ (0) σ x). (3)
Here, we have introduced the lattice spacing a, while A m and B m are non-universal coefficients that depend on the microscopic details of the model. For integrable models, these coefficients have been determined from Bethe Ansatz calculations [66][67][68] while for non-integrable models, they can be determined from numerical calculations of correlation functions. [START_REF] Hikihara | Correlation amplitudes for the spin-1/2 xxz chain in a magnetic field[END_REF][START_REF] Bouillot | [END_REF] Introducing the canonically conjugate linear combinations:

φ c = 1 √ 2 (φ ↑ + φ ↓ ) Π c = 1 √ 2 (Π ↑ + Π ↓ ), (4)
φ s = 1 √ 2 (φ ↑ -φ ↓ ) Π s = 1 √ 2 (Π ↑ -Π ↓ ), (5)
the bosonized Hamiltonian can be rewritten as H = H c + H s , where

H c = dx 2π u c K c (πΠ c ) 2 + u c K c (∂ x φ c ) 2 (6)
describes the total density fluctuations for incommensurate filling when umklapp terms are irrelevant, and

H s = dx 2π u s K s πΠ s + λ a √ 2 2 + u s K s (∂ x φ s ) 2 -2ΩA 2 0 dx cos √ 2θ s + U ⊥ aB 2 1 2 dx cos √ 8φ s (7)
describes the antisymmetric density fluctuations. In Eq. ( 7) and (6), u s and u c are respectively the velocity of antisymmetric and total density excitations, A 0 and B 1 are non universal coefficients [START_REF] Giamarchi | Quantum Physics in One Dimension[END_REF] while K s and K c are the corresponding Tomonaga-Luttinger (TL) exponents, which are perturbative with respect to U ⊥ . They can be expressed as a function of the velocity of excitations u, and Tomonaga-Luttinger liquid exponent K of the isolated chain as:

u c = u 1 + U ⊥ Ka πu 1/2 (8)
u s = u 1 - U ⊥ Ka πu 1/2 (9)
K c = K 1 + U ⊥ Ka πu -1/2 (10)
K s = K 1 - U ⊥ Ka πu -1/2 (11)
For an isolated chain of hard core bosons, we have u = 2t sin(πρ 0 σ ) and K = 1. Physical observables can also be represented in bosonization. The rung current, or the flow of bosons from the upper leg to the lower leg, is:

J ⊥ (j) = -iΩ(b † j,↑ b j ↓ -b † j,↓ b j ↑ ). = 2ΩA 2 0 sin √ 2θ s + . . . (12)
The chiral current, i.e. the difference between the currents of upper and lower leg, is defined as

J (j, λ) = -it σ σ(b † j,σ e iλσ b j+1,σ -b † j+1,σ e -iλσ b j,σ ), ( 13
) = u s K s π √ 2 ∂ x θ s + λ a √ 2 . ( 14
)
The density difference between the chains S z j = n j↑ -n j↓ , is written in bosonization as:

S z j = - √ 2 π ∂ x φ s -2B 1 sin( √ 2φ c -πρx) sin √ 2φ s , (15)
while the density of particles per rung is:

n j = - √ 2 π ∂ x φ c -2B 1 cos( √ 2φ c -πρx) cos √ 2φ s . ( 16
)
Let us discuss some simple limits of the Hamiltonian (7). When Ω = 0, U ⊥ = 0, and λ → 0, the antisymmetric modes Hamiltonian Eq. ( 7) reduces to a quantum sine-Gordon Hamiltonian. For K s > 1/4, the spectrum of H s is gapped and the system is in the so-called Meissner state 12,13 characterized by θ s = 0. In such state, the chiral current increases linearly with the applied flux at small λ, while the average rung current J ⊥ = 0 and its correlations J ⊥ (j)J ⊥ (0) decay exponentially with distance. The antisymmetric density correlations also decay exponentially with distance, while the symmetric ones behave as:

n i n j = - 2K c π 2 (i -j) 2 + e -|i-j|/ξ cos πn(i -j) |i -j| Kc , ( 17
)
where ξ is the correlation length resulting from the spectral gap of H s . With Ω = 0, U ⊥ = 0, the antisymmetric density fluctuations Hamiltonian (7) becomes again a quantum sine-Gordon model that can be related to the previous one by the duality transformation θ s → 2φ s , φ s → θ s /2, K s → 1/(4K s ). For K s < 1, the Hamiltonian H s has a gapped spectrum and φ s = π √ 8 for U ⊥ > 0 yielding a zig-zag density wave ground state and φ s = 0 for U ⊥ < 0 yielding a rung density wave ground state. 61,63,[START_REF] Orignac | [END_REF][73][74][75] In both density wave states, the expectation values of the spin and conversion current vanish, and their correlations decay exponentially. The Green's functions of the bosons also decay exponentially, so that the momentum distribution only has a Lorentzian shaped maximum at k = 0. However, in the zig-zag density wave state (U ⊥ > 0), we have:

S z j S z k ∼ C 1 e -|j-k|/ξ + C 2 cos πn(j -k) |j -k| Kc , (18)
n j n k ∼ - 2K c π 2 (j -k) 2 + C 3 cos πn(j -k) |j -k| Kc e -|j-k|/ξ , (19)
while in the rung density wave (U ⊥ < 0),

n j n k ∼ - 2K c π 2 (j -k) 2 + C 3 cos πn(j -k) |j -k| Kc , (20)
S z j S z k ∼ C 1 e -|j-k|/ξ + C 2 cos πn(j -k) |j -k| Kc e -|j-k|/ξ (21)
where K c depends on the interleg interaction, increasing when it is attractive and decreasing when it is repulsive as indicated in Eq. ( 8). The behavior of density correlations in real space is reflected in the corresponding static structure factors:

S c (q) = j e -iqj n j n 0 , (22)
S s (q) = j S z j S z 0 . (23)
In all phases, S c (q → 0) = 2Kc π |q| + o(q), while S s (q) ∼ S s (0) + Aq 2 + o(q 2 ) indicating that symmetric excitations are always gapless while antisymmetric excitations are always gapped. However, in the rung density wave, S c (q → πn) has a power law divergence ∼ |q-πn| Kc-1 (if

K c < 1) or a cusp ∼ C + C |q -πn| Kc-1 (if 1 < K c < 2)
and S s (q → πn) has only a Lorentzian-shaped maximum while in the zig-zag density-wave, S s (q → πn) shows a cusp or singularity while S c (q → πn) has a Lorentzianshaped maximum. The case of U αα = +∞ is peculiar as K s → 1. The Hamiltonian (1) can then be mapped to the Fermi-Hubbard model (see Sec. A 2). Bosonization of the Fermi-Hubbard model [START_REF] Giamarchi | Quantum Physics in One Dimension[END_REF] shows that the operator cos √ 8φ s is marginal in the renormalization group sense. On the attractive side, [START_REF] Giamarchi | Quantum Physics in One Dimension[END_REF] it is marginally relevant, and the density wave exists for all U ⊥ < 0. However, on the repulsive side, cos √ 8φ s is marginally irrelevant and the staggered density wave is absent.

With both Ω and U ⊥ nonzero and λ = 0, the Hamiltonian H s becomes the self-dual sine-Gordon model. 76,77 When both cosines are relevant (i. e. 1/4 < K s < 1) the Meissner phase (stable for |Ω| |U ⊥ |) is competing with the density wave phases (stable in the opposite limit). The competing phases are separated by an Ising critical point. 76,77 In the case of U αα = +∞, since the density wave is absent for U ⊥ > 0, one only has the Meissner state for all U ⊥ > 0. By contrast, for U ⊥ < 0, the charge density wave exists at Ω = 0 and an Ising critical point is present. Thus, phase diagrams for U ⊥ > 0 and U ⊥ < 0 are very different.

In the presence of flux (λ = 0), the density wave phases are stable. However,for U ⊥ = 0, in the Meissner phase, 12,13 when the flux λ exceeds the threshold λ c the commensurate-incommensurate transition takes place: [78][79][80] the ground state of H s then presents a nonzero density of sine-Gordon solitons forming a Tomonaga-Luttinger liquid. 12,13 The low energy properties of the incommensurate phase are described by the effective Hamiltonian: 13 in the ladder language, J (j) decreases and eventually vanishes for large flux values. Meanwhile, the conversion current correlations, density correlations and the Green's functions of the bosons decay with distance as a power law damped sinusoids. The effect of the interaction between identical spins on the commensurate-incommensurate transition has been largely investigated both numerically and theoretically. 23,27,29,30,38 Since U ⊥ can give rise a phase competing with the Meissner state in the absence of flux, its effect on the commensurate incommensurate transition induced by λ needs to be considered. Indeed, near the transition the scaling dimension of the field cos √ 8φ s is 2K * s (λ) 1, thus the cos √ 8φ s term in Eq. ( 7) is relevant and causes a gap opening. 54,56 A fermionization approach 52,53 allows to show that the flux induced transition remains in the Ising universality class. Moreover, this approach also predicts the existence of a disorder point 57,58 where incom-mensuration develops in some correlation functions even though the gap and the density wave phase persist. For instance, the bosonic Green function reads:

H * = dx 2π u * s (λ)K * s (λ)(πΠ * s ) 2 + u * s (λ) K * s (λ) (∂ x φ)
b j,σ b † k,σ = e i θc (ja) √ 2 e -i θc (ka) √ 2 e iσ θ * s (ja) √ 2 e -iσ θ * s (ka) √ 2 ∼ 1 |j -k| 1 4Kc
e iσq(λ)(j-k) e -|j-k|/ξ , (25) where q = π Π s /(a √ 2) and consequently the momentum distribution

n σ (k) = j b † j,σ b 0,σ e -ikj , (26)
instead of showing power-law divergences 37 at momentum ±q as in the vortex state, presents Lorentzianshaped maximas. In the bosonization picture, the disorder point can be understood as the superposition of the incommensuration induced by λΠ s and the gap opened by cos √ 8φ s . As λ further increases, the dimension K * s (λ) recovers the value K s . In the case of K s > 1, there is a second critical point, λ = λ BKT where K * s (λ BKT ) = 1 and the cos √ 8φ s operator becomes marginal. At that point, a Berezinskii-Kosterlitz-Thouless 59,60 takes place, 54 from the density wave phase to the gapless vortex state. 82 This allows to interpret the density wave state with incommensuration as a melted vortex state. By contrast, if K s < 1, the ground state remains in a gapped density wave for all values of λ > λ c .

Thus from simple low-energy description, we can derive the phases listed in the Table I

(k) C(k) Ss(k ∼ 0) Ss(k ∼ πn) Meissner power law |k| 1 4Kc -1 Lorentzian peak at k = 0 Quadratic Lorentzian peak at k = πn Meissner-DW single Lorentzian peak at k = 0 Lorentzian peak at k = 0 Quadratic ∼ |k -πn| Kc -1 Melted vortex
two Lorentzian peaks at ±q(λ) two Lorentzian peaks at ±2q(λ)

Quadratic ∼ |k -πn| Kc -1 Vortex power law peaks ∼ |k -q(λ)| 1/4Kc +1/4K * s -1 power law peaks ∼ |k -2q(λ)| 1/K * s -1 Linear ∼ |k -πn| Kc +K * s -1 TABLE I.
Fourier transform of the correlation functions of the different observables in the phases as predicted in bosonization. When a power law behavior is indicated, if the exponent is negative the singularity is a divergence. When the exponent is positive, the singularity is a cusp and the power law is simply a leading correction to a constant value. Lorentzian at kmax indicates that in the vicinity of its maximum kmax, the correlation function behaves as 1/(a + b(k -kmax) 2 ) with kmax = πn.

In the fourth column, "Quadratic" indicates that

Ss(k) = Ss(0) + Ss"(0)k 2 /2 + o(k 2 ) while "Linear" indicates a behavior Ss(k) = K * s |k|/π + o(k),

III. ISING TRANSITION AND DISORDER POINT

As discussed above in Sec. II A the application of the flux gives rise to an Ising transition point followed by a disorder point both of which can be described using a Majorana fermion representation.

A. Majorana Fermions representation and Quantum Ising transition

Let us now consider a value of the flux close at the commensurate-incommensurate transition, when K s = 1/2, fermionization 52,53 leads to a a detailed picture of the transition between the Meissner state and the density wave states. The fermionized Hamiltonian reads 51 :

H = -i u s 2 dx 2 j=1 (ζ R,j ∂ x ζ R,j -ζ L,j ∂ x ζ L,j ) -i j=1,2 m j dxζ R,j ζ L,j -ih dx(ζ R,1 ζ R,2 + ζ L,1 ζ L,2 ) + dx h 2 2πu s (27)
where

m j = m + (-) j-1 ∆ with, h = - λu s K s a , (28)
m = 2πΩA 2 0 a, (29)
∆ = π 2 U ⊥ (B 1 a) 2 , (30)
and {ζ ν,j (x), ζ ν ,j (x )} = δ nu,ν δ j,j δ(x -x ) are Majorana fermion field operators.

Hamiltonians of the form ( 27) have previously been studied in the context of spin-1 chains in magnetic field [83][START_REF] Wang | Field-induced ising criticality and incommensurability in anisotropic spin-1 chains[END_REF][START_REF] Essler | [END_REF] or spin-1/2 ladders 86,87 with anisotropic interactions. 88 The eigenvalues of ( 27) are:

E ± (k) 2 = (u s k) 2 + m 2 + h 2 + ∆ 2 ±2 h 2 (u s k) 2 + h 2 m 2 + ∆ 2 m 2 . ( 31
) For m = √ h 2 + ∆ 2 , E -(k) = u ∆ m |k| + O(k 2 )
, and a single Majorana fermion mode becomes massless at the transition 52 between the Meissner and the density wave state as expected at an Ising 89 transition. As a consequence, at the transition, the Von Neumann entanglement entropy

S vN = 1 3 (c c + c Ising ) ln L = 1 2 (1 + 1/2) ln L, while away from the transition it is S vN = 1 3 c c ln L = 1
3 ln L since the total density modes φ c are always gapless. A more detailed discussion of finite size scaling of entanglement entropies is found in Ref. 90.

B. Ising order and disorder parameters

At the point K s = 1/2, the bosonization operators cos θ s / √ 2, sin θ s / √ 2, cos √ 2φ s and sin √ 2φ s can be expressed in terms of the Ising order and disorder operators associated with the Majorana fermions operators of Eq. ( 27) as: [91][92][93][START_REF] Nersesyan | New theoretical approaches to strongly correlated systems[END_REF]

cos θ s √ 2 = µ 1 µ 2 sin θ s √ 2 = σ 1 σ 2 , (32) cos
√ 2φ s = σ 1 µ 2 sin √ 2φ s = µ 1 σ 2 . ( 33
)
With our conventions, for m j > 0 we have µ j = 0, σ j = 0 while m j < 0 we have µ j = 0, σ j = 0.

In terms of the Ising order and disorder fields,

b j,σ = e i θc √ 2 (µ 1 µ 2 + isign(σ)σ 1 σ 2 ) (34) S z j = i(ζ R,1 ζ R,2 -ζ L,1 ζ L,2 ) -2B 1 sin( √ 2φ c -πρx)µ 1 σ 2 , (35)
n j = - √ 2 π ∂ x φ c -2B 1 sin( √ 2φ c -πρx)µ 2 σ 1 , (36)
Let's consider first the case of h = 0, Ω > 0. For U ⊥ = 0 the system is in the Meissner phase with µ 1 µ 2 = 0. As U ⊥ > 0 increases, m 2 = m -∆ changes sign, so that σ 2 = 0 while m 1 remains positive and µ 1 = 0. As a result, sin √ 2φ s = 0, and we recover the zig-zag density wave phase. 50 With U ⊥ < 0, m 2 remains positive, while m 1 is changing sign. As a result, for large |U ⊥ |, σ 1 = 0 giving a nonzero cos √ 2φ s and a rung density wave sets in.

Instead as a function of h, we stress that in case of fixed U ⊥ , Ω and variable h, a phase transition is possible only if m 2 -∆ 2 > 0, i. e. only when for h = 0 we have

µ 1 µ 2 = 0. Then, for h > √ m 2 -∆ 2 , we will have µ 1 σ 2 = 0 (for U ⊥ > 0) or µ 2 σ 1 = 0 (for U ⊥ < 0).
Therefore, as in the case of the transition as a function of U ⊥ , one of the pairs of dual Ising variable is becoming critical at the transition while the other remains spectator.

C. Disorder point

The correlators of the Majorana fermion operators ζ ν,j (x)ζ ν ,j (x ) can be obtained from just two integrals 51 :

I 1 (x) = dk 2π e ikx E + (k)E -(k)(E + (k) + E -(k)) , (37)
I 2 (x) = dk 2π e ikx (E + (k) + E -(k)) (38)
by taking the appropriate number of derivatives with respect to x.

To estimate the asymptotic behavior of the Green's functions, one can apply a contour integral method [START_REF] Bender | Advanced mathematical methods for scientists and engineers[END_REF] as detailed in the Appendix B. The long distance behavior is determined by the branch cut singularities of the denominators in the upper half plane. For I 2 , the cut is obtained for uk = ±i m 2 (1 + ∆ 2 /h 2 ) cosh φ, so

I 2 (x) = O(e -|x| √ m 2 (1+∆ 2 /h 2 )/u
). As a result, the long distance behavior is dominated by I 1 (x). For h < m, its branch cut extends along the imaginary axis from

i|∆ - √ m 2 -h 2 |/u < k < i(∆ + √ m 2 -h 2 ), giving I 1 (x) ∼ e -|∆- √ m 2 -h 2 ||x| u
. This recovers the correlation length diverging as ∼ |m -√ h 2 + ∆ 2 | -1 near the Ising transition.

For h > m, the denominator in I 1 has two branch cuts that terminate into two branch points. The long distance behavior of Ī1 is determined by these two branch points as:

Ī1 (x) ∼ e -∆|x| u e i √ h 2 -m 2 |x| u ϕ 1 (x) + e -i √ h 2 -m 2 |x| u ϕ 1 (x) * , (39) with |ϕ 1 (x)| = O(x -1/2
), so that oscillations of wavevector √ h 2 -m 2 /u appear in the real space Majorana fermion correlators for h > m. The point h = m is called a disorder point. 57,58 If we calculate equal time correlation functions of the conversion current using Wick's theorem, the result depends on products of two Green's functions. The conversion current thus shows exponentially damped oscillations with wavevector 2 √ h 2 -m 2 /u s and correlation length u s /(2∆).

Moreover, the correlation functions of the Ising order and disorder fields are expressed in terms of Pfaffians of antisymmetric matrices whose elements are expressed in terms of the Majorana fermion Green's functions. [START_REF] Mccoy | Statistical Mechanics and Field Theory[END_REF] The presence of exponentially damped oscillations in the Majorana fermions Green's function thus also affects correlation functions of Ising order and disorder operators. [START_REF] Wang | Field-induced ising criticality and incommensurability in anisotropic spin-1 chains[END_REF] More precisely, when the large flux ground state is the CDW, for long distances:

σ 1 (x)µ 2 (x)σ 1 (0)µ 2 (0) ∼ e -2∆|x| u r , (40)
µ 1 (x)σ 2 (x)µ 1 (0)σ 2 (0) ∼ ( µ 1 σ 2 ) 2 = 0, (41)
µ 1 (x)µ 2 (x)µ 1 (0)µ 2 (0) ∼ e -∆|x| u √ r × × cos √ h 2 -m 2 x u , (42)
and when the ground state is the zig-zag density wave the long distance correlations of σ 2 µ 1 and σ 1 µ 2 are exchanged.

D. Effect of finite temperature

From the eigenenergies (31), we find the free energy per unit length as:

f = F L = h 2 2πu s -k B T r=± Λ 0 dk π ln 2 cosh E r (k) 2k B T (43)
The spin current is

J s = -usKs a ∂ h f with: ∂f ∂h = h πu s - r=± Λ 0 dk 2π tanh E r (k) 2k B T ∂E r (k) ∂h (44)
The integral (44) is convergent in the limit Λ → +∞. We can split (44) into a ground state contribution and a thermal contribution:

∂f ∂h = ∂e GS ∂h + r=± Λ 0 dk π 2 e Er (k) k B T + 1 ∂E r (k) ∂h , (45)
and we see that away from the critical point, the latter [START_REF] Sondhi | [END_REF][START_REF] Sachdev | Quantum Phase Transitions[END_REF] is represented on Fig. 1. The region where the corrections are linear in temperature is the quantum critical region. At fixed temperature, varying the applied flux, two regimes are possible. For k B T min(∆, m), only a narrow region of flux around the critical flux is inside the quantum critical region, and the current versus flux curve is barely modified. For k B T min(∆, m), the current versus flux curve is showing a broadened maximum that shifts progressively to higher flux. This behavior is shown on Fig. 2.

contribution is O(e -E-(0)/(k B T ) ) when E -(0) k B T . For E -(0) k B T the thermal contribution becomes O(k B T ). A crossover diagram
If we turn to the current susceptibility, which has a logarithmic divergence at the critical flux in the ground state, its positive temperature expression is:

∂ 2 f ∂ 2 h = 1 πu - r=± +∞ 0 dk 2π ∂ 2 E r (k) ∂h 2 tanh E r (k) 2k B T + ∂E r (k) ∂h 2 1 2k B T cosh 2 Er(k) 2k B T   , (46)
Exactly at the critical point h = √ m 2 -∆ 2 , we find that:

∂ 2 E -(k) ∂h 2 = h 2 2m∆u|k| + O(|k|), (47)
so that ∂J ∂h ∼ -2 m T=10 -1 m T=0.5 m FIG. 2. The current versus flux curves for ∆ = 0.2m at varying temperature. For temperatures small (T = 0.01m) compared with ∆, m the curve is indistinguishable from the zero temperature curve. As temperature become comparable with ∆ the maximum of the current becomes broader and shifts to higher flux. For temperature comparable with m, the maximum becomes very broad.

divergence of ∂ h J is controlled by the integral:

Λ 0 dk E -(k) tanh E -(k) 2k B T . (48)
If we take T = 0, the integral will have a logarithmic divergence in the limit of h → √ m 2 -∆ 2 , indicating the Ising transition. However, for any finite T , the hyperbolic tangent will cutoff the divergence for E -(k) k B T , and give instead a maximum scaling as ∼ ln(1/T ). Therefore, one expects that

∂ h J ∼ -ln[(m- √ h 2 + ∆ 2 ) 2 +(k B T ) 2 ]
. Thus, for very low temperature, the slope of the curve J versus h presents a maximum at h = √ m 2 -∆ 2 indicating the presence of an inflection point instead of the vertical tangent obtained at T = 0. If we turn to correlation functions, since our system is one dimensional, at any nonzero temperature its correlation functions always decay exponentially. [START_REF] Landau | Statistical Physics[END_REF] However, in the quantum Ising chain, the correlation length of operators that are long range ordered at zero temperature has been found [START_REF] Sachdev | [END_REF]100 to behave as ∼ u s (T M ) -1/2 e M/T where M is the gap at zero temperature. By contrast, operators with shortrange ordered correlations in the ground state still have a correlation length ∼ u s /M . The difference between the two classes of operators thus remain distinguishable until T ∼ M . Therefore, in the "renormalized classical" region, the zero temperature power law peaks in n σ (k → 0) turns into a narrow Lorentzian maximum, while the Lorentzian maximas in S c/s (k) and C(k) remain broad. The distinction between CDW and Meissner phase is lost only at a temperature k B T ∼ E -(0).

IV. SECOND INCOMMENSURATION WITH REPULSIVE INTERACTION

In previous investigations 37,38 a second incommensuration (2IC) was obtained when the flux λ = πn for a two-leg ladder of hard core bosons. The 2IC can be associated to the interchain hopping and manifests in the periodic oscillations of the correlation functions at wavevectors formed by a linear combinations of λ and Ω. A very simple picture of the second incommensuration can be obtained in the limit U ⊥ t where one can use a Jordan-Wigner representation for the bosons. Using a gauge transformation, the Hamiltonian (1) can be rewritten:

H = -t j,σ (b † j+1,σ b j,σ + b † j,σ b j+1,σ ) + U j n j↑ n j↓ + Ω 2 j (e iλj b † j,↑ b j,↓ + e -iλj b † j,↓ b j,↑ ). ( 49
)
In terms of the Jordan-Wigner fermions given in Appendix (A4) the interchain hopping has, in general, a complicated non-local expressions:

b † j,↑ b j,↓ = c † j,↑ η j↑ η j↓ c † j,↓ e iπ k<j (n k↑ +n k↓ ) , (50)
However, at half-filling, the charge is gapped so that one can approximate,

e iπ k<j (n k↑ +n k↓ ) (-) j , (51)
and the remaining gapless spin mode described by an effective spin chain model:

H = 4t 2 U n S n • S n+1 . (52)
The antihermitian operator η j↑ η j↓ commutes with the Hamiltonian, and can be replaced by one of its eigenvalues ±i. Then, the interchain hopping reduces to:

Ω 2 j (e i(λ-π)j ic † j ↑ c j ↓ + H.c.), (53)
and, having in mind the Jordan-Wigner transformation (A4), it reduces to Ω j S y j when λ = π. Therefore, it acts on the spin chain (52) as a uniform magnetic field, and induces a magnetization along the y axis. Such magnetization also gives rise to incommensuration 71 in the correlation functions of the spin components x and z. This treatment represents the simplest way to understand the origin of a second-incommensuration in the correlation functions. However, in the case away from half-filling, the second incommensuration could not be deduced as straightforwardly 38 and one had to resort to a modified mean-field theory.

Here, we want to present another approach, using a canonical transformation that avoids some of the shortcomings of the mean-field theory. If we bosonize the Jordan-Wigner fermionic version of the Hamiltonian (49) we obtain:

H = ν=c,s dx 2π u ν K ν (πΠ ν ) 2 + u ν K ν (∂ x φ ν ) 2 + Ω 2πa dx e i √ 2φc (e -i √ 2(θs+φs) + e -i √ 2(θs-φs) ) + H.c. - 2g 1⊥ (2πa) 2 dx cos √ 8φ s . (54)
The Hamiltonian (54) contains relevant perturbations of conformal spin ±1 that break Lorentz invariance and can in particular shift the minimum of the dispersion 101 of the low energy excitations away from k = 0. In such case, incommensurate modulation of the correlation functions can be observed. [102][103][104][105] On the basis of a modified mean-field theory, we have argued 38 previously that indeed the rung current correlator would develop a incommensuration in its correlation functions. Here, we will follow a different approach using a canonical transformation. Indeed, in the limit K c → 0, the field φ c becomes classical, and the incommensurate character of the correlations can be readily obtained. If we introduce the rescaling φ c = √ K c φc , θ c θc / √ K c , and use the unitary transformation (55)

U = exp -i dx π K c φc (x)∂ x φ s , (55)
we will have:

U † HU = dx 2π u c (π Πc ) 2 + u s K s (πΠ s ) 2 + (u s K s K c + u c ) (∂ x φc ) 2 + u c K c + u s K s (∂ x φ s ) 2 + K c dx(u s K s Π s ∂ x φc -u c Πc ∂ x φ s ) + Ω 2πa dx cos √ 2θ s cos √ 2φ s - 2g 1⊥ (2πa) 2 dx cos √ 8φ s (56)
The spin-charge interaction in the second line of the Hamiltonian is proportional to √ K c 1, and is an exactly marginal perturbation in the renormalization group sense. In a first approximation, we can neglect it. We then obtain the Hamiltonian of an XXZ chain in a uniform transverse field. 106,107 Using a rotation (see App. C) one can find the ground state of that Hamiltonian 106,107 and obtain its correlation functions. In the gapless phase, one finds:

ρ(j)ρ(j ) ∼ σ z (j)σ z (j ) ∼ (-1) j-j |j -j | + 1 2π 2 (j -j ) 2 cos h s (j -j ) u s ± λ(j -j ) (57) j ⊥ (j)j ⊥ (j ) ∼ (-1) j-j |j -j | + 1 2π 2 (j -j ) 2 cos h s (j -j ) u s ± λ(j -j ) , (58)
with h s = O(Ω). The correlation functions will therefore present periodic oscillations of wavevector formed of linear combinations of λ and h s /u s with integer coefficients, i. e. besides the incommensuration resulting from the flux, a second incommensuration resulting from interchain hopping is obtained. At large Ω a gapped phase can form in which either the spin-spin or the rung current correlation will show a quasi long range order. In such case, the oscillations associated with the second incommensuration become exponentially damped, but give rise to Lorentzian-like peaks in the structure factors. When Ω is low, a charge density wave can be stabilized. Such situation is possible in the case of attractive interaction, and making attraction between opposite spins stronger is detrimental to the observation of the second incommensuration. This explains why, in Ref. 51, we were not observing a competition of Ising and second incommensuration in the attractive case. At odds, in the repulsive case, the second incommensuration is very robust. Let is finally note that if we consider the effect of the marginal operator that we have neglected, its main effect will be to modify scaling dimension of the other operators. It may thus affect the location of the phase boundaries, but not the presence of an incommensuration in the correlation functions.

V. THE HARD-CORE LIMIT

In this section we report numerical results on the effect of the interaction between opposite spins when the repulsion between bosons of the same spin is infinite (hard-core case). Here we focus on the repulsive case, since results obtained in the attractive case have been discussed in Ref. 51, where we found that instead of having a single flux-driven Meissner to Vortex transition, the commensurate Meissner phase and the incommensurate Vortex phase leave space to a Meissner charge-density wave and to a melted vortex phase with short range order. The transition from the Meissner to the charge density wave phase was in the Ising universality class, as predicted by fermionization. With a repulsive interaction we find that the observation of the Ising transition becomes difficult even though signatures of a vortex melting remain visible.

We show results from DMRG simulations for the filling ρ = 0.5 per rung. We fix interchain hopping Ω/t and con-sider different values of the applied flux λ with varying the interaction strength U ⊥ . Simulations are performed in Periodic Boundary Conditions (PBC) for L = 32 and up to L = 64 in some selected cases, keeping up to M = 841 states during the renormalization procedure. The truncation error, that is the weight of the discarded states, is at most of order 10 -5 , while the error on the ground-state energy is of order 10 -4 at most. 108,109 At variance with attractive case, for filling different from unity, in the absence of an applied field we do not expect the transition from the superfluid Meissner phase to the density wave phase 63 since repulsion only gives rise to a marginally irrelevant perturbation. Thus, the phase diagram in the presence of flux is expected to be qualitatively different from the one with attraction.

In Fig. 3 we show the response functions, S c (k) for the attractive case and S s (k) for the repulsive case, as they evolve upon increasing the strength of interchain interaction, when λ = 0. As already discussed in Sec. II A,in the hard-core case, attractive interchain interaction is expected to give rise to charge density wave, while in the repulsive case there is not a spin density wave. In the left panel of Fig. 3, peaks in S c (k) (U ⊥ /t < 0) develop at k = π/2 and k = 3π/2 as attraction increases and the system enters the in-phase density wave phase. Meanwhile in the right panel of Fig. 3, S s (k) (U ⊥ /t > 0) never develops peaks and in fact becomes almost flat as the bosons become more localized, as repulsion is increased. Hence, as expected from marginal irrelevance of interchain repulsion, the spin density wave phase is unfavored.

In order to detect the density wave phases we choose a value of Ω/t sufficiently large and an applied flux close to the value at which the commensurate-incommensurate transition between the Meissner and the Vortex phases occurs in the absence of interchain interaction. Let us note that the Luttinger parameters K c and K s have a different dependence on the interchain interaction. In the attractive case K c is enhanced and K s is reduced, thus the region of stability of the Meissner phase is reduced and the system is more prone to reach the in-phase density wave and vortex regime. On the contrary, in the repulsive case K c is reduced and K s is increased and, as a consequence, the Meissner phase becomes more stable at the expense of the Vortex and density wave ones.

We consider the following case: Ω/t = 0.125 at different applied fluxes. At λ λ c (U ⊥ = 0), .i.e just before the C-IC transition occurs, the system never develops a density wave. In Fig. 4 we show the behavior of the spin and the charge response functions S s (k) and S c (k) respectively, for small and large interaction interaction strength. On increasing the strength the spin static structure factor develops shoulders at k = 2k F = ±π/2 signaling the incipient transition towards a density wave phase, while the static structure factor for low momentum show the expected linear behavior S c (k) 2Kc π |k| for gapless charge excitations. K c smoothly decreases as a function of U ⊥ , going from one, as for a non-interacting hard-core Bose system, towards 1/2 as shown from the slope of the low momentum linear behavior of S c (k) (see Fig. 4). The asymmetry between the attractive and the repulsive case persists in the presence of an applied flux, as shown in Fig. 5, where we follow the response functions change when we increase the interaction strength at fixed λ and Ω/t = 0.5. At small |U ⊥ |, panel (c) and (b) we start from Meissner phase where the momentum distribution has a single peak at k = 0, but for larger interaction strength, while in the attractive case we are in a melted Vortex phase, panel (a), in the repulsive case the system is still in the Meissner phase and S s (k) shows only shoulders at k = 2k F (panel (d)). In the following we investigate the system for fixed Ω/t = 0.125 and at a fixed applied flux for which the system is in the Vortex state in the absence of interaction between the chains (U ⊥ = 0). In the absence of the interaction the spin response function S s (k) displays the expected linear behavior at small momentum and a discontinuity in the derivative at k = 2k F 38 . As we increase the interaction strength (panel A in Fig. 6) the spin structure factor develops peaks at k = π/2 and k = 3π/2 and an almost quadratic behavior at small wavevector. The quadratic behavior indicates that spin excitations remain gapped, while the presence of peaks at k = π/2, 3π/2 is the signature of a zig-zag charge density wave (in the ladder language) or a spin density wave (in the spin-orbit language). The momentum distribution as well the rung-rung response function C(k) develop two separate peaks, that show negligible size effects, indicating the presence of an incommensuration. Thus, we can identify the phase to the so-called melted Vortex phase. 51 For large value of interaction, panel C of Fig. 6, the system is in strongly correlated Meissner phase, indeed momentum distribution shows only one peak at k = 0 whose height scales with size as n(k) = L 1-1/(4Kc) f (kL) where f (kL) can be expressed as a ratio of Gamma functions. At the same time S s (k) develops two peaks, whose shape is well fitted by a Lorentzian form, at k = π/2 and k = 3π/2, signalling the incipient transition towards the CDW-Meissner phase.

In panel B we have an intermediate situation. The sig- nature of the zig-zag DW can be detected from the peaks in the spin response function. The absence of a shortrange incommensuration is visible from the position at k = 0 of the peak in n(k) and from rung-rung correlation function, which show negligible size effects. We conjecture that this corresponds to the so-called charge-density Meissner phase.

In the upper panel of Fig. 6 we show the spin current J s as a function of the strength of interchain interaction when the system goes from the Vortex state to the Meissner one: there is no cusp indicating a square root threshold singularity typical of the C-IC transition, instead the spin current only shows at most a vertical tangent indicating a possible logarithmic divergence of its derivative. To summarize, under application of interleg repulsion, the Vortex phase becomes first a melted vortex phase via a BKT transition, then past the disorder point a DW-Meissner is formed, and finally the Meissner state is stabilized at large repulsion.

As discussed in the previous section, in the presence of the so-called second incommensuration, 37,38 the picture becomes more complicated. Indeed, in such a case, nearby λ πn there is a new incommensurate wavevector which gives, in the various structure factors, extra peaks whose magnitude of which is controlled by Ω. In order to illustrate that situation we have made simulations for a larger interchain hopping, namely Ω/t = 0.5, so that we have the C-CI transition nearby λ = nπ in absence of interchain interaction and therefore near the occurrence of the second incommensuration.

In Fig. 7 situation at fixed U ⊥ /t = 1.5 and Ω/t = 0.5 is shown. In the upper panel we follow the spin current as a function of the applied field. It shows the typical behavior of the Meissner phase when it increases as a function of λ, then it rapidly decreases when entering the Vortex phase which is however short-ranged ordered and finally for λ 0.75π enters the quasi long range ordered Vortex phase, as it can also seen from the typical finite size induced oscillations in this quantity. 110 The Meissner phase is shown in panel A, while the melted-Vortex phase is shown in panel B, where the spin response function has the expected peaks at k = π/2 and 3π/2, yet it has the low momentum behavior observed in the presence of a second incommensuration. 38 In this case, in the momentum distribution is possible to see besides the primary peaks also the secondary peaks related to the second incommensuration. These peaks can be seen also in the rung-rung correlation function C(k). However, both of these functions do not show appreciable size effects attesting the short range of the incommensurate order. In panel C we recover the quasi-long range ordered Vortex phase.

As a last comment we want to stress the fact that in the rung-rung current correlation function in the Meissner phase, see panel C of Fig. 6 and panel A of Fig. 7, shows respectively a Lorentzian-like peak and a cusp centered at k = 4k F = π as the result of higher order term in the Haldane expansion when we derive the rung current. This cusp is present since the exponent K c is decreasing with repulsion, thus enhancing the contribution of the contribution at π compared with the attractive case.

VI. CONCLUSIONS

To conclude, we have analyzed the phase diagram of boson ladder in the presence of an artificial gauge field, when a repulsive interchain interaction is switched on. We have shown, using bosonization, fermionization and DMRG approach, that the the commensurateincommensurate transition between the Meissner phase and the QLRO vortex phase is replaced by an Isinglike transition towards a commensurate zig-zag density wave phase. The fermionization approach has allowed to predict the existence of a disorder point after which the bosonic Green's functions and the rung current correlation function develop exponentially damped oscillations in real space while zig-zag density wave phase persists. This phase is recognized as a melted vortex phase. Differently from the attractive interaction, a secondincommensuration, i.e. an extra periodic oscillation of the correlation functions at wavevectors formed by a linear combinations of the flux and the interchain interaction, dominates even away from half-filling. As numerically shown, the hard core limit in the chains favors the zig-zag density wave phase. Our predictions on the melting of vortices in Bose-Einstein condensates and on the second incommensuration in optical lattices can be traced in current experiments by the measuring the static structure factors and momentum distributions, together with the rung current.

which can be deduced easily from the Holstein-Primakoff representation 113 of spin-1/2 operators. With such mapping, we can rewrite the Hamiltonian (1) as a two-leg ladder Hamiltonian in which the upper and the lower leg have uniform Dzyaloshinskii Moriya interaction. In the two leg ladder representation, Ω and U ↑↓ become the rung exchange interaction, t cos(λ/2) and U ↑↑ , U ↑↑ become the leg exchange interaction, t sin(λ/2) becomes the Dzyaloshinskii-Moriya term.

Mapping to spin-1/2 fermions

Another possible mapping in the case of hard core bosons U ↑↑ , U ↓↓ → ∞ is to the Hubbard model. This mapping is only valid when Ω = 0, but it allows to take advantage of the integrability of the Hubbard model. [114][115][116][117] The mapping, is obtained from the Jordan-Wigner transformation 118 given in AppendixA 2, that maps hard core bosons operators b jσ to fermion operators c jσ :

b jσ = η σ c j,σ e iπ k<j c † k,σ c k,σ , (A4) b † jσ b jσ = c † j,σ c j,σ , (A5)
where {η σ , η σ } + = δ σσ . The Hamiltonian (1) with Ω = 0 is rewritten as:

H = -t j,σ (c † j+1,σ e -iλσ c j,σ + H.c.) + U j n j,↑ n j,↓ (A6)
The gauge transformation 119 c j,σ = e -iλσj a j,σ reduces the Hamiltonian (A6) to the Hubbard form. The Hubbard model presents a spin-charge separation. When interactions are repulsive, and away from half-filling, charge and spin modes are gapless, whereas with attractive interactions charge modes are always gapless but spin modes are gapped. In terms of the original bosons, total density modes are always gapless away from halffilling, but the chain antisymmetric density fluctuations are gapped with attractive interaction giving rise to a symmetric density wave phase, gapless with repulsive interaction.

Appendix B: Asymptotic behavior of the Green's functions

To estimate the asymptotic behavior of the Green's functions, we apply a contour integral method [START_REF] Bender | Advanced mathematical methods for scientists and engineers[END_REF] to the integral

Ī1 (x) = ∞ -∞ dk 2π e ik|x| E -(k) . ( B1
)
The function E -(k) has only branch cut singularities in the upper half plane. The branch cuts arise either from h 2 (uk) 2 + h 2 m 2 + m 2 ∆ 2 < 0 or (uk) 2 + m 2 + ∆ 2 + h 2 -2 h 2 (uk) 2 + h 2 m 2 + m 2 ∆ 2 < 0. The first branch cut, obtained for u 2 k 2 < -m 2 (1+∆ 2 /h 2 ) gives a contribution decaying as e -m √ 1+∆ 2 /h 2 |x|/u , that can be ignored for |x| m/u. The contribution of the cuts of the second type depends whether h < m or h > m. For h < m, there is a single branch cut extending along the imaginary axis from i|∆ -√ m 2 -h 2 |/u < k < i(∆ + √ m 2 -h 2 ). We can rewrite the integral (B1) as: For h > m, there are two branch cuts given by:

(uk) 2 + m 2 + m 2 ∆ 2 h 2 = h ± i∆ 1 - m 2 h 2 cosh α, (B3)
and α real. The integration path in the complex plane is represented on Fig. 8. The branch cuts terminate at the branch points k

(±) d = i ∆ u ± √ h 2 -m 2 u
such that E -(k (B4) so that oscillations of wavevector √ h 2 -m 2 /u appear in the real space correlation functions for h > m. The point h = m is called a disorder point 57,58 . Disorder points are known to occur in frustrated quantum Ising chains in transverse field, 120 bilinear-biquadratic spin-1 chains, 121,122 frustrated spin-1/2 123,124 and spin-1 125,126 chains. They can be classified 58 into disorder points of the first kind (with parameter dependent incommensuration) and disorder point of the second kind (with parameter independent incommensuration). In our model, the disorder point is of the first kind.

Appendix C: Second incommensuration and canonical transformation

In this section, we give some details on the rotation 106,107 used to diagonalize the Hamiltonian obtained after the unitary transformation of Eq. (55). First, we rewrite our Hamiltonian (56) using nonabelian bosonization: 127 with g 1 = g 1⊥ . Using a π 2 rotation around the x axis 37 we can rewrite: and perform here the shift of the field φs → φ s -hsx √ 2us x and using a rotation of the SU (2) 1 primary fields 37 , we reexpress the observables as: In the gapless case, taking the expectation value gives the second incommensuration.

H s = 2πv 3 dx( J R • J R + J L • J L ) + g 1 dxJ z R J z L (C1) +g 1⊥ dx(J x R J x L + J y R J y L ) + Ω dx(J y R + J y L ) (C2)
H s = 2πv 3
U † ρ(x)U = ρ 0 - √ 2 π ∂ x φ c + cos √ 2(φ c -2πρ 0 x) cos √ 2 φs + h s u s x (C7) U † σ z (x)U = - 1 πa r,

FIG. 1 .

 1 FIG.1. Crossover diagram for the current. Below the dashed line, the low temperature region with h < hc is the "renormalized classical" regime, while the low temperature region with h > hc is the disordered regime. In both of these regions, the finite temperature correction to the zero temperature current is exponentially small. Above the dashed line, in the quantum critical region, thermal corrections are O(kBT ).

h 2 4πm∆u

 2 ln(1/T ). In the general case, the

FIG. 3 .

 3 FIG. 3. Left panel: S c (k) for the attractive case. Right panel: S s (k) for the repulsive case. Interaction strength is |U ⊥ |/t = ±1.0, 1.5, 2.0, 3.0 and 6.0 from bottom to top curves. Solid black curves indicate a Meissner state, while red solid curves indicate a CDW where the peaks at k = 2kF develop. Data from L = 64 DMRG simulations in PBC at λ = 0 for ρ = 0.5, at Ω/t = 0.5. Color online.

FIG. 4 .

 4 FIG.4. S c (k)and S s (k) for the repulsive case, respectively solid and dotted lines, for U ⊥ /t = 0.5 and 5.5, red and black curves respectively. Data from L = 64 DMRG simulations in PBC at λ = 0.0625π, for Ω/t = 0.125 and ρ = 0.5. Color online.

8 FIG. 5 .

 85 FIG. 5. S c (k) and S s (k) are respectively shown as dark-green and red solid lines in all panels. Blue dotted lines are for spin resolved momentum distribution nσ(k) whose argument has been shifted shifted of π and black dashed lines are for rung-rung correlation function C(k). Panels a, b, c and d are respectively for U ⊥ = -3, -1, 1 and 3. Data from L = 64 DMRG simulations in PBC at λ = π/8, for Ω/t = 0.5 and ρ = 0.5. Color online.

FIG. 6 .

 6 FIG.6. Upper panel shows the spin current Js as a function of strength of interchain interaction, solid red line is only a guide to eye. Panels below show 8S s (k) (red solid line), 4C(k) (black dashed line), and n(k) (blue dotted line) where the argument of this last quantity has been shifted of π. Left, center and right panel shows these quantities for the cases indicated by the point A, B and C in the upper panel, respectively corresponding to cases where the system is in the melted Vortex phase, in the CDW-Meissner phase and in the Meissner phase. Data shown are from DMRG simulations in PBC for L=32. Color online.

FIG. 7 .

 7 FIG. 7. Upper panel shows the spin current Js as a function of the applied flux, red line is only a guide to eye. Panels below show S s (k) (red solid line), C(k) (black line), nσ(k) (blue line) where the argument of this last quantity has been shifted of π. Left, center and right panels show these quantities for the cases indicated by the point A, B and C in the upper panel, respectively corresponding to cases where the system is in the Meissner phase, in the melted-Vortex phase and in the Vortex phase. Data shown are from DMRG simulations in PBC for L=64. Color online.

Ī1

  that I 1 (x) ∼ e -|∆- √ m 2 -h 2 ||x| u. This gives a correlation length diverging as ∼ |m-√ h 2 + ∆ 2 | -1 near the Ising transition.

d ) 2 √h 2 -m 2 |x| u ϕ 1 h 2 -m 2 |x| u ϕ 1

 22121 = 0. The long distance behavior of Ī1 is determined by these two branch points as:Ī1 (x) ∼ e -∆|x| u e i (x) + e -i √ (x) * ,

FIG. 8 .

 8 FIG.8. The integration path in complex k plane for h > m. The red curves are such that E-(k) 2 < 0. The red dots at extremities of the curve are the points where E-(k) = 0.

  at increasing λ:

				zig-zag CDW	zig-zag CDW
	Phase	momentum distribution	rung current structure factor	structure factor	structure factor
		n

  dx( JR • JR + JL • JL ) + g 1 dx JyAs we can see, either we obtain a fixed point with θs long range ordered or a gapless fixed point. In both cases since we may eliminate the ∂ x φs by a shift of the φs field, one has φs = hs √ 2us x. When θs is gapless, this gives rise to the second incommensuration of Ref. 38. To be more precise, if we consider the bosonized expression for the observables:U † ρ(x)U = ρ 0 -

										U † σ z (x)U = -	√ π 2	√ π 2 ∂ x φ s + cos ∂ x φ c + cos √ 2(φ c -2πρ 0 x) sin √ 2(φ c -2πρ 0 x) cos √ 2φ s √	2φ s (C5) (C6)
										U † j ⊥ (x)U = + cos Ω πa + cos	sin 2(θ s + rφ s ) sin √ 2(θ s + rφ s ) cos r=±1 √ √ 2φ c + sin √ √ 2φ c 2θ s cos( √ 2θ s sin( √ 2φ c + λx)	√	2φ c + λx)
	Finally, returning to abelian bosonization, we
	obtain 106,107
	H s =	dx 2π	uK(π Πs ) 2 +	u K	(∂ x φs ) 2 -	π	Ω √	2	dx∂ x φs
	+	2(g 1⊥ + g 1 ) (2πa) 2	dx cos	√	8 φs
	+	2(g 1⊥ -g 1 ) (2πa) 2	dx cos	√	8 θs		(C4)

R

Jy

L +g 1⊥ dx( Jx R Jx L + Jz R Jz L ) + Ω dx( Jz R + Jz L ) (C3)

  r =±

												e ir	√	2( θs+r φs)+irr	hsx us
					-cos	√	2(φ c -2πρ 0 x) cos	√	2 θs	(C8)
	U † j ⊥ (x)U =	Ω πa	r=±1	sin	√	2 θs + r φs +	h s u s	x cos	√	2φ c
	-	π	1 √	2	∂ x φs sin	√	2φ c + sin	√	2 θs cos(	√	2φ c + λx)
					+ sin	√	2 φs +	h s u s	x sin( √	2φ c + λx) (C9)

VII. ACKNOWLEDGMENTS

We acknowledge A. Celi, M. Calvanese Strinati and E. Tirrito for fruitful discussions. Simulations were performed at Università di Salerno, Università di Trieste and Democritos local computing facilities. M. Di Dio and S. De Palo thank F. Ortolani for the DMRG code. E. Orignac acknowledges hospitality from Università of Salerno.

Appendix A: Hard core boson limit and mappings That limit corresponds to U ↑↑ = U ↓↓ → +∞. In that limit, the bosonic ladder can be mapped to an anisotropic two-leg ladder model with Dzyaloshinskii-Moriya 111,112 interaction, and to the Hubbard model.

Mapping to a spin ladder

If we consider hard core bosons, we can use the mapping of hard core bosons to spins 1/2: