# Superfluidity vs. prethermalisation in a nonlinear Floquet system

Sen Mu, Nicolas Macé, Jiangbin Gong, Christian Miniatura, Gabriel Lemarié, Mathias Albert

► To cite this version:

Sen Mu, Nicolas Macé, Jiangbin Gong, Christian Miniatura, Gabriel Lemarié, et al. Superfluidity vs. prethermalisation in a nonlinear Floquet system. EPL - Europhysics Letters, 2022, 140 (5), pp.50001. 10.1209/0295-5075/aca4f0. hal-03892288

HAL Id: hal-03892288

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

Submitted on 7 Aug 2025

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
Superfluidity vs prethermalisation in a nonlinear Floquet system

S. MU¹, N. MACÉ³, J. GONG¹,², C. MINIATURA²,⁴, G. LEMARIÉ²,³,⁴ and M. ALBERT⁵

¹ Department of Physics, National University of Singapore, Singapore 117542, Singapore

² Centre for Quantum Technologies, National University of Singapore, Singapore 117543, Singapore

³ Laboratoire de Physique Théorique, IRSAMC, Université de Toulouse, CNRS, UPS, France

⁴ MajuLab, International Joint Research Unit UMI 3654, CNRS, Université Côte d'Azur, Sorbonne Université, National University of Singapore, Nanyang Technological University, Singapore

⁵ Université Côte d'Azur, CNRS, INPHYNI, France

**Abstract** – We show that superfluidity can be used to prevent thermalisation in a nonlinear Floquet system. Generically, periodic driving boils an interacting system to a featureless infinite temperature state. Fast driving is a known strategy to postpone Floquet heating with a large but always finite boiling time. In contrast, using a nonlinear periodically-driven system on a lattice, we show the existence of a continuous class of initial states which do not thermalise at all. This absence of thermalisation is associated to the existence and persistence of a stable superflow motion.

**Introduction.** – Thermalisation of isolated quantum systems has been the subject of intensive research during the past decade. In this context, it is important to distinguish ergodic systems, who dynamically explore all possible energetically accessible states, from non-ergodic ones who can only explore a limited fraction of such states. Beside fine-tuned integrable systems whose dynamics is constrained by constants of motion, a seminal example of a nontrivial non-ergodic behavior is many-body localisation happening in strongly disordered quantum interacting systems [1, 2]. Other notable ergodicity breaking phenomena include quantum many-body scars and Hilbert space fragmentation [3, 4].

In conservative systems, the dynamics can only take place within the energy shell since energy is conserved. In time-dependent systems, this is no longer true and time-periodic driving has been used to realize new nonequilibrium phases of matter, e.g. Floquet topological insulator and discrete time crystal [5, 6]. Generically, periodically-driven systems subject to interaction tend to evolve toward a featureless state akin to infinite temperature at long times [7–9]. Different strategies have been devised to escape or hamper this heating mechanism due to energy exchanges. A first example is achieving many-body localisation through strong disorder in quantum systems [10–13]. Another example is to bring the system into a long-lived prethermalization state where Floquet energy transfers are largely suppressed and heating is

postponed to exponentially large times. This prethermal state can be achieved by driving the system sufficiently fast. This mechanism is interesting because it is universal as it applies to both quantum [14–20] and classical systems [21–24] and does not rely on strong disorder. Nevertheless fast driving does not prevent ultimately the system to boil to an infinite temperature state.

In this Letter, we investigate the robustness of superfluidity in interacting Floquet systems as it could offer an escape route against prethermalisation. Indeed, in conservative systems, a superfluid at zero temperature [25] is immune to small enough perturbations. As long as its velocity remains smaller than a well defined critical velocity [26], scattering by impurities is suppressed and internal excitations cannot be activated, preventing energy redistribution and thus thermalisation. It is only above this threshold that the system may enter a route toward wave thermalisation [27] and eventually reach a statistical equilibrium with an effective temperature. First discovered in liquid Helium [28, 29], superfluidity was later shown to be more universal and was observed in various quantum fluids [30–33].

Recent studies show that quenched superfluids generally evolve to a thermal equilibrium state at large times, see e.g. [34, 35]. In this Letter, we show that, under suitable conditions, superfluidity can be maintained at all times in Floquet systems and escape prethermalisation and heating. Using a lattice model with periodically-kicked non-
linear interaction and onsite potentials, we build a driving protocol allowing the system to remain superfluid despite energy injection from the kicks, and show that the time $\tau_{boil}$ needed to boil the system diverges at the prethermal to superfluid transition.

Our paper is organized as follows. We first introduce our periodically-driven system and the new driving protocol allowing for superfluidity to emerge dynamically. We discuss its mapping to, and difference with, conservative superfluids. Superfluidity and heating features are then characterized by gradually increasing the complexity of the model. In the clean case, we perform a Bogoliubov analysis to unveil two important characteristic scales, the sound velocity and the healing length. We also describe the Floquet heating instability due to the presence of a nonlinearity. We then study the effect of a single impurity. Two distinct phases are found, a superfluid and a prethermal phase, exhibiting dramatically different dynamical behaviours. A superfluid flow is observed up to a critical velocity given by a Landau criterion [26, 36, 37]. Above that, the system reaches a prethermal phase with a large but finite boiling time. Strikingly, we find that the boiling time diverges at the transition between the superfluid and prethermal regions (up to the longest times accessible numerically). Thus, thermalisation is absent in the superfluid phase. Finally, we show that this phenomenology remains valid in the presence of disorder. We describe the statistical properties of the critical velocity using extreme value statistics and we find a very good agreement with numerical results. Technical details are given in the supplementary material (SM).

**Model.** — Driven systems with periodically-kicked onsite potentials have attracted extensive attention to study the dynamical localization transition and engineer exotic topological phases of matter [38–41]. There has also been a rising interest with temporally modulated interactions to design synthetic gauge fields or modify the transport properties [42–46]. In the same spirit, we consider here bosonic particles hopping in a one-dimensional lattice (with unit lattice constant $a = 1$) comprising $N$ sites (with periodic boundary conditions) where both the onsite potential and the mean-field (repulsive) interaction terms are periodically-kicked with the same time sequence (time period $T = 1$ set to unity, driving frequency $\Omega = 2\pi$). Setting $\hbar = 1$, the dynamics of our system is then governed by the following time-dependent Gross-Pitaevskii (GP) lattice equation

$$i\partial_t\psi_x = -\frac{J}{2}(\psi_{x+1} + \psi_{x-1}) + \text{III}(t)(V_x + \tilde{g}N_a|\psi_x|^2)\psi_x, \quad (1)$$

where $\text{III}(t) = \sum_{n \in \mathbb{Z}} \delta(t - n)$ is the Dirac comb, $x$ labels the lattice sites, $J > 0$ is the nearest-neighbor hopping amplitude, $V_x$ the onsite potential, $N_a$ the number of particles and $\tilde{g} > 0$ the two-body interaction strength. The wave function is normalized to $\sum_{x=1}^N |\psi_x(t)|^2 = 1$. In the following, we will consider $V_x = 0$ (clean case), $V_x = W\delta_{x,x_0}$

(single-site impurity with strength $W$ located at some given site $x_0$), and $V_x \in [-W/2, W/2]$ (site-uncorrelated uniformly distributed random potential). For future purposes, we define the renormalised nonlinear interaction strength $g = \tilde{g}N_a/N$. This model system could possibly be realized with cold atoms hopping in an optical lattice in the tight-binding regime. An additional optical potential (disordered or not) would be periodically flashed on the atoms. The interaction could also possibly be repeatedly switched on and off using a Feshbach resonance [47].

Noteworthy, the system described by Eq. (1) is formally equivalent to a nonlinear quantum kicked rotor (QKR) [48]. The QKR is a driven system exhibiting a rich phenomenology, from quantum chaos [49, 50] to topological effects [51, 52]. In particular, the QKR displays dynamical localisation [53], a phenomenon analogous to Anderson localisation [54] but in momentum space [55]. This system has been successfully implemented with cold atoms, see e.g. [56–58]. A variant of the QKR includes mean field GP interactions [48, 59, 60] and has been used to study the breakdown of dynamical localization [48, 61–63] and more recently prethermalization and wave-condensation [64].

