Exploring helical phases of matter in bosonic ladders

Ladder models of ultracold atoms offer a versatile platform for the experimental and theoretical study of different phenomena and phases of matter linked to the interplay between artificial gauge fields and interactions. Strongly correlated helical states are known to appear for specific ratios of the particle and magnetic flux densities and they can often be interpreted as a one-dimensional limit of fractional quantum Hall states, thus being called pretopological. Their signatures, however, are typically hard to observe due to the small gaps characterizing these states. Here we investigate bosonic ladder models at filling factor 1. Based on bosonization, renormalization group and matrix product state simulations we pinpoint two strongly correlated helical phases appearing at this resonance. We show that one of them can be accessed in systems with two-species hardcore bosons and on-site repulsions only, thus amenable for optical lattice experiments. Its signatures are sizable and stable over a broad range of parameters for realistic system sizes.

Introduction.-The experimental investigation of ultracold atomic gases has developed in the last years with remarkable stride. Ultracold atom simulators have been adopted to examine several many-body quantum problems, with the opportunity of tuning the ratio between their interactions and kinetic energies and to trap atomic gases in different geometries. In this context, the study of ladder geometries generated through optical lattices offer the possibility of testing the behavior of quantum systems at the border between one-and two-dimensional systems [1]. Their kinetic energy is firmly one-dimensional, but the presence of plaquettes around which the atoms move allows for the introduction of artificial gauge potentials with non-trivial effects.
In the last years, the effects of these artificial magnetic fluxes have been tested for ladder geometries with a transverse direction defined either in real or "synthetic" dimension [2]. In the latter case, an inner degree of freedom of the atoms is used to represent the transverse coordinate in the ladder. Synthetic dimensions offer, in this way, the possibility of defining both sharp edges and artificial gauge fluxes with suitable Raman couplings [3].
The introduction of artificial gauge potentials in general breaks time-reversal symmetry and allows for the onset of helical many-body states, characterized by the appearance of a net chiral current running in counterpropagating directions at the edges of the ladder. Such phenomena have been investigated for both fermionic [4][5][6][7] and bosonic [8,9] non-interacting gases.
Ultracold gases trapped in ladder systems are currently at the focus of great attention, due to their non-trivial behavior reminiscent of both superconducting [10] and quantum Hall systems [11][12][13][14], and the possibility of engineering interacting topological phases of matter [15][16][17][18][19][20]. In particular, these systems are characterized by several commensuration effects. The first kind is related to the physics of Mott insulators and appears for commensurate ratios between the number of particles and sites in the ladder [21][22][23][24][25][26]. The second kind concerns the ratio between the total artificial magnetic flux enclosed by the ladder and the number of plaquettes in the geometry: several non-trivial phases corresponding to crystals of magnetic vortices alternate when increasing the flux per plaquette at fixed particle density [27][28][29][30]. A third, and more elusive kind of commensuration has been at the center of an intensive study in the last years and it is related to the ratio between the number of particles and magnetic fluxes. For several fractional values of this ratio, suitable repulsive interactions among the atoms cause the onset of partially gapped states, characterized by a helical current. For both fermions and bosons, it has been shown that many of these helical many-body states are indeed the one-dimensional ladder limit of fractional quantum Hall states [31][32][33][34][35][36][37][38][39]. They can be considered as pretopological one-dimensional chiral states.
The physics of these many-body states can be understood in terms of specific resonances between the modulation length of the density of particles and the modulation of the hopping phases determining the external fluxes. In particular, bosonization techniques allow us to understand through semiclassical approximations the fermionic states appearing at filling factor 1/3 and the bosonic states at filling factor 1/2 as Laughlin-like pretopological states.
In this work, we focus on the case of 2-leg ladders of bosonic atoms at filling factor ν = 1. Similarly to their fermionic counterpart at ν = 1/2 [35], these states cannot be discussed in terms of a simple semiclassical approximation due to the competition of several operators becoming resonant. The physics at this resonance is in-deed dictated by competing interactions that concur in forming two different phases. The one dominating for hardcore bosons with rung repulsions is originated by the same effective interaction responsible for the formation of paradigmatic examples of intrinsically gapless symmetryprotected one-dimensional phases of matter [40]. In this respect, resonant ladder models offer a mechanism to form partially gapped phases from interactions that are usually irrelevant, thus providing a platform to study the features of novel strongly correlated helical systems that are intrinsically gapless.
Two-dimensional spinless bosonic systems in Hofstadter-Hubbard models in ladder geometries at ν = 1 have already been numerically studied [41,42], and it was pointed out that for this specific ratio and certain ranges of the hopping parameters, signatures of a second incommensurability appear in the correlation functions [43][44][45]. A systematic analysis of the appearance and characterization of the energy gaps in this system, however, is still missing.
In this work we aim to fill this gap, and we discuss in detail the physics of this system through bosonization and matrix-product-state simulation. We show, in particular, that helical states appear in a ladder of hardcore bosons. With suitable choices of the interleg interactions, the signatures of these states are considerably stronger than their fermionic counterpart at filling ν = 1/2. This opens the path for an experimental study of these strongly correlated and helical phases of matter: our simulations and analysis confirm indeed that one of the two helical phases under investigation can be accessed in bosonic ladders with contact interactions only, which is the most realistic scenario for atoms trapped in optical lattices.
In the following, we will discuss the possible phases of matter characterizing this system in the limit of weak tunneling along the rungs, and we will investigate in detail its main observables and correlation functions. Sec. I defines the lattice model and its interactions, which are described through an effective low-energy field-theory in Sec. II. The renormalization group analysis of the system at the ν = 1 resonance is presented in Sec. III and the main features of the model are analyzed in Sec. IV through matrix product state simulations. Finally, our conclusions are summarized in Sec. V. The appendices present the detail of our bosonization and renormalization group analysis.

I. THE MODEL
Ladder models of ultracold bosons trapped in optical lattices have been experimentally realized with rubidium gases both in spacial ladders [1], by using optical superlattices to isolate single two-leg ladders, and in synthetic dimensions [8,9]. In both cases, suitable pairs of Raman lasers allow us to engineer an effective long-time dynamics of the system described by a low-energy tight-binding Hamiltonian characterized by position-dependent phases in the hopping amplitudes of the atoms (see, for example, Refs. [46][47][48][49]). The resulting artificial gauge flux χ in each plaquette depends on the recoil momentum of the Raman lasers and can be varied by tuning the angle between them (see, for example, the experiments [50][51][52]).
In this work, we focus on the case of ultracold bosons trapped in a two-leg ladder geometry (see Fig. 1a). Their long-time dynamics is described by an interacting Hamiltonian of the form: where H 0 represents the kinetic energy of the atoms, whereas H U and H ⊥ define respectively the intraleg and interleg repulsive interactions. We consider a ladder geometry with L rungs and a gas of N atoms.
A. The single-particle physics We begin our analysis from the single-particle features of the model. The kinetic energy H 0 is given by two hopping terms describing a tight-binding model of spinless bosonic particles hopping in the ladder geometry. Our model is characterized by an intraleg tunneling amplitude t and an interleg hopping Ω (see Fig. 1): Hereafter we set the lattice spacing a = 1 for simplicity; the pseudospin index y = ±1, or, equivalently, y =↑, ↓, labels the transverse direction distinguishing two legs of the ladder, and the bosonic annihilation/creation operators b x,y /b † x,y satisfy the bosonic commutation algebra. We have chosen a gauge which preserves the translational invariance along the x direction, and the hopping phases along this direction define an artificial magnetic flux per plaquette χ = Ad l defined by the counterclockwise Aharonov-Bohm phase acquired by an atom moving along any of the ladder plaquettes.
In a translational invariant case, the single-particle Hamiltonian H 0 can be rewritten in momentum space: The spectrum of the Hamiltonian (4) is depicted in Fig. 1 (b) for a small value of Ω/t; it matches the typical spectrum of one-dimensional systems characterized by a spin-orbit coupling (whose role is played by χ in H t ) and the Zeeman splitting Ω. H t defines indeed two cosine dispersions for pseudospin y = ↑ and y = ↓, displaced in momentum space by the flux χ. In the case of bosons, the particles condense around the two minima of these dispersions and their momentum, which can be measured through time-of-flight experiments. As such, the momentum is locked to their spins and therefore with the leg degree of freedom (see, for example, [8,9]).
The role of H Ω is to open a gap 2Ω at k = 0. This gap can be interpreted as the bulk gap of a quantum Hall system in the one-dimensional ladder limit [11]. When Ω exceeds a threshold Ω c (χ), the minimum of the energy dispersion becomes non-degenerate at k = 0 and the system becomes a superfluid whose atoms are aligned in the x direction of the pseudospin. For non-interacting bosonic systems, this corresponds to the Meissner phase characterized by a chiral current proportional to the flux χ [1].
The non-interacting system is indeed characterized by two main phases: the Meissner phase appearing at small fluxes and sufficiently large interleg tunneling Ω, and a vortex phase appearing at larger fluxes or smaller values of Ω [1,10].
The introduction of repulsive interactions strengthen the Meissner phase and allows for the appearance of several partially gapped and strongly correlated states in the vortex phase.

