Constraints on anharmonic corrections of Fuzzy Dark Matter

The cold dark matter (CDM) scenario has proved successful in cosmology. However, we lack a fundamental understanding of its microscopic nature. Moreover, the apparent disagreement between CDM predictions and subgalactic-structure observations has prompted the debate about its behaviour at small scales. These problems could be alleviated if the dark matter is composed of ultralight fields $m \sim 10^{-22}\ \text{eV}$, usually known as fuzzy dark matter (FDM). Some specific models, with axion-like potentials, have been thoroughly studied and are collectively referred to as ultralight axions (ULAs) or axion-like particles (ALPs). In this work we consider anharmonic corrections to the mass term coming from a repulsive quartic self-interaction. Whenever this anharmonic term dominates, the field behaves as radiation instead of cold matter, modifying the time of matter-radiation equality. Additionally, even for high masses, i.e. masses that reproduce the cold matter behaviour, the presence of anharmonic terms introduce a cut-off in the matter power spectrum through its contribution to the sound speed. We analyze the model and derive constraints using a modified version of CLASS and comparing with CMB and large-scale structure data.


I. INTRODUCTION
The evidence collected over the last decades suggests that most of the matter in the universe exists in the form of dark matter (DM), whose effects have only been detected through its gravitational interaction. In particular, the assumption that dark matter is composed of non-relativistic particles, the so-called cold dark matter (CDM), has produced a remarkable concordance with the observational data over a wide range of scales and evolution epochs. It is one of the foundations of the succesful standard cosmological model ΛCDM.
Notwithstanding agreement with observations, several ingredients are lacking in our understanding of DM. In the first place, we have been unable to detect any non-gravitational interaction of DM. Most of the work in the field is currently devoted to direct, indirect detection and production searches. Owing to this effort it has been possible to tighten the parameter space of the most popular models. This lack of additional interactions makes it more difficult to discriminate between different models. There are many candidates that behave like CDM on cosmological scales, with masses ranging from the meV of the QCD axion [1] to the TeV [2,3] and going up to the 100 M of the primordial black holes [4]. The other ingredient missing is a precise understanding of the DM behaviour on small, i.e. galactic, scales. Even though most DM models mimic CDM on cosmological scales, their predictions usually differ on smaller scales [5] so they could be discriminated based only on their gravitational effects. In fact, there exist three long-standing debates, questioning the agreement between observations and the CDM theoretical predictions [6,7], the so-called 'too big to fail ' [8], 'missing satellites' [9] and especially the 'core-cusp' problem [10]. The 'core-cusp' problem refers to the discrepancy between the density profiles of CDM halos obtained in N -body simulations, that tend to be cuspy in the center, and the ones inferred from observations, that point to the existence of a central core. Although these problems are sometimes attributed to baryonic effects unaccounted for in the simulations [11], they remain one of the main challenges of the CDM model.
An interesting alternative that neatly solves the 'core-cusp' problem is fuzzy dark matter (FDM) [12]. In this picture, dark matter is composed of ultralight particles with m ∼ 10 −22 eV, so that its Compton wavelength (m −1 ) reaches astrophysical scales. Then, the formation of cusps is prevented [13]. The wave nature of the particles on the smallest scales makes them impossible to localize. While solving this problem, FDM behaves as a rapidly oscillating coherent scalar field, thus recovering the CDM behaviour on cosmological scales. In his groundbreaking work [14], Turner analyzed a homogeneous oscillating scalar field in an expanding universe. He showed that a rapidly oscillating scalar field with a power-law potential V (φ) ∝ φ n behaves as a perfect fluid with an effective equation of state w = (n − 2)/(n + 2). More general expressions can be obtained from a version of the virial theorem [15]. The results of Turner show that a massive scalar field, i.e. harmonic potential, oscillating coherently with a frequency much higher than the expansion rate behaves as CDM, at least at the background level. Afterwards, ultralight scalar fields have been thoroughly studied at the perturbation level [15][16][17][18][19], proving that the same conclusion holds. Perturbations of coherent oscillating scalar fields admit an effective fluid description with an effective sound speed nearly zero, like CDM. The main cosmological signature of these models is the supression of growth at small scales. Below some Jeans scale k −1 J the modes do not grow appreciably, translating into a cut-off in the matter power spectrum [18]. Although the work on ultralight fields has been mainly concerned with scalar fields, there are recent results on higher spin fields. It has been shown that abelian vectors at the background [20] and perturbation level [21], non-abelian vectors [22] and arbitrary-spin fields [23] behave in a similar way. Interestingly, the results of [23] show that it is possible to achieve an isotropic model of higher-spin dark matter as long as it is rapidly oscillating.
These ideas have been applied to the axion, a particularly well-motivated DM candidate. The standard QCD axion was initially proposed to solve the strong CP problem [24][25][26] in particle physics. Likewise, the appereance of many light scalar fields seems to be a generic feature of different string-theory scenarios. Some of these fields have a similar origin as the QCD axion, arising from the breaking of an approximate shift symmetry, and are usually known as axion-like particles (ALPs) or ultralight axions (ULAs) [27,28]. ALPs present similar periodic potentials but with a mass much smaller than the QCD axion that could lie in the range of ultralight fields m ∼ 10 −22 eV. While behaving like FDM, ALPs have a rich phenomenology based on their assumed interaction with matter. Aside from the standard searches for axions, there is a wealth of dedicated searches and projected experiments on the lookout for ultralight axions. These include studies of the neutral hydrogen distribution in the universe [29,30], laboratory constraints based on nuclear interactions [31], astrophysical bounds [32][33][34], gravitational wave searches [35,36] and analysis of CMB spectral distortions [37,38]. A prominent feature of the model is the presence of anharmonic corrections over the mass term in FDM. These corrections arise, to first order, as quartic corrections in the potential with the opposite sign of the mass term, i.e. attractive self-interactions. These effects have been studied, as well as the effect of the full axion potential [39,40] and their effect on the linear matter power spectrum seems to be negligible. However, self-interactions could modify non-linear structures in a significant way [41] Another possibility involves introducing a positive quartic correction, i.e. repulsive self-interactions. It is more difficult to find particle-physics models in this case [42], but the model is nonetheless well motivated as the simplest modification leading to a stable potential. This modification has been previously analyzed in some works [43][44][45][46][47]. The additional source of pressure from the repulsive self-interactions helps to solve the 'core-cusp' problem with larger masses [42]. Additionaly, unlike the axion case, it could explain the formation of vortices in galaxies [48].
In this work we will consider a fuzzy dark matter model with an additional quartic self-interaction. Using a modified version of the cosmological Boltzmann code class [49] and parameter-extraction code MontePython [50] we will constrain the parameters of the model with CMB [51] and large-scale structure (LSS) [52] data. Section II presents the model and the relevant equations for background and perturbation evolution. In section III, we review the averaging procedure when the field is rapidly oscillating and the effective fluid equations in this case. Section IV discusses a simplified model and estimates analytic bounds on the parameters, highlighting the main physical effects and the origin of the constraints on the model. In section V we present the result of the full numerical analysis and the final constraints on the model, as well as a discussion of its physical origin. Section VI summarizes the conclusions and prospects for future work.

