# Analytical approach to core-halo structure of fuzzy dark matter

Atsushi Taruya<EMAIL_ADDRESS>Center for Gravitational Physics
and Quantum Information, Yukawa Institute for Theoretical Physics, Kyoto
University, Kyoto 606-8502, Japan  Kavli Institute for the Physics and
Mathematics of the Universe (WPI), The University of Tokyo Institutes for
Advanced Study, The University of Tokyo, 5-1-5 Kashiwanoha, Kashiwa, Chiba
277-8583, Japan  Shohei Saga<EMAIL_ADDRESS>Sorbonne Université, CNRS, UMR7095,
Institut d’Astrophysique de Paris, 98bis boulevard Arago, F-75014 Paris,
France  Laboratoire Univers et Théories, Observatoire de Paris, Université
PSL, CNRS, Université de Paris, 5 place Jules Janssen 92190 Meudon, France

###### Abstract

Ultralight bosonic dark matter called fuzzy dark matter (FDM) has attracted
much attention as an alternative to the cold dark matter. An intriguing
feature of the FDM model is the presence of a soliton core, a stable dense
core formed at the center of halos. In this paper, we analytically study the
dependence of the soliton core properties on the halo characteristics by
solving approximately the Schrödinger-Poisson equation. Focusing on the
ground-state eigenfunction, we derive a key expression for the soliton core
radius, from which we obtain the core-halo mass relations similar to those
found in the early numerical work, but involving the factor dependent
crucially on the halo concentration and cosmological parameters. Based on the
new relations, we find that for a given cosmology, (i) there exist a
theoretical bound on the radius and mass of soliton core for each halo mass
(ii) incorporating the concentration-halo mass (C-M) relation into the
predictions, the core-halo relations generally exhibit a non-power-law
behavior, and with the C-M relation suppressed at the low-mass scales,
relevant to the FDM model, predictions tend to match the simulations well
(iii) the scatter in the C-M relation produces a sizable dispersion in the
core-halo relations, and can explain the results obtained from cosmological
simulations. Finally, the validity of our analytical treatment are critically
examined. A perturbative estimation suggests that the prediction of the core-
halo relations is valid over a wide range of parameter space, and the impact
of the approximation invoked in the analytical calculations is small, although
it is not entirely negligible.

††preprint: YITP-22-82

##  I Introduction

The dark matter, invisible matter component, is an essential constituent
accounting for $\sim 30\%$ of the energy contents of the Universe. Although
the origin of dark matter is not yet clarified, the structure formation driven
by the gravitational instability, consistent with large-scale structure
observations, suggests that the dark matter is nonbaryonic, and have
negligible velocity dispersion, referred to as the cold dark matter (CDM).
Since its foundation in early 1980s [1, 2, 3], the CDM has been playing an
important role to establish a concordant picture of the Universe, and it
successfully describes a wide range of observations on scales larger than
galaxies. However, on the subgalactic scales, recent observations suggest that
deviation from CDM predictions is visible, and remains significant especially
for the observations of the local Universe, highlighted as the small-scale
crises (see e.g., [4, 5] for review).

To relieve possible flaws in the CDM paradigm, ultralight particle of masses,
$m_{\phi}\sim 10^{-22}$-$10^{-20}$eV, called fuzzy dark matter (FDM) [6], has
recently attracted much attention as an alternative to the CDM (Refs. [7, 8,
9, 10] for review, see also Ref. [11]). While the predictions for the
evolution and formation of large-scale structure remain essentially the same
as in CDM model, and thus consistent with large-scale observations (e.g., Ref.
[12, 13] for recent works), many prominent features appear at small scales due
to a very long de Broglie wavelength that extends over astrophysical scales,
$\lambda_{\rm dB}=h/(m_{\phi}v)\simeq 0.48\,{\rm kpc}\,(10^{-22}\,{\rm
eV}/m_{\phi})(250{\rm km}/v)$, with $v$ being the particle velocity. The FDM
model can therefore modify small-scale structure formation, and thus the
small-scale observations can place a tighter bound on the FDM model (e.g.,
Refs. [14, 15, 16] for constraints from the Lyman-$\alpha$ forest).
Importantly, modifications of the small-scale structures appears manifest not
only through the small-scale cutoff of the linear power spectrum [6] (see also
Refs. [17, 18] for early works), but also by the quantum nature of dark matter
particles at a macroscopic level. The latter aspect is particularly
interesting, and can bring unique signature of FDM models to observationally
test and constrain their model parameters.

Numerical studies on the FDM structure formation show that the dark matter
halos commonly host a dense flat core called soliton with a universal density
profile at the central region, accompanied by granular structures with a
typical size of the de Broglie wavelength that extends over the outskirt of
halos [19]. A notable feature of the soliton was its relation to the host
halos. In early works, the mass and radius of FDM solitons, $M_{\rm c}$ and
$r_{\rm c}$, are found to be correlated tightly with halo masses in a power-
law form as $M_{\rm c}\propto M_{\rm h}^{\alpha}$ and $r_{\rm c}\propto M_{\rm
h}^{\beta}$, and Ref. [20] obtained $\alpha\sim 1/3$ and $\beta\sim-1/3$.
Later, these properties have been examined with different setups by several
groups [21, 22, 23, 24], and obtained a similar result but with slightly
different slopes (see also Ref. [25] for systematic study based on a particle-
based method). More recent simulations covering a large cosmological volume
suggest that not tight scaling relation exists between soliton structure and
halo mass, and rather indicate a large scatter in the correlation [26, 27].
Another intriguing property of the soliton core comes from the results of
high-resolution simulations [28, 29], showing that the soliton core is not
strictly stable, but rather moves around a halo center in a random fashion,
accompanying a rapid oscillation with the order-of-unity variation of its
amplitude. These phenomena are considered as a consequence of the wave
interference [30, 31].

So far, most of these properties has been investigated by numerical
simulations of the Schrödinger-Poisson equation, starting with various setups,
and little analytical work has been done (but see [32, 33]). While the soliton
core-halo mass relation has led to a number of works that have tried to test
and constrain the FDM model (e.g., Refs. [34, 35, 36, 37, 38]), a major
difficulty in characterizing such a relation in numerical simulations is a
severe requirement of the dynamic range. That is, a realistic simulation must
describe, in a relevant cosmological setup, both the halo merger processes and
wave nature of soliton dynamics, typical scales of that, respectively, appear
above Mpc and below kpc, thus requiring both a large simulation box and a
sufficient spatial resolution. Even the state-of-the-art simulation is still
challenging to trace an entire history of the formation and evolution of
cosmic structures while fulfilling both the conditions. Aim of this paper is
therefore to investigate analytically the soliton core structure, and to
address the core-halo relations in a cosmological setup.

In doing so, a crucial step is how to solve the Schrödinger-Poisson equation.
Coupled with the Poisson equation sourced by the mass density of FDM, the
governing equation for the wave function of FDM becomes nonlinear. In
particular, if one wants to address the soliton core structure, the
contributions of both the soliton and surrounding halo structures have to be
taken into account consistently. As we will describe in detail, for an
analytically tractable calculation, our approach is to consider the
contribution of halo as a smooth background, and to describe the soliton state
ignoring its self gravity. Although the flat core structure is realized under
a balance between quantum pressure and gravity and one thus expects the self-
gravity of the soliton to play a certain role, this setup greatly simplifies
the situation, and the problem is reduced to a linear eigenvalue problem.
Still, there exists no exact solution, and some approximations need to be
invoked to construct analytically the eigenfunctions of the soliton state. In
this paper, we apply the uniform asymptotic approximation invented by Langer
[39, 40, 41] to the description of soliton state. In the presence of a large
system parameter, this approximation is known to significantly outperform the
familiar Wenzel-Kramers-Brillouin (WKB) approximation, and provides a very
accurate way to estimate both eigenvalues and eigenfunctions. In particular,
the analytically constructed eigenfunctions by this method is globally valid
even around the so-called turning points, where the WKB approximation is
broken down. Making use of the Langer’s method, in this paper, we present
analytical solutions of the Schrödinger-Poisson equation under a spherically
symmetric halo potential. Focusing particularly on the ground-state wave
function for the zero angular momentum $(\ell=0)$, we derive the analytical
expressions for the core-halo relations. The resultant expressions apparently
resemble those found in numerical simulations, but importantly, they include
an additional factor dependent on cosmology and halo properties, with which
the predicted core-halo relations can differ from a power-law form found in
small-box simulations. We will discuss it in detail in comparison with
simulation results. Further, we estimate the impact of the soliton self-
gravity which we ignore on the core-halo relations, and the validity of our
analytical treatment is discussed.

This paper is organized as follows. In Sec. II, we begin by describing the
setup of the problem to characterize the soliton core structure in a halo.
Section III then presents the analytical construction of the soliton
eigenstate, focusing on the cases with zero angular momentum. The construction
for the nonzero angular momentum is presented in Appendix B. Section IV
presents the main results of this paper, and based on the treatment in
previous section, we derive the analytical expressions of core-halo relations,
and provide a general prediction for the core radius vs halo mass and core
mass vs halo mass relations, treating the halo concentration as a free
parameter. Further, imposing models of the concentration-halo mass relation
(C-M relation), the predicted relations are compared with simulations. In Sec.
V, the validity of our treatment ignoring the soliton self-gravity is
critically examined, and employing a perturbative calculation, we present an
estimate of its impact on the soliton core radius. Finally, Sec. VI is devoted
to summary of important findings and conclusion.

Throughout the paper, we work with the natural units, setting the speed of
light $c$ and Planck constant $\hbar$ to unity.

##  II Setup

###  II.1 Schrödinger-Poisson equation

Our starting point is to write down the basic equation for FDM in a
cosmological background. Ignoring the self-interaction terms, the scalar field
dark matter is generically described by the Klein-Gordon equation:

$\displaystyle\Box\phi-m_{\phi}^{2}\phi=0$ (1)

with $m_{\phi}$ being the mass of FDM. Consider the spacetime metric for the
longitudinal gauge. Assuming the flat geometry, we have

$\displaystyle
ds^{2}=-(1-2\Psi)dt^{2}+a^{2}(1-2\Psi)\,\delta_{ij}dx^{i}dx^{j}.$ (2)

Throughout the paper, we work with the weak-field limit, $|\Psi|\ll 1$, and
take the nonrelativistic limit of scalar field, introducing the following
functional form:

$\displaystyle\phi=\frac{1}{\sqrt{2\,m_{\phi}}\,a^{3/2}}\bigl{[}\varphi\,e^{-i\,m_{\phi}\,t}+\mbox{c.c}\bigr{]}.$
(3)

Since we are interested in the perturbations deep inside the Hubble horizon,
we can make use of the approximation such that $H/m_{\phi}\sim\epsilon\sim
k/m_{\phi}$ and $\Psi\sim\epsilon^{2}$ with $\epsilon$ being a small
parameter, i.e., $\epsilon\ll 1$. Then, the governing equation for the scalar
field at Eq. (1) is shown to be reduced to the Schrödinger-Poisson system
(e.g., Refs. [7, 10, 9]). Adopting the natural units$(\hbar=c=1)$, we have

$\displaystyle i\frac{\partial}{\partial
t}\varphi=\Bigl{[}-\frac{1}{2m_{\phi}\,a^{2}}\nabla_{x}^{2}+m_{\phi}\,\Psi\Bigr{]}\varphi.$
(4)

This has to be solved with the Poisson equation below:

$\displaystyle\frac{1}{a^{2}}\nabla^{2}_{x}\Psi=4\pi\,G\,\delta\rho_{\phi},$
(5)

where the density fluctuation $\delta\rho_{\phi}$ is given by

$\displaystyle\delta\rho_{\phi}\equiv\frac{m_{\phi}}{a^{3}}\big{\\{}|\varphi|^{2}-\langle|\varphi|^{2}\rangle\\}.$
(6)

Equations (4)$\sim$(6) are the basic equations to describe the cosmic
structure in FDM model. Note that taking the large $m_{\phi}$ limit, these
equations have been also used to describe the structure formation in the CDM
model [42], in particular, as an alternative way to solve the Vlasov-Poisson
equation (e.g., Refs. [43, 44, 45]).

###  II.2 FDM in a spherical halo

A primary goal of this paper is to derive analytically the soliton core
structure and to address its relation to the halo properties. Numerical
simulations show that the soliton core forms at the center of halos, and the
size and mass of the solitons are small enough compared to those of halos.
Thus, at a first-order approximation for analytical study, we shall ignore the
self-gravity of the soliton core, and treat the surrounding halo as a
background. The impact of the soliton self-gravity will be later discussed
(see Sec. V). It is to be noted that halos found in simulations have large
fluctuations in density, forming granular structures. While these fluctuations
affect the soliton core and induce the motion and oscillations through the
interaction with them, these nonstationary behaviors can be understood as the
results of wave interference and are well described by a superposition of the
ground (soliton) and excited states in a fixed potential [30, 31]. In this
respect,constructing the eigenstates under a smooth halo background would be
an important first step to describe the dynamical features of the soliton
core.

Based on these considerations, we consider a spherically symmetric halo with
the density profile described by the Navarro-Frenk-While (NFW) profile [46,
47]111We are interested in the region where the density is much larger than
the mean density, and the contribution of the term
$\langle|\varphi|^{2}\rangle$ in Eq. (6) is ignored. :

$\displaystyle m_{\phi}\,|\varphi|^{2}$ $\displaystyle\simeq\rho_{\rm
NFW}(r)\equiv\frac{\rho_{\rm s}}{(r/r_{\rm s})(1+r/r_{\rm s})^{2}}.$ (7)

Note that the density $\rho_{\rm NFW}$ defined above differs from the one
frequently used in the literature, and a factor of $1/a^{3}$ is factored out.
Plugging Eq. (7) into Eq. (6), the Poisson equation at Eq. (5) is analytically
solved, and the potential is expressed as (e.g., see Appendix D of Ref.
[48])222Equation (8) is obtained under the assumption that the density
fluctuation $\delta\rho_{\phi}$ follows the NFW profile even outside the
virial radius. While this is not in reality true in both the observations and
simulations, the modification of the outer profile hardly changes the inner
potential structure, and hence the analysis based on Eq. (8) is relevant to
describe the soliton core structure.

$\displaystyle\Psi=-4\pi\,G\,\frac{\rho_{\rm s}\,r_{\rm
s}^{2}}{a}\,\frac{\log(1+r/r_{\rm s})}{r/r_{\rm s}}.$ (8)

Treating this halo potential as a background, solving the Schrödinger-Poisson
(S-P) equation at Eq. (5) is reduced to a linear eigenvalue problem. Writing
the wave function in the form as

$\displaystyle\varphi=\sum_{n\ell m}u_{n\ell}(r)\,Y_{\ell
m}(\theta,\phi)\,e^{-i\,\epsilon\tau},$ (9)

where the quantity $\tau$ is the new time variable defined by
$\tau\equiv\int^{t}dt^{\prime}/\\{a(t^{\prime})\\}^{2}$, Eq. (4) with the
potential (8) yields

$\displaystyle-\frac{1}{2m_{\phi}}\frac{1}{r^{2}}\frac{d}{dr}\Bigl{(}r^{2}\frac{du_{n\ell}(r)}{dr}\Bigr{)}+\Bigl{[}-4\pi\,G\,m_{\phi}\,\rho_{\rm
s}\,r_{\rm s}^{2}\,a\,\frac{\log(1+r/r_{\rm s})}{r/r_{\rm
s}}+\frac{\ell(\ell+1)}{2m_{\phi}\,r^{2}}\Bigr{]}u_{n\ell}(r)=\epsilon\,u_{n\ell}(r).$
(10)

Equation (10) is further rewritten in a simplified form by introducing
dimensionless quantities:

$\displaystyle x\equiv\frac{r}{r_{\rm s}},\qquad\mathcal{E}\equiv
2m_{\phi}\,r_{\rm s}^{2}\epsilon,\qquad\alpha\equiv
8\pi\,G\,m_{\phi}^{2}\,\rho_{\rm s}\,r_{\rm s}^{4}\,a.$ (11)

We have

$\displaystyle-\frac{1}{x^{2}}\frac{d}{dx}\Bigl{(}x^{2}\frac{du_{n\ell}(x)}{dx}\Bigr{)}+\Bigl{[}-\alpha\,\frac{\log(1+x)}{x}+\frac{\ell(\ell+1)}{x^{2}}\Bigr{]}u_{n\ell}(x)=\mathcal{E}\,u_{n\ell}(x),$
(12)

Here, the quantity $\alpha$ characterizes the depth of the potential well, and
it typically takes a large value. Using the quantities characterizing the
properties of halos in Appendix A, we obtain

$\displaystyle\alpha$ $\displaystyle=7.21\times
10^{2}\,\Bigl{(}\frac{m_{\phi}}{10^{-22}{\rm
eV}}\Bigr{)}^{2}\,\Bigl{(}\frac{M_{\rm
h}}{10^{9}\,M_{\odot}}\Bigr{)}^{4/3}\Bigl{(}\frac{\Delta_{\rm
vir}(z)}{200}\Bigr{)}^{-1/3}\Bigl{(}\frac{\Omega_{\rm
m,0}h^{2}}{0.147}\Bigr{)}^{-1/3}\,\frac{f(c_{\rm vir})}{1+z}$ (13)

with the function $f$ defined by

$\displaystyle f(c_{\rm vir})\equiv\frac{1}{c_{\rm vir}\\{\ln(1+c_{\rm
vir})-c_{\rm vir}/(1+c_{\rm vir})\\}}.$ (14)