B. Introducing the interactions
In ultracold atom experiments, repulsive contact interactions among atoms are common. For a realization of the ladder model in synthetic dimension, atoms whose position differs only by the transverse y coordinate are actually located on the same physical position. Therefore contact repulsions are represented by the following density-density interactions with n x,y = b † x,y b x,y . Such interactions have the potential to drive the system from the vortex to different helical phases, if the total density n is commensurate to the flux χ [11,34,36].
For ultracold atoms trapped in optical lattices, the ratio between interaction and kinetic energy can be varied by changing the amplitude of the trapping lasers which define the depth of the trapping lattice. Concerning the ratio between U and V ⊥ , tunable intraspecies and interspecies interactions can be achieved through the introduction of suitable magnetic fields via Fano-Feshbach resonances [53]. While rubidium gases might be difficult to tune in a regime with considerable differences between intraspecies and interspecies interactions, a more promising platform might be offered by potassium [54][55][56][57].
For the remainder of this work, and, especially, in our numerical simulations, we will mostly consider the ladder hard-core bosonic (HCB) scenario, in which U/t → ∞ is the dominating interaction which excludes any double occupation of a single site (x, y), i.e. n x,y ≤ 1. We will consider instead a tunable interspecies interaction V ⊥ .

II. EFFECTIVE LOW-ENERGY DESCRIPTION OF THE MODEL
To study the physics of the interacting gas at zero temperature we set up a low-energy description of the system in the continuum limit based on bosonization [58,59]. In order to examine the chiral states appearing at commensurate values of the filling factor ν = πN/Lχ, we adopt the following bosonization identity for bosons [60], accounting for higher harmonics of the density modulations: here k 0 = πN/(2L) is the wavevector associated to the density modulations and we introduced two pairs of dual fields that satisfy the commutation algebra θ y (x ), ϕ y (x) = iπδ yy Θ(x − x). The fields θ y and ϕ y are associated to the charge and current fluctuations of the particles in the two legs and the density operator can be approximated with n(x, y) ≈ (k 0 − ∂ x θ y (x))/π, whereas the current density j y is proportional to ∂ x ϕ y . Finally, β p is a set of non-universal parameters, and, for our purposes, we will set β 0 = 1, β ±2 = 1/2, whereas all the other β's are set to 0. This choice is suitable to approximate hardcore bosons in the limit Ω, V ⊥ → 0 (see Appendix A for more detail). The general effective interacting Hamiltonian is mapped into a two-species Luttinger liquid, which can be described in terms of two decoupled charge c and spin s sectors. In the following we adopt the standard notation ϕ c/s = (ϕ ↑ ± ϕ ↓ )/ √ 2 and θ c/s = (θ ↑ ± θ ↓ )/ √ 2. We also introduce the parameters u q and K q (with q ∈ {c, s}) which represent the velocities and Luttinger parameters of the bosonic modes. The latter encode the effects of the intraleg interactions in the system. The interleg tunneling Ω and repulsion V ⊥ , instead, determine the presence of additional interactions, which can be cast into a generalized sine-Gordon form. The effective bosonized Hamiltonian is derived in Appendix A and assumes the form: with: See also [43,45] for an equivalent derivation. The operator O g is, in general, fast-oscillating due to the phase (2k 0 − χ)x; therefore, it averages to a non-zero value only if its phase becomes position independent, thus the system approaches the resonance 2k 0 = χ, corresponding to the commensurate filling factor ν = 1. The operator O h matches the interaction adopted in [40] to design intrinsically gapless symmetry-protected topological phases of matter. Analogously to these topological states of matter, h > 0 for the repulsive bosonic model, as well as for its fermionic counterpart studied in [35].
In the Hamiltonian (7) we neglected other fastoscillating terms that are responsible for the onset of different pre-topological states (see, for example, [11,31,32]) but are irrelevant in proximity of ν = 1. We also considered systems with densities far from the Mott commensurate states [21,24], therefore we also neglected further oscillating sine-Gordon interactions stemming from the onsite repulsion U . The main role of the onsite repulsion in our model is to affect the Luttinger parameters K q . Additional terms, less relevant in the renormalization group (RG) sense, have been neglected as well.
Finally, the values of the (bare) coupling constants h and g depend on V ⊥ and Ω respectively and can be approximated as: see Appendix A for more detail.
Differently from the pretopological Laughlin-like states [11,31,32], the physics of bosons at ν = 1 is dictated by the interplay of the non-commuting operators appearing in O g and O h . The g-term, in particular, require special attention: both the operators in Eq. (9) have the same amplitudes and scaling dimension, such that a simple semiclassical approach may not suffice for a clear understanding of the behavior of this system in the thermodynamic limit. To overcome this difficulty, we resort to a second-order RG analysis, in analogy with the fermionic ladder models at filling ν = 1/2 [35], and numerical simulations based on matrix product states (MPS).