Importantly, we consider here a driving protocol where the potential $V_x$ is turned on adiabatically instead of abruptly. Sudden quench effects have been recently studied with the nonlinear QKR model in the regime of strong nonlinearities. It has been found that the system shows interesting prethermal properties before reaching an infinite-temperature thermal behavior at large times [65]. By contrast, we focus on an adiabatic driving protocol where an additional controlling parameter $a_t$ is coupled to $V_x$, whose strength is slowly ramped up as $a_t = \tanh(t/\tau)$ with $t \in \mathbb{Z}$ and $\tau = 10^3$ unless specified otherwise. We will consider the subsequent dynamics at times $t \gg \tau$ ($a_t \sim 1$) where the onsite potential has reached its stationary value. Defining $\psi_x^{n\pm} = \psi_x(t = n+0^{\pm})$ and $\psi_k^{n\pm} = \psi_k(t = n+0^{\pm})$ with $\psi_k^n = \sum_x \psi_x^n e^{ikx}$ where $k \in [-\pi, \pi)$ is the quasi-momentum, the dynamics of our system is obtained by iterating the following nonlinear map,

$$\begin{aligned} \psi_k^{n+1-} &= e^{iJ \cos k} \psi_k^{n+}, \\ \psi_x^{n+1+} &= e^{-ia_n V_x} e^{-ig|\psi_x^{n+1-}|^2} \psi_x^{n+1-}. \end{aligned} \quad (2)$$

At this stage, it is interesting to comment on the connection of our model with the well-known conservative case. Our model Eq. (2) is explicitly time-dependent, and we study its dynamics at stroboscopic times. This gives rise to the existence of Floquet quasi-energy bands. Like in [64], we will consider a situation where these quasi-energy bands are well separated from each other so that their coupling, due to interactions, is weak. This implies constraints on the parameters of the model which we will discuss later. In this regime, it is also possible to work in the low quasi-energy sector, i.e. at the edge of a quasi-energy band. This is achieved by considering initial plane wave states with sufficiently small quasi-momenta, see [64]. Moreover, Floquet heating due to the interplay
between periodic driving and interactions must be considered in our model. Indeed, since energy is not conserved, the system is generally expected to evolve toward a featureless state which maximizes the entropy akin to infinite temperature at long times.

**Superfluid properties in the clean case.** — We will now show that our model shares some important formal features with the usual conservative superfluids, in particular the concepts of sound velocity and healing length. Adapting Bogoliubov theory to our periodically-driven system, we derive a low-energy excitation spectrum with a linear dispersion relation at low momenta that supports sound waves. Along this analysis, a threshold on the hopping amplitude $J$ is obtained to avoid fast Floquet heating [64, 66].

Superfluidity can be characterized by a stability analysis of the initial plane-wave mode, $\psi_x(t=0) = \frac{1}{\sqrt{N}}e^{-ik_0x}$. In the clean situation, the onsite potential is zero and the plane wave $\psi_x^0(t) = \frac{1}{\sqrt{N}}e^{-i\phi(t)}e^{-ik_0x}$ is an exact solution when $\dot{\phi}(t) = g\text{III}(t) - J\cos k_0$. In order to obtain the low-energy excitation spectrum, we perform a linear stability analysis with the ansatz $\psi_x(t) = \psi_x^0(t)(1 + \delta\psi_x(t))$. We decompose the perturbation $\delta\psi_x(t) = \sum_q u(q, t)e^{-iqx} + v^*(q, t)e^{iqx}$ into different plane-wave modes labeled by $q \in [-\pi, \pi)$. Linearizing Eq. 1, we obtain the time-dependent Bogoliubov-de Gennes equation for the perturbation:

$$i\partial_t \begin{pmatrix} u \\ v \end{pmatrix} = \mathcal{M}(q, t) \begin{pmatrix} u \\ v \end{pmatrix}. \quad (3)$$

Introducing $\lambda = 2J \sin(\frac{q}{2}) \sin(\frac{q}{2} + k_0)$ and the Pauli matrices, we have $\mathcal{M}(q, t) = \lambda \sigma_z + g\text{III}(t) (\sigma_z + i\sigma_y)$. Since the operator $\mathcal{M}(q, t)$ is periodic in time, we consider the one-period evolution operator associated with it, aka the Floquet operator, given by the time-ordered integration $U(q) = \mathcal{T}e^{-i\int_0^1 dt \mathcal{M}(q,t)}$, that is

$$U(q) = e^{-i\lambda\sigma_z} e^{-ig(\sigma_z + i\sigma_y)}. \quad (4)$$

We next employ the Baker-Campbell-Hausdorff formula [67] to approximate the effective Floquet Hamiltonian for $U(q)$ order by order in both $g$ and $\lambda$, i.e. $H_F(q) = -i\ln U(q) \approx H_F^{(1)}(q) + H_F^{(2)}(q) + \dots$. At first order in $g$ and $\lambda$, we find

$$H_F^{(1)}(q) = (\lambda + g) \sigma_z + ig \sigma_y \quad (5)$$

with eigenvalues $\omega_{\pm}^{(1)}(q) = \pm\omega(q)$ where $\omega(q) = \sqrt{(\lambda + g)^2 - g^2}$. At $k_0 = 0$, we get:

$$\omega(q) = 2\sqrt{gJ} \left| \sin \frac{q}{2} \right| \sqrt{\frac{J}{g} \sin^2 \frac{q}{2} + 1}. \quad (6)$$

At small momenta, $\sin(q/2) \sim q/2$ and this excitation spectrum adopts a form similar to the Bogoliubov spectrum for time-independent systems with a quadratic kinetic energy term with particle mass $m \sim 1/J$. In the limit

<figure><img src="image_2.png" /><figcaption>Fig. 1: Stroboscopic temporal dynamics of the population $\rho(k_0, t)$ of the initial mode for different initial momenta $k_0$ (Mach number $\mathcal{M} = k_0\xi$ where $\xi$ is the healing length). System parameters values are $J = 1.3$, $g = 0.1$ ($\xi \approx 3.6$, $c \approx 0.36$), $\Lambda = W/c = 0.41$ and $N = 1024$. This population remains almost unchanged when $k_0$ is below a certain critical value while a rapid drop of the population takes place when $k_0$ is above that critical value.</figcaption></figure>

$q \to 0$, we recover a linear dispersion relation $\omega(q) \sim cq$ where $c = \sqrt{gJ}$ plays the role of a sound velocity. We also see from the square-root term in Eq. (6) that $\xi = \sqrt{J/g}$ is a length scale playing the same role as the healing length in usual superfluid systems.

At this level of approximation, the excitation dynamics of our periodically driven system is identical to that for a time-independent system described by $H_F^{(1)}(q)$. In particular, Floquet heating is totally absent. It will arise when taking into account higher order terms in the expansion of $H_F(q)$. As shown in the SM, it is possible to exactly diagonalise $U(q)$ in Eq. (4) and look for dynamical instabilities (eigenvalues of the Floquet operator outside the unit disk in the complex plane). At $k_0 = 0$, the first unstable mode is $q = \pi$ and appears when $J + g > \pi/2$ which defines a threshold for Floquet heating. Actually, this criterion has a simple interpretation when considering the relation between the driving frequency $\Omega = 2\pi$ and the effective bandwidth $E_b$ of our model [15, 64, 68]. The single-particle bandwidth is given by $E_b = 2J$, and the correction from non-linearity results in $E_b = 2(J+g)$ for small $g$. An estimate for suppressing the direct inter Floquet band transition is achieved by setting $E_b = 2(J+g) < \pi$ to avoid fast Floquet heating. In the following, we will set the hopping amplitude $J < (\pi/2 - g)$ and $g = 0.1$.

**Superfluid flow across a single impurity.** —

Here, we investigate further the analogy between our periodically-driven system and a conservative superfluid by studying the fate of an initial plane wave in the presence of a single impurity of strength $W$ located at some given site $x_0$, $V_x = W\delta_{xx_0}$, that is ramped up adiabatically to avoid sudden quench effects. In the following, since we work with periodic boundary conditions, we set $x_0 = N/2$ without any loss of generality. We show that, below a certain critical velocity threshold, the system maintains its superfluid properties, i.e. can flow through the impurity
without any dissipation or back scattering.

In traditional superfluids, the nature of the flow interacting with a single localized impurity is governed by two independent dimensionless parameters (see [36, 37] and SM). The first one is the ratio between the velocity of the fluid and the sound velocity, aka the Mach number $M$. The second parameter is the dimensionless impurity strength $\Lambda$ (see SM). In this ($M, \Lambda$) parameter space, there is a critical line that separates the superfluid phase existing at low Mach numbers from a non-superfluid phase existing at higher Mach numbers. In the superfluid regime, the impurity cannot induce any excitation in the fluid and the flow is stationary. It therefore lasts forever. Above the critical line, the flow is slowed down by energy transfers from the coherent motion to internal excitations and enters a weak turbulent regime characterized by erratic and complex dynamics.