The quantity $c_{\rm vir}$ is the concentration parameter, $c_{\rm vir}\equiv
r_{\rm vir}/r_{\rm s}$, with $r_{\rm vir}$ being the virial radius of the halo
having the mass of $M_{\rm h}$. For typical values of $c_{\rm vir}\sim 2-10$,
the function $f$ varies from $f\sim 1$ to $0.07$ (see right panel of Fig. 1).
Left panel of Fig. 1 shows the effective potential normalized by the parameter
$\alpha$, i.e., $V(x)/\alpha\equiv-\log(1+x)/x+\ell(\ell+1)/x^{2}/\alpha$,
where the parameter $\alpha$ is set to $10^{3}$. Note that the parameter
$\alpha$ is not strictly a constant value, but it depends explicitly on the
redshift. In this respect, Eqs. (10) or (12) cannot be reduced to a stationary
problem. Nevertheless, estimating a typical timescale for the wave function,
defined by $\tau_{\rm osc}\equiv 2\pi/|\epsilon|$, gives

$\displaystyle\tau$ $\displaystyle=\frac{1}{4\pi\,G\,m_{\phi}\,\rho_{\rm
s}\,r_{\rm s}^{2}}\,\Bigl{(}\frac{|\mathcal{E}|}{\alpha}\Bigr{)}^{-1}$
$\displaystyle=1.34\times 10^{7}\,\,[{\rm
year}]\,\,\Bigl{(}\frac{|\mathcal{E}|}{\alpha}\Bigr{)}^{-1}\Bigl{(}\frac{m_{\phi}}{10^{-22}\,{\rm
eV}}\Bigr{)}^{-1}\,\Bigl{(}\frac{M_{\rm
h}}{10^{9}\,M_{\odot}}\Bigr{)}^{-2/3}\Bigl{(}\frac{\Delta_{\rm
vir}}{200}\Bigr{)}^{-1/3}\Bigl{(}\frac{\Omega_{\rm
m,0}h^{2}}{0.147}\Bigr{)}^{-1/3}\,\Bigl{\\{}\frac{g(c_{\rm
vir})}{0.1}\Bigr{\\}}$ (15)

with the function $g(c_{\rm vir})$ defined by

$\displaystyle g(c_{\rm vir})\equiv\frac{1}{c_{\rm vir}^{2}\,f(c_{\rm vir})},$
(16)

which typically has $g\sim 0.15-0.2$ for a broad range of $c_{\rm vir}$ (see
right panel of Fig. 1). As we will see later, the typical value of the energy
level is given by $|\mathcal{E}|\sim\alpha$ (Sec. III.2). Hence, the
oscillation period of the wave function is sufficiently shorter than the
cosmological timescales, and thus we can ignore the time dependence of the
parameter $\alpha$.

Figure 1: Left: effective potential of the Schrödinger-Poisson equation in a
NFW halo [see Eq. (12)] normalized by $\alpha$, i.e.,
$V(x)/\alpha\equiv-\log(1+x)/x+\ell(\ell+1)/x^{2}/\alpha$. The results are
plotted against the dimensionless radius, $x\equiv r/r_{\rm s}$. Blue, orange,
and green lines are the cases with $\ell=0$, $1$ and $5$, respectively. Here,
we set the parameter $\alpha$ to $10^{3}$ for illustration. Right: functions
$f(c_{\rm vir})$ and $g(c_{\rm vir})$, defined at Eqs.(14) and (16), as a
function of concentration parameter $c_{\rm vir}$.

##  III Analytical solutions for $\ell=0$

In this section, in order to describe the soliton core structure and its
properties, we focus on the cases with zero angular momentum ($\ell=0$), and
obtain an approximate solution for the radial part of the S-P equation. In
Appendix B, we also construct the approximate eigenstates for the nonzero
angular momentum, and the results are compared with numerical solutions.

###  III.1 Constructing analytical eigenfunctions and eigenvalues

Hereafter, dropping the subscript $\ell$, we denote the radial wave function
for $\ell=0$ by $u_{n}$. Introducing the new function $\tilde{u}_{n}\equiv
x\,u_{n}$ and setting the angular momentum $\ell$ to zero, we rewrite Eq. (12)
with the normal form below:

$\displaystyle\frac{d^{2}\tilde{u}_{n}(x)}{dx^{2}}+\alpha\,g_{0}(x)\tilde{u}_{n}(x)=0\,;\quad
g_{0}(x)\equiv\frac{\log(1+x)}{x}+\frac{\mathcal{E}}{\alpha}.$ (17)

We shall construct the bound-state solutions satisfying the boundary
condition333The inner boundary condition $\tilde{u}(0)=0$ ensures that the
mass profile, $M(x)\propto\int_{0}^{x}|u(x^{\prime})|^{2}\,x^{\prime
2}dx^{\prime}=\int_{0}^{x}|\tilde{u}(x^{\prime})|^{2}dx^{\prime}$, is regular
at $x=0$, and we have $M(0)=0$.:

$\displaystyle\tilde{u}_{n}(0)=0,\qquad\tilde{u}_{n}(\infty)=0.$ (18)

To obtain an approximate solution for the bound-state FDM distribution, we
first notice that the function $g_{0}(x)$ in Eq. (17) is a monotonically
decreasing function of $x$, and we have $g_{0}(0)=1+\mathcal{E}/\alpha$ and
$g_{0}(\infty)=\mathcal{E}/\alpha$. Since the eigenvalue $\mathcal{E}$ for the
bound-state solution must be negative ($\mathcal{E}<0$), the function $g_{0}$
should have a single zero-crossing (or turning) point. Denoting its position
by $x_{c}$, we have

$\displaystyle
g_{0}(x_{c})=\frac{\log(1+x_{c})}{x_{c}}+\frac{\mathcal{E}}{\alpha}=0.$ (19)

Note that the above equation implies that the position $x_{c}$ depends on the
eigenvalue $\mathcal{E}$. The eigenvalues are given as a set of discrete
values associated to the eigenfunction $\tilde{u}$, and will be later
determined by the boundary condition at Eq. (18).

Keeping in mind the fact that Eq. (17) has a single turning point, we consider
the following transformation, known as the Liouville-Green transformation (see
e.g., [49]):

$\displaystyle z=p(x),\qquad v_{n}(z)=\sqrt{p^{\prime}(x)}\,\tilde{u}_{n}(x).$
(20)

The prime denotes the derivative with respect to $x$. With this
transformation, Eq. (17) becomes

$\displaystyle\frac{d^{2}v_{n}(z)}{dz^{2}}+\Bigl{[}\alpha\frac{g_{0}(x)}{\\{p^{\prime}(x)\\}^{2}}+\delta\Bigr{]}v_{n}(z)=0,$
(21)

where the quantity $\delta$ is defined by

$\displaystyle\delta\equiv\frac{q}{\\{p^{\prime}\\}^{2}}\Bigl{(}\frac{1}{q}\Bigr{)}^{\prime\prime},\,\,\quad
q=\sqrt{p^{\prime}}.$ (22)

Here, we choose the function $p$ so as to satisfy444This choice is known as
the Langer’s transformation [39]. See Chapter 14.6 of Ref. [49]

$\displaystyle\alpha\,\frac{g_{0}(x)}{\\{p^{\prime}(x)\\}^{2}}=-z.$ (23)

Recall from Eq. (20) that $z$ is equal to $p$, the above equation is solved to
give an explicit form of $p$ as follows. First we take the square root of both
sides and obtain

$\displaystyle p^{1/2}p^{\prime}=\pm\alpha^{1/2}\,\sqrt{-g_{0}(x)}$ (24)

for $x>x_{c}$. Upon separation of variables, this becomes

$\displaystyle p^{1/2}dp=\pm\alpha^{1/2}\,\sqrt{-g_{0}(x)}\,dx.$ (25)

Integrating once, we have

$\displaystyle\frac{2}{3}z^{3/2}=\frac{2}{3}p^{3/2}=\pm\alpha^{1/2}\,\int_{x_{c}}^{x}\,\sqrt{-g_{0}(x^{\prime})}\,dx^{\prime}.$
(26)

Here we set the lower value of the integral to $x_{c}$. Taking the positive
sign of Eq. (26), we obtain

$\displaystyle
z=p(x)=\alpha^{1/3}\,\Bigl{[}\frac{3}{2}\int_{x_{c}}^{x}\sqrt{-g_{0}(x^{\prime})}\,dx^{\prime}\Bigr{]}^{2/3},$
(27)

which is valid at $x>x_{c}$ (i.e., $z>0$). Similarly, we obtain the relation
valid at $x\leq x_{c}$ (or $z\leq 0$):

$\displaystyle
z=p(x)=-\alpha^{1/3}\,\Bigl{[}\frac{3}{2}\int^{x_{c}}_{x}\sqrt{g_{0}(x^{\prime})}\,dx^{\prime}\Bigr{]}^{2/3}.$
(28)

With Eqs. (27) and (28), the function $z$ yields a smooth and monotonically
increasing function of $x$, and its sign changes from negative to positive at
$x_{c}$.

Using Eq. (23) and the relations given above, the differential equation given
at Eq. (21) is further recast as

$\displaystyle\frac{d^{2}v_{n}(z)}{dz^{2}}-z\,v_{n}(z)=-\delta\,v_{n}(z).$
(29)

Equations (27) and (28) ensure that the quantity $\delta$ defined at Eq. (22)
is shown to be of the order of $\mathcal{O}(\alpha^{-2/3})$ for all values of
$x$. Thus, in the limit of $\alpha\gg 1$, we can safely ignore the right-hand
side involving the factor $\delta$, and the above equation is reduced to

$\displaystyle\frac{d^{2}v_{n}(z)}{dz^{2}}-z\,v_{n}(z)\simeq 0.$ (30)

The general solution of this equation is known to be expressed in terms of the
Airy functions:

$\displaystyle
v_{n}(z)=\tilde{c}_{1}\,\mbox{Ai}(z)+\,\tilde{c}_{2}\,\mbox{Bi}(z)$ (31)

with $\tilde{c}_{1}$ and $\tilde{c}_{2}$ being the integration constants. Note
that the asymptotic behaviors of the Airy functions are (e.g., [50])

