Backreaction of axion coherent oscillations

# Backreaction of axion coherent oscillations

Mor Rozner Physics department and Asher Space Science Institute, Technion, Haifa 3200003, Israel    Vincent Desjacques Physics department and Asher Space Science Institute, Technion, Haifa 3200003, Israel
July 26, 2019
###### Abstract

We investigate how coherent oscillations backreact on the evolution of the condensate wave function of ultra-light axions in the non-relativistic regime appropriate to cosmic structure formation. The coherent oscillations induce higher harmonics beyond the fundamental mode considered so far when a self-interaction is present, and imprint oscillations in the gravitational potential. We emphasize that the effective self-interaction felt by the slowly-varying envelop of the wave function always differs from the bare Lagrangian interaction potential. We also point out that, in the hydrodynamical formulation of the Gross-Pitaevskii equation, oscillations in the gravitational potential result in an attractive force that counteracts the effect of the quantum pressure arising from the strong delocalization of the particles. Since these effects become significant on physical scales less than the (large) Compton wavelength of the particle, they are presumably not very relevant on the mildly nonlinear scales traced by intergalactic neutral hydrogen for axion masses consistent with the bounds from the Lyman- forest. However, they might affect the formation of virialized cosmological structures and their stability.

###### pacs:
98.80.-k, 98.65.-r, 95.35.+d, 98.80.Es

## I Introduction

Since the virial equilibrium considerations of Zwicky and his inference of a “missing mass” in galaxy clusters Zwicky:1933 (), dark matter has become a inescapable ingredient of nearly all viable cosmological models PlanckCollaboration;2016 (). Nonetheless, its nature has remained thus far elusive, and there is a plethora of models trying to explain it. Beyond the popular weakly interacting massive particles (WIMPs) Steigman/Turner;1985 (), primordial black holes Chapline/Frampton:2016 () or a modification to gravity such as MOND Milgrom;1983 () could also explain at least part of the observations.

In another class of models, dark matter is made of light bosons such as axions Abbott:1982af (); Preskill:1982cy (); Dine:1981rt (); Baldeschi:1983mq (); Sin:1992bg (); HBG:2000 (); Riotto:2000kh (); Amendola/Barbieri:2006 (); Sikivie/Yank:2009 (); Kim/Marsh:2016 (). Such particles could be produced in the early Universe from a symmetry-breaking event conjectured to solve the strong CP problem Weinberg:1977ma (); Wilczek:1977pj (); PQ1 (); PQ2 (). From a cosmological point of view, ultra-light axions with a mass much smaller than the mass of QCD axions are particularly interesting because they could solve some of the small-scale problems associated with standard cold dark matter (CDM) like cusps vs. cores inside dark matter halos (see Hui/Ostriker/etal:2017 (); marsh:2016 () and references therein). However, observations of the Lyman- forest appear to rule out ultra-light axions as dark matter if their mass is less than Irsic:2017yje (); Armengaud:2017nkf (), although uncertainties in the thermal history of the intergalactic medium somewhat weaken this constraint Zhang:2017chj ().

“Fuzzy” dark matter (FDM) such as ultra-light axions is in the form of a Bose-Einstein condensate due to the very large de Broglie wavelength of the particles (the details of the condensation process are still uncertain, see for instance Kolb:1993hw (); Semikoz:1995rd (); Sikivie/Yang:2009 (); Davidson:2013aba (); Guth:2014hsa (); Levkov:2018kau ()). In the non-relativistic limit, the (scalar) wave function of the mean field evolves according to the Gross-Pitaevskii (GP) equation. On applying the Madelung transformation Madelung:1927 (); Bohm:1952 (), this equation can be re-expressed as standard hydrodynamical equations, in which gravity and axion self-interactions are balanced by the “quantum” pressure. This hydrodynamical representation can be easily implemented numerically to study the formation of cosmic structures in FDM models woo/chiueh:2008 (); Schive/etal:2014 (); veltmaat/niemeyer:2016 (); Zhang:2016uiy (); Mocz:2017wlg (); schive/chiueh:2018 (), for which existing Boltzmann code can be modified to provide the suitable initial conditions (e.g., Hlozek/Grin/etal:2015, ).

