# Simulation of time crystal behavior for a few boson
chiral soliton model in a ring

Patrik Öhberg  SUPA, Institute of Photonics and Quantum Sciences, Heriot-Watt
University, Edinburgh EH14 4AS, United Kingdom  Ewan M. Wright  SUPA,
Institute of Photonics and Quantum Sciences, Heriot-Watt University, Edinburgh
EH14 4AS, United Kingdom  James C. Wyant College of Optical Sciences,
University of Arizona, Tucson, Arizona 85721, USA

###### Abstract

We present numerical simulations for a chiral soliton model with $N=2,3$
bosons in a ring, this being a few-particle version of our previous mean-field
model for a quantum time crystal OW19 . Following Syrwid, Kosior, and Sacha
(SKS), the notion is that a precise position measurement of one particle can
lead to spontaneous formation of a bright soliton that in a time crystal
should rotate intact for at least a few revolutions around the ring SKS_EPL .
In their work SKS find spontaneous formation of a soliton due to the position
measurement, but quantum fluctuations cause the soliton to subsequently decay
before it has a chance to perform even one revolution of the ring. Based on
this they conclude that time crystal dynamics are impossible. In contrast, for
our few boson chiral soliton model, and allowing for imprecise (weak)
measurements of the particle position, we show that time crystal behavior is
possible allowing for several revolutions of the spontaneously formed soliton
around the ring. We therefore argue that our chiral soliton model can realize
a quantum time crystal when weak position measurements are allowed for.

##  I Introduction

Since Wilczek’s seminal paper initiated the field of quantum time crystals in
2012 Wil12 , it has attracted intense interest, both fundamental and applied.
His model comprised $N$ bosons with attractive interactions in a ring that was
pierced by a flux tube, the key idea being that this model could yield a
ground state composed of a soliton that was rotating around the ring. Although
this model was subsequently shown to always yield a non-rotating ground state,
the time crystal idea had taken root. Now there is active research in both
continuous and discrete time crystals and for a variety of platforms including
optics, ultracold atoms, metamaterials, spin systems, and both closed and open
systems. For recent reviews of the area see SacZak18 ; Sacha20 ; ZalLukMon23 .

The issue undermining Wilczek’s initial proposal was nicely illustrated in a
paper by Syrwid, Kosior, and Sacha (SKS) in which they simulated N-bosons in a
ring that is pierced by a flux tube, and which was initiated from the
N-particle ground state which shows no signatures of rotation of localization
SKS_EPL . What they found was that precise position measurement of one
particle could lead to the spontaneous formation of a $(N-1)$-particle soliton
that is rotating, but that quantum fluctuations cause this soliton to
subsequently decay before it has a chance to perform even one revolution of
the ring. SKS presented results for particle numbers in the range $N=10-60$,
and the lifetime $t_{c}$ of the soliton was found to decrease with $N$. We
observed the same decay for a non-chiral version of our model with $N=3$, both
for precise and imprecise position measurements (with an uncertainty in the
position measurement), the lifetime being longer for the imprecise
measurement. The conclusion of SKS remains intact that Wilczek’s model cannot
show time crystal behavior SKS_EPL .

In 2019 we proposed a mean-field chiral soliton model for a quantum time
crystal OW19 . This paper was followed by discussion in the literature with
SKS SKS_comment -OW_lack , culminating in their paper highlighted in the
previous paragraph SKS_EPL . The goal of the present paper is to examine a few
boson $(N=2,3)$ limit of our previous chiral soliton model to numerically
assess whether time crystal behavior is possible. In particular, we find the
$N$-boson ground state, perform a position measurement of one particle, and
see if spontaneous formation of the resulting $(N-1)$-boson soliton can
persist over several revolutions of the ring as required for time crystal
behavior. We find that for an imprecise or weak position measurement time
crystal behavior is possible in our chiral soliton model, whereas for a
precise measurement quantum fluctuations cause the soliton to decay in
accordance with SKS SKS_EPL . We therefore argue that our chiral model can
realize a quantum time crystal when weak position measurements are given due
consideration.

The remainder of this paper is organized as follows: Section II describes the
few boson version of our chiral soliton model along with a description of the
parameters involved. The case of $N=2$ particles is covered in Sec. III and
shows the ideas involved in a transparent manner. Section IV contains the main
results and covers the three particle case that demonstrates that our few
boson chiral soliton model can display time crystal behavior in which a two-
boson soliton state can execute several revolutions of the ring. Summary and
conclusions are given in Sec. V.

##  II Chiral soliton model

###  II.1 N-particle Schrödinger equation

We consider a system of $N$ scalar bosons of mass $m$ located
circumferentially in a ring of radius $R$. The Schrödinger equation for the
N-particle wavefunction $\psi(\theta_{1},\theta_{2},\ldots\theta_{N},t)$ is in
scaled units Jackiw

