Towards an Analytic Construction of the Wavefunction of Boson Stars

TL;DR

This paper introduces a semi-analytic, highly accurate method for constructing the wavefunction profile of non-interacting boson stars using the Gross-Pitaevskii-Poisson equations in the Newtonian limit.

Contribution

It provides a novel semi-analytic solution to the equations governing boson stars, improving stability and accuracy over previous approaches.

Findings

Solution is stable and numerically precise with better than 10^{-6} accuracy.

Applicable to non-interacting scalar fields forming boson stars.

Enhances understanding of boson star profiles in the Newtonian regime.

Abstract

Light scalar fields can form gravitationally bound compact objects called boson stars. The profile of boson stars in the Newtonian limit is described by the Gross-Pitaevskii-Poisson equations. We present a semi-analytic solution to these equations and construct the profile of boson stars formed by a non-interacting scalar field. Our solution is stable with respect to numerical errors and has accuracy better than $10^{-6}$ over the entire range.

Equations

$$\mathcal{L}=g^{\mu \nu }({\partial }_{\mu }{\varphi }^{\ast })({\partial }_{\nu }\varphi )-m^2({\varphi }^{\ast }\varphi )-\frac{\lambda }{2}({\varphi }^{\ast }\varphi {)}^2$$

$$\varphi (\vec{r},t)=\sum _{nlm}{R}_{nl}(r)Y_{lm}(\theta ,\varphi )e^{-i{E}_{nlm}t}.$$

$$(1+2\Phi {)}^{-1}{\partial }_t^2\psi -{\nabla }^2\psi +m^2\psi +\frac{N\lambda }{2m}{\psi }^3=0.$$

$$e\psi =-\frac{1}{2m}{\nabla }^2\psi +m\Phi \psi +\frac{N\lambda }{4m^2}{\psi }^3$$

$$\rho =|{\partial }_t\varphi {|}^2+|\nabla \varphi {|}^2+m^2|\varphi {|}^2+\lambda |\varphi {|}^4\approx Nm{\psi }^2$$

$${\nabla }^2\Phi =4\pi GNm{\psi }^2$$

$$x=2rm,V=\frac{e}{2m}-\frac{\Phi }{2},\text{and}S=\sqrt{\frac{\pi GN}{2m}}\psi$$

$${\nabla }^{2}V=-S^2\text{and}{\nabla }^{2}S=-VS.$$

$${\int }_0^{\infty }{x}^2{S}^{2}dx=GMm,$$

$$x\to \frac{x}{f},S\to f^{2}S,V\to f^{2}V,M\to fM.$$

$$z=kx,S=k^{2}s,V=k^{2}v,GMm=2k\beta .$$

$${\nabla }^{2}s=-sv\text{and}{\nabla }^{2}v=-s^2\text{with}v(\infty )=-1.$$

$$s=\sum _{n=0}^{\infty }{s}_n{z}^n\text{and}v=\sum _{n=0}^{\infty }{v}_n{z}^n.$$

$$\nabla s=s^{\prime \prime }+\frac{2}{z}{s}^{\prime }=\sum _{n=0}^{\infty }(n+2)(n+3)s_{n+2}{z}^n.$$

$$\displaystyle sv=\Bigg{[}\sum_{i=0}^{\infty}s_{i}z^{i}\Bigg{]}\Bigg{[}\sum_{j=0}^{\infty}v_{j}z^{j}\Bigg{]}=\sum_{n=0}^{\infty}\sum_{m=0}^{n}s_{m}v_{n-m}z^{n}.$$

$$(n+2)(n+3)s_{n+2}$$

$$(n+2)(n+3)v_{n+2}$$

$$s=\sum _{n,m=0,0}^{\infty ,\infty }{s}_m^n{\left(\frac{e^{-z}}{z^{\sigma }}\right)}^n{z}^{-m},v=\sum _{n,m=0,0}^{\infty ,\infty }{v}_m^n{\left(\frac{e^{-z}}{z^{\sigma }}\right)}^n{z}^{-m}$$

