# Interaction Induced Anderson Transition in a Kicked One Dimensional Bose Gas

H. Olsen, P. Devillard, G. Aupetit-Diallo, P. Vignolo, M. Albert

► To cite this version:

H. Olsen, P. Devillard, G. Aupetit-Diallo, P. Vignolo, M. Albert. Interaction Induced Anderson Transition in a Kicked One Dimensional Bose Gas. Physical Review Letters, 2025, 135 (17), pp.173403. 10.1103/dyxq-qkmd. hal-05479942

HAL Id: hal-05479942

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

Submitted on 27 Jan 2026

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 CC BY 4.0 - Attribution - International License
Interaction induced Anderson transition in a kicked one dimensional Bose gas

H. Olsen,¹ P. Devillard,² G. Aupetit-Diallo,³ P. Vignolo,¹,⁴ and M. Albert¹,⁴

¹Université Côte d'Azur, CNRS, Institut de Physique de Nice, France

²Aix Marseille Univ, Université de Toulon, CNRS, CPT, Marseille, France

³SISSA, Via Bonomea 265, I-34136 Trieste, Italy

⁴Institut Universitaire de France

We investigate the Lieb-Liniger model of one-dimensional bosons subjected to periodic kicks. In both the non-interacting and strongly interacting limits, the system undergoes dynamical localization, leading to energy saturation at long times. However, for finite interactions, we reveal an interaction-driven transition from an insulating to a metallic phase at a critical kicking strength, provided the number of particles is three or more. Using the Bethe Ansatz solution of the Lieb-Liniger gas, we establish a formal correspondence between its dynamical evolution and an Anderson model in *N* spatial dimensions, where *N* is the number of particles. This theoretical prediction is supported by extensive numerical simulations for three particles, complemented by finite-time scaling analysis, demonstrating that this transition belongs to the orthogonal Anderson universality class.

The phenomenon of Anderson localization, originally introduced in the context of electron transport in disordered media [1], is now recognized as a universal feature of classical and quantum wave dynamics. In non-interacting systems, interference effects lead to the exponential localization of eigenstates and the suppression of diffusion. A key aspect of Anderson localization is the existence of a disorder-induced metal-insulator transition in three dimensions (3D), where a critical level of disorder separates a diffusive metallic phase at small disorder from a localized insulating phase at strong disorder [2]. This transition, known as the Anderson transition, shares key properties with second-order phase transitions and can be described within the framework of the one-parameter scaling theory [3]. According to this theory, the localization length diverges at criticality as $\ell \sim (W - W_c)^{-\nu}$, where *W* is the disorder strength, $W_c$ the critical disorder, and $\nu$ the critical exponent. On the other side of the transition this is the diffusion constant *D* which vanishes as $(W_c - W)^s$. Importantly, Anderson localization is fundamentally dependent on dimensionality: while all states are localized in one dimension (1D) for any amount of disorder, in 3D (and for larger dimensions), a mobility edge separates localized and extended states [3]. A major open question is how localization is affected by interactions, particularly in low-dimensional and driven systems. While interactions generally promote delocalization [4], strong disorder can stabilize a many-body localized (MBL) phase, which exhibits ergodicity breaking and prevents thermalization [5–9].

Dynamical localization, the analog of Anderson localization for driven systems, occurs in periodically kicked quantum systems, where momentum states play the role of spatial sites in the Anderson model [10–12]. The quantum kicked rotor (QKR) is a paradigmatic model in this context: in the absence of interactions, dynamical localization manifests itself as an energy saturation at long times due to destructive quantum interference. The QKR

has been extensively studied in atomic experiments, leading to key insights into Anderson localization [13], including the observation of Anderson transitions [14–18] and the effects of symmetries on weak localization [19, 20].

Ultracold atomic gases provide an ideal platform [21–23] to investigate the interplay between interactions and localization in driven quantum systems. One-dimensional (1D) Bose gases, described by the Lieb-Liniger model [24], are particularly well suited for exploring these effects, as they allow for precise control over interactions and can be routinely realized in experiments [25–28]. Despite recent progress, the fate of dynamical localization in the presence of interactions remains a subject of debate. In the non-interacting and strongly interacting (Tonks-Girardeau) limits [29, 30], dynamical localization persists for any number of particles [31, 32]. For two interacting bosons, localization is also always preserved [33]. However, for finite interactions and more than two particles, theoretical and experimental results are conflicting. Mean-field studies suggest that interactions induce subdiffusive energy growth, destroying dynamical localization [34–39], whereas, Rylands et al. argued that dynamical localization persists for an arbitrary number of particles in the weak-kicking regime [32]. Experimental studies have also provided contradictory conclusions: two experiments found that interactions destroy localization [40, 41], while another observed its robustness [42].

In this letter, we resolve this open problem by employing exact methods. We study the kicked Lieb-Liniger model, describing a 1D Bose gas subjected to periodic kicks, and establish a formal correspondence between its dynamical evolution and an Anderson model in *N* spatial dimensions, where *N* is the number of particles. This mapping enables us to predict the existence of an Anderson-type metal-insulator transition at a critical kicking strength, dependent on interactions for *N* ≥ 3. Our findings are supported by extensive numerical simu-
lations and finite-time scaling analysis for $N = 3$, providing new insights into the fundamental interplay between interactions and localization in driven quantum systems.