$i{\partial\Psi\over\partial
t}=\sum_{j=1}^{N}\left(-i{\partial\over\partial\theta_{j}}-A(\theta_{1},\theta_{2},\ldots\theta_{N})\right)^{2}\Psi+{g\over
2}\rho(\theta_{1},\theta_{2},\ldots\theta_{N})\Psi.$ (1)

Here $\theta_{j}=[-\pi,\pi],j=1,2,\ldots N$ are angular coordinates around the
ring, and time is in units of $2mR^{2}/\hbar$. The scaled gauge potential is
given by

$A(\theta_{1},\theta_{2},\ldots\theta_{N})=A^{(0)}+\kappa\rho(\theta_{1},\theta_{2},\ldots\theta_{N}),$
(2)

where $A^{(0)}$ represents the single-particle contribution to the gauge
potential, and $\rho(\theta_{1},\theta_{2},\ldots\theta_{N})$ describes the
density-dependent contribution to the gauge potential. The dimensionless
parameter $\kappa$ characterizes the strength of the density-dependent chiral
gauge potential, and the parameter $g$ controls the density-dependent and non-
chiral contribution to the system energy. The density-dependent gauge
potential is built from the combinations of the two-particle interactions
between the composite bosonic particles, the two-particle interaction being
denoted $\eta(x)$, with $x$ the separation between the two particles. In
particular, we set

$\rho(\theta_{1},\theta_{2},\ldots\theta_{N})=\sum_{i=1}^{N}\sum_{j>i}^{N}\eta(\theta_{i}-\theta_{j}),$
(3)

where $\rho(\theta_{1},\theta_{2},\ldots\theta_{N})$ is symmetric under
particle exchange by virtue of the symmetry $\eta(x)=\eta(-x)$. Since the
density and corresponding density-dependent gauge potential are symmetric
under particle exchange, which breaks Galilean invariance, we expect chiral
dynamics to appear in this model.

###  II.2 Model parameters

In our previous mean-field model the chiral soliton model was derived starting
from the spinor Schrödinger equation for a system of two-level bosonic atoms
trapped in a ring that are dipole coupled using a laser beam carrying orbital
angular momentum (OAM) characterized by the winding number $\ell$. By
preparing the atoms in the appropriate dressed state a scalar field analysis
may be obtained in which the atoms are subject to a vector potential of the
form in Eq. (2) with $A^{(0)}={\ell\over 2}$. Moreover, the density-dependent
gauge potential results from a collision-induced detuning, with the
consequence that the dimensionless parameter $\kappa\propto\ell$. This means
that if the winding number of the laser beam $\ell=0$, then $\kappa=0$ and
there will be no chiral dynamics and no possibility for time crystal behavior.
The term proportional to the dimensionless parameter $g$ describes energy-
shifts due to many-body interactions.

For numerical purposes it is impractical to consider s-wave many-body
interactions with $\eta(x)=\delta(x)$. Instead we here use the non-local model
PriValWri23

$\eta(x)=Q\cos^{2q}(x/2),\quad q=1,2,3,\ldots,$ (4)

where $Q$ is a normalization constant such that
$\int_{-\pi}^{\pi}\eta(x)dx=1$. This allows for numerical simulations to be
performed in conjunction with checking that the results obtained are not
particularly sensitive to the choice of $q$. In the simulations presented we
set $q=50$ but increasing this to $q=100$ made little difference, meaning that
the non-locality is not a key ingredient for the results we obtain. In this
sense our model is local for all practical purposes.

##  III Two-particle case

As a first example we consider the case of two bosons in a ring. This example
is over-simplified in that if one particle is measured this leaves only one
particle on the ring, so that the subsequent dynamics is linear. This example
is nonetheless quite transparent and illustrates the ideas involved.

###  III.1 Two-particle Schrödinger equation

The density for the two bosonic particles may be written as

$\displaystyle\rho(\theta_{1},\theta_{2})$ $\displaystyle\equiv$
$\displaystyle\eta(\theta_{1}-\theta_{2}),$ (5)

and the Schrödinger equation for $\Psi(\theta_{1},\theta_{2},t)$ then becomes
PriValWri23

$i{\partial\Psi\over\partial
t}=\sum_{j=1}^{2}\left(-i{\partial\over\partial\theta_{j}}-{\ell\over
2}-\kappa\eta(\theta_{1}-\theta_{2}))\right)^{2}\Psi+{g\over
2}\eta(\theta_{1}-\theta_{2})\Psi.$ (6)

Changing to center-of-mass (COM) and relative coordinates using
$s=(\theta_{1}+\theta_{2})/2,x=(\theta_{1}-\theta_{2})$, the transformed
Schrödinger equation for $\Psi(x,s,t)$ becomes

