# Absolutely Stable Discrete Time Crystals

Krzysztof Giergiel  Optical Sciences Centre, Swinburne University of
Technology, Melbourne 3122, Australia  Instytut Fizyki Teoretycznej,
Uniwersytet Jagielloński, ulica Profesora Stanislawa Lojasiewicza 11,
PL-30-348 Kraków, Poland  Jia Wang  Centre for Quantum Technology Theory,
Swinburne University of Technology, Melbourne 3122, Australia  Bryan J. Dalton
Centre for Quantum Technology Theory, Swinburne University of Technology,
Melbourne 3122, Australia  Peter Hannaford  Optical Sciences Centre, Swinburne
University of Technology, Melbourne 3122, Australia  Krzysztof Sacha  Instytut
Fizyki Teoretycznej, Uniwersytet Jagielloński, ulica Profesora Stanislawa
Lojasiewicza 11, PL-30-348 Kraków, Poland

###### Abstract

We show that interacting bosons on a ring which are driven periodically by a
rotating lattice potential can support absolutely stable discrete time
crystals. The absolute stability is demonstrated by an exact mapping of
discrete time crystal states to low-lying eigenstates of a time-independent
model that reveals spontaneous breaking of space translation symmetry. The
mapping ensures that there are no residual time-dependent terms that could
lead to heating of the system and destruction of discrete time crystals. With
the help of the Bethe ansatz solutions we also analyze periodically kicked
bosons where the mapping is approximate only and cannot guarantee the absolute
stability of discrete time crystals. However, the kicked boson model shows a
richer interplay between time and space symmetry breaking.

Ordinary space crystals correspond to periodic distributions of particles in
space which form despite the fact that the Hamiltonian is invariant under a
translation of all the particles by an arbitrary vector. This is the
phenomenon of spontaneous breaking of the continuous space translation
symmetry into a discrete space translation symmetry in a time-independent
many-body system. In 2012, research on time crystals was initiated where
spontaneous breaking of time translation symmetry is responsible for the
formation of a crystalline structure Wilczek (2012). Periodically driven many-
body closed systems can form discrete time crystals as a non-equilibrium
quantum state in which the long-time stable system response occurs at integer
multiples of the drive period Sacha (2015a); Khemani et al. (2016); Else et
al. (2016); Yao et al. (2017). The discrete time translation symmetry of a
time-periodic Hamiltonian is spontaneously broken and new periodic motion
emerges — a novel crystalline structure appears in time Sacha (2015a).

The stability of discrete time crystals is much less obvious than for space
crystals because the quasi-energy spectrum of the Floquet Hamiltonian has no
lower bound Sacha (2020). The key question arises whether discrete time
crystals really break ergodicity and are absolutely stable or they gradually
absorb energy from the drive and eventually heat up to a structureless
infinite temperature state as expected for a generic periodically driven many-
body system D’Alessio and Rigol (2014); Lazarides et al. (2014); Ponte et al.
(2015a).

Systems which do not thermalize are typically integrable systems but
integrability is fragile and usually requires fine-tuning of system parameters
Rigol et al. (2007); Santos and Rigol (2010a); Rigol and Santos (2010); Santos
and Rigol (2010b); Cassidy et al. (2011); Caux and Essler (2013). Time-
independent many-body systems in the presence of disorder can reveal many-body
localization (MBL) and emergent integrability Nandkishore and Huse (2015);
Altman and Vosk (2015); D’Alessio et al. (2016); Mierzejewski et al. (2018).
Therefore, periodically driven MBL systems are potential candidates for
realization of absolutely stable discrete time crystals Else et al. (2016);
Khemani et al. (2016); Else et al. (2017). While in the time-independent case,
the MBL can be proven Imbrie (2016), the absolute stability of discrete time
crystals in driven MBL systems Ponte et al. (2015b); Lazarides et al. (2015);
Abanin et al. (2016) is far from being obvious. Also, the absolute stability
of the first proposed discrete time crystal in ultra-cold bosonic atoms
bouncing resonantly on an oscillating atom mirror Sacha (2015a) has not been
proven. While previous studies using the truncated Wigner approximation,
Bogoliubov approach and two-mode theory have demonstrated the absence of
heating and ensured the stability of discrete time crystals for evolution
times longer than those of realistic experiments Wang et al. (2021a); Kuroś et
al. (2020); Wang et al. (2021b), their validity in the context of an infinite
number of bosons and infinite evolution time remains unproven. It should be
stressed though that even if theoretical analyses of the discrete time
crystals in both ultra-cold atoms bouncing on a mirror and MBL systems can
only guarantee the lack of heating and the stability of the discrete time
crystals for a finite evolution time, this time is much longer than the
duration of the experiments Zhang et al. (2017); Choi et al. (2017); Pal et
al. (2018); Rovny et al. (2018); Smits et al. (2018); Mi et al. (2022);
Randall et al. (2021); Frey and Rachel (2022).

Here, we show that resonantly driven bosons by a rotating lattice on a ring in
the Lieb-Liniger (LL) model can reveal absolutely stable discrete time
crystals. Eigenstates of the Floquet Hamiltonian that break the discrete time
translation symmetry can be mapped to low energy eigenstates of a time-
independent Hamiltonian in the rotating frame. The latter eigenstates are
stable and reveal spontaneous breaking of space translation symmetry in the
rotating frame, which corresponds to spontaneous breaking of the discrete time
translation symmetry of the original model. This model concretely provides an
example of the absence of quantum thermalization and the existence of stable
time crystals in a closed quantum many-body system without the need for either
disorder or integrability. Furthermore, an experimental realisation of this
model should be possible, since ultra-cold atoms on a ring have been created
and applying a periodic drive potential should be straightforward Gupta et al.
(2005).

Let us consider $N$ bosons with contact interactions on a ring with
circumference $2\pi$ which are periodically driven by a rotating lattice
potential $\cos(sx-\omega t)$. In this Letter we focus on $s=2$ but it can be
any integer number and allows realization of big discrete time crystals, see
Giergiel et al. (2018); Hannaford and Sacha (2022). The system reduces to the
LL Hamiltonian $H_{\rm LL}$ with an additional time-periodic drive

$\displaystyle H$ $\displaystyle=$ $\displaystyle H_{\rm
LL}+\lambda\sum_{i=1}^{N}\cos(2x_{i}-\omega t),$ (1) $\displaystyle H_{\rm
LL}$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{N}\frac{p_{i}^{2}}{2}+g_{0}\sum_{i<j}^{N}\delta(x_{i}-x_{j}),$
(2)

where we use $R$ and $\hbar^{2}/mR^{2}$ for the length and energy units,
respectively, where $R$ is the ring radius and $m$ the mass of the bosons,
$\lambda$ is the amplitude of the external driving potential and $g_{0}$
stands for the strength of the contact interactions 111The interaction
strength $g_{0}=2mR\omega_{\perp}a_{s}/\hbar$, where $\omega_{\perp}$ is the
frequency of the harmonic transverse confiniment and $a_{s}$ is the atomic
s-wave scattering length.. The system possesses discrete time and space
translation symmetries, i.e., the Hamiltonian (1) does not change if
$t\rightarrow t+T$ (where $T=2\pi/\omega$) or all $x_{i}\rightarrow
x_{i}+\pi$. We will see that for sufficiently strong interactions, these
symmetries can spontaneously be broken and the system starts evolving with a
period longer than $T$ forming a discrete time crystal.

Let us investigate the system in the moving frame of the rotating lattice,
where we first perform a time-dependent unitary transformation
$U_{t}=\exp(i\sum_{j}p_{j}\omega t/2)$, leading to a shift in the positions,
$x_{i}\rightarrow x_{i}+\omega t/2$, and next a second time-independent
unitary transformation $U_{p}=\exp\left(-i\sum_{j}x_{j}\omega/2\right)$,
leading to a shift in the momenta, $p_{i}\rightarrow p_{i}+\omega/2$ 222Note,
that if $\omega/2$ is not an integer number, the original periodic boundary
conditions become twisted boundary conditions after the transformation, see SM
. For simplicity, we always choose integer values of $\omega/2$.. Under these
transformations, we end up with the following exact time-independent
Hamiltonian

$\tilde{H}=\sum_{i=1}^{N}\left[\frac{p_{i}^{2}}{2}+\lambda\cos(2x_{i})\right]+g_{0}\sum_{i<j}^{N}\delta(x_{i}-x_{j}),$
(3)

where a constant term has been omitted. The state vector in the moving frame
is related to that in the laboratory frame via $\tilde{\psi}=U_{p}U_{t}\psi$
SM .

Suppose that the discrete space translation symmetry of the Hamiltonian (3),
i.e., the invariance under the shift $x_{i}\rightarrow x_{i}+\pi$, is
spontaneously broken in the ground state and low-energy eigenstates and only
the symmetry related to the periodic boundary conditions on a ring remains.
The corresponding time-independent single-particle probability densities in
the moving frame, $\tilde{P}(x)=\int dx_{2}\dots
dx_{N}|\tilde{\psi}(x,x_{2},\dots,.x_{N})|^{2}$, fulfill the periodic boundary
conditions, $\tilde{P}(x+2\pi)=\tilde{P}(x)$, but
$\tilde{P}(x+\pi)\neq\tilde{P}(x)$. When we return to the original laboratory
frame by means of the inverse $U_{p}$ and $U_{t}$ transformations, the lab
probability densities read $P(x,t)=\tilde{P}(x-\omega t/2)$ SM  and thus, due
to the spontaneous breaking of the space translation symmetry of the
Hamiltonian (3), they now also reveal spontaneous breaking of the discrete
time translation symmetry of the Hamiltonian (1), since $P(x,t)$ is periodic
with the period $2T$ but not with the period $T$. Hence, the system
spontaneously starts evolving with a period which is an integer multiple of
the period dictated by the drive and forms a discrete time crystal. Since the
equilibrium symmetry breaking states in the moving frame are assumed to have
an infinite lifetime in the thermodynamic limit (i.e., for
$N\rightarrow\infty$, $g_{0}\rightarrow 0$ but $g_{0}N=\rm const$), the
corresponding lab-frame states also exist for infinitely long time, and hence
can be rigorously recognized as an absolutely stable time crystal.

