\section{SMBF states and quantum fluctuations}

We show here the emergence of SMBF states as well as the Floquet character of its spectrum of quantum fluctuations for several canonical examples of variational ans\"atze and many-body Hamiltonians.

\subsection{MultiConfiguration Time-Dependent Hartree method}

We consider the following \textit{time-independent} general many-body Hamiltonian in second quantization, valid for both bosons and fermions:
\begin{equation}\label{eq:HamiltonianManyBodyGeneral}
    \hat{H}=\hat{h}+\hat{V}=\sum_{lk} h_{lk}\hat{c}_l^{\dagger}\hat{c}_k+
    \frac{1}{2}\sum_{lkjm} V_{lk,jm}\hat{c}_l^{\dagger}\hat{c}_j^{\dagger}\hat{c}_m\hat{c}_k
\end{equation}
Here, $\hat{h}$ is the single-particle part of the Hamiltonian (typically kinetic energy plus some time-independent external field) and $\hat{V}$ the many-body interacting part, $h_{lk}$ and $V_{lk,jm}$ being their respective matrix elements:
\begin{eqnarray}
    h_{lk}&=&\int\mathrm{d}\mathbf{x}~\left[\phi^*_l(\mathbf{x})\right]_\alpha h_{\alpha\beta}(\mathbf{x})\left[\phi_k(\mathbf{x})\right]_\beta \\
    \nonumber V_{lk,jm}&=&\sum_{\alpha\beta\lambda\mu}\int\mathrm{d}\mathbf{x}\int\mathrm{d}\mathbf{x}'~\left[\phi^*_l(\mathbf{x})\right]_\alpha \left[\phi_k(\mathbf{x})\right]_\beta V_{\alpha\beta,\lambda\mu}(\mathbf{x}-\mathbf{x}'')\left[\phi^*_j(\mathbf{x}')\right]_\lambda \left[\phi_m(\mathbf{x}')\right]_\mu
\end{eqnarray}
where Latin indices label single-particle states and Greek indices label possible discrete (spin, pseudospin, etc.) degrees of freedom of the wave functions.

The MCTDH method [37,38] restricts the dynamics of $N$ particles (bosons or fermions) to a fixed subspace of $M$ time-dependent single-particle states $\{\phi_i(\mathbf{x},t)\}^M_{i=1}$. This results in the following many-body trial wave function:
\begin{equation}\label{eq:MBAnsatz}
    \ket{\Psi(t)}=\sum_{\mathbf{n}}C_\mathbf{n}(t)\ket{\mathbf{n}~t},~\ket{\mathbf{n}~t}=\left(\prod^{M}_{i=1} \frac{\left[\hat{c}^{\dagger}_i(t)\right]^{n_i}}{\sqrt{n_i!}}\right)\ket{\mathbf{0}},~\sum^M_{i=1}n_i=N
\end{equation}
with $\ket{\mathbf{n}~t}$ a time-dependent Fock state, and $\mathbf{n}$ a $M$-dimensional vector with components $n_i$ that represent the occupation numbers of each orbital $\phi_i(\mathbf{x},t)$, whose associated creation operator is $\hat{c}^{\dagger}_i(t)$. For fermions, $n_i=0,1$, which implies $M\geq N$. We note that, in the limit $M\rightarrow\infty$, the MCTDH approximation becomes exact as it spans the complete single-particle Hilbert space.

By inserting this ansatz in the Dirac-Frenkel variational principle (\ref{eq:DiracFrenkel}), one arrives at the MCTDH equations, which take the same form for both bosons and fermions [37,38]
\begin{eqnarray}\label{eq:MCTDH}
i\hbar\frac{dC_\mathbf{n}}{dt}&=&\sum_{\mathbf{n}'}\mathcal{H}_{\mathbf{n}\mathbf{n}'}(t)C_{\mathbf{n}'},~\mathcal{H}_{\mathbf{n}\mathbf{n}'}(t)=\bra{\mathbf{n}~t}\hat{H}-\hat{f}\ket{\mathbf{n}'~t}\\
\nonumber i\hbar\partial_t \phi_{k}&=&f\phi_{k}+Q(t)\left[(h-f)\phi_{k}+\sum^{M}_{i,j,l,m=1}[\rho^{-1}]_{kj}(t)\rho_{ji,lm}(t)W_{lm}(t)\phi_i\right],~Q=1-P, ~P=\sum^M_{k=1}\ket{\phi_k(t)}\bra{\phi_k(t)}
\end{eqnarray}
In this equation, $P$ is the projector onto the restricted single-particle subspace, $\rho_{lk}(t)=\bra{\Psi(t)}\hat{c}_l^{\dagger}(t)\hat{c}_k(t)\ket{\Psi(t)}$ and $\rho_{lk,jm}(t)=\bra{\Psi(t)}\hat{c}_l^{\dagger}(t)\hat{c}_j^{\dagger}(t)\hat{c}_m(t)\hat{c}_k(t)\ket{\Psi(t)}$ are the projected one-body and two-body density matrices, respectively, and $W_{lm}(t)$ is the non-local time-dependent potential
\begin{equation}
    \left[W_{lm}(\mathbf{x},t)\right]_{\alpha\beta}=\sum_{\lambda\mu}\int\mathrm{d}\mathbf{x}'~ V_{\alpha\beta,\lambda\mu}(\mathbf{x}-\mathbf{x}')\left[\phi^*_l(\mathbf{x}',t)\right]_\lambda \left[\phi_m(\mathbf{x}',t)\right]_\mu
\end{equation}
The operators $h,f$ are the single-particle versions of the second-quantized operators $\hat{h},\hat{f}$ acting on the wave functions $\phi_i$. The operator $f$ is characterized by the matrix elements $f_{lk}(t)=\braket{\phi_{l}(t)|i\hbar\partial_t \phi_{k}(t)}$ and can be actually chosen to be any Hermitian operator in the restricted single-particle subspace since the ansatz of Eq. (\ref{eq:MBAnsatz}) is invariant under unitary transformations within the single-particle basis. In particular, we can choose $f$ as the projection of the single-particle Hamiltonian, $f=PhP$, and hence Eqs. (\ref{eq:MCTDH}) are reduced to
\begin{eqnarray}\label{eq:MCTDHSP}
i\hbar\frac{dC_\mathbf{n}}{dt}&=&\sum_{\mathbf{n}'}\mathcal{H}_{\mathbf{n}\mathbf{n}'}(t)C_{\mathbf{n}'},~\mathcal{H}_{\mathbf{n}\mathbf{n}'}(t)=\bra{\mathbf{n}~t}\hat{V}\ket{\mathbf{n}'~t}\\
\nonumber i\hbar\partial_t \phi_{k}&=&h\phi_{k}+Q(t)\sum^{M}_{i,j,l,m=1}[\rho^{-1}]_{kj}(t)\rho_{ji,lm}(t)W_{lm}(t)\phi_i
\end{eqnarray}
We remark that the matrix $\mathcal{H}(t)$ only depends on the time-dependent single-particle orbitals, while the one-body and two-body density matrices only depend on the coefficients $C_\mathbf{n}(t)$.

We now assume that all single-particle wave functions behave as a Floquet wave,
$\phi_{k}(\mathbf{x},t)=e^{-i\varepsilon_kt/\hbar}u_{k}(\mathbf{x},t),~u_{k}(\mathbf{x},t+T)=u_{k}(\mathbf{x},t)$. This implies $\ket{\mathbf{n}~t}=e^{-i\varepsilon_\mathbf{n}t/\hbar}\ket{\bar{\mathbf{n}}~t}$, with $\varepsilon_\mathbf{n}=\sum^{M}_{i=1} n_i\varepsilon_i$ and $\ket{\bar{\mathbf{n}}~t}$ a periodic Fock state. In turn, the time-dependent phase can be absorbed as $C_\mathbf{n}=e^{i\varepsilon_\mathbf{n}t/\hbar}\tilde{C}_\mathbf{n}$, which yields
\begin{equation}
    i\hbar\frac{d\tilde{C}_\mathbf{n}}{dt}=\sum_{\mathbf{n}'}\bar{\mathcal{H}}_{\mathbf{n}\mathbf{n}'}(t)\tilde{C}_{\mathbf{n}'},~\bar{\mathcal{H}}_{\mathbf{n}\mathbf{n}'}(t)=\bra{\bar{\mathbf{n}}~t}\hat{V}\ket{\bar{\mathbf{n}}'~t}+\varepsilon_\mathbf{n}\delta_{\mathbf{n}\mathbf{n}'}
\end{equation}
Since $\bar{\mathcal{H}}(t)$ is a linear periodic operator with respect to $\tilde{C}_\mathbf{n}$, there are solutions of the form $\tilde{C}_\mathbf{n}=e^{-i\varepsilon t/\hbar}\bar{C}_\mathbf{n}(t)$, with $\bar{C}_\mathbf{n}(t)$ periodic. In this case, the total many-body wave function is indeed an SMBF state,
\begin{equation}
        \ket{\Psi(t)}=e^{-i\varepsilon t/\hbar}\ket{u(t)},~\ket{u(t)}=\sum_{\mathbf{n}}\bar{C}_\mathbf{n}(t)\ket{\bar{\mathbf{n}}~t}
\end{equation}
$\ket{u(t)}=\ket{u(t+T)}$ being periodic. This also implies that the one-body and two-body density matrices behave as
\begin{equation}
    \rho_{lk}(t)=e^{-i(\varepsilon_l-\varepsilon_k)t/\hbar}\bar{\rho}_{lk}(t),~
    \rho_{lk,jm}(t)=e^{-i(\varepsilon_l+\varepsilon_j-\varepsilon_m-\varepsilon_k)t/\hbar}\bar{\rho}_{lk,jm}(t)
\end{equation}
In addition, $Q(t)$ is periodic, and $W_{lm}(t)=e^{-i(\varepsilon_m-\varepsilon_l)t/\hbar}\bar{W}_{lm}(t)$, where the functions with a bar are periodic. Therefore, in the equation for the $u_k$,
\begin{equation}
    i\hbar\partial_t u_{k}=(h-\varepsilon_k)u_{k}+Q(t)\sum^{M}_{i,j,l,m=1}[\bar{\rho}^{-1}]_{kj}(t)\bar{\rho}_{ji,lm}(t)\bar{W}_{lm}(t)u_i
\end{equation}
the rightmost term is self-consistently periodic. After expanding in Floquet components each $u_{k}(\mathbf{x},t)$ as in the main text, a complicated system of nonlinear equations can be derived. For a better understanding of the emergence of SMBF states and the behavior of their quantum fluctuations, we consider some well-known limits of the MCTDH equations for both bosons and fermions.

\subsection{Gross-Pitaevskii and Bogoliubov-de Gennes equations}

The GP equation (\ref{eq:HamiltonianEffective}) is obtained by considering the MCTDH equations for $N$ spin-0 bosons condensed in a single orbital $\phi_0(\mathbf{x},t)$. First, we note that the Hamiltonian (\ref{eq:HamiltonianManyBody}) is in fact equivalent to that in Eq. (\ref{eq:HamiltonianManyBodyGeneral}) as the field operator of spin-0 bosons reads
\begin{equation}
    \hat{\Psi}(\mathbf{x})=\sum_k \phi_k(\mathbf{x})\hat{c}_k
\end{equation}
where $k$ here runs over the complete single-particle basis of the Hilbert space.

For $M=1$, all density matrices are just scalars, $\rho_{00}=N$ and $\rho_{00,00}=N(N-1)$, and there is only one coefficient $C_\mathbf{n}$ in the expansion of Eq. (\ref{eq:MBAnsatz}), whose time evolution is just a trivial time-dependent phase. We define a single-particle mean-field operator
\begin{equation}
    h^{\rm{MF}}\equiv  h+(N-1)W_{00},~W_{00}(\mathbf{x},t)=\int\mathrm{d}\mathbf{x}'~ V(\mathbf{x}-\mathbf{x}')|\phi_0(\mathbf{x}',t )|^2
\end{equation}
By taking $f=Ph^{\rm{MF}}P$ in Eq. (\ref{eq:MCTDH}), we get $ i\partial_t\phi_0=h^{\rm{MF}}\phi_0$. The usual form of the GP equation $h^{\rm{MF}}=H_{\rm{GP}}$ is retrieved after defining $\Psi(\mathbf{x},t)\equiv \sqrt{N}\phi_0(\mathbf{x},t)$, where $N\simeq N-1$ for large $N$, and by replacing the actual interaction by a contact pseudopotential $g\delta(\mathbf{x}-\mathbf{x}')$.

A simple derivation of the BdG equations for the quantum fluctuations $\hat{\phi}(\mathbf{x},t)$ is obtained by direct substitution $\Psi(\mathbf{x},t)\rightarrow\Psi(\mathbf{x},t)+\hat{\phi}(\mathbf{x},t)$ in the GP equation and expansion up to linear order. In the specific case of an SMBF state discussed in the main text, one uses a decomposition of the type $\Psi(\mathbf{x},t)\rightarrow\left[u(\mathbf{x},t)+\hat{\phi}(\mathbf{x},t)\right]e^{-i\tilde{\mu}t/\hbar}$,
which yields Eq. (\ref{eq:BdGfieldequation}).

\subsection{Hartree-Fock equations and time-dependent Hartree-Fock approximation}

The Hartree-Fock (HF) equations for fermions are recovered from the MCTDH equations by considering $N$ fermions in $M=N$ orbitals, which means that there is only one Slater determinant present in the expansion of Eq. (\ref{eq:MBAnsatz}). Thus, the one-body and two-body density matrices of Eq. (\ref{eq:MCTDH}) are $\rho_{lk}(t)=\delta_{lk}$ and $\rho_{lk,jm}(t)=\delta_{lk}\delta_{jm}-\delta_{lm}\delta_{jk}$. We define a single-particle mean-field operator, the HF Hamiltonian
\begin{equation}
    h^{\rm{HF}}\equiv h+V^{\rm{HF}},~V^{\rm{HF}}\phi_k=\sum^M_{i=1}W_{ii}\phi_k-W_{ik}\phi_i
\end{equation}
After taking $f=Ph^{\rm{HF}}P$, we recover the usual form of the self-consistent time-dependent Hartree-Fock (TDHF) equations
\begin{equation}\label{eq:TDHFSC}
    i\hbar\partial_t\left[\phi_{k}(\mathbf{x},t)\right]_\alpha=\sum_\beta \int\mathrm{d}\mathbf{x}'~h_{\alpha\beta}^{\rm{HF}}(\mathbf{x},\mathbf{x}',t)\left[\phi_{k}(\mathbf{x}',t)\right]_\beta=\sum_\beta h_{\alpha\beta}(\mathbf{x})\left[\phi_{a}(\mathbf{x},t)\right]_\beta+\int\mathrm{d}\mathbf{x}'~V_{\alpha\beta}^{\rm{HF}}(\mathbf{x},\mathbf{x}',t)\left[\phi_{k}(\mathbf{x}',t)\right]_\beta
\end{equation}
The explicit expression of the nonlocal time-dependent HF potential $V^{\rm{HF}}_{\alpha\beta}(\mathbf{x},\mathbf{x}',t)$
is given by
\begin{equation}\label{eq:TDHFSCPotential}
    V^{\rm{HF}}_{\alpha\beta}(\mathbf{x},\mathbf{x}',t)= \sum^M_{k=1}\sum_{\lambda\mu}\delta(\mathbf{x}-\mathbf{x}')\int\mathrm{d}\mathbf{x}''~V_{\alpha\beta,\lambda\mu}(\mathbf{x}-\mathbf{x}'')\left[\phi^*_k(\mathbf{x}'',t)\right]_\lambda \left[\phi_k(\mathbf{x}'',t)\right]_\mu- V_{\alpha\mu,\lambda\beta}(\mathbf{x}-\mathbf{x}')\left[\phi_k(\mathbf{x},t)\right]_\mu \left[\phi^*_k(\mathbf{x}',t)\right]_\lambda
\end{equation}
If we assume Floquet solutions $\phi_k(\mathbf{x},t)=e^{-i
\varepsilon_k t/\hbar}u_k(\mathbf{x},t)$, with $u_k(\mathbf{x},t+T)=u_k(\mathbf{x},t)$, the HF Hamiltonian then becomes periodic and the Floquet states are self-consistent solutions of the TDHF equations. This is an explicit example of a fermionic SMBF state.

Moreover, the spectrum of quantum fluctuations, computed from the TDHF approximation (TDHFA), also takes the form of a periodic Floquet equation. First, we move to the full single-particle Floquet basis of the periodic HF Hamiltonian $\phi_k(\mathbf{x},t)=e^{-i
\varepsilon_k t/\hbar}u_k(\mathbf{x},t)$, where in the following indices $k=a,b=1\ldots M$ label occupied states and $k=c,d$ unoccupied ones. Due to the formal analogy between the BdG and the TDHFA equations (see for instance Appendix B of Ref. \citeSupp{deNova2017}), quantum fluctuations can be computed by considering small particle-hole excitations of the form:
\begin{equation}\label{eq:TDHFASMBFLinear}
    \phi_{a}(\mathbf{x},t)\rightarrow e^{-i
\varepsilon_a t/\hbar}\left[u_{a}(\mathbf{x},t)+\sum_c \hat{C}_{ca}(t)u_{c}(\mathbf{x},t)\right],~\hat{C}_{ca}(t)\equiv\hat{c}^\dagger_{a}(t)\hat{c}_{c}(t)
\end{equation}
By expanding the TDHF equations (\ref{eq:TDHFSC}) up to linear order and projecting onto the unoccupied states $u_{c}(\mathbf{x},t)$, we obtain the TDHFA equations
\begin{equation}\label{eq:TDHFAFloquetElectronHole}
      i\hbar\partial_t\hat{C}_{ca}=(\varepsilon_c-\varepsilon_a) \hat{C}_{ca}+\sum_{bd}\left[\bar{V}_{ca,bd}(t)-\bar{V}_{cd,ba}(t)\right]\hat{C}_{db}+\left[\bar{V}_{ca,db}(t)-\bar{V}_{cb,da}(t)\right]\hat{C}^{\dagger}_{db}\equiv \sum_{bd} N_{ac,bd}(t)\hat{C}_{db}+A_{ac,db}(t)\hat{C}^{\dagger}_{db}
\end{equation}
or, in matrix notation,
\begin{equation}\label{eq:TDHFAMatrix}
i\hbar\frac{d\hat{\mathbf{C}}}{dt}=\mathbf{X}(t)\hat{\mathbf{C}},~\hat{\mathbf{C}}=\left[\begin{array}{c}
\hat{C}_{ca}(t)\\
\hat{C}^{\dagger}_{ca}(t)
\end{array}\right],~\mathbf{X}=\left[\begin{array}{c|c}
\mathbf{N} & \mathbf{A} \\
\hline
-\mathbf{A}^* & -\mathbf{N}^*
\end{array}\right]
\end{equation}
where $\bar{V}_{lk,jm}(t)$ are the matrix elements of the interacting potential in terms of the $u_k(\mathbf{x},t)$. Since they are periodic, the matrix $\mathbf{X}(t)$ of the TDHFA equations is periodic and thus, the spectrum is also described in terms of Floquet bands. The BdG/TDHFA correspondence implies that a similar expansion to that of Eq. (\ref{eq:BdGmodesexpansion}) can be made \citeSupp{deNova2017}. A more rigorous derivation of the TDHFA equations involves the Heisenberg equation of motion for the operators $\hat{C}_{ca}(t)$, where pairs of operators are contracted in the spirit of the HF approximation to obtain a linear equation, namely, Eq. (\ref{eq:TDHFAFloquetElectronHole}).

\subsection{Bose-Hubbard model}
We consider a general Bose-Hubbard Hamiltonian of the type
\begin{equation}
    \hat{H}=-\sum_{ij} t_{ij}\hat{c}^\dagger_i\hat{c}_j+\sum_i(\epsilon_i-\mu)\hat{c}^\dagger_i\hat{c}_i+\frac{U}{2}\hat{c}^\dagger_i\hat{c}^\dagger_i\hat{c}_i\hat{c}_i
\end{equation}
with $i,j$ labeling lattice sites, $t_{ij}$ the hopping matrix and $\epsilon_i$ some on-site energy. We first consider the discrete GP equation, obtained by an ansatz of the form
\begin{equation}\label{eq:BHGPAnsatz}
    \ket{\Psi(t)}= \frac{\left[\sum_i\phi_i(t)\hat{c}^{\dagger}_i\right]^{N}}{\sqrt{N!}}\ket{\mathbf{0}}
\end{equation}
After inserted into the Dirac-Frenkel variational principle, it gives
\begin{equation}\label{eq:discreteGP}
i\hbar\partial_t\psi_i=\sum_{j}H_{ij}\psi_j,~H_{ij}=-t_{ij}+\delta_{ij}\left[(\epsilon_i-\mu)+U|\psi_i|^2\right]
\end{equation}
with $\psi_i\equiv\sqrt{N}\phi_i$. SMBF states emerge once more when one considers Floquet states $\psi_i(t)=e^{-i \tilde{\mu} t/\hbar}u_i(t)$, $u_i(t+T)=u_i(t)$, since the resulting nonlinear operator $H_{ij}$ is self-consistently periodic. Similarly, the associated discrete BdG equations for the quantum fluctuations (obtained after expanding $\psi_i\rightarrow[u_i(t)+\hat{b}_i]e^{-i \tilde{\mu} t/\hbar}$) are periodic:
\begin{equation}
    i\hbar\partial_t\hat{b}_i=\sum_{j}N_{ij}\hat{b}_j+Uu^2_i(t)\hat{b}^{\dagger}_i,~N_{ij}=-t_{ij}+\delta_{ij}\left[(\epsilon_i-\mu)+2U|u_i(t)|^2-\tilde{\mu}\right]
\end{equation}
Thus, the BdG spectrum is also described in terms of quasi-energy bands.

More interestingly, SMBF states can also arise within a Gutzwiller ansatz, which goes beyond the MCTDH framework and is able to describe the superfluid-Mott insulator phase transition [39],
\begin{equation}
    \ket{\Psi(t)}=\prod_i \ket{\phi_i(t)},~\ket{\phi_i(t)}=\sum^{\infty}_{n=0} c_{i,n}(t)\ket{n}_i,~\sum^{\infty}_{n=0} |c_{i,n}(t)|^2=1
\end{equation}
$\ket{n}_i$ being a Fock state of the lattice site $i$. When inserted into the Dirac-Frenkel variational principle, this ansatz yields the following equations of motion
\begin{equation}\label{eq:Gutzwiller}
    i\hbar\frac{dc_{i,n}}{dt}=\sum^{\infty}_{m=0}\mathcal{H}^{(i)}_{nm}(t)c_{i,m},~\mathcal{H}^{(i)}_{nm}=\delta_{nm}\left[n(\epsilon_i-\mu)+U\frac{n(n-1)}{2}\right]-\sum_{j}t_{ij}\left[\delta_{n,m+1}\sqrt{n}\Psi_j+\delta_{n+1,m}\sqrt{m}\Psi^*_j\right]
\end{equation}
with
\begin{equation}
    \Psi_i(t)=\braket{\hat{c}_i}=\sum^{\infty}_{n=0} \sqrt{n+1}c^*_{i,n}(t)c_{i,n+1}(t)
\end{equation}
Once more, if $\mathcal{H}^{(i)}(t)$ is periodic, one can look for Floquet solutions of Eq. (\ref{eq:Gutzwiller}), $c_{i,n}(t)=e^{-i\varepsilon_it/\hbar}a_{i,n}(t)$, with $a_{i,n}(t)$ periodic. Self-consistently, $\Psi_i(t)$ becomes periodic and, consequently, $\mathcal{H}^{(i)}(t)$ does. The on-site wave functions behave as Floquet states $\ket{\phi_i(t)}=e^{-i\varepsilon_it/\hbar}\ket{u_i(t)}$, and hence the total wave function $\ket{\Psi(t)}$ becomes an SMBF state.

Regarding quantum fluctuations, they are also computed by expanding to linear order as in the previous cases, $c_{i,n}(t)\rightarrow e^{-i\varepsilon_it/\hbar}\left[a_{i,n}(t)+\hat{b}_{i,n}(t)\right]$, which yields periodic BdG equations of the form \citeSupp{Caleffi2020}
\begin{equation}
    i\hbar\partial_t\hat{b}_{i,n}=\sum^{\infty}_{m=0}\sum_{j}N^{ij}_{nm}(t)\hat{b}_{j,m}+A^{ij}_{nm}(t)\hat{b}^{\dagger}_{j,m}
\end{equation}
or, in matrix notation,
\begin{equation}
    i\hbar\frac{d\hat{\mathbf{b}}}{dt}=\mathbf{M}(t)\hat{\mathbf{b}},~\hat{\mathbf{b}}=\left[\begin{array}{c}
\hat{b}_{i,n}\\
\hat{b}^{\dagger}_{i,n}
\end{array}\right],~\mathbf{M}(t)=\left[\begin{array}{c|c}
\mathbf{N} & \mathbf{A} \\
\hline
-\mathbf{A}^* & -\mathbf{N}^*
\end{array}\right]
\end{equation}
where

\begin{eqnarray}
    \nonumber N^{ij}_{nm}(t)&\equiv&{H}^{(i)}_{nm}(t)\delta_{ij}-t_{ij}\left[\sqrt{n}\sqrt{m}a_{n-1,i}(t)a^*_{m-1,j}(t)+\sqrt{n+1}\sqrt{m+1}a_{n+1,i}(t)a^*_{m+1,j}(t)\right]\\
    A^{ij}_{nm}(t)&\equiv&-t_{ij}\left[\sqrt{n}\sqrt{m+1}a_{n-1,i}(t)a_{m+1,j}(t)+\sqrt{n+1}\sqrt{m}a_{n+1,i}(t)a_{m-1,j}(t)\right]
\end{eqnarray}
are periodic operators, and thus, $\mathbf{M}(t)$ is periodic.

\begin{figure}[t]
\begin{tabular}{@{}ccc@{}}

    \includegraphics[width=0.33\columnwidth]{SM/IHFCSupFlecha} &
    \includegraphics[width=0.33\columnwidth]{SM/BHLSup} &
    \includegraphics[width=0.33\columnwidth]{SM/SolitonSupFlecha} \\ \\ \\
    \includegraphics[width=0.33\columnwidth]{SM/IHFC2D} &
    \includegraphics[width=0.33\columnwidth]{SM/BHL2D} &
    \includegraphics[width=0.33\columnwidth]{SM/SOL2D}
\end{tabular}
\caption{Upper row: Spatial profile of sound (solid blue) and flow (dashed red) velocities for different types of initial condition, where the shaded area represents the region in which the attractive constant potential $V(x)=-V_0$ is placed. Lower row: 2D plot of $|\Psi(x,t)|^2$ for the time evolution of each initial condition above, where in all cases $v=0.65,~V_0=1,~X=2$. Left column: IHFC. Central column: BHL. Right column: SGS.}
\label{fig:Model}
\end{figure}

\begin{figure}[t]\includegraphics[width=0.5\columnwidth]{SM/FloquetPartialCompa}
\caption{Comparison of the Floquet components $u_0(x)$ (solid line), $u_{-1}(x)$ (dashed line) and $u_{-2}(x)$ (dashed-dotted line) of the SMBF state for the three simulations of Fig. \ref{fig:Model}. Color code: IHFC (blue), BHL (red), SGS (black).}
\label{fig:FloquetComparison}
\end{figure}

\section{Robustness of the CES time crystal}

We provide here the technical details of the analysis of the robustness of the CES time crystal.

\subsection{Independence of the initial state and transient}

We examine the independence of the SMBF state with respect to the initial condition and the transient details. In particular, we consider three different types of initial condition at $t=0$ that asymptotically match (upstream and downstream) a subsonic plane wave $\Psi(x,0)\sim e^{ivx}$, and whose evolution is described by the same Hamiltonian (including the same attractive square-well potential). Therefore, their evolution is characterized by the same set of parameters $(v,X,V_0)$. A schematic representation of these initial conditions is provided in Fig. \ref{fig:Model}.

The first type of initial condition is the main model considered in this work, starting from an IHFC in which we introduce an attractive square well potential at $t=0$ (left column in Fig. \ref{fig:Model}). In order to analyze the robustness against the transient details, the square well is not suddenly quenched but instead introduced within a time scale $\tau$.

The second type of initial condition is an unstable stationary black-hole laser (BHL) solution within an attractive square well potential [51] (central column in Fig. \ref{fig:Model}). At $t=0$, some small noise is placed on top of it, triggering the dynamical instabilities that eventually will grow up to the nonlinear saturation regime, where the system will reach either the CES state or the GS [40,41].

Finally, we consider an initial condition in which we directly start from the GS, localized around the attractive square well. However, in the upstream region, we introduce a solitonic defect that travels with velocity $v_{\rm{sol}}$ towards the GS, eventually destabilizing it (right column in Fig. \ref{fig:Model}). We refer to this model as the soliton-ground state (SGS) model. We recall that the IHFC and BHL states were already considered in Ref. [41], while the SGS model is novel of this work.

In all cases, for a fixed set of parameters $(v,V_0,X)$ that lies within the CES region of the phase diagram, the same CES state is eventually reached. A comparison of their time evolution is shown in the bottom row of Fig. \ref{fig:Model}, where we observe that, after some transient that strongly depends on the initial condition, they all enter the CES state. Within numerical accuracy, the three CES states oscillate with the same period. Moreover, we have extracted the Floquet components $u_n(x)$ of each CES state and compared them in Fig. \ref{fig:FloquetComparison}, finding an excellent agreement.

When reached, the CES state is also quite insensitive to the details of the transient. Specifically, it is independent 1) for an IHFC, from the time scale $\tau$ at which the potential is introduced (except for sufficiently adiabatic rates, in whose case the system always evolves towards GS); 2) for a BHL solution, from the stochastic initial noise; 3) for the SGS model, from the velocity of the launched soliton $v_{\rm{sol}}$ (except for sufficiently fast solitons, whose short passage through the localized GS is not enough to push it to the CES state).

\begin{figure}[t]
\begin{tabular}{@{}ccc@{}}
    \includegraphics[width=0.33\columnwidth]{SM/Gaussian} &
    \includegraphics[width=0.33\columnwidth]{SM/DeltaArrow} &
    \includegraphics[width=0.33\columnwidth]{SM/FP} \\ \\ \\
    \includegraphics[width=0.33\columnwidth]{SM/Gaussian2D} &
    \includegraphics[width=0.33\columnwidth]{SM/Delta2D} &
    \includegraphics[width=0.33\columnwidth]{SM/FP2D}
\end{tabular}
\caption{Upper row: Spatial profile of sound (solid blue) and flow (dashed red) velocities for different configurations at $t=0$. Left: Attractive Gaussian well of Eq. (\ref{eq:GaussianPotential}), where the shaded area represents the potential profile. The velocity is $v=0.8$ and the potential parameters are $V_0=1, \sigma=1$. Center: Repulsive delta barrier of Eq. (\ref{eq:DeltaPotential}). The parameters are $Z=0.5$ and $v=0.9$. Right: Flat-profile BHL of Eq. (\ref{eq:FlatProfile}). Background parameters are $v=0.85,~c_2=0.4,~X=2$. Lower row: 2D plot of $|\Psi(x,t)|^2$ for the time evolution of the initial conditions of the upper row.}
\label{fig:Universality}
\end{figure}

\begin{figure}[t]
\begin{tabular}{@{}cc@{}}
    \includegraphics[width=0.5\columnwidth]{SM/Disorder1} &
    \includegraphics[width=0.5\columnwidth]{SM/Disorder7} \\ \\ \\
    \includegraphics[width=0.5\columnwidth]{SM/Disorier1} &
    \includegraphics[width=0.5\columnwidth]{SM/Blurier}

\end{tabular}
\caption{Time evolution of an IHFC in which the quenched attractive well is stochastically modulated in time through Eq. (\ref{eq:Disordered}). The background parameters of the flow are  $v=0.95,~V_0=1,~X=2$. Upper row: 2D plot of $|\Psi(x,t)|^2$. Lower row: Fourier spectrum of upper row,  $|\Psi(x,\omega)|^2$. Left column: $\epsilon=0.1$. Right column: $\epsilon=0.7$.}
\label{fig:Disorder}
\end{figure}

\begin{figure}[t]
\begin{tabular}{@{}cc@{}}
    \includegraphics[width=0.5\columnwidth]{SM/Lambda62D} &
    \includegraphics[width=0.5\columnwidth]{SM/Lambda42D} \\ \\ \\
    \includegraphics[width=0.5\columnwidth]{SM/Lambda6CV} &
    \includegraphics[width=0.5\columnwidth]{SM/Lambda4CV}
\end{tabular}
\caption{Analysis of the role of interactions as a function of the parameter $\lambda$ for an IHFC with background parameters $v=0.95,~V_0=1,~X=2$. Upper row: 2D plot of $|\Psi(x,t)|^2$. Lower row: 1D snapshot of the sound (solid blue) and flow (dashed red) velocities for the latest time of upper row.
Left: $\lambda=0.6$. Right: $\lambda=0.4$.}
\label{fig:Interactions}
\end{figure}

\begin{figure}[t]
\begin{tabular}{@{}cc@{}}
    \includegraphics[width=0.5\columnwidth]{SM/Lambda22D} &
    \includegraphics[width=0.5\columnwidth]{SM/Lambda02D} \\ \\ \\
    \includegraphics[width=0.5\columnwidth]{SM/Mu22D} &
    \includegraphics[width=0.5\columnwidth]{SM/Mu02D}
\end{tabular}
\caption{Analysis of the role of interactions as a function of the parameter $\lambda$ for an IHFC with background parameters $v=0.95,~V_0=1,~X=2$. Upper row: 2D plot of $|\Psi(x,t)|^2$. Lower row: 2D plot of $|\Psi(x,\omega)|^2$ for the simulations of upper row, once in the stationary regime. Left: $\lambda=0.2$. Right: $\lambda=0$.}
\label{fig:Scrodinger}
\end{figure}

\begin{figure}[t]
\begin{tabular}{@{}cc@{}}
    \includegraphics[width=0.5\columnwidth]{SM/Lambda222D} &
    \includegraphics[width=0.5\columnwidth]{SM/Lambda202D} \\ \\ \\
    \includegraphics[width=0.5\columnwidth]{SM/Mu222D} &
    \includegraphics[width=0.5\columnwidth]{SM/Mu202D}
\end{tabular}
\caption{Same as Fig. \ref{fig:Scrodinger} but for $V_0=2$.}
\label{fig:Scrodinger2}
\end{figure}

\subsection{Universality}

We examine the universality of the proposed SMBF state by considering other obstacles instead of the idealized attractive square well, finding in all cases that an SMBF state is reached provided that the background flow velocity is sufficiently high. Here we describe in detail three of those models, represented in upper row of Fig. \ref{fig:Universality}.

First, we consider that the potential quenched at $t=0$ is a more realistic Gaussian attractive well
\begin{equation}\label{eq:GaussianPotential}
    V(x)=-V_0e^{-\frac{x^2}{\sigma^2}}
\end{equation}

We also consider a repulsive delta barrier, similarly to Ref. [44],
\begin{equation}\label{eq:DeltaPotential}
    V(x)=Z\delta(x)
\end{equation}

Finally, we consider the flat-profile BHL model of Ref. [40], where the quench involves both the external potential and the coupling constant in such a way that
\begin{equation}\label{eq:FlatProfile}
    g(x)n_0+V(x)=E_b
\end{equation}
with $E_b$ some constant. In this way, an IHFC of density $n_0$ remains as a stationary solution after the quench. In particular, $g(x)$ is chosen as a piecewise function so that the speed of sound satisfies $c(x)=1$ for $|x|>X/2$, and $c(x)=c_2<v$ for $|x|<X/2$. As a result, the central finite region of size $X$ becomes supersonic and, consequently, the IHFC after the quench is an unstable stationary BHL solution. Once more, some small noise is added to trigger the dynamical instabilities.

For all the described models, we have found that an SMBF state can be eventually reached. Examples are displayed in lower row of Fig. \ref{fig:Universality}.

\subsection{Temporal disorder}

We analyze the stability of the time crystal against time-dependent disorder. Specifically, we consider stochastic perturbations resulting from driving the square well potential with a time-dependent modulation
\begin{equation}\label{eq:Disordered}
    V(x,t)=V(x)\left[1+\epsilon h(t)\right]
\end{equation}
with $V(x)$ the square well potential, $h(t)$ normally distributed white noise, and $\epsilon$ the strength of the modulation. We have found that the CES state is quite robust against stochastic perturbations, as can be seen in Fig. \ref{fig:Disorder}. In the left column, we analyze the case of an appreciable modulation with $\epsilon=0.1$, once the system is well in the CES state. We observe that the CES state keeps its periodicity (upper left panel) as well as its Floquet character, as seen from its Fourier spectrum (lower left panel).

Only in the presence of strong disorder is the periodicity spoiled, as seen in right column of Fig. \ref{fig:Disorder}, where we examine an extreme case $\epsilon=0.7$. The distortion of the periodicity (upper right panel) is translated into a blurring of the Floquet components (lower right panel), especially in the downstream region (blue circle).

We have also considered a stochastic driving in the coupling constant $g(t)=g\left[1+\epsilon h(t)\right]$, finding qualitatively similar results.

\subsection{Role of interactions}

Interactions are an essential ingredient for the formation of an SMBF state. In our model, they are hidden behind the system of units, $gn_0=1$. In order to keep track the role of interactions, we start from a regular subsonic IHFC with $v<1$, and introduce a dimensionless parameter $0\leq\lambda\leq1$ that reduces the strength of the interactions as $\lambda|\Psi(x,t)|^2$ (note that an IHFC is a stationary solution of the homogeneous GP equation for $t<0$ regardless the strength of the interaction). The sound speed then becomes $c(x,t)=\sqrt{\lambda|\Psi(x,t)|^2}$ and, in particular, the initial value of the homogeneous speed of sound is $c_0=\sqrt{\lambda}$. Therefore, whenever $\sqrt{\lambda}<v$, the IHFC is supersonic, which is known to be an energetically unstable flow.

We study the formation of an SMBF state starting from a supersonic IHFC in Fig. \ref{fig:Interactions}. For moderate values of the initial supersonic Mach number $M\equiv v/\sqrt{\lambda}$, the same type of SMBF state (upper left panel) is still reached. In the process, the system increases its upstream average density until it becomes subsonic (lower left panel). Since particle number must be conserved, this increase is reached at the expense of the downstream region, which reduces its average density, keeping its initially supersonic character. However, for increasing values of $M$ (right column), more upstream density increase is required to become subsonic. This is translated into an enlarged transient, as seen by the long-lasting presence of the upstream emitted train of solitons (black fringes in the upstream region of upper right panel).

The limit of small $\lambda$ is analyzed in detail Fig. \ref{fig:Scrodinger}, in which the system cannot become subsonic and the upstream black fringe pattern remains stationary, at least within the time scope of the simulations (upper left panel). Nevertheless, even in that regime, the system still behaves as an SMBF state, as seen by its Fourier spectrum (lower left panel). Finally, the case of $\lambda=0$ is well understood since it corresponds to the usual Schr\"odinger equation (upper right panel). Here, what we have for long times is simply the scattering of the incident wave packet $e^{ivx}$ from the left, giving rise to the interference between the incident and the reflected wave (upstream black fringe pattern) and to the transmission of the wave packet (homogeneous downstream density), plus the oscillation of the localized bound state of the well. This is clearly observed in the Fourier spectrum (lower right panel), which only displays two lines at the frequency of the kinetic energy of the scattering state associated to the incident wave, $\omega=v^2/2\approx 0.4512$, and at the negative energy $E_0$ of the only bound state of the well, $\omega=E_0\approx -0.6039$. We recall that the binding energies $E_n<0$ of an attractive square well are computed from the zeros of the transcendental equations
\begin{equation}
    \sqrt{-\frac{E_n+V_0}{E_n}}=\frac{1}{\tan\sqrt{\frac{V_0+E_n}{2}}X},~\sqrt{-\frac{E_n+V_0}{E_n}}=-\tan\sqrt{\frac{V_0+E_n}{2}}X
\end{equation}
where the number of bound states is
\begin{equation}
    N=\left\lfloor{\frac{\sqrt{2V_0}X}{\pi}}\right\rfloor+1
\end{equation}
with $\left\lfloor x \right\rfloor$ the floor function. The trivial linear combination of two eigenstates with different frequencies indeed results in a periodic oscillation of the density, but cannot be regarded as a Floquet state.

An ever more spectacular result is presented in Fig. \ref{fig:Scrodinger2}, where we analyze the case of a square well with larger amplitude $V_0$, resulting in the presence of two bound states with energies $E_0\approx-0.2035$ and $E_1\approx -1.4697$ in the Schr\"odinger case $\lambda=0$. As a result, the non-interacting case (upper right panel) has not even well-defined frequency for the density oscillations since the differences between the frequencies of the spectrum are incommensurate (lower right panel). However, when switching on interactions, we find that, after some larger transient, the system reaches again an SMBF state (upper left panel), where interactions have stabilized the frequency of oscillation and the Floquet character of the Fourier spectrum is recovered (lower left panel). This demonstrates the key role played by interactions in the formation of the SMBF state.

\section{Long-range order in time}

We study the presence of long-range order in time for the one-body correlation function of a condensate. The calculation proceeds along the same lines as the usual study of long-range order in space [52,53], explicitly accounting now for the time dependence. We first consider the case of a 1D condensate. In order to study this problem, a phase-density decomposition of the field operator $\hat{\Psi}(x,t)$ is typically used due to the infrared divergence of the phase fluctuations,
\begin{equation}\label{eq:phasedensity}
 \hat{\Psi}(x,t)=\Psi(x,t)+\hat{\phi}(x,t)=e^{i[\theta(x,t)+\delta\hat{\theta}(x,t)]}\sqrt{n(x,t)+\delta\hat{n}(x,t)}
\end{equation}
with $\hat{\phi}(x,t)$ the quantum fluctuations around the mean-field GP wave function $\Psi(x,t)=\sqrt{n(x,t)}e^{i\theta(x,t)}$. By absorbing the phase of the condensate into the field operator as $\hat{\phi}(x,t)\equiv e^{i\theta(x,t)}\hat{\chi}(x,t)$, we obtain after expanding to linear order in Eq. (\ref{eq:phasedensity}) that
\begin{eqnarray}
\delta\hat{\theta}(x,t)&=&\frac{\hat{\chi}(x,t)-\hat{\chi}^{\dagger}(x,t)}{2i\sqrt{n(x,t)}}\\
\nonumber \delta\hat{n}(x,t)&=&\sqrt{n(x,t)}\left[\hat{\chi}(x,t)+\hat{\chi}^{\dagger}(x,t)\right]
\end{eqnarray}
from where we can express the density and phase fluctuations in terms of the field fluctuations, whose dynamics is in turn described by the BdG equations. It is immediate to see that density and phase fluctuations obey canonical commutation relations
\begin{equation}
    [\delta\hat{n}(x,t),\delta\hat{\theta}(x',t)]=i\delta(x-x')
\end{equation}
since $[\hat{\chi}(x,t),\hat{\chi}^{\dagger}(x',t)]=\delta(x-x')$.

In order to study the long-range order of the one-body correlation function $G(x,x',t,t')=\braket{\hat{\Psi}^{\dagger}(x,t)\hat{\Psi}(x',t')}$, we focus on the simple case of a homogeneous condensate, $\Psi(x,t)=\sqrt{n_0}e^{i\theta_0}e^{-i\mu t}$ where, through the rest of this section, we set $\hbar=m=c_0=k_B=1$ unless otherwise specified, with $c_0=\sqrt{gn_0/m}$ the condensate speed of sound. The fluctuations of the field operator, $\hat{\phi}(x,t)=\hat{\chi}(x,t)e^{i\theta_0}e^{-i\mu t}$, are described by the usual Bogoliubov expansion:
\begin{equation}\label{eq:fieldfluctuations}
\hat{\chi}(x,t)=\frac{1}{\sqrt{L}}\sum_k \hat{\alpha}_{k} u_ke^{i(kx-\Omega_kt)}+\hat{\alpha}^{\dagger}_k v^*_ke^{-i(kx-\Omega_kt)}
\end{equation}
with $\hat{\alpha}_{k}$ the bosonic operator describing a quasiparticle, $L$ the size of the condensate, $u_k,v_k$ the usual Bogoliubov components
\begin{equation}
u_k=\frac{\frac{k^2}{2}+\Omega_k}{\sqrt{2k^2\Omega_k}},~
v_k=\frac{\frac{k^2}{2}-\Omega_k}{\sqrt{2k^2\Omega_k}}
\end{equation}
and $\Omega_k$ the Bogoliubov dispersion relation
\begin{equation}
\Omega_k=\sqrt{k^2+\frac{k^4}{4}}
\end{equation}
where the healing length satisfies $\xi=\hbar/mc_0=1$ in these units. The density and phase fluctuations admit similar expansions
\begin{eqnarray}\label{eq:phasedensityfluctuations}
\nonumber \delta\hat{n}(x,t)&=&\sqrt{\frac{n_0}{L}}\sum_k\hat{\alpha}_{k} \rho_ke^{i(kx-\Omega_kt)}+\hat{\alpha}^{\dagger}_k \rho^*_ke^{-i(kx-\Omega_kt)},~\rho_k=u_k+v_k=\sqrt{\frac{k^2}{2\Omega_k}}\\
\delta\hat{\theta}(x,t)&=&\frac{1}{\sqrt{n_0L}}\sum_k\hat{\alpha}_{k} \theta_ke^{i(kx-\Omega_kt)}+\hat{\alpha}^{\dagger}_k \theta^*_ke^{-i(kx-\Omega_kt)},~\theta_k=\frac{u_k-v_k}{2i}=-i\sqrt{\frac{\Omega_k}{2k^2}}
\end{eqnarray}
Since $\Omega_k\sim |k|$ for low $k$, then $\rho_k\sim |k|^{\frac{1}{2}}$ and $\theta_k\sim |k|^{-\frac{1}{2}}$, which explicitly shows that phase fluctuations are dominant in the infrared limit. However, in the ultraviolet limit of large $k$, both fluctuations are of the same order, since $\Omega_k\simeq k^2/2$ and then $|\rho_k|\simeq 1$,~$|\theta_k|\simeq 1/2$.

Thus, when computing the one-body correlation function, we concentrate on the critical infrared contribution of phase fluctuations and neglect density fluctuations
\begin{equation}\label{eq:OneBodyDensityPhaseFluctuations}
    G(x,x',t,t')\simeq \Psi^*(x,t)\Psi(x',t')\braket{e^{-i\delta\hat{\theta}(x,t)}e^{i\delta\hat{\theta}(x',t')}}
\end{equation}
We recall that the commutator between phase fluctuations at different times is a finite c-number, $\left[\delta\hat{\theta}(x,t),\delta\hat{\theta}(x',t')\right]=C(x-x',t-t')$, with
\begin{equation}\label{eq:PhaseCommutator}
    C(x,t)=-\frac{2i}{n_0L}\sum_k |\theta_k|^2e^{ikx}\sin\Omega_kt=-\frac{i}{\pi n_0}\int^{\infty}_{-\infty}\mathrm{d}k~ |\theta_k|^2e^{ikx}\sin\Omega_kt
\end{equation}
By employing the usual relation for two operators $A,B$ that satisfy $[A,[A,B]]=[B,[A,B]]=0$,
\begin{equation}
    e^{A}e^{B}=e^{A+B}e^{\frac{[A,B]}{2}}
\end{equation}
we get
\begin{equation}\label{eq:OneBodyExpectationValuePhase}
    G(x,x',t,t')=\Psi^*(x,t)\Psi(x',t')e^{\frac{C(x-x',t-t')}{2}}\braket{e^{i\left[\delta\hat{\theta}(x',t')-\delta\hat{\theta}(x,t)\right]}}
\end{equation}
The expectation value is computed from the general case of the expectation value in a Gaussian state of the imaginary exponential of a Hermitian operator $\hat{A}$ that is linear in annihilation and creation operators, and with vanishing expectation value $\braket{\hat{A}}=0$. Expanding the exponential yields
\begin{equation}
    \braket{e^{i\hat{A}}}=\sum^{\infty}_{n=0}\frac{\braket{i^n\hat{A}^n}}{n!}=
    \sum^{\infty}_{n=0}(-1)^n\frac{\braket{\hat{A}^{2n}}}{2n!}
\end{equation}
while Wick's theorem gives
\begin{equation}
    \braket{\hat{A}^{2n}}=\frac{2n!}{2^nn!}\braket{\hat{A}^2}
\end{equation}
where the prefactor takes into account the number of different pairs that can be taken from $2n$ elements. divided by $n!$ to avoid double-counting. Therefore,
\begin{equation}\label{eq:exponentialwick}
    \braket{e^{i\hat{A}}}=
    \sum^{\infty}_{n=0}(-1)^n\frac{\braket{\hat{A}^{2}}}{2^n n!}=e^{-\frac{\braket{\hat{A}^{2}}}{2}}
\end{equation}

Since phase fluctuations satisfy
\begin{eqnarray}\label{eq:PreGaussian}
\braket{\delta\hat{\theta}(x',t')-\delta\hat{\theta}(x,t)}&=&0 \\
\nonumber \braket{\left[\delta\hat{\theta}(x',t')-\delta\hat{\theta}(x,t)\right]^2}&\equiv&\Theta(x,x',t,t')=\braket{\delta\hat{\theta}^2(x',t')+\delta\hat{\theta}^2(x,t)-\delta\hat{\theta}(x',t')\delta\hat{\theta}(x,t)-\delta\hat{\theta}(x,t)\delta\hat{\theta}(x',t')}
\end{eqnarray}
and our state is Gaussian (either a $T=0$ ground state or a $T>0$ thermal state), we can make use of Eq. (\ref{eq:exponentialwick}), finding
\begin{equation}\label{eq:OneBodyCorrelationFunctionPhaseCorrelations}
G(x,x',t,t')=\Psi^*(x,t)\Psi(x',t')e^{\frac{C(x-x',t-t')}{2}}e^{-\frac{\Theta(x,x',t,t')}{2}}
\end{equation}
Thus, we only need to compute the phase correlations. Proceeding in a similar way to Eq. (\ref{eq:PhaseCommutator}) yields
\begin{eqnarray}\label{eq:phasecorrelation}
 \braket{\delta\hat{\theta}(x,t)\delta\hat{\theta}(x',t')}=\frac{1}{2\pi n_0}\int^{\infty}_{-\infty}\mathrm{d}k~ |\theta_k|^2\left[e^{ik(x-x')}e^{-i\Omega_k(t-t')}+2n_k\cos\left[k(x-x')-\Omega_k(t-t')\right]\right]
\end{eqnarray}
with $n_k$ the Bose occupation factor
\begin{equation}\label{eq:OccupationNumber}
    n_k=\frac{1}{e^{\frac{\Omega_k}{T}}-1}
\end{equation}
which results from considering a thermal Gaussian state with temperature $T$. By combining Eq. (\ref{eq:phasecorrelation}) with Eq. (\ref{eq:PreGaussian}), it is immediate to see that $\Theta(x,x',t,t')=\Theta(x-x',t-t')$, with $\Theta(x,t)=\Theta_0(x,t)+\Theta_T(x,t)$, where $\Theta_0(x,t)$ is the zero-temperature contribution
\begin{eqnarray}\label{eq:CorrelationZeroTemperature}
\Theta_0(x,t)=\frac{1}{\pi n_0}\int^{\infty}_{-\infty}\mathrm{d}k~ |\theta_k|^2\left[1-\cos\left(kx-\Omega_kt\right)\right]
\end{eqnarray}
and $\Theta_T(x,t)$ is the thermal contribution
\begin{eqnarray}\label{eq:CorrelationFiniteTemperature}
\Theta_T(x,t)=\frac{2}{\pi n_0}\int^{\infty}_{-\infty}\mathrm{d}k~ |\theta_k|^2n_k\left[1-\cos\left(kx-\Omega_kt\right)\right]
\end{eqnarray}
In the following, as we are examining the existence of long-range order, we focus on the asymptotic behavior for large $x$ and $t$. Due to dispersive effects, using a saddle-point approximation for large $t$ gives $C(x,t)\sim |t|^{-1/2}\simeq 0$, so phase fluctuations commute. Moreover, we note that for deriving Eq. (\ref{eq:OneBodyCorrelationFunctionPhaseCorrelations}) we have neglected density fluctuations. This implicitly assumes a long-wavelength approximation, since at high momentum the magnitude of both fluctuations is similar. Thus, we must cut all the integrals at values $|k|\simeq 1$. In the same spirit, we approximate all quantities by their low $k$ expansions in order to study possible divergences:
\begin{equation}
    \Omega_k\simeq |k|,~|\theta_k|^2\simeq\frac{1}{2|k|},~n_k\simeq\frac{1}{e^{\frac{|k|}{T}}-1}\simeq \frac{T}{|k|}
\end{equation}
We consider first the $T=0$ case,
\begin{eqnarray}\label{eq:CorrelationZeroTemperatureCut}
\Theta_0(x,t)&\simeq&\frac{1}{\pi n_0}\int^{1}_{-1}\mathrm{d}k~ |\theta_k|^2\left[1-\cos\left(kx-\Omega_kt\right)\right]
\simeq\frac{1}{2\pi n_0}\int^{1}_{-1}\mathrm{d}k~ \frac{1-\cos\left(kx-\Omega_kt\right)}{|k|}\\
\nonumber &=&\frac{1}{2\pi n_0}\left[\int^{1}_{0}\mathrm{d}k~ \frac{1-\cos k(x-t)}{k}+\int^{1}_{0}\mathrm{d}k~ \frac{1-\cos k(x+t)}{k}\right]
=\frac{1}{2\pi n_0}\left[f(x-t)+f(x+t)\right]
\end{eqnarray}
where we have defined the function
\begin{equation}
    f(x)\equiv\int^{1}_{0}\mathrm{d}k~ \frac{1-\cos kx}{k}=
    \int^{x}_{0}\mathrm{d}z~ \frac{1-\cos z}{z}
\end{equation}
This function has a logarithmic divergence for large $x$, $f(x)\simeq \ln |x|$. As a result, long-range phase-fluctuations diverge logarithmically as
\begin{eqnarray}\label{eq:CorrelationZeroTemperatureDivergence}
\Theta_0(x,t)&=&\frac{1}{2\pi n_0}\left[f(x-t)+f(x+t)\right]
\simeq\frac{1}{2\pi n_0}\left[\ln|x-t|+\ln|x+t|\right]
=\frac{\ln \left|x^2-t^2\right|}{2\pi n_0}
\end{eqnarray}

Gathering all together, and after momentarily restoring units, we find that the one-body correlation function at $T=0$ decays algebraically as
\begin{eqnarray}\label{eq:CorrelationFunctionZeroTemperatureLongRange}
    G(x,x',t,t')&\simeq&\Psi^*(x,t)\Psi(x',t')\left|\frac{(x-x')^2-c^2_0(t-t')^2}{\xi^2}\right|^{-\frac{1}{4\pi n_0\xi}}
\end{eqnarray}

The $T>0$ contribution for the phase correlations gives
\begin{equation}
\Theta_T(x,t)\simeq\frac{2}{\pi n_0}\int^{1}_{-1}\mathrm{d}k~ |\theta_k|^2n_k\left[1-\cos\left(kx-\Omega_kt\right)\right]\simeq
\frac{T}{\pi n_0}\int^{1}_{-1}\mathrm{d}k~ \frac{1-\cos\left(kx-|k|t\right)}{k^2}=\frac{T}{\pi n_0}\left[g(x-t)+g(x+t)\right]
\end{equation}
where
\begin{equation}
    g(x)\equiv\int^{1}_{0}\mathrm{d}k~ \frac{1-\cos kx}{k^2}=x\int^{x}_{0}\mathrm{d}z~ \frac{1-\cos z}{z^2}
\end{equation}
In the limit of large $x$, it behaves as
\begin{equation}
    g(x)\simeq |x|\int^{\infty}_{0}\mathrm{d}z~ \frac{1-\cos z}{z^2}=\frac{\pi}{2}|x|
\end{equation}
Hence, phase fluctuations diverge linearly for $T>0$:
\begin{eqnarray}\label{eq:CorrelationFiniteTemperatureLinear}
\Theta_T(x,t)&=&\frac{T}{\pi n_0}\left[g(x-t)+g(x+t)\right]
\simeq\frac{T}{2 n_0}\left[|x-t|+|x+t|\right]
\end{eqnarray}
After restoring units once more, this results in a exponential decay of the long-range order of the one-body correlation function:
\begin{eqnarray}\label{eq:CorrelationFunctionZeroTemperatureLongRange}
    G(x,x',t,t')&\simeq&\Psi^*(x,t)\Psi(x',t')e^{-\frac{mk_BT}{4 n_0\hbar^2 }\left[|(x-x')-c_0(t-t')|+|(x-x')+c_0(t-t')|\right]}
\end{eqnarray}
The phase-coherence length of the condensate is thus
\begin{equation}
    L_{\theta}=\frac{4n_0\hbar^2}{mk_BT}=\frac{4\hbar^2}{m\xi^2k_BT}(n_0\xi)\xi=
    4\frac{\mu_0}{k_BT}(n_0\xi)\xi
\end{equation}
with $\mu_0=gn_0$ the chemical potential.

In summary, long-range order in time is suppressed in the same manner as long-range order in space in 1D:
\begin{eqnarray}
    G(x,x',t,t')&\sim& \left|\frac{(x-x')^2-c^2_0(t-t')^2}{\xi^2}\right|^{-\frac{1}{4\pi n_0\xi}},~T=0\\
\nonumber G(x,x',t,t')&\sim& e^{-\frac{mk_BT}{4 n_0\hbar^2 }\left[|(x-x')-c_0(t-t')|+|(x-x')+c_0(t-t')|\right]},~T>0
\end{eqnarray}
i.e., algebraically at $T=0$ and exponentially for $T>0$.

We can repeat the same calculations for 2D. We have that phase correlations then read
\begin{eqnarray}\label{eq:phasecorrelation2D}
 \braket{\delta\hat{\theta}(\mathbf{x},t)\delta\hat{\theta}(\mathbf{x}',t')}=\frac{1}{(2\pi)^2 n_0}\int\mathrm{d}^2\mathbf{k}~ |\theta_\mathbf{k}|^2\left[e^{i\mathbf{k}(\mathbf{x}-\mathbf{x}')}e^{-i\Omega_\mathbf{k}(t-t')}+2n_\mathbf{k}\cos\left[\mathbf{k}(\mathbf{x}-\mathbf{x}')-\Omega_\mathbf{k}(t-t')\right]\right]
\end{eqnarray}
After cutting the integrals for values $k=|\mathbf{k}|<1$, for $T=0$ we have
\begin{eqnarray}\label{eq:CorrelationZeroTemperature2D}
\Theta_0(\mathbf{x},t)=\frac{1}{4\pi^2 n_0}\int^1_0\mathrm{d}k\int^{2\pi}_0\mathrm{d}\phi~ \left[1-\cos\left(kr\cos\phi-\Omega_kt\right)\right]
\end{eqnarray}
while for $T>0$,
\begin{eqnarray}\label{eq:CorrelationFiniteTemperature2D}
\Theta_T(\mathbf{x},t)=\frac{T}{2\pi^2 n_0}\int^1_0\mathrm{d}k\int^{2\pi}_0\mathrm{d}\phi~ \frac{1-\cos\left(kr\cos\phi-\Omega_kt\right)}{k}
\end{eqnarray}
$\phi$ being the polar angle between $\mathbf{k}\mathbf{x}=kr\cos\phi$ and $r=|\mathbf{x}|$.

\textit{Mutatis mutandis}, we find that at $T=0$ there is indeed long-range order in time for the one-body correlation function, while for $T>0$ it decays algebraically with the same exponent as in space,
\begin{equation}
    \eta=\frac{k_BT}{\hbar c_02\pi n_0\xi}=\frac{1}{2\pi n_0\xi^2}\frac{k_BT}{\mu_0}
\end{equation}

In summary, long-range order in time for the one-body correlation function behaves in exactly the same way as in space for the corresponding spatial dimension. Naively, one could think that, since time alone is 1D, it should always behave like 1D in space. However, since the divergence is controlled by the infrared behavior of the density of states, set by the number of spatial dimensions, and that in this limit $\Omega_k\sim k$, long-range order in time presents the same dependence as in space.

\section{Truncated Wigner calculation of quantum fluctuations}

We explore the robustness of the CES time crystal against quantum fluctuations via the Truncated Wigner method [54,55]. The Truncated Wigner approximation computes symmetric-ordered expectation values from ensemble averages of integrations of the GP equation
\begin{equation}
    i\hbar\partial_t\Psi_W(x,t)=\left[-\frac{\hbar^2}{2m}\partial_x^2+V(x)+g|\Psi_W(x,t)|^2\right]\Psi_W(x,t)
\end{equation}
where the stochastic initial condition is
\begin{equation}
\Psi_W(x,0)=\sqrt{n_0(x)+\delta n_W(x)}e^{i\theta_0(x)}e^{i\delta \theta_W(x)}
\end{equation}
with $\Psi_0(x)=\sqrt{n_0(x)}e^{i\theta_0(x)}$ the GP wave function describing the initial condensate at equilibrium and
\begin{eqnarray}
\delta\theta_W(x)&=&\frac{\chi_W(x)-\chi^*_W(x)}{2i\sqrt{n_0(x)}}\\
\nonumber \delta n_W(x)&=&\sqrt{n_0(x)}\left[\chi_W(x)+\chi^*_W(x)\right]
\end{eqnarray}
where $\chi_W(x)e^{i\theta_0(x)}=\phi_W(x,0)$ describes the quantum fluctuations of the field operator, sampled from the Wigner distribution describing the condensate. Due to the critical role of phase fluctuations, we are employing a quasi-condensate description to describe the present physics \citeSupp{Martin2010PRL,Martin2010NJP}.

In the particular case of this work, $\Psi_0(x)=e^{ivx}$, where we are recovering the system of units of the main text. Then, we have that
\begin{equation}\label{eq:TWExpansionNotConserving}
\Psi_W(x,0)=\sqrt{1+\frac{\delta n_W(x)}{n_0}}e^{ivx}e^{i\delta \theta_W(x)}
\end{equation}
The field fluctuations are written similarly to Eq. (\ref{eq:fieldfluctuations}),
\begin{equation}
\chi_W(x)=\frac{1}{\sqrt{L}}\sum_k \alpha_{k} u_ke^{ikx}+\alpha^*_k v^*_ke^{-ikx}
\end{equation}
as well as the density and phase fluctuations,
\begin{eqnarray}
\nonumber \delta n_W(x)&=&\sqrt{\frac{n_0}{L}}\sum_k \alpha_{k} \rho_ke^{ikx}+\alpha^*_k \rho^*_ke^{-ikx}\\
\delta \theta_W(x)&=&\frac{1}{\sqrt{n_0L}}\sum_k \alpha_{k} \theta_ke^{ikx}+\alpha^*_k \theta^*_ke^{-ikx}
\end{eqnarray}
The amplitudes $\alpha_k$ of each mode are stochastic variables that are sampled from the Wigner distribution of the initial equilibrium state. Our initial quantum state is the $T=0$ ground state in the comoving frame of the condensate, where it is at rest. Hence, our stochastic initial condition is sampled from a Gaussian distribution characterized by
\begin{equation}
\braket{\alpha_k}=0,~\braket{\alpha_{k'}\alpha_k}=\braket{\alpha^*_{k'}\alpha^*_k}=0,~\braket{\alpha^*_{k'}\alpha_k}=\frac{\delta_{kk'}}{2}
\end{equation}

Since our goal is the study of long-range order in time and space of the one-body correlation function $G(x,x',t,t')$, we restrict to large values of $|x-x'|,|t-t'|$, so
\begin{equation}
\braket{\Psi^*_W(x,t)\Psi_W(x',t')}=\frac{\braket{\hat{\Psi}^{\dagger}(x,t)\hat{\Psi}(x',t')}+\braket{\hat{\Psi}(x',t')\hat{\Psi}^{\dagger}(x,t)}}{2}\simeq\braket{\hat{\Psi}^{\dagger}(x,t)\hat{\Psi}(x',t')}=G(x,x',t,t')
\end{equation}
where the commutator for sufficiently separated points in space-time
\begin{equation}\label{eq:fieldcommutator}
 \left[\hat{\Psi}(x',t'),\hat{\Psi}^{\dagger}(x,t)\right]=\left[\hat{\phi}(x',t'),\hat{\phi}^{\dagger}(x,t)\right]=
\frac{1}{2\pi}\int^{\infty}_{-\infty}\mathrm{d}k~ e^{ik(x-x')}e^{-i\Omega_k(t-t')}-\frac{i}{\pi }\int^{\infty}_{-\infty}\mathrm{d}k~ e^{ik(x-x')}\sin\Omega_k(t-t')|v_k|^2
\end{equation}
is neglected on similar grounds to those invoked in the case of phase fluctuations [see discussion after Eq. (\ref{eq:CorrelationFiniteTemperature})].

\begin{figure}[t]
\begin{tabular}{@{}cc@{}}
    \includegraphics[width=0.5\columnwidth]{SM/Gtt} &
    \includegraphics[width=0.5\columnwidth]{SM/Gxt}
\end{tabular}
\caption{Absolute value of the correlation function $|G(x,x',t,t')|$. Left: Mean-field value of $|G(x,x',t,t')|$ for $x=-30$ and $x'=30$ as a function of $(t,t')$. Right: Comparison of $|G(x,x',t,0)|$ for $-x=x'=100$ as a function of $t$ between a Truncated Wigner simulation that includes quantum fluctuations (solid blue) and the mean-field prediction (dashed-black).}
\label{fig:TWA}
\end{figure}

For the simulations, we have chosen as mean-field parameters those of lower Fig. 2 of the main text, $v=0.95,~V_0=1,~X=2$. With respect to the specific parameters of the Truncated Wigner simulation, we have used $L\simeq 1885$ and $N=10^8$, with $n_0=N/L$ the condensate density. The number of modes is $N_m=3000$, which corresponds to a cut-off in $k$ space of $|k|<5$. Expectation values are evaluated after ensembles of $1000$ simulations.

Left Fig. \ref{fig:TWA}  shows the time dependence of the mean-field value of $|G(x,x',t,t')|$ for fixed spatial points far away, one upstream and one downstream. The one-body correlation function clearly exhibits long-range time-periodicity in both $t,t'$, revealing the time crystal character of the CES state. Quantum fluctuations, included via the Truncated Wigner method, do not destroy the time-periodic long-range order for typical times of experiments, thus allowing for a potential experimental observation. This can be seen in right Fig. \ref{fig:TWA}, where we compare the mean-field result for $|G(x,x',t,0)|$ with a Truncated Wigner simulation for long times.

Technically, one should use a number conserving approximation in the computations, taking the number of particles in the condensate as $N_0=N-N_{\rm{NC}}$, with $N_{\rm{NC}}$ the number of non-condensed particles
\begin{equation}
    N_{\rm{NC}}=\sum_k \left(|u_k|^2+|v_k|^2\right)\braket{\hat{\alpha}_k^{\dagger}\hat{\alpha}_k}+|v_k|^2=\sum_k \left(|u_k|^2+|v_k|^2\right)\left(|\alpha_k|^2-\frac{1}{2}\right)+|v_k|^2
\end{equation}
This is translated into a modified version of Eq. (\ref{eq:TWExpansionNotConserving}),
\begin{equation}\label{eq:TWExpansionConserving}
\Psi_W(x,0)=\sqrt{1-\frac{N_{\rm{NC}}}{N}+\frac{\delta n_W(x)}{n_0}}e^{ivx}e^{i\delta \theta_W(x)}
\end{equation}
For the parameters considered, we have not found any significant effect resulting from the use of a non-conserving approximation.

\bibliographystyleSupp{apsrev4-1}
\bibliographySupp{Hawking}

\end{document}