# Dynamics of spontaneous symmetry breaking in a space-time crystal

J.N. Stehouwer  Institute for Theoretical Physics and Center for Extreme
Matter and Emergent Phenomena, Utrecht University  J. Smits  Debye Institute
for Nanomaterials and Center for Extreme Matter and Emergent Phenomena,
Utrecht University  P. van der Straten  Debye Institute for Nanomaterials and
Center for Extreme Matter and Emergent Phenomena, Utrecht University  H.T.C.
Stoof  Institute for Theoretical Physics and Center for Extreme Matter and
Emergent Phenomena, Utrecht University

###### Abstract

We present the theory of spontaneous symmetry breaking (SSB) of discrete time
translations as recently realized in the space-time crystals of an atomic
Bose-Einstein condensate. The non-equilibrium physics related to such a
driven-dissipative system is discussed in both the Langevin as well as the
Fokker-Planck formulation. We consider a semi-classical and a fully quantum
approach, depending on the dissipation being either frequency independent or
linearly dependent on frequency, respectively. For both cases, the Langevin
equation and Fokker-Planck equation are derived, and the resulting equilibrium
distribution is studied. We also study the time evolution of the space-time
crystal and focus in particular on its formation and the associated dynamics
of the spontaneous breaking of a $Z_{2}$ symmetry out of the symmetry unbroken
phase, i.e., the equilibrium Bose-Einstein condensate before the periodic
drive is turned on. Finally, we compare our results with experiments and
conclude that our theory provides a solid foundation for the observations.

Diagram

##  I Introduction

Statistical physics [1, 2] has been extremely successful at describing both
the equilibrium properties as well as the non-equilibrium dynamics of many-
body systems. In the latter case this is not only true for small disturbances
of a system in equilibrium, that can be treated by linear-response theory and
therefore still probe equilibrium fluctuations of the system, but also for far
out-of-equilibrium dynamics such as quenches across a phase transition and the
ensuing relaxation towards a new equilibrium. At present there is also much
interest in studying systems which are constantly out of thermodynamic
equilibrium as many examples of this can be found in our mostly classical
everyday life, such as in fluids [3], cells [4] and suspensions [5].

In addition, the investigation of the non-equilibrium dynamics of quantum
many-body systems has obtained a large boost in the last decade from the
readily availability of ultra-cold atomic gases in research labs, see for
instance Refs. [6, 7, 8, 9, 10]. Ultra-cold atomic gases are especially suited
for non-equilibrium dynamics as essentially all their characteristics can be
manipulated and tuned almost at will. For instance, the magnetic or optical
trap in which the atoms are stored can be precisely engineered, allowing even
for a control over the dimensionality of the atomic gas that typically affects
the dynamics in a dramatic way [11]. In this manner also a periodic drive can
be easily applied such that the system cannot behave as it would do in the
standard equilibrium situation as the continuous time-translation symmetry is
broken into a discrete symmetry only. Nevertheless a non-equilibrium steady
state can occur, in which all the energy gain that is provided by the drive is
exactly dissipated away from the system [12].

An especially interesting example of such driven-dissipative systems are so-
called (discrete) time crystals: systems that have spontaneously broken the
discrete time-translational symmetry of the periodic drive [13, 14, 15]. Time
crystals have now been observed experimentally by multiple groups around the
world [16, 17, 18, 19, 20, 21, 22], but our focus in this paper lies on the
experimental realization of a space-time crystal in a Bose-Einstein condensate
(BEC) of sodium atoms as described in Refs. [23, 24]. Notice that only in this
case a space-time crystal is formed, which means that on top of the periodic
structure in time, with twice the period of the drive, it also possesses a
crystalline structure in space. Hence, both the discrete time-translation
symmetry as well as the continuous translation symmetry are broken
spontaneously at the same time.

In the experiment of interest one starts from a cigar-shaped, and therefore
axially symmetric gas, which is evaporatively cooled to such low temperatures
that an almost complete Bose-Einstein condensation has occurred. By suddenly
changing the radial trap frequency one excites a long-lived radial breathing
mode that effectively implements a periodic drive for the axial collective
modes of the Bose-Einstein condensate. As a response to this drive exactly one
crystalline mode, i.e., a standing wave, is excited in the axial direction. As
already mentioned, this mode also has a frequency of precisely half the drive
frequency and on the basis of its spontaneously broken symmetries it thus
corresponds precisely to a space-time crystal. Most interestingly for our
purposes, there exists a (hidden) Ising-like symmetry related to the time
evolution of the amplitude of the space-time crystal being roughly speaking
either in phase or out of phase with the drive. More precisely, the temporal
phase difference of a single realization of the experiment can be either
$\phi$ or $\phi+\pi$, where the phase $\phi$ depends on the microscopic
details of the experiment and is therefore not universal. This is the
spontaneous symmetry breaking of the $Z_{2}$ symmetry that has been observed
recently in Ref. [25] and that will be investigated theoretically in great
detail in this paper.

The layout of the rest of the paper is as follows. In Sec. II we explain the
theoretical methods that are needed to describe the spontaneous symmetry
breaking observed in the experiments. In particular, we will be using a
functional formulation of the Schwinger-Keldysh formalism that is very
convenient for our purposes. In the end it will lead us to the exact Langevin
and Fokker-Planck equations that accurately describe the full quantum dynamics
of both the amplitude and the phase of the space-time crystal, including the
nonlinear dissipation from the small thermal cloud in the atomic gas. In Sec.
III, we then split our discussion into two parts, depending on the assumptions
we make about the frequency dependence of the dissipation. We consider two
cases that are both allowed by the generalized fluctuation-dissipation theorem
that we derived in Sec. II, namely frequency-independent dissipation and a
dissipation rate that depends linearly on the frequency. More physically, the
validity of these two cases depends on the time scale of interest in the
experiment, describing either prethermalization or thermalization,
respectively. We investigate the Langevin and Fokker-Planck equations in each
of the two approximations. In particular, we will look at the equilibrium
solutions, but also discuss numerical solutions of the Fokker-Planck equation
that explicitly show the dynamics of the spontaneous symmetry breaking taking
place in the Bose-Einstein condensate after the drive is turned on and the
space-time crystal is formed out of quantum and thermal fluctuations. After
that, we are in a position to also compare our results with experiments in
Sec. IV. In Sec. V we end the paper with conclusions and an outlook.

##  II Schwinger-Keldysh set-up

The starting point for this paper comes from the results of Ref. [24]. In that
paper, the average dynamics of the space-time crystals were found to be
described very accurately by the Hamiltonian

$\hat{H}=-\hbar\delta\hat{a}^{\dagger}\hat{a}+\frac{\hbar\omega_{D}A_{D}}{8}(\hat{a}^{\dagger}\hat{a}^{\dagger}+\hat{a}\hat{a})+\frac{\hbar
g}{2}\hat{a}^{\dagger}\hat{a}^{\dagger}\hat{a}\hat{a},$ (1)

where $\delta$ is the detuning from resonance, $\omega_{D}$ is the driving
frequency, $A_{D}$ is the relative driving amplitude and $\hat{a}^{(\dagger)}$
is the annihilation (creation) operator of a quantum in the resonant axial
mode. Note that the detuning enters this Hamiltonian, as it describes the slow
dynamics in the frame rotating along with the drive. It is in this rotating
frame that the Hamiltonian becomes in a good approximation time independent
and that prethermalization can therefore occur. Furthermore, there is the
complex interaction parameter
$g=g^{\prime}+ig^{\prime\prime}=|g|e^{i\phi_{g}}$, which was introduced
phenomenologically to counter the exponential growth of the space-time crystal
due to the presence of the drive and to allow for a relaxation into a long-
lived prethermal state as observed experimentally. Phenomenologically, this
nonlinear interaction term leads to an excellent agreement with the (average)
growth and relaxation dynamics observed experimentally. In particular, the
nonlinear dissipation due to $g^{\prime\prime}<0$ is required, and the more
common linear dissipation does not work in this case. Since the Hamiltonian
has an imaginary part and is thus non-hermitian, we expect that also noise
will arise because of the fluctuation-dissipation theorem. Therefore, we will
expand and support the model of Ref. [24] by microscopically deriving this
imaginary piece and its associated noise. We will do this via a non-linear
coupling with a heat bath, resulting in both the imaginary part of the
Hamiltonian, as well as the noise that satisfies a generalized fluctuation-
dissipation theorem.

An in our opinion elegant way to develop our theory is by using the functional
techniques of the Schwinger-Keldysh formalism. With that, we will derive an
effective action that precisely corresponds to the Hamiltonian of Eq. (1) and
that automatically also incorporates the fluctuations generated by the non-
hermitian term in the Hamiltonian. As the dynamics of the drive is fully
coherent, we focus for simplicity first on the derivation of the nonlinear
interaction term by considering the case of $A_{D}=0$ and $\omega_{D}=0$,
i.e., without the drive. After that has been achieved, the effect of the drive
can still be easily incorporated exactly into the resulting Langevin and
Fokker-Planck equations. We hope that this makes the following discussion more
clear by avoiding unnecessary long expressions.

Let us first derive the microscopic action corresponding to the non-linear
coupling with the heat bath. The probability distribution $P[a^{*},a;t]$ for
the space-time crystal to be in the coherent state $|a\rangle$ of the
annihilation operator $\hat{a}$ at time $t$ can be expressed in terms of the
matrix elements $\langle a;t|a_{0};t_{0}\rangle$ and its complex conjugate
$\langle a;t|a_{0};t_{0}\rangle^{*}=\langle a_{0};t_{0}|a;t\rangle$, where
$|a;t\rangle$ and $|a_{0};t_{0}\rangle$ are the coherent states at the present
time $t$ and at the initial time $t_{0}$, respectively. For these matrix
elements, the ‘path’-integral expressions going through all possible field
configurations are calculated, e.g. in Refs. [26, 27]. Because of the fact
that the complex conjugate has an opposite time ordering, the actions
corresponding to these elements will have opposite evolutions in time. This
leads to the introduction of the Schwinger-Keldysh contour $\mathcal{C}$,
which obeys
$\int_{\mathcal{C}}dt^{\prime}=\int_{t_{0}}^{t}dt_{+}+\int_{t}^{t_{0}}dt_{-}$.
For our purposes of calculation equal-time correlation functions, we are
allowed to deform the contour such that it encircles the complete time axis by
taking the limits $t_{0}\to-\infty$ and $t\to\infty$. For more details, we
refer to Ref. [28].

