[1,2]Kuldeep J. Purohit

These authors contributed equally to this work.

These authors contributed equally to this work.

These authors contributed equally to this work.

[1]Physics Department, Faculty of Science, The Maharaja Sayajirao University
of Baroda, Vadodara, 390 002, Gujarat, India

2]Department of Physics, St. Xavier’s College, Ahmedabad, 380 009, Gujarat,
India

3]Theoretical Physics Division, Physical Research Laboratory, Ahmedabad, 380
009, Gujarat, India

4]Department of Physics, Indian Institute of Technology Gandhinagar, Palaj,
382 355, Gujarat, India

5]School of Fundamental Physics and Mathematical Sciences, Hangzhou Institute
for Advanced Study, UCAS, Hangzhou, 310 024, China

6]University of Chinese Academy of Sciences, Beijing, 100 190, China

# Formation of a Bose Star in a Rotating Cloud

<EMAIL_ADDRESS>Pravin Kumar Natwariya<EMAIL_ADDRESS>Jitesh
R. Bhatt<EMAIL_ADDRESS>Prashant K. Mehta<EMAIL_ADDRESS>*  [
[  [  [  [

###### Abstract

In this paper, we study the evolutions of a self-gravitating cloud of bosonic
dark matter with finite angular momentum and self-interaction. This is
achieved by using the sixth-order pseudospectral operator splitting method to
solve the system of nonlinear Schrödinger and Poisson equations. The initial
cloud is assumed to have mass density randomly distributed throughout three-
dimensional space. The dark matter particles in the initial cloud are in the
kinetic regime, i.e., their de Broglie wavelength is much smaller than the
halo size.

It is shown that Bose stars are indeed formed in the numerical simulation
presented here. The presence of angular momentum and self-interaction in the
initial cloud can significantly influence the star formation time in a non-
trivial fashion.

Furthermore, the plots of vorticity magnitude profile after the star formation
time indicates that the formed star may not have any intrinsic angular
momentum for the cases when the self-interaction among the particles is either
negligible or attractive. These results are in agreement with the earlier
analytical studies of an isolated rotating Bose star. However, for the case of
repulsive self-interaction, the vorticity magnitude analysis show a
possibility that the star formed in the numerical simulations may possess
intrinsic angular momentum. It is also shown that the average mass and radius
diagrams of the star are strongly influenced by the presence of angular
momentum in the initial cloud.

###### keywords:

Bose Stars, Star Formation, Ultra Light Dark Matter, Bose-Einstein
Condensation

##  1 Introduction

Although the idea of Bose star dates back to several decades, the recent
advances in the observational astronomy can provide an unique opportunity for
detecting them. For example, the detection of the gravitational-wave signal
GW190521 by the Advanced LIGO/Virgo collaboration [1] could be due to the
merger of two boson stars [3].The idea of a Bose star is closely related to
the concept of geons, first introduced by Wheeler in 1955 [31]. Geons are
considered to be macroscopic confined states of electromagnetic or
gravitational radiation, held together by the gravitational field created by
their own energy density. The notion of a boson star, formed by the scalar
degree of freedom present in fundamental theories of physics undergoing
gravitational collapse, was first introduced by Kaup in 1968 [11] and Ruffini
and Bonazzola in 1969 [19]. Such a star can be regarded as a macroscopic state
created by the balance between the quantum pressure of the scalar field due to
the Heisenberg uncertainty principle and the force of self-gravity. There has
been significant progress in finding solutions of the Einstein-Klein-Gordon
equations for various situations in the literature [29]. According to the
standard model of cosmology, known as the $\Lambda$CDM model, cold dark matter
constitutes approximately 83% of the total matter in the Universe. Various
candidates for dark matter, such as the QCD-axion [16, 12, 27, 9] and fuzzy
cold dark matter [10], are described by the scalar degree of freedom. Hence,
it is plausible that the Universe contains numerous Bose star-like objects.
Moreover, gravitational bound state structures formed by bosonic dark matter
may help in explaining the missing satellite problem faced by $\Lambda$CDM
model[for example see [15]]. There has been several studies that explore the
connection between Bose stars and phenomena like fast radio bursts or the
detection of axion stars [14, 17, 18, 2].

It should be noted that in the literature, not much attention has been paid to
the issue of boson star formation under realistic scenarios. However, some
interesting works have been published in recent years that focus on this issue
[15, 5, 28]. In [15], the authors study the gravitational collapse of a
randomly distributed, non-interacting cloud of scalar bosons in the kinetic
(low-density) regime. The kinetic regime is characterized by the de Broglie
wavelength $(mv)^{-1}$ and time scale $mv^{2}$ of the particles, which satisfy
the following conditions:

$mvR\gg 1\quad\text{and}\quad mv^{2}\tau_{\text{gr}}\gg 1$ (1)

Here, $m$ is the mass of the dark matter particle and $v$ is its typical
velocity. $R$ represents the halo size, and $\tau_{\text{gr}}$ is the
condensation time at which the star formation process begins. The authors
estimate the value of $\tau_{\text{gr}}$ using a heuristic argument. After the
time $\tau_{\text{gr}}$, numerical simulations show that the time evolution of
the scalar field exhibits strong fluctuations during star formation, and the
mass of the star continues to grow on average. The authors argue that their
study can also help in solving the missing satellite problem [15].

The role of self-interaction among the bosonic dark matter particles in the
kinetic regime was analyzed in [5]. This study shows that self-interaction in
the dark matter cloud can influence the mass and compactness of the Bose star.
Simulation of ultralight dark matter halo formation has been carried out in
[28]. The study demonstrates a solitonic core and confirms the core-halo mass
relations.

However, it should be emphasized that gravitational shearing can naturally
induce rotation in star-forming clouds. For instance, velocity gradients were
observed across the molecular clouds of the M33 galaxy, revealing their
rotation [13, 24]. Thus, considering the ubiquity of rotation in star-forming
clouds, we investigate the formation and evolution of a Bose star when the
initial dark matter distribution possesses finite angular momentum.
Additionally, we explore the role of self-interaction in the star formation
process while the initial dark matter density is randomly distributed
throughout space, satisfying the conditions associated with the kinetic
regime. In this work, we consider the self-interacting scalar dark matter
field described by a Gross-Pitaevskii-type equation. First, we would like to
note that there are existing studies in the literature that discuss related
issues. In the reference [20], numerical solutions of the Einstein-Klein-
Gordon and Einstein-Proca equations are considered to investigate the
formation process of a Bose star. The initial field density profiles are
examined with and without axisymmetry. It is demonstrated that in such cases,
the scalar star is always unstable and loses all its angular momentum,
resulting in the formation of a Saturn-like ring around the scalar star.
However, the authors do not take into account any self-interaction. Moreover,
in the Proca case, it was found that it might be possible to have a rotating
stable star.

In the references [7, 25], the authors consider the initial configuration as
an isolated boson star with finite angular momentum and study their stability.
The authors show that rotating Bose stars are always unstable when the self-
interaction between their constituent particles is either negligible or
attractive, leading to the loss of all angular momentum. Interestingly, when
the self-interaction is repulsive, it is possible to have a rotating star. In
reference [7], an analytical expression for the growth rate of the instability
was obtained. The dark matter field in this reference is described by a Gross-
Pitaevskii-type equation. On the other hand, the reference [25] considers the
Klein-Gordon field.

Therefore, we believe that the approach chosen in our present paper will
provide useful insights into the formation process of Bose stars.
Additionally, it will allow us to compare our results with the related works
discussed here.

This work is divided into the following sections: In Section 2, we discuss the
time evolution of boson particles for various cases, including
attractive/repulsive self-interaction and angular momentum. In Section 3, we
present the results for the aforementioned cases. Finally, we summarize and
conclude the results in Section 4. Throughout the paper, we use natural units
where $c$, $k_{B}$, and $\hbar$ are set to 1. Here, $c$ represents the speed
of light, $k_{B}$ is the Boltzmann constant, and $\hbar$ is the reduced Planck
constant. Boldface quantities indicate three-dimensional vectors.

##  2 Evolution of the Bose star

In the non-relativistic limit and for a high mean occupation number, the
system can be represented by a classical scalar field, $\psi(\bm{r},t)$, whose
time evolution can be described by the Gross-Pitaevskii-Poisson equations [5,
4, 8] :

$\displaystyle i\frac{\partial\psi}{\partial t}$
$\displaystyle=-\frac{1}{2m}\nabla^{2}\psi+m\Phi\psi+g|\psi|^{2}\psi,$ (2)
$\displaystyle\nabla^{2}\Phi$ $\displaystyle=4\pi Gm(|\psi|^{2}-n),$ (3)

Here, $\Phi\equiv\Phi(\bm{r},t)$ is the gravitational potential, $G$ is the
universal gravitational constant, and $n$ is the mean particle density of the
boson gas. The parameter $g$ represents the self-interaction coupling
constant, where positive values correspond to repulsive self-interaction and
negative values correspond to attractive self-interaction. The case $g=0$
corresponds to a non-interacting scenario. For a box with a size $L$ and a
total number of particles $N$, the number density is given by $n=N/L^{3}$. The
above equations (2) and (3) can be written in dimensionless form using the
following scaling: $r=(1/mv_{0})\tilde{r}$, $t=(1/mv_{0}^{2})\tilde{t}$,
$\Phi=v_{0}^{2}\tilde{\Phi}$, $\psi=(v_{0}^{2}\sqrt{m/4\pi G})\tilde{\psi}$,
and $g=(4\pi G/v_{0}^{2})\tilde{g}$:

$\displaystyle i\frac{\partial}{\partial\tilde{t}}\tilde{\psi}$
$\displaystyle=-\frac{1}{2}\tilde{\nabla}^{2}\tilde{\psi}+\tilde{\Phi}\tilde{\psi}+\tilde{g}|\tilde{\psi}|^{2}\tilde{\psi},$
(4) $\displaystyle\tilde{\nabla}^{2}\tilde{\Phi}$
$\displaystyle=|\tilde{\psi}|^{2}-\tilde{n},$ (5)

Here, $\tilde{n}=\dfrac{4\pi G}{mv_{0}^{4}}n$, and $N=\int d^{3}r|\psi|^{2}$
can be expressed as $N=\dfrac{v_{0}}{4\pi Gm^{2}}\tilde{N}$, where $v_{0}$ is
the reference velocity [5].

Next, numerical solutions of the above equations are obtained using the
pseudospectral operator splitting method. This is done by utilizing the code
AxioNyx111https://github.com/axionyx [23]. However, we needed to modify the
code to incorporate the effect of self-interaction in the method.

To proceed further, we need to specify the initial condition for
$\tilde{\psi}$ at time $\tilde{t}=0$. First, consider the case when there is
no angular momentum carried by the function $\tilde{\psi}$. In this case, the
wave function is assumed to be ’Gaussian’ and can be written in Fourier space
as:

$\tilde{\psi}_{p}=\frac{\tilde{N}^{1/2}}{\pi^{3/4}}e^{-\tilde{p}^{2}/2}e^{i\alpha_{p}},$
(6)

Here, $\alpha_{p}$ represents a random phase distributed over the Fourier
space between the numbers 0 and $2\pi$. To generate a randomly distributed
matter density distribution in space, we transform to the coordinate space to
obtain the profile of amplitude
$\tilde{\psi}(\tilde{x},\tilde{y},\tilde{z},\tilde{t}=0)$. This prescription
is the same as that considered in reference [15].

The angular momentum in the initial wavefunction is incorporated using the
following prescription: First, consider the function

$\displaystyle\tilde{f}(\tilde{\bm{r}},\tilde{t}=0)=e^{-\tilde{\bm{r}}^{2}/2}e^{il\tan^{-1}(\tilde{y}/\tilde{x})},$
(7)

This function is then transformed to Fourier space to obtain
$\tilde{f}(\bm{p})$. The initial wavefunction in momentum space can be
obtained as $\tilde{\psi}(\bm{p})=\tilde{f}(\bm{p})\times e^{i\alpha_{p}}$,
where $\alpha_{p}$ is the random phase as defined above. The input
wavefunction $\tilde{\psi}(\tilde{x},\tilde{y},\tilde{z},\tilde{t}=0)$ is
generated by transforming $\tilde{\psi}(\bm{p})$ back to coordinate space. For
all the numerical results presented in this paper, we use the aforementioned
prescription. It is possible to modify the prescription described by equation
(7); in fact, we have tried a couple of different prescriptions to introduce
angular momentum in $\tilde{\psi}(\tilde{x},\tilde{y},\tilde{z},\tilde{t}=0)$,
but it does not change the value of $\tilde{\tau}_{\text{gr}}$ for a given
value of angular momentum. The qualitative behavior of the numerical results
remains the same. Fig.(1) shows the initial density profile
$|\tilde{\psi}(\tilde{x},\tilde{y},\tilde{z},\tilde{t}=0)|^{2}$ for different
values of angular momentum in different planes.

---

---

Figure 1: The initial snapshots of density profiles of $|\tilde{\psi}|^{2}$ at
the initial time ($\tilde{t}=0$) in $\tilde{y}\tilde{z}$, $\tilde{z}\tilde{x}$
and $\tilde{x}\tilde{y}$ planes. The top panel from left to right describes
the initial cloud density for $\mathcal{\tilde{L}_{\rm tot}}=0$. The bottom
panel respectively correspond to $\mathcal{\tilde{L}_{\rm tot}}=5$ case.
Figure 2: Snapshots of $|\tilde{\psi}|^{2}$ at different values of time are
shown in $\tilde{z}\tilde{y}$-plane. All the plots represent the case when the
initial cloud has total-angular momentum $\mathcal{\tilde{L}_{\rm tot}}=5$ and
no self-interaction $\tilde{g}=0$. The gravitational condensation time
$\tilde{\tau}_{gr}$ for this case is around 15600. The plots in the upper
panel from the left to right respectively correspond to time
$\tilde{t}=0,\,0.5\tau_{gr},\,\&1\tau_{gr}$. Three plots in the bottom
correspond to
$1.15\tilde{\tau}_{gr},\,1.5\tilde{\tau}_{gr}\,\&1.8\tilde{\tau}_{gr}$.

To calculate the total angular momentum carried by the initial scalar field,
we consider the stress-energy tensor [30]:

$\displaystyle\widetilde{T}^{0i}=\frac{i}{2}\left(\tilde{\psi}\tilde{\partial}_{i}\tilde{\psi}^{*}-\tilde{\psi}^{*}\tilde{\partial}_{i}\tilde{\psi}\right).$
(8)

The total angular momentum in the $\tilde{x}$, $\tilde{y}$, and $\tilde{z}$
directions, in terms of the stress-energy tensor, can be defined as [30]:

$\displaystyle\tilde{J}^{ij}=\int(\tilde{x}^{i}\widetilde{T}^{j0}-\tilde{x}^{j}\widetilde{T}^{i0})d^{3}\tilde{x}.$
(9)

For the considered initial distribution in equation (7), the total angular
momentum in the $\tilde{x}$ and $\tilde{y}$ directions ($\tilde{J}^{2,3}$ and
$\tilde{J}^{3,1}$, respectively) is zero. Therefore, the total angular
momentum of the system is
$\mathcal{\tilde{L}_{\text{tot}}}=\tilde{J}^{1,2}=l\tilde{N}$. It is important
to note that once a star is nucleated, it will be performing random motion.
Therefore, even in the case with $l=0$, the star will have some orbital
angular momentum. In cases when the initial cloud has finite angular momentum,
it is hard to study angular momentum transfer from the cloud to the star
because of the presence of orbital and also the possibility of having
intrinsic angular momentum. In this case, we find it useful to study how the
vorticity magnitude is distributed outside the star formation region. However,
this approach does not allow us to quantify angular momentum transfer from the
surrounding cloud to the star. To determine the vorticity of the system, we
introduce the ’velocity’ ${\bf v}=\bm{\mathcal{J}}/\rho$, where
$\bm{\mathcal{J}}$ represents the current density and $\rho=|\psi|^{2}$. The
current density $\bm{\mathcal{J}}$ is defined from equation (1) as
$\bm{\mathcal{J}}=\frac{i}{2}(\psi{\bm{\nabla}}\psi^{*}-\psi^{*}{\bm{\nabla}}\psi)$.
From this, one can define the vorticity as
${\bf\Omega}={\bm{\nabla}}\times{\bf v}$, and the vorticity magnitude as
$\Omega=\sqrt{\Omega_{x}^{2}+\Omega_{y}^{2}+\Omega_{z}^{2}}$.

In the subsequent discussion, the magnitude of vorticity in different planes
has been used to understand how the angular momentum is distributed throughout
the space during the process of star formation.

It is important to emphasize that the presence of attractive self-interaction
can also allow for the possibility of forming bound states without the self-
gravity playing a significant role [5]. Thus, for a gravitational bound state
to be produced in the numerical solution, it is expected that the self-
coupling term remains small or comparable to the self-gravity term.
Furthermore, the self-interaction term can influence the typical size of the
structure produced by self-gravity, even when there is no angular momentum in
the rotating cloud. In the kinetic regime, one may not be able to use the
hydrodynamic formalism to estimate the Jeans length. Therefore, in the present
case, the following arguments can be used to estimate the typical size of the
gravitational structure and provide a bound on the coupling constant.

If we replace all spatial derivative terms in the GPP equations (i.e.,
equations (3) and (4)) with $1/\tilde{L}$, where $\tilde{L}$ is the length
scale, we get:

$i\frac{\partial\tilde{\psi}}{\partial
t}=\left[\frac{1}{2\tilde{L}^{4}}-\left(|\tilde{\psi}|^{2}-\tilde{n}\right)+\tilde{g}\frac{|\tilde{\psi}|^{2}}{\tilde{L}^{2}}\right]\tilde{\psi}\tilde{L}^{2}.$
(10)

The system would collapse under gravity. If there are no particles in the
surroundings, then a stationary state can be reached, i.e.,
$i\frac{\partial\tilde{\psi}}{\partial t}\sim\omega\tilde{\psi}$. By
relabeling $\tilde{L}$ with the symbol $\tilde{L}_{B}$, where the subscript
$B$ indicates a bound state, we have:

$\omega\tilde{\psi}=\left[\frac{1}{2\tilde{L}_{B}^{4}}-\left(|\tilde{\psi}|^{2}-\tilde{n}\right)+\tilde{g}\frac{|\tilde{\psi}|^{2}}{\tilde{L}_{B}^{2}}\right]\tilde{\psi}\tilde{L}_{B}^{2}.$
(11)

Here, the term with $\tilde{n}$ is neglected in comparison with the
$|\tilde{\psi}|^{2}$ term. The first and second terms in the square bracket
correspond to the quantum pressure and self-gravity effects, respectively,
while the third term is the contribution from self-interaction. For a
gravitationally bound state to be formed in the collapse, we require that the
second term in the square bracket remain larger than the term with self-
coupling, leading to the following condition:

$|\tilde{g}\tilde{L}_{B}^{-2}|<1.$ (12)

Thus, we can estimate the size of the bound structure by solving the
biquadratic equation:

$\tilde{L}_{B}^{2}=-\frac{\omega-\tilde{g}|\tilde{\psi}|^{2}}{2|\tilde{\psi}|^{2}}+\sqrt{\frac{(\omega-\tilde{g}|\tilde{\psi}|^{2})^{2}+2|\tilde{\psi}|^{2}}{4|\tilde{\psi}|^{4}}},$
(13)

where we have ignored the negative square root solution as it could result in
$\tilde{L}_{B}^{2}<0$. Since $\omega<0$ for gravitationally bound states, from
the above equation, it can be seen that for attractive self-interaction (i.e.,
$\tilde{g}<0$), the formed structure is smaller in size compared to the cases
with $\tilde{g}>0$ and $\tilde{g}=0$.

The restriction on the values of the coupling constant can also be obtained
using a different set of arguments [6]. Here, the timescale
$\tau_{\text{self}}$ for obtaining a bound state without self-gravity is given
by, $\tau_{\text{self}}=\dfrac{4\sqrt{2}mv^{2}}{3n^{2}g^{2}\pi}$, when there
is an attractive interaction between the constituent particles of the initial
cloud. For the case when both self-interaction and self-gravity are present,
their combined effect can give rise to a new timescale
$\tilde{\tau}_{\text{eff}}$ for structure formation, which is given by:

$\tilde{\tau}_{\text{eff}}=\frac{{\tilde{\tau}}_{\text{gr}}\tilde{\tau}_{\text{self}}}{\tilde{\tau}_{\text{gr}}+\tilde{\tau}_{\text{self}}}.$
(14)

Thus, when the condition $\tilde{\tau}_{\text{gr}}<\tilde{\tau}_{\text{self}}$
is satisfied, self-gravity can play a dominant role in structure formation.

##  3 Results and Discussion

The density profiles of the initial wave-functions in the $yz$, $zx$, and $xy$
planes for various values of angular momentum $\mathcal{L}_{\text{tot}}=0$ are
shown in Figure (1). In the absence of angular momentum or self-coupling in
the initial boson cloud, the only important timescale is the gravitational
condensation time described by the formula given in [15]:

$\tau_{\text{gr}}=\frac{b\sqrt{2}}{12\pi^{3}}\frac{mv_{0}^{6}}{G^{2}n^{2}\Lambda}.$
(15)

where $n$ is the boson number density and $\Lambda=\log(mv_{0}L)$, with
$\tilde{L}=mv_{0}L$ being the box length used in the numerical simulation. It
should be noted that the complete determination of the parameter $b$ is not
possible, but it can be of the order of unity [15]. In the subsequent section,
we will discuss the influence of angular momentum and self-interaction on the
parameter $b$. Using dimensionless units, we can write the gravitational
condensation time as:

$\tilde{\tau}_{\text{gr}}=\frac{b\sqrt{2}}{12\pi^{3}}\frac{1}{\tilde{n}^{2}\log(\tilde{L})}..$
(16)

To ensure the validity of the numerical results obtained using the modified
code, we have employed the following checks:

1) In the absence of angular momentum and self-interaction, the code produces
the results from [15] and [23], i.e., the estimated value of
$\tilde{\tau}_{\text{gr}}$ is in agreement with the numerical values obtained
by the code. Also, the power spectrum of the wave-function $F(\omega,t)$
(defined below) develops a peak at the negative frequency value after
$\tilde{t}>\tilde{\tau}_{\text{gr}}$.

2) The total energy of the system is conserved for the system of equations (4)
& (5). To check energy conservation, the total energy is calculated for the
initial data and compared with the total energy at each subsequent step. If at
any time the total energy of the system exceeds the initial energy by more
than 5%, we do not run the code any further.

Next, keeping the available computational resources in mind, we choose
$\tilde{N}=10$, $\tilde{L}=30$, and a resolution of $\tilde{L}^{3}/128^{3}$.
The typical value of $\tilde{\tau}_{\text{gr}}$ is around $8\times 10^{3}$,
and further increasing the resolution can significantly increase computational
time. Also, including repulsive interactions or higher angular momentum can
further increase the condensation time $\tilde{\tau}_{\text{gr}}$. It should
be noted that for this resolution, the continuum limit can be achieved, as
shown in [15](also see [5]). Here, the authors use the exceptionally stable
sixth order pseudospectral operator splitting method which is suitable for a
long statistical simulations. In this work we employ the same technique to
study time evolution of a self-gravitating cloud of a dark-matter in the
kinetic regime.

For the chosen values of $\tilde{N}$ and $\tilde{L}$, the highest allowed
value of $|\tilde{g}|$ for the formation of a gravitationally bound state is
4.56, as given by the bound condition (12). We have checked if the numerical
results show the process of star formation for $|\tilde{g}|>4.56$ ; however,
we do not find star formation occurring for either attractive or repulsive
self-interaction within the time allowed by our energy conservation test.
Moreover, if we reduce the highest values of $|\tilde{g}|$ by an order of
magnitude, the numerical solutions are very similar to the case with
$\tilde{g}=0$. Thus, in this work, we present the results for the following
three values of the coupling constant: $\tilde{g}=-4.56$, 0.0, and +4.56. The
data generated by the code is analyzed using The yt Project [26].

Here, we would like to note that recently in [15] and [25], the authors
studied the stability of an isolated rotating Bose star. It was demonstrated
that when the self-interaction among the constituent particles is attractive
or negligibly small, the star becomes unstable for any value of angular
momentum. The instability results in the shedding of particles (including
angular momentum) in the region outside the star, leading to the formation of
a Saturn ring-like structure around the star. The time scale associated with
the instability of a rotating boson star [7] was estimated using the growth
rate given by:

$\tau_{\text{ins}}=\frac{100}{2.2}\frac{l^{2}}{\alpha_{l}}\frac{1}{m^{3}G^{2}M_{s}^{2}}..$
(17)

where $l$, $M_{s}$, and $\alpha_{l}$ represent the values of angular momentum,
the mass of the star, and an angular-momentum dependent parameter,
respectively. However, the authors demonstrate that when the self-interaction
between the constituent particles is repulsive, it is possible to have a
stable Bose star configuration.

Thus, it would be rather interesting to see whether the Bose stars formed in
our numerical results are rotating or not. First, we would like to emphasize
that in the stability analysis by [7], the initial configuration is assumed to
be an isolated boson star. In our work, star formation occurs after a time
$\tilde{\tau}_{\text{gr}}$ from the initial state with a randomly distributed
density over the entire space. To understand how the instability influences
the presented numerical results, let us consider the ratio
${\tilde{\tau}_{\text{ins}}}/{\tilde{\tau}_{\text{gr}}}$. Using the values of
$\tilde{n}$ and $\tilde{L}$, the ratio
${\tilde{\tau}_{\text{ins}}}/{\tilde{\tau}_{\text{gr}}}\sim 10^{-2}$.
Furthermore, it should be noted that we have used the value of
$\tilde{\tau}_{\text{gr}}\sim 8000$. The typical range of
$\tilde{\tau}_{\text{gr}}$ is around 8000-25000 for $\tilde{g}\geq 0$. Thus,
one can only decrease ${\tilde{\tau}_{\text{ins}}}/{\tilde{\tau}_{\text{gr}}}$
by increasing $\tilde{\tau}_{\text{gr}}$. Since
${\tilde{\tau}_{\text{ins}}}/{\tilde{\tau}_{\text{gr}}}\leq 10^{-2}$, we
believe that any possible instability would have time to saturate.

To check when the star condensation begins, we employ the following three
tests:

  * •

The first test involves studying the time evolution of the maximum amplitude
$|\tilde{\psi}(\tilde{x},\tilde{t})|$. The typical value for the initial
maximum amplitude is $|\tilde{\psi}(\tilde{x},\tilde{t})|\sim 2\times
10^{-2}$. The sudden rise in the maximum amplitude could be due to gravity or
the self-interaction producing a bound state. This provides us with an
estimate of $\tilde{\tau}_{\text{eff}}$ for $\tilde{g}\leq 0$ or
$\tilde{\tau}_{\text{gr}}$ when $\tilde{g}>0$ for repulsive interaction.

  * •

The second test involves studying the power spectrum $F(t,\omega)$ of
$\psi(x,t)$ [15], defined as the Fourier transform of the correlator:

$F(\omega,t)=\int\frac{dt_{1}}{2\pi}d^{3}\bm{x}\
\psi^{*}(t,\bm{x}),\psi(t+t_{1},\bm{x}),e^{i\omega
t_{1}-t_{1}^{2}/\tau_{1}^{2}}$ (18)

Around the time when gravitational condensation occurs, $F(\omega,t)$ starts
developing a peak at $\omega_{s}<0$, indicating the formation of a
gravitationally bound state.

  * •

The third test involves checking if the formed structure is gravitationally
bound by satisfying the condition $|\tilde{g}\tilde{L}_{B}^{-2}|<1$.

These tests allow us to determine the onset of star condensation and verify
the formation of a gravitationally bound structure.

Figure (3(a)) shows the evolution of the maximum density,
$\max|\tilde{\psi}|^{2}$, with time $\tilde{t}$. The plots from top to bottom
correspond to cases with $\mathcal{\tilde{L}_{\text{tot}}}=0$, $3$, and $5$.
Initially, $\max|\tilde{\psi}|$ remains constant and small for a long time
$\tilde{t}\sim\tilde{\tau}_{\text{gr}}$, and then there is a sudden rise in
its value. It should be noted that the mean value of density continues to rise
monotonically for $\tilde{t}>\tilde{\tau}_{\text{gr}}$ until the convergence
test stops the code around $\tilde{t}\sim 21000$. Ideally, the growth in
density should stop at some point, but due to computational limitations, it is
difficult to observe this in the numerical computation. Previous studies
related to Bose star formation, such as [15] and [5], have also not observed
saturation in density in the kinetic regime. Additionally, the inclusion of
angular momentum in the initial amplitude does not help in achieving
saturation.

Figure (3(b)) describes the behavior of the spectrum
$\tilde{F}(\tilde{\omega},\tilde{t})$ as a function of frequency at
$\tilde{t}\sim 21000$, which is larger than the condensation time
$\tilde{\tau}{\text{gr}}$ for all three cases. The figure clearly shows the
presence of a prominent peak in the region where the frequency has negative
values for all different values of $\mathcal{\tilde{L}}_{\text{tot}}$. It
should be emphasized that this peak in $\tilde{F}(\tilde{\omega},\tilde{t})$
occurs soon after $\tilde{\tau}_{\text{gr}}$ and persists as long as the
convergence test remains valid. This peak indicates the presence of a
gravitationally bound state or a Bose star.

(a)

(b)

Figure 3: Figure 3(a) represents the time evolution of the maximum density
$|\psi|^{2}$ for different values of initial angular momentum of the bosonic
field for $\tilde{g}=0$. Figure 3(b) represents the energy spectrum,
$\tilde{F}(\tilde{\omega})$, as a function of $\tilde{\omega}$. In both
figures, $\mathcal{\tilde{L}_{\rm tot}}$ varies from top to bottom as 0, 3,
and 5. In figure 3(b), all plots are obtained at $\tilde{t}\sim 21000$.

To gain insight into the fate of the initial cloud’s angular momentum, we have
shown three snapshots of the density profile $|\tilde{\psi}|^{2}$ along with
the vorticity magnitude in the $\tilde{x}\tilde{y}$-plane at a time around
21000 in Figure (4). The vorticity magnitude is represented by the black
contours. The leftmost plot corresponds to the case with
$\mathcal{\tilde{L}}_{\text{tot}}=0$, where there is no initial vorticity. In
the other two plots with non-zero angular momentum values, most of the
vorticity magnitude remains outside the star-forming region. This suggests
that the stars formed in the cases with $\mathcal{\tilde{L}}_{\text{tot}}=3$
and $\mathcal{\tilde{L}}_{\text{tot}}=5$ may not possess intrinsic angular
momentum. These findings align with the stability analysis presented in [7].

The typical time interval for which the star persists after the condensation
time is approximately (1.5-2.0)$\tilde{\tau}_{\text{gr}}$, with the specific
value depending on $\mathcal{\tilde{L}}_{\text{tot}}$ and type of self-
interaction. The instability of a rotating boson star in the absence of self-
interaction is known, but the characteristic time for the instability to grow,
$\dfrac{\tilde{\tau}_{\text{ins}}}{\tilde{\tau}_{\text{gr}}}\sim 10^{-2}$, is
much shorter than (1.5-2.0)$\tilde{\tau}_{\text{gr}}$. Therefore, it is
expected that the instability saturates within the observed time interval.
Additionally, the presence of a peak in $\tilde{F}(\tilde{\omega},\tilde{t})$
at negative frequencies throughout the interval
(1.5-2.0)$\tilde{\tau}_{\text{gr}}$ suggests that the star formed in the case
with $\mathcal{\tilde{L}}_{\text{tot}}=0$ is stable against rotational
instability.

Figure 4:  Three snapshots of density profile $|\tilde{\psi}|^{2}$ together
with vorticity magnitude in $\tilde{x}\tilde{y}-$plane are shown at time
$\tilde{t}\sim$ 21000\. This the case when there is no self-interaction. The
first plot on the left corresponds to the case with no angular momentum
($\mathcal{\tilde{L}_{\rm tot}}=0$), while the second and the the third plots
respectively describe the cases when $\mathcal{\tilde{L}_{\rm tot}}=3$ and 5.
The black contours represent vorticity magnitude. The high density region
inside the redspots corresponds to the Bose star. The second and third plots
shows the presence of vorticity outside the star formation region. In
$\tilde{x}\tilde{z}-$plane & $\tilde{y}\tilde{z}-$planes also vorticity found
to be in the region outside the star.

In the left panel of Figure (3(a)), it can be seen that the maximum density
starts to increase and fluctuate strongly after the condensation time.
Therefore, the plots of total mass of the star $M(\tilde{t})$ versus time
$\tilde{t}$ can exhibit growth, but the data also shows large fluctuations.
Similarly, one would expect to observe fluctuations in the radius versus time
plots as well. The theoretical calculations based on stationary boson stars
predict that as the mass of the star increases, its radius shrinks [4, 22].
This behavior has also been observed in numerical simulations [15].
Furthermore, it has been shown that if the mass and radius of a non-rotating
Bose star increase simultaneously, they can be unstable [22]. In order to
smooth out the fluctuations in the data, we employ a digital filter developed
by Savitzky-Golay [21].Compared to other smoothing techniques(such as moving
average, etc), the Savitzky-Golay filter stands out for its ability to
simultaneously reduce noise and maintain the integrity of fine details present
in the data. By fitting a local polynomial function to a sliding window of
data points, the filter effectively captures the underlying trend while
suppressing high-frequency noise components.

Figure (5) shows the $M-R$ relation, $R_{99}$ measures distance from the
centre at which 99% mass of the star is located. It has been observed in our
simulation that $2R_{99}\,>\,L_{B}$. The red curve represents the case with no
angular momentum in the initial cloud. It can be clearly seen that the star
becomes more and more compact as its mass increases. Similar behavior is
observed for the blue curve representing the case with
$\mathcal{\tilde{L}}_{\text{tot}}=3$. However, for the case with
$\mathcal{\tilde{L}}_{\text{tot}}=5$, there are two branches of solutions. In
the lower part of the green curve, both the mass and radius of the star
increase together. This branch may be unstable [22], but the time scale of the
instability may be much longer than $\tilde{\tau}_{\text{gr}}$, allowing the
star to survive during the simulation time. The lowest mass on this diagram
occurs around the time when the condensation just begins, and with the passage
of time, the star moves towards the stable branch.

Figure 5: Mass–radius relation for Bose stars without any self-interaction for
different values of $\mathcal{\tilde{L}_{\rm tot}}$

###  3.1 Bose star formation in presence of self-interaction

Next, we investigate Bose star formation in the presence of finite self-
interaction. We examine the gravitational collapse of a self-interacting boson
field for various values of the coupling constant $\tilde{g}$. As discussed
earlier, star formation may not be possible for arbitrarily large values of
the coupling strength. For $\tilde{N}=10$ and $\tilde{L}=30$, we determine
that the maximum allowed value of $|\tilde{g}|$ is approximately 4.56.
However, if we decrease the coupling constant by an order of magnitude, the
results approach the case with no self-interaction. Therefore, in this
section, we focus our analysis on a high value of the coupling constant.

####  3.1.1 Attractive self-interaction

In this case, it is possible to form bound structures due to the balance
between the attractive force provided by the self-interaction and the quantum
pressure. Hence, it is crucial to examine the role of self-gravity in the
observed bound state. For this analysis, we consider $\tilde{g}=-4.56$.

We begin by analyzing the plots of maximum density $\max|\tilde{\psi}|^{2}$
versus time for different values of $\mathcal{\tilde{L}}_{\rm tot}$, as shown
in figure (6(a)). All the plots indicate that for a significant duration of
time, $\max|\tilde{\psi}|^{2}$ remains very small and constant. However, after
a certain time when condensation occurs, $\max|\tilde{\psi}|^{2}$ starts to
rise monotonically until the convergence test related to the total energy
stops the process. The specific time at which $\max|\tilde{\psi}|^{2}$ begins
to increase depends on the value of $\mathcal{\tilde{L}}_{\rm tot}$ present in
the initial cloud.

(a)

(b)

Figure 6: The caption is same as in figure (3), except here, we have
considered attractive self-interaction ($\tilde{g}=-4.56$) and all plots are
obtained at $\tilde{t}=16000$.

Next, Figure (6(b)) presents the plots of the spectrum
$\tilde{F}(\tilde{t},\tilde{\omega})$ versus frequency $\tilde{\omega}$ at
time $\tilde{t}\simeq 16000$ for three different values of
$\mathcal{\tilde{L}}_{\rm tot}$. All three plots in Figure (6(b)) demonstrate
the formation of a gravitationally bound state, indicated by the presence of a
peak in the negative frequency regime. As discussed earlier, when the
interaction is attractive, there exists an upper bound on the coupling
constant $\tilde{g}$ given by the condition (12). This estimation requires
determining the value of $\tilde{L}_{B}$. In the case of self-interaction, if
the Bose star is formed with net angular momentum, it could be unstable with
an instability timescale smaller than $\tilde{\tau}_{gr}$.

Figure (7) displays three snapshots of $|\psi|^{2}$ together with the
vorticity magnitude at time $\tilde{t}\sim 16000$ in the $xy$-plane for three
different values of $\mathcal{\tilde{L}}_{\rm tot}$. The vorticity magnitude
is depicted by the black contours, while the bright red spots indicate the
region where the Bose star has formed. The plots reveal that the black
contours are concentrated outside the star-forming region. Thus, in this case,
similar to the $\tilde{g}=0$ case, we believe that the formed star does not
have any intrinsic angular momentum, and our results demonstrate stable stars.
These findings align with the stability analysis presented in [7].

Figure 7: The caption is same as in figure (4), except here, we have
considered attractive self-interaction ($\tilde{g}=-4.56$) and all plots are
obtained at $\tilde{t}=16000$.

Figure (8) illustrates $M-R$ diagrams for different values of the total
angular momentum of the star. The plot is generated using the digital filter
discussed earlier. The red curve correspond to the case when there is no
angular momentum in the cloud. The lowest value of mass corresponds to the
time when the condensation begins. The curve indicates that as the star mass
grows, it becomes more compact. The blue curve exhibits a similar behavior.
The green curve corresponds to the case when the total angular momentum is 5.
In this case, for the lower branch, the mass of the star increases with the
radius, indicating instability. However, for the upper branch, we observe that
increasing the mass of the star leads to a more compact configuration.

Figure 8: Mass–radius relation for Bose stars with Attractive self-interaction
for different values of $\mathcal{\tilde{L}_{\rm tot}}$

####  3.1.2 Repulsive self-interaction

Next, we consider the case of repulsive self-interaction. In this case, it is
reasonable to expect that the bound state formed in the numerical results is
due to gravity. If the coupling constant is increased to values higher than
4.56, bound states can still be formed. However, studying such bound states
would require more computational time than allowed by the convergence test
implemented to validate the numerical results. Therefore, for now, we restrict
our analysis to $\tilde{g}=+4.56$.

First, we analyze the plots of the maximum density $\max|\tilde{\psi}|^{2}$
vs. time, as shown in Figure (9(a)). The topmost plot corresponds to the case
when the total angular momentum of the collapsing cloud is zero. The middle
and bottom plots correspond to the cases with $\mathcal{\tilde{L}}_{\rm
tot}=3$ and $\mathcal{\tilde{L}}_{\rm tot}=5$, respectively. All plots show
that $\max|\tilde{\psi}|^{2}$ remains small and constant until time
$\tilde{t}\sim\tilde{\tau}_{gr}$. After that, $\max|\tilde{\psi}|^{2}$ starts
to increase with time, until the convergence test limits the program runtime.

(a)

(b)

Figure 9: The caption is same as in figure (3), except here, we have
considered repulsive self-interaction ($\tilde{g}=+4.56$) and all plots are
obtained at $\tilde{t}=25500$.

Next, Figure (9(b)) illustrates the plots of the spectrum
$\tilde{F}(\tilde{t},\tilde{\omega})$ versus frequency $\tilde{\omega}$ at
time $\tilde{t}\simeq,25500$ for different values of $\mathcal{\tilde{L}}{\rm
tot}$. The presence of a peak at negative frequency in
$\tilde{F}(\tilde{t},\tilde{\omega})$ indicates the existence of a
gravitationally bound state. The stability analysis in [7] suggests that an
isolated Bose star with angular momentum might be stable. Therefore, it is of
interest to investigate what happens to the net angular momentum of the
collapsing cloud after the star formation.

Figure 10: Vorticity magnitude plotted in $\tilde{x}\tilde{y}$ plane at times
$\tilde{t}=22000,23300\ \textrm{and}\ 27400$ from left to right together with
the Density for $\tilde{g}=+4.56$ (repulsive self-interaction) and
$\mathcal{\tilde{L}_{\rm tot}}\,=\,5.0$. With the passage of time the net
vorticity vanishes from the surrounding of the star formation region. The
plots does not have any vorticity magnitude outside the star forming region.
The absence of vorticity is seen in the other two planes also. This may be
indicative of angular momentum entering the star.

Figure (10) presents three snapshots of density profiles together with the
vorticity magnitude when $\mathcal{\tilde{L}}_{\rm tot}=5$. The condensation
time for $\mathcal{\tilde{L}}_{\rm tot}=5$ is approximately 18000. The
leftmost snapshot, taken at $\tilde{t}\sim 22000$, is closer to
$\tilde{\tau}_{gr}\sim 18000$. In this plot, the black contours represent the
presence of finite vorticity magnitude outside the star-forming region,
represented by the red spot. The middle snapshot, taken at time 23300, shows a
reduction in the black contours around the star-forming region. This reduction
is not limited to the $xy$-plane alone; a significant decrease in vorticity is
also observed outside the star-forming region in the other two planes. The
rightmost snapshot, taken at time 27400, exhibits a near absence of vorticity
magnitude outside the star-forming region. Similarly, vorticity is not
observed outside the star-forming region in the other planes. A similar
behavior is observed for the $\mathcal{\tilde{L}}_{\rm tot}=3$ case as well,
where the vorticity magnitude outside the star-forming region starts
decreasing soon after the gravitational condensation.

Figure 11: Mass–radius relation for Bose stars with Repulsive self-interaction
for different values of $\mathcal{\tilde{L}_{\rm tot}}$

Figure (11) shows the $M-R$ diagram when the self-interaction is repulsive for
different values of angular momentum. In the absence of angular momentum in
the cloud, the red curve clearly indicates that just after the condensation,
the star begins with the lowest mass and, as it becomes more massive with the
passage of time, it becomes more compact. Interestingly, the vorticity plots
after the condensation begins clearly demonstrate that the vorticity vanishes
outside the star-forming region. This suggests that the angular momentum might
have been transferred to the star. In this case, both the blue and green
curves show two branches, where the mass of the star grows with increasing
radius, and vice versa. However, there is currently no analytical or numerical
work, to the best of our knowledge, that can determine whether a rotating star
with repulsive interaction can have a branch of solutions where the mass and
radius of the star grow together. These cases are new and require further
investigation. However, from the timescales relevant to our simulations, such
stars are formed just after the time $\tilde{\tau}_{gr}$ and they survive
until the convergence test allows the code to run. In that limited sense,
stars represented by the blue and green curves are stable.

In Figure (12(a)), we observe the relationship between condensation time
($\tau_{\text{gr}}$) and angular momentum per particle ($l$). Notably,
$\tau_{\text{gr}}$ increases in a nonlinear manner with increasing $l$. This
trend is consistent across cases involving attractive, repulsive, and no self-
interaction. Furthermore, it’s evident that, for a specific value of $l$,
$\tau_{\text{gr}}$ is shorter when there’s an attractive self-interaction
compared to the scenario with no self-interaction. Conversely, when there’s a
repulsive self-interaction, $\tau_{\text{gr}}$ is longer compared to the case
with no self-interaction. This underscores the significant influence of both
angular momentum in the initial field and the type of self-interaction between
particles on the behavior of $\tau_{\text{gr}}$. Figure (12(b)) provides
insights into estimating the parameter $b$ in equation (16) for various values
of $l$. Notably, when $l$ is equal to zero and there’s no self-interaction,
the parameter $b$ is approximately one. However, as we increase the angular
momentum ($l$), the parameter $b$ exhibits a nonlinear increase. Additionally,
the type of self-interaction also plays a role in influencing the behavior of
$b$.

(a) Variation of Condensation Time, ($\tau_{\text{gr}}$), with Angular
Momentum per Particle ($l$) for Different Types of Self-Interaction

(b) Variation of parameter $b$ with Angular Momentum per Particle ($l$) for
Different Types of Self-Interaction

Figure 12:

##  4 Conclusions

In conclusion, this work studies the numerical solutions of the nonlinear
Schrödinger and Poisson equations, which describe the formation and evolution
of a Bose star when the initial cloud has finite angular momentum. The star
formation occurs in all the numerical solutions presented here. The results
show that once a Bose star is formed in our numerical simulation, it persists
until our code violates the energy conservation test.

Since the work is carried out in the kinetic regime, with an initially
randomly distributed dark matter density distribution over the entire box, the
typical condensation time is on the order of $10^{3}$ in dimensionless units.
The numerical results demonstrate that the presence of angular momentum and
self-interaction can significantly affect the star condensation time.
Furthermore, the study shows that the presence of angular momentum and self-
interaction in the star-forming cloud strongly influences the mass and radius
of the Bose star.

Further we would like to note here that the star formed when the initial cloud
does not carry any intrinsic angular momentum performs a random motion during
its growth [15]. This allows the formed star to have ”orbital” angular
momentum. This in turn complicates the important question related with net
angular momentum transferred between the rotating initial cloud and the star
formed in the numerical simulation. To gain insight into this problem we have
studied how vorticity magnitude profile along with the star. The vorticity
magnitude provide a qualitative idea of angular momentum transfer between the
initial cloud and the star. More quantitative analysis of angular momentum
transfer between the initial cloud and the star shall be done in a future
work. The vorticity magnitude plotted along with the density profile for the
cases related with negligible or attractive interaction indicates that the
formed star may not possess any intrinsic angular momentum. But for the case
of repulsive interaction, our study shows the absence of vorticity outside the
star-forming region. This implies that the Bose star formed in this case might
be having an intrinsic angular momentum transferred from the initial cloud.
Thus we believe that our results about the angular momentum transfer are in
agreement with the analytical work carried out in Ref. [7] in which the
stability of a rotating Bose star with different interaction has been
analysed.

##  5 Acknowledgments

We are grateful to Bodo Schwabe, Mateja Gosenca, and Dilip Angom for helpful
discussions. And for their help with understanding and modifying the code
AxioNyx. We also extend our appreciation to the anonymous referee for their
constructive comments that significantly improved the quality of this research
paper. All the computations were accomplished on the Vikram-100 and
Vikram-1000 HPC cluster at the Physical Research Laboratory, Ahmedabad.

## References

  *   * Abbott et al [2020] Abbott R, et al (2020) GW190521: A Binary Black Hole Merger with a Total Mass of $150M_{\odot}$. Phys Rev Lett 125(10):101102. 10.1103/PhysRevLett.125.101102, arXiv:2009.01075 [gr-qc]
  * Buckley et al [2021] Buckley JH, Dev PSB, Ferrer F, et al (2021) Fast radio bursts from axion stars moving through pulsar magnetospheres. Phys Rev D 103:043015. 10.1103/PhysRevD.103.043015, URL https://link.aps.org/doi/10.1103/PhysRevD.103.043015
  * Bustillo et al [2021] Bustillo JC, Sanchis-Gual N, Torres-Forné A, et al (2021) Gw190521 as a merger of proca stars: A potential new vector boson of $8.7\times{}{10}^{-13}\text{ }\mathrm{eV}$. Phys Rev Lett 126:081101. 10.1103/PhysRevLett.126.081101, URL https://link.aps.org/doi/10.1103/PhysRevLett.126.081101
  * Chavanis [2011] Chavanis PH (2011) Mass-radius relation of newtonian self-gravitating bose-einstein condensates with short-range interactions. i. analytical results. Phys Rev D 84:043531. 10.1103/PhysRevD.84.043531, URL https://link.aps.org/doi/10.1103/PhysRevD.84.043531
  * Chen et al [2021] Chen J, Du X, Lentz EW, et al (2021) New insights into the formation and growth of boson stars in dark matter halos. Phys Rev D 104:083022. 10.1103/PhysRevD.104.083022, URL https://link.aps.org/doi/10.1103/PhysRevD.104.083022
  * Chen et al [2022] Chen J, Du X, Lentz EW, et al (2022) Relaxation times for bose-einstein condensation by self-interaction and gravity. Phys Rev D 106:023009. 10.1103/PhysRevD.106.023009, URL https://link.aps.org/doi/10.1103/PhysRevD.106.023009
  * Dmitriev et al [2021] Dmitriev AS, Levkov DG, Panin AG, et al (2021) Instability of rotating bose stars. Phys Rev D 104:023504. 10.1103/PhysRevD.104.023504, URL https://link.aps.org/doi/10.1103/PhysRevD.104.023504
  * Eby et al [2016] Eby J, Kouvaris C, Nielsen NG, et al (2016) Boson Stars from Self-Interacting Dark Matter. JHEP 02:028. 10.1007/JHEP02(2016)028, arXiv:1511.04474 [hep-ph]
  * Hogan and Rees [1988] Hogan CJ, Rees MJ (1988) AXION MINICLUSTERS. Phys Lett B 205:228–230. 10.1016/0370-2693(88)91655-3
  * Hu et al [2000] Hu W, Barkana R, Gruzinov A (2000) Fuzzy cold dark matter: The wave properties of ultralight particles. Phys Rev Lett 85:1158–1161. 10.1103/PhysRevLett.85.1158, URL https://link.aps.org/doi/10.1103/PhysRevLett.85.1158
  * Kaup [1968] Kaup DJ (1968) Klein-gordon geon. Phys Rev 172:1331–1342. 10.1103/PhysRev.172.1331, URL https://link.aps.org/doi/10.1103/PhysRev.172.1331
  * Kolb and Tkachev [1993] Kolb EW, Tkachev II (1993) Axion miniclusters and bose stars. Phys Rev Lett 71:3051–3054. 10.1103/PhysRevLett.71.3051, URL https://link.aps.org/doi/10.1103/PhysRevLett.71.3051
  * Larson [1984] Larson RB (1984) Gravitational torques and star formation. mnras 206:197–207. 10.1093/mnras/206.1.197
  * Levkov et al [2017] Levkov DG, Panin AG, Tkachev II (2017) Relativistic axions from collapsing bose stars. Phys Rev Lett 118:011301. 10.1103/PhysRevLett.118.011301, URL https://link.aps.org/doi/10.1103/PhysRevLett.118.011301
  * Levkov et al [2018] Levkov DG, Panin AG, Tkachev II (2018) Gravitational bose-einstein condensation in the kinetic regime. Phys Rev Lett 121:151301. 10.1103/PhysRevLett.121.151301, URL https://link.aps.org/doi/10.1103/PhysRevLett.121.151301
  * Preskill et al [1983] Preskill J, Wise MB, Wilczek F (1983) Cosmology of the invisible axion. Physics Letters B 120(1):127–132. https://doi.org/10.1016/0370-2693(83)90637-8, URL https://www.sciencedirect.com/science/article/pii/0370269383906378
  * Pshirkov [2017] Pshirkov MS (2017) May axion clusters be sources of fast radio bursts? Int J Mod Phys D 26(07):1750068. 10.1142/S0218271817500687, arXiv:1609.09658 [astro-ph.HE]
  * Raby [2016] Raby S (2016) Axion star collisions with Neutron stars and Fast Radio Bursts. Phys Rev D 94(10):103004. 10.1103/PhysRevD.94.103004, arXiv:1609.01694 [hep-ph]
  * RUFFINI and BONAZZOLA [1969] RUFFINI R, BONAZZOLA S (1969) Systems of self-gravitating particles in general relativity and the concept of an equation of state. Phys Rev 187:1767–1783. 10.1103/PhysRev.187.1767, URL https://link.aps.org/doi/10.1103/PhysRev.187.1767
  * Sanchis-Gual et al [2019] Sanchis-Gual N, Di Giovanni F, Zilhão M, et al (2019) Nonlinear dynamics of spinning bosonic stars: Formation and stability. Phys Rev Lett 123:221101. 10.1103/PhysRevLett.123.221101, URL https://link.aps.org/doi/10.1103/PhysRevLett.123.221101
  * Savitzky and Golay [1964] Savitzky A, Golay MJE (1964) Smoothing and differentiation of data by simplified least squares procedures. Analytical Chemistry 36:1627–1639. 10.1021/ac60214a047
  * Schiappacasse and Hertzberg [2018] Schiappacasse ED, Hertzberg MP (2018) Analysis of Dark Matter Axion Clumps with Spherical Symmetry. JCAP 01:037. 10.1088/1475-7516/2018/01/037, [Erratum: JCAP 03, E01 (2018)], arXiv:1710.04729 [hep-ph]
  * Schwabe et al [2020] Schwabe B, Gosenca M, Behrens C, et al (2020) Simulating mixed fuzzy and cold dark matter. Phys Rev D 102:083518. 10.1103/PhysRevD.102.083518, URL https://link.aps.org/doi/10.1103/PhysRevD.102.083518
  * Seigar [2005] Seigar MS (2005) The connection between shear and star formation in spiral galaxies. Monthly Notices of the Royal Astronomical Society: Letters 361(1):L20–L24. 10.1111/j.1745-3933.2005.00056.x, URL https://doi.org/10.1111/j.1745-3933.2005.00056.x, https://academic.oup.com/mnrasl/article-pdf/361/1/L20/2995591/361-1-L20.pdf
  * Siemonsen and East [2021] Siemonsen N, East WE (2021) Stability of rotating scalar boson stars with nonlinear interactions. Phys Rev D 103:044022. 10.1103/PhysRevD.103.044022, URL https://link.aps.org/doi/10.1103/PhysRevD.103.044022
  * Turk et al [2010] Turk MJ, Smith BD, Oishi JS, et al (2010) yt: A multi-code analysis toolkit for astrophysical simulation data. The Astrophysical Journal Supplement Series 192(1):9. 10.1088/0067-0049/192/1/9, URL https://dx.doi.org/10.1088/0067-0049/192/1/9
  * Turner [1986] Turner MS (1986) Cosmic and local mass density of “invisible” axions. Phys Rev D 33:889–896. 10.1103/PhysRevD.33.889, URL https://link.aps.org/doi/10.1103/PhysRevD.33.889
  * Veltmaat et al [2018] Veltmaat J, Niemeyer JC, Schwabe B (2018) Formation and structure of ultralight bosonic dark matter halos. Phys Rev D 98:043509. 10.1103/PhysRevD.98.043509, URL https://link.aps.org/doi/10.1103/PhysRevD.98.043509
  * Visinelli [2021] Visinelli L (2021) Boson stars and oscillatons: A review. Int J Mod Phys D 30(15):2130006. 10.1142/S0218271821300068, arXiv:2109.05481 [gr-qc]
  * Weinberg [1972] Weinberg S (1972) Gravitation and Cosmology: Principles and Applications of the General Theory of Relativity. John Wiley and Sons, New York
  * Wheeler [1955] Wheeler JA (1955) Geons. Phys Rev 97:511–536. 10.1103/PhysRev.97.511, URL https://link.aps.org/doi/10.1103/PhysRev.97.511