$i{\partial\Psi\over\partial t}=\left[-2{\partial^{2}\over\partial
x^{2}}-{1\over 2}{\partial^{2}\over\partial s^{2}}+2i\left({\ell\over
2}+\kappa\eta(x)\right){\partial\over\partial s}+2\left({\ell\over
2}+\kappa\eta(x)\right)^{2}+{g\over 2}\eta(x)\right]\Psi.$ (7)

###  III.2 Two-particle ground state

We seek a two-particle stationary state using the ansatz PriValWri23

$\Psi(x,s,t)=e^{-i\varepsilon t+ips}\varphi(x),$ (8)

with $\varepsilon$ the scaled energy of the two-particle system, and $p$ an
integer that is the winding number associated with the center-of-mass OAM of
the two particles. Substituting this ansatz into Eq. (7) yields

$\left[-{d^{2}\over dx^{2}}+\left({(p-\ell)\over
2}-\kappa\eta(x)\right)^{2}+{g\over
2}\eta(x)\right]\varphi(x)={\varepsilon\over 2}\varphi(x),$ (9)

where $\varepsilon$ is the scaled energy of the two-particle system.

Figure 1: (a) Ground state profile $|\varphi(x)|^{2}$ versus $x$ and $g$, and
(b) scaled energy $\varepsilon$ versus $g$. Parameters used are $N=2,\ell=2$,
$\kappa=0.2$, and $q=50$. Figure 2: (a) Scaled energy $\varepsilon$ versus COM
winding number $p$, and (b) ground state profile $|\varphi(x)|^{2}$ versus
$x$. Parameters are the same as Fig. 1 along with $g=-2$.

We have solved Eq. (9) numerically for a range of parameters using the
imaginary time method LehToiElo2007  applied to Eq. (7), in particular to
obtain the ground state $\varphi(x)$ and associated scaled energy
$\varepsilon$. Figure 1 shows numerical results for parameters
$\ell=2,\kappa=0.2,q=50$. Specifically, Fig. 1(a) shows a color coded plot of
the ground state profile $|\varphi(x)|^{2}$ versus $x$ for a variety of $g<0$
along the horizontal axis, and plot (b) shows the scaled energy $\varepsilon$
over the same range of $g$. The ground state profile $|\varphi(x)|^{2}$ is
normalized to unity on-axis so that the profile can be discerned for each
value of $g$. What is evident is that for $g>-0.5$ the ground state is uniform
to a high degree, whereas for $g<-0.5$ the ground state profile is localized
in terms of the relative coordinate $x$.

In the remainder of this section we choose $g=-2$ as a illustrative example,
though similar results could be obtained for different combinations of
parameters. Figure 2 shows some more detail of the solution for this choice.
In particular, Fig. 2(a) shows the scaled energy versus $p$ and reveals that
the ground state occurs for $p=\ell=2$, and (b) shows a line plot of the
ground state profile $|\varphi(x)|^{2}$ (properly normalized) versus $x$ for
$p=2$. Furthermore, the width (half-width at half-maximum) of the ground state
profile in Fig. 2(b) with respect the the relative coordinate $x$ of the two
particles is $w\sim 1$, the parameter $w$ characterizing the expected width of
any localized two-particle solutions that may arise.

###  III.3 Time crystal dynamics from measurement

To proceed we first recognize that the ground state represented by the
probability density in Fig. 2(b) does not correspond to a localized state on
the ring since $|\varphi(x)|^{2}$ depends only on the relative coordinate $x$,
and there is no information concerning the COM position $s$ of the particle.
We can probe further and reconstruct the ground state two-body wave function
$\Psi(\theta_{1},\theta_{2})$ (at $t=0$), where
$|\Psi(\theta_{1},\theta_{2})|^{2}$ is shown in Fig. 3(a). From this it is
straightforward to see that the probability density of finding the particle at
$\theta_{1}$ for any value of $\theta_{2}$,

$P(\theta_{1})=\int_{-\pi}^{\pi}~{}d\theta_{2}|\Psi(\theta_{1},\theta_{2})|^{2}={1\over
2\pi},$ (10)

is a flat distribution, once again verifying that for the ground state the
solution does not represent a localized solution that is pinned to a given COM
position. In their work SKS SKS_EPL  argued that a precise position
measurement of one particle could lead initially to spontaneous formation of a
localized solution, and we pursue and expand that concept here for our chiral
soliton model. In particular, we assume that the position of particle $2$ is
measured yielding a value $\theta_{2}^{(0)}$, but allow for the measurement to
be imprecise, so that the wave function for the remaining particle may be
written as

$\psi(\theta_{1},t=0)=\Psi(\theta_{1},\theta_{2}^{(0)})=\int_{-\pi}^{\pi}d\theta_{2}~{}G(\theta_{2}-\theta_{2}^{(0)})\Psi(\theta_{1},\theta_{2}),$
(11)