II. EXACT EVOLUTION
Let us assume a scalar field with Lagrangian and potential in a homogeneous and isotropic universe with a flat Robertson-Walker metric in conformal time η The equation of motion for a homogeneous scalar field in this background is where H =ȧ/a and˙≡ ∂/∂η . We choose initial conditions where φ 0 is chosen to match the desired energy density Ω φ today. These are the usual initial conditions when the axion-like particles are produced through a misalignment mechanism [53] and the field starts its evolution frozen. It is important to note that the choice of initial conditions has a deep impact in the subsequent evolution. In [44], the authors considered a case similar to ours, but with an initial velocityφ = 0. In this case, there is an initial phase of stiff-matter domination, absent in our case, constrained to be short enough not to spoil BBN. We now introduce scalar perturbations over a flat Robertson-Walker metric. Following the notation of [54], the general form of the perturbations is The equation of motion for the scalar field perturbation is where the gauge has not yet been fixed. We can introduce a different parameterization, reminiscent of a perfect fluid. The components of the perturbed energy-momentum tensor are [54] δT 0 We can rewrite (8) in terms of the fluid variables, introducing δ = δρ/ρ and u = (1 + w)(v − B). In the synchronous gauge, the metric variables read and the equations of motion arė where w = P/ρ is the equation of state and the adiabatic sound speed c 2 ad is Following the analysis of [18] we provide the system with initial conditions valid up to corrections of order (kη) 4 . The scalar field starts its evolution frozen in a value φ 0 with an equation of state w −1. As the universe expands the field starts rolling down the potential, when it reaches the minimum it undergoes rapid oscillations. These oscillations occur when the effective frequency ω eff ∼ V (φ) is bigger than the friction term H, so once the scalar field starts oscillating its frequency becomes much larger than the expansion rate, the inverse of the evolution time scale of the background. On the numerical side, this means that it becomes prohibitely expensive to compute the exact evolution of the field, following every oscillation. However, the huge difference between time scales allows us to average the equations of motion and turn to an effective description.

