Cosmological signatures of ultralight dark matter with an axionlike potential

TL;DR

This paper investigates how nonlinearities in axion-like dark matter potentials affect cosmological evolution, revealing delays in oscillation onset and small-scale instabilities that could influence structure formation.

Contribution

It provides the first detailed analysis of nonlinear potential effects on axion-like dark matter using the CLASS code, highlighting new dynamical features.

Findings

Delay in axion field oscillation onset due to nonlinearities

Identification of tachyonic instability in linear modes

Potential constraints on axion parameters from Ly-α observations

Abstract

Nonlinearities in a realistic axion field potential may play an important role in the cosmological dynamics. In this paper we use the Boltzmann code CLASS to solve the background and linear perturbations evolution of an axion field and contrast our results with those of CDM and the free axion case. We conclude that there is a slight delay in the onset of the axion field oscillations when nonlinearities in the axion potential are taken into account. Besides, we identify a tachyonic instability of linear modes resulting in the presence of a bump in the power spectrum at small scales. Some comments are in turn about the true source of the tachyonic instability, how the parameters of the axionlike potential can be constrained by Ly-$\alpha$ observations, and the consequences in the stability of self-gravitating objects made of axions.

Figures

Equations

$$V(\varphi )=m_{\varphi }^2{f}^2\left[1+\cos \left(\varphi /f\right)\right],$$

$$H^2$$

$$\dot{H}$$

$$\ddot{\varphi }$$

$$\frac{\kappa \dot{\varphi }}{\sqrt{6}H}\equiv {\Omega }_{\varphi }^{1/2}\sin (\theta /2),\frac{\kappa V^{1/2}}{\sqrt{3}H}\equiv {\Omega }_{\varphi }^{1/2}\cos (\theta /2),$$

$$y_1\equiv -2\sqrt{2}\frac{{\partial }_{\varphi }{V}^{1/2}}{H},$$

$${\theta }^{\prime }$$

$$y_1^{\prime }$$

$${\Omega }_{\varphi i}=a_{\mathrm{osc}}^{-3}{a}_i\frac{{\Omega }_{\varphi 0}}{{\Omega }_{r0}},y_{1i}=5{\theta }_i,\frac{m_{\varphi }^2}{H_i^2}=\frac{y_{1i}^2}{4}+\lambda {\Omega }_{\varphi i},$$

$$\ddot{\phi }=-3H\dot{\phi }-\left[\frac{k^2}{a^2}-m_{\varphi }^2\cos (\varphi /f)\right]\phi -\frac{1}{2}\dot{\varphi }\dot{\bar{h}},$$

$$\sqrt{\frac{2}{3}}\frac{\kappa \dot{\phi }}{H}\equiv -{\Omega }_{\varphi }^{1/2}{e}^{\alpha }\cos (\vartheta /2),\frac{\kappa y_1\phi }{\sqrt{6}}\equiv -{\Omega }_{\varphi }^{1/2}{e}^{\alpha }\sin (\vartheta /2),$$

$${\delta }_0^{\prime }$$

$${\delta }_1^{\prime }$$

Full text

I Introduction

Modern cosmological observations have brought about a large amount of data Ade et al. (2016); Albareti et al. (2016), making it possible to constrain, with high accuracy, theoretical models describing the Universe at large scales. The so-called Lambda Cold Dark Matter ($\Lambda$CDM) paradigm is very successful at reproducing cosmological observations Ade et al. (2016) but it requires a dark matter (DM) component ($\approx 26\%$), effectively described by collisionless particles that interacts mostly gravitationally with other matter componentsLiddle and Lyth (1993); Eisenstein and Hu (1997); Peebles (2017). However, there are longstanding discussions about how well the $\Lambda$CDM model describes the Universe at galactic and sub-galactic scalesWeinberg et al. (2014); Garrison-Kimmel et al. (2014). The solution to these problems may come from the specific properties of the DM, or from an interplay between the DM properties and kinematic effects with baryons, but still the incompleteness of galactic observations may impair our ability to infer the DM distribution properties from them. Given the current status, the further development of theoretical models still is very much desirable if one is to elucidate the properties of this matter component of the Universe.