where the function $G(\theta_{2}-\theta_{2}^{(0)})$ reflects the uncertainty
of the position measurement of particle $2$, and is a distribution which obeys
the periodicity of the ring, and with an adjustable width $\Delta\theta$. Here
for illustration we employ the following distribution

$G(\theta-\theta^{(0)})={\cal N}\cos^{2n}\left({\theta-\theta^{(0)}\over
2}\right),\quad n=1,2,3,\ldots$ (12)

with ${\cal N}$ a normalization constant. In the limit $n\gg 1$ this
distribution approaches a delta-function which corresponds to the precise
measurement model adopted in the paper of SKS, whereas for $n=1$ the
measurement uncertainty is around $\Delta\theta\sim 2$. After the measurement,
precise or imprecise, the wave function $\psi(\theta_{1},t=0)$ (suitably
normalized) for the remaining particle is given by Eq. (11), and subsequent
time development is governed by the single-particle Schrödinger equation

$i{\partial\psi\over\partial
t}=\left(-i{\partial\over\partial\theta_{1}}-{\ell\over 2}\right)^{2}\psi.$
(13)

Figure 3(b) shows a simulation of the quantum dynamics according to Eq. (13)
following a precise measurement for which we set $\theta_{2}^{(0)}=0$ and
$n=50$ so that $\Delta\theta\sim 0.2$. In this case the uncertainty is much
less than the expected width $w\sim 2$ for a two-particle localized state, so
the measurement is expected to be very disruptive. The figure shows a color
coded plot of the probability density $|\psi(\theta_{1},t)|^{2}$ versus $t$
and $\theta_{1}$. Then, for example, if we select
$\theta_{1}=\theta_{2}^{(0)}=0$, and trace the probability density in time
along the horizontal axis, we see what might be a periodic variation
indicative of a time crystal. But this is an illusion based on us choosing
$\theta_{1}=\theta_{2}^{(0)}$, a quantity we cannot know a priori, and if we
chose another $\theta_{1}$ at which to monitor the time variation of the
probability density we would generally find another much noisier variation
with time. Another way to see this from Fig. 3(b) is that the spatial density
profile $|\psi(\theta_{1},t=0)|^{2}$ becomes rapidly distorted as time
increases, not simply displaced as expected for a time crystal. This is in
perfect keeping with the findings of SKS that too precise a position
measurement of one particle can produce a localized state initially but this
will become rapidly dispersed before the initial localized solution can
exhibit a single rotation around the ring.

For time crystal behavior what is required is that the center of the initial
localized state basically translates in time and persists in shape for several
rotations around the ring. The speed and direction of the rotation are
intimately tied to the chiral gauge potential proportional to $\kappa$. Figure
4 shows an example of time crystal behavior using $n=1$ and $\Delta\theta\sim
2$, which means a much weaker measurement than in Figure 3(b). In this case
the uncertainty is greater than the expected width $w\sim 1$ for a two-
particle localized state, so the measurement is expected to be less
disruptive. This figure shows that the initial localized solution remains
largely intact while it rotates, with a velocity $v\sim-1$, and it persists
over several rotations around the ring. The initial position of the centroid
of the localized solution is set by the outcome of the initial measurement,
$\theta^{(0)}=0$ in Fig. 4(b). In addition, if one selects a value for
$\theta_{1}\neq 0$ in Fig. 4 the resulting time variation of the probability
density is then largely same for any other value of $\theta_{1}$, just
displaced in time. This example therefore exhibits time crystal behavior to a
good approximation. This may seem counter-intuitive at first, since one would
expect a localized initial state to expand and eventually interfere with
itself in the ring giving rise to complicated dynamics. However, it is
possible to construct special solutions to the single particle Schrödinger
equation for a ring, that are localised non-dispersive solutions with a non-
zero angular momentum, i.e. they rotate. Such solutions have been found by Qin
Qin  and are of the form

$\Phi(\theta,t)=C_{0}\cos(m\theta-nt)e^{i(m\theta-{1\over 2}(m^{2}+n^{2})t)},$
(14)

where $m$ and $n$ are integers. From Eq. (14) we notice that $\Phi(\theta,t)$
is a travelling localized solution with velocity $v=n/m$. For $m=-n=1$ we find
$v=-1$, in qualitative agreement with the results in Fig. 4. In addition,
because of the linearity of the single-particle Schröinger equation, we can
always construct superpositions of different states such that we have a non-
zero background. This is indeed the type of solution we generally find after
measuring one particle with sufficiently large uncertainty. We therefore
conclude that when extracting information about the two-particle ground state
by measuring the position of one particle with sufficiently large uncertainty,
the ground state is projected onto a single particle wave function which is
localised, shape preserving, and rotating in the ring.