III. RENORMALIZATION GROUP ANALYSIS
A first insight of the possible thermodynamical phases and properties of the model is given by a simple scaling analysis of the interactions in the effective Hamiltonian (7). The linear combination of fields appearing in the two terms of the O g operator (9) and in O h (8) do not commute. Therefore, in the ground states of the model, each of these three terms favors the minimization of different combinations of densities and currents which are not compatible with each other. Both the contributions in O g are characterized by a scaling dimension D g = K c + K s + K −1 s /2 and, having the same coupling constant, none of the two can dominate over the other. The scaling dimension of O h is instead D h = 2K s , such that this operator would be always irrelevant, in the RG sense, for repulsive contact interaction for which K s > 1. The form of the Hamiltonian (7) suggests the potential existence of three phases: The first is a Luttinger phase in which both interactions are irrelevant. It corresponds to the vortex phase of the ladder, in which both the charge and spin sectors are gapless and all the two-point correlation functions decay algebraically.
The second phase is the phase dominated by the interaction O h . Hereafter, we will call it h-phase. In this phase the spin sector is gapped and the charge sector is gapless. The h-phase, semiclassically, corresponds to the situation in which the field θ s is pinned to one of the minima of O h . As a consequence, the pseudospin fluctuations and rung current on the ladder are suppressed in the bulk of the system. However, it is important to remark that, for finite systems with open boundary conditions, a rung current appears at the edges of the system (as depicted in Fig. 6 d). We emphasize that this phase is qualitatively different from the Meissner phase appearing at small fluxes since, in the Meissner phase, it is the rung current being ordered rather than the pseudospin density, corresponding to the field ϕ s being pinned to a semiclassical minimum.
The third phase is the one in which the operator O g dominates, and we will refer to it as the g-phase. The study of this phase is less straightforward and its characteristics can be intuitively understood through a mean-field analysis based on the mapping into a Wess-Zumino-Witten (WZW) model (see, for example [59]) proposed in [43,45]. We sketch in the following the key elements to gain insight into this phase. At the meanfield level, we can separate the O g operator into the operators o c = e i √ 2θc , acting on the charge sector, and o s,± = e i √ 2(ϕs±θs) , acting on the spin sector. Assuming that the operators o s,± acquire a non-zero expectation value o s,± , the operator o c opens a gap in the charge sector and suppresses the fluctuations of the charge density. The situation is different in the spin sector: the operators o s,± can be mapped into the chiral current operators J y R and J y L of a WZW model describing the spin sector. With an additional mapping from this WZW model into a rotated Luttinger liquid Hamiltonian [43,45], it is possible to show that these operators do not open a gap in the spin sector for small values of their coupling constants. This can be understood by considering that the current operator J y R + J y L constitutes a perturbation proportional to ∂ xφs in the rotated Luttinger liquid: such perturbation shifts the expectation value ∂ xφs , thus it introduces incommensurability without opening a spin gap [43].
In conclusion, mean-field arguments suggest that the g-phase is characterized by a gapped charge sector and a gapless spin sector. Such mean-field analysis neglects the interactions mixing the spin and charge sectors, and may provide only an approximate description of the gphase. This study was performed in [43,45] to investigate the appearance of a second incommensurability effect in the correlation functions of the system for values of Ω comparable with t, with ladders typically close or within the Meissner phase. We believe, however, that a similar analysis can be extended also in our regime for Ω t and systems with well-separated Meissner and ν = 1 resonances. We mention, that for contact interactions only we do not observe the onset of the g-phase in this regime.
The scaling dimensions D h and D g allow us to obtain a simplified phase diagram as a function of the Luttinger parameters (dashed lines in Fig. 2 (b)): the system is gapless for large values of K c and K s (Luttinger/vortex phase); the h-phase appears for small values of K s , such that O h is more and more relevant; the g-phase occupies instead a region for small values of K c and intermediate values of K s . This suggests the possibility of accessing such a peculiar phase within our setup with suitable interactions.
Given the complexity of the g-interactions, we apply a Wilsonian RG study of the Hamiltonian (7) at second order in the interaction parameters g and h to obtain a more accurate phase diagram. In Appendix A 2 we derive the following RG differential equations: Note that in the above we only consider the flow of couplings and Luttinger parameters, fully neglecting the flow of velocities. However, we checked that the velocities flow very slowly and thus yield only minor corrections to the phase diagram, justifying our assumption of constant u q 's. Before analyzing in more detail the solution of these four equations, let us observe that some limiting case can be easily extracted. Setting g = 0 corresponds with the system away from the ν = 1 resonance; in this case only h and K s flow and they follow the standard (second-order) behavior of the sine-Gordon model [58,59] dictated by the Thouless equations. The operator O h is relevant for K s < 1, such that the system at g = 0 flows in the hdominated phase for any bare K s < 1; for bare K s > 1, the system will flow to the gapped h-phase for sufficiently large bare values of h.
This general behavior describes also what happens to the system for large values of the K c Luttinger parameter. For K c , thus D g sufficiently large, the g interaction is highly irrelevant and we expect to find the Luttinger phase for values of K s 1 (if h/t 1) and the gapped h-phase for K s < 1 (see the right side of Fig. 2 (b)).
Let us address in the following the main features of the phase diagram, that, in its full complexity, depends on the four bare parameters of g, h, K c and K s .
In Fig. 2 (a) we present stacks of a restricted twoparameter phase diagram obtained from the numerical solution of the RG equations, as a function of the bare initial values of K c and K s , and considering fixed values of the bare coupling constants in Eq. (10) obtained by setting the non-universal parameter β 2 = 1/2 (see Appendices A and B). The numerical results are obtained by setting Ω = 0.05t and varying values of V ⊥ /t = {1, 3, 5, 7, 9}. We observe that larger values of Ω would in principle enhance the extension of the g-phase. However, they would also considerably increase the extension of the Meissner phase, such that it is convenient to consider Ω/t 1 in order not to overshadow the properties of the system around ν = 1. This is also consistent with the validity of the perturbation theory approach to the renormalization group analysis with relies on the bare parameters in Eq. (10) being smaller than t. Concerning the non-universal constant β 2 , its value is expected to weakly depend on Ω and V ⊥ as in analogous systems [61,62]. The phase diagram displayed in Fig. 2 is qualitatively stable under small variations of this parameter. In general, a smaller value of β 2 increases the ratio of the bare coupling constants g/h and favors the g-phase over the h-phase.
Our strategy to determine the phase diagram is to set an upper threshold t * ≈ t and a lower thresholdτ ≈ t/1000 to the moduli of g and h. If, during the flow, any of the two coupling constants increases beyond the threshold, we consider that the system reached a strongcoupling (partially gapped) phase at a value of the flow parameter l * ∈ R such that g(l * ) or h(l * ) = t * ∈ R. Depending on which coupling constant is dominating, we characterize the gapped phase as being of the g or h kind. In these cases, the gap of the system can be estimated through ∆ ∼ t * e −l * , as depicted in the inset of Fig. 2 (b) for a selection of points representing the hardcore boson system. If instead, both the coupling constants drop below the lower threshold, we conclude that the system flows into the Luttinger/vortex phase.
This strategy, however, must take into account possible divergences of the Luttinger parameters: for certain ranges of the bare Luttinger parameters, K s diverges for a finite value of the flow parameter l (see Appendix B for further detail). This happens, in particular, in regions where g dominates over h, whereas h tends to flow belowτ . When such a divergence occurs, we label the resulting phase based on the dominant coupling in proximity of the divergence. This corresponds to lowering the upper threshold t * to a value around 0.2t for most of the g-phase and we represented this difficulty in Fig. 2 (b) by shading the regions where either divergences appears or the numerical solution converges slowly. The divergences of the Luttinger parameters suggest that the second order RG equations may not be sufficient to rigorously determine the phase diagram and the third order contributions should be considered as well.
Therefore, the resulting phase diagram in Fig. 2 (b) must be considered as a qualitative phase diagram valid for small values of the bare parameters in Eq. (10): we checked that its features are essentially stable for V ⊥ < 5t (thus for a bare value (10) h < 0.09t). Above this threshold, non-physical features appear for large values of K s indicating that higher orders of the perturbation theory are not negligible in this regime.
The second-order results show that the g-phase is considerably smaller than what expected by the simple first-order scaling analysis, whereas the h and Luttinger phases increase their extension. The behavior of the system close to the border between the h and g phase cannot be precisely determined: we cannot distinguish whether there is a direct phase transition between these two gapped phases or rather an extended Luttinger region that separates them.
The phase diagram in Fig. 2 , the phase diagram is very similar, and we resort to a detailed discussion of the projection V ⊥ = t presented in panel (b). The colored regions correspond to the numerical outcomes of the second order RG flow (see main text): the gray region depicts the Luttinger/vortex phase, the yellow region corresponds to the g-phase and the red region to the h-phase (shaded regions mark stiff or divergent solutions of the RG equations). The dotted lines signal the boundaries of the phases, obtained from a first-order scaling analysis: . From first to second order, the size of the g-phase shrinks significantly. Our numerical simulations of interacting systems are projected onto the the path formed by the circles changing color from green to blue as the amplitude V ⊥ /t is increased. Luttinger liquid parameters are extracted from fits of the density/spin fluctuations presented in Fig. 5. The inset shows an estimate of the gap ∆/t * ∝ e −l * computed from second order RG with initial values taken from the MPS simulations. The estimates of the gap are in agreement to the nontrivial change in the correlation functions and chiral current after a critical value of the interaction V ⊥ /t 6, which signal a transition from Luttinger liquid to h-dominated phase. In our simulations, we find no traces of a g-dominated phase.
represents the ladder of hardcore bosons with V ⊥ = 0, which lies in proximity of the phase transition between the Luttinger and the h-phase. In such a scenario, O h is indeed marginal and the role of the resonant O g interac-tion becomes crucial in determining the thermodynamic behavior of the system. This is due to the second order term, proportional to g 2 which enhances the value of h during the RG flow. Away from the resonance χ = 2k 0 , g averages to 0, and h is suppressed for repulsive interactions, which yield K s > 1 (see Eq. (12) below). At resonance, instead, the coupling g not only can give rise to the g-phase, but it also strengthen the h-phase.
To gain a more realistic insight of the system we restrict our discussion to the case of hard-core bosons on the ladder geometry. Based on the mapping of the single leg subsystems into fermions (see Appendix A) we can estimate the bare values of the Luttinger parameters as: It is important to stress that these relations provide an estimate of the bare parameters in the perturbative limit V ⊥ 2πt and Ω t, only. In this respect, we emphasize that there is an important distinction between K c and K s . The spin sector is considerably influenced by the interaction O h , which tends to reduce the value of K s during the RG flow. This implies that, even in the vortex phase, when O h is irrelevant, the physical value of K s displays major deviations from Eq. (12) which considerably overestimates it. The discrepancy between the bare value of Ks in (12) and the corresponding renormalized parameter is even more severe at the resonances. In the Meissner resonance, for example, the tunneling Ω gaps the spin sector and K s is supposed to diverge during the RG flow. At ν = 1, instead, the coupling constant h can be further increased by the g 2 contribution and the flow of K s is modified by g in a non-trivial way. For these reasons, Eq. (12) fails in predicting the physical value of K s for ν = 1 beyond the perturbative regime V ⊥ 2πt. The charge sector, instead, is not affected by O h and, considering the system outside the resonance or away from the g-phase, K c does not considerably change during the RG flow. Hence, differently from the spin sector, Eq. (12) provides a reasonable estimate for the physical value of K c and, indeed, it qualitatively agrees with the K c parameter extracted from the numerical tensor network simulation (see Fig. 5 and the next section). Nonetheless, we expect major deviations for very strong interactions: in the limit Ω → 0 the model can be mapped into a spinless fermionic Hubbard model [11] and K c → 1/2 for V ⊥ → ∞. The lower limit of K c for large V ⊥ , however, depends on Ω and it is not captured by our simple approximation, as showed by the comparison with the numerical results in Fig. 5 (d).
When we consider a system with fixed Ω, thus a fixed bare value of g, the phase diagram depends on the remaining three initial values of the RG flow. However, based on Eqs. (10) and (12) for the model of hard-core bosons, the bare parameters K c , K s , and h are determined from the value of V ⊥ only. This implies that when we vary V ⊥ , the system defines a trajectory in this threeparameter phase diagram (see Fig. 2 (a)).
With respect to Fig. 2 (b), we observe that for larger values of V ⊥ the g-phase is shifted towards larger values of K s . The colored dots in Fig. 2 (b) depict the projection of the spiral trajectory in Fig. 2 (a), determined by the MPS simulations in the next sections, over the K c , K s plane at V ⊥ = t.
Based on this trajectory, we predict that the system at V ⊥ = 0 lies in proximity of the edges between the Luttinger and the two gapped phases, it evolves close to the boundary between the h and g phases for intermediate values of V ⊥ and, finally, it develops a larger gap entering deep in the h-phase.
We emphasize, however, that this second-order RG analysis based on Eqs. (10,11) and (12) provides accurate results only for h < t, thus for sufficiently small V ⊥ k 2 0 . Therefore, in the next sections we complement the previous RG predictions with the MPS results spanning over a broad range of V ⊥ .

