Spin-gap spectroscopy in a bosonic flux ladder

Ultracold bosonic atoms trapped in a two-leg ladder pierced by a magnetic field provide a minimal and quasi-one-dimensional instance to study the interplay between orbital magnetism and interactions. Using time-dependent matrix-product-states simulations, we investigate the properties of the so-called"Meissner"and"vortex"phases which appear in such system, focusing on experimentally accessible observables. We discuss how to experimentally monitor the phase transition, and show that the response to a modulation of the density imbalance between the two legs of the ladder is qualitatively different in the two phases. We argue that this technique can be used as a tool for many-body spectroscopy, allowing to quantitatively measure the spin gap in the Meissner phase. We finally discuss its experimental implementation


Introduction
Orbital magnetism (OM) encompasses a host of phenomena that arise in systems of charged particles subject to an applied magnetic field. Because the Bohr-van Leeuwen theorem forbids its appearance in an ensemble of classical particles [1,2], OM is a trademark of quantum mechanics since its early days. In the case of electrons in solids, for instance, OM effects include Landau diamagnetism [3], and the integer and fractional quantum Hall effects [4,5].
Flux ladders (FL) composed of two (or more) coupled one-dimensional subparts with a magnetic field perpendicular to the ladder plane are among the simplest setups where OM can appear. FL are quasi-one-dimensional, and thus still amenable to an efficient theoretical treatment in the presence of interactions, either using bosonization [6] or numerical methods based on matrix-product states (MPS) [7,8]. Establishing the connection with two-dimensional physics for studying FL is one of the major motivations in this research field.
Bosonic two-leg FLs have been particularly studied, in part due to the simplicity of the model, and in part because of the recent experimental realization with ultracold atoms in suitably designed optical lattices [9]. Using the bosonization technique, the pioneering work of [10] predicts the appearance of vortex (V) and Meissner (M) phases paralleling the phenomenology of superconductors. A V phase is characterized by non-vanishing inter-leg ("transverse") current, and a M phase by vanishing transverse current. For strong interactions and commensurate densities, a phase transition between a Mott-insulator (MI) and a superfluid (SF) also appears [11]. According to these fieldtheory treatments of the low-energy part of the model, two-leg ladders feature generally two excitation branches, related to "charge" (or "density") degrees of freedom on the one hand, and to "spin" degrees of freedom on the other [10,11]. The MI phases then correspond to the opening of a charge gap, and the M phases to the opening of a spin gap. All the four situations obtained by combining these two classifications -V-SF, M-SF, V-MI and M-MI -are possible. Numerical studies of microscopic models of interacting bosonic FLs have confirmed the existence of these four phases and more, revealing an extraordinarily rich phenomenology . For instance, it has been proposed recently that precursors of the physics of the fractional quantum Hall effect, and in particular of Laughlin wave functions, might appear in experimentally-relevant bosonic FL [19,24,34,35].
Experimentally, the realization of bosonic FL belongs to a more general effort to realize effective gauge potentials coupling to ultracold atoms in spite of their electrical neutrality [36,37]. The experiment of [9] creates a one-dimensional array of isolated ladders with a total flux per plaquette Φ = π/2 induced by combining laser-assisted hopping with a periodic spatial modulation of the lattice. In this experiment, each site of the ladder is in reality a one-dimensional bosonic gas with many atoms, with the result that the interaction energy per atom was very weak compared to inter-and intra-leg tunneling energies. Recently, the role of interactions in bosonic FL was experimentally investigated for two particles [38].
The experiment of [39] exploits the concept of "synthetic dimension". Each leg of the ladder can be represented by internal (spin) states of the atom, and the magnetic flux is due to Raman transitions coupling the internal states. The idea of synthetic dimension has been recently generalized to momentum space lattices [40]. Importantly, in the synthetic dimension approach, the two legs are not separated in space, but fully overlapping. As a result, interactions are short-ranged in real space, but have almost infinite range along the synthetic dimension. This makes interacting models using the synthetic dimension approach quite different from models with short-range interactions [41,42].
Fermionic flux ladders can also be explored experimentally with ultracold atoms using similar approaches as in the bosonic case [43][44][45][46]. Theoretical studies have highlighted the presence of fractional charge excitations and predicted a host of novel phases of matter (such as charge-, bond-and density-waves or orbital antiferromagnets) leading to a more complex phenomenology than the V-M competition of the bosonic case [47][48][49]. Triggered by the interest in the quantum Hall effect, analogues of the chiral modes which characterize both integer and fractional phases have been discussed [24,34,42,[50][51][52][53][54][55]. Here, J is the tunneling amplitude between nearest-neighbor sites in the longitudinal direction j, J ⊥ is the tunneling amplitude in the transverse direction m, Φ is the gauge flux piercing each plaquette, and U is the on-site interaction strength, taken equal for both legs.
In this article, we propose an experimentally-feasible method to distinguish the M and V phases in the bosonic FL and to characterize their low-energy excitation spectrum. It is known that the M and V phases can be distinguished qualitatively by time-of-flight methods [9,31,32,45,56,57]. We show that they also respond differently to a periodic "spin" modulation, and we interpret our results as a measure of the spin gap in the M phase. We support our claims by presenting numerical simulations performed both in the dilute non-interacting limit and in the dense interacting case. This extends previous work studying dynamical protocols to probe bosonic or fermionic systems in one dimension [58][59][60][61][62][63]. Finally, we show how to adapt the proposal of [64], initially designed to realize two-dimensional systems with an effective magnetic flux, to the realization of FL with strong on-site interactions. This scheme is well suited to the spectroscopic method probing the spin gap, although we note that the method can also be used in other implementations of bosonic FL.
The article is organized as follows. In Section 2, we introduce the model, and in Section 3 we briefly discuss some aspects of its phase diagram. In Section 4, we present our theory for the spin-gap spectroscopy and the numerical simulations supporting our statements. In Section 5, we discuss a possible experimental implementation of bosonic FL using state-dependent lattices and laser-induced tunneling, and discuss how the proposed measurement could be carried out. We finally draw our conclusions in Section 6, and provide some technical details in the appendices.

Model and notations
We consider a gas of interacting bosonic atoms loaded into an optical lattice at zero temperature. The system is a FL composed of two coupled one-dimensional systems immersed in a (possibly synthetic) magnetic field. A sketch of the ladder is shown in Figure 1, where j and m identify the longitudinal and transverse directions of the ladder respectively. Such a system can be modeled by the following tight-binding Hamiltonian including interactions [10]: Here,b j,m (b † j,m ) annihilates (creates) a boson on site j and on the leg m,n j,m =b † j,mb j,m is the local density operator on the leg m, J and J ⊥ denote the tunneling amplitude between two nearest-neighbor (NN) sites in the longitudinal and transverse direction, respectively, and Φ is the magnetic flux per plaquette. The inter-particle interaction is taken into account by the Bose-Hubbard on-site interaction U , which we take equal for both legs. We denote by L the total number of rungs of the ladder, and consider open boundary conditions (OBC). The total number of particles N defines the particle density per rung n through n = N/L. In the Hamiltonian in Equation (1), the gauge flux is set in such a way that the tunneling matrix elements on the transverse links of the ladder are complex, and the longitudinal ones are real. We will refer to this choice as the experimental gauge (ex). It is convenient to make the Hamiltonian in Equation (1) wheren j,m =b † j,mb j,m =d † j,md j,m . The choice of the gauge as in Equation (2) will be referred to as the condensed-matter gauge (cm). In the following, if not explicit, we use J as reference energy scale.
For non-interacting bosons (U = 0), the Hamiltonian in Equation (2) can be diagonalized in momentum space by introducing the operatorsd k,m = L −1/2 j e ikjd j,m . The two energy bands are given by E ± (k) = −2J cos(k) cos(Φ/2) ± 4J 2 sin 2 (k) sin 2 (Φ/2) + J 2 ⊥ . The structure of the lower band E − (k) changes with J ⊥ or Φ. When J ⊥ exceeds a critical value J ⊥,c = 2J sin(Φ/2) tan(Φ/2), the lower energy band has one minimum at k = 0. When J ⊥ < J ⊥,c , the lower band features two symmetric minima at k = ±k M (Φ, J ⊥ ). In the former case, the system is in the Meissner phase (M), whereas it is in the vortex phase (V) in the latter. By tuning J ⊥ and/or Φ, the system can undergo the M-V phase transition [18]. This transition persists for non-zero repulsive interactions, but the critical value J ⊥,c (that depends on U, n, Φ in general) can be strongly modified by interactions [28,31].

Momentum distribution functions and phase diagram of interacting bosonic flux ladders
We numerically study the properties of bosonic FL with repulsive interactions (U > 0) using a MPS-based algorithm [8]. The ground state (GS) of the system is found after a local variational search in the MPS space. At finite U , we keep d loc = 3 states for the local Hilbert space (see Appendix A for details and a critical discussion).
According to bosonization, the M phase is distinguished from the V phase by the presence of a gap appearing in the spin sector of the low-energy theory [10,21] (hereafter denoted as "spin gap"). As a consequence, the two phases differ also in the so-called central charge c that, in this context, roughly speaking gives half the number of gapless modes [6]. When the particle density is less than unity, n < 1 (which is the situation that will be studied in this article), the charge sector is always gapless (no MI phase). The spin sector is gapped in the M phase (thus c = 1 if n < 1), and gapless in the V phase (thus c = 2 if n < 1). Monitoring the change of c with variations of parameters J ⊥ , Φ, n allows one to track the M-V phase transition (see Appendix A). MPS methods are well-suited to extract the entanglement entropy from which the central charge is deduced [65].
Another possibility would be the direct numerical computation of the spin gap that distinguishes the two phases. Such measurement is typically performed in ladder or more general models with two decoupled species, where the number of particles for each species is a conserved quantity [66,67]. In our situation, however, the spin gap can not be accessed directly: only the total number of particles is a conserved quantity when J ⊥ and Φ are both non-zero. As a result, there is no quantum number associated with the spin sector (outside of the low-energy sector). This makes the computation of the spin gap unfeasible in practice. We propose in the next section a spectroscopic method that can be used to estimate the spin gap.
We begin by reviewing a method to study the phase diagram [20,23,28,31,32], which can be easily implemented in experiments [9]. We focus on the momentum distribution functions (MDF), both leg-resolved and total. Time-of-flight measurements readily give access to the total MDF; in some experimental schemes, such as the one discussed in Section 5, is even possible to measure it only for a specific leg. The leg-resolved MDF in the experimental gauge is defined as where the expectation value is computed over the GS of the Hamiltonian in Equation (1). Since the MDF is periodic with period 2π, we restrict the momentum variable to k ∈ [−π : π). By using the unitary transformation introduced before, the MDF in the experimental and in the condensed-matter gauge are simply related by a momentum shift, i.e., n In the condensed-matter gauge, the MDF displays one peak centered at k = 0 for the M phase, and two symmetric peaks at k = ±k M for the V phase, reminiscent of the single or double minimum of the lower energy band when U = 0 [32]. We report in Figure 2 the MDF for J ⊥ /J = 1.50 (panel (c)) and J ⊥ /J = 1.75 (panel (d)) for several values of U and n = 1/2. For sufficiently low values of J ⊥ , the two-peak structure of the MDF is observed for all U . For large enough J ⊥ , the V-M phase transition occurs when U is increased beyond a critical value; in this case, we see the emergence of a third peak at k = 0, which eventually dominates the MDF when one enters the M phase.
To go beyond these qualitative features and to quantitatively distinguish M and V phases, we define the imbalance ratio (IR) The IR takes the values 0 < δn < 1 in the V phase and −1 < δn < 0 in the M phase. We propose to find the transition points by imposing the condition δn = 0. The IR provides a simple and experimentally accessible observable to distinguish V and M phases, although it is not an order parameter in the sense of Landau theory. A more rigorous numerical characterization of the two phases is provided in Appendix A, where we show that, for n = 1/2, the transition point identified by δn = 0 is very close to the point where the central charge introduced earlier changes from c = 2 to c = 1 [21]. By monitoring the variations of the IR with a control parameter, for instance J ⊥ , we can obtain a qualitative phase diagram for the Hamiltonian in Equation (1), and analyze how the presence of interactions affects the critical point at which the V-M phase transition occurs. A similar analysis was discussed in [31]. The phase diagram in the U -J ⊥ plane for a fixed flux per plaquette of Φ = 0.8 π is shown in Figure 2(a) for n = 1/4 and (b) for n = 1/2. Red points correspond to δn < 0 (M phase), green points correspond to δn > 0 (V phase), and the yellow line represents the critical line separating the two phases. We first focus on the case with n = 1/4 ( Figure 2(a)). Repulsive interactions U > 0 shift the critical value of J ⊥ with respect to the non-interacting case J ⊥,c (U = 0) 5.9 J [18]. We find that J ⊥,c (U, n) is a monotonous and decreasing function of U , with J ⊥,c (∞, n) 3.8 J for hard-core bosons (U → ∞) and Φ = 0.8π. For a larger particle density, the shift of J ⊥,c (U, n) is expected to be enhanced further with respect to the n = 1/4 case. The numerical simulation confirms this expectation, as we show in Figure 2(b) for n = 1/2.

Spin gap spectroscopy
In the previous Section, we characterized the M and V phases by looking at the MDF. In this Section, we study the response of the bosonic ladder to a periodic imbalance of the particle number on the two legs, and we show that the system displays different responses in the M and V phases. We interpret our method as a spectroscopic tool that detects and measures the presence of the spin gap in the M phase predicted by bosonization.

Model and observables
We consider the HamiltonianĤ (ex) 0 in Equation (1), and add a time-periodic perturbationV = F (t)N s proportional to the difference of populations between the two legs (hereafter denoted as spin imbalance), withN m = jn j,m the particle number per leg, with F (t) = δ 1 sin(ωt) Θ(t), and with Θ(t) the unit step function. Here, we denote by δ 1 and ω the amplitude and frequency of the modulation, respectively. The total Hamiltonian is thuŝ In what follows, to ease the notation, the superscripts denoting the experimental gauge in the Hamiltonian in Equation (6) will be omitted. We consider the time evolution of the mean energy, where |Ψ(t) is the time-evolved state starting from the GS of the bosonic ladder. We define the energy absorption rate (EAR) aṡ Within linear response theory, the EAR per unit frequency probes the imaginary part of the response function, i.e.ε(ω)/ω ∝ Im[χ Ns−Ns (ω)], where χ Ns−Ns (ω) = ∞ 0 dt e iωt χ Ns−Ns (t) and where is the response function in real time. HereN s (t) = e iĤ 0 t/ N s e −iĤ 0 t/ is the number imbalance expressed in the interaction picture with respect toĤ 0 . Notice that, by means of Equation (7), the EAR can be experimentally accessed by measuring the total spin imbalance in time Ψ(t)|N s |Ψ(t) (see Section 5).
If we denote by ∆E s the value of the spin gap, a spectroscopic method that identifies it should consist of a periodic modulation of the system that is sensitive to its presence, so that the system does not absorb energy as long as ω < ∆E s , and energy absorption can occur only for ω > ∆E s . We thus expect Im[χ Ns−Ns (ω)] = 0 if ω < ∆E s and Im[χ Ns−Ns (ω)] > 0 otherwise.
To compute the response in time to the modulation in Equation (6), we first compute the GS by means of the variational MPS-based algorithm discussed in Section 2. The time-evolved state, |Ψ(t) is computed by using the time-evolving-block-decimation (TEBD) algorithm [8,68,69] with a fourth-order Trotter expansion [70,71] with time step dt (during the time evolution, we fix the maximum bond link D max,t used to describe the MPS state at time t).

Results for dilute gases
We first analyze a very dilute gas (n 1) where interaction effects are weak and the physics is expected to be close to the free case, for which the critical line is analytically known [18]. We focus on the limit of hard-core bosons (U → ∞). Differently from the equilibrium results presented in Section 4.1, the increased numerical complexity of simulating the time evolution forces us to restrict ourselves to smaller values of the system size, namely L = 24.
The results are shown in Figure 3, for L = 24, n = 1/12 and different values of Φ. The amplitude of the density modulation is δ 1 = 0.4 J, and we use different values of ω, depending on the value of Φ. During the time evolution, the time-dependent Hamiltonian in Equation (6) is taken constant within each Trotter step. Therefore, the time step dt has to be chosen small enough to ensure the reliability of this approximation for all the values of ω that we consider. We have verified that we can choose the time step dt = 10 −2 /J in the Trotter expansion (see also Appendix B for a deeper discussion). In Figure 3(a), we show the relative energy variation, ∆E(t) = E(t) − E(0) for the two typical cases. In the M phase (J ⊥ /J = 7.0), there is no net energy absorption for sufficiently small ω, whereas the system absorbs energy for all ω in the V phase    (J ⊥ /J = 2.0). The EAR is extracted from the slope of ∆E(t) represented by the black dashed line. To remove the fast oscillations of ∆E(t) and extract the longtimes linear trend, we perform M linear fits to ∆E(t) using different ranges of t. Accordingly, we obtain a set of values for the EAR per unit frequency, {ε q (ω)/ω} M q=1 , from which we compute the mean valueε(ω)/ω = M −1 M q=1ε q (ω)/ω, and the standard deviation, We take the latter as a measure of the uncertainty on the determined slope.
In Figure 3 (8). This is consistent with the curve plotted in Figure 3(e), which tends to 0 for low values of J ⊥ . A similar behaviour is also expected for the values of Φ used in Figure 3(b)-(d), but the considered value of J ⊥ was not small enough to highlight it.
As we previously stated, in the presence of a spin gap, the system is expected to absorb energy only if ω > ∆E s . In Figure 3(f ), we show the EAR per unit frequency as function of ω for Φ/π = 0.24, both in the V phase (J ⊥ /J = 0.1) and in the M phase (J ⊥ /J = 0.5). For low modulation frequencies, we observe that the system can absorb energy in the V phase for values of the modulation frequency down to ω ∼ 10 −2 J/ . In contrast, in the M phase, energy absorption starts from a finite frequency threshold, the value of which can be considered as a qualitative estimate of the spin gap. For high frequencies ω, one observes a drop of the response in both phases, as expected from the general behaviour of the susceptibility χ(ω) [6].

Results for strongly interacting gases
We now move to the discussion of the strongly correlated case. To approach this regime, we consider hard-core bosons (U → ∞) and higher density with respect to the case in Section 4.2. Previously, in Section 2, we showed how the presence of interactions shifts the critical point for the V-M phase transition. We here demonstrate that this shift is also detected by the periodic modulation of the density imbalance. For concreteness, we focus on n = 1/4 and Φ = 0.8 π, as for the data in Figure 2(a).
The numerical results are shown in Figure 4. In Figure 4(a), we plot ∆E(t) for different values of J ⊥ /J using δ 1 = 0.4 J and ω = 10 −2 J/ . In Figure 4(b), we display the EAR per unit frequency as a function of J ⊥ /J for the same set of data. We use the same smoothing procedure as in the previous Section 4.2. The EAR per unit frequency  vanishes for J ⊥ ≥ 5.0J, and becomes nonzero when J ⊥ = 4.0J and below. In Section 2, we estimated the critical value for the V-M transition, J ⊥,c /J 3.8, slightly lower than the observed threshold for energy absorption. This quantitative discrepancy may be due both to finite-size effects, and to the fact that ω is possibly larger than the spin gap for J ⊥ /J = 4.0.
In panels (c) and (d), we show the same analysis for a larger value of ω, namely ω = 10 −1 J/ . For J ⊥ /J 4.0, the system absorbs energy until saturation starts to take place. Instead, for J ⊥ /J 4.0, energy absorption is suppressed. Differently from the data in panels (a) and (b), we see a nonzero energy absorption also for J ⊥ /J = 5.0, 6.0. We ascribe this fact to the larger value of ω, possibly overcoming the value of the spin gap. The numerical complexity of the problem prevents us to use lower values of ω, as the required simulation times t are beyond our numerical possibilities. For a more critical discussion of the numerical data, see Appendix B.
Concluding, our results are compatible with the opening of a spin gap around J ⊥ /J 4.0, which is in qualitative agreement with the phase diagram presented in Section 2 for hard-core bosons and n = 1/4. We thus conclude that the protocol we propose provides an experimentally accessible way to detect and measure the spin gap in the bosonic ladder all the way from the weakly to the strongly interacting regime.

Discussion
We conclude this Section with a discussion of the choice of the perturbation used to probe the system. Modulating the spin imbalanceN s is a natural choice to probe the properties of the system in the spin sector from an experimental perspective (see Section 5). As pointed out in Section 4.1, both the M and V phases are gapless, and thus the choice of the modulation is crucial to distinguish them, since a generic one will in principle be sensitive to the presence of the gapless excitations and thus lead to absorption in both cases.
To display a counter-example, we show an additional calculation where the perturbation leads to energy absorption irrespective of whether the system is in the M or V phase. Instead of the spin density [Equation (6)], we perturb the system using the perturbationV = F (t)Ĵ s , where the longitudinal spin-current operatorĴ s is defined asĴ In Equation (9),Ĵ j,m is the current operator on the link between site j and j + 1, and on the leg m:Ĵ The results of the simulation are shown in Figure 5. We use the same system parameters as in Figure 3(c). The blue points correspond to the data forε(ω)/ω, as a function of J ⊥ , when the system is modulated by usingN s , whereas the red point correspond tȯ ε(ω)/ω whenĴ s is instead used. As we show in the figure, when we perturb the system usingĴ s , energy absorption takes place both in the V and in the M phase. Thus, the choice of using the spin current as a perturbation does not allow us to probe the spin gap, differently from the case when the spin density is used.

Experimental realization using laser-induced tunneling
As discussed in the Introduction, most experimental realizations of bosonic flux ladders with cold atoms do not strictly realize the situation described by the Hamiltonian in Equation (1) due to different interaction terms. In the approach of [9], the interaction energy per atom is very weak due to the large number of atoms per site, and in [39], interactions are long-ranged in the synthetic (spin) dimension. The bosonic FL with strong, short-range interactions, but only for two particles, has been also investigated [38]. Here, we discuss an alternative experimental realization that follows from the proposal of [64] for realizing the Harper-Hofstadter Hamiltonian in a square optical lattice. This scheme naturally realizes a bosonic FL with short-range (on-site) interactions and low filling around or below one atom per site. We first review the scheme described in [64]. We consider an atomic species with two long-lived internal states connected by an ultra-narrow optical transition as used in optical atomic clocks [72]. This can be realized, e.g. using the singlet 1 S 0 = g GS and a metastable 3 P 0 = e state in group-II or Ytterbium atoms. The atoms are trapped in two dimensions by a strong confining potential along z, and in the x − y plane by a state-dependent square optical lattice trapping atoms in different sublattices depending on their internal state (see Figure 6 and [64,73]). The y lattice of period d y is chosen to trap atoms in both internal states identically. The x potential is formed by the sum of a short lattice with spacing d x , V x,µ (x) = µ V 0,x cos 2 (πx/d x + φ SL ), with µ = g, e and with µ = +1 for g and −1 for e, and of a long lattice with spacing 2d x , W µ (x) = W µ cos 2 (πx/2d x + φ W ), with a well-controlled relative phase φ W [74]. By suitably choosing the depths of the x lattices, one can suppress standard tunneling along x within each sublattice g or e.
A laser of wavevector k L is then used to coherently couple states g and e, thereby inducing hopping between the g and e sublattices. This laser-assisted tunneling process [75,76] is described by a tight-binding Hamiltonian of the form (1) with with d y = d y e y . For Ytterbium atoms, for instance, d y 380 nm and 2π/k L 578 nm, leading to a maximum value of Φ max 0.66 when the coupling laser propagates along y. The value of Φ can be tuned between 0 and Φ max by changing the direction of propagation of the laser. A calculation of the band structure leads to laser-induced tunneling energies of J ⊥ /h ∼ 100 Hz for V 0,x = 8 E r,x and W µ V 0,x , where E r,x /h 3 kHz is the recoil energy associated with the period-d x lattice [64]. Note that J ⊥ is proportional to the power of the coupling laser, and that the intra-leg tunneling J is tuneable independently by changing the depth of the y lattice.
The simultaneous presence of the superlattice and laser coupling enlarges the unit cell to 2d x , with in general four non-equivalent sites per unit cell (two associated with g and two with e). This corresponds to four different types of g − e "links" and to four different transition frequencies, which are non-degenerate for a generic φ W . By a suitable choice of φ W , two of these links can be made degenerate [64]. Connecting all neighboring lattice sites with resonant laser-assisted tunneling then requires three different transition frequencies ω 1 , ω 1 ± W/ (where W is related to the amplitudes W e , W g ). Choosing W large enough compared to the laser-induced tunneling energies J ⊥ ensures that a given laser frequency only enables tunneling for the links where it is resonant (typically one can choose W/h ∼ 8 kHz and W/J ⊥ ∼ 80). This setup leads to a two-dimensional Hofstadter optical lattice with a uniform flux Φ through each unit cell. This fully connected Hofstadter lattice can be reduced in a straightforward manner to an array of two-leg ladders by removing every other frequency ω 1 ± W/ (see Figure 6(b)). Similarly, three-leg ladders could be realized by removing only one frequency, for instance ω 1 + W/ .
Focusing now on the two-leg ladder geometry, each leg of the ladders is associated with a different internal state g ≡ +1/2 or e ≡ −1/2. In this situation, time-offlight and state-dependent imaging (see, e.g., [77]) gives access to the leg-resolved MDF. Furthermore, a non-zero detuning δ 1 = ω 1 − ω eg of the coupling laser from the atomic resonance ω eg generates a term ∝N s , as desired for the spectroscopy protocol presented in Section 4. Frequency modulation of ω 1 is straightforward to implement using acoustoor electro-optical modulators, and energy absorption can be detected by monitoring the changes of the MDF.

Conclusions
In this article, we have investigated the properties of bosonic flux ladders from the dilute to the strongly correlated regime. For particle densities n < 1, the phase transition from a Meissner to a vortex phase is qualitatively unchanged, but quantitatively strongly affected by interactions. With the help of numerical simulations, we have shown that this phase transition can be observed by recording the momentum distribution, and that its precise location is well identified by the "imbalance ratio" characterizing the multior single-peak character of the momentum distribution.
Moreover, we have discussed a spectroscopic method that employs a periodic modulation of the spin imbalance between the two legs as a probe of the excitation spectrum. Gapped spin-like excitations in the Meissner phase prevent energy absorption below a certain frequency threshold, that we identified with the spin gap; in contrast, energy absorption occurs at all frequencies in the vortex phase. As such, monitoring the energy absorbed versus the modulation frequency allows one to measure not only the location of the phase transition, but also the value of the spin gap.
The characterization of the low-energy properties of a quantum many-body system is as important as the characterization of the state itself. Since we have shown that the protocols discussed in this article are within the reach of state-of-the art experiments, we believe that our work will motivate further interest in the study of the low-energy properties of complex quantum phases by indicating an effective procedure to be applied in the non-trivial cases where gapped and gapless excitations of different nature coexist.
where s 1 is a non-universal value and c is the central charge, which gives the number of gapless modes in the system. Thus, for n < 1, one predicts c = 2 in the V phase, and c = 1 in the M phase, where the spin sector is gapped [21]. The analysis of the EE and of the central charge for the data of the phase diagram in Figure 2(b), with n = 1/2, is reported in Figure A1. In panels (a) and (b), we show the EE for different values of U/(2J) as in the legends, across the V-M phase transition (see Figure 2(b), magenta line). We perform a fit with Equation (A.1) (black dashed lines in Figure 2(b)) to extract the central charge. Close to the V-M phase transition, Equation (A.1) fails to describe the behaviour of the EE, but sufficiently far away from the transition point the fit agrees well with the numerical data. Such behaviour of the EE has been observed also in other models [34,55,[78][79][80], and ascribed to the fact that, in the vicinity of the phase transition, the low-energy excitations become massive because of the presence of a gapped low-energy spectrum, and the leading order of S( ) is not described by Equation (A.1) any more.
We We now discuss the numerical estimation of the critical point for hard-core bosons at n = 1/4 and Φ/π = 0.8. We have chosen U → ∞ to simulate longer chains (L = 96) and reduce finite-size effects while having a sufficiently low numerical complexity. We compute the total MDF in the condensed-matter gauge and the EE, from which we extract the central charge. The result is shown in Figure A2. Deep in the V phase, the MDF displays two symmetric peaks with respect to k = 0. As the V-M phase transition is approached, additional peaks around k = 0 start to appear, and one peak eventually dominates when one enters the M phase. The phase transition is also signaled by the jump of the central charge, which drops from c = 2 in the V phase to c = 1 in the M phase. The EE and central charge display the same behaviour as in the previous case, and are analyzed in the same way. We finally estimate J ⊥,c (∞, n) 3.8 J from the behaviour of the central charge, which agrees with the value we estimate by measuring the IR (Figure 2(a)).

Appendix B. Details on the time-dependent numerical calculations
In this appendix, we discuss the effect of a finite value of the bond link D max,t and of the time step dt in the numerical calculation using the TEBD algorithm. In order to ensure the reliability of our data for long times, the value of D max,t must be large enough to take into account the increasing amount of entanglement in the system, which is particularly important for the deep V phase. We first focus on the data in panels (a) and (b), which are taken using ω = 10 −2 J/ , with D max,t = 300 (for J ⊥ /J = 1.0, 2.0) and D max,t = 200 (for J ⊥ /J ≥ 3.0). In our simulations, we see that the bond link D t starts to saturate to D max,t at the sites around L/2 after a time which is smaller than the total simulation time.
In order to see how this fact affects our data of ∆E(t), we compare the results for ∆E(t) by using D max,t = 200 and D max,t = 300. The result is shown in Figure B1. In particular, we separately show ∆E(t) in the V phase ( Figure B1(a)) and in the M phase ( Figure B1(b)). The data at D max,t = 200 and D max,t = 300 are shown using solid and dashed lines respectively. As evident from the figure, the curves with D max,t = 200 become significantly different in the deep V phase (J ⊥ /J = 1.0, 2.0) from the curves computed using D max,t = 300 for times which are between t = 500 /J and t = 750 /J, i.e., after D t has saturated to D max,t = 200 almost on all sites of the chain. Indeed, in the V phase, where we have c = 2, we see that the saturation of the bond link to D max,t = 200 starts after t 160 /J, for J ⊥ /J = 1.0, 2.0, and after t = 370 /J for J ⊥ /J = 3.0. Instead, in the M phase, where we have c = 1, the bond link increases in time with a smaller rate with respect to the data in the V phase: for the data at J ⊥ /J = 4.0, we start to see saturation of the bond link to D max,t = 200 after t 1500 /J, whereas D t never saturates for J ⊥ /J > 4.0. Thus, from this analysis, we see that we need to use at least D max,t = 300 for the data at J ⊥ /J = 1.0, 2.0, whereas we can use D max,t = 200 for the others.
We perform the same analysis for the data in Figure 4(c), which are taken at ω = 10 −1 J/ . In this case, we can simulate up to shorter times with respect to the case in Figure 4(a). This allows us to use larger values of the bond link, which we choose D max,t = 350. In order to see the effect of the finite value of D max,t , we then compare these data of ∆E(t) with the data computed using D max,t = 300. The result is shown in Figure B2. In this case, we observe that D t starts to saturate to D max,t = 350 already at t = 25 /J in the V phase. The fact the D t grows in time with a larger rate with respect to the case in Figure B1 is due to the larger value of ω that we use.
As in Figure B1, the data in the M phase ( Figure B2(b)) are less sensitive to the bond link difference with respect to the data in the V phase because of the smaller amount of entanglement. As we found for the data in Figure B1, we here see that the different values of D max,t during the TEBD algorithm do not drastically affect the qualitative behaviour of ∆E(t), and thus of the EAR, at least for the times considered for the fits.
As we pointed out in Section 4.2, also the time step dt has to be properly chosen in order to ensure the correct convergence of the TEBD algorithm. In our algorithm, during the time evolution and for each Trotter step, the Hamiltonian is taken constant within the time interval dt. Since the Hamiltonian depends explicitly on time through the function F (t) [see Equation (6)], it is important to check the validity of this approximation for the choice dt = 10 −2 /J, specifically in the large-ω limit considered in the data in Figure 3(f ). To do so, we compared the results of the simulations in Figure 3(f ) with the results of a simulation with the same parameters but using dt = 10 −3 /J. We found that, |ε(ω, dt = 10 −2 /J)/ω −ε(ω, dt = 10 −3 /J)/ω|/J 10 −4 even for the largest values of ω that we consider (not shown), whereε(ω, dt)/ω indicates the data series of the EAR per unit frequency taken using the time step dt. Therefore, the choice of dt = 10 −2 /J in the fourth-order Trotter expansion is sufficient to ensure the correct convergence of the TEBD algorithm.

Appendix C. Additional data in the interacting regime
We now extend the discussion carried out in Figure 3(e) and Figure 4(d). In the former case (L = 24 and N = 2), we showed that the energy absorption starts to be suppressed when J ⊥ /J 6.0, in agreement with the analytical result for free bosons J ⊥,c (U = 0) 5.9 J [18], whereas in the latter case (L = 24 and N = 6), the suppression of the energy absorption was observed approximatively from J ⊥ /J 4.0. This result is in agreement with the critical point J ⊥,c 3.8 J that we numerically estimated in Figure A2, in the case of a long chain (L = 96) at filling n = 1/4. In order to pinpoint the reliability of these numerical data and in order to ensure that the shift of the critical point, estimated from the spectroscopic method, that we observe from  is not an artifact due to finite-size effect, we here show additional numerical data, simulating hard-core bosons with the same parameters as in Figure 3(e) and Figure 4(d), but using a different value of the density, n = 1/6 (i.e., N = 4 with L = 24). The data of the EAR per unit frequency as a function of J ⊥ /J are shown in Figure C1(a): the energy absorption is nonzero for J ⊥ /J 4.0, and it starts to be suppressed between J ⊥ 4.0 J and J ⊥ 5.0 J, suggesting that the spin gap opens between these two value of J ⊥ /J. As we did for the n = 1/4, we compare this result with the bahaviour of the central charge ( Figure C1(b)), computed simulating hard-core bosons at n = 1/6 and L = 96. As evident from the figure, the central charge drops from values which are close to c = 2 (vortex phase) to values close to c = 1 (Meissner phase) between J ⊥ = 4.5 J and J ⊥ = 5.0 J, in agreement with the value estimated by looking at the EAR per unit frequency in panel (a). In the light of these results together with the results discussed in the previous appendices, we are confident about the reliability of the computed energy change ∆E(t) and EARε(ω).