# Nonclassicality of axion-like dark matter through gravitational self-
interactions

Michael Kopp<EMAIL_ADDRESS>Department of Physics, Stockholm University,
SE-106 91 Stockholm, Sweden  Nordita, KTH Royal Institute of Technology and
Stockholm University, Hannes Alfvéns väg 12, SE-106 91 Stockholm, Sweden
Vasileios Fragkos<EMAIL_ADDRESS>Department of Physics,
Stockholm University, SE-106 91 Stockholm, Sweden  Igor Pikovski
<EMAIL_ADDRESS>Department of Physics, Stockholm University, SE-106
91 Stockholm, Sweden  Department of Physics, Stevens Institute of Technology,
Hoboken, NJ 07030, USA

###### Abstract

Axion-like particles (ALPs) are promising dark matter candidates. They are
typically described by a classical field, motivated by large phase space
occupation numbers. Here we show that such a description is accompanied by a
quantum effect: squeezing due to gravitational self-interactions. For a
typical QCD axion today, the onset of squeezing is reached on $\mathrm{\mu
s}$-scales and grows over millennia. Thus within the usual models based on the
classical Schrödinger-Poisson equation, a type of Gross-Pitaevskii equation,
any viable ALP is nonclassical. We also show that squeezing may be relevant on
the scales of other self-gravitating systems such as galactic haloes, or
solitonic cores. Conversely, our results highlight the incompleteness and
limitations of the classical single field description of ALPs.

##  I Introduction

The search for dark matter is one of the main challenges of modern physics. A
viable and popular dark matter candidate is a light (mass $m\ll 1\,\rm{eV}$)
self-gravitating quantum scalar field, such as the QCD axion or any other
axion-like particle (ALP). Extensive experimental searches focus on its
distinct wavelike signatures Schive  _et al._ (2014); Pozo  _et al._ (2020);
Arvanitaki  _et al._ (2020), and its theoretically expected coupling to matter
Marsh (2016); Sikivie (2021); Backes  _et al._ (2021). This type of ALP dark
matter, also including ultra-light or “fuzzy dark matter,” is usually
described by a classical scalar field rather than particles Abbott and Sikivie
(1983); Preskill  _et al._ (1983); Dine and Fischler (1983); Widrow and Kaiser
(1993); Peebles (2000); Hu  _et al._ (2000); Guth  _et al._ (2015); Hertzberg
(2016); Hui  _et al._ (2017); Hui (2021). The classical field description is
motivated by the de Broglie wave length $\lambdabar_{\rm deB}=\hbar/(mv_{\rm
rms})$ being significantly larger than the inter-particle distance $d$, and
thus the mean phase space occupation number $(\lambdabar_{\rm deB}/d)^{3}$ is
large. If a misalignment mechanism produces the ALP then initially it would be
in a coherent state Abbott and Sikivie (1983); Preskill  _et al._ (1983), and
thus a classical field.

An open question is whether dark matter can exhibit observable quantum
features. The difference between quantum and classical dynamics can be subtle,
and there are several non-classicality measures in quantum information science
Goldberg  _et al._ (2020). One such non-classicality witness is the inability
to represent the system in terms of classical mixtures of coherent states. The
most prominent quantum effect of this form is squeezing Drummond and Ficek
(2013). In the context of quantum optics, squeezing of bosonic systems is a
well established benchmark for non-classicality that has been observed in a
range of different systems Tanaś and Kielich (1983); Breitenbach  _et al._
(1997); Dodonov (2002); Wollman  _et al._ (2015), with many applications in
quantum-enhanced metrology Gross  _et al._ (2010) and quantum information
processing Zhong  _et al._ (2020); Hosten  _et al._ (2016). The
electromagnetic case is also instructive to understand the difference between
the classical and quantum states of the ALP. The starting point of quantum
optics and our understanding of the difference between classical and quantum
coherence was the Hanburry-Brown-Twiss experiment Brown and Twiss (1954). It
focused on the correlation between intensities from a beam of light that is
received at two spatially separated detectors. While the observed correlations
can be described classically, Glauber Glauber (1963) and Sudrashan Sudarshan
(1963) used it to develop a quantum description of light that showed how one
could distinguish classical from quantum effects: photon number correlations
can exhibit classical or quantum features – the latter typically referred to
as anti-bunching. Squeezed light is a particular quantum state that may
exhibit such quantum signatures that can be verified experimentally. As we
will show in this work, ALPs may well be in such a state.

In this paper we consider the dynamics of a self-gravitating quantum field and
show that it undergoes rapid self-squeezing due to self-interactions, a
quantum feature that cannot be captured by a classical field description. Our
results hold in the standard framework used for describing a classical scalar
field, the Gross-Pitaevskii equation (GPE) in the form of the Schrödinger-
Poisson equation. The squeezing happens on $\mu s$-timescales and is
continuously reproduced. We show formation of squeezing in three different
physical scenarios:  on cosmological scales, galactic scales, and in solitonic
cores that form in a galactic centers. These scenarios have the common feature
that they are approximately self-gravitating. On even smaller scales, such as
haloscopes (ground based experiments searching for axions) Sikivie (1983) we
find that squeezing becomes negligible even if the galactic halo undergoes
strong squeezing. A QCD axion, for example, reaches $9\mathrm{dB}$ squeezing –
reduction of the vacuum noise for one quadrature (and simultaneous increase
for the opposite quadrature) by an order of magnitude – after less than
$100\mu s$ for the cosmological volume down to the solitonic core volume.
Thus, under standard assumptions about the dynamical description of the ALP,
we arrive at the result that an initial “classical” quantum state quickly
evolves into a non-classical squeezed state. In contrast to previous works
Hertzberg (2016); Dvali and Zell (2018), here we find an example of quantum
behavior that forms very rapidly. The results suggest that either ALPs exhibit
non-negligible quantum features or that the usually adopted classical
description by a single classical GP scalar field is incomplete.

##  II The ALP model

We model the axion-like particle using a non-relativistic scalar quantum field
$\hat{\psi}(\bm{x},t)$. In a cosmological setting it is embedded in a
homogeneously and isotropically expanding Friedmann-Robertson-Walker (FRW)
universe. The resulting Heisenberg equation of motion of
$\hat{\psi}(\bm{x},t)$ is, see for instance Kopp  _et al._ (2017),

$\displaystyle i\hbar\partial_{t}\hat{\psi}(\bm{x},t)$
$\displaystyle=-\frac{\hbar^{2}}{2mA(t)^{2}}\bm{\nabla}^{2}\hat{\psi}(\bm{x},t)+m\hat{\Phi}(\bm{x},t)\hat{\psi}(\bm{x},t)$
(1) $\displaystyle\bm{\nabla}^{2}\hat{\Phi}(\bm{x},t)$
$\displaystyle=\frac{4\pi
Gm}{A(t)}\Big{(}\hat{\psi}^{\dagger}(\bm{x},t)\hat{\psi}(\bm{x},t)-\overline{\hat{\psi}^{\dagger}(\bm{x},t)\hat{\psi}(\bm{x},t)}\Big{)}\,,$

where $A(t)$ is the scale factor, which appears due to our use of comoving
coordinates $\bm{x}$. We remove the homogeneous mode, denoted by the overbar,
from the source of the gravitational potential in the Poisson equation. Our
approach therefore does not quantise the homogeneous mode of the ALP. This
allows us to avoid an IR divergence and conceptual issues of such a Newtonian
quantum cosmology. The homogeneous mode obeys the Friedmann equation
$H(t)^{2}=\frac{8\pi G}{3}\rho_{0}\,A(t)^{-3}$. Here, $A(t)$ is the scale
factor of the universe normalised to $A(t_{0})=1$, where $t_{0}$ is present
age of the universe. $H=\dot{A}(t)/A(t)$ is the Hubble expansion rate, and
$\rho_{0}=mn_{0}$ the present-day mean energy density of the ALP, with $n_{0}$
the present ALP mean number density and $m$ the mass of the ALP. We set the
speed of light $c=1$. For simplicity, we neglect here other types of matter
and the cosmological constant.111Inclusion of a cosmological constant is
straightforward and would not affect our approach based on (1), the order of
magnitude of our results and our conclusions. On the other hand, inclusion of
additional degrees of freedom will affect the quantum coherence, as discussed
in section III. To be consistent with the gravitational potential being
entirely due to the ALP in (1) we focus on a purely ALP dominated background
cosmology. For our purposes, this is a sufficiently accurate description of
the late universe. While a cosmological constant could be trivially included
(and will not change the order of magnitude of our results), inclusion of
other types of clustering matter would require a quantum mechanical modelling
which is beyond the scope of this paper, but will likely play a crucial role
in the quantum-classical transition of the ALP.

###  II.1 Hartree ansatz

Our approach to solve the Heisenberg equation (1) is to impose the Hartree
ansatz for the wave function of the ALP quantum field, and – crucially – the
assumption that this ansatz remains dynamically valid, which guarantees the
validity of the GPE Parkins and Walls (1998); L. Pitaevskii (2003); Alon  _et
al._ (2007); Erdos  _et al._ (2010). By doing so, we treat the ALP in a
conservative way since the GPE remains valid, but at the same time obtain an
analytically tractable nonlinear quantum description. The large squeezing at
very short timescales that we find here is precisely due to guaranteeing the
validity of the GPE by assuming the Hartree ansatz. Thus our results show that
squeezing should arise whenever this ansatz is justified.

The Hartree ansatz for the ALP quantum state implies that only a single mode
with mode function $\psi(\bm{x},t)$ and operator $\hat{a}(t)$ is relevant such
that approximately

$\hat{\psi}(\bm{x},t)=\frac{\psi(\bm{x},t)}{\sqrt{N}}\hat{a}(t)\,,$ (2)

with all other modes $\widehat{\delta\psi}(\bm{x},t)$ contributing to
$\hat{\psi}(\bm{x},t)$ neglected. This means that negligence of these other
modes defines $\psi(\bm{x},t)$ and $\hat{a}(t)$ through (2). The normalisation
$\sqrt{N}$ is for convenience and implies that the resulting GPE, for quantum
states $|\Psi(t_{i})\rangle$ containing $N$ (or approximately $N$) particles,
will not depend on $N$ in the limit $N\rightarrow\infty$.222The GPE field is
normalised to $\int_{V}d^{3}\\!x|\psi(\bm{x},t)|^{2}=N$, so that the spatial
average is $\overline{|\psi|^{2}}=N/V\equiv n_{0}$. Here $V$ is the comoving
volume, assumed to be very large, roughly the size of the Hubble volume
$H^{-3}_{0}$. Examples of such Hartree states are “non-classical” Fock states
$|N\rangle=\frac{1}{\sqrt{N!}}(\hat{a}^{\dagger}(t_{i}))^{N}|0\rangle$, or
“classical” coherent states
$|\alpha\rangle=e^{-N/2}\sum^{\infty}_{n=0}\frac{\alpha^{n}}{\sqrt{n!}}|n\rangle$,
with $\alpha=\sqrt{N}$ so that
$\hat{\psi}(x,t_{0})|\alpha\rangle=\psi(x,t_{0})|\alpha\rangle$. The state
$|\alpha\rangle$ describes a classical field configuration $\psi(x,t_{0})$
with negligible quantum fluctuations
$\mathrm{Var}(\hat{\psi})^{1/2}/\psi=N^{-1/2}$.

###  II.2 System of Kerr and Gross-Pitaevski equations

Projecting the Heisenberg equation (1) onto the mode function, that is
$\frac{1}{\sqrt{N}\hbar}\int_{V}d^{3}\\!x(\psi^{*}(\bm{x},t)\,\mathrm{Eq.}(1))$,
we obtain the Kerr oscillator equation for $\hat{a}(t)$:

$\displaystyle i\partial_{t}\hat{a}(t)$
$\displaystyle=\omega(t)\hat{a}(t)+2\chi(t)\hat{a}^{\dagger}(t)\hat{a}(t)\hat{a}(t)$
(3a) with coefficients (see Appendix E for a more rigorous derivation),
$\displaystyle\omega(t)$ $\displaystyle\equiv\tilde{\omega}(t)-\mu(t)$ (3b)
$\displaystyle\mu(t)$
$\displaystyle\equiv\frac{1}{N}i\int_{V}d^{3}\\!x\psi^{*}(\bm{x},t)\partial_{t}\psi(\bm{x},t)$
(3c) $\displaystyle\tilde{\omega}(t)$
$\displaystyle\equiv\frac{1}{N}\frac{\hbar}{2mA^{2}(t)}\int_{V}d^{3}\\!x|\bm{\nabla}\psi(\bm{x},t)|^{2}$
(3d) $\displaystyle\chi(t)$
$\displaystyle\equiv\frac{1}{2N^{2}}\frac{m}{\hbar}\int_{V}d^{3}\\!x\,\Phi(\bm{x},t)|\psi(\bm{x},t)|^{2}$
(3e) $\displaystyle=-\frac{1}{2N}\omega(t)\,.$ (3f)

The last equality (3f) follows from the Gross-Pitaevskii equation (GPE) for
$\psi(\bm{x},t)$, which is identical in form to the Heisenberg equation (1)
with all the “hats” removed, see Appendix E for a derivation.

Note that the common approach in the ALP literature is to set
$\hat{\psi}(\bm{x},t)\rightarrow\psi(\bm{x},t)$, and thus
$\hat{a}(t)\rightarrow\sqrt{N}$. While this indeed solves (3), it is not
strictly a valid solution since it does not satisfy
$[\hat{a},\hat{a}^{\dagger}]=1$. It neglects the quantum noise which undergoes
strong squeezing and can therefore grow in time, as we will now show.

The relevant dynamics of $\hat{a}$ is equivalent to a Kerr-Hamiltonian which
is well-known in quantum optics Tanaś and Kielich (1983); Yurke and Stoler
(1988) and in the context of BECs Parkins and Walls (1998). An initially
coherent state $|\alpha=\sqrt{N}\rangle$ evolves unitarily according to the
Hamiltonian
$\hat{H}=\hbar\omega\hat{a}^{\dagger}\hat{a}+\hbar\chi\hat{a}^{\dagger}\hat{a}^{\dagger}\hat{a}\hat{a}\,,$
with time-independent parameters $\omega$ and $\chi$. The Heisenberg equation
of motion $i\hbar\dot{\hat{a}}=[\hat{a},\hat{H}]$ is equivalent to (3) with
time-independent coefficients and with solution
$\hat{a}(t)=e^{-it(\omega+2\hat{a}^{\dagger}(t_{i})\hat{a}(t_{i})\chi)}\hat{a}(t_{i})$.
As we show below, assuming these parameters to be constant is a good
approximation for the time-scales we consider. In the following we set
$t_{i}=0$ to simplify expressions.

###  II.3 Quadrature squeezing

Defining the quadrature operator
$\hat{X}_{\theta}(t)=\hat{a}(t)e^{-i\theta}+\hat{a}^{\dagger}(t)e^{i\theta}\,,$
where $\theta$ is the quadrature angle, squeezing is present at
$\theta=\theta_{-}(t)$ if the minimal variance
$V_{-}(t)\equiv\mathrm{Var}(\hat{X}_{\theta_{-}(t)}(t))<1$ is below the vacuum
level (see Appendix A for more details). One effect of the non-linearity is to
generate such squeezing, which is still described in the Gaussian
approximation and takes place on much shorter timescales than non-Gaussian
effects such as generation of superposition states Yurke and Stoler (1988). In
the limit of large $N$, the minimal variance is given explicitly in Appendix
B, eq. (21). We define the squeezing timescale as the time on which $V_{-}(t)$
crosses below $e^{-2}\simeq 0.135$, i.e. the vacuum noise is reduced by a
factor 0.135 (or about $9\,\mathrm{dB}$ squeezing). Or in terms of the
squeezing parameter

$r(t)\equiv-\tfrac{1}{2}\ln V_{-}(t)\,,$ (4)

the squeezing timescale is $r(t_{\rm sqz})=1$. This yields

$N|\chi|t_{\rm sqz}\simeq\tfrac{1}{2}\sinh(1)\simeq 0.6\,.$ (5)

For large $N$, $V_{-}(t)$ approaches $0$ very closely, with its mininum
corresponding to a maximum of $r(t)$ given

$r_{\rm
max}\simeq\ln\big{(}3^{-1/2}2^{5/6}N^{1/6}\big{)}\simeq\ln\big{(}N^{1/6}\big{)}\,,$
(6)

with $r(t_{\rm max})=r_{\rm max}$ and

$N|\chi|t_{\rm max}\simeq 2^{-5/3}N^{1/6}\simeq 0.3N^{1/6}\,.$ (7)

The squeezing angle defining the orientation of the squeezed quadrature can be
approximated by

$\theta_{-}(t)\simeq\mathrm{sgn}(\chi t)\frac{\pi}{4}-\frac{1}{2}\arctan(2\chi
Nt)-(\omega+2\chi N)t\,,$ (8)

which for $t>t_{\rm sqz}$ approximately co-rotates with the classical solution
$a_{c}(t)=\alpha e^{-it(\omega+2\chi N)}$ of the Kerr model which is a good
approximation to the mean $\langle\hat{a}(t)\rangle$. The conjugate quadrature
with $\theta_{+}=\theta_{-}+\pi/2$ is exactly anti-squeezed. The $N$-dependent
phase shift provided by the arc tangent approximately cancels for $t>t_{\rm
sqz}$ with the remaining offset at $t_{\rm max}$ vanishing as
$\mathrm{sgn}(\chi t_{\rm max})N^{-1/6}$, thus producing an approximately
amplitude-squeezed state.333Eq. (8) corrects previous results, eqs. 26-28 in
Tanaś  _et al._ (1991). We also note that our expressions for $r_{\rm max}$
and $t_{\rm max}$ are more accurate than those presented in eq. 8 of Bajer
_et al._ (2002).

###  II.4 Estimation of the Kerr model parameters

We can express the parameters (3f) that govern the dynamics and the squeezing
in terms of known cosmological parameters and the ALP mass. We estimate
$\tilde{\omega}(a)$ by treating $\psi(\bm{x},t)$ in perturbation theory (see
Appendix D), or by solving for it numerically using the full GPE. In both
cases one approximately finds

$\displaystyle\frac{\hbar}{m}\tilde{\omega}(A)$
$\displaystyle=K(A)\simeq\tfrac{1}{2}v^{2}_{\rm rms}=\tfrac{1}{2}A\,v^{2}_{\rm
rms}(A=1)\,,$ (9)

where $mK$ is the average ALP kinetic energy, and $v^{2}_{\rm rms}$ is the
density weighted average of the ALP velocity field squared, see Kopp  _et al._
(2017), or in the context of cold dark matter (CDM) eq. 106 in Hertzberg
(2014). Numerical simulations show $v^{2}_{\rm rms}(A=1)\simeq 10^{-6}$ in the
present-day universe. This is the result of nonlinear structure formation
through gravitational instability of $\psi$ from primordial density
perturbations. On scales larger than $\lambdabar_{\rm deB}=\hbar/(mv_{\rm
rms})$, $|\psi(\bm{x},t)|^{2}$ and the density of of cold dark matter are
virtually the same Widrow and Kaiser (1993); Schive  _et al._ (2014); Mocz
_et al._ (2018); Kopp  _et al._ (2017). Therefore estimating $K(t)$ can be
done with any conventional CDM N-body simulation or CDM perturbation theory.
The squeezing timescales are determined by $\chi$ which is related to the
average potential energy $W=\frac{\hbar}{m}N\,\chi$ Kopp  _et al._ (2017).
Using the Layzer-Irvine equation $\frac{d}{dA}[A(W(A)+K(A))]=-K(A)$, eq. 47 of
Kopp  _et al._ (2017), we get with $K(A)\propto A$ from (9),
$W=-\frac{3}{2}K$. Thus combined

$\omega(t)=3\tilde{\omega}(t)\,,\quad\chi(t)=-\frac{3}{2}\frac{\tilde{\omega}(t)}{N}\,,\quad\mu(t)=-2\tilde{\omega}(t)\,.$
(10)

We have therefore tied all parameters appearing in the Kerr oscillator to a
simple expression for $\tilde{\omega}(t)$, eq. (9), entirely fixed by known
cosmological parameters and the ALP mass $m$.

Figure 1: Squeezing generated through unitary time evolution of a dark matter
coherent state $|\alpha\rangle$ in the Hartree approximation visualised by the
Wigner function
$W(z,t)=(\pi\hbar)^{-1}\langle\alpha|\hat{D}(z)(-1)^{\hat{a}^{\dagger}\hat{a}}\hat{D}^{\dagger}(z)|\alpha\rangle$,
see Cahill and Glauber (1969), where
$\hat{D}(z)=e^{z\hat{a}^{\dagger}-z^{*}\hat{a}}$. Darker hues of blue
correspond to larger values of $W$. Top panel: Initial coherent state. Middle
panel: Onset of squeezing. Bottom panel: Maximal achievable squeezing for
$\sqrt{N}=2000$. The initial time $t=0$ is an arbitrary moment in the late
universe (i.e. close to today). The red (large) dot indicates the classical
solution $a_{\rm cl}(t)=\alpha=2000$, and the black (small) dot the quantum
mean field $\langle\hat{a}(t)\rangle$. The classical solution is time
independent since (3f) holds for the ALP.

The number $N$ of axions in the volume $V$, with $V^{1/3}\lesssim H^{-1}_{0}$
is approximately

$N\simeq
H^{-3}_{0}\frac{\rho_{0}}{m}=10^{93}\frac{10^{-5}\,\mathrm{eV}}{m}\,.$ (11)

Note that this number deviates from ALP mode occupation numbers used elsewhere
in the literature, e.g. $10^{61}$ in Sikivie and Yang (2009) and $10^{26}$ in
Guth  _et al._ (2015) which refer to phase space occupation numbers in the
early and late universe, respectively. Our $N$ is the occupation number of the
single GP mode.

###  II.5 Onset of squeezing

We find the squeezing timescale of the ALP in the present-day universe from
(5), (9) and (10) to be

$t_{\rm sqz}\simeq\frac{1}{\tilde{\omega}}\simeq\frac{\hbar}{mv^{2}_{\rm
rms}}=\frac{\lambdabar_{\rm deB}}{v_{\rm rms}}\stackrel{{\scriptstyle
A=1}}{{\simeq}}\frac{10^{-5}\,\mathrm{eV}}{m}\,66\mu\mathrm{s}\,.$ (12)

We have thus shown that squeezing is rapidly generated for an ALP, much faster
than the age of the universe $\simeq H_{0}^{-1}\simeq 10^{18}\,\mathrm{s}$ for
all $m>10^{-22}\,\mathrm{eV}$. In other words, the single mode treatment of
the ALP predicts quantum effects for the full mass range of viable ALP models
on cosmologically short timescales. This also justifies our approximation of
the Kerr model parameters by the their present-day values, since
$|\dot{\omega}/\omega|\simeq|\dot{\chi}/\chi|\simeq H$. Note that the
(quantum) squeezing time scale $\lambdabar_{\rm deB}/v_{\rm rms}$ of the
operator $\hat{a}(t)$ essentially coincides with the (classical) coherence
time scale of the GP field $\psi(\bm{x},t)$, see for instance Graham and
Rajendran (2013).