$$\displaystyle\nabla^{2}s=\sum_{n,m=0,0}^{\infty,\infty}\Big{(}s^{n}_{m}n^{2}+2s^{n}_{m-1}n(n\sigma+m-2)$$

$$\displaystyle+s^{n}_{m-2}(n\sigma+m-3)(n\sigma+m-2)\Big{)}\left(\frac{e^{-z}}{z^{\sigma}}\right)^{n}z^{-m}$$

$$sv$$

$$=\sum _{n,m=0,0}^{\infty ,\infty }\left(\sum _{p,q=0,0}^{n,m}{s}_q^p{v}_{m-q}^{n-p}\right){\left(\frac{e^{-z}}{z^{\sigma }}\right)}^n{z}^{-m}$$

$$\sum _{p,q=0,0}^{n,m}$$

$$s_1^1+2(\sigma -1)s_0^1=-s_1^1{v}_0^0-s_0^1{v}_1^0$$

$$s_{M+1}^1$$

$$+(\sigma +M-1)(\sigma +M-2)s_{M-1}^1=0$$

$$s_M^1=-\frac{(\sigma +M-1)(\sigma +M-2)}{2M}{s}_{M-1}^1.$$

$$s_M^1=s_0^1\frac{(\sigma {)}_M(\sigma -1{)}_M}{M!}\cdot (-2{)}^{-M}$$

$$s^{(1)}=\alpha \frac{e^{-z}}{z^{1-\beta }}\sum _{m=0}^{\infty }\left(\frac{(1-\beta {)}_m(-\beta {)}_m}{m!}(-2z{)}^{-m}\right)$$

$$s^{(1)}=\frac{\alpha }{2^{\beta }z}{W}_{\beta ,-\frac{1}{2}}(2z)=\frac{\alpha }{2^{\beta }z}{e}^{-z}U(-\beta ,0,2z)$$

$$s^{(1)}=s_{(M)}^{(1)}+R_M^1\text{with}{s}_{(M)}^{(1)}=\sum _{m=0}^M{s}_m^1\frac{e^{-z}}{z^{\beta -1}}{z}^{-m}$$

$$s=s_{(M)}^{(N)}+R_M^N\text{with}{s}_{(M)}^{(N)}=\sum _{n,m=0,0}^{N,M}{s}_m^n{\left(\frac{e^{-z}}{z^{\beta -1}}\right)}^n{z}^{-m}$$

Full text

††preprint: UCI-HEP-TR-2017-04

I Introduction

Axions are an attractive solution 1977PhRvL..38.1440P ; 1977PhRvD..16.1791P ; 1978PhRvL..40..223W ; 1978PhRvL..40..279W ; 1981PhLB..104..199D to the strong CP problem of QCD 1976PhLB…63..334C ; 1976PhRvL..37..172J , and also provide an attractive and natural dark matter candidate (for reviews see Kim:2008hd ; Cheng:1987gp ). This has motivated multiple search strategies for axions and axion-like particles 1983PhRvL..51.1415S ; 1985PhRvD..32.2988S ; 2010PhRvL.104d1301A ; 2011PhRvD..84l1302H ; Rybka:2014cya ; Sikivie:2013laa ; 2012PhRvD..85c5018B ; Horns:2012jf ; Budker:2013hfa ; Graham:2013gfa ; 1999NIMPA.425..480Z ; Stadnik:2013raa . It has also been argued that ultra-light axions Hui:2016ltb can solve problems encountered by the usual CDM models 2011AJ….141..193O ; 1993MNRAS.264..201K ; 1999ApJ…522…82K ; 2000PhRvL..84.4525K .