Such a critical line also exists in our periodically-driven system. In our case, since $J \sim 1/m$ and $\hbar = 1$, the velocity corresponds to $Jk_0$ and the Mach number simply reads $M = Jk_0/c = k_0\xi$ whereas the effective impurity strength is given by $\Lambda = W\xi/J = W/c$, see Table 1 in the SM. Concretely, we initialize the system in the single plane wave mode $\psi_x(t=0) = \frac{1}{\sqrt{N}}e^{-ik_0x}$ with quasi-momentum $k_0$. We then iterate the quantum map Eq. 2 with $V_x = W\delta_{x,x_0}$ for a given $W$. The nature of the flow is then established by looking at the evolution of the population of the initial mode $\rho(k_0, t) = |\psi_{k_0}(t)|^2$ over a very long time. Figure 1 illustrates our numerical experiment for a given set of parameters. Remarkably, $\rho(k_0, t \gg \tau)$ shows two different behaviors. For initial momenta $k_0$ below a certain critical threshold $k_c$, the population of the initial plane wave remains very close to unity, meaning that a superflow is preserved. Above $k_c$, $\rho(k_0, t)$ drops to zero after a certain time, meaning that the initial state is only marginally populated compared to other modes and the superflow is lost.

Next, we repeat the numerical simulation for different impurity strength $W$ to quantitatively extract the critical line in the ($M, \Lambda$) plane, see Fig.2. The agreement between our data, obtained for a Floquet system, and the analytical prediction for this critical line in a conservative quantum fluid [36, 37], namely

$$ \Lambda_c(M) = \frac{\sqrt{2}}{4M} \sqrt{-8M^4 - 20M^2 + 1 + (1 + 8M^2)^{3/2}}, \quad (7) $$

is truly remarkable.

**Superfluid to prethermal transition.** — Above the critical velocity, the initial state is unstable and the flow enters a regime analogous to the weak turbulent regime in time-independent systems [36, 37]. Following [64], we expect our system to enter a long-lived prethermal state before Floquet heating drives the system into an infinite temperature state after some boiling time $\tau_{\text{boil}}$. A natural question that then arises is whether, below the critical line, $\tau_{\text{boil}}$ is finite (albeit larger than the observation time

<figure><img src="image_3.png" /><figcaption>Fig. 2: Phase diagram of the fluid in the presence of a single-site impurity. In our case, the Mach number and renormalized impurity strength are $M = k_0\xi$ and $\Lambda = W/c$, where $\xi = \sqrt{J/g}$ is the healing length and $c = \sqrt{gJ}$ is the sound velocity. The critical line (green triangles) delineating the superfluid state from a prethermal state is numerically determined from the time evolution of the population in the initial mode. Our data match remarkably well the analytical prediction Eq. (7), obtained for conservative quantum fluids (red solid line). In each region, the color map quantifies the boiling time $\tau_{\text{boil}}$. Parameters values are set to $J = 1.3$, $g = 0.1$ and $N = 1024$.</figcaption></figure>

considered in Figure 1), in which case, we would have a metastable superfluid state, or infinite, in which case we would have a true superfluid state. In this Section, we show that superfluidity prevents thermalisation at all times in our periodically-driven nonlinear system.

The boiling process induced by Floquet heating can be conveniently characterised by the temporal evolution of the variance $\sigma_k^2(t) = \langle k^2 \rangle - \langle k \rangle^2$ of the quasi-momentum distribution [64]. In the prethermal regime $t \ll \tau_{\text{boil}}$, $\sigma_k^2(t)$ shows a plateau while at $t \gg \tau_{\text{boil}}$ it saturates to the value $\pi^2/3$ obtained for a uniform momentum distribution, a signature of an infinite temperature state. The transition from these two behaviors is sharp and the time at which this happens defines the boiling time. Figures 3(b) and (c) present two examples of simulations in, respectively, the superfluid regime and the prethermal regime. In the first case, the variance remains small and no boiling time could be defined up to the maximal integration time, i.e. $10^8$ periods. On the contrary, inset (c) demonstrates the existence of a finite boiling time although it is exponentially large. Note that, though the variance remains small before $\tau_{\text{boil}}$ in both insets (b) and (c), the underlying dynamical behaviors are different, see SM.

We now present one of our most important results. In Fig.3a, we plot the evolution of $\tau_{\text{boil}}$ as a function of the Mach number $M = k_0\xi$ at fixed dimensionless impurity
<figure><img src="image_4.png" /><figcaption>Fig. 3: Temporal behavior of the variance $\sigma_k^2(t)$ of the momentum distribution and boiling time obtained by iterating the nonlinear map Eq. 2 for different initial momenta $k_0$. Parameter values are $J = 1.3$, $g = 0.1$, $\Lambda = W/c = 0.39$ and $N = 1024$. (a) Boiling time $\tau_{\text{boil}}$ as a function of the Mach number $M = k_0\xi$. The vertical black line marks the critical Mach number $M_c$ inferred from the theoretical prediction Eq. (7). (b) The variance for $k_0 < k_c$ remains small and does not change (up to $10^8$ periods). (c) The variance for $k_0 > k_c$ remains small for an exponentially long time before it abruptly jumps and saturates to $\pi^2/3$. The blue vertical dashed line defines $\tau_{\text{boil}}$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4: Probability distribution function (PDF) $P_N(M_c)$ of the critical Mach number $M_c = Jk_c/c$ for 4 different disorder strengths $\Lambda = W/c$ at fixed system size $N = 1024$, interaction strength $g = 0.1$ and hopping parameter $J = 1$. Symbols are the histograms extracted from numerical simulations using $10^4$ disorder configurations: $\Lambda = 0.05$ (black dots), $\Lambda = 0.15$ (blue dots), $\Lambda = 0.3$ (green dots) and $\Lambda = 0.6$ (red dots). The shaded curves are the corresponding analytical predictions from Eq. (8).</figcaption></figure>

strength $\Lambda = W/c$. Strikingly, we find that $\tau_{\text{boil}}$ jumps discontinuously when $M$ crosses the critical value $M_c$ for the superfluid regime previously determined in Figure 2. In the prethermal phase, $\tau_{\text{boil}}$ is always finite while in the superfluid phase, thermalisation of the system does not take place up to the longest times considered ($t = 10^8$). In Fig.3a, the jump in size of $\tau_{\text{boil}}$ just below and above $M_c$ is of nearly two orders of magnitude. This sharp transition at $M_c$ strongly suggests that the superfluid is immune to Floquet heating. This is reflected in Fig.2 where the color map in the phase diagram represents $\tau_{\text{boil}}$ as a function of $M$ and $\Lambda$, as obtained from iterating the Floquet quantum map up to $10^7$ periods. The excellent agreement between the critical line below which superfluidity is maintained and the critical line below which $\tau_{\text{boil}}$ diverges leads us to conclude that a superfluid state below its critical velocity is stable in our periodically-driven system, i.e. robust against thermalisation due to either interaction with a localized impurity or to Floquet heating, or to their combined effects.

**Disordered case.** — We now address the dynamical behavior of our system in the presence of a site-uncorrelated disordered potential with uniform distribution $V_x \in [-W/2, W/2]$. As expected from [69], our numerical investigations (data not shown) indicate that the phenomenology observed in Fig.2 and Fig.3 with a single impurity remains entirely correct in the presence of a disordered potential. The only difference is that the critical line obtained in Fig.2 now depends on the configuration of disorder considered and $M_c$ becomes a random vari-

able. We therefore discuss here the statistical properties of this random critical momentum $M_c$ as a function of the different parameters.

To compute the distribution of $M_c$, we use arguments based on screening and extreme value statistics [69]. First, we recall that $\xi$ is the minimal length scale for density fluctuations. Therefore, spatial details of the disorder potential on scales below $\xi$ cannot be resolved by the system. We therefore renormalise the disorder through a coarse graining procedure over a scale of a few $\xi$. Then, we isolate the maximum value of the renormalised disorder and apply criterion Eq.(7) to obtain the critical velocity. Using extreme value statistics, it is possible to compute the full distribution of the critical velocity as it is explained in the SM. This approach gives a rather accurate description of the statistical properties of $M_c$ as demonstrated below. Our central result is the cumulative distribution of the critical $M_c$ which is given by the following formula

$$ \Phi_N(M_c) = 1 - \left[ \frac{1 + \text{erf} \left( \frac{\Lambda_c(M_c)}{\Lambda_0} \right)}{2} \right]^{\frac{N}{\alpha\xi}} \quad (8) $$