To reproduce the non-hermitian interaction term in the Hamiltonian we first of
all make use of the fact that for ultra-cold atoms a possible frequency
dependence of the interaction comes mostly from the dependence of the
interatomic collisions on the center-of-mass energy. Therefore, we imagine a
decoupling of the interaction in the particle-particle channel and consider a
linear interaction of the heat bath with $\psi(t)\psi(t)$, where $\psi(t)$
refers to the field of the space-time crystal. This implies that two quanta
from the resonant axial mode can combine into a single quantum of the heat
bath and, because of time-reversal symmetry of the microscopic action, that
one quantum of the heat bath can also decay into two quanta of the resonant
axial mode. Since the quadratic piece of the Hamiltonian is not affected by
these scattering processes with the heat bath, we can in more detail write for
the microscopic action that is the starting point of our discussion

$\displaystyle S[\psi^{*},\psi,\psi^{*}_{R},\psi_{R}]$
$\displaystyle=-\frac{1}{\sqrt{V}}\sum_{\mathbf{k}}\int_{\mathcal{C}}dt\bigg{\\{}t(\mathbf{k})\psi(t)\psi(t)\psi^{*}_{\mathbf{k}}(t)+t^{*}(\mathbf{k})\psi_{\mathbf{k}}(t)\psi^{*}(t)\psi^{*}(t)\bigg{\\}}$
(2)
$\displaystyle+\int_{\mathcal{C}}dt\psi^{*}(t)\left(i\hbar\frac{\partial}{\partial
t}+\hbar\delta\right)\psi(t)+\sum_{\mathbf{k}}\int_{\mathcal{C}}dt\psi^{*}_{\mathbf{k}}(t)\left(i\hbar\frac{\partial}{\partial
t}-\epsilon(\mathbf{k})\right)\psi_{\mathbf{k}}(t),$

where
$\psi_{R}(\mathbf{x},t)=\sum_{\mathbf{k}}\psi_{\mathbf{k}}e^{i\mathbf{k}\cdot\mathbf{x}}/\sqrt{V}$
is the reservoir field describing the homogeneous heat bath with volume $V$,
$t({\bf k})$ are the ‘tunneling’ matrix elements or scattering amplitudes
between the space-time crystal and the heat bath, $\epsilon({\bf k})$ denotes
the energy of a reservoir state with wave vector ${\bf k}$, and $\mathcal{C}$
is the deformed Schwinger-Keldysh contour that was introduced above. Note that
since we are dealing in Eq. (1) with a single-mode Hamiltonian for the space-
time crystal, there is no position dependence in the field $\psi(t)$ in
contrast to the reservoir field $\psi_{R}({\bf x},t)$, which in the
thermodynamic limit $V\rightarrow\infty$ describes the continuum required for
the occurrence of dissipation and the resulting ‘arrow of time’.

The next step is to integrate out the heat bath since the action is quadratic
in $\psi_{\mathbf{k}}$. This gives the desired Schwinger-Keldysh action for
the full quantum dynamics of the space-time crystal as

$\displaystyle S[\psi^{*},\psi]$ $\displaystyle=$
$\displaystyle\int_{\mathcal{C}}dt\psi^{*}(t)\left(i\hbar\frac{\partial}{\partial
t}+\hbar\delta\right)\psi(t)$ $\displaystyle-$
$\displaystyle\int_{\mathcal{C}}dt\int_{\mathcal{C}}dt^{\prime}\psi^{*}(t)\psi^{*}(t)\frac{\hbar
g(t,t^{\prime})}{2}\psi(t^{\prime})\psi(t^{\prime}).$

Here we introduced the very important two-point function on the Schwinger-
Keldysh contour $\hbar g(t,t^{\prime})=(2/\hbar
V)\sum_{\mathbf{k}}t^{*}(\mathbf{k})G(\mathbf{k};t,t^{\prime})t(\mathbf{k})$,
which determines the coupling strength of the effective interaction.
Furthermore, the reservoir Green’s function $G(\mathbf{k};t,t^{\prime})$
satisfies $(i\hbar\partial/\partial
t-\epsilon(\mathbf{k}))G(\mathbf{k};t,t^{\prime})=\hbar\delta(t,t^{\prime})$.
A Feynman diagram corresponding to the effective interaction is shown in Fig.
1. Note in particular the appearance of a fourth-order term in $\psi$ in our
effective action, which would become a quadratic selfenergy term [28] if we
would have used a linear coupling with the heat bath as is common in the
Caldeira-Leggett-like models of macroscopic quantum mechanics [29].

(30,30) thick i1,i2 o1,o2 fermioni1,v1 fermioni2,v1fermionv2,o1 fermionv2,o2
photonv1,v2 v1,v2 $t^{\prime}$v1 $t$v2

Figure 1: Feynman diagram of the effective interaction in the space-time
crystal. The wiggly line corresponds to the Green’s function or propagator of
the reservoir.

We have almost completed the derivation of the effective non-equilibrium
action. The last step is projecting the two (forward and backward) branches of
the Schwinger-Keldysh contour onto the real axis to obtain the real-time
dynamics. This is achieved by making the transformation
$\psi(t_{\pm})=a(t)\pm\xi(t)/2$, where $a(t)$ represent the average or
‘classical’ part of the $\psi$ field and $\xi(t)$ the fluctuations, as we will
see explicitly in a moment. Because of the fourth-order nature of the
effective interaction, we can in principle get fluctuations up to fourth order
in $\xi$, instead of only up to quadratic order as in the Caldeira-Leggett-
like models with a linear coupling to the heat bath [28]. Nonetheless, we can
still work out the action. Up to second order in the fluctuations, we find

$\displaystyle S^{\rm eff}[a^{*},a;\xi^{*},\xi]$ $\displaystyle=\int
dt\bigg{\\{}a^{*}(t)\bigg{(}i\hbar\frac{\partial}{\partial
t}+\hbar\delta\bigg{)}\xi(t)+\xi^{*}(t)\bigg{(}i\hbar\frac{\partial}{\partial
t}+\hbar\delta\bigg{)}a(t)\bigg{\\}}$ (4) $\displaystyle-\int dt\int
dt^{\prime}\bigg{\\{}a^{*}(t)a^{*}(t)\hbar
g^{(-)}(t-t^{\prime})a(t^{\prime})\xi(t^{\prime})+\xi^{*}(t)a^{*}(t)\hbar
g^{(+)}(t-t^{\prime})a(t^{\prime})a(t^{\prime})\bigg{\\}}$ $\displaystyle-\int
dt\int dt^{\prime}\xi^{*}(t)a^{*}(t)\hbar
g^{K}(t-t^{\prime})a(t^{\prime})\xi(t^{\prime}),$

where the retarded, advanced, and Keldysh functions, $g^{(\pm)}(t-t^{\prime})$
and $g^{K}(t-t^{\prime})$ respectively, are related to the analytic (known as
bigger and lesser) pieces of $g(t,t^{\prime})$ via

$g^{(\pm)}(t-t^{\prime})=\pm\theta(\pm(t-t^{\prime}))\left(g^{>}(t-t^{\prime})-g^{<}(t-t^{\prime})\right),$
(5)

and

$g^{K}(t-t^{\prime})=g^{>}(t-t^{\prime})+g^{<}(t-t^{\prime}),$ (6)

with $\theta(t)$ the usual Heaviside step function. Using these definitions
and our explicit expression for the two-point function on the Schwinger-
Keldysh contour $g(t,t^{\prime})$ we can now prove that the Fourier transforms
of the above retarded, advanced, and Keldysh functions are related by the
famous fluctuation-dissipation theorem

$g^{K}(\omega)=\left(g^{(+)}(\omega)-g^{(-)}(\omega)\right)(1+2N(\omega)),$
(7)

where we introduced the Bose or Planck distribution
$N(\omega)=1/(e^{\hbar\omega/k_{B}T}-1)$ and the temperature $T$ of the heat
bath. This crucial relation will ultimately ensure that our probability
distribution $P[a^{*},a;t]$ relaxes to the correct physical equilibrium.

The effective action also has a third-order contribution, namely

$\displaystyle-\frac{1}{8}\int dt\int dt^{\prime}$
$\displaystyle\bigg{(}\xi^{*}(t)\xi^{*}(t)\hbar
g^{(-)}(t-t^{\prime})a(t^{\prime})\xi(t^{\prime})$
$\displaystyle+\xi^{*}(t)a^{*}(t)\hbar
g^{(+)}(t-t^{\prime})\xi(t^{\prime})\xi(t^{\prime})\bigg{)},$

whereas the fourth-order term turns out to be exactly zero. In the remainder
of this paper, we will ignore these third-order terms in the fluctuations,
because they lead to quantum corrections to the equations of motion resulting
from the Hamiltonian in Eq. (1), that we already know to be unimportant for an
understanding of the experiments of interest to us here. In the literature the
neglect of these higher-order fluctuations is sometimes also known as the
truncated Wigner approximation [30, 31].

There are now two ways to proceed and extract the non-equilibrium dynamics
from the effective action in Eq. (4). First, we can integrate out the
fluctuations and arrive at an effective action for the field $a(t)$ alone.
Quantizing this action and writing down the associated Schrödinger equation
leads then to the Fokker-Planck equation for the probability distribution
$P[a^{*},a;t]$. Second, we can first derive the Langevin equation for $a(t)$
and then obtain the desired Fokker-Planck equation from that. We will follow
the latter approach here, because it avoids some operator-ordering problems
that arise in the former approach. But for consistency reasons and as a check
on our calculations, we illustrate also the direct derivation of the Fokker-
Planck equation in the following section.