**Model.**— We consider $N$ identical bosons of mass $m$ and coordinates $y_i$, living on a ring of size $L$ with point-like repulsive interactions of strength $g$, and periodic boundary conditions. In addition, the atoms are subjected to a kick potential of the form $K \cos(2\pi y_i/L) \sum_n \delta(t' - n\tau)$ with period $\tau$ and amplitude $K$. After rescaling time with $\tau$, positions with $\ell = L/2\pi$ and energy with $\varepsilon = m\ell^2/\tau^2$ the two contributions of the Hamiltonian at which the atoms are subjected read

$$H_0 = \sum_{i=1}^{N} \frac{p_i^2}{2} + c \sum_{i>j} \delta(x_i - x_j), \quad (1)$$

$$H_{\text{kick}} = H_K \sum_{n=-\infty}^{+\infty} \delta(t-n), \quad H_K = K \sum_{i=1}^{N} \cos(x_i), \quad (2)$$

where $c = g/(\varepsilon\ell)$, $K = K/(\tau\varepsilon)$, $t = t'/\tau$, $x_i = y_i/\ell$, $p_i = -i\hbar e \partial/\partial x_i$ the momenta of the particles satisfying the commutation relation $[x_j, p_j] = i\hbar e$, where $\hbar e = \hbar/(\varepsilon\tau)$ is the effective Planck's constant. $H_0$ is the Lieb-Liniger Hamiltonian which can be exactly solved with Bethe Ansatz [24]. The eigenstates are labeled by a set of $N$ rapidities $\lambda_i$ and can be written, in the fundamental sector $x_1 \le x_2 \le \dots \le x_N$, as

$$\Psi_{\{\lambda_j\}}(\{x_j\}) = \sum_{P \in S_N} A_P \exp\left(i \sum_{k=1}^{N} \lambda_{P(k)} x_k\right), \quad (3)$$

where $S_N$ is the permutation group of $N$ elements. The $A_P$ are coefficients given by

$$A_P = \mathcal{N}(-1)^{\operatorname{sgn}(P)} \prod_{1 \le k < j \le N} [\lambda_{P(j)} - \lambda_{P(k)} - ic] \quad (4)$$

with $\mathcal{N}$ the normalization constant. In addition, the rapidities are all real and fulfill the $N$ coupled Bethe equations

$$\lambda_j = I_j - \frac{1}{\pi} \sum_{k=1}^{N} \arctan \left( \frac{\lambda_j - \lambda_k}{c} \right), \quad (5)$$

where the Bethe numbers $I_j$ are integers for odd $N$ and half-integers for even $N$, and must all be distinct. The energy of a Bethe state is given by $E_{\vec{\lambda}} = \frac{\hbar^2}{2} \sum_{i=1}^{N} \lambda_i^2$ and its momentum $P_{\vec{\lambda}} = \hbar e \sum_{i=1}^{N} \lambda_i$. The ground state corresponds to the lowest available set of Bethe numbers. For instance, for $N = 3$, it is given by $\vec{I} = (-1, 0, +1)$. Due to the bosonic statistics, these states are symmetric under permutations of spatial coordinates and rapidities, allowing us to restrict our calculations to the fundamental rapidity sector $I_1 < I_2 < \dots < I_N$. It is important

to note that the rapidities $\lambda_i$ differ from the actual momenta of the particles except in the non-interacting limit. In the limit of infinite interactions ($c \to +\infty$), the rapidities coincide with the momenta of the corresponding free fermions via the Bose-Fermi mapping [30]. In general, however, they do not reside on a regular lattice (which would have spacing $2\pi\ell/L = 1$ in our units), but are instead coupled through the Bethe equations (5).

Before proceeding further, we summarize a few key results relevant to our study. In the non-interacting limit ($c = 0$), the system reduces to $N$ independent quantum kicked rotors, each of which undergoes dynamical localization, characterized by energy saturation at long times [11]. At short times, the momentum-space dynamics of each rotor is diffusive, in agreement with the classical description, with a diffusion constant approximately given by $D \simeq K^2/4$. However, beyond the localization time $\tau_{\text{loc}} \simeq K^2/4\hbar_e^2$, quantum interference effects halt diffusion, leading to energy saturation at a characteristic momentum scale $p_{\text{loc}} \simeq K^2/4\hbar_e$, the localisation momentum [12]. In the opposite limit of infinite interactions ($c \to +\infty$), the system maps onto a gas of free fermions via the Bose-Fermi correspondence [30]. In this regime, the localization properties remain qualitatively similar to the non-interacting case, as previously demonstrated in theoretical and numerical studies [31, 32]. Additionally, the case of two interacting bosons ($N = 2$) has been investigated for arbitrary interaction strength [33], revealing that dynamical localization persists regardless of the interaction strength.

**General mapping to an Anderson-like model.**— Since the Hamiltonian is time-periodic, the system's evolution can be described by iterating the time evolution over one period (here and in the following $\tau = 1$)

$$U = e^{-iH_0/\hbar e} e^{-iH_K/\hbar e}. \quad (6)$$

For $N = 1$ it is well known that the eigenstates of this Floquet operator are dynamically localized in momentum space [11]. This localization can be understood through a formal analogy between the Floquet eigenvalue equation and a one-dimensional Anderson model. A natural question then arises: to what extent do these results generalize to a larger number $N$ of interacting particles? Here, we demonstrate that a similar formal mapping can be established in the $N$-dimensional rapidity space. Instead of relying on the tangent map approach [10, 11], we follow the strategy outlined in Ref. [43]. Our basis of states consists of the Bethe states $|\vec{\lambda}\rangle$, rather than the usual momentum eigenstates. The Floquet eigenstates of the evolution operator (6) satisfy $U|\phi_\alpha\rangle = e^{-i\mu_\alpha/\hbar e}|\phi_\alpha\rangle$, where $\mu_\alpha$ are the Floquet quasi-energies. Moving to the interaction representation by defining $|\tilde{\phi}_\alpha\rangle = e^{iH_0/2\hbar e}|\phi_\alpha\rangle$, we note that this unitary transformation does not alter the state structure in rapidity space, since $H_0$ is diagonal in
this basis. The Floquet equation then takes the form

$$ \sum_{\vec{\mu} \neq \vec{\lambda}} W_{\vec{\lambda}, \vec{\mu}} \langle \vec{\mu} | \tilde{\phi}_{\alpha} \rangle + W_{\vec{\lambda}, \bar{\lambda}} \langle \bar{\lambda} | \tilde{\phi}_{\alpha} \rangle = 2 \cos \left( \frac{\mu_{\alpha}}{\hbar_{e}} \right) \langle \bar{\lambda} | \tilde{\phi}_{\alpha} \rangle, \quad (7) $$

with the matrix elements

$$ W_{\vec{\lambda}, \vec{\mu}} = 2 \operatorname{Re} \left\{ \exp \left( -i \frac{E_{\vec{\lambda}} + E_{\vec{\mu}}}{2\hbar_e} \right) M_{\vec{\lambda}, \vec{\mu}} \right\}, \quad (8) $$

$$ M_{\vec{\lambda}, \vec{\mu}} = \langle \vec{\lambda} | e^{-iH_K / \hbar_e} | \vec{\mu} \rangle, \quad (9) $$

and $\operatorname{Re}\{X\}$ denotes the real part of a complex number $X$. This equation closely resembles an Anderson model with long-range hopping in an $N$-dimensional rapidity space. However, for this analogy to hold rigorously, two key conditions must be met. First, the on-site energy $W_{\vec{\lambda}, \vec{\lambda}}$ must exhibit pseudo-random behavior, which is ensured if $\hbar_e$ is incommensurate with $4\pi$ [10]. Second, the hopping terms must decay faster than $1/||\vec{\lambda} - \vec{\mu}||$ [44]. This requirement is satisfied in our model since the matrix elements $M_{\vec{\lambda}, \vec{\mu}}$ decay exponentially [45].

This formal correspondence leads to a highly nontrivial prediction: an Anderson-type insulator-to-metal transition should emerge for $N \ge 3$. This, in particular, explains why no transition was observed in the two-particle case studied in Ref. [33]. However, several subtleties need to be addressed. Unlike a standard Anderson model describing a particle on a $d$-dimensional lattice, our system prohibits double occupancy of a site, as two rapidities cannot be equal. Additionally, the many-body wave function is invariant under rapidity permutations, allowing us to restrict the problem to the fundamental sector $\lambda_1 < \lambda_2 < \dots < \lambda_N$. Moreover, disorder in our model is genuinely $N$-dimensional only when interactions are finite. In both the non-interacting limit ($c \to 0$) and the Tonks-Girardeau limit ($c \to \infty$), the Bethe equations (5) decouple, meaning that all rapidities evolve independently [45]. In these limits, the system effectively reduces to $N$ independent copies of the one-dimensional kicked rotor, which explains why dynamical localization persists. By contrast, for finite $c$, the rapidities are coupled, and the disorder acquires an $N$-dimensional structure due to the non-factorization of the evolution operator.

Under these conditions, we expect an Anderson transition at a critical stochasticity parameter $K_c$, which depends on both the interaction strength $c$ and the effective Planck constant $\hbar_e$. Furthermore, we anticipate that $K_c$ diverges in both the weakly interacting ($c \to 0$) and strongly interacting ($c \to \infty$) limits, reinforcing the absence of a transition in these regimes.

*Results for $N = 3$ bosons.*— To illustrate our findings, we consider the case of $N = 3$ bosons, where extensive numerical calculations can be performed. To do so, we truncate the Hilbert space by restricting the Bethe numbers to the range $-N_s \le I_j \le N_s$, where $N_s$ is an integer

<figure><img src="image_2.png" /><figcaption>FIG. 1: Upper panel: time evolution of the total energy of the three particles for different values of the stochasticity parameter $K$ at finite interaction strength $c = 10$. Close to the critical point $K_c/\hbar_e \simeq 2.1$ (green curve) anomalous diffusion is observed with exponent $2/3$ (black dashed line). For $K = 1.5\hbar_e < K_c$ (blue curve) the energy growth saturates as expected in the localized regime. Above the critical point $K = 3.2\hbar_e > K_c$ (red curve) the dynamics is diffusive up to some time where the evolution is sensitive to finite size effects. Lower panel: phase diagram of the dynamical phases in the $c-K$ plane. The circles correspond to numerical determination of the critical line while the color code is just a guide to the eyes based on a fit of the data. Here $\hbar_e = 2.89$, $N = 3$ and $N_s = 29$. The white arrow indicates the path we follow to analyse the transition at $c = 10$.</figcaption></figure>

(since $N = 3$ is odd). The computation of the Floquet operator matrix elements (6) then reduces to evaluating $M_{\vec{\lambda}, \vec{\mu}} = \langle \vec{\lambda} | e^{-iH_K / \hbar_e} | \vec{\mu} \rangle$ since $e^{-iH_0 / \hbar_e}$ is diagonal in this basis. These matrix elements can be computed analytically using (3) [45]. The total matrix size in the truncated basis is given by $N_M = \binom{2N_s+1}{3}$. We initialize the system in the ground state, characterized by Bethe numbers $\vec{I} = (-1, 0, 1)$, and evolve it numerically by iterating the Floquet operator in the truncated Hilbert space with $N_s = 29$ for selected values of $K$, $c$, and $\hbar_e$. The total energy is then computed as $\langle E(t) \rangle = \sum_{\vec{\lambda}} |\alpha_{\vec{\lambda}}(t)|^2 E_{\vec{\lambda}}$, where $\alpha_{\vec{\lambda}}(t)$ represents the amplitude of the many-body wave function in the Bethe state $|\vec{\lambda}\rangle$.

This approach successfully reproduces the expected behavior in both the non-interacting and strongly interacting limits. It also enables us to compute the phase diagram shown in Fig. 1, which demonstrates the existence of a metal-insulator phase transition at a critical value of $K$ that depends on the interaction strength and diverges
<figure><img src="image_3.png" /><figcaption>FIG. 2: Finite time scaling applied to the numerical results with $c = 10$ and $N = 3$. The time evolution of $\langle E \rangle$ is computed as a function of time, from 20 to 200 kicks, for several values of $K/\hbar_e$ between $K/\hbar_e = 0.4$ to $K/\hbar_e = 2.9$. The scaling function $\ln \Lambda$ (panel a), with $\Lambda = \langle E \rangle / t^{2/3}$, displays a lower branch (blue) associated with the localized regime and an upper branch (red) associated with the diffusive regime. The continuous curve is a fit using a Taylor expansion of the scaling function up to fourth order and a critical exponent $\nu = 1.56$. The dependence of the scaling parameter $\xi(K)$ is displayed on panel b) and shows a divergent behavior around the critical point $K_c/\hbar_e = 2.13$. The continuous light-blue curve is a fit using a critical exponent $\nu = 1.56$ (see text).</figcaption></figure>