IV. NUMERICAL RESULTS
To obtain more quantitative predictions about the behavior of the model, we simulated the ground state of Hamiltonian (1) for hard-core bosons and system sizes on the order of typical experiments, i.e. L = 64. The simulations were performed through density matrix renormalization group, adopting a matrix product state (MPS) ansatz with open boundary conditions and bond dimensions varying up to 512 different states being kept in the approximation of the reduced density matrix. This leads to a discarded probability of ∆ρ < 10 −8 , enough to assume converged values roughly up to the eighth digit.
In the following we present the main features of the system, including the estimates of the chiral current its main fluctuations, correlation functions and an overview of some dynamical properties.

A. Chiral current
We focus on values of the density of hard-core bosons and flux per plaquette in proximity to the ν = 1 (for example N = 48 and χ = 3π/4 for L = 64) and we vary the interactions V ⊥ and the flux χ. In particular, we vary the flux in units of 2π/L because all the observables in the open system display strong and regular oscillations when continuously varying χ. The chosen discretization corresponds to the standard unit quantum of momentum for translationally invariant chains. We observe however that other choices are equally well justified: similar results are obtained by slightly larger/smaller flux variations 2π/(L ± 1) which respectively correspond to the quantum of momentum and the flux quantum in ladders with Dirichlet boundary conditions.
The chiral current j c is a crucial observable to detect the onset of chiral phases of the ladder. In the case of non-interacting bosonic systems, it shows the typical linear dependence on the flux in the one-dimensional analog of the Meissner phase followed by a shark-fin transition to a vortex phase [1,10], of which we still spot the decaying tail on the left sides of Figs. 3 and 4. The corresponding operator reads: Here the chiral current is defined based on the gauge choice in Eq. (3) and we consider its average over the whole chain length. Fig. 3 depicts j c for N = 48 as a function of the flux displacement around the resonance at filling factor ν = 1. For values of V ⊥ < 6t (not shown) the chiral current does not display any discontinuity with respect to the vortex phase. This is compatible with the system being in a Luttinger liquid phase, consistently with the second-order RG results close to K c = K s = 1 (see Fig. 2). When increasing V ⊥ in the approximate range (6t, 8t) the system develops a high peak in the chiral current exactly at the resonance with filling factor ν = 1 (blue and green curves). We interpret this peak as a feature due to the system evolving towards the critical region at the edge between the h and g phases. By increasing further the interaction, at V ⊥ = 10t, the current signature flattens again. Finally, for 10t < V ⊥ < 18t the chiral current develops a double cusp pattern, which is a hall mark of the two commensurate-incommensurate phase transitions [58,63,64] that separate a (partially) gapped commensurate and chiral phase at the resonance (the h phase) from the truly gapless vortex phase that dominates when the flux χ is sufficiently displaced from the resonant point.
We observe that, for the system sizes we analyzed (L = 64), the chiral current cusps are considerably more evident than the analogous cusps appearing for the pretopological ν = 1/2 Laughlin-like states at χ = 4k 0 and χ = 2π − 4k 0 (the latter being the particle-hole symmetric state of the Laughlin-like state). This is in agreement with previous works [11,32,34], in which the fully resolved Laughlin-like gaps are presented for both larger interactions and larger system sizes. Contrary to the weak and ambiguous signal in the proximity of ν = 1/2, we consistently observe three aligned j c points across the resonance (see Figures 3 and 4), which indicate a sizable gap. This allows us to conclude that the main signatures of the ν = 1 states are more feasible from an experimental viewpoint.
The behavior of the chiral current in the h phase can be deduced by the relation: in which E GS denotes the ground state energy. If the interaction associated with amplitude h is the dominating term, a gap in the spin sector is formed and the chiral current follows a linear behavior with respect the flux χ between the two cusps and is predicted to cross 0 exactly at the resonance. This is true provided there are no background effects such as the long decaying tail from the strong shark-fin signal at small fluxes. In reality, the resonance at ν = 1 is thus established on top of a slowly decaying background current which is a remnant of the vortex phase. We verified that, for different particle numbers, the position of the two cusps of the chiral current shift accordingly to the resonance following the relation χ = 2k 0 (see Fig. 4): the signature of a chiral phase at integer resonance appears indeed for a broad range of particle numbers and its features are visible as long as the ν = 1 res-onance is sufficiently separated from the Meissner phase of the system.
We point out that, in our simulations, the typical double cusp pattern of the resonant state is mostly evident for rung interactions below V ⊥ = 18t when Ω = 0.05t. For larger values of V ⊥ the Meissner and melted-vortex states dominate more than half of the phase diagram accessible through tuning the flux, ultimately resulting in a poor resolution of the ν = 1 signal.

