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

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\"odinger-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 1 eV) selfgravitating quantum scalar field, such as the QCD axion or any other axion-like particle (ALP). Extensive experimental searches focus on its distinct wavelike signatures [1][2][3], and its theoretically expected coupling to matter [4][5][6]. This type of ALP dark matter, also including ultralight or "fuzzy dark matter," is usually described by a classical scalar field rather than particles [7][8][9][10][11][12][13][14][15][16]. The classical field description is motivated by the de Broglie wave length λ deB = /(mv rms ) being significantly larger than the inter-particle distance d, and thus the mean phase space occupation number (λ deB /d) 3 is large. If a misalignment mechanism produces the ALP then initially it would be in a coherent state [7,8], 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 [17]. 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 [18]. 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 [19][20][21][22], with many applications in quantum-enhanced metrology [23] and quantum information processing [24,25]. The electromagnetic case is also instructive to understand the difference between the classical and quantum states of the * michael.kopp@su.se † vasileios.fragkos@fysik.su.se ‡ igor.pikovski@fysik.su.se ALP. The starting point of quantum optics and our understanding of the difference between classical and quantum coherence was the Hanburry-Brown-Twiss experiment [26]. 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 [27] and Sudrashan [28] 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 selfgravitating 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 µ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) [29] we find that squeezing becomes negligible even if the galactic halo undergoes strong squeezing. A QCD axion, for example, reaches 9dB 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µ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 [14,30], 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 nonrelativistic scalar quantum fieldψ(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ψ(x, t) is, see for instance [31], where A(t) is the scale factor, which appears due to our use of comoving coordinates 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 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 =Ȧ(t)/A(t) is the Hubble expansion rate, and ρ 0 = mn 0 the presentday 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. 1 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. 1 Inclusion 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.

A. 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 [32][33][34][35]. 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.

