# Field moment expansion method for interacting Bosonic systems

Andrew Eberhardt  Corresponding author.
<EMAIL_ADDRESS>Kavli Institute for Particle Astrophysics and Cosmology,
Menlo Park, 94025, California, USA  Physics Department, Stanford University,
Stanford, California, USA  SLAC National Accelerator Laboratory  Michael Kopp
Nordita, KTH Royal Institute of Technology and Stockholm University, Hannes
Alfvéns väg 12, SE-106 91 Stockholm, Sweden  Alvaro Zamora  Kavli Institute
for Particle Astrophysics and Cosmology, Menlo Park, 94025, California, USA
Physics Department, Stanford University, Stanford, California, USA  SLAC
National Accelerator Laboratory  Tom Abel  Kavli Institute for Particle
Astrophysics and Cosmology, Menlo Park, 94025, California, USA  Physics
Department, Stanford University, Stanford, California, USA  SLAC National
Accelerator Laboratory

###### Abstract

We introduce a numerical method and python package, CHiMES, that simulates
quantum systems initially well approximated by mean field theory using a
second order extension of the classical field approach. We call this the field
moment expansion method. In this way, we can accurately approximate the
evolution of first and second field moments beyond where the mean field theory
breaks down. This allows us to estimate the quantum breaktime of a classical
approximation without any calculations external to the theory. We investigate
the accuracy of the field moment expansion using a number of well studied
quantum test problems. Interacting Bosonic systems similar to scalar field
dark matter are chosen as test problems. We find that successful application
of this method depends on two conditions: the quantum system must initially be
well described by the classical theory, and that the growth of the higher
order moments be hierarchical.

##  I Introduction

Interacting many body bosonic systems describe a wide array of interesting
phenomena. This includes Bose-Einstein condensates (BEC) Anderson  _et al._
(1995); Davis  _et al._ (1995), electromagnetic radiation Bialynicki-Birula
(1977), and scalar field dark matter (SFDM) Hu  _et al._ (2000); Mocz  _et
al._ (2019); Guth  _et al._ (2015). Their dynamical properties are often
explored using a classical mean field thoery (MFT) approximation, the Gross-
Pitaevskii equations, or Schrödinger-Poisson equations in the case of SFDM
Hertzberg (2016); Hu  _et al._ (2000); Guth  _et al._ (2015); Leggett (2001a);
Carr  _et al._ (2000); Baer (2000); Minguzzi  _et al._ (2004).

Numerically, MFT is preferable to an exact quantum field description which,
for a system with $M$ interacting modes and total number of particles
$n_{tot}$, would involve simulating a Hilbert space of dimensional $\sim
n_{tot}^{M-2}$ Sikivie and Todarello (2017). For large $M$ or $n_{tot}$ an
exact quantum treatment is infeasible. Rather than try and implement an exact
quantum solver it is simpler to extend the classical theory using correction
terms that capture quantum effects on the classical physics Alon  _et al._
(2007); Prezhdo and Pereverzev (2000); Boßmann  _et al._ (2020); Heller
(1976); Royer (1992); Sreedharan  _et al._ (2020); Lentz  _et al._ (2019).

When occupation numbers are large and interactions weak the MFT is known to
accurately describe the dynamics of these systems Leggett (2001a); Bialynicki-
Birula (1977); Gross (1961); Pitaevskii (1960); Leggett (2001b). However, any
interacting system with a nonlinearity will exhibit quantum corrections on
some time scale Yurke and Stoler (1986); Sreedharan  _et al._ (2020);
Lewenstein and You (1996); Caballero-Benitez  _et al._ (2008); Chakrabarty
_et al._ (2018); Chakrabarty (2021); Kopp  _et al._ (2021). The MFT, tracking
only the mean value of the field operator, cannot account for these quantum
terms Alon  _et al._ (2007). This means that effects like phase diffusion,
quantum squeezing, and fragmentation inherently require beyond MFT approach
Alon  _et al._ (2007); Sreedharan  _et al._ (2020); Yurke and Stoler (1986);
Minguzzi  _et al._ (2004). The effect of these corrections on MFT is a current
topic of interest Dvali and Zell (2018); Hertzberg (2016); Sikivie and
Todarello (2017); Dvali  _et al._ (2017); Sreedharan  _et al._ (2020);
Chakrabarty  _et al._ (2018); Kopp  _et al._ (2021); Alon  _et al._ (2007);
Sreedharan  _et al._ (2020); Lentz  _et al._ (2019, 2020). This motivates the
development of numerical methods which can capture beyond MFT physics Alon
_et al._ (2007); Prezhdo and Pereverzev (2000); Boßmann  _et al._ (2020);
Heller (1976); Royer (1992); Sreedharan  _et al._ (2020); Lentz  _et al._
(2019). We will refer to the time at which MFT can no longer accurately
approximate the evolution of the underlying system as the “quantum breaktime”.

The classical theory is generally achieved as a limit of the quantum field
theory. When the expectation value of the field operator is large compared to
the variance of the field operator it is sensible to replace the field
operators in the quantum field theory with their expectation values Ballentine
. The expectation value of the quantum field is then called the classical
field. Quantum coherent states, with parameters large compared to unity,
satisfy this approximation criterion Ballentine . Tracking higher order
moments and their effects on the evolution of observables has been studied
with success for position and momentum operators Heller (1976); Prezhdo and
Pereverzev (2000). Likewise, conceptually similar methods of expanding the
field about its mean value have proven useful as correction MFT Royer (1992);
Georges and Yedidia (1991); Lentz  _et al._ (2019) as well as helping explore
interesting physics Glimm  _et al._ (1976); Boßmann  _et al._ (2020); Tanaś
and Kielich (1983); Leggett (2001a).

We apply these techniques to extend MFT approximation to include terms
proportional to higher order central moments and then integrate coupled
differential equations governing the evolution of the the mean and these
higher order moments. We will refer to this method as the “field moment
expansion” (FME). The main focus of this work will be to introduce a publicly
available code repository which implements a solver which tracks both the mean
field values and the second moments of the field. In this way we can simulate
systems initially well described by MFT into regimes where quantum corrections
become important.

There are two main benefits to this approach. First, prior to the quantum
breaktime, FME produces a more accurate approximation of the quantum field
expectation values than MFT. Secondly, FME provides an internal estimation of
the breaktime, and can therefore estimate its own regime of validity. This
internal assessment means no calculations external to the theory are necessary
to estimate the quantum breaktime. This is in contrast to approximations of
this timescale calculated by methods external to MFT, see for example Dvali
and Zell (2018); Dvali  _et al._ (2013, 2017). Additionally, when compared
with other mean field extensions Alon  _et al._ (2007), the FME scales as
$\mathcal{O}(M^{2}\,\log M)$, depending only on the grid size as opposed to
the particle occupation numbers $n_{p}$.

The method is applicable to any interacting scalar field system assuming that
a number of criteria are met. First, the initial correction terms must be
subleading order. The method predictions cannot be trusted past the time when
the quantum corrections become large. However, we will show that the FME is
able to approximate the evolution of the quantum system longer than MFT.
Second, the correction terms must grow hierarchically. For a term, $F$ of
order $m$ written as a function of moments of order less than and equal to
$m$, given $F^{m}(\text{moments of order $\leq m$})$, the evolution the terms
must satisfy $F^{1}>F^{2}>F^{3}\dots$. Generally, these criteria will be met
if the initial conditions are a coherent state with large mode occupation
numbers.

In this work we test the field moment expansion using two test problems that
have been well studied in the literature, for which exact quantum solutions
are possible and that exhibit a breakdown of the MFT on some timescale Yurke
and Stoler (1986); Hertzberg (2016); Sikivie and Todarello (2017); Erken  _et
al._ (2012). For each we show that the FME provides a more accurate solution
until the quantum breaktime. Most importantly, we demonstrate that the method
can be used in this case to accurately predict the quantum breaktime.

The paper is organized as follows. In Section II, we discuss background on
interacting Bosonic systems, the MFT, and FME approximations. Section III
explains our numerical implementation. In Section IV, we demonstrate that FME
is accurate for a number of quantum test problems. Conclusions regarding the
overall utility of these methods and future directions are presented in
Section V.

##  II Background

###  II.1 Quantum description

We start from the following Hamiltonian, used to describe non-relativistic
scalar fields Hertzberg (2016); Sikivie and Todarello (2017); Erken  _et al._
(2012).

$\hat{H}=\sum_{j}^{M}\omega_{j}\hat{a}_{j}^{\dagger}\hat{a}_{j}+\sum_{ijkl}^{M}\frac{\Lambda_{kl}^{ij}}{2}\hat{a}_{k}^{\dagger}\hat{a}_{l}^{\dagger}\hat{a}_{i}\hat{a}_{j}\,.$
(1)

Where the sums are performed over the $M$ system modes. For appropriately
chosen $\Lambda$, $\omega$, and $M$ this Hamiltonian describes a wide range of
physical systems. The first sum describes the kinetic energy of the system and
the second the self interactions. $\hat{a}_{i}$ is the annihilation operator
on mode $i$, which is defined by its commutation and action on number
eigenstates. A number eigenstate is written as

$\ket{\set{n}}=\ket{n_{1},n_{2},\dots,n_{M}}\,,$ (2)

where $n_{i}$ describes the number of particles occupying the $i$th mode. Here
we will take the modes to represent momentum eigenstates with momentum
$p_{i}$. The number eigenstates form an orthonormal basis such that
$\braket{n^{\prime}_{j}}{n_{i}}=\delta_{i\,j}\delta_{n\,n^{\prime}}$. We can
now describe the $\hat{a}$ operators as follows

$\displaystyle[\hat{a}_{i},\hat{a}_{j}]$ $\displaystyle=0\,,$ (3)
$\displaystyle[\hat{a}_{i},\hat{a}^{\dagger}_{j}]$
$\displaystyle=\delta_{ij}\,,$ (4)
$\displaystyle\hat{a}_{j}^{\dagger}\ket{n_{j}}$
$\displaystyle=(n_{j}+1)^{1/2}\ket{n_{j}+1}\,,$ (5)
$\displaystyle\hat{a}_{j}\ket{n_{j}}$
$\displaystyle=n_{j}^{1/2}\ket{n_{j}-1}\,,$ (6)
$\displaystyle\hat{N}_{j}\ket{n_{j}}$
$\displaystyle\equiv\hat{a}_{j}^{\dagger}\hat{a}_{j}\ket{n_{j}}=n_{j}\ket{n_{j}}\,\,.$
(7)