B. Fluctuations and estimate of the Luttinger parameters
To better compare the numerical results with the RG predictions, we must locate the state of the system as a function of V ⊥ in the phase diagram obtained by the RG equations (see the colored dots in Fig. 2). To this purpose we must estimate the Luttinger parameters K c and K s from the DMRG results. The Luttinger parameters determine most of the features of the ground state of the system, including its spin and charge fluctuations, the decay of its two-point correlation functions and even the dynamics of its excitations.
The first estimates we derive for the Luttinger parameters are obtained by fitting the bipartite charge and spin fluctuations [65]. In particular, we define where we introduced the total number of particles / the total magnetization N c/s in a bipartition of size . In bosonization, the leading order of these fluctuations read and depend thus on the correlations of the density fields θ c/s . If the corresponding sector is gapless, bipartite fluctuations follow a logarithmic dependence (see Appendix C): In the above, d denotes the chord distance of open boundary systems which is defined as The function f c/s accounts for the total charge/spin fluctuation of the system, which is related to the U (1) gauge transformations in the corresponding sectors: if the total charge/magnetization is conserved, f c/s = ε (small constant); otherwise, f c/s can be approximated by a linear contribution in the bipartition length, f c/s = γ + f 0 , which equally distributes the total fluctuations on each site. From the fits of the DMRG results we derive that the value of K c remains essentially constant as a function of the flux χ [see Fig. 5 (a)]. This is the expected behavior across the Meissner-vortex phase, since the Meissner phase presents a gap in the spin sector only. When considering the integer resonance at χ = 2k 0 , the independence of K c from the flux χ suggests that also in this case the physics at the resonance is not determined by the charge sector. The situation is the opposite for the spin fluctuations. In the Meissner phase the spin fluctuations are maximized, the linear term given by f s in Eq. (17) is completely dominating over the logarithmic contribution and K s ceases to be meaningful in this gapped phase (its estimate based on Eq. (17) is systematically wrong). For large values of the fluxes in the vortex phase, K s 1, consistently with having repulsive rung interaction. Finally, exactly at the resonance, we can observe a dip in the fitted value of K s [see Fig. 5 (b)], which clearly shows that the spin fluctuations are suppressed at the resonance, compatible with the formation of a new gap in the spin sector. This is the first indication that, at ν = 1 and sufficiently large V ⊥ , the system enters the h-phase, whereas the g-phase, expected to have a gapped charge sector, is not reached. Interestingly, we find also a signature of the pinned spin sector by fitting the constant f 0 in the fluctuations of the total magnetization: it is large in the Meissner and small, but significantly increasing with respect to the vortex environment [see the inset in Fig. 5(c)], inside the resonant state at filling n = χ/π, thus indicating relevant corrections to the logarithmic law.
Considering the system exactly at the resonance, we can therefore fit the values of the Luttinger parameters as a function of V ⊥ based on Eq. (17). The results are used to define the dots of the physical trajectory presented in Fig. 2. The values of V ⊥ and of the fitted Luttinger parameters can then be adopted as initial conditions of the RG flow determined by the equations (11). The numerical solutions of these flow equations allows for an estimate of the gap of the system based on the final value of the flow parameter: ∆ ≈ t * e −l * (see the previous section). The so-obtained gap is represented in Fig. 2b as a function of V ⊥ and it is in striking agreement with the previous analysis of the measured chiral current: the gap is negligible for V ⊥ < 6t and it opens in correspondence of the first peak of the chiral current at V ⊥ ∼ 6t, 8t. This suggests that the system lies in the Luttinger/vortex phase (or inside the h-phase but very close to the boundary with the gapless region) until it reaches the region between the h and g phases for V ⊥ ∼ 6t, 8t. For stronger interactions, the predicted gap ∆ remains sizable, consistently with the system entering deep within the h-phase.
We emphasize that our numerical results and the estimate of the spin correlation functions (see the next subsection) suggest a tiny physical gap, such that the related correlation length ξ s is comparable with a large portion of the system size (L = 64). This implies a small nonuniversal constant which should enter as a multiplicative factor in the estimate of the gap presented in the inset of Fig. 2 (b). As such, the inset describes the general trend of the gap in arbitrary units.
We conclude our analysis of the suppression of the spin fluctuations by considering its implication on the rung current. When the θ s field becomes semiclassically pinned due to the opening of the gap in the h-phase, the rung current must be suppressed as well. For systems with open boundary conditions, however, the rung current has a non-zero value at the edges and it is expected to decay in its bulk. This is what we observe in our results [see Fig. 6 (d)] where the rung current is slightly suppressed in the bulk, compared to the vortex phase.

C. Correlations
In order to better characterize the gapped phase at the resonance, we study the decay of several two-point correlation functions. Indeed, the correlation functions allow us to easily distinguish the Meissner and vortex phases, and provide further indications that the system develops a gap in the spin sector at the integer resonance. In particular, the following set of observables is studied in Fig. 6: We rewrite them in bosonized form and find (23) in which f s/c and g y,y denote the corrections to the leading order. In particular, the leading order is obtained by separating the p = 0 contribution from higher harmonics p = 0 in Eq. (6). It is most convenient to plot the above correlation functions as a function of the renormalized chord distancẽ which highlights the leading order decay caused by the momentum field expectation values. In particular, the connected part of C s shows the expected behavior: it follows an algebraic decay (C s ∝ d −1/Ks ) in the vortex phase; it decays exponentially in the Meissner phase, in which the related unconnected correlation rapidly saturates to a constant due to the pinning of ϕ s (the exponential correction depicted in blue in Fig. 6 (c) is due to the subleading terms f s ); finally, it displays an exponential tail in the h-phase at filling ν = 1 due to the gap in the spin sector and the pinning of θ s [see the bending in the red curve in Fig. 6 (c)].
In contrast to C s , C c decays algebraically everywhere (C c ∝d −1/Kc ), with no signature of any exponential tail [ Fig. 6 (a)]. This is expected as the charge sector remains always gapless.
Finally, the intra-wire Green's function G y,y shows a clear exponential tail in the h-phase [ Fig. 6 (b)], in agreement to our statements above. In the Meissner phase instead, we find G y,y ≈ G y,−y ∝d −1/(4Kc) which is predicted by the leading order of Eq.(23).

