# Twist and Turn Squeezing in a Multi-Mode Bose-Einstein Condensate

Junbang Liu  Department of Quantum Science, Research School of Physics,
Australian National University, Canberra, Ngunnawal Country, Australia.
Thomas Bartlett  Department of Quantum Science, Research School of Physics,
Australian National University, Canberra, Ngunnawal Country, Australia.
Joseph Hope  Department of Quantum Science, Research School of Physics,
Australian National University, Canberra, Ngunnawal Country, Australia.  Simon
Haine  Department of Quantum Science, Research School of Physics, Australian
National University, Canberra, Ngunnawal Country, Australia.
<EMAIL_ADDRESS>
###### Abstract

Here we examine the generation of Twist and Turn (TNT) Squeezing in a large
atom-number Bose-Einstein Condensate for the purposes of generating quantum-
enhanced states for atom interferometry. Unlike previous analysis, we examine
situations where the multi-mode dynamics is significant, and cannot be
captured by a simple single-mode model. We find that in some regimes, with
careful choice of the rotation parameter, we can still obtain squeezing much
more rapidly than via one-axis twisting (OAT).

##  I Introduction

There is currently considerable interest in the production of entangled states
in Bose-Einstein condensates with the motivation of enhancing the sensitivity
of atom interferometers and atomic clocks Pezzè  _et al._ (2018); Szigeti  _et
al._ (2021). Without many-particle entanglement, the phase sensitivity of such
experiments is fundamentally constrained to the shot-noise limit (SNL)
$\Delta\phi=1/\sqrt{N}$ Giovannetti  _et al._ (2006); Pezzé and Smerzi (2009).
In recent years, experiments in atomic systems based on the one-axis twisting
(OAT) squeezing scheme of Kitagawa and Ueda Kitagawa and Ueda (1993); Mølmer
and Sørensen (1999) have demonstrated metrologically useful spin-squeezing
Esteve  _et al._ (2008); Leroux  _et al._ (2010), and sub-shot-noise phase
detection Gross  _et al._ (2010); Riedel  _et al._ (2010); Muessel  _et al._
(2014) in proof-of-principle experiments. However, typical experiments are
limited to only moderate quantum enhancement due to constraints on the state
preparation time imposed by dephasing Li  _et al._ (2008, 2009); Haine
(2018a), multi-mode dynamics Haine and Johnsson (2009); Haine  _et al._
(2014), or in the case of atomic gravimetry, expansion of the freely
propagating wave-packets Szigeti  _et al._ (2020). This leads to a degree of
quantum enhancement that is considerably less than the theoretical optimum.
Recently, a related method known as Twist and Turn (TNT) squeezing Micheli
_et al._ (2003); Sorelli  _et al._ (2019); Mirkhalaf  _et al._ (2018) has been
demonstrated Strobel  _et al._ (2014); Muessel  _et al._ (2015). The TNT
Hamiltonian, which is a specific case of the Lipkin-Meshkov-Glick Hamiltonian
Lipkin  _et al._ (1965) using the same nonlinear interactions as OAT with an
additional linear rotation, typically reaches larger degrees of quantum
enhancement for the same interaction time. As TNT is based on the same
interactions that leads to OAT squeezing, it can be implemented in the same
experimental set-ups.

As well as spin-squeezing, it has been shown that TNT dynamics is capable of
generating highly entangled states beyond the Gaussian regime Strobel  _et
al._ (2014); Muessel  _et al._ (2015). In this case, the metrologically useful
entanglement can be quantified via the quantum Fisher information (QFI)
Demkowicz-Dobrzanski  _et al._ (2015); Braunstein and Caves (1994); Paris
(2009); Tóth and Apellaniz (2014), and can be accessed via interaction based
readouts Davis  _et al._ (2016); Hosten  _et al._ (2016); Fröwis  _et al._
(2016); Macrì  _et al._ (2016); Linnemann  _et al._ (2016); Nolan  _et al._
(2017); Szigeti  _et al._ (2017); Anders  _et al._ (2018); Mirkhalaf  _et al._
(2018); Haine (2018b); Hayes  _et al._ (2018); Haine (2021).

So far, experimental demonstrations of TNT dynamics have been restricted to
samples of a few hundred atoms. One reason for this is that larger atom
numbers introduce multimode dynamics Li  _et al._ (2009); Haine and Johnsson
(2009); Opanchuk  _et al._ (2012); Haine  _et al._ (2014). These dynamics can
cause the two spin components to spatially separate. In the case of OAT, this
separation can significantly increase the rate of entangling dynamics Riedel
_et al._ (2010); Li  _et al._ (2008); Nolan and Haine (2018). However, for
efficient TNT dynamics, we require continuous spatial overlap to enable
continuous coupling between the two spin components. In this paper, we
investigate how this multimode dynamics effects the implementation of TNT
dynamics.

The structure of this paper is as follows: In section II we introduce the
theoretical model used. In section III we derive an effective single-mode
model from our multi-mode model, and recap ideal single-mode behaviour of both
the OAT and TNT hamiltonians. In section IV we identify three parameter
regimes, and discuss how multi-mode dynamics effects the ability to generate
spin-squeezing and entanglement in each of these regimes, before summarising
our findings in section V.

##  II Model

We consider an atomic BEC with two relevant atomic states $|a\rangle$ and
$|b\rangle$, confined in a potential $V(x)$. For example, these levels could
be the two hyperfine ground states of 87Rb, in which case
$|a\rangle\equiv|F=1,m_{F}=0\rangle$ and $|b\rangle\equiv|F=2,m_{F}=0\rangle$.
We also allow for continuous coupling between the two states via a radio
frequency or microwave transition, which provides the ‘turn’ in the TNT
dynamics. Introducing the bosonic field operators
$\hat{\psi}_{i}(\mathbf{r})$, which annihilate an atom of state $|i\rangle$
from position $\mathbf{r}$, and obey the usual bosonic commutation relations

$\displaystyle\left[\hat{\psi}_{i}(\mathbf{r})\,,\,\,\hat{\psi}^{\dagger}_{j}(\mathbf{r}^{\prime})\right]=\delta_{ij}\delta(\mathbf{r}-\mathbf{r}^{\prime})\,,$
(1)

the Hamiltonian describing the system is

$\displaystyle\hat{H}$
$\displaystyle=\sum_{j=a,b}\int\hat{\psi}^{\dagger}_{j}(\mathbf{r})\hat{H}_{0}\hat{\psi}_{j}(\mathbf{r})\,\,d^{3}\mathbf{r}$
$\displaystyle+\sum_{j,k=a,b}\frac{U_{ij}}{2}\int\hat{\psi}^{\dagger}_{j}(\mathbf{r})\hat{\psi}^{\dagger}_{k}(\mathbf{r})\hat{\psi}_{j}(\mathbf{r})\hat{\psi}_{k}(\mathbf{r})\,\,d^{3}\mathbf{r}$
$\displaystyle+\frac{\hbar\Omega}{2}\int\left(\hat{\psi}^{\dagger}_{a}(\mathbf{r})\hat{\psi}_{b}(\mathbf{r})+\hat{\psi}^{\dagger}_{b}(\mathbf{r})\hat{\psi}_{a}(\mathbf{r})\right)\,\,d^{3}\mathbf{r}\,,$
(2)

where

$\displaystyle\hat{H}_{0}$
$\displaystyle=\frac{-\hbar^{2}}{2m}\bm{\nabla}^{2}+V(x)$ (3)

is the single particle Hamiltonian, and

$U_{ij}=\frac{4\pi\hbar^{2}a_{ij}}{m}$ (4)

is the inter-particle interactions between atoms in states $|i\rangle$ and
$|j\rangle$, parameterised by the $s$-wave scattering length $a_{ij}$ Chin
_et al._ (2010). The final term in Eq. (2) describes continuous coupling
between states $|a\rangle$ and $|b\rangle$ with Rabi frequency $\Omega$. We
assume that all atoms are initially in state $|a\rangle$ in the ground
motional state before applying a $\frac{\pi}{2}$ pulse at $t=0$ to coherently
couple 50% of the population to state $|b\rangle$, such that the state of each
atom is $\frac{1}{\sqrt{2}}\left(|a\rangle+|b\rangle\right)$. The system then
evolves under Eq. (2).

The goal of the scheme is to produce an entangled state, which when used as
the input state of an atom interferometer, is capable of providing
sensitivities better than the SNL. Atom interferometery is best described via
the SU(2) collective pseudo-spin operators defined by

$\displaystyle\hat{J}_{x}$
$\displaystyle=\frac{1}{2}\int\left(\hat{\psi}^{\dagger}_{b}(\mathbf{r})\hat{\psi}_{a}(\mathbf{r})+\hat{\psi}^{\dagger}_{a}(\mathbf{r})\hat{\psi}_{b}(\mathbf{r})\right)\,d^{3}\mathbf{r}\,,$
(5a) $\displaystyle\hat{J}_{y}$
$\displaystyle=-\frac{i}{2}\int\left(\hat{\psi}^{\dagger}_{b}(\mathbf{r})\hat{\psi}_{a}(\mathbf{r})-\hat{\psi}^{\dagger}_{a}(\mathbf{r})\hat{\psi}_{b}(\mathbf{r})\right)\,d^{3}\mathbf{r}\,,$
(5b) $\displaystyle\hat{J}_{z}$
$\displaystyle=\frac{1}{2}\int\left(\hat{\psi}^{\dagger}_{a}(\mathbf{r})\hat{\psi}_{a}(\mathbf{r})-\hat{\psi}^{\dagger}_{b}(\mathbf{r})\hat{\psi}_{b}(\mathbf{r})\right)\,d^{3}\mathbf{r}\,.$
(5c)