We now proceed to investigate the spontaneous symmetry breaking of our model
in detail. In the presence of the external potential ($\lambda\neq 0$) and
when the interaction is either attractive ($g_{0}<0$) but small, or is
repulsive ($g_{0}>0$), the ground state of the system (3) is a Bose-Einstein
condensate (BEC) where all bosons occupy the single-particle ground state
which is a balanced superposition of two wave-packets localized in each well
of the external potential in (3) Pethick and Smith (2002). The width of the
wave-packets can be estimated by employing the harmonic oscillator
approximation for the potential wells and it reads $\sigma\approx
1/\sqrt{\Omega}$, where the frequency of the harmonic oscillator
$\Omega=2\sqrt{\lambda}$. In the presence of the attractive interactions
($g_{0}<0$) and increasing their strength, we enter a self-trapping regime for
bosons where the mean-field lowest energy solutions are degenerate and each of
them can be approximated by all bosons occupying a single wave-packet
localized in one well of the external potential Pethick and Smith (2002);
Jackson and Weinstein (2004). The space translation symmetry of the
Hamiltonian (3) is spontaneously broken and the self-trapped states live
forever in the thermodynamic limit. Increasing further the strength of the
attractive interactions, we enter the bright soliton regime Pethick and Smith
(2002), i.e., the wave-packets localized in the potential wells start
shrinking due to strong attractive interactions and resembling the bright
soliton solutions. Thus, we can observe three regimes: (i) weakly-interacting
regime with no spontaneous symmetry breaking, (ii) moderate-interaction regime
and spontaneous breaking of the space translation symmetry of (3), and (iii)
strong-interaction regime where the spontaneous breaking of the symmetry
corresponds to the formation of bright soliton wavepackets of width
$\xi=2/|g_{0}(N-1)|<\sigma$, the width of the single-particle ground state in
the potential well Pethick and Smith (2002). While in the regime (iii) low-
energy excitations are related to collective excitations of the entire soliton
Weiss and Castin (2009); Sacha et al. (2009), the low-energy excitations in
the regime (ii) correspond to quantum depletion of the self-trapped BEC and
transfer of bosons to the other potential well Pethick and Smith (2002); Ziń
et al. (2008); Oleś et al. (2010).

(a)
(b) (c)

Figure 1:  (a) Solid black lines show the excitation spectrum of (3), i.e.
eigenenergies minus the ground state energy, for $N=9$ and $\lambda=1.5$.
Vertical dashed line indicates critical value of $g_{0}(N-1)$ for the quantum
phase transition to the discrete time crystal regime obtained for
$N\rightarrow\infty$ from the results presented in (c). Green lines are
related to the ground state level and the first excited level of the center of
mass of the system, i.e., the eigenvalues of $H_{\rm cm}$, which are doubly
degenerate because the bosons can be located in one of the wells of the
double-well potential. Red circles are exact quasi-energies of the kicked LL
model (cf. (4) with $T=\pi/31$) which are relevant to time crystal states and
which are perfectly reproduced by the spectrum of (3). (b) is an enlargement
of (a) in the vicinity of the critical point. (c) Log-log plots of the
difference, $\Delta E$ between the lowest eigenenergies of (3) vs. $N$ for
different fixed values of $g_{0}(N-1)$ as indicated in the figure. The
critical value of $g_{0}(N-1)\approx-0.18$ corresponds to an algebraic
decrease of $\Delta E$ with $N$ and it agrees with the two-mode prediction for
$N\rightarrow\infty$ SM . For weaker interactions $\Delta E$ approaches a
constant value while for stronger interactions, $\Delta E$ decreases
exponentially with $N$.

In order to diagonalize the Hamiltonian (3), we employ the eigenbasis of the
undriven Lieb-Liniger model which can be obtained with the help of the Bethe
ansatz Gaudin (2014); SM . Diagonalization of the Hamiltonian matrix yields
the energy spectrum of (3) which is depicted in Fig. 1 as a function of
$g_{0}(N-1)$ for $N=9$ and $\lambda=1.5$. Let us first focus on the strongest
interactions presented in Fig. 1(a) which are nearly in the regime (iii) where
the width $\xi$ of the bright soliton becomes smaller than the width $\sigma$
of the single-particle ground state in the potential well, i.e., when
$|g_{0}(N-1)|\gtrsim 2\sqrt{2}\lambda^{1/4}=3.13$ for $\lambda=1.5$. Without
the external potential ($\lambda=0$), the lowest energy solution within the
mean-field approach would represent a bright soliton,
$\phi_{0}(x-q)=\cosh^{-1}[(x-q)/\xi]/\sqrt{2\xi}$, which is occupied by all
bosons and which can be localized at any point $q$ on the ring Castin (2001).
In the presence of the double well potential ($\lambda\neq 0$) and for
$\xi\lesssim\sigma$, there are two possible ground state locations of the
soliton at the bottoms of the potential wells, i.e., $q=\pi/2$ or $3\pi/2$.
Low-energy excitations of the bright soliton correspond to excitations of its
center of mass which can be described by the Hamiltonian of a single-body of
mass $N$, $H_{\rm cm}=-\partial_{q}^{2}/2N+\lambda N\int
dx|\phi_{0}(x-q)|^{2}\cos(2x)$ Weiss and Castin (2009); Sacha et al. (2009).
Solving $H_{\rm cm}\chi_{n}(q)=E_{n}\chi_{n}(q)$, we obtain the spectrum of
the center of mass excitations and an estimate for the corresponding single-
particle densities, $\tilde{P}(x)\approx\int
dq|\chi_{n}(q)|^{2}|\phi_{0}(x-q)|^{2}$. Only one excited energy level that
belongs to such a spectrum can be seen in Fig. 1(a) because the regime (iii)
has not been reached. Actually, in order to obtain the spectrum of the center
of mass excitations we have not assumed that $\phi_{0}$ in $H_{\rm cm}$ is
given by the bright soliton profile but applied the self-consistency method SM
.

In the regime (iii), the low-energy spectrum of (3) corresponds to center of
mass excitations of the bright soliton. When we decrease the strength of the
attractive interactions, we lose such a single-body character of the low
energy spectrum which takes place for $\xi\gtrsim\sigma$. Then, low-energy
excitations lead to quantum depletion of a BEC localized in one of the
potential wells and transfer of bosons to the other well Ziń et al. (2008);
Oleś et al. (2010). This moderate-interaction regime corresponds to the self-
trapping of a BEC where bosons prefer to localize in one of the potential
wells but do not form a bound state like in the bright soliton case. The self-
trapping properties are observed in the ground and excited eigenstates up to
the so-called symmetry breaking edge, i.e., the corresponding excitation
energy is proportional to $N$ SM .

One can ask how weak the attractive interactions should be in order to recover
the space translation symmetry of (3)? The critical interaction strength for
the phase transition between the symmetry-broken and symmetry-preserving
phases can be estimated by means of the two-mode approach because the
interactions in this regime are weak and not able to modify the shape of the
single-particle wave-packets localized in the potential wells which are used
in the two-mode approximation Ziń et al. (2008); Sacha (2015a); SM . Indeed,
the two-mode prediction agrees with the numerical results shown in Fig. 1(c)
where at the critical value of $g_{0}(N-1)\approx-0.18$, the energy gap
between the lowest energy eigenstates decreases algebraically with $N$. For
stronger interactions the gap decreases exponentially while for weaker
interactions it approaches a constant value Sacha (2015a). Thus, if the
interactions are sufficiently strong and $N\rightarrow\infty$, the symmetry-
preserving eigenstates are degenerate and their superpositions form symmetry-
broken eigenstates which live forever.

Figure 2:  Single-particle probability densities corresponding to
superpositions of the two lowest energy eigenstates of (3) plotted in the lab
frame for different moments of time as indicated in the panels for $N=9$. In
each panel, collections of the densities for different interaction strength
$g_{0}(N-1)$ are presented. Solid white lines indicate the critical
interaction strength $g_{0}(N-1)\approx-0.18$ for the transition to the
discrete time crystal regime when $N\rightarrow\infty$. Densities above these
lines reveal a decay of the $2T$-periodic evolution due to tunneling — for
non-interacting bosons, complete tunneling takes place at $t\approx 302T$.
Densities below the solid white lines show discrete time crystal evolution for
$N\rightarrow\infty$. The other parameters are the same as in Fig. 1, i.e.,
$\lambda=1.5$ and $T=\pi/31$.

When we return to the laboratory frame, the symmetry-broken eigenstates of (3)
will evolve with a period of $2T$, demonstrating absolutely stable discrete
time crystals. In Fig. 2 we present the time evolution of superpositions of
the lowest symmetry-preserving eigenstates of (3) for different interaction
strengths. At $g_{0}(N-1)\approx-0.18$ and for $N\rightarrow\infty$, the
quantum phase transition to the discrete time crystal regime occurs and the
superpositions evolve with the period $2T$ — for $N=9$ the ground state
becomes practically degenerate for stronger interactions, i.e.,
$g_{0}(N-1)\approx-0.3$. At $g_{0}(N-1)$ around $-3$, there is a crossover to
the bright soliton regime where low-energy excitations have single-body
character because bosons form a bound bright soliton state Weiss and Castin
(2009); Sacha et al. (2009).

Having analyzed the system (1), we now switch to a more general driving when
the mapping of discrete time crystal states to the low-lying eigenstates of a
time-independent Hamiltonian is approximate only but which shows richer
interplay between different symmetries. Let us exchange the time-periodic
perturbation in the Hamiltonian (1) with

$H_{1}=\lambda
T\sum_{i=1}^{N}\cos(2x_{i})\sum_{m=-\infty}^{+\infty}\delta(t-mT),$ (4)

which describes periodic kicking of the particles with the period
$T=2\pi/\omega$. When we perform the same unitary transformation to the moving
frame as previously, i.e., $U_{t}$, and a similar shift of the momenta
$U_{p}$, then using the rotating-wave approximation or the Magnus expansion we
can obtain an effective Hamiltonian identical to (3), see SM .

We have already analyzed the system (3); thus, in the present case of the
time-periodic kicking (4), we have to demonstrate only that the low-energy
eigenstates of (3) reproduce well the relevant exact eigenstates of the
Floquet Hamiltonian which are also eigenstates of the Floquet evolution
operator. The Floquet evolution operator is the evolution operator of the
system over a single driving period which in the case of the time-periodic
kicking (4) reads