For a QCD axion with $m\simeq 10^{-5}\,\mathrm{eV}$, this timescale is about
$66\mathrm{\mu s}$. One can compare this to other quantum mechanical
timescales appearing in the Kerr model, such as the Ehrenfest timescale
$t_{\rm Ehr}=N^{1/2}t_{\rm sqz}$, or Schrödinger “cat” creation time $t_{\rm
cat}=Nt_{\rm sqz}$ Yurke and Stoler (1988). Both of these timescales far
exceed the age of the universe for the ALP mass range.

|  Cosmology |  Solitonic core |  Milkyway
---|---|---|---
$t_{\mathrm{sqz}}[\mathrm{\mu s}]$ |  $66\,\bigl{(}\frac{10^{-5}\mathrm{eV}}{m}\bigr{)}$ |  $1400\,\bigl{(}\frac{10^{-5}\mathrm{eV}}{m}\bigr{)}$ |  $33\,\bigl{(}\frac{10^{-5}\mathrm{eV}}{m}\bigr{)}$
$t_{\mathrm{max}}$ [yr] |  $3500\,(\frac{10^{-5}\mathrm{eV}}{m})^{\frac{7}{6}}$ |  $0.5\,\bigl{(}\frac{10^{-5}\mathrm{eV}}{m}\bigr{)}^{\frac{4}{3}}$ |  $33000\,(\frac{10^{-5}\mathrm{eV}}{m})^{\frac{7}{6}}$
$r_{\mathrm{max}}$ |  $36+\frac{1}{6}\ln{\bigl{(}\frac{10^{-5}\mathrm{eV}}{m}\bigr{)}}$ |  $24+\frac{1}{3}\ln{\bigl{(}\frac{10^{-5}\mathrm{eV}}{m}\bigr{)}}$ |  $32+\frac{1}{6}\ln{\bigl{(}\frac{10^{-5}\mathrm{eV}}{m}\bigr{)}}$
$t_{\mathrm{Ehr}}$ [yr] |  $10^{35}(\frac{10^{-5}\mathrm{eV}}{m})^{\frac{3}{2}}$ |  $10^{21}\bigl{(}\frac{10^{-5}\mathrm{eV}}{m}\bigr{)}^{2}$ |  $10^{32}(\frac{10^{-5}\mathrm{eV}}{m})^{\frac{3}{2}}$
Table 1: Timescales and size of squeezing due to gravitational self-
interactions of the ALP, for various physical scenarios. $t_{\mathrm{sqz}}$
refers to the timescale on which $9\mathrm{dB}$ squeezing ($r=1$) is formed,
while $t_{\mathrm{max}}$ is the time for which the maximal squeezing
$r_{\mathrm{max}}$ is reached. In comparison, $t_{\mathrm{Ehr}}$ is the
Ehrenfest timescale which is often considered as a benchmark for the breakdown
of the classical field description. Our results show that quantum effects
arise on much shorter timescales.

###  II.6 Maximal squeezing

The maximum squeezing that can be produced in ALPs is obtained by inserting
(11) into (6):

$r_{\rm max}\simeq 36+\frac{1}{6}\ln\frac{10^{-5}\,\mathrm{eV}}{m}\,.$ (13)

This immense squeezing is similar in size to the squeezing of inflaton
perturbations produced during inflation, which is known to lead to observably
large quantum signatures Grishchuk and Sidorov (1990); Albrecht  _et al._
(1994); Polarski and Starobinsky (1996); Nelson and Jess Riedel (2017). The
time $t_{\rm max}$, eq. (7), at which maximal squeezing is reached is

$t_{\rm max}=0.5\,t_{\rm
sqz}N^{1/6}\simeq\Big{(}\frac{10^{-5}\,\mathrm{eV}}{m}\Big{)}^{7/6}\,3500\,\mathrm{yr}\,.$
(14)

Therefore even the extreme squeezing $r_{\rm max}$ would be easily reached in
the present universe, justifying our use of constant Kerr model parameters,
for $m\geq 10^{-12}\mathrm{eV}$.

Since $\omega+2N\chi=0$ the squeezing angle evolves slowly between
$\theta_{-}\simeq-\frac{\pi}{4}$ and $\theta_{-}\simeq-\frac{\pi}{4}N^{-1/6}$
in the interval $0<t\lesssim t_{\rm max}$. This implies that there is an
approximately fixed quadrature for which squeezing remains strongest and grows
as $r(t)$ over an extended period $t_{\rm max}$.

As an example we show in Fig. 1 the evolution of an initial coherent state
with $\alpha=2000$ into a squeezed coherent state using the Wigner
representation. Over time, the quadrature squeezing grows as the quadrature
rotates, while the number fluctuation remains constant (see the Appendices B
and C for more details).

We also consider two other physical situations with drastically smaller
volumes for which the GPE dynamics, or equivalently the Hartree ansatz, may be
better justified: galactic haloes and solitonic cores in dark matter haloes.
The calculations are analogous to the cosmological case, but without the scale
factor in eqs. (1). The details are presented in the Appendices H and I, and
the results are summarized in Table I. Importantly, the onset of squeezing is
independent of $N$, thus similar timescales $t_{\mathrm{sqz}}$ are found. In
Appendix I we also show that squeezing of axions contained in a volume that
isn’t gravitationally bound, such as an axion haloscope. In this scenario the
haloscope volume is part of a larger gravitationally bound volume, the
galactic halo which is modelled using the Hartree state ansatz adopted in this
paper, and exhibit negligible squeezing.

The simple approximation used here is common in the study of ALPs; it assumes
a single GP mode (the Hartree ansatz which assures the validity of the GPE).
This GP mode may well approximate some physical situations such as the
solitonic core.

##  III Discussion

Some important conclusions can be drawn from the results based on this simple
model. One is that the quantum effect described here accompanies the dynamics
whenever the GP assumption is established through the common Hartree ansatz.
Thus even for large occupancies a single classical field description that is
commonly assumed in ALP cosmology Widrow and Kaiser (1993); Peebles (2000); Hu
_et al._ (2000); Guth  _et al._ (2015); Hertzberg (2016); Hui  _et al._
(2017); Hui (2021) is incomplete. Our results therefore also highlight the
need to scrutinize the range of validity of the GP-ansatz in cosmology, which
is at the core of many predictions of ALP behavior such as interference
fringes Schive  _et al._ (2014). Another conclusion is that the quantum effect
we describe would be continuously regenerated on very short timescales even if
the ALP state were to collapse onto a coherent state through environmental
decoherence.

While a similar non-relativistic QFT description for ALPs was previously used
Hertzberg (2016); Erken  _et al._ (2012); Banik  _et al._ (2015); Sikivie and
Todarello (2017); Chakrabarty  _et al._ (2018), the specific quantum effect we
isolate here is novel in the context of ALPs and takes place on much shorter
timescales than other expected quantum phenomena. This squeezing timescale is
related to some previous findings. It matches the “classical break time”
introduced in Dvali and Zell (2018), which characterizes the onset of non-
linearities in the classical description. It can also be related to the
thermalisation timescale found in Sikivie and Yang (2009); Erken  _et al._
(2012), where a plane wave mode expansion was used to study thermalisation of
macroscopically occupied modes. The relaxation of these modes forms a Bose-
Einstein condensate (BEC) with an extension of the entire Hubble patch Erken
_et al._ (2012). This thermalisation timescale can be related to our squeezing
timescale as $Ht_{\rm th}\simeq\frac{1}{Ht_{\rm sqz}}$. Thus mode
thermalisation and self-squeezing are efficient at distinct periods with
thermalisation happening before self-squeezing.

It is instructive to compare our results to quantum effects in laboratory
BECs. Quantum revivals due to the Kerr effect have been demonstrated Greiner
_et al._ (2002), as well as spin squeezing of internal states Gross (2012).
The ground state of a trapped BEC has been shown to be a squeezed displaced
state Dunningham  _et al._ (1998), as found here, see the Appendix G. However,
in lab experiments the trapping typically dominates self-interactions and the
free-fall time is very short. Thus dynamical generation of self-squeezing as
described here has not been observed yet, but related theory for cold atoms
has been developed Johnsson and Haine (2007); Wüster  _et al._ (2008); Haine
and Johnsson (2009) and experiments with self-squeezed BECs might become
feasible in the near future.

From a purely theoretical perspective, our results show that quantum effects
are present in ALPs even in the limit where one typically assumes a purely
classical description – the GP equation from the Hartree ansatz is in fact
accompanied by a quantum effect that cannot be avoided. But experimentally, it
would be a challenging task to show non-classicality, more formidable than
verification of ALPs in the first place. To understand this better, let us
consider the verification of squeezing of light which has been developed
theoretically and experimentally in the 80’s using homodyne detection Yuen and
Chan (1983); Abas  _et al._ (1983), which measures the squeezed quadratures
directly but which also requires control of a reference coherent state.
Axions, if they exist, would be very difficult to control due to their feeble
interactions such that a homodyne detection scheme and other forms of quantum
state tomography, see e.g. Lvovsky and Raymer (2009), that require control of
a reference quantum state seem unrealisable for the axion. Thus the proposals
to detect squeezing in condensed cold atoms mentioned in the previous
paragraph, e.g. Johnsson and Haine (2007); Wüster  _et al._ (2008); Haine and
Johnsson (2009), might not be applicable to axions. Another idea is to attempt
to transfer the squeezing of the axion field into a squeezing of the
haloscope’s cavity mode and then infer the squeezing of that electromagnetic
mode using the methods for detecting squeezed light. This approach might not
be able to detect the Kerr-squeezing of the Hartree mode because the squeezing
of the axion field on scales of the haloscope volume is severely reduced, see
Appendix I. Thus, the indirect tests through interferometeric correlation
functions as mentioned in the Introduction might be better suited for
squeezing verification in axions. This could be achieved by measuring the
intensity correlation function
$G_{I}^{(2)}(t-t^{\prime})=\langle\hat{a}^{\dagger}(t)\hat{a}^{\dagger}(t^{\prime})\hat{a}(t)\hat{a}(t^{\prime})\rangle$
Walls (1983); Leuchs (1986). At $t=t^{\prime}$ squeezed states of light can
have $0<G_{I}^{(2)}/\langle\hat{a}^{\dagger}\hat{a}\rangle<1$ so that detected
photons are “anti-bunched”, whereas for a coherent state this is equal to 1 as
a consequence of the Poissonian statistics. Any other classical phase space
distribution of the electromagnetic radiation field can only exceed this
value, which is why squeezed light is called non-classical. In the case of
Kerr-squeezing discussed in this paper, however, the squeezing anti-bunching
does not occur, see the dashed line in Fig. 3. This is related to the fact
that the fundamental dynamics preserves the number state. It is thus expected
that detecting the squeezing of the axion, given a detection of the axion in
the first place, will be a difficult task and will require new approaches or
the study of indirect consequences on other fields.