Assuming a simple Mach-Zehnder interferometer composed of a
$\frac{\pi}{2}-\pi-\frac{\pi}{2}$ pulse sequence, if measurements capable of
resolving the mean and variance of the number difference at the output of the
interferometer are made, the phase sensitivity is given by

$\displaystyle\Delta\phi$ $\displaystyle=\frac{\xi}{\sqrt{N}}\,,$ (6)

where $\xi$ is the Wineland spin-squeezing parameter Wineland  _et al._ (1992)
defined by

$\displaystyle\xi$
$\displaystyle=\sqrt{\frac{N\mathrm{Var}(\hat{J}_{\theta})}{\langle\hat{J}_{x}\rangle^{2}}}\,,$
(7)

with

$\displaystyle\hat{J}_{\theta}$
$\displaystyle=\hat{J}_{z}\cos\theta+\hat{J}_{y}\sin\theta\,,$ (8)

where the angle $\theta$ is chosen to minimise
$\mathrm{Var}(\hat{J}_{\theta})$. This state is then converted into a state
with reduced fluctuations in number difference by a rotation by $-\theta$
around the $\hat{J}_{x}$, before entering the interferometer. In practice,
this is done via an additional coherent Rabi pulse.

##  III Ideal single-mode dynamics

We begin by revising the ideal _single-mode_ dynamics of both the OAT and TNT
Hamiltonians. We do this by expanding the field operator into an orthonormal
single-particle basis, and then make the assumption that only one mode is
occupied, ie

$\displaystyle\hat{\psi}_{a}(\mathbf{r})$
$\displaystyle=\sum_{j=0}^{\infty}\hat{a}_{j}u_{a,j}(\mathbf{r})\approx
u_{a,0}(\mathbf{r})\hat{a}_{0}$ (9a) $\displaystyle\hat{\psi}_{b}(\mathbf{r})$
$\displaystyle=\sum_{j=0}^{\infty}\hat{b}_{j}u_{b,j}(\mathbf{r})\approx
u_{b,0}(\mathbf{r})\hat{b}_{0}\,.$ (9b)

This approximation has been shown to be reasonably valid in tight confining
potentials with small numbers of atoms Gross  _et al._ (2010); Riedel  _et
al._ (2010), but breaks down when the population increases to the point where
the energy due to atomic interactions dominate Haine  _et al._ (2014). Making
this approximation in Eq. (2) gives

$\displaystyle\hat{H}$
$\displaystyle=\hbar\chi_{aa}\hat{a}^{\dagger}\hat{a}^{\dagger}\hat{a}\hat{a}+\hbar\chi_{bb}\hat{b}^{\dagger}\hat{b}^{\dagger}\hat{b}\hat{b}+2\hbar\chi_{ab}\hat{a}^{\dagger}\hat{a}\hat{b}^{\dagger}\hat{b}$
$\displaystyle+\frac{\hbar\Omega}{2}\left(\eta\hat{a}^{\dagger}\hat{b}+\eta^{*}\hat{b}^{\dagger}\hat{a}\right)\,,$
(10)

where

$\displaystyle\chi_{ij}=\frac{U_{ij}}{2\hbar}\int\left|u_{i}(\mathbf{r})\right|^{2}\left|u_{j}(\mathbf{r})\right|^{2}\,d^{3}\mathbf{r}\,,$
(11)

and

$\displaystyle\eta=\int
u_{a}^{*}(\mathbf{r})u_{b}(\mathbf{r})\,d^{3}\mathbf{r}\,,$ (12)

and we have made the replacements $\hat{a}_{0}\rightarrow\hat{a}$,
$\hat{b}_{0}\rightarrow\hat{b}$, $u_{a(b),0}(\mathbf{r})\rightarrow
u_{a(b)}(\mathbf{r})$ for notational simplicity. Introducing the single-mode
version of the pseudo-spin operators (Eq. (5))

$\displaystyle\hat{J}_{x}$
$\displaystyle=\frac{1}{2}\left(\hat{a}^{\dagger}\hat{b}+\hat{b}^{\dagger}\hat{a}\right)\,,$
(13a) $\displaystyle\hat{J}_{y}$
$\displaystyle=\frac{i}{2}\left(\hat{a}^{\dagger}\hat{b}-\hat{b}^{\dagger}\hat{a}\right)\,,$
(13b) $\displaystyle\hat{J}_{z}$
$\displaystyle=\frac{1}{2}\left(\hat{a}^{\dagger}\hat{a}-\hat{b}^{\dagger}\hat{b}\right)\,,$
(13c)

and assuming perfect overlap ($\eta=1$), the Hamiltonian (Eq. (10)) can be
written as

$\displaystyle\hat{H}$
$\displaystyle=\hbar\chi\hat{J}_{z}^{2}+\hbar\chi_{-}(\hat{N}-1)\hat{J}_{z}+\hbar\Omega\hat{J}_{x}\,,$
(14)

where $\chi=\chi_{aa}+\chi_{bb}-2\chi_{ab}$, $\chi_{-}=\chi_{aa}-\chi_{bb}$,
and $\hat{N}=\hat{a}^{\dagger}\hat{a}+\hat{b}^{\dagger}\hat{b}$ is the total
number of atoms. For a fixed number of atoms $N$, when $\chi_{-}=0$, and
$\Omega=\chi N/2$, this is the well studied TNT Hamiltonian, and when
$\Omega=0$ we recover the OAT Hamiltonian. Specifically, we define

$\displaystyle\hat{H}_{\mathrm{OAT}}$
$\displaystyle=\hbar\chi\hat{J}_{z}^{2}\,,$ (15a)
$\displaystyle\hat{H}_{\mathrm{TNT}}$
$\displaystyle=\hbar\chi\hat{J}_{z}^{2}+\Omega\hat{J}_{x}\,.$ (15b)

In writing Eq. (15a) and Eq. (15b), we have neglected terms that only depend
on $\hat{N}$, as they do not result in any observable effects.

Figures 1 and 2 show the evolution of an initial coherent spin state (CSS)
under these Hamiltonians. We visualise the state using the Husimi $Q$-function
Arecchi  _et al._ (1972); Agarwal (1998). The most notable difference between
OAT and TNT dynamics, is the speed at which the state evolves away from a CSS.
Figure (3) shows the variance of the three pseuedo-spins ($\hat{J}_{x}$,
$\hat{J}_{y}$, $\hat{J}_{z}$) for both cases. We see that the timescale for
the variance to increase from $\sim N/4$ (that of a CSS) to $\sim N^{2}/8$
(that of a highly non-classical state, such as a twin Fock-state, for example)
is approximately 40 times faster for the TNT dynamics.

Figure 1: The Husimi $Q$-function $Q(\theta,\varphi)$, for the state evolving
under Eq. (15a) for $N=100$. The $Q$-function is defined as
$Q(\theta,\varphi)=\left|\langle\theta,\varphi|\Psi(t)\rangle\right|^{2}$,
where
$|\theta,\varphi\rangle=e^{i\varphi\hat{J}_{z}}e^{i\theta\hat{J}_{y}}|J_{z}=N/2\rangle$
represents the spin coherent state along $\theta$ and $\varphi$ directions
corresponding to rotating the maximal $\hat{J}_{z}$ eigenstate around
azimuthal and polar angles $\\{\theta,\varphi\\}$ Arecchi  _et al._ (1972).
(a): $\chi t=0$ and (b): $\chi t=0.05$, (c): $\chi t=0.1$, (d): $\chi t=0.3$
Figure 2: The Husimi $Q$-function $Q(\theta,\varphi)$, for the state evolving
under Eq. (15b) for $N=100$. The $Q$-function is defined as
$Q(\theta,\varphi)=\left|\langle\theta,\varphi|\Psi(t)\rangle\right|^{2}$,
where
$|\theta,\varphi\rangle=e^{i\varphi\hat{J}_{z}}e^{i\theta\hat{J}_{y}}|J_{z}=N/2\rangle$
represents the spin coherent state along $\theta$ and $\varphi$ directions
corresponding to rotating the maximal $\hat{J}_{z}$ eigenstate around
azimuthal and polar angles $\\{\theta,\varphi\\}$ Arecchi  _et al._ (1972).
(a): $\chi t=0$ and (b): $\chi t=0.02$, (c): $\chi t=0.04$, (d): $\chi t=0.06$

Two quantities of significant interest for sensing applications are the spin-
squeezing parameter $\xi$, and the quantum Fisher information. Figure (4a)
shows the spin-squeezing parameter for both cases. We see that in the case of
TNT, $\xi$ decreases significantly faster than for OAT, and OAT reaches a
lower minimum. Specifically, it takes OAT more than twice as long to reach the
minimum $\xi$ from TNT dynamics. Alternatively, if the state preparation time
is constrained, TNT achieves a value of $\xi$ 2.3 times smaller for the same
state preparation time. Beyond this minimum, ‘ _oversqueezing_ ’ prevents any
further reduction in $\xi$ despite the presence of metrologically useful
entanglement Tóth (2012). When measurements are made that can resolve the full
probability distribution rather than just the mean and variance of the
collective spin, as is required for resolving the full metrological potential
of highly entangled quantum states Haine (2018b), a more relevant metric for
quantifying the sensitivity is the quantum Fisher information (QFI), which
relates to the sensitivity via the quantum Cramer-Rao bound
$\Delta\phi=1/\sqrt{F_{Q}}$. For an initial pure state entering a MZ
interferometer, the QFI is $F_{Q}=4\mathrm{Var}(\hat{J}_{y})$. However, if we
allow for an additional $\hat{J}_{x}$ rotation before the MZ interferometer,
the QFI can be increased further. Specifically, we use the definition