with $M_c = k_c\xi$, $\Lambda_0 = \sqrt{\frac{\alpha\xi}{6}}$, $\Lambda = \sqrt{\frac{\alpha\xi}{6}}\frac{W}{c}$ where $\alpha = 2.78$ (see SM). Then the probability distribution function (PDF) of the critical $M_c$ is simply the derivative of $\Phi_N$, $P_N(M) = \Phi'_N(M)$. Note that this simple expression contains nontrivial scaling properties of $M_c$ with respect to the system size $N$, the disorder strength $W$, the hopping amplitude $J$ and the non-linearity strength $g$ via $c$ and $\xi$.
<figure><img src="image_6.png" /><figcaption>Fig. 5: Median critical Mach number $\overline{M_c} = \overline{k_c}\xi$ as a function of $\Lambda = W/c$ for system size $N = 1024$, interaction strength $g = 0.1$ and hopping parameter $J = 1$. Black dots: Data obtained from numerical simulations using $10^4$ disorder configurations. Thick blue line: Extreme value statistics prediction. Dashed line: Theory prediction for a single impurity of strength $W$.</figcaption></figure>

Fig. 4 shows the comparison between numerically extracted histograms of the critical $M_c$ and our analytical predictions for different values of the disorder strength $W$. Our statistical model remarkably reproduces not only the typical critical $M_c$ and its fluctuations but also the full probability distribution. In order to quantify the typical value of the critical $M_c$, we have computed its median $\overline{M_c}$, defined as $\Phi_N(\overline{M_c}) = 1/2$, as a function of the disorder strength $W$ and the results are shown in Fig. 5. One can see that it is substantially below the critical Mach number predicted for a single impurity but very well captured by the effective disorder computed from extreme value arguments. In addition, we have also checked the scaling with the system size $N$ (see SM) and interaction strength $g$ (data not shown) which give the same level of agreement. Altogether, these results demonstrate that the superfluid to prethermal transition in the disordered case can be understood and accurately described by a coarse-graining procedure. In other words, the single-site impurity case can be seen as an effective description of the uniformly distributed disorder case provided its strength is appropriately renormalised.

Last, we have also addressed the Floquet heating instability with disorder. We find that the distribution of the critical $M_c$ determined above predicts very well the onset of a diverging $\tau_{\text{boil}}$ (see SM). For a superfluid moving through disorder with initial $M$ below its renormalised critical $M_c$, thermalisation is not seen to occur under periodic driving.

**Conclusion.** — In this letter, we have shown that superfluidity can be used to prevent thermalisation in a non-linear Floquet system. While fast driving is a well-known mechanism to postpone Floquet heating to exponentially

long times, we have demonstrated a model where the boiling time diverges at the prethermal to superfluid transition. This promotes superfluidity as a new mechanism for ergodicity breaking in Floquet systems.

This opens interesting avenues to study Floquet topological superfluids [70] and dynamical phase transitions [71] such as the recently proposed dynamical Berezinskii-Kosterlitz-Thouless phase transition [72].

***

We thank P. E. Larré for fruitful discussions. This study has been supported by the French National Research Agency (ANR) under projects COCOA ANR-17-CE30-0024, MANYLOK ANR-18-CE30-0017 and GLADYS ANR-19-CE30-0013, the EUR grant NanoX No. ANR-17-EURE-0009 in the framework of the “Programme des Investissements d’Avenir”, and by the Singapore Ministry of Education Academic Research Fund Tier I (WBS No. R-144-000-437-114). Computational resources were provided by the facilities of Calcul en Midi-Pyrénées (CALMIP) and the National Supercomputing Centre (NSCC), Singapore.

## REFERENCES

[1] NANDKISHORE R. and HUSE D. A., *Annual Review of Condensed Matter Physics*, **6** (2015) 15.

[2] ABANIN D. A., ALTMAN E., BLOCH I. and SERBYN M., *Rev. Mod. Phys.*, **91** (2019) 021001.

[3] SERBYN M., ABANIN D. A. and PAPIĆ Z., *Nature Physics*, **17** (2021) 675.

[4] MOUDGALYA S., BERNEVIG B. A. and REGNAULT N., *Rep. Prog. Phys.*, **85** (2022) 086501.

[5] WEITENBERG C. and SIMONET J., *Nature Physics*, **17** (2021) 1342.

[6] ELSE D. V., BAUER B. and NAYAK C., *Phys. Rev. X*, **7** (2017) 011026.

[7] D’ALESSIO L. and RIGOL M., *Phys. Rev. X*, **4** (2014) 041048.

[8] PONTE P., CHANDRAN A., PAPIĆ Z. and ABANIN D. A., *Annals of Physics*, **353** (2015) 196.

[9] REITTER M., NÄGER J., WINTERSPERGER K., STRÄTER C., BLOCH I., ECKARDT A. and SCHNEIDER U., *Phys. Rev. Lett.*, **119** (2017) 200402.

[10] PONTE P., PAPIĆ Z., HUVENEERS F. M. C. and ABANIN D. A., *Phys. Rev. Lett.*, **114** (2015) 140401.

[11] LAZARIDES A., DAS A. and MOESSNER R., *Phys. Rev. Lett.*, **115** (2015) 030402.

[12] KHEMANI V., LAZARIDES A., MOESSNER R. and SONDHI S. L., *Phys. Rev. Lett.*, **116** (2016) 250401.

[13] ELSE D. V., BAUER B. and NAYAK C., *Phys. Rev. Lett.*, **117** (2016) 090402.

[14] ABANIN D. A., DE ROECK W. and HUVENEERS F. M. C., *Phys. Rev. Lett.*, **115** (2015) 256803.

[15] BUKOV M., GOPALAKRISHNAN S., KNAP M. and DEMLER E., *Phys. Rev. Lett.*, **115** (2015) 205301.

[16] MORI T., KUWAHARA T. and SAITO K., *Phys. Rev. Lett.*, **116** (2016) 120401.

[17] MALLAYYA K., RIGOL M. and DE ROECK W., *Phys. Rev. X*, **9** (2019) 021027.
[18] RUBIO-ABADAL A et al., Phys. Rev. X, **10** (2020) 021044.

[19] PENG P., YIN C., HUANG X., RAMANATHAN C. and CAPPELLARO P., *Nature Physics*, **17** (2021) 444.

[20] SHKEDROV C., MENASHES M., NESS G., VAINBAUM A., ALTMAN E. and SAGI Y., Phys. Rev. X, **12** (2022) 011041.

[21] MORI T., Phys. Rev. B, **98** (2018) 104303.

[22] RAJAK A., DANA I. and DALLA TORRE E. G., Phys. Rev. B, **100** (2019) 100302.

[23] HOWELL O., WEINBERG P., SELS D., POLKOVNIKOV A. and BUKOV M., Phys. Rev. Lett., **122** (2019) 010602.

[24] HODSON W. and JARZYNSKI C., Phys. Rev. Research, **3** (2021) 013219.

[25] LEGGETT A. J., Rev. Mod. Phys., **71** (1999) S318.

[26] LANDAU L. D., J. Phys. U.S.S.R., **5** (1941) 71.

[27] NAZARENKO S., *Wave turbulence* Vol. 825 (Springer Science & Business Media) 2011.

[28] ALLEN J. F. and MISENER A. D., Nature, **141** (1938) 3558.

[29] KAPITZA P., Nature, **141** (1938) 3558.

[30] OSHEROFF D. D., RICHARDSON R. C. and LEE D. M., Phys. Rev. Lett., **28** (1972) 885.

[31] RAMAN C. et al., Phys. Rev. Lett., **83** (1999) 2502.

[32] AMO A. et al., Nature Physics, **5** (2009) 805.

[33] MICHEL C. et al., Nat. Commun., **9** (2018) 2108.

[34] BUCHHOLD M., HEYL M. and DIEHL, S., Phys. Rev. A, **94** (2016) 013601.

[35] SCOQUART T., DELANDE D. and CHERRORET N., Phys. Rev. A, **106** (2022) L021301.

[36] HAKIM V., Phys. Rev. E, **55** (1997) 2835.

[37] LEBOEUF P. and PAVLOFF N., Phys. Rev. A, **64** (2001) 033602.

[38] ARTUSO R., BORGONOVI F., GUARNERI I., REBUZZINI L. and CASATI G., Phys. Rev. Lett., **69** (1992) 3302.

[39] PINQUAN Q., CHUANHAO Y. and SHU C., Phys. Rev. B, **90** (2014) 054303.

[40] CADEZ T., MONDAINI R. and SACRAMENTO P. D., Phys. Rev. B, **96** (2017) 144301.