in both previously mentioned limits. In the following, we explain how this phase diagram was constructed and focus on a specific interaction strength, $c = 10$, where the critical value of $K$ remains small enough to allow for accurate numerical calculations despite our limited system size.

As a first result, we present the energy dynamics in Fig. 1 for three different values of the stochasticity parameter $K$. This provides an initial indication of the phase transition, as we identify three distinct behaviors. For $K < K_c$, the energy saturates at long times, a hallmark of dynamical localization. Conversely, for $K > K_c$, the energy grows unboundedly, exhibiting a linear increase, albeit with corrections due to finite-size effects. Between these two regimes, we identify a critical behavior at $K = K_c$, characterized by anomalous diffusion, with the energy scaling as $\langle E(t) \rangle \sim t^{2/3}$. This scaling is consistent with the predictions of the one-parameter scaling theory in three-dimensional Anderson transitions [3, 15].

We now proceed to a more precise characterization of the transition. Identifying a phase transition is challenging due to finite-time and finite-size effects. This difficulty arises because, at the transition, the localization momentum $p_{\text{loc}}$ diverges and can significantly exceed the system size, which then acts as an upper bound for the effectively observed localization momentum. To circumvent this issue, Lemarié et al. [15] proposed a finite-time scaling analysis inspired by the finite-size scaling method developed for the Anderson model [46, 47]. We closely follow this approach to demonstrate that the observed transition indeed belongs to the orthogonal Anderson universality class.

To apply this method, we first recall some important theoretical results. In the localized regime ($K < K_c$),

<figure><img src="image_4.png" /><figcaption>FIG. 3: (a) The rescaled quantity $\ln \Lambda = \ln[\langle E(t) \rangle t^{-2/3}]$ as a function of $K$ for different times between $t = 5$ (light blue) and $t = 150$ (dark blue). All curves intersect approximately at the critical point ($K_c/\hbar_e \simeq 2.1$, $\ln \Lambda_c \simeq 1.54$) demonstrating the existence of a metal-insulator transition. (b) Determination of the critical exponent $\nu$ by fitting $(\ln \Lambda)'(K_c, t) \sim t^{1/3\nu}$ in log-log scale. Parameters are the same as in Fig. 2.</figcaption></figure>

the energy saturates at long times as $\langle E \rangle \sim p_{\text{loc}}^2$. If an Anderson transition exists, we expect $p_{\text{loc}}$ to diverge near the critical point as $p_{\text{loc}} \sim (K_c - K)^{-\nu}$. On the other side of the transition ($K > K_c$), the dynamics is diffusive, with the energy growing as $\langle E \rangle \sim D(K)t$ where the diffusion constant $D(K)$ vanishes as $(K - K_c)^s$. In three dimensions, it is known that $s = \nu$ [48]. At the critical point ($K = K_c$), anomalous diffusion is expected, characterized by an exponent $2/3$ [15]. The one-parameter scaling theory postulates that the quantity

$$ \Lambda(K,t) = \langle E(t) \rangle t^{-2/3} = f(\xi(K)t^{-1/3}), \quad (10) $$

where $f$ is an arbitrary function, and the scaling parameter $\xi(K)$ depends only on $K$. To test this hypothesis, we present in Fig. 2a our numerical results obtained for $c = 10$, $\hbar_e = 2.89$, and several values of $K$. Each dataset is shifted along the x-axis by a quantity $\ln \xi(K)$ to align all the data points onto the same universal scaling curve. The values of $\xi(K)$ are determined by minimizing the distance between the corresponding values of $\ln(\xi(K)t^{-1/3})$ for each value of $\ln \Lambda$. Remarkably, all data collapse onto a single universal function, with a lower branch corresponding to the localized regime and an upper branch corresponding to the diffusive regime. Figure 2b displays the scaling parameter $\xi(K)$, which diverges at $K_c$ following a power law consistent with the prediction of the orthogonal Anderson universality class, namely $\nu = 1.56$. To model this divergence and account for finite-size effects, we fit our data using the function $\xi(K)^{-1} = \alpha |K - K_c|^{\nu} + \beta$.