$F_{Q}=4\mathrm{V}(\hat{J}_{\theta})$ (16)

where $\hat{J}_{\theta}=\cos\theta\hat{J}_{y}+\sin\theta\hat{J}_{z}$, and
$\theta$ is chosen to maximise the QFI. We note that a full probability-
resolving measurement is not required if an interaction-based readout is used
after the interferometer Nolan  _et al._ (2017); Mirkhalaf  _et al._ (2018);
Haine (2018b). Figure (4b) shows the QFI for both OAT and TNT for $N=10^{5}$
particles. We see that TNT reaches a higher value, and achieves this maximum
significantly more quickly. Specifically, TNT reaches the threshold
$F_{Q}\approx N^{2}/2$, which is the QFI of the highly non-classical twin-Fock
state, $\sim 40$ times faster than OAT. Alternatively, at the time when TNT
reaches this threshold, the QFI is $\sim 337$ times larger than for OAT.

Figure 3: Variance of the three components of the pseudo-spin for (a) OAT and
(b) TNT dynamics. The black dotted line represents the threshold $N^{2}/8$,
which is the approximate value that OAT dynamics plateaus at for long time.
The number of atoms was $N=10^{5}$.  Figure 4: (a) Evolution of the spin
squeezing parameter for OAT (blue solid line) and TNT dynamics (red dashed
line). (b) The QFI for OAT (blue solid line) and TNT (red dashed line). The
three horizontal black dotted lines represent $F_{Q}=N$ (the limit for
unentangled particles), $N^{2}/2$, (the long-time plateu for OAT dynamics) and
$N^{2}$ (the Heisenberg limit). The number of particles was $N=10^{5}$.

##  IV Multi-Mode Dynamics

We now investigate the effects of multi-mode dynamics. Assuming a cigar shaped
trapping potential where the two transverse dimensions ($y$ & $z$) are much
tighter than the longitudinal dimension ($x$), the dynamics in the $x$
direction is well approximated by the Heisenberg equations of motion

$\displaystyle i\hbar\frac{d}{dt}\hat{\psi}_{a}$
$\displaystyle=\left(\hat{H}_{1D}+\tilde{U}_{aa}\hat{\psi}^{\dagger}_{a}\hat{\psi}_{a}+\tilde{U}_{ab}\hat{\psi}^{\dagger}_{b}\hat{\psi}_{b}\right)\hat{\psi}_{a}+\frac{\hbar\Omega}{2}\hat{\psi}_{b}$
(17a) $\displaystyle i\hbar\frac{d}{dt}\hat{\psi}_{b}$
$\displaystyle=\left(\hat{H}_{1D}+\tilde{U}_{ab}\hat{\psi}^{\dagger}_{a}\hat{\psi}_{a}+\tilde{U}_{bb}\hat{\psi}^{\dagger}_{b}\hat{\psi}_{b}\right)\hat{\psi}_{b}+\frac{\hbar\Omega}{2}\hat{\psi}_{a}$
(17b)

where

$\displaystyle\hat{H}_{1D}$
$\displaystyle=\frac{-\hbar^{2}}{2m}\frac{\partial^{2}}{\partial
x^{2}}+\frac{1}{2}m\omega_{x}^{2}x^{2}\,,$ (18)

and $\tilde{U}_{ij}=U_{ij}/A_{\perp}$ is the dimensionally reduced effective
one dimensional interaction strength obtained by dividing the three
dimensional interaction strength by a parameter characterising the transverse
area of the system. Throughout this work, we use $m=87$ amu, and
$A_{\perp}=10^{-10}$ m2.

We model the situation where initially all the atoms in are in state
$|a\rangle$ in the ground motional state, and then coherently transfer half
the population to state $|b\rangle$. Unless $U_{aa}=U_{bb}=U_{ab}$, this new
state will not be the ground state, and motional dynamics will occur. As
$\chi\propto U_{aa}+U_{bb}-2U_{ab}\neq 0$, we cannot obtain entangling
dynamics without also exciting motional dynamics. We consider three distinct
cases that provide qualitatively different dynamics:

  * •

Case I: $U_{aa}=U_{bb}>U_{ab}$. In this case, the two components will undergo
breathing oscillations, but will tend to breathe-together. Specifically, we
choose $a_{aa}=a_{bb}=100.0a_{0}$, $a_{ab}=97.0a_{0}$, where $a_{0}=5.29\times
10^{-11}$ m.

  * •

Case II: $a_{bb}>a_{aa}>a_{ab}$. In this case, the components tend to
separate, as one component breathes inwards while the other breathes outwards,
such that the overlap of the two components is significantly decreased.
However, by adjusting the relative atom numbers, a _breathe-together_ solution
exists Sinatra, A. and Castin, Y. (2000). Specifically, we chose
$a_{aa}=95.0a_{0}$, $a_{bb}=100.0a_{0}$, $a_{ab}=90.0a_{0}$.

  * •

Case III: $U_{aa}>U_{ab}>U_{bb}$. The two components will tend to separate,
and no breathe-together solution exists. Specifically, we choose
$a_{aa}=100.0a_{0}$, $a_{bb}=95.0a_{0}$, $a_{ab}=97.0a_{0}$. This represents
the scattering parameters of the two hyperfine ground-states of 87Rb.

To illustrate the three cases, we calculate the density distribution under the
mean-field approximation by solving the Gross-Pitaevskii equation Dalfovo  _et
al._ (1999), obtained by making the substitution
$\hat{\psi}_{j}(x)\rightarrow\psi_{j}(x)$ in equations (17a,17b). Figures (5,
6, 7) shows the evolution of the density distribution for these three cases.
In particular, note that in cases I the two components evolve together, while
in case II and III, the two components spatially separate. This separation
will hinder the ability to implement TNT dynamics, as the varying spatial
overlap will complicate the coherent coupling required for the $\hat{J}_{x}$
rotation. Additionally, the spin-squeezing parameter requires a high degree of
overlap between the modes.

Figure 5: Evolution of the density of each component for Case I.
$|\psi_{a}(x,t)|^{2}$ (blue solid line), $|\psi_{b}(x,t)|^{2}$ (red dashed
line), compared to the initial state
$|\psi_{a}(x,0)|^{2}=|\psi_{b}(x,0)|^{2}=|\psi_{0}(x)|^{2}$ (black dotted
line). Both components vary only slightly from the initial condition, but
remain identical to each other.  Figure 6: Evolution of the density of each
component for Case II when a 50/50 beamsplitter is implemented.
$|\psi_{a}(x,t)|^{2}$ (blue solid line), $|\psi_{b}(x,t)|^{2}$ (red dashed
line), compared to the initial state
$|\psi_{a}(x,0)|^{2}=|\psi_{b}(x,0)|^{2}=|\psi_{0}(x)|^{2}$ (black dotted
line). Figure 7: Evolution of the density of each component for Case III.
$|\psi_{a}(x,t)|^{2}$ (blue solid line), $|\psi_{b}(x,t)|^{2}$ (red dashed
line), compared to the initial state
$|\psi_{a}(x,0)|^{2}=|\psi_{b}(x,0)|^{2}=|\psi_{0}(x)|^{2}$ (black dotted
line).

The Gross-Pitaeveskii equation is incapable of capturing the evolution of the
quantum statistics, which is required to investigate spin-squeezing and
entanglement. To investigate this effect, we simulate the dynamics of the
system using the Truncated Wigner (TW) method, which has previously been used
to model the dynamics of quantum gases Steel  _et al._ (1998); Sinatra  _et
al._ (1995); Norrie  _et al._ (2006); Drummond and Opanchuk (2017), and unlike
the GPE, can be used to model non-classical particle correlations Haine and
Ferris (2011); Ruostekoski and Martin (2013); Haine  _et al._ (2014); Szigeti
_et al._ (2017); Haine (2018b); Szigeti  _et al._ (2020). The derivation of
the TW method has been described in detail elsewhere Drummond and Hardman
(1993); Steel  _et al._ (1998); Blakie  _et al._ (2008). Briefly, the equation
of motion for the Wigner function of the system can be found from the von-
Neumann equation by using correspondences between differential operators on
the Wigner function and the original quantum operators Gardiner and Zoller
(2004). By truncating third- and higher-order derivatives (the TW
approximation), a Fokker–Planck equation (FPE) is obtained. The FPE is then
mapped to a set of stochastic partial differential equations for complex
fields $\psi_{j}(x,t)$, which loosely correspond to the original field
operators $\hat{\psi}_{j}(x,t)$, with initial conditions stochastically
sampled from the appropriate Wigner distribution Blakie  _et al._ (2008);
Olsen and Bradley (2009). The complex fields obey the partial differential
equations

$\displaystyle i\hbar\frac{d}{dt}\psi_{a}$
$\displaystyle=\hat{H}_{1D}\psi_{a}+\frac{\hbar\Omega}{2}\psi_{b}$
$\displaystyle+\left(\tilde{U}_{aa}\left(|\psi_{a}|^{2}-\frac{1}{\Delta
x}\right)+\tilde{U}_{ab}\left(|\psi_{b}|^{2}-\frac{1}{\Delta
x}\right)\right)\psi_{a}$ (19a) $\displaystyle i\hbar\frac{d}{dt}\psi_{b}$
$\displaystyle=\hat{H}_{1D}\psi_{b}+\frac{\hbar\Omega}{2}\psi_{a}$
$\displaystyle+\left(\tilde{U}_{ab}\left(|\psi_{a}|^{2}-\frac{1}{\Delta
x}\right)+\tilde{U}_{bb}\left(|\psi_{b}|^{2}-\frac{1}{\Delta
x}\right)\right)\psi_{b}$ (19b)