The annihilation operator can also be used to define the complex quantum field
$\hat{\psi}(x)$, which is related to $\hat{a}$ by Fourier transform.

$\hat{\psi}(x)=\sum_{i}\hat{a}_{i}u^{\dagger}_{i}(x)$ (8)

where $u^{\dagger}_{i}(x)$ is the eigenstate of the momentum operator with
eigenvalue $p_{i}$. The Heisenberg equation describes the dynamics of these
operators. For an arbitrary operator $\hat{A}$ with time independent
Hamiltonian the equation of motion is written

$\partial_{t}\hat{A}=\frac{i}{\hbar}[\hat{H},\hat{A}]\,.$ (9)

Hereafter we set $\hbar\equiv 1$. We can now solve for the evolution of our
field operators,

$\displaystyle\partial_{t}\hat{a}_{p}$
$\displaystyle=i[\hat{H},\hat{a}_{p}]=-i\left[\omega_{p}\hat{a}_{p}+\sum_{ijl}\Lambda^{ij}_{pl}\hat{a}_{l}^{\dagger}\hat{a}_{i}\hat{a}_{j}\right]\,.$
(10)

In this work we will be taking the constants $\Lambda^{ij}_{pl}$ to be of the
following form

$\Lambda^{ij}_{pl}=\left(\frac{C}{2(p_{p}-p_{i})^{2}}+\frac{C}{2(p_{p}-p_{j})^{2}}+\Lambda_{0}\right)\delta^{ij}_{pl}\,,$
(11)

where the constant $C$ describes a long range $r^{-1}$ potential and
$\Lambda_{0}$ characterizes the strength of contact interactions.
$\delta^{ij}_{pl}$ is the Kroneker delta. If $C=\frac{-4\pi Gm^{2}}{L}$ and
$\Lambda_{0}=0$, then taking a Fourier transform of equation (33) yields the
familiar second quantized Schrödinger–Poisson equations. Where $G$ is the
gravitational constant, $L$ is the volume of the box for which the quantum
field is periodic, and $m$ is the mass of the field. Here we are working in
1-D.

$\displaystyle\partial_{t}\hat{\psi}(x)$
$\displaystyle=-i\left[\frac{-\nabla^{2}}{2m}+m\hat{V}(x)\right]\hat{\psi}(x)\,,$
(12) $\displaystyle\nabla^{2}\hat{V}$ $\displaystyle=4\pi
Gm\,\hat{\psi}^{\dagger}(x)\hat{\psi}(x)\,.$ (13)

Here we have started with a complex quantum field operator. However, it is
possible to derive this set of equations as a non-relativistic and weak
gravity limit of the real scalar Klein Gordon field. This can be done
following the derivation in Suárez and Chavanis (2015); Lentz  _et al._
(2019). These limits need to be kept in mind when determining where this set
of equations and approximations of it are valid.

While the above analysis is true for an arbitrary quantum state, within the
stated limits, it is useful to define an initial quantum state for which the
mean field theory starts as an accurate approximation of the quantum field
theory. The “most classical” state is the coherent state, for which MFT is
initially exact, parameterized by the complex vector
$\vec{z}\in\mathbb{C}^{M}$, which can be written as a sum of number
eigenstates as

$\ket{\vec{z}}=\bigotimes_{i=1}^{M}\exp\left[-\frac{|z_{i}|^{2}}{2}\right]\sum_{n_{i}=0}^{\infty}\frac{z_{i}^{n_{i}}}{\sqrt{n_{i}!}}\ket{n_{i}}\,.$
(14)

When representing a coherent state numerically we truncate the above sum when
the square norm of $\braket{\vec{z}}{\vec{z}}\geq.995$.

A coherent state is thought to describe the initial state of the axion field
if produced via the misalignment mechanism Abbott and Sikivie (1983); Preskill
_et al._ (1983).

This state implies that a measurement of the particle number in the mode $i$
would be Poisson distributed with expectation value $|z_{i}|^{2}$. This state
is special because it has the property that the expectation value any normally
ordered operator composed of $\hat{a}$ and $\hat{a}^{\dagger}$ with respect to
this state is given by simply replacing the operators with the parameter $z$,
i.e.

$\braket{\vec{z}}{\,f(\set{\hat{a}^{\dagger}})\,g(\set{\hat{a}})\,}{\vec{z}}=f(\vec{z}^{\dagger})\,g(\vec{z})\,.$
(15)

This will be important when deriving the mean field theory.

###  II.2 Mean field approximation

The mean field is simply the expectation value of the field operator, we
define a mean field in position and momentum space respectively as

$\displaystyle\psi(x)$ $\displaystyle\equiv\braket{\hat{\psi}(x)}\,,$ (16)
$\displaystyle a_{i}$ $\displaystyle\equiv\braket{\hat{a}_{i}}\,.$ (17)

The higher order moments can then be calculated from the mean field operators.
Occupation numbers, $N^{cl}$ are given as the amplitude of the field
operators, e.g.

$N_{i}^{cl}=|a_{i}|^{2}$ (18)

The mean field theory is attained simply by taking an expectation value of the
equations of motion and then approximating the operators by their expectation
value. Let us say that the operator $\hat{A}$ evolves according to the
following equation of motion

$\partial_{t}\hat{A}=f(\hat{A})\,.$ (19)

And that $\hat{A}$ corresponds to some dynamic observable of the system. If
the expectation value of $\hat{A}$ is large compared to its root variance than
we can make the following approximation

$\partial_{t}\braket{\hat{A}}=\braket{f(\hat{A})}\approx
f(\braket{\hat{A}})\,.$ (20)

This approximation is one way to transition to a mean field theory
description. It is identical to Ehrenfest theorem if we replaced the position
and momentum operators with field operators.

This statement that the mean field theory is accurate at a time $T$ implies
the following conditions

$\mathrm{E}[\hat{A}(T)]\gg\sqrt{\mathrm{Var}[\hat{A}(T)]}\,.$ (21)

Note this is a condition both on the evolution of the operator $\hat{A}$ and
the quantum state that the expectation is taken with respect to.

For the mean field approximation to hold we need only that the approximation
in equation (20) remain accurate on the time scale of the evolution. We will
make the requirements more precise in the next section, but from here we can
see qualitatively why large occupation number tends to motivate the mean field
theory approximation.

Let us assume that $\hat{A}=\set{\hat{a},\hat{a}^{\dagger},\dots}$ is the set
of operators generated by field operators $\hat{a}$ and $\hat{a}^{\dagger}$,
as will be the case in the next sections. It is always possible to write the
right hand side of equation (19) in terms of normally ordered operators. This
means that if we are in a coherent state by equation (15) the mean field
approximation is an equality so long as we remain in a coherent state. This is
true regardless of expected occupation number, however, if we are only in an
approximately coherent state with $|z|^{2}=n\gg 1$ then we can expect the
variance in our field operators to be approximately governed by a Poisson
distribution; i.e.
$\mathrm{E}[\hat{N}]\sim\mathrm{Var}[\hat{N}]\sim|z|^{2}=n$, this means that
the fractional deviation in the expectation value will go as
$\sqrt{\mathrm{Var}[\hat{N}]}/\mathrm{E}[\hat{N}]\sim 1/\sqrt{N}$. This is
easily made into an estimate of the fractional field variance by recalling
that $\hat{N}=\hat{a}^{\dagger}\hat{a}$. So the fractional deviation in the
field values is small for large occupation numbers $n\gg 1$.

Large occupation numbers is not enough on its own to ensure an accurate MFT
description. The mean field theory will be accurate if the occupation numbers
be large but also that the quantum state remained approximately coherent, or
equivalently, that the distribution of the number operator has Poisson
distributed expectation and variance proportional to some power of $n$. But
the actual condition that needs to be satisfied is equation (21).

It is easy to imagine a state for which this assumption is not met. A number
eigenstate, $\ket{\set{n}}$, for example has field expectation
$\braket{\set{n}}{\,\hat{a}\,}{\set{n}}=0$ even when $n\gg 1$. Therefore, it
does not satisfy the condition in equation (21) and will not be well described
by a single classical field even for large $n$. It was demonstrated in Sikivie
and Todarello (2017); Chakrabarty (2021) that a number eigenstate did not
approach a single classical field description even at large occupation number.
It was then shown in Hertzberg (2016) that this state could be approximated by
an ensemble of classical fields with ensembled expectation $0$ and amplitude
$n$.

It is important to note that large occupation number is only a proxy for the
accuracy of the mean field theory. For example, a coherent state evolved by
the free particle Hamiltonian will always be perfectly described by classical
field theory even as $\mathrm{E}[\hat{N}]\ll 1$. Conversely, a number
eigenstate will never be well described by a single classical field even as
$\mathrm{E}[\hat{N}]\gg 1$. Going forward we will phrase our estimation of the
accuracy of the classical field theory in terms of the condition described in
equation (21) for the field operator $\hat{a}$.

It should also be noted that it is possible to reproduce the quantum evolution
of the number operator without reproducing the evolution of the field itself.
For example, a so-called field number state in the large $N$ limit approaches
the classical evolution of the mode occupations, however, it has vanishing
field expectation regardless of $N$. For states like these, equation (21)
should be expressed in terms of the number operator.

###  II.3 Quantum corrections

Systems initially well described by the mean field theory will eventually
diverge from this description on some time scale if the Hamiltonian is non-
Harmonic. The specific causes of this are of interest in the literature
Lewenstein and You (1996); Sreedharan  _et al._ (2020); Caballero-Benitez  _et
al._ (2008); Yurke and Stoler (1986); Kopp  _et al._ (2021) but for our
purposes we can think of them generally as a delocalization in phase space. In
this section we will show an example of how deviation from the classical field
theory occurs, and discuss a way to paramterize it.

Specifically that the variance in the field operators becomes of order the
expectation value violating the condition in equation (21). We will clarify
this point using the following example.

Consider the toy Hamiltonian on one mode

$\hat{H}=\frac{\hat{p}^{2}}{2m}+\lambda_{L}\hat{q}^{2}+\lambda_{NL}\hat{q}^{4}\,.$
(22)

Where $p\equiv-i\hbar\nabla_{q}$ The first two terms define a harmonic
oscillator and the last term some non-linearity. The set of operators that we
are interested in are the position and momentum operators
$\hat{A}=\set{\hat{q},\hat{p}}$. Note that this Hamiltonian can be recast in
terms of the $\hat{a}$ operator using the relation