While the GP equation encodes the leading contributions to the hydrodynamic gradient expansion, higher-order effects, which are generally neglected, might play an important role at smaller scales (i.e., in the formation of “Bose stars” Ruffini:1969qy (); Tkachev:1986tr ().) and possibly affect the cosmic structure formation in FDM models. In particular, axion coherent oscillations will generate higher harmonics beyond the fundamental mode considered so far when a self-interaction is present (see e.g. visinelli/etal:2018, ; Eby:2017teq, , for a recent study in the context of “dense” Bose stars), and imprint oscillations in the gravitational potential Khmelnitsky/Rubakov:2013 (). The goal of this paper is to investigate how the coherent oscillations backreact on the evolution of the condensate wave function – in the non-relativistic regime appropriate to the large scale structure of the Universe – and how this backreaction manifests itself in the hydrodynamical formulation of the GP equation. Obviously, since the GP equation captures the leading contributions at large scales, this backreaction can only be significant on scales comparable to, or smaller than, the (huge) Compton wavelength of the particle.

Our paper is organized as follows. After a review of the ultra-light axion model of cold dark matter (Sec. §II), we pursue with a discussion of the backreaction of axion coherent oscillations at the homogeneous level (Sec. §III), before taking into account some of the backreaction effects in the derivation of the Gross-Pitaevskii system (Sec. §IV). We discuss the implications of our findings in Sec. §V and conclude in Sec. §VI.

We shall use the natural units throughout, with a gravitational coupling constant given in terms of the Planck mass . In these units, the Hubble constant is . Finally, Greek indices , etc. run over the 4 spacetime dimensions, whereas Latin indices , etc. run over the spatial dimensions solely.

## Ii Axion as a cold dark matter candidate

We begin with a brief overview of the theoretical results relevant to the description of axions in the large scale structure of the Universe (see marsh:2016, , for a review).

### ii.1 Lagrangian and characteristic scales

The starting point is the action for the axion (real scalar) field ,

 S[ϕ]=∫d4x√−g[12(∂ϕ)2−Λ4(1−cosϕf)], (1)

where is a sort of condensation scale, is the decay constant and is the metric determinant.

Expanding for and including the quartic coupling, one obtains a potential of the form

 V(ϕ)=12m2aϕ2−14!λϕ4, (2) wherem2a=Λ4f2% andλ=m2af2.

Our sign convention is such that when the quartic interaction is attractive. While our analysis is valid for any light boson in the form of a condensate, we shall adopt a fiducial axion mass and decay constant in all illustrations. This gives . Higher order terms ( and higher) are negligible as long as and this remains true as well when taking into account the high phase-space density. Notice that has dimension of energy, and that the sign of the self-interaction coupling leads to an attractive force. This will be relevant for all our considerations.

The relative importance of time-derivatives and gradients of the field can be estimated from a few characteristics ratios. Firstly, , where is the Hubble rate, represents the importance of Hubble friction relative to the coherent oscillations of the condensate. In matter domination, we have

 H0ma≃2.1×10−11ha−3/2(ma10−22 eV)−1. (3)

We have long after at the onset of oscillations, which occurs at a scale factor for the fiducial parameters. Secondly, , where is the comoving wavenumber and is the scale factor normalized to unity at the present epoch, characterizes the importance of kinetic terms. We have

 kmaa≃6.4×10−5a−1(kKpc−1)(ma10−22 eV)−1. (4)

This ratio involves the (reduced) Compton wavelength of the particle rather than its de Broglie wavelength , where is a typical velocity. After matter-radiation equality, never reaches unity on scales unless the particle mass is much less than . Finally, the dimensionless coupling constant quantifies the importance of the self-interaction. Note that always appears multiplied by a factor of , where is the axion energy density.