For all these purposes, it is important to understand the spatial structure of axion-like particles if they constitute dark matter. Most importantly, it is crucial to know whether or not these particles clump into compact objects (see Liebling:2012fv for a review). The Jeans instability 1985MNRAS.215..575K ; Bianchi:1990mha ; 1994PhRvL..72.2516S ; Mielke:2000mh ; 1986PhRvL..57.2485C indicates that a uniform density of axions is unstable, indicating the formation of large compact objects (sometimes called boson stars 1968PhRv..172.1331K ; Ruffini:1969qy ). The boson stars are prevented from completely collapsing by a bosonic analogue of the Fermi pressure in white dwarfs 1926MNRAS..87..114F .

Once the boson stars form, further cosmological evolution can occur by scattering of these boson stars off other boson stars, as well as baryonic matter Eby:2017xaw . Such scatterings may either enhance or decrease the number of these objects. To analyze these scatterings, we must have a detailed understanding of the bound states, including their energies and profiles.

At least for fairly dilute111There is a different set of solutions with large density: dense boson stars Braaten:2016kzc ; Braaten:2015eeu ; Eby:2016cnq ; Bai:2016wpg . systems, the compact objects are bound states of a nonlinear Gross-Pitaevskii-Poisson equations, which we re-derive below. Numerical solutions to these equation have found the bound state energies and mass-radius relation, both in the cases with no self-interactions, as well as including self interactions, either attractive or repulsive Ruffini:1969qy ; Barranco:2010ib ; 2007JCAP…06..025B ; 1989PhRvA..39.4207M ; 1996PhRvD..53.2236L ; 2000NewA….5..103G ; 2003PhRvD..68b3511A ; Chavanis:2011zm ; Eby:2014fya ; Eby:2015hsq ; Yang:2015paa ; Kan:2017uhj . However, the profiles are only available numerically; they are computationally expensive to find and are difficult to extend to perturbations of the boson stars. Furthermore, they tend to have numerical instabilities in the tails of the profile. The author of Chavanis:2011zi follows a different approach: using a simple Gaussian ansatz for the density profile they were able to obtain the mass-radius relation and ground-state energy within a 10% deviation of the numerical solution. The Gross-Pitaevskii-Poisson equations have also been studied in the context of quantum state reduction by Penrose1998 ; Moroz:1998dh ; Tod1992 .

In this paper, we introduce a combination of analytical and numerical methods to find approximate solutions to the Gross-Pitaevskii-Poisson equations. We illustrate this for the case of no self-interactions among the bosons (the interacting case will be treated in an upcoming paper Kling:2017hjm ). We also perform a detailed analysis of the accuracy of our results. We show that our method is much less computationally expensive than previous approaches, nevertheless we demonstrate that we find excellent agreement with a full numerical solution over the entire parameter range. Furthermore, our method is numerically stable and does not diverge far away from the star, which can occur for a purely numerical approach.

II The Gross-Pitaevskii-Poisson Equations

II.1 Derivation

Let us consider a complex scalar field $\phi(\vec{r},t)$ described by the Lagrangian

[TABLE]

We can then expand the field $\phi$ in spherical harmonics222Note that this expansion allows us to choose all $R_{nl}(r)$ to be real.

[TABLE]

We assume that only the ground state $(n,l,m)=(1,0,0)$ is populated. In this case the field takes the simple form $\phi(\vec{r},t)=(2E/N)^{-\frac{1}{2}}\psi(r)e^{-iEt}$, where we denote the ground-state energy as $E$. The real function $\psi(r)$ describes the radial profile of the star and is sometimes called the wavefunction. We have chosen a normalization $\int\psi^{2}dV=1$, which allows us to identify $\psi^{2}$ with the probability density. Let us further make two simplifying assumptions

Newtonian Gravity

The field is only weakly coupled to gravity such that we can use a Newtonian approximation. This allows us to introduce the Newtonian potential $\Phi$ in the metric $g_{\mu\nu}=\text{diag}(1+2\Phi,-1,-1,-1)$. The Newtonian potential is related to the energy density via the Poisson equation $\nabla^{2}\Phi=4\pi G\rho$.

Non-Relativistic

The ground state is non-relativistic. In this case we can write $E=m+e$ with binding energy $e\ll m$. This implies $e\psi,\Phi\psi,\nabla\psi\ll m\psi$.