$\hat{a}=\frac{1}{\sqrt{2}}(\hat{q}+i\hat{p})\,,$ (23)

but this does not change the analysis.

The classical equations of motion can be found using equations (9) and (22) as
well as the commutation relation $[q,p]=i\hbar$ and then applying Ehrenfest’s
theorem. For the classical variables $p$, $q$, related to the classical field
by $a=\frac{1}{\sqrt{2}}(q+ip)$ (recall there is only a single mode), we
obtain the following equations of motion

$\displaystyle\partial_{t}p$
$\displaystyle=-2\,\lambda_{L}\,q-4\,\lambda_{NL}\,q^{3}$ (24)
$\displaystyle\partial_{t}q$ $\displaystyle=p$ (25)

With $(q,p)|_{t=0}=(q_{i},0)$. We use a symplectic leap-frog integrator to
solve the classical equations of motion.

Our initial wavfunction will be Gaussian, which is initially well localized in
$q$,

$\braket{q}{\phi}=\sqrt{\frac{1}{\sqrt{2\pi\sigma}}}\exp\left(-\frac{(q-q_{i})^{2}}{4\sigma^{2}}\right)$
(26)

The quantum equations can be solved by integrating Schrödinger’s equation

$\partial_{t}\ket{\phi}=-i\hat{H}\,\ket{\phi}\,.$ (27)

We use a symplectic spectral leap-frog integrator to perform the integration.

Let us also parameterize the inequality in equation (21) by defining the
following quantity $Q_{q}$

$Q_{q}\equiv\frac{\braket{\hat{q}^{2}}-\braket{\hat{q}}^{2}}{q_{i}^{2}}\,.$
(28)

This can be used to parameterize the quantum theories deviation from the
classical theory because the nonlinearity is spatial. It should be noted that
this is not the only parameter that can be constructed with this property.

Figure 1: Here we show how quantum corrections cause deviations from the mean
field theory. Plotted are the Husimi functions for two different non-linear
models at three different times. The top row corresponds to the weakly
nonlinear case with $\lambda_{NL}=0.073$, and the bottom row to the strongly
nonlinear case with $\lambda_{NL}=73$. Each column corresponds to a different
time, $T$ in the evolution. The red dot in each panel indicates the solution
solved obtained using the classical mean field thoery, $a$. The cyan triangle
in each panel represents the actual mean of the $\hat{a}$ operator
$\braket{\hat{a}}$. At all times the weakly nonlinear model closely adheres to
the classical solution. The strongly nonlinear case has the classical field
theory start as a good approximation but strongly diverge over the course of
the evolution.

In figure 1 we track a quantum phase space analog, the Husimi function
Ballentine , of the wavefunction, the classical mean field theory
approximation of $\mathrm{E}[\hat{a}]$, the exact quantum value of
$E[\hat{a}]$ for two different strength nonlinearities. In both cases the
initial spread of the wavefunction is $\sigma=0.025$, the initial location
$q_{i}=0.25$, $\hbar=0.01$, $\lambda_{L}=100$. This means that our parameter
$Q_{q}=0.01$. Note also that equation (23) implies that the occupation number
here is $n=q_{i}^{2}=0.0625$. Clearly we are not in the high occupation
regime. The top row and bottoms rows have $\lambda_{NL}=0.073$ and
$\lambda_{NL}=73$, representing the weakly nonlinear and strongly nonlinear
respectively.

For the nonlinear case we see in the second panel that squeezing and phase
diffusion have caused the Husimi function to be poorly localized around the
classical solution, however, at this time the Husimi function is still well
approximated by a squeezed Gaussian. At this time the classical and quantum
solutions begin to diverge. This is approximately the quantum breaktime, where
the $Q_{q}$ parameter is starting to approach $\mathcal{O}(1)$. In the
rightmost panel of the nonlinear evolution we see that phase diffusion has
removed most of the information about the phase of the field. The classical
and quantum solutions now deviate by an $\mathcal{O}(1)$ fraction.

Throughout this work we will be interested in functions of the field operators
and their Fourier transforms. We will define the quantum breaktime as we have
done in this section, by first defining a quantity that paramterizes the
deviation from the classical theory. The following definition will be used
moving forward

$Q\equiv\sum_{i}^{M}\frac{\braket{\delta\hat{a}_{i}^{\dagger}\delta\hat{a}_{i}}}{n_{tot}}\,.$
(29)

Where $n_{tot}\equiv\sum_{i}^{M}n_{i}$ and
$\delta\hat{a}\equiv\hat{a}-\braket{\hat{a}}$.

When $Q$ is small a single classical field can accurately capture both the
first and second moments of the field operator and the classical approximation
in equation (18) is valid. When this parameter ceases being small, such a
description inaccurately capture both the field and occupation number
expectations. Therefore, we will use $Q$ to define a quantum breaktime,
$t_{br}$ condition as when $Q=.15$, i.e.

$Q(t_{br})\equiv 0.15\,.$ (30)

Note that this is not intended to be the unique usable definition of the
quantum breaktime. It is only intended to allow us to quantify when the
assumption in equation (18) breaks down. The specific value $0.15$ is
arbitrary and the scaling of the breaktime with occupation number should be
relatively insensitive to the specific choice of $Q(t_{br})$. We choose this
specific value because it indicates that the correction terms are becoming the
same order as the leading order classical terms and analysis of our test
problems indicate that the field moment solver reliably reproduces quantum
results when $Q<0.15$.

###  II.4 Field moment expansion

For well behaved probability distributions with well defined moments, we can
uniquely identify the distribution by its moments. This is true for functions
of quantum mechanical operators. Consider a set of operators $\hat{A}$ and
function of this set $f(\hat{A})$. Let us assume we can write $f$ as a sum of
products of normally ordered generators of $\hat{A}$, and that there exists
some integer $R$ for which every term contains $R$ or fewer elements of the
generator of $\hat{A}$. Let us also assume that the $i$th moment of each
element of $\hat{A}$ is well defined for $i\leq R$. We can then write the
expectation value of $f$ as a sum of terms weighted by central moments as

$\braket{f(\hat{A})}=f(\braket{\hat{A}})+\sum_{j=2}^{R}\frac{1}{j!}\left[\prod_{k=1}^{j}\left(\sum_{\hat{a}\in\hat{A}}\hat{\delta}_{\hat{a}}\right)\right]f(\braket{\hat{A}})\,.$
(31)

Where $\braket{\hat{A}}\equiv\set{\braket{a_{1}},\braket{a_{2}},\dots}$ is the
set of the expectation values of the elements of the set $\hat{A}$, and the
product of the $\hat{\delta}_{\hat{a}}$ operators is defined as

$\hat{\delta}_{\hat{a}_{1}}\hat{\delta}_{\hat{a}_{2}}\dots\equiv\braket{\delta\hat{a}_{1}\,\delta\hat{a}_{2}\,\dots}\frac{\partial}{\partial\braket{\hat{a}_{1}}}\frac{\partial}{\partial\braket{\hat{a}_{2}}}\dots\,.$
(32)

Where the operators $\delta\hat{a}$ are normally ordered.

Consider the Hamiltonian in equation (1). We can solve for the equation of
motion for the $\hat{a}$ operator using equation (9) giving

$\displaystyle\partial_{t}\hat{a}_{p}$
$\displaystyle=i[\hat{H},\hat{a}_{p}]=-i\left[\omega_{p}\hat{a}_{p}+\sum_{ijl}\Lambda^{ij}_{pl}\hat{a}_{l}^{\dagger}\hat{a}_{i}\hat{a}_{j}\right]\,.$
(33)

We see that the highest order operator in the equation of motion is third
order in $\hat{a}$ and $\hat{a}^{\dagger}$. Replacing the $\hat{a}$ and
$\hat{a}^{\dagger}$ operators with their expectation values in this equation
gives the classical field theory given by the following inequality

$\displaystyle\partial_{t}\braket{\hat{a}_{p}}$
$\displaystyle=i\braket{[\hat{H},\hat{a}_{p}]}$
$\displaystyle=-i\left[\omega_{p}\braket{\hat{a}_{p}}+\sum_{ijl}\Lambda^{ij}_{pl}\braket{\hat{a}_{l}^{\dagger}\hat{a}_{i}\hat{a}_{j}}\right]$
(34)
$\displaystyle\approx-i\left[\omega_{p}\braket{\hat{a}_{p}}+\sum_{ijl}\Lambda^{ij}_{pl}\braket{\hat{a}_{l}^{\dagger}}\braket{\hat{a}_{i}}\braket{\hat{a}_{j}}\right]\,.$
(35)

We can rewrite equation (34) the form of equation (31).

$\displaystyle\partial_{t}$
$\displaystyle\braket{\hat{a}_{p}}=i\braket{[\hat{H},\hat{a}_{p}]}$
$\displaystyle=-i\left[\omega_{p}\braket{\hat{a}_{p}}+\sum_{ijl}\Lambda^{ij}_{pl}\braket{\hat{a}_{l}^{\dagger}\hat{a}_{i}\hat{a}_{j}}\right]$
$\displaystyle=-i\left[\omega_{p}\braket{\hat{a}_{p}}+\sum_{ijl}\Lambda^{ij}_{pl}\left(\braket{\hat{a}_{l}^{\dagger}}\braket{\hat{a}_{i}}\braket{\hat{a}_{j}}\right.\right.$
(36a)
$\displaystyle+\braket{\delta\hat{a}_{i}\delta\hat{a}_{j}}\braket{\hat{a}^{\dagger}_{l}}+\braket{\delta\hat{a}^{\dagger}_{l}\delta\hat{a}_{i}}\braket{\hat{a}_{j}}+\braket{\delta\hat{a}^{\dagger}_{l}\delta\hat{a}_{j}}\braket{\hat{a}_{i}}$
(36b)
$\displaystyle+\left.\left.\braket{\delta\hat{a}^{\dagger}_{l}\delta\hat{a}_{i}\delta\hat{a}_{j}}\right)\right]\,.$
(36c)
$\displaystyle\approx-i\left[\omega_{p}\braket{\hat{a}_{p}}+\sum_{ijl}\Lambda^{ij}_{pl}\left(\braket{\hat{a}_{l}^{\dagger}}\braket{\hat{a}_{i}}\braket{\hat{a}_{j}}+\braket{\delta\hat{a}_{i}\delta\hat{a}_{j}}\braket{\hat{a}^{\dagger}_{l}}+\braket{\delta\hat{a}^{\dagger}_{l}\delta\hat{a}_{i}}\braket{\hat{a}_{j}}+\braket{\delta\hat{a}^{\dagger}_{l}\delta\hat{a}_{j}}\braket{\hat{a}_{i}}\,\right)\right].$
(36d)