Since our focus is on the late-time, sub-horizon evolution of the axion field relevant to the formation of cosmic structures, we will consider a regime in which and . Furthermore, since we are interested in scales at which the density does not considerably exceed its background value, we also have . As a result, powers of , and are small and become rapidly negligible in a perturbative expansion (see, e.g., ratra:1991, ; Hwang/Noh:2009, ).

### ii.2 The Gross-Pitaevskii equation

The classical Euler-Lagrange equation leads to the Klein-Gordon (KG) equation

 □ϕ≡gμν∇μ(∂νϕ)=∂ϕV, (5)

where is the d’Alembertian operator. We choose the conformal Newtonian gauge, in which the metric of the perturbed flat FRW spacetime takes the form

 ds2=a2[−(1+2Ψ)dη2+(1−2Φ)δijdxidxj] (6)

and the Bardeen potentials are of order on the scales we are interested in.

Considering the weak field limit, where is small w.r.t. , the KG equation reduces to

 (7)

Here, is the conformal Hubble rate and a dot designates a derivative w.r.t. the conformal time . This equation simplifies drastically in the absence of any anisotropic stress. In our regime of interest, this is justified because, in the non-relativistic limit, the (second-order) anisotropic stress of the axion field (also known as “quantum pressure”) is

 πij=14m2aa2ρa(∂iρa∂jρa−13δij(∂kρa)2), (8)

where . Therefore,

 ∣∣πij∣∣∼k2m2aa2ρa. (9)

Einstein equations then imply

 ∣∣Φ−Ψ∣∣∼a2Gk2∣∣πij∣∣∼ρa(mamP)2. (10)

For an axion energy density equal to the present-day critical density, , we find that is negligible. This shows that the impact of the axion anisotropic stress on the Bardeen potentials can be safely neglected in the regime under consideration. Note that, in a realistic cosmological model including other particles such as massive neutrinos, there would be a (first order) anisotropic stress. With (which we shall refer to as the gravitational potential), the KG equation simplifies to

 ¨ϕ+2H˙ϕ−4˙Ψ˙ϕ+a2(1+2Ψ)∂ϕV=0. (11)

We retained the term for the following reason: when oscillations in the gravitational potential are taken into account, it is only suppressed by a factor of rather than relative to leading-order contribution . Conversely, we ignored the term proportional to because it contributes only at second order.

The Gross-Pitaevskii equation is derived under the assumption that the axion field undergoes coherent oscillations. Therefore, one makes the ansatz

 ϕ(η,\em x) =√2Re[ψ(η,\em x)e−imat] (12) t =∫η0dη′a(η′),

where the slowly-varying, complex envelop changes on a timescale of order . The frequency appears because we consider the non-relativistic limit, in which the typical energy of a particle is . Further assuming that the gravitational potential varies slowly (so that can be neglected), one eventually obtains

 ia(∂ηψ+32Hψ) =−12ma∂i∂iψ (13) +maa2(Ψ−18f2|ψ|2)ψ +fast oscillating piece.

Factors of arise with the partial derivatives, so that this equation is truly quantum (See for instance suarez/chavanis:2015, ; mocz/lancaster/etal:2018, , for a discussion of the “classical” limit ). Furthermore, the gravitational energy can now be recast in the form of an interaction potential appropriate to the Newtonian limit. Finally, the fast oscillatory piece involves the harmonics , and varies on a timescale much shorter than . Therefore, it is usually neglected and one is left with a nonlinear Schrödinger equation. Note that the ansatz Eq.(12) properly describes the long-term behaviour of the axion oscillation only if the axion self-interaction can be neglected.

We will now investigate the backreaction of the self-interaction and of the rapidly varying gravitational potential on the equation of motion, which arises when the axion field undergoes coherent oscillations.

## Iii Coherent oscillations of the homogeneous background

