Axion Structure Formation I: The Co-motion Picture
MMNRAS , 1– ?? (2018) Preprint 23 October 2018 Compiled using MNRAS L A TEX style file v3.0
Axion Structure Formation I: The Co-motion Picture
Erik W. Lentz, , (cid:63) Thomas R. Quinn, Leslie J Rosenberg Institut f¨ur Astrophysik, Georg-August Universitat G¨ottingen, G¨ottingen, Deutschland 37707; Department of Astronomy, University of Washington, Seattle, WA, USA 98195-1580; Department of Physics, University of Washington, Seattle, WA, USA 98195-1580;
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
Axions as dark matter is an increasingly important subject in astrophysics andcosmology. Experimental and observational searches are mounting across the massspectrum of axion-like particles, many of which require detailed knowledge of axionstructure over a wide range of scales. Current understanding of axion structures isfar from complete, however, due largely to controversy in modeling the candidate’shighly-degenerate state. The series Axion Structure Formation seeks to develop a con-sistent model of QCD axion dark matter dynamics that follows their highly-degeneratenature to the present using novel modeling techniques and sophisticated simulations.This inaugural paper presents the problem of describing many non-relativistic axionswith minimal degrees of freedom and constructs a theory of axion infall for the limit ofcomplete condensation. The derived model is shown to contain axion-specific dynam-ics not unlike the exchange-correlation influences experienced by identical fermions.Perturbative calculations are performed to explore the potential for imprints in earlyuniverse structures.
Key words: cosmology: dark matter – cosmology: early Universe – galaxies: forma-tion – galaxies: haloes – galaxies: structure
The QCD (quantum chromodynamic) axion is a well-motivated candidate to solve problems in both particlephysics and cosmology. The theory of axions originated in1977 as a result of a scalar-induced axial symmetry over theQCD sector, used to solve the strong (QCD) CP problem ofparticle physics (Peccei and Quinn 1977). The axion particlewas proposed the following year from a spontaneous break-ing of that axial symmetry, solving the strong CP problemdynamically for energies about and below λ QCD (cid:38)
200 MeV(Weinberg 1978; Wilczek 1978). Experiments and astrophys-ical observations have replaced this Peccei-Quinn-Weinberg-Wilczek axion (Mimasu and Sanz 2015; Patrignani and Par-ticle Data Group 2016) with more complex and unified mod-els, pushing the breaking scale ever higher and the axionmass ever lower. These new theories also have increasinglyaddressed cosmology (Marsh 2016), where the dark matter(DM) problem has been gaining in prominence.As a DM candidate, the QCD axion is highly attrac-tive due to its well-bounded parameter space of mass andcouplings to the standard model (SM) of particle physics,Fig. 1. Starting at milli-eV masses, there is a bound above (cid:63)
E-mail: [email protected] which the axion would have been seen in various astrophys-ical processes (Raffelt 2008; Isern et al. et al. et al. et al. et al. et al. c (cid:13) a r X i v : . [ a s t r o - ph . C O ] O c t E. W. Lentz et al. [] -6 -5 -4 -3 -2 -1 mass [eV/c ] -17 -16 -15 -14 -13 -12 -11 -10 -9 -8 -7 C o u p li n g | g a γγ | [ G e V − ] K S V Z D F S Z S N A CASTADMX Canonical Galactic Degeneracy
Figure 1.
The parameter space of axion-like particles over theplausible QCD axion DM region, with benchmark theories (red),cosmological constraints (orange-to-brown gradient, in order ofseverity) and degeneracy (blue gradient). Experimental and ob-servational constraints are also shown from ADMX (Asztalos et al. et al. M (cid:12) / kpc for the local galactic DM density and 300 km/s for the velocitydispersion. ephemeral memory of gaseous and stellar motions, thoughsome progress has been made in recent years (Valluri et al. f a , and the inflation scale f I = T I / π where T I is the as-sociated cosmological radiation temperature. In brief, if thePQ symmetry is unbroken during inflation, meaning that f a < T I / π , then the axion’s parent scalar field has zerovacuum expectation value (VEV), leaving the angular ax-ion field randomly distributed over each Hubble patch whenthe f a is reached after inflation. Such a distribution resultsin the large average misalignment angle (cid:10) φ a,i (cid:11) = π /
3. Ifthe PQ symmetry is broken during inflation, φ PQ varies onscales above the Hubble volume post-inflation. This effec- tively makes φ PQ uniform in the visible universe and closeto a free parameter of the model.Issues occur with each inflation possibility. The random-ness of φ PQ in the pre-PQ breaking inflation case gives riseto topological defects and seeds other structures that mayprove difficult for local searches such as haloscope experi-ments, as much of the mass may exist in compact objectssuch as mini-clusters (Kolb and Tkachev 1993; Tinyakov et al. et al. O (100 GeV ) for micro-eVaxions, far above the QCD axion mass (Banik and Sikivie2015). It is reasonable to assume that the axions inhabitsome form of Bose-Einstein condensate (BEC) from the ex-pected degeneracy alone.Axion DM dynamics go through several stages after cre-ation. Radiation gravity and the potential of the axion fielddominate axion interactions early after creation, respectivelyequilibrating with the cosmological radiation perturbationsor inducing violent collapse in areas of large over-density(Redondo 2013; Gorghetto et al. et al.
MNRAS , 1– ?? (2018) xion Structure Formation I on relic axions is straightforward and attractive. In the limitof ultra-light axions (ULAs), m a ∼ ( few ) × − eV, the deBroglie wavelength scale enters the galactic scale (Arvanitaki et al. et al. et al. D − catastrophes (Sikivie 1999a,b;Duffy and Sikivie 2008; Chakrabarty and Sikivie 2017). Itis also claimed that our own solar system is close to one ofthese caustics of a fitted MW axion halo. Cold flows feed-ing into the outermost rings would contribute to the localaxion energy spectra, significantly influencing the results ofterrestrial search experiments. Proponents of caustic ringsappreciate the compactness of the MFT description, but be-lieve that it does not adequately maintain the quantum na-ture of the degenerate Bose system. The fault with the MFTsimplification is reported to lie in the mis-estimation of therelevant timescales of gravitational equilibration, leading toan imbalance in state occupation (Banik et al. et al. et al. Post-QCD phase transition axions are well-described by asingle quantum scalar field operator ( ˆ φ ) defined over a dy-namic 3+1 pseudo-Riemannian manifold ( M ) governedby classical Einstein’s equations. Such a full description isextremely accurate but excessive for common cosmologicalstructures, so several simplifications are made in the inter-ests of computability. On the largest scales, the geometry istaken to be nearly homogeneous, isotropic, while the smallerscales of interest are taken to satisfy weak-field conditions(Misner et al. et al. H − Λ c πG ρ, (1)2 ¨ aa + (cid:18) ˙ aa (cid:19) − Λ c = − πGc ¯ p, (2)where ¯ ρ is the total expected mean co-moving energy den-sity, ¯ p is the total expected mean pressure, a ( t ) is the FLRWscale factor, H = ˙ a/a is the Hubble flow, and the dots repre-sent derivatives in physical time as measured by a co-movingobserver. The next order terms obey a parabolic potential-like equation.∆Φ = 4 πGac ρ (1) (3)where ∆ = (cid:126) ∇ is the flat spatial Laplacian, ρ (1) is the devi-ation in the frame’s energy density from the mean, and the MNRAS , 1– ????
MNRAS , 1– ?? (2018) xion Structure Formation I on relic axions is straightforward and attractive. In the limitof ultra-light axions (ULAs), m a ∼ ( few ) × − eV, the deBroglie wavelength scale enters the galactic scale (Arvanitaki et al. et al. et al. D − catastrophes (Sikivie 1999a,b;Duffy and Sikivie 2008; Chakrabarty and Sikivie 2017). Itis also claimed that our own solar system is close to one ofthese caustics of a fitted MW axion halo. Cold flows feed-ing into the outermost rings would contribute to the localaxion energy spectra, significantly influencing the results ofterrestrial search experiments. Proponents of caustic ringsappreciate the compactness of the MFT description, but be-lieve that it does not adequately maintain the quantum na-ture of the degenerate Bose system. The fault with the MFTsimplification is reported to lie in the mis-estimation of therelevant timescales of gravitational equilibration, leading toan imbalance in state occupation (Banik et al. et al. et al. Post-QCD phase transition axions are well-described by asingle quantum scalar field operator ( ˆ φ ) defined over a dy-namic 3+1 pseudo-Riemannian manifold ( M ) governedby classical Einstein’s equations. Such a full description isextremely accurate but excessive for common cosmologicalstructures, so several simplifications are made in the inter-ests of computability. On the largest scales, the geometry istaken to be nearly homogeneous, isotropic, while the smallerscales of interest are taken to satisfy weak-field conditions(Misner et al. et al. H − Λ c πG ρ, (1)2 ¨ aa + (cid:18) ˙ aa (cid:19) − Λ c = − πGc ¯ p, (2)where ¯ ρ is the total expected mean co-moving energy den-sity, ¯ p is the total expected mean pressure, a ( t ) is the FLRWscale factor, H = ˙ a/a is the Hubble flow, and the dots repre-sent derivatives in physical time as measured by a co-movingobserver. The next order terms obey a parabolic potential-like equation.∆Φ = 4 πGac ρ (1) (3)where ∆ = (cid:126) ∇ is the flat spatial Laplacian, ρ (1) is the devi-ation in the frame’s energy density from the mean, and the MNRAS , 1– ???? (2018) E. W. Lentz et al.
Newtonian potential is related to the co-moving geometryvia Φ = − h / h is the diagonal time component ofthe co-moving metric perturbation from FLRW. Relativisticgravity as a classical theory can be interpreted as averagingover quantum length scales. Free fundamental particles areassociated with two quantum length scales, Compton andde Broglie, though only Compton is covariantly meaningful;the de Broglie lengths parameterize frame-dependent cor-rections from wave kinematics. As the fundamental lengthscale of a particle, the Compton length determines the limitby which a light-like vector field may distinguish space-likestructure. The Compton scale is on the order of centimetersfor QCD axion DM. This is the smoothing length that shallbe used to evaluate the axion gravitational potential.The effective axion equations of motion are also derivedfrom a fully covariant description. Describing the axion sys-tem as a second-quantized operator, the axion field’s govern-ing equations are characterized by a Lagrangian operator ofthe formˆ L = (cid:16) ∇ µ ˆ φ † (cid:17) g µν (cid:16) ∇ ν ˆ φ (cid:17) − V ( ˆ φ ) (4)where ∇ µ is a covariant derivative over M and V is theaxion interaction potential. The axion field operator is ofBose type, meaning that it commutes under particle ex-change. The field obeys the Klein-Gordon equation over ef-fective fields in the semi-classical limit ∇ µ ∇ µ ˆ φ + (cid:16) m a (cid:126) (cid:17) ˆ φ − λ | ˆ φ | ˆ φ = 0 (5)where the fourth order expansion of the axion potential isused instead of the full cosine axion potential V ( ˆ φ ) = m a (cid:126) | ˆ φ | − λ | ˆ φ | + O ( | ˆ φ | ) (6)and λ = ∝ f − a is the coupling strength of the four-fieldcoupling strength. Evaluating the covariant derivatives inthe linear approximation, the derivative term of the Klein-Gordon operator becomes g µν ∇ µ ∇ ν ˆ φ = a − (cid:3) ˆ φ + 3 a − H∂ t ˆ φ − a − (cid:0) h − (cid:1) µν ∂ µ ∂ ν ˆ φ − a − H (cid:0) h − (cid:1) ν ∂ ν ˆ φ − a − h∂ ˆ φ + 12 a − ∂ ν ( h ) ∂ ν ˆ φ + O ( h ) (7)where (cid:3) = − ∂ t + ∆ is the flat d’Alembert operator, and allindex contraction, raising and lowering, are performed usingthe Minkowski metric ( η ) in agreement with the previously-mentioned weak field formalism (Misner et al. φ = (cid:126) √ m (cid:16) ˆ ψe − im a t/ (cid:126) + ˆ ψ ∗ e im a t/ (cid:126) (cid:17) (8)Substituting this form of the field operator, and taking theNewtonian limit for gravity, h = − ψ ∗ (cid:18) ∂ t − H (cid:19) ˆ ψ + ˆ ψ ∗ (cid:126) (cid:126) ∇ a m a ˆ ψ + m a Φ | ˆ ψ | − λ (cid:126) m a | ˆ ψ | + { h.c. } + (quickly oscillating terms)= O ( h , m / (cid:126) , mh/ (cid:126) ) (9) The quickly oscillating terms are ignored as their actionsevaluate to 0 over first-order timescales. The terms of theeffective axion action can be easily related to the terms of thefield Hamiltonian operator as the action measure includes afactor of √− g = a + O ( h ):ˆ H f = ˆ K + ˆ U + ˆ V (10)ˆ K = − ˆ ψ ∗ (cid:126) (cid:126) ∇ a m a ˆ ψ (11)ˆ U = λ | ˆ ψ | (12)ˆ V = ˆ ψ ∗ m a Φ ˆ ψ (13)which confirms that Φ is indeed the canonical Newtoniangravitational potential.The many-body axion Schr¨odinger operator equationcan now be written in the comoving parameterization i (cid:126) ∂ t ˆ ψ = − (cid:126) (cid:126) ∇ a m a ˆ ψ + λ | ˆ ψ | ˆ ψ + m a Φ ˆ ψ (14)This is the operator equation governing the quantum axionfield in the presence of classical gravity, which serves as arobust starting point for much of axion structure-formationmodels after the QCD phase transition.Interactions during the remainder of the radiation eraare dominated by gravitational perturbations from relativis-tic species and the axion point interaction. Structure forma-tion of axions are linear so long as axion perturbations re-main small during appreciable self-interaction. After a timethe axions diffuse beyond appreciable self-interaction, themean free path exceeds the Hubble radius, and their num-ber becomes a conserved quantity.The expansion rate and diffuse gravitational potentialskeep the axions in the perturbative regime during the re-mainder of radiation domination. The field operator pictureis still sensible during this time of point-wise and separablepotentials. Even so, there is some sub-dominant but increas-ingly important physics at work. Inter-axion gravitation con-nects every axion to every other axion via the long-rangeun-screenable Coulomb potential V grav ( (cid:126)x i , (cid:126)x j ) ∼ = − πGm a | (cid:126)x i − (cid:126)x j | (15)for well-separated axions in this limit of gravity. Such ahighly-correlated interaction is thought to maintain theaxion condensate through “gravitational thermalization”(Erken et al. MNRAS , 1– ?? (2018) xion Structure Formation I describes a new model of axion SF that more fully exploresthe state space. This section follows similar path to Lentz et al. (2018). Em-bellishment on the finer points of Schr¨odinger and Wignermethods in the context of axion self-gravitation are given inAppxs. A, B, C.
While in principle correct, the second quantization formal-ism used by many axion MFTs is not well-suited for ax-ion structure formation. It encourages a separable decom-position among axions not supported by their dynamics,and supplies a superfluous freedom to create/annihilatesaid separated particles. The technique presented here in-stead restricts its scope by falling back to first quantiza-tion, made possible by conserved axion number. Rephrasingthe field theory into a quantum mechanical approach, theSchr¨odinger picture of Hilbert space uncovers a many-bodywave-tracing evolution equation. Evaluating the dynamicalequation over an arbitrary element of the many-body solu-tion space gives the form i (cid:126) ∂ t Ψ( (cid:126)x , ..., (cid:126)x m ; t ) = ˆ H Ψ( (cid:126)x , ..., (cid:126)x m ; t ) (16)where Ψ is the many-body axion wave-function and ˆ H is themany-body Hamiltonian operator on wave-function space.Note that the symbol for the first-quantized many-bodyHamiltonian operator ˆ H is similar but distinct from theHubble flow H . The wave-function also inherits an exact ex-change symmetry due to the commuting nature of the axionfield operator.Ψ ( (cid:126)x , ..., (cid:126)x i , ..., (cid:126)x j , ..., (cid:126)x N ; t ) = +Ψ ( (cid:126)x , ..., (cid:126)x j , ..., (cid:126)x i , ..., (cid:126)x N ; t ) ∀ i, j (17)The many-body Hamiltonian ˆ H may be broken down intorecognizable kinetic and potential energy partsˆ H = ˆ T + ˆ V = M (cid:88) i − (cid:126) ∇ i a m a + (cid:88) i Φ i (18)where Φ i is the gravitational potential on the i-th axion.The gravitational interaction requires further clarifica-tion. Strictly speaking the gravitational potential may comefrom multiple massive and radiative species. Consideringonly massive species like axion DM, stars, atomic gas, etc.,the potential may be broken up linearly, using the Coulombgauge choice where the potential vanishes when infinitely farfrom all sources: ρ = ρ a + ρ (cid:48) and (19) ∇ Φ = 4 πGρ = 4 πG (cid:0) ρ a + ρ (cid:48) (cid:1) (20)impliesΦ = Φ a + Φ (cid:48) , (21) where ∇ Φ a = 4 πGρ a (22) ∇ Φ (cid:48) = 4 πGρ (cid:48) (23)and all potential terms not sourced by axionic densitiesare consolidated into ρ (cid:48) and Φ (cid:48) and referred to hereafteras “baryons”. These objects are highly classical and requireno further gravitational consideration than has already beenmade in the literature. The inter-axion gravitational Poissonequation may be inverted using standard Greens functiontechniquesΦ a ( (cid:126)x ) = − (cid:90) d x (cid:48) πGρ a ( (cid:126)x (cid:48) ) | (cid:126)x − (cid:126)x (cid:48) | (24)From the previous discussion of stress-energy smoothing it isdetermined that the axion density should be sampled on theCompton scale. The potential of a collection of Compton-limited axions is thenΦ a ( (cid:126)x ) = − N (cid:88) i (cid:90) d x (cid:48) πGρ C ( (cid:126)x i ) | (cid:126)x − (cid:126)x i | (25)where the sum is over every axion and ρ C is the axion’sCompton-limited distribution profile. The result looks simi-lar to electrons interacting electro-statically in simple atomicand molecular systems. The gravitational potential for anaxion from other axions is thenΦ a,j = − N (cid:88) i,i (cid:54) = j (cid:90) d x (cid:48) πGρ C ( (cid:126)x i ) | (cid:126)x j − (cid:126)x i | (26)where the contribution of the axion onto itself is expresslyomitted. For the remainder of the paper we shall use theshorthand definition of the inter-axion gravitational poten-tial φ ij = φ ( (cid:126)x i , (cid:126)x j ) = − (cid:90) d x (cid:48) πGρ C ( (cid:126)x i ) | (cid:126)x j − (cid:126)x i | (27)The form of the many-body Hamiltonian may now be rewrit-tenˆ H = N (cid:88) i − (cid:126) ∇ a m a + 12 N (cid:88) i (cid:54) = j φ ij + N (cid:88) i Φ (cid:48) ( (cid:126)x i ) (28)where the inter-axion potential sum is over both i and j indexes such that i (cid:54) = j . A more detailed version of this section may be found inAppx. A, which outlines the context of the many-bodybosonic Schr¨odinger equation, determines its implicit andexplicit symmetries, and derives solutions. Presented here isa class of useful solutions to motivate later derivations.The axion DM Hamiltonian is far from separable acrosssingle particle lines, making a Fock space representation andreduction of states untenable. In the case of no external po-tentials, there exists a more concise form to solutions of theinterchange-symmetric Schr¨odinger equation, which centerson inter-particle correlators. The center-of-mass solutions of
MNRAS , 1– ????
MNRAS , 1– ???? (2018) E. W. Lentz et al. the system are then spanned by fully symmetric combina-tions of inter-particle correlatorsΨ ( x, t ) ∈ sp (cid:40) perm (cid:32)(cid:89) ij ψ α ij ( (cid:126)x i − (cid:126)x j , t ) (cid:33)(cid:41) α ∈ A (29)where sp () gives the linear span of a collection of functions, perm () is the permanent operation, A is the set of many-body solution indexes, ˆ ψ a is a solution to the single-bodyequation i (cid:126) ∂ t ψ a = − (cid:126) a µ ∇ ( (cid:126)x i − (cid:126)x j ) ψ α + φ ( (cid:126)x i − (cid:126)x j ) ψ a (30)and µ is a particle reduced mass and scales as 1 /N ,thereby giving a naturally-extended potential length scaleto the inter-particle correlations. The inter-particle poten-tial’s length scale, the number of particles N , and the globalreach of permutation (anti-)symmetry are therefore key tothe development of the system’s structure. This correlatorform explores the solution space more efficiently than thesingle-particle expansion natural to Fock representations,especially in the presence of highly-correlated interactions.The utility of this form lies in the implicit embedding ofexact exchange symmetry and the effects of that constrainton possible actions of the system; this is crucial in the ex-pression of physics beyond the mean field. Only condensedconfigurations are considered for the remainder of this pa-per, where condensed here is taken to mean occupation ofa single basis element of Exp. 29 with α ij = α . More de-tailed considerations of condensation, including a general-ized density-of-states derivation conforming to early axiondynamics, will appear in a later paper of this series. Condensing solutions of the many-body axion Schr¨odingerequation to the salient degrees of freedom is next. Apply-ing the Runge-Gross theorem Runge and Gross (1984a,b),which proves the existence of an injective mapping from thepotentials of a many-body quantum mechanical system to asingle-body density when given an initial many-body wave-function, reveals that the only important degree of freedomfor the identical system is its density, ρ . The Runge-Grosstheorem therefore reveals the possibility of reducing dimen-sionality without loss of generality of the axion state. Reduc-tion of a cosmological volume of axions to a single body’sdynamics is therefore possible. Our treatment aims for sucha description.The approach pursued here phrases the density in termsof distribution functions (DFs), which track the densityof system attributes through phase space. The many-bodywave-function can be reduced into a phase-space pseudo-DF via a many-body Wigner transform. At the classicallevel ( (cid:126) → ∂ t f ( N ) + N (cid:88) i (cid:126)p i a m a (cid:126) ∇ i f ( N ) − N (cid:88) i (cid:126) ∇ i Φ (cid:48) · (cid:126) ∇ p i f ( N ) − N (cid:88) i (cid:54) = j (cid:126) ∇ i φ ij · (cid:126) ∇ p i f ( N ) = O ( (cid:126) ) (31) where f ( N ) is the N-body DF. See Appx. B for more detailon Wigner transforms and their properties. At this pointneither the specific dynamics nor exchange conditions on thestate are enforced, implying that such a Boltzmann equationmay also by applied to analogous classes of fermionic ormixed systems.Reducing the degrees of freedom to a single body isnow a matter of integration. Projecting from the full N bod-ies to a two-body DF is straightforward, covered in detailin Appx. C, but reduction from two to one body is morecomplex due to the presence of two-body correlations. Thenormal “molecular chaos” approach to distribution theorywould lead one to take a near-uncorrelated form of the dis-tribution f (2) = f (1) ⊗ f (1) + g . The correlation function ofthis form is often restricted to 1 /N scaling by local consid-erations of classical two-particle scattering Fokker (1914);Kolmogoroff (1931). However, the condensed state solutionsfrom Sec. 3.2 implies no such scaling due to the global reachof the permutation condition. Instead, the extremal natureof the condensate is used to construct a correlation optimiza-tion problem. Defining the two-body correlation function bydecomposing the two-body DF as f (2) ( w , w , t ) = ˜ g × f (1) ( w , t ) f (1) ( w , t ) , (32)correlation functions of the extremal case of a condensedBose fluid are found with minimal computation˜ g = C − λ f + λ f + (33)where f + is the symmetrized single-body DF, and λ and λ are integrals of motion, and C is constrained by the correla-tion present in the initial configuration; details are availablein Appx. C. The resulting single-body equation of motionmay now be written0 = ∂ t f (1) + (cid:126)pa m a · (cid:126) ∇ f (1) − (cid:126) ∇ Φ (cid:48) · (cid:126) ∇ p f (1) − (cid:126) ∇ ¯Φ · (cid:126) ∇ p f (1) − N − N (cid:90) d w j (cid:126) ∇ φ j · (cid:126) ∇ p f (1) ( w , t ) f (1) ( w j , t ) × C − − ( λ + λ ) f + λ f + (34)where ¯Φ is the axion Newtonian gravitational potential, clas-sically averaged¯Φ = N (cid:90) d w φ f (1) (35)The extra-classical physics is a direct result of the exchangesymmetry shared by each axion pair and the non-trivial cor-relations between axions; this influence is referred to as theexchange-correlation (XC) interaction. The effect of the XCpseudo-forces is almost classical, in that the removal of anyone axion does little to change the result.While compact, the form of Eqn. 34 is not simple tocompute for many calculations in structure formation. Con-veniently expanded forms of the new Boltzmann-like equa-tion and other equations relevant for perturbative calcula-tion are provided in Appx. D. This section investigates the presence of unique structurein early axion structure formation. “Early” here is taken to
MNRAS , 1– ?? (2018) xion Structure Formation I be synonymous with linearly perturbative over a homoge-neous and isotropic background, such as in most radiation-dominated to matter-dominated recombination era modelsof SF.Perturbative dynamics in the hydrodynamic limit arewell recorded in the literature, including references Binneyand Tremaine (2008) and Mo et al. (2010). The observablesof fluid density and velocity are tracked in this demonstra-tion, though higher moments are straightforward to com-pute. The density takes the classical definition of densityperturbations in a co-moving frame ρ ( (cid:126)x ) = M (cid:90) d pf ( (cid:126)x, (cid:126)p, t ) = ¯ ρa ( t ) (1 + δ ( (cid:126)x, t )) (36)where M is the enclosed mass and ¯ ρ is the average densityover the Hubble volume. The fluid velocity takes on a similardefinition (cid:104) (cid:126)v ( (cid:126)x ) (cid:105) = 1 ρ (cid:90) d pf ( (cid:126)x, (cid:126)p, t ) (cid:126)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 ∂ t f (1) + (cid:126)pma · (cid:126) ∇ f (1) − (cid:126) ∇ Φ (cid:48) · (cid:126) ∇ p f (1) − (cid:126) ∇ ¯Φ · (cid:126) ∇ p f (1) = 0 (38)is expanded in orders of density and velocity perturbations.Expectations for the zeroth and first (cid:126)v moments of the dis-tribution are found to have dynamics0 th : ∂ t δ + 1 a (cid:126) ∇ · { (1 + δ ) (cid:104) (cid:126)v (cid:105)} = 0 (39)1 st : ∂ t { (1 + δ ) (cid:104) (cid:126)v (cid:105)} + H (1 + δ ) (cid:104) (cid:126)v (cid:105) + (1 + δ ) a (cid:126) ∇ (cid:0) Φ (cid:48) + Φ (cid:1) + 1 a (cid:126) ∇ · { (1 + δ ) (cid:104) (cid:126)v ⊗ (cid:126)v (cid:105)} = 0 (40)giving density and momentum conservation laws respec-tively. An Euler-like form may be derived through a com-bination of Eqns. 39 and 40: ∂ t (cid:104) (cid:126)v (cid:105) + H (cid:104) (cid:126)v (cid:105) + 1 a (cid:16) (cid:104) (cid:126)v (cid:105) · (cid:126) ∇ (cid:17) (cid:104) (cid:126)v (cid:105) + 1 a (cid:126) ∇ (cid:0) Φ (cid:48) + Φ (cid:1) + 1(1 + δ ) a (cid:126) ∇ · { (1 + δ ) ( (cid:104) (cid:126)v ⊗ (cid:126)v (cid:105) − (cid:104) (cid:126)v (cid:105) ⊗ (cid:104) (cid:126)v (cid:105) ) } = 0 , (41)where the last term disappears as cold dark matter has van-ishingly small velocity dispersion. The fluctuation equationof motion is derived by combining the Euler form and matterconservation again, and to first order looks like a dampedharmonic instability ∂ t δ + 2 H∂ t δ = 4 πG ¯ ρ (cid:0) δ + δ (cid:48) (cid:1) + O ( δ ) (42)where δ (cid:48) is the energy density fluctuation from complemen-tary species. This gives the expected collapse on all con-cerned scales, in the limit of no back reactions.Calculating perturbative Bose hydrodynamics followsvery similarly. The lowest order XC effects modify the per-ceived depth of the self-gravitational well, which is analo-gous to modifying the observable cosmological dark matterdensity, but only within the DM species. The fluctuationequation of motion is derived to be ∂ t δ + 2 H∂ t δ = 4 πG ¯ ρ (cid:0) Cδ + δ (cid:48) (cid:1) − πG ¯ ρ λ + λ W δ + O ( λ , δ ) , (43) [] P ( k ) | δC | = | δC | = | δC | = -3 -2 -1 k [Mpc − ]10 -5 -4 -3 -2 -1 F r a c . D i ff . Figure 2.
Perturbation power spectra evolved to z = 0. The up-per panel shows the perturbation-theory estimated power spectraevolved under both CDM and Bose dynamics. The CDM spec-tra is consistent with that found in Spergel et al. (2007). Bosespectra are calculated using the same relic density as CDM, plusthe gravity amplification derived in Eqns. 43 and 44. Deviationsin the Bose initial correlations δC vary over multiple orders inmagnitude. Spectra are shown well into the non-linear regime, asdemarcated by the vertical dashed line. The lower panel showsfractional power differences between the CDM and Bose spectra. where the λ s are evaluated at the initial spectrum and W represents the phase space volume over which the λ s arecalculated. The used expansion of the correlation functioncan be found in Appx. D. MFT is reached in the homo-geneous limit, meaning to zeroth order C = C = 1 and λ i = λ i, = 0. Strictly speaking, the XC has no impacton first order dynamics as the non-trivial components of C and the λ s enter at first order in perturbations. However,an inflated XC impact will be entertained here for the pur-poses of demonstration. A simple relation between the initialcorrelation and the constraints is found at first order in per-turbations C − C = λ , + λ , W . (44)Past and present probes of structure are capable of de-termining the early power spectra and concentration of grav-itating species with percent level precision (Spergel et al. et al. | δC | ∈ [0 . , . MNRAS , 1– ????
Perturbation power spectra evolved to z = 0. The up-per panel shows the perturbation-theory estimated power spectraevolved under both CDM and Bose dynamics. The CDM spec-tra is consistent with that found in Spergel et al. (2007). Bosespectra are calculated using the same relic density as CDM, plusthe gravity amplification derived in Eqns. 43 and 44. Deviationsin the Bose initial correlations δC vary over multiple orders inmagnitude. Spectra are shown well into the non-linear regime, asdemarcated by the vertical dashed line. The lower panel showsfractional power differences between the CDM and Bose spectra. where the λ s are evaluated at the initial spectrum and W represents the phase space volume over which the λ s arecalculated. The used expansion of the correlation functioncan be found in Appx. D. MFT is reached in the homo-geneous limit, meaning to zeroth order C = C = 1 and λ i = λ i, = 0. Strictly speaking, the XC has no impacton first order dynamics as the non-trivial components of C and the λ s enter at first order in perturbations. However,an inflated XC impact will be entertained here for the pur-poses of demonstration. A simple relation between the initialcorrelation and the constraints is found at first order in per-turbations C − C = λ , + λ , W . (44)Past and present probes of structure are capable of de-termining the early power spectra and concentration of grav-itating species with percent level precision (Spergel et al. et al. | δC | ∈ [0 . , . MNRAS , 1– ???? (2018) E. W. Lentz et al. [] T l ( l + ) [ µ K ] CDM Baseline | δC | = | δC | = | δC | = F r a c . D i ff . Figure 3.
Perturbative baryon angular power spectra at recom-bination as influenced by background DM. The top panel showsthe perturbative angular power spectra induced into the baryoniccomponent by both CDM and Bose background distributions.Deviations in the Bose initial correlations δC vary over multi-ple orders in magnitude. Calculations were performed in CAMB(Howlett et al. The model derived in Section 3 displays a number of interest-ing properties. The extremal nature of an exact condensatereveals the nature of inter-axion correlation and exchangein a most natural way, producing new physics on all con-sidered scales. The unique physics appears both in terms ofhigh-order dynamical moments and as novel gravity-inspiredpotentials, all stemming from the correlation function ˜ g . So-lutions to the conformal Lagrange multipliers in turn dependon global invariants of the distribution, perhaps due to theunbounded nature of the extremal correlations and exchangesymmetry. More realistic depictions of primordial axions arelikely to ruin this exceptional case by introduction of newlength scales and state transitions, though the current formmay be considered to represent the best-case scenario forseeing physics beyond the dissociated limit.The evolution of linear structures strictly reveals no de-parture of Bose DM from CDM infall. Deviations enter onlyat the second order in perturbation and reveal only milddivergences in dynamics and early structure when those lib-erties are taken; only marginal impacts to current isolatedsmall-scale structures should be expected. This lack of scaledependence is intriguing, and again held to be related to theunbounded nature of the correlation and the non-local ex-change symmetry within the region of interest. The marginalresult of new structure is not surprising as the non-trivialportion of the correlation function is highly non-linear.The final unresolved problem for the dynamics of con-densed QCD axions is specifying the value of initial corre-lation. A complete model tracking of quantum correlations from pre-inflation to the current era is needed to make defini-tive statements about the presence of axion physics beyondthe mean field. Such models, even for a subset of early cos-mological history, are currently lacking in the literature. The well-motivated axion DM candidate has properties thatallow it to simultaneously form a Bose-Einstein conden-sate and a level of inter-particle correlation unavailable tomost bosons studied in material science. This paper soughtto derive a robust and computationally effective model ofaxion infall from first principles, and has done so in thecondensed limit. The unique nature of such a maximallycondensed axion fluid displays distinct departures from thestandard collision-less fluid picture in the relevant case ofself-gravitation. These cross-correlation impacts do not ap-pear to significantly influence the outcome of linear struc-ture formation, but may impact structures entering the non-linear regime. The next entries to this series introduce nu-merical simulation techniques to Bose gravitational collapseand explore the nuances of imperfect and dynamical axioncondensation.
We would like to thank Jens Niemeyer, Katy Clough, DavidMarsh, Bodo Schwabe, Jan Veltmaat, and Xiaolong Du fortheir productive discussions in the refinement of this paper.We also gratefully acknowledge the support of the U.S. De-partment of Energy office of High Energy Physics and theNational Science Foundation. TQ was supported in part bythe NSF grant AST-1514868. EL and LR were supported inpart by the DOE grant DE-SC0011665.
REFERENCES
R. D. Peccei and H. R. Quinn, Physical Review Letters , 1440(1977).S. Weinberg, Phys. Rev. Lett. , 223 (1978).F. Wilczek, Phys. Rev. Lett. , 279 (1978).K. Mimasu and V. Sanz, Journal of High Energy Physics , 173(2015), arXiv:1409.4792 [hep-ph] .C. Patrignani and Particle Data Group, Chinese Physics C ,100001 (2016).D. J. E. Marsh, Phys. Rep. , 1 (2016), arXiv:1510.07633 .G. G. Raffelt, “Astrophysical axion bounds,” in Axions: Theory,Cosmology, and Experimental Searches , edited by M. Kuster,G. Raffelt, and B. Beltr´an (Springer Berlin Heidelberg,Berlin, Heidelberg, 2008) pp. 51–71.J. Isern, E. Garc´ıa-Berro, L. G. Althaus, and A. H. C´orsico, A&A , A86 (2010), arXiv:1001.5248 [astro-ph.SR] .A. H. C´orsico, 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. , 010 (2012),arXiv:1211.3389 [astro-ph.SR] .N. Viaux, M. Catelan, P. B. Stetson, G. G. Raffelt, J. Redondo,A. A. R. Valcarce, and A. Weiss, A&A , A12 (2013),arXiv:1308.4627 [astro-ph.SR] .L. F. Abbott and P. Sikivie, Physics Letters B , 133 (1983).M. Dine and W. Fischler, Physics Letters B , 137 (1983).MNRAS , 1– ?? (2018) xion Structure Formation I J. Preskill, M. B. Wise, and F. Wilczek, Physics Letters B ,127 (1983).J. E. Kim, Phys. Rev. Lett. , 103 (1979).M. A. Shifman, A. I. Vainshtein, and V. I. Zakharov, NuclearPhysics B , 493 (1980).M. Dine, W. Fischler, and M. Srednicki, Physics Letters B ,199 (1981).A. R. Zhitnitsky, Sov. J. Nucl. Phys. , 260 (1980), [Yad.Fiz.31,497(1980)].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 ADMXCollaboration, Physical Review Letters , 041301 (2010),arXiv:0910.5914 [astro-ph.CO] .ADMX Collaboration, In Preparation , (2017).M. Arik, S. Aune, K. Barth, A. Belov, S. Borghi, H. Br¨auninger,G. Cantatore, J. M. Carmona, S. A. Cetin, J. I. Col-lar, E. Da Riva, T. Dafni, M. Davenport, C. Eleftheri-adis, N. Elias, G. Fanourakis, E. Ferrer-Ribas, P. Friedrich,J. Gal´an, J. A. Garc´ıa, A. Gardikiotis, J. G. Garza, E. N.Gazis, T. Geralis, E. Georgiopoulou, I. Giomataris, S. Gni-nenko, H. G´omez, M. G´omez Marzoa, E. Gruber, T. Guth¨orl,R. Hartmann, S. Hauf, F. Haug, M. D. Hasinoff, D. H. H.Hoffmann, F. J. Iguaz, I. G. Irastorza, J. Jacoby, K. Jakovˇci´c,M. Karuza, K. K¨onigsmann, R. Kotthaus, M. Krˇcmar,M. Kuster, B. Laki´c, P. M. Lang, J. M. Laurent, A. Lio-lios, A. Ljubiˇci´c, G. Luz´on, 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´as, M. Tsagri, K. vanBibber, T. Vafeiadis, J. Villar, J. K. Vogel, S. C. Yildiz,and K. Zioutas (CAST Collaboration), Phys. Rev. Lett. ,091302 (2014).M. Valluri, S. R. Loebman, J. Bailin, A. Clarke, V. P. Debat-tista, and G. Stinson, in
The General Assembly of Galaxy Ha-los: Structure, Origin and Evolution , IAU Symposium, Vol.317, edited by A. Bragaglia, M. Arnaboldi, M. Rejkuba, andD. Romano (2016) pp. 358–359, arXiv:1510.06006 .J. Binney, ArXiv e-prints (2017), arXiv:1706.01374 .E. W. Kolb and I. I. Tkachev, Physical Review Letters , 3051(1993), hep-ph/9303313 .P. Tinyakov, I. Tkachev, and K. Zioutas, J. Cosmology Astropart.Phys. , 035 (2016), arXiv:1512.02884 .M. Fairbairn, D. J. E. Marsh, J. Quevillon, and S. Rozier, ArXive-prints (2017), arXiv:1707.03310 .N. Banik and P. Sikivie, ArXiv e-prints (2015), arXiv:1501.05913.J. Redondo, Journal of Cosmology and Astroparticle Physics , 008 (2013).M. Gorghetto, E. Hardy, and G. Villadoro, Journal of High En-ergy Physics , 151 (2018), arXiv:1806.04677 [hep-ph] .A. Vaquero, J. Redondo, and J. Stadler, ArXiv e-prints (2018),arXiv:1809.09241 .A. Arvanitaki, S. Dimopoulos, S. Dubovsky, N. Kaloper,and J. March-Russell, Phys. Rev. D , 123530 (2010),arXiv:0905.4720 [hep-th] .H.-Y. Schive, T. Chiueh, and T. Broadhurst, Nature Physics ,496 (2014), arXiv:1406.6586 .J. Veltmaat, J. C. Niemeyer, and B. Schwabe, ArXiv e-prints(2018), arXiv:1804.09647 .P. Sikivie, Phys. Rev. D , 063501 (1999a), astro-ph/9902210 .P. Sikivie, Nuclear Physics B Proceedings Supplements , 110(1999b), astro-ph/9810286 .L. D. Duffy and P. Sikivie, Phys. Rev. D , 063508 (2008),arXiv:0805.4556 .S. S. Chakrabarty and P. Sikivie, in APS April Meeting Abstracts (2017).N. Banik, A. J. Christopherson, P. Sikivie, and E. M. Todarello, Phys. Rev. D , 123540 (2015).O. Erken, P. Sikivie, H. Tam, and Q. Yang, Phys. Rev. D ,063520 (2012), arXiv:1111.1157 [astro-ph.CO] .J. Binney and S. Tremaine, Galactic Dynamics: Second Edi-tion, 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).C. W. Misner, K. S. Thorne, and J. A. Wheeler,
San Francisco:W.H. Freeman and Co., 1973 (1973).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. Cut-taia, P. de Bernardis, G. de Zotti, J. Delabrouille, J.-M.Delouis, E. Di Valentino, J. M. Diego, O. Dor´e, M. Dous-pis, A. Ducout, X. Dupac, S. Dusini, G. Efstathiou, F. El-sner, 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´enova-Santos, M. Gerbino, T. Ghosh,J. Gonz´alez-Nuevo, K. M. G´orski, 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¨anen, 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´opez-Caniego, P. M. Lubin, Y.-Z. Ma, J. F. Mac´ıas-P´erez, G. Maggio, D. Maino, N. Mandolesi, A. Mangilli,A. Marcos-Caballero, M. Maris, P. G. Martin, M. Martinelli,E. Mart´ınez-Gonz´alez, S. Matarrese, N. Mauri, J. D. McEwen,P. R. Meinhold, A. Melchiorri, A. Mennella, M. Migliaccio,M. Millea, S. Mitra, M.-A. Miville-Deschˆenes, 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. Po-lastri, G. Polenta, J.-L. Puget, J. P. Rachen, M. Reinecke,M. Remazeilles, A. Renzi, G. Rocha, C. Rosset, G. Roudier,J. A. Rubi˜no-Mart´ın, B. Ruiz-Granados, L. Salvati, M. San-dri, 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. Vib-ert, P. Vielva, F. Villa, N. Vittorio, B. D. Wandelt, I. K. We-hus, M. White, S. D. M. White, A. Zacchei, and A. Zonca,ArXiv e-prints (2018), arXiv:1807.06209 .E. W. Lentz, T. R. Quinn, and L. J. Rosenberg, ArXiv e-prints(2018), arXiv:1808.06378 .E. Runge and E. K. U. Gross, Phys. Rev. Lett. , 997 (1984a).E. Runge and E. K. U. Gross, Phys. Rev. Lett. , 997 (1984b).A. D. Fokker, Annalen der Physik , 810 (1914).A. Kolmogoroff, Mathematische Annalen , 415 (1931).H. Mo, F. C. van den Bosch, and S. White, Galaxy Formationand Evolution, by Houjun Mo , Frank van den Bosch , SimonWhite, Cambridge, UK: Cambridge University Press, 2010 (2010).D. N. Spergel, R. Bean, O. Dor´e, 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. Wei-land, E. Wollack, and E. L. Wright, ApJS , 377 (2007),astro-ph/0603449 .MNRAS , 1– ????
San Francisco:W.H. Freeman and Co., 1973 (1973).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. Cut-taia, P. de Bernardis, G. de Zotti, J. Delabrouille, J.-M.Delouis, E. Di Valentino, J. M. Diego, O. Dor´e, M. Dous-pis, A. Ducout, X. Dupac, S. Dusini, G. Efstathiou, F. El-sner, 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´enova-Santos, M. Gerbino, T. Ghosh,J. Gonz´alez-Nuevo, K. M. G´orski, 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¨anen, 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´opez-Caniego, P. M. Lubin, Y.-Z. Ma, J. F. Mac´ıas-P´erez, G. Maggio, D. Maino, N. Mandolesi, A. Mangilli,A. Marcos-Caballero, M. Maris, P. G. Martin, M. Martinelli,E. Mart´ınez-Gonz´alez, S. Matarrese, N. Mauri, J. D. McEwen,P. R. Meinhold, A. Melchiorri, A. Mennella, M. Migliaccio,M. Millea, S. Mitra, M.-A. Miville-Deschˆenes, 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. Po-lastri, G. Polenta, J.-L. Puget, J. P. Rachen, M. Reinecke,M. Remazeilles, A. Renzi, G. Rocha, C. Rosset, G. Roudier,J. A. Rubi˜no-Mart´ın, B. Ruiz-Granados, L. Salvati, M. San-dri, 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. Vib-ert, P. Vielva, F. Villa, N. Vittorio, B. D. Wandelt, I. K. We-hus, M. White, S. D. M. White, A. Zacchei, and A. Zonca,ArXiv e-prints (2018), arXiv:1807.06209 .E. W. Lentz, T. R. Quinn, and L. J. Rosenberg, ArXiv e-prints(2018), arXiv:1808.06378 .E. Runge and E. K. U. Gross, Phys. Rev. Lett. , 997 (1984a).E. Runge and E. K. U. Gross, Phys. Rev. Lett. , 997 (1984b).A. D. Fokker, Annalen der Physik , 810 (1914).A. Kolmogoroff, Mathematische Annalen , 415 (1931).H. Mo, F. C. van den Bosch, and S. White, Galaxy Formationand Evolution, by Houjun Mo , Frank van den Bosch , SimonWhite, Cambridge, UK: Cambridge University Press, 2010 (2010).D. N. Spergel, R. Bean, O. Dor´e, 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. Wei-land, E. Wollack, and E. L. Wright, ApJS , 377 (2007),astro-ph/0603449 .MNRAS , 1– ???? (2018) E. W. Lentz et al.
C. Howlett, A. Lewis, A. Hall, and A. Challinor, J. CosmologyAstropart. Phys. , 027 (2012), arXiv:1201.3654 [astro-ph.CO] .W. H. Zurek, eprint arXiv:quant-ph/0306072 (2003), quant-ph/0306072 .C. K. Zachos, D. B. Fairlie, and T. L. Curtright,
Quantum me-chanics in phase space , by Zachos, Cosmas. Quantum me-chanics in phase space, editors, Cosmas K. Zachos...[et al.].,Singapore; New Jersey, N.J.: World Scientific, c2005. (2005).D. I. Bondar, R. Cabrera, D. V. Zhdanov, and H. A. Rabitz,Phys. Rev. A , 263 (2013), arXiv:1202.3628 [quant-ph] . APPENDIX A: THE MANY-BODYSCHR ¨ODINGER EQUATION AND ONESPECIAL CLASS OF SYSTEMS
Let us consider the problem of axion infall as a solutionof the many-body Schr¨odinger equation under inter-particleNewtonian gravitation. i (cid:126) ∂ t Ψ ( (cid:126)x , ..., (cid:126)x N ; t ) = ˆ H Ψ ( (cid:126)x , ..., (cid:126)x N ; t ) (A1)with the Hamiltonian of the formˆ H = ˆ T + ˆ V int = − N (cid:88) i (cid:126) ∇ a m a + N (cid:88) i 1) ‘parti-cles’ with Hamiltonian H int defined over the function space L ( (cid:60) N ( N − ) ⊗ C ( (cid:60) ) with an analogous Bose constraint L ( (cid:60) N ( N − ) ⊗ C ( (cid:60) ) /S N ( N − (A14)and all the implications thereof. In this case the Hamiltonianis separable and the stationary states are spanned by singlepermanents of single-‘particle’ stationary statesΦ (cid:48) (∆ ij ) ∈ sp (cid:40) perm (cid:32)(cid:89) ij φ α ij ( (cid:126)x i − (cid:126)x j , t ) (cid:33)(cid:41) α ∈ A (A15)and the pure states asΨ (cid:48) (∆ ij , t ) ∈ sp (cid:40) perm (cid:32)(cid:89) ij ψ α ij ( (cid:126)x i − (cid:126)x j , t ) (cid:33)(cid:41) α ∈ A (A16)where sp () gives the linear span of a collection of functions, perm () is the permanent operation, A is the set of many-body solution indexes, ˆ ψ a is a solution to the single-bodyequation i (cid:126) ∂ t ψ a = − (cid:126) a µ ∇ ( (cid:126)x i − (cid:126)x j ) ψ α + φ ( (cid:126)x i − (cid:126)x j ) ψ a (A17)Solutions of this N ( N − 1) space may be used in the originaluse-case of N axions by applying an appropriate projection p : ∆ ij = (cid:126)x i − (cid:126)x j ∀ i, j (A18)on the Hamiltonian and function space pair. (cid:16) ˆ H (cid:16) (cid:60) N ( N − (cid:17) , C (cid:16) (cid:60) N ( N − (cid:17)(cid:17) → (cid:16) ˆ H (cid:16) (cid:60) N +1 (cid:17) , C (cid:16) (cid:60) N +1 (cid:17)(cid:17) (A19)It is apparent that this projection both respects the modu-lated many-body solution space and the N-body Schr¨odingerequation. It can also be shown that the span of solu-tions to the N -body problem can be represented by pro-jected solutions of the N ( N − p (cid:16) L ( (cid:60) N ( N − ) ⊗ C ( (cid:60) ) /S N ( N − (cid:17) is dense in L ( (cid:60) N ) ⊗ C ( (cid:60) ) /S N . APPENDIX B: THE MANY-BODY WIGNERFUNCTION AND ITS SUPER-DE BROGLIELIMIT 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. B1 Wigner Transform Phase space DFs are often far better suited for tracking sys-tem observables than wave-functions, due to the hyperbol-icity of their governing dynamics. The many-body Bose sys-tem can be reduced into a phase-space pseudo-DF via themany-body Wigner function f ( N ) ( (cid:126)x , (cid:126)p , ..., (cid:126)x N , (cid:126)p N , t )= (cid:90) d x (cid:48) · · · d x (cid:48) N e i (cid:80) Nj (cid:126)p j · (cid:126)x (cid:48) j / (cid:126) × Ψ † ( (cid:126)x + (cid:126)x (cid:48) / , ..., (cid:126)x N + (cid:126)x (cid:48) N / , t ) × Ψ( (cid:126)x − (cid:126)x (cid:48) / , ..., (cid:126)x N − (cid:126)x (cid:48) N / , t ) (B1)For brevity, the arguments of Ψ † and Ψ are understood tobe { (cid:126)x i + (cid:126)x (cid:48) i / } and { (cid:126)x i − (cid:126)x (cid:48) i / } respectively. While not atrue DF, in that f ( N ) can take on negative values, it doesshare a number of properties with true DFs, such as beingreal-valued f ( N ) ( (cid:126)x, (cid:126)p ) ∗ = (cid:90) d x (cid:48) e − i (cid:80) Nj (cid:126)p j · (cid:126)x (cid:48) j / (cid:126) (cid:16) Ψ † Ψ (cid:17) ∗ = (cid:90) d x (cid:48) e − i (cid:80) Nj (cid:126)p j · (cid:126)x (cid:48) j / (cid:126) ΨΨ † = (cid:90) d x (cid:48) e i (cid:80) Nj (cid:126)p j · (cid:126)x (cid:48) j / (cid:126) Ψ † Ψ, ( (cid:126)x (cid:48) → − (cid:126)x (cid:48) )= f ( N ) ( (cid:126)x, (cid:126)p ) , (B2) MNRAS , 1– ?? (2018) E. W. Lentz et al. and providing properly-normalized probability densities overprojections to 3D single-body spaces ρ ( (cid:126)x ) = (cid:90) d p d N − wf ( N ) , (B3) ρ ( (cid:126)p ) = (cid:90) d x d N − wf ( N ) . (B4)The Wigner pseudo-DF can also be expected to provide a re-liable DF in certain circumstances (Zurek 2003; Zachos et al. et al. B2 Super de-Broglie Equation of Motion The equation of motion for the distribution is found bystraightforward differentiation ∂ t f ( N ) ( (cid:126)x , (cid:126)p , ..., (cid:126)x N , (cid:126)p N , t )= (cid:90) d x (cid:48) · · · d x (cid:48) N e i (cid:80) Nj (cid:126)p j · (cid:126)x (cid:48) j / (cid:126) (cid:16) ∂ t Ψ † Ψ + Ψ † ∂ t Ψ (cid:17) (B5)Judicious substitution of the Schr¨odinger equation trans-forms the expression into an equation of motion ∂ t f ( N ) ( (cid:126)x , (cid:126)p , ..., (cid:126)x N , (cid:126)p N , t )= 1 i (cid:126) (cid:90) d x (cid:48) · · · d x (cid:48) N e i (cid:80) Nj (cid:126)p j · (cid:126)x (cid:48) j / (cid:126) (cid:16) − ˆ H Ψ † Ψ + Ψ † ˆ H Ψ (cid:17) (B6)where we use the Hamiltonian of Eqn.28.The kinetic contribution can be simplified through suc-cessive application of the product rule of differentiation andconsolidation of terms, and as expected provides transportterms K = 1 i (cid:126) (cid:126) a m N (cid:88) i (cid:90) d x (cid:48) · · · d x (cid:48) N e i (cid:80) Nj (cid:126)p j · (cid:126)x (cid:48) j / (cid:126) × (cid:16) ∇ (cid:126)x i + (cid:126)x (cid:48) i / Ψ † Ψ − Ψ † ∇ (cid:126)x i − (cid:126)x (cid:48) i / Ψ (cid:17) = − N (cid:88) i (cid:126)p i a m (cid:126) ∇ i f ( N ) + O ( (cid:126) ) (B7)where de Broglie wavelength contributions and higher or-der are inconsequential. The quantum diffusion scales of theQCD axion are many orders of magnitude smaller than thedesired scale of simulation, causing them to average out tozero over cosmological coarse graining.The external potentials give forcing contributions, cal-culated through Taylor expansion of the potentials about (cid:126)x i = 0: F ext = 1 i (cid:126) N (cid:88) i (cid:90) d x (cid:48) · · · d x (cid:48) N e i (cid:80) Nj (cid:126)p j · (cid:126)x (cid:48) j / (cid:126) × (cid:16) Φ (cid:48) ( (cid:126)x i + (cid:126)x (cid:48) i / † Ψ − Ψ † Φ (cid:48) ( (cid:126)x i − (cid:126)x (cid:48) i / (cid:17) = N (cid:88) i (cid:126) ∇ i Φ (cid:48) · (cid:126) ∇ p i f ( N ) + O ( (cid:126) ) (B8)Lastly, the two-body interaction potentials also contribute aforcing term F int = 12 i (cid:126) N (cid:88) i (cid:54) = j (cid:90) d x (cid:48) · · · d x (cid:48) N e i (cid:80) Nj (cid:126)p j · (cid:126)x (cid:48) j / (cid:126) × (cid:16) φ ij ( (cid:126)x i + (cid:126)x (cid:48) i / , (cid:126)x j + (cid:126)x (cid:48) j / † Ψ − Ψ † φ ij ( (cid:126)x i − (cid:126)x (cid:48) i / , (cid:126)x j − (cid:126)x (cid:48) j / (cid:17) = N (cid:88) i (cid:54) = j (cid:126) ∇ i φ ij · (cid:126) ∇ p i f ( N ) + O ( (cid:126) )recalling that φ ij is the two-body gravitational potential.The result is a Liouville-like equation of motion that appearsvery similar to the classical form for a collision-less fluid ∂ t f ( N ) = K + F ext + F int (B9)= − N (cid:88) i (cid:126)p i a m (cid:126) ∇ i f ( N ) + N (cid:88) i (cid:126) ∇ i Φ (cid:48) · (cid:126) ∇ p i f ( N ) + 12 N (cid:88) i (cid:54) = j (cid:126) ∇ i φ ij · (cid:126) ∇ p i f ( N ) + O ( (cid:126) ) (B10)We concern ourselves with the Liouville-like equation com-posed of only lowest order terms. Note that so far we havenot enforced specific dynamics or even exchange conditionson the state. Such a Boltzmann equation may also by appliedto analogous classes of fermionic or mixed systems. APPENDIX C: PROJECTION OF THEMANY-BODY DISTRIBUTION FUNCTION Projection of the many-body DF down to the total phasespace density is accomplished through integration of phasespaces 2 , , ..., N . For general quantum and classical systems,correlations in the (psuedo-)DF f ( N ) may be present be-tween any subset of particles. Fortunately, the topic systemprovides us with insight to the correlations in the form ofa complete basis of solutions to the governing Schr¨odingerequation. C1 Reduction to Single DF We follow the well-used approach of density projection byintegrating out the supurfluous degrees of freedom f ( M ) = (cid:90) d w M +1 · · · d w N f ( N ) (C1)where w i = ( (cid:126)x i , (cid:126)p i ) is the i-th six dimensional phase space.Projection of the governing super de Broglie Liouville equa- MNRAS , 1– ?? (2018) xion Structure Formation I tion ∂ t f ( N ) + N (cid:88) i (cid:126)p i a m a (cid:126) ∇ i f ( N ) − N (cid:88) i (cid:126) ∇ i Φ (cid:48) · (cid:126) ∇ p i f ( N ) − N (cid:88) i (cid:54) = j (cid:126) ∇ i φ ij · (cid:126) ∇ p i f ( N ) = O ( (cid:126) ) (C2)is expected to follow analogously.Performing the single-body projection on the lowest or-der of the Liouville equation0 = (cid:82) d w · · · d w N × (cid:104) ∂ t f ( N ) + (cid:80) Ni (cid:126)p i m (cid:126) ∇ i f ( N ) − (cid:80) Ni (cid:126) ∇ i Φ (cid:48) · (cid:126) ∇ p i f N − (cid:80) Ni (cid:54) = j (cid:126) ∇ i φ ij · (cid:126) ∇ p i f ( N ) (cid:105) (C3)reduces the explicit time component to T (1) ≡ (cid:90) d w · · · d w N ∂ t f ( N ) = ∂ t f (1) , (C4)The kinetic transport term reduces as expected, due to thedivergence theorem and axion conservation K (1) ≡ (cid:90) d w · · · d w N N (cid:88) i (cid:126)p i a m (cid:126) ∇ i f ( N ) = (cid:126)pa m a (cid:126) ∇ f (1) , (C5)The external potential forcing term reduces similarly F (1) ext ≡ − (cid:90) d w ··· d w N N (cid:88) i (cid:126) ∇ i Φ (cid:48) · (cid:126) ∇ p i f ( N ) = − (cid:126) ∇ Φ (cid:48) · (cid:126) ∇ p f (1) (C6)It is the self-gravity forcing term that is the most trouble-some, as the two-body gravitational interactions contain ameasure of two-body correlations F (1) int ≡ − (cid:90) d w · · · d w N N (cid:88) i (cid:54) = j (cid:126) ∇ i φ ij · (cid:126) ∇ p i f ( N ) = − (cid:90) d w j N (cid:88) j> (cid:126) ∇ φ j · (cid:126) ∇ p f (2) (C7) C2 Correlation Hints Let us take stock of correlations in the many-body distribu-tion function before attempting further simplification. Re-call that the considered many-body DF may be expressed interms of a single condensed state f ( N ) ( (cid:126)x , (cid:126)p , ..., (cid:126)x N , (cid:126)p N , t )= (cid:90) d x (cid:48) · · · d x (cid:48) N e i (cid:80) Nj (cid:126)p j · (cid:126)x (cid:48) j / (cid:126) (cid:89) i (cid:54) = j ˆ ψ † a (cid:0) (cid:126)x i − (cid:126)x j + (cid:126)x (cid:48) i / − (cid:126)x (cid:48) j / , t (cid:1) × (cid:89) i (cid:54) = j ˆ ψ a (cid:0) (cid:126)x i − (cid:126)x j − (cid:126)x (cid:48) i / (cid:126)x (cid:48) j / , t (cid:1) (C8)where the condensed state is of the form derived in Appx A.Condensate wave-functions comprised solely of two-axioncorrelators conspire with the reach of inter-axion gravity toexpose that two-body correlations exist between axions and that they are the only correlations relevant to the forma-tion of large scale structure. Finally, the long-range natureof exchange symmetry in this NR co-moving context indi-cates that the correlation function between axions should beintrinsically scale independent. C3 Correlation Function To derive the form of two-body correlations relevant for largescale structure, we first rewrite the two-body DF, f (2) , fromEqn. C7 in a form that makes the distinction between cor-relations and total density more obvious f (2) ( w , w , t ) = ˜ g × f (1) ( w , t ) f (1) ( w , t ) (C9)where the correlation function ˜ g is seen to be symmetric overboth phase spaces.This correlation form has some convenient properties,including converging with standard definitions of correlationin the condensed matter literature in the limit of N → ∞ .In addition to the standard normalization of the single-bodyDF, (cid:90) d w f (1) = 1 , (C10)there is also a normalization condition on the correlationfunction, (cid:90) d w ˜ gf (1) = 1 . (C11)This correlation function is not expected to be of greatcomputational utility as an explicitly wave-function expres-sion, though it can be shown that the totally condensedaxion system inhabits an extremal correlation state in thesense that the functional I = (cid:90) d w d w ˜ g (C12)has vanishing physical variation about the condensed state( δI = 0). A modified functional may be constructed I (cid:48) = I + λ (cid:18)(cid:90) d w f (1) − (cid:19) + λ (cid:18)(cid:90) d w ˜ gf (1) − (cid:19) (C13)using correlation-normalization Lagrange multipliers, con-straining the correlator to physicality. Beyond this point wecontinue the derivation using physical conjectures due to theintractability of evaluating the full many-body expressions.We conjecture that ˜ g ’s dependence on f (1) goes as f + = 12 (cid:16) f (1) ( w , t ) + f (1) ( w , t ) (cid:17) . (C14)The symmetric nature of ˜ g ( f + ) implies that the extremiza-tion of Eqn. C13 is equivalent to extremization of I (cid:48)(cid:48) = (cid:90) d w ˜ g + λ (cid:18)(cid:90) d w f (1) − (cid:19) + λ (cid:18)(cid:90) d w ˜ gf (1) − (cid:19) (C15)This functional can be simply extremized over f (1) to findthe condensed ˜ g . Variational principles show that the ex-tremals of I (cid:48)(cid:48) satisfy ∂ ˜ g∂f (1) + λ + λ (cid:18) ˜ g + f (1) ∂ ˜ g∂f (1) (cid:19) = 0 . (C16) MNRAS , 1– ???? To derive the form of two-body correlations relevant for largescale structure, we first rewrite the two-body DF, f (2) , fromEqn. C7 in a form that makes the distinction between cor-relations and total density more obvious f (2) ( w , w , t ) = ˜ g × f (1) ( w , t ) f (1) ( w , t ) (C9)where the correlation function ˜ g is seen to be symmetric overboth phase spaces.This correlation form has some convenient properties,including converging with standard definitions of correlationin the condensed matter literature in the limit of N → ∞ .In addition to the standard normalization of the single-bodyDF, (cid:90) d w f (1) = 1 , (C10)there is also a normalization condition on the correlationfunction, (cid:90) d w ˜ gf (1) = 1 . (C11)This correlation function is not expected to be of greatcomputational utility as an explicitly wave-function expres-sion, though it can be shown that the totally condensedaxion system inhabits an extremal correlation state in thesense that the functional I = (cid:90) d w d w ˜ g (C12)has vanishing physical variation about the condensed state( δI = 0). A modified functional may be constructed I (cid:48) = I + λ (cid:18)(cid:90) d w f (1) − (cid:19) + λ (cid:18)(cid:90) d w ˜ gf (1) − (cid:19) (C13)using correlation-normalization Lagrange multipliers, con-straining the correlator to physicality. Beyond this point wecontinue the derivation using physical conjectures due to theintractability of evaluating the full many-body expressions.We conjecture that ˜ g ’s dependence on f (1) goes as f + = 12 (cid:16) f (1) ( w , t ) + f (1) ( w , t ) (cid:17) . (C14)The symmetric nature of ˜ g ( f + ) implies that the extremiza-tion of Eqn. C13 is equivalent to extremization of I (cid:48)(cid:48) = (cid:90) d w ˜ g + λ (cid:18)(cid:90) d w f (1) − (cid:19) + λ (cid:18)(cid:90) d w ˜ gf (1) − (cid:19) (C15)This functional can be simply extremized over f (1) to findthe condensed ˜ g . Variational principles show that the ex-tremals of I (cid:48)(cid:48) satisfy ∂ ˜ g∂f (1) + λ + λ (cid:18) ˜ g + f (1) ∂ ˜ g∂f (1) (cid:19) = 0 . (C16) MNRAS , 1– ???? (2018) E. W. Lentz et al. [] f + /f + ,max ˜ g Figure C1. Correlation functions for several types of systems,both bosonic and fermionic. The lines that lie below unity at f + = f + ,max (red to green) are from non-interacting Fermion dis-tributions, organized in a spin-symmetric configuration, param-eterized by temperature decreasing with initial correlation like T = T f / ( C − , where T f is the Fermi temperature. The linesat and above unity (teal to violet) are from condensed scalar Bo-son distributions shaped into a spherical Gaussian in phase space,and also become colder with more extreme initial correlation viathe dispersion length relation ˜ σ p / ˜ σ x ∝ C . The symmetry of ˜ g over phase spaces also allows for sym-metrization of the extremal characteristic equation, resultingin the equivalent extremization condition ∂ ˜ g∂f + + λ + λ (cid:18) ˜ g + f + ∂ ˜ g∂f + (cid:19) = 0 . (C17)Solutions to the condition are found to be˜ g = C − λ f + λ f + (C18)where C is constrained by the correlation present in the ini-tial configuration. For instance, a separable Bose condensategives C = 1, with Lagrange multipliers also conforming tothe uncorrelated dynamics of standard MFTs.The extremal values of λ and λ may be found throughthe constraint equations. The correlation constraint, Eqn.C11, for the extremal solutions becomes1 = (cid:90) d w f (1) ( w , t ) C − λ f + λ f + . (C19)See Fig. C1 for an example of the behavior of this correlator.The forcing term of Eqn.C7 may now be written F int = − (cid:90) d w j N (cid:88) j> (cid:126) ∇ φ j · (cid:126) ∇ p f (2) == − (cid:90) d w j N (cid:88) j> (cid:126) ∇ φ j · (cid:126) ∇ p f (1) ( w , t ) f (1) ( w j , t ) × C − λ f + λ f + (C20) producing an new collision-less Boltzmann-like equation0 = ∂ t f (1) + (cid:126)pm · (cid:126) ∇ f (1) − (cid:126) ∇ Φ (cid:48) · (cid:126) ∇ p f (1) − NN − (cid:90) d w j (cid:126) ∇ φ j · (cid:126) ∇ p f (1) ( w , t ) f (1) ( w j , t ) × C − λ f + λ f + (C21)with the last term containing both the Newtonian and theexchange-correlation contributions of inter-axion gravita-tion. APPENDIX D: EXPANSIONS AND LIMITS OFTHE CONDENSED BOLTZMANN EQUATION The full form of the Boltzmann-like equation for condensedBose systems found in Appx. C is quite compact, thoughdifficult to use analytically due to the presence of non-polynomial non-linearities. Semi-analytic techniques, such asthe early universe structure formation calculations of Section4, would benefit from an expansion in orders. We provide thefirst few orders here.Expanding F int from Eqn. C20 in orders of λ f + andorganizing in terms of powers of f (1) produces F int = − (cid:90) d w j N (cid:88) j> (cid:126) ∇ φ j · (cid:126) ∇ p f (1) ( w , t ) f (1) ( w j , t ) × ( C − λ f + ) ∞ (cid:88) n =0 ( − λ f + ) n = − ( N − C (cid:126) ∇ ¯Φ · (cid:126) ∇ p f (1) + ( N − λ + λ (cid:32) (cid:126) ∇ ¯Φ · (cid:126) ∇ p (cid:16) f (1) (cid:17) + (cid:126) ∇ (cid:90) d w φ (cid:16) f (1) (cid:17) · (cid:126) ∇ p f (1) (cid:33) − ( N − λ ( λ + λ )4 (cid:32) (cid:126) ∇ ¯Φ · (cid:126) ∇ p (cid:16) f (1) (cid:17) + 2 (cid:126) ∇ (cid:90) d w φ (cid:16) f (1) (cid:17) · (cid:126) ∇ p (cid:16) f (1) (cid:17) + (cid:126) ∇ (cid:90) d w φ (cid:16) f (1) (cid:17) · (cid:126) ∇ p f (1) (cid:33) + ( N − λ ( λ + λ )8 (cid:32) (cid:126) ∇ ¯Φ · (cid:126) ∇ p (cid:16) f (1) (cid:17) + 3 (cid:126) ∇ (cid:90) d w φ (cid:16) f (1) (cid:17) · (cid:126) ∇ p (cid:16) f (1) (cid:17) + 3 (cid:126) ∇ (cid:90) d w φ (cid:16) f (1) (cid:17) · (cid:126) ∇ p (cid:16) f (1) (cid:17) + (cid:126) ∇ (cid:90) d w φ (cid:16) f (1) (cid:17) · (cid:126) ∇ p f (1) (cid:33) + O (cid:18)(cid:16) f (1) (cid:17) (cid:19) MNRAS , 1– ?? (2018) xion Structure Formation I The governing Boltzmann-like equation then expands to0 = ∂ t f (1) + (cid:126)pm · (cid:126) ∇ f (1) − (cid:126) ∇ Φ (cid:48) · (cid:126) ∇ p f (1) − (cid:126) ∇ ¯Φ · (cid:126) ∇ p f (1) + λ + λ (cid:32) (cid:126) ∇ ¯Φ · (cid:126) ∇ p (cid:16) f (1) (cid:17) + (cid:126) ∇ (cid:90) d w φ (cid:16) f (1) (cid:17) · (cid:126) ∇ p f (1) (cid:33) − λ ( λ + λ )4 (cid:32) (cid:126) ∇ ¯Φ · (cid:126) ∇ p (cid:16) f (1) (cid:17) + 2 (cid:126) ∇ (cid:90) d w φ (cid:16) f (1) (cid:17) · (cid:126) ∇ p (cid:16) f (1) (cid:17) + (cid:126) ∇ (cid:90) d w φ (cid:16) f (1) (cid:17) · (cid:126) ∇ p f (1) (cid:33) + λ ( λ + λ )8 (cid:32) (cid:126) ∇ ¯Φ · (cid:126) ∇ p (cid:16) f (1) (cid:17) + 3 (cid:126) ∇ (cid:90) d w φ (cid:16) f (1) (cid:17) · (cid:126) ∇ p (cid:16) f (1) (cid:17) + 3 (cid:126) ∇ (cid:90) d w φ (cid:16) f (1) (cid:17) · (cid:126) ∇ p (cid:16) f (1) (cid:17) + (cid:126) ∇ (cid:90) d w φ (cid:16) f (1) (cid:17) · (cid:126) ∇ p f (1) (cid:33) + O (cid:18)(cid:16) f (1) (cid:17) (cid:19) (D1)where we have taken N − N in cosmological contexts. Energy, momentum, and angularmomentum remain conserved at each level of expansion andin the full correlation form.Solutions to the Lagrange multipliers may also be ap-proximated in some cases, using a geometric expansion ofthe characteristic constraint equation1 = (cid:90) d w f (1) ( w , t ) C − λ f + λ f + assuming values of | λ f + | not exceeding unity over the re-gion of integration. The functional condition is then trans-lated into a sequence of quasi-algebraic conditions δ k = (cid:90) d w f (1) (cid:16) C − λ f (1) / (cid:17) ∞ (cid:88) n ≥ k (cid:32) nk (cid:33) (cid:18) − λ (cid:19) n (cid:16) f (1) (cid:17) n − k − (cid:90) d w f (1) λ f (1) / ∞ (cid:88) n ≥ k − (cid:32) nk − (cid:33) (cid:18) − λ (cid:19) n (cid:16) f (1) (cid:17) n − ( k − (D2)for k = 0 , , , ... . Only k = 0 and one other equation arelinearly independent. Here k = 0 , | λ f | ansatz. Order ofthe self-expectation values of f (1) ¯ f l ≡ (cid:90) d w (cid:16) f (1) (cid:17) l +1 (D3)are chosen as the expansion parameter, with the normal-ization constraint setting ¯ f = 1. The first few powers are analytically solvable, though in general a numerical schemeis needed to solve to arbitrary order. MNRAS , 1– ????