$U(T)=e^{-iH_{\rm LL}T}\;e^{-i\lambda T\sum_{i=1}^{N}\cos(2x_{i})}.$ (5)

This unitary operator can be diagonalized in the eigenbasis of the LL
Hamiltonian and knowing its eigenphases $\phi_{n}$, we can calculate the
quasi-energies of the system, $E_{n}=-\phi_{n}/T$. The obtained exact quasi-
energies, which are relevant to discrete time crystal states, are also shown
in Fig. 1 and they are perfectly reproduced by the low-lying spectrum of (3).
We have chosen the same $\lambda$ and $\omega$ as previously and consequently
we expect exactly the same spectrum as in the case of the Hamiltonian (1), and
indeed the obtained quasi-energies are indistinguishable in the plot.

Both the system (1) and the kicked LL model (4) can be described using the
Floquet formalism, where all physically relevant quasi-energies lie in a
single Floquet zone. For small values of $N$, Floquet states can be obtained
numerically. However, very high-order terms, which are neglected in numerical
calculations, may introduce tiny couplings between Floquet states with similar
quasi-energies that can lead to the decay of discrete time crystals as
$t\rightarrow\infty$, especially when we first take the $N\rightarrow\infty$
limit. In the case of the system (1), we know that quasi-energies
corresponding to the discrete time crystal can be unfolded and they correspond
to the low-energy spectrum of the time-independent Hamiltonian (3).
Furthermore, we have a guarantee that they are not coupled to any other
Floquet states in any order. In the case of the kicked LL model, such absolute
stability of the discrete time crystals can not be guaranteed. The absolutely
stable discrete time crystal has been designed so that the perturbation
contains only one harmonic responsible for the formation of the crystal. In
the case of the kicked LL model, there are many other harmonics whose
influence can be reduced but cannot be fully eliminated.

So far we have analyzed the Floquet states that break both discrete space and
time translational symmetry in the periodically kicked LL model by switching
to the frame moving with the frequency $\omega/2$ by means of the unitary
transformation $U_{t}$. Let us show that the periodically kicked LL model (4)
can also support states that spontaneously break the discrete space
translation symmetry without breaking the discrete time translation symmetry,
which does not exist in the corresponding LL model driven by the rotating
lattice. To investigate these states, we can also switch to the moving frame
but with the frequency $\omega$, and we still obtain the same effective
Hamiltonian (3) SM . Then, however, eigenstates of (3) that reveal spontaneous
breaking of the discrete space translation symmetry, do not break the discrete
time translation symmetry. Indeed, when we return to the laboratory frame, all
eigenstates of (3) evolve with the driving period $T$.

To summarize, we have analyzed interacting bosons on a ring with various
periodic perturbations which turns out to be a suitable system for realization
of discrete time crystals. Most importantly, this system can reveal absolutely
stable discrete time crystals which is demonstrated by an exact mapping of the
discrete time crystal states to low-lying eigenstates of a time-independent
Hamiltonian. The periodically driven bosons on a ring which we consider here
can also reveal big discrete time crystals which spontaneously evolve with a
period many times longer than the driving period Giergiel et al. (2018). These
kind of systems possess many sites in the time dimension and are suitable for
investigation of condensed matter phenomena in the time domain Sacha (2015b);
Hannaford and Sacha (2022), see also condenced matter in phase space crystals
Guo et al. (2013); Guo (2021).

This research was funded by the National Science Centre, Poland, Projects No.
2018/31/B/ST2/00349 (K.G.) and No. 2021/42/A/ST2/00017 (K.S.) and the
Australian Research Council, Project DP190100815. K.G. acknowledges the
support of the Polish National Agency for Academic Exchange Bekker Programme
(BPN/BEK/2021/1/00339).

## References

  * Wilczek (2012) F. Wilczek, Phys. Rev. Lett. 109, 160401 (2012).
  * Sacha (2015a) K. Sacha, Phys. Rev. A 91, 033617 (2015a), URL http://link.aps.org/doi/10.1103/PhysRevA.91.033617.
  * Khemani et al. (2016) V. Khemani, A. Lazarides, R. Moessner, and S. L. Sondhi, Phys. Rev. Lett. 116, 250401 (2016).
  * Else et al. (2016) D. V. Else, B. Bauer, and C. Nayak, Phys. Rev. Lett. 117, 090402 (2016).
  * Yao et al. (2017) N. Y. Yao, A. C. Potter, I.-D. Potirniche, and A. Vishwanath, Phys. Rev. Lett. 118, 030401 (2017).
  * Sacha (2020) K. Sacha, _Time Crystals_ (Springer International Publishing, Switzerland, Cham, 2020), ISBN 978-3-030-52523-1, URL https://doi.org/10.1007/978-3-030-52523-1.
  * D’Alessio and Rigol (2014) L. D’Alessio and M. Rigol, Phys. Rev. X 4, 041048 (2014).
  * Lazarides et al. (2014) A. Lazarides, A. Das, and R. Moessner, Phys. Rev. E 90, 012110 (2014).
  * Ponte et al. (2015a) P. Ponte, A. Chandran, Z. Papić, and D. A. Abanin, Ann. Phys. (Amsterdam) 353, 196 (2015a).
  * Rigol et al. (2007) M. Rigol, V. Dunjko, V. Yurovsky, and M. Olshanii, Phys. Rev. Lett. 98, 050405 (2007).
  * Santos and Rigol (2010a) L. F. Santos and M. Rigol, Phys. Rev. E 81, 036206 (2010a).
  * Rigol and Santos (2010) M. Rigol and L. F. Santos, Phys. Rev. A 82, 011604 (2010).
  * Santos and Rigol (2010b) L. F. Santos and M. Rigol, Phys. Rev. E 82, 031130 (2010b).
  * Cassidy et al. (2011) A. C. Cassidy, C. W. Clark, and M. Rigol, Phys. Rev. Lett. 106, 140405 (2011).
  * Caux and Essler (2013) J.-S. Caux and F. H. L. Essler, Phys. Rev. Lett. 110, 257203 (2013).
  * Nandkishore and Huse (2015) R. Nandkishore and D. A. Huse, Annu. Rev. Condens. Matter Phys. 6, 15 (2015).
  * Altman and Vosk (2015) E. Altman and R. Vosk, Annu. Rev. Condens. Matter Phys. 6, 383 (2015).
  * D’Alessio et al. (2016) L. D’Alessio, Y. Kafri, A. Polkovnikov, and M. Rigol, Adv. Phys. 65, 239 (2016).
  * Mierzejewski et al. (2018) M. Mierzejewski, M. Kozarzewski, and P. Prelovšek, Phys. Rev. B 97, 064204 (2018), URL https://link.aps.org/doi/10.1103/PhysRevB.97.064204.
  * Else et al. (2017) D. V. Else, B. Bauer, and C. Nayak, Phys. Rev. X 7, 011026 (2017).
  * Imbrie (2016) J. Z. Imbrie, Journal of Statistical Physics 163, 998 (2016), eprint 1403.7837.
  * Ponte et al. (2015b) P. Ponte, Z. Papić, F. M. C. Huveneers, and D. A. Abanin, Phys. Rev. Lett. 114, 140401 (2015b).
  * Lazarides et al. (2015) A. Lazarides, A. Das, and R. Moessner, Phys. Rev. Lett. 115, 030402 (2015).
  * Abanin et al. (2016) D. Abanin, W. D. Roeck, and F. Huveneers, Ann. Phys. 372, 1 (2016).
  * Wang et al. (2021a) J. Wang, P. Hannaford, and B. J. Dalton, New J. Phys. 23, 063012 (2021a).
  * Kuroś et al. (2020) A. Kuroś, R. Mukherjee, W. Golletz, F. Sauvage, K. Giergiel, F. Mintert, and K. Sacha, New Journal of Physics 22, 095001 (2020), URL https://doi.org/10.1088/1367-2630/abb03e.
  * Wang et al. (2021b) J. Wang, K. Sacha, P. Hannaford, and B. J. Dalton, Phys. Rev. A 104, 053327 (2021b), URL https://link.aps.org/doi/10.1103/PhysRevA.104.053327.
  * Zhang et al. (2017) J. Zhang, P. W. Hess, A. Kyprianidis, P. Becker, A. Lee, J. Smith, G. Pagano, I.-D. Potirniche, A. C. Potter, A. Vishwanath, et al., Nature 543, 217 (2017).
  * Choi et al. (2017) S. Choi, J. Choi, R. Landig, G. Kucsko, H. Zhou, J. Isoya, F. Jelezko, S. Onoda, H. Sumiya, V. Khemani, et al., Nature 543, 221 (2017).
  * Pal et al. (2018) S. Pal, N. Nishad, T. S. Mahesh, and G. J. Sreejith, Phys. Rev. Lett. 120, 180602 (2018), URL https://link.aps.org/doi/10.1103/PhysRevLett.120.180602.
  * Rovny et al. (2018) J. Rovny, R. L. Blum, and S. E. Barrett, Phys. Rev. Lett. 120, 180603 (2018), URL https://link.aps.org/doi/10.1103/PhysRevLett.120.180603.
  * Smits et al. (2018) J. Smits, L. Liao, H. T. C. Stoof, and P. van der Straten, Phys. Rev. Lett. 121, 185301 (2018), URL https://link.aps.org/doi/10.1103/PhysRevLett.121.185301.
  * Mi et al. (2022) X. Mi, M. Ippoliti, C. Quintana, A. Greene, Z. Chen, J. Gross, F. Arute, K. Arya, J. Atalaya, R. Babbush, et al., Nature 601, 531 (2022), ISSN 1476-4687, URL https://doi.org/10.1038/s41586-021-04257-w.
  * Randall et al. (2021) J. Randall, C. E. Bradley, F. V. van der Gronden, A. Galicia, M. H. Abobeih, M. Markham, D. J. Twitchen, F. Machado, N. Y. Yao, and T. H. Taminiau, Science 374, 1474 (2021), eprint 2107.00736.
  * Frey and Rachel (2022) P. Frey and S. Rachel, Science Advances 8, eabm7652 (2022), eprint https://www.science.org/doi/pdf/10.1126/sciadv.abm7652, URL https://www.science.org/doi/abs/10.1126/sciadv.abm7652.
  * Gupta et al. (2005) S. Gupta, K. W. Murch, K. L. Moore, T. P. Purdy, and D. M. Stamper-Kurn, Phys. Rev. Lett. 95, 143201 (2005).
  * Giergiel et al. (2018) K. Giergiel, A. Kosior, P. Hannaford, and K. Sacha, Phys. Rev. A 98, 013613 (2018), URL https://link.aps.org/doi/10.1103/PhysRevA.98.013613.
  * Hannaford and Sacha (2022) P. Hannaford and K. Sacha, Association of Asia Pacific Physical Societies Bulletin 32, 12 (2022), eprint 2202.05544.
  * (39) See Supplemental Material.
  * Pethick and Smith (2002) C. Pethick and H. Smith, _Bose-Eistein condensation in dilute gases_ (Cambridge University Press, Cambridge, England, 2002).
  * Jackson and Weinstein (2004) R. K. Jackson and M. I. Weinstein, Journal of Statistical Physics 116, 881 (2004), ISSN 1572-9613, URL https://doi.org/10.1023/B:JOSS.0000037238.94034.75.
  * Weiss and Castin (2009) C. Weiss and Y. Castin, Phys. Rev. Lett. 102, 010403 (2009), URL https://link.aps.org/doi/10.1103/PhysRevLett.102.010403.
  * Sacha et al. (2009) K. Sacha, C. A. Müller, D. Delande, and J. Zakrzewski, Phys. Rev. Lett. 103, 210402 (2009), URL https://link.aps.org/doi/10.1103/PhysRevLett.103.210402.
  * Ziń et al. (2008) P. Ziń, J. Chwedeńczuk, B. Oleś, K. Sacha, and M. Trippenbach, EPL (Europhysics Letters) 83, 64007 (2008), URL http://stacks.iop.org/0295-5075/83/i=6/a=64007.
  * Oleś et al. (2010) B. Oleś, P. Ziń, J. Chwedeńczuk, K. Sacha, and M. Trippenbach, Laser Physics 20, 671 (2010), ISSN 1555-6611, URL http://dx.doi.org/10.1134/S1054660X10050130.
  * Gaudin (2014) M. Gaudin, _The Bethe Wavefunction_ (Cambridge University Press, 2014).
  * Castin (2001) Y. Castin, in _Coherent atomic matter waves_ , edited by R. Kaiser, C. Westbrook, and F. David (Springer Berlin Heidelberg, Berlin, Heidelberg, 2001), pp. 1–136, ISBN 978-3-540-45338-3.
  * Sacha (2015b) K. Sacha, Sci. Rep. 5, 10787 (2015b), URL https://www.nature.com/articles/srep10787.
  * Guo et al. (2013) L. Guo, M. Marthaler, and G. Schön, Phys. Rev. Lett. 111, 205303 (2013), URL https://link.aps.org/doi/10.1103/PhysRevLett.111.205303.
  * Guo (2021) L. Guo, _Phase Space Crystals_ , 2053-2563 (IOP Publishing, 2021), ISBN 978-0-7503-3563-8, URL https://dx.doi.org/10.1088/978-0-7503-3563-8.
  * Korepin et al. (1993) V. E. Korepin, N. M. Bogoliubov, and A. G. Izergin, _Quantum Inverse Scattering Method and Correlation Functions_ (Cambridge University Press, Cambridge, UK, 1993).
  * Ribeiro et al. (2008) P. Ribeiro, J. Vidal, and R. Mosseri, Phys. Rev. E 78, 021106 (2008), URL https://link.aps.org/doi/10.1103/PhysRevE.78.021106.
  * Dziarmaga (2004) J. Dziarmaga, Phys. Rev. A 70, 063616 (2004), URL https://link.aps.org/doi/10.1103/PhysRevA.70.063616.
  * Blanes et al. (2010) S. Blanes, F. Casas, J. A. Oteo, and J. Ros, European Journal of Physics 31, 907 (2010), URL https://doi.org/10.1088%2F0143-0807%2F31%2F4%2F020.