To obtain the Langevin equation we decouple the quadratic term in the
fluctuations by a Hubbard-Stratonovich transformation, see e.g. Ref. [26]. In
this case this can most conveniently be achieved by multiplying the integrand
of our functional integral over $a(t)$ and $\xi(t)$ by the Gaussian integral

$1=\int d[\eta^{*}]d[\eta]\exp\left(\frac{iS^{\rm
eff}[\eta^{*},\eta]}{\hbar}\right),$ (8)

where we take

$\displaystyle S^{\rm eff}[\eta^{*},\eta]=\int dt\int
dt^{\prime}\bigg{(}\eta^{*}(t)-\int
dt^{\prime\prime}\xi^{*}(t^{\prime\prime})a^{*}(t^{\prime\prime})\hbar
g^{K}(t^{\prime\prime}-t)\bigg{)}\frac{{g^{K}}^{-1}(t-t^{\prime})}{\hbar}\bigg{(}\eta(t^{\prime})-\int
dt^{\prime\prime}\hbar
g^{K}(t^{\prime}-t^{\prime\prime})\xi(t^{\prime\prime})a(t^{\prime\prime})\bigg{)}.$
(9)

This will by construction precisely cancel out the quadratic term in $\xi$ of
our effective action, and we thus end up with an action that is completely
linear in $\xi$. Adding the above $S^{\rm eff}[\eta^{*},\eta]$ to our
effective action in Eq. (4) indeed leads to

$\displaystyle S^{\rm eff}[a^{*},a,\xi^{*},\xi,\eta^{*},\eta]$
$\displaystyle=\int dt\bigg{\\{}a^{*}(t)\left(i\hbar\frac{\partial}{\partial
t}+\hbar\delta\right)-\int dt^{\prime}a^{*}(t^{\prime})a^{*}(t^{\prime})\hbar
g^{(-)}(t^{\prime}-t)a(t)-a(t)\eta^{*}(t)\bigg{\\}}\xi(t)$ (10)
$\displaystyle+\int dt\xi^{*}(t)\bigg{\\{}\left(i\hbar\frac{\partial}{\partial
t}+\hbar\delta\right)a(t)-\int dt^{\prime}a^{*}(t)\hbar
g^{(+)}(t-t^{\prime})a(t^{\prime})a(t^{\prime})-\eta(t)a^{*}(t)\bigg{\\}}$
$\displaystyle+\int dt\int
dt^{\prime}\eta^{*}(t)\frac{{g^{K}}^{-1}(t-t^{\prime})}{\hbar}\eta(t^{\prime}).$

Because the effective action has become linear in the fluctuations after the
Hubbard-Stratonovich transformation, integrating out $\xi(t)$ now leads to
$\delta$ functionals in a similar manner as in the well-known identity $\int
dke^{ikx}=2\pi\delta(x)$. As a result we can immediately read off from Eq.
(10) that our Langevin equation becomes

$\displaystyle i\hbar\frac{\partial}{\partial t}a(t)$ $\displaystyle=$
$\displaystyle-\hbar\delta a(t)+\int dt^{\prime}a^{*}(t)\hbar
g^{(+)}(t-t^{\prime})a(t^{\prime})a(t^{\prime})$ (11) $\displaystyle+$
$\displaystyle\eta(t)a^{*}(t),$

together with the complex conjugate equation for $a^{*}(t)$. Furthermore, we
also conclude from our effective action in Eq. (10) that the noise
correlations obey

$\displaystyle\langle\eta(t)\eta(t^{\prime})\rangle=\langle\eta^{*}(t)\eta^{*}(t^{\prime})\rangle=0,$
(12)
$\displaystyle\langle\eta^{*}(t)\eta(t^{\prime})\rangle=i\hbar^{2}g^{K}(t-t^{\prime}).$
(13)

Note in particular that the Langevin equation contains multiplicative noise,
due to the fact that the dissipation occurs in the non-linear term in our
case. As always in physics, this multiplicative noise must be interpreted as
Stratonovich noise that leads also to noise-induced drift terms as we show
next.

##  III Frequency dependence of the dissipation

At this point, we will split our discussion up into two parts in order to
further analyze the general Langevin equations given by Eq. (11). We will do
this by calculating the equations of motion for the averages $\langle
a(t)\rangle$ and $\langle a^{*}(t)\rangle$ and the fluctuations $\langle
a^{*}(t)a(t)\rangle$, and then inferring the Fokker-Planck equation from these
equations of motion. We can, as a first approximation, assume that the
dissipation doesn’t have any frequency dependence, meaning that the
interactions act locally in time. Our second method will be the situation
where we allow for a linear frequency dependence. These two options will be
our two different approaches.

The goal of our first approximation, which is a semi-classical approach, is to
describe the relatively short time scales explored by the experiments in which
the prethermalization of the space-time crystal occurs. This will already
allow us to study the interesting physics associated with the spontaneous
symmetry breaking, as we will see explicitly in Sec. III.1.3, and to make a
direct comparison with experiments in Sec. IV. The second approach builds upon
the first one and uses its techniques to describe the full quantum dynamics in
the low-frequency limit that is needed to show that on the longest time scales
the system ultimately relaxes to the correct thermal equilibrium. These long
time scales have not been explored experimentally and are also expected to be
very difficult to explore in practice. However, we add this quantum approach
here for completeness and to show the consistency of our theory. In
particular, it shows the crucial role played by the fluctuation-dissipation
theorem that is automatically incorporated in the Schwinger-Keldysh formalism
as we have seen.

###  III.1 Semi-classical dynamics

So let us start with the semi-classical approach. We make the approximation
that the dynamics is dominated by a single frequency $\bar{\omega}$ and the
fluctuation-dissipation theorem can be evaluated at twice that typical
frequency. In the case without a drive, that we are still considering at the
moment, we have for instance $\bar{\omega}=-\delta$, where it should be noted
that with $\omega_{D}=0$ we have that $-\delta$ just equals the Bogoliubov
excitation frequency of the resonant axial mode. This means that our
fluctuation-dissipation theorem takes the form

$g^{K}=2ig^{\prime\prime}(1+2N(2\bar{\omega}))\simeq
2ig^{\prime\prime}\frac{k_{B}T}{\hbar\bar{\omega}},$ (14)

where $g^{K}\equiv g^{K}(2\bar{\omega})$, $g=g^{\prime}\pm
ig^{\prime\prime}\equiv g^{(\pm)}(2\bar{\omega})$, and for consistency we also
applied the semi-classical approximation to the Bose distribution. Stated in
the time domain, we thus have
$g^{(\pm,K)}(t-t^{\prime})=g^{(\pm,K)}\delta(t-t^{\prime})$ and all memory
effects have disappeared. Note that this local or Markovian approximation is
indeed only valid at shorter time scales as the Bose distribution diverges for
$\omega\rightarrow 0$. At the longest time scales a full quantum approach is
therefore needed, which we develop in Sec. III.2.

####  III.1.1 Langevin formulation

In the semi-classical approximation a general solution of the Langevin
equation in Eq. (11) reads

$a(t)=e^{i\delta
t}\bigg{\\{}a(0)-ig\int_{0}^{t}dt^{\prime}a^{*}(t^{\prime})a(t^{\prime})a(t^{\prime})e^{-i\delta
t^{\prime}}-\frac{i}{\hbar}\int_{0}^{t}dt^{\prime}\eta(t^{\prime})a^{*}(t^{\prime})e^{-i\delta
t^{\prime}}\bigg{\\}},$ (15)

and similarly for $a^{*}(t)$. With that, we can calculate the equal-time
correlation functions that contain $\eta(t)$, $a(t)$, and $a^{*}(t)$. The
results of interest to us become

$\displaystyle\langle\eta(t)a^{*}(t)\rangle=-\frac{\hbar g^{K}}{2}\langle
a(t)\rangle,$ (16) $\displaystyle\langle\eta(t)a^{*}(t)a^{*}(t)\rangle=-\hbar
g^{K}\langle a^{*}(t)a(t)\rangle,$ (17)

where we used that $g^{K}$ is purely imaginary as explicitly shown in Eq.
(14), and $\langle\eta(t)\rangle=\langle\eta^{*}(t)\rangle=0$. With these
results, we conclude that

$i\hbar\frac{\partial}{\partial t}\langle a(t)\rangle=-\hbar\delta\langle
a(t)\rangle+\hbar g\langle a^{*}(t)a(t)a(t)\rangle-\frac{\hbar
g^{K}}{2}\langle a(t)\rangle,$ (18)

and the complex conjugate equation for $\langle a^{*}(t)\rangle$. In addition
the fluctuations are determined by

$i\hbar\frac{\partial}{\partial t}\langle a^{*}(t)a(t)\rangle=2i\hbar
g^{\prime\prime}\langle a^{*}(t)a^{*}(t)a(t)a(t)\rangle-2\hbar g^{K}\langle
a^{*}(t)a(t)\rangle.$ (19)

This is the full set of equations of motion in the semi-classical approach,
where the term on the right-hand side of Eq. (18) proportional to $g^{K}$ is
the noise-induced drift term from the Stratonovich noise.

For the physical interpretation of Eq. (19) we now show that it corresponds to
the Boltzmann equation. The occupation numbers $N(t)$ of the space-time
crystal are given by [28]

$\langle a^{*}(t)a(t)\rangle=N(t)+\frac{1}{2},$ (20)

containing explicitly both thermal and quantum fluctuations, respectively.
Using Wick’s theorem [26, 27] to work out the four-point correlation function,
we are able to write

$\frac{\partial}{\partial
t}N(t)=-\Gamma\left(N(t)+\frac{1}{2}\right)\bigg{\\{}\left(N(t)+\frac{1}{2}\right)-\frac{k_{B}T}{\hbar\bar{\omega}}\bigg{\\}},$
(21)

where $\Gamma\equiv-4g^{\prime\prime}$ is a rate of decay that determines the
relaxation of the occupation numbers towards the correct semi-classical
equilibrium. It is important to remember that we are still discussing the
situation without the presence of a drive, and the equilibrium we have
obtained here will therefore later on in Sec. III.1.3 be used as an initial
state for the numerical solution of the Fokker-Planck equation that describes
most clearly the spontaneous breaking of the $Z_{2}$ symmetry after the drive
is turned on. Before we start discussing this Fokker-Planck treatment, we
would like to point out that we can write the above rate equation into a more
familiar Boltzmann form as

