# Lattice modulation spectroscopy of one-dimensional quantum gases: Universal scaling of the absorbed energy

Roberta Citro, Eugene Demler, Thierry Giamarchi, Michael Knap, Edmond Orignac

► To cite this version:

Roberta Citro, Eugene Demler, Thierry Giamarchi, Michael Knap, Edmond Orignac. Lattice modulation spectroscopy of one-dimensional quantum gases: Universal scaling of the absorbed energy. Physical Review Research, 2020, 2 (3), pp.033187. 10.1103/PhysRevResearch.2.033187. hal-02916922

HAL Id: hal-02916922

https://hal.science/hal-02916922v1

Submitted on 4 Sep 2024

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
Lattice modulation spectroscopy of one-dimensional quantum gases:
Universal scaling of the absorbed energy

R. Citro,¹,² E. Demler,³ T. Giamarchi,⁴ M. Knap,⁵,⁶ and E. Orignac ⁷

¹Dipartimento di Fisica “E.R. Caianiello,” Università degli Studi di Salerno and CNR-SPIN c/o University of Salerno,
Via Giovanni Paolo II 132, 84084 Fisciano, Salerno, Italy

²INFN, Sezione di Napoli, Gruppo Collegato di Salerno, 84084 Fisciano, Salerno, Italy

³Department of Physics, Harvard University, Cambridge, Massachusetts 02138, USA

⁴DQMP, University of Geneva, 24 Quai Ernest-Ansermet, 1211 Geneva, Switzerland

⁵Department of Physics and Institute for Advanced Study, Technical University of Munich, 85748 Garching, Germany

⁶Munich Center for Quantum Science and Technology, Schellingstraße 4, 80799 München, Germany

⁷University of Lyon, ENS de Lyon, Université Claude Bernard, CNRS, Laboratoire de Physique, 69342 Lyon, France

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

(Received 10 March 2020; accepted 25 June 2020; published 3 August 2020)

Lattice modulation spectroscopy is a powerful tool for probing low-energy excitations of interacting many-body systems. By means of bosonization we analyze the absorbed power in a one-dimensional interacting quantum gas of bosons or fermions, subjected to a periodic drive of the optical lattice. For these Tomonaga-Luttinger liquids we find a universal $\omega^3$ scaling of the absorbed power, which at very low-frequency turns into an $\omega^2$ scaling when scattering processes at the boundary of the system are taken into account. We confirm this behavior numerically by simulations based on time-dependent matrix product states. Furthermore, in the presence of impurities, the theory predicts an $\omega^2$ bulk scaling. While typical response functions of Tomonaga-Luttinger liquids are characterized by exponents that depend on the interaction strength, modulation spectroscopy of cold atoms leads to a universal power-law exponent of the absorbed power. Our findings can be readily demonstrated in ultracold atoms in optical lattices with current experimental technology.

DOI: 10.1103/PhysRevResearch.2.033187

I. INTRODUCTION

Cold atomic systems offer an unprecedented level of control of the properties of interacting quantum systems [1,2] and allow for the realization of a plethora of novel phases and phenomena that were previously inaccessible in other experiments. They have given access to the “experimental solution” of certain models and hence can be referred to as quantum simulators. Among those phenomena, one paradigmatic example with bosons on a lattice is the transition between a superfluid and Mott insulator state. Such a transition was successfully observed in three dimensions [3], two dimensions [4], and one dimension [5,6]. For the latter, the transition is found to be in the universality class of the Berezinskii-Kosterlitz-Thouless transitions. Cold atoms have thus provided a remarkable way of testing the universal properties of such models.

In order to analyze quantitatively the properties of the
correlated phases and the transitions between them, it is
important to develop a detailed understanding of different
experimental probes. Among them is lattice modulation spec-
troscopy [7]. This technique consists of modulating the ampli-
tude [8–10] or the phase [11,12] of the optical lattice in which

the atoms are trapped. After some time, the energy deposited
in the system or the number of doubly occupied states is
measured to characterize the underlying state [13,14]. In
general, driving of the hopping energy provides a novel form
of Floquet engineering, which enables atypical Hamiltonians
and exotic states of matter to be produced and controlled
[15,16].

For bosons this probe can determine energy gaps and
thus locate the Mott-to-superfluid transition [5-7]. Moreover,
specific modes of the superfluid such as the Higgs mode
could be excited [4], which is expected to occur in three- and
two-dimensional superfluids.

Despite this effort, several questions remain to be in-
vestigated for bosons and fermions in one dimension. For
these systems, a symmetry-broken state cannot exist because
of strong phase fluctuations even at zero temperature [17].
Hence, only quasi-long-range order can exist as characterized
by a power-law decay of certain correlation functions. This
result is part of the more general properties of Tomonaga-
Luttinger liquids (TLLs) that are expected to describe most
of the interacting one-dimensional (1D) quantum problems
[17,18]. Given the absence of true long-range superfluid order,
one may wonder whether the response to shaking in a one-
dimensional bosonic system would also show traces of a
Higgs mode as in higher dimensions [4]. More generally,
this prompts an analysis of the response to shaking of a
one-dimensional TLL.

In the present paper we perform such an analysis in both
the gapless (superfluid) and the gapped phase (Mott insulator).

Published by the American Physical Society under the terms of the
Creative Commons Attribution 4.0 International license. Further
distribution of this work must maintain attribution to the author(s)
and the published article's title, journal citation, and DOI.
Using a combination of field-theoretic and numerical matrix product state (MPS) techniques, we obtain the response of the system to the shaking of the optical lattice. We find that this response is a power law of the shaking frequency, with a universal exponent. This is quite remarkable in view of the fact that the TLL is normally characterized by responses which show nonuniversal power-law behavior, with exponents depending on the interaction strength. The choice of modulating the amplitude of the optical lattice is important. Would one modulate the phase instead, the conductivity would be obtained [11] (as periodic phase modulation translates to a periodic force), yielding nonuniversal exponents, while the amplitude modulation modifies the tunneling and consequently the kinetic energy of the system.

Our results also show that a well-formed Higgs mode does not exist in one dimension, as the response to shaking is monotonically increasing up to shaking frequencies of the order of the bandwidth of the system. In order to connect to experiments, we also analyze and discuss the effect of this response to open boundary conditions, as realized, for example, with a box potential, which constitute a relevant perturbation at low frequencies. In the presence of impurities similar effects are observed, provided the concentration of impurities is small and one can simply add energy absorption due to different impurities.

While the focus of this paper is on one-dimensional systems, we note that a similar approach can be used to describe modulation experiments in higher-dimensional systems that allow for a hydrodynamic description. Time-periodic modulation of the interaction strength or the transverse confinement potential will result in the resonant parametric generation of excitation pairs [19–22], analogous to the 1D case discussed in this paper. The analysis of the absorbed energy can be done following the same approach that we discuss here, with the main difference being the phase space for collective modes. In Sec. VI we comment on the relevance of our analysis beyond ensembles of cold atoms and suggest possible applications of our results to chains of Josephson junctions and pump and probe experiments in electron systems.

The paper is organized as follows. In Sec. II we introduce the models that we are investigating. In Sec. III we discuss the analytical calculation of the absorbed power of a one-dimensional gas subject to a lattice modulation within a Tomonaga-Luttinger liquid treatment and compare the low-energy behavior with results obtained from time-dependent matrix product states. In Sec. IV we discuss the edge effects treated through an effective boundary potential that couples to the density of the fluid and also analyze the effect of a single impurity in the bulk of the system. In Sec. V we consider modulation spectroscopy in the case of gapped systems such as would occur for the Bose- and the Fermi-Hubbard model with repulsive interactions and commensurate filling. In Sec. VI we present a summary and discuss our findings.

## II. MODELS

We consider fermionic or bosonic ultracold atoms confined to a 1D tube subjected to a deep lattice potential. For deep enough potentials, such a system can be described in a tight-binding approximation by a Hubbard model [23]. This leads

to the Bose-Hubbard model for spinless bosons [1]

$$H_b^0 = \sum_l \left[ -J_0(b_{l+1}^\dagger b_l + b_l^\dagger b_{l+1}) + \frac{U}{2} n_l(n_l - 1) \right], \quad (1)$$

where $b_l^\dagger (b_l)$ creates (annihilates) a particle on site $l$ and $n_l = b_l^\dagger b_l$ counts the particles on site $l$, and the Fermi-Hubbard model [2] for spin-1/2 fermions

$$H_f^0 = \sum_{l,\sigma} \left[ -J_0(c_{l+1,\sigma}^\dagger c_{l,\sigma} + c_{l,\sigma}^\dagger c_{l+1,\sigma}) + \frac{U}{2} n_{l,\sigma} n_{l,-\sigma} \right], \quad (2)$$

where again $c_{l,\sigma}^\dagger$ and $c_{l,\sigma}$ are the creation and annihilation operators, respectively.

In modulation spectroscopy, the trapped atoms are probed by modulating the strength of the longitudinal periodic potential. The modulation lowers or raises the potential barrier between two consecutive minima, and thus to leading order modifies the strength of the tunneling $J_0$ as well as the interaction within one well. Modulation of $J$ is expected to be much larger since it depends exponentially on the barrier height. We note that modulation of the Hamiltonian as a whole does not lead to energy absorption and what is important is the difference in the relative modulation strength of the two terms in the Hamiltonian. Therefore, we consider the time-dependent Hamiltonian, in which only the tunneling amplitude is modulated $J_0 + \delta J(t)$, giving rise to

$$H_v(t) = H_v^0 + \delta J(t) O_v, \quad (3)$$

with $v \in b, f$ and

$$O_b = \sum_l (b_{l+1}^\dagger b_l + b_l^\dagger b_{l+1}) \quad (4)$$

for bosons and

$$O_f = \sum_{l,\sigma} (c_{l+1,\sigma}^\dagger c_{l,\sigma} + c_{l,\sigma}^\dagger c_{l+1,\sigma}) \quad (5)$$

for fermions. The labels $f$ and $b$ refer to fermions and bosons, respectively.

We work in the linear-response limit, with $\delta J(t) = \delta J \cos(\omega t)$ and $\delta J \ll J_0, U$. In the linear response, the absorbed power is given by [24]

$$\mathcal{P} = \frac{\omega}{2} |\delta J|^2 \mathrm{Im}\chi_v(\omega), \quad (6)$$

where

$$\chi_v(\omega) = i \int_0^{+\infty} e^{i\omega t} \langle [O_v(t), O_v(0)] \rangle \quad (7)$$

is the response function. In Sec. III we calculate this response function in a low-energy and long-wavelength limit.

## III. BOSONIZATION

In the low-energy or long-wavelength limit, interacting bosons and fermions can be described within an effective continuum theory called the Tomonaga-Luttinger liquid [17,18,25–33]. The low-energy excitations are phononic collective excitations with a linear dispersion that describe density fluctuations (and when applicable spin fluctuations). In general, spin and density fluctuations propagate with different
velocities, a phenomenon known as spin-charge separation [34,35]. In that low-energy limit, the original particles appear as coherent states of the collective modes. As a result, all the observables of the original system are expressible in terms of the collective modes. We will thus use in this section the bosonized representation to calculate the response function (7) and hence the resulting absorbed power.