Supplemental Material

In Sec. I of this Supplemental Material we present: details of the
transformations used in the Letter to describe the LL model driven by a
rotating lattice potential, analysis of consequences of breaking of space
translation symmetry in the moving frame on time translation symmetry in the
laboratory frame, two-mode analysis of the formation of discrete time
crystals, a short description of the Bethe ansatz approach, and analysis of
the bright soliton regime in the LL model.

In Sec. II we present: a derivation of the effective Hamiltonian of the
resonantly kicked LL model and a demonstration that in the kicked LL model
space translation symmetry can be spontaneously broken without breaking of
time translation symmetry.

##  I Lieb-Liniger model driven by a rotating lattice potential

Let us consider $N$ bosons with contact interactions on a ring with
circumference $2\pi$ which are periodically driven by a rotating lattice
potential $\cos(sx-\omega t)$ where $s$ is an integer number. The system
reduces to the Lieb-Linger (LL) Hamiltonian $H_{\rm LL}$ with an additional
time-periodic drive

$\displaystyle H$ $\displaystyle=$ $\displaystyle H_{\rm
LL}+\lambda\sum_{i=1}^{N}\cos(sx_{i}-\omega t),$ (S1) $\displaystyle H_{\rm
LL}$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{N}\frac{p_{i}^{2}}{2}+g_{0}\sum_{i<j}^{N}\delta(x_{i}-x_{j}),$
(S2)

where we use $R$ and $\hbar^{2}/mR^{2}$ for the length and energy units,
respectively, where $R$ is the ring radius and $m$ the mass of the bosons. The
interaction strength $g_{0}=2mR\omega_{\perp}a_{s}/\hbar$, where
$\omega_{\perp}$ is the frequency of the harmonic transverse confinement and
$a_{s}$ is the atomic s-wave scattering length. The system possesses discrete
time and space translation symmetries, i.e., the Hamiltonian (S1) does not
change if $t\rightarrow t+T$ (where $T=2\pi/\omega$) or all $x_{i}\rightarrow
x_{i}+2\pi/s$.

We have chosen the potential of the form of $\cos(sx-\omega t)$ to drive the
LL model but one can choose many different potentials $V(sx-\omega t)$ and
similar phenomena that we consider in the Letter can be realized.

#### Note on conventions

In the Bethe ansatz literature the following convention is often used Korepin
et al. (1993); Gaudin (2014)

$\displaystyle H^{\prime}$ $\displaystyle=$ $\displaystyle H^{\prime}_{\rm
LL}+\lambda^{\prime}\sum_{i=1}^{N}\cos(sx_{i}-\omega t),$ (S3) $\displaystyle
H^{\prime}_{\rm LL}$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{N}p_{i}^{2}+2c\sum_{i<j}^{N}\delta(x_{i}-x_{j}).$
(S4)

To switch to the convention used in (S1)-(S2), the following scaling has to be
applied

$\displaystyle c=g_{0},$ (S5) $\displaystyle\lambda^{\prime}=2\lambda,$ (S6)
$\displaystyle T^{\prime}=T/2,$ (S7) $\displaystyle H^{\prime}=2H.$ (S8)

###  I.1 Transformation to moving frame: Hamiltonian and quantum state

The model Hamiltonian (S1) becomes easier to investigate in the moving frame
of the rotating lattice, where we first perform a time-dependent unitary
transformation $U_{t}=\exp(i\sum_{j}p_{j}\omega t/s)$, leading to a shift in
the positions, $x_{i}\rightarrow x_{i}+\omega t/s$, and next a second time-
independent unitary transformation
$U_{p}=\exp\left(-i\sum_{j}x_{j}\omega/s\right)$, leading to a shift in the
momenta, $p_{i}\rightarrow p_{i}+\omega/s$. Under these transformations, the
time-dependent quantum Liouville-von Neumann equation in the lab frame
$i\partial_{t}\rho=[H,\rho]$ becomes
$i\partial_{t}\tilde{\rho}=[\tilde{H},\tilde{\rho}]$ in the moving frame,
where $\tilde{\rho}=U_{p}U_{t}\rho U_{t}^{\dagger}U_{p}^{\dagger}$ and $\rho$
are notations for density operators in the moving-frame and lab frame,
respectively. As will be shown below, we end up with the following exact time-
independent Hamiltonian

$\tilde{H}=\sum_{i=1}^{N}\left[\frac{p_{i}^{2}}{2}+\lambda\cos(sx_{i})\right]+g_{0}\sum_{i<j}^{N}\delta(x_{i}-x_{j}),$
(S9)

where a constant term has been omitted. This new Hamiltonian possesses the
same discrete space translation symmetry as the Hamiltonian (S1).

The quantities $\tilde{\rho}(\vec{x},\vec{x}^{\prime};t)$ and
$\rho(\vec{x},\vec{x}^{\prime};t)$ are matrix elements of $\tilde{\rho}$ and
$\rho$ for position eigenstates. Thus,
$\rho(\vec{x},\vec{x}^{\prime};t)=\left\langle\vec{x}\right|\,\rho_{t}\,\left|\vec{x}^{\prime}\right\rangle$,
with $\vec{x}\equiv\\{x_{1},x_{2},...,x_{N}\\}$ and the subscript $t$
specifies the density operator at time $t$. The diagonal elements of
$\rho(\vec{x},\vec{x}^{\prime};t)$ are related to the position single-particle
probability density in the laboratory frame $P(x,t)$ via

$P(x,t)=\int dx_{2}dx_{3}...dx_{n}\,\rho(\vec{x},\vec{x},;t),$ (S10)

where $\vec{x}=\\{x,x_{2},x_{3},...x_{N}\\}$.

To show that the moving frame density operator $\tilde{\rho}$ satisfies the
equation $i\partial_{t}\tilde{\rho}=[\tilde{H},\tilde{\rho}]$ involving the
transformed Hamiltonian $\tilde{H}$ we begin with the time-dependent
Schrödinger equation $(H-i\partial_{t})\,\left|\Psi\right\rangle=0$ for a
quantum state $\left|\Psi\right\rangle$ in the laboratory frame, and $H$ given
by (S1). The quantum state satisfies periodicity conditions for the ring of
length $2\pi$, thus
$\left\langle\vec{x}+2\pi|\Psi\right\rangle=\left\langle\vec{x}|\Psi\right\rangle$,
where $\vec{x}+2\pi\equiv\\{x_{1}+2\pi,x_{2}+2\pi,\ldots,x_{N}+2\pi\\}$.