where $\Delta$ is discretisation size of the spatial $x$ grid. By averaging
over many trajectories with stochastically sampled initial conditions,
expectation values of quantities corresponding to symmetrically ordered
operators in the full quantum theory can be obtained via the correspondence
$\langle\\{f(\hat{\psi}^{\dagger}_{j},\hat{\psi}_{j})\\}_{\mathrm{sym}}\rangle=\overline{f[\psi_{j}^{*},\psi_{j}]}$,
where ‘sym’ denotes symmetric ordering and the overline denotes the mean over
many stochastic trajectories. The initial conditions for the simulations are
chosen as $\psi_{a}(x,0)=\Psi_{0}(x)+\eta_{a}(x)$,
$\psi_{b}(x,0)=\eta_{b}(x)$, where $\Psi_{0}(x)$ is the ground state of the
single-component time-independent Gross-Pitaevski equation, and
$\eta_{j}(\xi)$ are complex Gaussian noises satisfying
$\overline{\eta^{*}_{i}(x_{n})\eta_{j}(x_{m})}=\frac{1}{2}\delta_{m,n}\delta_{i,j}/\Delta$,
for spatial grid points $x_{m}$ and $x_{n}$. At $t=0$, the $\pi/2$ beam
splitting pulse which initiates the dynamics is implemented via the
transformation
$\psi_{a}(x)=\frac{1}{\sqrt{2}}\left(\psi_{a}(x,0)+\psi_{b}(x,0)\right)$,
$\psi_{a}(x)=\frac{1}{\sqrt{2}}\left(\psi_{b}(x,0)-\psi_{a}(x,0)\right)$.

###  IV.1 Case: I

In order to implement TNT dynamics, the rotation parameter $\Omega$ must be
set to the optimum value $\Omega=\chi N/2$. While this is simple to do in the
single-mode model, in the multimode model the effective $\chi$ depends on the
time-dependent density distributions. Furthermore, it has previously been
shown that when strong multimode dynamics are present, estimates of $\chi$
based on Eq. (11) are poor Haine  _et al._ (2014). We estimate the effective
$\chi$ by first setting $\Omega=0$ and calculating the variance of the three
pseudo-spin operators, and comparing to ideal single-mode OAT dynamics for the
same number of atoms. The single-mode dynamics resulting from Eq. (15a) is
also calculated via the TW method 111The quantum dynamics for the single-mode
model is obtained using the TW method by solving the ODEs
$i\dot{\alpha}=\frac{1}{2}\chi(|\alpha|^{2}-|\beta|^{2})\alpha+\Omega\beta$,
$i\dot{\beta}=-\frac{1}{2}\chi(|\alpha|^{2}-|\beta|^{2})\beta+\Omega\alpha$,
where $\alpha$ and $\beta$ are the stochastic variables corresponding to
$\hat{a}_{0}$ and $\hat{b}_{0}$ respectively.. Figure (8) shows a comparison
between the multi-mode dynamics and single-mode dynamics, with the parameter
$\chi$ in the single-mode mode adjusted to provide the best match to the
multimode dynamics. We use this as our estimate of $\chi$ when choosing a
value of $\Omega$ when implementing TNT dynamics in the multi-mode model.

Figure 8: Comparison between the multi-mode dynamics and ideal single-mode
behaviour for Case I, when only OAT dynamics is implemented. When we choose
$\chi=2.2\times 10^{-3}$ s-1, we have excellent agreement in the variance of
the pseudospins. $V(J_{y})$ as calculated from the single-mode model is
indicated via the red dotted line. Figure 9: Comparison between multi-mode OAT
dynamics and ideal OAT single-mode behaviour for Case I, for the spin-
squeezing parameter (a) and QFI (b). Figure 10: Comparison between the multi-
mode dynamics and ideal single-mode behavior when TNT dynamics is implemented.
The value of $\Omega$ is chosen as $\Omega=\chi N/2$, where the value
$\chi=2.2\times 10^{-3}$ s-1, was used. $V(J_{z})$ as calculated from the
single-mode model is indicated via the red dotted line.

Figure 11: Comparison between multi-mode TNT dynamics and ideal TNT single-
mode behaviour for Case I, for the spin-squeezing parameter (a) and QFI (b)

.

We see here that setting $\chi=2.2\times 10^{-3}$ s-1 provides excellent
agreement between the single-mode model and the multimode model for both the
pseudospin variances (fig. 8) spin-squeezing parameter, and QFI (fig.9). Using
this value in our choice of $\Omega=N\chi/2$ should therefore result in TNT
dynamics. Figure 10 shows the multimode dynamics compared to the single mode
dynamics. We see that there is good agreement between the single-mode and
multi-mode models, which is unsurprising given that the two components evolve
identically, ensuring that the overlap between the two components remains
constant, and that the relative phase is the same at each point in space. This
ensures the coupling term results in the pure $J_{x}$ rotation required for
TNT dynamics. In this parameter regime, the spin squeezing parameter and QFI
also display excellent agreement with the ideal single-mode dynamics (figure
11). Importantly, we see that development of entanglement (characterised by
large quantum Fisher information) occurs much more rapidly than with OAT
dynamics.

###  IV.2 Case II

When $U_{bb}\neq U_{aa}$, the two components will evolve differently,
effecting both the spatial overlap, and the relative phase. As before, we
first simulate the system with $\Omega=0$ to investigate the agreement between
the single-mode and multi-mode systems for OAT dynamics. However, as can be
seen in figure (6), the two components begin to separate, which will inhibit
the performance of the TNT dynamics. In order to prevent this, we exploit the
_breathe-together_ solution Li  _et al._ (2008). This is achieved by replacing
the initial 50/50 beam-splitter with an asymmetric beamsplitter such that each
component experiences the same interaction strength. Specifically, this is
achieved by choosing the beam-splitter angle such that the ratio of population
in each component, $N_{a}$ and $N_{b}$, satisfy

$\frac{N_{a}}{N_{b}}=\frac{U_{bb}-U_{ab}}{U_{aa}-U_{ab}}\,.$ (20)

As can be seen from figure 12, this choice of initial condition results in the
motional dynamics being frozen out. The dynamics of the quantum statistics,
however, is shown in figure 13, as well as the dynamics for the equivalent
single-mode system. The different initial conditions result in reduced final
spin variances for both the multi-mode and ideal single-mode dynamics. By
comparing both cases, we can infer that the effective interaction parameter
$\chi$ is $\chi=5.25\times 10^{-3}$ Hz. Again, we use this value to choose the
appropriate value of $\Omega$ for implementation of TNT dynamics. However, as
the breathe-together solution is asymmetric in population, the optimum value
of $\Omega$ will be slightly different. We account for this by trialling a
range of values, as shown in figure 15. We find that the most effective
entangling dynamics (that is, the dynamics that leads to the largest spin-
variance) is when $\Omega=0.85\chi N/2$.

Figure 12: Evolution of the density of each component for Case II when an
asymmetric beamspitter is implemented, in order to satisfy the breathe-
together solution. $|\psi_{a}(x,t)|^{2}$ (blue solid line),
$|\psi_{b}(x,t)|^{2}$ (red dashed line), compared to the initial conditions
$|\psi_{a}(x,0)|^{2}$ and $|\psi_{b}(x,0)|^{2}$ (black dotted lines). Both
components deviate only slightly from their initial conditions.  Figure 13:
Comparison between multi-mode OAT dynamics and ideal OAT single-mode behaviour
for Case II with the breathe-together initial condition. For the single-mode
simulation, $\chi=5.25\times 10^{-3}$ Hz was found to have the best agreement
with the multimode results. $V(J_{y})$ as calculated from the single-mode
model is indicated via the red dotted line.  Figure 14: Comparison between
multi-mode OAT dynamics and ideal OAT single-mode behaviour for Case II with
the breathe-together solution, for the spin-squeezing parameter (a) and QFI
(b). Figure 15: Comparison between single-mode TNT dynamics (a) and multimode
TNT dynamics ((b) and (c)). The initial condition was chosen to satisfy the
breathe-together solution. In (a) and (b), the optimal rotation rate
($\Omega=0.85\chi N/2$) was used, which gives slightly better performance than
the usual TNT solution ($\Omega=\chi N/2$), shown in (c).

In figure 16 we see good agreement between the ideal single-mode and multi-
mode TNT dynamics for the spin-squeezing parameter and QFI. Importantly, the
entangling dynamics occurs $\sim 40$ times faster than for OAT alone.

Figure 16: Comparison between multi-mode TNT dynamics and ideal TNT single-
mode behaviour for Case II with the breathe-together solution, for the spin-
squeezing parameter (a) and QFI (b). $\Omega=0.85\chi N/2$ was used for both
cases.

###  IV.3 Case III