[41] BOMANTARA R. W. and GONG J., Phys. Rev. Lett., **120** (2018) 230405.

[42] WANG J. and GONG J., Phys. Rev. Lett., **102** (2009) 244102.

[43] GONG J., MORALES-MOLINA L. and HĂNGGI P., Phys. Rev. Lett., **103** (2009) 133002.

[44] GRESCHNER S., SUN G., POLETTI D. and SANTOS L., Phys. Rev. Lett., **113** (2014) 215303.

[45] MEINERT F., MARK M. J., LAUBER K., DALEY A. J. and NAGERL H. C., Phys. Rev. Lett., **116** (2016) 205301.

[46] DUVAL C., DELANDE D. and CHERRORET N., Phys. Rev. A, **105** (2022) 033309.

[47] CHIN C., GRIMM R., JULIENNE P. and TIESINGA E., Rev. Mod. Phys., **82** (2010) 1225.

[48] SHEPELYANSKY D., Phys. Rev. Lett., **70** (1993) 1787.

[49] IZRAILEV F. M., Phys. Rep., **196** (1990) 299.

[50] HAAKE F., *Quantum Signatures of Chaos* Springer-Verlag, Berlin, Heidelberg 2006.

[51] DAHLHAUS J. P., EDGE J. M., TWORZYDŁO J. and BEENAKKER C. W., Phys. Rev. B, **84** (2011) 115133.

[52] HO D. Y. H. and GONG J., Phys. Rev. Lett., **109** (2012) 010601.

[53] CASATI G., CHIRIKOV B. V., IZRAELEV F. M. and FORD J., *Stochastic behavior in classical and quantum Hamiltonian systems* Springer 1979.

[54] ABRAHAMS E., ED., *50 Years of Anderson Localization* World Scientific 2010.

[55] GREMPEL D. R., PRANGE R. E. and FISHMAN S., Phys. Rev. A, **29** (1984) 1639.

[56] MOORE F.L., ROBINSON J.C., BHARUCHA C.F., WILLIAMS P.E. and RAIZEN M.G., Phys. Rev. Lett., **73** (1994) 2974.

[57] CHABÉ J., LEMARIÉ G., GRÉMAUD B., DELANDE D., SZRIFTGISER P. and GARREAU, J. C., Phys. Rev. Lett., **101** (2008) 255702.

[58] LEMARIÉ G., CHABÉ J., SZRIFTGISER P., GARREAU J. C., GRÉMAUD B. and DELANDE D., Phys. Rev. A, **80** (2009) 043626.

[59] GLIGORIĆ G., BODYFELT J. D. and FLACH S., Eur. Phys. Lett., **96** (2011) 30004.

[60] LELLOUCH S., RANÇON A., DE BIÈVRE S., DELANDE D. and GARREAU J. C., Phys. Rev. A, **101** (2020) 043624.

[61] CHERRORET N., VERMERSCH B., GARREAU J. C. and DELANDE D., Phys. Rev. Lett., **112** (2014) 170603.

[62] CAO A. et al., arXiv, (2021) 2106.09698.

[63] TOH J. H. S. et al., arXiv, (2021) 2106.13773.

[64] HALDAR P., MU S., GEORGEOT B., GONG J., MINIATURA C. and LEMARIÉ G., arXiv, (2021) 2109.14347.

[65] MARTINEZ M., LARRÉ P. E., DELANDE D. and CHERRORET N., Phys. Rev. A, **106** (2022) 043304.

[66] LELLOUCH S., BUKOV M., DEMLER E. and GOLDMAN N., Phys. Rev. X, **7** (2017) 021015.

[67] ACHILLES R. and BONFIGLIOLI A., *Archive for History of Exact Sciences*, **66** (2012) 295.

[68] GENSKE M. and ROSCH A., Phys. Rev. A, **92** (2015) 062108.

[69] ALBERT M., PAUL T., PAVLOFF N. and LEBOEUF P., Phys. Rev. A, **82** (2010) 011602.

[70] FOSTER M. S., GURARIE V., DZERO M. and YUZBASHYAN E. A., Phys. Rev. Lett., **113** (2014) 076403.

[71] HEYL M., Rep. Prog. Phys., **81** (2018) 054001.

[72] SCOQUART T., DELANDE D. and CHERRORET N., Phys. Rev. A, **106** (2022) L021301.
# Supplementary material for "Superfluidity vs prethermalisation in a nonlinear Floquet system"

S. MU¹, N. MACÉ³, J. GONG¹,², C. MINIATURA²,⁴, G. LEMARIÉ²,³,⁴ and M. ALBERT⁵

¹ Department of Physics, National University of Singapore, Singapore 117542, Singapore

² Centre for Quantum Technologies, National University of Singapore, Singapore 117543, Singapore

³ Laboratoire de Physique Théorique, IRSAMC, Université de Toulouse, CNRS, UPS, France

⁴ MajuLab, International Joint Research Unit UMI 3654, CNRS, Université Côte d'Azur, Sorbonne Université, National University of Singapore, Nanyang Technological University, Singapore

⁵ Université Côte d'Azur, CNRS, INPHYNI, France

## Abstract -*** Missing author ***

This supplementary contains the following material arranged by sections:

1. We present the correspondence between the superfluid states obtained in our periodically-kicked Gross-Pitaevskii model and in the usual time-independent Gross-Pitaevskii equation.

2. We detail the derivation of the Bogoliubov expansion for our model in the clean case ($V_x = 0$).

3. We detail the derivation of the effective impurity strength in the disordered Floquet model.

**Mapping to the Gross-Pitaevskii equation.** — The one-dimensional Gross-Pitaevskii equation describing the time evolution of the macroscopic wave function $\Phi(x, t)$ of a system of identical bosons of mass $m$ interacting with a single impurity of strength $W$ located at the origin reads:

$$i\hbar\partial_t\Phi(x,t) = -\frac{\hbar^2}{2m}\frac{\partial^2}{\partial x^2}\Phi(x,t) + g|\Phi(x,t)|^2\Phi(x,t) + W\delta(x)\Phi(x,t). \quad (1)$$

The normalisation condition reads $\int dx |\Phi(x,t)|^2 = N_0$, where $N_0$ is the total number of atoms, and we call $n_0$ the atomic number density.

Setting $\Phi(x, t) = \exp(-i\mu t/\hbar)\phi(x)$, we now look for a stationary solution with the following boundary condition:

$$\phi(x) \underset{|x| \to +\infty}{\to} \sqrt{n_0} e^{imv_0x/\hbar}. \quad (2)$$

It describes a macroscopic stationary wave function with asymptotic density $n_0$ and velocity $v_0$. Since $\phi(x)$ solves

$$\mu\phi(x) = -\frac{\hbar^2}{2m}\phi''(x) + g|\phi(x)|^2\phi(x) + W\delta(x)\phi(x), \quad (3)$$

where $\phi''$ denotes the second space derivative, we see that the boundary condition imposes

$$\mu = gn_0 + \frac{1}{2}mv_0^2. \quad (4)$$

We introduce the sound velocity $c$ and the healing length $\xi$ according to the standard formulae [1]

$$mc^2 = gn_0, \quad \xi = \frac{\hbar}{\sqrt{mgn_0}} = \frac{\hbar}{mc}. \quad (5)$$

Using the rescaled variables $z = x/\xi$ and $\psi = \phi/\sqrt{n_0}$, Eq. (3) takes on the following universal form:

$$\left[1 + \frac{1}{2}\mathcal{M}^2\right] \psi(z) = -\frac{1}{2}\psi''(z) + \Lambda\delta(z)\psi(z) + |\psi(z)|^2\psi(z) \quad (6)$$

which depends on two dimensionless parameters, the Mach number $\mathcal{M}$ and the effective impurity strength $\Lambda$ given by:

$$\mathcal{M} = \frac{v_0}{c}, \quad \Lambda = \frac{m}{\hbar^2} W \xi = \frac{W}{\hbar c}. \quad (7)$$

The existence of the superfluid solution $\psi(z)$ is therefore governed by $\mathcal{M}$ and $\Lambda$ only as demonstrated in [2, 3] for instance.

