Axion Structure Formation I

# Axion Structure Formation I: The Co-motion Picture

Erik W. Lentz, Thomas R. Quinn, Leslie J Rosenberg
Institut für Astrophysik, Georg-August Universitat Göttingen, Göttingen, Deutschland 37707;
Department of Astronomy, University of Washington, Seattle, WA, USA 98195-1580;
Department of Physics, University of Washington, Seattle, WA, USA 98195-1580;
E-mail: erik.lentz@uni-goettingen.de
Accepted XXX. Received YYY; in original form ZZZ
###### Abstract

Axions as dark matter is an increasingly important subject in astrophysics and cosmology. Experimental and observational searches are mounting across the mass spectrum of axion-like particles, many of which require detailed knowledge of axion structure over a wide range of scales. Current understanding of axion structures is far from complete, however, due largely to controversy in modeling the candidate’s highly-degenerate state. The series Axion Structure Formation seeks to develop a consistent model of QCD axion dark matter dynamics that follows their highly-degenerate nature to the present using novel modeling techniques and sophisticated simulations. This inaugural paper presents the problem of describing many non-relativistic axions with minimal degrees of freedom and constructs a theory of axion infall for the limit of complete condensation. The derived model is shown to contain axion-specific dynamics not unlike the exchange-correlation influences experienced by identical fermions. Perturbative calculations are performed to explore the potential for imprints in early universe structures.

###### keywords:
cosmology: dark matter – cosmology: early Universe – galaxies: formation – galaxies: haloes – galaxies: structure
pubyear: 2018pagerange: Axion Structure Formation I: The Co-motion PictureLABEL:lastpage

July 20, 2019

## 1 Introduction

The QCD (quantum chromodynamic) axion is a well-motivated candidate to solve problems in both particle physics and cosmology. The theory of axions originated in 1977 as a result of a scalar-induced axial symmetry over the QCD sector, used to solve the strong (QCD) CP problem of particle physics (Peccei and Quinn, 1977). The axion particle was proposed the following year from a spontaneous breaking of that axial symmetry, solving the strong CP problem dynamically for energies about and below MeV (Weinberg, 1978; Wilczek, 1978). Experiments and astrophysical observations have replaced this Peccei-Quinn-Weinberg-Wilczek axion (Mimasu and Sanz, 2015; Patrignani and Particle Data Group, 2016) with more complex and unified models, pushing the breaking scale ever higher and the axion mass ever lower. These new theories also have increasingly addressed cosmology (Marsh, 2016), where the dark matter (DM) problem has been gaining in prominence.

As a DM candidate, the QCD axion is highly attractive due to its well-bounded parameter space of mass and couplings to the standard model (SM) of particle physics, Fig. 1. Starting at milli-eV masses, there is a bound above which the axion would have been seen in various astrophysical processes (Raffelt, 2008; Isern et al., 2010; Córsico et al., 2012; Viaux et al., 2013; Patrignani and Particle Data Group, 2016; Marsh, 2016) such as anomalous energy transport in SN1987a. Approaching micro-eV masses, there is a bound below which the maximal (misalignment) creation mechanism would produce more axions than there is DM (Abbott and Sikivie, 1983; Dine and Fischler, 1983; Preskill et al., 1983; Patrignani and Particle Data Group, 2016). This lower bound is somewhat soft as axion creation mechanisms can be suppressed in the details of some axion theories (Marsh, 2016). The two diagonal lines of Fig. 1 represent benchmark axion models. KSVZ (Kim-Shifman-Vainshtein-Zakharov) represents a theory where the axion couples to hadrons only (Kim, 1979; Shifman et al., 1980), and DFSZ (Dine-Fischler-Srednicki-Zhitnitsky) couples to both hadrons and leptons (Dine et al., 1981; Zhitnitsky, 1980) as consistent with grand unified theories. The search window is then given by the region between KSVZ and DFSZ and the lower and upper mass bounds.

All direct DM searches have in common the need to understand the candidate abundance in the vicinity of their experimental apparatus. Astronomical observations are unable to provide accurate estimates of abundances due to difficulties in extracting local DM contributions from the ephemeral memory of gaseous and stellar motions, though some progress has been made in recent years (Valluri et al., 2016; Binney, 2017). Instead, simulations that capture the relevant physics of formation are performed to create Milky Way (MW) analogues, and are used to estimate DM populations. Methods for performing these simulations are already quite sophisticated, capable of detailing the formation of cosmological structure down to the inner workings of galaxies. Often starting from cosmic microwave background (CMB)-motivated initial conditions, simulations of structure formation must accurately describe the violent process of collapse, collision, and equilibration of matter into stars, galaxies, halos, and other forms.

There are several critical points in the early universe that significantly impact the modern relic axion state. The first centers on the relative placement of the spontaneous symmetry-breaking scale of the axion’s parent scalar field, , and the inflation scale where is the associated cosmological radiation temperature. In brief, if the PQ symmetry is unbroken during inflation, meaning that , then the axion’s parent scalar field has zero vacuum expectation value (VEV), leaving the angular axion field randomly distributed over each Hubble patch when the is reached after inflation. Such a distribution results in the large average misalignment angle . If the PQ symmetry is broken during inflation, varies on scales above the Hubble volume post-inflation. This effectively makes uniform in the visible universe and close to a free parameter of the model.

Issues occur with each inflation possibility. The randomness of in the pre-PQ breaking inflation case gives rise to topological defects and seeds other structures that may prove difficult for local searches such as haloscope experiments, as much of the mass may exist in compact objects such as mini-clusters (Kolb and Tkachev, 1993; Tinyakov et al., 2016; Fairbairn et al., 2017). Post-PQ breaking inflation, which produces relatively lighter axions, introduces a degree of anthropic reasoning to the state of the visible universe. A narrow range of misalignment angles are allowed to reproduce the observed DM density, severely constraining the parameter. Even with inflation-enforced flatness, striking the necessary balance of axion DM with other components makes even a semi-random axion phase precarious.

The second critical point occurs when the post-symmetry-breaking post-inflation misalignment angles leave the freeze-out state. This occurs for the QCD axion during the QCD phase transition. The condensation of the quark-gluon plasma creates an effective potential for the axions such that the potential minimum is located at a point of vanishing CP violation. There are several popular models for producing cosmological axions, including the topological defect decays from phase windings, decays from parent particles, thermal populations, and vacuum realignment. We restrict ourselves to the defect-free realignment mechanism at this time, though the other possibilities may be explored in future work.

The misalignment process is very clean, creating nearly homogeneous axions within a Hubble patch, and restricting the scale of spatial oscillation to be of order the horizon size or larger in the absence of topological defects. The initial axion state is incredibly degenerate, with the canonical critical temperature estimated to be of for micro-eV axions, far above the QCD axion mass (Banik and Sikivie, 2015). It is reasonable to assume that the axions inhabit some form of Bose-Einstein condensate (BEC) from the expected degeneracy alone.