The equation of motion is the Klein-Gordon equation, $\Box\phi+m^{2}\phi+\lambda|\phi|^{2}\phi=0$, which we can write in terms of the wavefunction as

[TABLE]

In the non-relativistic approximation we can write

[TABLE]

and recover the Schrödinger-type equation. For our non-relativistic approximation to be consistent, the last term should be sufficiently small i.e. $\frac{N\lambda}{4m^{2}}\psi^{2}\ll m$.

The energy density of the complex scalar field is

[TABLE]

where we used the non-relativistic approximation in the last step. Newtons equation therefore takes the simple form

[TABLE]

II.2 Scaling Symmetry of the Gross-Pitaevskii-Poisson System

As derived above, the ground state of the self-gravitating boson star in the non-relativistic limit can be described by a wavefunction $\psi(r)$ and a gravitational potential $\Phi(r)$ satisfying the Gross-Pitaevskii-Poisson equations given in (4) and (6). For the remainder of this paper, we will focus on the non-interacting case $\lambda=0$, and postpone the discussion of finite self-interactions $\lambda\neq 0$ to a separate study Kling:2017hjm . For simplification, let us introduce the following dimensionless variables,

[TABLE]

The equations (4) and (6) become

[TABLE]

The wavefunction normalization condition $\int\psi^{2}dV=1$ can be rewritten as

[TABLE]

where we introduced the star mass $M=Nm$. These equations describe the hydrostatic equilibrium of the boson star, in which the gravitational attraction caused by the potential $V$ is balanced by a repulsive quantum pressure. This quantum pressure arises from Heisenberg’s uncertainty principle and prevents the system from gravitational collapse.

Let us note that Eq. 8 and 9 are invariant under the scaling

[TABLE]

where $f$ is a scaling factor. This implies that solutions corresponding to different masses $M$ can be related to a unique solution by rescaling. To obtain the unique solution, we have to fix the scale by choosing a reference scale $k$. Although there are many different ways to fix $k$, a particularly useful choice for our discussion333For the numerical integration in appendix A, we will choose a different reference scale $k$ which sets $V(0)=1$. is to set $-k^{2}=V(\infty)=\frac{e}{2m}$ which transforms as $k\to fk$. We can then introduce the scale invariant coordinate $z$, wavefunction $s$, potential $v$ and mass $\beta$ via

[TABLE]

Note that $S$ and $V$ are functions of the scaling dependent coordinate $x$ while $s$ and $v$ are functions of the scale independent coordinate $z$. Using the scale independent variables, we can write the Gross-Pitaevskii-Poisson equations as

[TABLE]

To obtain the solution corresponding to a boson star with mass $M$, we then have to perform the rescaling given in Eq. 11 with $k=\frac{GMm}{2\beta}$. In the following section, we will obtain an approximate analytical form for $s,v$ and the mass parameter $\beta$.

III Series Expansions

It has been shown Tod1992 ; Choquard:2007aa that the Gross-Pitaevskii-Poisson system given by (12) has a unique square normalizable solution for $s,v$ with $s>0$ for all values of $z$. The authors of Moroz:1998dh also provide a numerical solution. However, the authors have also shown that this numerical solution will diverge at some finite value of $z$ and can therefore not be used to describe the profile at large radius. The author of Chavanis:2011zi follows a different approach and approximates the density profile by a Gaussian, which is able to approximately reproduce thermodynamical properties of the boson star but otherwise fails to describe the profile, in particular at large radius. We attempt to solve this problem by providing an analytical expression for $s$ and $v$ which describes the profile at all radii with high precision.

III.1 Expansion at Small Radius: $z=0$

Since both $s$ and $v$ are well behaved around $z=0$, we can expand them in a series expansion

[TABLE]

We can then write the Laplacian on the left hand side of Eq. 12 as

[TABLE]

and the right hand side as a Cauchy product

[TABLE]

By matching the coefficients in Eq. 12 we obtain the recursion relations

[TABLE]