D. Dynamics and velocities
To further probe the validity of our RG predictions, we extract the velocities of the Luttinger liquid through the time evolution of spin and charge excitations, which we can then pair with the outcomes of the Luttinger parameters extracted from the bipartite fluctuations.
The velocities can be approximated by tracking the propagation of a local excitation in time. For convenience, we choose to generate a hole at the central site of the lattice: at time t = 0 we apply b L/2,↑ to the (static) ground state |GS and evolve it in time, i.e. we simulate and measure the local density / magnetization at all times n c/s (x, t) = ψ(t)|n c/s (x, t = 0)|ψ(t) . By fitting the propagation of the wavefront in n c/s (see Fig. 7), we obtain the estimates of u c/s presented in Fig. 8. In the charge velocities, we do not expect to find a difference between vortex and integer fillings, and indeed the simulations overlap. On the contrary, we find two different curves distinguishing between non-Meissner and Meissner phase. This behavior is expected as the pinning of the ϕ s field in the Meissner phase results in a renormalization of the effective charge velocity.
An analog reasoning holds for the comparison of the spin velocities in the vortex phase and ν = 1 resonance [grey and red lines in Fig. 8(b)]: O h is irrelevant for small values of the rung repulsion V ⊥ such that the two u s velocities coincide for small interactions. However, for ν = 1 [red line in Fig. 8(b)], u s displays a clear reduction beyond a critical value of V ⊥ ≈ 6t at which it separates with the corresponding velocity in the vortex phase. We interpret this as a renormalization effect and it is a further evidence of a transition between the vortex/Luttinger phase and the helical h-phase for V ⊥ ∼ 6t.
The estimates of u q and K q allow us to evaluate the superfluid stiffness in response to an infinitesimal variation of the vector potential. Due to the Galilean invariance of the interaction V ⊥ and the small values of Ω, the bare values of the superfluid stiffness are roughly constant in both charge and spin sectors: The spin stiffness is expected to be heavily renormalized by the interactions (which primarily act in the spin sector), and we focus in the following on the extraction of the charge stiffness only. Panel (c) of Fig. 8 refers to the quasi-particle propagation in the charge sector for an interacting system. For weak interactions, V ⊥ t , the data display an approximately constant charge stiffness for the vortex and resonant phases, consistently with a Luttinger liquid phase [see panel (d)]. On the contrary, in the Meissner phase the stiffness is renormalized and reduced also for weak interactions. By increasing the interaction, we see that the vortex and resonant states behave differently starting from V ⊥ /t 6. This is a further confirmation that a gap in the spin sector is opened for larger interactions, consistently with the onset of the h-phase.

V. CONCLUSIONS AND PERSPECTIVES
Ultracold atoms hopping in ladder geometries and subject to artificial magnetic fluxes are known to generate rich phase diagrams. For commensurate values of the ratio between the number of fluxes and the number of atoms, helical states with a net chiral current may appear. The simplest and most evident example of these helical states are the non-interacting Meissner phase for bosons and the helical state at flux χ = 2k F for fermions. It is known, however, that additional strongly correlated helical states originate for suitable values of the filling factor ν and suitable interactions [31][32][33][34][35][36][37][38][39].
In this work, we argued that a two-leg ladder of hardcore bosons at ν = 1 is characterized by one of these strongly correlated helical phases of matter, which we called the h-phase. We studied its signatures in terms of correlation functions, fluctuations and dynamical evolution.

Interactions
Laughlin ν = 1/3 Our DMRG simulations show that the strongly correlated ν = 1 helical h-phase can be accessed in systems with contact interactions only, similarly to the bosonic Laughlin-like state at filling factor ν = 1/2 [32]. With respect to the pretopological Laughlin-like state, however, the chiral current and gap signatures we observe are considerably stronger for a broad range of parameters when comparing systems with the same particle density, repulsive interaction and interleg tunneling. This is particularly relevant for experiments based on bosonic ladders like the Rb gases studied in the experiments in Refs. [8,9].
Comparing our findings with the analogous fermionic systems, we also observe that the strongly-correlated phases of bosons can be reached through interactions with a shorter range than their fermionic counterpart, as common for several fractional quantum Hall states. This intuitively explains also why the signals we detect for bosons at filling ν = 1 are considerably larger than their fermionic counterpart (a pretopological K = 8 fractional quantum Hall state) at filling ν = 1/2 [35]. In this respect, we find that bosonic systems are more suitable for the experimental characterization of these helical and strongly correlated phases of matter (see the summary in Table I).
The onset of the bosonic h-phase is caused by a particular mechanism that allows the interaction O h acting on the spin sector of the theory to become relevant despite the presence of repulsive interactions only, which normally suppress it. This happens for bosons at the integer filling factor ν = 1 and it is analogous to the pairing mechanism determining the appearance of the pretopological K = 8 phase in fermionic systems at ν = 1/2 [35]. We stress that engineering (strong) pairing mechanisms analogous to the interaction O h is recently at the focus of several proposals, due to its potential for the development of topological and strongly-correlated paired phases of matter [40,66]. The combination we adopt of artificial gauge potentials, equivalent to a spin-orbit coupling, and a commensurate particle density, can be adopted also in this framework.