We see in this form that the classical equations of motion are given by
equations (36a), then (36b) and (36c) act as “quantum” corrections to the mean
field theory. We can also see now the manner in which our classicality
condition in equation (21) is technically imprecise. What is actually required
is that the terms in (36b) and (36c) remain small compared to the terms in
(36a). However, it is important to keep in mind that the accuracy of the mean
field theory is not intrinsically a function of occupation number but instead
a property of how the moments of the field operators compare.

The condition in equation (21) comes about by assuming that the moments grow
hierarchically, that is that the first order terms in (36a) start out largest
and that the terms in (36b) grow faster than the terms in (36c) and so on. If
we then also assume that the second order central moments in (36b) are all
approximately the same order, and that the dynamics are approximately number
conserving, we see that taking a ratio of the first and second order terms
gives us the parameter $Q$ and asserting that $Q\ll 1$ is equivalent then to
equation (21).

In order to integrate equation (36) we couple the evolution of the field
operator $\hat{a}$ to the evolution of the central moments. If we assume that
the central moments grow hierarchically and we are interested in evolution
only until the classical field theory breaks we can truncate equation (36) at
the second order terms, this has the benefit of better computational scaling.
The equations of motion for the higher order moments can be found using
equation (9) and then expanding to second order using equation (31). We will
start with the $\mathrm{Cov}[\hat{a}_{i},\hat{a}_{j}]$ operator which has the
following equation of motion.

$\displaystyle\partial_{t}\braket{\delta\hat{a}_{i}\delta\hat{a}_{j}}$
$\displaystyle\approx-i\left[(\omega_{i}+\omega_{j})\braket{\delta\hat{a}_{i}\delta\hat{a}_{j}}+\sum_{kp}\Lambda^{kp}_{ji}\braket{\hat{a}_{k}}\braket{\hat{a}_{p}}\right.$
(37)
$\displaystyle\left.+\sum_{kpl}\Lambda^{kp}_{ji}\left(\braket{\delta\hat{a}_{l}\delta\hat{a}_{j}}\braket{\hat{a}^{\dagger}_{k}}\braket{\hat{a}_{p}}+\braket{\delta\hat{a}_{j}\delta\hat{a}_{p}}\braket{\hat{a}^{\dagger}_{k}}\braket{\hat{a}_{l}}+\braket{\delta\hat{a}_{k}^{\dagger}\delta\hat{a}_{j}}\braket{\hat{a}_{l}}\braket{\hat{a}_{p}}\right)+(i\leftrightarrow
j)\right]\,.$

This equation can be broken down into three types of terms. The first term and
the second line of terms are the kinetic and potential terms respectively. The
terms are proportional to the covariance operator. The second term on the
second line is proportional both to the potential energy and
$[\hat{a},\hat{a}^{\dagger}]$. This term guarantees that the covariance
operator will grow even if initially zero so long as their is some
nonlinearity in the Hamiltonian.

We can solve for the $\mathrm{Cov}[\hat{a}^{\dagger}_{i},\hat{a}_{j}]$
operator in the same manner.

$\displaystyle\partial_{t}\braket{\delta\hat{a}_{i}^{\dagger}\delta\hat{a}_{j}}$
$\displaystyle\approx
i\left[(\omega_{i}-\omega_{j})\braket{\delta\hat{a}_{i}^{\dagger}a_{j}}+\right.$
(38)
$\displaystyle\left.+\sum_{kpl}\Lambda^{kp}_{ji}\left(\braket{\delta\hat{a}_{j}\delta\hat{a}_{k}}\braket{\hat{a}^{\dagger}_{p}}\braket{\hat{a}^{\dagger}_{l}}+\braket{\delta\hat{a}^{\dagger}_{p}\delta\hat{a}_{j}}\braket{\hat{a}^{\dagger}_{l}}\braket{\hat{a}_{k}}+\braket{\delta\hat{a}^{\dagger}_{l}\delta\hat{a}_{l}}\braket{\hat{a}^{\dagger}_{p}}\braket{\hat{a}_{k}}\right)+(c.c.,\,i\leftrightarrow
j)\right]\,.$

Which has the same structure as the previous equation without the term
proportional to the commutation operator. $c.c$ indicates complex conjugate.

The term in equation (37) proportional to the the commutation between
$\hat{a}$ and $\hat{a}^{\dagger}$ gives us a qualitative sense of how large
occupation number implies classicality. If the quantum state is approximately
coherent then these second order central moments are near $0$ by equation
(15). On a timescale $\sim\mathcal{O}(\Lambda^{-1})$ the second order central
moments will have grown by a factor $[\hat{a},\hat{a}^{\dagger}]=1$. In the
large occupation number limit $[\hat{a},\hat{a}^{\dagger}]=1\ll 1$ meaning our
lowest order quantum corrections contribute vanishingly to the evolution of
the mean field.

###  II.5 Penrose-Onsager criterion

When the Penrose-Onsager (PO) criterion Penrose and Onsager (1956) is
satisfied we can write

$\braket{\hat{a}^{\dagger}_{i}\hat{a}_{j}}=\vec{z}_{i}^{\dagger}\vec{z}_{j}\,.$
(39)

That is that the expectation values of the second field moment matrix
$M_{ij}\equiv\braket{\hat{a}^{\dagger}_{i}\hat{a}_{j}}$ can be written as an
outer product of a single vector $\vec{z}$ with its complex conjugate.

When the PO criterion is satisfied $\hat{M}_{ij}$ contains a single nonzero
eigenvalue, called the principal eigenvalue, equal to the square norm of
$\vec{z}$, i.e. $\lambda_{p}=\sum_{i}|z_{i}|^{2}$. Where
$\vec{z}/\sqrt{\sum_{i}|z_{i}|^{2}}$ is the corresponding principal
eigenvector, $\vec{\xi}_{p}$. This vector is not a priori equal to the
classical field but when the classical field adequately describes the system
we expect the PO criterion to be satisfied.

Note that $\mathrm{Tr}[\hat{M}_{ij}]=n_{tot}$ and, therefore, in number
preserving systems the trace of $\hat{M}_{ij}$ is a conserved quantity,. When
the system is well described by the classical theory we expect that the
principal eigenvalue is very close to $n_{tot}$, more specifically we expect
Leggett (2001a)

$\frac{\lambda_{p}}{n_{tot}}-1\ll 1\,.$ (40)

Because the FME tracks both second central moment
$\braket{\delta\hat{a}_{i}^{\dagger}\delta\hat{a}_{j}}$ and
$\braket{\hat{a}_{i}}$ we can use this method to approximate $\hat{M}_{ij}$ as

$\hat{M}_{ij}^{\mathrm{FME}}=\braket{\delta\hat{a}_{i}^{\dagger}\delta\hat{a}_{j}}^{\mathrm{FME}}+(\braket{\hat{a}_{i}}^{\dagger}\braket{\hat{a}_{j}})^{\mathrm{FME}}\,.$
(41)

It is important to note that while both $Q$ and equation (40) parameterize the
deviation from the classical theory the two are technically distinct in the
following way:

When $Q\not\ll 1$ it implies that the classical approximation in equation (18)
is breaking down, i.e. $Q\sim
1\,,\rightarrow|\braket{\hat{a}}|^{2}\neq\braket{\hat{a}^{\dagger}\hat{a}}$.
Note that this is not a useful parameterization for states which track the
evolution of the mode occupations but not the field values themselves. On the
other hand, $\lambda_{p}/n_{tot}-1\not\ll 1$ implies that $\hat{M}_{ij}$
cannot be described by a single eigenvector. Neither explicitly implies that
the classical field poorly approximates occupation numbers, but both can be
used to give an approximate sense of how closely a system is adhering to the
classical field theory.

##  III Numerical implementation

The full code repository for the simulation and data analyses of the classical
and expanded field theories performed here is publicly available at
https://github.com/andillio/CHiMES.

###  III.1 Mean field theory

The evolution of the mean field, $a^{cl}$, is solved by integrating the
classical field equations of motion given

$\displaystyle\partial_{t}a^{cl}_{p}=-i\left[\omega_{p}a^{cl}_{p}+\sum_{ijl}\Lambda^{kp}_{ji}a^{cl\dagger}_{l}a^{cl}_{i}a^{cl}_{j}\right]\,.$
(42)

The initial conditions are chosen such that the field values correspond to a
coherent state with parameter $\vec{z}=\vec{a}^{cl}$. Note that this implies
the initial square amplitudes of the classical field give the mode occupation
number expectations, i.e.

$\displaystyle|a^{cl}_{p}|^{2}\biggr{\rvert}_{t=0}=E[\hat{N}_{p}]$ (43)

We use a fourth order Runga-Kutta update scheme to update the field Runge
(1895); Kutta (1901). The update function for the field at mode $p$ is given

$\displaystyle F(a_{p})=a_{p}\left(1-i\,\omega_{p}\,\Delta
t\right)-i\,f(a)_{p}\,\Delta t\,,$ (44)

and the update scheme is then

  1. 1.

$k_{1}=F\left(a^{cl}_{p}(t)\right)$

  2. 2.

$k_{2}=F\left(a^{cl}_{p}(t)+k_{1}/2\right)$

  3. 3.

$k_{3}=F\left(a^{cl}_{p}(t)+k_{2}/2\right)$

  4. 4.

$k_{4}=F\left(a^{cl}_{p}(t)+k_{3}\right)$

  5. 5.

$a^{cl}_{p}(t+\Delta t)=\frac{1}{6}\left(k_{1}+2k_{2}+2k_{3}+k_{4}\right)$

which takes the field at a time $t$, $a^{cl}(t)$, to a time $t+\Delta t$,
$a^{cl}(t+\Delta t)$. The function $f(a)$ defines the potential term and is
given as follows

$\displaystyle f(a)_{p}=\mathcal{F}\left[V(x)\,\psi(x)\right]_{p}\,,$ (45)
$\displaystyle
V(x)=\mathcal{F}^{-1}\left[\mathcal{F}\left[\psi^{\dagger}(y)\psi(y)\right]_{i}\left(\frac{C}{k_{i}^{2}}+\Lambda_{0}\right)\right](x)\,.$