When $U_{aa}>U_{ab}>U_{bb}$, there is no breathe-together solution. As before,
we first simulate the system with $\Omega=0$ to investigate the behaviour in a
multimode system for OAT dynamics. However, as can be seen in figure (7), the
two components begin to separate, which will inhibit the performance of the
TNT dynamics. Additionally, as can be seen in figure 17, the asymmetry in
scattering lengths results in an additional rotation of the collective spin
around the $J_{z}$ axis, which will further inhibit TNT dynamics, as we
require a coherent rotation around the collective spin direction. We can
correct for this term by adding an additional rotation of angular frequency
$\omega_{r}$, either by adjusting the detuning between the two levels, or by
dynamically rotating the relevant rotation axis for our TNT dynamics. Figure
18 shows the spin dynamics of the system under OAT dynamics with this
additional correction. While not mimicking the single mode OAT dynamics
perfectly, it displays qualitatively similar behaviour.

Figure 17: Multimode OAT dynamics behaviour for Case III when 50/50
beamspliter is implemented. (a) the expectation values of the spin operators,
and (b) the variances of the spin operators Figure 18: Multimode OAT dynamics
behaviour for Case III when a rotation about the $\hat{J}_{z}$ axis is added
to compensate for the spin-precession. (a) the expectation values of the spin
operators, and (b) the variances of the spin operators.

We have modelled TNT dynamics with this additional rotation term (figure 19).
We see that, while the timescale of the entangling dynamics is significantly
faster than OAT alone, it does not behave as well as ideal TNT dynamics. The
reason for this is that while the additional rotation partially corrects for
the drifting phase difference between the two components, the multimode
dynamics ensures that there is a dynamic and spatially-varying phase
difference, so this cancellation is imperfect. This effect has an even more
pronounced effect on the spin-squeezing and QFI (figure 21); there is little
improvement in either the rate or depth of spin-squeezing achievable by
implementing TNT dynamics.

Figure 19: Multimode TNT dynamics behaviour for Case III, for the spin
variance (a) and three mean value of angular operators (b). Figure 20:
Multimode OAT dynamics behaviour for Case III, for the spin squeeze parameter
(a) and QFI (b). Figure 21: Multimode TNT dynamics behaviour for Case III, for
the spin squeeze parameter (a) and QFI (b).

##  V Discussion

Our modelling indicates that in case I and case II systems, TNT can be used to
significantly speed up spin-squeezing and entangling dynamics even in regimes
where significant multi-mode dynamics is present. Importantly, in case II
systems, use of the breathe-together solution can be used to ensure strong
mode-overlap between the two components. In both of these cases, the spin-
squeezing dynamics of the full system is well approximated by an effective
two-mode model.

However, in case III systems, the implementation of TNT dynamics provides
little-to-no benefit over conventional OAT dynamics. We attribute this to a
spatially varying phase-profile of the two components, which results in a
variation in the rotation axis of the effective TNT rotation.

Our results indicate that in some regimes, spin-squeezing and entanglement
generation via atomic self-interaction is achievable in BECs with large
numbers of atoms, even in the presence of multi-mode dynamics, and that TNT
dynamics can be used to decrease the state-preparation time required in order
to achieve this squeezing. Alternatively, faster entangling dynamics can
result in better spin-squeezing in cases where the state-preparation time is
limited Haine and Hope (2020); Szigeti  _et al._ (2020, 2021).

##  VI ACKNOWLEDGMENTS

This research was undertaken with the assistance of re-sources and services
from the National Computational Infrastructure (NCI), which is supported by
the Australian Government. The Australian National University is situated on
land traditionally owned by the Ngunnawal people.