Note that requiring the profile to be smooth at the origin implies $s_{1}=v_{1}=0$ and therefore also all odd coefficients $s_{2n+1},v_{2n+1}$ vanish. The profile at small radius $z$ can therefore be fully parametrized in terms of the wavefunction and potential at the origin: $s_{0}$ and $v_{0}$. We can obtain $s_{0}$ and $v_{0}$ from a fit to the numerical solution, as discussed in appendix A.

For practical purposes, we will truncate the infinite series of Eq. 13 at some $N$. This is shown in Fig. 1. The upper panel shows both the numerical solution as well as the truncated expansion $s_{(N)}$ for different values of $N$. The lower panel shows the deviations of the truncated expansion from the numerical solution. We can see that already a small number of terms in the series expansion provides a sub-percent level accuracy for the inner part of the profile. A better accuracy can be obtained by including more terms in the expansion. However, the accuracy of the series expansion is limited by the accuracy of $s_{0}$ and $v_{0}$ which in Fig. 1 is about $10^{-12}$. Note that the series expansion of Eq. 13 diverges at $z>4$ and a different parametrization has to be chosen.

III.2 Expansion at Large Radius: $z=\infty$

At large radius $z\to\infty$, we expect the wavefunction to decrease at least exponentially, $\psi\sim e^{-kr}$, and the potential to approach the Newtonian potential of a central point mass $\Phi\approx\frac{GM}{r}$. A series expansion at large radius must be able to reproduce these limits. Let us choose the following general ansatz for the form of the solution

[TABLE]

Here we assume the existence of a $\sigma\in\mathbb{R}$, whose meaning will become clear later. The Laplacian of $s$ can be written as

[TABLE]

Note that we introduced the short-hand notation $s_{n,-1}=s_{n,-2}=0$. The right hand side of Eq. 12 can be rewritten as Cauchy product:

[TABLE]

By matching the coefficients, Eq. 12 we obtain the recursion relations

[TABLE]

Let us note the following properties of $s$ and $v$: i) Normalizability requires $s^{0}_{0}=0$. Eq. 20 then implies that all coefficients $s^{0}_{m}$ vanish as well. This means that the wavefunction decays at least exponentially. ii) Eq. 21 then implies that all $v^{0}_{m}=0$ for $m>1$. This means that at large radius, the potential is described by the Newtonian potential $v^{(0)}=v^{0}_{0}+\frac{v^{0}_{1}}{z}$. All other terms in the expansion of $v$ are at least exponentially suppressed. iii) Eq. 20 and 21 further imply that the potential contains only non-vanishing components $v^{n}_{m}$ for even $n$ while the wavefunction only has non-vanishing components $s^{n}_{m}$ for odd $n$.

IV The Solution for the Wavefunction

IV.1 The Wavefunction at Leading Order

Using Eq. 20 and 21, we are able to recursively calculate all coefficients in the expansion of $s$ and $v$. Let us first analyze the $n=1$ expansion of $s$ which provides us both a leading order approximation and a deeper understanding about the form of the solution. We have seen before that at leading order the potential is given by $v^{(0)}=v^{0}_{0}+\frac{v^{0}_{1}}{z}$. Setting $m=0$, Eq. 20 reads $s^{1}_{0}=-s^{1}_{0}v^{0}_{0}$ which implies $v^{0}_{0}=-1$ as expected from our scale choice which fixes $v(\infty)=-1$. Setting $m=1$, Eq. 20 reads

[TABLE]

which implies $v^{0}_{1}=2(1-\sigma)$. This is a remarkable result: the exponent $\sigma$ in the series expansion is related to the the total mass of the system. In the notation of Eq. 11, we can write $\sigma=1-\frac{1}{2}v^{0}_{1}=1-\beta$. Let us now calculate the remaining coefficients by setting $m=M+1$ with $M\geq 1$. Then Eq. 20 can be written as

[TABLE]

We can therefore recursively compute the coefficients $s^{1}_{M}$ by

[TABLE]