We first apply the unitary transformation $U_{t}$. It is straight-forward
using the well-known identity
$\exp(S)\;\Omega\;\exp(-S)=\Omega+[S,\Omega]+(1/2)[S,[S,\Omega]]+\ldots$ and
the commutation rules $[x_{j^{\prime}},p_{j}]=i\delta_{j^{\prime}j}$ to show
that $U_{t}\,x_{i}\,U_{t}^{-1}=x_{i}+\omega t/s$ and
$U_{t}\,p_{i}\,U_{t}^{-1}=p_{i}$, so that

$U_{t}\,H\,U_{t}^{-1}=\sum_{i=1}^{N}\frac{p_{i}^{2}}{2}+g_{0}\sum_{i<j}^{N}\delta(x_{i}-x_{j})+\lambda\sum_{i=1}^{N}\cos(sx_{i}).$
(S11)

We then have

$U_{t}\,H\,U_{t}^{-1}\left(U_{t}\,\left|\Psi\right\rangle\right)-iU_{t}\,\frac{\partial}{\partial
t}\left(U_{t}^{-1}\;U_{t}\,\left|\Psi\right\rangle\right)=0.$ (S12)

Now

$\displaystyle\frac{\partial}{\partial
t}\left(U_{t}^{-1}U_{t}\,\left|\Psi\right\rangle\right)$ $\displaystyle=$
$\displaystyle\frac{\partial}{\partial
t}\left(U_{t}^{-1}\right)\times\left(U_{t}\,\left|\Psi\right\rangle\right)+\left(U_{t}^{-1}\right)\times\frac{\partial}{\partial
t}\left(U_{t}\,\left|\Psi\right\rangle\right)$ (S13) $\displaystyle=$
$\displaystyle(\frac{-i\omega}{s}\sum_{j}p_{j})\left(U_{t}^{-1}\right)\times\left(U_{t}\,\left|\Psi\right\rangle\right)+\left(U_{t}^{-1}\right)\times\frac{\partial}{\partial
t}\left(U_{t}\,\left|\Psi\right\rangle\right),$

since it is straight-forward to show that $\frac{\partial}{\partial
t}\left(U_{t}^{-1}\right)=(\frac{-i\omega}{s}\sum_{j}p_{j})\left(U_{t}^{-1}\right)$.
As $U_{t}$ commutes with $p_{j}$, substituting into (S12) then gives

$\displaystyle\left(\sum_{i=1}^{N}\frac{p_{i}^{2}}{2}+g_{0}\sum_{i<j}^{N}\delta(x_{i}-x_{j})+\lambda\sum_{i=1}^{N}\cos(sx_{i})-\frac{\omega}{s}\sum_{i}p_{i}-i\frac{\partial}{\partial
t}\right)\left(U_{t}\,\left|\Psi\right\rangle\right)$ $\displaystyle=$
$\displaystyle 0,$ (S14)

and consequently

$\displaystyle\left(\sum_{i=1}^{N}\frac{(p_{i}-\frac{\omega}{s})^{2}}{2}+g_{0}\sum_{i<j}^{N}\delta(x_{i}-x_{j})+\lambda\sum_{i=1}^{N}\cos(sx_{i})-\frac{N\omega^{2}}{2s^{2}}-i\frac{\partial}{\partial
t}\right)\left(U_{t}\,\left|\Psi\right\rangle\right)$ $\displaystyle=$
$\displaystyle 0.$ (S15)

Note the momentum off-set $\omega/s$ when only $U_{t}$ is applied.

Secondly, we apply the unitary transformation $U_{p}$. It is straight-forward
to show that $U_{p}\,x_{i}\,U_{p}^{-1}=x_{i}$ and
$U_{p}\,p_{i}\,U_{p}^{-1}=p_{i}+\omega/s$, so that from the last equation we
obtain

$\displaystyle
U_{p}\left(\sum_{i=1}^{N}\frac{(p_{i}-\frac{\omega}{s})^{2}}{2}+g_{0}\sum_{i<j}^{N}\delta(x_{i}-x_{j})+\lambda\sum_{i=1}^{N}\cos(sx_{i})-\frac{N\omega^{2}}{2s^{2}}\right)U_{p}^{-1}\left(U_{p}U_{t}\,\left|\Psi\right\rangle\right)-iU_{p}\frac{\partial}{\partial
t}\left(U_{p}^{-1}U_{p}U_{t}\,\left|\Psi\right\rangle\right)$ $\displaystyle=$
$\displaystyle 0,$ (S16)

and finally

$\displaystyle\left(\sum_{i=1}^{N}\frac{p_{i}{}^{2}}{2}+g_{0}\sum_{i<j}^{N}\delta(x_{i}-x_{j})+\lambda\sum_{i=1}^{N}\cos(sx_{i})-\frac{N\omega^{2}}{2s^{2}}-i\frac{\partial}{\partial
t}\right)\left(U_{p}U_{t}\,\left|\Psi\right\rangle\right)$ $\displaystyle=$
$\displaystyle 0,$ (S18)

noting that $U_{p}$ is not time dependent, so $\frac{\partial}{\partial
t}\left(U_{p}^{-1}\right)=0$. Apart from the constant
$-\frac{N\omega^{2}}{2s^{2}}$, the Hamiltonian in (S18) is the same as (S9).
The quantum state in the moving frame is
$|\tilde{\Psi}\rangle=U_{p}U_{t}\,\left|\Psi\right\rangle$.

The generalization to mixed states, which in the laboratory frame would be of
the form $\rho=\sum_{\Psi}P_{\Psi}|\Psi\rangle\langle\Psi|$, is obvious, and
results in the Liouville-von Neumann equation (described above) involving the
time-independent Hamiltonian (S9).

###  I.2 Position probability density in moving and laboratory frames

Since the Hamiltonian (S9) is time-independent, there exist well-defined
energy eigenstates and equilibrium states, so the density matrix elements in
the moving frame $\tilde{\rho}(\vec{x},\vec{x}^{\prime})$ can be time
independent. The corresponding matrix elements in the lab frame at time $t$
are given by

$\displaystyle\rho(\vec{x},\vec{x}^{\prime};t)$ $\displaystyle=$
$\displaystyle\left\langle\vec{x}\right|\,U_{t}^{\dagger}U_{p}^{\dagger}\;\tilde{\rho}\;U_{p}U_{t}\,\left|\vec{x}^{\prime}\right\rangle=\left\langle\vec{x}-\omega
t/s\right|U_{p}^{\dagger}\;\tilde{\rho}\;U_{p}\left|\vec{x}^{\prime}-\omega
t/s\right\rangle=\left\langle\vec{x}-\omega
t/s\right|\tilde{\rho}\left|\vec{x}^{\prime}-\omega t/s\right\rangle
e^{i\varphi(\vec{x},\vec{x}^{\prime})}$ (S19) $\displaystyle=$
$\displaystyle\tilde{\rho}\left(\vec{x}-\frac{\omega
t}{s},\vec{x}^{\prime}-\frac{\omega
t}{s}\right)e^{i\varphi(\vec{x},\vec{x}^{\prime})},$

where
$\varphi(\vec{x},\vec{x}^{\prime})=\sum_{i}(x_{i}-x_{i}^{\prime})\omega/s$ is
a phase that only depends on relative distances and originates from the
operator $U_{p}$. It does not depend on the quantum state. Using (S10) we then
see that the single-particle probability densities $P(x,t)$ and
$\tilde{P}(x,t)$ in the laboratory and moving frames are related via

$P(x,t)=\tilde{P}(x-\omega t/s).$ (S20)

Suppose in the moving frame a time-independent
$\tilde{\rho}(\vec{x},\vec{x}^{\prime})$ has a spatial periodicity $2\pi$,
i.e.,
$\tilde{\rho}\left(\vec{x}-2\pi,\vec{x}^{\prime}-2\pi\right)=\tilde{\rho}\left(\vec{x},\vec{x}^{\prime}\right)$.
Let us consider what happens in the laboratory frame at time $t+sT$

$\displaystyle\rho(\vec{x},\;\vec{x}^{\prime};\;t+sT)$ $\displaystyle=$
$\displaystyle\tilde{\rho}\left(\vec{x}-2\pi-\frac{\omega
t}{s},\;\vec{x}^{\prime}-2\pi-\frac{\omega
t}{s}\right)e^{i\varphi(\vec{x},\vec{x}^{\prime})}$ (S21) $\displaystyle=$
$\displaystyle\tilde{\rho}\left(\vec{x}-\frac{\omega
t}{s},\;\vec{x}^{\prime}-\frac{\omega
t}{s}\right)e^{i\varphi(\vec{x},\vec{x}^{\prime})}$ (S22) $\displaystyle=$
$\displaystyle\rho(\vec{x},\;\vec{x}^{\prime};\;t).$ (S23)

So, in the laboratory frame, the density matrix elements
$\rho(\vec{x},\vec{x}^{\prime};t)$ have time periodicity $sT$. From (S20) the
same applies to the single-particle probability density in the laboratory
frame.

There could of course be other possible space periodicities for the time
independent moving frame density matrix elements, such as $2\pi/s$, since
there are $s$ identical potential wells between $0$ and $2\pi$ and the
Hamiltonian (S9) is invariant under translation of all bosons in space by
$2\pi/s$. But suppose this is not the case, i.e.,

$\tilde{\rho}\left(\vec{x}-\frac{2\pi}{s},\;\vec{x}^{\prime}-\frac{2\pi}{s}\right)\neq\tilde{\rho}\left(\vec{x},\;\vec{x}^{\prime}\right).$
(S24)

Now consider what happens at time $t+T$ in the laboratory frame. Similarly to
before we would have

$\displaystyle\rho(\vec{x},\vec{x}^{\prime};\;t+T)$ $\displaystyle=$
$\displaystyle\tilde{\rho}\left(\vec{x}-\frac{\omega(t+T)}{s},\;\vec{x}^{\prime}-\frac{\omega(t+T)}{s}\right)e^{i\varphi(\vec{x},\vec{x}^{\prime})}$
(S25) $\displaystyle=$
$\displaystyle\tilde{\rho}\left(\vec{x}-\frac{2\pi}{s}-\frac{\omega
t}{s},\;\vec{x}^{\prime}-\frac{2\pi}{s}-\frac{\omega
t}{s}\right)e^{i\varphi(\vec{x},\vec{x}^{\prime})}$ (S26) $\displaystyle\neq$
$\displaystyle\tilde{\rho}\left(\vec{x}-\frac{\omega
t}{s},\;\vec{x}^{\prime}-\frac{\omega
t}{s}\right)e^{i\varphi(\vec{x},\vec{x}^{\prime})},$ (S27)