Finally, to precisely locate the critical point and extract the critical exponent $\nu$, we plot in Fig. 3a $\ln \Lambda(K, t)$ as a function of $K$ for different time values. This method is particularly convenient since, at the critical point, $\Lambda(K_c, t)$ is time-independent. We thus determine $K_c$ by identifying the crossing points of all curves. Furthermore, we extract the critical exponent by analyzing the behav-
ior of $\ln \Lambda(K, t)$ near $K_c$. Assuming $\xi(K) \sim |K - K_c|^{-\nu}$, the scaling function depends on $(K - K_c)t^{1/3\nu}$, leading to $\ln \Lambda(K, t) \sim \ln \Lambda_c + C_1(K - K_c)t^{1/3\nu} + \dots$. Thus, the slope of $\ln \Lambda$ near $K_c$ is proportional to $t^{1/3\nu}$. By plotting this slope as a function of time on a log-log scale, we determine the critical exponent via a linear fit. As shown in Fig. 3b, this yields $\nu = 1.56 \pm 0.11$, in agreement with the state-of-the-art value for the orthogonal Anderson transition [49], $\nu = 1.57 \pm 0.02$. For comparison, the exponent for the unitary Anderson transition is $\nu = 1.43 \pm 0.04$ [49].

**Conclusion.**— We have investigated the dynamics of a periodically kicked one-dimensional Bose gas with repulsive contact interactions. By establishing a formal correspondence between the system's evolution and an Anderson model in an $N$-dimensional rapidity space, we have demonstrated the existence of an interaction-driven metal-insulator transition for $N \ge 3$. Our extensive numerical simulations for three bosons confirm this transition and allow us to extract the critical stochastic parameter $K_c$ as a function of interactions as well as the associated critical exponent.

This work, grounded in exact calculations, opens new perspectives on the interplay between interactions and Anderson localization. A key open question remains the reconciliation of our findings with conflicting experimental results [40–42]. Additionally, the connection with mean-field approaches, which predict delocalization as soon as interactions are introduced, warrants further investigation. A deeper analysis of our effective model in the infinite-dimensional limit could provide crucial insights into this discrepancy.

We would like to acknowledge helpful discussions with L. Bellando De Castro, J. C. Garreau, B. Grémaud, F. Hébert, M. Landini, C. Miniatura, A. Rançon and F. Werner. This work has benefited from the financial support of Agence Nationale de la Recherche under Grant No ANR-21-CE47-0009 Quantum-SOPHA and Grant No ANR-23-PETQ-0001 Dyn1D France 2030.

Note.—During the publication process, we found that A. Yang et al [50] independently posted a preprint on the same topic at the same time.

[1] P. W. Anderson, Phys. Rev. **109**, 1492 (1958).

[2] D. J. Thouless, Physics Reports **13**, 93 (1974).

[3] E. Abrahams, P. W. Anderson, D. C. Licciardello, and T. V. Ramakrishnan, Phys. Rev. Lett. **42**, 673 (1979).

[4] T. Giamarchi and H. J. Schulz, Phys. Rev. B **37**, 325 (1988).

[5] D. M. Basko, I. L. Aleiner, and B. Altshuler, Ann. Phys. **321**, 1126 (2006).

[6] R. Nandkishore and D. Huse, Annu. Rev. Condens. Matter Phys. **6**, 15 (2015).

[7] F. Alet and N. Laflorencie, Comptes Rendus Physique **19**, 408 (2018).

[8] P. Ponte, Z. Papić, F. Huveneers, and D. A. Abanin, Phys. Rev. Lett. **114**, 140401 (2015).

[9] P. Ponte, A. Chandran, Z. Papić, and D. A. Abanin, Ann. Phys. **353**, 196 (2015).

[10] S. Fishman, D. R. Grempel, and R. E. Prange, Phys. Rev. Lett. **49**, 509 (1982).

[11] D. R. Grempel, R. E. Prange, and S. Fishman, Phys. Rev. A **29**, 1639 (1984).

[12] D. Shepelyansky, Phys. Rev. Lett. **56**, 677 (1986).

[13] F. L. Moore, J. C. Robinson, C. F. Bharucha, Bala Sundaram, and M. G. Raizen, Phys. Rev. Lett. **75**, 4598 (1995).

[14] J. Chabé, G. Lemarié, B. Grémaud, D. Delande, P. Szriftgiser, and J. C. Garreau, Phys. Rev. Lett. **101**, 255702 (2008).

[15] G. Lemarié, J. Chabé, P. Szriftgiser, J. C. Garreau, B. Grémaud, and D. Delande, Phys. Rev. A **80**, 043626 (2009).

[16] J. Wang, and A. M. Garcia-Garcia, Phys. Rev. E **79**, 036206 (2009).

[17] G. Lemarié, H. Lignier, D. Delande, P. Szriftgiser, and J. C. Garreau, Phys. Rev. Lett. **105**, 090601 (2010).

[18] M. Lopez, J. F. Clément, P. Szriftgiser, J. C. Garreau, and D. Delande, Phys. Rev. Lett. **108**, 095701 (2012).

[19] C. Hainaut, I. Manai, J. F. Clément, J. C. Garreau, P. Szriftgiser, G. Lemarié, N. Cherroret, D. Delande, and R. Chicireanu, Nat. Commun. **9**, 1382 (2018).

[20] G. Lemarié, C. A. Müller, D. Guéry-Odelin, and C. Miniatura, Phys. Rev. A **95**, 043626 (2017).

[21] M. Lewenstein, A. Sanpera, V. Ahufinger, B. Damski, A. S. De, and U. Sen, Adv. in Physics **56**, 243 (2007).

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

[23] M. Lewenstein, A. Sanpera, and V. Ahufinger, *Ultracold atoms in optical lattices: Simulating Quantum Many-Body Physics* (Oxford University Press, Oxford, U.K. 2012).

[24] E. H. Lieb, and W. Liniger, Phys. Rev. **130**, 1605 (1963).

[25] B. Paredes, A. Widera, V. Murg, O. Mandel, S. Fölling, I. Cirac, G. V. Shlyapnikov, T. W. Hänsch, and I. Bloch, Nature **429**, 227 (2004).

[26] T. Kinoshita, T. Wenger, and D. S. Weiss, Science **305**, 1125 (2004).

[27] T. Kinoshita, T. Wenger, and D. S. Weiss, Phys. Rev. Lett. **95**, 190406 (2005).

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

[29] L. Tonks, Phys. Rev. **50**, 955 (1936).

[30] M. Girardeau, J. Math. Phys. **1**, 516 (1960).

[31] V. Vuatelet, and A. Rançon, Phys. Rev. A **104**, 043302 (2021).

[32] C. Rylands, E. B. Rozenbaum, V. Galitski, and R. Konik, Phys. Rev. Lett. **124**, 155302 (2020).

[33] R. Chicireanu, and A. Rançon, Phys. Rev. A **103**, 043314 (2021).

[34] D. L. Shepelyansky, Phys. Rev. Lett. **70**, 1787 (1993).

[35] A. S. Pikovsky, and D. L. Shepelyansky, Phys. Rev. Lett. **100**, 094101 (2008).

