---
author:
- 'L. Arturo Ureña-López,'
bibliography:
- 'references.bib'
title: 'Scalar field dark matter with a cosh potential, revisited'
---

Introduction {#sec:intro}
============

There is no doubt that one of the most fascinating riddles of modern cosmology is the dark matter (DM) that seems to be an ubiquitous component in the universe, and specially one that is indispensable for the formation of cosmological structure. DM is an essential part of the successful model Lambda Cold Dark Matter ($\Lambda$CDM), which has become the standard paradigm to understand the cosmos and its evolution up to its present state. In this model, DM is simply described by collisionless particles that interacts mostly gravitationally with other matter components, and makes up about $26\%$ of the total matter budget. Quite amazingly, the theoretical predictions of the $\Lambda$CDM model agree well with a wide range of cosmological and astrophysical observations [@PlanckCollaboration2018; @Bull2016; @Bertone2016]. However, the physical properties of the DM component remain as evasive as ever, mostly because the particle of the standard model that describes it as Weakly Interacting Massive Particle (WIMP), seems to be absent from most of the so-called direct detection experiments (eg [@Akerib2014; @Aprile2017; @Akerib2018; @Baudis2018]).

Given the crisis of the WIMP hypothesis, a major trend in modern studies about DM is characterized by the ’no stone left unturned’ approach [@Bertone2018], which asks for a thorough search of different alternatives to the standard CDM model. One possibility that has shown a rich phenomenology is the so-called Scalar Field Dark Matter (SFDM) model, which has been studied relentlessly for almost two decades now (under different names: fuzzy dark matter, wave dark matter, ultra-light axion particles, etc), see [@Hui2017; @Magana2012; @Marsh2016; @Lee:2017qve]. The common characteristics in all these variants are, firstly, the presence of a scalar field (SF), whether complex or real, which is endowed with a potential that contains, explicitly or implicitly, a mass term of the form $m^2_a \phi^2$, and secondly the (bare) mass of the SF is very light, of the order of $m_a \sim 10^{-22} \, \mathrm{eV}/c^2$[^1]. The foregoing properties are enough for the rich phenomenology we mentioned before, that allows the comparison of the model with a wide range of cosmological and astrophysical data, including, among others, gravitational waves, black holes, $21$-cm constraints, etc, see [@Amendola2006; @Hlozek2015a; @Urena-Lopez2016; @Ikeda2019; @Barack2018; @Brito2017; @Nebrin2018; @Sarkar2016; @Baumann2018] for some selected examples.

The properties of SFDM can be extended if one includes the presence of a self-interacting term of the fourth order, in the form $V(\phi) =
(1/2)m^2_a \phi^2 + (g_4/4) \phi^4$, where $g_4 > 0$ is a dimensionless constant. Self-interacting SF models have also been studied and their signatures as DM model have been widely discussed in [@Fan2016; @Li2014; @RINDLER-DALLER2014; @Goodman2000]. One can even consider the inclusion of higher order terms in the SF potential. The most famous case is the axion-like potential, $$V(\phi) = m^2_a f^2_a \left[ 1-\cos\left( \phi/f_a \right) \right] = \frac{m^2_a}{2} \phi^2 - \frac{m^2_a}{24 f^2_a} \phi^4 + \ldots \, , \label{eq:00}$$ where $f_a$ is called, for historical reasons, the axion decay constant. This type of models are known as axion-like particles (ALP) [@Marsh2016a], and the anharmonic nature of the ALP potential  produces observable signatures. The more noticeable effect appears for the evolution of linear density perturbations: there is an overgrowth of the density contrast with respect to the CDM case, which is characterized as a bump in the mass power spectrum (MPS) [@Zhang2017c; @Zhang2017b; @Cedeno2017].

As for the non-linear formation of structure, one has to consider to the non-relativistic limit of the Einstein-Klein-Gordon (EKG) system, which is the so-called Schrodinger-Poisson (SP) system [@Ruffini1969; @Seidel1990; @Guzman2004a; @Guzman2006]. There was early evidence about the similarities between the CDM model and the solutions of the SP system [@Widrow1993; @Woo2009], but it was until the results in [@Schive2014c], see also [@Schive2016; @Mocz2017; @Amin2019; @Li2018], that there appeared a clear and separate picture for the formation of structure under the SFDM hypothesis, specially for the differences with respect to CDM at small scales.

The gravitationally bounded objects that one could identify as DM galaxy halos all have a common structure: one central soliton surrounded by a Navarro-Frenk-White-like envelope created by the interference of the Schrodinger wave functions[@Schive2014c], features that have been confirmed by dedicated numerical simulations [@Schwabe2016; @Veltmaat2016; @Veltmaat2018; @Du2017; @Du2018].

The presence of a central soliton in all of SFDM galaxy halos has motivated studies about galactic kinematics to infer, first, the presence of such soliton structure, and, second, to determine the mass scale of the underlying SF particle [@Schive2014c; @Marsh2015a; @Gonzalez-Morales2012; @Gonzalez-Morales2017; @Bernal2018; @DeMartino2018; @Calabrese2016; @Marsh2018; @Lora2015; @Lora2012; @Urena-Lopez2017; @Robles2015; @Robles2018; @Bar2018a; @Bar:2019bqz; @Broadhurst2019a; @Lee2019]. The results are not yet conclusive, and depending on the analysis one may argue for positive presence of a soliton object and a SF mass of around $m_{a 22} \simeq 1$, or just an upper bound for the latter, $m_{a 22} < 0.4$. More recently, the presence of a soliton structure has been tested using data from rotation curves in galaxies, and then it is inferred that $m_{a22} > 10$ [@Bar2018a; @Bar:2019bqz].

The foregoing results on the SF mass, given that small scale structure seems to require light masses, are in tension, to say it mildly, with Lyman-$\alpha$ observations, which in contrast seem to demand larger values of the SF mass [@Amendola2006; @Irsic2017; @Armengaud2017]. According to the latter, the lower bound on the SF mass is $m_{a 22} > 21.3$, but this result is obtained from $N$-body simulations that cannot yet capture the whole properties of the SFDM model [@Nori2018AX-GADGET:Models; @Nori2018a], see for instance [@Zhang2018; @Zhang2017; @Zhang2017a; @Zhang2017d; @Li2018] for critical comments.

Given the motivations above, and as an added contribution to the studies of SFDM, the main aim in this paper is to revisit the properties of the hyperbolic counterpart of the axion-like potential , the cosh potential that was first studied in [@Sahni2000; @Matos2000; @Matos2001; @Matos2004], $$V(\phi) = m^2_a f^2_a \left[\cosh\left( \phi/f_a \right) - 1\right] = \frac{m^2_a}{2} \phi^2 + \frac{m^2_a}{24 f^2_a} \phi^4 + \ldots\, . \label{eq:0}$$

Although we have made an expansion of the cosh potential  similarly to that of the axion-like one , there are differences that go beyond the fourth order. Firstly, the cosh potential resembles the exponential one $V(\phi) \simeq (m^2_a f^2/2) e^{\pm \phi/f_a}$ for $|\phi|/f_a \gg 1$, whereas it becomes the standard free one $V(\phi) \simeq (m^2_a/2) \phi^2$ for $|\phi|/f_a \ll 1$. The latter is the desired form for the cosh potential to work as CDM at late times, whereas the former is an additional advantage of the cosh potential that has been used before to avoid any fine tuning of the initial conditions within a cosmological setting [@Ferreira1998; @Ferreira1997; @Copeland1998; @Sahni2000; @Matos2001; @Matos:2009hf]. However, the full properties of the cosh potential  were left aside and remain quite understudied (although see [@Beyer2014] for a similar model), and we intend to cover this gap.

A summary of the paper is as follows. In Sec. \[sec:cosmological-se\] we revisit the cosmological solutions, in an expanding universe, of the cosh potential , for both the background and linear perturbations of the SF. In doing that we will use the same numerical approach as in [@Urena-Lopez2016; @Cedeno2017], based on an amended version of the Boltzmann code CLASS (Cosmic Linear Anisotropic Solving System, [@Lesgourgues:2011re; @Blas2011; @Lesgourgues2011b; @Lesgourgues2011]). We will use the numerical solutions to uncover the differences between the cosh potential and other SFDM models, in particular regarding the early evolution of the SF density and the final form of the MPS.

In Sec. \[sec:self-grav\], we study the properties of the self-gravitating objects that can be formed under the cosh potential , but only in the non-relativistic limit and under the quartic approximation. These approximations are justified, as our purpose is to use the numerical solutions to obtain models for small galaxy halos. We shall show that, because of the cosmological constraints, the self-gravitating objects show a different scaling symmetry with respect to the free case, and that its positive self-interaction, as opposed to the negative one for the axion-like potential , has some advantages for their formation and gravitational stability. Finally, we discuss the main results and conclusions in Sec. \[sec:discussion-and\].

Cosmological setup \[sec:cosmological-se\]
==========================================

Here we present the most relevant equations of motion of SFDM endowed with the cosh potential , for both the background and linear perturbations. As mentiones before, we deal with the stage of rapid oscillations following the work in [@Urena-Lopez2016; @Cedeno2017], which is based on the polar transformation of the Klein-Gordon equation originally presented in [@Urena-Lopez2014; @Reyes-Ibarra:2010jje] for inflationary models (see also [@Roy2018; @Roy:2013wqa] for applications to quintessence models of dark energy).

Background Dynamics \[sec:background-dy\]
-----------------------------------------

The equation of motion for a scalar field $\phi$ endowed with the potential , in a homogeneous and isotropic space-time with null spatial curvature, is $$\ddot{\phi} = -3 H \dot{\phi} - m_a^2 f_a \sinh(\phi/f_a)  \, , \label{eq:kge}$$ where a dot denotes derivative with respect to cosmic time $t$, and $H
= \dot{a}/a$ is the Hubble parameter with $a(t)$ the scale factor. The scalar field energy density and pressure are given, respectively, by the canonical expressions: $\rho_\phi = (1/2)\dot{\phi}^2+V(\phi)$ and $p_\phi = (1/2)\dot{\phi}^2-V(\phi)$, whereas other matter components are the same as in the standard cosmological model.

We define a new set of variables adapted specifically for the cosh potential,

\[eq:newvars\] $$\begin{aligned}
\frac{\kappa \dot{\phi}}{\sqrt{6} H} &\equiv&  \Omega^{1/2}_\phi \sin(\theta/2) \, , \label{eq:newvarsa} \\
\frac{\kappa V^{1/2}}{\sqrt{3} H} &\equiv& \Omega^{1/2}_\phi \cos(\theta/2) = -\sqrt{\frac{2}{3}} \frac{\kappa m_\phi f_\phi}{H} \sinh(\phi/2f_\phi) \, , \label{eq:newvarsb} \\
y_1 &\equiv& -2 \sqrt{2} \frac{\partial_\phi V^{1/2}}{H} = -2 \frac{m_\phi}{H} \cosh(\phi/2f_\phi)\, , \label{eq:newvarsc}\end{aligned}$$

with $\Omega_\phi = \kappa^2 \rho_\phi/3H^2$ the standard SF density parameter, and $\kappa^2 = 8\pi G$. Similarly to the trigonometric case [@Cedeno2017], there is a constraint equation that arises from the known identity $\cosh^2 x - \sinh^2 x = 1$, which in terms of the new variables  reads $$y^2_1 = 4 \frac{m^2_a}{H^2} - \lambda \Omega_\phi \left( 1+\cos \theta \right) \, , \label{eq:constraint}$$ where for convenience we have defined the parameter $\lambda =
-3/\kappa^2 f^2$. Being an intrinsic mathematical property of the cosh potential , Eq.  should be satisfied at all times during the cosmic evolution.

After some straightforward algebra, the KG equation  takes the form of the following dynamical system:

\[eq:new4\] $$\begin{aligned}
  \theta^\prime &=& -3 \sin \theta + y_1 \, , \label{eq:new4a} \\
    y^\prime_1 &=& \frac{3}{2}\left( 1 + w_{tot} \right) y_1 + \frac{\lambda}{2} \Omega_\phi \sin \theta \, , \label{eq:new4b} \\
  \Omega^\prime_\phi &=& 3 (w_{tot} + \cos\theta)
  \Omega_\phi \label{eq:new4c} \, ,\end{aligned}$$

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 (EoS) is $w_{tot} = p_{tot}/\rho_{tot}$. For the SF itself, the EoS is explicitly given by $w_\phi = - \cos \theta$ in Eq. .

Notice that the dynamical system  is of general applicability for the three most known examples of SFDM. In fact, for $\lambda=0$, the dynamical system becomes that of the free case [@Urena-Lopez2016; @Hlozek2015a; @Amendola2006], whereas for $\lambda > 0$ (which formally requires $f^2_a < 0$ so that we can change from hyperbolic to trigonometric functions) becomes that of the axion (trigonometric) case [@Cedeno2017; @Zhang2017c; @Zhang2017b]

One critical step in the numerical solution of the equations of motion  is to find the correct initial conditions of the dynamical variables. For that, we consider that well within radiation domination the SF amplitude is such that $|\phi|/f \gg 1$, and then the potential  is approximated by an exponential one. Under this assumption, the field rapidly approaches the known scaling solution of exponential potentials [@Ferreira1997; @Ferreira1998; @Copeland1998] under which the SF density evolves like the dominant component. For the present case, the SF must keep, initially, a constant ratio with respect to the radiation density, $\rho_\phi/\rho_r = \mathrm{const.}$, and its equation of state must then be $w_\phi = -1/3$. Hence, the natural expressions for the initial conditions are, \[eq:match\] $$\Omega_{\phi i} = -12/\lambda \, , \quad \cos\theta_i = -1/3 \, . \label{eq:scaling}$$

Eq.  fixes the initial values of two variables only, which leaves $y_{1i}$ as the available parameter to adjust the final contribution of $\Omega_\phi$ as the DM component. This requires to match the early exponential behavior of potential  to its quadratic one at late times. The precise matching is explained in [@Matos2001] (see their Eq. 22), and the resultant relation between the SF mass $m_a$ and $\lambda$ is written as, $$m_{a22} \simeq 2.57 \times 10^{-10} h \, \mathrm{eV} \, \left( \frac{\lambda}{12} \right)^2 \left( 1- \frac{12}{\lambda} \right)^{3/2} \frac{\Omega^2_{\phi 0}}{\Omega^{3/2}_{r0}} \, , \label{eq:matcha}$$ where $h$ is the reduced Hubble parameter, and $\Omega_{r0}$ ($\Omega_{\phi 0}$) are the present values of the density parameter of radiation (SF).

From here we can calculate an expression that can be useful for the estimation of the initial conditions, in terms of the initial value of the SF mass to Hubble ratio. After some straightforward algebra, we obtain from Eq.  that

$$\frac{m_a}{H_i} = 1.5 \left[ \left( \frac{\lambda}{3} - 4 \right) \frac{\Omega_{\phi0}}{\Omega_{r0}} a_i \right]^2 \, . \label{eq:matchb}$$

In consequence, the initial condition required for variable $y_1$ can be written, after combining Eqs.  and  together, as $$y_{1i} = \sqrt{8} \left( 1 + \frac{1}{2} \frac{m^2_a}{H^2_i} \right)^{1/2} \, . \label{eq:initialb}$$

In practice, we calculate the initial mass to Hubble ratio in the form $m_a/H_i = A \times$\[Eq. \], where $A$ is a numerical coefficient that is adjusted using a shooting method already within CLASS to obtain the desired value $\Omega_{\phi 0}$ at the present time. The system converges rapidly after a few iterations, and the resultant SF mass to Hubble ratio is substituted into Eq.  to complete the initial conditions.

In the left panel of Fig. \[fig:rho-omega\] we show the evolution of the SF density $\rho_\phi$, for different values of the field mass $m_\phi$, in comparison with that of the standard CDM and radiation components of the fiducial $\Lambda$CDM model in [@PlanckCollaboration2018]. It can be clearly seen that initially the SF density redshifts as radiation, and eventually matches the CDM line during the stage of rapid oscillations. The different curves are labeled in terms of the SF mass, to ease the comparison with other SFDM models. In Table \[tab:values\] we report the values of $\lambda$ and the corresponding values of the SF mass $m_a$; the latter were read out from the numerical solutions, but it can be verified that they agree very well with those estimated directly from Eq. .

       $\lambda$           $m_a$ \[eV\]           $\Lambda_g$      $M_{\rm max}$ \[$M_\odot$\]
  ------------------- ----------------------- ------------------- -----------------------------
   $1.2 \times 10^2$   $1.1 \times 10^{-25}$   $1.0 \times 10^1$      $7.9 \times 10^{14}$
   $1.2 \times 10^3$   $1.3 \times 10^{-23}$   $1.0 \times 10^2$      $2.2 \times 10^{13}$
   $3.6 \times 10^3$   $1.2 \times 10^{-22}$   $3.0 \times 10^2$      $4.1 \times 10^{12}$
   $1.2 \times 10^4$   $1.4 \times 10^{-21}$   $1.0 \times 10^3$      $6.8 \times 10^{11}$
   $1.2 \times 10^5$   $1.4 \times 10^{-19}$   $1.0 \times 10^4$      $2.1 \times 10^{10}$

  : Different values of the free parameter $\lambda$ and their corresponding mass values $m_a$. The latter are closed to those obtained from Eq.  for the fiducial values $\Omega_{\phi 0} = 0.26$, $\Omega_{r0} = 9.13 \times 10^{-5}$, and $h=0.67$. The self-interaction strength $\Lambda_g$ and the maximum mass $M_{\rm max}$ indicated here refer to those of self-gravitation objects in Eqs.  and , respectively. \[tab:values\]

We see that the matching of SFDM to CDM happens later (earlier) for smaller (larger) values of the SF mass, which is a common feature of SFDM models. For a quick comparison, we also plotted the results from the free case (dashed lines), which shows that both models start their rapid oscillations at the same time. This is not surprising, as the evolution of the SFDM with a cosh potential can be parametrized in terms of the SF mass only (recall Eq. ), so that the appearance of the rapid oscillations is solely dictated by the mass to Hubble ratio $m_a/H$, and then its late behavior should match that of the free model too.

As mentioned before, the SF density has a non-negligible contribution during radiation domination, but we note that such contribution diminishes for larger values of the SF mass. For a better illustration of this effect, we show the ratio of the different matter components with respect to the critical density in the right panel of Fig. \[fig:rho-omega\] (see also [@Matos2000; @Matos2001]), and make a comparison again with the results from $\Lambda$CDM. In this form, the contribution of the SF component is noticeable during radiation domination. Also, it must be noticed that the radiation, matter, $\Lambda$ domination eras do not suffer any alteration and proceed just like in the $\Lambda$CDM case.

![image](rhos){width="0.5\linewidth"} ![image](omegas){width="0.5\linewidth"}

One common concern about the presence of the SF at early times is its contribution as an extra relativistic degree of freedom, that increases the expansion rate at the time of nucleosynthesis. The relativistic degrees of freedom are given by $N_{\rm eff} = \rho_{\rm rad}/\rho_\nu$, where now $\rho_{\rm rad}$ includes the contribution from the SF density at early times, and $\rho_\nu$ is the neutrino density. The resultant effect is shown in Fig. \[fig:nucleo\], for the indicated values of the SF mass and in terms of the quantity $\Delta N_{\rm eff} \equiv N_{\rm eff} -3.046$.

The non-negligible contribution of the SF to $\Delta N_{\rm eff}$ can be clearly seen during radiation domination and at the time around nucleosynthesis (see also Fig. \[fig:rho-omega\]). The present constraints, from [@PlanckCollaboration2018], then indicate that $\Omega_\phi < 3.5 \times 10^{-2}$, which in turn translates, by means of Eq. , into $|\lambda| > 3.4 \times 10^2$. Considering the values shown in Table \[tab:values\], the cosh potential  then evades the nucleosynthesis constraints as long as its field mass is $m_{a22} > 10^{-2} \mathrm{eV}$.

![The extra relativistic degrees of freedom $\Delta N_{\rm eff}$ induced by the non-negligible contribution of the SF density at around the time of nucleosynthesis. The shaded region represents the reported constraint in [@PlanckCollaboration2018], $\Delta N_{\rm eff} < 0.3$ at $95\%$CL. \[fig:nucleo\]](nucleo){width="0.5\linewidth"}

As a final note, we contrast the nucleosynthesis constraint on the cosh potential to that applied on a SFDM with a quartic self-interaction (although with a complex SF) [@Li2014]: for the latter there is a kinetic dominate phase (which manifests as a stiff-fluid behavior) prior to the radiation-like behavior of the SF. Thus, one has to finely tune the values of the SF mass and the quartic interacting parameter to avoid any undesirable effects upon the nucleosynthesis process[^2]. The cosh potential  evades such tight constraints because the early radiation-like solution  is a quite stable attractor solution of the equations of motion .

Linear Density Perturbations \[sec:linear-den\]
-----------------------------------------------

We now turn our attention to the linear field perturbations $\varphi$ around the background value in the form $\phi(x,t) = \phi(t)+\varphi(x,t)$. As in previous works, we choose the synchronous gauge with the line element $ds^2 = -dt^2+a^2(t)(\delta_{ij}+\bar{h}_{ij})dx^idx^j$, where $\bar{h}_{ij}$ is the tensor of metric perturbations. The equation of motion for a given Fourier mode $\varphi(k,t)$ reads [@Ratra1988; @Ferreira1997; @Ferreira1998; @Perrota1999] $$\ddot{\varphi} = - 3H \dot{\varphi} - \left[\frac{k^2}{a^2} +m_a^2\cosh(\phi/f)\right] \varphi -
  \frac{1}{2} \dot{\phi} \dot{\bar{h}} \, , \label{eq:13}$$ where a dot means derivative with respect the cosmic time, $\bar{h} = {\bar{h}^j}_j$ and $k$ is a comoving wavenumber.

Following [@Cedeno2017; @Urena-Lopez2016], we can choose appropriate variables to transform Eq.  into the following dynamical system,

\[eq:dp\] $$\begin{aligned}
\delta^\prime_0 &=&  \left[-3\sin\theta-\frac{k^2}{k^2_J}(1 - \cos \theta) \right] \delta_1 + \frac{k^2}{k^2_J} \sin \theta \, \delta_0 - \frac{\bar{h}^\prime}{2}(1-\cos\theta) \, , \label{eq:dp0} \\
\delta^\prime_1 &=& \left[-3\cos \theta - \left( \frac{k^2}{k^2_J} - \frac{\lambda \Omega_\phi}{2y_1} \right) \sin\theta \right] \delta_1 + \left(\frac{k^2}{k^2_J} - \frac{\lambda \Omega_\phi}{2y_1} \right) \left(1 + \cos \theta \right) \, \delta_0 - \frac{\bar{h}^\prime}{2} \sin \theta \, , \label{eq:dp1}\end{aligned}$$

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$. In this approach $\delta_0 \equiv \delta \rho_\phi/\rho_\phi$ is the scalar field density contrast and $\delta_1$ is a second density constrast that arises naturally under the new variables. As before for Eqs. , Eqs.  have again the same structure for the cosh, axion and free potentials, and the only difference is the value of the interaction parameter $\lambda$ for each case.

For the initial conditions, we use the same solutions at early times considered for a quadratic potential (which allows us to use one single Boltzmann code for three different potentials). For the particular case of the exponential potential the true attractor solution is $\delta_0 = (4/15) \delta_{CDM}$[@Matos2001; @Ferreira1998], and such solution is quickly reached just after some $e$-folds (see Fig. \[fig:deltas\] below).

The numerical results illustrating both the cosh and quadratic cases, for $m_{a22} = 2.3 \times 10^{-2}$ and $\lambda = 5\times 10^2$, are shown in Fig. \[fig:deltas\]. The plots replicate those presented in [@Matos2000; @Matos2001], which were the first to show the evolution of density perturbations for the cosh potential obtained directly from a Boltzmann code (the now outdated CMBFAST).

In the case of large enough scales, we see that the SF density contrast reaches the attractor solution at early times mentioned above, but quickly joins the standard CDM evolution once the SF oscillations start. In contrast, for small enough scales, again the density contrast reaches the same attractor solution at early times but it ends up oscillating at late times with almost a constant amplitude. The results from the cosh potential match those of the quadratic one at late times, which is a common feature in our numerical results: as long as the SF mass is the same, the two models are almost indistinguishable.

![image](deltas){width="0.5\linewidth"} ![image](deltas1){width="0.5\linewidth"}

Finally, we show in Fig. \[fig:mps\] the MPS at the present time for a cosh potential, also in comparison with the results that are obtained for the quadratic one with the same values of the SF mass. As expected from the discussion above about the density contrast, the MPS is practically indistinguishable in the two cases, which means that the well-known cut-off in the MPS of the cosh potential (first reported in [@Matos2001] for any SFDM model) just depends upon the values of the SF mass. The cosh potential is then an example in which a more involved SF potential does not provide, within the cosmological context, different results from those of the free case. The reason for this is, of course, the close relation between the SF mass $m_a$ and the self-interaction parameter $\lambda$ in Eq. , which is required to obtain an appropriate radiation to matter transition for a corrrect contribution of SFDM at late times.

![MPS at the present time for SFDM with a cosh potential, for the values of the SF mass and $\lambda$ as indicated in the plot. Also shown is the quadratic potential with the same values of the SF mass (dashed lines in colour). The MPS appears to be the same for both models, see the text for more details. \[fig:mps\]](mpsk){width="0.5\linewidth"}

Self-gravitating objects \[sec:self-grav\]
==========================================

A separate question of physical interest is the type of self-gravitating objects that are formed under a cosh potential. This topic has been treated before in [@Urena-Lopez2002b], and here we just give a brief account of the main results reported throughout the specialized literature (see [@Liebling2017; @Schunck2003GeneralStars] for recent reviews).

The key parameter here is the quartic coupling that arises from the series expansion of potential  up to the fourth order: $V(\phi) \simeq (m^2_a/2) \phi^2 + (m^2_a f^{-2}_a/24) \phi^4$. Following standard nomenclature [@Liebling2017], although adapted for real scalar fields (see Eqs.  below), the self-interaction strength that is of physical interest is quantified by the combination[^3] $$\Lambda_g \equiv m^2_a f^{-2}_a/(32 \pi G m^2_a) = 3/(12\kappa^2 f^2_a) = - \lambda/12 \, , \label{eq:quartic}$$ and then determined solely by the interaction parameter $f_a$ in the potential . In passing by, notice also that the initial condition  can also be written as $\Omega_{\phi i} = \Lambda^{-1}_g$.

For the values of $\lambda$ suggested by the cosmological constraints discussed above, we find that the self-interaction strength must be of the order of $\Lambda_g \sim 10^2$ or larger. In any case, we find that cosmological constraints seem to suggest that self-gravitating objects belong to the so-called strong-regime, $\Lambda_g \gg 1$. In what follows we will revise the properties of such objects.

General stability and equilibrium configurations \[sec:general-stability\]
--------------------------------------------------------------------------

To begin with, there is a maximum mass for stable configurations (see [@Colpi1986; @Balakrishna1998; @Urena-Lopez2002b]), that is approximately given by $M_{\rm max} \simeq 0.22 \Lambda^{1/2}_g
m^2_{\rm Pl}/m_\phi$, and which can be a very large number for the field masses of cosmological interest[^4]. Moreover, as we have seen in Sec. \[sec:cosmological-se\] above, for the particular case of the cosh potential  there exists a close relation between the SF mass and the self-interaction strength, see Eq. , and then the maximum mass can be written in terms of $\Lambda_g$ alone. The two aforementioned relations, for the fiducial cosmological values assumed above, are

\[eq:masses\] $$\begin{aligned}
    m_{a22} &\simeq& 1.37 \times 10^{-5} \, \Lambda^2_g \left( 1- 1/\Lambda_g \right)^{3/2} \, , \label{eq:massesa} \\ 
    M_{\rm max} &\simeq& 2.15 \times 10^{16} \, M_\odot \, \left( \Lambda_g - 1 \right)^{-3/2} \, . \label{eq:massesb}\end{aligned}$$

Some examples of $M_{\rm max}$ are shown in Table \[tab:values\] for reference. Notice that the maximum mass could be as low as $10^{10} \, M_\odot$, which means that, for large enough values of the interaction parameter $\Lambda_g$, some configurations of astrophysical interest could be close to the point of gravitational instability. In the least troublesome case, unstable configuration just migrate to stable ones by means of gravitational cooling, but in other more involved cases they could collapse into black holes [@Seidel1991; @Balakrishna1998; @Alcubierre2003Numerical2-oscillatons; @Guzman2006; @Urena-Lopez2012; @Liebling2017].

However, the case of more astrophysical interest is the non-relativistic one which is obtained in the weak field limit of the Einstein-Klein-Gordon (EKG) system. Under the ansatz $\sqrt{8\pi G} \phi = e^{-im_at} \psi + \mathrm{c.c.}$, where $\psi$ is a non-relativistic, complex wave function, the EKG system can be recast as the Gross-Pitaevskii-Poisson (GPP) system of equations (eg [@Chavanis2011]). Using the field mass to define dimensionless time and distance, $t \to t/m_a$ and $\mathbf{x} \to \mathbf{x}/m_a$, the final form of the GPP equations is (for details see [@Guzman2004a; @Lee1996; @Chavanis2011]),

\[eq:GPP\] $$\begin{aligned}
 i \dot{\psi} &=& - \frac{1}{2} \nabla^2 \psi + \left( U + \Lambda_g |\psi|^2 \right) \psi \, , \label{eq:GPPa} \\
 \nabla^2 U &=& |\psi|^2 \, , \label{eq:GPPb}\end{aligned}$$

where $U$ is the Newtonian gravitational potential. If $\Lambda_g =0$, we recover from Eqs.  the well-known Schrodinger-Poisson (SP) system that has been the standard workhorse to study the formation of large scale structure in SFDM models (eg [@Schive2014; @Hwang2009; @Veltmaat2016; @Veltmaat2018; @Mocz2017]).

The GPP system  is invariant under the following scaling symmetry, $$\left\{ t,\mathbf{x}, \psi,U, \Lambda_g \right\} \; \to \; \left\{ \alpha^{-2} \hat{t}, \alpha^{-1} \hat{\mathbf{x}}, \alpha^2 \hat{\psi}, \alpha^2 \hat{U}, \alpha^{-2} \hat{\Lambda}_g \right\} \, , \label{eq:scaling-sym}$$ where $\alpha$ is an arbitrary parameter. This somehow eases the numerical effort to find the equilibrium configurations of the GPP system if we can choose an appropriate value for the scaling parameter $\alpha$. The usual choice in the non-interacting case ($\Lambda_g =0$) is $\alpha = |\psi_c|^{-1/2}$, so that the rescaled wavefunction can simply have a central value of order unity, $\hat{\psi}_c =1$ [@Guzman2004a]. This cannot be done anymore in the interacting case, and then the most sensible choice is to use instead the strength parameter to unity so that $\alpha = \Lambda^{-1/2}_g$. After this, there is not any further scaling available for Eqs.  to ease the numerical analysis.

Of particular interest here are the so-called equilibrium configurations of the GPP system, which are stationary, spherically symmetric, solutions of the wave function in the form $\psi(t,r) = \varphi(r) e^{-i\sigma t}$, where $\sigma$ is a constant parameter. Upon substitution of the foregoing ansatz in Eqs. , the resultant system of equations is

\[eq:GPP-stationary\] $$\begin{aligned}
 \partial_{rr} (r \varphi) &=& 2 \left( U - \sigma + \varphi^2 \right) \varphi \, , \label{eq:GPP-stationarya} \\
 \partial_{rr} (r U) &=& \varphi^2 \, . \label{eq:GPP-stationaryb}\end{aligned}$$

Notice that we have made use of the scaling symmetry  and then the re-scaled self-interaction parameter has been set to unity, $\hat{\Lambda}_g =1$. For purposes of simplicity in the notation, all quantities in Eqs.  are assumed to have been re-scaled according to . Apart from the numerical mass, other quantities were also calculated for each one of the solutions; their numerical values are shown in Table \[tab:equilibrium\].

   $\Lambda_g$   $\varphi_c$   $\sigma$   $M_T$   $r_{95}$
  ------------- ------------- ---------- ------- ----------
        0            1.0        -0.69     2.06      3.93
        1            0.1        -0.06     0.63     10.89
                     0.5        -0.46     2.23      5.83
                     1.0        -1.25     4.84      4.36
                     1.5        -2.47     8.78      3.78

  : Characteristic values of non-interacting and interacting equilibrium configurations. The quantities shown are the strength parameter $\Lambda_g$, the central field $\varphi_c$, the eigen-frequency $\sigma$, the total mass $M_T$, and the radius containing $95\%$ of the total mass $r_{95}$. \[tab:equilibrium\]

Weak and strong field limits
----------------------------

It is known that there is a universal density profile for non-interacting equilibrium configurations, which can be approximated by the formula: $\rho(r) =
\rho_c/(1+r^2/r^2_s)^8$ [@Marsh2015a; @Schive2014c], and whose central density $\rho_c$($= |\psi_c|^2$) and scale radius $r_s$ are related in the scale-independent form $\rho_c r^4_s = \mathrm{const.}$ (see Eqs. ). The analytic density profile has been widely used in diverse studies of SFDM and in their comparison with diverse astrophysical data [@Gonzalez-Morales2017; @Herrera-Martin2019; @Urena-Lopez2017; @Marsh2015a; @Bernal2018].

Things are not quite straightforward in the interacting case, as one cannot find, in general, a universal density profile that is just rescaled for all possible configurations. However, if the central field value is weak enough, $\hat{\varphi}_c \ll 1$, the self-interaction term in Eq.  must become negligible, and in this limit we should somehow recover the non-interacting case. This is exactly what we found in our numerical solutions, which is also shown in the left panel of Fig. \[fig:equilibriuma\] for the relation between the total mass $\hat{M}_T$ and the $95\%$ radius $\hat{r}_{95}$. Weaker configurations have larger radius, and it is at larger radii that we see the equivalence between the interacting and non-interacting cases; in terms of the central field, the equivalence is reached for $\hat{\varphi}_c \lesssim 0.1$.

![image](mass-radius){width="0.5\linewidth"} ![image](densities){width="0.5\linewidth"}

The main advantage here is that we can use, for the interacting case, all the known scaling symmetries of the SP system, avoiding all the hassle of calculating configurations case by case. Even more, despite the presence of the self-interacting term $\Lambda_g$, we can choose a different scaling parameter in the form $\alpha \to \alpha \, \Lambda^{1/2}_g$, with which the numerical solutions would rightly take the same form as in the non-interacting case.

On the other hand, there is the belief that interacting configurations have an analytic expression for their density in the form $\rho(r) = \rho_c \sin(\pi r/r_s)/(\pi r/r_s)$, for $r \leq r_s$ and zero otherwise. The scale radius of the configuration is given by $r_s = (\pi/\sqrt{2}) |\Lambda_g|^{1/2} m^{-1}_\phi$, whereas the central density would be a free parameter [@Diez-Tejedor2014b; @Chavanis2012; @Lee1996]. We have not found any evidence of such common profile in the numerical results of the GPP system . As shown in the right panel of Fig. \[fig:equilibriuma\], the normalized density profiles, $\rho(r)/\rho_c$, for some of the numerical solutions listed in Table \[tab:equilibrium\]. The profiles differ one from each other, as they correspond to different values of $\varphi_c$, and they all differ from the analytic profile presented before.

The belief in the analytic expression seems to have been originated from the Newtonian approximation used in [@Lee1996] of the relativistic equilibrium configurations obtained in [@Colpi1986] for the limit $\Lambda_g \to \infty$. It was assumed in [@Lee1996] that Newtonian limit of the fundamental frequency is $\omega/m \to 1$, and then the density must obey the equation $\nabla^2 \rho = - \rho$. However, the correct result is $\omega/m \to
1 - \sigma/2$, where $\sigma$ is the eigen-frequency in Eq. . Hence, the equilibrium configuration must be obtained again from an eigenvalue problem to determine the correct value of $\sigma$. Actually, one can see in Fig. 3 of [@Colpi1986] that the analytic profile is just an approximation and that the actual, numerically obtained, field profile deviates considerably from it at large values of $r$ [^5].

As a final note in the strong field limit, we must recall that there is a critical mass for equilibrium configurations, see Eq. , which also corresponds to a critical field value, estimated to be $\kappa_0 \phi \simeq
\Lambda^{-1/2}_g$ [@Urena-Lopez2002b; @Colpi1986]. For non-relativistic configurations this translates, in terms of re-scaled quantities, into the constraint $\hat{\varphi}_c \lesssim 1$. As shown in the right panel of Fig. \[fig:equilibriuma\], it may be possible to reach the analytic profile $\sin (\pi r)(\pi r)$ if $\hat{\varphi} \gg 1$, but this means that the resultant configuration is well within the unstable branch of relativistic equilibrium solutions (as shown in [@Guzman2013]), which for that very reason makes it of less interest for cosmological and astrophysical applications.

Discussion and conclusions \[sec:discussion-and\]
=================================================

We have revisited the general properties of one SFDM model with a cosh potential, which was proposed as an alternative model to CDM almost two decades ago. This is an example, similar to the axion case, in which one can include non-linear terms in the SF potential that open the door for distinguishable features with respect to the simple quadratic case, which has been used as the paradigm for SFDM models.

At the cosmological level, we reviewed the solutions of the cosh model for physical quantities in the background and their linear perturbations. The distinctive feature of the cosh model is its scaling solution during the epoch of radiation domination, under which the SF has a non-negligible contribution to the radiation component. Such scaling solution is a strong attractor solution, which helps the SFDM model to evade the fine tuning problem, to get the right DM contribution at the present time, that is unavoidable for the quadratic and axion cases. By means of the nucleosynthesis constraint for extra, non-thermal, relativistic degrees of freedom, we found an upper bound on the interaction parameter, $\kappa f_a \lesssim 9.4 \times 10^{-2}$ ($f_a \lesssim 1.9 \times 10^{-2} m_{\rm Pl}$), which is equivalent to a lower bound on the SF mass, $m_{a22} > 10^{-2}$.

As for the analysis of the linear density perturbations, we confirmed the existence of an attractor solution for them at early times, which also alleviates the fine tuning problem of SF linear perturbations. The resultant MPS summarizes well the general properties of the density perturbations with the presence of a sharp cut-off of power at small scales. Although the cut-off had already been shown to exist in [@Matos2001], we obtained here that the cut-off in the cosh case is quite similar to that of the quadratic one, and then for both cases the properties of the density perturbations are characterized by the SF mass only.

We also revisited the properties of the self-gravitating objects obtained from the cosh potential, although only for the Newtonian limit and using a quartic approximation to the potential. The self-interaction term for the cosh potential is positive definite, which then contributes with an extra repulsive force to prevent the gravitational collapse of the SF equilibrium configurations. This reinforces their gravitational stability, as the maximum mass they can attain is larger than for the non-interacting case.

However, we made estimations for the SF amplitude required to represent realistic galaxies by means of soliton structures and found that, in such limit, the self-interaction term can be neglected. Hence, the soliton configurations from SFDM with a cosh potential can be safely represented by the standard free ones, whose known properties are closely linked to the SF mass only [@Schive2014c].

Another question that usually arises within the context of SF DM models is the possible formation of a Bose-Einstein condensate (BEC). Much attention has been put forward in the formation of a BEC for the axion case ever since the results reported in [@Sikivie2009]. However, such results have been heavily contested because of the attractive (negative) self-interaction attributed to axion-like potentials [@Guth2015]. Actually, previous works had already shown that BEC formation is possible for repulsive (positive) self-interactions by means of an inverse particle cascade that populates the zero-mode in the distribution function of particles [@PineiroOrioli2015], whereas the opposite happens, eg the depopulation of the zero-mode by a direct particle cascade, for attractive (negative) self-interactions. The cosh potential  has a repulsive self-interaction, and then the formation of a BEC is an expected result that would clearly differentiate it from axion-like potentials.

A SFDM model with a cosh potential is a candidate as good as the free model to be an alternative to CDM, with the added advantages of being a non-linear potential with a positive definite self-interaction. As is currently known for such kind of DM models, there is a good agreement with observations from large (CMB, MPS) to small scales (eg dwarf spheroidals), but the question remains about its appropriateness for the full range of scales. Of particular interest are constraints from Lyman-$\alpha$ observations, which have since the study in [@Amendola2006] been a point of major concern for the preferred model with $m_{a22} \sim 1$ [@Zhang2017b; @Zhang2017a; @Irsic2017; @Armengaud2017; @Nori2018a; @Li2018], but we may have to wait for the development of better numerical simulations before we can decide firmly about the viability of the model.

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; and CONACyT México under Grant No. 17899.

[^1]: Herafter, and for purposes of simplicity in the notation of the mathematical expressions, we will use natural units with $c=1=\hbar$.

[^2]: A similar study of the quartic self-interaction model for a real scalar field can be found in [@Cembranos2018], although the initial conditions were seemingly chosen so that the early stiff and radiation like behaviors were avoided. In such case, the nucleosynthesis constraints are not important, and both the mass and self-interaction parameters can be varied more freely than in the case considered in [@Li2014].

[^3]: Notice that the interaction parameter in the case of boson stars, as originally defined in [@Colpi1986], is $\Lambda_g \equiv g_4/(4\pi G m^2_a)$, where $g_4$ is the quartic coupling constant as in $(g_4/4)\phi^4$. In the case of real scalar fields, a more appropriate definition is $\Lambda_g \equiv 3 g_4/(8\pi G m^2_a)$ [@Urena-Lopez2002b; @Urena-Lopez2012]. Once we consider the Newtonian limit of the Klein-Gordon equation for real scalar fields, it appears convenient to further take $\Lambda_g \equiv 3 g_4/(16\pi G m^2_a)$, see [@Lee1996; @Guzman2004a; @Guzman2006]. The latter is the definition used to write $\Lambda_g$ in Eq. .

[^4]: The corresponding maximum mass in the non-interacting case ($\Lambda_g =0$) is $M_{\rm max} \simeq 0.6 m^2_{\rm Pl}/m_a$, see for instance [@Seidel1990; @Seidel1991; @Alcubierre2003Numerical2-oscillatons]. As for ALP models with the axion potential , that has a negative self-interaction strength $\Lambda_g < 0$, the estimated maximum mass is $M_{\rm max} \simeq |\Lambda_g|^{-1/2} m^2_{\rm Pl}/m_a$ [@Guzman2006; @Chavanis2016; @Liebling2017; @Helfer:2016ljl]; then, the larger the self-interaction the less massive are stable configurations.

[^5]: The assumed profile with a large self-interaction also appears from the so-called Thomas-Fermi (TF) approximation, in which one neglects the spatial derivatives of the wavefunction $\psi$. The TF profile has been shown to be at variance with galactic observations [@Arbey2003; @Bohmer2007; @Diez-Tejedor2014b]. It would be, for the reasons explained in the text, incorrect to use those negative results against SFDM models.