Where $\mathcal{F}$ and $\mathcal{F}^{-1}$ define the Fourier transform and
inverse Fourier transform respectively.
$\psi(x)=\sum_{i}a^{cl}_{i}u^{\dagger}_{i}(x)$, as in equation (8), i.e.
$\psi$ is the inverse Fourier transform of $a^{cl}$.

###  III.2 Field moment expansion

The evolution of the field moments is solved by integrating the coupled
equations (36d), (37), and (38). The initial conditions are chosen such that
the initial values of the moments correspond to those of a coherent state with
parameter $\vec{z}$. Meaning

$\displaystyle\braket{\hat{a}_{p}}^{\mathrm{FME}}\biggr{\rvert}_{t=0}$
$\displaystyle=z_{p}\,,$ (46)
$\displaystyle\braket{\delta\hat{a}_{i}\delta\hat{a}_{j}}^{\mathrm{FME}}\biggr{\rvert}_{t=0}$
$\displaystyle=0\,,$ (47)
$\displaystyle\braket{\delta\hat{a}^{\dagger}_{i}\delta\hat{a}_{j}}^{\mathrm{FME}}\biggr{\rvert}_{t=0}$
$\displaystyle=0\,.$ (48)

This solver uses the following update functions

$\displaystyle F^{\mathrm{FME}}_{a}(a_{p},\,\delta A_{ij},\,\delta B_{ij})$
$\displaystyle=a_{p}\left(1-i\,\omega_{p}\,\Delta t\right)-i\,\Delta
t\,\left(\,f(a^{1})_{p}+g_{1}(a,\,\delta A_{ij})_{p}+g_{2}(a,\,\delta
B_{ij})_{p}\,\right)\,,$ (49) $\displaystyle
F^{\mathrm{FME}}_{aa}(a_{p},\,\delta A_{ij},\,\delta B_{ij})$
$\displaystyle=\delta A_{ij}\left(1-i\,(\omega_{i}+\omega_{j})\Delta
t\right)-i\Delta t\left(h(a)_{ij}+g_{3}(a,\,\delta
A_{ij})_{ij}+g_{4}(a,\,\delta B_{ij})_{ij}\right)$ (50) $\displaystyle
F^{\mathrm{FME}}_{ba}(a_{p},\,\delta A_{ij},\,\delta B_{ij})$
$\displaystyle=\delta B_{ij}\left(1+i\,(\omega_{i}-\omega_{j})\Delta
t\right)-i\Delta t\left(g_{5}(a,\,\delta A_{ij})_{ij}+g_{6}(a,\,\delta
B_{ij})_{ij}\right)$ (51)

We use the following Runga-Kutta integration scheme to update the field
moments

  1. 1.

$k_{1}^{a}=F^{\mathrm{FME}}_{a}\left(\braket{\hat{a}_{p}}^{\mathrm{FME}}(t),\,\braket{\delta\hat{a}_{i}\delta\hat{a}_{j}}^{\mathrm{FME}}(t),\,\braket{\delta\hat{a}^{\dagger}_{i}\delta\hat{a}_{j}}^{\mathrm{FME}}(t)\right)$

  2. 2.

$k_{1}^{aa}=F^{\mathrm{FME}}_{aa}\left(\braket{\hat{a}_{p}}^{\mathrm{FME}}(t),\,\braket{\delta\hat{a}_{i}\delta\hat{a}_{j}}^{\mathrm{FME}}(t),\,\braket{\delta\hat{a}^{\dagger}_{i}\delta\hat{a}_{j}}^{\mathrm{FME}}(t)\right)$

  3. 3.

$k_{1}^{ba}=F^{\mathrm{FME}}_{ba}\left(\braket{\hat{a}_{p}}^{\mathrm{FME}}(t),\,\braket{\delta\hat{a}_{i}\delta\hat{a}_{j}}^{\mathrm{FME}}(t),\,\braket{\delta\hat{a}^{\dagger}_{i}\delta\hat{a}_{j}}^{\mathrm{FME}}(t)\right)$

  4. 4.

$k_{2}^{a}=F^{\mathrm{FME}}_{a}\left(\braket{\hat{a}_{p}}^{\mathrm{FME}}(t)+k_{1}^{a}/2,\,\braket{\delta\hat{a}_{i}\delta\hat{a}_{j}}^{\mathrm{FME}}(t)+k_{1}^{aa}/2,\,\braket{\delta\hat{a}^{\dagger}_{i}\delta\hat{a}_{j}}^{\mathrm{FME}}(t)+k_{1}^{ba}/2\right)$

  5. 5.

$k_{2}^{aa}=F^{\mathrm{FME}}_{aa}\left(\braket{\hat{a}_{p}}^{\mathrm{FME}}(t)+k_{1}^{a}/2,\,\braket{\delta\hat{a}_{i}\delta\hat{a}_{j}}^{\mathrm{FME}}(t)+k_{1}^{aa}/2,\,\braket{\delta\hat{a}^{\dagger}_{i}\delta\hat{a}_{j}}^{\mathrm{FME}}(t)+k_{1}^{ba}/2\right)$

  6. 6.

$k_{2}^{ba}=F^{\mathrm{FME}}_{ba}\left(\braket{\hat{a}_{p}}^{\mathrm{FME}}(t)+k_{1}^{a}/2,\,\braket{\delta\hat{a}_{i}\delta\hat{a}_{j}}^{\mathrm{FME}}(t)+k_{1}^{aa}/2,\,\braket{\delta\hat{a}^{\dagger}_{i}\delta\hat{a}_{j}}^{\mathrm{FME}}(t)+k_{1}^{ba}/2\right)$

  7. 7.

$k_{3}^{a}=F^{\mathrm{FME}}_{a}\left(\braket{\hat{a}_{p}}^{\mathrm{FME}}(t)+k_{2}^{a}/2,\,\braket{\delta\hat{a}_{i}\delta\hat{a}_{j}}^{\mathrm{FME}}(t)+k_{2}^{aa}/2,\,\braket{\delta\hat{a}^{\dagger}_{i}\delta\hat{a}_{j}}^{\mathrm{FME}}(t)+k_{2}^{ba}/2\right)$

  8. 8.

$k_{3}^{aa}=F^{\mathrm{FME}}_{aa}\left(\braket{\hat{a}_{p}}^{\mathrm{FME}}(t)+k_{2}^{a}/2,\,\braket{\delta\hat{a}_{i}\delta\hat{a}_{j}}^{\mathrm{FME}}(t)+k_{2}^{aa}/2,\,\braket{\delta\hat{a}^{\dagger}_{i}\delta\hat{a}_{j}}^{\mathrm{FME}}(t)+k_{2}^{ba}/2\right)$

  9. 9.

$k_{3}^{ba}=F^{\mathrm{FME}}_{ba}\left(\braket{\hat{a}_{p}}^{\mathrm{FME}}(t)+k_{2}^{a}/2,\,\braket{\delta\hat{a}_{i}\delta\hat{a}_{j}}^{\mathrm{FME}}(t)+k_{2}^{aa}/2,\,\braket{\delta\hat{a}^{\dagger}_{i}\delta\hat{a}_{j}}^{\mathrm{FME}}(t)+k_{2}^{ba}/2\right)$

  10. 10.

$k_{4}^{a}=F^{\mathrm{FME}}_{a}\left(\braket{\hat{a}_{p}}^{\mathrm{FME}}(t)+k_{3}^{a},\,\braket{\delta\hat{a}_{i}\delta\hat{a}_{j}}^{\mathrm{FME}}(t)+k_{3}^{aa},\,\braket{\delta\hat{a}^{\dagger}_{i}\delta\hat{a}_{j}}^{\mathrm{FME}}(t)+k_{3}^{ba}\right)$

  11. 11.

$k_{4}^{aa}=F^{\mathrm{FME}}_{aa}\left(\braket{\hat{a}_{p}}^{\mathrm{FME}}(t)+k_{3}^{a},\,\braket{\delta\hat{a}_{i}\delta\hat{a}_{j}}^{\mathrm{FME}}(t)+k_{3}^{aa},\,\braket{\delta\hat{a}^{\dagger}_{i}\delta\hat{a}_{j}}^{\mathrm{FME}}(t)+k_{3}^{ba}\right)$

  12. 12.

$k_{4}^{ba}=F^{\mathrm{FME}}_{ba}\left(\braket{\hat{a}_{p}}^{\mathrm{FME}}(t)+k_{3}^{a},\,\braket{\delta\hat{a}_{i}\delta\hat{a}_{j}}^{\mathrm{FME}}(t)+k_{3}^{aa},\,\braket{\delta\hat{a}^{\dagger}_{i}\delta\hat{a}_{j}}^{\mathrm{FME}}(t)+k_{3}^{ba}\right)$

  13. 13.

$\braket{\hat{a}_{p}}^{\mathrm{FME}}(t+\Delta
t)=\frac{1}{6}\left(k^{a}_{1}+2k^{a}_{2}+2k^{a}_{3}+k^{a}_{4}\right)$

  14. 14.

$\braket{\delta\hat{a}_{i}\delta\hat{a}_{j}}^{\mathrm{FME}}(t+\Delta
t)=\frac{1}{6}\left(k^{aa}_{1}+2k^{aa}_{2}+2k^{aa}_{3}+k^{aa}_{4}\right)$

  15. 15.

$\braket{\delta\hat{a}_{i}^{\dagger}\delta\hat{a}_{j}}^{\mathrm{FME}}(t+\Delta
t)=\frac{1}{6}\left(k^{ba}_{1}+2k^{ba}_{2}+2k^{ba}_{3}+k^{ba}_{4}\right)$

which takes the field moments are time $t$,
$\braket{\hat{a}}^{\mathrm{FME}}(t)$,
$\braket{\delta\hat{a}_{i}\delta\hat{a}_{j}}^{\mathrm{FME}}(t)$,
$\braket{\delta\hat{a}^{\dagger}_{i}\delta\hat{a}_{j}}^{\mathrm{FME}}(t)$, to
a time $t+\Delta t$, $\braket{\hat{a}}^{\mathrm{FME}}(t+\Delta t)$,
$\braket{\delta\hat{a}_{i}\delta\hat{a}_{j}}^{\mathrm{FME}}(t+\Delta t)$,
$\braket{\delta\hat{a}^{\dagger}_{i}\delta\hat{a}_{j}}^{\mathrm{FME}}(t+\Delta
t)$.