Axion DM dynamics go through several stages after creation. Radiation gravity and the potential of the axion field dominate axion interactions early after creation, respectively equilibrating with the cosmological radiation perturbations or inducing violent collapse in areas of large over-density (Redondo, 2013; Gorghetto et al., 2018; Vaquero et al., 2018). Axion self-interaction reduces quickly with further cosmological expansion, leaving dynamics dominated by axion-conserving interactions from environmental and self gravity. QCD axion DM interactions are dominated by gravity by the time of the matter-radiation transition; this letter concentrates on structure formation after that transition.

The state of axion structure formation after self-interactions consists primarily of two viewpoints. The first finds sufficiency in a mean-field description of a condensed axion fluid while the second holds that the full quantum description is necessary to describe the axion field. Relic axions created under the misalignment mechanism easily satisfy the canonical requirements of a thermal BEC. It therefore seems reasonable to assume that every axion is in essentially the same state. Proponents of the mean field description motivate simplification along this single-state constraint over every axion. Implementing such a mean field theory (MFT) on relic axions is straightforward and attractive. In the limit of ultra-light axions (ULAs), eV, the de Broglie wavelength scale enters the galactic scale (Arvanitaki et al., 2010; Schive et al., 2014; Veltmaat et al., 2018). QCD axions are at least sixteen orders of magnitude heavier and will not see outright wave diffusion effects on galactic scales. This reduces the dynamics to a classical pressure-less fluid, the same used in the cosmological standard model CDM. The overlap of axion and standard cold DM (CDM) models would relegate the detection of axions to direct or otherwise non-gravitational methods.

The second portion of the axion structure formation community proposes that quantum effects extend far beyond the de Broglie scale. This invites the possibility that QCD axion DM may produce observable signatures unique to axions on galactic scales. Models built on requirements of coherent energy and angular momentum conservation at the quantum level produce halos where axions collapse as a series of cold phase-space sheets, drawn preferentially towards the galactic plane by their rotational properties. The resulting post-collapse halo is found to contain of a series of caustic rings, the effect of catastrophes (Sikivie, 1999a, b; Duffy and Sikivie, 2008; Chakrabarty and Sikivie, 2017). It is also claimed that our own solar system is close to one of these caustics of a fitted MW axion halo. Cold flows feeding into the outermost rings would contribute to the local axion energy spectra, significantly influencing the results of terrestrial search experiments. Proponents of caustic rings appreciate the compactness of the MFT description, but believe that it does not adequately maintain the quantum nature of the degenerate Bose system. The fault with the MFT simplification is reported to lie in the mis-estimation of the relevant timescales of gravitational equilibration, leading to an imbalance in state occupation (Banik et al., 2015; Erken et al., 2012). The more complete axion models used in the gravitational thermalization research are implemented at the level of individual axion quanta and unfortunately cannot be extended to cosmological scales.

Non-scalability notwithstanding, a more complete quantum description of axions than MFT is required for several reasons. First, current axion MFTs do not track inter-axion correlations. Standard condensate MFTs use a single-particle Fock decomposition of the axion solutions combined with an assumption of separability between axion states to produce a compact single Slater permanent representation of the system. However, gravity is an infinite range, un-screenable, and highly-correlated interaction that makes such a separation in the gravitational solution space unfeasible. It is likely that the constraint of axion infall models to such confined states introduce unphysical forces to the axions’ dynamics.

Further, the expected energy level of a cosmological system of axions is nowhere close to its ground state, nor is it expected to reach thermal equilibrium. Classically, self-gravitating systems collapse to a virialized quasi-equilibrium, which is known to be distinct from a thermalized fluid by phenomena such as the gravi-thermal catastrophe (Binney and Tremaine, 2008). Gravitationally-virialized fluids are also expected to maintain knowledge of their infall history, which may be valuable to astronomers and dark matter researchers (Valluri et al., 2016). Structures formed with non-local Newtonian gravity may conceptually be quite different than the canonical thermodynamic approach, requiring rigorous definitions of condensation and condensates. The ignoring of correlations in current axion MFTs also extends to the interactions, smoothing out potential two-body forces. There is therefore a need for a beyond-MFT model of axion dark matter, capable of tracking axions in and out of equilibrium.

It is the objective of this series to provide a model of axion structure formation sufficient to resolve axion-specific physics on astronomical scales and determine the extent of unique structures. This first letter derives a new model of QCD axion DM structure formation in the condensed limit and calculates its potential imprints on the early visible universe. The remainder of this paper is structured as follows: Section 2 describes the context for axion structure formation after the QCD phase transition and into modern times; Section 3 introduces the new continuum model of structure formation for a system of perfectly condensed self-gravitating bosons; Section 4 performs calculations of condensed axion structure formation in the small-perturbation limit, exploring the extent of unique axion structures in the early universe; Section 5 discusses the implications of the new model and the perturbative results; and Section 6 summarizes the results and outlines the scope of upcoming letters.

## 2 Axion Physics after the QCD Phase Transition

Post-QCD phase transition axions are well-described by a single quantum scalar field operator () defined over a dynamic 3+1 pseudo-Riemannian manifold () governed by classical Einstein’s equations. Such a full description is extremely accurate but excessive for common cosmological structures, so several simplifications are made in the interests of computability. On the largest scales, the geometry is taken to be nearly homogeneous, isotropic, while the smaller scales of interest are taken to satisfy weak-field conditions (Misner et al., 1973). The space is then found to evolve under a Friedmann-Lemaitre-Robertson-Walker (FLRW) cosmology with small geometric perturbations. Cosmological measurements favor the flat variety of FLRW cosmologies (Planck Collaboration et al., 2018), requiring a critical energy density among all gravitating species. Einstein’s equations then decouple in orders of metric perturbations, with the background flat FLRW cosmology evolving under the horizon-scale-averaged stress-energy

 H2−Λc23=8πG3¯ρ, (1) 2¨aa+(˙aa)2−Λc2=−8πGc2¯p, (2)

where is the total expected mean co-moving energy density, is the total expected mean pressure, is the FLRW scale factor, is the Hubble flow, and the dots represent derivatives in physical time as measured by a co-moving observer. The next order terms obey a parabolic potential-like equation.

 ΔΦ=4πGac4ρ(1) (3)

where is the flat spatial Laplacian, is the deviation in the frame’s energy density from the mean, and the Newtonian potential is related to the co-moving geometry via where is the diagonal time component of the co-moving metric perturbation from FLRW. Relativistic gravity as a classical theory can be interpreted as averaging over quantum length scales. Free fundamental particles are associated with two quantum length scales, Compton and de Broglie, though only Compton is covariantly meaningful; the de Broglie lengths parameterize frame-dependent corrections from wave kinematics. As the fundamental length scale of a particle, the Compton length determines the limit by which a light-like vector field may distinguish space-like structure. The Compton scale is on the order of centimeters for QCD axion DM. This is the smoothing length that shall be used to evaluate the axion gravitational potential.

The effective axion equations of motion are also derived from a fully covariant description. Describing the axion system as a second-quantized operator, the axion field’s governing equations are characterized by a Lagrangian operator of the form

 (4)