Using the rising factorials $(x)_{n}=x(x+1)(x+2)\cdots(x+n-1)$, $(x)_{0}=1$ we can write the coefficients explicitly as

[TABLE]

The leading order wavefunction $s^{(1)}$ can therefore be written as

[TABLE]

Here we have introduced the normalization parameter $\alpha=s^{1}_{0}$. The far-field solution is described by only two free parameters: the wavefunction normalization $\alpha$ and the total mass parameter $\beta$. At very large radius, the wavefunction approaches $s=\alpha e^{-z}z^{\beta-1}$.

Let us note that we can rewrite the far-field solution as

[TABLE]

where $W_{\beta,-\frac{1}{2}}(2z)$ is the Whittaker function which can also be expressed in terms of the confluent hypergeometric function $U$. This result is not surprising: when considering the leading order potential $v^{(0)}=-1+\frac{2\beta}{z}$, Eq. 12 turns into the Whittaker equations. For a more detailed discussion on different representations of the leading order wavefunction, see appendix B.

It is known that the series expansion of the confluent hypergeometric function $U(a,b,z)$, and therefore also the expansion of $s^{(1)}$, are not Cauchy convergent. However, the infinite sum has a finite value and behaves convergently for a finite number of terms. This is illustrated in the upper panel of Fig. 2 where we show the coefficients $s^{1}_{m}$. We can see that they converge for $m<6$ while they start to diverge again for $m>6$. Let us therefore split the function $s^{(1)}$ into a finite series $s^{(1)}_{(M)}$ truncated after $M$ terms and the corresponding remainder $R_{M}$:

[TABLE]

In the right panel of Fig. 2 we show the relative size of the remainder $R^{1}_{M}$ with respect to $s^{(1)}$ for different values of $z$. We can see that finite series first converges quickly, even for values of $z$ close the convergence radius $z>2$. At some $M$, which depends on the value of $z$, the remainder reaches a minimum and diverges for large $M$. Note that the coefficients $s^{1}_{m}$ and therefore also the remainders $R_{M}$ are oscillating which allows the infinite series $s^{(1)}$ to be finite. In this work, we will avoid the problem of divergence by truncating the series expansion of $s^{(1)}$ at $m=M$ and ignoring the remainder. This precision of this approximation should be sufficient for most applications444As shown in Daalhuis1992 , it is possible to perform a hyperasymptotic expansions in which we truncate the series at $m=M$ and perform another series expansion for the remainder. This procedure can be repeated until the desired precision is reached..

Let us now compare the leading order approximation with the numerical solution. To obtain the expansion parameters $\alpha$ and $\beta$, we perform a fit of the leading order expansion to the numerical solution, as explained in appendix A. In the upper panel of Fig. 3 we compare the full leading order approximation $s^{(1)}$ and the truncated series $s^{(1)}_{(M)}$ for different $M$ to the numerical solution. We can see that the full leading order series expansion $s^{(1)}$ converges to the numerical profile for $z>2$ and is already well described by the truncated series with $M=1$.

The lower panel shows the normalized differences between the truncated series expansion and the full leading order solution, $(s^{(1)}_{(M)}-s^{(1)})/s^{(1)}$, and deviation of the full leading order solution compared to the numerical solution. For increasing $z$, the differences between the numerical solution and the leading order solution $s^{(1)}$ become exponentially small. For $z\gtrsim 12$ we can see that the precision of the leading order approximation $s^{(1)}$ starts to be limited by the precision of the expansion parameters $\alpha$ and $\beta$, which we estimated to be of the order $10^{-9}$. The colored dashed lines show the remainder $R^{1}_{M}$ of the truncated series expansion $s^{(1)}_{(M)}$. We can see the for $M=4$ the uncertainty due to truncation is comparable to the uncertainty of the input parameters.

At small $z$, we can see that difference between the numerical solution an $s^{(1)}$ increase and higher order terms $N>1$ start to be important.

IV.2 Next to Leading Order Contributions