[36] S. Flach, D. O. Krimer, and Ch. Skokos, Phys. Rev. Lett. **102**, 024101 (2009).

[37] G. Gligorić, J. D. Bodyfelt, and S. Flach, Europhys. Lett. **96**, 30004 (2011).

[38] N. Cherroret, B. Vermersch, J. C. Garreau, and D. De-
lande, Phys. Rev. Lett. **112**, 170603 (2014).

[39] S. Lellouch, A. Rançon, S. De Bièvre, D. Delande, and J. C. Garreau, Phys. Rev. A **101**, 043624 (2020).

[40] A. Cao, R. Sajjad, H. Mas, E. Q. Simmons, J. L. Tanlimco, E. Nolasco-Martinez, T. Shimasaki, H. E. Kondakci, V. Galitski, and D. M. Weld, Nat. Phys. **18**, 1302 (2022).

[41] J. H. S. Toh, K. C. McCormick, X. Tang, Y. Su, X.-W. Luo, C. Zhang, S. Gupta, Nature Physics **18**, 1297 (2022).

[42] Y. Guo, S. Dhar, A. Yang, Z. Chen, H. Yao, M. Horvath, L. Ying, M. Landini, and H.-C. Nägerl, Science **389**, 716 (2025).

[43] S. Notarnicola, F. Iemini, D. Rossini, R. Fazio, A. Silva, and A. Russomanno, Phys. Rev. E **97**, 022202 (2018).

[44] A. D. Mirlin, Y. V. Fyodorov, F.-M. Dittes, J. Quezada, and T. H. Seligman, Phys. Rev. E **54**, 3221 (1996).

[45] See Supplemental Material [url] for additional information about the derivation of the equations, which includes

Refs. [51–55].

[46] A. MacKinnon, and B. Kramer, Phys. Rev. Lett. **47**, 1546 (1981).

[47] J. L. Pichard, and G. Sarma, J. Phys. C **14**, L127 (1981).

[48] F. Wegner, Z. Phys. **B25**(1976).

[49] K. Slevin and T. Ohtsuki, Phys. Rev. Lett. **78**, 4083 (1997).

[50] A. Yang, Z. Chen, Y. Guo, M. Landini, H-C Nägerl, and L. Ying, arXiv:2503.04683 (2025).

[51] A. Lamacraft, https://austen.uk/courses/tqm/lieb-liniger-model.

[52] V. E. Korepin, N. M. Bogoliubov, and A. G. Izergin, *Quantum Inverse Scattering Method and Correlation Functions* (Cambridge University Press, Cambridge, 1983).

[53] N. A. Slavnov, Theoret. and Math. Phys. **79**, 502 (1989).

[54] J. De Nardis and M. Panfil, J. Stat. Mech. P02019 (2015).

[55] E. Granet and F. Essler, Sci. Post. Phys. **9**, 082 (2020).
Supplemental material: Interaction induced Anderson transition in a kicked one dimensional Bose gas

H. Olsen,<fnref n="1" /> P. Devillard,<fnref n="2" /> G. Aupetit-Diallo,<fnref n="3" /> P. Vignolo,<fnref n="1" />,<fnref n="4" /> and M. Albert<fnref n="1" />,<fnref n="4" />

BETHE ANSATZ SOLUTION FOR THREE PARTICLES

In this section, we recall the expressions of the eigenstates of Hamiltonian (1) from the main text for the case of $N = 3$ identical bosons [1, 2]. Following the Bethe Ansatz, the wavefunction in each position sector is expressed as a superposition of plane waves with fixed energy. In the fundamental position sector $x_1 < x_2 < x_3$, the Bethe wavefunction takes the form

$$ \Psi(x_1, x_2, x_3) = A_{123}e^{i(\lambda_1 x_1 + \lambda_2 x_2 + \lambda_3 x_3)} + A_{132}e^{i(\lambda_1 x_1 + \lambda_3 x_2 + \lambda_2 x_3)} \\ + A_{213}e^{i(\lambda_2 x_1 + \lambda_1 x_2 + \lambda_3 x_3)} + A_{231}e^{i(\lambda_2 x_1 + \lambda_3 x_2 + \lambda_1 x_3)} \\ + A_{312}e^{i(\lambda_3 x_1 + \lambda_1 x_2 + \lambda_2 x_3)} + A_{321}e^{i(\lambda_3 x_1 + \lambda_2 x_2 + \lambda_1 x_3)}, \quad (1) $$

where the positions $x_j$ are expressed in units of $\ell = L/(2\pi)$. The non-normalized amplitudes $A_{ijk}$, which depend on the rapidities $\lambda_j$ and the interaction strength $c$, are determined by enforcing the cusp condition at $x_i = x_j$. This accounts for the discontinuity in the derivative due to the contact interactions [1, 2]. The amplitudes are given by

$$ A_{123}(\lambda_1, \lambda_2, \lambda_3, c) = (\lambda_1 - \lambda_2 + ic)(\lambda_1 - \lambda_3 + ic)(\lambda_2 - \lambda_3 + ic) $$

$$ A_{132}(\lambda_1, \lambda_2, \lambda_3, c) = (\lambda_1 - \lambda_2 + ic)(\lambda_1 - \lambda_3 + ic)(\lambda_2 - \lambda_3 - ic) $$

$$ A_{213}(\lambda_1, \lambda_2, \lambda_3, c) = (\lambda_1 - \lambda_2 - ic)(\lambda_1 - \lambda_3 + ic)(\lambda_2 - \lambda_3 + ic) $$

$$ A_{312}(\lambda_1, \lambda_2, \lambda_3, c) = (\lambda_1 - \lambda_2 + ic)(\lambda_1 - \lambda_3 - ic)(\lambda_2 - \lambda_3 - ic) $$

$$ A_{321}(\lambda_1, \lambda_2, \lambda_3, c) = (\lambda_1 - \lambda_2 - ic)(\lambda_1 - \lambda_3 - ic)(\lambda_2 - \lambda_3 - ic) $$

$$ A_{231}(\lambda_1, \lambda_2, \lambda_3, c) = (\lambda_1 - \lambda_2 - ic)(\lambda_1 - \lambda_3 - ic)(\lambda_2 - \lambda_3 + ic). $$

Finally the rapidities $\lambda_j$ are obtained by imposing periodic boundary conditions, leading to the following system of coupled equations

$$ \begin{aligned} \lambda_1 &= J_1 + \frac{1}{2\pi}\theta(\lambda_1 - \lambda_2) + \frac{1}{2\pi}\theta(\lambda_1 - \lambda_3), \\ \lambda_2 &= J_2 + \frac{1}{2\pi}\theta(\lambda_2 - \lambda_1) + \frac{1}{2\pi}\theta(\lambda_2 - \lambda_3), \\ \lambda_3 &= J_3 + \frac{1}{2\pi}\theta(\lambda_3 - \lambda_1) + \frac{1}{2\pi}\theta(\lambda_3 - \lambda_2), \end{aligned} \quad (2) $$

where $\theta(x) = -2 \arctan(x/c)$, and $J_1, J_2$ and $J_3$ are distinct relative integers.

CALCULATION OF THE MATRIX ELEMENTS [EQ.(9)]

We want to compute

$$ M_{\bar{\lambda}, \bar{\mu}} = \langle \lambda_1, \dots, \lambda_N | e^{-iK/\hbar_e \sum_j \cos x_j} | \mu_1, \dots, \mu_N \rangle, \quad (3) $$