Realizing that the discrete version of the kinetic energy term $-\frac{\hbar^2}{2m}\frac{\partial^2}{\partial x^2}\Phi(x,t)$ writes $-\frac{\hbar^2}{2ma^2}[\Phi_{x+a}(t) + \Phi_{x-a}(t)]$
<table><caption>Table 1: Mapping obtained at low energy between the continuous conservative Gross-Pitaevskii model given by Eq.(1) and the non-linear Floquet lattice model given by Eq.(1) in the main text (with lattice constant $a = 1$ and $\hbar = 1$).</caption><thead><tr><th>Continuous GP<br>Model Eq.(1)</th><th>Kicked GP lattice<br>Model Eq.(1) in the main text</th></tr></thead><tbody><tr><td><i>p</i></td><td><i>k</i></td></tr><tr><td><i>m</i></td><td>1/<i>J</i></td></tr><tr><td><i>gn</i><sub>0</sub></td><td><i>g</i></td></tr><tr><td><i>W</i></td><td><i>W</i></td></tr><tr><td><i>v</i><sub>0</sub></td><td><i>J</i><i>k</i><sub>0</sub></td></tr><tr><td><i>c</i> = &radic;<i>gn</i><sub>0</sub>/<i>m</i></td><td><i>c</i> = &radic;<i>Jg</i></td></tr><tr><td><i>&xi;</i> = <i>&hbar;</i>/<i>mc</i></td><td><i>&xi;</i> = &radic;<i>J</i>/<i>g</i></td></tr><tr><td><i>M</i> = <i>v</i><sub>0</sub>/<i>c</i></td><td><i>M</i> = <i>Jk</i><sub>0</sub>/<i>c</i> = <i>k</i><sub>0</sub><i>&xi;</i></td></tr><tr><td><i>&Lambda;</i> = <i>mW</i><i>&xi;</i>/<i>&hbar;</i><sup>2</sup></td><td><i>&Lambda;</i> = <i>W</i><i>&xi;</i>/<i>J</i> = <i>W</i>/<i>c</i></td></tr></tbody></table>

$2\Phi_x(t)$ on a lattice of spatial period $a$, we see that the hopping amplitude in this discrete version is simply $\hbar^2/(ma^2)$. The same scaling is obviously found from the dispersion relation of $-J\cos(ka)$ obtained in our lattice model in the absence of kicks. For $ka \ll 1$, we get $-J\cos(ka) \approx -J(1-(ka)^2/2)$. With $p = \hbar k$, the effective mass of our discrete model would simply be $\hbar^2/(Ja^2)$. Setting $\hbar = 1$ and $a = 1$, we thus get the mapping $J \leftrightarrow 1/m$ between the continuous GP model given above by Eq.(1) and our kicked GP lattice model given by Eq.(1) in the main text. With this in mind, the formal mapping given in Table 1 immediately follows. Along with the quantum map associated to Eq.(1) of the main text, we choose to spotlight the momentum variables here. As a consequence, hereafter, we express the Mach number as $M = k_0\xi$ in our kicked system. We also choose to express the effective impurity strength as $\Lambda = W/c$. One of the major results of this Letter is that a superfluid state can persist under the Floquet dynamics just like it can persist in continuous conservative systems.

To complement the results on the critical Mach number $M_c$ shown in Fig.2 of the main text, we present here additional results showing the scaling behavior of the critical curve when plotted in the $(M, \Lambda)$ plane. In Fig.1a, we plot the numerically-extracted critical momenta $k_c$ as a function of $W$ for different hopping amplitudes $J$. As one can see, we obtain different non-overlapping curves. However when $k_c$ and $W$ are properly rescaled by $\xi$ and $c$, all these curves collapse onto a single universal critical line, see Fig.1b. The agreement with the analytical prediction $W/c = \Lambda_c(M_c)$, where

$$ \Lambda_c(y) = \frac{\sqrt{2}}{4y} \sqrt{-8y^4 - 20y^2 + 1 + (1 + 8y^2)^{3/2}}, \quad (8) $$

proves just remarkable.

Bogoliubov expansion in the periodically kicked system without impurities. — In the main text, we have seen that a sound velocity $c = \sqrt{gJ}$ and a healing length $\xi = \sqrt{J/g}$ can be defined in our periodically-kicked

<figure><img src="image_7.png" /><figcaption>Fig. 1: (a) Numerically-extracted critical momentum $k_c$ (symbols) as a function of the bare impurity strength $W$ for different hopping amplitudes $J$ at interaction strength $g = 0.1$ and system size $N = 1024$. All curves are different. (b) When $k_c$ and $W$ are rescaled to $M_c = k_0\xi$ and $\Lambda = W/c$, where $\xi$ is the healing length and $c$ the sound velocity, all previous curves collapse onto a single universal critical line. The analytical prediction (solid line) is given by $W/c = \Lambda_c(M_c)$ through Eq.(8). As one can see, the agreement between the data (symbols) and the theory is very good.</figcaption></figure>

system by looking at low-energy excitations on top of a stationary plane-wave solution. Here we provide more details about this derivation and give additional numerical results on the linear dispersion relation and on the dynamical instabilities.

Let us consider a clean superfluid ($V_x = 0$) in a one dimensional lattice under periodic driving. Eq.(1) in the main text then reduces to

$$ i\partial_t \psi_x = -\frac{J}{2}(\psi_{x+1} + \psi_{x-1}) + \tilde{g}N_a \Pi(t) |\psi_x|^2 \psi_x, \quad (9) $$

where $\Pi(t) = \sum_{n \in \mathbb{Z}} \delta(t-n)$. We can readily check that any plane wave state $\psi_x^0(t) = \frac{1}{\sqrt{N}}e^{-i\phi(t)}e^{-ik_0x}$ is an exact solution of Eq.(9) provided $\dot{\phi}(t) = g\Pi(t) - J\cos k_0$ where $g = \tilde{g}N_a/N$. Given this time-dependent solution, the basic idea of the Bogoliubov expansion is to consider a small perturbation on top of it. There are several ways to parameterise the perturbation, and we choose the following form for technical convenience [4], $\psi_x(t) = \psi_x^0(t)(1+\delta\psi_x(t))$. Using this parametrisation and after linearising the nonlinear equation Eq. (9), we can obtain the time-dependent Bogoliubov-de Gennes equations for the perturbation, which take the general form,

$$ i\partial_t \begin{pmatrix} \delta\psi_x(t) \\ \delta\psi_x^*(t) \end{pmatrix} = \mathcal{L}(t) \begin{pmatrix} \delta\psi_x(t) \\ \delta\psi_x^*(t) \end{pmatrix}, \quad (10) $$

where

$$ \mathcal{L}(t) = \begin{pmatrix} \mathcal{F}_x(t) & g\Pi(t)|\psi_x^0|^2 \\ -g\Pi(t)|\psi_x^0|^2 & -\mathcal{F}_x^*(t) \end{pmatrix}. \quad (11) $$

Here, the linear operator $\mathcal{F}_x(t)$ acting on $\delta\psi_x(t)$ gives

$$ \begin{aligned} \mathcal{F}_x(t)\delta\psi_x(t) = & -\frac{J}{2} \frac{\psi_{x+1}^0(t)\delta\psi_{x+1}(t) + \psi_{x-1}^0(t)\delta\psi_{x-1}(t)}{\psi_x^0(t)} \\ & + 2g\Pi(t)\delta\psi_x(t) - \frac{i\partial_t\psi_x^0(t)}{\psi_x^0(t)}\delta\psi_x(t). \end{aligned} \quad (12) $$
For the chosen exact solution, we have $\frac{i\partial_t\psi_x^0(t)}{\psi_x^0(t)} = \dot{\phi}(t)$ and $\psi_{x\pm 1}^0(t)/\psi_x^0(t) = \exp(\mp ik_0)$. Next, we decompose the perturbation into different Fourier modes labeled by $q$,

$$\delta\psi_x(t) = \sum_q u(q,t)e^{-iqx} + v^*(q,t)e^{iqx} \quad (13)$$

with $q \in [-\pi, \pi]$. Substituting the decomposition into Eq. 10, we have

$$i\partial_t \begin{pmatrix} u \\ v \end{pmatrix} = \mathcal{M}(q,t) \begin{pmatrix} u \\ v \end{pmatrix}, \quad (14)$$

where

$$\begin{aligned} \mathcal{M}(q,t) &= \lambda(q) \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix} + g \text{III}(t) \begin{pmatrix} 1 & 1 \\ -1 & -1 \end{pmatrix} \\ &= \lambda(q) \sigma_z + g \text{III}(t) (\sigma_z + i\sigma_y) \end{aligned} \quad (15)$$

with $\lambda(q) = 2J \sin(q/2) \sin(q/2 + k_0)$ and the Pauli matrices $\sigma_y$ and $\sigma_z$.

The 1-period Floquet operator associated to Eq. (14) reads:

$$U(q) = \mathcal{T} e^{-i \int_0^1 dt \mathcal{M}(q,t)} = e^{-i\lambda \sigma_z} e^{-ig(\sigma_z + i\sigma_y)} \quad (16)$$