## References

  * Pezzè  _et al._ (2018) Luca Pezzè, Augusto Smerzi, Markus K. Oberthaler, Roman Schmied, and Philipp Treutlein, “Quantum metrology with nonclassical states of atomic ensembles,” Rev. Mod. Phys. 90, 035005 (2018).
  * Szigeti  _et al._ (2021) Stuart S. Szigeti, Onur Hosten, and Simon A. Haine, “Improving cold-atom sensors with quantum entanglement: Prospects and challenges,” Applied Physics Letters 118, 140501 (2021), https://doi.org/10.1063/5.0050235 .
  * Giovannetti  _et al._ (2006) Vittorio Giovannetti, Seth Lloyd, and Lorenzo Maccone, “Quantum metrology,” Phys. Rev. Lett. 96, 010401 (2006).
  * Pezzé and Smerzi (2009) Luca Pezzé and Augusto Smerzi, “Entanglement, nonlinear dynamics, and the Heisenberg limit,” Phys. Rev. Lett. 102, 100401 (2009).
  * Kitagawa and Ueda (1993) Masahiro Kitagawa and Masahito Ueda, “Squeezed spin states,” Phys. Rev. A 47, 5138–5143 (1993).
  * Mølmer and Sørensen (1999) Klaus Mølmer and Anders Sørensen, “Multiparticle entanglement of hot trapped ions,” Phys. Rev. Lett. 82, 1835–1838 (1999).
  * Esteve  _et al._ (2008) J. Esteve, C. Gross, A. Weller, S. Giovanazzi, and M. K. Oberthaler, “Squeezing and entanglement in a Bose-Einstein condensate,” Nature 455, 1216 (2008).
  * Leroux  _et al._ (2010) Ian D. Leroux, Monika H. Schleier-Smith, and Vladan Vuletić, “Implementation of cavity squeezing of a collective atomic spin,” Phys. Rev. Lett. 104, 073602 (2010).
  * Gross  _et al._ (2010) C. Gross, T. Zibold, E. Nicklas, J. Esteve, and M. K. Oberthaler, “Nonlinear atom interferometer surpasses classical precision limit,” Nature 464, 1165–1169 (2010).
  * Riedel  _et al._ (2010) Max F. Riedel, Pascal Böhi, Yun Li, Theodor W. Hänsch, Alice Sinatra, and Philipp Treutlein, “Atom-chip-based generation of entanglement for quantum metrology,” Nature 464, 1170–1173 (2010).
  * Muessel  _et al._ (2014) W. Muessel, H. Strobel, D. Linnemann, D. B. Hume, and M. K. Oberthaler, “Scalable spin squeezing for quantum-enhanced magnetometry with bose-einstein condensates,” Phys. Rev. Lett. 113, 103004 (2014).
  * Li  _et al._ (2008) Yun Li, Y. Castin, and A. Sinatra, “Optimum spin squeezing in bose-einstein condensates with particle losses,” Phys. Rev. Lett. 100, 210401 (2008).
  * Li  _et al._ (2009) Yun Li, P. Treutlein, J. Reichel, and A. Sinatra, “Spin squeezing in a bimodal condensate: spatial dynamics and particle losses,” The European Physical Journal B 68, 365–381 (2009).
  * Haine (2018a) Simon A Haine, “Quantum noise in bright soliton matterwave interferometry,” New Journal of Physics 20, 033009 (2018a).
  * Haine and Johnsson (2009) Simon A. Haine and Mattias T. Johnsson, “Dynamic scheme for generating number squeezing in Bose-Einstein condensates through nonlinear interactions,” Phys. Rev. A 80, 023611 (2009).
  * Haine  _et al._ (2014) S. A. Haine, J. Lau, R. P. Anderson, and M. T. Johnsson, “Self-induced spatial dynamics to enhance spin squeezing via one-axis twisting in a two-component Bose-Einstein condensate,” Phys. Rev. A 90, 023613 (2014).
  * Szigeti  _et al._ (2020) Stuart S. Szigeti, Samuel P. Nolan, John D. Close, and Simon A. Haine, “High-precision quantum-enhanced gravimetry with a bose-einstein condensate,” Phys. Rev. Lett. 125, 100402 (2020).
  * Micheli  _et al._ (2003) A. Micheli, D. Jaksch, J. I. Cirac, and P. Zoller, “Many-particle entanglement in two-component bose-einstein condensates,” Phys. Rev. A 67, 013607 (2003).
  * Sorelli  _et al._ (2019) Giacomo Sorelli, Manuel Gessner, Augusto Smerzi, and Luca Pezzè, “Fast and optimal generation of entanglement in bosonic josephson junctions,” Phys. Rev. A 99, 022329 (2019).
  * Mirkhalaf  _et al._ (2018) Safoura S. Mirkhalaf, Samuel P. Nolan, and Simon A. Haine, “Robustifying twist-and-turn entanglement with interaction-based readout,” Phys. Rev. A 97, 053618 (2018).
  * Strobel  _et al._ (2014) Helmut Strobel, Wolfgang Muessel, Daniel Linnemann, Tilman Zibold, David B. Hume, Luca Pezzè, Augusto Smerzi, and Markus K. Oberthaler, “Fisher information and entanglement of non-gaussian spin states,” Science 345, 424–427 (2014).
  * Muessel  _et al._ (2015) W. Muessel, H. Strobel, D. Linnemann, T. Zibold, B. Julia-Diaz, and M. K. Oberthaler, “Twist-and-turn spin squeezing in Bose-Einstein condensates,” arXiv:1507.02930 (2015).
  * Lipkin  _et al._ (1965) H.J. Lipkin, N. Meshkov, and A.J. Glick, “Validity of many-body approximation methods for a solvable model: (i). exact solutions and perturbation theory,” Nuclear Physics 62, 188 – 198 (1965).
  * Demkowicz-Dobrzanski  _et al._ (2015) R. Demkowicz-Dobrzanski, M. Jarzyna, and J. Kolodynski, “Quantum limits in optical interferometry,” Progress in Optics 60, 345 (2015).
  * Braunstein and Caves (1994) Samuel L. Braunstein and Carlton M. Caves, “Statistical distance and the geometry of quantum states,” Phys. Rev. Lett. 72, 3439–3443 (1994).
  * Paris (2009) Matteo G. A. Paris, “Quantum estimation for quantum technology,” International Journal of Quantum Information 07, 125–137 (2009).
  * Tóth and Apellaniz (2014) Géza Tóth and Iagoba Apellaniz, “Quantum metrology from a quantum information science perspective,” Journal of Physics A: Mathematical and Theoretical 47, 424006 (2014).
  * Davis  _et al._ (2016) E. Davis, G. Bentsen, and M. Schleier-Smith, “Approaching the Heisenberg limit without single-particle detection,” Phys. Rev. Lett. 116, 053601 (2016).
  * Hosten  _et al._ (2016) O. Hosten, R. Krishnakumar, N. J. Engelsen, and M. A. Kasevich, “Quantum phase magnification,” Science 352, 1552–1555 (2016).
  * Fröwis  _et al._ (2016) Florian Fröwis, Pavel Sekatski, and Wolfgang Dür, “Detecting large quantum Fisher information with finite measurement precision,” Phys. Rev. Lett. 116, 090801 (2016).
  * Macrì  _et al._ (2016) Tommaso Macrì, Augusto Smerzi, and Luca Pezzè, “Loschmidt echo for quantum metrology,” Phys. Rev. A 94, 010102 (2016).
  * Linnemann  _et al._ (2016) D. Linnemann, H. Strobel, W. Muessel, J. Schulz, R. J. Lewis-Swan, K. V. Kheruntsyan, and M. K. Oberthaler, “Quantum-enhanced sensing based on time reversal of nonlinear dynamics,” Phys. Rev. Lett. 117, 013001 (2016).
  * Nolan  _et al._ (2017) Samuel P. Nolan, Stuart S. Szigeti, and Simon A. Haine, “Optimal and robust quantum metrology using interaction-based readouts,” Phys. Rev. Lett. 119, 193601 (2017).
  * Szigeti  _et al._ (2017) Stuart S. Szigeti, Robert J. Lewis-Swan, and Simon A. Haine, “Pumped-up SU(1,1) interferometry,” Phys. Rev. Lett. 118, 150401 (2017).
  * Anders  _et al._ (2018) Fabian Anders, Luca Pezzè, Augusto Smerzi, and Carsten Klempt, “Phase magnification by two-axis countertwisting for detection-noise robust interferometry,” Phys. Rev. A 97, 043813 (2018).
  * Haine (2018b) Simon A. Haine, “Using interaction-based readouts to approach the ultimate limit of detection-noise robustness for quantum-enhanced metrology in collective spin systems,” Phys. Rev. A 98, 030303(R) (2018b).
  * Hayes  _et al._ (2018) Anthony J Hayes, Shane Dooley, William J Munro, Kae Nemoto, and Jacob Dunningham, “Making the most of time in quantum metrology: concurrent state preparation and sensing,” Quantum Science and Technology 3, 035007 (2018).
  * Haine (2021) Simon A Haine, “Searching for signatures of quantum gravity in quantum gases,” New Journal of Physics 23, 033020 (2021).
  * Opanchuk  _et al._ (2012) B. Opanchuk, M. Egorov, S. Hoffmann, A. I. Sidorov, and P. D. Drummond, “Quantum noise in three-dimensional BEC interferometry,” EPL (Europhysics Letters) 97, 50003 (2012).
  * Nolan and Haine (2018) Samuel P. Nolan and Simon A. Haine, “Generating macroscopic superpositions with interacting bose-einstein condensates: Multimode speedups and speed limits,” Phys. Rev. A 98, 063606 (2018).
  * Chin  _et al._ (2010) Cheng Chin, Rudolf Grimm, Paul Julienne, and Eite Tiesinga, “Feshbach resonances in ultracold gases,” Rev. Mod. Phys. 82, 1225–1286 (2010).
  * Wineland  _et al._ (1992) D. J. Wineland, J. J. Bollinger, W. M. Itano, F. L. Moore, and D. J. Heinzen, “Spin squeezing and reduced quantum noise in spectroscopy,” Phys. Rev. A 46, R6797–R6800 (1992).
  * Arecchi  _et al._ (1972) F. T. Arecchi, Eric Courtens, Robert Gilmore, and Harry Thomas, “Atomic coherent states in quantum optics,” Phys. Rev. A 6, 2211–2237 (1972).
  * Agarwal (1998) G. S. Agarwal, “State reconstruction for a collection of two-level systems,” Phys. Rev. A 57, 671–673 (1998).
  * Tóth (2012) Géza Tóth, “Multipartite entanglement and high-precision metrology,” Phys. Rev. A 85, 022322 (2012).
  * Sinatra, A. and Castin, Y. (2000) Sinatra, A. and Castin, Y., “Binary mixtures of bose-einstein condensates: Phase dynamics and spatial dynamics,” Eur. Phys. J. D 8, 319–332 (2000).
  * Dalfovo  _et al._ (1999) Franco Dalfovo, Stefano Giorgini, Lev P. Pitaevskii, and Sandro Stringari, “Theory of Bose-Einstein condensation in trapped gases,” Rev. Mod. Phys. 71, 463–512 (1999).
  * Steel  _et al._ (1998) M. J. Steel, M. K. Olsen, L. I. Plimak, P. D. Drummond, S. M. Tan, M. J. Collett, D. F. Walls, and R. Graham, “Dynamical quantum noise in trapped Bose-Einstein condensates,” Phys. Rev. A 58, 4824–4835 (1998).
  * Sinatra  _et al._ (1995) A Sinatra, F Castelli, L A Lugiato, P Grangier, and J P Poizat, “Effective two-level model versus three-level model,” Quantum and Semiclassical Optics: Journal of the European Optical Society Part B 7, 405 (1995).
  * Norrie  _et al._ (2006) A. A. Norrie, R. J. Ballagh, and C. W. Gardiner, “Quantum turbulence and correlations in bose-einstein condensate collisions,” Phys. Rev. A 73, 043617 (2006).
  * Drummond and Opanchuk (2017) Peter D. Drummond and Bogdan Opanchuk, “Truncated wigner dynamics and conservation laws,” Phys. Rev. A 96, 043616 (2017).
  * Haine and Ferris (2011) S. A. Haine and A. J. Ferris, “Surpassing the standard quantum limit in an atom interferometer with four-mode entanglement produced from four-wave mixing,” Phys. Rev. A 84, 043624 (2011).
  * Ruostekoski and Martin (2013) J. Ruostekoski and A. D. Martin, “The truncated wigner method for bose gases,” in  _Quantum Gasses_ , edited by The editor (Imperial College Press, 2013) pp. 203–214.
  * Drummond and Hardman (1993) P. D. Drummond and A. D. Hardman, “Simulation of quantum effects in raman-active waveguides,” EPL (Europhysics Letters) 21, 279 (1993).
  * Blakie  _et al._ (2008) P. B. Blakie, A. S. Bradley, M. J. Davis, R. J. Ballagh, and C. W. Gardiner, “Dynamics and statistical mechanics of ultra-cold Bose gases using c-field techniques,”  _Advances in Physics_ , Advances in Physics 57, 363–455 (2008).
  * Gardiner and Zoller (2004) C. W. Gardiner and P. Zoller,  _Quantum Noise: A Handbook of Markovian and Non-Markovian Quantum Stochastic Methods with Applications to Quantum Optics_ , 3rd ed. (Springer, Berlin and Heidelberg, 2004).
  * Olsen and Bradley (2009) M.K. Olsen and A.S. Bradley, “Numerical representation of quantum states in the positive-P and wigner representations,” Optics Communications 282, 3924 – 3929 (2009).
  * Note (1) The quantum dynamics for the single-mode model is obtained using the TW method by solving the ODEs $i\mathaccentV{dot}05F{\alpha}=\frac{1}{2}\chi(|\alpha|^{2}-|\beta|^{2})\alpha+\Omega\beta$, $i\mathaccentV{dot}05F{\beta}=-\frac{1}{2}\chi(|\alpha|^{2}-|\beta|^{2})\beta+\Omega\alpha$, where $\alpha$ and $\beta$ are the stochastic variables corresponding to $\mathaccentV{hat}05E{a}_{0}$ and $\mathaccentV{hat}05E{b}_{0}$ respectively.
  * Haine and Hope (2020) Simon A. Haine and Joseph J. Hope, “Machine-designed sensor to make optimal use of entanglement-generating dynamics for quantum sensing,” Phys. Rev. Lett. 124, 060402 (2020).


# Resolution to Sutner’s Conjecture

William Boyles Department of Mathematics, North Carolina State University,
Raleigh, NC 27695<EMAIL_ADDRESS>
##  1 Introduction

Consider a game played on a simple graph $G=(V,E)$ where each vertex consists
of a clickable light. Clicking any vertex $v$ toggles on on/off state of $v$
and its neighbors. One wins the game by finding a sequence of clicks that
turns off all the lights. When $G$ is a $5\times 5$ grid, this game was
commercially available from Tiger Electronics as Lights Out.

Sutner was one of the first to study these games mathematically. He showed
that for any $G$ the initial configuration of all lights on is solvable [3].
He also found that when $d(G)=\text{dim}\left(\ker{(A+I)}\right)$ over the
field $\mathbb{Z}_{2}$, where $A$ is the adjacency matrix of $G$, is 0 all
initial configurations are solvable. In particular, 1 out of every $2^{d(G)}$
initial configurations are solvable, while each solvable configuration has
$2^{d(G)}$ distinct solutions [3]. When investigating $n\times n$ grid graphs,
Sutner conjectured the following relationship:

$\displaystyle d_{2n+1}$ $\displaystyle=2d_{n}+\delta_{n}\text{,
}\delta_{n}\in\\{0,2\\}$ $\displaystyle\delta_{2n+1}$
$\displaystyle=\delta_{n},$

where $d_{n}=d(G)$ for $G$ an $n\times n$ grid graph [3].

We resolve this conjecture in the affirmative. We use results from Sutner that
give the nullity of a $n\times n$ board as the GCD of two polynomials in the
ring $\mathbb{Z}_{2}[x]$ [4]. We then apply identities from Hunziker,
Machiavelo, and Park that relate the polynomials $(2n+1)\times(2n+1)$ grids
and $n\times n$ grids [2]. We then apply a result from Ore about the GCD of
two products [6]. Together, these results allow us to prove Sutner’s
conjecture. We then go further and show for exactly which values of $n$
$\delta_{n}$ is 0 or 2.

##  2 Fibonacci Polynomials

######  Definition 1.

Let $f_{n}(x)$ be the polynomial in the ring $\mathbb{Z}_{2}[x]$ defined
recursively by

$f_{n}(x)=\begin{cases}0&n=0\\\ 1&n=1\\\
xf_{n-1}(x)+f_{n-2}(x)&\text{otherwise }.\end{cases}$

These polynomials are often referred to as Fibonacci polynomials because when
defined over the reals, evaluating $f_{n}(x)$ at $x=1$ gives the $n$th
Fibonacci number. Since coefficients of $f_{n}(x)$ are either 1 or 0, one can
visualize them by coloring squares black or white to represent the
coefficients. For example,

$f(6,x)=1x^{5}+0x^{4}+0x^{3}+0x^{2}+1x+0=\blacksquare\square\square\square\square\blacksquare\square.$

Plotting $f(1,x)$ through $f(128,x)$ in the same way, aligning terms of the
same degree, we see a Sierpinski triangle rotated 90 degrees.

Sutner, using a well-known connection between the Sierpinski triangle and
parity of binomial coefficients notes [4]

$f_{n}(x)=\sum_{i}{\binom{n+i}{2i+1}x^{i}\mod 2}.$

Sutner showed how to calculate $d_{n}$ as the degree of the GCD of two
polynomials in $\mathbb{Z}_{2}[x]$ [4].

######  Theorem 2.1 (Sutner).

For all $n\in\mathbb{N}$.

$d_{n}=\deg{\gcd\left(f_{n+1}(x),f_{n+1}(x+1)\right)}.$

The recursive definition given in Theorem 2.1 and Sutner’s formula in terns of
the parity of binomial coefficients provide brute force ways to calculate
$f_{n}(x)$. Hunziker, Machiavelo, and Park show the following identity that
makes calculating $f_{n}(x)$ easier when $n$ is divisible by powers of 2 [2].

######  Theorem 2.2 (Hunziker, Machiavelo, and Park).

Let $n=b\cdot 2^{k}$ where $b$ and $k$ are non-negative integers. Then

$f_{n}(x)=x^{2^{k}-1}f_{b}^{2^{k}}(x).$

In particular, we will use this result to relate $f_{2n+2}(x)$ and
$f_{4n+4}(x)$ to $f_{n+1}(x)$.

######  Corollary 2.2.1.

The following identities hold:

$\displaystyle f_{2n+2}(x)$ $\displaystyle=xf_{n+1}^{2}(x)$ $\displaystyle
f_{4n+4}(x)$ $\displaystyle=x^{3}f_{n+1}^{4}(x).$

###### Proof.

Notice that $2n+2=(n+1)2^{1}$ and $4n+4=(n+1)2^{2}$. Thus, our desired
identities follow from Theorem 2.2. ∎

Now that we have a way to express $f_{2n+2}(x)$ and $f_{4n+4}(x)$ as a product
of $f_{n+1}(x)$ and a power of $x$, we need a way to express the GCD of
products so we can relate $d_{2n+1}$ and $d_{n}$. This is exactly what a
number-theoretic result from Ore provides [6].

######  Theorem 2.3 (Ore).

Let $a$, $b$, $c$, and $d$ be integers. Let $(a,b)$ denote $\gcd{(a,b)}$. Then

$(ab,cd)=(a,c)(b,d)\left(\frac{a}{(a,c)},\frac{d}{(b,d)}\right)\left(\frac{c}{(a,c)},\frac{b}{(b,d)}\right).$

Although Ore’s result deals specifically with integers, both the integers and
$\mathbb{Z}_{2}[x]$ are Euclidean domains, so the result will still hold.

Hunziker, Machiavelo, and Park also showed the following identity [2].

######  Theorem 2.4 (Hunziker, Machiavelo, and Park).

A polynomial $\tau(x)$ in $\mathbb{Z}_{2}[x]$ divides both $f_{n}(x)$ and
$f_{m}(x)$ if and only if it divides $f_{\gcd{(m,n)}}$. In particular,

$\gcd{\left(f_{m}(x),f_{n}(x)\right)}=f_{\gcd{(m,n)}}(x).$

We specifically will use the following corollary:

######  Corollary 2.4.1.

For some polynomial $\tau(x)$ in $\mathbb{Z}_{2}[x]$, let $n\geq 0$ be the
smallest integer such that $\tau(x)$ divides $f_{n}(x)$. Then for all $m\geq
0$, $\tau(x)$ divides $f_{m}(x)$ if and only if $n$ divides $m$.

###### Proof.

Let $\tau(x)$ be some polynomial in $\mathbb{Z}_{2}[x]$. Let $f_{n}(x)$ be the
smallest Fibonacci polynomial such that $\tau(x)$ divides $f_{n}(x)$.

Assume that $\tau(x)$ divides $f_{m}(x)$ for some number $m$. Then $\tau(x)$
is a common factor of $f_{m}(x)$ and $f_{n}(x)$, so Theorem 2.4 tells us that
$\tau(x)$ divides $f_{\gcd{(m,n)}}(x)$. Since $f_{n}(x)$ is the smallest
Fibonacci polynomial that is divisible by $\tau(x)$,

$\gcd{(m,n)}\geq n.$

This inequality only holds if $\gcd{(m,n)}=n$. Thus, $m$ must be a multiple of
$n$ as desired.

Now assume that $m$ is a multiple of $n$. Then $\gcd{(m,n)}=n$. Theorem 2.4
tells us

$\gcd{\left(f_{m}(x),f_{n}(x)\right)}=f_{\gcd{(m,n)}}(x)=f_{n}(x).$

Since $\tau(x)$ divides $f_{n}(x)$, and $f_{n}(x)$ is the GCD of $f_{m}(x)$
and $f_{n}(x)$, $\tau(x)$ must also divide $f_{m}(x)$, as desired. ∎

In particular, we will use the following instances of Corollary 2.4.1 to
determine when $\delta_{n}$ is 0 or 2.

######  Corollary 2.4.2.

The following are true:

  1. (i)

$x$ divides $f_{n}(x)$ if and only if $n\equiv 0\mod 2$.

  2. (ii)

$x+1$ divides $f_{n}(x+1)$ if and only if $n\equiv 0\mod 2$.

  3. (iii)

$x+1$ divides $f_{n}(x)$ if and only if $n\equiv 0\mod 3$.

  4. (iv)

$x$ divides $f_{n}(x+1)$ if and only if $n\equiv 0\mod 3$.

###### Proof.

Notice,

  1. (i)

We find that $f_{2}(x)=x$ is the smallest Fibonacci polynomial divisible by
$x$, so we apply Corollary 2.4.1 to get the desired result.

  2. (ii)

Follows from (i) by substituting $x+1$ for $x$.

  3. (iii)

We find that $f_{3}(x)=(x+1)^{2}$ is the smallest Fibonacci polynomial
divisible by $x+1$, so we apply Corollary 2.4.1 to get the desired result.

  4. (iv)

Follows from (iii) by substituting $x+1$ for $x$.

∎

##  3 Proof of Sutner’s Conjecture

Finally, we are ready to prove Sutner’s conjecture [3].

######  Theorem 3.1.

For all $n\in\mathbb{N}$,

$d_{2n+1}=2d_{n}+\delta_{n},$

where $\delta_{n}\in\\{0,2\\}$, and $\delta_{2n+1}=\delta_{n}$.

###### Proof.

Let $(a,b)$ denote $\gcd{(a,b)}$. Applying the results from Theorems 2.1, 2.2,
and 2.3,

$\displaystyle d_{2n+1}$
$\displaystyle=\deg\left(f_{2n+2}(x),f_{2n+2}(x+1)\right)$
$\displaystyle=\deg\left(xf^{2}_{n+1}(x),(x+1)f^{2}_{n+1}(x+1)\right)$
$\displaystyle=\deg(x,x+1)\left(f^{2}_{n+1}(x),f^{2}_{n+1}(x+1)\right)\left(\frac{x+1}{(x,x+1)},\frac{f^{2}_{n+1}(x)}{(f^{2}_{n+1}(x),f^{2}_{n+1}(x+1))}\right)\left(\frac{x}{(x,x+1)},\frac{f^{2}_{n+1}(x+1)}{(f^{2}_{n+1}(x),f^{2}_{n+1}(x+1))}\right)$
$\displaystyle=\deg\left(f_{n+1}(x),f_{n+1}(x+1)\right)^{2}\left(x+1,\frac{f^{2}_{n+1}(x)}{(f_{n+1}(x),f_{n+1}(x+1))^{2}}\right)\left(x,\frac{f^{2}_{n+1}(x+1)}{(f_{n+1}(x),f_{n+1}(x+1))^{2}}\right)$
$\displaystyle=\deg\left(f_{n+1}(x),f_{n+1}(x+1)\right)^{2}\left(x+1,\frac{f_{n+1}(x)}{(f_{n+1}(x),f_{n+1}(x+1))}\right)\left(x,\frac{f_{n+1}(x+1)}{(f_{n+1}(x),f_{n+1}(x+1))}\right)$
$\displaystyle=2d_{n}+\deg\left(x+1,\frac{f_{n+1}(x)}{(f_{n+1}(x),f_{n+1}(x+1))}\right)\left(x,\frac{f_{n+1}(x+1)}{(f_{n+1}(x),f_{n+1}(x+1))}\right).$