$\frac{\partial}{\partial t}N(t)=-\Gamma\langle
a^{*}(t)a(t)\rangle\bigg{\\{}N(t)(1+N_{k})-N_{k}(1+N(t))\bigg{\\}},$ (22)

with $N_{k}$ describing the equilibrium occupation numbers of the reservoir,
the momentum ${\bf k}$ being fixed by the energy conservation $\epsilon({\bf
k})=\hbar\bar{\omega}$. Note in particular the overall factor $\langle
a^{*}(t)a(t)\rangle$ on the right-hand side that is due to the multiplicative
noise and effectively renormalizes the rate for quanta of the space-time
crystal to scatter into and out of the heat bath. A Feynman diagram of the
process involved is shown in Fig. 2.

(25,25) thin i o fermioni,v1 fermionv2,o fermion, right, tension=0.05v2,v1
photon, right, tension =0.05v1,v2 v1,v2

Figure 2: Feynman diagram of the scattering process involved in the right-hand
side or ‘collision integral’ of the Boltzmann equation. A straight line
represents $\langle a^{*}(t)a(t)\rangle$ and a wiggly line $\langle
a^{*}_{k}(t)a_{k}(t)\rangle$.

####  III.1.2 Fokker-Planck formulation

We now turn to the Fokker-Planck description of the obtained Langevin
dynamics. We can show that the term $2i\hbar g^{\prime\prime}\langle
a^{*}a^{*}aa\rangle$ on the right-hand side of Eq. (19) is taken into account
by the streaming terms of the Fokker-Planck equation that follow from Eq.
(18). It therefore does not lead to a contribution in the diffusive part of
the Fokker-Planck equation. As a result the diffusion term only has to
reproduce the two terms proportional to $g^{K}$, i.e., the last term in the
right-hand side of Eq. (18) and of Eq. (19) that both arise from the
multiplicative noise. This completely fixes the diffusion term of the Fokker-
Planck equation to be

$\displaystyle-\frac{\hbar g^{K}}{2}\bigg{[}\frac{\partial}{\partial
a^{*}}\bigg{\\{}|a|^{2}\frac{\partial}{\partial a}P[a^{*},a;t]\bigg{\\}}$
$\displaystyle\hskip 28.45274pt+\frac{\partial}{\partial
a}\bigg{\\{}|a|^{2}\frac{\partial}{\partial
a^{*}}P[a^{*},a;t]\bigg{\\}}\bigg{]}.$

Specifically for our case, we had to choose a symmetric diffusion operator.
Together with the streaming terms, which follow directly from Eq. (18), we
find that the Fokker-Planck equation becomes

$\displaystyle i\hbar\frac{\partial}{\partial t}P[a^{*},a;t]$ $\displaystyle=$
$\displaystyle-\frac{\partial}{\partial a}\bigg{\\{}(-\hbar\delta+\hbar
g|a|^{2})aP[a^{*},a;t]\bigg{\\}}$ (23) $\displaystyle+$
$\displaystyle\frac{\partial}{\partial a^{*}}\bigg{\\{}(-\hbar\delta+\hbar
g^{*}|a|^{2})a^{*}P[a^{*},a;t]\bigg{\\}}$ $\displaystyle-$
$\displaystyle\frac{\hbar g^{K}}{2}\frac{\partial}{\partial
a^{*}}\bigg{\\{}|a|^{2}\frac{\partial}{\partial a}P[a^{*},a;t]\bigg{\\}}$
$\displaystyle-$ $\displaystyle\frac{\hbar g^{K}}{2}\frac{\partial}{\partial
a}\bigg{\\{}|a|^{2}\frac{\partial}{\partial a^{*}}P[a^{*},a;t]\bigg{\\}}.$

From the Langevin equation and the resulting equations of motion for the
average moments $\langle a(t)\rangle$ and $\langle a^{*}(t)a(t)\rangle$ we
have found the Fokker-Planck description that is fully equivalent to the
Langevin dynamics. This is the most common procedure in non-equilibrium
statistical physics. However, as already mentioned, we can also derive the
Fokker-Planck equation directly from the effective action, which is a
convenient consistency-check for our theory. This is what we briefly present
now.

Performing the Gaussian integration over the fluctuation field $\xi(t)$, we
obtain from Eq. (4) the desired effective action of the $a(t)$ field alone

$\displaystyle S^{\rm eff}[a^{*},a]\equiv\int dtL(t)$ (24) $\displaystyle=\int
dt\frac{1}{\hbar
g^{K}|a(t)|^{2}}\bigg{\lvert}\bigg{(}i\hbar\frac{\partial}{\partial
t}+\hbar\delta-\hbar g|a(t)|^{2}\bigg{)}a(t)\bigg{\rvert}^{2}.$

Since we still have to perform a functional integral over $a(t)$ to obtain the
partition function, this effective action corresponds to the quantum theory of
a ‘classical’ system with the Lagrangian $L(t)$. The Schrödinger equation of
this quantum theory is exactly the Fokker-Planck equation we are interested
in. To obtain it we thus need to quantize the theory of Eq. (24). In
particular, the canonical momentum $\pi(t)$ of $a(t)$ satisfies

$\pi(t)=\frac{i\hbar}{\hbar
g^{K}|a(t)|^{2}}\bigg{(}-i\hbar\frac{\partial}{\partial t}+\hbar\delta-\hbar
g^{*}|a(t)|^{2}\bigg{)}a^{*}(t),$ (25)

and the complex conjugate for the canonical momentum $\pi^{*}(t)$ of
$a^{*}(t)$. The rest of this derivation goes alone the same way as the
discussion of Ref. [28], i.e., we can now construct the Hamiltonian $H(t)$
from the Lagrangian $L(t)$ to obtain

$\displaystyle H(t)$ $\displaystyle=$
$\displaystyle\frac{\pi(t)}{i\hbar}(-\hbar\delta+\hbar g|a(t)|^{2})a(t)$
$\displaystyle-$ $\displaystyle\frac{\pi^{*}(t)}{i\hbar}(-\hbar\delta+\hbar
g^{*}|a(t)|^{2})a^{*}(t)+\frac{g^{K}|a(t)|^{2}}{\hbar}\lvert\pi(t)\rvert^{2}.$

After requiring the usual commutation relations between coordinates and
canonical momenta, the Schrödinger equation

$i\hbar\frac{\partial}{\partial t}P[a^{*},a;t]=HP[a^{*},a;t],$ (27)

then indeed exactly reproduces the Fokker-Planck equation of Eq. (23). As it
stands the Hamiltonian in Eq. (III.1.2) has clearly operator-ordering problems
but these have fortunately already been resolved from the Langevin derivation
that we carried out first. Interestingly, normal ordering is not appropriate
for the diffusion term in the Fokker-Planck equation, which is a consequence
of the Stratonovich nature of the multiplicative noise that gives rise to a
noise-induced drift term. Normal ordering would lead to the more mathematical
Itô interpretation of the multiplicative noise.

Of course, we can also rederive the Langevin equations from the Fokker-Planck
equation. Taking moments of the Fokker-Planck equation and integrating by
parts nicely reproduces the desired equations of motion in Eqs. (18) and (19).
From this we can conclude that the Langevin equation and the Fokker-Planck
equation are really two sides of the same coin. As a result the equilibrium
probability distribution of the Fokker-Planck equation in Eq. (23) should
agree with the equilibrium occupation numbers obtained from the Boltzmann
equation in Eq. (19). It can therefore be proven that the equilibrium
probability distribution is equal to the semi-classical ideal-gas solution

$P[a^{*},a]=\frac{\hbar\bar{\omega}}{\pi
k_{B}T}\exp\bigg{(}-\frac{\hbar\bar{\omega}}{k_{B}T}|a|^{2}\bigg{)},$ (28)

a fact that is shown in App. A in more detail.

All the above was achieved without the drive terms
$\hbar\omega_{D}A_{D}(\hat{a}^{\dagger}\hat{a}^{\dagger}+\hat{a}\hat{a})/8$ in
our starting Hamiltonian in Eq. (1). Adding the drive to the system is
actually the key ingredient to observe the symmetry breaking we are interested
in, since this is the only term that has explicit phase dependence. More
precisely, the drive breaks the $U(1)$ phase symmetry down to a $Z_{2}$
symmetry $\hat{a}\rightarrow-\hat{a}$ only. Because the time evolution due to
the drive is unitary it is straightforward to add these terms to our model.
Our Langevin equation in Eq. (11) gets just an additional term on the right-
hand side, namely $\hbar\omega_{D}A_{D}a^{*}(t)/4$, and
$\hbar\omega_{D}A_{D}a(t)/4$ for the complex conjugate equation for
$a^{*}(t)$. These will then also have to enter into the streaming terms of the
Fokker-Planck equation, leading to

$\displaystyle i\hbar\frac{\partial}{\partial t}P[a^{*},a;t]=$
$\displaystyle-\frac{\partial}{\partial a}\bigg{\\{}(-\hbar\delta+\hbar
g|a|^{2})aP[a^{*},a;t]+\frac{\hbar\omega_{D}A_{D}}{4}a^{*}P[a^{*},a;t]\bigg{\\}}$
(29) $\displaystyle+\frac{\partial}{\partial
a^{*}}\bigg{\\{}(-\hbar\delta+\hbar
g^{*}|a|^{2})a^{*}P[a^{*},a;t]+\frac{\hbar\omega_{D}A_{D}}{4}aP[a^{*},a;t]\bigg{\\}}$
$\displaystyle-\frac{\hbar g^{K}}{2}\frac{\partial}{\partial
a^{*}}\bigg{\\{}|a|^{2}\frac{\partial}{\partial
a}P[a^{*},a;t]\bigg{\\}}-\frac{\hbar g^{K}}{2}\frac{\partial}{\partial
a}\bigg{\\{}|a|^{2}\frac{\partial}{\partial a^{*}}P[a^{*},a;t]\bigg{\\}}.$