Where the functions are given by

$\displaystyle g_{1}(a,\,\braket{\delta a_{i}\delta a_{j}})_{p}$
$\displaystyle=\mathcal{F}\left[\mathcal{F}^{-1}\left[\mathcal{F}_{y}\left[\braket{\delta\psi(x)\delta\psi(y)}\psi^{\dagger}(y)\right]_{i}\left(\frac{C}{k_{i}^{2}}+\Lambda_{0}\right)\right](x,x)\right]_{p}$
(52) $\displaystyle g_{2}(a,\,\braket{\delta a_{i}^{\dagger}\delta
a_{j}})_{p}$
$\displaystyle=\mathcal{F}\left[\mathcal{F}^{-1}\left[\mathcal{F}_{x}\left[\braket{\delta\psi^{\dagger}(x)\delta\psi(y)}\psi(x)\right]_{i}\left(\frac{C}{k_{i}^{2}}+\Lambda_{0}\right)\right](x,x)\right]_{p}$
$\displaystyle+\mathcal{F}\left[\mathcal{F}^{-1}\left[\mathcal{F}\left[\braket{\delta\psi^{\dagger}(x)\delta\psi(x)}\right]_{i}\left(\frac{C}{k_{i}^{2}}+\Lambda_{0}\right)\right](x)\,\psi(x)\right]_{p}$
(53) $\displaystyle h(a)_{ij}$
$\displaystyle=\mathcal{F}_{xy}\left[K(x,y)\,\psi(x)\,\psi(y)\right]_{ij}$
(54) $\displaystyle g_{3}(a,\,\braket{\delta a_{i}\delta a_{j}})_{ij}$
$\displaystyle=\mathcal{F}_{xy}\left[\mathcal{F}^{-1}\left[\mathcal{F}_{y}\left[\braket{\delta\psi(x)\delta\psi(y)}\psi^{\dagger}(x)\right]_{i}\left(\frac{C}{k_{i}^{2}}+\Lambda_{0}\right)\right](x,y)\psi(x)\right]_{ij}$
$\displaystyle+\mathcal{F}_{xy}\left[\mathcal{F}^{-1}\left[\mathcal{F}\left[|\psi(x)|^{2}\right]_{i}\left(\frac{C}{k_{i}^{2}}+\Lambda_{0}\right)\right](x)\,\braket{\delta\psi(x)\delta\psi(y)}\right]_{ij}+(i\leftrightarrow
j)$ (55) $\displaystyle g_{4}(a,\,\braket{\delta a^{\dagger}_{i}\delta
a_{j}})_{ij}$
$\displaystyle=\mathcal{F}_{xy}\left[\mathcal{F}^{-1}\left[\mathcal{F}_{x}\left[\braket{\delta\psi^{\dagger}(x)\delta\psi(y)}\psi(x)\right]_{i}\left(\frac{C}{k_{i}^{2}}+\Lambda_{0}\right)\right](x,y)\,\psi(x)\right]_{ij}+(i\leftrightarrow
j)$ (56) $\displaystyle g_{5}(a,\,\braket{\delta a_{i}\delta a_{j}})_{ij}$
$\displaystyle=\mathcal{F}_{xy}\left[\mathcal{F}^{-1}\left[\mathcal{F}_{x}\left[\braket{\delta\psi(x)\delta\psi(y)}\,\psi^{\dagger}(x)\right]_{i}\left(\frac{C}{k_{i}^{2}}+\Lambda_{0}\right)\right](x,y)\,\psi^{\dagger}(x)\right]_{ij}-(c.c.,\,i\leftrightarrow
j)$ (57) $\displaystyle g_{6}(a,\braket{\delta a^{\dagger}_{i}\delta
a_{j}})_{ij}$
$\displaystyle=\mathcal{F}_{xy}\left[\mathcal{F}^{-1}\left[\mathcal{F}_{x}\left[\braket{\delta\psi^{\dagger}(x)\delta\psi(y)}\,\psi(x)\right]_{i}\left(\frac{C}{k_{i}^{2}}+\Lambda_{0}\right)\right](x,y)\,\psi^{\dagger}(x)\right]_{ij}$
$\displaystyle+\mathcal{F}_{xy}\left[\mathcal{F}^{-1}\left[\mathcal{F}\left[|\psi(x)|^{2}\right]_{i}\left(\frac{C}{k_{i}^{2}}+\Lambda_{0}\right)\right](x)\,\braket{\delta\psi^{\dagger}(x)\delta\psi(y)}\right]_{ij}-(c.c,\,i\leftrightarrow
j)\,.$ (58)

The position space fields are related to the momentum space arguments by

$\displaystyle\psi(x)$ $\displaystyle=\sum_{i}a_{i}u^{\dagger}_{i}(x)\,,$ (59)
$\displaystyle\braket{\delta\psi^{\dagger}(x)\psi(y)}$
$\displaystyle=\sum_{ij}\braket{\delta\hat{a}^{\dagger}_{i}\delta\hat{a}_{j}}u_{i}(x)u_{j}^{\dagger}(y)\,,$
(60) $\displaystyle\braket{\delta\psi(x)\psi(y)}$
$\displaystyle=\sum_{ij}\braket{\delta\hat{a}_{i}\delta\hat{a}_{j}}u^{\dagger}_{i}(x)u_{j}^{\dagger}(y)\,.$
(61)

###  III.3 Quantum field theory

The evolution of the quantum system is solved by integrating Schrödinger’s
equation

$\partial_{t}\ket{\vec{z}}=-i\hat{H}\,\ket{\vec{z}}\,.$ (62)

We integrate this equation using the QIBS repository available at
https://github.com/andillio/QIBS.

##  IV Test problems

The purpose of this section is to compare the FME and MFT approximations of
quantum systems using solutions which can be evaluated exactly, either
analytically or numerically. In general, we will look to demonstrate that the
field moment expansion is successful based on the following criteria

  1. 1.

Provides a more accurate approximation of the expectation value of the field
operator, $\braket{\hat{a}}$ at least until the quantum breaktime as defined
in equation (30).

  2. 2.

An accurate approximation of when the quantum breaktime occurs.

Note that the expectation values of the occupation numbers can be found using
$Q$, $n_{tot}$, and $\braket{\hat{a}}$. Therefore, achieving the two goals
listed above also implies that the field moment expansion can accurately
approximate the expectation values of the occupation numbers.

###  IV.1 Kerr nonlinearity

Figure 2:  Here we show the evolution of mean field for a given more for each
method in each of our test problems. The Kerr nonlinearity is shown in the top
left. The repulsive and attractive contact interactions in the top and bottom
right respectively. Attractive long range interactions are shown in the bottom
left. Nonzero coupling constants and mode numbers are given in each panel. On
the vertical axis we plot the absolute value of the expectation of one of the
mode field operators as a fraction of the initial value. On the horizontal
axis we plot the time as a fraction of the quantum breaktime. Therefore, $t=1$
corresponds to the quantum brekatime in all plotted systems, shown in dashed
lighted gray. In all cases, the classical field solution, shown in dotted red,
has diverged from the exact quantum solution, shown in solid blue. However,
the field moment expansion solution, shown in dashed green, remains an
accurate approximation of the exact quantum solution at least until the
breaktime. Figure 3:  Here we show the fractional error and approximation of
$Q$ for each method in each of our test problems. In the top row we show the
error in each method’s approximation of the magnitude of the mean field as a
fraction of the exact value of the mean field. In the bottom row we show each
method’s approximation of $Q$. The Kerr nonlinearity is shown on the far left
column. The repulsive and attractive contact interactions in the middle left
and right columns respectively. Attractive long range interactions are shown
in the far right column. Nonzero coupling constants and mode numbers are given
above each column. On the horizontal axis we plot the time as a fraction of
the quantum breaktime. Therefore, $t=1$ corresponds to the quantum brekatime
in all plotted systems, shown in dashed lighted gray. In all cases, the
classical field, shown in dotted red, has a relatively large fractional error
by the breaktime. However, the field moment expansion solution, shown in
dashed green, remains an accurate approximation of the exact quantum solution,
shown in blue, until at least the breaktime. Likewise, the field moment
expansion accurately approximates $Q$ until the breaktime.

In this section we examine the Kerr nonlinearity which can be described by the
following Hamiltonian

$\displaystyle\hat{H}$
$\displaystyle=\omega\,\hat{a}^{\dagger}\hat{a}+\frac{\Lambda_{0}}{2}(\hat{a}^{\dagger}\hat{a})^{2}$
$\displaystyle=\left(\omega+\frac{\Lambda_{0}}{2}\right)\,\hat{a}^{\dagger}\hat{a}+\frac{\Lambda_{0}}{2}\hat{a}^{\dagger}\hat{a}^{\dagger}\hat{a}\hat{a}$
$\displaystyle=\omega_{0}\,\hat{a}^{\dagger}\hat{a}+\frac{\Lambda_{0}}{2}\hat{a}^{\dagger}\hat{a}^{\dagger}\hat{a}\hat{a}\,.$
(63)

This is a special case of our Hamiltonian in equation (1) with $C=0,\,M=1$.

This problem is interesting because it admits an exact solution, and so the
time scales on which it diverges from the classical solution can be found
analytically Yurke and Stoler (1986); Tanaś and Kielich (1983). Our initial
condition will be a coherent state, see equation (14). The exact wavefunction
can be given as a function of time as follows

$\ket{\phi(t)}=\exp\left[-\frac{|z|^{2}}{2}\right]\sum_{n}e^{-it(\Lambda_{0}n^{2}/2+\omega
n)}\frac{z^{n}}{\sqrt{n!}}\ket{n}\,.$ (64)

Figure 4:  Here we show the evolution of the central field moments, normalized
by $|z|^{3}$, over time. We see all three moments become $\sim\mathcal{O}(1)$
on a time scale set by the nonlinearity. The second moments becoming
relatively large by the quantum breaktime, shown in dashed light gray. The
moments growth is hierarchical, i.e. the second moments become large before
the third moment. Here we set $z=5$, $\Lambda=1\times 10^{-3}$, and
$\omega=1$.

Given this, it is straightforward to calculate the exact evolution of the
normally ordered central moments and the field expectation. The quantum
evolution of the expectation of the field is characterized by a decaying
amplitude which is not captured in the classical theory, as shown in the top
left panel of figure 2. Here we set $z=5$, $\Lambda=1\times 10^{-3}$, and
$\omega=1$.