namely

$$ M_{\bar{\lambda}, \bar{\mu}} = \int \prod_j dx_j e^{-iK/\hbar_e \cos x_j} \Psi_{\lambda_1, \dots, \lambda_N}^*(x_1, \dots, x_N) \Psi_{\mu_1, \dots, \mu_N}(x_1, \dots, x_N), \quad (4) $$

<footnote n="1">¹Université Côte d'Azur, CNRS, Institut de Physique de Nice, France</footnote>

<footnote n="2">²Aix Marseille Univ, Université de Toulon, CNRS, CPT, Marseille, France</footnote>

<footnote n="3">³SISSA, Via Bonomea 265, I-34136 Trieste, Italy</footnote>

<footnote n="4">⁴Institut Universitaire de France</footnote>
with

$$ \Psi_{\lambda_1, \dots, \lambda_N}(x_1, \dots, x_N) = \sum_Q A_Q(\{\lambda_i\}, c) \sum_P \theta_P(x_1, \dots, x_N) e^{i \sum_j x_{P(j)} \lambda_{Q(j)}}, \quad (5) $$

being the many-body wavefunction for N Lieb-Liniger bosons restricted to the fundamental rapidity sector $\lambda_1 < \lambda_2 < \dots < \lambda_N$ where Q and P are permutations of $S_N$. Putting expression (5) in Eq. (4), one obtains

$$ M_{\vec{\lambda}, \vec{\mu}} = \sum_{Q, Q'} A_Q^* A_{Q'} \int \prod_j dx_j e^{-iK/\hbar_e \cos x_j} \sum_{P, P'} \theta_P(x_1, \dots, x_N) \theta_{P'}(x_1, \dots, x_N) e^{-i \sum_q (x_{P(q)} \lambda_{Q(q)} - x_{P'(q)} \mu_{Q'(q)})}, \quad (6) $$

with $A_Q = A_Q(\{\lambda_i\}, c)$ and $A_{Q'} = A_{Q'}(\{\mu_i\}, c)$ which are defined in Eq. (4) of the main text. This is only different from 0 when $P = P'$, which leads to