Figure 3: (a) Reconstructed ground state probability density
$|\Psi(\theta_{1},\theta_{2})|^{2}$, and (b) simulation for a measurement with
$n=50$ that localizes particle $2$ with an uncertainty $\Delta\theta\sim 0.2$.
Picking a given value of $\theta_{1}$ one can trace the time evolution of the
system and it does not reflect time crystal behavior. Parameters are the same
as as in Fig. 1 along with $g=-2$. Figure 4: Simulation for an imprecise
measurement with $n=1$ and an uncertainty $\Delta\theta\sim 2$ that couples
the remaining particle after the measurement of particle 2 to a non-dispersing
solution for the ring. Picking a given value of $\theta_{1}$ one can trace the
time evolution of the system and it is quite periodic, consistent with time
crystal behavior over several periods. Parameters are the same as as in Fig. 1
along with $g=-2$.

The basic message from our two-particle simulations is that using a weak
measurement of particle position can lead to time crystal behavior for the
remaining particle on the ring. This is the case since the weak measurement
can dominantly excite a non-dispersive and rotating solution on the ring,
whereas the more precise position measurement excites a broader wavepacket of
the ring modes that leads to the observed decay and lack of time crystal
behavior. The benefit of the two-particle case is that it reveals the role
played by the position measurement uncertainty in a transparent manner.

##  IV Three-particle case

Although the two-particle case provides some insight it does not speak to the
bigger issue of whether the spontaneously formed and rotating solution after
the measurement can persist even for a multi-particle state, so the rotating
solution may be viewed even approximately as a rotating soliton. For this
reason we next consider the three-particle case so that two particles are left
on the ring after the measurement, with the possibility of exciting the two-
particle analogue of a soliton.

###  IV.1 Three-particle Schrödinger equation

The density for the $N=3$ bosonic particles may be written as

$\rho(\theta_{1},\theta_{2},\theta_{3})\equiv\eta(\theta_{1}-\theta_{2})+\eta(\theta_{1}-\theta_{3})+\eta(\theta_{2}-\theta_{3}),$
(15)

which by construction is symmetric under exchange of any pair of particle
coordinates. The Schrödinger equation for
$\Psi(\theta_{1},\theta_{2},\theta_{3},t)$ may then be written as

$i{\partial\Psi\over\partial
t}=\sum_{j=1}^{3}\left(-i{\partial\over\partial\theta_{j}}-{\ell\over
2}-\kappa\rho(\theta_{1},\theta_{2},\theta_{3}))\right)^{2}\Psi+{g\over
2}\rho(\theta_{1},\theta_{2},\theta_{3})\Psi.$ (16)

To proceed it is useful to use the following Jacobi coordinates Jacobi
appropriate to this three-body problem

$s={1\over 3}(\theta_{1}+\theta_{2}+\theta_{3}),\quad
x=(\theta_{1}-\theta_{2}),\quad y=(\theta_{1}+\theta_{2}-2\theta_{3}),$ (17)

and inverting the Jacobi coordinates we obtain

$\displaystyle\theta_{1}$ $\displaystyle=$ $\displaystyle s+{1\over
6}y+{1\over 2}x,$ $\displaystyle\theta_{2}$ $\displaystyle=$ $\displaystyle
s+{1\over 6}y-{1\over 2}x,$ $\displaystyle\theta_{3}$ $\displaystyle=$
$\displaystyle s-{1\over 3}y.$ (18)

From these results we find

$\displaystyle(\theta_{1}-\theta_{2})$ $\displaystyle=$ $\displaystyle x$
$\displaystyle(\theta_{1}-\theta_{3})$ $\displaystyle=$ $\displaystyle{1\over
2}y+{1\over 2}x,$ $\displaystyle(\theta_{2}-\theta_{3})$ $\displaystyle=$
$\displaystyle{1\over 2}y-{1\over 2}x,$ (19)

so that the density $\rho$ in these Jacobi coordinates becomes

$\rho(x,y)=\eta\left({x\over 2}\right)+\eta\left({y\over 4}+{x\over
4}\right)+\eta\left({y\over 4}-{x\over 4}\right),$ (20)

which is independent of the COM coordinate $s$. The transformed Schrödinger
equation for $\Psi(x,y,s,t)$ then becomes

$i{\partial\Psi\over\partial t}=-\left(2{\partial^{2}\over\partial
x^{2}}+6{\partial^{2}\over\partial y^{2}}+{1\over 3}{\partial^{2}\over\partial
s^{2}}\right)\Psi-2i\kappa\rho(x,y){\partial\Psi\over\partial
s}+3[\kappa\rho(x,y)]^{2}\Psi+{g\over 2}\rho(x,y)\Psi.$ (21)

###  IV.2 Three-particle ground state

We seek a three-particle stationary state using the ansatz

$\Psi(x,y,s,t)=e^{-i\varepsilon t+ips}\varphi(x,y),$ (22)

where $p$ is the integer valued winding number associated with the center-of-
mass OAM of the three particles. Substituting the above ansatz into Eq. (21)
yields

