# Spin-gap spectroscopy in a bosonic flux ladder

Marcello Calvanese Strinati, Fabrice Gerbier, Leonardo Mazza

► To cite this version:

Marcello Calvanese Strinati, Fabrice Gerbier, Leonardo Mazza. Spin-gap spectroscopy in a bosonic flux ladder. New Journal of Physics, 2018, 20, pp.015004. 10.1088/1367-2630/aa9ca2. hal-01696270

HAL Id: hal-01696270

https://hal.sorbonne-universite.fr/hal-01696270

Submitted on 30 Jan 2018

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
PAPER • OPEN ACCESS

# Spin-gap spectroscopy in a bosonic flux ladder

To cite this article: Marcello Calvanese Strinati et al 2018 *New J. Phys.* **20** 015004

View the [article online](#) for updates and enhancements.
**New Journal of Physics**

The open access journal at the forefront of physics

Deutsche Physikalische Gesellschaft DPG

IOP Institute of Physics

Published in partnership
with: Deutsche Physikalische
Gesellschaft and the Institute
of Physics

<figure><img src="image_2.png" /></figure>

CrossMark

OPEN ACCESS

RECEIVED

9 August 2017

REVISED

21 October 2017

ACCEPTED FOR PUBLICATION

23 November 2017

PUBLISHED

8 January 2018

Original content from this
work may be used under
the terms of the *Creative
Commons Attribution 3.0
licence*.

Any further distribution of
this work must maintain
attribution to the
author(s) and the title of
the work, journal citation
and DOI.

<figure><img src="image_3.png" /></figure>

PAPER

Spin-gap spectroscopy in a bosonic flux ladder

Marcello Calvanese Strinati<fnref n="1" />,<fnref n="2" />,<fnref n="4" />, Fabrice Gerbier<fnref n="3" /> and Leonardo Mazza<fnref n="2" />

E-mail: marcello.calvanesestrinati@gmail.com

**Keywords:** flux ladders, synthetic gauge fields, Meissner and vortex phases, matrix-product states

Abstract

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-state simulations, we investigate the properties of the so-called ‘Meissner’ and ‘vortex’ phases which appear in such a system, focusing on the experimentally accessible observables. We discuss how to experimentally monitor the phase transition, and show that the response to the 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 us to quantitatively measure the spin gap in the Meissner phase. We finally discuss its experimental implementation.

**1. Introduction**

Orbital magnetism (OM) encompasses a host of phenomena that arise in the 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 has been a trademark of quantum mechanics since its early days. In the case of electrons in solids, for instance, the OM effects include Landau diamagnetism [3] and the integer and fractional quantum Hall effects [4, 5].

Flux ladders (FLs) composed of two (or more) coupled one-dimensional subparts with a magnetic field perpendicular to the ladder plane are among the simplest setups in which OM can appear. FLs are quasi-one-dimensional, and thus still amenable to 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 FLs is one of the major motivations of this research field.

Bosonic two-leg FLs have been particularly well studied, in part due to the simplicity of the model, and in part because of the recent experimental realization of 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. The V phase is characterized by a nonvanishing inter-leg (‘transverse’) current, and the M phase is characterized by the 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 field-theory treatments of the low-energy part of the model, two-leg ladders generally feature 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 four situations obtained by combining these two classifications—V-SF, M-SF, V-MI and M-MI—are possible. Numerical studies of the microscopic models of interacting bosonic FLs have confirmed the existence of these four phases and more, revealing an extraordinarily rich phenomenology [12–33]. For instance, it has been proposed recently that the

<footnote n="1">¹ NEST, Scuola Normale Superiore & Istituto Nanoscienze-CNR, I-56126 Pisa, Italy</footnote>

<footnote n="2">² Département de Physique, Ecole Normale Supérieure / PSL Research University, CNRS, 24 rue Lhomond, F-75005 Paris, France</footnote>

<footnote n="3">³ Laboratoire Kastler Brosser, Collège de France, CNRS, ENS-PSL Research University, UPMC-Sorbonne Universités, 11 place Marcelin Berthelot, F-75005 Paris, France</footnote>

<footnote n="4">⁴ Author to whom any correspondence should be addressed.</footnote>
precursors of the physics of the fractional quantum Hall effect—and in particular of the Laughlin wave functions—might appear in experimentally relevant bosonic FLs [19, 24, 34, 35].

Experimentally, the realization of a bosonic FL belongs to a more general effort to realize effective gauge potentials coupled to ultracold atoms in spite of their electrical neutrality [36, 37]. The experiment in [9] creates a one-dimensional array of isolated ladders with a total flux per plaquette of $\Phi = \pi/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 is very weak compared to the inter- and intra-leg tunneling energies. Recently, the role of interactions in bosonic FLs 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, the interactions are short-ranging in real space, but have an almost infinite range along the synthetic dimension. This makes the 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 to the bosonic case [43–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 V-M competition in the bosonic case [47–49]. Triggered by 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–55].

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 noninteracting limit and in the dense interacting case. This extends previous work studying dynamical protocols to probe bosonic or fermionic systems in one dimension [58–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 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.

## 2. 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 the interactions [10]:

$$ \hat{H}_0^{(\text{ex})} = -J \sum_{j=1}^{L-1} \sum_{m=\pm 1/2} (\hat{b}_{j,m}^{\dagger} \hat{b}_{j+1,m} + \text{H.c.}) + J_{\perp} \sum_{j=1}^{L} (\hat{b}_{j-1/2}^{\dagger} \hat{b}_{j+1/2} e^{-i\Phi j} + \text{H.c.}) \\ + \frac{U}{2} \sum_{j=1}^{L} \sum_{m=\pm 1/2} \hat{n}_{j,m} (\hat{n}_{j,m} - 1). \qquad (1) $$

Here, $\hat{b}_{j,m}$ ($\hat{b}_{j,m}^{\dagger}$) annihilates (creates) a boson on site j and on leg m, $\hat{n}_{j,m} = \hat{b}_{j,m}^{\dagger} \hat{b}_{j,m}$ is the local density operator on leg m, J and $J_{\perp}$ denote the tunneling amplitude between the two nearest-neighbor (NN) sites in the longitudinal and transverse direction, respectively, and $\Phi$ 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 to be 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$.
<figure><img src="image_4.png" /><figcaption>Figure 1. A schematic representation of the two-leg ladder. Here, *J* is the tunneling amplitude between nearest-neighbor sites in the longitudinal direction *j*, *J*<sub>⊥</sub> is the tunneling amplitude in the transverse direction *m*, *Φ* is the gauge flux piercing each plaquette, and *U* is the on-site interaction strength, which is taken to be equal for both legs.</figcaption></figure>

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) translationally invariant, swapping the gauge flux with the longitudinal links, by using the unitary transformation $\hat{d}_{j,m} = e^{-i\Phi_{jm}} \hat{b}_{j,m}$. The transformed Hamiltonian reads

$$ \hat{H}_0^{(\text{cm})} = -J \sum_j \sum_{m=\pm 1/2} (\hat{d}_{j,m}^\dagger \hat{d}_{j+1,m} e^{i\Phi_m} + \text{H.c.}) + J_\perp \sum_j (\hat{d}_{j-\frac{1}{2}}^\dagger \hat{d}_{j+\frac{1}{2}} + \text{H.c.}) \\ + \frac{U}{2} \sum_j \sum_{m=\pm 1/2} \hat{n}_{j,m}(\hat{n}_{j,m} - 1), \qquad (2) $$

where $\hat{n}_{j,m} = \hat{b}_{j,m}^\dagger \hat{b}_{j,m} = \hat{d}_{j,m}^\dagger \hat{d}_{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 the reference energy scale.

For noninteracting bosons ($U=0$), the Hamiltonian in equation (2) can be diagonalized in momentum space by introducing the operators $\hat{d}_{k,m} = L^{-1/2} \sum_j e^{ikj} \hat{d}_{j,m}$. The two energy bands are given by $E_{\pm}(k) = -2J \cos(k) \cos(\Phi/2) \pm \sqrt{4J^2 \sin^2(k) \sin^2(\Phi/2) + J_{\perp}^2}$. The structure of the lower band $E_{-}(k)$ changes with $J_{\perp}$ or $\Phi$. When $J_{\perp}$ exceeds a critical value $J_{\perp,c} = 2J \sin(\Phi/2) \tan(\Phi/2)$, the lower energy band has one minimum at $k=0$. When $J_{\perp} < J_{\perp,c}$, the lower band features two symmetric minima at $k = \pm k_M(\Phi, J_{\perp})$. 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_{\perp}$ and/or $\Phi$, the system can undergo an M-V phase transition [18]. This transition persists for nonzero repulsive interactions, but the critical value $J_{\perp,c}$ (which depends on $U$, $n$, $\Phi$ in general) can be strongly modified by interactions [28, 31].

### 3. Momentum distribution functions and the phase diagram of interacting bosonic flux ladders

We numerically study the properties of a bosonic FL with repulsive interactions ($U>0$) using an 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_{\text{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 the 'spin gap'). As a consequence, the two phases also differ in the so-called central charge *c*, which 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 the variations of parameters $J_{\perp}$, $\Phi$, *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 a measurement is typically performed in the 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 cannot be accessed directly: only the total number of particles is a conserved quantity when $J_{\perp}$ and $\Phi$ are both nonzero. 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. In the next section, we propose a spectroscopic method that can be used to estimate the spin gap.
<figure><img src="image_5.png" /><figcaption>Figure 2. The phase diagram in the U/(2J) versus J$_{⊥}$/J plane. (a) We use n = 1/4, L = 72 and Φ = 0.8 π. The points on the critical line (yellow line) are found by imposing δn = 0. The green points identify the V phase, whereas the red points identify the M phase. (b) The same analysis as in panel (a) but using n = 1/2 and L = 24. In the latter case, we have to use smaller values of L because of the higher numerical complexity. The critical line is shifted towards smaller values of J$_{⊥}$. To increase the numerical accuracy, we compute the critical line using L = 48 (magenta line), which overlaps the one obtained using L = 24 (yellow line). (c)–(d) The data for n (cm) (k) for the phase diagram in panel (b), using L = 48, for (c) J$_{⊥}$/J = 1.50, (d) J$_{⊥}$/J = 1.75.</figcaption></figure>

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 the total. Time-of-flight measurements readily give access to the total MDF; in some experimental schemes, such as the one discussed in section 5, it is even possible to measure it for a specific leg. The leg-resolved MDF in the experimental gauge is defined as

$$n_m^{(\text{ex})}(k) = \langle \hat{b}_{k,m}^{\dagger} \hat{b}_{k,m} \rangle = \frac{1}{L} \sum_{j,h=1}^{L} e^{-ik(j-h)} \langle \hat{b}_{j,m}^{\dagger} \hat{b}_{h,m} \rangle, \quad (3)$$

where the expectation value is computed over the GS of the Hamiltonian in equation (1). Since the MDF is periodic with the period $2\pi$, we restrict the momentum variable to $k \in [-\pi : \pi)$. 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_m^{(\text{cm})}(k) = n_m^{(\text{ex})}(k - m\Phi)$. The total MDFs are accordingly $n_m^{(\text{ex})}(k) = \sum_{m=\pm 1/2} n_m^{(\text{ex})}(k)$ and $n_m^{(\text{cm})}(k) = \sum_{m=\pm 1/2} n_m^{(\text{cm})}(k)$.

In the condensed-matter gauge, the MDF displays one peak centered at $k = 0$ for the M phase, and two symmetric peaks at $k = \pm k_M$ for the V phase, reminiscent of the single or double minimum of the lower energy band when $U = 0$ [32]. In figure 2 we report 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)

$$\delta n(\Phi, J_{\perp}, U, n) := \frac{n^{(\text{cm})}(k_M) - n^{(\text{cm})}(0)}{n^{(\text{cm})}(k_M) + n^{(\text{cm})}(0)}. \quad (4)$$

The IR takes the values 0 < δn < 1 in the V phase and −1 < δn < 0 in the M phase. We propose finding 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 the 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 noninteracting 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.

# 4. 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.

## 4.1. Model and observables

We consider the Hamiltonian $\hat{H}_0^{(\text{ex})}$ in equation (1), and add a time-periodic perturbation $\hat{V} = F(t)\hat{N}_s$ proportional to the population difference between the two legs (hereafter denoted as spin imbalance),

$$ \hat{N}_s = \hat{N}_{+\frac{1}{2}} - \hat{N}_{-\frac{1}{2}}, \qquad (5) $$

where $\hat{N}_m = \sum_j \hat{n}_{j,m}$ is the particle number per leg, with $F(t) = \delta_1 \sin(\omega t) \Theta(t)$, and $\Theta(t)$ is the unit step function. Here, we denote by $\delta_1$ and $\omega$ the amplitude and frequency of the modulation, respectively. The total Hamiltonian is thus

$$ \hat{H}^{(\text{ex})}(t) = \hat{H}_{0}^{(\text{ex})} + \hat{F}(t) \hat{N}_{s}, \qquad (6) $$

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, $E(t) = \langle \Psi(t) | \hat{H}(t) | \Psi(t) \rangle$, where $|\Psi(t)\rangle$ is the time-evolved state starting from the GS of the bosonic ladder. We define the energy absorption rate (EAR) as

$$ \dot{\varepsilon}(\omega) = \lim_{T \to \infty} \frac{1}{T} \int_0^T dt \frac{\partial E}{\partial t} = \lim_{T \to \infty} \frac{1}{T} \int_0^T dt \frac{\partial F}{\partial t} \langle \Psi(t) | \hat{N}_s | \Psi(t) \rangle. \qquad (7) $$

Within linear response theory, the EAR per unit frequency probes the imaginary part of the response function, i.e. $\dot{\varepsilon}(\omega)/\omega \propto \text{Im}[\chi_{N_s-N_s}(\omega)]$, where $\chi_{N_s-N_s}(\omega) = \int_0^\infty dt e^{i\omega t} \chi_{N_s-N_s}(t)$ and where

$$ \chi_{N_s-N_s}(t) = -\frac{i}{\hbar} \langle \Psi_0 | [\hat{N}_s(t), \hat{N}_s(0)] | \Psi_0 \rangle \qquad (8) $$

is the response function in real time. Here $\hat{N}_s(t) = e^{i\hat{H}_0 t/\hbar} \hat{N}_s$ is the number imbalance expressed in the interaction picture with respect to $\hat{H}_0$. Notice that by means of equation (7), the EAR can be experimentally accessed by measuring the total spin imbalance in time $\langle \Psi(t) | \hat{N}_s | \Psi(t) \rangle$ (see section 5).

If we denote by $\Delta E_s$ the value of the spin gap, the spectroscopic method that identifies it should consist of a periodic modulation of the system which is sensitive to its presence, so that the system does not absorb energy as long as $\hbar\omega < \Delta E_s$, and energy absorption can only occur for $\hbar\omega > \Delta E_s$. We thus expect $\text{Im}[\chi_{N_s-N_s}(\omega)] = 0$ if $\hbar\omega < \Delta E_s$ and $\text{Im}[\chi_{N_s-N_s}(\omega)] > 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, $|\Psi(t)\rangle$ is computed by using the time-evolving block decimation (TEBD) algorithm [8, 68, 69] with a fourth-order Trotter expansion [70, 71] and a time step $dt$ (during the time evolution, we fix the maximum bond link $D_{\text{max},t}$ used to describe the MPS state at time $t$).

## 4.2. Results for dilute gases

We first analyze a very dilute gas (n ≪ 1) where the 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
<figure><img src="image_6.png" /><figcaption>**Figure 3.** The energy absorption for dilute bosons ($n = 1/12 \ll 1$). This observable monitors the response to a modulation of the spin imbalance with an amplitude of $\delta_1 = 0.4$ J and a frequency of $\omega$ as shown in the legends. (a) Relative energy variation, $\Delta E(t) = E(t) - E(0)$ in the V phase ($J_\perp/J = 2.0$, magenta full line) and in the M phase ($J_\perp/J = 7.0$, blue full line) for $\Phi/\pi = 0.80$. The energy absorption rate (EAR) per unit frequency and its error bars are extracted from a linear fit to $\Delta E(t)$ (black dashed line, see text for details). (b)–(e) EAR per unit frequency $\dot{\epsilon}(\omega)/\omega$ as a function of $J_\perp$ for (b) $\Phi/\pi = 0.16$, (c) $\Phi/\pi = 0.24$, (d) $\Phi/\pi = 0.64$ and (e) $\Phi/\pi = 0.80$. The insets in (b)–(f) show the noninteracting phase diagram ($U = 0$) for reference. The V–M phase transition for noninteracting bosons is indicated by the magenta arrows in the main plots. (f) The EAR per unit frequency as a function of $\omega$ for $\Phi/\pi = 0.24$ in the V (magenta data) and in the M phase (blue data). For all plots reported here, the simulations were done for hard-core bosons using $L = 24$ and $N = 2$.</figcaption></figure>

bosons ($U \to \infty$). Unlike 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 $\Phi$. The amplitude of the density modulation is $\delta_1 = 0.4$ J, and we use different values of $\omega$, depending on the value of $\Phi$. During the time evolution, the time-dependent Hamiltonian in equation (6) is taken as 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 values of $\omega$ that we consider. We have verified that we can choose the time step $dt = 10^{-2} \hbar/J$ in the Trotter expansion (see also appendix B for a deeper discussion). In figure 3(a), we show the relative energy variation, $\Delta E(t) = E(t) - E(0)$ for the two typical cases. In the M phase ($J_\perp/J = 7.0$), there is no net energy absorption for sufficiently small $\omega$, whereas the system absorbs energy for all $\omega$ in the V phase ($J_\perp/J = 2.0$). The EAR is extracted from the slope of $\Delta E(t)$ represented by the black dashed line. To remove the fast oscillations of $\Delta E(t)$ and extract the long-time linear trend, we perform M linear fits to $\Delta E(t)$ using different ranges of $t$. Accordingly,
we obtain a set of values for the EAR per unit frequency, $\{\dot{\epsilon}_q(\omega)/\omega\}_{q=1}^M$, from which we compute the mean value $\dot{\epsilon}(\omega)/\omega = M^{-1}\sum_{q=1}^M \dot{\epsilon}_q(\omega)/\omega$, and the standard deviation, $\sigma_{\dot{\epsilon}} = \sqrt{M^{-1}\sum_{q=1}^M [\dot{\epsilon}_q(\omega)/\omega - \dot{\epsilon}(\omega)/\omega]^2}$. We take the latter as a measure of the uncertainty on the determined slope.

In figures 3(b)–(e), we show the EAR per unit frequency as a function of $J_{\perp}$ for four different values of $\Phi$. In the insets, we show the noninteracting phase diagram ($U=0$), where the red solid line corresponds to the critical line $J_{\perp,c}$ ($U=0$), and the black dashed line indicates the line at which we cut the phase diagram. The behavior of the EAR per unit frequency is in agreement with the opening of the spin gap at the expected value of $J_{\perp,c}$. The system absorbs energy in the V phase ($J_{\perp} \lesssim J_{\perp,c}$), whereas the energy absorption ceases as the V-M phase transition takes place ($J_{\perp} \gtrsim J_{\perp,c}$).

When $J_{\perp} = 0$, the relation $[\hat{N}_s, \hat{H}_0] = 0$ holds, implying that $[\hat{N}_s(t), \hat{N}_s] = 0$ and $\chi_{N_s-N_c}(\omega) = 0$ from equation (8). This is consistent with the curve plotted in figure 3(e), which tends to 0 for low values of $J_{\perp}$. Similar behavior is also expected for the values of $\Phi$ used in figures 3(b)–(d), but the considered value of $J_{\perp}$ was not small enough to highlight it.

As we previously stated, in the presence of a spin gap, the system is only expected to absorb energy if $\hbar\omega > \Delta E_s$. In figure 3(f), we show the EAR per unit frequency as a function of $\omega$ for $\Phi/\pi = 0.24$, both in the V phase ($J_{\perp}/J = 0.1$) and in the M phase ($J_{\perp}/J = 0.5$). For low modulation frequencies, we observe that the system can absorb energy in the V phase for values of modulation frequency down to $\omega \sim 10^{-2} J/\hbar$. 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 $\omega$, one observes a drop in the response in both phases, as expected from the general behavior of the susceptibility $\chi(\omega)$ [6].

### 4.3. 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 \to \infty$) and a 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. Here we demonstrate that this shift is also detected by the periodic modulation of the density imbalance. For concreteness, we focus on $n = 1/4$ and $\Phi = 0.8\pi$, as for the data in figure 2(a).

The numerical results are shown in figure 4. In figure 4(a), we plot $\Delta E(t)$ for different values of $J_{\perp}/J$ using $\delta_1 = 0.4J$ and $\omega = 10^{-2}J/\hbar$. In figure 4(b), we display the EAR per unit frequency as a function of $J_{\perp}/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_{\perp} \ge 5.0J$, and becomes nonzero when $J_{\perp} = 4.0J$ and below. In section 2, we estimated the critical value for the V-M transition, $J_{\perp,c} \simeq 3.8$, which is slightly lower than the observed threshold for energy absorption. This quantitative discrepancy may be due both to the finite-size effects, and to the fact that $\hbar\omega$ is possibly larger than the spin gap for $J_{\perp}/J = 4.0$.

In panels (c) and (d), we show the same analysis for a larger value of $\omega$, namely $\omega = 10^{-1}J/\hbar$. For $J_{\perp}/J \lesssim 4.0$, the system absorbs energy until saturation starts to take place. Instead, for $J_{\perp}/J \gtrsim 4.0$, energy absorption is suppressed. Unlike the data in panels (a) and (b), we also see nonzero energy absorption for $J_{\perp}/J = 5.0, 6.0$. We ascribe this fact to the larger value of $\hbar\omega$, possibly overcoming the value of the spin gap. The numerical complexity of the problem prevents us from using lower values of $\omega$, as the required simulation times $t$ are beyond our numerical possibilities. For a more critical discussion of the numerical data, see appendix B.

To conclude, our results are compatible with the opening of a spin gap around $J_{\perp}/J \simeq 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 of detecting and measuring the spin gap in the bosonic ladder all the way from the weakly to the strongly interacting regime.

### 4.4. Discussion

We conclude this section with a discussion of the choice of the perturbation used to probe the system. Modulating the spin imbalance $\hat{N}_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 for distinguishing them, since a generic one will, in principle, be sensitive to the presence of 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 perturbation $\hat{V} = F(t)\hat{J}_s$, where the longitudinal spin-current operator $\hat{J}_s$ is defined as
<figure><img src="image_7.png" /><figcaption>Figure 4. Energy absorption for hard-core bosons with $n = 1/4$. (a) Energy absorption $\Delta E(t)$, and (b) the EAR per unit frequency for a modulation frequency $\omega = 10^{-2} J/\hbar$. The magenta arrow indicates the estimated position of the M-V phase transition for hard-core bosons, $J_{\perp,c}/J \simeq 3.8$ (see appendix B). (c), (d) The energy absorption and EAR per unit frequency for $\omega = 10^{-1} J/\hbar$. In this case, we see a less sharp decrease in the energy absorption for $J_{\perp}/J \gtrsim 4.0$. For all plots, the spin imbalance is modulated with an amplitude of $\delta_1 = 0.4 J$ and we use $L = 24$ and $N = 6$. The data at $J_{\perp}/J = 1.0, 2.0$ in panels (a) and (b) is obtained using a bond length $D_{max,t} = 300$, whereas the other data is taken using $D_{max,t} = 200$. For panels (c) and (d), we use $D_{max,t} = 350$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 5. The EAR per unit frequency $\dot{\epsilon}(\omega)/\omega$ as a function of $J_{\perp}$ for a dilute gas ($n = 1/12$) at $\Phi = 0.24 \pi$; other parameters are as in figure 3(c). The blue or red points are obtained by modulating the spin imbalance $\hat{N}_s$ or the spin current $\hat{J}_s$, respectively. The magenta vertical dashed line marks the M-V transition point $J_{\perp,c}$ for noninteracting bosons. The spin current is insensitive to the opening of the spin gap, and energy absorption takes place with roughly equal strength in both M and V phases.</figcaption></figure>

$$ \hat{\mathcal{J}}_s \equiv \sum_j \hat{\mathcal{J}}_{s,j} = \sum_j (\hat{\mathcal{J}}_{j,\frac{1}{2}} - \hat{\mathcal{J}}_{j,-\frac{1}{2}}). \qquad (9) $$

In equation (9), $\hat{\mathcal{J}}_{j,m}$ is the current operator on the link between site j and j + 1, and on leg m:

$$ \hat{\mathcal{J}}_{j,m} = -iJ(\hat{b}_{j,m}^{\dagger}\hat{b}_{j+1,m} - \text{H.c.}). \qquad (10) $$

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 $\dot{\epsilon}(\omega)/\omega$, as a function of $J_{\perp}$, when the system is modulated by using $\hat{N}_s$,
<figure><img src="image_9.png" /><figcaption>Figure 6. A possible experimental realization following [64]. (a) A sketch of the laser arrangement; (b) a spin-dependent lattice realizing a collection of disconnected two-leg ladders.</figcaption></figure>

whereas the red point corresponds to $\dot{\epsilon}(\omega)/\omega$ when $\hat{J}_s$ is used instead. As we show in the figure, when we perturb the system using $\hat{J}_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, unlike the case in which the spin density is used.

## 5. 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-range in the synthetic (spin) dimension. The bosonic FL with strong, short-range interactions, but for only two particles, has also been 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 $^1S_0 = g$ GS and a metastable $^3P_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 that traps 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,\mu}(x) = \epsilon_\mu V_{0,x} \cos^2(\pi x/d_x + \phi_{SL})$, with $\mu = g, e$ and with $\epsilon_\mu = +1$ for $g$ and $-1$ for $e$, and of a long lattice with spacing $2d_x$, $W_\mu(x) = W_\mu \cos^2(\pi x/2d_x + \phi_W)$, with a well-controlled relative phase $\phi_W$ [74]. By suitably choosing the depths of the $x$ lattices, one can suppress the 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

$$ \Phi = \frac{k_L \cdot d_y}{2\pi}, \qquad (11) $$

where $d_y = d_y e_y$. For Ytterbium atoms, for instance, $d_y \approx 380$ nm and $2\pi/k_L \approx 578$ nm, leading to a maximum value of $\Phi_{max} \approx 0.66$ when the coupling laser propagates along $y$. The value of $\Phi$ can be tuned between 0 and $\Phi_{max}$ by changing the propagation direction of the laser. A calculation of the band structure leads to laser-induced tunneling energies of $J_\perp/h \sim 100$ Hz for $V_{0,x} = 8 E_{r,x}$ and $W_\mu \ll V_{0,x}$, where $E_{r,x}/h \approx 3$ kHz is the recoil energy associated with the period-$d_x$ lattice [64]. Note that $J_\perp$ is proportional to the power of the coupling laser, and that the intra-leg tunneling $J$ is tunable 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$, in general, with four nonequivalent 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 nondegenerate for a generic
$\phi_W$. By a suitable choice of $\phi_W$, two of these links can be made degenerate [64]. Connecting all neighboring lattice sites to resonant laser-assisted tunneling then requires three different transition frequencies $\omega_1$, $\omega_1 \pm W/\hbar$ (where W is related to the amplitudes $W_c$, $W_g$). Choosing a W that is large enough compared to the laser-induced tunneling energies $J_\perp$ ensures that a given laser frequency only enables tunneling for the links where it is resonant (typically one can choose $W/\hbar \sim 8$ kHz and $W/J_\perp \sim 80$). This setup leads to a two-dimensional Hofstadter optical lattice with a uniform flux $\Phi$ 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 $\omega_1 \pm W/\hbar$ (see figure 6(b)). Similarly, three-leg ladders can be realized by removing only one frequency, for instance $\omega_1 + W/\hbar$.

Focusing now on the two-leg ladder geometry, each leg of the ladders is associated with a different internal state $g \equiv +1/2$ or $e \equiv -1/2$. In this situation, time-of-flight and state-dependent imaging (see e.g. [77]) gives access to the leg-resolved MDF. Furthermore, nonzero detuning $\delta_1 = \omega_1 - \omega_{eg}$ of the coupling laser from the atomic resonance $\omega_{eg}$ generates a term $\propto \hat{N}_s$, as desired for the spectroscopy protocol presented in section 4. Frequency modulation of $\omega_1$ is straightforward to implement using acousto- or electro-optical modulators, and energy absorption can be detected by monitoring the changes of the MDF.

## 6. 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 multi- or 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, which 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 nontrivial cases where gapped and gapless excitations of a different nature coexist.

## Acknowledgments

We thank J Beugnon, M Bosch, Aguilera, R Bouganne, S De Palo and R Fazio for fruitful discussions. We are also grateful to D Rossini for providing the MPS code and for support. L M was supported by LabEX ENS-ICFP: ANR-10-LABX-0010/ANR-10-IDEX-0001-02 PSL*. This work was granted access to the HPC resources of MesoPSL financed by the Region Ile de France and the project Equip@Meso (reference ANR-10-EQPX-29-01) of the program Investissements d'Avenir supervised by the Agence Nationale pour la Recherche. We also acknowledge the CINECA award under the ISCRA initiative, for the availability of high performance computing resources and support.

## Appendix A. Numerical analysis of the phase diagram

In this appendix, we discuss our results on the phase diagram obtained in figure 2. For sufficiently large $L$, the phase transition from the V to the M phase can be numerically detected by computing the central charge, which is extracted from the entanglement entropy (EE). The EE is defined by $S(\ell) = -\text{Tr}[\hat{\rho}_\ell \log(\hat{\rho}_\ell)]$, $\hat{\rho}_\ell$ being the reduced density matrix of a bipartition of a chain of length $\ell$. In the case of OBC, the leading behavior of the EE computed on the GS is predicted to be [65]

$$S(\ell) = s_1 + \frac{c}{6} \log \left[ \left( \frac{2L}{\pi} \right) \sin \left( \frac{\pi \ell}{L} \right) \right], \quad (A.1)$$
<figure><img src="image_10.png" /></figure>

**Figure A1.** Analysis of the EE for the phase diagram in figure 2 for *n* = 1/2. We show (a) the EE for *J*<sub>⊥</sub>/J = 1.75 and (b) for *J*<sub>⊥</sub>/J = 2.00, for different values of *U*/*J* across the V-M phase transition. The phase transition is detected from the sudden change of the EE. Such a change is well reflected by the central charge, which is shown as a function of *U*/*(2J)* for (c) *J*<sub>⊥</sub>/J = 1.75 and (d) *J*<sub>⊥</sub>/J = 2.00. The uncertainties are estimated as explained in the text. Sufficiently far away from the transition point, the values of *c* that we fit are in agreement with the expected ones. The transition points estimated from the jump of the central charge are in agreement with what we found in figure 2(b) by looking at the IR (magenta dotted line).

where *s*₁ is a nonuniversal 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].

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*/(2*J*) 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 behavior of the EE, but sufficiently far away from the transition point the fit agrees well with the numerical data. Such EE behavior has also been observed in other models [34, 55, 78–80], and is 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 show the central charge as a function of $U/(2J)$, for the same set of data, in panels (c) and (d). We ascribe the fact that we do not fit $c=1$ or $c=2$ exactly to finite-size effects. Because of the oscillatory behavior of the EE and of the choice of OBC, to fit the EE and compute the values of $c$ for each value of $U/(2J)$, we repeat the fit $N_c$ times, introducing a cut-off $L_c$ which we vary from $L_c=1$ to $L_c=N_c$. For each repetition of the fit, we only include points in the range $\ell \in [L_c : L - L_c]$. We accordingly obtain a set of values for the central charge, $\{c_{L_c}^{N} \}_{L_c=1}^{N_c}$, from which we estimate the mean value as $\bar{c} = N_c^{-1} \sum_{L_c=1}^{N_c} c_{L_c}$, and the uncertainty by means of the standard deviation $\sigma_c = \sqrt{N_c^{-1} \sum_{L_c=1}^{N_c} (c_{L_c} - \bar{c})^2}$. As we see from the figure, the phase transition from the V phase to the M phase is identified by the jump of the central charge. Furthermore, sufficiently far away from the V-M transition point, the fitted values of $c$ are in agreement with the expected values predicted by bosonization. We conclude by noting that the transition points estimated from the IR in figure 2 are in agreement with the one estimated from the numerically determined central charge, the latter being known to signal the M-V phase transition [21].

The numerical simulations at finite *U* are performed by truncating the local Hilbert space on each site *j* and leg *m*, which we denote by *H*<sub>*j*,m*</sub>. We define by $|r\rangle_{j,m}$ the local Fock space such that *H*<sub>*j*,m*</sub> = span $\{|r\rangle_{j,m}\}_{r=0}^{d_{\text{loc}}-1}$, where $d_{\text{loc}} = \text{dim}(*H*<sub>*j*,m*</sub>). Let $\hat{\mathbb{P}}_r^{(j,m)} = |r\rangle_{j,m}\langle r|$ be the local projector over the state $|r\rangle_{j,m}$. The local density operator is then $\hat{n}_{j,m} = \sum_{r=0}^{d_{\text{loc}}-1} r \hat{\mathbb{P}}_r^{(j,m)}$. The suitable choice for $d_{\text{loc}}$ depends on the values of $U/(2J)$; we can
<figure><img src="image_11.png" /><figcaption>**Figure A2.** Analysis of the V-M phase transition on the $U/J \to \infty$ line of the phase diagram in figure 2(a). We simulate HCBs using $n = 1/4$ and $L = 96$. We show (a) the total MDF in the condensed-matter gauge for different values of $J_\perp/J$ across the V-M phase transition, and (b) the central charge as a function of $J_\perp/J$, computed as explained in the text. We see that the total MDF becomes single-peaked around $J_\perp/J \simeq 3.8$ (see also figure 2(a)). Accordingly, the central charge drops down to $c = 1$ around the same value of $J_\perp/J$, where the magenta dotted line indicates the value of $J_{\perp,c}/J$ found by looking at the IR. As in figure A1, the critical value agrees with the one estimated via the IR. Also, the behavior of the EE fails to be described by equation (A.1) (not shown) close to the transition point, and thus the values of $c$ that we fit deviate from the expected ones. Instead, sufficiently far away from the phase transition, the values of $c$ that we fit are in agreement with the expected values ($c = 2$ in the V phase and $c = 1$ in the M phase).</figcaption></figure>

keep up to $d_{\text{loc}}$ states for $H_{j,m}$ if the probability of finding $d_{\text{loc}} - 1$ particles on the site $j$ of leg $m$ is small with respect to the local density, i.e. $\langle \hat{\mathbf{P}}_{d_{\text{loc}}-1}^{(j,m)} \rangle \ll \langle \hat{n}_{j,m} \rangle$, for all $j$ and $m$, where the expectation value is computed on the GS of the system. Thus, we choose to verify that $L^{-1}\sum_j \langle \hat{\mathbf{P}}_{d_{\text{loc}}-1}^{(j,m)} \rangle \ll L^{-1}\sum_j \langle \hat{n}_{j,m} \rangle$ for all $m$. In our numerical simulations, we see that this condition is fulfilled using $d_{\text{loc}} = 3$, for all the values of $U/(2J)$ that we use, since we verify that $L^{-1}\sum_j \langle \hat{\mathbf{P}}_{d_{\text{loc}}-1}^{(j,m)} \rangle \lesssim \epsilon L^{-1}\sum_j \langle \hat{n}_{j,m} \rangle$, for all $m$, where $\epsilon = 3 \times 10^{-2}$ is a small numerical factor.

We now discuss the numerical estimation of the critical point for hard-core bosons at $n = 1/4$ and $\Phi/\pi = 0.8$. We have chosen $U \to \infty$ 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 behavior as in the previous case, and are analyzed in the same way. We finally estimate $J_{\perp,c}(\infty, n) \simeq 3.8 J$ from the behavior of the central charge, which agrees with the value we estimate by measuring the IR (figure 2(a)).

## Appendix B. Details of the time-dependent numerical calculations

In this appendix, we discuss the effect of a finite value of the bond link $D_{\text{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_{\text{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 is taken using $\omega = 10^{-2} J/\hbar$, with $D_{\text{max},t} = 300$ (for $J_\perp/J = 1.0, 2.0$) and $D_{\text{max},t} = 200$ (for $J_\perp/J \ge 3.0$). In our simulations, we see that the bond link $D_t$ starts to saturate to $D_{\text{max},t}$ in 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 $\Delta E(t)$ data, we compare the results for $\Delta E(t)$ by using $D_{\text{max},t} = 200$ and $D_{\text{max},t} = 300$. The result is shown in figure B1. In particular, we separately show $\Delta E(t)$ in the V phase (figure B1(a)) and in the M phase (figure B1(b)). The data at $D_{\text{max},t} = 200$ and $D_{\text{max},t} = 300$ are shown using solid and dashed lines respectively. As is evident from the figure, the curves with $D_{\text{max},t} = 200$ become significantly different in the deep V phase ($J_\perp/J = 1.0, 2.0$) from the curves computed using $D_{\text{max},t} = 300$ for times which are between $t = 500\hbar/J$ and $t = 750\hbar/J$, i.e. after $D_t$ has saturated to $D_{\text{max},t} = 200$ on almost 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_{\text{max},t} = 200$ starts after $t \simeq 160\hbar/J$, for $J_\perp/J = 1.0, 2.0$, and after $t = 370\hbar/J$ for $J_\perp/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_\perp/J = 4.0$, we start to see the saturation of the bond link to $D_{\text{max},t} = 200$ after
<figure><img src="image_12.png" /></figure>

**Figure B1.** The numerical data of $\Delta E(t)$ for the data in figure 4(a). We show $\Delta E(t)$ for (a) $J_\perp/J = 1.0$ (red lines), $J_\perp/J = 2.0$ (blue lines) and $J_\perp/J = 3.0$ (magenta lines), and (b) $J_\perp/J = 4.0$ (green lines), $J_\perp/J = 5.0$ (grey lines) and $J_\perp/J = 6.0$ (black lines). Solid lines are taken using $D_{\text{max},t} = 200$, whereas dashed lines are taken with $D_{\text{max},t} = 300$.

<figure><img src="image_13.png" /></figure>

**Figure B2.** The numerical data of $\Delta E(t)$ for the data in figure 4(c). We show $\Delta E(t)$ for (a) $J_\perp/J = 1.0$ (red lines), $J_\perp/J = 2.0$ (blue lines) and $J_\perp/J = 3.0$ (magenta lines), and (b) $J_\perp/J = 4.0$ (green lines), $J_\perp/J = 5.0$ (grey lines) and $J_\perp/J = 6.0$ (black lines). Solid lines are taken using $D_{\text{max},t} = 350$, whereas dashed lines are taken with $D_{\text{max},t} = 300$.

$t \approx 1500 \, \hbar/J$, whereas $D_t$ never saturates for $J_\perp/J > 4.0$. Thus, from this analysis, we see that we need to use at least $D_{\text{max},t} = 300$ for the data at $J_\perp/J = 1.0, 2.0$, whereas we can use $D_{\text{max},t} = 200$ for the others.

We perform the same analysis for the data in figure 4(c), which is taken at $\omega = 10^{-1} J/\hbar$. In this case, we can simulate 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 as $D_{\text{max},t} = 350$. In order to see the effect of the finite value of $D_{\text{max},t}$, we then compare the data of $\Delta E(t)$ with the data computed using $D_{\text{max},t} = 300$. The result is shown in figure B2. In this case, we observe that $D_t$ starts to saturate to $D_{\text{max},t} = 350$ already at $t = 25 \, \hbar/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 $\omega$ that we use.

As in figure B1, the data in the M phase (figure B2(b)) is 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 see here that the different values of $D_{\text{max},t}$ during the TEBD algorithm do not drastically affect the qualitative behavior of $\Delta E(t)$, and thus of the EAR, at least for the times considered for the fits.

As we pointed out in section 4.2, the time step $dt$ also 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 to be 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} \, \hbar/J$, specifically in the large-$\omega$ 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} \, \hbar/J$. We found that $|\dot{\varepsilon}(\omega, dt = 10^{-2} \, \hbar/J)/\omega - \dot{\varepsilon}(\omega, dt = 10^{-3} \, \hbar/J)/\omega| \lesssim 10^{-4}$ even for the largest values of $\omega$ that we consider (not shown), where $\dot{\varepsilon}(\omega, dt)/\omega$ indicates the data series of the EAR per unit frequency taken using the time step $dt$. Therefore, the choice of $dt = 10^{-2} \, \hbar/J$ in the fourth-order Trotter expansion is sufficient to ensure the correct convergence of the TEBD algorithm.
<figure><img src="image_14.png" /></figure>

**Figure C1.** (a) EAR per unit frequency for $n = 1/6$ and $L = 24$. Other simulation parameters as in figures 3(e) and 4(d). Here, the energy absorption starts to be suppressed between $J_⊥ = 4.0$ J and $J_⊥ = 5.0$ J, i.e. at an intermediate value between the critical point at $n = 1/12$ ($J_{⊥,c} \simeq 5.9$ J) and the one at $n = 1/4$ ($J_{⊥,c} \simeq 3.8$ J). (b) The central charge $c$ for hard-core bosons at $n = 1/6$ as in panel (a), but using $L = 96$. We see that $c$ drops from $c = 2$ (vortex phase) to $c = 1$ (Meissner phase) between $J_⊥ = 4.5$ J and $J_⊥ = 5.0$ J (magenta line), in agreement with the value estimated by looking at the EAR per unit frequency in panel (a). We conclude that the behavior of the EAR per unit frequency in panel (a) also correctly signals the opening of the spin gap for this value of density.

## Appendix C. Additional data in the interacting regime

We now extend the discussion carried out in figures 3(e) and 4(d). In the former case ($L = 24$ and $N = 2$), we showed that energy absorption starts to be suppressed when $J_⊥/J \simeq 6.0$, in agreement with the analytical result for free bosons $J_{⊥,c}(U = 0) \simeq 5.9$ J [18], whereas in the latter case ($L = 24$ and $N = 6$), the suppression of energy absorption was observed approximately from $J_⊥/J \simeq 4.0$. This result is in agreement with the critical point $J_{⊥,c} \simeq 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 this numerical data and ensure that the shift of the critical point estimated by the spectroscopic method, which we observe in figure 3(e) to figure 4(d), is not an artifact due to the finite-size effect, here we show additional numerical data, simulating hard-core bosons with the same parameters as in figures 3(e) and 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$ is shown in figure C1(a): the energy absorption is nonzero for $J_⊥/J \lesssim 4.0$, and it starts to be suppressed between $J_⊥ \simeq 4.0$ J and $J_⊥ \simeq 5.0$ J, suggesting that the spin gap opens between these two values of $J_⊥/J$. As we did for $n = 1/4$, we compare this result with the behavior of the central charge (figure C1(b)), computing simulated hard-core bosons at $n = 1/6$ and $L = 96$. As is 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 light of these results, together with the results discussed in the previous appendices, we are confident about the reliability of the computed energy change $\Delta E(t)$ and EAR $\dot{\varepsilon}(\omega)$.

## References

[1] Bohr N 1972 *Niels Bohr Collected Works* vol 1 (Amsterdam: Elsevier)

[2] Van Leeuwen H J 1921 *J. Phys. Radium* **2** 361–77

[3] Landau L and Lifshitz E 2013 *Statistical Physics* vol 5 (Amsterdam: Elsevier)

[4] Klitzing K V, Dorda G and Pepper M 1980 *Phys. Rev. Lett.* **45** 494–7

[5] Tsui D C, Stormer H L and Gossard A C 1982 *Phys. Rev. Lett.* **48** 1559–62

[6] Giamarchi T 2003 *Quantum Physics in One Dimension* International Series of Monographs on Physics (Oxford: Clarendon)

[7] Schollwöck U 2005 *Rev. Mod. Phys.* **77** 259–315

[8] Schollwöck U 2011 *Ann. Phys.* **326** 96–192

[9] Atala M, Aidelsburger M, Lohse M, Barreiro J T, Paredes B and Bloch I 2014 *Nat. Phys.* **10** 588

[10] Orignac E and Giamarchi T 2001 *Phys. Rev. B* **64** 144515

[11] Donohue P and Giamarchi T 2001 *Phys. Rev. B* **63** 180508

[12] Granato E 2005 *Phys. Rev. B* **72** 104521

[13] Rizzi M, Cataudella V and Fazio R 2006 *Phys. Rev. B* **73** 100502

[14] Dhar A, Maji M, Mishra T, Pai R V, Mukerjee S and Paramekanti A 2012 *Phys. Rev. A* **85** 041602

[15] Petrescu A and Le Hur K 2013 *Phys. Rev. Lett.* **111** 150601

[16] Dhar A, Mishra T, Maji M, Pai R V, Mukerjee S and Paramekanti A 2013 *Phys. Rev. B* **87** 174501

[17] Wei R and Mueller E J 2014 *Phys. Rev. A* **89** 063617

[18] Tokuno A and Georges A 2014 *New J. Phys.* **16** 073005

[19] Petrescu A and Le Hur K 2015 *Phys. Rev. B* **91** 054520
[20] Di Dio M, Citro R, De Palo S, Orignac E and Chiofalo M L 2015 *Eur. Phys. J.* **224** 525–31

[21] Piraud M, Heidrich-Meisner F, McCulloch I P, Greschner S, Vekua T and Schollwöck U 2015 *Phys. Rev.* **B** **91** 140406

[22] Uchino S and Tokuno A 2015 *Phys. Rev.* **A** **92** 013625

[23] Di Dio M, De Palo S, Orignac E, Citro R and Chiofalo M L 2015 *Phys. Rev.* **B** **92** 060506

[24] Cornfeld E and Sela E 2015 *Phys. Rev.* **B** **92** 115446

[25] Kolley F, Piraud M, McCulloch I P, Schollwöck U and Heidrich-Meisner F 2015 *New J. Phys.* **17** 092001

[26] Greschner S, Huerga D, Sun G, Poletti D and Santos L 2015 *Phys. Rev.* **B** **92** 115120

[27] Natu S S 2015 *Phys. Rev.* **A** **92** 053623

[28] Greschner S, Piraud M, Heidrich-Meisner F, McCulloch I P, Schollwöck U and Vekua T 2015 *Phys. Rev. Lett.* **115** 190402

[29] Uchino S 2016 *Phys. Rev.* **A** **93** 053629

[30] Bilitewski T and Cooper N R 2016 *Phys. Rev.* **A** **94** 023630

[31] Greschner S, Piraud M, Heidrich-Meisner F, McCulloch I P, Schollwöck U and Vekua T 2016 *Phys. Rev.* **A** **94** 063628

[32] Orignac E, Citro R, Di Dio M, De Palo S and Chiofalo M L 2016 *New J. Phys.* **18** 055017

[33] Anisimovas E, Račiūnas M, Sträter C, Eckardt A, Spielman I B and Juzeliūnas G 2016 *Phys. Rev.* **A** **94** 063632

[34] Calvanese Strinati M, Cornfeld E, Rossini D, Barbarino S, Dalmonte M, Fazio R, Sela E and Mazza L 2017 *Phys. Rev.* **X** **7** 021033

[35] Petrescu A, Piraud M, Roux G, McCulloch I P and Le Hur K 2017 *Phys. Rev.* **B** **96** 014524

[36] Dalibard J, Gerbier F, Juzeliūnas G and Öhberg P 2011 *Rev. Mod. Phys.* **83** 1523–43

[37] Goldman N, Juzeliūnas G, Öhberg P and Spielman I B 2014 *Rep. Prog. Phys.* **77** 126401

[38] Tai M E, Lukin A, Rispoli M, Schittko R, Menke T, Borgnia D, Preiss P M, Grusdt F, Kaufman A M and Greiner M 2017 *Nature* **546** 519–23

[39] Stuhl B K, Lu H I, Aycock L M, Genkina D and Spielman I B 2015 *Science* **349** 1514–8

[40] An F A, Meier E J and Gadway B 2017 *Science Advances* **3** e1602685

[41] Celi A, Massignan P, Ruseckas J, Goldman N, Spielman I B, Juzeliūnas G and Lewenstein M 2014 *Phys. Rev. Lett.* **112** 043001

[42] Barbarino S, Taddia L, Rossini D, Mazza L and Fazio R 2016 *New J. Phys.* **18** 035010

[43] Zhang X, Bishof M, Bromley S L, Kraus C V, Safronova M S, Zoller P, Rey A M and Ye J 2014 *Science* **345** 1467–73

[44] Pagano G et al 2014 *Nat. Phys.* **10** 198–201

[45] Mancini M et al 2015 *Science* **349** 1510–3

[46] Livi L F et al 2016 *Phys. Rev. Lett.* **117** 220401

[47] Narozhny B N, Carr S T and Nersesyan A A 2005 *Phys. Rev.* **B** **71** 161101

[48] Carr S T, Narozhny B N and Nersesyan A A 2006 *Phys. Rev.* **B** **73** 195114

[49] Roux G, Orignac E, White S R and Poilblanc D 2007 *Phys. Rev.* **B** **76** 195105

[50] Mazza L, Aidelsburger M, Tu H H, Goldman N and Burrello M 2015 *New J. Phys.* **17** 105001

[51] Barbarino S, Taddia L, Rossini D, Mazza L and Fazio R 2015 *Nat. Comm.* **6** 8134

[52] Budich J C, Laflamme C, Tschirsich F, Montangero S and Zoller P 2015 *Phys. Rev.* **B** **92** 245121

[53] Lacki M, Pichler H, Sterdyniak A, Lyras A, Lembessis V E, Al-Dossary O, Budich J C and Zoller P 2016 *Phys. Rev.* **A** **93** 013604

[54] Ghosh S K, Greschner S, Yadav U K, Mishra T, Rizzi M and Shenoy V B 2017 *Phys. Rev.* **A** **95** 063612

[55] Haller A, Rizzi M and Burrello M 2017 arXiv:1707.05715

[56] Bloch I, Dalibard J and Zwerger W 2008 *Rev. Mod. Phys.* **80** 885–964

[57] Aidelsburger M, Atala M, Nascimbène S, Trotzky S, Chen Y A and Bloch I 2011 *Phys. Rev. Lett.* **107** 255301

[58] Kollath C, Iucci A, McCulloch I P and Giamarchi T 2006 *Phys. Rev.* **A** **74** 041604

[59] Kollath C, Iucci A, Giamarchi T, Hofstetter W and Schollwöck U 2006 *Phys. Rev. Lett.* **97** 050402

[60] Iucci A, Cazalilla M A, Ho A F and Giamarchi T 2006 *Phys. Rev.* **A** **73** 041608

[61] Dalla Torre E G, Berg E and Altman E 2006 *Phys. Rev. Lett.* **97** 260401

[62] Berg E, Dalla Torre E G, Giamarchi T and Altman E 2008 *Phys. Rev.* **B** **77** 245119

[63] Dalla Torre E G 2013 *J. Phys. B: At. Mol. Opt. Phys.* **46** 085303

[64] Gerbier F and Dalibard J 2010 *New J. Phys.* **12** 033007

[65] Calabrese P and Cardy J 2004 *J. Stat. Mech.* **2004** P06002

[66] Arita R, Kuroki K, Aoki H and Fabrizio M 1998 *Phys. Rev.* **B** **57** 10324–7

[67] Moreno A, Muramatsu A and Manmana S R 2011 *Phys. Rev.* **B** **83** 205113

[68] Vidal G 2003 *Phys. Rev. Lett.* **91** 147902

[69] Vidal G 2004 *Phys. Rev. Lett.* **93** 040502

[70] Suzuki M 1976 *Prog. Theor. Phys.* **56** 1454

[71] Suzuki M 1991 *J. Math. Phys.* **32** 400

[72] Ludlow A D, Boyd M M, Ye J, Peik E and Schmidt P O 2015 *Rev. Mod. Phys.* **87** 637–701

[73] Yi W, Daley A J, Pupillo G and Zoller P 2008 *New J. Phys.* **10** 073015

[74] Folling S, Trotzky S, Cheinet P, Feld M, Saers R, Widera A, Muller T and Bloch I 2007 *Nature* **448** 1029–32

[75] Ruostekoski J, Dunne G V and Javanainen J 2002 *Phys. Rev. Lett.* **88** 180401

[76] Jaksch D and Zoller P 2003 *New J. Phys.* **5** 56

[77] Goldman N, Beugnon J and Gerbier F 2012 *Phys. Rev. Lett.* **108** 255303

[78] Lee J, Sachdev S and White S R 2016 *Phys. Rev.* **B** **94** 115112

[79] Rodney M, Song H F, Lee S S, Le Hur K and Sørensen E S 2013 *Phys. Rev.* **B** **87** 115132

[80] Zhuang Y, Changlani H J, Tubman N M and Hughes T L 2015 *Phys. Rev.* **B** **92** 035154