We begin with a discussion of the homogeneous KG equation. We shall see that also exhibits oscillations, which are largest at the onset of the axion coherent oscillations. However, for a realistic cosmological model including a component of relativistic particles, this is a small effect for our fiducial axion mass. This calculation will also provide us with an ingredient useful to the discussion of the inhomogeneous case, that is, the value of classical action for the background solution.

### iii.1 Free oscillations

We begin with the homogeneous KG equation

 ¨¯ϕ+2H˙¯ϕ+m2aa2¯ϕ=0. (14)

Following ratra:1991 (); Hwang/Noh:2009 (), one substitutes a solution of the form

 ¯ϕ(η)∼√2[ϕ1c(η)cos(mat)+ϕ1s(η)sin(mat)], (15)

where is the cosmic time. The assumption is that both and vary on the long timescale . The dependence of the sine and cosine on (rather than ) enable use to easily handle the mass term . At order (which quickly drops below unity after the onset of axion oscillations), both satisfy the differential equation Hwang/Noh:2009 ()

 ˙ϕ+32Hϕ=0. (16)

The axion background energy density and pressure are given by

 ¯ρa =¯ρ0+¯ρ2ccos(2mat)+¯ρ2ssin(2mat) (17) ¯Pa =¯P0+¯P2ccos(2mat)+¯P2ssin(2mat),

with

 ¯ρ0 =m2a[(ϕ21c+ϕ21s)+1maa(˙ϕ1cϕ1s−˙ϕ1sϕ1c)] (18) ¯ρ2c =m2a[1maa(˙ϕ1cϕ1s+˙ϕ1sϕ1c)] ¯ρ2s =m2a[1maa(˙ϕ1sϕ1s−˙ϕ1cϕ1c)],

and

 ¯P0 =m2a[1maa(˙ϕ1cϕ1s−˙ϕ1sϕ1c)] (19) ¯P2c =m2a[(ϕ21s−ϕ21c)+1maa(˙ϕ1cϕ1s+˙ϕ1sϕ1c)] ¯P2s =m2a[−2ϕ1cϕ1s+1maa(˙ϕ1sϕ1s−˙ϕ1cϕ1c)]

up to order .

Although this suggests that , a closer look at the Wronskian reveals that this is not the case. Namely, let and be the two independent solutions of the homogeneous equation. The Wronskian is given by

 W(y1,y2) =∣∣∣y1y2˙y1˙y2∣∣∣ (20) =maaϕ1cϕ1s+(˙ϕ% 1sϕ1c−˙ϕ1cϕ1s) ×cos(mat)sin(mat).

From a famous theorem of Abell, it must satisfy

 W(y1,y2)=~ce−2∫Hdη=~ca2, (21)

where does not depend on time for and to be linearly independent. Since and in radiation domination (RD era), while and in matter domination (MD era), the solution to Eq.(16) is (e.g., HBG:2000 (); Hwang/Noh:2009 (); Turner:1983he ())

 ϕ∝a−3/2 (22)

in both eras. This implies

 ϕ1s(η)=const⋅ϕ1c(η)+O(Hma). (23)

Explicit solutions for and are given in Appendix §A.

Eq. (23) leads to the cancellation of the term at order in and . Therefore, the slowly varying part of the pressure is only of order . By contrast, the amplitude and of the second harmonics is of order as was first pointed out by Khmelnitsky/Rubakov:2013 (). We will discuss the implications of this in Sec. §V.

### iii.2 Hubble rate

We turn to the evolution of the background FRW universe and assume that the latter is filled by the axions and by a relativistic component with homogeneous density and pressure . Hence, the Friedmann equations read

 H2 =8πG3(¯ρa+¯ρr)a2 (24) ˙H =−4πG3(¯ρa+¯Pa+2¯ρr)a2.

In light of the oscillatory pieces in and , must be of the form

 H=H0+H2ccos(2mat)+H2ssin(2mat)+… (25)