$\varepsilon\varphi(x,y)=-\left(2{\partial^{2}\over\partial
x^{2}}+6{\partial^{2}\over\partial y^{2}}\right)\varphi+3\left[{p\over
3}+\kappa\rho(x,y)\right]^{2}\varphi+{g\over 2}\rho(x,y)\varphi,$ (23)

where $\varepsilon$ is the scaled energy of the three-particle system.

We have solved Eq. (23) numerically for a range of parameters using the
imaginary time method LehToiElo2007 , in particular to obtain the ground state
$\varphi(x,y)$ and associated scaled energy $\varepsilon$. Figure 5 shows
numerical results for parameters $\ell=3,g=-2,\kappa=0.2,q=50$. Specifically,
Fig. 5(a) shows the scaled energy versus $p$ and reveals that the ground state
occurs for $p=\ell=3$, and (b) shows a line plot of the ground state profile
$|\varphi(x,y)|^{2}$ versus $x$ and $y$ for $p=3$.

Figure 5: (a) Scaled energy $\varepsilon$ versus COM winding number $p$, and
(b) ground state profile $|\varphi(x,y)|^{2}$ for $p=\ell=3$. Parameters are
$\ell=3,g=-2,\kappa=0.2,q=50$.

###  IV.3 Time crystal dynamics from measurement

To proceed we first recognize that the ground state represented by the
probability density in Fig. 5(b) does not correspond to a localized state on
the ring since $|\varphi(x,y)|^{2}$ depends only on the Jacobi coordinates $x$
and $y$, and there is no information concerning the COM position $s$ of the
particle. We can also reconstruct the ground state three-body wave function
for one coordinate fixed, say $\theta_{3}=0$, and
$|\Psi(\theta_{1},\theta_{2},\theta_{3}=0)|^{2}$ (at $t=0$) is shown in Fig.
6(a): We note that the width of this ground state density profile is around
$w\sim 1$, that provides a characteristic scale for localized three-particle
states that may appear in the system.

A few words are in order about transforming from $\varphi(x,y)$ to
$\Psi(\theta_{1},\theta_{2},\theta_{3})$. First we note that the Jacobi
coordinates $x$ and $y$ can be re-expressed as

$x=\theta_{1}^{\prime}-\theta_{2}^{\prime},\quad
y=\theta_{1}^{\prime}+\theta_{2}^{\prime},$ (24)

where $\theta_{1,2}^{\prime}=\theta_{1,2}-\theta_{3}$. This means that, for a
given $\theta_{3}$, the transformation between the $\varphi(x,y)$ in the
Jacobi coordinates $x$ and $y$, and
$\Psi(\theta_{1}^{\prime},\theta_{2}^{\prime})$ in the coordinates
$\theta_{1}^{\prime}$ and $\theta_{2}^{\prime}$ can be performed using a
$45^{\circ}$ rotation. We note that at this point the wave function
$\Psi(\theta_{1}^{\prime},\theta_{2}^{\prime})$ should be symmetrized if
needed. One can then transform from
$\Psi(\theta_{1}^{\prime},\theta_{2}^{\prime})$ to
$\Psi(\theta_{1},\theta_{2},\theta_{3})$ using simple linear displacements
$\theta_{1,2}=\theta_{1,2}^{\prime}+\theta_{3}$.

In their work SKS argued that a precise position measurement of one particle
could lead initially to spontaneous formation of a localized solution SKS_EPL
. We assume that the position of particle $3$ is measured yielding a value
$\theta_{3}^{(0)}$, but allow for the measurement to be imprecise, so that the
wave function for the remaining particles may be written as

$\psi(\theta_{1},\theta_{2},t=0)=\Psi(\theta_{1},\theta_{2},\theta_{3}^{(0)})=\int_{-\pi}^{\pi}d\theta_{3}~{}G(\theta_{3}-\theta_{3}^{(0)})\Psi(\theta_{1},\theta_{2},\theta_{3}),$
(25)

where the function $G(\theta_{3}-\theta_{3}^{(0)})$ reflects the uncertainty
of the position measurement, and we again use the model in Eq. (12). After the
measurement the wave function $\psi(\theta_{1},\theta_{2},t=0)$ (suitably
normalized) for the remaining particles is given by Eq. (25), and subsequent
time development is governed by the two-particle Schrödinger equation

$i{\partial\psi\over\partial
t}=-\left({\partial^{2}\over\partial\theta_{1}^{2}}+{\partial^{2}\over\partial\theta_{2}^{2}}\right)\psi+2i\left({\ell\over
2}+\kappa\rho\right)\left({\partial\over\partial\theta_{1}}+{\partial\over\partial\theta_{2}}\right)\psi+2\left({\ell\over
2}+\kappa\rho\right)^{2}\psi+{g\over 2}\rho\psi,$ (26)