where is a covariant derivative over and is the axion interaction potential. The axion field operator is of Bose type, meaning that it commutes under particle exchange. The field obeys the Klein-Gordon equation over effective fields in the semi-classical limit

 ∇μ∇μ^ϕ+(maℏ)2^ϕ−λ3!|^ϕ|2^ϕ=0 (5)

where the fourth order expansion of the axion potential is used instead of the full cosine axion potential

 V(^ϕ)=m2aℏ2|^ϕ|2−λ4!|^ϕ|4+O(|^ϕ|6) (6)

and is the coupling strength of the four-field coupling strength. Evaluating the covariant derivatives in the linear approximation, the derivative term of the Klein-Gordon operator becomes

 gμν∇μ∇ν^ϕ =a−2□^ϕ+3a−2H∂t^ϕ −a−2(h−1)μν∂μ∂ν^ϕ−3a−2H(h−1)0ν∂ν^ϕ −a−2h∂0^ϕ+12a−2∂ν(h)∂ν^ϕ+O(h2) (7)

where is the flat d’Alembert operator, and all index contraction, raising and lowering, are performed using the Minkowski metric () in agreement with the previously-mentioned weak field formalism (Misner et al., 1973).

The effective Newtonian limit eliminates the lowest order terms and is consistent with non-relativistic fields in the co-moving frame, allowing for the non-relativistic field expansion

 ^ϕ=ℏ√2m(^ψe−imat/ℏ+^ψ∗eimat/ℏ) (8)

Substituting this form of the field operator, and taking the Newtonian limit for gravity, , the real time Klein-Gordon operator equation reduces to

 ^ψ∗(∂t−32H)^ψ+^ψ∗ℏ2→∇22a2ma^ψ+maΦ|^ψ|2 −λℏ43!m4a|^ψ|4+{h.c.}+(quickly oscillating terms) =O(h2,m2/ℏ2,mh/ℏ) (9)

The quickly oscillating terms are ignored as their actions evaluate to over first-order timescales. The terms of the effective axion action can be easily related to the terms of the field Hamiltonian operator as the action measure includes a factor of :

 ^Hf =^K+^U+^V (10) ^K =−^ψ∗ℏ2→∇22a2ma^ψ (11) ^U =λ4!|^ψ|4 (12) ^V =^ψ∗maΦ^ψ (13)

which confirms that is indeed the canonical Newtonian gravitational potential.

The many-body axion Schrödinger operator equation can now be written in the comoving parameterization

 iℏ∂t^ψ=−ℏ2→∇22a2ma^ψ+λ4!|^ψ|2^ψ+maΦ^ψ (14)

This is the operator equation governing the quantum axion field in the presence of classical gravity, which serves as a robust starting point for much of axion structure-formation models after the QCD phase transition.

Interactions during the remainder of the radiation era are dominated by gravitational perturbations from relativistic species and the axion point interaction. Structure formation of axions are linear so long as axion perturbations remain small during appreciable self-interaction. After a time the axions diffuse beyond appreciable self-interaction, the mean free path exceeds the Hubble radius, and their number becomes a conserved quantity.

The expansion rate and diffuse gravitational potentials keep the axions in the perturbative regime during the remainder of radiation domination. The field operator picture is still sensible during this time of point-wise and separable potentials. Even so, there is some sub-dominant but increasingly important physics at work. Inter-axion gravitation connects every axion to every other axion via the long-range un-screenable Coulomb potential

 Vgrav(→xi,→xj)≅−4πGm2a|→xi−→xj| (15)

for well-separated axions in this limit of gravity. Such a highly-correlated interaction is thought to maintain the axion condensate through “gravitational thermalization” (Erken et al., 2012), though this may be somewhat irrelevant as normal dispersion processes are not anticipated to push the DM out of the canonical condensate phase.

As matter energy density establishes equality and then dominance over radiation, axion self-gravity becomes the primary source of structure-formation (SF) dynamics. Maintaining the field operator form of the system becomes more difficult during this process as the Fock representation of the many-axion system becomes more tenuous. Degradation in the separability and perturb-ability of the Hamiltonian at best requires the use of a longer chain of Fock states. Regardless, one can expect computations over a separable set of axions to become untenable as we enter the DM and dark energy (DE) dominated eras. MFTs are then seen to fail in describing axion SF in the late universe as they are constrained from tracing the correlated states. The next section describes a new model of axion SF that more fully explores the state space.

## 3 A New Model of DM Axions

This section follows similar path to Lentz et al. (2018). Embellishment on the finer points of Schrödinger and Wigner methods in the context of axion self-gravitation are given in Appxs. A, B, C.

### 3.1 Many-Body Quantum Axions as a Many-Body Schrödinger Equation

While in principle correct, the second quantization formalism used by many axion MFTs is not well-suited for axion structure formation. It encourages a separable decomposition among axions not supported by their dynamics, and supplies a superfluous freedom to create/annihilate said separated particles. The technique presented here instead restricts its scope by falling back to first quantization, made possible by conserved axion number. Rephrasing the field theory into a quantum mechanical approach, the Schrödinger picture of Hilbert space uncovers a many-body wave-tracing evolution equation. Evaluating the dynamical equation over an arbitrary element of the many-body solution space gives the form

 iℏ∂tΨ(→x1,...,→xm;t)=^HΨ(→x1,...,→xm;t) (16)

where is the many-body axion wave-function and is the many-body Hamiltonian operator on wave-function space. Note that the symbol for the first-quantized many-body Hamiltonian operator is similar but distinct from the Hubble flow . The wave-function also inherits an exact exchange symmetry due to the commuting nature of the axion field operator.

 Ψ(→x1,...,→xi,...,→xj,...,→xN;t) =+Ψ(→x1,...,→xj,...,→xi,...,→xN;t) ∀i,j (17)

The many-body Hamiltonian may be broken down into recognizable kinetic and potential energy parts

 ^H=^T+^V=M∑i−ℏ2∇2i2a2ma+∑iΦi (18)

where is the gravitational potential on the i-th axion.

The gravitational interaction requires further clarification. Strictly speaking the gravitational potential may come from multiple massive and radiative species. Considering only massive species like axion DM, stars, atomic gas, etc., the potential may be broken up linearly, using the Coulomb gauge choice where the potential vanishes when infinitely far from all sources:

 ρ =ρa+ρ′ and (19) ∇2Φ =4πGρ=4πG(ρa+ρ′) (20)

implies

 Φ =Φa+Φ′, (21)

where

 ∇2Φa =4πGρa (22) ∇2Φ′ =4πGρ′ (23)

and all potential terms not sourced by axionic densities are consolidated into and and referred to hereafter as “baryons”. These objects are highly classical and require no further gravitational consideration than has already been made in the literature. The inter-axion gravitational Poisson equation may be inverted using standard Greens function techniques

 Φa(→x)=−∫d3x′4πGρa(→x′)|→x−→x′| (24)

From the previous discussion of stress-energy smoothing it is determined that the axion density should be sampled on the Compton scale. The potential of a collection of Compton-limited axions is then

 Φa(→x)=−N∑i∫d3x′4πGρC(→xi)|→x−→xi| (25)