due to the previous inequality and consequently
$\rho(\vec{x},\vec{x}^{\prime};t+T)\neq\rho(\vec{x},\vec{x}^{\prime};t)$.

Thus, the failure of space periodicity $2\pi/s$ in the moving frame results in
the failure to have time periodicity $T$ in the laboratory frame. Conversely,
if in the moving frame $\tilde{\rho}(\vec{x},\vec{x}^{\prime})$ has a spatial
periodicity $2\pi/s$, then
$\rho(\vec{x},\vec{x}^{\prime};t+T)=\rho(\vec{x},\vec{x}^{\prime};t)$ and the
density matrix elements $\rho(\vec{x},\vec{x}^{\prime};t)$ then have time
periodicity $T$, the same as that of the drive in (S1). We have therefore
shown that spontaneous breaking of the discrete space translation symmetry of
(S9) corresponds to spontaneous breaking of the discrete time translation
symmetry of the original Hamiltonian (S1). That is, in the moving frame the
symmetry-broken ground state of (S9) corresponds to a state which evolves with
the period $sT$ in the laboratory frame and consequently describes a discrete
time crystal. For a space crystal we have a large number of regular
repetitions in space of an observable measured at any given time, whereas for
the time crystal with $s\gg 1$ we have a large number of regular repetitions
in time of an observable measured at any given position. For $s\gg 1$ we deal
with a big discrete time crystal Giergiel et al. (2018); Hannaford and Sacha
(2022).

###  I.3 Two-mode approximation

Let us focus on the $s=2$ case which corresponds to the double-well potential
in the Hamiltonian (S9). If we consider only a single particle, then the
lowest energy eigenstates of (S9) are symmetric and anti-symmetric
superpositions of two wavepackets, $w_{1,2}(x)$, localized in the potential
wells, i.e., $w_{1}(x)\pm w_{2}(x)$. The corresponding eigenenergies are
slightly split by

$J=-2\int_{0}^{2\pi}dx\;w_{2}^{*}(x)\left[\frac{p^{2}}{2}+\lambda\cos(2x)\right]w_{1}(x),$
(S28)

which is the tunneling amplitude — if a particle is initially prepared in
$w_{1}$ (or $w_{2}$) state, it tunels to the other potential well after time
$\pi/J$. If there are $N$ non-interacting bosons ($g_{0}=0$), then in the
Hilbert subspace spanned by the two modes $w_{1,2}$, the Hamiltonian (S9)
reduces to
$\tilde{H}=-J(\hat{a}_{1}^{\dagger}\hat{a}_{2}+\hat{a}_{2}^{\dagger}\hat{a}_{1})/2$,
where a constant term is omitted and $\hat{a}_{1,2}$ are the standard bosonic
anihilation operators. If the contact interactions are present ($g_{0}\neq
0$), then in the same Hilbert subsapce, the Hamiltonian reads Pethick and
Smith (2002); Sacha (2020)

$\tilde{H}\approx-\frac{J}{2}(\hat{a}_{1}^{\dagger}\hat{a}_{2}+\hat{a}_{2}^{\dagger}\hat{a}_{1})+\frac{U}{2}\left(\hat{a}_{1}^{\dagger}\hat{a}_{1}^{\dagger}\hat{a}_{1}\hat{a}_{1}+\hat{a}_{2}^{\dagger}\hat{a}_{2}^{\dagger}\hat{a}_{2}\hat{a}_{2}\right),$
(S29)

with $U=U_{11}-2U_{12}$ where

$U_{ij}=g_{0}\int_{0}^{2\pi}dx\;|w_{i}(x)|^{2}|w_{j}(x)|^{2}.$ (S30)

The two-mode approximation (S29) is valid provided the interaction energy per
particle is much smaller than the energy gap to the next single-particle
excited eigenstate which is not included in the Hilbert subspace spanned by
the two modes $w_{1,2}$, i.e., provided $NU\ll 2\sqrt{\lambda}$.

The Hamiltonian (S29) can be reduced to the form of the Lipkin-Meshkov-Glick
model by defining the spin operators,
$\hat{S}_{x}=(\hat{a}_{1}^{\dagger}\hat{a}_{2}+\hat{a}_{2}^{\dagger}\hat{a}_{1})/2$
and
$\hat{S}_{z}=(\hat{a}_{2}^{\dagger}\hat{a}_{2}-\hat{a}_{1}^{\dagger}\hat{a}_{1})/2$.
If we use $\hat{S}_{x,z}=\sum_{i=1}^{N}\sigma_{i}^{x,z}/2$, where
$\sigma_{i}^{x,z}$ are the Pauli matrices, the Lipkin-Meshokov-Glick model
takes the form of a system of $N$ spin-1/2 particles with all-to-all
interactions Ribeiro et al. (2008); Sacha (2020). In the Lipkin-Meshkov-Glick
model, there is a quantum phase transition between the paramagnetic and
ferromagnetic phases when the ferromagnetic interactions ($g_{0}<0$) between
the spins are sufficiently strong. In the two-mode Hamiltonian (S29) this
phase transition corresponds to self-trapping of bosons in one of the modes
$w_{1,2}$ (i.e., in one of the potential wells) and consequently spontaneous
breaking of the space translation symmetry of the Hamiltonian (S9) if the
attractive interactions ($g_{0}<0$) between bosons are sufficiently strong
Pethick and Smith (2002); Jackson and Weinstein (2004).

The critical interaction strength for the spontaneous breaking of the space
translation symmetry of the Hamiltonian (S29) can be obtained by applying the
mean-field theory where one looks for the ground state of the system in the
form of a product state (i.e., a Bose-Einstein condensate state)

$\tilde{\psi}(\vec{x})=\prod_{i=1}^{N}\left[a_{1}w_{1}(x_{i})+a_{2}w_{2}(x_{i})\right].$
(S31)

The mean-field ground state is determined by the minimal value of the energy
functional

$E=N\left[-\frac{J}{2}(a_{1}^{*}a_{2}+a_{2}^{*}a_{1})+\frac{U(N-1)}{2}\left(|a_{1}|^{4}+|a_{2}|^{4}\right)\right],$
(S32)

with the constraint $|a_{1}|^{2}+|a_{2}|^{2}=1$. For $\gamma>-1$, where

$\gamma=\frac{U(N-1)}{J},$ (S33)

the lowest energy state corresponds to the symmetric superposition of the
modes, i.e., $a_{1}=a_{2}=1/\sqrt{2}$ in (S31). However, at $\gamma=-1$, there
is a bifurcation and two degenerate symmetry-broken states are born which
possess the lowest energy when $\gamma<-1$

$a_{1,2}=\left(\frac{1\pm\sqrt{1-1/\gamma^{2}}}{2}\right)^{1/2}.$ (S34)

Thus, in the mean-field limit (i.e., when $N\rightarrow\infty$ but $\gamma=\rm
constant$ or equivalently $g_{0}N=\rm constant$), there are two degenerate
ground states which for $\gamma\ll-1$ describe a Bose-Einstein condensate
located in one of the wells of the double-well potential in (S9) with $s=2$.

For finite $N$, one can diagonalize the two-mode Hamiltonian (S29) and the
lowest energy eigenstates are not exactly degenerate and they obey the space
translation symmetry of the Hamiltonian (S9) — for $\gamma\ll-1$ they
correspond to

$|\tilde{\psi}_{\pm}\rangle\approx\frac{1}{\sqrt{2}}\left(|N,0\rangle\pm|0,N\rangle\right),$
(S35)

where $N$ bosons occupy either the mode $w_{1}$ (the state $|N,0\rangle$) or
the mode $w_{2}$ (the state $|0,N\rangle$). However, the energy splitting
between the corresponding eigenenergies, $E_{-}-E_{+}$, goes exponentially
quickly to zero with an increase of $N$ but fixed $\gamma$ Ziń et al. (2008);
Sacha (2015a). Thus, in the mean-field limit if the system is initially
prepared in one of the symmetry broken states, $|N,0\rangle$ or $|0,N\rangle$,
it lives in such a state infinitely long. Even if one was able to prepare in
the experiment one of the states (S35), any decoherence (e.g., measurement of
the position of a single boson) would lead to the collapse to one of the
symmetry-broken states, $|N,0\rangle$ or $|0,N\rangle$, because
$|\tilde{\psi}_{\pm}\rangle$ are actually Schrödinger cat-like states. The
collapse indicates the spontaneous symmetry breaking phenomenon Sacha (2015a).

We would like to stress that not only the two lowest eigenstates of the
Hamiltonian (S29) reveal spontaneous symmetry breaking but all eigenstates up
to the so-called symmetry broken edge, $E_{\rm edge}$, which can be estimated
by means of the mean-field theory Sacha (2020)

$E_{\rm
edge}-E_{\pm}=N\frac{J}{2}\left(|\gamma|+\frac{1}{|\gamma|}-2\right)>0.$ (S36)

It is also worth mentioning that low-energy states of the Hamiltonian (S9)
correspond to the $s:1$ resonant driving of the particles in the laboratory
frame. That is, low-energy states of (S9) are related to
$p_{i}\approx\omega/s$ in the laboratory frame and thus, in the classical
description, to the particles which are moving along the ring with a period
$s$ times longer than the driving period $T$. For $s=2$ considered in this
section this means the $2:1$ resonance.

###  I.4 Bethe ansatz approach

For $\lambda=0$, eigenstates of (S9) can be found analytically by means of the
Bethe ansatz approach Gaudin (2014). The basic idea of the approach is the
observation that between collisions of the particles, where certain
$x_{i}=x_{j}$, the particles are described by a product of momentum
eigenstates $e^{ik_{1}x_{1}}e^{ik_{2}x_{2}}\dots e^{ik_{N}x_{N}}$ which has to
be symmetrized if the particles are bosons. In order to fulfill boundary
conditions imposed by the contact interaction potential and the periodic
boundary conditions on a ring, the parameters $k_{i}$ (called quasi-momenta)
of the eigenstates have to satisfy the following Bethe equations