where $\mathcal{T}$ denotes the time-ordered integration operator. It is a product of 2 non-commuting unitary operators. Next, we note that $(\sigma_z + i\sigma_y)$ is a nilpotent matrix as it squares to zero. As a consequence,

$$e^{-ig(\sigma_z + i\sigma_y)} = 1 - ig(\sigma_z + i\sigma_y), \quad (17)$$

and we easily get:

$$U(q) = \begin{pmatrix} e^{-i\lambda(1-ig)} & -ie^{i\lambda g} \\ ie^{-i\lambda g} & e^{i\lambda(1+ig)}. \end{pmatrix} \quad (18)$$

From the perspective of dynamical instability, we can extract the "Lyapunov" exponents by solving the eigenvalue equation

$$U(q)|\phi_{\pm}(q)\rangle = e^{i\omega_{\pm}(q)}|\phi_{\pm}(q)\rangle \quad (19)$$

with eigenvalues given by

$$e^{i\omega_{\pm}(q)} = (\cos \lambda - g \sin \lambda) \pm \sqrt{g^2 - (g \cos \lambda + \sin \lambda)^2}. \quad (20)$$

The appearance of negative imaginary parts in $\omega_{\pm}(q)$ indicates a dynamical instability, i.e. an exponential growth of the associated mode $q$ at a rate $s_q = -\text{Im}\omega(q)$. Considering different hopping amplitude $J$ with a fixed weak nonlinearity strength $g \ll 1$, we find that the initial solution with $k_0 = 0$ is unstable against perturbations when $J + g > \pi/2$, as presented in Fig.2b. This instability will grow and dominate the dynamics as soon as the unstable $q$ modes get populated. This explains the reason why we have chosen hopping amplitudes satisfying $J + g < \pi/2$ in the main text, for example $J = 1$ or $J = 1.3$ for $g = 0.1$. It has also been shown that in the prethermal plateau

reached by the system, the boiling time is shorter for larger hopping amplitude [5]. It is thus important that our numerical simulation time is larger than the boiling times observed in the prethermal phase. Otherwise, it would be difficult to relate the observed stability of the superfluid state to a genuine robustness against Floquet heating or to a boiling time well beyond the simulation time.

When $\omega(q)$ is real, its associated mode is norm-conserving and stationary. These stroboscopic states are similar to the eigenstates of a time-independent Hamiltonian. To illustrate this point further, we define the Floquet Hamiltonian associated to the first Brillouin zone of the excitation spectrum through:

$$H_F(q) = -i \ln U(q). \quad (21)$$

It effectively describes the dynamics of the excitations in our periodically-driven system. Generally, this operator can be highly non-local since $U(q)$ features the concatenation of 2 non-commuting unitary operators. There are several methods to approximate $H_F(q)$ in the fast-driving frequency regime [6, 7]. Here, we employ instead the Baker-Campbell-Hausdorff (BCH) formula [8]. We focus here on the dynamics of low-energy excitations around an initial plane wave state at $k_0 = 0$, i.e. for $q$ small. In this case, $\lambda \approx Jq^2/2 \ll 1$. We further assume weak interaction strengths $g \ll 1$.

Given two possibly non commuting operators $X$ and $Y$, the BCH formula gives the operator $Z$ satisfying

$$e^Z = e^X e^Y \quad (22)$$

as a formal series in terms of $X$, $Y$ and their iterated commutators:

$$\begin{aligned} Z = & X + Y + \frac{1}{2}[X, Y] \\ & + \frac{1}{12}[X, [X, Y]] - \frac{1}{12}[Y, [X, Y]] + (\cdots) \end{aligned} \quad (23)$$

where $(\cdots)$ indicates terms involving higher commutators of $X$ and $Y$. Note that the series is not necessarily convergent for generic cases. Plugging $X = -i\lambda\sigma_z$ and $Y = -ig(\sigma_z + i\sigma_y)$ in Eq. (23), we can write the first few terms for the effective Floquet Hamiltonian,

$$\begin{aligned} H_F(q) &= H_F^{(1)}(q) + H_F^{(2)}(q) + (\cdots) \\ H_F^{(1)}(q) &= ig\sigma_y + (\lambda + g)\sigma_z \end{aligned} \quad (24)$$

$$H_F^{(2)}(q) = -ig\lambda\sigma_x$$

where $(\cdots)$ denotes higher order terms in $g$ and $\lambda$. In the main text, we presented the effective quasi-energies by considering $H_F^{(1)}(q)$ only. Here we consider the quasi-energies $\pm\omega(q)$ at next order, that is for $H_F^{(1)}(q) + H_F^{(2)}(q)$. We find:

$$\omega(q) = \sqrt{gJ} |q| \sqrt{\frac{J(1-g^2)q^2}{4g} + 1}. \quad (25)$$
<figure><img src="image_8.png" /><figcaption>Fig. 2: (a) Quasi-energy $\omega(q)$ (positive branch) of the Floquet operator Eq.(16) as a function of $q$ (orange line). The linear dispersion relation with sound velocity $c = \sqrt{gJ}$ is indicated by the black dashed line. The green vertical dashed lines mark $1/\xi$ with $\xi = \sqrt{J/g}$. Parameters values are $k_0 = 0, g = 0.1$ and $J = 1$. (b) Dynamical instability diagram in the $(J, q)$ plane. The instability rate $s_q = -\text{Im}\omega(q)$ is quantified with the color map. Instability occurs at $q = \pi$ when $J + g \approx \pi/2$. The black dashed line denotes the instability threshold $J = \pi/2 - g$. Parameters values are $k_0 = 0$ and $g = 0.1$.</figcaption></figure>

We observe that the linear spectrum observed at $q \to 0$ for $H_F^{(1)}(q)$ still persists for $H_F^{(1)}(q) + H_F^{(2)}(q)$ with the same sound velocity $c = \sqrt{gJ}$. We also see that the healing length is slightly renormalised to $\xi = \sqrt{J(1-g^2)/g} \approx \sqrt{J/g}$ since $g \ll 1$. Alternatively, we can numerically diagonalise Eq.(19) to obtain $\omega(q)$ to all orders and extract the sound velocity from the $q \to 0$ linear behavior.

In Fig.2a, we plot $\omega(q)$ as a function of $q$ and a linear dispersion relation is clearly observed for small $q$ up to momenta of order $\xi^{-1}$ with sound velocity $c = \sqrt{gJ}$. Note that the dynamical instability for large hopping amplitude $J$ does not appear in the approximation presented above, since it happens in a regime where the assumption $\lambda \ll 1$ does not hold anymore and originates from higher order terms in the series expansion. As demonstrated in Fig.2b, the instability indeed occurs at $q \approx \pi$, i.e. at $\lambda \approx 2J$. An intuitive picture would be that the excitations, coming in pairs $\pm q$ of momenta, proliferate when the excitation gap is on resonance with the driving frequency.

**Differences between the superfluid and prethermal phases close to the critical line for the single-site impurity case.** — In the phase diagram of Fig.2 of the main text, we concluded that there is a critical line separating a prethermal regime where the Floquet system eventually heats up to infinite temperature from a superfluid regime where no heating ever arises. Here we will present more details on the different dynamical behaviors of the system when it is slightly above or below the critical line.

As seen in Fig.3a, the respective variances $\sigma_k^2(t)$ obtained for an initial mode $k_0$ started slightly above or slightly below the critical line remain identically small for a very long time. As a consequence, they cannot accurately discriminate the two phases, at least until the

<figure><img src="image_9.png" /><figcaption>Fig. 3: Dynamical behavior of our periodically-kicked system in the presence of a single-site impurity. (a) Momentum variance as a function of time for an initial plane wave state $k_0$ just above (prethermal state) and just below (superfluid state) the critical line. (b) Momentum distribution of the system in the prethermal phase at time $t=0$ and at time $t=354813$ slightly before before $\tau_{\text{boil}}$, see dashed vertical line in Panel (a). (c) Same as (b) but in the superfluid phase. (d) Time evolution of the population in mode $k_0$ just above (prethermal state) and just below (superfluid state) the critical line. In all panels, the blue color denotes the superfluid phase and the orange color the prethermal phase. Other parameters values are the same as in Fig. 3 of the main text.</figcaption></figure>