where the sum is over every axion and is the axion’s Compton-limited distribution profile. The result looks similar to electrons interacting electro-statically in simple atomic and molecular systems. The gravitational potential for an axion from other axions is then

 Φa,j=−N∑i,i≠j∫d3x′4πGρC(→xi)|→xj−→xi| (26)

where the contribution of the axion onto itself is expressly omitted. For the remainder of the paper we shall use the shorthand definition of the inter-axion gravitational potential

 ϕij=ϕ(→xi,→xj)=−∫d3x′4πGρC(→xi)|→xj−→xi| (27)

The form of the many-body Hamiltonian may now be rewritten

 ^H=N∑i−ℏ2∇22a2ma+12N∑i≠jϕij+N∑iΦ′(→xi) (28)

where the inter-axion potential sum is over both and indexes such that .

### 3.2 Solutions of the Many-Body Equation

A more detailed version of this section may be found in Appx. A, which outlines the context of the many-body bosonic Schrödinger equation, determines its implicit and explicit symmetries, and derives solutions. Presented here is a class of useful solutions to motivate later derivations.

The axion DM Hamiltonian is far from separable across single particle lines, making a Fock space representation and reduction of states untenable. In the case of no external potentials, there exists a more concise form to solutions of the interchange-symmetric Schrödinger equation, which centers on inter-particle correlators. The center-of-mass solutions of the system are then spanned by fully symmetric combinations of inter-particle correlators

 Ψ(x,t)∈sp⎛⎝{perm(∏ijψαij(→xi−→xj,t))}α∈A⎞⎠ (29)

where gives the linear span of a collection of functions, is the permanent operation, is the set of many-body solution indexes, is a solution to the single-body equation

 iℏ∂tψa=−ℏ22a2μ∇2(→xi−→xj)ψα+ϕ(→xi−→xj)ψa (30)

and is a particle reduced mass and scales as , thereby giving a naturally-extended potential length scale to the inter-particle correlations. The inter-particle potential’s length scale, the number of particles , and the global reach of permutation (anti-)symmetry are therefore key to the development of the system’s structure. This correlator form explores the solution space more efficiently than the single-particle expansion natural to Fock representations, especially in the presence of highly-correlated interactions. The utility of this form lies in the implicit embedding of exact exchange symmetry and the effects of that constraint on possible actions of the system; this is crucial in the expression of physics beyond the mean field. Only condensed configurations are considered for the remainder of this paper, where condensed here is taken to mean occupation of a single basis element of Exp. 29 with . More detailed considerations of condensation, including a generalized density-of-states derivation conforming to early axion dynamics, will appear in a later paper of this series.

### 3.3 Distillation to a Condensed Distribution Function

Condensing solutions of the many-body axion Schrödinger equation to the salient degrees of freedom is next. Applying the Runge-Gross theorem Runge and Gross (1984a, b), which proves the existence of an injective mapping from the potentials of a many-body quantum mechanical system to a single-body density when given an initial many-body wave-function, reveals that the only important degree of freedom for the identical system is its density, . The Runge-Gross theorem therefore reveals the possibility of reducing dimensionality without loss of generality of the axion state. Reduction of a cosmological volume of axions to a single body’s dynamics is therefore possible. Our treatment aims for such a description.