III. AVERAGED EVOLUTION
The study of the cosmological evolution of a fast oscillating scalar was first performed in [14]. Basically, if the oscillation frequency of the scalar field is much higher than the expansion rate of the universe, the cosmological evolution becomes independent of the periodic phase of the field at leading order. Consequently, the Einstein equations can be approximately solved averaging in time the energy-momentum tensor where If the field is periodic, we can consider an integer number of periods as the integration interval. However, similar results can be reached for fast-evolving bounded solutions averaging over time spans much bigger than the inverse of its frequency but much smaller than the inverse of the expansion rate, ω −1 T H −1 . The averaging error in both cases results O(HT ).
To leading order we can drop the averages of total time derivatives, so it can be proved [19] that and with this result the effective equation of state can be written as for power-law potentials V (φ) ∝ φ n . As it can be seen, a massive scalar field, V = m 2 φ 2 /2, would behave as CDM.
For this particular case we can solve the equation of motion through a WKB expansion. Thanks to this adiabatic expansion in the parameter O(H/ma) we can perform the averages explicitly, isolating the fast evolving factor and integrating by parts, as explained in [19]. For our model, we must compute the first correction in λ. If the mass term is dominant, via a WKB expansion we can write so the first anharmonic correction to the equation of state is In this effective description the background evolution of the field is described through its density and its effective equation of state, using the conservation equationρ where for the equation of state w we will use the formula that smoothly interpolates between the radiation-like w 1/3 and matter-like w 0 behaviour. Now we can apply the same trick to the evolution of the perturbations. The equations of motion for the fluid variables arė where δ ≡ δρ / ρ , u ≡ (1 + w) v stand for averaged quantities and w, c 2 s are the effective equation of state and sound speed. To complete the system there only remains to compute the effective sound speed In contrast with the adiabatic sound speed, the sound speed is gauge-dependent. But, as we will show now, the gauge ambiguities remain of order O(H/ω eff ) so our effective sound speed turns out to be gauge-independent. In fact, identical expressions have previously been obtained working in the comoving gauge [44] and in the Newtonian gauge [19]. To leading order we have Then, using (8) and (29) we can obtain the result and finally compute the effective sound speed for a generic gauge As we anticipated, the gauge ambiguities in the metric perturbations remain of order O(H/ω eff ), so the final expression holds in any gauge. Moreover, it can be rewritten in a manifestly gauge-invariant form substituting δφ by its gaugeinvariant perturbation [54] and using the relations we obtain This expression agrees with the result obtained in [19] working in the Newtonian gauge, so the same conclusions apply. In particular, a generic feature of this kind of models is a suppression of growth c 2 s 1 for small scales k ω eff . In the case of a power-law potential V (φ) = C n φ n , for large scales k ω eff we have For a harmonic potential n = 2, the zero-order term drops out and we must calculate the first-order corrections in k.
Our potential of interest is a polynomial V (φ) = 1 2 m 2 φ 2 + λ 4 φ 4 , a mass term plus an anharmonic correction, in this case we have [19] where ρ is the energy density of the scalar field and the anharmonic correction is assumed to be small. In our numerical solution we will use an effective sound speed suggested by the form of (32) and that smoothly interpolates between all the regimes of interest.