This is now the complete Fokker-Planck equation, containing all the non-
equilibrium dynamics of the space-time crystal in the semiclassical
approximation. As this equation accurately describes the prethermalization of
the space-time crystal observed experimentally, it is the most important
result of our paper.

(a) After 0 oscillations.

(b) After 10 oscillations.

(c) After 20 oscillations.

(d) After 30 oscillations.

(e) After 40 oscillations.

(f) After 50 oscillations.

Figure 3: Dynamical behavior of the probability distribution $P[a^{*},a;t]$ in
the complex $(a^{\prime},a^{\prime\prime})$ plane. Numerically, fifty
oscillations, i.e., fifty periods of the drive, are simulated, and a picture
was taken after every ten oscillations. The number of oscillations past is
indicated below each figure. For clarity, we have suppressed the center part
of the distribution after 20 oscillations, since otherwise the range of the
color scale is suppressed by one order of magnitude.

####  III.1.3 Numerical solution

We have numerically solved the Fokker-Planck equation in Eq. (29) by using an
implicit finite-difference method, see App. B for details. We have used the
same parameters as obtained experimentally in Ref. [24], namely $\delta/2\pi$
= 2 Hz, $|g|/2\pi$ = 3$\times 10^{-4}$ Hz, $\phi_{g}$ = -0.4$\pi$, and $A_{D}$
= 0.091, except that we have increased $|g|$ by a factor of 10 in order to
decrease $|\langle a\rangle|^{2}$ by the same factor for numerical
convenience. In particular, we simulate what happens when we start with the
thermal Gaussian initial distribution in Eq. (28) and then let the system
evolve in time under influence of the drive. Note that our initial condition
is $U(1)$ invariant, i.e., rotationally invariant in the complex
$(a^{\prime},a^{\prime\prime})$ plane, and therefore does not contain a
preferred phase. It is also invariant under inversion with respect to the
origin $|a|=0$, which is the $Z_{2}$ symmetry we are most interested in. So
the initial condition does not break any symmetries as desired for a
discussion of spontaneous symmetry breaking. Our numerical findings of the
spontaneous symmetry breaking in our time crystal is shown in Fig. 3. We
clearly see that initially the drive is only squeezing the probability
distribution, which breaks the initial $U(1)$ symmetry down to a $Z_{2}$
symmetry, as only the inversion symmetry in the complex plane remains. This
symmetry breaking is not spontaneous as the drive is explicitly breaking the
$U(1)/Z_{2}$ symmetry. In this stage of the evolution the probability
distribution is squeezed, but still has a maximum in the origin of the complex
plane. In terms of order parameters this phase is thus characterized by
$\langle aa\rangle\neq 0$ and $\langle a\rangle=0$, so the $Z_{2}$ symmetry is
unbroken. Only at a later stage this single maximum in the origin falls apart
into two maxima away from the origin that are related by the $Z_{2}$ inversion
symmetry.

Figure 4: Time evolution of the reduced probability distribution on the
$\phi=3\pi/4$ (mod $\pi$) diagonal, obtained by performing the integral of the
probability distribution $P[a^{*},a;t]$ over all lines perpendicular to this
diagonal.

To emphasize this effect, we have integrated the probability distribution over
all lines perpendicular to the $\phi=3\pi/4$ (mod $\pi$) diagonal and show the
results for this reduced probability distribution in Fig. 4. Initially, the
reduced distribution is Gaussian with its center at $|a|$ = 0. After the drive
has been turned on, the central part of the distribution spreads out to the
left and right and after 20 oscillations two distinct peaks appear. However,
note that even for a large number of oscillations there remains a small part
of the distribution near $a$ = 0. This comes about because the strength of the
driving $\hbar\omega_{D}A_{D}/4$ is multiplied with $a$ and $a^{*}$ in the
Fokker-Planck equation and so for $|a|\simeq 0$ the effect of the drive
becomes negligible. Moreover the noise is multiplicative and vanishes in this
limit.

The appearance of the two symmetry-related peaks corresponds to the true
spontaneous breaking of the $Z_{2}$ symmetry and the formation of the space-
time crystal as now both $\langle aa\rangle\neq 0$ and $\langle a\rangle\neq
0$. Note that in principle the probability distribution is inversion symmetric
at all times, as the Hamiltonian does not break this symmetric explicitly.
This implies that if one performs many experiments with the same thermal
initial condition for the Bose-Einstein condensate, the observed space-time
crystal will have a temporal phase of $\phi$ or $\phi+\pi$ with exactly equal
probability [25].

###  III.2 Quantum dynamics

In classical non-equilibrium physics we usually deal with white noise, and a
frequency-independent fluctuation-dissipation theorem. Perhaps the most famous
example of white noise is Brownian motion with the Einstein relation between
the diffusion constant and temperature. In the above semi-classical approach
for the experimentally relevant (prethermalization) time scales, the same
approximation was made by evaluating the exact fluctuation-dissipation theorem
in Eq. (7) at twice the typical frequency $\bar{\omega}$. Quantum noise,
however, is in general colored noise and this complicates the analytical or
numerical analysis of the relevant Langevin problem considerable. Fortunately,
if we are interested in the longest time scales of the problem in which
ultimately relaxation towards thermal equilibrium occurs, we are allowed to
make another approximation, i.e., we can take the long-wavelength or in this
case the low-frequency limit. In that limit $g^{K}(\omega)$ goes to the
nonzero constant $g^{K}(0)$ and we have

$g^{K}(t-t^{\prime})\simeq g^{K}(0)\delta(t-t^{\prime}),$ (30)

which again corresponds to white noise in frequency space. We would like to
take the same limit in all quantities that appear in the fluctuation-
dissipation theorem, but this is not possible because we have

$(1+2N(\hbar\omega))\simeq\frac{2k_{B}T}{\hbar\omega},$ (31)

which in fact diverges for $\omega\rightarrow 0$. In order to compensate for
this divergence in the exact fluctuation-dissipation theorem in Eq. (7), we
therefore obtain that in frequency space we must have that
$g^{\prime\prime}(\omega)=(\hbar\omega/4ik_{B}T)g^{K}(0)$ or in the time
domain that

$g^{\prime\prime}(t-t^{\prime})=\frac{g^{K}(0)}{4k_{B}T}\hbar\frac{\partial}{\partial
t}\delta(t-t^{\prime}).$ (32)

####  III.2.1 Langevin formulation

With that, our Langevin equation in Eq. (11) becomes

$\displaystyle\bigg{(}1$ $\displaystyle-$ $\displaystyle\frac{\hbar
g^{K}}{2k_{B}T}a^{*}(t)a(t)\bigg{)}i\hbar\frac{\partial}{\partial t}a(t)$
$\displaystyle=$ $\displaystyle-\hbar\delta a(t)+\hbar
g^{\prime}a^{*}(t)a(t)a(t)+\eta(t)a^{*}(t),$

and the complex conjugate equation for $a^{*}(t)$. Notice that for
$g^{\prime}$ we are allowed to take the naive zero-frequency limit and that
from now on we use $g^{K}\equiv g^{K}(\omega=0)$. Moreover, for the same
reasons as before we have briefly returned to the case without drive, as the
unitary evolution of the drive is easily incorporated later on as we have
seen. We also introduce the (stochastic) quantities

$I(t)=1-\frac{\hbar g^{K}}{2k_{B}T}|a(t)|^{2}$ (34)

and thus

$I^{*}(t)=1+\frac{\hbar g^{K}}{2k_{B}T}|a(t)|^{2},$ (35)

as $g^{K}$ is purely imaginary. These will play a key role in our following
discussion.

We now want to work out again the relevant correlation functions for the noise
$\eta(t)$ with certain functions of $a(t)$ and $a^{*}(t)$, in order to find
the equations of motion for $\langle a(t)\rangle$ and $\langle
a^{*}(t)a(t)\rangle$, similar to Eqs. (18) and (19). By dividing the Langevin
equation in Eq. (III.2.1) by $I(t)$, we can see that our task now is to
calculate quantities like $\langle\eta(t)a^{*}(t)/I(t)\rangle$ and
$\langle\eta(t)a^{*}(t)a^{*}(t)/I(t)\rangle$. Let us therefore also define

$A^{*}(t)\equiv\frac{a^{*}(t)}{I(t)}.$ (36)

We can expand this function as

$A^{*}(t)=A^{*}(0)+\frac{\partial A^{*}(0)}{\partial
a^{*}}(a^{*}(t)-a^{*}(0))+\frac{\partial A^{*}(0)}{\partial a}(a(t)-a(0)),$
(37)

neglecting all higher-order contribution. So we always need to remember that
these stochastic quantities we introduced are explicit functions of $a^{*}(t)$
and $a(t)$.

Using again that
$\langle\eta(t)\rangle=\langle\eta^{*}(t)\rangle=\langle\eta(t)\eta(t^{\prime})\rangle=\langle\eta^{*}(t)\eta^{*}(t^{\prime})\rangle=0$,
and the general solution for $a(t)$ in this case given by

$a(t)=a(0)-\frac{i}{\hbar}\int_{0}^{t}dt^{\prime}\bigg{\\{}-\frac{\hbar\delta}{I(t^{\prime})}a(t^{\prime})+\hbar
g^{\prime}A^{*}(t^{\prime})a(t^{\prime})a(t^{\prime})+\eta(t^{\prime})A^{*}(t^{\prime})\bigg{\\}},$
(38)

we are in the position to calculate the desired correlation functions. The
results become

$\langle\eta(t)A^{*}(t)\rangle=-\frac{\hbar
g^{K}}{2}\left\langle\frac{\partial A^{*}(t)}{\partial
a^{*}}A(t)\right\rangle,$ (39)

and

$\displaystyle\langle\eta(t)a^{*}(t)A^{*}(t)\rangle$ $\displaystyle=$
$\displaystyle-\frac{\hbar g^{K}}{2}\left\langle A(t)A^{*}(t)\right\rangle$
(40) $\displaystyle-$ $\displaystyle\frac{\hbar g^{K}}{2}\left\langle
a^{*}(t)\frac{\partial A^{*}(t)}{\partial a^{*}}A(t)\right\rangle.$