The approach pursued here phrases the density in terms of distribution functions (DFs), which track the density of system attributes through phase space. The many-body wave-function can be reduced into a phase-space pseudo-DF via a many-body Wigner transform. At the classical level (, the transform produces a DF governed by a collision-less Boltzmann equation

 ∂tf(N)+N∑i→pia2ma→∇if(N)−N∑i→∇iΦ′⋅→∇pif(N) −12N∑i≠j→∇iϕij⋅→∇pif(N)=O(ℏ) (31)

where is the N-body DF. See Appx. B for more detail on Wigner transforms and their properties. At this point neither the specific dynamics nor exchange conditions on the state are enforced, implying that such a Boltzmann equation may also by applied to analogous classes of fermionic or mixed systems.

Reducing the degrees of freedom to a single body is now a matter of integration. Projecting from the full N bodies to a two-body DF is straightforward, covered in detail in Appx. C, but reduction from two to one body is more complex due to the presence of two-body correlations. The normal “molecular chaos” approach to distribution theory would lead one to take a near-uncorrelated form of the distribution . The correlation function of this form is often restricted to scaling by local considerations of classical two-particle scattering Fokker (1914); Kolmogoroff (1931). However, the condensed state solutions from Sec. 3.2 implies no such scaling due to the global reach of the permutation condition. Instead, the extremal nature of the condensate is used to construct a correlation optimization problem. Defining the two-body correlation function by decomposing the two-body DF as

 f(2)(w1,w2,t)=~g×f(1)(w1,t)f(1)(w2,t), (32)

correlation functions of the extremal case of a condensed Bose fluid are found with minimal computation

 ~g=C−λ1f+1+λ2f+ (33)

where is the symmetrized single-body DF, and and are integrals of motion, and is constrained by the correlation present in the initial configuration; details are available in Appx. C. The resulting single-body equation of motion may now be written

 0=∂tf(1)+→pa2ma⋅→∇f(1)−→∇Φ′⋅→∇pf(1)−→∇¯Φ⋅→∇pf(1) −N−1N∫d6wj→∇ϕ1j⋅→∇pf(1)(w1,t)f(1)(wj,t) ×C−1−(λ1+λ2)f+1+λ2f+ (34)

where is the axion Newtonian gravitational potential, classically averaged

 ¯Φ=N∫d6w2ϕ12f(1) (35)

The extra-classical physics is a direct result of the exchange symmetry shared by each axion pair and the non-trivial correlations between axions; this influence is referred to as the exchange-correlation (XC) interaction. The effect of the XC pseudo-forces is almost classical, in that the removal of any one axion does little to change the result.

While compact, the form of Eqn. 34 is not simple to compute for many calculations in structure formation. Conveniently expanded forms of the new Boltzmann-like equation and other equations relevant for perturbative calculation are provided in Appx. D.

## 4 Early Axion Structures

This section investigates the presence of unique structure in early axion structure formation. “Early” here is taken to be synonymous with linearly perturbative over a homogeneous and isotropic background, such as in most radiation-dominated to matter-dominated recombination era models of SF.

Perturbative dynamics in the hydrodynamic limit are well recorded in the literature, including references Binney and Tremaine (2008) and Mo et al. (2010). The observables of fluid density and velocity are tracked in this demonstration, though higher moments are straightforward to compute. The density takes the classical definition of density perturbations in a co-moving frame

 ρ(→x)=M∫d3pf(→x,→p,t)=¯ρa(t)3(1+δ(→x,t)) (36)

where is the enclosed mass and is the average density over the Hubble volume. The fluid velocity takes on a similar definition

 ⟨→v(→x)⟩=1ρ∫d3pf(→x,→p,t)→p (37)

The fluid is also considered to be nearly static in the co-moving frame, with flows of order .

The linearized hydrodynamics of CDM are derived first. The co-moving CDM Boltzmann equation

 ∂tf(1)+→pma2⋅→∇f(1)−→∇Φ′⋅→∇pf(1)−→∇¯Φ⋅→∇pf(1)=0 (38)

is expanded in orders of density and velocity perturbations. Expectations for the zeroth and first moments of the distribution are found to have dynamics

 0th :∂tδ+1a→∇⋅{(1+δ)⟨→v⟩}=0 (39) 1st +1a→∇⋅{(1+δ)⟨→v⊗→v⟩}=0 (40)

giving density and momentum conservation laws respectively. An Euler-like form may be derived through a combination of Eqns. 39 and 40:

 ∂t⟨→v⟩+H⟨→v⟩+1a(⟨→v⟩⋅→∇)⟨→v⟩+1a→∇(Φ′+Φ) +1(1+δ)a→∇⋅{(1+δ)(⟨→v⊗→v⟩−⟨→v⟩⊗⟨→v⟩)}=0, (41)

where the last term disappears as cold dark matter has vanishingly small velocity dispersion. The fluctuation equation of motion is derived by combining the Euler form and matter conservation again, and to first order looks like a damped harmonic instability

 ∂2tδ+2H∂tδ=4πG¯ρ(δ+δ′)+O(δ2) (42)

where is the energy density fluctuation from complementary species. This gives the expected collapse on all concerned scales, in the limit of no back reactions.

Calculating perturbative Bose hydrodynamics follows very similarly. The lowest order XC effects modify the perceived depth of the self-gravitational well, which is analogous to modifying the observable cosmological dark matter density, but only within the DM species. The fluctuation equation of motion is derived to be

 −4πG¯ρλ1+λ22W3δ+O(λ22,δ2), (43)

where the s are evaluated at the initial spectrum and represents the phase space volume over which the s are calculated. The used expansion of the correlation function can be found in Appx. D. MFT is reached in the homogeneous limit, meaning to zeroth order and . Strictly speaking, the XC has no impact on first order dynamics as the non-trivial components of and the s enter at first order in perturbations. However, an inflated XC impact will be entertained here for the purposes of demonstration. A simple relation between the initial correlation and the constraints is found at first order in perturbations

 C−1=C1=λ1,1+λ2,1W. (44)

Past and present probes of structure are capable of determining the early power spectra and concentration of gravitating species with percent level precision (Spergel et al., 2007; Planck Collaboration et al., 2018). Bose dynamics with minor density fluctuations have nearly indistinguishable dynamics from CDM near the radiation-matter transition. The linear evolution of large scale structure formation from a Harrison-Zeldovich primordial power spectrum produces only weak deviations from CDM for our inflated condition on the initial correlation, Figs. 2, 3. Linear Bose structures are found to be consistent with the CDM structures to the one percent level for initial correlations in the range . Structure in the non-linear regime require a more study, which will be applied in the next entries of this series.

## 5 Discussion

The model derived in Section 3 displays a number of interesting properties. The extremal nature of an exact condensate reveals the nature of inter-axion correlation and exchange in a most natural way, producing new physics on all considered scales. The unique physics appears both in terms of high-order dynamical moments and as novel gravity-inspired potentials, all stemming from the correlation function . Solutions to the conformal Lagrange multipliers in turn depend on global invariants of the distribution, perhaps due to the unbounded nature of the extremal correlations and exchange symmetry. More realistic depictions of primordial axions are likely to ruin this exceptional case by introduction of new length scales and state transitions, though the current form may be considered to represent the best-case scenario for seeing physics beyond the dissociated limit.

The evolution of linear structures strictly reveals no departure of Bose DM from CDM infall. Deviations enter only at the second order in perturbation and reveal only mild divergences in dynamics and early structure when those liberties are taken; only marginal impacts to current isolated small-scale structures should be expected. This lack of scale dependence is intriguing, and again held to be related to the unbounded nature of the correlation and the non-local exchange symmetry within the region of interest. The marginal result of new structure is not surprising as the non-trivial portion of the correlation function is highly non-linear.

The final unresolved problem for the dynamics of condensed QCD axions is specifying the value of initial correlation. A complete model tracking of quantum correlations from pre-inflation to the current era is needed to make definitive statements about the presence of axion physics beyond the mean field. Such models, even for a subset of early cosmological history, are currently lacking in the literature.

## 6 Summary

The well-motivated axion DM candidate has properties that allow it to simultaneously form a Bose-Einstein condensate and a level of inter-particle correlation unavailable to most bosons studied in material science. This paper sought to derive a robust and computationally effective model of axion infall from first principles, and has done so in the condensed limit. The unique nature of such a maximally condensed axion fluid displays distinct departures from the standard collision-less fluid picture in the relevant case of self-gravitation. These cross-correlation impacts do not appear to significantly influence the outcome of linear structure formation, but may impact structures entering the non-linear regime. The next entries to this series introduce numerical simulation techniques to Bose gravitational collapse and explore the nuances of imperfect and dynamical axion condensation.

## 7 Acknowledgments

We would like to thank Jens Niemeyer, Katy Clough, David Marsh, Bodo Schwabe, Jan Veltmaat, and Xiaolong Du for their productive discussions in the refinement of this paper. We also gratefully acknowledge the support of the U.S. Department of Energy office of High Energy Physics and the National Science Foundation. TQ was supported in part by the NSF grant AST-1514868. EL and LR were supported in part by the DOE grant DE-SC0011665.

## References

• Peccei and Quinn (1977) R. D. Peccei and H. R. Quinn, Physical Review Letters 38, 1440 (1977).
• Weinberg (1978) S. Weinberg, Phys. Rev. Lett. 40, 223 (1978).
• Wilczek (1978) F. Wilczek, Phys. Rev. Lett. 40, 279 (1978).
• Mimasu and Sanz (2015) K. Mimasu and V. Sanz, Journal of High Energy Physics 6, 173 (2015)arXiv:1409.4792 [hep-ph] .
• Patrignani and Particle Data Group (2016) C. Patrignani and Particle Data Group, Chinese Physics C 40, 100001 (2016).
• Marsh (2016) D. J. E. Marsh, Phys. Rep. 643, 1 (2016)arXiv:1510.07633 .
• Raffelt (2008) G. G. Raffelt, “Astrophysical axion bounds,” in Axions: Theory, Cosmology, and Experimental Searches, edited by M. Kuster, G. Raffelt,  and B. Beltrán (Springer Berlin Heidelberg, Berlin, Heidelberg, 2008) pp. 51–71.
• Isern et al. (2010) J. Isern, E. García-Berro, L. G. Althaus,  and A. H. Córsico, A&A 512, A86 (2010)arXiv:1001.5248 [astro-ph.SR] .
• Córsico et al. (2012) A. H. Córsico, L. G. Althaus, A. D. Romero, A. S. Mukadam, E. García-Berro, J. Isern, S. O. Kepler,  and M. A. Corti, J. Cosmology Astropart. Phys. 12, 010 (2012)arXiv:1211.3389 [astro-ph.SR] .
• Viaux et al. (2013) N. Viaux, M. Catelan, P. B. Stetson, G. G. Raffelt, J. Redondo, A. A. R. Valcarce,  and A. Weiss, A&A 558, A12 (2013)arXiv:1308.4627 [astro-ph.SR] .
• Abbott and Sikivie (1983) L. F. Abbott and P. Sikivie, Physics Letters B 120, 133 (1983).
• Dine and Fischler (1983) M. Dine and W. Fischler, Physics Letters B 120, 137 (1983).
• Preskill et al. (1983) J. Preskill, M. B. Wise,  and F. Wilczek, Physics Letters B 120, 127 (1983).
• Kim (1979) J. E. Kim, Phys. Rev. Lett. 43, 103 (1979).
• Shifman et al. (1980) M. A. Shifman, A. I. Vainshtein,  and V. I. Zakharov, Nuclear Physics B 166, 493 (1980).
• Dine et al. (1981) M. Dine, W. Fischler,  and M. Srednicki, Physics Letters B 104, 199 (1981).
• Zhitnitsky (1980) A. R. Zhitnitsky, Sov. J. Nucl. Phys. 31, 260 (1980), [Yad. Fiz.31,497(1980)].
• 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, J. Clarke,  and ADMX Collaboration, Physical Review Letters 104, 041301 (2010)arXiv:0910.5914 [astro-ph.CO] .
• Arik et al. (2014) M. Arik, S. Aune, K. Barth, A. Belov, S. Borghi, H. Bräuninger, G. Cantatore, J. M. Carmona, S. A. Cetin, J. I. Collar, E. Da Riva, T. Dafni, M. Davenport, C. Eleftheriadis, N. Elias, G. Fanourakis, E. Ferrer-Ribas, P. Friedrich, J. Galán, J. A. García, A. Gardikiotis, J. G. Garza, E. N. Gazis, T. Geralis, E. Georgiopoulou, I. Giomataris, S. Gninenko, H. Gómez, M. Gómez Marzoa, E. Gruber, T. Guthörl, R. Hartmann, S. Hauf, F. Haug, M. D. Hasinoff, D. H. H. Hoffmann, F. J. Iguaz, I. G. Irastorza, J. Jacoby, K. Jakovčić, M. Karuza, K. Königsmann, R. Kotthaus, M. Krčmar, M. Kuster, B. Lakić, P. M. Lang, J. M. Laurent, A. Liolios, A. Ljubičić, G. Luzón, S. Neff, T. Niinikoski, A. Nordt, T. Papaevangelou, M. J. Pivovaroff, G. Raffelt, H. Riege, A. Rodríguez, M. Rosu, J. Ruz, I. Savvidis, I. Shilon, P. S. Silva, S. K. Solanki, L. Stewart, A. Tomás, M. Tsagri, K. van Bibber, T. Vafeiadis, J. Villar, J. K. Vogel, S. C. Yildiz,  and K. Zioutas (CAST Collaboration), Phys. Rev. Lett. 112, 091302 (2014).
• Valluri et al. (2016) M. Valluri, S. R. Loebman, J. Bailin, A. Clarke, V. P. Debattista,  and G. Stinson, in The General Assembly of Galaxy Halos: Structure, Origin and Evolution, IAU Symposium, Vol. 317, edited by A. Bragaglia, M. Arnaboldi, M. Rejkuba,  and D. Romano (2016) pp. 358–359, arXiv:1510.06006 .
• Binney (2017) J. Binney, ArXiv e-prints  (2017), arXiv:1706.01374 .
• Kolb and Tkachev (1993) E. W. Kolb and I. I. Tkachev, Physical Review Letters 71, 3051 (1993)hep-ph/9303313 .
• Tinyakov et al. (2016) P. Tinyakov, I. Tkachev,  and K. Zioutas, J. Cosmology Astropart. Phys. 1, 035 (2016)arXiv:1512.02884 .
• Fairbairn et al. (2017) M. Fairbairn, D. J. E. Marsh, J. Quevillon,  and S. Rozier, ArXiv e-prints  (2017), arXiv:1707.03310 .
• Banik and Sikivie (2015) N. Banik and P. Sikivie, ArXiv e-prints  (2015), arXiv:1501.05913 .
• Redondo (2013)
• Gorghetto et al. (2018) M. Gorghetto, E. Hardy,  and G. Villadoro, Journal of High Energy Physics 7, 151 (2018)arXiv:1806.04677 [hep-ph] .
• Vaquero et al. (2018) A. Vaquero, J. Redondo,  and J. Stadler, ArXiv e-prints  (2018), arXiv:1809.09241 .
• Arvanitaki et al. (2010) A. Arvanitaki, S. Dimopoulos, S. Dubovsky, N. Kaloper,  and J. March-Russell, Phys. Rev. D 81, 123530 (2010)arXiv:0905.4720 [hep-th] .
• Schive et al. (2014) H.-Y. Schive, T. Chiueh,  and T. Broadhurst, Nature Physics 10, 496 (2014)arXiv:1406.6586 .
• Veltmaat et al. (2018) J. Veltmaat, J. C. Niemeyer,  and B. Schwabe, ArXiv e-prints  (2018), arXiv:1804.09647 .
• Sikivie (1999a) P. Sikivie, Phys. Rev. D 60, 063501 (1999a)astro-ph/9902210 .
• Sikivie (1999b)
• Duffy and Sikivie (2008) L. D. Duffy and P. Sikivie, Phys. Rev. D 78, 063508 (2008)arXiv:0805.4556 .
• Chakrabarty and Sikivie (2017) S. S. Chakrabarty and P. Sikivie, in APS April Meeting Abstracts (2017).
• Banik et al. (2015) N. Banik, A. J. Christopherson, P. Sikivie,  and E. M. Todarello, Phys. Rev. D 91, 123540 (2015).
• 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] .
• Binney and Tremaine (2008) J. Binney and S. Tremaine, Galactic Dynamics: Second Edition, by James Binney and Scott Tremaine. ISBN 978-0-691-13026-2 (HB). Published by Princeton University Press, Princeton, NJ USA, 2008. (Princeton University Press, 2008).
• Misner et al. (1973) C. W. Misner, K. S. Thorne,  and J. A. Wheeler, San Francisco: W.H. Freeman and Co., 1973 (1973).
• Planck Collaboration et al. (2018) Planck Collaboration, N. Aghanim, Y. Akrami, M. Ashdown, J. Aumont, C. Baccigalupi, M. Ballardini, A. J. Banday, R. B. Barreiro, N. Bartolo, S. Basak, R. Battye, K. Benabed, J.-P. Bernard, M. Bersanelli, P. Bielewicz, J. J. Bock, J. R. Bond, J. Borrill, F. R. Bouchet, F. Boulanger, M. Bucher, C. Burigana, R. C. Butler, E. Calabrese, J.-F. Cardoso, J. Carron, A. Challinor, H. C. Chiang, J. Chluba, L. P. L. Colombo, C. Combet, D. Contreras, B. P. Crill, F. Cuttaia, P. de Bernardis, G. de Zotti, J. Delabrouille, J.-M. Delouis, E. Di Valentino, J. M. Diego, O. Doré, M. Douspis, A. Ducout, X. Dupac, S. Dusini, G. Efstathiou, F. Elsner, T. A. Enßlin, H. K. Eriksen, Y. Fantaye, M. Farhang, J. Fergusson, R. Fernandez-Cobos, F. Finelli, F. Forastieri, M. Frailis, E. Franceschi, A. Frolov, S. Galeotta, S. Galli, K. Ganga, R. T. Génova-Santos, M. Gerbino, T. Ghosh, J. González-Nuevo, K. M. Górski, S. Gratton, A. Gruppuso, J. E. Gudmundsson, J. Hamann, W. Handley, D. Herranz, E. Hivon, Z. Huang, A. H. Jaffe, W. C. Jones, A. Karakci, E. Keihänen, R. Keskitalo, K. Kiiveri, J. Kim, T. S. Kisner, L. Knox, N. Krachmalnicoff, M. Kunz, H. Kurki-Suonio, G. Lagache, J.-M. Lamarre, A. Lasenby, M. Lattanzi, C. R. Lawrence, M. Le Jeune, P. Lemos, J. Lesgourgues, F. Levrier, A. Lewis, M. Liguori, P. B. Lilje, M. Lilley, V. Lindholm, M. López-Caniego, P. M. Lubin, Y.-Z. Ma, J. F. Macías-Pérez, G. Maggio, D. Maino, N. Mandolesi, A. Mangilli, A. Marcos-Caballero, M. Maris, P. G. Martin, M. Martinelli, E. Martínez-González, S. Matarrese, N. Mauri, J. D. McEwen, P. R. Meinhold, A. Melchiorri, A. Mennella, M. Migliaccio, M. Millea, S. Mitra, M.-A. Miville-Deschênes, D. Molinari, L. Montier, G. Morgante, A. Moss, P. Natoli, H. U. Nørgaard-Nielsen, L. Pagano, D. Paoletti, B. Partridge, G. Patanchon, H. V. Peiris, F. Perrotta, V. Pettorino, F. Piacentini, L. Polastri, G. Polenta, J.-L. Puget, J. P. Rachen, M. Reinecke, M. Remazeilles, A. Renzi, G. Rocha, C. Rosset, G. Roudier, J. A. Rubiño-Martín, B. Ruiz-Granados, L. Salvati, M. Sandri, M. Savelainen, D. Scott, E. P. S. Shellard, C. Sirignano, G. Sirri, L. D. Spencer, R. Sunyaev, A.-S. Suur-Uski, J. A. Tauber, D. Tavagnacco, M. Tenti, L. Toffolatti, M. Tomasi, T. Trombetti, L. Valenziano, J. Valiviita, B. Van Tent, L. Vibert, P. Vielva, F. Villa, N. Vittorio, B. D. Wandelt, I. K. Wehus, M. White, S. D. M. White, A. Zacchei,  and A. Zonca, ArXiv e-prints  (2018), arXiv:1807.06209 .
• Lentz et al. (2018) E. W. Lentz, T. R. Quinn,  and L. J. Rosenberg, ArXiv e-prints  (2018), arXiv:1808.06378 .
• Runge and Gross (1984a) E. Runge and E. K. U. Gross, Phys. Rev. Lett. 52, 997 (1984a).
• Runge and Gross (1984b) E. Runge and E. K. U. Gross, Phys. Rev. Lett. 52, 997 (1984b).
• Fokker (1914) A. D. Fokker, Annalen der Physik 348, 810 (1914).
• Kolmogoroff (1931) A. Kolmogoroff, Mathematische Annalen 104, 415 (1931).
• Mo et al. (2010) H. Mo, F. C. van den Bosch,  and S. White, Galaxy Formation and Evolution, by Houjun Mo , Frank van den Bosch , Simon White, Cambridge, UK: Cambridge University Press, 2010 (2010).
• Spergel et al. (2007) D. N. Spergel, R. Bean, O. Doré, M. R. Nolta, C. L. Bennett, J. Dunkley, G. Hinshaw, N. Jarosik, E. Komatsu, L. Page, H. V. Peiris, L. Verde, M. Halpern, R. S. Hill, A. Kogut, M. Limon, S. S. Meyer, N. Odegard, G. S. Tucker, J. L. Weiland, E. Wollack,  and E. L. Wright, ApJS 170, 377 (2007)astro-ph/0603449 .
• Howlett et al. (2012) C. Howlett, A. Lewis, A. Hall,  and A. Challinor, J. Cosmology Astropart. Phys. 1204, 027 (2012)arXiv:1201.3654 [astro-ph.CO] .
• Zurek (2003) W. H. Zurek, eprint arXiv:quant-ph/0306072  (2003), quant-ph/0306072 .
• Zachos et al. (2005) C. K. Zachos, D. B. Fairlie,  and T. L. Curtright, Quantum mechanics in phase space , by Zachos, Cosmas.  Quantum mechanics in phase space, editors, Cosmas K. Zachos…[et al.]., Singapore; New Jersey, N.J.: World Scientific, c2005. (2005).
• Bondar et al. (2013) D. I. Bondar, R. Cabrera, D. V. Zhdanov,  and H. A. Rabitz, Phys. Rev. A 88, 263 (2013)arXiv:1202.3628 [quant-ph] .