In the previous section we have analyzed the leading order $N=1$ contribution of the series expansion in Eq. 20 to the wavefunction. We found that at large $z\gtrsim 12$, the contribution from next-to-leading order terms $N>1$ is smaller than the uncertainty induced by the uncertainty of the parameter $\alpha$ and $\beta$. We concluded that in this range the $N>1$ terms can be safely ignored. However, at intermediate $z$ in the range $2<z<12$, the next to leading order terms become important, as we have seen in Fig. 3 and terms of higher order in $N$ have to be included. Let us introduce the truncated solution $s^{(N)}_{(M)}$ and the corresponding remainder $R^{N}_{M}$ via

[TABLE]

In the upper panel of Fig. 4 we show the truncated expansion $s^{(N)}_{(M)}$ for the wavefunction for different choices of $N$, fixing $M=4$. The lower panel shows the corresponding deviations of the truncated expansion from the numerical solution. We can see that the series always diverges for $z<1$. Including the $N=3$ and $N=5$ terms will increase the accuracy of our series expansion for an intermediate radius $z\approx 3$ to an $\mathcal{O}(10^{-5})$ level. Taking into account additional terms $N>5$ does not significantly increase the accuracy of the expansion. In this case, the dominant contribution to the remainder comes from neglected terms with $m>M$, which become more important at small $z$.

IV.3 Combined Result and Matching

We can now combine the two truncated solutions, $s_{(N)}$ and $s^{(N)}_{(M)}$, obtained for both small and large values of $z$ by matching them at a matching point $z^{*}$. This is shown for $s_{(N)}$ with $N=10$ and $s^{(N)}_{(M)}$ with $N=1$, $M=1$ in the upper panel of Fig. 5. Here the truncated solution takes the simple form

[TABLE]

We can see that already such few terms in the series expansion are sufficient to describe the wavefunction well. The accuracy is at the few percent level at the matching point $z=2.5$ and orders of magnitude better at small and large $z$, as shown in the lower panel. A precision of $10^{-5}$ at the matching point can be achieved using $s_{(N)}$ with $N=50$ and $s^{(N)}_{(M)}$ with $N=5$, $M=4$, where the precision at small and large $z$ is again limited by the precision of the expansion parameters $s_{0}$, $v_{0}$, $\alpha$ and $\beta$ at the $10^{-12}$ level.

We can also determine the expansion parameters by matching the small and large radius wavefunction and their derivatives at a matching point $z^{*}$. We have performed such a matching using $s_{(N)}$ with $N=250$ and $s^{(N)}_{(M)}$ with $N=10$ and $M=6$ and obtain

[TABLE]

To estimate the uncertainty associated with the matching procedure, we performed multiple matchings for $3<z^{*}<3.5$. The obtained values for expansion parameters are consistent with those in Eq. 32 obtained through the fitting of the large radius solution but have a significantly worse precision. This is not surprising, since the precision of our series expansion is expected to be worst at the matching point $z^{*}$.

V Conclusion

We have found a semi-analytic solution to the coupled Poisson-Newton equations describing dilute boson stars. We have shown that our solution is stable to numerical errors, and that it reproduces the numerical results with accuracy better than $10^{-5}$ over the entire range. Improvements in accuracy can easily be attained for a small expense of numerical work.

There are many possible applications of our methods. The simplest one is to consider interacting bosons, when a potential for the bosons is added. Such a potential can significantly modify the solution, because the interactions can be much stronger than gravitational. Rotation can also modify the solution. In all these case, the large number of parameters makes it impractical to find purely numerical solutions; our semi-analytical method would be better suited for these problems.

Another open question is the stability of these solutions. For example, it is not known how the boson stars are affected by external perturbations e.g. by another star nearby. Once again, an accurate knowledge of the profiles is a requirement for the stability analysis. We hope to return to these and other questions in future work.

Acknowledgements.

This work is supported by NSF under Grant PHY-1620638.

Appendix A Numerical Integration and Fitting