where (to be distinguished from the present-day Hubble rate ) is the dominant, slowly-varying contribution, while and both vary on a timescale . Using Eq.(16), this yields the following equations for the fast oscillating contributions to the Hubble rate,

 ˙H2c+2maaH2s=−4πGm2aa2(ϕ21s−ϕ21c)[1+O(Hma)] (26)

and

 ˙H2s−2maaH2c=8πGm2aa2ϕ1cϕ1s[1+O(Hma)] (27)

In the left-hand side, the term (resp. ) is of order relative to (resp. ). Hence, the leading order contribution to and are given by

 H2s ≈−2πGmaa(ϕ21s−ϕ21c) (28) H2c ≈−4πGmaaϕ1cϕ1s.

Using the first Friedmann equation and setting , this implies

 H2s (29) H2c ∼−32(Hma)m2aϕ1cϕ1s¯ρ0+¯ρrH0.

Since at leading order (see Sec. §III.1), we arrive at

 H2s =34(Hma)cos(2ϑ)1+ΩrΩa(a0a)H0 (30) H2c =−34(Hma)sin(2ϑ)1+ΩrΩa(a0a)H0.

In particular, in the simple adiabatic evolution outlined in Appendix §A. We will see in Sec. §V that is, in fact, equal to the value of the (real) classical action for the homogeneous solution.

In radiation domination, for which , and are suppressed by a factor of relative to . In matter domination, this suppression relative to scales as . Therefore, axion oscillations can have a significant impact on the expansion rate only if their onset occurs around the matter-radiation equality. This would happen for an axion mass . However, this range of axion mass is, as of today, ruled out by Lyman- forest measurements if all the dark matter is in the form of axions Irsic:2017yje (). Therefore, viable axion models (consistent with Lyman- forest data) have a negligible impact on the expansion rate. Hence, and can be safely ignored, except when one also considers fluctuations in the gravitational potential (see Sec. §IV.1).

### iii.3 Including a quartic interaction

While the background KG equation of a free scalar field involves only the fundamental mode in the long-time asymptotics regime considered here, interactions will generate higher-order harmonics as is well-known from anharmonic oscillators. Namely, when a quartic self-interaction is included, the KG equation for the background scalar field takes the textbook form of a driven, damped harmonic oscillator:

 ¨¯ϕ+2H˙¯ϕ+m2aa2¯ϕ=a2λ3!¯ϕ3. (31)

As is well known, substituting a solution of the form into this equation generates contributions proportional to owing to the term in the right-hand side. Therefore, one should consider instead

 ¯ϕ(η)=√2∑n=1,3(ϕnccos(nmt)+ϕnssin(nmt)), (32)

in which generally because the anharmonicity induced by the self-interaction shifts frequencies relative to the harmonic motion. Furthermore, the third harmonic is suppressed by a factor of relative to the fundamental one.

Substituting Eq. (32) into Eq. (31), retaining contributions up to linear in and discarding terms of order (which decays rapidly below unity after the onset of the axion oscillations), we arrive at the following set of algebraic (rather than differential) equations:

 (m2a−m2)a2ϕ1c =a2λ4(ϕ31c+ϕ1c% ϕ21s) (33) (m2a−m2)a2ϕ1s =a2λ4(ϕ31s+ϕ1s% ϕ21c)

for the fundamental mode, and

 (m2a−9m2)a2ϕ3c =a2λ12(ϕ31c−3ϕ% 1cϕ21s) (34) (m2a−9m2)a2ϕ3s =a2λ12(3ϕ21cϕ1% s−ϕ31s)

for the third harmonic. Hence, the relative frequency shift is given by

 (35)

while the amplitude of the third harmonic reads

 ϕ3c =−196f2(ϕ31c−3ϕ1cϕ21s) (36) ϕ3s =196f2(ϕ31s−3ϕ1sϕ21c).

The scaling implies that the frequency shifts decays as , whereas the amplitude of the third harmonic behaves like , that is, it decays quickly with the expansion of the Universe. Furthermore, both the squared frequency shift and the ratio have a present-day amplitude . They become of order unity only when approaches (signaling the breakdown of our perturbative approach). Notwithstanding, we will see in the forthcoming Section that the third and, more generally, higher order harmonics should be taken into account because they generate contributions to the self-interaction potential at all order in the coupling .