The far left column of figure 3 shows the result of applying the field moment
expansion to this system. We can see that $Q$ effectively parameterizes the
time when the fractional error in the classical theory is no longer small.
Until this point the field moment expansion provides a solution with a much
lower fractional deviation. Likewise, until the breaktime the field moment
expansion estimate of $Q$ remains accurate. Therefore, the field moment
expansion provides both a more accurate solution until the breaktime, and
successfully provides an accurate estimation of this time.

The evolution central moments are shown in figure 4. As expected for an
initially coherent state the central moments all start out at $0$ and then
grow on a time scale set by the nonlinearity. We see also that the moments
grow hierarchically, with the second central moments, normalized by $|z|^{3}$,
becoming $\sim\mathcal{O}(1)$ faster than the third central moment. We see
both the initial accuracy of the mean field theory and the hierarchically
growth conditions are met for this system.

The field moment expansion remains more accurate, past the breaktime, however,
we can see in figure 4 that past this time the third moment begins to become
relevant. Since we have truncated our expansion at second order past this time
is where we expect our solver to fail. Therefore, even assuming hierarchical
growth, the field moment expansion is not reliable past the time when the
highest moment in its truncation becomes large.

###  IV.2 Contact interactions

Unlike the previous test problem, scalar field dark matter systems involve a
large number of modes. Therefore, it is prudent to test the accuracy of the
field moment expansion on a system with multiple modes. We select the system
given by the Hamiltonian in equation (1) with $M=5$. This Hamiltonian has been
used as a test problem in Hertzberg (2016); Sikivie and Todarello (2017);
Erken  _et al._ (2012) and so will serve as a good benchmark to test the field
moment expansion. It also contains a self-interaction term which is present in
many models of scalar field dark matter.

This system models a contact interaction and linear dispersion, which implies
the following

$\displaystyle\Lambda^{ij}_{pl}$
$\displaystyle=\Lambda_{0}\,\delta^{ij}_{pl}\,,$ (65)
$\displaystyle\omega_{j}$ $\displaystyle=j\,\omega_{0}$ (66)

Where $\Lambda_{0}<0$ defines an attractive interaction and $\Lambda_{0}>0$ a
repulsive one. We evolve a coherent state defined by parameter
$\vec{z}\in\mathbb{C}^{5}$. In order to test how the solution behaves with
scaled occupation number we will simulate a benchmark coherent state
$\ket{\vec{z};r}$ where
$\vec{z}=(0,\,\sqrt{2r}\,e^{i\theta_{1}},\,\sqrt{2r}\,e^{i\theta_{2}},\,\sqrt{1r}\,e^{i\theta_{3}},\,0)$
and the phases are drawn from a uniform random distribution, $\theta_{i}\sim
U[0,2\pi)$ with fixed random seed.

Figure 5:  Here we plot the evolution the expectation of the number and field
operators. In the top row we show the expectation value of the occupation of
each mode. In the bottom row we show the amplitude of the expectation of the
field operator. In general this amplitude decays over time due to quantum
effects like phase diffusion. The middle and left columns show the results of
simulating contact interactions. The left column shows the repulsive case with
$\Lambda_{0}>0$ and the middle column the attractive case with
$\Lambda_{0}<0$. The right most column shows evolution of the long range
attractive interactions. For reference, the quantum breaktime for each system
is shown in dashed light gray. Nonzero nonlinear paramters are given above
each column. In all cases $\omega_{0}=1$, $M=5$, and $r=3$.

For this system, in general, the occupations of the modes will thermalize and
the expectation of the field itself will decay. This is shown in the right two
columns of figure 5. Here we set $\Lambda_{0}=\pm 0.1$, $\omega_{0}=1$, $M=5$,
and $r=3$.

We now test the field moment expansion to assess its accuracy and ability to
approximate the quantum break time. As in the previous problem we can see that
the field moment expansion solution remains close to the quantum solution past
where the deviation between the classical field theory and quantum field
theory becomes large. This is show for mode 3 in middle two columns of figure
2.

Figure 6: Here we show FME, in green, estimate of the breaktime compared to
the exact quantum result, in blue, for a number of $r$ values. We show both
the results of an analysis using $Q$ and the PO condition, dashed and solid,
respectively. We see that FME provides a close estimate of the breaktime. Here
we set $\omega_{0}=1$, $M=5$.

The field moment expansion also successfully predicts the quantum breaktime.
This is shown for $r=3$ in the middle two columns of figure 3. We can see that
the field moment expansion approximation of the field as well as the $Q$
parameter remains accurate until past the breaktime. We use the field moment
expansion to estimate the breaktime for a number of different values of $r$.
This is shown for the repulsive and attractive potentials in the top and
middle panels, respectively, of figure 6. There we also show an approximation
of the breaktime using the PO condition. The results of the two breaktime
definitions approximately agree. We see that the field moment expansion
closely approximates the breaktime in all cases.

###  IV.3 Long range interactions

SFDM can include self-interactions like those in the previous section.
However, given that we expect dark matter to be nearly collisionless, long
range interactions like those found in gravity are going to govern much of the
evolution. Therefore, we now turn towards modeling a system with a $1/r$
potential.

We still start with the Hamiltonian described in equation (1). Like in contact
interaction test problems we will use $M=5$. Long range interactions and a
quadratic dispersion relation can be modeled using the following constants

$\displaystyle\Lambda_{pl}^{ij}$
$\displaystyle=\left(\frac{C}{2(p_{p}-p_{i})^{2}}+\frac{C}{2(p_{p}-p_{j})^{2}}\right)\delta^{ij}_{pl}\,$
(67) $\displaystyle\omega_{j}$ $\displaystyle=\frac{j^{2}}{2}\omega_{0}\,.$
(68)

Where again $C<0$ gives an attractive potential, and $C>0$ a repulsive one.
Here we evolve the same benchmark coherent state, $\ket{\vec{z};r}$, as the
last section, where
$\vec{z}=(0,\,2r\,e^{i\theta_{1}},\,2r\,e^{i\theta_{2}},\,1r\,e^{i\theta_{3}},\,0)$
and the phases are drawn from a uniform random distribution, $\theta_{i}\sim
U[0,2\pi)$ with fixed random seed.

The density of dark matter is well measured Aghanim  _et al._ (2020). And
therefore the quantity $n_{tot}C$ should be fixed as we vary the occupation
number. This means that larger occupation results in a lower value of $C$. As
we scale our reference state then we will also scale the coupling constant
sending $C\rightarrow C/r$.

The evolution of this system is shown in the right column of figure 5 for
$r=3$. The field moment expansion produces an accurate estimation of the field
until the breaktime, see the right most column of figure 2. We can see that
the fractional error in the field moment expansion estimation of the field is
close to zero up to and past the breaktime. Additionally, the field moment
expansion successfully predicts $Q$ until the breaktime and consequently the
breaktime itself, shown in the right most column of figure 3.

##  V Conclusions

In all the test problems the FME successfully approximated the first and
second order moments of the exact quantum evolution when the correction terms
were subleading order, i.e. when $Q\not\sim 1$. Therefore, we can say that the
FME provides

  1. 1.

A more accurate approximation of the expectation value of the field operator,
$\braket{\hat{a}}$ at least until the quantum breaktime as defined in equation
(30).

  2. 2.

An accurate approximation of when the quantum breaktime occurs.

This is not terribly surprising because we solved the FME to second order.
Intuitively, a second order approximation should remain accurate for longer
than a first order approximation like MFT. Likewise, because the benchmarks of
classicality we used, the $Q$ parameter and the PO condition, are themselves
based on second order moments of the evolution, the fact that the FME can
approximate the breaktime is not surprising.

However, it is important to note that the method makes a number of assumptions
about the system. Specifically, we make an assumption about the initial
conditions and the evolution of the system as follows

  1. 1.

Initial conditions assumption: The initial conditions should be well
approximated by MFT. That is initially, $Q\ll 1$ and $1-\lambda_{p}/n_{tot}\ll
1$.

  2. 2.

Evolution assumption: The evolution of the field moments should be
hierarchical. That is, for a term in the equations of motion of the field
operators, $F$ of order $m$ written as a function of moments of order less
than and equal to $m$, given $F^{m}(\text{moments of order $\leq m$})$, over
the evolution the terms must satisfy $F^{1}>F^{2}>F^{3}\dots$.

For coherent state initial conditions and the Hamiltonian in equation (1),
these assumptions are generally satisfied. However, for initial conditions
such as number eigenstates, which do not have have hierarchically ordered
moment terms, neither FME or MFT will accurately approximate the quantum
evolution. The specific manner in which the quantum solutions approach a
classical description will be explored in a later paper.

If it is known that the initial conditions of a system are well described by
MFT, as in the case of coherent states, then the FME can provide a reliable
check of the timescales on which the MFT approximation remains valid. This
method could be applied to the evolution of ultra light scalar field dark
matter created via the misalignment mechanism. Because the misalignment
mechanism creates a system initially well described by the classical theory
FME could be used to approximate the quantum breaktime of this system. This
application will be explored in a later publication.

In this work we have focused on the results of solvers that integrate
equations expanded to second order in the field moments. The equations of
motion are obtained via a truncation of an infinite series of coupled
differential equations. However, it may be numerically feasible to instead
create a closure relation using Wick’s theorem, or as is done in Prezhdo and
Pereverzev (2000), which allows moments beyond second order to be approximated
using lower order moments. Such a solver may remain accurate for longer than
the current implementation.

###### Acknowledgements.

A.E., A.Z., and T.A. are supported by the U.S. Department of Energy under
contract number DE-AC02-76SF00515. .