IV. HEURISTIC CONSTRAINTS ON THE NON-HARMONIC CONTRIBUTION
In this section we will discuss the simplest limits that constrain the model. With this objective let us assume a simple cosmology composed of radiation, cosmological constant and our scalar field where Ω i = 8πGρ i /(3H 2 0 ) are the abundances with i ∈ {φ, rad, Λ} which correspond to scalar field, radiation and cosmological constant respectively.
• Limits on λ from background evolution. The position of the peaks in the CMB temperature spectrum, especially the first one, is very sensitive to the amount of matter and the redshift of equality z eq . We can assume that to have a viable model of dark matter this quantities remain essentially the same as in ΛCDM. In this case, to have a dark matter behaviour that resemble CDM the anharmonic corrections at this time should be small This imposes an upper limit on λ, namely excluding the orange region in Figure 1.
• Limits on m from perturbation evolution. If λ is small enough, the background evolution of the effective fluid is identical to ΛCDM. In this case we can obtain limits from the behaviour of the perturbations. From (26) and (27) it can be seen that if we neglect the expansion rate, c 2 s k 2 H 2 , density perturbations evolve according tö producing an oscillatory behaviour instead of the standard growth. To avoid a clear disagreement with observations, the effect of a non-negligible sound speed must be small Translating into a lower bound in the allowed masses as before, we assume that z eq corresponds to the standard value. To obtain a conservative limit we choose k = 0.2 Mpc −1 , the highest mode observed in LSS at the linear level, so that excluding the blue region in Figure 1.
• Observable effects of anharmonic corrections. Finally, there is a region in the parameter space that we cannot yet exclude and where the effects of anharmonic corrections to the sound speed may be important Imposing that the second term dominates over the harmonic contribution yields an upper bound on λ corresponding to a region that we cannot exclude right away, but where effects of the anharmonic corrections to the sound speed are to be expected.

An additional result that can be obtained from (44) is the Jeans wavenumber
Sub-Hubble modes below this Jeans wavenumber, k < k J , grow while modes with k > k J are suppressed. In the massive case with λ = 0 we obtain Now, since we have seen that the quartic correction affects the sound speed, it will also affect the Jeans scale. It is natural to ask what combination of parameters (m, λ) can have a similar impact on structure formation as the case (m,λ = 0). To this end, we look for the combination that gives the same Jeans scale at the matter-radiation equality. Since its scaling in time is not significantly modified, this simple estimate should capture the essential features of structure formation in both models. Equating both sound speeds and inserting the result (50) we have an estimate for λ λ = 4.96 × 10 −100 m 10 −24 eV This simple result suggests for instance that, at the linear level, structure formation should be similar in the models (m = 10 −26 eV,λ = 0) and (m = 10 −24 eV, λ 4.96 × 10 −98 ), a result that we will check with the full numerical solution. This estimate is represented in Figure 1 for two different massesm. After discussing some approximate bounds on our model and its physical origin, we will devote the next section to the full numerical solution.