According to recent studies, axion DM has become a compelling candidate to replace CDMHui et al. (2017); Marsh (2016); Magana and Matos (2012), and even some experiments have been already set up to have a direct detection of this elusive particle Asztalos et al. (2010); Avignone et al. (1998); Bernabei et al. (2001); Morales et al. (2002); Zioutas et al. (2005); Bradley et al. (2003). In particular, there are several proposals for detection of ultralight axions (ULA) using laser interferometers Aoki and Soda (2016), analyzing the frequency and dynamics of pulsars Khmelnitsky and Rubakov (2014); Blas et al. (2017), and also in gravitational wave detectors Aoki and Soda (2017). Nonetheless, there are still many open questions such as what is the right axion mass limits one can place by using, for instance, galactic kinematicsUreña López et al. (2017); Gonzalez-Morales et al. (2016); Chen et al. (2016); Calabrese and Spergel (2016); Marsh and Pop (2015), and Lyman-$\alpha$ observationsIr$\check{{\rm s}}$i$\check{{\rm c}}$ et al. (2017); Armengaud et al. (2017). At the cosmological level, axion models have been studied considering it as a free scalar field, i.e., as a scalar field endowed with a quadratic potential $V(\phi)=m^{2}\phi^{2}/2$ Ureña López and González-Morales (2016); Marsh and Ferreira (2010); Marsh et al. (2012a); Marsh and Silk (2013); Hlozek et al. (2015); Marsh et al. (2012b). However, a more realistic form of the axion potential is the trigonometric one,

[TABLE]

where $m_{\phi}$ is the axion mass and $f$ is the decay constant of the axion. The axion potential (1) originally arose in QCD with the aim to solve the strong CP problem Peccei and Quinn (1977a, b); Weinberg (1978); Wilczek (1978), and the potential of such field arises from non-perturbative effects which generate a periodic behavior after the breaking of the Peccei-Quinn symmetry $U(1)_{PQ}$ due to instantonsSikivie (2008); Cheng and Kaplan (2001); Dine and Fischler (1983). More recently, it has been argued that axions emerge in string theories from the compactification of extra dimensionsSvrcek and Witten (2006); Olsson (2007); Arvanitaki et al. (2010).

Given the motivations above, our aim in this paper is to study the axion field as source of DM with its corresponding trigonometric potential (1). For that purpose, we present, for the first time, an analysis of the cosmological evolution, from radiation domination up to the present day, of an axion field taking into account the whole properties of the potential (1). This is accomplished by: $\left.1\right)$ transforming the standard cosmological equations for both, the background and the linear perturbations into a dynamical system, and $\left.2\right)$ using an amended version of the Boltzmann code CLASSLesgourgues (2011) to obtain accurate numerical solutions. We analyze the differences in the linear process of structure formation of the axion field with respect to the free (quadratic potential) and the CDM cases. For the sake of concreteness we present all the results for a fiducial model with axion mass $m_{\phi}=10^{-22}$ eV, but we have verified that the qualitative features hold for other masses in the range, $10^{-26}<m_{\phi}/$eV$<10^{-20}$.

II Background Dynamics

The equations of motion for a scalar field $\phi$ endowed with the potential (1), in a homogeneous and isotropic space-time with null spatial curvature, are given by

[TABLE]

where $\kappa^{2}=8\pi G$, $\rho_{j}$ and $p_{j}$ are the energy and pressure density of ordinary matter, a dot denotes derivative with respect to cosmic time $t$, and $H=\dot{a}/a$ is the Hubble parameter. The scalar field energy density and pressure are given by the canonical expressions $\rho_{\phi}=(1/2)\dot{\phi}^{2}+V(\phi)$ and $p_{\phi}=(1/2)\dot{\phi}^{2}-V(\phi)$.

We define a new set of polar coordinates as in Copeland et al. (1998); Ureña López and González-Morales (2016); Ureña López (2016),

[TABLE]

with which the Klein-Gordon equation (2c) takes the form of the following dynamical system:

[TABLE]