## Appendix A The Many-Body Schrödinger Equation and One Special Class of Systems

Let us consider the problem of axion infall as a solution of the many-body Schrödinger equation under inter-particle Newtonian gravitation.

 iℏ∂tΨ(→x1,...,→xN;t)=^HΨ(→x1,...,→xN;t) (45)

with the Hamiltonian of the form

 ^H =^T+^Vint =−N∑iℏ2∇22a2ma+N∑i

where the external potential from the main text has been ignored. The domain of the Hamiltonian operator is contained in the Hilbert space , the time operator is spanned by the set of smooth complex functions, and together spans the domain of the Schrödinger operator. The product space is isomorphic to the -fold tensor product

 L2(R3)⊗...⊗L2(R3)⊗C10(R) (47)

as well as the N-body Fock domain

 (L2(R3)⊗C10(R))⊗...⊗(L2(R3)⊗C10(R))/T (48)

where is the equal-time equivalence relation among all elements of the product. Being contained in the Fock domain implies any solution of the many-body Schrödinger equation may be written as combinations of single-body products to arbitrary accuracy

 Ψ=∏i∈{1,...,N}α∈Aψα(→xi)+... (49)

Finding these combinations can be difficult. With such a strongly correlated force as gravity, these linear combinations are often of infinite length and slowly converging, making computation difficult when faced with finite resources. As a goal of the axion structure formation research project is to succinctly describe the dynamics of the axion system, an additional tool is needed.