$\displaystyle 2\pi k_{i}+2\sum_{j\neq
i}^{N}\arctan\left(\frac{k_{i}-k_{j}}{g_{0}}\right)=2\pi I_{i},$ (S37)

where the $I_{i}$’s are arbitrary integers (half integers) for odd (even) $N$
which label the eigenstates of the system. The sums

$\displaystyle P_{\rm total}$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{N}k_{i},$ (S38) $\displaystyle E$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{N}\frac{k_{i}^{2}}{2}$ (S39)

are eigenvalues of the total momentum and energy of the unperturbed LL model,
respectively. The corresponding eigenstates read

$\tilde{\psi}_{\vec{k}}(\vec{x})\propto\prod_{n<m}\left[\frac{\partial}{\partial
x_{n}}-\frac{\partial}{\partial x_{m}}+g_{0}\;{\rm
sgn}(x_{n}-x_{m})\right]\det[e^{ik_{j}x_{i}}].$ (S40)

When the external potential is on, diagonalization of the Hamiltonian (S9) can
be performed in the eigenbasis of the unperturbed LL model (S40). The only
non-diagonal matrix elements correspond to the external potential

$\lambda\sum_{i=1}^{N}\langle\tilde{\psi}_{\vec{k}^{\prime}}|\cos(sx_{i})|\tilde{\psi}_{\vec{k}}\rangle.$
(S41)

In the remaining part of Sec. I we will focus on the case of the double-well
potential in (S9), i.e., $s=2$.

###  I.5 Bright soliton regime

Let us start with the unperturbed LL model ($\lambda=0$) and apply the mean-
field approximation in order to find the ground state of the Hamiltonian (S9).
That is, assume that in the moving frame we look for the lowest energy state
in the Hilbert space of product states

$\tilde{\psi}(\vec{x})=\prod_{i=1}^{N}\phi_{0}(x_{i}).$ (S42)

To find the minimal value of
$\langle\tilde{\psi}|\tilde{H}|\tilde{\psi}\rangle$ we have to minimize the
following energy functional

$E=N\int_{0}^{2\pi}dx\left[\frac{1}{2}|\partial_{x}\phi_{0}(x)|^{2}+\frac{g_{0}(N-1)}{2}|\phi_{0}(x)|^{4}\right],$
(S43)

with the constraint $\langle\phi_{0}|\phi_{0}\rangle=1$, which reduces to the
solution of the Gross-Pitaevskii equation

$-\frac{1}{2}\partial_{x}^{2}\phi_{0}(x)+g_{0}(N-1)|\phi_{0}(x)|^{2}\phi_{0}(x)=\mu\phi_{0}(x),$
(S44)

where $\mu$ is the chemical potential. For attractively interacting bosons on
a ring, there is a critical interaction strength $g_{0}(N-1)=\pi/2$ when the
uniform solution, $\phi_{0}=1/\sqrt{2\pi}$, becomes unstable and a bright
soliton is born. If $2/[g_{0}(N-1)]\ll 2\pi$, the soliton solution is well
approximated by the free space bright soliton wavefunction

$\phi_{0}(x-q)=\sqrt{\frac{1}{2\xi}}\cosh^{-1}\left[\frac{x-q}{\xi}\right],$
(S45)

where

$\xi=\frac{2}{g_{0}(N-1)},$ (S46)

and $q$ is the location of the solitonic center which can be arbitrary. The
corresponding chemical potential

$\mu=-\frac{g_{0}^{2}(N-1)^{2}}{8}.$ (S47)

The Bethe-ansatz approach (see Sec. I.4) allows one to obtain the exact many-
body ground state which, due to the continuous space translation symmetry of
the unperturbed LL Hamiltonian, corresponds to the uniform single-particle
probability density $\tilde{P}(x)$. However, the continuous space translation
symmetry is spontaneously broken (e.g., when we measure particles’ positions)
and the bright soliton probability profile $|\phi_{0}(x-q)|^{2}$ emerges in
the density of the detected bosons. The location $q$ of the soliton is
randomly determined in the detection process and can be any random position on
a ring, i.e., the space translation symmetry is broken spontaneously.

If the bosons are initially prepared in the bright soliton state, i.e., the
Bose-Einstein condensate state (S42) with $\phi_{0}$ given in (S45), the mean-
field approach predicts it is a stable stationary state. This prediction is
true only in the limit of $N\rightarrow\infty$ and $g_{0}\rightarrow 0$ with
$g_{0}N=\rm constant$. In the full many-body description, the center of mass
of the bosons is described by the free-particle Hamiltonian, $H_{\rm
cm}=P_{\rm cm}^{2}/(2N)$ Castin (2001), and even if initially the center of
mass is much better localized than the bright soliton width $\xi$, it starts
spreading leading to delocalization of the soliton. However, for a finite but
large $N$, the delocalization takes a very long time and has not been observed
experimentally yet.

In the presence of the double-well potential [$\lambda\neq 0$ and $s=2$ in
(S9)] we are interested in low-energy eigenstates of the system which can be
obtained numerically as described in Sec. I.4. In the bright soliton regime,
the character of low-energy excitations can be analyzed by employing a quantum
description to the center of mass of the soliton. If the double-well potential
is weak and the width $\xi$ of the bright soliton smaller than the width
$\sigma$ of the single-particle wave-packets $w_{1}$ or $w_{2}$ (see Sec.
I.3), the profile of the bright soliton is not much affected by the presence
of the potential but the center of mass behavior is. For $\xi<\sigma\approx
1/\sqrt{\Omega}$, where $\Omega=2\sqrt{\lambda}$ is the frequency of the
harmonic approximation of the potential wells, the center of mass of the
soliton can be described by the quantum Hamiltonian

$H_{\rm cm}=-\frac{1}{2N}\partial_{q}^{2}+\lambda
N\int_{0}^{2\pi}dx|\phi_{0}(x-q)|^{2}\cos(2x),$ (S48)

with $\phi_{0}$ given in (S45). The Hamiltonian (S48) can be derived by means
of the many-body approach or by an extended Bogoliubov approach where the
soliton position is treated in a non-perturbative way Dziarmaga (2004); Weiss
and Castin (2009); Sacha et al. (2009). The eigenequation $H_{\rm
cm}\chi_{n}(q)=E_{n}\chi_{n}(q)$ describes the ground state of the center of
mass of the system and its excitations. The latter are the low-energy
excitations of the system in the bright soliton regime where bosons form a
bound state and taking a single particle from such a bound state costs energy
given by modulus of the chemical potential (S47) which is larger than the
excitations of the center of mass of the system. Thus, in the bright soliton
regime the excitations of the system have the character of a single-body of
mass $N$ described by the Hamiltonian (S48).

Figure S1: The lowest peak is the exact single-particle probability density
$\tilde{P}(x)$ corresponding to one of the symmetry-broken ground states of
(S9) for $N=9$, $\lambda=1.5$ and $g_{0}(N-1)=-2.66$, cf. Fig. 1 in the
Letter. The middle-height curves are the deconvoluted $|\phi_{0}|^{2}$ profile
(solid line) and the analytical bright soliton profile (dash line)
corresponding to $g_{0}(N-1)=-2.66$. The highest peak is the ground-state
probability density, $|\chi_{0}|^{2}$, of $H_{\rm cm}$.

The strongest interactions analyzed in the Letter are not yet in the bright
soliton regime although they are close to this regime. Therefore, in order to
decribe excitations of the center of mass of bosons we do not assume that
$|\phi_{0}|^{2}$ in (S48) is given by the bright soliton solution (S45) but
find $|\phi_{0}|^{2}$ in a self-consistent way. That is, having the exact
single-particle probability density $\tilde{P}(x)$ correpsonding to one of the
symmetry-broken ground states of (S9) and starting with the ground state
$\chi_{0}(q)$ of the harmonic Hamiltonian, $H_{\rm
cm}\approx-\partial_{q}^{2}/(2N)+N\Omega^{2}q^{2}/2$, we can obtain
$|\phi_{0}|^{2}$ by solving the equation

$\tilde{P}(x)\approx\int_{0}^{2\pi}dx\;|\phi_{0}(x-q)|^{2}|\chi_{0}(q)|^{2},$
(S49)

which is an approximation of the single-particle probability density of the
system in the bright soliton state. By means of the Fourier transform we get
$|\phi_{0}|^{2}$ and substitute it in (S48) and find a new $\chi_{0}(q)$ and
continue the procedure until $|\phi_{0}|^{2}$ converges. A few iterations is
sufficient to obtain the final $|\phi_{0}|^{2}$ which well reproduces the
bright soliton profile (S45), see Fig. S1. Having $|\phi_{0}|^{2}$ we can find
the ground state energy and excited energy levels of (S48) which describe
excitations of the center of mass of the system.

##  II Periodically kicked Lieb-Liniger model

###  II.1 Dicrete time crystals

Let us consider the LL model which is periodically kicked in time,

$H=H_{\rm LL}+\lambda
T\sum_{i=1}^{N}\cos(sx_{i})\sum_{m=-\infty}^{+\infty}\delta(t-mT),$ (S50)

where $T=2\pi/\omega$, and perform the same unitary transformation to the
moving frame and the same shift of the momenta as in Sec. I.1

$\displaystyle x_{i}$ $\displaystyle\rightarrow$ $\displaystyle
x_{i}+\frac{\omega}{s}t,$ (S51) $\displaystyle p_{i}$
$\displaystyle\rightarrow$ $\displaystyle p_{i}+\frac{\omega}{s}.$ (S52)

The resulting Hamiltonian, with a constant term omitted, reads

$\displaystyle\tilde{H}$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{N}\left[\frac{p_{i}^{2}}{2}+\lambda\cos(sx_{i}+\omega
t)\sum_{q=-\infty}^{+\infty}e^{iq\omega t}\right]$ (S54)
$\displaystyle+g_{0}\sum_{i<j}^{N}\delta(x_{i}-x_{j}),$

where we have used

$\sum_{m=-\infty}^{+\infty}\delta(t-mT)=\frac{1}{T}\sum_{q=-\infty}^{+\infty}e^{iq\omega
t}.$ (S55)