With that, the full set of equations of motion are now obtained as

$i\hbar\frac{\partial}{\partial t}\langle
a(t)\rangle=\bigg{\langle}(-\hbar\delta+\hbar
g^{\prime}a^{*}(t)a(t))\frac{a(t)}{I(t)}\bigg{\rangle}-\frac{\hbar
g^{K}}{2}\bigg{\langle}\frac{\partial A^{*}(t)}{\partial
a^{*}}A(t)\bigg{\rangle},$ (41)

the complex conjugate equation for $\langle a^{*}(t)\rangle$, and

$\displaystyle i\hbar\frac{\partial}{\partial t}\langle a^{*}(t)a(t)\rangle$
(42) $\displaystyle=\bigg{\langle}(-\hbar\delta+\hbar
g^{\prime}a^{*}(t)a(t))\bigg{(}\frac{a(t)}{I(t)}-\frac{a^{*}(t)}{I^{*}(t)}\bigg{)}\bigg{\rangle}$
$\displaystyle-\frac{\hbar g^{K}}{2}\left\langle a^{*}(t)\frac{\partial
A^{*}(t)}{\partial a^{*}}A(t)+a(t)\frac{\partial A(t)}{\partial
a}A^{*}(t)\right\rangle$ $\displaystyle-\hbar g^{K}\langle
A^{*}(t)A(t)\rangle,$

which may be directly compared to Eqs. (18) and (19).

####  III.2.2 Fokker-Planck formulation

Following the same line of thought as in the semi-classical case, we obtain
for the case without drive the Fokker-Planck equation

$\displaystyle i\hbar\frac{\partial}{\partial t}P[a^{*},a;t]$ $\displaystyle=$
$\displaystyle-\frac{\partial}{\partial a}\bigg{\\{}(-\hbar\delta+\hbar
g^{\prime}|a|^{2})\frac{a}{I}P[a^{*},a;t]\bigg{\\}}$ (43) $\displaystyle+$
$\displaystyle\frac{\partial}{\partial a^{*}}\bigg{\\{}(-\hbar\delta+\hbar
g^{\prime}|a|^{2})\frac{a^{*}}{I^{*}}P[a^{*},a;t]\bigg{\\}}$ $\displaystyle-$
$\displaystyle\frac{\hbar g^{K}}{2}\frac{\partial}{\partial
a^{*}}\bigg{\\{}\frac{a}{I^{*}}\frac{\partial}{\partial
a}\frac{a^{*}}{I}P[a^{*},a;t]\bigg{\\}}$ $\displaystyle-$
$\displaystyle\frac{\hbar g^{K}}{2}\frac{\partial}{\partial
a}\bigg{\\{}\frac{a^{*}}{I}\frac{\partial}{\partial
a^{*}}\frac{a}{I^{*}}P[a^{*},a,t]\bigg{\\}},$

where we want to emphasize the fact that $g$ and $g^{*}$ in Eq. (23) have now
both become equal to $g^{\prime}$. As a result the dissipation in the
equations of motion effectively does no longer enter via the nonlinear
coupling constant, but through the imaginary parts of the quantities $I$ and
$I^{*}$ instead.

This Fokker-Planck equation can be derived directly from the effective action
in Eq. (4) in the same manner as described before. Instead of a constant
$g^{\prime\prime}$ in Eq. (24), we now have a contribution containing a time
derivative. For this situation, the canonical momentum $\pi(t)$ therefore
reads

$\pi(t)=\frac{i\hbar I(t)}{\hbar
g^{K}|a(t)|^{2}}\bigg{(}-i\hbar\frac{\partial}{\partial t}+\hbar\delta-\hbar
g^{\prime}|a(t)|^{2}\bigg{)}a^{*}(t),$ (44)

and the Hamiltonian becomes after performing the usual Legendre transformation

$\displaystyle H(t)$ $\displaystyle=$ $\displaystyle\frac{\pi(t)}{i\hbar
I(t)}(-\hbar\delta+\hbar g^{\prime}|a(t)|^{2})a(t)$ $\displaystyle-$
$\displaystyle\frac{\pi^{*}(t)}{i\hbar I^{*}(t)}(-\hbar\delta+\hbar
g^{\prime}|a(t)|^{2})a^{*}(t)-\frac{g^{K}|a(t)|^{2}}{\hbar|I(t)|^{2}}\lvert\pi(t)\rvert^{2},$

which after quantization, and resolution of the usual operator-ordering
problems, indeed exactly reproduces the Fokker-Planck equation in Eq. (43).

If we also include the drive terms again, similar to the approach in the semi-
classical case, the full Fokker-Planck equation reads now

$\displaystyle i\hbar\frac{\partial}{\partial t}P[a^{*},a;t]=$
$\displaystyle-\frac{\partial}{\partial a}\bigg{\\{}(-\hbar\delta+\hbar
g^{\prime}|a|^{2})\frac{a}{I}P[a^{*},a;t]+\frac{\hbar\omega_{D}A_{D}}{4}\frac{a^{*}}{I}P[a^{*},a;t]\bigg{\\}}$
(46) $\displaystyle+\frac{\partial}{\partial
a^{*}}\bigg{\\{}(-\hbar\delta+\hbar
g^{\prime}|a|^{2})\frac{a^{*}}{I^{*}}P[a^{*},a;t]+\frac{\hbar\omega_{D}A_{D}}{4}\frac{a}{I^{*}}P[a^{*},a;t]\bigg{\\}}$
$\displaystyle-\frac{\hbar g^{K}}{2}\frac{\partial}{\partial
a^{*}}\bigg{\\{}\frac{a}{I^{*}}\frac{\partial}{\partial
a}\frac{a^{*}}{I}P[a^{*},a;t]\bigg{\\}}-\frac{\hbar
g^{K}}{2}\frac{\partial}{\partial
a}\bigg{\\{}\frac{a^{*}}{I}\frac{\partial}{\partial
a^{*}}\frac{a}{I^{*}}P[a^{*},a;t]\bigg{\\}}.$

Although this appears to be more complicated then in the semi-classical case,
the equilibrium distribution can now be obtained analytically and ultimately
reads

$P[a^{*},a]\propto|I(|a|^{2})|\exp\bigg{(}-\frac{1}{k_{B}T}\bigg{\\{}-\hbar\delta|a|^{2}+\frac{\hbar
g^{\prime}}{2}|a|^{4}+\frac{\hbar\omega_{D}A_{D}}{8}(a^{*}a^{*}+aa)\bigg{\\}}\bigg{)}.$
(47)

Notice that the equilibrium distribution scales with the Boltzmann factor
$e^{-E/k_{B}T}$, where $E[a^{*},a]$ is the classical energy associated with
the quantum Hamiltonian in Eq. (1), which was the starting point of our
developments. Apart from this Boltzmann factor, we find also a prefactor equal
to $|I(|a|^{2})|$. This precisely corresponds to the expected change in phase-
space volume [32, 33] as our Hamiltonian dynamics is modified into

$i\hbar I\frac{\partial a}{\partial t}=\frac{\partial E[a^{*},a]}{\partial
a^{*}},~{}-i\hbar I^{*}\frac{\partial a^{*}}{\partial t}=\frac{\partial
E[a^{*},a]}{\partial a}.$ (48)

In the overdamped limit $\hbar|g^{K}|\gg k_{B}T$ this prefactor becomes
$\hbar|g^{K}||a|^{2}/2k_{B}T$ and therefore the stationary solution will
always be zero at $|a|=0$. On the other hand, in the underdamped limit, we
have $I\simeq 1$ and the prefactor doesn’t depend on $a$ or $a^{*}$ and thus
it can be absorbed in the normalisation constant. For the experimental
parameters of Ref. [24], the underdamped limit is indeed a very good
approximation, which confirms that dissipation is only playing a minor role in
this case as expected for a superfluid droplet that is almost completely Bose-
Einstein condensed.

##  IV Comparison with experiments

In this section, our goal is to show the result for the stationary solution of
the Fokker-Planck equation for the quantum dynamics in Eq. (47) and compare it
with the semi-classical dynamics, either using the Langevin equation or the
Fokker-Planck equation. The semi-classical dynamics within the Langevin
formulation, after including also technical noise, is found to agree very well
with experiments in Ref. [25]. In Fig. 5 we therefore compare the result for
the equilibrium distribution from the Langevin equation with the result from
the Fokker-Planck equation, both for the semi-classical approach, and obtain,
as expected, excellent agreement. The results for the Langevin equation are
obtained by calculating 100,000 trajectories using random values of $a$ and
$a^{*}$ from the initial distribution of Eq. (28) and integrating the
evolution of the individual trajectories over 100 oscillations, as discussed
in Ref. [25]. The results are then binned to obtain the probability
distribution $P[a^{*},a]$. The results for the Fokker-Planck distribution are
obtained by taking an initial Gaussian distribution and letting it evolve to
an equilibrium using Eq. (29). Not only do the two equations lead to the same
equilibrium, also the dynamics towards the equilibrium is very similar in both
cases and this numerically validates the equivalence between the two
theoretical frameworks.

(a) Equilibrium distribution for the semi-classical approach using the
Langevin equation.

(b) Equilibrium distribution for the semi-classical approach using the Fokker-
Planck equation.

(c) Equilibrium distribution for the quantum approach using the Langevin
equation.

(d) Equilibrium distribution for the quantum approach using the Fokker-Planck
equation.

Figure 5: Comparison between the two numerical simulations of the equilibrium
probability distributions $P[a^{*},a]$ for both the semi-classical and quantum
cases, respectively. All the results are calculated for the experimental
parameters of Ref. [24].