## Iv Coherent oscillations in the presence of perturbations

In this Section, we investigate how the axion coherent oscillations affect the equation of motion of the slowly-varying envelop of the axion condensate when perturbations around the background are included. We focus on the backreaction of the time-dependent gravitational potential, and on the third harmonic induced by the quartic coupling.

### iv.1 Time dependence of the gravitational potential

Let us explore first how oscillations in the axion density and pressure backreact on the gravitational potential .

To carry out this calculation, we start from the (scalar) Einstein equations

 −3H2Ψ−3H˙Ψ−∂i∂iΨ =4πGa2δρa (37) (2¨aa−H2)Ψ+3H˙Ψ+¨Ψ =4πGa2δPa,

in the longitudinal gauge, in which we set (this is justified so long as there is no significant anisotropic stress). The density and pressure perturbations are given by and , in which the background values and are evaluated as in Eqs.(18) and (19). As first recognized by Khmelnitsky/Rubakov:2013 (), oscillations of the background pressure on a timescale , see Eq.(19), leave a similar imprint on the gravitational potential. On decomposing the latter as

 Ψ =Ψ0+Ψ2ccos(2mat)+Ψ2ssin(2mat) (38) =Ψ0+Ψ2e−2imat+Ψ∗2e2imat,

we can solve for the slow- () and fast-oscillating ( and ) contributions to the potential. Separating out the different harmonics and taking into account the dominant terms solely, we obtain

 Δ\em xΨ0 =4πGa2(ρ0−¯ρ0) (39) −m2aΨ2s =πG(P2s−¯P2s) −m2aΨ2c =πG(P2c−¯P2c).

The last two relations follow from the Einstein equation for the isotropic stress. The average pressures and can generally be expressed as

 ¯P2s =−¯ρ0sin(2ϑ) (40) ¯P2c =−¯ρ0cos(2ϑ),

where, again, in the simplified model discussed in Appendix §A.

These equations show that the amplitude of and is

 (41)

Unlike the gravitational slip which is independent of scale, the suppression of the fast-oscillating potentials strongly depends on wavenumber. These become of order on the Compton scale , which is at recombination.

### iv.2 Including a quartic interaction

We are now in a position to derive the non-relativistic limit of the KG equation, Eq.(11), taking into account the feedback from both the oscillations in the gravitational potential and the quartic interaction. When the latter is included, one should consider the ansatz

 ϕ(η,\em x)=√2Re[ψ1(η,\em x% )e−imat+ψ3(η,\em x)e−3imat] (42)

as emphasized in the homogeneous case, see Sec. §III.

Because the leading order term cancels out in the harmonic, we shall retain terms up to order for , which are of the form and . By contrast, the mass term does not vanish from the harmonics. Therefore, the time-dependence of can be ignored when . We expect that the higher order harmonics neglected here can also be treated in the steady-state approximation. Regarding the gravitational potential, we retains term up to order times the leading contribution such as to include the feedback from oscillations in the gravitational potentials. The harmonics and yield the equations

 ia(∂ηψ1+32Hψ1) =−12ma∂i∂iψ1 (43) +maa2[(Ψ0−18f2|ψ1|2)ψ1 −3Ψ2ψ∗1−18f2(ψ∗1)2ψ3] ∂i∂iψ3+m2aa2 (8ψ3+112f2ψ31)=0. (44)

Owing to the third harmonic, the fast oscillatory terms usually neglected in the original Gross-Pitaevskii equation have now vanished from Eq.(43) (see also Namjoo:2017nia, ). Note also that, while the second equation is linear in the dimensionless coupling , the first includes terms up to order to emphasize that the backreaction of on is only second-order. The harmonics and furnish equations for and , which are precisely the complex conjugate of Eqs.(43) – (44). In particular, the complex potential is replaced by in those equations. Finally, the minus sign in front of is, as we shall see shortly, important for the physical interpretation of this term.