boiling time $\tau_{\text{boil}}$ for the prethermal phase. However, the underlying dynamical behaviors in these two phases are markedly different. This is exemplified in Fig.3b and Fig.3c, where we show the momentum density profiles of our system in the prethermal and superfluid phases at the initial time and at a later time below $\tau_{\text{boil}}$ shown by the vertical dashed line in Panel (a). In the prethermal phase close to the critical line, the momentum distribution obtained for initial modes at $k_0 \neq 0$ shifts towards the "lowest" quasi-energy mode $k=0$ where it grows a momentum component. This phenomenon has been studied in [5] for the disordered case and signals a (classical) wave condensation process at work: the momentum distribution develops a macroscopically occupied mode on top of a thermal distribution described by the Rayleigh-Jeans (RJ) thermal distribution. Importantly, the macroscopically occupied mode is no longer the initial mode. In marked contrast, the momentum distribution obtained in the superfluid phase close to the critical line for initial modes at $k_0 \neq 0$ remains unaffected. This is further illustrated in Fig.3d where we show the population in the initial mode for both phases as a function of time. The initial mode $k_0$ always remains macroscopically occupied in the superfluid phase while it drops much earlier than the boiling time $\tau_{\text{boil}}$ in the prethermal phase.

**Derivation of the effective impurity strength in the disordered case.** — We explain here the renormali-
<figure><img src="image_10.png" /><figcaption>Fig. 4: Probability distribution function (PDF) $P_N(M_c)$ of the critical Mach number $M_c = k_c \xi$ for different system sizes $N$ at constant effective disorder strength $\Lambda = W/c = 0.3$ ($J = 1$ and $g = 0.1$): $N = 512$ (grey dots), $N = 1024$ (blue dots), and $N = 2048$ (red dots). The full curves are analytical predictions extracted from Eq. (30) and symbols are numerical simulations of the disordered Floquet model as explained in the main text. We used $10^4$ disorder realisations.</figcaption></figure>

sation procedure that we use to effectively replace the full disorder potential by a single-site impurity model. The first step is a screening argument. Since the system cannot resolve spatial potential scales smaller than the healing length, we first construct a coarse-grained potential by summing up all the impurities in a window of size $\alpha\xi$, with $\alpha$ a parameter of order unity that will be fixed by fitting a single set of data. Starting from $N$ uncorrelated impurities of strength $V_x \in [-\frac{W}{2}, \frac{W}{2}]$ with box uniform distribution $(\langle V_x \rangle = 0, \langle V_x V_{x'} \rangle = \frac{W^2}{12} \delta_{xx'})$, we define a new renormalised potential

$$U_j = \sum_{x \in \text{box } j} V_x. \qquad (26)$$

Each box contains $M_\alpha = \alpha\xi/a$ impurities, where $a$ is the lattice spacing ($a = 1$ in our work) and the total number of boxes is $B = N/M_\alpha$. The idea is then to look for the dominant contribution that comes from the maximum $U_m = \max\{U_j\}$ of this new disorder coarse-grained potential. At this stage, everything boils down to a problem of extreme value statistics.

We first look at the statistics of the random potential $U_j$. Note that $U_j$ is a sum of independent random variables with the same uniform distribution. It rapidly converges ($M_\alpha \ge 5$ is enough for our purpose) to a normal distribution with zero mean and variance $\sigma^2 = M_\alpha \frac{W^2}{12}$. The statistics of the maximum value $U_m = \max\{\tilde{u}_j\}$ is a standard problem in extreme value statistics, and we briefly explain it here. It is easy to find the cumulative distribution of the maximum, i. e. the probability for $U_m$ to be smaller than a given value $U$. Indeed, for the maximum to be smaller than $U$, all of the $V_x$ contributing to $U_j$ have

to be smaller than $U$. Since the $B$ random variables $U_j$ are independent, the result is just the product of their respective cumulative distributions $F(U)$. We therefore end up with a simple equation for the cumulative distribution $\mathcal{F}_B$ of $U_m$,

$$\mathcal{F}_B(U_m) = [F(U_m)]^B. \qquad (27)$$

Note that for a Gaussian variable $z$ with distribution $p(z) = \exp[-z^2/2\sigma^2]/\sqrt{2\pi\sigma^2}$, we have

$$F(z) = \int_{-\infty}^{z} dyp(y) = \frac{1 + \operatorname{erf}\left(\frac{z}{\sqrt{2\sigma^2}}\right)}{2}. \qquad (28)$$

We now use the crude approximation of replacing the full complexity of the microscopic disorder with a single-site impurity of effective strength $U_m$ and apply the criterion in Eq. (8) above to compute the critical Mach number $M_c = k_c \xi$, namely

$$\frac{U_m}{c} = \Lambda_c(M_c). \qquad (29)$$

This relation being one-to-one we finally obtain the statistical distribution of $M_c$ through the following relation between the cumulative distributions

$$\Phi_N(M_c) = 1 - \mathcal{F}_B(U_m = c\Lambda_c(M_c)). \qquad (30)$$

This yields Eq. (8) in the main text. This simple expression contains however nontrivial scaling properties of $M_c$ (and thus $k_c$) with the system size $N$, the disorder strength $W$, and the hopping amplitude $J$ and interaction strength $g$ (via $c$ and $\xi$). Finally, the parameter $\alpha$ is determined by fitting the theory to a single set of data. We obtained $\alpha = 2.78$ which is indeed a parameter of order unity.

In addition to the data shown in the main text, we demonstrate in Fig. 4 that the scaling of the critical momentum $k_c$ with respect to the system size is also well reproduced by this model. As the system size $N$ is increased, the distribution is shifted to lower momenta. This can be easily understood: By increasing $N$, we increase the number $B$ of boxes and thus also increase the probability to find a large $U_j$ value. Note however, that the Gaussian approximation for the random variable $U_j$ predicts that the critical momentum $k_c$ vanishes in the thermodynamic limit $N \to \infty$ at fixed number of particles per site $N_a/N$. This is actually an artefact of this Gaussian approximation since $U_j$ is bounded by $M_\alpha W/2$. This pathology can be corrected by taking the true distribution of $U_j$ that can be exactly computed. However, we chose to use the Gaussian approximation in order to avoid cumbersome expressions and calculations that are not really important for general illustration purposes.

Last, in Fig. 5, we present additional results on the distribution of the critical Mach number separating a finite boiling time and an "infinite" one for the periodically-kicked nonlinear disordered lattice. Note that by "infiniti" we mean that we were not able to observe any boiling process for an initial momentum $k_0$ slightly below the
<figure><img src="image_11.png" /><figcaption>Fig. 5: Probability distribution function (PDF) $P_N(M_c)$ of the critical Mach number $M_c = k_c \xi$ determined from both the dynamical behaviors of the population $\rho(k_0, t)$ (red dots) and of the variance $\sigma_k^2(t)$ (blue triangles) at constant disorder strength $\Lambda = W/c = 0.3$ ($J = 1.3$ and $g = 0.1$). The blue area is the analytical prediction. We have used here 360 disorder realizations to obtain the PDF.</figcaption></figure>

critical line up to the longest time explored in our nu-
merical simulation ($10^7$ periods). This longest simulation
time was however enough to observe the finite boiling time
$\tau_{boil} \sim 10^6$ for $k_0$ slightly above the critical line. We find
an excellent agreement between the probability distribu-
tion functions (PDF) obtained from the dynamical behav-
iors of the population $\rho(k_0, t)$ and of the variance $\sigma_k^2(t)$.
Specifically, the critical $M_c$ is obtained as the largest value
of $k_0\xi$ for which $\rho(k_0, t)$ never drops below 0.5 and $\sigma_k^2(t)$
never increases above 0.5 within the simulation time $10^7$
in each disorder realization. This confirms again that the
superfluid phase we report in this Letter for our weakly
disordered Floquet model is robust against heating when
$J + g < \pi/2$.

REFERENCES

[1] PETHICK C. J. and SMITH H., *Bose-Einstein condensation in dilute gases* (Cambridge university press) 2008.

[2] HAKIM V., *Phys. Rev. E*, **55** (1997) 2835.

[3] LEBOEUF P. and PAVLOFF N., *Phys. Rev. A*, **64** (2001) 033602.

[4] BUKOV M., GOPALAKRISHNAN S., KNAP M. and DEMLER E., *Phys. Rev. Lett.*, **115** (2015) 205301.

[5] HALDAR P., MU S., GEORGEOT B., GONG J., MINIATURA C. and LEMARIÉ G., arXiv e-prints, (2021) arXiv:2109.14347.

[6] GOLDMAN N. and DALIBARD J., *Phys. Rev. X*, **4** (2014) 031027.

[7] BUKOV M., D’ALESSIO L. and POLKOVNIKOV A., *Advances in Physics*, **64** (2015) 139.

[8] ACHILLES R. and BONFIGLIOLI A., *Archive for History of Exact Sciences*, **66** (2012) 295.