ACKNOWLEDGMENTS
The authors wish to thank R. Citro  The field theoretical description of the hard-core boson ladder model has been developed based on bosonization techniques. In particular, the analysis of the behavior at the resonance χ = 2k 0 requires a careful approach since there are two commonly used approximations that cannot be applied in this case. The first is related to the use of several harmonics to map the creation and annihilation operators of the lattice model into the low-energy description in the continuum: the resonance at χ = 2k 0 appears evident only when taking into account higher harmonics, as in the case of the second-incommensurability effects studied in [43][44][45]. This is well-known in the analysis of the one-dimensional limits of fractional quantum Hall states (see, for example, [11,[31][32][33][34][35][36]), where higher harmonics can also be interpreted as multi-particle interaction processes appearing in higher orders of perturbation theory [67][68][69]. The second important characteristics is that, for bosons at ν = 1, there are two of these multi-particle processes that resonate and compose O g in Eq. (9). These have exactly the same scaling behavior under the RG flow and hinder the possibility of adopting a simple semiclassical analysis. This situation is analogous to the case of fermionic ladder models at ν = 1/2 [35] and requires a second-order renormalization group analysis to be examined. In this context the role of O h in Eq. (7) emerges and plays a crucial role, since it mixes and competes with the mentioned g-terms.
In the following we will first present some of the details related to the derivation of the effective Hamiltonian (7), then we will discuss the main steps to derive the RG equation in (11) and how we deal with them numerically.

Bosonization of the hard-core boson model
We mentioned that the analysis of the resonant states appearing at specific values of the filling factor ν requires a bosonization of the lattice operators in terms of the series expansion of vertex operators in Eq. (6) (see, for example, [60]). This series expansion relies on the non-universal coefficients β p which are difficult to evaluate in non-integrable models like ours. Here we first obtain an estimate of the most relevant of these nonuniversal parameters in our hard-core bosonic model, β ±2 , via a Jordan-Wigner transformation and the standard bosonization of fermions. This is rigorous when considering separate chains, i.e., in the limit Ω → 0 and V ⊥ → 0. Then we briefly discuss our choice of setting them to such value, β ±2 = 1/2, irrespective of Hamiltonian parameters.
Hard-core boson operators b ( †) can be related to onedimensional fermionic operators c ( †) via the Jordan- Since the total density of particles, n 0 = N tot /L, coincides for both bosons and fermions, the Fermi momentum is given by k 0 = πn 0 /2. The standard bosonization of the fermionic operators (see, for example, [58,60]) reads: where θ y and ϕ y are two pairs of dual fields obeying commutation relations (Θ being the Heaviside function): and we did not introduce any Klein factor, due to the bosonic nature of our particles which imposes operators in the two different legs to commute. The density of each pseudo-spin species can be approximated at first order in the harmonics expansion as: The symmetrized form of the Jordan-Wigner string of Eq. (A1) reads then: (−1) j<x n(j,y) = e iπ j<x n(j,y) + e −iπ j<x n(j,y) 2 Substituting the expressions (A2) and (A5) into Eq. (A1), we obtain: b x,y ≈ k 0 2π 2 e iϕy(x) + e −2ik0x 2 e i(ϕy(x)+2θy(x)) + e 2ik0x 2 e i(ϕy(x)−2θy(x)) .
In the limit Ω → 0, space-inversion symmetry implies indeed that β p = β −p . When considering Ω > 0, however, it relates the β p and β −p coefficients of the bosons of different legs. Since we focus in a regime with small values of Ω/t, we assume hereafter that these coefficients are the same for both the bosonic species, such that β 2 = β −2 . Even in case of more general coefficients, however, the analysis we present in the following is not qualitatively affected.
Concerning a precise evaluation of these coeffcients, the estimate of the β p parameters for general values of the interactions, the transverse hopping Ω and the density of the system remains an open problem. Analytical solutions and numerical approximation are known only for (several) integrable models (see, for example, [61,62]).
Our model, in particular, can be mapped into a fermionic Hubbard model with interaction V ⊥ > 0 in the limit Ω → 0 (see also a similar analysis in [11]). For the fermionic Hubbard model in the limit V ⊥ → ∞, the density modulations with momentum 2k 0 are suppressed [58], such that we expect β ±2 to be equally suppressed for large V ⊥ . As we will discuss below, setting the value of β 2 might have important implications in the determination of the phase diagram through the RG equations (11), due to the dependence of the bare values in (10) on β 2 . However, since a precise estimate of β 2 for large interactions is difficult to derive and beyond the scope of our work, we decide to approximate it with the constant 1/2 irrespectively of the Hamiltonian parameters. This provides results compatible with the numerical simulations and we checked that the phase diagram in not affected by small variations of this parameter.
At this point, we are set to derive the g− and h-terms of the effective model in Eq. (7). This is most conveniently done by using the charge c and spin s sectors of the model, as introduced in Sec. II: We begin our analysis from a gauge-equivalent formulation of the non-interacting Hamiltonian of Eq. (3): Although this gauge choice breaks the translational invariance, it turns out to be more convenient to make resonant terms evident. By plugging in Eq. (6), the interleg hopping term reads indeed: +H.c. (A9) where we did not explicitly write additional fast oscillating terms (F.O.) which are not resonant for χ = 2k 0 , i.e., ν = 1. A direct comparison with Eq. (7) yield the initial value of the coupling constant g = β 2 k 0 /2π 2 in the RG flow, reported in Eq. (10).
Coming to the interacting part of the model, we need to supplement Eq. (A4) with additional higher-order oscillating terms [60]. The first of such terms can be estimated from a point-splitting procedure, i.e., by evalu- In this way, we obtain the correction : (A10) The interspecies interaction becomes then: where we neglected fast oscillating addends, boundary terms (such as ∂ x θ c ) and less relevant ones (such as (∂ x θ c ) cos(2 √ 2θ s )). The initial value h = V ⊥ β 2 2 k 2 0 /2π 2 for the RG flow, reported in Eq. (10), is then readily identified. The first line of Eq. (A11) is responsible for giving the initial condition for the Luttinger parameters K q of Eq. (12).

Renormalization group equations
The second-order renormalization group analysis of the Hamiltonian (7) follows the analogous fermionic case at filling ν = 1/2 [35]. We apply, in particular, the Wilsonian RG in momentum space at second order. Here we derive the corresponding flow equations (11) in detail for the resonant regime, i.e., χ = 2k 0 .
It is crucial to anticipate here that we consider the interaction terms h and g in Eq. (7) as a perturbation of the gapless Luttinger liquid. Therefore, we expect our results to be valid when the bare values of g and h are much smaller than 1. In particular, by assuming Ω = 0.05 and N/L = 48/64 as in the numerical simulations, we obtain that the general features of the phase diagram presented in Fig. 2 (b) are stable until V ⊥ < 5t (cf. Fig. 2 (a)). Beyond this threshold non-physical features emerge which are not compatible with the numerical simulations, thus indicating that a second-order perturbation approach fails beyond this limit.
For each of the bosonic fields we distinguish fast and slow modes, separated by an effective cutoff in momentum space that we labelΛ. Furthermore, we introduce an ultraviolet momentum cutoff Λ >Λ. The fast oscillating modes are characterized byΛ < k < Λ and we are interested in the limit Λ/Λ = 1+dl, with dl infinitesimal. The bosonic fields (q = c, s) can thus be decomposed into: Moreover, we resort to Euclidean space, i.e., we use coordinates z = (x, τ ≡ it) and we denote d 2 z = dx dτ , such that the following duality relations hold: We separate the action of the resonant model into a Gaussian and an interacting part, including both the g− and h−terms of Eq. (7), and we consider the latter as a perturbation: Our aim is to obtain an effective action for the slow modes only, by integrating out the fast degrees of freedom: where the expectation values are taken on the (fast) Gaussian action only, and we identified the effective action at the second order of perturbation theory. In the following we make extensive use of the following wellknown key property of Gaussian integrals (G): with φ the fields over which the integration is carried out. The Gaussian correlation functions for the fields in first order of ln(Λ/Λ) read: Here the logarithm captures the scaling behavior, and C q (r) is a function of r = u 2 q (τ 1 − τ 2 ) 2 + (x 1 − x 2 ) 2 , such that C q (0) = 1. In the following we will consider C q (r) to be suitably short-ranged; in the case of a sharp cutoff, C q (r) = J 0 (Λr) and the Bessel function J 0 does not satisfactorily fulfil this assumption, but C q (r) can be made sufficiently short-ranged with more refined cutoffs.
In particular, such an optimization can be achieved by a deformation of the integration Λ Λ dp → ∞ 0 dpf n (p, Λ) with f n = Λ n /(p n + Λ n ), in which the sharp cutoff is realized for n → ∞. An analytic result can be obtained in case of n = 2, i.e. C q (r) = ΛrK 1 (Λr) with the modified Bessel function K 1 of the second kind. This function has an exponentially decaying asymptotic form zK 1 (z) ≈ zπ/2e −z which is thus sufficiently short-ranged. To be more precise, the numerical values of zK 1 (z) above machine precision are confined to the interval 0 ≤ z < 13π, and moreover K 1 (2π) ≈ 0.0062 is already pretty small for short distances of r ∼ 2π/Λ. For this reason, we will later truncate the space-time integrations of a product containing the function C q to the interval (0, α) in which the cutoff parameter α ∼ 2π/Λ = a is on the order of the lattice spacing.
For the sake of convenience, let us rewrite the interacting part of the action as: where we introduced the shorthand notation: and O ν g,µ = e iν √ 2(θc+ϕs+µθs) . (A21) Thus we obtain for the first-order contribution: where the operators on the right hand sides of the equations are intended to be constructed with slow fields only, and we dropped extra-labels for simplicity. Together with the factor (Λ/Λ) 2 from the rescaling of the integration domains for the new action after a RG (infinitesimal) step, and keeping in mind that Λ/Λ = 1+dl, Eqs. (A22)-(A23) return S I with the usual first-order dependence of the couplings g and h on the scaling dimensions, D h = 2K s and D g = K c + K s + K −1 s /2, in the RG Eqs. (11)). The many addends of the second-order contribution could be compactly written as: where we introduced the connected part of the quadratic expectation values, i.e., It is straightforward to obtain which RG correction each addend gives rise to (hereafter, we drop the integration symbol over the Euclidean space, for the sake of space): where we indicated the dependence on the Euclidean coordinates z 1 and z 2 by simple labels, and we took into account the additional Baker-Campbell-Hausdorff phase factor arising in Eq. (A28).
We can now split the integrals into the integral of the center of mass coordinate (z 1 + z 2 )/2, and the relative coordinate z 12 = (z 1 − z 2 ), and recall the assumed short-range character of the C q functions -see the pre-vious discussion, just below Eqs. (A18)-(A19). This constrains the two points to be in a neighbourhood of size δz (α, α/u q ) from each other, where α is of the order of the lattice spacing a. The effect is two-fold: first, this allows us to Taylor-expand the field operators in Eqs. (A25)-(A28); second, a factor α 2 /u q is generated by the integration. Henceforth we group the addends of Eq. (A24) according to their RG behaviour: • The terms with ν = −ν in Eq. (A25) contribute to the quadratic part of the action: as well as those with ν = −ν and µ = +µ in Eq. (A28): Kc uc + 1 where we dropped two terms mixing ϕ s and θ c , thus the spin and change sectors. Such terms break the possibility of rigorously splitting the system into these separate sectors and would require instead to introduce a unitary matrix, dependent on l, to diagonalize the Gaussian part of the action at each scale. For the sake of simplicity, however, we neglect this mixing and we consider only the (standard) correction to the Luttinger parameters K s and K c ; • the terms with ν = −ν and µ = −µ in Eq. (A28) generate a correction to O h : • finally, the terms with ν = −νµ in Eqs. (A26)-(A27) give rise to a contribution proportional to O g + O † g : • all other terms are considerably more irrelevant.
In our numerical solutions of the RG equations, we checked that small variations of the non-universal parameter α around the lattice spacing a do not qualitatively affect the phase diagram in Fig. 2 (a), consistently with having a sufficiently sharp cutoff function C q .
Finally, we point out that, differently from the standard sine-Gordon interactions, the terms in Eq. (A30) generate also a non-trivial flow of the velocities. In particular, there are second-order corrections to u c and u s proportional to dlg 2 u 2 c − u 2 s /u 2 c . In the regime with a bare value g t and small interactions, the flow of the velocities is expected to have only negligible effects and we do not take in consideration their flow. To observe a quantitative difference between the bare and physical values of the velocity, we can compare the gray (vortex phase) and red (ν = 1 resonance) curves in Fig. (8).
In the vortex phase the system is described by a standard sine-Gordon model which fulfills Lorentz invariance (thus the velocities do not flow). Only a small difference in u s can be observed when comparing the red and gray curves, thus confirming that the flow of the velocities is practically negligible in our regime of interest. In the following, we provide additional information about the numerical solutions of the renormalization group equations presented in Eqs. (11). Fig. 2 is obtained by explicitly solving the differential equations using an implict Runge-Kutta method, in particular by fixing the initial/boundary conditions of Eqs. (11). In this four-dimensional phase space, we put on the axes the Luttinger parameters before the RG flow K q ≡ K q (l = 0), while fixing g/h(l = 0). We thus manage to present two-dimensional cuts of the general phase diagram which display in which regions of these bare parameters the operators are relevant. Representatives of the solutions of the RG equations for Ω = 0.05t and V ⊥ = t are presented in Fig. 9. The flow of the Luttinger parameters in the MPS explored region is in general slow, and the bare values roughly correspond to the renormalized ones at l = l * (see panels (c) and (d) presenting the flow in the h-phase). This way we can conveniently combine the phase diagram and the path of our MPS simulations. In particular, feeding the extracted Luttinger parameters from the MPS simulations as initial conditions, we confirm by estimating the energy gap ∆ presented in Fig. 2 (b) that the RG predicts a flow inside the h phase in the thermodynamic limit, after a critical value of V ⊥ ≈ 6t. It must be stressed that the estimate of the RG predicted energy gap are given in arbitrary units, and, being unable to resolve it with the MPS simulations, its quantitative size remains unknown.
By observing the different points of the phase diagram presented in Fig. 9, we observe a Luttinger liquid in which the couplings flow very fast towards the lower threshold The h-dominated phase is slightly refined: we note the expected exponential growth of the coupling while the Luttinger liquid parameter Ks flows towards zero. On the contrary, the g-dominated phase is drastically refined, as the Luttinger parameters show sudden divergences (see panels (e), (f), and an explanation in the text), causing the couplings to diverge accordingly. These divergences are followed by a non-monotonicity in the growth of the couplings which is better visualized by gradually moving the initial values of K c/s from the Luttinger towards the g-dominated phase, a path we present in Fig. 10.
τ and the Luttinger parameters remain constant [panels (a) and (b)]. In the h-dominated phase [panels (c) and (d)], we always find a fast exponential growth of the couplings, hitting the upper threshold at l = l * < 6. On the contrary, the entire g-phase suffers from divergences due to a divergent K s [see Fig. 9 panels (e) and (f)]. These divergent solutions of the RG equations can be qualitatively understood by assuming that K s reaches a large value K s 1, K c during its flow: in this limit, we can approximate the dominant part of its RG equations with dKs dl = aK 3 s , which upon integration becomes K −1 s = √ c − 3al, diverging at l = c/3a. Interestingly, the second order RG equations result in a non-monotonic curvature of g(l) and h(l) in parts of the Luttinger liquid phase [see 10]. In this cases, the Luttinger parameters reach asymptotic values such that the coupling constants h and g decay to zero after a nonmonotonic flow. This happens in particular in the portion of the Luttinger liquid phase where g is predicted to be relevant according to a simple first order analysis. By examining the flow in proximity of the g-dominated phase, we observe that h flows to negative values before decaying towards zero. Naturally, this non-monotonicity of the coupling constants connects smoothly to the reported divergences of the g-dominated phase in which the asymptotic value of K s appears to be infinite. In conclu-sion, the divergence is not caused by numerical errors, but it is rather a unique feature of the second order RG equations in proximity of the g-dominated phase. This also suggests that a more thorough study of the intricate g-phase must include third order corrections.
Concerning the velocities of the system, their value is considered to be approximately constant throughout the RG flow and has been estimated from the relation of Eq. (26), resulting in v q = v 0 /K q . This relation, however, is supposed to be affected by large deviations for large values of the interaction V ⊥ , as suggested by the behavior in Fig. 8 (b) which shows u s asymptotically reaching a minimum value ∼ u 0 /4 for large interactions. The failure of the approximation v s = v 0 /K s for large interactions concurs to increase the errors in the RG flow predictions for this strong interaction regime. In particular, we see that in the off-resonant case g = 0, the RG flow erroneously predicts a new onset of the h phase for V ⊥ > 6t at large values of K s . This is a combined effect of the limitations of the second-order perturbation theory and the approximation v s = v 0 /K s , which lead to dK s /dl ∝ −h 2 K 5 s at the beginning of the flow. In this specific regime of large interactions, this decrease of the Luttinger parameter K s appears to be too fast and yields a non-physical reappearance of the h-phase for large V ⊥ which is not observed anywhere outside of the resonance . By extending the first order RG approximation to second order, we note that the Luttinger liquid phase grows drastically larger, reaching far inside the region predicting the relevance of g at first order. This Luttinger liquid phase is not trivial, as the Luttinger parameters flow to asymptotic values significantly different from the initial ones, followed by a rapid change of the coupling amplitudes. We observe thus a subtle non-monotonic behavior of the coupling constants, which smoothly extends towards the g-dominated region. Towards the g-dominated region the asymptotic value of Ks approaches infinite. We report thus a sudden divergence of Ks inside the g phase (see Fig. 9 (f)), demonstrating the breakdown of the second order RG approach as indicated by the shaded regions in the phase diagrams.
in our numerical simulations. Based on the numerical solution of the RG equation, we conclude that our RG results are reliable only in the range V ⊥ < 5t when addressing large values of the Luttinger parameter K s .