A multimode treatment will not necessarily inhibit squeezing, which is an
inherent consequence of the self-interactions. For instance number-squeezing
was experimentally observed in a multi-well-trapped BEC Orzel  _et al._
(2001), and in the context of scalar field dark matter quadrature squeezing of
the GP mode has been observed to persist in a toy model where the Hartree
ansatz was replaced by five Fourier modes Eberhardt  _et al._ (2021).
Multiple localisation sites, analogs of ALP DM haloes, are known to lead to
the emergence of several macroscopically occupied GP-like fields Cederbaum and
Streltsov (2003). A systematic approach that reveals the necessity for and
determines the dynamics of additional mode functions has been established in
Alon  _et al._ (2007, 2008), and applied to what could be considered analog
situations of ALP cosmology Katsimiga  _et al._ (2017); Sreedharan  _et al._
(2020). These latter studies revealed that collision events between BEC
solitons composed of attractively interacting atoms can necessitate a
dynamical increase of required c-number fields, even if the initial state was
accurately described by a single GP field. Large squeezing is also prone to
decoherence when additional interactions are taken into account Zurek (2001);
Nelson and Jess Riedel (2017). For example, adding baryons would a priori
prevent the pure state Hartree ansatz (2) for the ALP that is used throughout
this paper. Other extensions of relevance will be the addition of unresolved
modes of the ALP or even more exotic coupling mechanisms such as time-dilation
induced effects Pikovski  _et al._ (2015). These will add to the decoherence
of the system and likely limit the attainable squeezing. But since it is
constantly regenerated we expect squeezing to persist at least on some time-
scales. The multi-mode treatment, inclusion of environments and resulting
decoherence will be considered in subsequent work.

Note added: After submission of this article a related work appeared Kuß and
Marsh (2021), in which two-mode squeezing of the axion is discussed. The
mechanism and the type of squeezing studied in Kuß and Marsh (2021) is
different from that discussed here. In our mechanism squeezing is due to self-
gravity, i.e. the nonlinear self-interaction of inhomogeneities of the scalar
field in the late universe, whereas in Kuß and Marsh (2021) squeezing of the
inhomogeneous modes of the axion is a linear process due to the strong time-
dependence of the expansion of space in the very early inflationary universe
Grishchuk and Sidorov (1990); Albrecht  _et al._ (1994); Polarski and
Starobinsky (1996); Nelson and Jess Riedel (2017).

##  IV Conclusion

In this paper we have shown quantum squeezing in ALP models of dark matter,
using the Hartree approximation which guarantees the validity of the Gross-
Pitaevskii equation. Our results highlight the quantum nature of the models,
even when they are expected to yield purely classical results. Non-classical
squeezed states are continuously formed at very short timescales in different
physical scenarios. On the one hand, our results motivate searches for
observable signatures of quantum effects of ALP dark matter. On the other
hand, the results challenge the validity of the simple models that are
routinely employed in ALP cosmology.

We thank Tom Abel, Alexander Balatsky, Andrew Eberhardt, Benjo Fraser, David
Marsh and Frank Wilczek for helpful discussions. This work was supported by
the Swedish Research Council under grant no. 2019-05615, the European Research
Council under Grant No. 742104 and by the Branco Weiss Fellowship – Society in
Science.

## References

  * Schive  _et al._ (2014) H.-Y. Schive, T. Chiueh, and T. Broadhurst, Nature Physics 10, 496 (2014), arXiv:1406.6586 [astro-ph.GA] .
  * Pozo  _et al._ (2020) A. Pozo, T. Broadhurst, I. de Martino, T. Chiueh, G. F. Smoot, S. Bonoli, and R. Angulo, arXiv e-prints , arXiv:2010.10337 (2020), arXiv:2010.10337 [astro-ph.GA] .
  * Arvanitaki  _et al._ (2020) A. Arvanitaki, S. Dimopoulos, M. Galanis, L. Lehner, J. O. Thompson, and K. Van Tilburg, Phys. Rev. D 101, 083014 (2020), arXiv:1909.11665 [astro-ph.CO] .
  * Marsh (2016) D. J. E. Marsh, Phys. Rep. 643, 1 (2016), arXiv:1510.07633 [astro-ph.CO] .
  * Sikivie (2021) P. Sikivie, Reviews of Modern Physics 93, 015004 (2021), arXiv:2003.02206 [hep-ph] .
  * Backes  _et al._ (2021) K. M. Backes, D. A. Palken, S. A. Kenany, B. M. Brubaker, S. B. Cahn, A. Droster, G. C. Hilton, S. Ghosh, H. Jackson, S. K. Lamoreaux, A. F. Leder, K. W. Lehnert, S. M. Lewis, M. Malnou, R. H. Maruyama, N. M. Rapidis, M. Simanovskaia, S. Singh, D. H. Speller, I. Urdinaran, L. R. Vale, E. C. van Assendelft, K. van Bibber, and H. Wang, Nature 590, 238 (2021), arXiv:2008.01853 [quant-ph] .
  * Abbott and Sikivie (1983) L. F. Abbott and P. Sikivie, Physics Letters B 120, 133 (1983).
  * Preskill  _et al._ (1983) J. Preskill, M. B. Wise, and F. Wilczek, Physics Letters B 120, 127 (1983).
  * Dine and Fischler (1983) M. Dine and W. Fischler, Physics Letters B 120, 137 (1983).
  * Widrow and Kaiser (1993) L. M. Widrow and N. Kaiser, ApJ 416, L71 (1993).
  * Peebles (2000) P. J. E. Peebles, ApJ 534, L127 (2000), arXiv:astro-ph/0002495 [astro-ph] .
  * Hu  _et al._ (2000) W. Hu, R. Barkana, and A. Gruzinov, Phys. Rev. Lett. 85, 1158 (2000), arXiv:astro-ph/0003365 [astro-ph] .
  * Guth  _et al._ (2015) A. H. Guth, M. P. Hertzberg, and C. Prescod-Weinstein, Phys. Rev. D 92, 103513 (2015), arXiv:1412.5930 [astro-ph.CO] .
  * Hertzberg (2016) M. P. Hertzberg, JCAP 2016, 037 (2016), arXiv:1609.01342 [hep-ph] .
  * Hui  _et al._ (2017) L. Hui, J. P. Ostriker, S. Tremaine, and E. Witten, Phys. Rev. D 95, 043541 (2017), arXiv:1610.08297 [astro-ph.CO] .
  * Hui (2021) L. Hui, arXiv e-prints , arXiv:2101.11735 (2021), arXiv:2101.11735 [astro-ph.CO] .
  * Goldberg  _et al._ (2020) A. Z. Goldberg, A. B. Klimov, M. Grassl, G. Leuchs, and L. L. Sánchez-Soto, AVS Quantum Science 2, 044701 (2020).
  * Drummond and Ficek (2013) P. D. Drummond and Z. Ficek,  _Quantum squeezing_ , Vol. 27 (Springer Science & Business Media, 2013).
  * Tanaś and Kielich (1983) R. Tanaś and S. Kielich, Optics Communications 45, 351 (1983).
  * Breitenbach  _et al._ (1997) G. Breitenbach, S. Schiller, and J. Mlynek, Nature 387, 471 (1997).
  * Dodonov (2002) V. V. Dodonov, Journal of Optics B: Quantum and Semiclassical Optics 4, R1 (2002).
  * Wollman  _et al._ (2015) E. E. Wollman, C. Lei, A. Weinstein, J. Suh, A. Kronwald, F. Marquardt, A. A. Clerk, and K. Schwab, Science 349, 952 (2015).
  * Gross  _et al._ (2010) C. Gross, T. Zibold, E. Nicklas, J. Estève, and M. K. Oberthaler, Nature 464, 1165 (2010), arXiv:1009.2374 [cond-mat.quant-gas] .
  * Zhong  _et al._ (2020) H.-S. Zhong, H. Wang, Y.-H. Deng, M.-C. Chen, L.-C. Peng, Y.-H. Luo, J. Qin, D. Wu, X. Ding, Y. Hu, _et al._ , Science 370, 1460 (2020).
  * Hosten  _et al._ (2016) O. Hosten, N. J. Engelsen, R. Krishnakumar, and M. A. Kasevich, Nature 529, 505 (2016).
  * Brown and Twiss (1954) R. H. Brown and R. Q. Twiss, The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 45, 663 (1954).
  * Glauber (1963) R. J. Glauber, Physical Review 131, 2766 (1963).
  * Sudarshan (1963) E. Sudarshan, Physical Review Letters 10, 277 (1963).
  * Sikivie (1983) P. Sikivie, Phys. Rev. Lett. 51, 1415 (1983).
  * Dvali and Zell (2018) G. Dvali and S. Zell, JCAP 2018, 064 (2018), arXiv:1710.00835 [hep-ph] .
  * Kopp  _et al._ (2017) M. Kopp, K. Vattis, and C. Skordis, Phys. Rev. D 96, 123532 (2017), arXiv:1711.00140 .
  * Parkins and Walls (1998) A. S. Parkins and D. F. Walls, Phys. Rep. 303, 1 (1998).
  * L. Pitaevskii (2003) S. S. L. Pitaevskii,  _Bose-Einstein condensation_, The International Series of Monographs on Physics (Oxford University Press, USA, 2003).
  * Alon  _et al._ (2007) O. E. Alon, A. I. Streltsov, and L. S. Cederbaum, Physics Letters A 362, 453 (2007), arXiv:cond-mat/0607490 [cond-mat.other] .
  * Erdos  _et al._ (2010) L. Erdos, B. Schlein, and H.-T. Yau, Annals of Mathematics 172, 291 (2010).
  * Yurke and Stoler (1988) B. Yurke and D. Stoler, Physica B+C 151, 298 (1988).
  * Tanaś  _et al._ (1991) R. Tanaś, A. Miranowicz, and S. Kielich, Phys. Rev. A 43, 4014 (1991).
  * Bajer  _et al._ (2002) J. Bajer, A. Miranowicz, and R. Tanaś, Czechoslovak Journal of Physics 52, 1313 (2002), arXiv:quant-ph/0212158 [quant-ph] .
  * Hertzberg (2014) M. P. Hertzberg, Phys. Rev. D 89, 043521 (2014), arXiv:1208.0839 [astro-ph.CO] .
  * Mocz  _et al._ (2018) P. Mocz, L. Lancaster, A. Fialkov, F. Becerra, and P.-H. Chavanis, Phys. Rev. D 97, 083519 (2018), arXiv:1801.03507 [astro-ph.CO] .
  * Cahill and Glauber (1969) K. E. Cahill and R. J. Glauber, Physical Review 177, 1882 (1969).
  * Sikivie and Yang (2009) P. Sikivie and Q. Yang, Phys. Rev. Lett. 103, 111301 (2009), arXiv:0901.1106 [hep-ph] .
  * Graham and Rajendran (2013) P. W. Graham and S. Rajendran, Phys. Rev. D 88, 035023 (2013), arXiv:1306.6088 [hep-ph] .
  * Grishchuk and Sidorov (1990) L. P. Grishchuk and Y. V. Sidorov, Phys. Rev. D 42, 3413 (1990).
  * Albrecht  _et al._ (1994) A. Albrecht, P. Ferreira, M. Joyce, and T. Prokopec, Phys. Rev. D 50, 4807 (1994), arXiv:astro-ph/9303001 [astro-ph] .
  * Polarski and Starobinsky (1996) D. Polarski and A. A. Starobinsky, Classical and Quantum Gravity 13, 377 (1996), arXiv:gr-qc/9504030 [gr-qc] .
  * Nelson and Jess Riedel (2017) E. Nelson and C. Jess Riedel, arXiv e-prints , arXiv:1711.05719 (2017), arXiv:1711.05719 [quant-ph] .
  * Erken  _et al._ (2012) O. Erken, P. Sikivie, H. Tam, and Q. Yang, Phys. Rev. D 85, 063520 (2012), arXiv:1111.1157 [astro-ph.CO] .
  * Banik  _et al._ (2015) N. Banik, A. Christopherson, P. Sikivie, and E. M. Todarello, arXiv e-prints , arXiv:1509.02081 (2015), arXiv:1509.02081 [astro-ph.CO] .
  * Sikivie and Todarello (2017) P. Sikivie and E. M. Todarello, Physics Letters B 770, 331 (2017), arXiv:1607.00949 [hep-ph] .
  * 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] .
  * Greiner  _et al._ (2002) M. Greiner, O. Mandel, T. W. Hänsch, and I. Bloch, Nature 419, 51 (2002), arXiv:cond-mat/0207196 [cond-mat.soft] .
  * Gross (2012) C. Gross, Journal of Physics B: Atomic, Molecular and Optical Physics 45, 103001 (2012).
  * Dunningham  _et al._ (1998) J. A. Dunningham, M. J. Collett, and D. F. Walls, Physics Letters A 245, 49 (1998).
  * Johnsson and Haine (2007) M. T. Johnsson and S. A. Haine, Phys. Rev. Lett. 99, 010401 (2007), arXiv:quant-ph/0609225 [quant-ph] .
  * Wüster  _et al._ (2008) S. Wüster, B. Dabrowska-Wüster, S. M. Scott, J. Close, and C. Savage, Physical Review A 77, 023619 (2008).
  * Haine and Johnsson (2009) S. A. Haine and M. T. Johnsson, Phys. Rev. A 80, 023611 (2009), arXiv:0812.2087 [quant-ph] .
  * Yuen and Chan (1983) H. P. Yuen and V. W. S. Chan, Optics Letters 8, 177 (1983).
  * Abas  _et al._ (1983) G. L. Abas, V. W. S. Chan, and T. K. Yee, Optics Letters 8, 419 (1983).
  * Lvovsky and Raymer (2009) A. I. Lvovsky and M. G. Raymer, Reviews of Modern Physics 81, 299 (2009), arXiv:quant-ph/0511044 [quant-ph] .
  * Walls (1983) D. F. Walls, Nature 306, 141 (1983).
  * Leuchs (1986) G. Leuchs, in  _Frontiers of nonequilibrium statistical physics_ (Springer, 1986) pp. 329–360.
  * Orzel  _et al._ (2001) C. Orzel, A. K. Tuchman, M. L. Fenselau, M. Yasuda, and M. A. Kasevich, Science 291, 2386 (2001).
  * Eberhardt  _et al._ (2021) A. Eberhardt, A. Zamora, M. Kopp, and T. Abel, arXiv e-prints , arXiv:2111.00050 (2021), arXiv:2111.00050 [hep-ph] .
  * Cederbaum and Streltsov (2003) L. S. Cederbaum and A. I. Streltsov, Physics Letters A 318, 564 (2003), arXiv:cond-mat/0310697 [cond-mat] .
  * Alon  _et al._ (2008) O. E. Alon, A. I. Streltsov, and L. S. Cederbaum, Phys. Rev. A 77, 033613 (2008), arXiv:cond-mat/0703237 [cond-mat.other] .
  * Katsimiga  _et al._ (2017) G. C. Katsimiga, G. M. Koutentakis, S. I. Mistakidis, P. G. Kevrekidis, and P. Schmelcher, New Journal of Physics 19, 073004 (2017), arXiv:1612.09151 [quant-ph] .
  * 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] .
  * Zurek (2001) W. H. Zurek, Nature 412, 712 (2001), arXiv:quant-ph/0201118 [quant-ph] .
  * Pikovski  _et al._ (2015) I. Pikovski, M. Zych, F. Costa, and Č. Brukner, Nature Physics 11, 668 (2015).
  * Kuß and Marsh (2021) J. L. J. Kuß and D. J. E. Marsh, The Open Journal of Astrophysics 4, 11 (2021), arXiv:2106.03528 [astro-ph.CO] .
  * Schive  _et al._ (2014) H.-Y. Schive, M.-H. Liao, T.-P. Woo, S.-K. Wong, T. Chiueh, T. Broadhurst, and W.-Y. P. Hwang, Phys. Rev. Lett. 113, 261302 (2014).
  * Simon (2002) C. Simon, Phys. Rev. A 66, 052323 (2002), arXiv:quant-ph/0110114 [quant-ph] .
  * Lee (2018) J.-W. Lee, J. Korean Phys. Soc. 73, 1596 (2018), arXiv:1510.07968 [hep-th] .
  * Navarro  _et al._ (1996) J. F. Navarro, C. S. Frenk, and S. D. M. White, ApJ 462, 563 (1996), arXiv:astro-ph/9508025 [astro-ph] .
  * Asztalos  _et al._ (2010) S. J. Asztalos, G. Carosi, C. Hagmann, D. Kinion, K. van Bibber, M. Hotz, L. J. Rosenberg, G. Rybka, J. Hoskins, J. Hwang, P. Sikivie, D. B. Tanner, R. Bradley, and J. Clarke, Phys. Rev. Lett. 104, 041301 (2010).