Notice that if we substitute $x+1$ for $x$,

$\left(x+1,\frac{f_{n+1}(x)}{(f_{n+1}(x+1),f_{n+1}(x))}\right)\text{ becomes
}\left(x,\frac{f_{n+1}(x+1)}{(f_{n+1}(x),f_{n+1}(x+1))}\right).$

Thus, we see that these two remaining GCD terms in our expression for
$d_{2n+1}$ are either both 1 or not 1 simultaneously. This means we can
further simplify to

$d_{2n+1}=2d_{n}+2\deg\left(x,\frac{f_{n+1}(x+1)}{(f_{n+1}(x),f_{n+1}(x+1))}\right).$

So, we see that

$d_{2n+1}=2d_{n}+\delta_{n}\text{, where
}\delta_{n}=2\deg\left(x,\frac{f_{n+1}(x+1)}{(f_{n+1}(x),f_{n+1}(x+1))}\right).$

Thus, $\delta_{n}\in\\{0,2\\}$ as desired.

What remains is to show that $\delta_{n}=\delta_{2n+1}$. Applying Corollary
2.2.1,

$\displaystyle d_{4n+3}$
$\displaystyle=\deg\left(x^{3}f^{4}_{n+1}(x),(x+1)^{3}f^{4}_{n+1}(x+1)\right)$
$\displaystyle=\deg\left(x^{3},(x+1)^{3}\right)\left(f^{4}_{n+1}(x),f^{4}_{n+1}(x+1)\right)\left(x^{3},\frac{f^{4}_{n+1}(x+1)}{(f^{4}_{n+1}(x),f^{4}_{n+1}(x+1))}\right)\left((x+1)^{3},\frac{f^{4}_{n+1}(x)}{(f^{4}_{n+1}(x),f^{4}_{n+1}(x+1))}\right)$
$\displaystyle=\deg\left(f_{n+1}(x),f_{n+1}(x+1)\right)^{4}\left(x^{3},\frac{f^{4}_{n+1}(x+1)}{(f_{n+1}(x),f_{n+1}(x+1))^{4}}\right)\left((x+1)^{3},\frac{f^{4}_{n+1}(x)}{(f_{n+1}(x),f_{n+1}(x+1))^{4}}\right)$
$\displaystyle=\deg\left(f_{n+1}(x),f_{n+1}(x+1)\right)^{4}\left(x^{3},\frac{f^{3}_{n+1}(x+1)}{(f_{n+1}(x),f_{n+1}(x+1))^{3}}\right)\left((x+1)^{3},\frac{f^{3}_{n+1}(x)}{(f_{n+1}(x),f_{n+1}(x+1))^{3}}\right)$
$\displaystyle=\deg\left(f_{n+1}(x),f_{n+1}(x+1)\right)^{4}\left(x,\frac{f_{n+1}(x+1)}{(f_{n+1}(x),f_{n+1}(x+1))}\right)^{3}\left(x+1,\frac{f_{n+1}(x)}{(f_{n+1}(x),f_{n+1}(x+1))}\right)^{3}$
$\displaystyle=4d_{n}+3\delta_{n}.$

Also, from our work previously in this proof,

$\displaystyle d_{4n+3}$ $\displaystyle=d_{2(2n+1)+1}$
$\displaystyle=2d_{2n+1}+\delta_{2n+1}$
$\displaystyle=2\left(2d_{n}+\delta_{n}\right)+\delta_{2n+1}$
$\displaystyle=4d_{n}+2\delta_{n}+\delta_{2n+1}.$

For these two expressions for $d_{4n+3}$ to be equal, we must have
$\delta_{2n+1}=\delta_{n}$, as desired. ∎

This result may have been proven prior by Yamagishi [5]. However, Yamagishi
does not mention the connection to Sutner’s conjecture, and the proof provided
is not as direct as the one we provide because Yamagishi’s work is concerned
with tori rather than grids.

##  4 When does $\delta_{n}=2$?

Theorem 3.1 proves Sutner’s conjecture as stated and even gives a formula for
finding $\delta_{n}$. However, this formula is somewhat messy, containing one
polynomial division and two polynomial GCDs. We can improve this formula to
just a modulo operation on $n$. We’ll do so by using the divisibility
properties established in Corollary 2.4.2.

######  Theorem 4.1.

The value of $\delta_{n}$ is 2 if and only if $n+1$ is divisible by 3.

###### Proof.

From our work in Theorem 3.1, we know that

$\delta_{n}=2\deg\left(x+1,\frac{f_{n+1}(x)}{(f_{n+1}(x),f_{n+1}(x+1))}\right).$

So we see that $\delta_{n}$ is 2 exactly when $f_{n+1}(x)$ can be divided
without remainder by $x+1$ more times than $f_{n+1}(x+1)$.

For $n+1$ is not divisible by 3, Corollary 2.4.2 tells us that $f_{n+1}(x)$ is
not divisible by $x+1$. So in this case, $\delta_{n}=0$, as desired.

For $n+1$ divisible by 3, let $n+1=b\cdot 2^{k}$ for some integers $b,k\geq 0$
where $b$ is odd. Since $n+1$ is divisible by 3, $b$ must also be divisible by
3. Applying Corollary 2.2.1,

$f_{n+1}(x)=x^{2^{k}-1}f_{b}^{2^{k}}(x)\text{ and
}f_{n+1}(x+1)=(x+1)^{2^{k}-1}f_{b}^{2^{k}}(x+1).$

Since $b$ is an odd multiple of 3, Corollary 2.4.2 tell us that $x+1$ divides
$f_{b}(x)$, but $x+1$ does not divide $f_{b}(x+1)$. So,

$f_{n+1}(x)=x^{2^{k}-1}(x+1)^{2^{k}}g^{2^{k}}(x)\text{ and
}f_{n+1}(x+1)=(x+1)^{2^{k}-1}x^{2^{k}}g^{2^{k}}(x+1),$

for some $g(x)\in\mathbb{Z}_{2}[x]$, where $g(x)$ and $g(x+1)$ are both
divisible by neither $x$ nor $x+1$. So, we see that $f_{n+1}(x)$ can be
divided without remainder by $x+1$ one more time than $f_{n+1}(x+1)$. So,
$\delta_{n}=2$, as desired. ∎

##  5 Future Work

There are many other relationships with $d_{n}$, many yet to be proven. For
example, Sutner mentions that for all $k\in\mathbb{N}$, $d_{2^{k}-1}=0$ [3].
We believe that the following relationships hold, but are unaware of a proof.

######  Conjecture 5.1.

There are infinitely many $n$ such that $d_{n}=2$. In particular, for all
$k\in\mathbb{N}$, $d_{2\cdot 3^{k}-1}=2$.

This conjecture is similar to Sutner’s result that shows there are infinitely
many $n$ such that $d_{n}=0$.

######  Conjecture 5.2.

Let $a$ be an odd natural number. If $a$ is not divisible by 21, then for all
$k\in\mathbb{N}$,

$d_{a^{k}-1}=d_{a-1}.$

Goshima and Yamagishi conjectured a similar statement on tori instead of grids
and for $a$ prime [1].

######  Theorem 5.1.

The case of $a=3$ for Conjecture 5.2 and 5.1 are equivalent.

###### Proof.

For $a=3$, Conjecture 5.2 says that for all $k\in\mathbb{N}$,

$d_{3^{k}-1}=d_{3-1}=0.$

Since $3^{k}$ is divisible by 3, Theorem 4.1 tells us that
$\delta_{3^{k}-1}=2$. So, applying Theorem 3.1,

$d_{2\cdot 3^{k}-1}=2d_{3^{k}-1}+\delta_{3^{k}-1}=2,$

exactly what Conjecture 5.1 states. One can apply all same steps in reverse to
shows that Conjecture 5.1 implies the $a=3$ case of Conjecture 5.2. ∎

## References

  * [1] Masato Goshima and Masakazu Yamagishi, _On the dimension of the space of harmonic functions on a discrete torus_ , Experimental Mathematics 19 (2010), no. 4, 421–429.
  * [2] Markus Hunziker, António Machiavelo, and Jihun Park, _Chebyshev polynomials over finite fields and reversibility of $\sigma$-automata on square grids_, Theoretical Computer Science 320 (2004), no. 2, 465–483.
  * [3] Klaus Sutner, _Linear cellular automata and the Garden-of-Eden_ , The Mathematical Intelligencer 11 (1989), no. 2, 49–53.
  * [4] , _sigma-automata and chebyshev-polynomials_ , Theoretical Computer Science 230 (1996), 49–73.
  * [5] Masakazu Yamagishi, _Periodic harmonic functions on lattices and chebyshev polynomials_ , Linear Algebra and its Applications 476 (2015), 1–15.
  * [6] Øystein Ore, _Number theory and its history_ , p. 109, McGraw-Hill, 1948, Section 5-4, Problem 2.