Following Desjacques/etal:2018 (), we rescale the coordinates and the fields according to

 η→1maη=~η,\em x=~\em x,ψi→mafψi=~ψi, (45) ρ→1f2ρ=~ρ,Ψi→m2aΨi=~Ψi.

Note that this is not the only possible coordinate redefinition which removes (at least part of) the characteristic scales of the problem (see, for instance, Levkov:2016rkk, , for a different choice). The prescription Eq.(45) has the advantage that corrections to the dimensionless Gross-Pitaevskii-Poisson system arise with factors of .

Upon applying this coordinate and field transformation to Eqs. (43) and (44), these equations can be recast into the form

 ia(∂ηψ1+32Hψ1) =−12∂i∂iψ1+a2[(Ψ0−18|ψ1|2)ψ1 −3Ψ2ψ∗1−18(ψ∗1)2ψ3] (46) ∂i∂iψ3 +a2(8m2aψ3+112ψ31)=0 (47)

We omitted the tildes of the new coordinates and fields to avoid clutter.

These equations must be supplemented by “Poisson” equations for the gravitational potentials , and . Assuming that both the energy density and pressure are dominated by the mass term (so that can be interpreted as a wave function), these equations read

 Δ\em xΨ0 =4π~Ga2(ρ0−¯ρ0) (48) −m2aΨ2s =π~G(P2s−¯P2s) −m2aΨ2c =π~G(P2c−¯P2c)

in the transformed coordinates, where the gravitational constant now has dimension of energy squared. Furthermore, the general expression for the energy and pressure of a scalar field yields

 ρ0 =|ψ1|2 (49) P2s =i2[ψ21−(ψ∗1)2] P2c =−12[ψ21+(ψ∗1)2].

Finally, the background density and pressure must be rescaled according to , .

Note that it is not possible to absorb, through any coordinate and field redefinition, the factor of which appears in the equation for and for the oscillatory pieces and of the gravitational potential. In other words, the relative importance of the higher harmonics and the oscillatory potential is an absolute scale that cannot be transformed away.

## V Discussion

We will now discuss the implications of the system Eqs. (46) – (48). While we focus on the evolution of the fundamental mode , we emphasize that the considerations drawn in this Section also apply to the higher-order harmonics.

### v.1 Backreaction of the higher harmonics

As emphasized above, is evaluated in steady-state approximation because, unlike , the mass term does not vanish. When the kinetic term can also be neglected  111This is similar to the Thomas-Fermi approximation (see Dalfovo:1999zz, , in the context of Bose-Einstein Condensates) which, however, is applied to the fluid representation of the GP equation, cf. Boehmer:2007um ()., which is true so long as , Eq.(47) furnishes a simple algebraic relation between the fundamental mode and the third harmonic,

 ψ3=−196m2aψ31. (50)

Once this solution is substituted into the equation for the fundamental mode, it changes the content of the square brackets in Eq.(46) into

 Ψ0ψ1−18|ψ1|2(1−196m2a|ψ1|2)ψ1−3Ψ2ψ∗1. (51)

This shows that, at order , the third harmonic “renormalizes” the self-interaction potential experienced by the fundamental mode into

 −18f2|ψ1|2ψ1→−18f2|ψ1|2ψ1+1768f4|ψ1|4ψ1, (52)

where we have momentarily switched back to dimensionfull coordinates and fields.