##  Appendix A Quadrature variance and squeezing

A quadrature is defined as

$\hat{X}_{\theta}(t)=\hat{a}(t)e^{-i\theta}+\hat{a}^{\dagger}(t)e^{i\theta}\,,$
(15)

and has variance

$\displaystyle\mathrm{Var}(\hat{X}_{\theta}(t))$
$\displaystyle=1+2(\langle\hat{a}^{\dagger}(t)\hat{a}(t)\rangle-|\langle\hat{a}(t)\rangle|^{2})$
$\displaystyle\quad+e^{-2i\theta}\mathrm{Var}(\hat{a}(t))+e^{2i\theta}\mathrm{Var}(\hat{a}^{\dagger}(t))$
(16)

where
$\mathrm{Var}(\hat{O})\equiv\langle\hat{O}^{2}\rangle-\langle\hat{O}\rangle^{2}$.
The minimising angle $\theta_{-}$ is

$e^{2i\theta_{-}(t)}=-\sqrt{\frac{\mathrm{Var}(\hat{a}(t))}{\mathrm{Var}(\hat{a}^{\dagger}(t))}}$
(17)

leading to the minimal variance
$V_{-}(t)\equiv\mathrm{Var}(\hat{X}_{\theta_{-}}(t))$,

$V_{-}(t)=1+2(\langle\hat{a}^{\dagger}(t)\hat{a}(t)\rangle-|\langle\hat{a}(t)\rangle|^{2})-2|\mathrm{Var}(\hat{a}(t))|\,.$
(18)

Since $\mathrm{Var}(\hat{a}^{\dagger}(t))=(\mathrm{Var}(\hat{a}(t)))^{*}$ the
only necessary ingredients are $\langle\hat{a}(t)\rangle$,
$\langle\hat{a}^{2}(t)\rangle$ and
$\langle\hat{a}^{\dagger}(t)\hat{a}(t)\rangle$.

##  Appendix B Details of the Kerr model

It is convenient to write
$\hat{H}=\hbar(\omega-\chi)\hat{n}+\hbar\chi\hat{n}^{2}$ because the linear
and nonlinear parts of $\hat{H}$ commute in this form so that
$U=e^{-i\hat{H}t/\hbar}=e^{-i(\omega-\chi)t\hat{n}}e^{-i\chi t\hat{n}^{2}}$.
For a coherent state $|\alpha\rangle$ this implies
$\langle\hat{a}^{\dagger}\hat{a}\rangle=|\alpha|^{2}=N$. We assume again
w.l.o.g that $\alpha=\sqrt{N}$.

To evaluate $V_{-}$ and $\theta_{-}$ in the Kerr model, the only additional
quantities we need are $\langle\hat{a}\rangle$ and
$\langle\hat{a}^{2}\rangle$. We calculate these expectation values using the
result
$e^{\hat{A}}\hat{B}e^{-\hat{A}}=\sum^{\infty}_{k=0}\tfrac{1}{k!}[\hat{A},\hat{B}]_{k}$,
where $[\hat{A},\hat{B}]_{k}=[\hat{A},[\hat{A},\hat{B}]_{k-1}]$,
$[\hat{A},\hat{B}]_{0}=\hat{B}$. Assigning $\hat{A}=i\hat{H}t/\hbar$ and
$\hat{B}=\hat{a}$, or $\hat{B}=\hat{a}^{2}$, gives
$\hat{U}^{-1}\hat{a}\hat{U}$ and $\hat{U}^{-1}\hat{a}^{2}\hat{U}$ which can be
evaluated using $[\hat{a},\hat{a}^{\dagger}]=1$, results in

$\displaystyle\langle\hat{a}\rangle$
$\displaystyle=\langle\alpha|\hat{U}^{-1}\hat{a}\hat{U}|\alpha\rangle$
$\displaystyle=\alpha e^{-i\omega t}e^{-2N\sin^{2}(\chi t)}e^{-iN\sin(2\chi
t)}$ (19) $\displaystyle\langle\hat{a}^{2}\rangle$
$\displaystyle=\langle\alpha|\hat{U}^{-1}\hat{a}^{2}\hat{U}|\alpha\rangle$
$\displaystyle=\alpha^{2}e^{-i2(\omega+\chi)t}e^{-2N\sin^{2}(2\chi
t)}e^{-iN\sin(4\chi t)}$ (20)

Inserting this into (17) and (18) gives the squeezing angle $\theta_{-}(t)$
and size $V_{-}(t)$ for the Kerr oscillator. The resulting expressions are
difficult to analyse further analytically due to their non-algebraic
structure. In Bajer  _et al._ (2002) a very accurate algebraic approximation
for $V_{-}(t)$ has been obtained by substituting $N\rightarrow\tau/(\chi t)$
and performing a Taylor expansion in $t$ (with $\tau$ fixed) which gives

$\displaystyle V_{-}(t)$ $\displaystyle\simeq 1-4\tau
s+8\tau^{2}+\frac{8\tau^{3}(5+12\tau^{2})}{Ns}-\frac{16\tau^{4}}{N}$ (21)
$\displaystyle s$ $\displaystyle\equiv\sqrt{1+4\tau^{2}}$ $\displaystyle\tau$
$\displaystyle\equiv N|\chi|t\,.$

Note that $V_{-}(t)$ is independent of the sign of $\chi$. To obtain the
squeezing timescale we hold $\tau$ fixed while taking the limit
$N\rightarrow\infty$ in (21). Then we solve $-\tfrac{1}{2}\ln V_{-}(t_{\rm
sqz})=1$ for $t_{\rm sqz}$ to obtain (12). To obtain the time of maximal
squeezing we replace $|\chi|t\rightarrow\tilde{\tau}N^{-5/6}$ in (21), and
while holding $\tilde{\tau}$ fixed, we expand in $1/N$. The leading term is

$\displaystyle V_{-}(\tilde{\tau})\simeq
N^{-1/3}\frac{2^{9}\tilde{\tau}^{6}+1}{2^{4}\tilde{\tau}^{2}}\,.$ (22)

Taking a $\tilde{\tau}$ derivative, we find the minimum at $\tilde{\tau}_{\rm
max}=2^{-5/3}$, which gives (7). Inserting this back into
$V_{-}(\tilde{\tau}_{\rm max})$ corresponding to $r_{\rm
max}=-\frac{1}{2}\,\ln V_{-}(\tilde{\tau}_{\rm max})$ and results in (7). Fig.
2 shows $t_{\rm sqz}$, $t_{\rm max}$ and $r_{\rm max}$ as a function of $N$,
and compares our approximations (dot-dashed) for (5), (7) and (6) to their
exact result (full lines).