As the axions are bosons, the commutation relations of the quantum field requires the system to only occupy fully symmetric states under particle exchange

 Ψ(→x1,...,→xi,...,→xj,...,→xN;t) =Ψ(→x1,...,→xj,...,→xi,...,→xN;t) ∀i,j (50)

which further refines the solution domain to

 (L2(R3)⊗C10(R))⊗...⊗(L2(R3)⊗C10(R))/(T∘SN) (51)

where is the equivalence among the symmetric group of axion permutations.

To find a more useful and concise form of the many-body solutions, let us rewrite the Hamiltonian in a form that showcases the unique shape of the potential. The kinetic term describes the energy due to motion of axions relative to a prescribed frame of reference, but motion of the collection of particles may just as well be parameterized by the motion of the system’s center of mass and the relative particle motions. An appropriate gradient re-parameterization makes this clear

 →∇i =1NN∑j(→∇j+(→∇i−→∇j)) =1N→∇com+1N∑i≠j→∇Δij (52)

where , , and , making

 ∇2i =1N2∇2com+1N2∑i≠j∇2Δij (53) +2N2∑i≠j→∇com⋅→∇Δij+1N2∑k≠j→∇Δij⋅→∇Δik

The alternate form of can then be written as

 ^T=−ℏ22a2Mtot∇2com−∑i