B. System of Kerr and Gross-Pitaevski equations
Projecting the Heisenberg equation (1) onto the mode function, that is 1 (1)), we obtain the Kerr oscillator equation forâ(t): The GPE field is normalised to V d 3 x|ψ(x, t)| 2 = N , so that the spatial average is |ψ| 2 = N/V ≡ n 0 . Here V is the comoving volume, assumed to be very large, roughly the size of the Hubble volume H −3 0 .
with coefficients (see Appendix E for a more rigorous derivation), The last equality (3f) follows from the Gross-Pitaevskii equation (GPE) for ψ(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ψ(x, t) → ψ(x, t), and thusâ(t) → √ N . While this indeed solves (3), it is not strictly a valid solution since it does not satisfy [â,â † ] = 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â is equivalent to a Kerr-Hamiltonian which is well-known in quantum optics [19,36] and in the context of BECs [32]. An initially coherent state |α = √ N evolves unitarily according to the HamiltonianĤ = ωâ †â + χâ †â †ââ , with timeindependent parameters ω and χ. The Heisenberg equation of motion i ȧ = [â,Ĥ] is equivalent to (3) with time-independent coefficients and with solutionâ(t) = e −it(ω+2â † (ti)â(ti)χ)â (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.

C. Quadrature squeezing
Defining the quadrature operatorX θ (t) =â(t)e −iθ + a † (t)e iθ , where θ is the quadrature angle, squeezing is present at θ = θ − (t) if the minimal variance V − (t) ≡ Var(X θ−(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 [36]. In the limit of large N , the minimal variance is given explicitly in Appendix B, eq. (B3). We define the squeezing timescale as the time on which V − (t) crosses below e −2 0.135, i.e. the vacuum noise is reduced by a factor 0.135 (or about 9 dB squeezing). Or in terms of the squeezing parameter the squeezing timescale is r(t sqz ) = 1. This yields For large N , V − (t) approaches 0 very closely, with its mininum corresponding to a maximum of r(t) given with r(t max ) = r max and The squeezing angle defining the orientation of the squeezed quadrature can be approximated by which for t > t sqz approximately co-rotates with the classical solution a c (t) = αe −it(ω+2χN ) of the Kerr model which is a good approximation to the mean â(t) . The conjugate quadrature with θ + = θ − + π/2 is exactly antisqueezed. The N -dependent phase shift provided by the arc tangent approximately cancels for t > t sqz with the remaining offset at t max vanishing as sgn(χt max )N −1/6 , thus producing an approximately amplitude-squeezed state. 3

D. 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ω(a) by treating ψ(x, t) in perturbation theory (see Appendix D), or by solving for it numerically using the full GPE. In both cases one approximately finds where mK is the average ALP kinetic energy, and v 2 rms is the density weighted average of the ALP velocity field squared, see [31], or in the context of cold dark matter (CDM) eq. 106 in [39]. Numerical simulations show v 2 rms (A = 1) 10 −6 in the present-day universe. This is the result of nonlinear structure formation through gravitational instability of ψ from primordial density perturbations. On scales larger than λ deB = /(mv rms ), |ψ(x, t)| 2 and the density of of cold dark matter are virtually the same [1,10,31,40]. Therefore estimating K(t) can be done with any conventional CDM N-body simulation or CDM perturbation theory. The squeezing timescales are determined by χ which is related to the average potential energy W = m N χ [31]. Using the Layzer-Irvine equation eq. 47 of [31], we get with K(A) ∝ A from (9), (10) We have therefore tied all parameters appearing in the Kerr oscillator to a simple expression forω(t), eq. (9), entirely fixed by known cosmological parameters and the ALP mass m.
The number N of axions in the volume Note that this number deviates from ALP mode occupation numbers used elsewhere in the literature, e.g. 10 61 in [42] and 10 26 in [13] 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.

E. Onset of squeezing
We find the squeezing timescale of the ALP in the present-day universe from (5), (9) and (10)   We have thus shown that squeezing is rapidly generated for an ALP, much faster than the age of the universe H −1 0 10 18 s for all m > 10 −22 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 |ω/ω| |χ/χ| H. Note that the (quantum) squeezing time scale λ deB /v rms of the operatorâ(t) essentially coincides with the (classical) coherence time scale of the GP field ψ(x, t), see for instance [43].
For a QCD axion with m 10 −5 eV, this timescale is about 66µs. One can compare this to other quantum mechanical timescales appearing in the Kerr model, such as the Ehrenfest timescale t Ehr = N 1/2 t sqz , or Schrödinger "cat" creation time t cat = N t sqz [36]. Both of these timescales far exceed the age of the universe for the ALP mass range.

F. Maximal squeezing
The maximum squeezing that can be produced in ALPs is obtained by inserting (11) into (6): 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 [44][45][46][47]. The time t max , eq. (7), at which maximal squeezing is reached is Therefore even the extreme squeezing r max would be easily reached in the present universe, justifying our use of constant Kerr model parameters, for m ≥ 10 −12 eV.
Since ω + 2N χ = 0 the squeezing angle evolves slowly between θ − − π 4 and θ − − π 4 N −1/6 in the interval 0 < t t 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 max .
As an example we show in Fig. 1 the evolution of an initial coherent state with α = 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 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 [10][11][12][13][14][15][16] 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 [1]. 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 [14,[48][49][50][51], 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 [30], which characterizes the onset of non-linearities in the classical description. It can also be related to the thermalisation timescale found in [42,48], 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 [48]. This thermalisation timescale can be related to our squeezing timescale as Ht th 1 Htsqz . 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 [52], as well as spin squeezing of internal states [53]. The ground state of a trapped BEC has been shown to be a squeezed displaced state [54], as found here, see the Appendix G. However, in lab experiments the trapping typically dominates selfinteractions 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 [55][56][57] 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 [58,59], 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. [60], 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. [55][56][57], 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  [61,62]. At t = t squeezed states of light can have 0 < G (2) I / â †â < 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 selfinteractions. For instance number-squeezing was experimentally observed in a multi-well-trapped BEC [63], 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 [64]. Multiple localisation sites, analogs of ALP DM haloes, are known to lead to the emergence of several macroscopically occupied GP-like fields [65]. A systematic approach that reveals the necessity for and determines the dynamics of additional mode functions has been established in [34,66], and applied to what could be considered analog situations of ALP cosmology [67,68]. 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 [47,69]. 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 ex-otic coupling mechanisms such as time-dilation induced effects [70]. 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 [71], in which two-mode squeezing of the axion is discussed. The mechanism and the type of squeezing studied in [71] 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 [71] 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 [44][45][46][47].

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.

Appendix B: Details of the Kerr model
It is convenient to writeĤ = (ω −χ)n+ χn 2 because the linear and nonlinear parts ofĤ commute in this form so that U = e −iĤt/ = e −i(ω−χ)tn e −iχtn 2 . For a coherent state |α this implies â †â = |α| 2 = N . We assume again w.l.o.g that α = √ N . To evaluate V − and θ − in the Kerr model, the only additional quantities we need are â and â 2 . We calculate these expectation values using the result eÂBe −Â = Inserting this into (A3) and (A4) gives the squeezing angle θ − (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 [38] a very accurate algebraic approximation for V − (t) has been obtained by substituting N → τ /(χt) and performing a Taylor expansion in t (with τ fixed) which gives Note that V − (t) is independent of the sign of χ. To obtain the squeezing timescale we hold τ fixed while taking the limit N → ∞ in (B3). Then we solve − 1 2 ln V − (t sqz ) = 1 for t sqz to obtain (12). To obtain the time of maximal squeezing we replace |χ|t →τ N −5/6 in (B3), and while holdingτ fixed, we expand in 1/N . The leading term is Taking aτ derivative, we find the minimum atτ max = 2 −5/3 , which gives (7). Inserting this back into V − (τ max ) corresponding to r max = − 1 2 ln V − (τ max ) and results in (7). Fig. 2 shows t sqz , t max and r 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 θ − . We proceed in a similar fashion as we did for V − . We insert (B1) and (B2) into Var(â)/(Var(â)) * and we replace N → T /(χt). Then we take the log and Taylor expand the result in χt with the result This can be simplified into the expression (8). The signum function in (8) arises due a branch cut in the log.
The squeezed coherent state approximation breaks down due to non-Gaussianities around the time t max . In order to verify that a squeezed coherent state is a good approximation to the Kerr-evolved coherent state until t max we evaluate Var(n) = 1 4 cosh(4ρ)−1 +N (cosh(2ρ)−cos(2φ) sinh(2ρ)) (C10) for the squeezed coherent state | √ N , r(t)e 2iφ with φ = θ − (t) and φ = 0 in Fig. 3. We see that inclusion of the squeezing orientation is crucial to match the value Var(n) = N of the Kerr-evolved coherent state. This also explains the "coincidence" of θ − (t max ) π 4 N −1/6 and V 1/2 + 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 . 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 µ andω 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(x, t) and velocity field u(x, t) satisfying the Euler-Poisson equation. Since approximately 2 m 2 |∇ψ(x, t)| 2 = n(x, t)|u(x, t)| 2 , see e.g. [10], [31], the mean kinetic energy per mass is where v = u/A is the peculiar velocity field, and u is the canonical velocity field. In linear perturbation theory the comoving density field can be written as n(x, A) = n 0 (1+ δ(x, A)) with the density perturbation δ(x, A) related to the velocity field via ∇ · v = −AHδ(x, A). During the assumed matter dominated expansionȦ/A = H = H 0 A −3/2 we thus get to leading order in perturbation theory In Fourier space, taking the limit V → ∞, this becomes where we used the ergodic theorem to replace |δ k | 2 → P (k = |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 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 ∝ (AHf ) 2 . Here f Ω 6/11 m is the growth rate of linear matter density perturbations, Ω m = Ω m,0 A −3 H 2 0 /H 2 the fraction of matter and H = H 0 (Ω m,0 A −3 + 1 − Ω m,0 ) 1/2 the expansion rate. w.r.t. φ * (x, t) and c * n (t) gives Neglecting the time dependence of ω and χ we find which is the solution of the Schrödinger equation for a Kerr oscillator. For an initial coherent state |Ψ(t = 0) = |α with α = √ N we have c n (t = 0) = e −N/2 N n/2 √ n! so that n(t) = N and n 2 (t) = N 2 +N . The φ evolution equation, (E5) and (E6), then simplifies to which is the GPE used in the main text, after redefining ψ ≡ φ √ N and restoring the scale factor A(t). Equations (E15), (E16) then establish the validity of the constraints d 3 x |φ(x, t)| 2 = 1 and n |c n (t)| 2 = 1, which we could have enforced via Lagrange multipliers [34]. Noticing that n, t|Ĥ − i ∂ t |n, t = n|Ĥ|n , with |n = |n, t = 0 andĤ = (ω − χ)n + χn 2 , we verify (3) forâ(t) = e iĤt/ b (t = 0)e −iĤt/ in the Heisenberg picture.
This equals to Ψ K (t = 0)|e iĤt/ b (t = 0)e −iĤt/ |Ψ K (t = 0) = α|â(t)|α , 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 G: Ground state in a single-GP-mode approximation If thermalisation as described in [48] indeed keeps the axion at its momentary ground state |Ψ 0 (t) , then it might be natural to expect that this ground state is related to our squeezed state. [48] 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 [32,54]. To leading order in a perturbative solutionâ = α + δâ, this ground state is found to be a squeezed coherent state matching our approximate solution (F1) at t 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 max . Interestingly, the relation to the BEC ground state might persist even when non-Gaussianities are included: it was found in [32,54] that the ground state Wigner function of the Kerr Hamiltonian is bent into a "banana" producing a so-called number-squeezed state.
relation betweenâ andb operators we find using a calculation similar to that presented in Appendix B S † (ζ)bŜ(ζ) = √ N h √ N â cosh(r) −â † e 2iθ− sinh(r) , (I7) with α = √ N and ζ = re 2iθ− . Using these results (I5) can be simplified most easily in the "Kerr picture", see Appendix F, to give where V − (t) = e −2r(t) is the minimal variance ofâ. Note that although this result looks quite intuitive it is nontrivial since theâ-squeezed coherent state is not a product state of ab-squeezed coherent state and aĉ-squeezed coherent state, such that the haloscope subvolume is entangled with the rest of the halo. This means that at t sqz , defined by V − (t sqz ) = e −2 , when the galacticâ mode gets significantly squeezed, thê b mode squeezing is or assuming N h /N 1 r h (t sqz ) 0.43 N h N r(t sqz ) = 1 .
Thus, the squeezing of a subvolume (b mode) is significantly reduced compared to squeezing of the full volume (â 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 theb mode is drastically reduced. Since r max < ∞, one finds For the galactic halo we assume a stationary Navarro-Frenk-White (NFW) density profile [75] ρ(r) = ρ 0 A sufficiently accurate approximation to the total gravitational energy is given by Furthermore, we need to find the approximate number of axions in the axion haloscope volume [76] N h = M m = ρ(r = 8kpc)V h m 6 × 10 18 10 −5 eV m .
For comparison, we find that the Ehrenfest timescale is approximately t Ehr = √ N t sqz 10 32 ( 10 −5 eV m ) 3 2 yr. Finally, let us evaluate the numerical value of the squeezing of the haloscope mode at t sqz and t max . We find from (I10) and (I11) that