Appendix C: Further details about the correlation functions
The calculation of the correlation functions of the θ and ϕ fields for finite systems depends on the boundary conditions. In principle, we should separately consider the charge and spin sectors. For the charge sector indeed, no particle can tunnel in or out of the system from the boundaries, and the charge current must exactly vanish at the boundaries. Therefore ∂ x ϕ c = 0 for both x = 0 and x = L. For the spin sector, instead, the total spin density is not conserved and, although the chiral current must vanish in average at the boundaries, the previous Neumann boundary conditions for ϕ s is too restrictive, as demonstrated by the non-vanishing rung current at the edges of the system (Fig. 6 (d)). Hence, one should consider more general kinds of boundary conditions [70][71][72][73]. We expect, however, the deviation from the Neumann boundary conditions to be small (proportional to Ω/t), and, also for the spin case, we will consider ∂ x ϕ s = 0 at the edges. Based on the construction in [60] (see also [70,71]), we can define the charge and spin fields (at time t = 0) in the form: where labels the sum over the even integers k, with k = 0 from −∞ to +∞. Here b q,k labels two sets of bosonic annihilation and creation operators related to the spin and charge sectors respectively. In order to fulfil the Neumann boundary conditions, we must impose: From the previous relations we derive: Considering the case x → 0, and neglecting functions of x only, we obtain the behavior of Eq. (17) for the spin and charge fluctuations. Given the deviation from the Neumann boundary conditions of the spin sector, however, we must account for an additional weak space dependence represented by the function f s ( ).