Let us switch for a moment to the interaction picture,
$e^{i\tilde{H}_{0}t}\tilde{H}e^{-i\tilde{H}_{0}t}$, where
$\tilde{H}_{0}=\sum_{i=1}^{N}p_{i}^{2}/2$. Then, the part of the Hamiltonian
corresponding to the driving potential reads

$\displaystyle\sum_{i=1}^{N}\sum_{q=-\infty}^{+\infty}e^{i\tilde{H}_{0}t}\lambda\cos(sx_{i}+\omega
t)e^{iq\omega t}e^{-i\tilde{H}_{0}t}$ $\displaystyle=$
$\displaystyle\frac{\lambda}{2}\sum_{i=1}^{N}\sum_{q=-\infty}^{+\infty}e^{ip_{i}^{2}t/2}\left(e^{isx_{i}}e^{i(1+q)\omega
t}+e^{-isx_{i}}e^{-i(1-q)\omega t}\right)e^{-ip_{i}^{2}t/2}.$ (S56)

When we calculate matrix elements of (LABEL:intpic) in the momentum basis of
bosons, $|\vec{k}\rangle=|k_{1},\dots,k_{N}\rangle$, we obtain a sum of the
terms in the following form

$\frac{\lambda}{2}\left(e^{isx_{i}}e^{i(1+q)\omega
t}+e^{-isx_{i}}e^{-i(1-q)\omega t}\right)e^{i(k_{j}^{2}-k_{l}^{2})t/2}.$ (S58)

Low-energy states of the Hamiltonian $\tilde{H}_{0}$, i.e., states with
$k_{j}\approx 0$, correspond to resonant states in the laboratory frame, i.e.,
to states with $p_{i}\approx\omega/s$, cf. (S52). If the periodic driving
potential and interactions between bosons are weak, to describe the resonant
driving we may restrict to the Hilbert subspace with $k_{j}^{2}\ll\omega$.
Then, all time-dependent terms (S58) are quickly oscillating and can be
neglected except those with $q=\pm 1$. Such a rotating-wave approximation
leads to an effective Hamiltonian of the resonantly kicked LL model which in
the Schrödinger picture takes the form

$\displaystyle\tilde{H}\approx\sum_{i=1}^{N}\left[\frac{p_{i}^{2}}{2}+\lambda\cos(sx_{i})\right]+g_{0}\sum_{i<j}^{N}\delta(x_{i}-x_{j}),$
(S59)

identical to the Hamiltonian (S9) of the LL model driven by the rotating
lattice potential. Note, however, that in the present case it is an
approximate description only and not the exact description as in the case of
the rotating lattice potential. The neglected quickly oscillating terms do not
allow us to claim that the discrete time crystals in the kicked LL model,
which we describe in the Letter, are absolutely stable.

Another way to obtain the effective Hamiltonian (S59) is to apply the Magnus
expansion Blanes et al. (2010). Kicked models are special as the Magnus
expansion can be reduced to the Baker-Campbell-Hausdorff formula. The problem
of finding the Floquet Hamiltonian is that of finding a logarithm of the
unitary operator of the evolution of the system over a single period $T$. In
the moving frame we obtain

$\tilde{H}=\frac{i}{T}\log\left[e^{-i\tilde{H}_{LL}T}\;e^{-i\lambda
T\tilde{H}_{kick}}\right],$ (S60)

where

$\tilde{H}_{kick}=\sum_{i=1}^{N}\cos(sx_{i}).$ (S61)

Employing the Baker-Campbell-Hausdorff formula one gets the following leading
terms

$\tilde{H}\approx\tilde{H}_{LL}+\lambda\tilde{H}_{kick}-\frac{i\pi\lambda}{\omega}\;[\tilde{H}_{LL},\;\tilde{H}_{kick}],$
(S62)

where the last term can be omitted if the driving frequency $\omega$ is high,
i.e., if the resonant kicking of the bosons corresponds to highly excited
momentum states in the laboratory frame, $p_{i}\approx\omega/s$, cf. (S15).
This way we arrive at the same effective Hamiltonian as obtained within the
rotating-wave approximation, see (S59). We should stress that it is not
guaranted that the Magnus expansion converges and therefore there is no
guarantee that discrete time crystals in the kicked LL model are absolutely
stable.

###  II.2 Spontaneous breaking of the space translation symmetry without
breaking the time translation symmetry

If, instead of the transformations (S51)-(S52), we apply

$\displaystyle x_{i}$ $\displaystyle\rightarrow$ $\displaystyle x_{i}+\omega
t,$ (S63) $\displaystyle p_{i}$ $\displaystyle\rightarrow$ $\displaystyle
p_{i}+\omega,$ (S64)

then the kicked LL Hamiltonian takes the form

$\displaystyle\tilde{H}$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{N}\left[\frac{p_{i}^{2}}{2}+\lambda\cos(sx_{i}+s\omega
t)\sum_{q=-\infty}^{+\infty}e^{iq\omega t}\right]$ (S66)
$\displaystyle+g_{0}\sum_{i<j}^{N}\delta(x_{i}-x_{j}).$

In the interaction picture, the terms analogous to (S58) read

$\frac{\lambda}{2}\left(e^{isx_{i}}e^{i(s+q)\omega
t}+e^{-isx_{i}}e^{-i(s-q)\omega t}\right)e^{i(k_{j}^{2}-k_{l}^{2})t/2},$ (S67)

and only those with $q=\pm s$ are slowly varying in time and form the
effective Hamiltonian of the system within the rotating-wave approximation
which is identical to (S59). However, spontaneous breaking of the space
translation symmetry in the moving frame does not imply breaking of the time
translation symmetry in the laboratory frame because the relation between the
single-particle probability density in the laboratory and moving frames is the
following, cf. Sec. I.2

$P(x,t)=\tilde{P}(x-\omega t).$ (S68)

Thus, in the kicked LL model we can observe a richer interplay between time
and space translation symmetry breaking. This is possible because the time-
dependent perturbation in (S66) contains many different harmonics in time and
many different resonant drivings are possible Sacha (2015a); Giergiel et al.
(2018). This is contrary to the time-dependent perturbation in (S1) where
there is only one harmonic in time in the time-dependent perturbation.
Resonant driving and formation of resonant islands in classical single-
particle phase space, which support discrete time crystal states in the
quantum many-body description, are described in the next subsection.

###  II.3 Single-particle classical phase space

Figure S2: Stroboscopic maps (i.e. postion $x$ and momentum $p$ of the
particle after each period $T$ of the time evolution of the particle) in the
laboratory frame for the kicked single particle (left column) and the single
particle driven by the rotating lattice potential (right column) for $s=2$ and
$T=\pi/31$. Top row corresponds to $\lambda=1.5$ while bottom row to
$\lambda=40$.

Let us consider the single-particle counterpart of the kicked LL model

$H_{\rm single}=\frac{p^{2}}{2}+\lambda
T\cos(sx)\sum_{m=-\infty}^{+\infty}\delta(t-mT).$ (S69)

Classical motion of the particle can be analyzed with the help of a
stroboscopic map where after each evolution period, i.e., at $t_{n}=nT$,
points $(x_{n},p_{n})$ in classical phase space are plotted. In the kicked
system, $(x_{n},p_{n})$ are solutions of the following Hamilton equations

$\displaystyle x_{n+1}$ $\displaystyle=$ $\displaystyle x_{n}+p_{n}T,$ (S70)
$\displaystyle p_{n+1}$ $\displaystyle=$ $\displaystyle p_{n}+s\lambda
T\sin(sx_{n}).$ (S71)

The particle is moving on the ring, so we restrict to $x_{n}$ and $x_{n+1}$
modulo $2\pi$. Moreover, if we shift $p_{n}$ and $p_{n+1}$ by $2\pi/T$, then
(S71) does not change but also (S70) does not change because $x_{n}$ and
$x_{n+1}$ are defined modulo $2\pi$. Thus, the phase space structure is
periodic along the momentum axis with the period $2\pi/T$.

In Fig. S2 examples of the stroboscopic maps are shown for $s=2$, $T=\pi/31$
and two different values of the strength of the kicking, i.e., $\lambda=1.5$
(used in the Letter) and $\lambda=40$. In the latter case the kicking is so
strong that chaotic sea appears between the clearly visible resonance islands.

Figure S2 also presents similar stroboscopic maps but for the case when the
particle is driven by the rotating lattice potential, i.e., when the single-
particle Hamiltonian reads

$H_{\rm single}=\frac{p^{2}}{2}+\lambda\cos(2x-\omega t).$ (S72)

Then, there are only two islands corresponding to the $2:1$ resonant driving
of the particle, i.e., a particle starting at the center of one of the islands
at $t=0$ jumps to the other island at $t=T$ and returns to the initial island
at $t=2T$. In the quantum description, if the particle is initially prepared
in a wave-packet localized in one of the islands, it jumps between the islands
every period $T$ for some time but such a $2T$-periodic evolution is broken
due to tunneling process. Only if we switch to the many-body case and the
particles sufficiently strongly attract each other, the tunneling process
stops and a discrete time crystal is formed Sacha (2015a).

In the case of the rotating lattice potential (S72), there are only $2:1$
resonance islands and regardless of how strong the driving is, we do not see
irregular dynamics. Switching to the moving frame, $x\rightarrow x+\omega
t/2$, and performing the shift in the momentum, $p\rightarrow p+\omega/2$, we
obtain the time-independent Hamiltonian of the particle in the double well
potential. Low-energy eigenstates of the quantum many-body counterpart of such
a system correspond to the absolutely stable discrete time crystals which we
consider in the Letter.

In the case of the kicked particle (S69), we can see a few different
resonances in Fig. S2, i.e., co-rotating ($p=\pi/T$) and counter-rotating
($p=-\pi/T$) $2:1$ resonances and co-rotating ($p=2\pi/T$) and counter-
rotating ($p=-2\pi/T$) $1:1$ resonances. The $1:1$ resonances allow us to
realize spontaneous breaking of the discrete space translation symmetry (i.e.
the symmetry under the shift $x\rightarrow x+\pi$) without breaking the
discrete time translation symmetry if the many-body system is appropriately
driven as described in Sec. II.2. The presence of different resonances offers
richer interplay between time and space symmetry breaking but on the other
hand it does not allow one to claim that discrete time crystals in the kicked
LL model are absolutely stable. Indeed, apart from a resonant term there are
always non-resonant ones which cannot be eliminated and can lead to decay of
the crystal on a long time scale.