where $\lambda=3/\kappa^{2}f^{2}$ and $\Omega_{\phi}=\kappa^{2}\rho_{\phi}/3H^{2}$ is the standard scalar field density parameter. Here, a prime denotes derivative with respect to the number of $e$-foldings $N\equiv\ln(a/a_{i})$, with $a$ the scale factor of the Universe and $a_{i}$ its initial value, and the total equation of state $w_{tot}=p_{tot}/\rho_{tot}$. For $\lambda=0$ in Eq. (4) the dynamical system for the free case is recovered; see Ref. Ureña López and González-Morales (2016).

One critical step in the numerical solution of Eqs. (2) and (4) is to find the correct initial conditions of the dynamical variables. For the axion case, it can be shown that we must satisfy the following constraints,

[TABLE]

where $a_{\rm osc}$ is the value of the scale factor at the onset of the oscillations of the field $\phi$ around the minimum of the potential (1). The solution of Eqs. (5) provides appropriate seed values that the CLASS code adjusts through a shooting procedure to obtain the correct value of the axion density parameter $\Omega_{\phi 0}$ at the present time.

In Fig. 1 we show the evolution of the axion energy density $\rho_{\phi}$ in comparison with that of CDM (all other cosmological quantities are the same as in the fiducial $\Lambda$CDM modelAde et al. (2016)). The numerical examples correspond to $\lambda=0,10,10^{2},10^{3},10^{4},10^{5}$. We can clearly see that $\rho_{\phi}$ evolves just like CDM after the onset of the field oscillations. The latter are delayed by the presence of the decay parameter $\lambda$, and also the transition to the CDM behavior occurs more abruptly for larger values of $\lambda$. This is just a consequence of the increase in the steepness of the potential (1) for $\lambda\gg 1$, which in turn makes it more difficult to find a reliable numerical solution of Eqs. (4). The largest value considered for the decay parameter was $\lambda=10^{5}$. Although larger values would be desirable, we are already close to the expected upper bound on $\lambda$. As estimated in Ref. Diez-Tejedor and Marsh (2017), the axion field can provide the whole of the DM budget as long as $m_{\phi}/\sqrt{\lambda}>6\times 10^{-27}$ eV. In particular, a conservative estimate is that $\lambda\lesssim 10^{8}$ if $m_{\phi}=10^{-22}$ eV, although other considerations can provide stronger constraints Diez-Tejedor and Marsh (2017); Visinelli (2017); Arias et al. (2012).

III Linear Density Perturbations and Mass Power Spectrum

Let us now consider the case of linear perturbations $\varphi$ of the axion field in the form $\phi(x,t)=\phi(t)+\varphi(x,t)$. As for the metric, we choose the synchronous gauge with the line element $ds^{2}=-dt^{2}+a^{2}(t)(\delta_{ij}+h_{ij})dx^{i}dx^{j}$, where $h_{ij}$ is the tensor of metric perturbations. The linearized Klein-Gordon equation for a given Fourier mode $\varphi(k,t)$ reads Ratra (1991); Ferreira and Joyce (1997, 1998); Perrotta and Baccigalupi (1999):

[TABLE]

where a dot means derivative with respect the cosmic time, $\bar{h}={h^{j}}_{j}$ and $k$ is a comoving wavenumber.

As shown in Ref.Ureña López and González-Morales (2016), we can transform Eq. (6) into a dynamical system by means of the following (generalized) change of variables,

[TABLE]

with $\alpha$ and $\vartheta$ the new variables needed for the evolution of the scalar field perturbations. But if we further define $\delta_{0}=-e^{\alpha}\sin(\theta/2-\vartheta/2)$ and $\delta_{1}=-e^{\alpha}\cos(\theta/2-\vartheta/2)$, then Eq. (6) takes on a more manageable form,

[TABLE]