### A. Bosonized representation for fermions

In the fermionic case, away from commensurate filling, the Hamiltonian (2) has the bosonized representation [17,29]

$$H_f = H_\rho + H_\sigma, \qquad (8)$$

$$H_\rho = \int \frac{dx}{2\pi} \left[ u_\rho K_\rho (\pi \Pi_\rho)^2 + \frac{u_\rho}{K_\rho} (\partial_x \phi_\rho)^2 \right], \qquad (9)$$

$$H_\sigma = \int \frac{dx}{2\pi} \left[ u_\sigma K_\sigma (\pi \Pi_\sigma)^2 + \frac{u_\sigma}{K_\sigma} (\partial_x \phi_\sigma)^2 \right] - \frac{2g_{1\perp}}{(2\pi\alpha)^2} \int dx \cos \sqrt{8}\phi_\sigma, \qquad (10)$$

where $u_\rho$ and $u_\sigma$ are the density and spin velocity, respectively, and $K_\rho$ and $K_\sigma$ are the density and spin Tomonaga-Luttinger exponents, respectively, where for small U one has $u_\rho = v_F \sqrt{1+U/\pi v_F}$, $u_\sigma = v_F \sqrt{1-U/\pi v_F}$, $K_\rho = \frac{1}{\sqrt{1+U/\pi v_F}}$, $K_\sigma = \frac{1}{\sqrt{1-U/\pi v_F}}$, $g_{1\perp} = U$, and $[\phi_\nu(x), \Pi_\nu(x')] = i\delta(x-x')\delta_{\nu,\nu'}$. A general expression can be obtained for arbitrary U as discussed in [17,29]. For repulsive interactions, $H_\sigma$ is renormalized to a fixed-point Hamiltonian $H_\sigma^*$ with $K_\sigma^* = 1$ and $g_{1\perp}^* = 0$, yielding gapless excitations with linear dispersion $\omega = u_\sigma |k|$. For attractive interactions without external magnetic field the spin Hamiltonian $H_\sigma$ is gapped while the density Hamiltonian $H_\rho$ remains gapless [17]. At half filling, umklapp processes are present [17]. They contribute to the bosonized Hamiltonian (8) a term

$$H_{\text{umk}} = -\frac{4g_3}{(2\pi\alpha)^2} \int dx \cos \sqrt{8n}\phi_{\rho}, \qquad (11)$$

but, as shown in Appendix A, when umklapp processes are irrelevant in the renormalization group sense [36], they add only a subdominant contribution to the absorbed power at low frequency. When the umklapp processes are relevant, they open a gap in the spectrum.

In the perturbative limit, since $O_f$ is proportional to the kinetic energy in Eq. (2), its bosonized form is simply the bosonized Hamiltonian of noninteracting spinless fermions divided by $J_0$. Changing to the spin and charge fields, we find [17,29]

$$O_f = \sum_{\nu=\rho,\sigma} O_{f,\nu} = O_{f,\rho} + O_{f,\sigma},$$

$$O_{f,\nu} = 2a \sin(k_F a) \int \frac{dx}{2\pi} [(\pi \Pi_\nu)^2 + (\partial_x \phi_\nu)^2]. \qquad (12)$$

Due to spin-charge separation, the absorbed power is the sum of a spin and a density contribution. To find an expression of $O_f$ applicable away from the perturbative limit, we note that $O_f$ is obtained by differentiating the Fermi-Hubbard Hamiltonian (2) with respect to $J_0$. Assuming that the identity

carries over to the bosonized description, we have

$$O_{f,\rho} = \int \left[ \frac{\partial(u_\rho K_\rho)}{\partial J_0} (\pi \Pi_\rho)^2 + \frac{\partial}{\partial J_0} \left( \frac{u_\rho}{K_\rho} \right) (\partial_x \phi_\rho)^2 \right] \frac{dx}{2\pi} \qquad (13)$$

and a similar approximation for $O_{f,\sigma}$. As a further approximation, in the repulsive case, we take the fixed-point values in $H_\sigma$ and write

$$O_{f,\sigma} = \int \frac{\partial u_\sigma}{\partial J_0} [(\pi \Pi_\sigma)^2 + (\partial_x \phi_\sigma)^2] \frac{dx}{2\pi}. \qquad (14)$$

Applying Eqs. (13) and (14) in the perturbative case, Eq. (12) is recovered. It is important to note that the full expression of the fermion kinetic energy contains, besides the linear dispersion valid near the Fermi points, corrections coming from band curvature. So the expressions (13) and (14) are really the most relevant terms in an expansion of the operator $O_f$ in a series of operators of increasing scaling dimensions. The contributions of operators of higher scaling dimensions are subdominant at low frequency, as shown in Appendix A.

### B. Bosonized representation for bosons

In this section we turn to bosons. The Bose-Hubbard Hamiltonian with $U > 0$ has the bosonized representation

$$H = \int \frac{dx}{2\pi} \left[ uK(\pi\Pi)^2 + \frac{u}{K}(\partial_x\phi)^2 \right], \qquad (15)$$

where $\Pi(x)$ and $\phi(x)$ are conjugate operators that describe the boson density fluctuations, $u$ is their velocity, and $K$ is the Tomonaga-Luttinger exponent [25] that to the lowest-order approximation is determined by $uK = \frac{\pi\rho_0}{m}$ and $\frac{u}{K} = \frac{U}{\pi}$, where $\rho_0$ is the boson density, while its dependence on general values of the interaction can be found in [18,37–39]. In Appendix B the Luttinger parameters of the Bose-Hubbard model as a function of the interaction U and system size L are shown. At integer filling, umklapp processes contribute a term proportional to $\cos 2\phi$ to the Hamiltonian (15). Similarly to the fermionic case, their contribution is subdominant as long as the system remains in a Tomonaga-Luttinger liquid ground state. The limit $U \to 0$ of the Hamiltonian (15) is singular, with the velocity $u$ vanishing to recover the quadratic dispersion of noninteracting bosons above a condensate and the Tomonaga-Luttinger exponent going to $+\infty$. Thus, in contrast to the fermionic case of Sec. III A, it is impossible to derive a bosonized representation of (4) by considering the noninteracting limit. However, assuming as in Sec. III A that the identity $O_b = \frac{\partial H_b}{\partial J_0}$ is applicable to the bosonized Hamiltonian (15), we find

$$O_b = \int \frac{dx}{2\pi} \left[ \frac{\partial(uK)}{\partial J_0} (\pi\Pi)^2 + \frac{\partial}{\partial J_0} \left( \frac{u}{K} \right) (\partial_x\phi)^2 \right]. \qquad (16)$$

This expression is similar to (12). Moreover, in the hard-core limit $U \to +\infty$, bosons can be mapped to noninteracting spinless fermions [40] and the fermionic expression (12) yields an explicit form of $O_b$ which fully agrees with (16). As we discussed in the fermionic case, the expression (16) is only the first term in a series of operators of increasing scaling dimension that represent the various band curvature terms coming from the dispersion of the lattice model.
### C. Response function in an infinite system

With repulsive interactions, both for fermions and for bosons, the calculation of the response function (7) reduces to the calculation of the response function of an operator of the form $\int dx[A\Pi^2 + B(\partial_x\phi)^2]$ for a Hamiltonian quadratic in $\Pi$ and $\partial_x\phi$. That calculation is further simplified by rewriting the bosonized form of the operator $O_{b,f}$ as a linear combination of the Hamiltonian and an operator proportional to $\int (\partial_x\phi)^2$. In the bosonic case,

$$O_b = \frac{1}{uK} \frac{\partial(uK)}{\partial J_0} H - \int \frac{dx}{\pi} \frac{u}{K^2} \frac{\partial K}{\partial J_0} (\partial_x\phi)^2, \quad (17)$$

and in the fermionic case, for the perturbative limit,

$$O_{f,v} = 2a \sin(k_F a) \left[ \frac{H_v}{u_v K_v} + \left(1 - \frac{1}{K_v^2}\right) \int dx (\partial_x \phi_v)^2 \right], \quad (18)$$

while in the nonperturbative limit,

$$O_{f,v} = \frac{\partial}{\partial J_0}(u_v K_v) \left[ \frac{H_v}{u_v K_v} - \int \frac{u_v}{K_v^2} \frac{\partial K_v}{\partial J_0} (\partial_x \phi_v)^2 \right]. \quad (19)$$

Since the Hamiltonian is time independent, the response function reduces, up to a proportionality factor, to the one of $\int dx(\partial_x\phi)^2$. We note that this is the same response function as in the case where the on-site interaction is modulated. Furthermore, according to Eq. (19), the response function (7) vanishes for noninteracting fermions since $K_\rho = K_\sigma = 1$ for any $J_0$ in that case. This can be established more directly from the lattice Hamiltonian by noting that for $U = 0$, $O_f$ is proportional to the Hamiltonian. More importantly, Eq. (19) also shows that the contribution of the spin excitations calculated at the fixed point $K_\sigma^* = 1$ is vanishing. This indicates that for interacting fermions the dominant contribution comes from the density response. Due to the fact that the drive is coupling only to the density and not to the spin, this is expected to be the case on general grounds. Similarly, in the bosonic case, in the limit $U \to +\infty$, where $K = 1$ for all $J$, the response function (7) is also vanishing. Again, this is more directly established by noting that $O_b$ is directly proportional to the hard-core boson Hamiltonian in that limit.

We calculate $\chi(\omega)$ by taking the analytic continuation $\chi(\omega) = \chi_M(i\omega_n \to \omega + i0_+)$ of the Matsubara correlation function

$$\chi_M(i\omega_n) = \int d\tau e^{i\omega_n \tau} \langle T_\tau O_v(\tau) O_v(0) \rangle. \quad (20)$$

For the sake of definiteness, we perform the calculation for bosons. Using translational invariance, we find that

$$\frac{1}{L} \chi_M(i\omega_n) = \left( \frac{u}{\pi K^2} \frac{\partial K}{\partial J_0} \right)^2 \int dx d\tau \times e^{i\omega_n \tau} \langle T_\tau (\partial_x \phi)^2 (x, \tau) (\partial_x \phi)^2 (0, 0) \rangle. \quad (21)$$

Details on the evaluation of $\chi(\omega)$ can be found in Appendix C and for zero temperature the final result is

$$\frac{1}{L} \mathrm{Im} \chi(\omega) = \mathfrak{F}(K) \omega^2 e^{-|\omega|\alpha/u} \mathrm{sgn}(\omega), \quad (22)$$

where $\alpha$ is a short distance cutoff (of the order of the lattice spacing) and $\mathfrak{F}(K) = \frac{1}{16u} \frac{1}{K^2} (\frac{\partial K}{\partial J_0})^2$. Only the behavior

<figure><img src="image_3.png" /><figcaption>FIG. 1. Response function $\chi(\omega)$ to a modulation of the lattice at the frequency $\omega$ and temperature $T$ evaluated from Eq. (20) for $F(K) = 1$ and $\alpha/u = 1$.</figcaption></figure>

for $|\omega| \ll u/a \sim \delta J$ is reliably predicted by bosonization. For frequencies of order of the bandwidth, the linearized approximation for the dispersion certainly breaks down and high-energy excited states not described by bosonization can contribute as well to the energy absorption.

At finite temperature, Eq. (22) becomes

$$\frac{1}{L} \mathrm{Im} \chi(\omega) = \mathfrak{F}(K) \omega^2 e^{-|\omega|\alpha/u} \coth\left(\frac{\omega}{4T}\right), \quad (23)$$

so the response function behaves as $\sim \omega T$ when $\omega \ll T$ and as $\omega^2$ when $T \ll \omega$ (see Fig. 1). Thus the absorbed power is

$$\mathcal{P}_b = \frac{L}{32u} |\omega|^3 \left( \frac{\delta J}{K} \frac{\partial K}{\partial J_0} \right)^2 \coth \left( \frac{\omega}{4T} \right) e^{-|\omega|\alpha/u} \quad (24)$$

for bosons and

$$\mathcal{P}_f = \frac{L}{32u\rho} |\omega|^3 \left( \frac{\delta J}{K_\rho} \frac{\partial K_\rho}{\partial J_0} \right)^2 \coth \left( \frac{\omega}{4T} \right) e^{-|\omega|\alpha/u_\rho} \quad (25)$$

for fermions.

It has a universal power-law dependence on frequency, with an exponent independent of interactions. This universal behavior has to be contrasted with the conductivity [17,41] where the power-law exponent varies with the Tomonaga-Luttinger parameter and thus depends on the microscopic interaction strength. Here only the prefactor depends on the logarithmic derivative of the Tomonaga-Luttinger parameter with respect to the hopping amplitude. In Appendix B the dependence of this prefactor on system size and interaction $U$ is reported.

### D. Numerical results

In order to elucidate the universal frequency exponent of the absorbed power density $\mathcal{P}_b/L \sim \delta J^2 \omega^3$ predicted from the Tomonaga-Luttinger theory, we numerically evaluate the energy absorption in the Bose-Hubbard model (1) using matrix product states [42,43]. In particular, we consider systems with 120 sites and noninteger boson density $\rho = 1.2$, to fully avoid
<figure><img src="image_4.png" /><figcaption>FIG. 2. Absorbed power density. Using matrix product states, we have evaluated the absorbed power in a periodically driven Bose-Hubbard model on an open chain of L = 120 sites and density ρ = 1.2 for two values of the interaction strength (see the legend). The absorbed power density, renormalized by the drive strength δJ $^{2}$, shows a universal ω 3 scaling, as predicted from the Tomonaga-Luttinger theory. Deviations at low frequencies are decreasing with increasing system size and are expected to arise from the residual contributions of system edges which add an ω $^{2}$/L contribution that is leading in frequency but vanishing in the thermodynamic limit.</figcaption></figure>

umklapp processes. We have checked the convergence of our results with the bond dimension of the matrix product state which ranges from χ = 400 to 800.

Our objective is to simulate an experimental protocol to measure the absorption. To this end, we first compute the ground state of our model and then apply a periodic modulation of the kinetic energy of the form J ( t ) = J$_{0}$ + δJ sin( ωt ). We choose the driving strength δJ = 0 . 1 J$_{0}$ , small enough such that the absorbed energy increases linearly in time, as required from the linear-response theory. We evolve the system for a few drive periods and extract the absorbed power density for a range of modulation frequencies (see Fig. 2). The power density scales as ω $^{3}$, in agreement with the Tomonaga- Luttinger liquid prediction. At low frequencies there are small deviations from the predicted scaling, as expected from a contribution from boundaries (see Sec. IV A).

IV. BROKEN TRANSLATIONAL SYMMETRY

So far, our considerations have been restricted to an infinite system without defects. This section focuses instead on the case of systems with broken translational symmetry caused by either boundaries or impurities.

A. Effects of boundaries

Since trapped atoms are systems of finite length, the ef- fect of boundaries on their response must in principle be considered. We examine in this section the effect of edge potentials that pin the density and this can potentially modify the response to shaking.

The bosonized Hamiltonian in the presence of forward- scattering edge potentials becomes [17,44–47]

$$
H = \int_0^L \frac{dx}{2\pi} \left[ uK(\pi\Pi)^2 + \frac{u}{K}(\partial_x\phi)^2 \right] \\
- \frac{V}{\pi} [\partial_x\phi(0) + \partial_x\phi(L)], \qquad (26)
$$

with the Dirichlet boundary conditions φ (0) = 0 and φ ( L ) = − πN . Those boundary conditions ensure that no current can leak through the edges of the system. As ρ ( x ) = − ∂$_{x}$φ /π , the terms V simply represent a forward scattering in the vicinity of the system edges. Note that with the Dirichlet boundary conditions, a backscattering term − V$_{b}$ cos 2 φ (0) can be re- duced to a forward-scattering term [17], so there is no loss of generality in Eq. (26). Since in bosonization the particle- hole symmetry is φ → − φ and Π → − Π, V vanishes in a particle-hole symmetric system [46]. In the absence of such symmetry however, those terms can be nonzero. When one considers only the static properties, the edge potential can be eliminated by modifying the Dirichlet boundary conditions [47]. However, when we modulate the lattice, the edge poten- tial can be time dependent V = V ( J ( t )), and for that reason, it is better to retain the original boundary conditions. When we differentiate the Hamiltonian (26) with respect to J$_{0}$ , as in Eq. (16), the edge potential in (26) gives an extra edge contribution proportional to ∂$_{J}$$_{0}$V to the operator O$_{b}$ ,

$$
O_b = \int_0^L \frac{dx}{2\pi} \left[ \frac{\partial}{\partial J_0} (uK)(\pi\Pi)^2 + \frac{\partial}{\partial J_0} \left(\frac{u}{K}\right) (\partial_x\phi)^2 \right] \\ - \frac{1}{\pi} \frac{\partial V}{\partial J_0} [\partial_x\phi(0) + \partial_x\phi(L)]. \qquad (27)
$$

The response coming from the edge potential is calculated in Appendix D. It contributes

$$
\mathcal{P}_{\text{edge}} = \left(\frac{V}{u}\right)^2 \left[ \delta J \frac{\partial}{\partial J_0} \left( \ln \frac{VK}{u} \right) \right]^2 \frac{\omega^2}{\pi} \quad (28)
$$

to the absorbed power. The total absorbed power is therefore P$_{tot}$ = P$_{edge}$ + P$_{bulk}$ ∼ ω 2 + Lω $^{3}$. The edge response domi- nates below a crossover frequency ω ∗ ∼ 1 /L . The boundary potential V remains to be determined. A possible approach is to see how the Friedel oscillations are affected by these scattering potentials at the boundary.

B. Friedel oscillations and determination of the edge potentials

The edge potential V in Eq. (26) can be deduced from the Friedel oscillations [48–50] in the density profile of the ground state. The explicit calculation of the density profile in Appendix E leads us to an expression valid sufficiently far from edges

$$
\langle \rho(x) \rangle \sim \frac{\cos(2k'_F x - \varphi)}{\left[ \sin\left(\frac{\pi x}{L}\right) \right]^K}, \qquad (29)
$$

where 2 k ′ F = 2 k$_{F}$ + 4 KV uL and ϕ = 2 KV u , with k$_{F}$ = πN$_{0}$/L the nominal Fermi wave vector of the Friedel oscillations in a sys- tem of length L containing N$_{0}$ bosons. From (29), consecutive zeros of the Friedel oscillations are separated by the distance π 2 k ′ F = L 2 N$_{0}$ − L N 2 0 πu + O ( L/N 3 $_{0}$), instead of L 2 N$_{0}$ , revealing the
presence of the edge potential. In the thermodynamic limit, $2k'_F$ reduces to $2k_F$. However, the phase shift $\varphi = 2KV/u$ persists and the fitted expression of the Friedel oscillations obtained by MPS techniques reveals the presence or absence of a potential near the edge.

### C. Effects of a single impurity

Let us finally consider a single impurity located at $x_0$ whose potential energy is given by $H_{imp} = V\rho(x)\delta(x - x_0)$. Within the bosonization approach this term gives rise to two terms in the Luttinger liquid Hamiltonian: a term $-\frac{1}{\pi}\partial_x\phi(x_0)$, which corresponds to a forward-scattering process, and a term proportional to $\cos[2\phi(x_0)]$, which corresponds to backscattering. Using the same treatment as above, the first term will be leading to a dominant $\omega^2$ scaling function for the absorbed power, while the backscattering will contribute a term proportional to $\omega^{2(K-1)}$ with $K > 1$ and thus less relevant at low frequency. Thus the presence of a single impurity would lead to a dominant $\omega^2$ contribution to the absorbed power.

## V. GAPPED SYSTEMS

In the case of fermions with attractive interactions or in the case of fermions or bosons with repulsive interactions at commensurate filling, the spectrum can become gapped. The response in that gapped regime can be calculated either in the Luther-Emery limit [17,51] or in the more general case using the form-factor expansion [52]. Both methods predict a threshold in absorption power at the gap.

For the Bose-Hubbard model, below that gap the power absorption will be zero. In the Fermi-Hubbard case with repulsive interactions, we have seen that the spin response was suppressed, so only the density response contributed. In the gapped state, the density response also does not contribute at frequency lower than the gap, making the threshold observable as well. In the case of the Fermi-Hubbard model with attractive interaction, since the density modes are gapless, the response at low frequency will be the $\sim |\omega|^3$ contribution. The threshold at the gap then appears as a cusp-shaped rapid increase of absorption. The physical interpretation of such threshold is quite simple. At frequencies lower than the binding energy of two fermions of opposite spins, the pairs of fermions behave as an interacting boson gas [27], yielding the $\sim |\omega|^3$ contribution to the absorbed power. As the frequency becomes comparable to the binding energy of the pair, another absorption channel from dissociation of the pairs becomes available, leading to the rapid increase of absorption.

For concreteness, let us first consider the Fermi-Hubbard model in the Mott insulating phase for the particular case of the Luther-Emery limit where $K_\rho = 1/2$. In that limit, the resulting absorbed power is given by (see Appendix F for a detailed calculation)

$$ \mathcal{P}_{\text{LE}} = \frac{L\Delta^2}{u} \left( \frac{5a\delta J}{2u} \sin(k_F a) \right)^2 \sqrt{\left(\frac{\omega}{2}\right)^2 - \Delta^2}, \quad (30) $$

where $2\Delta < \omega < 4\Delta$, leading to a cusp singularity at $\omega = 2\Delta$.

This analysis can be extended away from the Luther-Emery point to any value of the Luttinger parameter by the

form-factor expansion for the sine-Gordon model [52–55], as detailed in Appendix G (see also [56]). The main result for the absorbed power (with $1/2 < K_\rho < 1$) is

$$ \mathcal{P}_{\text{FF}} = L \left( \frac{2K_{\rho}}{K_{\rho}^2 + 1} \right)^2 \frac{2M_s^2 |2a \sin(k_F a) \delta J|^2}{\pi u_{\rho}^3 \nu^2} \sqrt{\omega^2 - 4M_s^2}, \quad (31) $$

where we have found a threshold at twice the mass of the soliton $M_s$ and $\nu = K_\rho/(1 - K_\rho)$ for the Fermi-Hubbard model. For $K_\rho < 1/2$, besides the threshold behavior (31), discrete peaks coming from resonant absorption by soliton-antisoliton bound states become possible. This behavior could be readily observed in current experiments with cold atoms in the Mott insulating regime [6,57].

The same threshold behavior as in Eqs. (30) and (31) was also obtained in the opposite limit of a weak lattice [56] in which the depth of the periodic potential was modulated. One may thus speculate whether such threshold behavior is also observed for intermediate lattice strengths.

## VI. SUMMARY AND OUTLOOK

We have analyzed in the linear response the power absorbed by one-dimensional fermions and bosons in the Tomonaga-Luttinger liquid [25] or Luther-Emery liquid [51] phase, for the amplitude modulation of an optical lattice. In the Tomonaga-Luttinger liquid, we have found that the absorbed power possesses a universal $\omega^3$ power-law onset, which has been confirmed by numerical simulations based on matrix product states. We have also shown that this power law crosses over to $\sim \omega^2$ at low frequency in finite systems, when edge effects are taken into account. A similar $\omega^2$ behavior is found for systems with a single impurity located in the bulk.

Such universal behavior is surprising since in the Tomonaga-Luttinger liquid theory, response functions usually show nonuniversal exponents determined by the interaction strength [58]. The universal $\omega^3$ scaling of the absorbed power can be readily measured for ultracold atoms in optical lattices confined to one dimension by measuring the energy change over time. In Luther-Emery liquid phases, which can be obtained for commensurate densities or with spin-1/2 fermions having an attractive interaction, the absorbed power vanishes below a gap and shows a marked onset above, thus making it possible to identify this energy scale.

The discussion in this paper focused on experiments with spinless ultracold atoms. Before concluding this section we briefly review other systems in which ideas developed in this paper can be tested experimentally.

Bosonic spin mixtures in optical lattices can be used to realize lattice spin Hamiltonians and spinor condensates [59–61]. Recent experiments by Jepsen et al. [62] used the magnetic field dependence of the interspecies scattering length to realize XXZ spin chains with tunable anisotropy of interactions. In the regime of easy plane anisotropy XXZ chains are in the gapless regime, while the easy axis case corresponds to the gapped regime. Periodic modulation of $J_z/J_\perp$ can be achieved in this system through periodic modulation of the magnetic field and should have an effect equivalent to modulation of the interaction strength for spinless
bosons. These experiments have high local resolution, which will allow one to spatially resolve spin patterns induced by modulation of the interaction anisotropy. Hence predictions of our paper for both gapless and gapped regimes can be checked experimentally.

Recent progress in superconducting nanotechnology makes it possible to engineer arrays of coupled Josephson junctions whose parameters can be controlled dynamically. Lähteenmäki et al. [63] demonstrated a dynamical Casimir effect in a one-dimensional chain of Josephson junctions, in which the Josephson energy of the junctions was modulated by periodically changing the background magnetic flux. Parametric generation of photons at half the modulation frequency observed in these experiments is the direct analog of energy absorption in the Luttinger liquid discussed in our paper. Recent experiments by Kuzmin et al. [64] demonstrated the possibility of tuning a chain of Josephson junctions through the superconductor to insulator transition and explored evolution of the collective phase mode across the transition. Hence 1D superconducting metamaterials make it possible to study modulation spectroscopy of 1D systems in both gapless and gapped phases.

Although the focus of this paper has been on one-dimensional systems, a similar analysis can be applied to study periodic driving of higher-dimensional systems provided their lower-energy excitations allow a field-theoretic description. Modulation of the kinetic energy of bosons in optical lattices has been considered in the context of the Higgs mode in systems with broken U(1) symmetry [4,65,66]. In the superfluid phase in $d = 2, 3$ close to the critical point the imaginary part of the response function of the operator of kinetic energy develops a broad peak at the energy equal to the Higgs mode frequency and has a universal scaling form proportional to $\omega^{d+1}$ at smaller frequencies. The latter is determined by the process of resonant excitation of pairs of Goldstone modes with opposite momenta mediated by the virtual excitation of the Higgs mode. This process is equivalent to the mechanism of exciting pairs of Luttinger liquid phonons considered in our paper for one-dimensional systems. Thus energy absorption rate at low frequencies has a general scaling form $\omega^{d+2}$. We also note that our formalism should be useful for analyzing pump and probe experiments in interacting electron systems [67–71]. Recent experiments by

von Hoegen et al. [72] have observed parametric excitation of Josephson plasmons in YBCO superconductors following resonant excitation of apical oxygen phonons. The microscopic mechanism of phonon-plasmon coupling is modulation of the superfluid density in copper-oxide planes by the phonon-induced motion of oxygen atoms. Analogously to what we have discussed in this paper, resonant parametric excitation of plasmon pairs has been a crucial component of experiments by von Hoegen et al. One important difference, however, is that three-wave mixing between phonons and plasmons involves two different types of plasmons, the so-called lower and upper Josephson plasmons. The formalism developed in our paper can be extended to the case of parametric instabilities involving different types of collective excitations. We expect that resonant parametric interactions between phonons and collective excitations of many-body electron systems should be a ubiquitous phenomenon. Excited phonons can modulate several parameters of electron systems, including effective mass, interactions, and carrier density. Thus pump and probe experiments can be used to achieve parametric driving of a broad range of collective modes, including plasmons in superconductors, spin waves in magnets, and phasons in incommensurate charge-density-wave systems.

## ACKNOWLEDGMENTS

This work was supported in part by the Swiss National Science Foundation under Division II. M.K. acknowledges support from the Technical University of Munich, Institute for Advanced Study, funded by the German Excellence Initiative and the European Union FP7 under Grant Agreement No. 291763, the Deutsche Forschungsgemeinschaft (DFG) (German Research Foundation) under Germany’s Excellence Strategy No. EXC-2111–390814868, the European Research Council under the European Union’s Horizon 2020 research and innovation program (Grant Agreement No. 851161), and DFG Grants No. KN1254/1-1, No. KN1254/1-2, and No. DFG TRR80 (Project F8). E.D. acknowledges support from Harvard-MIT CUA, AFOSR-MURI: Photonic Quantum Matter (Award No. FA95501610323), and DARPA DRINQS program (Award No. D18AC00014). R.C. acknowledges hospitality at DQMP, Geneve (CH) and at Harvard-MIT CUA. E.O. acknowledges hospitality at DQMP, Geneve (CH).

## APPENDIX A: IRRELEVANT PERTURBATIONS

We want now evaluate the corrections to the response function (21) in the presence of an irrelevant perturbation $H_p = (g/2\pi) \int dx \int d\tau \cos[2\phi(x, \tau)]$. To do this we apply second-order perturbation theory and thus we need to evaluate the average of the time-ordered product of operators

$$ \Phi(x, \tau) = \langle T_{\tau} \partial_{x}^{2} \phi(x, \tau) \cos[2\phi(1)] \cos[2\phi(2)] \partial_{x}^{2} \phi(0, 0) \rangle, \quad (A1) $$

where we used the compact notation $1 \equiv (x_1, \tau_1)$ and similarly for 2. This average can be estimated from the correlator

$$ \langle T_{\tau} e^{i\lambda\partial_{x}\phi(x,\tau)} e^{i\mu\partial_{x}\phi(0,0)} \cos[2\phi(1)] \cos[2\phi(2)] \rangle, \quad (A2) $$

taking the second derivative with respect to $\lambda$ and $\mu$ in the limit $\lambda, \mu \to 0$ and keeping in mind the identity

$$ \langle T_{\tau} \Pi_j e^{iq_j f(\phi(x_j, \tau_j))} \rangle = \exp \left( - \sum_{i>j} q_i \langle T_{\tau} f(\phi(x_i, \tau_i)) f(\phi(x_j, \tau_j)) \rangle \right), \quad (A3) $$
where $f(\phi(x, \tau)) = \partial_x\phi(x, \tau)$ or $\phi(x, \tau)$. This leads to the expression

$$ \Phi(x, \tau) = 64 \int d1 \int d2 e^{-2\langle T_{\tau} [\phi(1) - \phi(2)]^2 \rangle} \langle T_{\tau} \partial_x \phi(x, \tau) \partial_x \phi(0, 0) \rangle \langle T_{\tau} \partial_x \phi(x, \tau) \partial_x \phi(0, 0) \rangle \\ \times [\langle T_{\tau} \partial_x \phi(x, \tau) \phi(1) \rangle \langle T_{\tau} \phi(1) \partial_x \phi(0) \rangle - \langle T_{\tau} \partial_x \phi(x, \tau) \phi(1) \rangle \langle T_{\tau} \phi(2) \partial_x \phi(0) \rangle]. \quad (A4) $$

This expression correspond to Hartree and Fock diagrams. Explicitly one has

$$ e^{-2\langle T_{\tau} [\phi(1) - \phi(2)]^2 \rangle} = \left( \frac{\alpha^2}{(\alpha + u|\tau_1 - \tau_2|)^2 + (x_1 - x_2)^2} \right)^K, \quad (A5) $$

while

$$ \langle T_{\tau} \partial_x \phi(x, \tau) \phi(x_1, \tau_1) \rangle = \partial_x G(x - x_1, \tau - \tau_1) = -\frac{K}{2} \frac{(x - x_1)}{(x - x_1)^2 + (u|\tau - \tau_1| + \alpha)^2}, \quad (A6) $$

$$ \langle T_{\tau} \partial_x \phi(0, 0) \phi(x_1, \tau_1) \rangle = \partial_x G(x_1, \tau_1) = \frac{K}{2} \frac{x_1}{x_1^2 + (u|\tau_1| + \alpha)^2}, $$

and similarly for the term with $1 \to 2$. In (A4) one can factorize the term

$$ \int dx \int d\tau \left( \frac{\alpha^2}{(\alpha + u|\tau|)^2 + x^2} \right)^K = \int d\tau \frac{\alpha^{2K}}{(\alpha + u|\tau|)^{2K-1}} \int \frac{dy}{(1+y^2)^K} = \frac{2\alpha^2}{u} I(K), \quad (A7) $$

which for $K > 1$ gives no power-law correction, and recognize the convolution integral

$$ \int d\tau_1 \int dx_1 \partial_x G(x - x_1, \tau - \tau_1) \partial_x G(x_1, \tau_1) = \int \frac{dq}{2\pi} \int \frac{d\omega}{2\pi} q^2 |G(q, \omega)|^2 e^{i(qx - \omega\tau)}. \quad (A8) $$

Thus the Hartree correction of $\Phi(x, \tau)$ reduces to

$$ \begin{aligned} \Phi_H(x, \tau) &= \frac{\alpha^2}{u} I(K) \int \frac{dq}{2\pi} \int \frac{d\omega}{2\pi} q^2 |G(q, \omega)|^2 e^{i(qx - \omega\tau)} \partial_x^2 G(x, \tau) \\ &= \frac{\alpha^2}{u} I(K) \int \frac{dq}{2\pi} \int \frac{d\nu}{2\pi} q^4 G(q, \nu) G(-q, -\nu) G(q, \omega - \nu), \end{aligned} \quad (A9) $$

where $G(q, \nu) = \frac{uK}{\nu^2+(uq)^2}$. One can easily show that this integral does not increase with $\omega$, so the Hartree correction can be neglected.

The Fock correction instead reduces to the integral

$$ \int d1 \int d2 \left( \frac{\alpha^2}{(\alpha + u|\tau_1 - \tau_2|)^2 + (x_1 - x_2)^2} \right)^K \partial_x G(x - x_1, \tau - \tau_1) \partial_x G(x_2, \tau_2) \partial_x^2 G(x, \tau). \quad (A10) $$

Turning to the Fourier transform representation, we find

$$ \int dx \int d\tau \frac{e^{-iqx+i\omega\tau}}{[x^2+(u|\tau|+\alpha)^2]} = \left[ \sqrt{\left(q^2+\frac{\omega^2}{u^2}\right)\alpha^2} \right]^{K-1} \frac{1}{2^{K-1}\Gamma(K+1)} K^{K-1} \left[ \sqrt{\left(q^2+\frac{\omega^2}{u^2}\right)\alpha^2} \right], \quad (A11) $$

where $K$ is the Bessel function of the second kind. It can be expanded for small $q, \omega$, and $K-1$ a noninteger: While the zeroth-order term is $q$ and $\omega$ independent, the first nonanalytic correction will be of the order $(q^2 + \frac{\omega^2}{u^2})^{K-1}$, which is subdominant compared to $q^2 + \frac{\omega^2}{u^2}$ when $K > 1$. In the integral (A11) we recognize (up to a factor $g^2$) the self-energy correction from a diagrammatic point of view and thus we can write that qualitatively

$$ \Sigma(q, \omega) \approx \Sigma(0, 0) + C \left( q^2 + \frac{\omega^2}{u^2} \right)^{K-1}, \quad (A12) $$

neglecting holomorphic terms of order $(q^2 + \frac{\omega^2}{u^2})\alpha^2$ and higher. We can thus evaluate the Fock correction (A10), which is

$$ \begin{aligned} & \int \frac{dv}{2\pi} \int \frac{dq}{2\pi} q^4 [\Sigma(q, \nu) G^2(q, \nu) G(-q, \omega - \nu) + \Sigma(-q, \omega - \nu) G^2(-q, \omega - \nu) G(q, \nu)] \\ &= g^2 \int \frac{dv}{2\pi} \int \frac{dq}{2\pi} q^4 \left(q^2 + \frac{\nu^2}{u^2} \alpha^2\right)^{K-1} \left(q^2 + \frac{\nu^2}{u^2}\right)^{-2} \left(q^2 + \frac{(\omega - \nu)^2}{u^2}\right)^{-1}. \end{aligned} \quad (A13) $$
<figure><img src="image_5.png" /><figcaption>FIG. 3. Luttinger liquid parameters. We evaluate (a) the Luttinger parameter K, (b) the Luttinger velocity u, and (c) the derivative of the Luttinger parameter with respect to the kinetic energy $dK/dJ_0$ as a function of the inverse system size L for various values of the interaction strength U and fixed density $\rho_0 = 1.2$. Lines (dashed and solid) are guides to the eyes.</figcaption></figure>

By simple power counting this integral behaves as $\omega^{2(K-1)}$ and for $K > 2$ it is subdominant compared to the term $\omega^2$ as $\omega \to 0$. So when $K > 2$ and the umklapp scattering is irrelevant, the intensity of the modulation spectroscopy behaves as

$$ \mathrm{Im}\chi(\omega) \sim \omega^2 + g^2\omega^{2(K-1)} + O(\omega^{2K-2}). \qquad (A14) $$

## APPENDIX B: LUTTINGER PARAMETERS OF THE BOSE-HUBBARD MODEL

We evaluate the Luttinger parameters for the Bose-Hubbard model at filling $\rho = 1.2$ using matrix product states. The Luttinger parameter K is obtained from the density-density correlation function and the Luttinger velocity u from the compressibility $\kappa^{-1} = \frac{\partial^2 E}{\partial N^2} = \frac{K}{\pi u}$. The Luttinger parameters are shown as a function of the inverse system size in Fig. 3. From the Luttinger parameters we can evaluate the prefactor of the absorbed power, as stated in Eq. (22) (see Fig. 4). We compare the analytically predicted prefactor obtained from Luttinger liquid theory (closed circles) with the prefactor of the absorbed power in the time-evolved many-body state (stars), which is on the same order of magnitude.

## APPENDIX C: EVALUATION OF THE RETARDED CORRELATION FUNCTION

For the sake of definiteness, we present the calculation in the case of bosons. The fermionic case proceeds along the

<figure><img src="image_6.png" /><figcaption>FIG. 4. Prefactor of the absorbed power. The analytically predicted prefactor $F(K)$ obtained in Eq. (22) (closed circles) is compared to the numerically obtained prefactor from the full time evolution (stars). The data are shown for various values of the interaction strength U and for a fixed density $\rho_0 = 1.2$ as a function of the inverse system size $1/L$.</figcaption></figure>

same line, with a simple change of prefactor. Using Wick's theorem, the correlator in Eq. (21) is rewritten

$$ \begin{aligned} \langle T_{\tau}(\partial_x \phi)^2(x, \tau) (\partial_x \phi)^2(0, 0) \rangle &= 2[\langle T_{\tau} \partial_x \phi(x, \tau) \partial_x \phi(0, 0) \rangle]^2 = \frac{K^2}{2} \frac{[(u|\tau| + \alpha)^2 - x^2]^2}{[(u|\tau| + \alpha)^2 + x^2]^4}, \end{aligned} \qquad (C1) $$

leading to

$$ \frac{1}{L} \chi_M(\omega) = \frac{\pi}{8} \left( \frac{u}{\pi K} \frac{\partial K}{\partial J_0} \right)^2 \int_{-\infty}^{\infty} \frac{e^{i\omega\tau}}{(u|\tau| + \alpha)^3} d\tau. \qquad (C2) $$

We need the integral

$$ I(\omega) = \int_{-\infty}^{\infty} \frac{e^{i\omega\tau}}{(u\tau + \alpha)^3} d\tau \qquad (C3) $$

and its analytic continuation. We write

$$ \frac{1}{(u|\tau| + \alpha)^3} = \frac{1}{2} \int_0^\infty k^2 e^{-k(u|\tau| + \alpha)} dk \qquad (C4) $$

and obtain

$$ I(\omega) = \frac{1}{2} \int_{0}^{+\infty} dk k^2 e^{-k\alpha} \left[ \frac{1}{i\omega + uk} - \frac{1}{i\omega - uk} \right]. \qquad (C5) $$

We find the analytic continuation $i\omega \to \omega + i0_+$ of Eq. (C5) using the identity

$$ \frac{1}{x + i0_+} = P\left(\frac{1}{x}\right) - i\pi\delta(x), \qquad (C6) $$

which gives

$$ \begin{aligned} \mathrm{Im}\chi(i\omega \to \omega + i0_+) &= \frac{L}{16} \left(\frac{u}{K} \frac{\partial K}{\partial J_0}\right)^2 \int_0^{+\infty} dk k^2 \\ & \times e^{-k\alpha} [\delta(\omega - uk) - \delta(\omega + uk)] \end{aligned} \qquad (C7) $$

$$ = \frac{L}{16u} \left( \frac{1}{K} \frac{\partial K}{\partial J_0} \right)^2 \omega^2 e^{-|\omega|\alpha/u} \mathrm{sgn}(\omega). \qquad (C8) $$

This leads to Eq. (19).
## APPENDIX D: CALCULATION OF THE RESPONSE FUNCTION IN THE CASE OF A SYSTEM WITH BOUNDARIES

In the case of a system with boundaries described by the Hamiltonian (23) with the operator $O_b$ given by (24), we first rewrite

$$O_b = \frac{\partial}{\partial J_0}(uK)\frac{H}{uK} - \int \frac{dx}{\pi} \frac{u}{K^2} \frac{\partial K}{\partial J_0}(\partial_x \phi)^2 \\ + \frac{V}{\pi} \frac{\partial}{\partial J_0} \left[ \ln \left( \frac{uK}{V} \right) \right] [\partial_x \phi(0) + \partial_x \phi(L)] \quad (D1)$$

and as before we only have to calculate the response function of the bulk term proportional to $(\partial_x \phi)^2$ and the edge term proportional to $\partial_x \phi(0) + \partial_x \phi(L)$. To perform the calculation, we first rescale the fields $\phi = \sqrt{K}\tilde{\phi}$ and $\Pi = \tilde{\Pi}/\sqrt{K}$ and introduce the Fourier decomposition (E7) and (E8) to rewrite the Hamiltonian (26) in terms of shifted harmonic oscillators

$$H = \frac{u}{2} \sum_{n=1}^{+\infty} \left[ \pi \tilde{\Pi}_n^2 + \left( \frac{\pi n}{L} \right)^2 \frac{\tilde{\phi}_n^2}{\pi} \right] \\ - \frac{V}{\pi} \sqrt{\frac{2K}{L}} \sum_{n=1}^{+\infty} [1 + (-1)^n] \frac{\pi n}{L} \tilde{\phi}_n \quad (D2)$$

and the operator $O_b$, without the contribution proportional to the Hamiltonian,

$$O_b = -\frac{u}{\pi K} \frac{\partial K}{\partial J_0} \sum_n \left(\frac{\pi n}{L}\right)^2 \tilde{\phi}_n^2 \\ + \frac{V}{\pi} \frac{\partial}{\partial J_0} \left[ \ln \left( \frac{uK}{V} \right) \right] \sqrt{\frac{2K}{L}} \sum_{n=1}^{+\infty} [1 + (-1)^n] \frac{\pi n}{L} \tilde{\phi}_n. \quad (D3)$$

We now introduce $\tilde{\phi}_n$ such that

$$\tilde{\phi}_n = \tilde{\phi}_n - \frac{L}{\pi n} \sqrt{\frac{2K}{L}} [1 + (-1)^n] \frac{V}{u} \quad (D4)$$

to have a Hamiltonian purely quadratic in $\tilde{\phi}_n$. In terms of the new operators,

$$O_b = -\frac{u}{\pi K} \frac{\partial K}{\partial J_0} \sum_n \left(\frac{\pi n}{L}\right)^2 \tilde{\phi}_n^2 \\ + \frac{V}{\pi} \frac{\partial}{\partial J_0} \left[ \ln \left( \frac{uK}{V} \right) \right] \sqrt{\frac{2K}{L}} \sum_{n=1}^{+\infty} [1 + (-1)^n] \frac{\pi n}{L} \tilde{\phi}_n. \quad (D5)$$

The first line gives the contribution calculated in Appendix C. The second gives the contribution coming from the edge potential. The necessary Matsubara correlator is

$$\left\{ \frac{V}{\pi} \frac{\partial}{\partial J_0} \left[ \ln \left( \frac{uK}{V} \right) \right] \right\}^2 \frac{2K}{L} \sum_{n=1}^{+\infty} \left( \frac{\pi n}{L} \right)^2 [1 + (-1)^n]^2 \\ \times \langle T_{\tau} \tilde{\phi}_n(\tau) \tilde{\phi}_n(0) \rangle. \quad (D6)$$

After taking the Fourier transform and making the analytic continuation, we find

$$\mathrm{Im} \chi_{\mathrm{edge}}(\omega + i0_+) = \left\{ \frac{2V}{\pi} \frac{\partial}{\partial J_0} \left[ \ln \left( \frac{uK}{V} \right) \right] \right\}^2 \frac{\pi K}{L} \sum_{n=1}^{+\infty} \sum_{n=1}^{+\infty} \frac{2n\pi}{L} \pi \\ \times \left[ \delta\left(\omega - \frac{2\pi nu}{L}\right) - \delta\left(\omega + \frac{2\pi nu}{L}\right) \right]. \quad (D7)$$

In the limit of $L \to +\infty$, we end up with

$$\mathrm{Im} \chi_{\mathrm{edge}}(\omega + i0_+) = \left\{ \frac{2V}{u} \frac{\partial}{\partial J_0} \left[ \ln \left( \frac{u}{V K} \right) \right] \right\}^2 \frac{\omega K}{2\pi}, \quad (D8)$$

yielding the edge contribution (25), to be added to the bulk contribution.

## APPENDIX E: FRIEDEL OSCILLATIONS

We consider a Bose-Hubbard chain of $M$ sites with open boundary conditions. Its Hamiltonian is

$$H = -J_0 \sum_{l=1}^{M-1} (b_l^\dagger b_{l+1} + b_{l+1}^\dagger b_l) + \frac{U}{2} \sum_{l=1}^{M} n_l(n_l - 1). \quad (E1)$$

We introduce the fictitious sites 0 and $M+1$ to write

$$H = -J_0 \sum_{l=0}^{M} (b_l^\dagger b_{l+1} + b_{l+1}^\dagger b_l) + U \sum_{j=1}^{M} n_l(n_l - 1), \quad (E2)$$

and $b_0 = b_{M+1} = 0$. The bosonized Hamiltonian reads

$$H = \int_0^L \frac{dx}{2\pi} \left[ uK(\pi\Pi)^2 + \frac{u}{K}(\partial_x\phi)^2 \right] \\ - \frac{V_1}{\pi} \partial_x\phi(0) - \frac{V_2}{\pi} \partial_x\phi(L), \quad (E3)$$

with $L = (M+1)a$, and we have included some forward-scattering potentials $V_1$ and $V_2$ at the edges. Our original boson Hamiltonian is symmetric under the reflection $b_l \to b_{M+1-l}$. Using the bosonized expressions of the boson annihilation operators [17,25], we find that under reflection

$$P\phi(x)P^\dagger = -\phi(L-x) - \pi N, \quad (E4)$$

$$P\Pi(x)P^\dagger = -\Pi(L-x), \quad (E5)$$

so that $V_1 = V_2$. The boundary conditions are derived in the fermion case from consideration of the noninteracting limit [45]. In the boson case, we have to consider the expression of the density

$$\rho(x) = \rho_0 - \frac{1}{\pi} \partial_x \phi + A \cos(2\phi(x) - 2\pi \rho_0 x), \quad (E6)$$

which implies through the continuity equation that $j = \partial_t \phi/\pi$. Since no current can leak through the edges of the system, we must have $\partial_t \phi(0) = \partial_t \phi(L) = 0$. So we must impose the Dirichlet boundary conditions $\phi(0, t) = \phi_0$ and $\phi(L, t) = \phi_1$. Moreover, since the number of particles in the system is an integer, by integrating (E6) we find that $(\phi_1 - \phi_0)\pi$ is an integer. We can choose, for instance, $\phi_0 = 0$ and $\phi_1 = -\pi N$, where $N$ is the number of particles added to the initial number of particles in the ground state $N_{GS}$ with $\rho_0 = N_{GS}/L$. We note
that $\partial_x \phi$ can still be nonvanishing as an operator, so we can a priori have edge scattering potentials $V_1$ and $V_2$ in (E3).

Now we introduce the Fourier decomposition

$$ \phi(x) = -\frac{\pi N x}{L} + \sum_{n=1}^{+\infty} \sqrt{\frac{2}{L}} \sin\left(\frac{n\pi x}{L}\right) \phi_n e^{-n\epsilon/2}, \quad (E7) $$

$$ \Pi(x) = \sum_{n=1}^{+\infty} \sqrt{\frac{2}{L}} \sin\left(\frac{n\pi x}{L}\right) \Pi_n e^{-n\epsilon/2}, \quad (E8) $$

which allows us to rewrite

$$ H = \frac{uL}{2\pi K} \left(\frac{\pi N}{L}\right)^2 + \frac{\pi N}{L} \frac{V_1 + V_2}{\pi} \\ + \frac{u}{2} \sum_{n=1}^{+\infty} \pi K \Pi_n^2 + \frac{1}{\pi K} \left(\frac{\pi n}{L}\right)^2 \phi_n^2 \\ - \sum_{n=1}^{+\infty} \sqrt{\frac{2}{L}} \frac{n}{L} [V_1 + (-)^n V_2] \phi_n. \quad (E9) $$

Until now, we have made no assumption concerning the symmetry of our bosonized Hamiltonian under parity. Using the Fourier expansion (E7), we can show that under a parity transformation $P\phi_n P^\dagger = (-1)^n \phi_n$. In the Hamiltonian (E9) $V_1$ and $V_2$ are exchanged by the parity transformation. So we recover $V_1 = V_2$ for a parity invariant Hamiltonian. To find the ground state, we have to minimize the first line with respect to $N$ and determine the shift of oscillators imposed by the edge potentials. The minimization with respect to $N$ yields

$$ N = E\left(\frac{1}{2} - \frac{K}{\pi u}(V_1 + V_2)\right), \quad (E10) $$

while the shift of oscillators is

$$ \langle \phi_n \rangle = \frac{K}{\pi u} \sqrt{2L} \frac{V_1 + (-)^n V_2}{n}. \quad (E11) $$

The expectation value of $\phi(x)$ in the ground state is then

$$ \langle \phi(x) \rangle = -\frac{\pi x}{L} E\left(\frac{1}{2} - \frac{K}{\pi u}(V_1 + V_2)\right) \\ + \frac{2K}{\pi u} \sum_{n=1}^{+\infty} \frac{V_1 + (-)^n V_2}{n} \sin\left(\frac{n\pi x}{L}\right) e^{-n\epsilon/2}. \quad (E12) $$

We thus have

$$ \langle \phi(x) \rangle = \frac{2K}{\pi u} \left[ V_1 \arctan \left( \frac{\sin\left(\frac{\pi x}{L}\right)}{e^{\epsilon/2} - \cos\left(\frac{\pi x}{L}\right)} \right) \right. \\ \left. - V_2 \arctan \left( \frac{\sin\left(\frac{\pi x}{L}\right)}{e^{\epsilon/2} + \cos\left(\frac{\pi x}{L}\right)} \right) \right] \\ - \frac{\pi x}{L} E\left(\frac{1}{2} - \frac{K}{\pi u}(V_1 + V_2)\right). \quad (E13) $$

Taking the limit of $\epsilon \to 0$, for $x$ far enough from an edge, we find the simplified expression

$$ \langle \phi(x) \rangle = \frac{KV_1}{u} - \frac{K(V_1 + V_2)}{uL} \frac{x}{L} - \frac{\pi x}{L} E\left(\frac{1}{2} - \frac{K}{\pi u}(V_1 + V_2)\right), \quad (E14) $$

which is a periodic function of $K(V_1 + V_2)/\pi u$ of period 1. So we can restrict ourselves to $|K(V_1 + V_2)/(\pi u)| < 1/2$ and drop the integer part in Eq. (E14). Using Luttinger liquid theory and the expression (E6), we derive

$$ \langle \rho(x) \rangle = \frac{1}{L} \left[ N_{\text{tot}} + \frac{K(V_1 + V_2)}{\pi u} \right] + A \cos \left[ \frac{2KV_1}{u} \right. \\ \left. - \frac{2\pi x}{L} \left( N_{\text{tot}} + \frac{K(V_1 + V_2)}{\pi u} \right) \right] \left( \frac{\pi a}{L \sin\left(\frac{\pi x}{L}\right)} \right)^K, \quad (E15) $$

and we see that far from the edges, the Friedel oscillations behaves as if the number of particles was $N'_{tot} = N_{tot} + \frac{K(V_1+V_2)}{\pi u}$. The expression (E15) applies only when $\alpha \ll x$ and $\alpha \ll L - x$. It corresponds to the effective Dirichlet boundary conditions $\phi(0) = \frac{KV_1}{u}$ and $\phi(L) = -\frac{KV_2}{u}$ that result from the phase shift on $\phi(x)$ imposed by the edge potentials. When $0 < x < \alpha = \epsilon L$, we cannot take the limit $\epsilon \to 0$ in Eq. (E13). There $\langle \phi(x) \rangle = O(x/\alpha) \to 0$, ensuring that the original Dirichlet boundary conditions are satisfied.

## APPENDIX F: LUTHER-EMERY LIMIT

Let us consider the case of the Fermi-Hubbard model in the Mott insulating phase. When looking at the bosonized expression of the operator $O_f$ Eq. (12), we can rescale the fields $\Pi_\rho \to \Pi_\rho/\sqrt{K_\rho}$ and $\phi_\rho \to \sqrt{K_\rho}\phi_\rho$ such that the operator $O_{f,\rho}$ becomes

$$ O_{f,\rho} = 2a \sin(k_F a) \int \frac{dx}{2\pi} \left[ \frac{(\pi \Pi_\rho)^2}{K_\rho} + K_\rho (\partial_x \phi_\rho)^2 \right]. \quad (F1) $$

The bosonized Hamiltonian for the fermions can also be written in terms of the rescaled fields and at the Luther-Emery point [51] $K_\rho = \frac{1}{2}$ it becomes

$$ H_\rho = \int \frac{dx}{2\pi} u_\rho [(\pi \Pi_\rho)^2 + (\partial_x \phi_\rho)^2] - \frac{2g_3}{(2\pi\alpha)^2} \int dx \cos 2\phi_\rho. \quad (F2) $$

That Hamiltonian is rewritten by introducing the pseudofermions

$$ \Psi_R(x) = \frac{e^{i[\theta(x) - \phi(x)]}}{\sqrt{2\pi\alpha}}, \quad (F3) $$

$$ \Psi_L(x) = \frac{e^{i[\theta(x) + \phi(x)]}}{\sqrt{2\pi\alpha}} \quad (F4) $$

in the form of a gapped noninteracting Hamiltonian

$$ H_\rho = -iu_\rho \int dx (\Psi_R^\dagger \partial_x \Psi_R - \Psi_L \partial_x \Psi_L) - \Delta (\Psi_R^\dagger \Psi_L + \Psi_L^\dagger \Psi_R), \quad (F5) $$

where $\Delta = \frac{2g_3}{4\pi\alpha}$. In terms of the pseudofermions we can rewrite the operator $O_\rho$ as

$$ O_\rho = O_1 + O_2, \quad (F6) $$

$$ O_1 = -\frac{5a}{2} i \sin(k_F a) \int dx (\Psi_R^\dagger \partial_x \Psi_R - \Psi_L \partial_x \Psi_L), \quad (F7) $$

$$ O_2 = 3\pi a \sin(k_F a) \int dx \rho_R \rho_L, \quad (F8) $$
where $\rho_{L,R} = \Psi_{L,R}^{\dagger}\Psi_{L,R}$ is the density operator of the L, R fermions. We have to evaluate the Matsubara correlator

$$ \chi(\tau) = \sum_{i=1,2} \sum_{j=1,2} \chi_{ij}(\tau), \qquad (F9) $$

with $\chi_{ij}(\tau) = \langle T_{\tau} O_i(\tau) O_j(0) \rangle$. This correlator can be expressed in terms of the creation and annihilation operators through the representation $\Psi_v = \frac{1}{\sqrt{L}} \sum_k e^{ikx} c_{k,v}$, in terms of which the Hamiltonian (F5) is written as

$$ H_{\rho} = \sum_{k} u_{\rho} k (c_{kR}^{\dagger} c_{kR} - c_{kL}^{\dagger} c_{kL}) - \Delta (c_{kR}^{\dagger} c_{kL} + \text{H.c.}). \quad (F10) $$

This Hamiltonian can be diagonalized by standard Bogoliubov transformations and expressed in the

form

$$ H'_{\rho} = \sum_{k} E_{k} (c_{k+}^{\dagger} c_{k+} - c_{k-}^{\dagger} c_{k-}), \qquad (F11) $$

with $E_k = \sqrt{(u_\rho k)^2 + \Delta^2}$, $c_{kR} = \cos \varphi_k c_{k+} - \sin \varphi_k c_{k-}$, and $c_{kL} = \sin \varphi_k c_{k+} + \cos \varphi_k c_{k-}$. Then the calculations of the correlators (F9) proceeds by applying Wick's theorem once the single-particle Green's function are known:

$$ \begin{aligned} \langle T_{\tau} c_{kR(L)}(\tau) c_{kR(L)}^{\dagger}(0) \rangle &= \frac{1}{2} \left[ \operatorname{sgn}(\tau) \pm \frac{u_{\rho} k}{E(k)} \right] e^{-|\tau|E(k)}, \\ \langle T_{\tau} c_{kR(L)}(\tau) c_{kL(R)}^{\dagger}(0) \rangle &= \frac{\Delta}{2E(k)} e^{-|\tau|E(k)}. \end{aligned} \qquad (F12) $$

The results for the correlators are

$$ \begin{aligned} \chi_{11}(\tau) &= \frac{25}{4} a^2 \sin^2(k_F a) \Delta^2 \sum_k \frac{k^2}{E(k)^2} e^{-2|\tau|E(k)}, \\ \chi_{22}(\tau) &= \frac{9\pi^2 a^2 \sin^2(k_F a)^2}{L^2} \sum_{k_1, \dots, k_4} \delta_{k_1+k_2, k_3+k_4} \exp \left( -|\tau| \sum_{j=1}^4 E(k_j) \right) \\ & \quad \times \frac{1}{16} \left[ \left( \operatorname{sgn}(\tau) + \frac{u_{k_4}}{E(k_4)} \right) \left( \operatorname{sgn}(\tau) + \frac{u_{k_3}}{E(k_3)} \right) - \frac{\Delta^2}{E(k_3)E(k_4)} \right] \\ & \quad \times \left[ \left( \operatorname{sgn}(-\tau) + \frac{u_{k_1}}{E(k_1)} \right) \left( \operatorname{sgn}(-\tau) + \frac{u_{k_2}}{E(k_2)} \right) - \frac{\Delta^2}{E(k_1)E(k_2)} \right], \\ \chi_{12}(\tau) &= \chi_{21}(\tau) = 0. \end{aligned} \qquad (F13) $$

The correlator $\chi_{22}$ can be simplified close to the threshold where an expansion up to order $O(k^2)$ can be performed such that $\chi_{22}(\tau) \sim (\frac{\Delta}{|\tau|})^{7/2} e^{-4|\tau|\Delta}$. In the complex frequency plane the correlator $\chi_{ij}(i\omega) = \int_{-\infty}^{\infty} d\tau \langle T_{\tau} O_i(\tau) O_j(0) \rangle$ can be analytically extended to evaluate the imaginary part. The result of the calculation gives

$$ \operatorname{Im}\chi_{11}(i\omega \rightarrow \omega + i0) = \frac{25a^2 \sin(k_F a)^2}{2} L \frac{\Delta^2}{\omega u^3} \sqrt{\frac{\omega^2}{4} - \Delta^2}, \qquad (F14) $$

which shows a threshold at $\omega = 2\Delta$ while $\operatorname{Im}\chi_{22}$ has a threshold at $\omega = 4\Delta$. The resulting absorbed power in the Luther-Emery limit is

$$ P_{\text{LE}} = \frac{L \Delta^2}{u} \left( \frac{5a\delta J}{2u} \sin(k_F a) \right)^2 \sqrt{\left( \frac{\omega}{2} \right)^2 - \Delta^2} \quad (\text{F15}) $$

for $2\Delta < \omega < 4\Delta$. This analysis can be extended away from the Luther-Emery point to any value of $K = K_\rho$ by using the form-factor expansion for the sine-Gordon model [52–55], as detailed in Appendix G (see also [56]).

## APPENDIX G: FORM-FACTOR APPROACH

In the present Appendix we want to extend the results derived using the Luther-Emery limit to any value of $K_\rho$. For $K_\rho > 1/2$ the excitations are massive solitons and antisolitons of mass $M_s$, while for $K < 1/2$ we also have breathers of mass $M_n = 2M_s \sin(n\frac{\pi}{2}\frac{K}{1-K})$ with $1 \le n < \frac{1}{K_\rho} - 1$ an integer. Working in the vicinity of the Luther-Emery point, the low-energy Hamiltonian is

$$ \begin{aligned} H = & -i \frac{u_{\rho}}{2} \left( 2K_{\rho} + \frac{1}{2K_{\rho}} \right) \int dx (\Psi_{R}^{\dagger} \partial_{x} \Psi_{R} - \Psi_{L} \partial_{x} \Psi_{L}) \\ & - \frac{g}{2\pi\alpha} (\Psi_{R}^{\dagger} \Psi_{L} + \Psi_{L}^{\dagger} \Psi_{R}) \\ & - \pi u_{\rho} \left( 2K_{\rho} - \frac{1}{2K_{\rho}} \right) \int dx \hat{\rho}_{R} \hat{\rho}_{L}, \end{aligned} \qquad (G1) $$

where the operator $\hat{\rho}_a = \Psi_a^\dagger \Psi_a$, $a = L, R$. It is the Hamiltonian of a massive Thirring model [73],

$$ \begin{aligned} H = & -i\bar{v} \int dx (\Psi_R^\dagger \partial_x \Psi_R - \Psi_L \partial_x \Psi_L) \\ & - M \int dx (\Psi_R^\dagger \Psi_L + \Psi_L^\dagger \Psi_R) - \bar{g} \int dx \hat{\rho}_R \hat{\rho}_L, \end{aligned} \qquad (G2) $$

where $\bar{v} = \frac{u}{2}(2K + \frac{1}{2K})$, $M = \frac{g}{2\pi\alpha}$, and $\bar{g} = \pi u_\rho(2K - \frac{1}{2K})$.
The kinetic energy operator is related to the component $T^{11}$ of the momentum-energy tensor (see [74]), so

$$ \frac{u_{\rho}}{2a \sin(k_{F}a)} \mathcal{O} = \frac{2K_{\rho}}{2K_{\rho}^{2} + 1} \int dx T^{11}(x) + 2\pi u_{\rho} \frac{K_{\rho}^{2} - 1}{K_{\rho}^{2} + 1} \int dx \hat{\rho}_{R} \hat{\rho}_{L}. \quad (G3) $$

Since $\hat{\rho}_R \hat{\rho}_L = \hat{\psi}_R^\dagger \hat{\psi}_L^\dagger \hat{\psi}_R \hat{\psi}_L$, that operator can only have matrix elements between the ground state of the massive Thirring model and a state containing two solitons and two antisolitons, i.e., a state with energy at least $4M_s$. So that term will not contribute for frequencies $\omega < 4M_s$ and we will have

$$ \mathrm{Im} \chi_{\rho}(\omega) = L \left( \frac{4K_{\rho}a \sin(k_{F}a)}{u_{\rho}(K_{\rho}^{2} + 1)} \right)^{2} \mathrm{Im} \chi^{ss}(\omega), \quad (G4) $$

where the contribution $\mathrm{Im} \chi^{ss}(\omega)$ of the $T^{11}$ component of the momentum-energy tensor is obtained from the form-factor expansion [53–55]. For the lowest excited state formed of a single soliton-antisoliton pair, we have

$$ \mathrm{Im} \chi^{ss}(\omega) = 2\pi \int \frac{d\theta}{(2\pi)^2} |\langle 0 | T^{11} | \theta, \bar{\theta} \rangle_{ss}|^2 \\ \times \delta\left(\frac{M_s}{u_\rho}(\sinh\theta + \sinh\bar{\theta})\right) \\ \times \delta(\omega - M_s(\cosh\theta + \cosh\bar{\theta})), \quad (G5) $$

where, according to [55], the form factor of the energy-momentum tensor is

$$ \langle 0 | T^{11} | \theta_1, \theta_2 \rangle_{ss} = -2i \frac{M^2}{u} \cosh^2 \left( \frac{\theta_1 + \theta_2}{2} \right) \\ \times \sinh \left( \frac{\theta_1 - \theta_2}{2} \right) \frac{F_+(\theta_1 - \theta_2)}{\nu}, \quad (G6) $$

$$ F_+(\theta) = \frac{i \cosh(\theta/2)}{\sinh\left(\frac{i\pi - \theta}{2\nu}\right)} F_{\min}(\theta), \quad (G7) $$

$$ F_{\min}(\theta) = \exp \left[ \int_{0}^{+\infty} \frac{dt}{t} \frac{\sinh \frac{1}{2}(1-\nu)t}{\sinh \frac{\nu t}{2} \cosh \frac{t}{2}} \frac{1 - \cosh \left(1 - \frac{\theta}{i\pi}\right)t}{2 \sinh t} \right]. \quad (G8) $$

Here $\nu = K_\rho/(1-K_\rho)$ for the Fermi-Hubbard model and $\nu = K/(2-K)$ for the Bose-Hubbard model. When $\theta \to 0$, $F_+(\theta) \to 1$, giving, for $\omega \to 2M^+$,

$$ \mathrm{Im} \chi^{ss}(i\omega \to \omega + i0) = \frac{4M_s^3}{\pi \omega u_\rho v^2} \sqrt{\left(\frac{\omega}{2M_s}\right)^2 - 1}. \quad (G9) $$

Now we have a threshold at twice the mass of the soliton. A similar threshold behavior was also obtained in the case of modulation of a weak optical lattice [56]. Technically, this can be understood as follows. We can always subtract an operator proportional to the Hamiltonian (G2) from the operator $\mathcal{O}$. So we would obtain an equivalent result if $T^{11}$ were replaced by a term proportional to $\Psi_R^\dagger \Psi_L + \Psi_L^\dagger \Psi_R$, which is precisely the perturbing term in [56].

[1] I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. 80, 885 (2008).

[2] T. Esslinger, Annu. Rev. Condens. Matter Phys. 1, 129 (2010).

[3] M. Greiner, O. Mandel, T. Esslinger, T. W. Hänsch, and I. Bloch, Nature (London) 415, 39 (2002).

[4] M. Endres, T. Fukuhara, D. Pekker, M. Cheneau, P. Schauss, C. Gross, E. Demler, S. Kuhr, and I. Bloch, Nature (London) 487, 454 (2012).

[5] E. Haller, R. Hart, M. J. Mark, J. G. Danzl, L. Reichsöllner, M. Gustavsson, M. Dalmonte, G. Pupillo, and H.-C. Nägerl, Nature (London) 466, 597 (2010).

[6] G. Boéris, L. Gori, M. D. Hoogerland, A. Kumar, E. Lucioni, L. Tanzi, M. Inguscio, T. Giamarchi, C. D'Errico, G. Carleo, G. Modugno, and L. Sanchez-Palencia, Phys. Rev. A 93, 011601(R) (2016).

[7] T. Stöferle, H. Moritz, C. Schori, M. Köhl, and T. Esslinger, Phys. Rev. Lett. 92, 130403 (2004).

[8] C. Kollath, A. Iucci, T. Giamarchi, W. Hofstetter, and U. Schollwöck, Phys. Rev. Lett. 97, 50402 (2006).

[9] P. Bordia, H. Lüschen, U. Schneider, M. Knap, and I. Bloch, Nat. Phys. 13, 460 (2017).

[10] A. Rubio-Abadal, M. Ippoliti, S. Hollerith, D. Wei, J. Rui, S. L. Sondhi, V. Khemani, C. Gross, and I. Bloch, Phys. Rev. X 10, 021044 (2020).

[11] A. Tokuno and T. Giamarchi, Phys. Rev. Lett. 106, 205301 (2011).

[12] R. Citro, N. Andrei, and Q. Niu, Phys. Rev. B 68, 165312 (2003).

[13] A. Tokuno and T. Giamarchi, Phys. Rev. A 85, 061603(R) (2012).

[14] D. Greif, T. Uehlinger, G. Jotzu, L. Tarruell, and T. Esslinger, Science 340, 1307 (2013).

[15] G. Pieplow, F. Sols, and C. E. Creffield, New J. Phys. 20, 073045 (2018).

[16] C. Schweizer, F. Grusdt, M. Berngruber, L. Barbiero, E. Demler, N. Goldman, I. Bloch, and M. Aidelsburger, Nat. Phys. 15, 1168 (2019).

[17] T. Giamarchi, *Quantum Physics in One Dimension*, International Series of Monographs on Physics Vol. 121 (Clarendon, Oxford, 2004).

[18] M. A. Cazalilla, R. Citro, T. Giamarchi, E. Orignac, and M. Rigol, Rev. Mod. Phys. 83, 1405 (2011).

[19] K. Staliunas, S. Longhi, and G. J. de Valcarcel, Phys. Rev. Lett. 89, 210406 (2002).

[20] K. Staliunas, S. Longhi, and G. J. de Valcarcel, Phys. Rev. A 70, 011601(R) (2004).

[21] P. Engels, C. Atherton, and M. A. Hoefer, Phys. Rev. Lett. 98, 095301 (2007).

[22] R. Nath and L. Santos, Phys. Rev. A 81, 033626 (2010).

[23] A. Georges and T. Giamarchi, in *Many-Body Physics with Ultracold Gases, Proceedings of the Les Houches Summer School of Theoretical Physics, XCIV, 2010*, edited by C. Salomon, G. V. Shlyapnikov, and L. F. Cugliandolo (Oxford University Press, Oxford, 2013).

[24] L. D. Landau and E. M. Lifshitz, *Statistical Physics* (Pergamon, New York, 1959).
[25] F. D. M. Haldane, Phys. Rev. Lett. 47, 1840 (1981).

[26] V. J. Emery, in *Highly Conducting One-Dimensional Solids*, edited by J. T. Devreese, R. P. Evrard, and V. E. van Doren (Plenum, New York, 1979), p. 247.

[27] J. Sólyom, Adv. Phys. 28, 201 (1979).

[28] H. Fukuyama and H. Takayama, in *Electronic Properties of Inorganic Quasi-One-Dimensional Compounds: Part I—Theoretical*, edited by P. Monceau (Reidel, New York, 1985), p. 41.

[29] H. J. Schulz, in *Mesoscopic Quantum Physics, Proceedings of the Les Houches Summer School of Theoretical Physics, LXI, 1995*, edited by E. Akkermans, G. Montambaux, J. L. Pichard, and J. Zinn-Justin (Elsevier, Amsterdam, 1995), p. 533.

[30] J. Voit, Rep. Prog. Phys. 58, 977 (1995).

[31] K. Schönhammer, Luttinger liquids: the basic concepts, in *Strong Interactions in Low Dimensions*, Physics and Chemistry of Materials with Low-Dimens, Vol. 25, edited by D. Baeriswyl and L. Degiorgi (Springer, Dordrecht, 2004), pp. 93–136.

[32] A. O. Gogolin, A. A. Nersesyan, and A. M. Tsvelik, *Bosonization and Strongly Correlated Systems* (Cambridge University Press, Cambridge, 1999).

[33] N. Nagaosa, *Quantum Field Theory in Strongly Correlated Electronic Systems* (Springer, Heidelberg, 1999).

[34] Y. Jompol, C. J. B. Ford, J. P. Griffiths, I. Farrer, G. A. C. Jones, D. Anderson, D. A. Ritchie, T. W. Silk, and A. J. Schofield, *Science* 325, 597 (2009).

[35] G. Salomon, J. Koepsell, J. Vijayan, T. A. Hilker, J. Nespolo, L. Pollet, I. Bloch, and C. Gross, *Nature (London)* 565, 56 (2019).

[36] R. Shankar, Int. J. Mod. Phys. B 4, 2371 (1990).

[37] M. A. Cazalilla, J. Phys. B 37, S1 (2004).

[38] S. Rachel, N. Laflorencie, H. F. Song, and K. Le Hur, Phys. Rev. Lett. 108, 116401 (2012).

[39] A. M. Läuchli, arXiv:1303.0741.

[40] P. Jordan and E. Wigner, Z. Phys. 47, 631 (1928).

[41] T. Giamarchi, Phys. Rev. B 44, 2905 (1991).

[42] F. Verstraete, V. Murg, and J. I. Cirac, Adv. Phys. 57, 143 (2008).

[43] U. Schollwöck, Ann. Phys. (NY) 326, 96 (2011).

[44] S. Eggert and I. Affleck, Phys. Rev. B 46, 10866 (1992).

[45] M. Fabrizio and A. O. Gogolin, Phys. Rev. B 51, 17827 (1995).

[46] V. Brunel, M. Bocquet, and T. Jolicoeur, Phys. Rev. Lett. 83, 2821 (1999).

[47] I. Affleck, J. Phys. A: Math. Gen. 31, 2761 (1998).

[48] J. Friedel, Philos. Mag. 43, 153 (1952).

[49] R. Egger and H. Grabert, Phys. Rev. Lett. 75, 3505 (1995).

[50] S. Rommer and S. Eggert, Phys. Rev. B 62, 4370 (2000).

[51] A. Luther and V. J. Emery, Phys. Rev. Lett. 33, 589 (1974).

[52] F. H. L. Essler and R. M. Konik, in *From Fields to Strings: Circumnavigating Theoretical Physics: Ian Kogan Memorial Collection*, edited by M. Shifman, A. Vainshtein, and J. Wheater (World Scientific, Singapore, 2004), Vol. 1, Pt. 2, p. 684.

[53] M. Karowski and P. Wiesz, Nucl. Phys. B 139, 455 (1978).

[54] H. Babujian, A. Fring, M. Karowski, and A. Zapletal, Nucl. Phys. B 538, 535 (1999).

[55] H. Babujian and M. Karowski, Nucl. Phys. B 620, 407 (2002).

[56] A. Iucci, M. A. Cazalilla, A. F. Ho, and T. Giamarchi, Phys. Rev. A 73, 041608(R) (2006).

[57] E. Haller, M. J. Mark, R. Hart, J. G. Danzl, L. Reichsöllner, V. Melezhik, P. Schmelcher, and H.-C. Nägerl, Phys. Rev. Lett. 104, 153203 (2010).

[58] C. Kollath, A. Iucci, I. McCulloch, and T. Giamarchi, Phys. Rev. A 74, 041604(R) (2006).

[59] L.-M. Duan, E. Demler, and M. D. Lukin, Phys. Rev. Lett. 91, 090402 (2003).

[60] A. Widera, S. Trotzky, P. Cheinet, S. Fölling, F. Gerbier, I. Bloch, V. Gritsev, M. D. Lukin, and E. Demler, Phys. Rev. Lett. 100, 140401 (2008).

[61] D. M. Stamper-Kurn and M. Ueda, Rev. Mod. Phys. 85, 1191 (2013).

[62] N. Jepsen et al. (unpublished).

[63] P. Lähteenmäki, G. S. Paraoanu, J. Hassel, and P. J. Hakonen, Proc. Natl. Acad. Sci. USA 110, 4234 (2013).

[64] R. Kuzmin, R. Mencia, N. Grabon, N. Mehta, Y. H. Lin, and V. E. Manucharyan, Nat. Phys. 15, 930 (2019).

[65] D. Podolsky, A. Auerbach, and D. P. Arovas, Phys. Rev. B 84, 174522 (2011).

[66] D. Pekker and C. Varma, Annu. Rev. Condens. Matter Phys. 6, 269 (2015).

[67] T. Kampfrath, K. Tanaka, and K. A. Nelson, Nat. Photon. 7, 680 (2013).

[68] D. N. Basov, R. D. Averitt, D. van der Marel, M. Dressel, and K. Haule, Rev. Mod. Phys. 83, 471 (2011).

[69] C. Giannetti, M. Capone, D. Fausti, M. Fabrizio, F. Parmigiani, and D. Mihailovic, Adv. Phys. 65, 58 (2016).

[70] D. N. Basov, R. D. Averitt, and D. Hsieh, Nat. Mater. 16, 1077 (2017).

[71] Y. Cao, D. G. Mazzone, D. Meyers, J. P. Hill, X. Liu, S. Wall, and M. P. M. Dean, Philos. Trans. A 377, 20170480 (2018).

[72] A. von Hoegen, M. Fechner, M. Först, J. Porras, B. Keimer, M. Michael, E. Demler, and A. Cavalleri, arXiv:1911.08284.

[73] H. Thacker, Rev. Mod. Phys. 53, 253 (1982).

[74] C. Itzykson and J. B. Zuber, *Quantum Field Theory* (McGraw-Hill, New York, 1980), p. 143.