To obtain a numerical solution, it is convenient to use the variables $S$, $V$ and coordinate $x$ given in Eq. 11 with a reference scale $k$ chosen such that $V(0)=1$. As shown in Tod1992 , the solutions of Eq. 8 can then be parametrized by the central value of the wavefunction, $S_{0}=S(0)$ and categorized into three distinct classes: for $S_{0}>S_{0}^{*}$ the wavefunction diverges for at large radius towards positive infinity, for $S_{0}=S_{0}^{*}$ the wavefunction converges to zero, is positive definite and square integrable, while for $S_{0}<S_{0}^{*}$ the wavefunction diverges for at large radius towards negative infinity.

Using a Runge-Kutta 4 method with constant step size $\Delta_{x}$, we perform the numerical integration until the wavefunction starts to diverge and iteratively optimize the central value of the wavefunction $S_{0}$ to find $S_{0}^{*}$. The precision of the wavefunction needed for the numerical solution to stay finite until a large value of $x$, which is needed to fit the large range solution, increases exponentially with the radial coordinate $x$. In this study we use a precision of 35 significant figures for $S_{0}$, providing a converging numerical solution for $x<35$. The accuracy of the numerical solution is limited by the step size $\Delta_{x}$. In this study, we use $\Delta_{x}=10^{-3}$, providing an accuracy of the solution of order $\mathcal{O}\left(\Delta_{x}^{4}\right)\approx 10^{-12}$.

To obtain the expansion parameters $\alpha,\beta$, we fit the Newtonian potential $V(x)=-k^{2}+\frac{2k\beta}{x}$ and the Whittaker solution $S(x)=\frac{k\alpha}{2^{\beta}x}W_{\beta,-\frac{1}{2}}(2kx)$ to the numerical solution for $V$ and $S$ at large $x$. To avoid systematic effects due to the truncation of subleading terms $n>1$ of the series expansion in Eq. 17, we restrict the fitting range to $x>x^{*}$, where the fraction of mass outside radius $x^{*}$ contributed less than $10^{-12}$ to the total mass of the boson star. The expansion parameters at small radius can then be obtained through $s_{0}=k^{-2}S_{0}$ and $v_{0}=k^{-2}$. We find that

[TABLE]

The uncertainty of $s_{0},v_{0}$ and $\beta$ were estimated by the difference of the best fit values to potential and wavefunction. For $\alpha$ we state the uncertainty of wavefunction fit.

Appendix B Representations of the Leading Order Wavefunction

When considering only the potential $v^{(0)}=-1+\frac{2\beta}{z}$, Eq. 12 can be written as

[TABLE]

or

[TABLE]

Here we performed a change of variables to $y=2z$ and $w=2z\,s^{(1)}$. This is the Whittaker equation Abramowitz:1974 , which has the solution555Note that in principle there are three additional solutions: (1) $w\sim M_{\beta,\frac{1}{2}}(z)$ which diverges at $z\to\infty$, (2) $w\sim M_{\beta,-\frac{1}{2}}(z)$ which is undefined for every value of $\beta,z$ and (3) $w\sim W_{\beta,\frac{1}{2}}(z)=W_{\beta,-\frac{1}{2}}(z)$. $w=C\cdot W_{\beta,-\frac{1}{2}}(y)$ and therefore implies $s^{(1)}=C\cdot(2z)^{-1}W_{\beta,-\frac{1}{2}}(2z)$. Note that the Whittaker function is related to the confluent hypergeometric function $U$ through $W_{\kappa,\mu}(z)=e^{-\frac{1}{2}z}z^{\frac{1}{2}+\mu}U(\frac{1}{2}+\mu-\kappa,1+2\mu,z)$ which was used in Eq. 27. We can expand the Whittaker function in an infinite series

[TABLE]

This is the same solution as Eq. 26 and we can identify $\alpha=C2^{\beta-1}$.

The far-field wavefunction can also be written in terms of the Bateman function $k_{\nu}(z)$

[TABLE]

The Bateman function Batemann1931 can be expressed via the integral form

[TABLE]

which is finite and convergent.