To compare the semi-classical dynamics to the quantum dynamics, the Langevin
equation of Eq. (III.2.1) is used to obtain the stochastic trajectories in the
quantum approach. Since in this approach the imaginary part $g^{\prime\prime}$
of the interaction parameter has effectively been removed from the equations
of motion in favor of $g^{K}$, we use $g^{\prime}=|g|$ with $|g|$ the
interaction parameter as previously used in the semi-classical approach. In
detail the probability distributions are very different as implicitly also
found in Ref. [24]. One important difference is immediately clear from the
location of the peaks. In the semi-classical simulations, the peaks are
located almost on the diagonal with $\phi\simeq 3\pi/4~{}({\rm mod}~{}\pi)$,
whereas in the ultimate equilibrium state the peaks are located on the
imaginary axis. In the experiment we find the peak to locate nearly on the
real axis, but in the experiment we have determined the amplitude of the
highly-excited axial mode, and this amplitude is proportional to the time
derivative of $a$ in the laboratory frame and this leads to a phase shift of
$\pi/4$. Taking this subtle point into account, we obtain excellent agreement
between the experiment and the simulations of the semi-classical Langevin
equation as shown in Ref. [24]. The location of the maxima in the semi-
classical equilibrium also depends sensitively on the detuning used, whereas
in the true equilibrium the peaks are always on the imaginary axis independent
of detuning. For instance, in our numerical solution of the Fokker-Planck
equation in Sec. III.1.3 we used $g^{\prime}\ll g^{\prime\prime}$. In that
case we obtain from the equations of motion from $\langle a(t)\rangle$ that in
the prethermal state

$\cos(2\phi)=\frac{4\delta}{\omega_{D}A_{D}}~{},$ (49)

with $\sin(2\phi)<0$, and

$|\langle
a\rangle|^{2}=\frac{k_{B}T}{\hbar\bar{\omega}}+\frac{1}{|g|}\sqrt{\bigg{(}\frac{\omega_{D}A_{D}}{4}\bigg{)}^{2}-\delta^{2}}~{}.$
(50)

The latter agrees with Ref. [24] except for the first term in the right-hand
side that represents the contribution from the noise-induced drift and that is
small compared to the second term under the conditions of interest. In that
regime we can actually take the limit $\omega_{D}A_{D}\gg{\rm
max}(4|\delta|,|g|k_{B}T/\hbar\bar{\omega})$ to show, next to $|\langle
a\rangle|\simeq\sqrt{\omega_{D}A_{D}/4|g|}$, that the phase of the space-time
crystal becomes in general

$\phi\simeq\frac{\pi-\phi_{g}}{2}~{}({\rm mod}~{}\pi),$ (51)

where $-\pi<\phi_{g}<0$ is the phase of $g$ [25].

We thus confirm here that our semi-classical approach can very accurately
describe the experiments and that the experiments have explored up to now only
the relatively short time scales in which prethermalization occurs, but not
yet the longer time scales on which the final thermalization would take place.
The latter may be very hard to realize in practice for two reasons. First, the
dissipation in the experiment also affects the radial breathing mode that acts
as a drive for the axial collective modes. As a result the amplitude of the
drive is not constant and slowly decays at those long time scales. Second, our
analysis is based on the Hamiltonian in Eq. (1), which is valid in the
rotating frame of the drive after neglecting the non-resonant terms. On the
long time scales of interest for the thermalization, these non-resonant, and
therefore oscillating terms will lead to heating that requires a further
analysis outside the scope of our paper.

For completeness, we reiterate that when we look at the width of the peaks, we
notice that our theoretical results are always narrower than the experimental
data. This can be explained by the fact that in the experiments additional
broadening occurs from the fluctuations in the initial conditions, most
importantly the number of atoms in the Bose-Einstein condensate. As is shown
in Ref. [25] this can be most easily accounted for in our theoretical modeling
by allowing for some fluctuations in the detuning $\delta$.

##  V Conclusion and outlook

The main achievements of our work can be summarized as follows. We developed
the general theory to describe the non-equilibrium physics of the space-time
crystal in an atomic Bose-Einstein condensate. Starting from a non-linear
(cubic) coupling of the space-time crystal to a heat bath, we applied the
tools of the Schwinger-Keldysh formalism to derive the general form of the
Langevin equation. At that point, we distinguished two cases: a semi-classical
approach, where the frequency dependence of the dissipation was neglected and
a fully quantum approach, where this frequency dependency was included in the
long wave-length limit. For both cases, we presented the full system of
equations describing the dynamics: both in the Langevin formulation as well as
in the Fokker-Planck formulation. We also showed the consistency of the two
formulations. Furthermore, we were able to derive stationary solutions for
both of these Fokker-Planck equations, which, as expected, turned out to scale
with the Boltzmann factor $e^{-E/k_{B}T}$, where $E[a^{*},a]$ is the
appropriate energy for that specific situation.

By solving the semi-classical Fokker-Planck equation numerically, we were able
to visualize very explicitly the dynamical features of the spontaneous
breaking of the $Z_{2}$ symmetry and the associated formation of the space-
time crystal. Since the semi-classical dynamics describes the
prethermalization at relative short time scales, it can be directly compared
with the experimental data. Excellent agreement between theory and experiment
is obtained in this manner. The quantum method on the other hand makes it
possible to discuss the ultimate thermalization of our theoretical model, but
such long time scales are presently not observable experimentally and may be a
topic of further investigations.

An open end in our discussion is the origin of the cubic coupling term between
the space-time crystal and the heat bath. Although its theoretical
implications work out nicely as these are mostly based on universal features
such as the fluctuation-dissipation theorem, its precise (experimental) origin
remains unclear. In particular, we would like to understand better which modes
of the thermal cloud effectively represent the heat bath in our approach. In
this manner it might be possible to derive also the strength of the
dissipation and the noise from first principles. Nonetheless, our results
shows that we have already obtained a solid theoretical description of the
experiments, with minimal approximations, that underpins the spontaneous
symmetry breaking occurring in these kind of time crystals.

An interesting area for further research will be to use our Fokker-Planck
equation in order to study the phenomenon of ‘tunneling’ from one energy
minimum to the other. Although we did not look into this yet in our
discussion, the non-equilibrium physics of this problem is also fully captured
by our model and may also be in reach in future experiments with space-time
crystals in Bose-Einstein condensates.

##  Appendix A Semi-classical equilibrium without drive

The equilibrium in the semi-classical approach has to satisfy the stationary
Fokker-Planck equation for that case, which reads

$\displaystyle 0$ $\displaystyle=$ $\displaystyle-\frac{\partial}{\partial
a}\bigg{\\{}(-\hbar\delta+\hbar g|a|^{2})aP[a^{*},a]\bigg{\\}}$ (52)
$\displaystyle+$ $\displaystyle\frac{\partial}{\partial
a^{*}}\bigg{\\{}(-\hbar\delta+\hbar g^{*}|a|^{2})a^{*}P[a^{*},a]\bigg{\\}}$
$\displaystyle-$ $\displaystyle\frac{\hbar g^{K}}{2}\frac{\partial}{\partial
a^{*}}\bigg{\\{}|a|^{2}\frac{\partial}{\partial a}P[a^{*},a]\bigg{\\}}$
$\displaystyle-$ $\displaystyle\frac{\hbar g^{K}}{2}\frac{\partial}{\partial
a}\bigg{\\{}|a|^{2}\frac{\partial}{\partial a^{*}}P[a^{*},a]\bigg{\\}}.$

Furthermore, the fluctuation-dissipation theorem in the semi-classical
approximation tells us that we have the relation
$g^{K}=2ig^{\prime\prime}(k_{B}T/\hbar\bar{\omega})$. Using then a phase
independent probability distribution $P[a^{*},a]=P[|a|^{2}]$, we can easily
see that both the terms with $-\hbar\delta$ as well as the terms with $\hbar
g^{\prime}$ (the real part of $\hbar g$) drop out. What is left is

$\displaystyle 0$ $\displaystyle=$
$\displaystyle\bigg{(}\frac{k_{B}T}{\hbar\bar{\omega}}\bigg{)}\bigg{[}\frac{\partial}{\partial
a^{*}}\bigg{\\{}|a|^{2}\frac{\partial}{\partial
a}P\bigg{\\}}+\frac{\partial}{\partial
a}\bigg{\\{}|a|^{2}\frac{\partial}{\partial a^{*}}P\bigg{\\}}\bigg{]}$ (53)
$\displaystyle+$ $\displaystyle 4|a|^{2}P+2|a|^{4}\frac{\partial
P}{\partial|a|^{2}}.$

Working this out further, we notice that

$\frac{\partial}{\partial a^{*}}\bigg{\\{}|a|^{2}\frac{\partial}{\partial
a}P\bigg{\\}}+\frac{\partial}{\partial
a}\bigg{\\{}|a|^{2}\frac{\partial}{\partial a^{*}}P\bigg{\\}}\\\
=4|a|^{2}\frac{\partial
P}{\partial|a|^{2}}+2|a|^{4}\frac{\partial^{2}P}{(\partial|a|^{2})^{2}},$ (54)

which means that the stationary Fokker-Planck equation becomes

$0=2|a|^{4}\bigg{\\{}\bigg{(}\frac{k_{B}T}{\hbar\bar{\omega}}\bigg{)}\frac{\partial^{2}P}{(\partial|a|^{2})^{2}}+\frac{\partial
P}{\partial|a|^{2}}\bigg{\\}}+4|a|^{2}\bigg{\\{}\bigg{(}\frac{k_{B}T}{\hbar\bar{\omega}}\bigg{)}\frac{\partial
P}{\partial|a|^{2}}+P\bigg{\\}}.$ (55)

We therefore conclude that the probability distribution obeys simply

$\frac{\partial P}{\partial|a|^{2}}=-\frac{\hbar\bar{\omega}}{k_{B}T}P,$ (56)

which is the defining equation for the Gaussian ideal-gas solution

$P[a^{*},a]\propto\exp\bigg{(}-\frac{\hbar\bar{\omega}}{k_{B}T}|a|^{2}\bigg{)},$
(57)

proportional to the expected classical Boltzmann factor. Note that there is no
change in the phase-space volume in this case and there is no additional
prefactor.

##  Appendix B Numerical simulation of semi-classical dynamics