$\displaystyle\mbox{Ai}(z)\sim\left\\{\begin{array}[]{ll}{\displaystyle\frac{e^{-(2/3)z^{3/2}}}{2\sqrt{\pi}\,z^{1/4}}}&(z\to\infty)\\\
\\\
{\displaystyle\frac{\sin\Bigl{[}(2/3)(-z)^{3/2}+\pi/4\Bigr{]}}{\sqrt{\pi}\,(-z)^{1/4}}}&(z\to-\infty)\end{array}\right.,\quad\mbox{Bi}(z)\sim\left\\{\begin{array}[]{ll}{\displaystyle\frac{e^{(2/3)z^{3/2}}}{\sqrt{\pi}\,z^{1/4}}}&(z\to\infty)\\\
\\\
{\displaystyle\frac{\cos\Bigl{[}(2/3)(-z)^{3/2}+\pi/4\Bigr{]}}{\sqrt{\pi}\,(-z)^{1/4}}}&(z\to-\infty)\end{array}\right..$
(38)

Thus, imposing the boundary condition $\tilde{u}_{n}(\infty)=0$, the constant
$\tilde{c}_{2}$ must be zero. Another boundary condition,
$\tilde{u}_{n}(0)=0$, leads to the following condition,

$\displaystyle\mbox{Ai}\bigl{(}z(0)\bigr{)}=0,$ (39)

where we used the fact that both of the functions $z(x)$ and $g_{0}(x)$ are
nonvanishing at $x=0$.

Equation (39) determines the eigenvalues $\mathcal{E}$ through the function
$z(0)$ given at Eq. (28), in which the integrand $\sqrt{g_{0}(x)}$ manifestly
depends on $\mathcal{E}$ [see Eq. (17] for the definition of $g_{0}(x)$). It
is the transcendental equation and has to be evaluated numerically to get the
eigenvalues. Although this can be easily and quickly done with the standard
numerical techniques, we can derive an analytical expression for energy
eigenvalues by further employing an approximation. Since the quantity $z$ is
of the order of $\mathcal{O}(\alpha^{1/3})$ and $z<0$ at $x=0$, we can make
use of the asymptotic form of the Airy function Ai by taking the limit of
$\alpha\gg 1$. Then, Eq. (39) is rewritten with [see Eq. (38)]

$\displaystyle\frac{1}{\sqrt{\pi}\,\\{-z(0)\\}^{1/4}}\sin\Bigl{(}\frac{2}{3}\\{-z(0)\\}^{3/2}+\frac{\pi}{4}\Bigr{)}=0.$
(40)

Note that Eq. (28) implies $z(0)\neq 0$. Thus, the condition given above is
recast as [using Eq. (28)]

$\displaystyle\frac{2}{3}\\{-z(0)\\}^{2/3}+\frac{\pi}{4}=n\,\pi\,\,$
$\displaystyle\qquad\Longleftrightarrow\,\,\alpha^{1/2}\int_{0}^{x_{c}}\sqrt{g_{0}(x^{\prime})}\,dx^{\prime}=\Bigl{(}n-\frac{1}{4}\Bigr{)}\pi$
(41)

with the variable $n$ being the positive integer (i.e., $n=1,\,2,\,\cdots$).
Assuming further $x_{c}\ll 1$ and expanding the function $g_{0}(x)$ as
$g_{0}(x)\simeq 1-x/2+\mathcal{E}/\alpha$, the turning point is easily
evaluated to give $x_{c}=2(1+\mathcal{E}/\alpha)$. Substituting the relations
into the above, the integration is analytically performed, and Eq. (41)
becomes

$\displaystyle\alpha^{1/2}\,\frac{4}{3}\Bigl{(}1+\frac{\mathcal{E}}{\alpha}\Bigr{)}^{3/2}=\Bigl{(}n-\frac{1}{4}\Bigr{)}\pi.$
(42)

Solving the above equation with respect to $\mathcal{E}$, we finally get an
analytical expression for the energy eigenvalues:

$\displaystyle\Bigl{(}\frac{\mathcal{E}}{\alpha}\Bigr{)}_{\rm
approx}=-1+\alpha^{-1/3}\Bigl{\\{}\frac{3}{4}\Bigl{(}n-\frac{1}{4}\Bigr{)}\pi\Bigr{\\}}^{2/3}.$
(43)

Note that the $n=1$ implies the ground state. The above analytical expression
is in general valid as long as we consider $\alpha\gg 1$ and for a small
integer $n$.

In Table 1, normalizing the energy eigenvalues by $\alpha$, we list their
numerical values for the lowest five eigenstates. Here, we specifically show
the results of $\alpha=10^{3}$. Table 1 demonstrates how well the analytical
estimation at Eq. (43) can work. Compared to the results obtained from
numerical calculation and Eq. (39), the eigenvalues from (43) start to deviate
as increasing $n$. On the other hand, results obtained by solving the
transcendental equation (39) are surprisingly accurate even for $n=5$. We will
see in next subsection how this semianalytical estimation works well in more
general cases with various values of $\alpha$.

Table 1: Comparison of analytically and numerically estimated eigenvalues for the zero angular-momentum case. Setting the parameter $\alpha$ to $10^{3}$, numerical values of $\mathcal{E}$ are summarized by normalizing them by $\alpha$, particularly for the lowest five eigenstates.  $\mathcal{E}/\alpha$ | Numerical | Eq. (39) | Eq. (43)
---|---|---|---
$n=1$ | -0.86680 | -0.86687 | -0.85383
$n=2$ | -0.78318 | -0.78323 | -0.74286
$n=3$ | -0.72307 | -0.72311 | -0.65244
$n=4$ | -0.67531 | -0.67535 | -0.57261
$n=5$ | -0.63553 | -0.63557 | -0.49965

To sum up, the analytical solution for the bound-state radial wave function
which satisfies the boundary conditions $\tilde{u}(0)=0=\tilde{u}(\infty)$ is
given by

$\displaystyle\tilde{u}_{n}(x)=x\,u_{n}(x)=\left\\{\begin{array}[]{ll}{\displaystyle\frac{\\{z(x)\\}^{1/4}}{\\{-g_{0}(x)\\}^{1/4}}\,\mbox{Ai}(z(x))}\,;\quad
z(x)=\alpha^{1/3}\Bigl{[}\frac{3}{2}\int_{x_{c}}^{x}\sqrt{-g_{0}(x^{\prime})}\,dx^{\prime}\Bigr{]}^{2/3},&(x>x_{c})\\\
\\\
{\displaystyle\frac{\\{-z(x)\\}^{1/4}}{\\{g_{0}(x)\\}^{1/4}}\,\mbox{Ai}(z(x))}\,;\quad
z(x)=-\alpha^{1/3}\Bigl{[}\frac{3}{2}\int_{x}^{x_{c}}\sqrt{g_{0}(x^{\prime})}\,dx^{\prime}\Bigr{]}^{2/3},&(0\leq
x\leq x_{c})\end{array}\right.$ (47)

with the function $g_{0}$ defined at Eq. (17). The eigenvalue $\mathcal{E}$ is
obtained by solving Eq. (39) for a given $\alpha$, and the turning point
$x_{c}$ is determined by $g_{0}(x_{c})=0$. Note that the analytical
expressions given above is not normalized, and for a proper definition, it has
to be divided by the normalization constant $\mathcal{N}$, defined by
$\mathcal{N}\equiv\\{\int
dx^{\prime}|\tilde{u}_{n}(x^{\prime})|^{2}\\}^{1/2}$.

Figure 2: Eigenfunctions for radial wave function $\tilde{u}_{n}(x)$ for
$\alpha=10$ (left), $10^{2}$ (middle), and $10^{3}$ (right). Results of the
lowest five eigenstates are plotted as a function of dimensionless radius
$x=r/r_{\rm s}$. shown. In the upper panels, while the thin solid lines are
the analytical results based on Eq. (47), thick dashed lines represent the
numerical solutions obtained by solving the matrix eigenvalue problem. Note
that all the results are normalized. On the other hand, the bottom panels show
the fractional difference between the analytical and numerical results defined
by $\Delta\tilde{u}_{n}/\tilde{u}_{n,{\rm numerical}}$ with
$\Delta\tilde{u}_{n}=\tilde{u}_{n,{\rm analytical}}-\tilde{u}_{n,{\rm
numerical}}$, where the functions $\tilde{u}_{n,{\rm analytical}}$ and
$\tilde{u}_{n,{\rm numerical}}$ denote the analytical and numerical
eigenstates, respectively. Here, to avoid the division by zero, we stop
plotting the results when the amplitude of the wave functions becomes smaller
than $10^{-9}$ at the outer part.

Figure 3: Same as Fig. 3, but we here plot the square of the wave function
$u_{n}(x)=x\,\tilde{u}_{n}(x)$, normalized by the one evaluated at the origin,
i.e., $|u_{n}(x)|^{2}/|u_{n}(0)|^{2}$. Note that the plotted results
correspond to the density profile normalized by the central density,
$\rho(x)/\rho(0)$.

###  III.2 Comparison with numerical solutions

Having obtained the analytical expressions for eigenvalues and eigenfunctions,
we now compare them with numerical solutions. Based on the standard technique
to solve the stationary problem of Schrödinger equation, we discretize Eq.
(17). Then, solving differential equation under the boundary conditions at Eq.
(18) is reduced to a matrix eigenvalue problem in the following form (e.g.,
Ref. [51]):

$\displaystyle\left(\begin{array}[]{ccccc}2/\Delta^{2}+V(x_{1})&-1/\Delta^{2}&\cdots&0&0\\\
-1/\Delta^{2}&2/\Delta^{2}+V(x_{2})&-1/\Delta^{2}&\cdots&0\\\
\vdots&\ddots&\ddots&\ddots&\vdots\\\
0&\ddots&-1/\Delta^{2}&2/\Delta^{2}+V(x_{n-1})&-1/\Delta^{2}\\\
0&\cdots&0&-1/\Delta^{2}&2/\Delta^{2}+V(x_{n})\\\
\end{array}\right)\left(\begin{array}[]{c}\tilde{u}(x_{1})\\\
\tilde{u}(x_{2})\\\ \vdots\\\ \tilde{u}(x_{n-1})\\\
\tilde{u}(x_{n})\end{array}\right)=\mathcal{E}\left(\begin{array}[]{c}\tilde{u}(x_{1})\\\
\tilde{u}(x_{2})\\\ \vdots\\\ \tilde{u}(x_{n-1})\\\
\tilde{u}(x_{n})\end{array}\right)$ (63)

with the quantity $\Delta$ and the function $V(x)$, respectively, defined by
$\Delta\equiv x_{i+1}-x_{i}$ and $V(x)=-\alpha\,\log(1+x)/x$.

Figures 3 and 3 show the wave functions of the lowest five eigenstates (i.e.,
$n=1,\cdots,5$) for the parameters $\alpha=10$ (left), $10^{2}$ (middle), and
$10^{3}$ (right). In upper panels of Fig. 3, thick dashed lines are the
numerical results of the function $\tilde{u}_{n}$, which are obtained by
setting the inner and outer boundaries to $x_{1}=0$ and $x_{n}=50$ for
$\alpha=10$ and $20$ for $\alpha=10^{2}$ and $10^{3}$ with the number of grids
$n=10^{4}$. These results are compared to the analytical results depicted as
thin solid lines, with the amplitude of each eigenfunction properly
normalized. Bottom panels of Fig. 3 plot the fractional difference between the
analytical and numerical results, defined by $(\tilde{u}_{n,{\rm
analytical}}-\tilde{u}_{n,{\rm numerical}})/\tilde{u}_{n,{\rm numerical}}$. On
the other hand, Fig. 3 shows the square of the wave function
$u_{n}=x\,\tilde{u}_{n}$, normalized it by the one evaluated at the origin,
i.e., $|u_{n}(x)|^{2}/|u_{n}(0)|^{2}$, which corresponds to the density
profile normalized by the central density, $\rho(x)/\rho(0)$. The meaning of
line types and colors are the same as those shown in the upper panels of Fig.
3.

Clearly the agreement between analytical numerical results is excellent even
for a small value of $\alpha$, with which we naively expect the analytical
prediction to be inaccurate. A closer look at the fractional difference
reveals a discrepancy at larger radii $x$ for small $\alpha$, where spiky
features arising from the zero-crossing points also become prominent.
Nevertheless, these behaviors appear manifest only when the wave functions
fall off and approach zero, thus giving no serious impact even from the
quantitative point of view. Indeed, looking at the density profiles plotted in
logarithmic scales (Fig. 3), we hardly see a difference between analytical and
numerical results.

Next look at the eigenvalues. Table 1 summarizes the results in the case of
$\alpha=10^{3}$. Further, we plot in Fig. 4 the eigenvalues for the lowest
five eigenstates as a function of $\alpha$. Here, the eigenvalues computed
from Eqs. (39) and (43) are shown, depicted as solid and dotted lines,
respectively. Again, we see an excellent agreement between numerical and
analytical estimations. It is surprising that the analytical estimation with
Eq. (39) remains accurate even at $\alpha\lesssim 1$, and reproduce numerical
results. Thus, we conclude that the solutions constructed analytically in Sec.
III.1 provide a fast and reliable estimate of the eigenfunctions and
eigenvalues, which can be used to study quantitatively the soliton core
properties.

Figure 4: Left: energy eigenvalue plotted as a function of $\alpha$. To be
precise, the plotted results are $1+\mathcal{E}/\alpha$, and results for the
lowest five eigenstates are shown. Lines from bottom to top indicate the
results from $n=1$ to $5$. Right: fractional error of the analytically
estimated values of the eigenvalue $\mathcal{E}$,
$\Delta\mathcal{E}/\mathcal{E}_{\rm numerical}$. Again, results for the lowest
five eigenstates are shown.

###  III.3 Analytical estimation of soliton core

In this subsection, before addressing the core-halo relations, we shall
compare the ground-state wave function ($n=1$) with the soliton density
profile found in numerical simulations.

It has been found in numerical simulations that the central core structure of
FDM halos is well described by the following fitting form [52, 53]:

$\displaystyle\rho_{\rm soliton}(r)=\frac{\rho_{\rm c}}{\\{1+\gamma(r/r_{\rm
c})^{2}\\}^{8}},$ (64)

with the central density $\rho_{\rm c}$ given by555Here, the quantities
$\rho_{\rm c}$ and $r_{\rm c}$ are defined as the comoving density and
comoving radius, respectively.

$\displaystyle\rho_{\rm
c}=\frac{0.019}{a}\,\Bigl{(}\frac{m_{\phi}}{10^{-22}}\Bigr{)}^{-2}\Bigl{(}\frac{r_{\rm
c}}{1\,\mbox{kpc}}\Bigr{)}^{-4}\,\,[M_{\odot}\mbox{pc}^{-3}].$ (65)

In the above, the core size $x_{\rm c}$ is defined to be the radius at which
the density drops to one-half of the central density, i.e., $\rho_{\rm
soliton}(x_{\rm c})=\rho_{\rm c}/2$. This gives the constant
$\gamma=2^{1/8}-1\simeq 0.091$.

Similarly, as we have seen in Fig. 3, the ground-state eigenfunctions ($n=1$),
depicted as black curves, commonly have a flat core followed by a sharp drop,
irrespective of the parameter $\alpha$. Although our treatment ignoring the
self-gravity cannot precisely predict the soliton core density $\rho_{\rm c}$,
making use of the asymptotic properties of the Airy function, one can
analytically express the size of the soliton core structure near the origin.
From Eq. (47), we have

$\displaystyle u_{n}(x)$ $\displaystyle=\frac{\tilde{u}_{n}(x)}{x}$
$\displaystyle\stackrel{{\scriptstyle x\ll
1}}{{\longrightarrow}}\frac{1}{\sqrt{\pi}\,\bigl{\\{}g_{0}(x)\bigr{\\}}^{1/4}}\,\sin\Bigl{[}\frac{2}{3}\bigl{\\{}-z(x)\bigr{\\}}^{3/2}+\frac{\pi}{4}\Bigr{]}$
$\displaystyle\simeq(-1)^{n+1}\frac{\alpha^{1/2}\\{g_{0}(0)\\}^{1/4}}{\pi^{1/2}}\,\Bigl{(}1-\beta\,x^{2}+\cdots\Bigr{)},$
(66)

where the function $z$ is Taylor expanded in the last line. Rewriting further
the derivatives of $z$ in terms of those of the function $g_{0}$, the
coefficient $\beta$ is expressed with a help of Eqs. (23) and (28) as

$\displaystyle\beta$
$\displaystyle=\frac{\alpha}{6}\,g_{0}(0)+\frac{5}{96}\Bigl{\\{}\frac{g_{0}^{\prime}(0)}{g_{0}(0)}\Bigr{\\}}^{2}-\frac{1}{24}\frac{g_{0}^{\prime\prime}(0)}{g_{0}(0)}.$
(67)

Note that in deriving Eq. (67), we do not use any functional form of $g_{0}$.
To get a more explicit expression, we substitute the function $g_{0}(x)$ given
at Eq. (17). We then obtain

$\displaystyle\beta=\frac{\alpha}{6}\Bigl{(}1+\frac{\mathcal{E}}{\alpha}\Bigr{)}+\frac{17+32\mathcal{E}/\alpha}{1152(1+\mathcal{E}/\alpha)^{2}}.$
(68)

Equating especially the ground-state eigenfunction
$|u_{1}(x)|^{2}/|u_{1}(0)|^{2}$ with $\rho_{\rm soliton}(x)/\rho_{\rm c}$ at
Eq. (64), we find that the dimensionless core radius $x_{\rm c}$, defined by
$x_{\rm c}\equiv r_{\rm c}/r_{\rm s}$, is expressed as

$\displaystyle x_{\rm
c}=p\,\beta^{-1/2}=p\,\sqrt{\frac{6}{\alpha(1+\mathcal{E}/\alpha)}},$ (69)

with the constant $p$ chosen to be $p=0.65$. Here, we have ignored the second
term in Eq. (68), which merely gives a subdominant contribution. To be
precise, in Eq. (69), a relevant numerical factor is $p=\sqrt{4\gamma}\simeq
0.60$, but we here adopt the valpue of $0.65$, which in fact gives a better
agreement with Eq. (69) when comparing the analytically estimated density
profile with the fitting formula below.

Figure 5: Comparison of the ground-state $(n=1)$ density profile (black solid)
with the fitting function of the soliton density profile (red dashed). Plotted
results are the normalized density profile, $\rho(x)/\rho(0)$, so that it
approaches unity at $x\to 0$. From left to right panels, the results for
$\alpha=5\times 10^{3}$, $10^{4}$, and $5\times 10^{4}$ are shown. In plotting
the fitting function in Eq. (64), we adopt the core radius $x_{\rm c}$ given
at Eq. (69). Note that the ground-state density profiles are identical to
those shown in Fig. 3 (see thin black solid lines).

Using the expression of the core radius in Eq. (69), Fig. 5 compares the
normalized density profile obtained from the analytical expression of the
ground-state wave function with Eq. (64). Plugging Eq. (69) into Eq. (64), the
analytical expression depicted as solid lines describes remarkably well the
fitting function (red dashed) over a wide range of the parameter $\alpha$ even
at the outskirt of the profile. Thus, our analytical result can give a good
description for the simulated soliton profile. In particular, Eq. (69) is the
key to derive analytically the core-halo relation, which we will discuss in
more detail.

##  IV Predicting soliton core-halo relations

In this section, on the basis of the analytical expression at Eq. (69), we
investigate in detail properties of the soliton core, focusing particularly on
its relation to the halo mass and other parameters.

###  IV.1 Core-halo mass-concentration relations

Let us recast the expression at Eq. (69) with the comoving core radius $r_{\rm
c}$, showing explicitly the parameter dependence. Through the relation $r_{\rm
c}=x_{\rm c}\,r_{\rm s}$ with the scale radius given by $r_{\rm s}=r_{\rm
vir}/c_{\rm vir}$, using the expression of $\alpha$ at Eq. (13) and the
definition of $r_{\rm vir}$ at Eq. (88) gives

$\displaystyle r_{\rm c}$
$\displaystyle=p\sqrt{\frac{6}{1+\mathcal{E}/\alpha}}\frac{r_{\rm vir}/c_{\rm
vir}}{\alpha^{1/2}}$ $\displaystyle\simeq
1.83\,\,\mbox{[kpc]}\,\,a^{-1/2}\,\Bigl{(}\frac{p}{0.65}\Bigr{)}\,\Bigl{(}\frac{m_{\phi}}{10^{-22}{\rm
eV}}\Bigr{)}^{-1}\Bigl{(}\frac{M_{\rm
h}}{10^{9}\,M_{\odot}}\Bigr{)}^{-1/3}\Bigl{(}\frac{\Delta_{\rm
vir}}{200}\Bigr{)}^{-1/6}\Bigl{(}\frac{\Omega_{\rm
m,0}h^{2}}{0.147}\Bigr{)}^{-1/6}\,\Bigl{\\{}\frac{g(c_{\rm
vir})}{1+\mathcal{E}/\alpha}\Bigr{\\}}^{1/2},$ (70)

where the function $g(c_{\rm vir})$ is defined in Eq. (16).

The above relation is compared to Eq. (7) of Ref. [53], who derived it based
on the kinematic argument with a help of numerical experiments. Apart from the
last factor, the dependencies on FDM mass, halo mass, and cosmological
parameters exactly coincide with each other666In Eq. (7) of Ref. [53], the
cosmological parameters are set to typical values, and their expression of the
core radius apparently misses the dependence on $\Omega_{\rm m,0}h^{2}$. Note
also that Eq. (7) is given as a physical radius, and thus the scale factor
dependence differs from ours by the factor of $a$. . Also, the numerical
prefactor of $1.83$ in our expression is rather close to theirs.

However, an important difference is the last factor, $\\{g(c_{\rm
vir})/(1+\mathcal{E}/\alpha)\\}^{1/2}$. While the numerator of this factor is
given as a function of $c_{\rm vir}$, the denominator, $1+\mathcal{E}/\alpha$,
depends on the parameter $\alpha$. Thus, through Eq. (13), the actual
parameter dependence of $r_{\rm c}$ can differ significantly from the one in
Ref. [53].

To see the explicit dependence of the expression (70) on the halo properties,
we hereafter fix the cosmological parameters and FDM mass, respectively, to
$\Omega_{\rm m,0}=0.276$, $h=0.677$, $\Omega_{\Lambda}=0.724$, and
$m_{\phi}=8\times 10^{-23}$ eV. Then, for a given redshift $z$, the soliton
core size $r_{\rm c}$ is described as a function of the halo mass and
concentration parameter, $M_{\rm h}$ and $c_{\rm vir}$. While the
concentration parameter is also given as a function of the halo mass and
redshift, its dependence is known to be sensitive to the cosmology, mainly
through the initial power spectrum. Hence, as a generic prediction of the
core-halo relation, we first consider $c_{\rm vir}$ to be an independent
parameter of the halo mass. We shall then introduce the models of
concentration-halo mass relation, and discuss their core-halo relations (Sec.
IV.2).

Figure 6: Soliton core-halo relation given as function of halo mass $M_{\rm
h}$ and concentration parameter $c_{\rm vir}$, with the latter dependence
shown in the color scales. The analytical predictions of the core radius vs
halo mass (left) and core mass vs halo mass (right) relations, given
respectively at Eqs. (70) and (72), are plotted at $z=0$. Here, we adopt the
mass of the FDM and the cosmological parameters as follows: $m_{\phi}=8\times
10^{-23}$ eV, $\Omega_{\rm m,0}=0.276$, $\Omega_{\Lambda}=0.724$ and
$h=0.677$. For comparison, we also plot the results measured from the
numerical simulations of S-P equation. The dark green and purple points
respectively represent the results from soliton merger simulations [27] and
cosmological simulations [26]. Note that in each panel, regions with no color
indicate that soliton core is not theoretically allowed to form. The boundary
on the allowed core radius and/or core mass is described by a single power-law
function of halo mass, whose scaling is depicted as gray dashed lines (see the
text in detail).

Left panel of Fig. 6 plots the predicted core radius given as a function of
$M_{\rm h}$ and $c_{\rm vir}$, with the latter dependence shown in the color
scales. In similar manner, right panel plots the predicted core mass. Here, we
adopt the fitting form of the soliton density profile at Eq. (64) with the
soliton core density $\rho_{\rm c}$ given by Eq. (65). The core mass is then
estimated from (e.g., [25]):

$\displaystyle M_{\rm c}\simeq 4\pi(0.2225)\rho_{\rm c}\,r_{\rm c}^{3},$ (71)

which gives

$\displaystyle M_{\rm c}=2.91\times
10^{7}\,\,[M_{\odot}]\,\,a^{-1/2}\,\Bigl{(}\frac{p}{0.65}\Bigr{)}^{-1}\,\Bigl{(}\frac{m_{\phi}}{10^{-22}{\rm
eV}}\Bigr{)}\Bigl{(}\frac{M_{\rm
vir}}{10^{9}\,M_{\odot}}\Bigr{)}^{1/3}\Bigl{(}\frac{\Delta_{\rm
vir}}{200}\Bigr{)}^{1/6}\Bigl{(}\frac{\Omega_{\rm
m,0}h^{2}}{0.147}\Bigr{)}^{1/6}\Bigl{\\{}\frac{g(c_{\rm
vir})}{1+\mathcal{E}_{n}/\alpha}\Bigr{\\}}^{-1/2}.$ (72)

Since the core mass shown here partly uses the relation determined by the
numerical simulations [in particular, the core density at Eq. (65)], Eq. (72)
is not strictly the first-principle prediction, compared to the core radius
shown in the left panel, where we only used the analytical results given in
Sec. III. Nevertheless, the predicted core mass can be used as an independent
cross check, and it provides some insights into the one obtained from
numerical simulations, as we will discuss below.

In both panels of Fig. 6, we also plot the results measured from numerical
simulations, depicted as filled circles. These data are taken from Ref. [27].
Two different colors indicate the results obtained from either soliton merger
simulations (dark green) or simulations started from the cosmological initial
condition (purple, Ref. [26])777To be precise, simulation data of Ref. [26]
have started from the CDM initial condition. Thus, unlike the FDM initial
condition, no small-scale cutoff was imposed in the initial conditions.
Nevertheless, the evolved power spectra measured at later time are shown to
look quantitatively similar to those expected from the FDM initial condition..
Multiplying by the scale factor $a^{1/2}$, the results are all scaled to those
at $z=0$.

In Fig. 6, a cautious remark is that the parameter $\alpha$ becomes smaller
than unity at the region satisfying the condition $1\lesssim\,(M_{\rm
h}/10^{7}M_{\odot})^{3/4}\,f(c_{\rm vir})$, which roughly corresponds to
reddish triangular regions. Although our analytical description of the core-
halo relations become inaccurate there, these parameter regions are somewhat
extreme and there are in fact no simulation data points. Apart from these
points, important findings from Fig. 6 are summarized below:

  * •

All of the simulation results lie at the regions allowed by predictions. To be
precise, measured core radii and core masses are consistent with predictions
with moderately small values of concentration parameter, $c_{\rm
vir}\sim\mathcal{O}(1-10)$.

  * •

In predictions, there exist clear boundaries on the allowed core radius and
core mass for each halo mass. Denoting, respectively, the upper and lower
bounds on the predicted core radius and core mass by $r_{\rm c}^{\rm limit}$
and $M_{\rm c}^{\rm limit}$, one finds that they gradually changes with halo
mass, scaled as $r_{\rm c}^{\rm limit}\propto M_{\rm h}^{-1/9}$ and $M_{\rm
c}^{\rm limit}\propto M_{\rm h}^{1/9}$ (see gray dashed lines). Both of the
bounds appear in the limit of $c_{\rm vir}\to 0$.

Making use of the analytical properties of the wave functions in Sec. III, we
can derive explicit expressions for the critical radius $r_{\rm c}^{\rm
limit}$ and mass $M_{\rm c}^{\rm limit}$ as follows. We first notice that
taking the limit $c_{\rm vir}\to 0$, the parameter $\alpha$ yields
$\alpha\to\infty$. This implies that we can safely use the approximate
expression at Eq. (43) to evaluate the ground-state $(n=1)$ eigenvalue. Then,
writing the parameter $\alpha$ given in Eq. (13) as
$\alpha=\tilde{\alpha}\,f(c_{\rm vir})$, the factor of $\\{g(c_{\rm
vir})/(1+\mathcal{E}/\alpha)\\}^{1/2}$ can be computed explicitly in the limit
of $c_{\rm vir}\to 0$. We obtain

$\displaystyle\Bigl{\\{}\frac{g(c_{\rm
vir})}{1+\mathcal{E}_{n}/\alpha}\Bigr{\\}}^{1/2}\,\,\stackrel{{\scriptstyle
c_{\rm vir}\to
0}}{{\longrightarrow}}\,\,\tilde{\alpha}^{1/6}\,\Bigl{(}\frac{8}{9\pi}\Bigr{)}^{1/3}.$
(73)

Recalling that the quantity $\tilde{\alpha}$ depends on various model
parameters including the FDM mass and halo mass [Eq. (13)], plugging Eq. (73)
into Eq. (70) finally yields the following relation:

$\displaystyle r_{\rm c}^{\rm
limit}=3.59\,[\mbox{kpc}]\,\,a^{-1/3}\,\Bigl{(}\frac{p}{0.65}\Bigr{)}\,\Bigl{(}\frac{m_{\phi}}{10^{-22}}\Bigr{)}^{-2/3}\Bigl{(}\frac{M_{\rm
h}}{10^{9}\,M_{\odot}}\Bigr{)}^{-1/9}\Bigl{(}\frac{\Delta_{\rm
vir}}{200}\Bigr{)}^{-2/9}\Bigl{(}\frac{\Omega_{\rm
m,0}h^{2}}{0.147}\Bigr{)}^{-2/9},$ (74)

which gives the upper bound on the core radius. Similarly, the lower bound on
the soliton core mass is obtained by substituting Eq. (73) into Eq. (72):

$\displaystyle M_{\rm c}^{\rm limit}=1.48\times
10^{7}\,[M_{\odot}]\,\,a^{-2/3}\,\Bigl{(}\frac{m_{\phi}}{10^{-22}}\Bigr{)}^{2/3}\Bigl{(}\frac{M_{\rm
h}}{10^{9}\,M_{\odot}}\Bigr{)}^{1/9}\Bigl{(}\frac{\Delta_{\rm
vir}}{200}\Bigr{)}^{2/9}\Bigl{(}\frac{\Omega_{\rm
m,0}h^{2}}{0.147}\Bigr{)}^{2/9}.$ (75)

Equations (74) and (75) derived above precisely reproduce the boundary curves
shown in Fig. 6. The presence of upper (lower) bound on core radius (mass) are
related to the structure of halo potential near the center. In other words,
the inner structure of halo profile determines the allowable soliton core size
and mass, and they can be indeed changed depending on the inner slope of halo
density profile. In Appendix C, we examine this in detail in a class of
generalized halo profiles, and show that the upper bound on the core radius
$r_{\rm c}^{\rm limit}$ generally scales with respect to the inner slope of
the density profile $s$, defined by $\rho\propto r^{-s}$ $(0\leq s<2)$, as
[Eq. (130)]:

$\displaystyle r_{\rm c}^{\rm
limit}=3.59\,\,\mbox{[kpc]}\,\,\frac{\mathcal{T}(s)}{\mathcal{T}(1)}a^{-1/(4-s)}\,\Bigl{(}\frac{p}{0.65}\Bigr{)}\,\Bigl{(}\frac{m_{\phi}}{10^{-22}\,\mbox{eV}}\Bigr{)}^{-2/(4-s)}\Bigl{(}\frac{M_{\rm
h}}{10^{9}\,M_{\odot}}\Bigr{)}^{-s/3/(4-s)}\Bigl{(}\frac{\Delta_{\rm
vir}}{200}\,\,\frac{\Omega_{\rm m,0}h^{2}}{0.147}\Bigr{)}^{-(1-s/3)/(4-s)},$

where the function $\mathcal{T}$, defined at Eq. (131), is a monotonically
decreasing function of the slope $s$, and it becomes zero when the slope
approaches $2$. Although the above expression is valid at $0\leq s<2$, the
behavior of the function $\mathcal{T}$ suggests that the FDM soliton, if
exists, does not form a flat core for a halo having a steeper inner slope of
$s\geq 2$.

###  IV.2 Core-halo mass relations in CDM and FDM cosmologies

In this subsection, incorporating the concentration-mass (C-M) relation
$c_{\rm vir}(M_{\rm h})$ into the results shown Fig. 6, we compare the
predicted core-halo relations with simulation results.

Figure 7: Concentration-mass (C-M) relation for the FDM and CDM models.
Adopting the analytical models presented in Refs. [54, 55] (see Appendix A.2
for their analytical expressions), the concentration parameter $c_{\rm vir}$
is plotted as a function of halo mass in blue and red lines for FDM and CDM
models, respectively. The faint and dark shaded area indicate the $1\sigma$
and $2\sigma$ errors, assuming the log-normal distribution of $c_{\rm vir}$
with dispersion given by $0.16$ dex.

It is known that the C-M relation acquires various dependencies on cosmology
and redshift through the formation and merger history of halos. As a result,
at low-mass halos of our interest, it becomes rather sensitive to the cutoff
of the initial power spectrum. Here, as representative examples, we examine
the two C-M relations given for different cosmological models. One is the C-M
relation of the FDM model. The linear power spectrum of this model exhibits a
sharp cutoff at the wave number $k_{1/2}\sim
4.5(m_{\phi}/10^{-22}\,\mbox{eV})^{4/9}$ Mpc-1 [6, 7], and this cutoff leads
to the formation of low-concentration halos with a suppressed abundance,
analogous to the case in the warm dark matter model (e.g., Refs. [56, 57]). At
present, there is a little work to numerical study the C-M relation in a
relevant cosmological setup. We adopt the analytical C-M relation proposed by
Ref. [54], who applied it to predict the two-point statistics of galaxies and
weak lensing based on the halo model prescription. Another C-M relation we
examine is the one of the CDM model. In contrast to the FDM model, the power
spectrum of the CDM model does not have any typical cutoff at relevant scales
of structure formation. Thus, a sizable amount of halos is formed, and this
results in the high concentration halos at small halo masses. We adopt the C-M
relation given by Ref. [55], which has been calibrated by cosmological
$N$-body simulations to quantitatively match the measured C-M relation (see
also Refs. [58, 59, 60] for recent improved modeling).

The analytical expressions for the C-M relation of both models are summarized
in Appendix A.2. Based on these, Fig. 7 shows the C-M relations of the FDM
(blue) and CDM (red) models. Since the C-M relation measured in simulations is
known to have a large scatter (e.g., Refs. [61, 55, 62, 63, 64, 65]), we also
show in Fig. 7 the $1$ and 2$\sigma$ errors around the mean C-M relation,
assuming the log-normal distribution with the scatter of $0.16$ dex [65].

Using these C-M relations and their scatter in Fig. 7, we compute the core-
halo relations, and the predictions are plotted in Figs. 9 and 9 for FDM and
CDM models, respectively. Together with the fitting function found numerically
by Ref. [53] (gray dashed), simulation results are also shown as small filled
circles, for which we further estimate the median values and dispersions,
depicted, respectively, as large filled circles and error bars.

Overall, the predicted core-halo relations exhibit a non-power-law behavior
for both the FDM and CDM models. At the small halo masses of $M_{\rm
h}\lesssim 10^{11}\,M_{\odot}$, differences between the predictions become
manifest. Obviously, the results from the C-M relation of the CDM model, which
predicts high-concentration halos, fails to reproduce the trend seen in
numerical simulations. On the other hand, adopting the C-M relation of the FDM
model, the predicted core-halo relations gives a close agreement with the
results obtained from numerical simulations. Interestingly, at $M_{\rm h}\sim
10^{8}-10^{9}\,M_{\odot}$, the predicted core-halo relations get closer to the
scaling relation found by Ref. [53] (gray dashed), which predicts the core-
halo relation in a power-law form of $r_{\rm c}\propto M_{\rm h}^{-1/3}$ and
$M_{\rm c}\propto M_{\rm h}^{1/3}$. Also, scatters in the predicted core-halo
relations resemble those shown in the simulations.

These results suggest that provided the C-M relation well-calibrated with
numerical simulations, our analytical formulas given at Eqs. (70) and (72) can
successfully describe the core-halo relations found in numerical simulations.
Nevertheless, one caveat to be noted is that the data points shown in Figs. 9
and 9 are obtained from the simulations with a small box size, $L=10\,h^{-1}$
Mpc for the cosmological simulations [26], and $L=300$ kpc for the soliton
merger simulations. Further, while the latter simulations are not strictly
made with a cosmological setup, the initial conditions of the former
simulations is not precisely consistent with the FDM model having a small-
scale cutoff. In this respect, in both simulations, the evolved halos may not
necessarily trace the C-M relations expected from those obtained from a
relevant cosmological setup. Indeed, we see a small discrepancy with
predictions, which appears manifest around the halo mass of $M_{\rm h}\sim
10^{9}\,M_{\odot}$. In addition, the analytical C-M relation for the FDM model
is not designed to account for the low-mass halos considered here. In this
respect, the predicted core-halo relations adopting the C-M relation of Ref.
[54] might not be accurate at $M_{\rm h}\lesssim 10^{11}M_{\odot}$. Still, the
results shown in Fig. 9 is very promising, and a more quantitative
investigation along the direction would shed light on clarifying the origin
and diversity of core-halo relations.

Finally, a closer look at the predictions in the FDM model reveals that the
scatter in core-halo relations gets smaller as decreasing the halo mass. This
is because the FDM model prefers low-concentration halos at small mass scales.
As we have seen in Sec. IV.1, the allowable size and mass of the soliton core
become converged to a certain value for each halo mass in the limit of $c_{\rm
vir}\to 0$ [see Eqs. (74) and (75)]. As a result, the scatter in the predicted
core-halo relation gets reduced for low-concentration halos. Indeed, similar
trend can be clearly seen in numerical simulations, and data points from the
two different simulations converge to a similar core radius and core mass at
$M_{\rm h}\sim 10^{8}\,M_{\odot}$, giving consistent results with predictions
of the FDM model.

Figure 8: The soliton core-halo mass relation adopting the C-M relation for
the FDM model. Incorporating the C-M relation of Ref. [54] into the
predictions shown in Fig. 6, the results of the core radius vs halo mass
(left) and the core mass vs halo mass (right) relations are plotted, together
with measured results from numerical simulations, for which the median values
and dispersions are also evaluated in each halo mass bin and are plotted as
large filled circles and error bars, respectively . In plotting the
predictions, we assume, for each halo mass, the log-normal distribution of
$c_{\rm vir}$, and evaluate the median and the scatter in the core-halo
relation. In each panel, the thick solid line is the median relation, while
the faint and dark shaded areas respectively show the $1\sigma$ and $2\sigma$
scatter arising from the scatter in $c_{\rm vir}$. Note that the median
relations shown here are hardly distinguishable from the predictions computed
with the mean C-M relation. For reference, the gray dashed lines are the
scaling relations numerically found by Ref. [53].

Figure 9: Same as Fig. 9, but the predictions adopting the C-M relation of the
CDM model by Ref. [55] are shown.

##  V Discussion: impact of soliton self-gravity

So far, we have ignored the self gravity of soliton core, and constructed
analytically the eigenfunctions of the S-P system under the potential of a
spherically symmetric halo. Focusing particularly on the ground-state
eigenfunction, we derived the analytical expression of the soliton core, which
is then used to predict the soliton core-halo relations. In this section, we
discuss the effect of soliton self-gravity, and estimate in particular its
impact on the soliton core radius.

As we have seen in Sec. II.1, a proper account of the soliton self-gravity is
to add the contribution of the soliton core to the gravitational potential, on
top of the potential of the background halo. The effective potential of Eq.
(12) is then changed to

$\displaystyle-\alpha\frac{\log(1+x)}{x}+\frac{\ell(\ell+1)}{x^{2}}$
$\displaystyle\qquad\longrightarrow\,\,-\alpha\frac{\log(1+x)}{x}+\frac{\ell(\ell+1)}{x^{2}}+\Delta
V(x),$ (76)

where the correction term $\Delta V$ is related to the gravitational potential
of the soliton, $\Psi_{\rm soliton}$:

$\displaystyle\Delta V(x)=\frac{2m_{\phi}^{2}r_{\rm
s}^{2}a^{2}}{\hbar^{2}}\Psi_{\rm soliton}(x).$ (77)

Since the potential $\Psi_{\rm soliton}$ is determined by the Poisson equation
sourced by the square of the wave function itself [see Eqs. (5) with (6)], the
S-P equation having the effective potential at Eq. (76) is no longer reduced
to a linear eigenvalue problem. Nevertheless, if the potential of the self-
gravitating soliton is small enough compared to the background halo potential,
a perturbative treatment can be applied, and the impact of soliton self-
gravity is estimated based on our treatment in Sec. III.

Let us suppose that the soliton core is described by the ground-state
eigenfunction in Sec. III.1, with its density profile approximately given by
the fitting function at Eq. (64). Adopting this fitting form, the Poisson
equation is analytically solved, and the gravitational potential induced by
the soliton core is expressed as [66]

$\displaystyle\Psi_{\rm soliton}(x)$ $\displaystyle=-\pi\,G\,\frac{\rho_{\rm
c}\,r_{\rm s}^{2}}{a}\,\frac{x_{\rm
c}^{2}}{53760\,\gamma^{3/2}}\,\Bigl{[}\frac{\gamma^{1/2}\,h(x/x_{\rm
c})}{[1+\gamma(x/x_{\rm c})^{2}]^{6}}+\frac{3465}{x/x_{\rm
c}}\tan^{-1}\Bigl{(}\frac{\gamma^{1/2}x}{x_{\rm c}}\Bigr{)}\Bigr{]}$ (78)

with $\gamma\simeq 0.091$ and the function $h(y)$ given by

$\displaystyle h(y)$
$\displaystyle=3465\gamma^{5}y^{10}+19635\gamma^{4}y^{8}+45738\gamma^{3}y^{6}$
$\displaystyle+55638\gamma^{2}y^{4}+36685\gamma y^{2}+11895.$ (79)

Note that the quantity $x$ is the dimensionless radius normalized by $r_{\rm
s}$, $x=r/r_{\rm s}$.

Plugging the explicit functional form into Eq. (76), we can repeat the same
analysis as we have done in Sec. III as follows. First notice that the
function $g_{0}(x)$ introduced at Eq. (17) is changed to

$\displaystyle
g_{0}(x)\longrightarrow\frac{\log(1+x)}{x}+\frac{\mathcal{E}}{\alpha}-\frac{\rho_{\rm
c}}{\rho_{\rm s}}\,\tilde{\Psi}_{\rm soliton}(x),$ (80)

where the function $\tilde{\Psi}$ is related to the potential $\Psi_{\rm
soliton}$ through $\Psi_{\rm soliton}=4\pi\,G\,\rho_{\rm c}r_{\rm
s}^{2}\,\tilde{\Psi}_{\rm soliton}$. One important point to note is that the
modification due to soliton potential does not affect the monotonicity of the
function $g_{0}$. That is, the function $g_{0}$ still has a single turning
point, though its location is perturbatively shifted to an inner radius. Thus,
even if the soliton self-gravity is taken into account, the ground-state
eigenfunction which describes the soliton core structure is still given by Eq.
(47), with the function $g_{0}$ replaced with Eq.(80). Hence, taking the limit
of $x\ll 1$, we have the same asymptotic form of the eigenfunction as we
obtained in Sec. III.3 [see Eq. (67)], but the explicit form of the
coefficient $\beta$ is now changed to

$\displaystyle\beta$
$\displaystyle\longrightarrow\frac{\alpha}{6}\Bigl{\\{}1+\frac{\mathcal{E}}{\alpha}-\frac{\rho_{\rm
c}}{\rho_{\rm s}}\tilde{\Psi}_{\rm soliton}(0)\Bigr{\\}},$ (81)

where we have only taken into account the leading-order term proportional to
the parameter $\alpha$. This modification implies that the soliton self-
gravity can change the core size $x_{\rm c}$ from Eq. (69) to

$\displaystyle x_{\rm c}\longrightarrow
p\,\Bigl{[}\frac{6}{\alpha\,\bigl{\\{}1+\mathcal{E}_{1}/\alpha+(\rho_{\rm
c}/\rho_{\rm s})(x_{\rm c}^{2}/14/\gamma)\bigr{\\}}}\Bigr{]}^{1/2},$ (82)

with $p$ given by $0.65$. Here we used the fact that $\tilde{\Psi}_{\rm
soliton}(0)=-x_{\rm c}^{2}/(14\,\gamma)$. In the above, the energy eigenvalue
$\mathcal{E}$ has to be consistently evaluated, taking the soliton potential
into account. Employing the perturbative calculation familiar in quantum
mechanics, this is estimated to give

$\displaystyle\mathcal{E}_{1}/\alpha$
$\displaystyle\simeq\mathcal{E}_{1}^{(0)}/\alpha+\Bigl{(}\frac{\rho_{\rm
c}}{\rho_{\rm s}}\Bigr{)}\int
dx\,\tilde{u}_{1}^{(0)}(x)\,\widetilde{\Psi}_{\rm
soliton}(x)\,\tilde{u}_{1}^{(0)}(x),$ (83)

where the quantity $\mathcal{E}_{1}^{(0)}$ stands for the unperturbed energy
eigenvalue of the ground state, which is obtained analytically from Eq. (39).
The function $\tilde{u}_{1}^{(0)}$ represents the unperturbed ground-state
eigenfunction that is properly normalized.

The new expression of the soliton core radius at Eq. (82) involves the soliton
core radius $x_{\rm c}$, which we evaluate with the unperturbed result at Eq.
(69). As a result, the ratio of the new core radius to the unperturbed
counterpart ignoring the soliton self-gravity, which we denote by
$\mathcal{R}$, is expressed as

$\displaystyle\mathcal{R}\equiv\frac{x_{\rm c}({\rm
w/\,self\,\,gravity})}{x_{\rm c}({\rm
w/o\,\,self\,\,gravity})}=\Biggl{[}1+\frac{\rho_{\rm c}/\rho_{\rm
s}}{1+\mathcal{E}_{1}/\alpha}\,\Bigl{\\{}\frac{x_{\rm c}^{2}}{14\,\gamma}+\int
dx\,\bigl{\\{}\tilde{u}_{1}^{(0)}(x)\bigr{\\}}^{2}\tilde{\Psi}_{\rm
soliton}(x)\Bigr{\\}}\Biggr{]}^{-1/2}.$ (84)

Using the explicit expressions for the quantities $\alpha$, $x_{\rm c}$,
$\rho_{\rm s}$ and adopting the relation for $\rho_{\rm c}$, respectively
given at Eqs. (13), (69), (87) and (65), the above expression is simplified to
give

$\displaystyle\mathcal{R}=\Bigl{[}1+1.73\Bigl{\\{}1-\int
dx\,\bigl{\\{}\tilde{u}_{1}^{(0)}(x)\bigr{\\}}^{2}V_{\rm
soliton}(x)\Bigr{\\}}\Bigr{]}^{-1/2},$ (85)

where the function $V_{\rm soliton}$ is the reduced (minus) soliton potential
defined by $V_{\rm soliton}(x)\equiv-\tilde{\Psi}_{\rm
soliton}(x)/\tilde{\Psi}(0)$. Note that the reduced potential $V_{\rm
soliton}$ is a monotonically decreasing function of $x$, and it goes to $1$ at
$x\to 0$ and $0$ at $x\to\infty$. Thus the integral given above generally
takes a value in between $0$ and $1$. Accordingly the ratio $\mathcal{R}$ is
expected to be smaller than $1$, and is around $\mathcal{R}\sim 0.6-1$888In
Eq. (85), we adopt the soliton core density $\rho_{\rm c}$ found in numerical
simulations, and the depth of the soliton potential relative to the halo
potential is fixed. In this sense, there is no control parameter to change the
soliton self-gravity itself. Rather, in our case, the soliton self-gravity is
determined by the quantities characterizing the halo properties. .

Equation (85) implies that the self gravity of soliton always makes the core
radius small. How it is small actually depends on the parameters, which are
wholly encapsulated in the integral involving the ground-state eigenfunction
and reduced potential. Note that these functions are solely characterized by
the parameter $\alpha$. We thus first show the behavior of the ratio
$\mathcal{R}$ as a function of $\alpha$. The result is plotted in the left
panel of Fig. 10. We then found that the ratio $\mathcal{R}$ is nearly
constant over a wide range of the parameter $\alpha$. As a result, the soliton
self-gravity changes the core size by $20-30$%. Here, for illustrative
purpose, the plotted range of the parameter $\alpha$ was chosen to be
sufficiently wide, $0.1\leq\alpha\leq 10^{5}$. While this range partly
includes the region where our analytical estimation of the soliton core size
become inaccurate (i.e., $\alpha\lesssim 1$), the resultant change in
$\mathcal{R}$ is fairly small, and no notable (sudden) transition appears. In
the right panel, fixing the cosmology and FDM mass to those in the fiducial
setup (see Sec. IV.1), the dependence on the halo mass and concentration
parameter is shown for the case of $z=0$. The result depicted as color scales
shows that the impact of soliton self gravity becomes gradually significant as
decreasing the halo mass and increasing the halo concentration. As a result,
the ratio $\mathcal{R}$ reaches $\sim 0.7$ as shown in the bluish upper
triangular region. Note that at this region, the parameter $\alpha$ also
becomes small and eventually below unity. However, the parameters of this
region are rather extreme, and halos having such a high concentration
parameter can never be realized in numerical simulations with a relevant
cosmological setup.

Hence, we conclude that the soliton self-gravity can be non-negligible for the
setup considered in Sec. IV, but its impact would not be large, suggesting
that our analytical treatment ignoring the soliton self gravity is still
valid, and can give a solid prediction for soliton core structure even from a
quantitative viewpoint. Nevertheless, we note that our estimation given here
relies on the perturbative treatment, and the contribution of the soliton
potential is assumed to be small enough, compared to the halo potential. In
this sense, the actual size of the impact may be quantitatively changed. For
more accurate estimation, a nonperturbative approach by numerically solving
the S-P equation is crucial, including self-consistently the soliton
potential. This is left to our future study.

Figure 10: Left: impact of soliton self-gravity on the core radius. Taking the
backreaction of the soliton self-gravity into account, the estimated core
radius is normalized by the one ignored the self-gravity, which we denote by
$\mathcal{R}$ (see Eq. 85), and is plotted as a function of the dimensionless
parameter $\alpha$. Right: same as left panel, but the result is plotted as
function of the halo mass $M_{\rm h}$ and the concentration parameter $c_{\rm
vir}$, with the ratio $\mathcal{R}$ depicted as color scales. Note that the
mass of FDM is fixed to $m_{\phi}=3\times 10^{-23}$ eV and we follow the same
parameter setup as adopted in Figs. 6–9 (see Sec. IV.1).

##  VI Conclusion

In this paper, we have investigated analytically the structural properties of
the soliton core, the stable dense core formed at the center of virialized
halos in the fuzzy dark matter (FDM) model. Early numerical studies suggest
that the size and mass of such a structure are tightly related to the host
halo properties, given in a power-law form. On the other hand, more recent
studies with high-resolution simulations reveal rich dynamical properties of
the soliton core, also indicating that no tight relation between the soliton
core structure and host halo mass exists, but instead there is a large
scatter. To address these issues, we have presented an analytical description
of the soliton core in a cosmological setup.

A crucial step to derive the soliton core properties is to solve the
Schrödinger-Poisson (S-P) equation in an analytically tractable manner. In our
treatment, we considered the gravitational potential of halo as a smooth
background potential, and ignored the self-gravity of the soliton. With this
setup, the problem to solve the S-P equation is reduced to a linear eigenvalue
problem. Still, no exact solution is known, and to obtain analytically the
eigenfunctions of soliton core, we have applied the Langer’s method of uniform
asymptotic approximation. We found that this method provides an excellent
approximation which agrees well with the numerical solution of the linear
eigenvalue problem, and thus gives an accurate description of the soliton
eigenstates over a wide range of model parameters. Then, from the ground-state
eigenfunction, the key expression characterizing the soliton core radius was
derived [Eqs. (69) or (70) ], and combining the fitting formula for the core
density $\rho_{\rm c}$ estimated from simulation [Eq. (65)], the expression
for the soliton core mass was further obtained [Eq. (72)]. While these
expressions apparently resemble those found in an early work [53], the newly
derived formulas contain an additional factor, which depends not only on the
cosmology and FDM mass but also on quantities characterizing the halo
properties. Based on our new analytical formulas, a quantitative prediction of
the soliton core-halo relations was made, and these were compared with results
of numerical simulations of S-P equation.

Our important findings are summarized as follows:

  * •

The soliton core-halo relations, i.e., $r_{\rm c}$-$M_{\rm h}$ and $M_{\rm
core}$-$M_{\rm h}$ relations, depends sensitively not only on the cosmology
and FDM mass but also on the halo concentration characterized by $c_{\rm
vir}$. In general, a lower concentration halo can host the soliton having a
larger core radius and smaller core mass (see Fig. 6). Comparison of the
predicted core-halo relations with simulation data suggests that the soliton
core prefers halos with moderately small values of concentration parameter,
$c_{\rm vir}\sim\mathcal{O}(1-10)$.

  * •

The predicted core-halo relations plotted as a function of halo concentration
parameter (Fig. 6) shows that for each halo mass, there exist the allowable
core size and core mass for the soliton states, respectively given as the
upper and lower bounds. For fixed cosmological parameters and FDM mass, the
theoretical upper limit on the core radius, $r_{\rm c}^{\rm limit}$, and lower
limit on the core mass, $M_{\rm c}^{\rm limit}$, are simply given as a
function of halo mass, and are scaled as $r_{\rm c}^{\rm limit}\propto M_{\rm
h}^{-1/9}$ and $M_{\rm c}^{\rm limit}\propto M_{\rm h}^{1/9}$ [see Eqs. (74)
and (75) for more explicit expressions]. A general analysis in Appendix C
shows that these limits sensitively depend on the inner slope of halo density
profile [see Eq. (130)].

  * •

Incorporating explicitly the concentration-halo mass (C-M) relation into the
predictions, the resultant core-halo relations generally deviate from a power-
law form. In particular, the C-M relation monotonically decreasing with the
halo mass, as typically seen in the CDM model, makes the core-halo relations
flatter, apparently inconsistent with measured results in simulations (Fig.
9). On the other hand, the nonmonotonic C-M relation like the one indicated in
the FDM model [54] eventually exhibit a power-law feature at the small halo
masses of $M_{\rm h}\sim 10^{8}-10^{9}\,M_{\odot}$ (Fig. 9), leading to a
close agreement with simulation results.

  * •

Overall, the scatter around C-M relation produces a non-negligible amount of
scatter in the soliton core-halo relations. This can explain the diversity
advocated in Ref. [27]. An interesting trend seen in the simulation results
may be that the scatter in the core-halo relations gets reduced as we go to
small halo mass. This is rather consistent with predictions, in which a
suppression of scatter is explained by the fact that the allowable size and
mass of the soliton core are limited.

We have also examined the validity of our treatment ignoring the soliton self-
gravity. A perturbative estimation of its impact on the soliton core radius
indicates that including the soliton self-gravity makes the effective
potential of the S-P equation deep, leading to a shrinking of the core radius.
The degree to which the core radius shrinks is characterized only by a single
parameter, $\alpha$, and adopting the core density $\rho_{\rm c}$ found in
numerical simulations [Eq. (65)], we found that the soliton self-gravity may
change the results by $\sim 20\%$, irrespective of the parameters we
considered. Hence, the soliton self-gravity cannot be completely ignored, but
its impact may not be large, and the predicted core-halo relations ignoring
the self gravity would be still relevant, at least qualitatively, hence
worthwhile for a further detailed comparison with new numerical simulations.
Nevertheless, the validity of this estimation as well as a more proper
description of soliton core have to be investigated, and we leave it to future
work.

Finally, the analytical description presented in this paper would provide a
useful and powerful route to characterize the dynamical properties of soliton
core. As briefly mentioned in Secs. I and II.2, the soliton core seen in
simulations is not strictly stable. Through the wave interference with
granules distributed in a halo, the soliton core undergoes oscillations and
random walks. These dynamical features can be indeed described as a
superposition of the ground and excited states, including those having nonzero
angular momenta [30, 31]. Although previous studies have numerically
constructed such a dynamical state, the analytical construction of the
eigenfunctions, described in Sec. III for $\ell=0$ and Appendix B for
$\ell\neq 0$, would be advantageous to fast realize a dynamical soliton with
less computational cost. Importantly, the present treatment provides a fairly
general way to construct the eigenstates, and can be applied not only to the
soliton in the NFW halo, but also to the one in a generalized halo profile
expressed in an analytical form. Following proposals by Refs. [67, 68], one
can even use our analytical treatment to generate the FDM halo structure
consisting of many granules, each of which resembles the soliton eigenstates
(see also Refs. [69, 70]). Analytically describing the dynamical soliton
states including halos is a very interesting and important subject and would
offer a useful tool to observationally test the FDM model (e.g., Refs. [71,
72] for related works). We will address it near future.

###### Acknowledgements.

We are grateful to Kohei Hayashi and Hei Yin Jowett Chan for providing us
their simulation data of soliton core structures and useful comments, Hiroki
Kawai, Yusuke Manita, Takahiro Nishimichi, Masamune Oguri and Takuya Takahashi
for discussions. We would like to also thank Mona Dentler and David J. Marsh
for correcting the concentration-mass relation of fuzzy dark matter model used
in this paper. This work was supported by Grant-in-Aid for JSPS Fellows No.
17J10553 (S. S.) and in part by MEXT/JSPS KAKENHI Grant Numbers JP20H05861 and
JP21H01081 (A. T. ). A. T. acknowledges the support from JST AIP Acceleration
Research Grant Number JP20317829. S. S. is supported by JSPS Overseas Research
Fellowships. Numerical computation was carried out partly at the Yukawa
Institute Computer Facility.

##  Appendix A Halo density profile and mass concentration

In this Appendix, we summarize important characteristics of the dark matter
halos used in the main text.

###  A.1 Characteristics of halo density profile

The dark matter halo, a virialized self-gravitating dark matter clump, is
known to have a characteristic structure in density. In particular, a
spherical mean of the density profile measured in cosmological $N$-body
simulations is known to be quantitatively described by the Navarro-Frenk-White
(NFW) profile [46, 47]. The comoving density of this profile is expressed as
[Eq. (7)]

$\displaystyle\rho_{\rm NFW}(r)=\frac{\rho_{\rm s}}{(r/r_{\rm s})(1+r/r_{\rm
s})^{2}},$ (86)

where the radius $r_{\rm s}$ describes the transition scale where the
logarithmic slope of the density profile changes from $-3$ to $-1$ as we
approach the halo center. The density $\rho_{\rm s}$ represents the
characteristic density at $r=r_{\rm s}$, and it is also expressed as
follows999In the definition of $\rho_{\rm NFW}$ at Eq. (7), a factor of
$1/a^{3}$ is factored out from the density profile. As a result, the
dependence of $\rho_{\rm s}$ on the scale factor is somewhat different from
the one in the literature.:

$\displaystyle\rho_{\rm s}=\frac{\Delta_{\rm vir}\,\rho_{\rm
m,0}}{3}\,\frac{c_{\rm vir}^{3}}{\ln(1+c_{\rm vir})-c_{\rm vir}/(1+c_{\rm
vir})}.$ (87)

Here, we introduced the concentration parameter $c_{\rm vir}$ defined by
$c_{\rm vir}\equiv r_{\rm vir}/r_{\rm s}$ with the radius $r_{\rm vir}$ being
the virial radius. The virial radius is linked to the halo mass through

$\displaystyle r_{\rm vir}=\Bigl{(}\frac{3\,M_{\rm h}}{4\pi\,\Delta_{\rm
vir}\rho_{\rm m,0}}\Bigr{)}^{1/3}.$ (88)

The quantity $\Delta_{\rm vir}$ is the virial overdensity, which characterizes
the mean overdensity of halos within the virial radius. While this is
precisely described by the spherical collapse model [73, 74, 75], a simple and
accurate fitting formula, which can be used in a flat universe with
cosmological constant, is known [76] (see also Ref. [77]):

$\displaystyle\Delta_{\rm vir}(z)$
$\displaystyle=\frac{18\pi^{2}+82\,\\{\Omega_{\rm m}(z)-1\\}-39\\{\Omega_{\rm
m}(z)-1\\}^{2}}{\Omega_{\rm m}(z)}.$ (89)

Here the quantity $\Omega_{\rm m}(z)$ is the matter density parameter at
redshift $z$, given by

$\displaystyle\Omega_{\rm m}(z)$ $\displaystyle=\frac{(1+z)^{3}\,\Omega_{\rm
m,0}}{(1+z)^{3}\Omega_{\rm m,0}+\Omega_{\Lambda,0}}$ (90)

with the subscript $0$ indicating the quantities given at the present time.

###  A.2 Concentration-mass relations

For a quantitative characterization of the halo density profile, the key
parameter is the concentration parameter $c_{\rm vir}$, which characterizes
the size and mass concentration of each halo. In the context of the CDM
cosmology, there have been numerous works to quantify its cosmological and
halo mass dependence based on the cosmological $N$-body simulations, and
accurately calibrated models have been exploited. One of such representative
models is given by Ref. [55]. The concentration parameter of this model, which
we denote by $c_{\rm vir}^{\rm B}$, is expressed as follows:

$\displaystyle c_{\rm vir}^{\rm B}(M_{\rm h},P_{\rm CDM})=A\frac{1+z_{\rm
coll}(M_{\rm h},\,P_{\rm CDM})}{1+z},$ (91)

with the constant $A$ chosen as $3.13$ [78]. The redshift $z_{\rm coll}$
characterizes the collapse time for halos with mass $M_{\rm h}$ in the CDM
cosmological model with the linear power spectrum $P_{\rm CDM}$. It is
determined through

$\displaystyle\frac{D(z_{\rm coll})}{D(z)}\,\sigma(f_{\rm coll}\,M_{\rm
h},\,P_{\rm CDM})=\delta_{\rm c},$ (92)

where the function $D$ is the linear growth factor, and $\sigma$ is the root-
mean-square amplitude of the linear density field smoothed with the top-hat
filter function over the radius determined by the halo mass $M_{\rm h}$, i.e.,
$\\{3\,M_{\rm h}/(4\pi\,\rho_{\rm m,0}\\}^{1/3}$. The quantity $\delta_{\rm
c}$ is the critical threshold for the linear density contrast, which we adopt
the fitting form given by Ref. [78]: $\delta_{\rm
c}=1.59+0.0314\ln\sigma_{8}(z)$ with $\sigma_{8}$ being the same rms
fluctuation as we defined above, but smoothed at the specific radius of
$8\,h^{-1}$ Mpc. Here, the constant $f_{\rm coll}$ is set to $0.01$.

On the other hand, in the FDM cosmology, there is a little numerical work to
calibrate the C-M relation based on the simulations of S-P equation. Still,
there are several analytical works to model it in order to assess the small-
scale problem as well as to constrain the FDM mass from weak lensing
observations [79, 80, 54]. Here, we adopt the model presented by Ref. [54] and
examine the core-halo relation based on our analytical treatment. The
concentration parameter of this model is expressed as a correction to the
model of Ref. [55] in the following form:

$\displaystyle c_{\rm vir}^{\rm F}(M_{\rm h})=c_{\rm vir}^{\rm B}(M_{\rm
h},P_{\rm FDM})\Delta_{\rm c}^{\rm FDM},$ (93)

where the concentration parameter $c_{\rm vir}^{\rm B}$ is obtained from Eq.
(91) with the linear power spectrum $P_{\rm FDM}$ given in the FDM
cosmological model [see e.g., Eqs. (8) with (9) in Ref. [6]]. The correction
term $\Delta_{\rm c}^{\rm FDM}$ is defined by

$\displaystyle\Delta_{\rm c}^{\rm
FDM}=\Bigl{\\{}1+\frac{M_{0}^{\rm(c)}}{M_{\rm
h}}\Bigr{\\}}^{-\gamma_{0}}\,\Bigl{\\{}1+\gamma_{1}\,\frac{f_{\rm
coll}\,M_{0}^{\rm(c)}}{M_{\rm h}}\Bigr{\\}}^{-\gamma_{2}},$ (94)

with the coefficient $\gamma_{1}$ defined by $\gamma_{1}=15$ and indices
$\gamma_{0}$ and $\gamma_{2}$, respectively, given by $\gamma_{0}=d\ln c_{\rm
vir}^{\rm B}/d\ln M_{\rm h}|_{4\,M_{0}^{\rm(c)}}$ and $\gamma_{2}=0.3$. Here,
the characteristic mass $M_{0}^{\rm(c)}$ is related to the scaling mass
$M_{0}^{\rm(n)}$ through $M_{0}^{\rm(c)}=M_{0}^{\rm(n)}/f_{\rm coll}$, and
$M_{0}^{\rm(n)}$ is defined by

$\displaystyle M_{0}^{\rm(n)}=1.6\times
10^{10}\,\,M_{\odot}\,\Bigl{(}\frac{m_{\phi}}{10^{-22}\,{\rm
eV}}\Bigr{)}^{-4/3}.$ (95)

##  Appendix B Analytical solutions $\ell\neq 0$

In this Appendix, extending the treatment in Sec. III.1, we construct
analytically the approximate solutions of the S-P equation given at Eq. (12)
in the case of $\ell\neq 0$.

###  B.1 Constructing analytical Eigenstates

Similarly to the $\ell=0$ case in Sec. III.1, we rewrite Eq. (12) with the
normal form, keeping the angular momentum nonzero. Again introducing the new
function $\tilde{u}_{n\ell}\equiv x\,u_{n\ell}$, the equation for the radial
wave function is given in the following form:

$\displaystyle\frac{d^{2}\tilde{u}_{n\ell}(x)}{dx^{2}}+\alpha\,g_{\ell}(x)\tilde{u}_{n\ell}(x)=0\,;$
$\displaystyle\qquad\quad
g_{\ell}(x)\equiv-\frac{\ell(\ell+1)}{\alpha\,x^{2}}+\frac{\log(1+x)}{x}+\frac{\mathcal{E}}{\alpha}.$
(96)

Imposing the same boundary conditions as given at Eq. (18), we construct
analytically the eigenfunction and eigenvalues of this system. Following the
same steps as we examined in the $\ell=0$ case, we apply the Liouville-Green
transformation $(x,\tilde{u}_{n\ell})\to(z,v_{n\ell})$ in Eq. (20). Equation
(96) is then rewritten with

$\displaystyle\frac{d^{2}v_{n\ell}(z)}{dz^{2}}+\Bigl{[}\alpha\frac{g_{\ell}(x)}{\\{p^{\prime}(x)\\}^{2}}+\delta\Bigr{]}v_{n\ell}(z)=0,$
(97)

where the quantity $\delta$ is defined by Eq. (22).

Comparing the above equation with Eq. (21) in the $\ell=0$ case, an important
difference appears in the function $g_{\ell}$, which has now the two turning
points. That is, the function becomes zero at $x=x_{1}$ and $x_{2}$, and is
positive and negative, respectively, at $x_{1}\leq x\leq x_{2}$ and $x<x_{1}$
or $x>x_{2}$. In order to obtain a uniformly valid approximation, the function
$z=p(x)$ has to be chosen in such a way that it behaves like the function
$g_{\ell}$, having two zero-crossing points. Here, we specify the functional
form of $z$ in the following way [40, 41]:

$\displaystyle\frac{g_{\ell}(x)}{\\{p^{\prime}(x)\\}^{2}}$
$\displaystyle=\beta^{2}-z^{2}$ $\displaystyle=\beta^{2}-\\{p(x)\\}^{2},$ (98)

where the quantity $\beta$ is a constant which will be given later. Equation
(98) is the first-order differential equation for the function $p$, which is
analytically integrated to give the following relations:

$\displaystyle\int_{x_{1}}^{x}\sqrt{g_{\ell}(x^{\prime})}\,dx^{\prime}$
$\displaystyle=\int_{-\beta}^{p}\sqrt{\beta^{2}-\tau^{2}}\,d\tau$
$\displaystyle=\frac{1}{2}\beta^{2}\,\arccos\Bigl{(}-\frac{p}{\beta}\Bigr{)}+\frac{p}{2}\,\sqrt{\beta^{2}-p^{2}}$
(99)

for $x_{1}<x<x_{2}$,

$\displaystyle\int_{x}^{x_{1}}\sqrt{-g_{\ell}(x^{\prime})}\,dx^{\prime}$
$\displaystyle=\int_{p}^{-\beta}\sqrt{\tau^{2}-\beta^{2}}\,d\tau$
$\displaystyle=-\frac{1}{2}\Bigl{\\{}\beta^{2}\,\cosh^{-1}\Bigl{(}-\frac{p}{\beta}\Bigr{)}+p\,\sqrt{p^{2}-\beta^{2}}\Bigr{\\}}$
(100)

for $x\leq x_{1}$, and

$\displaystyle\int_{x_{2}}^{x}\sqrt{-g_{\ell}(x^{\prime})}\,dx^{\prime}$
$\displaystyle=\int_{\beta}^{p}\sqrt{\tau^{2}-\beta^{2}}\,d\tau$
$\displaystyle=-\frac{1}{2}\Bigl{\\{}\beta^{2}\,\cosh^{-1}\Bigl{(}\frac{p}{\beta}\Bigr{)}-p\,\sqrt{p^{2}-\beta^{2}}\Bigr{\\}}$
(101)

for $x\geq x_{2}$.

From Eqs. (99)-(101) given above, one can construct the function $p$ expressed
in terms of the variable $x$ as follows. Provided the functional form of
$g_{\ell}$, the integrals at the left-hand side of these equations are
numerically performed, and the results are given as a function of $x$. On the
other hand, the functions at the right-hand side of Eqs. (99)-(101) are
evaluated as a function of $p$. That is, Eqs. (99)-(101) are symbolically
written in the form as $A(x)=B(p)$. Hence, solving numerically this implicit
equation for $p$, i.e., $p=B^{-1}(A(x))$, the functional form of $p(x)$ is
determined. In the above, the constant $\beta$ is arbitrary, but we shall
specify it in such a way that with the Liouville-Green transformation, the
turning points $x_{1}$ and $x_{2}$ are mapped to $p=-\beta$ and $p=\beta$,
respectively. From Eq. (99), this is equivalent to imposing the following
condition:

$\displaystyle\beta^{2}=\frac{2}{\pi}\int_{x_{1}}^{x_{2}}\sqrt{g_{\ell}(x^{\prime})}\,dx^{\prime}.$
(102)

With the specific choice at Eq. (98), Eq. (97) is now recast as

$\displaystyle\frac{d^{2}v_{n\ell}(z)}{dz^{2}}+\alpha\bigl{(}\beta^{2}-z^{2}\bigr{)}v_{n\ell}(z)=-\delta(z)v_{n\ell}(z).$
(103)

Ignoring the term at right-hand side, the general solution of this equation is
described by the parabolic cylinder functions $U$ and $\overline{U}$ [41,
40]101010The function $U$ given here is exactly the same one as defined in
Ref. [50]. Also, it is related to the Weber’s parabolic cylinder function
$D_{\lambda}$ through $D_{\lambda}(x)=U(-\lambda-1/2,x)$. On the other hand,
the function $\overline{U}$ slightly differ from the $V$ in Ref. [50] by the
factor of $\Gamma(1/2-b)$, i.e., $U(b,x)=\Gamma(1/2-b)\,V(a,x)$.:

$\displaystyle v_{n\ell}(z)$
$\displaystyle=\tilde{c}_{1}\,U\Bigl{(}-\frac{1}{2}\alpha^{1/2}\beta^{2},\,\sqrt{2\alpha^{1/2}}\,z\Bigr{)}$
$\displaystyle\qquad+\tilde{c}_{2}\,\overline{U}\Bigl{(}-\frac{1}{2}\alpha^{1/2}\beta^{2},\,\sqrt{2\alpha^{1/2}}\,z\Bigr{)}$
(104)

or

$\displaystyle\tilde{u}_{n\ell}(x)$
$\displaystyle=\Biggl{[}\frac{|\beta^{2}-\\{z(x)\\}^{2}|}{|g_{\ell}(x)|}\Biggr{]}^{1/4}$
$\displaystyle\times\Biggl{\\{}c_{1}\,U\Bigl{(}-\frac{1}{2}\alpha^{1/2}\beta^{2},\,\sqrt{2\alpha^{1/2}}\,z(x)\Bigr{)}$
$\displaystyle\quad\quad+c_{2}\,\overline{U}\Bigl{(}-\frac{1}{2}\alpha^{1/2}\beta^{2},\,\sqrt{2\alpha^{1/2}}\,z(x)\Bigr{)}\Biggr{\\}},$
(105)

where the coefficients $\tilde{c}_{i}$ and $c_{i}$ are the integration
constants.

Note that the parabolic cylinder functions have the following asymptotic
behaviors:

$\displaystyle U(b,\zeta)$
$\displaystyle\stackrel{{\scriptstyle\zeta\to\infty}}{{\sim}}\frac{1}{\zeta^{b+1/2}}\,e^{-\zeta^{2}/4},$
$\displaystyle\overline{U}(b,\zeta)$
$\displaystyle\stackrel{{\scriptstyle\zeta\to\infty}}{{\sim}}\frac{2^{1/2}\,\zeta^{b-1/2}}{\pi^{1/2}}\,\Gamma\Bigl{(}\frac{1}{2}-b\Bigr{)}\,e^{\zeta^{2}/4},$
(106)

and

$\displaystyle
U(b,\zeta)\stackrel{{\scriptstyle\zeta\to-\infty}}{{\sim}}\frac{(2\pi)^{1/2}(-\zeta)^{b-1/2}}{\Gamma(1/2+b)}\,e^{\zeta^{2}/4}.$
(107)

Let us recall from Eqs. (99)–(102) that the function $p(x)$ or $z(x)$ is a
monotonically increasing function of $x$, and varies from negative to
positive. Thus, the quantity $p$ approaches a negative constant value at $x=0$
and goes to infinity in the limit of $x\to\infty$. From Eq. (106), this
indicates that the second term in Eqs. (104) or (105) diverge in the limit of
$x\to+\infty$. Hence, to fulfill the boundary condition $\tilde{u}(\infty)=0$,
the coefficients $\overline{c}_{2}$ and $c_{2}$ must be zero.

On the other hand, another boundary condition $\tilde{u}(0)=0$ yields

$\displaystyle
U\Bigl{(}-\frac{1}{2}\alpha^{1/2}\beta^{2},\,\sqrt{2\alpha^{1/2}}\,z(0)\Bigr{)}=0,$
(108)

where we used the fact that the functions $g_{\ell}(x)$ and
$|\beta^{2}-\\{z(x)\\}^{2}|$ are nonzero at $x=0$. Equation (108) is the
transcendental equation for the eigenvalue $\varepsilon/\alpha$. That is,
given the parameters $\alpha$ and $\ell$, the turning points $x_{1,2}$ are
determined, and the constant $\beta$ and quantity $z(0)$ are computed, with
both of them given as a function of $\mathcal{E}/\alpha$. Then, substituting
them into Eq. (108), solving it yields a discrete set of eigenvalues
$\mathcal{E}/\alpha$.

Once the eigenvalues are determined from Eq. (108), the functional form of the
eigenfunctions for $\ell\neq 0$ are specified in an unambiguous way except for
the normalization factor. To sum up, the (unnormalized) eigenfunctions for
$\ell\neq 0$ that are constructed analytically are expressed as

$\displaystyle\tilde{u}_{n\ell}(x)$
$\displaystyle=x\,u_{n\ell}(x)=\Bigl{[}\frac{|\beta^{2}-\\{z(x)\\}^{2}|}{|g_{\ell}(x)|}\Bigr{]}^{1/4}\,U\Bigl{(}-\frac{1}{2}\alpha^{1/2}\beta^{2},\,\sqrt{2\alpha^{1/2}}\,z(x)\Bigr{)},\,\,\,$
$\displaystyle\beta^{2}=\frac{2}{\pi}\int_{x_{1}}^{x_{2}}\sqrt{g_{\ell}(x^{\prime})}\,dx^{\prime}$
(109)

with the function $g_{\ell}$ given at Eq. (96). The function $p=z(x)$ is
obtained by Eqs. (99)–(101).

Table 2: Lowest five eigenvalues normalized by $\alpha$ in the cases of the nonzero angular momenta $\ell=2$, $5$ and $10$, setting the parameter $\alpha$ to $\alpha=10^{3}$ (see Fig. 12). Analytical results obtained from Eq. (108) are summarized. For comparison, parentheses show the results by solving numerically the matrix eigenvalue problem, adopting the boundaries $x_{1}=0$ and $x_{n}=5$ with the number of grids $n=10^{4}$.  $\mathcal{E}_{n\ell}/\alpha$ | $\ell=2$ | $\ell=5$ | $\ell=10$
---|---|---|---
$n=1$ | -0.77454 (-0.77278) | -0.67810 (-0.67756) | -0.57018 (-0.56999)
$n=2$ | -0.71757 (-0.71620) | -0.63909 (-0.63863) | -0.54483 (-0.54466)
$n=3$ | -0.67142 (-0.67029) | -0.60531 (-0.60491) | -0.52189 (-0.52174)
$n=4$ | -0.63259 (-0.63163) | -0.57561 (-0.57526) | -0.50100 (-0.50087)
$n=5$ | -0.59911 (-0.59828) | -0.54916 (-0.54885) | -0.48187 (-0.48175)

Figure 11: Eigenfunctions for the radial wave function having nonzero angular
momentum, $\tilde{u}_{n\ell}(x)$. From left to right, fixing the parameter
$\alpha$ to $10^{3}$, eigenfunctions for the lowest five eigenstates are shown
for $\ell=2$, $5$, and $10$, respectively. In the upper panels, thin solid
lines are the eigenfunctions constructed analytically, while the thick dashed
lines are the results obtained by solving the matrix eigenvalue problem
numerically. Both of them are properly normalized. Bottom panels show the
fractional difference between analytical and numerical results, defined by
$(\tilde{u}_{n\ell,{\rm analytical}}-\tilde{u}_{n\ell,{\rm
numerical}})/\tilde{u}_{n\ell,{\rm numerical}}$ (see Fig. 3). Note that to
avoid the division by zero, we stop plotting the results when the amplitude of
the wave functions becomes smaller than $10^{-6}$ in the inner and outer
regions.

Figure 12:  Same as Fig. 12, but we plot the square of the normalized wave
function, $u_{n\ell}(x)=x\,\tilde{u}_{n\ell}(x)$, i.e., $|u_{n\ell}(x)|^{2}$.
Line and color types shown hare are the same as those used in Fig. 12. In each
panel, the gray dashed lines are the asymptotic slope of the wave function at
$x\lesssim 1$ predicted by the WKB approximation [see Eq. (115)].

###  B.2 Comparison with numerical solutions

Let us compare the analytical eigenstates for $\ell\neq 0$ to the numerical
results obtained by solving the matrix eigenvalue problem. For illustration,
we below focus on the specific angular momenta of $\ell=2$, $5$, and $10$.

Table 2 presents the lowest five eigenvalues for non-zero angular momenta,
fixing the dimensionless parameter $\alpha$ to $\alpha=10^{3}$. In Figs. 12
and 12, the corresponding (normalized) eigenfunctions $\tilde{u}_{n\ell}$ and
the square of eigenfunctions, $|u_{n\ell}|^{2}$, are, respectively, plotted as
a function of $x=r/r_{s}$. Line types and colors are the same as in Fig. 3.
Note that in obtaining the numerical results, we solve Eq. (63), with the
potential $V(x)$ being replaced with
$V(x)=-\alpha\,\log(1+x)/x+\ell(\ell+1)/x^{2}$.

Overall, we see an excellent agreement between analytical and numerical
results, similarly to the $\ell=0$ case. Although a closer look at the
fractional difference shown in the bottom panels of Fig. 12 reveals that the
eigenfunctions exhibit discrepancies, these become non-negligible only when
the eigenfunction becomes close to zero, and hence their net impact is
considered to be small. Hence, we conclude that the analytically constructed
eigenstates are accurate enough even for the nonzero angular momenta, and
would be useful to characterize quantitatively the dynamical properties of the
soliton core.

Finally, one notable behavior in the square of the eigenfunction,
corresponding to the soliton density profile, is the characteristic peak
structure shown in Fig. 12. As increasing $\ell$, the density peak generically
appears at a larger radius, with its width narrower. In particular, the inner
structure exhibits a power-law behavior, which we find to be scaled as
$|u_{n\ell}|^{2}\propto x^{2\sqrt{\ell(\ell+1)}-1}$ at $x\ll 1$. This is
regardless of whether it is in the ground or excited state (see gray dashed
lines). In Appendix B.3, we derive this power-law form explicitly based on the
WKB approximation.

###  B.3 Asymptotic behaviors near origin

In this Appendix, we discuss the asymptotic behavior of the eigenfunction
$\tilde{u}_{n\ell}$ near the halo center, and derive analytically the power-
law form of the density profile as shown in Fig. 12. In doing so, rather than
using Eq. (109), we apply the WKB approximation and obtain an analytical
expression for the eigenfunction valid near the origin $(x\approx 0)$, where
the function $g_{\ell}(x)$ is supposed to be negative. The (unnormalized) WKB
solution is then given in the form as

$\displaystyle\tilde{u}_{n\ell}^{\rm
WKB}(x)\simeq\frac{1}{\\{-g_{\ell}(x)\\}^{1/4}}\,$
$\displaystyle\quad\times\Biggl{[}d_{1}\,\exp\Bigl{\\{}\alpha^{1/2}\int_{x}^{x_{1}}\sqrt{-g_{\ell}(x^{\prime})}dx^{\prime}\Bigr{\\}}$
$\displaystyle\quad+d_{2}\,\exp\Bigl{\\{}-\alpha^{1/2}\int_{x}^{x_{1}}\sqrt{-g_{\ell}(x^{\prime})}dx^{\prime}\Bigr{\\}}\Biggr{]}$
(110)

with the upper bound of the integral, $x_{1}$, being the first turning point,
$g_{\ell}(x_{1})=0$. In the above, the coefficients $d_{1}$ and $d_{2}$ are
not arbitrary, but are related with each other through the inner boundary
condition, $\tilde{u}_{n\ell}(0)=0$. This relation is rewritten in the
following form:

$\displaystyle\frac{d_{1}}{d_{2}}=-\exp\Bigl{[}-2\alpha^{1/2}\int^{x_{1}}_{0}\,\sqrt{-g_{\ell}(x^{\prime})}\,dx^{\prime}\Bigr{]}.$
(111)

Recalling that near the origin, the function $g_{\ell}$ is approximated as
$g_{\ell}\simeq\ell(\ell+1)/x^{2}/\alpha$, the integral in the exponent above,
$\int^{x_{1}}_{0}\,\sqrt{-g_{\ell}(x^{\prime})}\,dx^{\prime}$, is shown to
exhibit a logarithmic divergence. To be precise, we may replace the lower
bound of this integral with a nonzero small and positive value of $\eta$, and
perform the integration. Assuming $x_{1}\ll 1$, this leads to
$\sqrt{\ell(\ell+1)/\alpha}\log(x_{1}/\eta)$. Taking the limit $\eta\to 0$,
the integral thus diverges positively. Note that the logarithmic divergence
generically appears regardless of whether $x_{1}$ is small or not. In other
words, the right-hand side of Eq. (111) becomes vanishing, and hence the
coefficient $d_{1}$ must be zero. Equation (110) is recast as

$\displaystyle\tilde{u}_{n\ell}^{\rm
WKB}(x)=\frac{d_{2}}{\\{-g_{\ell}(x)\\}^{1/4}}\,\,\exp\Bigl{\\{}-\alpha^{1/2}\int_{x}^{x_{1}}\sqrt{-g_{\ell}(x^{\prime})}dx^{\prime}\Bigr{\\}}.$
(112)

With the wave function in Eq. (112), we derive an asymptotic power-law form at
$x\ll 1$. To obtain an expression without assuming a specific value of
$x_{1}$, we consider a sufficiently small but still finite value of
$x_{0}\,(\ll x_{1})$, below which the function $g_{\ell}(x)$ is approximately
described by $g_{\ell}(x)\simeq\ell(\ell+1)/x^{2}/\alpha$. Focusing on the
region of $x\lesssim x_{0}$, the integral in the exponent of Eq. (112) is
evaluated to give

$\displaystyle\int_{x}^{x_{1}}\sqrt{-g_{\ell}(x^{\prime})}\,dx^{\prime}$
$\displaystyle\quad=\int_{x_{0}}^{x_{1}}\sqrt{-g_{\ell}(x^{\prime})}\,dx^{\prime}+\frac{\sqrt{\ell(\ell+1)}}{\alpha^{1/2}}\,\log\Bigl{(}\frac{x_{0}}{x}\Bigr{)}.$
(113)

Plugging this into Eq. (112), the WKB solution near the origin $(x\lesssim
x_{0})$ is rewritten in the following form:

$\displaystyle\tilde{u}_{n\ell}^{\rm WKB}(x)$ $\displaystyle\simeq
d_{2}\,\frac{\alpha^{1/4}\,}{\\{\ell(\ell+1)\\}^{1/4}}\,x^{1/2}\,\Bigl{(}\frac{x}{x_{0}}\Bigr{)}^{\sqrt{\ell(\ell+1)}}\,\exp\Bigl{\\{}-\alpha^{1/2}\int_{x_{0}}^{x_{1}}\sqrt{-g_{\ell}(x^{\prime})}dx^{\prime}\Bigr{\\}}.$
(114)

This shows that the density profile for the soliton state, proportional to
$|u_{n\ell}|^{2}$, exhibits a power-law behavior near the origin and it is
proportional to

$\displaystyle|u_{n\ell}^{\rm
WKB}(x)|^{2}=\Bigl{|}\frac{\tilde{u}_{n\ell}^{\rm
WKB}(x)}{x}\Bigr{|}^{2}\propto x^{2\sqrt{\ell(\ell+1)}-1},\quad(x\ll x_{0})$
(115)

Since the above expression has no explicit dependence on $n$, this scaling
generally holds for both ground and excited states. Note cautiously that in
deriving Eq. (115), we assumed nonzero $\ell$ and used the asymptotic behavior
of $g_{\ell}$. In this respect, Eqs. (114) and (115) are valid only for
$\ell\neq 0$ and cannot be applied to the $\ell=0$ case.

##  Appendix C Soliton core size in the limit of $c_{\rm vir}\to 0$

In this Appendix, we discuss the dependence of the soliton core size on the
halo density profile. To be specific, we investigate how the upper bound on
the core radius, which is obtained by taking the limit $c_{\rm vir}\to 0$,
varies with the inner slope of density profile.

Consider a class of spherically symmetric halos having a characteristic scale
$r_{\rm s}$, below which the radial profile asymptotically follows a power-law
form. One can write such a profile formally as

$\displaystyle\rho(x)=\rho_{\rm s}\,D(x),\,\,\,x=r/r_{\rm s}$ (116)

with the dimensionless function $D$ asymptotically behaving as

$\displaystyle D(x)\stackrel{{\scriptstyle x\ll 1}}{{\simeq}}x^{-s}.$ (117)

Note that Eq. (116) is the comoving density, and a factor of $1/a^{3}$ is
factorized in the same way as defined in Eq. (7). Thus, for the NFW profile,
$\rho_{\rm s}$ is given by Eq. (87) and $D$ becomes $D(x)=1/x/(1+x)^{2}$. The
mass and potential profiles of this model are, respectively, given by

$\displaystyle M(x)$ $\displaystyle=4\pi\,\rho_{\rm s}\,r_{\rm
s}^{3}\int_{0}^{x}dy\,y^{2}D(y),$ (118) $\displaystyle\Psi(x)$
$\displaystyle=-\frac{4\pi\,G\,\rho_{\rm s}\,r_{\rm
s}^{2}}{a}\int_{x}^{\infty}\frac{dy}{y^{2}}\int_{0}^{y}dw\,w^{2}\,D(w)$
$\displaystyle\equiv-\frac{4\pi\,G\,\rho_{\rm s}\,r_{\rm
s}^{2}}{a}\,\tilde{\Psi}(x).$ (119)

Below, we shall consider the case that both the mass and potential are finite
at the center. This implies the inner slope $s$ has to be less than $2$. In
such a case, the dimensionless potential $\tilde{\Psi}$ at $x\ll 1$ is
generally written in the following form:

$\displaystyle\tilde{\Psi}(x)\simeq\tilde{\Psi}(0)-\frac{x^{2-s}}{(3-s)(2-s)},\quad(s<2).$
(120)

With the halo mass defined by $M_{\rm h}=(4\pi/3)\rho_{\rm m,0}\Delta_{\rm
vir}r_{\rm vir}^{3}$ at $r=r_{\rm vir}$, the characteristic density $\rho_{\rm
s}$ and the scale radius $r_{\rm s}$ are expressed in terms of the
concentration parameter $c_{\rm vir}=r_{\rm vir}/r_{\rm s}$ as follows:

$\displaystyle\rho_{\rm s}$ $\displaystyle=\frac{\rho_{\rm m,0}\Delta_{\rm
vir}}{3}\,c_{\rm vir}^{3}\Bigl{\\{}\int_{0}^{c_{\rm
vir}}y^{2}D(y)dy\Bigr{\\}}^{-1},$ $\displaystyle\simeq\frac{\rho_{\rm
m,0}\Delta_{\rm vir}}{3}(3-s)\,c_{\rm vir}^{s},$ (121) $\displaystyle r_{\rm
s}$ $\displaystyle=c_{\rm vir}^{-1}\,\Bigl{(}\frac{3}{4\pi}\frac{M_{\rm
h}}{\rho_{\rm m,0}\Delta_{\rm vir}}\Bigr{)}^{1/3},$ (122)

where the last line of Eq. (121) is valid for $c_{\rm vir}\ll 1$.

Provided basic quantities to characterize the halos, we estimate the soliton
core size in the limit of $c_{\rm vir}\to 0$. As we have discussed in Sec.
III, the S-P equation ignoring the soliton self-gravity is reduced to a linear
eigenvalue problem, and focusing on $\ell=0$, the governing equation is recast
in the form similar to Eq. (17), but the halo potential is replaced with Eq.
(119). Nevertheless, this modification does not affect the monotonicity of the
halo potential, and we can follow the same steps as in the case of the NFW
profile to analytically construct the eigenstates. At the end, we obtain the
same analytical expression as given at Eq. (47), from which the expression of
soliton core size can be derived based on the asymptotic properties of the
Airy function (see Sec. III.3). We have111111In deriving Eq. (123), we
consider only the leading-order term proportional to $\alpha$, and ignore
other contributions. :

$\displaystyle x_{\rm c}=p\sqrt{\frac{6}{\alpha\,g_{0}(0)}}\,;\quad p=0.65,$
(123)

with the function $g_{0}(x)$ given by

$\displaystyle g_{0}(x)=\tilde{\Psi}(x)+\mathcal{E}/\alpha.$ (124)

The above expression is regarded as a generalized formula for the core radius.
From this, the comoving core radius $r_{\rm c}$ is written as

$\displaystyle r_{\rm c}$ $\displaystyle=x_{\rm c}\,r_{\rm s}$
$\displaystyle=p\sqrt{\frac{6}{\bigl{\\{}\tilde{\Psi}(0)+\mathcal{E}/\alpha\bigr{\\}}}}\,\frac{r_{\rm
s}}{\alpha^{1/2}}$ (125)

In the above, the parameter $\alpha$ is defined similarly as in the case of
the NFW profile, i.e., $\alpha=8\pi\,G\,m_{\phi}^{2}\,\rho_{\rm s}\,r_{\rm
s}^{4}\,a$ [see Eq. (11)]. Importantly, however, its dependence on the
concentration parameter differs generally from the one given in Eq. (13). To
explicitly see this, we use Eqs. (121) and (122) to rewrite the expression of
$\alpha$ with

$\displaystyle\alpha$ $\displaystyle=8\pi\,G\,m_{\phi}^{2}\,a\frac{\rho_{\rm
m,0}\Delta_{\rm vir}}{3}\Bigl{(}\frac{3}{4\pi}\frac{M_{\rm h}}{\rho_{\rm
m,0}\Delta_{\rm vir}}\Bigr{)}^{4/3}$
$\displaystyle\quad\qquad\times\Bigl{\\{}c_{\rm vir}\int_{0}^{c_{\rm
vir}}y^{2}\,D(y)dy\Bigr{\\}}^{-1}$

Thus, taking the limit $c_{\rm vir}\ll 1$, we obtain

$\displaystyle\alpha\simeq\Bigl{(}\frac{6}{\pi}\Bigr{)}^{1/3}G\,m_{\phi}^{2}\,a\Bigl{(}\frac{M_{\rm
h}^{4}}{\rho_{\rm m,0}\Delta_{\rm vir}}\Bigr{)}^{1/3}(3-s)\,c_{\rm
vir}^{s-4}.$ (126)

Figure 13: Behavior of the function $\mathcal{T}(s)$. The result normalized by
the one evaluated at $s=1$ is plotted as a function of the inner slope of the
density profile, $s$.

In Eq. (125), the expression of the soliton core size also involves the factor
$\tilde{\Psi}(0)+\mathcal{E}/\alpha$, which is given as a function of
$\alpha$. To express its dependence on $\alpha$ analytically, we make use of
the approximation in the limit of $\alpha\gg 1$. As we saw in Eq. (126), this
limit is ensured if we also take the limit $c_{\rm vir}\to 0$ for the inner
slope of the halo profile, $s<3$. Then, the key expression is Eq. (41).
Similarly to what we did in Sec. III.1, we use the expression for
$\tilde{\Psi}$ valid near the center, i.e., Eq. (120). Recalling that the
function $g_{0}$ vanishes at the turning point $x_{\rm
c}\simeq[(3-s)(2-s)\\{\tilde{\Psi}(0)+\mathcal{E}/\alpha\\}]^{1/(2-s)}$ and
assuming $x_{\rm c}\ll 1$, we obtain

$\displaystyle\alpha^{1/2}\int_{0}^{x_{c}}\sqrt{g_{0}(x^{\prime})}\,dx^{\prime}=\Bigl{(}n-\frac{1}{4}\Bigr{)}\pi$
$\displaystyle\Longrightarrow\,\,\alpha^{1/2}\frac{\\{(3-s)(2-s)\\}^{1/(2-s)}}{2-s}\,B\Bigl{(}\frac{1}{2-s},\,\frac{3}{2}\Bigr{)}$
$\displaystyle\qquad\times\bigl{\\{}\tilde{\Psi}(0)+\mathcal{E}/\alpha\bigr{\\}}^{1/2+1/(2-s)}=\Bigl{(}n-\frac{1}{4}\Bigr{)}\pi,$
(127)

where the function $B$ is the beta function. Setting the integer $n$ to unity
(this corresponds to the ground state), the above equation is recast as

$\displaystyle\tilde{\Psi}(0)+\mathcal{E}/\alpha=\alpha^{-(2-s)/(4-s)}\Bigl{(}\frac{3\pi}{4}\Bigr{)}^{2(2-s)/(4-s)}$
$\displaystyle\times\Biggl{[}\frac{2-s}{\bigl{\\{}(3-s)(2-s)\bigr{\\}}^{1/(2-s)}}\frac{1}{B(1/(2-s),\,3/2)}\Biggr{]}^{2(2-s)/(4-s)}.$
(128)

Figure 14: Dependence of the core radius-halo mass relation on the inner slope
of halo density profile, $s$. The upper bound on the soliton core radius in
the limit of $c_{\rm vir}\to 0$, given at Eq. (130), is plotted as a function
of halo mass, fixing model parameters for the cosmology and FDM mass to the
default setup (see Sec. IV.1). Note that the results for $s=1$ are identical
to the lower and upper limits in the core radius and mass shown in Fig. 6.

Substituting Eq. (128) into Eq. (125), the comoving core radius valid in the
limit $\alpha\gg 1$ leads to

$\displaystyle r_{\rm c}$ $\displaystyle=p\,\frac{6^{1/2}\,r_{\rm
s}}{\alpha^{1/(4-s)}}$
$\displaystyle\times\Biggl{[}\frac{4}{3\pi}\,\frac{\bigl{\\{}(3-s)(2-s)\bigr{\\}}^{1/(2-s)}}{2-s}\,B\Bigl{(}\frac{1}{2-s},\,\frac{3}{2}\Bigr{)}\Biggr{]}^{(2-s)/(4-s)}$
(129)

Substituting further Eqs. (126) and (125) into the above, we arrive at the
final expression for the core radius in the limit of $c_{\rm vir}\to 0$:

$\displaystyle r_{\rm c}$
$\displaystyle=3.59\,\,\mbox{[kpc]}\,\,\frac{\mathcal{T}(s)}{\mathcal{T}(1)}a^{-1/(4-s)}\,\Bigl{(}\frac{p}{0.65}\Bigr{)}\Bigl{(}\frac{m_{\phi}}{10^{-22}\,\mbox{eV}}\Bigr{)}^{-2/(4-s)}\Bigl{(}\frac{M_{\rm
h}}{10^{9}\,M_{\odot}}\Bigr{)}^{-s/3/(4-s)}\Bigl{(}\frac{\Delta_{\rm
vir}}{200}\,\,\frac{\Omega_{\rm m,0}h^{2}}{0.147}\Bigr{)}^{-(1-s/3)/(4-s)},$
(130)

where the function $\mathcal{T}(s)$ is defined by

$\displaystyle\mathcal{T}(s)\equiv
49.064\,\bigl{\\{}721.305\,(3-s)\bigr{\\}}^{-1/(4-s)}$
$\displaystyle\times\Biggl{[}\frac{4}{3\pi}\,\frac{\bigl{\\{}(3-s)(2-s)\bigr{\\}}^{1/(2-s)}}{2-s}\,B\Bigl{(}\frac{1}{2-s},\,\frac{3}{2}\Bigr{)}\Biggr{]}^{(2-s)/(4-s)}.$
(131)

Notice that the dependence on the concentration parameter $c_{\rm vir}$
disappears at Eq. (130). This means that for a fixed cosmological model and a
given halo mass, Eq. (130) gives a bound on the core radius. Note that setting
the slope $s$ to unity, the above expression reproduces the upper bound given
at Eq. (74).

Figure 13 plots the functions $\mathcal{T}$ normalized by the one evaluated at
$s=1$. This shows that the function $\mathcal{T}$ is a monotonically
decreasing functions of $s$, and approaches zero at $s\to 2$, indicating that
for halos with the inner slope of $s=2$, the soliton core disappears. Since
the potential at the halo center generally diverges for $s\geq 2$, this
implies that for $s\geq 2$, the fuzzy dark matter soliton, if exists, does not
form a flat core, but rather can have a cusp at the center.

Figure 14 shows the core-halo relations in the $c_{\rm vir}\to 0$ limit,
focusing specifically on the soliton core radius. Setting the model parameters
to those adopted in Figs. 6-9, we vary the inner slope of halo profiles from
$s=0$ to $1.5$. Here, Steepening the inner slope substantially reduces
allowable core sizes, and the predicted curve for $s=1.5$ is inconsistent with
simulations.

## References

  * Peebles [1982] P. J. E. Peebles, ApJ 263, L1 (1982).
  * Bond  _et al._ [1982] J. R. Bond, A. S. Szalay, and M. S. Turner, Phys. Rev. Lett. 48, 1636 (1982).
  * Blumenthal  _et al._ [1984] G. R. Blumenthal, S. M. Faber, J. R. Primack, and M. J. Rees, Nature 311, 517 (1984).
  * Bullock and Boylan-Kolchin [2017] J. S. Bullock and M. Boylan-Kolchin, ARA&A 55, 343 (2017), arXiv:1707.04256 [astro-ph.CO] .
  * Tulin and Yu [2018] S. Tulin and H.-B. Yu, Phys. Rep. 730, 1 (2018), arXiv:1705.02358 [hep-ph] .
  * Hu  _et al._ [2000] W. Hu, R. Barkana, and A. Gruzinov, Phys. Rev. Lett. 85, 1158 (2000), arXiv:astro-ph/0003365 [astro-ph] .
  * Marsh [2016a] D. J. E. Marsh, Phys. Rep. 643, 1 (2016a), arXiv:1510.07633 [astro-ph.CO] .
  * Niemeyer [2020] J. C. Niemeyer, Progress in Particle and Nuclear Physics 113, 103787 (2020), arXiv:1912.07064 [astro-ph.CO] .
  * Hui [2021] L. Hui, ARA&A 59, 10.1146/annurev-astro-120920-010024 (2021), arXiv:2101.11735 [astro-ph.CO] .
  * Ferreira [2021] E. G. M. Ferreira, A&A Rev. 29, 7 (2021), arXiv:2005.03254 [astro-ph.CO] .
  * Hui  _et al._ [2017] L. Hui, J. P. Ostriker, S. Tremaine, and E. Witten, Phys. Rev. D 95, 043541 (2017), arXiv:1610.08297 [astro-ph.CO] .
  * Hložek  _et al._ [2018] R. Hložek, D. J. E. Marsh, and D. Grin, MNRAS 476, 3063 (2018), arXiv:1708.05681 [astro-ph.CO] .
  * Laguë  _et al._ [2022] A. Laguë, J. R. Bond, R. Hložek, K. K. Rogers, D. J. E. Marsh, and D. Grin, J. Cosmology Astropart. Phys 2022, 049 (2022), arXiv:2104.07802 [astro-ph.CO] .
  * Armengaud  _et al._ [2017] E. Armengaud, N. Palanque-Delabrouille, C. Yèche, D. J. E. Marsh, and J. Baur, MNRAS 471, 4606 (2017), arXiv:1703.09126 [astro-ph.CO] .
  * Kobayashi  _et al._ [2017] T. Kobayashi, R. Murgia, A. De Simone, V. Iršič, and M. Viel, Phys. Rev. D 96, 123514 (2017), arXiv:1708.00015 [astro-ph.CO] .
  * Rogers and Peiris [2021] K. K. Rogers and H. V. Peiris, Phys. Rev. Lett. 126, 071302 (2021), arXiv:2007.12705 [astro-ph.CO] .
  * Khlopov  _et al._ [1985] M. I. Khlopov, B. A. Malomed, and I. B. Zeldovich, MNRAS 215, 575 (1985).
  * Nambu and Sasaki [1990] Y. Nambu and M. Sasaki, Phys. Rev. D 42, 3918 (1990).
  * Schive  _et al._ [2014a] H.-Y. Schive, T. Chiueh, and T. Broadhurst, Nature Physics 10, 496 (2014a), arXiv:1406.6586 [astro-ph.GA] .
  * Schive  _et al._ [2014b] H.-Y. Schive, M.-H. Liao, T.-P. Woo, S.-K. Wong, T. Chiueh, T. Broadhurst, and W. Y. P. Hwang, Phys. Rev. Lett. 113, 261302 (2014b), arXiv:1407.7762 [astro-ph.GA] .
  * Schwabe  _et al._ [2016] B. Schwabe, J. C. Niemeyer, and J. F. Engels, Phys. Rev. D 94, 043513 (2016), arXiv:1606.05151 [astro-ph.CO] .
  * Du  _et al._ [2017] X. Du, C. Behrens, J. C. Niemeyer, and B. Schwabe, Phys. Rev. D 95, 043519 (2017), arXiv:1609.09414 [astro-ph.GA] .
  * Mocz  _et al._ [2017] P. Mocz, M. Vogelsberger, V. H. Robles, J. Zavala, M. Boylan-Kolchin, A. Fialkov, and L. Hernquist, MNRAS 471, 4559 (2017), arXiv:1705.05845 [astro-ph.CO] .
  * Mina  _et al._ [2022] M. Mina, D. F. Mota, and H. A. Winther, Astron. Astrophys. 662, A29 (2022), arXiv:2007.04119 [astro-ph.CO] .
  * Nori and Baldi [2021] M. Nori and M. Baldi, MNRAS 501, 1539 (2021), arXiv:2007.01316 [astro-ph.CO] .
  * May and Springel [2021] S. May and V. Springel, MNRAS 506, 2603 (2021), arXiv:2101.01828 [astro-ph.CO] .
  * Chan  _et al._ [2022] H. Y. J. Chan, E. G. M. Ferreira, S. May, K. Hayashi, and M. Chiba, MNRAS 511, 943 (2022), arXiv:2110.11882 [astro-ph.CO] .
  * Veltmaat  _et al._ [2018] J. Veltmaat, J. C. Niemeyer, and B. Schwabe, Phys. Rev. D 98, 043509 (2018), arXiv:1804.09647 [astro-ph.CO] .
  * Schive  _et al._ [2020] H.-Y. Schive, T. Chiueh, and T. Broadhurst, Phys. Rev. Lett. 124, 201301 (2020), arXiv:1912.09483 [astro-ph.GA] .
  * Li  _et al._ [2021] X. Li, L. Hui, and T. D. Yavetz, Phys. Rev. D 103, 023508 (2021), arXiv:2011.11416 [astro-ph.CO] .
  * Zagorac  _et al._ [2022] J. L. Zagorac, I. Sands, N. Padmanabhan, and R. Easther, Phys. Rev. D 105, 103506 (2022), arXiv:2109.01920 [astro-ph.CO] .
  * Chavanis [2019] P.-H. Chavanis, Phys. Rev. D 100, 123506 (2019), arXiv:1905.08137 [astro-ph.CO] .
  * Desjacques and Nusser [2019] V. Desjacques and A. Nusser, MNRAS 488, 4497 (2019), arXiv:1905.03450 [astro-ph.CO] .
  * Safarzadeh and Spergel [2020] M. Safarzadeh and D. N. Spergel, ApJ 893, 21 (2020), arXiv:1906.11848 [astro-ph.CO] .
  * Hayashi and Obata [2020] K. Hayashi and I. Obata, MNRAS 491, 615 (2020), arXiv:1902.03054 [astro-ph.CO] .
  * Broadhurst  _et al._ [2020] T. Broadhurst, I. De Martino, H. N. Luu, G. F. Smoot, and S. H. H. Tye, Phys. Rev. D 101, 083012 (2020), arXiv:1902.10488 [astro-ph.CO] .
  * Burkert [2020] A. Burkert, ApJ 904, 161 (2020), arXiv:2006.11111 [astro-ph.GA] .
  * Hayashi  _et al._ [2021] K. Hayashi, E. G. M. Ferreira, and H. Y. J. Chan, ApJ 912, L3 (2021), arXiv:2102.05300 [astro-ph.CO] .
  * Langer [1932] R. E. Langer, Trans. Amer. Math. Soc. 34, 447 (1932).
  * Langer [1959] R. E. Langer, Trans. Amer. Math. Soc. 90, 113 (1959).
  * Olver [1975] F. W. J. Olver, Philosophical Transactions of the Royal Society of London Series A 278, 137 (1975).
  * Widrow and Kaiser [1993] L. M. Widrow and N. Kaiser, ApJ 416, L71 (1993).
  * Uhlemann  _et al._ [2014] C. Uhlemann, M. Kopp, and T. Haugg, Phys. Rev. D 90, 023517 (2014), arXiv:1403.5567 [astro-ph.CO] .
  * Uhlemann  _et al._ [2019] C. Uhlemann, C. Rampf, M. Gosenca, and O. Hahn, Phys. Rev. D 99, 083524 (2019), arXiv:1812.05633 [astro-ph.CO] .
  * Garny  _et al._ [2020] M. Garny, T. Konstandin, and H. Rubira, J. Cosmology Astropart. Phys 2020, 003 (2020), arXiv:1911.04505 [astro-ph.CO] .
  * Navarro  _et al._ [1996] J. F. Navarro, C. S. Frenk, and S. D. M. White, ApJ 462, 563 (1996), arXiv:astro-ph/9508025 [astro-ph] .
  * Navarro  _et al._ [1997] J. F. Navarro, C. S. Frenk, and S. D. M. White, ApJ 490, 493 (1997), arXiv:astro-ph/9611107 [astro-ph] .
  * Saga  _et al._ [2020] S. Saga, A. Taruya, M.-A. Breton, and Y. Rasera, MNRAS 498, 981 (2020), arXiv:2004.03772 [astro-ph.CO] .
  * Nayfeh [1993] A. H. Nayfeh,  _Introduction to perturbation techniques_ (Wiley-VCH, 1993).
  * Abramowitz and Stegun [1965] M. Abramowitz and I. A. Stegun,  _Handbook of mathematical functions with formulas, graphs, and mathematical tables_ (Dover, 1965).
  * Luna Zagorac  _et al._ [2021] J. Luna Zagorac, I. Sands, N. Padmanabhan, and R. Easther, arXiv e-prints , arXiv:2109.01920 (2021), arXiv:2109.01920 [astro-ph.CO] .
  * Schive  _et al._ [2014c] H.-Y. Schive, T. Chiueh, and T. Broadhurst, Nature Physics 10, 496 (2014c), arXiv:1406.6586 [astro-ph.GA] .
  * Schive  _et al._ [2014d] H.-Y. Schive, M.-H. Liao, T.-P. Woo, S.-K. Wong, T. Chiueh, T. Broadhurst, and W. Y. P. Hwang, Phys. Rev. Lett. 113, 261302 (2014d), arXiv:1407.7762 [astro-ph.GA] .
  * Dentler  _et al._ [2022] M. Dentler, D. J. E. Marsh, R. Hložek, A. Laguë, K. K. Rogers, and D. Grin, MNRAS 515, 5646 (2022), arXiv:2111.01199 [astro-ph.CO] .
  * Bullock  _et al._ [2001] J. S. Bullock, T. S. Kolatt, Y. Sigad, R. S. Somerville, A. V. Kravtsov, A. A. Klypin, J. R. Primack, and A. Dekel, MNRAS 321, 559 (2001), arXiv:astro-ph/9908159 [astro-ph] .
  * Schneider  _et al._ [2012] A. Schneider, R. E. Smith, A. V. Macciò, and B. Moore, MNRAS 424, 684 (2012), arXiv:1112.0330 [astro-ph.CO] .
  * Ludlow  _et al._ [2016] A. D. Ludlow, S. Bose, R. E. Angulo, L. Wang, W. A. Hellwing, J. F. Navarro, S. Cole, and C. S. Frenk, MNRAS 460, 1214 (2016), arXiv:1601.02624 [astro-ph.CO] .
  * Diemer and Joyce [2019] B. Diemer and M. Joyce, ApJ 871, 168 (2019), arXiv:1809.07326 [astro-ph.CO] .
  * Wang  _et al._ [2020] J. Wang, S. Bose, C. S. Frenk, L. Gao, A. Jenkins, V. Springel, and S. D. M. White, Nature 585, 39 (2020), arXiv:1911.09720 [astro-ph.CO] .
  * Ishiyama  _et al._ [2021] T. Ishiyama, F. Prada, A. A. Klypin, M. Sinha, R. B. Metcalf, E. Jullo, B. Altieri, S. A. Cora, D. Croton, S. de la Torre, D. E. Millán-Calero, T. Oogi, J. Ruedas, and C. A. Vega-Martínez, MNRAS 506, 4210 (2021), arXiv:2007.14720 [astro-ph.CO] .
  * Jing [2000] Y. P. Jing, ApJ 535, 30 (2000), arXiv:astro-ph/9901340 [astro-ph] .
  * Neto  _et al._ [2007] A. F. Neto, L. Gao, P. Bett, S. Cole, J. F. Navarro, C. S. Frenk, S. D. M. White, V. Springel, and A. Jenkins, MNRAS 381, 1450 (2007), arXiv:0706.2919 [astro-ph] .
  * Reed  _et al._ [2011] D. S. Reed, S. M. Koushiappas, and L. Gao, MNRAS 415, 3177 (2011), arXiv:1008.1579 [astro-ph.CO] .
  * Bhattacharya  _et al._ [2013] S. Bhattacharya, S. Habib, K. Heitmann, and A. Vikhlinin, ApJ 766, 32 (2013), arXiv:1112.5479 [astro-ph.CO] .
  * Diemer and Kravtsov [2015] B. Diemer and A. V. Kravtsov, ApJ 799, 108 (2015), arXiv:1407.4730 [astro-ph.CO] .
  * Chiang  _et al._ [2021] B. T. Chiang, H.-Y. Schive, and T. Chiueh, Phys. Rev. D 103, 103019 (2021), arXiv:2104.13359 [astro-ph.CO] .
  * Lin  _et al._ [2018] S.-C. Lin, H.-Y. Schive, S.-K. Wong, and T. Chiueh, Phys. Rev. D 97, 103523 (2018), arXiv:1801.02320 [astro-ph.CO] .
  * Yavetz  _et al._ [2022] T. D. Yavetz, X. Li, and L. Hui, Phys. Rev. D 105, 023512 (2022), arXiv:2109.06125 [astro-ph.CO] .
  * Dalal  _et al._ [2021] N. Dalal, J. Bovy, L. Hui, and X. Li, J. Cosmology Astropart. Phys 2021, 076 (2021), arXiv:2011.13141 [astro-ph.CO] .
  * Dandoy  _et al._ [2022] V. Dandoy, T. Schwetz, and E. Todarello, J. Cosmology Astropart. Phys 2022, 081 (2022), arXiv:2206.04619 [astro-ph.CO] .
  * Kawai  _et al._ [2022] H. Kawai, M. Oguri, A. Amruth, T. Broadhurst, and J. Lim, ApJ 925, 61 (2022), arXiv:2109.04704 [astro-ph.CO] .
  * Dalal and Kravtsov [2022] N. Dalal and A. Kravtsov, arXiv e-prints , arXiv:2203.05750 (2022), arXiv:2203.05750 [astro-ph.CO] .
  * Gunn and Gott [1972] J. E. Gunn and I. Gott, J. Richard, ApJ 176, 1 (1972).
  * Gunn [1977] J. E. Gunn, ApJ 218, 592 (1977).
  * Peebles [1980] P. J. E. Peebles,  _The large-scale structure of the universe_ (Princeton University Press, 1980).
  * Bryan and Norman [1998] G. L. Bryan and M. L. Norman, ApJ 495, 80 (1998), arXiv:astro-ph/9710107 [astro-ph] .
  * Nakamura and Suto [1997] T. T. Nakamura and Y. Suto, Progress of Theoretical Physics 97, 49 (1997), arXiv:astro-ph/9612074 [astro-ph] .
  * Mead  _et al._ [2015] A. J. Mead, J. A. Peacock, C. Heymans, S. Joudaki, and A. F. Heavens, MNRAS 454, 1958 (2015), arXiv:1505.07833 [astro-ph.CO] .
  * Marsh and Silk [2014] D. J. E. Marsh and J. Silk, MNRAS 437, 2652 (2014), arXiv:1307.1705 [astro-ph.CO] .
  * Marsh [2016b] D. J. E. Marsh, arXiv e-prints , arXiv:1605.05973 (2016b), arXiv:1605.05973 [astro-ph.CO] .