where $\rho=\eta(\theta_{1}-\theta_{2})$ is the scaled density for the two
remaining particles.

Figure 6(b) shows a simulation of the quantum dynamics according to Eq. (26)
following a precise measurement for which we set $\theta_{3}^{(0)}=0$ and
$n=50$ so that $\Delta\theta\sim 0.2$. In this case the uncertainty is much
less than the expected width $w\sim 1$ for a three-particle localized state,
so the measurement is expected to be very disruptive. The figure shows a color
coded plot of the probability density for measuring a second particle

$P(\theta_{1},t)=\int_{-\pi}^{\pi}d\theta_{2}|\psi(\theta_{1},\theta_{2},t)|^{2},$
(27)

versus $t$ and $\theta_{1}$. Then, for example, if we select
$\theta_{1}=\theta_{3}^{(0)}=0$, and track the probability density in time
along the horizontal axis, we see what might be a periodic variation
indicative of a time crystal. But this is an illusion based on us choosing
$\theta_{1}=\theta_{3}^{(0)}$, a quantity we cannot know a priori, and if we
chose another $\theta_{1}$ at which to monitor the time variation of the
probability density we would generally find another variation with time.
Another way to see this from Fig. 6(b) is that the spatial density profile
$|\psi(\theta_{1},t=0)|^{2}$ becomes rapidly distorted as time increases. This
is in perfect keeping with the findings of SKS that too precise a position
measurement of one particle can produce an localized state initially but this
will become rapidly dispersed before the initial localized solution can
exhibit a single rotation around the ring.

For time crystal behavior what is required is that the center of the initial
localized state basically translates in time and persists in shape for several
rotations around the ring. The speed and direction of the rotation are
intimately tied to the chiral gauge potential proportional to $\kappa$. Figure
7 shows an example of time crystal behavior using $n=1$ and $\Delta\theta\sim
2$, which means a much weaker measurement than in Figure 6(b): In this case
the uncertainty is greater than the expected width $w\sim 1$ for a three-
particle localized state, so the measurement is expected to be less
disruptive. This figure shows that the initial localized solution remains
largely intact while it rotates, and it persists over several rotations around
the ring. In addition, if one selects a value for $\theta_{1}$ in Fig. 7 the
resulting time variation of the probability density is then largely same for
any other value of $\theta_{1}$, just displaced in time. This example
therefore exhibits time crystal behavior to a good approximation.

To summarise, we deduce that the initial measurement of particle $3$ leads to
spontaneous formation of a localized state that is a two-particle analogue of
a mean-field chiral soliton Jackiw , and that this two-particle soliton can
continue rotating whilst intact over several rotation periods. This can happen
here since we considered a weak measurement as opposed to the precise
measurement that leads to decay of the rotating soliton, see Fig. 6(b).

Figure 6: (a) Reconstructed ground state probability density
$|\Psi(\theta_{1},\theta_{2},\theta_{3}=0)|^{2}$, and (b) simulation for a
measurement with $n=50$ that localizes particle $3$ with an uncertainty
$\Delta\theta\sim 0.2$. Picking a given value of $\theta_{1}$ one can trace
the time evolution of the system and it does not reflect time crystal
behavior. Parameters are the same as as in Fig. 5 along with $g=-2$. Figure 7:
Simulation for an imprecise measurement with $n=1$ and an uncertainty
$\Delta\theta\sim 2$ that couples the remaining two particles after the
measurement of particle $3$ to the two-particle limit of the chiral soliton.
Picking a given value of $\theta_{1}$ one can trace the time evolution of the
system and it is quite periodic, consistent with time crystal behavior over
several periods. Parameters are the same as as in Fig. 5 along with $g=-2$.

##  V Summary and conclusions

In 2019 we proposed a mean-field chiral soliton model for a quantum time
crystal OW19 , and this proposal met with some criticism in the literature
SKS_EPL ,SKS_comment -OW_lack . The goal of the present paper was to examine a
few boson $(N=2,3)$ limit of our previous chiral soliton model to numerically
assess whether time crystal behavior is possible in the context of the type of
calculation in Ref. SKS_EPL . In particular, we initialize the system in the
$N$-boson ground state, perform a position measurement of one particle, and
see if spontaneous formation of the resulting $(N-1)$-boson localized state
(soliton for $N=2$) can persist over several revolutions of the ring as
required for time crystal behavior. We find that for an imprecise or weak
position measurements time crystal behavior is possible in our chiral soliton
model, whereas for a precise measurement quantum fluctuations cause the
soliton to decay in accordance with SKS_EPL .