The Fokker-Planck equation of Eq. (29) is given in terms of the eigenvalue $a$
of the annihilation operator and its complex conjugate $a^{*}$. The equation
looks deceptively simple with the convection terms proportional to $\delta$,
$g$, and $\omega_{D}$, and the diffusive term proportional to $g^{K}$. The
probability distribution $P[a^{*},a;t]$ can be represented on a grid and the
spatial derivatives are evaluated using finite differences. The time-evolution
of the probability distribution is then integrated forward in time using
simple first-order finite difference in time. However, since $P$ is a real
function, only terms on the right-hand side of Eq. (29) that are imaginary
contribute, whereas the real terms cancel. Furthermore, some of the terms are
proportional to $|a|^{2}$ and since $a$ becomes large, these terms become
large at the edges of the grid. At the same time the probability distribution
can become narrow in some directions and thus the step size $\Delta a$ has to
be chosen rather small. This makes the number of points in the grid large and
precludes the use of implicit finite-difference methods. In the case of
explicit finite-difference methods, the step size $\Delta t$ in time is bound
by $\Delta a$ through the Courant condition, which requires small step sizes
in time.

To avoid the cancellation of terms in the evolution of $P$, we use real
coordinates $x=(a+a^{*})/2$ and $y=(a-a^{*})/2i$ and rewrite the Fokker-Planck
equation in terms of these coordinates using

$\frac{\partial}{\partial a}=\frac{1}{2}\left(\frac{\partial}{\partial
x}-i\frac{\partial}{\partial y}\right),\qquad\frac{\partial}{\partial
a^{*}}=\frac{1}{2}\left(\frac{\partial}{\partial x}+i\frac{\partial}{\partial
y}\right).$ (58)

For the numerical integration, it is beneficial to write the evolution
equation in the form of a conservation law and after some algabra we find

$\frac{\partial}{\partial t}P=-\frac{\partial}{\partial
x}j_{x}^{P}-\frac{\partial}{\partial y}j_{y}^{P},$ (59)

with the probability current densities

$j_{x}^{P}=-\delta
yP+\left(g^{\prime}y+g^{\prime\prime}x\right)\left(x^{2}+y^{2}\right)P-\frac{\omega_{D}A_{D}}{4}yP+\frac{g^{K}}{4i}\left(x^{2}+y^{2}\right)\frac{\partial}{\partial
x}P,$ (60)

and

$j_{y}^{P}=\delta
xP+\left(-g^{\prime}x+g^{\prime\prime}y\right)\left(x^{2}+y^{2}\right)P-\frac{\omega_{D}A_{D}}{4}xP+\frac{g^{K}}{4i}\left(x^{2}+y^{2}\right)\frac{\partial}{\partial
y}P,$ (61)

where we remind ourselves that $g^{K}$ is purely imaginary. These equations
are used for the simulations shown in Fig. 3.

There are many finite-difference methods proposed for this convection-
diffusion problem. After some trial and error, we have selected the MacCormack
method [34], where the evolution is obtained after first making a predictor
step, displacing the grid over half a grid spacing $\Delta a$ to the right,
followed by a corrector step displacing the grid back by shifting it half a
grid spacing to the left. This makes the method conditionally stable, compared
to the simple forward time, centered space (FTCS) method. The stability is
determined by the Courant-Friedrich-Lewy criterion, or Courant condition,
given in our case by

$\frac{|j_{x}|\Delta t}{P\Delta a}\leq 1,\qquad\frac{|j_{y}|\Delta t}{P\Delta
a}\leq 1.$ (62)

Since the width of $P$ can become small in the evolution, we have to choose
$\Delta a$ small and thus require a large number of time-steps for the
evolution.

The Fokker-Planck can also be rewritten in polar coordinates $\rho$ and $\phi$
and again after some algebra one arrives at

$\frac{\partial}{\partial
t}P=-\frac{1}{\rho}\frac{\partial}{\partial\rho}\left(\rho
j^{P}_{\rho}\right)-\frac{1}{\rho}\frac{\partial}{\partial\phi}j_{\phi}^{P},$
(63)

with

$j_{\rho}^{P}=-\frac{\omega_{D}A_{D}}{4}\rho\sin 2\phi
P+g^{\prime\prime}\rho^{3}P+\frac{g^{K}}{4i}\rho^{2}\frac{\partial}{\partial\rho}P,$
(64)

and

$j_{\phi}^{P}=\delta\rho P-\frac{\omega_{D}A_{D}}{4}\rho\cos 2\phi
P-g^{\prime}\rho^{3}P+\frac{g^{K}}{4i}\rho\frac{\partial}{\partial\phi}P,$
(65)

where we have expressed the right-hand side of Eq. (63) in the form of a
conservation law in polar coordinates. Note that the division by $\rho$ for
each of the terms can be avoided by performing the integration in time using
$\rho P(\rho,\phi)$ instead of $P(\rho,\phi)$.

These equations allow us to find the location of the maximum of the
equilibrium distribution using $j^{P}_{\rho}=0$ and $j^{P}_{\phi}=0$. The
maximum is located at

$\langle\rho\rangle^{2}=\frac{1}{|g|}\left(\delta\cos\phi_{g}+\sqrt{\left(\frac{\omega_{D}A_{D}}{4}\right)^{2}-\delta^{2}\sin^{2}\phi_{g}}\right),$
(66)

and

$\tan
2\phi=\frac{g^{\prime\prime}\langle\rho\rangle^{2}}{\delta-g^{\prime}\langle\rho\rangle^{2}}.$
(67)

The first result has already been derived by Smits et al. [24]. In the case of
small $g^{\prime}$, these equations correspond to Eqs. (49) and (50) without
the noise-induced drift that gives the small difference between the location
of the maximum of the probability distribution and the average
$\langle\rho\rangle$.

## References

  * [1] L. Landau and E. Lifshitz, Statistical Physics, (Elsevier Butterworth-Heinemann, London, 1984)
  * [2] K. Huang, Statistical Mechanics. (Wiley, Hoboken, 1987).
  * [3] M. C. Cross and P. C. Hohenberg, Rev. Mod. Phys. 65, 851 (1993).
  * [4] T. E. Angelini, E. Hannezo, X. Trepat, M. Marquez, J. J. Fredberge, and D. A. Weitz, Proc. Natl. Acad. Sci. USA 108, 4714 (2011).
  * [5] I. Theurkauff, C. Cottin-Bizonne, J. Palacci, C. Ybert, and L. Bocquet, Phys. Rev. Lett. 108, 268303 (2012).
  * [6] M. Greiner, O. Mandel, T. Hänsch, and I. Bloch, Nature 419, 51 (2002)
  * [7] T. Kinoshita, T. Wenger, and D. Weiss, Nature 440, 900 (2006).
  * [8] C. Weitenberg, M. Endres, J.F. Sherson, M. Cheneau, P. Schausz, T. Fukuhara, I. Bloch, and S. Kuhr, Nature 471, 319 (2011).
  * [9] J. Eisert, M. Friesdorf, and C. Gogolin, Nat. Phys. 11, 124 (2015).
  * [10] T. Langen, R. Geiger, and J. Schmiedmayer, Annu. Rev. Condens. Matter Phys. 6, 201 (2015)
  * [11] H. T. C. Stoof, Nature 440, 877 (2006).
  * [12] J. Li, arXiv:2009.12865.
  * [13] A. Shapere and F. Wilczek, Phys. Rev. Lett. 109, 160402 (2012).
  * [14] F. Wilczek, Phys. Rev. Lett. 109, 160401 (2012).
  * [15] K. Sacha, Time Crystals, (Springer, Cham, 2020).
  * [16] J. Zhang, P. W. Hess, A. Kyprianidis, P. Becker, A. Lee, J. Smith, G. Pagano, I.-D. Potirniche, A. C. Potter, A. Vishwanath, N. Y. Yao, and C. Monroe, Nature 543, 217 (2017).
  * [17] S. Choi, J. Choi, R. Landig, G. Kucsko, H. Zhou, J. Isoya, F. Jelezko, S. Onoda, H. Sumiya, V. Khemani, C. von Keyserlingk, N. Y. Yao, E. Demler, and M. D. Lukin, Nature 543, 221 (2017).
  * [18] S. Autti, V. B. Eltsov, and G. E. Volovik, Phys. Rev. Lett. 120, 215301 (2018).
  * [19] J. Rovny, R. L. Blum, and S. E. Barrett, Phys. Rev. Lett. 120, 180603 (2018).
  * [20] S. Pal, N. Nishad, T. S. Mahesh, and G. J. Sreejith, Phys. Rev. Lett. 120, 180602 (2018).
  * [21] N. Träger et al., arXiv:1911.13192.
  * [22] H. Kessler, P. Kongkhambut, C. Georges, L. Mathey, J. G. Cosme, and A. Hemmerich, arXiv:2012.08885.
  * [23] J. Smits, L. Liao, H. T. C. Stoof, and P. van der Straten, Phys. Rev. Lett. 121, 185301 (2018).
  * [24] J. Smits, H. T. C. Stoof, and P. van der Straten, New J. Phys. 22, 105001 (2020).
  * [25] J. Smits, H. T. C. Stoof, and P. van der Straten, arXiv:2103.03748
  * [26] H. T. C. Stoof, K. B. Gubbels, and D. B. M Dickerscheid, Ultracold Quantum Fields, (Springer, Berlin, 2009).
  * [27] J. W. Negele and H. Orland, Quantum Many-Particle Systems, (Addison-Wesley, New York, 1988).
  * [28] H. T. C. Stoof, J. Low Temp. Phys. 114, 11 (1999); J. Low Temp. Phys. 124, 431 (2001).
  * [29] A.O. Caldeira and A.J. Leggett, Phys. A 121, 587 (1983).
  * [30] E. P. Wigner, Phys. Rev. 40, 749 (1932).
  * [31] P. D. Drummond and A. D. Hardman, Europhys. Lett. 21, 279 (1993).
  * [32] D. Xiao, J. Shi, and Q. Niu, Phys. Rev. Lett. 95, 137204 (2005).
  * [33] F. Freimuth, R. Bamler, Y. Mokrousov, A. Rosch, Phys. Rev. B 88, 214409 (2013).
  * [34] R.W. MacCormack, 4th Aerodynamic Testing Conference, AIAA Paper 1969-354 (1969).