Finally we calculate $\theta_{-}$. We proceed in a similar fashion as we did
for $V_{-}$. We insert (19) and (20) into
$\mathrm{Var}(\hat{a})/(\mathrm{Var}(\hat{a}))^{*}$ and we replace
$N\rightarrow T/(\chi t)$. Then we take the log and Taylor expand the result
in $\chi t$ with the result

$\displaystyle\theta_{-}$ $\displaystyle\simeq-\frac{i}{4}\ln\Big{(}\
e^{-i4\omega t}\frac{e^{-8iT}(2T+i)}{2T-i}\Big{)}$ (23)
$\displaystyle=-\frac{i}{4}\ln\Big{(}\frac{2T+i}{2T-i}\Big{)}-\omega t-2T\,.$
(24)

This can be simplified into the expression (8). The signum function in (8)
arises due a branch cut in the log.

Figure 2: Squeezing timescales $t_{\rm sqz}$ and $t_{\rm max}$ (upper panel)
and maximal squeezing (lower panel), and their scaling with $N$, for the Kerr
model. Full lines are exact, dashed lines approximations used in the main
text.

##  Appendix C Squeezed coherent states

Displacement operator

$\hat{D}(\beta)=e^{\beta\hat{a}^{\dagger}-\beta^{*}\hat{a}}$ (25)

and squeeze operator

$\hat{S}(\zeta)=e^{\frac{1}{2}(\zeta^{*}\hat{a}^{2}-\zeta(\hat{a}^{\dagger})^{2})}$
(26)

have the properties

$\displaystyle\hat{D}^{\dagger}(\beta)\hat{a}\hat{D}(\beta)$
$\displaystyle=\hat{a}+\beta$ (27)
$\displaystyle\hat{S}^{\dagger}(\zeta)\hat{a}\hat{S}(\zeta)$
$\displaystyle=\hat{a}\cosh(\rho)-\hat{a}^{\dagger}e^{2i\phi}\sinh(\rho)\,,$
(28)

where

$\zeta=\rho e^{2i\phi}$ (29)

with $\rho\geq 0$. A squeezed coherent state is defined as

$|\beta,\zeta\rangle=\hat{D}(\beta)\hat{S}(\zeta)|0\rangle\,.$ (30)

Evaluating (18), (17) we find

$\displaystyle\langle\hat{a}\rangle$ $\displaystyle=\beta$ (31) $\displaystyle
V_{-}$ $\displaystyle=e^{-2\rho}$ (32) $\displaystyle\theta_{-}$
$\displaystyle=\phi\,.$ (33)

Since we know $\langle\hat{a}(t)\rangle$, $V_{-}(t)=:e^{-2r(t)}$ and
$\theta_{-}(t)$ for the Kerr oscillator, see (19), (21), (8), we can construct
an approximate solution to the Kerr oscillator quantum state in the
Schrödinger picture as $|\Psi^{\rm
SC}(t)\rangle\simeq|\langle\hat{a}(t)\rangle,r(t)e^{2i\theta_{-}(t)}\rangle$.
In the main text we have approximated $\langle\hat{a}\rangle\simeq a_{c}(t)$.
This is justified because we consider quantum evolution only up to the time
$t_{\rm max}\simeq 0.5N^{1/6}t_{\rm sqz}$ which is much smaller than the
Ehrenfest time $t_{\rm Ehr}=N^{1/2}t_{\rm sqz}$.

The squeezed coherent state approximation breaks down due to non-Gaussianities
around the time $t_{\rm max}$.

Figure 3: $\mathrm{Var}(\hat{n})$ for a squeezed coherent state $|\beta,\rho
e^{2i\phi}\rangle$ with $\beta=\sqrt{N}=10^{6}$ and $\rho=r(t)$ with $0\leq
r(t)\leq r_{\rm max}$ of the Kerr oscillator (21). Inserting $\phi=0$
corresponds to amplitude squeezing, accompanied by $\hat{n}$-squeezing,
whereas inserting $\phi=\theta_{-}(t)$ of the Kerr evolution prevents
$\hat{n}$-squeezing. The dashed line shows $\mathrm{Var}(\hat{n})$ for a
reference Kerr state.

In order to verify that a squeezed coherent state is a good approximation to
the Kerr-evolved coherent state until $t_{\rm max}$ we evaluate

$\mathrm{Var}(\hat{n})=\tfrac{1}{4}\big{(}\cosh(4\rho)-1\big{)}+N(\cosh(2\rho)-\cos(2\phi)\sinh(2\rho))$
(34)

for the squeezed coherent state $|\sqrt{N},r(t)e^{2i\phi}\rangle$ with
$\phi=\theta_{-}(t)$ and $\phi=0$ in Fig. 3. We see that inclusion of the
squeezing orientation is crucial to match the value $\mathrm{Var}(\hat{n})=N$
of the Kerr-evolved coherent state. This also explains the “coincidence” of
$\theta_{-}(t_{\rm max})\simeq\tfrac{\pi}{4}N^{-1/6}$ and $V^{1/2}_{+}\simeq
N^{1/6}$. As can be seen in the lower panel of Fig. 1, this combination of
angle and vertical extend of the phase space distribution leads to an order
unity range in radial direction that is independent of $N$.

##  Appendix D Perturbation theory estimate of $K(A)$

The goal of this Appendix is to show the physical origin, the approximate
order of magnitude, as well as the $A$-scaling of the kinetic energy $K(A)$.
This result was used derive a relation (10) between $\mu$ and $\tilde{\omega}$
appearing in the Kerr oscillator equation (3). A more precise estimate of
$K(A)$ could be extracted from a cosmological N-body or Schrödinger-Poisson
simulation if needed.

The mean kinetic energy $mK(A)$ and the mean potential energy $mW(A)$ of an
ALP can be evaluated in perturbation theory. Here we use linear perturbation
theory of cold dark matter modelled by a pressureless perfect fluid, with
particle density $n(\bm{x},t)$ and velocity field $\bm{u}(\bm{x},t)$
satisfying the Euler-Poisson equation. Since approximately
$\tfrac{\hbar^{2}}{m^{2}}|\bm{\nabla}\psi(\bm{x},t)|^{2}=n(\bm{x},t)|\bm{u}(\bm{x},t)|^{2}$,
see e.g. Widrow and Kaiser (1993), Kopp  _et al._ (2017), the mean kinetic
energy per mass is

$\displaystyle K(A)$ $\displaystyle=\frac{1}{2A^{2}N}\int_{V}d^{3}\\!x\
n(\bm{x},A)\bm{u}^{2}(\bm{x},A)$ (35)
$\displaystyle=\frac{1}{2N}\int_{V}d^{3}\\!x\ n(\bm{x},A)\bm{v}^{2}(\bm{x},A)$
(36)

where $\bm{v}=\bm{u}/A$ is the peculiar velocity field, and $\bm{u}$ is the
canonical velocity field. In linear perturbation theory the comoving density
field can be written as $n(\bm{x},A)=n_{0}(1+\delta(\bm{x},A))$ with the
density perturbation $\delta(\bm{x},A)$ related to the velocity field via
$\bm{\nabla}\cdot\bm{v}=-AH\delta(\bm{x},A)$. During the assumed matter
dominated expansion $\dot{A}/A=H=H_{0}A^{-3/2}$ we thus get to leading order
in perturbation theory

$\displaystyle K(A)$ $\displaystyle\simeq\frac{1}{2V}\int_{V}d^{3}\\!x\
\bm{v}^{2}(\bm{x},A)$ (37) $\displaystyle=\frac{A}{2V}\int_{V}d^{3}\\!x\
\bm{v}^{2}(\bm{x},A=1)\,.$ (38)

In Fourier space, taking the limit $V\rightarrow\infty$, this becomes

$\displaystyle K(A)$
$\displaystyle\simeq\frac{A}{2}H_{0}^{2}\int\frac{d^{3}k}{(2\pi)^{3}}\frac{|\delta_{\bm{k}}(A=1)|^{2}}{k^{2}}$
(39) $\displaystyle=\tfrac{1}{2}A\,v^{2}_{\rm rms}(A=1)$ (40) $\displaystyle
v^{2}_{\rm rms}(A=1)$
$\displaystyle:=H_{0}^{2}\int\frac{dk}{2\pi^{2}}P(k,A=1)\,,$ (41)

where we used the ergodic theorem to replace $|\delta_{\bm{k}}|^{2}\rightarrow
P(k=|\bm{k}|)$, and performed the angular integrals. For a standard matter
power spectrum $P(k)$ this evaluates to $10^{-5}$. The actual value of
$v^{2}_{\rm rms}$ used in the main text is about a ten times smaller due to
our negligence of the cosmological constant. Inclusion of a cosmological
constant in linear perturbation would give rise to an overall pre-factor
$f^{2}=0.27$ and slightly shallower scaling than $A$, namely
$K\propto(AHf)^{2}$. Here $f\simeq\Omega_{m}^{6/11}$ is the growth rate of
linear matter density perturbations,
$\Omega_{m}=\Omega_{m,0}A^{-3}H_{0}^{2}/H^{2}$ the fraction of matter and
$H=H_{0}(\Omega_{m,0}A^{-3}+1-\Omega_{m,0})^{1/2}$ the expansion rate.

##  Appendix E Derivation of the coupled GPE and Kerr oscillator

Throughout this section we work in the Schrödinger picture, so that the
quantum state $|\Psi(t)\rangle$ of the non-relativistic ALP is time dependent.
For simplicity we discard the cosmological evolution, so that the Hamiltonian
is

$\hat{\mathcal{H}}=\int
d^{3}\\!x\,\frac{\hbar^{2}}{2m}\nabla\hat{\psi}^{\dagger}(\bm{x})\nabla\hat{\psi}(\bm{x})-\\\
\frac{Gm^{2}}{2}\int
d^{3}\\!x\,d^{3}\\!x^{\prime}\,\frac{\hat{\psi}^{\dagger}(\bm{x})\hat{\psi}^{\dagger}(\bm{x}^{\prime})\hat{\psi}(\bm{x})\hat{\psi}(\bm{x}^{\prime})}{|\bm{x}-\bm{x}^{\prime}|}\,.$
(42)

$|\Psi(t)\rangle$ satisfies the Schrödinger equation

$i\hbar\partial_{t}|\Psi(t)\rangle=\hat{\mathcal{H}}|\Psi(t)\rangle\,.$ (43)

The Hartree ansatz for a $n$-particle wave function in the Schrödinger picture
is

$\Psi_{n}(\bm{x}_{1},...,\bm{x}_{n},t)=\phi(\bm{x}_{1},t)...\phi(\bm{x}_{n},t)\,.$
(44)

Here we introduced the normalised one-particle wave function $\phi(\bm{x},t)$.
This $n$-particle state can be written basis-independent as
$|n,t\rangle=\frac{1}{\sqrt{n!}}(\hat{b}^{\dagger}(t))^{n}|0\rangle$, where
$\hat{b}^{\dagger}(t)\equiv\int
d^{3}\\!x\,\phi(\bm{x},t)\hat{\psi}^{\dagger}(\bm{x})$ and
$\hat{\psi}^{\dagger}(\bm{x})$ is the time-independent field operator in the
Schrödinger picture. A generic (e.g. initially coherent) Hartree state
$|\Psi(t)\rangle$ then takes the form

$|\Psi(t)\rangle=\sum_{n=0}^{\infty}c_{n}(t)|n,t\rangle\,.$ (45)

Our goal is to derive evolution equations for $\phi(\bm{x},t)$ and $c_{n}(t)$.
Varying the action $S=\int
dt\langle\Psi(t)|\hat{\mathcal{H}}-i\hbar\partial_{t}|\Psi(t)\rangle$ w.r.t.
$\phi^{*}(\bm{x},t)$ and $c_{n}^{*}(t)$ gives