where $k_{J}^{2}\equiv a^{2}H^{2}y_{1}$ is the (squared) Jeans wavenumber and a prime again denotes derivative with respect to the number of $e$-folds, $N$. Notice that the new dynamical variable $\delta_{0}$ is the axion density contrast, as a straightforward calculation using Eqs. (3) and (7) shows that $\delta\rho_{\phi}/\rho_{\phi}=(\dot{\phi}\dot{\varphi}+\partial_{\phi}V\varphi)/\rho_{\phi}=\delta_{0}$. This implies that Eq. (8a) is the closest expression one can find to a fluid equation for the evolution of the axion density contrast. The physical interpretation of $\delta_{1}$ is by no means as direct as that of $\delta_{0}$, and then Eq. (8b) tells us of the difficulties to match the equations of motion of scalar field linear perturbations to those of a standard fluidHu (1998). For the initial conditions, we use the attractor solutions at early timesUreña López and González-Morales (2016) given by $\delta_{0i}=-\bar{h}_{i}\theta^{2}_{i}/84$ and $\delta_{1i}=-\bar{h}_{i}\theta_{i}/7$, where $\bar{h}_{i}$ and $\theta_{i}$ are, respectively, the initial values of the trace of metric perturbations $\bar{h}$ and the background angular variable $\theta$.

The solution of Eqs. (8) are useful to build up cosmological observables such as the mass power spectrum (MPS), which we show for the axion field and CDM in Fig. 2. It is well known that there is a characteristic cut-off in the MPS of a free field, and this feature is also present for the axion case, although the cut-off is shifted towards smaller scales (larger wavenumbers). But more prominently, the axion MPS presents an excess of power, even compared to CDM, at scales close to the cut-off if $\lambda\gg 1$. Such excess was reported before in Refs.Zhang and Chiueh (2017a, b) (see also Ref. Suarez and Matos (2011) for an early indication of such power excess in scalar field models) and attributed to the so-called extreme initial conditions of the background field (under our approach, this means $\phi/f\to 0$). As we shall explain below, the excess should be rightfully attributed to the extreme value of $\lambda\gg 1$ [which in turn has an effect on the initial conditions via Eqs. (5)], and then ultimately to the decay constant $f$.

Also shown in Fig. 2 are the free (with mass $m_{\phi}^{*}=3.635\times 10^{-22}$eV and $\lambda=0$) and extreme cases (with $\lambda=1.3\times 10^{5}$) whose MPS differs by $50\%$ from that of CDM at the same wavenumber, namely $k_{50\%}=11.218h/$Mpc. However, it is important to highlight that both cases have a very different behavior at smaller and larger values of $k$, which means that the axion MPS is non-degenerate with respect to that of the free case. Moreover, this also shows that the axion case $(m_{\phi},\lambda\neq 0)$ cannot be exactly matched to a free case $(m_{\phi}^{*},\lambda=0)$.

As for the excess of power at some scales in the MPS, we first note that the presence of $\lambda$ in Eq. (8b) defines an effective wavenumber in the form $k^{2}_{eff}=k^{2}-\lambda a^{2}H^{2}\Omega_{\phi}/2$, which, in contrast to the ratio $k^{2}/k^{2}_{J}$ that appears in the free case, it could be positive as well as negative. Taking advantage of the similarities with the free caseUreña López and González-Morales (2016), we will study the homogeneous solutions of Eqs. (8) (without the driving terms) after the onset of the rapid oscillations of the axion field. We first discard all the trigonometric terms, and then Eq. (8) can be written as: $\delta^{\prime}_{0}\simeq-(k^{2}/k^{2}_{J})\delta_{1}$ and $\delta^{\prime}_{1}\simeq(k^{2}_{eff}/k^{2}_{J})\delta_{0}$. Under the assumption that both functions $k^{2}_{eff}$ and $k^{2}_{J}$ are approximately constant, we obtain that the density contrast has a harmonic solution of the form $\delta_{0}\sim C_{0}\cos(\omega N)$, where the (squared) fundamental frequency is $\omega^{2}=k^{2}k^{2}_{eff}/k^{2}_{J}$ and $C_{0}$ is an integration constant. Just like in the free case, it can be seen that if $0<|\omega^{2}|\ll 1$ the homogeneous solution of Eqs. (8) becomes irrelevant and then the axion density contrast can grow similarly to that of CDM. Similarly, if $1\ll\omega^{2}$ the growth of the density contrast is strongly suppressed and there appears a sharp cut-off in the MPS at large $k$ (small scales). But now we must also consider the possibility that $\omega^{2}\ll-1$, for which the homogeneous solution changes to $\delta_{0}\sim C_{0}\cosh(|\omega|N)$, and then the growth of the given mode $k$ is even enhanced beyond the CDM case. We dub the latter effect as the tachyonic instability of linear perturbations.