Clearly the small particle numbers and parameters used here are not compatible
with experiments, but are rather intended to show that time crystal behavior
is not ruled out in our model. We remark that the same behavior could be seen
for other combinations of parameters as long as a weak measurement is
employed, with $\Delta\theta>w$. Restricting to a few bosons allowed for
numerical simulations to be performed for our chiral soliton model without
undue approximations, but our conclusions may still be of relevance in the
more general case with large $N$. According to the Wilczek model analyzed in
Ref. SKS_EPL  the soliton lifetime $t_{c}$ is expected to scale with $N$.
Assuming that this scaling also applies to our chiral model, if $N=3$ can show
time crystal behavior then it should be expected for larger $N$ as the
lifetime increases. Therefore we argue that the chiral model can realize a
quantum time crystal when weak position measurements are incorporated into the
scheme.

Figure 8: Simulation for an imprecise measurement with $n=1$ and an
uncertainty $\Delta\theta\sim 2$ that couples the remaining particle after the
measurement of particle $2$ to a linear non-dispersive wave. Picking a given
value of $\theta_{1}$ one can trace the time evolution of the system and it is
quite periodic, consistent with time crystal behavior over several periods.
Parameters are the same as in Fig. 5 along with $g=-2$.

Finally, we discuss the issue of whether our simulations point to genuine time
crystal behavior. This would require that the periodic motion seen in Fig. 7
would persist for all times for $N=3$, and we cannot be assured that this is
the case based on numerical solutions. What we can say for sure is the time
crystal behavior can persist over several revolutions of the ring if imprecise
measurements are allowed for. Specifically, the first position measurement
sets the time crystal in motion and establishes a time reference by setting
the initial position of the $N=2$ soliton peak, but further measurements will
be required to monitor if the time crystal continues to exhibit periodic
motion. Indeed, for $N=3$ we can go one step further and ask whether time
crystal behavior would persist if we performed a second imprecise measurement,
of say particle $2$, acting on the two-particle state in Fig. 7. Figure 8
shows a simulation for an imprecise measurement with $n=1$ and an uncertainty
$\Delta\theta\sim 2$, along with $\theta_{2}^{(0)}$ chosen consistent with the
probability density for such a measurement between $[-\pi,\pi]$, that couples
the remaining particle at $t=20$ to a linear non-dispersing state after the
measurement of particle $2$. Then picking a given value of $\theta_{1}$ one
can trace the time evolution of the system and it is quite periodic,
consistent with time crystal behavior over several periods. For our $N=3$ case
this validates that the time crystal behavior can survive more than one
measurement. In future publications we plan to extend this initial numerical
study to larger particle numbers, both numerically and analytically, and
thereby move the project towards experimentally accessible parameter ranges.

##  VI Acknowledgements

We acknowledge helpful discussions with Brian Anderson, Erika Andersson,
Stewart Lang, Joel Priestley, and Gerard Valentí-Rojas.

## References

  * (1) F. Wilczek, “Quantum time crystals,” Phys. Rev. Lett. 109, 160401 (2012).
  * (2) K. Sacha and J. Zakrzewski, “Time crystals: A review,” Rep. Prog. Phys. bf 81, (2018).
  * (3) K. Sacha, Time Crystals (Springer International Publishing, 2020).
  * (4) M. P. Zaletel, M. Lukin, C. Monroe, C. Nayak, F. Wilczek, and N. Y. Yao,“Colloquium: Quantum and classical discrete time crystals,” Rev. Mod. Phys. 95, 031001 (2023).
  * (5) A. Syrwid, A. Kosior, and K. Sacha, EPL 134, 66001 (2021).
  * (6) P. Öhberg and E. M. Wright, Phys. Rev. Lett. 123, 250402 (2019).
  * (7) A. Syrwid, A. Kosior, and K. Sacha, Phys. Rev. Lett. 124, 178901 (2020).
  * (8) P. Öhberg and E. M. Wright, Phys. Rev. Lett. 124, 178902 (2020).
  * (9) A. Syrwid, A. Kosior, and Sacha, Phys. Rev. Res. 2, 032038(R) (2020).
  * (10) P. Öhberg and E. M. Wright, Comment on “Lack of a time crystal in a chiral soliton model” by Syrwid, Kosior, and Sacha, archiv:2008.10940 [cond-mat.quant-gas].
  * (11) J. Priestley, G. Valenti-Rojas, E. M. Wright, and P. Öhberg, “The gauge coupled two-body problem in a ring,” J. Phys. A: Math. Theor. 56, 015305 (2023).
  * (12) See, for example, L. Lehtovaara, J. Toivanen, and J. Eloranta, J. of Comp. Phys. 221, 148 (2007).
  * (13) J. Qin, “Non-spreading matter-wave packets in a ring,” Phys. Scr. 94, 115402 (2019).
  * (14) F. Calogero, “Solution of a three body problem in one dimension,” J. Math Phys. 10, 2191 (1969).
  * (15) U. Aglietti, L. Griguolo, R. Jackiw, S.-Y. Pi, and D. Seminara , “Anyons and chiral solitons on a line,” Phys. Rev. Lett. 77, 4406 (1996).