V. NUMERICAL EVOLUTION AND CONSTRAINTS
We modify the publicly available Boltzmann code class [49] and include this ultralight scalar field as a new species, that will assume the role of dark matter. Now, we summarize the key changes in the code and the evolution scheme chosen for the scalar field.
• At the background level, we start solving the equation (4) with initial conditionsφ = 0 and φ = φ 0 . The initial value φ 0 is chosen internally with a built-in shooting algorithm such as to match the energy density Ω φ (today) required. As a technical aside, it is critical to start with a sensible initial guess for φ 0 , so that the shooting algorithm converges quickly. In [18] the authors provide analytical formulae for the initial guess in the harmonic case, that works as well if the anharmonic corrections are small. If the quadratic and quartic terms are comparable it is more difficult to find analytical expressions that fit our purposes. In our case, we precompute an interpolation table for different values of m, λ and φ 0 yielding some value Ω φ (m, λ, φ 0 ). We only compute a coarse table, so that we still use the shooting algorithm to adjust φ 0 and achieve the desired precision in Ω φ .
With the initial conditions provided, the field starts its evolution frozen, slowly rolling down the potential until its natural frequency term in (4) dominates and it undergoes rapid oscillations. In this case it is computationally expensive to follow every oscillation so we turn to the averaged equations when V (φ) > 3H.
In the averaged regime, we solve (24), matching continuously with the solution in the exact regime, and compute the pressure using the effective equation of state (25).
• At the perturbation level, we first solve (12) and (13) with adiabatic initial conditions δ = u = 0. For each mode k we start the integration early enough to ensure that we start well within the exact regime, V (φ) 3H. In the averaged regime, V (φ) > 3H, we solve the equations (26) and (27) with the sound speed given by (39).
Some results for temperature and matter power spectra are shown in Figures 2 and 3. They show the impact of different choices of m and λ, while the other cosmological parameters are fixed to their Planck [51] best-fit values. As anticipated, the main cosmological signature is the appearance of a cut-off in the matter power spectrum. This cut-off has already been discussed in the harmonic case [18]. In our case, we see that the anharmonic terms produce a similar effect. On the left, results for a massive scalar field without self-interaction. On the right, results for different self-interaction strengths for a mass that is indistinguishable from CDM with λ = 0.

A. Physical effects
The main physical effect responsible for the appearance of a cut-off in the matter power spectrum has already been discussed. In the averaged regime, the scalar field that supplies the dark matter component behaves like a fluid with a non-negligible sound speed. On small scales, above a certain Jeans scale k J , the density perturbations oscillate and the growth is suppressed. This effect is illustrated in the Figure 4 for modes above and below k J .
In the case of the CMB temperature power spectrum, it is far more difficult to disentangle the physical effect responsible for each feature. We split the effects in two categories, those coming from the modified background evolution and those coming from the perturbations. Furthermore, we will refer to two extreme cases (m = 10 −27 , λ = 0) and (m = 10 −24 , λ = 10 −97 ) as m-case and λ-case respectively. To gain some insight into the CMB spectrum structure, we will rely on simplified, analytical estimates [55][56][57] and work in the Newtonian gauge.

Background evolution
The modified equation of state (25) changes the background evolution, modifying in particular the redshift of matter-radiation equality z eq and in general the expansion history a(τ ). In the m-case, the field transitions directly from the frozen value w −1 to a matter-like phase, while in the λ-case there is an intermediate radiation-like phase. There are two key effects • First peak position. The position peak of the first peak can be estimated as where the angular diameter d a distance is defined as d s | dec is the sound horizon of the photon-baryon plasma evaluated at decoupling and the sound speed for the baryon-photon plasma is The angular diameter distance is almost unaffected but the sound horizon is slightly modified. Compared to ΛCDM we obtain relative deviations on peak of about +2%, shift to the left, in the m-case and −0.7%, shift to the right, in the λ-case. Both are compatible with the tiny deviations observed in Figure 2.
• Damping envelope. Another physical scale that is modified is the diffusion length where is the Thomson scattering rate. The diffusion length governs the damping envelope, e −( / D ) 2 , through the relation For a reference multipole = 820, corresponding to the third acoustic peak, in the λ-case we obtain a modified damping envelope that produces an enhancement of 6% compared to ΛCDM, that can explain the overall increase of power in Figure 2. For the m-case, we obtain the puzzling result of a suppression of 0.7%, in clear disagreement with the observed effect. However, we will shortly see how a novel effect in the perturbation evolution can account for this overall amplification.