## References

  * Anderson  _et al._ (1995) M. H. Anderson, J. R. Ensher, M. R. Matthews, C. E. Wieman, and E. A. Cornell, Science 269, 198 (1995).
  * Davis  _et al._ (1995) K. B. Davis, M.-O. Mewes, M. R. Andrews, N. J. Van Druten, D. S. Durfee, D. M. Kurn, and . Ketterle,  _Bose-Einstein Condensation in a Gas of Sodium Atoms_ , Tech. Rep. (1995).
  * Bialynicki-Birula (1977) I. Bialynicki-Birula, “Classical limit of quantum electrodynamics. [review],” (1977).
  * Hu  _et al._ (2000) W. Hu, R. Barkana, and A. Gruzinov, Phys. Rev. Lett. 85, 1158 (2000).
  * Mocz  _et al._ (2019) P. Mocz, A. Fialkov, M. Vogelsberger, F. Becerra, M. A. Amin, S. Bose, M. Boylan-Kolchin, P.-H. Chavanis, L. Hernquist, L. Lancaster, F. Marinacci, V. H. Robles, and J. Zavala, Phys. Rev. Lett. 123, 141301 (2019).
  * Guth  _et al._ (2015) A. H. Guth, M. P. Hertzberg, and C. Prescod-Weinstein, Phys. Rev. D 92, 103513 (2015).
  * Hertzberg (2016) M. P. Hertzberg, Journal of Cosmology and Astroparticle Physics 2016, 037 (2016), arXiv:1609.01342 [hep-ph] .
  * Leggett (2001a) A. J. Leggett, Reviews of Modern Physics 73, 307 (2001a).
  * Carr  _et al._ (2000) L. D. Carr, C. W. Clark, and W. P. Reinhardt, Physical Review A - Atomic, Molecular, and Optical Physics 62, 063610 (2000), arXiv:9911177 [cond-mat] .
  * Baer (2000) R. Baer, Physical Review A - Atomic, Molecular, and Optical Physics 62, 063810 (2000).
  * Minguzzi  _et al._ (2004) A. Minguzzi, S. Succi, F. Toschi, M. P. Tosi, and P. Vignolo, “Numerical methods for atomic quantum gases with applications to Bose-Einstein condensates and to ultracold fermions,”  (2004).
  * Sikivie and Todarello (2017) P. Sikivie and E. M. Todarello, Phys. Lett. B 770, 331 (2017), arXiv:1607.00949 [hep-ph] .
  * Alon  _et al._ (2007) O. Alon, A. Streltsov, and L. Cederbaum, Phys. Rev. A 77 (2007), 10.1103/PhysRevA.77.033613.
  * Prezhdo and Pereverzev (2000) O. V. Prezhdo and Y. V. Pereverzev, The Journal of Chemical Physics 113, 6557 (2000), https://doi.org/10.1063/1.1290288 .
  * Boßmann  _et al._ (2020) L. Boßmann, N. Pavlović, P. Pickl, and A. Soffer, Journal of Statistical Physics 178, 1362 (2020), arXiv:1905.06164 [math-ph] .
  * Heller (1976) E. Heller, Journal of Chemical Physics 64, 63 (1976).
  * Royer (1992) A. Royer, Foundations of Physics 22, 727 (1992).
  * Sreedharan  _et al._ (2020) A. Sreedharan, S. Choudhury, R. Mukherjee, A. Streltsov, and S. Wüster, Phys. Rev. A 101, 043604 (2020), arXiv:1904.11878 [cond-mat.quant-gas] .
  * Lentz  _et al._ (2019) E. W. Lentz, T. R. Quinn, and L. J. Rosenberg, Monthly Notices of the Royal Astronomical Society 485, 1809 (2019), arXiv:1810.09226 .
  * Gross (1961) E. P. Gross, Il Nuovo Cimento 20, 454 (1961).
  * Pitaevskii (1960) L. P. Pitaevskii, J. Exptl. Theoret. Phys. 13, 646 (1960).
  * Leggett (2001b) A. J. Leggett, Rev. Mod. Phys. 73, 307 (2001b).
  * Yurke and Stoler (1986) B. Yurke and D. Stoler, Phys. Rev. Lett. 57, 13 (1986).
  * Lewenstein and You (1996) M. Lewenstein and L. You, Phys. Rev. Lett. 77, 3489 (1996).
  * Caballero-Benitez  _et al._ (2008) S. Caballero-Benitez, E. Ostrovskaya, M. Gulacsi, and Y. Kivshar, (2008).
  * Chakrabarty  _et al._ (2018) S. S. Chakrabarty, S. Enomoto, Y. Han, P. Sikivie, and E. M. Todarello, Phys. Rev. D 97, 043531 (2018), arXiv:1710.02195 [hep-ph] .
  * Chakrabarty (2021) S. S. Chakrabarty, “Density perturbations in axion-like particles: classical vs quantum field treatment,” (2021), arXiv:2105.09749 [hep-ph] .
  * Kopp  _et al._ (2021) M. Kopp, V. Fragkos, and I. Pikovski, “Nonclassicality of axion-like dark matter through gravitational self-interactions,” (2021), arXiv:2105.13451 [astro-ph.CO] .
  * Dvali and Zell (2018) G. Dvali and S. Zell, Journal of Cosmology and Astroparticle Physics 2018, 064 (2018).
  * Dvali  _et al._ (2017) G. Dvali, C. Gomez, and S. Zell, JCAP 06, 028 (2017), arXiv:1701.08776 [hep-th] .
  * Lentz  _et al._ (2020) E. W. Lentz, T. R. Quinn, and L. J. Rosenberg, Nuclear Physics B 952, 114937 (2020).
  * (32) L. E. Ballentine,  _Quantum Mechanics A Modern Developement_.
  * Georges and Yedidia (1991) A. Georges and J. S. Yedidia, Journal of Physics A: Mathematical and General 24, 2173 (1991).
  * Glimm  _et al._ (1976) J. Glimm, A. Jaffe, and T. Spencer, Annals of Physics 101, 610 (1976).
  * Tanaś and Kielich (1983) R. Tanaś and S. Kielich, Optics Communications 45, 351 (1983).
  * Dvali  _et al._ (2013) G. Dvali, D. Flassig, C. Gomez, A. Pritzel, and N. Wintergerst, Phys. Rev. D 88, 124041 (2013), arXiv:1307.3458 [hep-th] .
  * Erken  _et al._ (2012) O. Erken, P. Sikivie, H. Tam, and Q. Yang, Phys. Rev. D 85, 063520 (2012).
  * Suárez and Chavanis (2015) A. Suárez and P.-H. Chavanis, Phys. Rev. D 92, 023510 (2015), arXiv:1503.07437 [gr-qc] .
  * Abbott and Sikivie (1983) L. Abbott and P. Sikivie, Physics Letters B 120, 133 (1983).
  * Preskill  _et al._ (1983) J. Preskill, M. B. Wise, and F. Wilczek, Phys. Lett. B 120, 127 (1983).
  * Penrose and Onsager (1956) O. Penrose and L. Onsager, Phys. Rev. 104, 576 (1956).
  * Runge (1895) C. Runge, Mathematische Annalen 46, 167 (1895).
  * Kutta (1901) W. Kutta, Zeit. Math. Phys. 46, 435 (1901).
  * Aghanim  _et al._ (2020) N. Aghanim, Y. Akrami, M. Ashdown, J. Aumont, C. Baccigalupi, M. Ballardini, A. J. Banday, R. B. Barreiro, N. Bartolo, S. Basak, R. Battye, K. Benabed, J. P. Bernard, M. Bersanelli, P. Bielewicz, J. J. Bock, J. R. Bond, J. Borrill, F. R. Bouchet, F. Boulanger, M. Bucher, C. Burigana, R. C. Butler, E. Calabrese, J. F. Cardoso, J. Carron, A. Challinor, H. C. Chiang, J. Chluba, L. P. Colombo, C. Combet, D. Contreras, B. P. Crill, F. Cuttaia, P. De Bernardis, G. De Zotti, J. Delabrouille, J. M. Delouis, E. Di Valentino, J. M. Diego, O. Doré, M. Douspis, A. Ducout, X. Dupac, S. Dusini, G. Efstathiou, F. Elsner, T. A. Enßlin, H. K. Eriksen, Y. Fantaye, M. Farhang, J. Fergusson, R. Fernandez-Cobos, F. Finelli, F. Forastieri, M. Frailis, A. A. Fraisse, E. Franceschi, A. Frolov, S. Galeotta, S. Galli, K. Ganga, R. T. Génova-Santos, M. Gerbino, T. Ghosh, J. González-Nuevo, K. M. Górski, S. Gratton, A. Gruppuso, J. E. Gudmundsson, J. Hamann, W. Handley, F. K. Hansen, D. Herranz, S. R. Hildebrandt, E. Hivon, Z. Huang, A. H. Jaffe, W. C. Jones, A. Karakci, E. Keihänen, R. Keskitalo, K. Kiiveri, J. Kim, T. S. Kisner, L. Knox, N. Krachmalnicoff, M. Kunz, H. Kurki-Suonio, G. Lagache, J. M. Lamarre, A. Lasenby, M. Lattanzi, C. R. Lawrence, M. Le Jeune, P. Lemos, J. Lesgourgues, F. Levrier, A. Lewis, M. Liguori, P. B. Lilje, M. Lilley, V. Lindholm, M. López-Caniego, P. M. Lubin, Y. Z. Ma, J. F. Maciás-Pérez, G. Maggio, D. Maino, N. Mandolesi, A. Mangilli, A. Marcos-Caballero, M. Maris, P. G. Martin, M. Martinelli, E. Martínez-González, S. Matarrese, N. Mauri, J. D. McEwen, P. R. Meinhold, A. Melchiorri, A. Mennella, M. Migliaccio, M. Millea, S. Mitra, M. A. Miville-Deschênes, D. Molinari, L. Montier, G. Morgante, A. Moss, P. Natoli, H. U. Nørgaard-Nielsen, L. Pagano, D. Paoletti, B. Partridge, G. Patanchon, H. V. Peiris, F. Perrotta, V. Pettorino, F. Piacentini, L. Polastri, G. Polenta, J. L. Puget, J. P. Rachen, M. Reinecke, M. Remazeilles, A. Renzi, G. Rocha, C. Rosset, G. Roudier, J. A. Rubiño-Martín, B. Ruiz-Granados, L. Salvati, M. Sandri, M. Savelainen, D. Scott, E. P. Shellard, C. Sirignano, G. Sirri, L. D. Spencer, R. Sunyaev, A. S. Suur-Uski, J. A. Tauber, D. Tavagnacco, M. Tenti, L. Toffolatti, M. Tomasi, T. Trombetti, L. Valenziano, J. Valiviita, B. Van Tent, L. Vibert, P. Vielva, F. Villa, N. Vittorio, B. D. Wandelt, I. K. Wehus, M. White, S. D. White, A. Zacchei, and A. Zonca, Astronomy and Astrophysics 641 (2020), 10.1051/0004-6361/201833910, arXiv:1807.06209 .