$\displaystyle i\hbar\partial_{t}\phi(\bm{x},t)$
$\displaystyle=-\frac{\hbar^{2}}{2m}\nabla^{2}\phi(\bm{x},t)+$ (46)
$\displaystyle\quad\
+\frac{\overline{n^{2}}(t)-\bar{n}(t)}{\bar{n}^{2}(t)}m\Phi(\bm{x},t)\phi(\bm{x},t)$
$\displaystyle\nabla^{2}\Phi(\bm{x},t)$ $\displaystyle\equiv 4\pi
Gm\bar{n}(t)|\phi(\bm{x},t)|^{2}$ (47) $\displaystyle
i\hbar\partial_{t}c_{n}(t)$ $\displaystyle=\langle
n,t|\hat{\mathcal{H}}-i\hbar\partial_{t}|n,t\rangle c_{n}(t)$ (48)
$\displaystyle\langle
n,t|\hbar^{-1}\hat{\mathcal{H}}-i\partial_{t}|n,t\rangle$
$\displaystyle=n(\omega(t)-\chi(t))+n^{2}\chi(t)$ (49)
$\displaystyle\bar{n}(t)$ $\displaystyle\equiv\sum_{n}|c_{n}(t)|^{2}n$ (50)
$\displaystyle\overline{n^{2}}(t)$
$\displaystyle\equiv\sum_{n}|c_{n}(t)|^{2}n^{2}$ (51) $\displaystyle\omega(t)$
$\displaystyle\equiv\tilde{\omega}(t)-\mu(t)$ (52)
$\displaystyle\tilde{\omega}(t)$ $\displaystyle\equiv\int
d^{3}\\!x\,\frac{\hbar}{2m}|\nabla\phi(\bm{x},t)|^{2}$ (53)
$\displaystyle\mu(t)$ $\displaystyle\equiv-i\int
d^{3}\\!x\,\phi^{*}(\bm{x},t)\partial_{t}\phi(\bm{x},t)$ (54)
$\displaystyle\chi(t)$ $\displaystyle\equiv\frac{m}{2\hbar\bar{n}(t)}\int
d^{3}\\!x\,\Phi(\bm{x},t)|\phi(\bm{x},t)|^{2}\,.$ (55)

Neglecting the time dependence of $\omega$ and $\chi$ we find

$c_{n}(t)=c_{n}(t\\!=\\!0)e^{-it((\omega-\chi)n+\chi n^{2})}\,,$ (56)

which is the solution of the Schrödinger equation for a Kerr oscillator. For
an initial coherent state $|\Psi(t=0)\rangle=|\alpha\rangle$ with
$\alpha=\sqrt{N}$ we have $c_{n}(t\\!=\\!0)=e^{-N/2}\frac{N^{n/2}}{\sqrt{n!}}$
so that $\bar{n}(t)=N$ and $\overline{n^{2}}(t)=N^{2}+N$. The $\phi$ evolution
equation, (46) and (47), then simplifies to

$\displaystyle i\hbar\partial_{t}\phi(\bm{x},t)$
$\displaystyle=-\frac{\hbar^{2}}{2m}\nabla^{2}\phi(\bm{x},t)+m\Phi(\bm{x},t)\phi(\bm{x},t)$
(57a) $\displaystyle\nabla^{2}\Phi(\bm{x},t)$ $\displaystyle\equiv 4\pi
GmN|\phi(\bm{x},t)|^{2}$ (57b)

which is the GPE used in the main text, after re-defining
$\psi\equiv\phi\sqrt{N}$ and restoring the scale factor $A(t)$. Equations
(56), (57) then establish the validity of the constraints $\int
d^{3}\\!x\,|\phi(\bm{x},t)|^{2}=1$ and $\sum_{n}|c_{n}(t)|^{2}=1$, which we
could have enforced via Lagrange multipliers Alon  _et al._ (2007). Noticing
that $\langle n,t|\hat{\mathcal{H}}-i\hbar\partial_{t}|n,t\rangle=\langle
n|\hat{H}|n\rangle$, with $|n\rangle=|n,t\\!=\\!0\rangle$ and
$\hat{H}=\hbar(\omega-\chi)\hat{n}+\hbar\chi\hat{n}^{2}$, we verify (3) for
$\hat{a}(t)=e^{i\hat{H}t/\hbar}\hat{b}(t=0)e^{-i\hat{H}t/\hbar}$ in the
Heisenberg picture.

In the main text we have worked in the Heisenberg picture where
$\hat{\psi}(\bm{x},t)=\phi(\bm{x},t)\hat{a}(t)$. As a check let us compare the
result of the mean field
$\langle\alpha|\hat{\psi}(\bm{x},t)|\alpha\rangle=\phi(\bm{x},t)\langle\alpha|\hat{a}(t)|\alpha\rangle$
in the Heisenberg picture to that in Schrödinger picture
$\langle\Psi(t)|\hat{\psi}(\bm{x})|\Psi(t)\rangle$. For this we assume
$\hat{\psi}(\bm{x})=\phi(\bm{x},t)\hat{b}(t)$. This seems somewhat odd since
$\hat{\psi}(\bm{x})$ should be time-independent, but is not. This is merely an
artifact of restricting the field operator to contain only the single GP mode,
and does not pose a problem. We thus have
$\langle\Psi(t)|\hat{\psi}(\bm{x})|\Psi(t)\rangle=\phi(\bm{x},t)\langle\Psi(t)|\hat{b}(t)|\Psi(t)\rangle$.
Let us focus on $\langle\Psi(t)|\hat{b}(t)|\Psi(t)\rangle$ and switch to the
“Kerr picture” that we define through $|\Psi_{\rm
K}(t)\rangle\equiv\sum_{n}c_{n}(t)|n,t\\!=\\!0\rangle$. Clearly,
$\langle\Psi(t)|\hat{b}(t)|\Psi(t)\rangle=\langle\Psi_{\rm
K}(t)|\hat{b}(t\\!=\\!0)|\Psi_{\rm K}(t)\rangle$. This equals to
$\langle\Psi_{\rm
K}(t\\!=\\!0)|e^{i\hat{H}t/\hbar}\hat{b}(t\\!=\\!0)e^{-i\hat{H}t/\hbar}|\Psi_{\rm
K}(t\\!=\\!0)\rangle=\langle\alpha|\hat{a}(t)|\alpha\rangle$, and we have
established the equivalence to the Heisenberg picture calculation. Similar
calculations establish the equivalence of all correlators and thus in
particular the squeezing.

In summary, using a variational ansatz for a quantum state in the Hartree
approximation we have derived the coupled GP and Kerr oscillator equations for
an initial coherent state of a non-relativistic ALP field.

##  Appendix F Squeezed state approximation

In the “Kerr picture” where $\hat{\psi}_{\rm
K}(\bm{x},t)=\psi(\bm{x},t)\tfrac{\hat{a}(t_{i})}{\sqrt{N}}$ and $|\Psi_{\rm
K}(t)\rangle=e^{-i\hat{H}t/\hbar}|\alpha\rangle$, the quantum state can be
approximated by a squeezed coherent state

$|\Psi^{\rm SC}_{\rm
K}(t)\rangle=\hat{D}\big{(}a_{c}(t)\big{)}\hat{S}\Big{(}r(t)e^{2i\theta_{-}(t)}\Big{)}|0\rangle\,,$
(58)

where
$\hat{D}(\beta)=e^{\beta\hat{a}^{\dagger}(t_{i})-\beta^{*}\hat{a}(t_{i})}$,
and
$\hat{S}(\zeta)=e^{\frac{1}{2}(\zeta^{*}\hat{a}^{2}(t_{i})-\zeta(\hat{a}^{\dagger}(t_{i}))^{2})}$,
are the displacement and squeezing operators, respectively Walls (1983).
Properties of this state are shown in Fig.1 and Fig. 3 (the curve labeled by
$\theta_{-}$).

A “classical” wave function $|\Psi^{\rm cl}_{\rm K}(t)\rangle$ would then be a
state with $\hat{a}(t_{i})|\Psi^{\rm cl}_{\rm K}(t)\rangle=a_{c}(t)|\Psi^{\rm
cl}(t)\rangle$ at all times, which is only true for a coherent state
$|\Psi^{\rm cl}_{\rm K}(t)\rangle=\hat{D}(a_{c}(t))|0\rangle$. The fidelity
$|\langle\Psi^{\rm cl}_{\rm K}(t)|\Psi^{\rm SC}_{\rm
K}(t)\rangle|^{2}=|\langle
0|\hat{S}\big{(}r(t)e^{2i\theta_{-}(t)}\big{)}|0\rangle|^{2}$ decays at the
squeezing time $t_{\rm sqz}$ indicating deviation from the classical
approximation.

##  Appendix G Ground state in a single-GP-mode approximation

If thermalisation as described in Erken  _et al._ (2012) indeed keeps the
axion at its momentary ground state $|\Psi_{0}(t)\rangle$, then it might be
natural to expect that this ground state is related to our squeezed state.
Erken  _et al._ (2012) does not explore the nature of this ground state.
However, the ground state of a self-interacting BEC in the single mode
approximation is discussed in Dunningham  _et al._ (1998); Parkins and Walls
(1998). To leading order in a perturbative solution
$\hat{a}=\alpha+\delta\hat{a}$, this ground state is found to be a squeezed
coherent state

$|\Psi_{0}(t)\rangle\simeq e^{-i\mu_{0}t}|\Psi(t_{\rm max})\rangle\,,$ (59)

matching our approximate solution (58) at $t_{\rm max}$ apart from some
irrelevant global phase. Squeezing in the Kerr model reaches a maximum only
because the quantum state ceases to be well approximated by a Gaussian after
$t_{\rm max}$. Interestingly, the relation to the BEC ground state might
persist even when non-Gaussianities are included: it was found in Dunningham
_et al._ (1998); Parkins and Walls (1998) that the ground state Wigner
function of the Kerr Hamiltonian is bent into a “banana” producing a so-called
number-squeezed state.

##  Appendix H Solitonic cores

For the solitonic core case, the relevant parameter to compute in order to
obtain the squeezing timescale is

$\displaystyle N\chi$
$\displaystyle\equiv\frac{1}{2N}\frac{m}{\hbar}\int_{V}d^{3}\\!x\,\Phi(\bm{x},t)|\psi(\bm{x},t)|^{2}\,,$
(60)

where we have assumed that $\chi$ is constant in time for the timescales we
are interested. Assuming a spherically symmetric dark matter halo and making
use of the radial core density profile in Schive  _et al._ (2014), it is
straightforward to find the Newtonian potential by solving the Poisson
equation. Inserting (60) into (5), we get

$t_{\mathrm{sqz}}\simeq 2\times
10^{8}\biggl{(}\frac{m}{10^{-22}\mathrm{eV}}\biggr{)}\biggl{(}\frac{r_{c}}{\mathrm{kpc}}\biggr{)}^{2}\mathrm{yr}\,,$
(61)

where $r_{\mathrm{c}}$ is the solitonic core radius. It is convenient to
express the squeezing timescale in terms of the host halo mass $M_{h}$. To
achieve that, we make use of the relation between the core radius $r_{c}$ and
the core mass $M_{c}$ Schive  _et al._ (2014)

$\frac{r_{c}}{\mathrm{kpc}}=\frac{5.5\times
10^{7}}{\bigl{(}\frac{m}{10^{-22}\mathrm{eV}}\bigr{)}^{2}}\frac{M_{\odot}}{M_{c}}$
(62)