We expect that higher-order harmonics will introduce corrections at all orders , as discussed in Appendix §C . In other words, the effective self-interaction felt by the fundamental mode is not the quartic interaction of the “bare” Lagrangian. Therefore, the heuristic procedure advocated by Eby:2016cnq (), in which in the bare Lagrangian is replaced by in the GP equation, must be applied with care. Although these higher-order contributions to the self-interaction are suppressed by powers of , these may have implications for the stability of equilibrium configurations as investigated in, e.g., Ruffini:1969qy (); vv (); Boehmer:2007um (); Chavanis:2011zi (); Desjacques/etal:2018 (). In these studies, one either solves a differential reflecting the hydrostatic equilibrium, or attempts to minimize the energy functional using an educated guess for the solution. In all cases, the stability sensitively depends on the behavior at short distances or, equivalently, at high densities where corrections to the self-interaction become important. Obviously, higher-order bare interactions, if present, will become relevant at high densities and generate additional corrections proportional to Eby:2016cnq (); Schiappacasse:2017ham (); Chavanis:2017loo (); visinelli/etal:2018 ().

When the kinetic term in Eq.(47) cannot be neglected, the latter takes the form of the inhomogeneous Helmholtz (elliptic) equation. Namely, on defining

 θ=8m2aa2ψ3,ξ=√8maa% \em x, (53)

we can write Eq.(47) as

 Δξθ+θ=−a212ψ31. (54)

This equation can be solved using Green’s functions. For the real part of for instance, with , the general solution is

 ~θ(ξ)=~θ0(ξ)−a212∫d3ξ0G(ξ,ξ0)Re[ψ31(ξ0)], (55)

where is the solution to the homogeneous Helmholtz equation, , and the stationary wave Green’s function

 G(ξ,ξ0)=−cos(|ξ−ξ0|)4π|ξ−ξ0| (56)

is suitable for the description of the steady-state solutions we are interested in. In the case of spherical symmetry, the homogeneous Helmholtz equation reduces to the Lane-Emden equation with .

In practice, in a numerical cosmological simulation of the large scale structure, the coupled system could be solved using a relaxation approach, in which the solution to Eq.(46) with provides a good initial guess for . This trial solution is then substituted into Eq.(55) assuming the homogeneous solution vanishes, (This follows from the fact that is sourced only by a non-vanishing ), etc., until the required convergence is achieved. Note that the solution will differ from the simple scaling , so that the product takes the general form of a complex function times . This brings us to the backreaction of the gravitational potential and the presence of dissipative terms.

### v.2 Backreaction of the gravitational potential

While, in Eq.(46), the purely imaginary coefficient encodes the dilution of the (conserved) number of particles owing to the expansion, and purely real terms such as or correspond to the leading gravitational and self-interaction contributions to the energy, the term implies a generally complex coefficient and, as such, could provide a source of dissipation.

In order to further investigate this issue, we write in Eq.(46) as

 3Ψ2ψ∗1 =α[|ψ1|2ψ1+(¯P2c+i¯P2s)ψ∗1] (57) =α[|ψ1|2ψ1−¯ρ0e2iϑψ∗1] =α[|ψ1|2ψ1−¯ρ0e2i(ϑ−φ)ψ1].

Here, is a real positive constant, and is the argument or phase of , i.e. . In the last equality, the amplitude of first term relative to the dominant contribution is . Therefore, it leads to a significant correction to the usual Newtonian self-interaction on scales smaller than the Compton length. Furthermore, since it can be written as a real coefficient multiplying , it contributes to the momentum conservation equation solely once the GP system is expressed as fluid equations. Conversely, the second term cannot generally be recast in the form of a real coefficient multiplying and, thus, contributes to both the continuity and momentum conservation equations (see below). Finally, since the coupling is absent at the homogeneous level, the phase of the homogeneous solution must satisfy . This is no longer the case when inhomogeneities are present.

The Madelung transformation leads to a reformulation of Eq. (46) as (real) hydrodynamical equations Madelung:1927 (); Bohm:1952 (). The imaginary part of , proportional to , manifests itself as a source term in the continuity equation, i.e.

 ∂ηρ+3Hρ+∇\em x(ρ\em u)=2αaρ¯ρ0sin(2(⟨φ⟩−φ)). (58)

Such a term arises because the Lagrangian leading to the equation of motion Eqs.(46) is no longer invariant under a global phase transformation owing to the presence of a term