To determine the linear modes that suffer a tachyonic instability we proceed as follows. We note that both the Jeans wavenumber $k_{J}$ and the effective wavenumber $k_{eff}$ are functions of background quantities only, and then their evolution can be easily calculated for different values of $\lambda$ and $k$. This is shown in Fig. 3, where we see that only a limited range of $k$, and for a finite lapse of time after the onset of the rapid oscillations of the axion field, will be affected by the tachyonic instability $\omega^{2}\ll-1$. The wavenumbers shown in Fig. 3 constitute a representative set of modes that allow us to get a better comprehension of the bump in the MPS for the extreme case $\lambda=10^{5}$. Large scales with $k<10^{-2}$Mpc*-1* are not affected by the tachyonic instability because for them the condition $\omega^{2}\ll-1$ is never satisfied. Modes with $k>10^{-2}$Mpc*-1* start to be affected as the condition $\omega^{2}\ll-1$ is satisfied at the onset of the oscillations of the axion field, but the time lapse of the tachyonic effect is reduced as the wavenumber increases (in other words, it takes less and less time for $\omega^{2}$ to change from negative to positive again), so that for small scales with $k>30$Mpc*-1* the tachyonic instability never happens. In fact, $\omega^{2}\gg 1$ at all times for those latter modes and the result is that the amplitude of their perturbations must be highly suppressed. Thus, we infer from Fig. 3 that the tachyonic instability only affects the modes with $10^{-2}<k/\textrm{Mpc}^{-1}<30$, which includes the range of wavenumbers responsible for the bump in the MPS in Fig. 2, i.e., approximately $1<k/\textrm{Mpc}^{-1}<10$.

IV Discussion and conclusions

We fully computed the MPS for the axion potential and, within the range of physical parameters that we were able to explore, showed that its features do change significantly in the case $f/m_{Pl}\ll 1$ ($\lambda\gg 1$), for which linear perturbations in a certain range of wavenumbers suffer a tachyonic instability and are able to grow more than their CDM counterparts. This causes the appearance of a bump in the MPS which is close to the cut-off scale, which in turn is also displaced towards larger wavenumbers in comparison to the free case. Our results are in agreement with the semi-analytical studies of the axion field in Refs. Zhang and Chiueh (2017a, b); Suarez and Matos (2011) (see also Diez-Tejedor and Marsh (2017)), which were the first to suggest the existence of a bump in the MPS. However, we were able to show that such effect results from the condition $\lambda\gg 1$Suarez and Matos (2011), rather than from extreme initial condition $\phi_{i}/f\ll 1$ as suggested in Refs. Zhang and Chiueh (2017a, b).

Just recently a set of new constraints on the axion mass based on the analysis of Lyman-$\alpha$ forest had been presented. The strongest constraint comes from high resolution spectra, implying $m_{\phi}>37.5\times 10^{-22}$eV  Ir$\check{{\rm s}}$i$\check{{\rm c}}$ et al. (2017) and $m_{\phi}>29\times 10^{-22}$eV Armengaud et al. (2017), at the 2-$\sigma$ confidence level. To extrapolate such constraints to the full axion potential is not straightforward. In Fig. 4 we show the 1-dimensional MPS ($P^{1{\rm D}}$) for the full axion potential, relative to that of the $\Lambda$CDM. The $P^{1{\rm D}}$ is closely related to the flux power spectrum, $P_{\rm F}$, that is the actual observable in surveys like BOSS Palanque-Delabrouille et al. (2013), HIRES/MIKESViel et al. (2013) and XQ100López et al. (2016). For reference we have included horizontal lines at the approximate precision at which such experiments can prove the $P^{1{\rm D}}$, and the $k$-range they cover. This can be read as indication that such observations would be able to set constraints in the axion mass vs decay parameter ($m_{\phi}$ vs $\lambda$) plane, provided that the actual quantity one have to predict is $P_{\rm F}$ following an analysis similar to that in Ir$\check{{\rm s}}$i$\check{{\rm c}}$ et al. (2017); Armengaud et al. (2017). It will also be interesting to know whether future surveys as the Dark Energy Spectroscopic Instrument (DESI)Levi et al. (2013) and the Large Synoptic Survey Telescope (LSST)Ivezic et al. (2008), in case a bump and a cut-off in the MPS, or the $P^{1{\rm D}}$ are detected, will also be able to spot the differences between the free case and the full axion one and in turn put constraints on the decay parameter $\lambda$.