and then express the core mass in terms of the host halo mass $M_{h}$ Schive
_et al._ (2014)

$M_{c}=\frac{1}{4\sqrt{\alpha}}\biggl{(}\frac{\zeta{(z)}}{\zeta{(0)}}\biggr{)}^{\frac{1}{6}}\frac{\bigr{(}4.4\times
10^{7}\bigl{)}^{\frac{2}{3}}}{\frac{m}{10^{-22}\mathrm{eV}}}M_{\odot}^{\frac{2}{3}}M_{\mathrm{h}}^{\frac{1}{3}}\,,$
(63)

where $z=A^{-1}-1$. For further details, the reader is referred to Schive  _et
al._ (2014). The squeezing timescale as a function of the host halo mass is
given by

$t_{\mathrm{sqz}}=7\times
10^{14}\biggl{(}\frac{10^{-22}\mathrm{eV}}{m}\biggr{)}\biggl{(}\frac{M_{\odot}}{M_{\mathrm{h}}}\biggr{)}^{\frac{2}{3}}\,\mathrm{yr}\,.$
(64)

For a dark matter halo of mass $M_{\mathrm{h}}=2\times 10^{12}M_{\odot}$ we
find

$t_{\mathrm{sqz}}\simeq
1400\biggl{(}\frac{10^{-5}\mathrm{eV}}{m}\biggr{)}\mathrm{\mu s}$ (65)

The number of axions inside the solitonic core is

$N=\frac{M_{c}}{m}\simeq 3\times
10^{62}\biggl{(}\frac{10^{-5}\mathrm{eV}}{m}\biggr{)}^{2}$ (66)

Having obtained the squeezing timescale and the total number of axions inside
a solitonic core, it is straightforward to compute $t_{\mathrm{max}}$,
$r_{\mathrm{max}}$ and $t_{\mathrm{Ehr}}$. The values can be found in the main
text (Table I).

##  Appendix I Axion haloscope and Milkyway

The Hartree, or single mode, ansatz can only be applied to sufficiently
isolated systems. In the present context of a purely dark matter dominated
universe this means that the volume to which we can apply the Hartree ansatz
must be to a good approximation self-gravitating. This was the case for the
entire universe and the solitonic core. The axions within a haloscope volume
are not self-gravitating, in the sense that their motion is dominated by the
external galactic gravitational potential and not their self-gravity. To
remain consistent with the Hartree ansatz, the premise of this work, we assume
that the whole galactic halo – plausibly the smallest self-gravitating volume
containing the haloscope – is described by a Hartree state under time
evolution and then focus on the squeezing within a subvolume comprising the
haloscope as the observable. This provides an estimate of the squeezing
timescale for axions within a haloscope.

The entire galactic halo contains $N$ axions and the assumed Hartree ansatz
implies that the only relevant operator is

$\hat{a}=\int d^{3}x\,\phi^{*}(\bm{x},t)\hat{\psi}(\bm{x},t)\,,$ (67)

where $\phi(\bm{x},t)=\psi(\bm{x},t)/\sqrt{N}$ is the 1-particle wave function
of a galactic axion. We now decompose $\phi(\bm{x},t)$ into two orthogonal
functions
$\phi_{h}:=\frac{\sqrt{N}}{\sqrt{N_{h}}}\,\phi(\bm{x},t)\theta_{h}(\bm{x})$
and
$\phi_{\bar{h}}:=\frac{\sqrt{N}}{\sqrt{N-N_{h}}}\,\phi(\bm{x},t)(1-\theta_{h}(\bm{x}))$,
where $\theta_{h}(\bm{x})$ equals 1 within the haloscope and vanishes outside
of it, and $N_{h}$ is the number of axions in the haloscope. We then have

$\hat{a}=\frac{\sqrt{N_{h}}}{\sqrt{N}}\hat{b}+\frac{\sqrt{N-N_{h}}}{\sqrt{N}}\hat{c}$
(68)

with
$[\hat{a},\hat{a}^{\dagger}]=[\hat{b},\hat{b}^{\dagger}]=[\hat{c},\hat{c}^{\dagger}]=1$,
$[\hat{b},\hat{c}^{\dagger}]=[\hat{c},\hat{b}^{\dagger}]=0$ and
$[\hat{b},\hat{a}^{\dagger}]=\frac{\sqrt{N_{h}}}{\sqrt{N}}$, and

$\displaystyle\hat{b}=$ $\displaystyle\int
d^{3}x\,\phi_{h}^{*}(\bm{x},t)\hat{\psi}(\bm{x},t)$ (69)
$\displaystyle\hat{c}=$ $\displaystyle\int
d^{3}x\,\phi_{\bar{h}}^{*}(\bm{x},t)\hat{\psi}(\bm{x},t)\,.$ (70)

A similar decomposition has been employed in Simon (2002); Lee (2018) to
investigate the entanglement between a subvolume and the remaining part of a
BEC. The squeezing of the haloscope mode $\hat{b}$ is given by

$V^{h}_{-}(t)=1+2(\langle\hat{b}^{\dagger}(t)\hat{b}(t)\rangle-|\langle\hat{b}(t)\rangle|^{2})-2|\mathrm{Var}(\hat{b}(t))|\,,$
(71)

where we assume as before that the initial quantum state (the quantum state in
the Heisenberg picture) of the galactic halo is in an $\hat{a}(t_{i})$-mode
squeezed coherent state $|\alpha\rangle$. We have shown that unitary evolution
can be approximated by
$\hat{U}(t)=\hat{D}(\sqrt{N})\hat{S}(r(t)e^{2i\theta_{-}(t)})\hat{D}^{\dagger}(\sqrt{N})$,
see eqs. (25), (26) and (58). Using the commutation relation between $\hat{a}$
and $\hat{b}$ operators we find using a calculation similar to that presented
in Appendix B

$\displaystyle\hat{D}^{\dagger}(\alpha)\hat{b}\hat{D}(\alpha)$
$\displaystyle=\hat{b}+\frac{\sqrt{N_{h}}}{\sqrt{N}}\alpha$ (72)
$\displaystyle\hat{S}^{\dagger}(\zeta)\hat{b}\hat{S}(\zeta)$
$\displaystyle=\frac{\sqrt{N_{h}}}{\sqrt{N}}\big{(}\hat{a}\cosh(r)-\hat{a}^{\dagger}e^{2i\theta_{-}}\sinh(r)\big{)}\,,$
(73)

with $\alpha=\sqrt{N}$ and $\zeta=re^{2i\theta_{-}}$. Using these results (71)
can be simplified most easily in the “Kerr picture”, see Appendix F, to give

$V^{h}_{-}(t)=1+\frac{N_{h}}{N}\Big{(}V_{-}(t)-1\Big{)}\,,$ (74)

where $V_{-}(t)=e^{-2r(t)}$ is the minimal variance of $\hat{a}$. Note that
although this result looks quite intuitive it is nontrivial since the
$\hat{a}$-squeezed coherent state is not a product state of a
$\hat{b}$-squeezed coherent state and a $\hat{c}$-squeezed coherent state,
such that the haloscope subvolume is entangled with the rest of the halo.

This means that at $t_{\rm sqz}$, defined by $V_{-}(t_{\rm sqz})=e^{-2}$, when
the galactic $\hat{a}$ mode gets significantly squeezed, the $\hat{b}$ mode
squeezing is

$V^{h}_{-}(t_{\rm sqz})=1+\frac{N_{h}}{N}(e^{-2}-1)\,,$ (75)

or assuming $N_{h}/N\ll 1$

$r^{h}(t_{\rm sqz})\simeq 0.43\frac{N_{h}}{N}\ll r(t_{\rm sqz})=1\,.$ (76)

Thus, the squeezing of a subvolume ($\hat{b}$ mode) is significantly reduced
compared to squeezing of the full volume ($\hat{a}$ mode), and bounded by
$N_{h}/N$, the ratio of the number of axions in the subvolume and the total
number of axions. Similarly, the maximum squeezing of the $\hat{b}$ mode is
drastically reduced. Since $r_{\rm max}<\infty$, one finds

$r^{h}(t_{\rm max})\lesssim 0.5\frac{N_{h}}{N}\ll r_{\rm max}\,.$ (77)

For the galactic halo we assume a stationary Navarro–Frenk–White (NFW) density
profile Navarro  _et al._ (1996)

$\rho(r)=\frac{\rho_{0}}{\frac{r}{R_{s}}\biggl{(}1+\frac{r}{R_{s}}\biggr{)}^{2}}$
(78)

to calculate $\chi$. The parameters $\rho_{0}$ and $R_{s}$ vary from halo to
halo. Solving the Poisson equation, the gravitational potential is

$\Phi(r)=-\frac{4\pi
G\rho_{0}R_{s}^{2}}{\frac{r}{R_{s}}}\ln\Bigl{(}1+\frac{r}{R_{s}}\Bigr{)}$ (79)

so that

$N\chi=\frac{1}{2}\frac{1}{N\hbar}\int_{V}d^{3}\\!x\,\Phi(\bm{x})\rho(\bm{x})\,.$
(80)

A sufficiently accurate approximation to the total gravitational energy is
given by

$N\chi\simeq\frac{1}{2}\frac{M_{\rm
MW}}{N\hbar}\bar{\Phi}\simeq-\frac{1}{2}\frac{m}{\hbar}10^{-6}$ (81)

where $\bar{\Phi}\simeq 10^{-6}$ is the average Galactic potential and $M_{\rm
MW}\simeq 10^{12}M_{\odot}$ the Galactic mass. Substitution of $N\chi$ into
(5) gives

$t_{\mathrm{sqz}}\simeq
30\biggl{(}\frac{10^{-5}\mathrm{eV}}{m}\biggr{)}\,\mathrm{\mu}\mathrm{s}\,.$
(82)

To evaluate the timescale of maximal squeezing (and its magnitude expressed by
$r_{\mathrm{max}}$) we need

$N=\frac{M_{\rm
MW}}{m}=10^{83}\biggl{(}\frac{10^{-5}\mathrm{eV}}{m}\biggr{)}\,.$ (83)

Furthermore, we need to find the approximate number of axions in the axion
haloscope volume Asztalos  _et al._ (2010)

$N_{h}=\frac{M}{m}=\frac{\rho(r=8\mathrm{kpc})V_{h}}{m}\simeq 6\times
10^{18}\biggl{(}\frac{10^{-5}\mathrm{eV}}{m}\biggr{)}\,.$ (84)

Inserting this into (6) and (7) we get

$r_{\mathrm{max}}\simeq\ln(N^{1/6})\simeq
31.8+\frac{1}{6}\ln\biggl{(}\frac{10^{-5}\mathrm{eV}}{m}\biggr{)}$ (85)

and

$t_{\mathrm{max}}\simeq 0.5t_{\mathrm{sqz}}N^{1/6}\simeq
33000\biggl{(}\frac{10^{-5}\mathrm{eV}}{m}\biggr{)}^{\frac{7}{6}}\mathrm{yr}\,.$
(86)

For comparison, we find that the Ehrenfest timescale is approximately
$t_{\mathrm{Ehr}}=\sqrt{N}t_{\mathrm{sqz}}\simeq
10^{32}(\frac{10^{-5}\mathrm{eV}}{m})^{\frac{3}{2}}\mathrm{yr}$.

Finally, let us evaluate the numerical value of the squeezing of the haloscope
mode at $t_{\rm sqz}$ and $t_{\rm max}$. We find from (76) and (77) that

$r^{h}(t_{\rm sqz})\simeq r^{h}(t_{\rm max})\lesssim
0.5\frac{N_{h}}{N}=10^{-65}\,.$ (87)