$$ M_{\vec{\lambda}, \vec{\mu}} = \sum_{Q, Q'} A_Q^* A_{Q'} \sum_P \int_{\Gamma_P} \prod_j dx_j e^{-iK/\hbar_e \cos x_j} e^{-i \sum_q x_{P(q)} (\lambda_{Q(q)} - \mu_{Q'(q)})}, \quad (7) $$

with $\Gamma_P$ delimiting the integration domain $x_{P(1)} < \dots < x_{P(N)}$. As every sector is equivalent, we can write the matrix element in a factorized form as

$$ M_{\vec{\lambda}, \vec{\mu}} = \sum_{Q, Q'} A_{Q, Q'} S_{Q, Q'}, \quad (8) $$

with $A_{Q, Q'} = A_Q^* A_{Q'}$ and

$$ S_{Q, Q'} = N! \int_{x_1 < \dots < x_N} \prod_j dx_j e^{-iK/\hbar_e \cos x_j} e^{-ix_j(\lambda_{Q(j)} - \mu_{Q'(j)})}. \quad (9) $$

While dealing with periodic boundary conditions (PBC), this formulation of the matrix element becomes ill-defined. The boundary condition imposes for every ordering to loop, which results in the integration domain having the form $x_1 < \dots < x_N < x_1$. Writing the sector this way may suggest that the integration domain has no dimension. However, we can overcome this issue by treating one of the variables (in this case, $x_1$) as a moving boundary for the other ($N-1$) variables. The corresponding domain is actually $N$ times larger than $x_1 < \dots < x_N$ and leads to the ensuing expression

$$ S_{Q, Q'} = N! \int_0^{L/\ell} dx_1 e^{-iK/\hbar_e \cos x_1} e^{-ix_1(\lambda_{Q(1)} - \mu_{Q'(1)})} \prod_{j=2}^{N} \int_{x_{j-1}}^{L/\ell} dx_j e^{-iK/\hbar_e \cos x_j} e^{-ix_j(\lambda_{Q(j)} - \mu_{Q'(j)})}. \quad (10) $$

By using the following identity

$$ e^{-iK/\hbar_e \cos x} = \sum_{n=-\infty}^{\infty} (-i)^n J_n(K/\hbar_e) e^{-inx}, \quad (11) $$

the S function takes the form

$$ S_{Q, Q'} = N! \sum_{m_1, \dots, m_N} (-i)^{\sum_j m_j} J_{m_1}(K/\hbar_e) \cdots J_{m_N}(K/\hbar_e) s_{Q, Q', \vec{l}}, \quad (12) $$

with

$$ s_{Q, Q', \vec{l}} = \int_0^{2\pi} dx_1 e^{-il_1 x_1} \prod_{j=2}^{N} \int_{x_{j-1}}^{2\pi} dx_j e^{-il_j x_j}, \quad (13) $$

where

$$ l_j = (\lambda_{Q(j)} - \mu_{Q'(j)}) + m_j. \quad (14) $$
<figure><img src="image_5.png" /><figcaption>FIG. 1: Matrix elements $|M_{\tilde{\lambda}, \tilde{\mu}}|$ for the case $c = 10$, $I_1 = J_1 = -1$, $I_2 = J_2 = 0$, as functions of $|J_3 - I_3|$, for different values of the stochasticity parameter $K/\hbar_e$.</figcaption></figure>

By performing the integration in Eq. (13) for $N = 3$, one obtains

$$s_{Q,Q',\vec{l}} = \frac{ie^{-i2\pi(l_1+l_2+l_3)}}{l_1 l_2 (l_1+l_2) l_3 (l_2+l_3) (l_1+l_2+l_3)} \left( e^{i2\pi(l_1+l_2+l_3)} l_1 l_2 (l_1+l_2) - l_2 l_3 (l_2+l_3) + e^{i2\pi l_1} (l_1+l_2) l_3 (l_1+l_2+l_3) - e^{i2\pi(l_1+l_2)} l_1 (l_2+l_3) (l_1+l_2+l_3) \right). \quad (15)$$

Expression (15) has some artificial singularities that can be eliminated by considering separately all the cases where the denominator vanishes as follows:

$$s_{Q,Q',\vec{l}} = \begin{cases} -2\pi e^{-i2\pi(l_1+l_2+l_3)} \left( \frac{1}{l_2+l_3} e^{i\pi l_1} \Phi_1 + \frac{1}{l_2+l_3} e^{i\pi(l_1+l_2)} \Phi_2 \right) & \text{if } l_3 \neq -l_2 \\ e^{-i2\pi l_1} \left( \frac{ie^{i2\pi(l_1+l_2)}}{l_2^2(l_1+l_2)} + \frac{e^{i2\pi l_1}(il_2+l_1(-i+2\pi l_2))}{l_1^2 l_2^2} - \frac{i}{l_1^2(l_1+l_2)} \right) & \text{if } l_3 = -l_2, l_1 \neq -l_2, l_1 \neq 0, l_2 \neq 0 \\ \frac{-i+ie^{i2\pi l_2}+2\pi l_2(1+i\pi l_2)}{l_2^3} & \text{if } l_3 = -l_2, l_1 = 0, l_2 \neq 0 \\ \frac{i-ie^{-i2\pi l_1}+2\pi l_1-i2\pi^2 l_1^2}{l_1^3} & \text{if } l_2 = 0, l_3 = 0, l_1 \neq 0 \\ \frac{-2(i+\pi l_1+e^{-i2\pi l_1}(-i+\pi l_1))}{l_1^3} & \text{if } l_3 = -l_2, l_1 = -l_2, l_1 \neq 0 \\ \frac{4\pi^3}{3} & \text{if } l_1 = 0, l_2 = 0, l_3 = 0 \end{cases} \quad (16)$$

where

$$\Phi_1 = \begin{cases} \frac{\text{Sinc}(\pi l_1) - e^{i\pi l_2} \text{Sinc}(\pi(l_1+l_2))}{-\pi l_1 \cos(\pi l_1) + (1-i\pi l_1) \sin(\pi l_1)} & \text{if } l_2 \neq 0 \\ \frac{l_2}{\pi l_1^2} & \text{if } l_1 \neq 0 \\ -i\pi & \text{otherwise} \end{cases} \quad (17)$$

and

$$\Phi_2 = \begin{cases} \frac{-\text{Sinc}(\pi(l_1+l_2))+e^{i\pi l_3}\text{Sinc}(\pi(l_1+l_2+l_3))}{\pi(l_1+l_2)\cos(\pi(l_1+l_2)) + (-1+i\pi(l_1+l_2))\sin(\pi(l_1+l_2))} & \text{if } l_3 \neq 0 \\ \frac{\pi(l_1+l_2)\cos(\pi(l_1+l_2)) + (-1+i\pi(l_1+l_2))\sin(\pi(l_1+l_2))}{\pi(l_1+l_2)^2} & \text{if } l_1 \neq -l_2 \\ i\pi & \text{otherwise,} \end{cases} \quad (18)$$

with $\text{Sinc}(z) = \sin(z)/z$. We show in Fig. 1 the behaviour of $|M_{\tilde{\lambda}, \tilde{\mu}}|$, where $\tilde{\lambda}$ and $\tilde{\mu}$ are determined by the sets of relative integers $\vec{I}$ and $\vec{J}$ respectively (see Eq. (2)), for the case $c = 10$, $I_1 = J_1 = -1$, $I_2 = J_2 = 0$, as functions of $|J_3 - I_3|$. We can observe that for the range of values of the stochasticity parameter $K/\hbar_e$ that are relevant for this work, the matrix elements $|M_{\tilde{\lambda}, \tilde{\mu}}|$ decrease exponentially with respect the rapidities distance.
Remark that it would be possible to express the matrix elements $M_{\tilde{\chi}, \tilde{\mu}}$ as functions of form factors [3–5]. This would allow to have an analytical formula for the limiting case $1/c \ll 1$ and an arbitrary number of particles $N$ [6], but the method proposed here is more efficient for the numerical calculation at any value of the interaction strength and for the specific case of $N = 3$. We leave the larger $N$ case for future investigations.

## Condition for dimensionality reduction

Here we formally demonstrate that a system of $N$ kicked Lieb-Liniger bosons reduces to $N$ kicked independent particles ($N$ independent 1D systems) only in the non-interacting regime ($c = 0$) and in the fermionized regime ($c \to \infty$). When this occurs, the matrix elements

$$M_{\tilde{\chi}, \tilde{\mu}} = \langle \lambda_1, \dots, \lambda_N | e^{-iK/\hbar_e \sum_j \cos x_j} | \mu_1, \dots, \mu_N \rangle, \quad (19)$$

can be rewritten as

$$M_{\tilde{\chi}, \tilde{\mu}} = \sum_{Q, Q'} \beta_{Q, Q'} \prod_{j=1}^{N} \langle \lambda_{Q(j)} | e^{-iK/\hbar_e \cos x_j} | \mu_{Q'(j)} \rangle, \quad (20)$$

which explicitly factorizes in terms of single-particle contributions. The term $\beta_{Q, Q'}$ and the sum over $Q$ and $Q'$ account for the symmetry properties of the many-body wavefunction. Starting from expression (7) for the matrix elements, we have

$$M_{\tilde{\chi}, \tilde{\mu}} = \sum_{Q, Q'} A_Q^* A_{Q'} \sum_P \int_{\Gamma_P} \prod_j dx_j e^{-iK/\hbar_e \cos x_j} e^{-i \sum_q x_{P(q)} (\lambda_{Q(q)} - \mu_{Q'(q)})}. \quad (21)$$

Without loss of generality, we can set $j = P(q)$, leading to

$$M_{\tilde{\chi}, \tilde{\mu}} = \sum_{Q, Q'} A_Q^* A_{Q'} \sum_P \int_{\Gamma_P} \prod_j dx_j e^{-iK/\hbar_e \cos x_j} e^{-ix_j (\lambda_{Q-P-1}(q) - \mu_{Q'-P-1}(q))}. \quad (22)$$

Eq. (22) simplifies to

$$M_{\tilde{\chi}, \tilde{\mu}} = \sum_{\tilde{Q}, \tilde{Q}'} A_{\tilde{Q}}^* A_{\tilde{Q}'} \prod_j \int_0^{2\pi} dx_j e^{-iK/\hbar_e \cos x_j} e^{-ix_j (\lambda_{\tilde{Q}(q)} - \mu_{\tilde{Q}'(q)})}. \quad (23)$$

where we define $\tilde{Q} = QP^{-1}$ and $\tilde{Q}' = Q'P^{-1}$. This factorization holds only if $A_{\tilde{Q}}^* A_{\tilde{Q}'} = A_Q^* A_{Q'}$ for all $Q, Q'$, and $P$, which translates into the condition

$$\prod_{1 \le j, j' < k, k' \le N} \left( 1 - \frac{ic}{\lambda_{Q(j)} - \lambda_{Q(k)}} \right) \left( 1 + \frac{ic}{\mu_{Q'(j')} - \mu_{Q'(k')}} \right) = \prod_{1 \le j, j' < k, k' \le N} \left( 1 - \frac{ic}{\lambda_{\tilde{Q}(j)} - \lambda_{\tilde{Q}(k)}} \right) \left( 1 + \frac{ic}{\mu_{\tilde{Q}'(j')} - \mu_{\tilde{Q}'(k')}} \right). \quad (24)$$

This equality holds for all $Q, Q'$, and $P$ only if $c = 0$ or $c \to \infty$.

## MAPPING TO THE ANDERSON MODEL

In this section we detail the procedure to get Eq. (7) of the main text from the Floquet equation $U|\phi_\alpha\rangle = e^{-i\frac{\mu_\alpha}{\hbar_e}}|\phi_\alpha\rangle$, with $U = e^{-iH_0/\hbar_e}e^{-iH_K/\hbar_e}$. By defining $|\tilde{\phi}_\alpha\rangle = e^{iH_0/2\hbar_e}|\phi_\alpha\rangle$, one obtains the time-reversal symmetrized Floquet equation

$$\tilde{U}|\tilde{\phi}_{\alpha}\rangle = e^{-i\frac{\mu_{\alpha}}{\hbar_e}}|\tilde{\phi}_{\alpha}\rangle, \quad (25)$$
with

$$ \tilde{U} = e^{-iH_0/(2\hbar_e)} e^{-iH_K/\hbar_e} e^{-iH_0/(2\hbar_e)}. \quad (26) $$

Furthermore, using the identity $\tilde{U}^\dagger U |\tilde{\phi}_\alpha\rangle = |\tilde{\phi}_\alpha\rangle$, the following equation

$$ \tilde{U}^\dagger |\tilde{\phi}_\alpha\rangle = e^{i \frac{\mu_\alpha}{\hbar_e}} |\tilde{\phi}_\alpha\rangle \quad (27) $$

holds too. By projecting Eqs. (25) and (27) onto the basis state $\langle \vec{\lambda} |$, one gets

$$ \sum_{\vec{\nu}, \vec{\mu}, \vec{\beta}} \langle \vec{\lambda} | e^{\mp i H_0 / (2\hbar_e)} | \vec{\nu} \rangle \langle \vec{\nu} | e^{\mp i H_K / \hbar_e} | \vec{\mu} \rangle \langle \vec{\mu} | e^{\mp i H_0 / (2\hbar_e)} | \vec{\beta} \rangle \langle \vec{\beta} | \tilde{\phi}_\alpha \rangle = e^{-i \frac{\mu_\alpha}{\hbar_e}} \langle \vec{\lambda} | \tilde{\phi}_\alpha \rangle, \quad (28) $$

namely

$$ \sum_{\vec{\mu}} e^{\mp i(E_{\vec{\lambda}} + E_{\vec{\mu}})/(2\hbar_e)} \langle \vec{\lambda} | e^{\mp i H_K / \hbar_e} | \vec{\mu} \rangle \langle \vec{\mu} | \tilde{\phi}_\alpha \rangle = e^{\mp i \frac{\mu_\alpha}{\hbar_e}} \langle \vec{\lambda} | \tilde{\phi}_\alpha \rangle. \quad (29) $$

By defining the matrix elements $M_{\vec{\lambda},\vec{\mu}} = \langle \vec{\lambda} | e^{-iH_K/\hbar_e} | \vec{\mu} \rangle$ and by setting

$$ W_{\vec{\lambda}, \vec{\mu}} = 2 \operatorname{Re} \left\{ \exp\left(-i \frac{E_{\vec{\lambda}} + E_{\vec{\mu}}}{2\hbar_e}\right) M_{\vec{\lambda}, \vec{\mu}} \right\}, \quad (30) $$

one gets straightforwardly Eq. (7) of the main text:

$$ \sum_{\vec{\mu} \neq \vec{\lambda}} W_{\vec{\lambda}, \vec{\mu}} \langle \vec{\mu} | \tilde{\phi}_{\alpha} \rangle + W_{\vec{\lambda}, \vec{\lambda}} \langle \vec{\lambda} | \tilde{\phi}_{\alpha} \rangle = 2 \cos\left(\frac{\mu_{\alpha}}{\hbar_e}\right) \langle \vec{\lambda} | \tilde{\phi}_{\alpha} \rangle. \quad (31) $$

## NUMERICAL PROCEDURE

Once the matrix elements are obtained, the main computational effort is completed. The evolution operator over one period is given by:

$$ \hat{U} = e^{-i \frac{\hat{H}_0}{\hbar_e}} e^{-i \frac{\hat{H}_K}{\hbar_e}}. \quad (32) $$

The matrix elements of the evolution operator in the Bethe basis $|\vec{\lambda}\rangle$ are:

$$ U_{\{\lambda_i\}, \{\mu_i\}} = \langle \lambda_1, \dots, \lambda_N | \hat{U} | \mu_1, \dots, \mu_N \rangle, \quad (33) $$

which can be written as:

$$ U_{\{\lambda_i\}, \{\mu_i\}} = e^{-i E_{\lambda} / \hbar_e} \langle \lambda_1, \dots, \lambda_N | e^{-i \frac{\hat{H}_K}{\hbar_e}} | \mu_1, \dots, \mu_N \rangle. \quad (34) $$

From our previous calculations, the matrix elements of the evolution operator are given by:

$$ U_{\{\lambda_i\}, \{\mu_i\}} = e^{-i \frac{\hbar_e}{2} (\lambda_1^2 + \dots + \lambda_N^2)} M_{\vec{\lambda}, \vec{\mu}}. \quad (35) $$

We recall that $M_{\vec{\lambda},\vec{\mu}}$ is given by:

$$ M_{\vec{\lambda}, \vec{\mu}} = \sum_{Q, Q'} A_{Q, Q'} N! \sum_{|m_1|, \dots, |m_N| \le M} (-i)^{\sum_j m_j} J_{m_1}(\kappa) \cdots J_{m_N}(\kappa) s_{Q, Q'}(\{l_i\}). \quad (36) $$

Here $M$ acts as a cut-off in the sum over Bessel functions. According to asymptotic properties of Bessel functions we use $M > [2K/\hbar_e]$ where $[\cdot]$ is the integer part. Since the full Hilbert space is infinite-dimensional, we truncate it to make numerical computations feasible. As explained in the main text, we restrict Bethe numbers to the range $-N_s \le \lambda_i \le N_s$ with $N_s$ up to 29. Furthermore, due to the symmetry of the wave function, it is sufficient to consider only ordered sets of Bethe numbers. This results in a matrix of size $N_m = \binom{2N_s+1}{3}$. Finally, we have to ensure that
the matrix elements are properly normalized. Within the definition of the Bethe amplitude given in this appendix this is not the case. Proper normalization is restored by dividing each matrix elements by $\sqrt{\langle\vec{\lambda}|\vec{\lambda}\rangle\langle\vec{\mu}|\vec{\mu}\rangle}$. These norms can be directly evaluated using $M_{\vec{\lambda},\vec{\lambda}}$ at $K = 0$.

To compute the time evolution of the system, we expand its wave function in the Lieb-Liniger basis as:

$$|\psi(t)\rangle = \sum_{\{\lambda_i\}} \alpha_{\{\lambda_i\}}(t) |\lambda_1, \dots, \lambda_N\rangle, \qquad (37)$$

where the coefficients $\alpha_{\{\lambda_i\}}(t)$ evolve stroboscopically as:

$$\alpha_{\{\lambda_i\}}(t+1) = \sum_{\{\mu_i\}} U_{\{\lambda_i\},\{\mu_i\}}(\kappa) \alpha_{\{\mu_i\}}(t), \qquad (38)$$

with $\alpha_{\{\lambda_i\}}(0) = 1$ in the ground state and zero otherwise.

The total energy of the system at time $t$ is given by:

$$\langle E(t) \rangle = \sum_{\{\lambda_i\}} |\alpha_{\{\lambda_i\}}(t)|^2 E_{\{\lambda_i\}}, \qquad (39)$$

where $E_{\{\lambda_i\}} = \frac{\hbar^2}{2}(\lambda_1^2 + \dots + \lambda_N^2)$ is the energy associated with the quantum numbers.

[1] E. H. Lieb, and W. Liniger, Phys. Rev. 130, 1605 (1963).

[2] A. Lamacraft, https://austen.uk/courses/tqm/lieb-liniger-model.

[3] V. E. Korepin, N. M. Bogoliubov, and A. G. Izergin, Quantum Inverse Scattering Method and Correlation Functions (Cambridge University Press, Cambridge, 1983).

[4] N. A. Slavnov, Theoret. and Math. Phys. 79, 502 (1989).

[5] J. De Nardis and M. Panfil, J. Stat. Mech. P02019 (2015).

[6] E. Granet and F. Essler, Sci. Post. Phys. 9, 082 (2020).