Some discussion about the consequences of a so-called extreme axion case in the formation of cosmological structure has been put forward in Ref. Schive and Chiueh (2017), by considering the equivalence between N-body simulations and the Schrodinger-Poisson system first hinted at in Ref. Schive et al. (2014). From our perspective, the formation of structure under the extreme axion case is better captured by a Gross-Pitaevskii-type of equation with a negative-definite quartic self-interaction given by $g_{4}=-m^{2}_{\phi}f^{-2}/6$ (the coefficient $g_{4}$ can be read off from the series expansion of potential (1) up to the fourth order: $V\simeq(m^{2}_{\phi}/2)\phi^{2}-(m^{2}_{\phi}f^{-2}/4!)\phi^{4}$. According to diverse studies in Refs. Guzman and Urena-Lopez (2006); Chavanis (2011); Levkov et al. (2017), the parameter that quantifies the strength of the self-interaction is the combination $g_{4}/(4\pi Gm^{2}_{\phi})=-\lambda/3$. In the extreme case $\lambda\gg 1$, equilibrium gravitational configurations of the coupled Gross-Pitaevskii-Poisson system present stable and unstable branches (see also Refs. Barranco et al. (2015); Helfer et al. (2017) for the relativistic axion case), and the critical quantities at the transition point between the two branches have been found to be $\phi_{c}/f\simeq\lambda^{-1/2}$ (for the central field value) and $M_{c}\simeq\lambda^{-1/2}\,m^{2}_{\rm Pl}/m_{\phi}$ (for the total mass), where $m_{\rm Pl}$ is the Planck mass Guzman and Urena-Lopez (2006); Chavanis (2011); Levkov et al. (2017). On one hand, stable configurations then correspond to field values $\phi/f\leq\lambda^{-1/2}$, and then the gravitational stability of an axion configuration requires a more diluted field for larger $\lambda$. At the same time, the critical total mass $M_{c}$ also decreases for larger $\lambda$, and for the fiducial model considered throughout we find $M_{c}\sim 10^{9}\,M_{\odot}$ if $\lambda=10^{5}$. As already noted in Refs. Levkov et al. (2017), this means that even the less massive halo objects in a typical simulation (see for instance Schive et al. (2014)) would be in risk to collapse into black holes.

All of the above lead us to wonder about the possibility of having $\lambda<0$, so that the quartic self-interaction $g_{4}$ is strictly positive definite. In such a case, the gravitational stability of bounded objects is instead enhanced by the presence of $\lambda$ and then the difficulties of the extreme axion case are easily avoided Colpi et al. (1986); Urena-Lopez (2002); Liebling and Palenzuela (2012); Herdeiro et al. (2016); Urena-Lopez et al. (2012); Li et al. (2014); Rindler-Daller and Shapiro (2012, 2014); Diez-Tejedor et al. (2014). This requires, at least formally, that $(f/m_{\rm Pl})^{2}<0$ and then the trigonometric potential 1 is replaced by the hyperbolic one studied in Refs. Matos and Urena-Lopez (2000, 2001); Sahni and Wang (2000). The study of such case is part of ongoing work that will be presented elsewhere.

Acknowledgements

We would like to thank Alberto Diez-Tejedor for useful discussions and comments. FXLC acknowledges CONACYT for financial support. AXG-M acknowledges support from Cátedras CONACYT and UCMEXUS-CONACYT collaborative project funding. This work was partially supported by Programa para el Desarrollo Profesional Docente; Dirección de Apoyo a la Investigación y al Posgrado, Universidad de Guanajuato, research Grants No. 732/2017 y 878/2017; CONACyT México under Grants No. 167335, No. 179881, No. 269652 and Fronteras 281; and the Fundación Marcos Moshinsky.