Perturbation evolution
In the tightly coupled regime, the photon temperature fluctuation evolves according tö In the standard scenario, ignoring slow changes in R, φ and ψ from the expansion, we havë This produces an oscillatory pattern with frequency ω = kc s γ and zero-point displaced by an amount −(1 + R)ψ. The main part of the temperature Sachs-Wolfe effect comes from the contribution |Θ 0 + ψ| 2 | dec , so the displacement of the zero-point of the oscillations gives the characteristic asymmetry between odd and even peaks in the CMB temperature spectrum. Our modification of dark matter produces two interrelated effects, oscillation and suppression of growth at small scales.
• Effects of suppression of growth at small scales. The suppression of dark matter density perturbations at small scales also suppresses the gravitational wells ψ, shifting the zero-point of the oscillation back to zero. This effect, alone, reduces the asymmetry among the peaks, decreasing the odd and increasing the even peaks. This explains the characteristic enhancement of the second peak with respect to the third one in Figure 2.
• Effects of oscillatory behaviour. There only remains to explain one effect, the striking gain in peak amplitude in the m-case. According to the modification in the damping envelope, the peaks should be slightly suppressed and their enhancement is actually related to a resonance effect. In the standard scenario, the term ψ behaves like a constant external force, shifting the equilibrium position of the photon oscillations. In our case, it is not constant anymore, but oscillates with a frequency kc s given by the sound speed of the dark matter perturbations (39). These two frequencies, kc s and kc s γ , are comparable for a range of k values, as shown in Figure 5, producing a resonant effect that increases the height of the peaks, as shown in Figure 6. Moreover, since according to (39) the scale of the crossover in Figure 5 evolves ∝ a, as we go from decoupling back in time it moves to smaller k. That is to say, although the crossover at decoupling is located around k 0.1 Mpc −1 , smaller k have also fulfilled the resonance condition at previous times, so they have also got amplified.

B. Observational constraints
To compare this model with CMB and LSS observations and refine the heuristic constraints obtained in section IV, we use the public parameter-extraction code MontePython [50]. We will compare our results with two different data sets: CMB measurements by Planck and large-scale structure information by WiggleZ [52]. We perform two analysis, Planck only and Planck+WiggleZ. In each case we vary the six ΛCDM base model parameters, in addition to the foreground parameters, plus m and λ, the mass and anharmonic parameter. We choose logarithmic priors in our model parameters, as shown in Table I.
It is important to note that to perform an accurate comparison with LSS data we must restrict our analysis to linear scales k < ∼ 0.2 h/Mpc. The non-linear module in class includes HaloFit [58], but since it has not been calibrated for our model we restrict our analysis to linear scales without non-linear corrections. It is to be expected that, in the future, as more N -body simulations with ultralight fields become available, non-linear information will allow us to tighten the constraints.
We do not observe any significant degeneracy between m, λ and the rest of cosmological parameters. Best-fit results are shown in Table II, while the marginalized countour for our model parameters is represented in Figure 7. I: Prior ranges on the base ΛCDM parameters and the model parameters m and λ. A symbol − means that there is no prior. Additionally, the fixed parameters include the neutrino properties. In our case, two massless neutrinos plus a massive one with m = 0.06 eV, such that N eff = 3.046 and mν /Ων = 93.14 eV.
Parameter minimum maximum

VI. CONCLUSIONS
The presence of self-interactions in the ultralight field potential can lead to the appearance of new backgroundevolution phases, like the radiation-like due to our quartic potential. This modified background evolution, and especially its critical effect on the sound speed of dark matter perturbations, can lead to significant differences from observations. The observational signatures of the anharmonic contribution are similar to the mass term, the most prominent being the appearance of a cut-off in the matter power spectrum. This produces constraints for masses that would be otherwise indistinguishable from CDM, i.e. m > ∼ 10 −24 eV. Our constraints on λ complement other bounds present in the literature, e.g. [46]. This bounds on λ follow a scaling law with m 4 according to (42). We can extrapolate the results to higher masses using the 2σ region of Figure 7, obtaining an approximate constraint on λ log 10 (λ) < −91.86 + 4 log 10 m 10 −22 eV , for masses m > 10 −24 eV. So far, we have only analyzed linear observables, but in fact larger effects on non-linear scales are expected. The available parameter space could be further constrained in the future using cosmological information with non-linear observables, as more simulations with ultralight fields become available. Even without non-linear information, using the formula (51) one could put forward the proposal that similar results on structure formation could be obtained for higher masses with a positive λ. For instance, results form = 10 −22 eV might be reproduced with masses m 10 −5 eV adding a self-interaction of the order of λ 10 −24 , very close to the limit that can be obtained from (61). Nevertheless, a definitive answer to this suggestive proposal require a fully non-linear analysis.