where . The second to last term from Eqn. 53 is found to vanish from as each may be matched with a , and the sum of the last term in Eqn. 53 is found to contribute half of the inter-axion energy. The many-body Hamiltonian can then be rewritten

 ^H =−ℏ22a2Mtot∇2com−∑i

and can be separated into a center-of-mass Hamiltonian and an interacting Hamiltonian

 ^Hcom =−12Mtot∇2com (56) ^Hint =−∑i

Note that appears to be separable over the , though the s are non-independent over the many-body solution space. To circumvent this, let use consider a slightly different domain: a system of ‘particles’ with Hamiltonian defined over the function space with an analogous Bose constraint

 L2(R3N(N−1))⊗C10(R)/SN(N−1) (58)

and all the implications thereof. In this case the Hamiltonian is separable and the stationary states are spanned by single permanents of single-‘particle’ stationary states

 Φ′(Δij)∈sp⎛⎝{perm(∏ijϕαij(→xi−→xj,t))}α∈A⎞⎠ (59)

and the pure states as

 Ψ′(Δij,t)∈sp⎛⎝{perm(∏ijψαij(→xi−→xj,t))}α∈A⎞⎠ (60)

where gives the linear span of a collection of functions, is the permanent operation, is the set of many-body solution indexes, is a solution to the single-body equation

 iℏ∂tψa=−ℏ22a2μ∇2(→xi−→xj)ψα+ϕ(→xi−→xj)ψa (61)

Solutions of this space may be used in the original use-case of axions by applying an appropriate projection

 p: Δij=→xi−→xj ∀i,j (62)

on the Hamiltonian and function space pair.

 (^H(R3N(N−1)+1),C(R3N(N−1)+1)) →(^H(R3N+1),C(R3N+1)) (63)

It is apparent that this projection both respects the modulated many-body solution space and the N-body Schrödinger equation. It can also be shown that the span of solutions to the -body problem can be represented by projected solutions of the -body problem, implying is dense in .

## Appendix B The Many-Body Wigner Function and its Super-de Broglie Limit

This is a short review of the many-body Wigner transform, some of its qualities, and its equation of motion in the super-de Broglie limit.

### b.1 Wigner Transform

Phase space DFs are often far better suited for tracking system observables than wave-functions, due to the hyperbolicity of their governing dynamics. The many-body Bose system can be reduced into a phase-space pseudo-DF via the many-body Wigner function

 f(N)(→x1,→p1,...,→xN,→pN,t) =∫d3x′1⋅⋅⋅d3x′Nei∑Nj→pj⋅→x′j/ℏ× Ψ†(→x1+→x′1/2,...,→xN+→x′N/2,t)× Ψ(→x1−→x′1/2,...,→xN−→x′N/2,t) (64)

For brevity, the arguments of and are understood to be and respectively. While not a true DF, in that can take on negative values, it does share a number of properties with true DFs, such as being real-valued

 f(N)(→x,→p)∗ =∫d3x′e−i∑Nj→pj⋅→x′j/ℏ(Ψ†Ψ)∗ =∫d3x′e−i∑Nj→pj⋅→x′j/ℏΨΨ† =∫d3x′ei∑Nj→pj⋅→x′j/ℏΨ†Ψ, (→x′→−→x′) =f(N)(→x,→p), (65)

and providing properly-normalized probability densities over projections to 3D single-body spaces

 ρ(→x) =∫d3p1d6(N−1)wf(N), (66) ρ(→p) =∫d3x1d6(N−1)wf(N). (67)

The Wigner pseudo-DF can also be expected to provide a reliable DF in certain circumstances (Zurek, 2003; Zachos et al., 2005; Bondar et al., 2013), namely on scales much larger than the de Broglie wavelength. The differences between a Wigner function and a true DF stem from the uncertainty principles present in the underlying wave-function/quantum-state representation. Complementary parameters such as position and momentum cannot be fully articulated simultaneously when quantum limited, but we see soon that the two become independent at the effective structure formation scales, far above the de Broglie length. We choose the DF form as it fully captures the dynamics of diffusive quantum systems and is most convenient for N-Body simulation, a numerical technique approached later in this series.

### b.2 Super de-Broglie Equation of Motion

The equation of motion for the distribution is found by straightforward differentiation

 ∂tf(N)(→x1,→p1,...,→xN,→pN,t) =∫d3x′1⋅⋅⋅d3x′Nei∑Nj→pj⋅→x′j/ℏ(∂tΨ†Ψ+Ψ†∂tΨ) (68)

Judicious substitution of the Schrödinger equation transforms the expression into an equation of motion

 ∂tf(N)(→x1,→p1,...,→xN,→pN,t) =1iℏ∫d3x′1⋅⋅⋅d3x′Nei∑Nj→pj⋅→x′j/ℏ(−^HΨ†Ψ+Ψ†^HΨ) (69)

where we use the Hamiltonian of Eqn.28.

The kinetic contribution can be simplified through successive application of the product rule of differentiation and consolidation of terms, and as expected provides transport terms

 K=1iℏℏ22a2mN∑i∫d3x′1⋅⋅⋅d3x′Nei∑Nj→pj⋅→x′j/ℏ× (∇2(→xi+→x′i/2)Ψ†Ψ−Ψ†∇2(→xi−→x′i/2)Ψ) =−N∑i→pia2m→∇if(N)+O(ℏ) (70)

where de Broglie wavelength contributions and higher order are inconsequential. The quantum diffusion scales of the QCD axion are many orders of magnitude smaller than the desired scale of simulation, causing them to average out to zero over cosmological coarse graining.

The external potentials give forcing contributions, calculated through Taylor expansion of the potentials about :

 Fext =1iℏN∑i∫d3x′1⋅⋅⋅d3x′Nei∑Nj→pj⋅→x′j/ℏ× (Φ′(→xi+→x′i/2)Ψ†Ψ−Ψ†Φ′(→xi−→x′i/2)Ψ) =N∑i→∇iΦ′⋅→∇pif(N)+O(ℏ2) (71)

Lastly, the two-body interaction potentials also contribute a forcing term

 Fint=12iℏN∑i≠j∫d3x′1⋅⋅⋅d3x′Nei∑Nj→pj⋅→x′j/ℏ× (ϕij(→xi+→x′i/2,→xj+→x′j/2)Ψ†Ψ −Ψ†ϕij(→xi−→x′i/2,→xj−→x′j/2)Ψ) =N∑i≠j→∇iϕij⋅→∇