Relativistic collapse of axion stars
Florent Michel Ian G. Moss



 

florent.c.michel@durham.ac.uk
Centre for Particle Theory, Durham University,
South Road, Durham, DH1 3LE, UK

ian.moss@newcastle.ac.uk
School of Mathematics, Statistics and Physics, Newcastle University, 
Newcastle Upon Tyne, NE1 7RU, UK

 
We study the gravitational collapse of an axion field in null coordinates, assuming spherical symmetry. 
Compared with previous studies, we use a simpler numerical scheme which can run, for relevant parameters, 
in a few minutes or less on a desktop computer. We use it to accurately determine the domains of parameter 
space in which the axion field forms a black hole, an axion star or a relativistic Bosenova.


Relativistic collapse of axion stars
    Ian G. Moss
    December 30, 2023
====================================





§ INTRODUCTION


Amongst the possible dark matter candidates, a coherent scalar field
with very low mass is an enticing possibility. The idea originated with the 
QCD axion <cit.>, but the concept has since been extended to a class of
axion-like particles (ALP's) with ultra-light masses <cit.>. In ALP scenarios, 
the dark matter forms gravitationally bound objects which may form into
galaxy cores <cit.>, or for larger masses into axion mini-clusters 
<cit.>. 
These objects are often stable only for a particular mass range, leaving the possibility of 
detectable cosmological signatures from the axion bound structures or from 
the remnants of their collapse.

ALP's are characterised by their mass m and decay constant (or symmetry
breaking scale) f.
Coherent ALP dark matter scenarios envision the dark matter energy density 
in the form of large-scale coherent axion oscillations of frequency ∼ m,
with density parameter <cit.>

    Ω_ ALP∼ 0.1(f 10^17  GeV)^2
    (m 10^-22  eV)^1/2,

although this is rather dependent on initial conditions. 
Spatial gradients
in the oscillating axion field induce “quantum” pressure forces
which are capable of supporting structures on the Kpc scale
for axion masses  around m∼ 10^-22  eV, or galaxy
Halo scales for m∼ 10^-24  eV <cit.>.


We follow the recent trend of referring to stable axion structures as axion stars
(though the term Bose star is also frequently used in this context). So far
three distinct scenarios of gravitational collapse for APL's have been identified
<cit.>: they can settle down quietly to an axion star; 
they can radiate away energy in bursts of relativistic axions or they can 
collapse to a black hole.
The second outcome is a relativistic analogue of the Bosenova phenomena 
in cold-atom physics <cit.>. Like the cold atoms in a Bosenova, the axions have an 
attractive self-interaction force which can overcome the quantum pressure. We 
will use the term Bosenova in this paper to refer to the axion collapse and radiation
phenomenon.

The fate of an axion clump can be represented on phase
diagrams labelled by parameters describing the axion properties and the initial conditions. Recently, Helfer et al. 
<cit.> have  produced a phase diagram for spherically symmetric 
collapse with axion decay constant f and the initial mass of the 
axion clump, and they have speculated that there 
is a tricritical point joining phase boundaries between the three outcomes. 
The aim of this paper is to provide convincing numerical evidence for the tricritical 
point using a particularly amenable form of the field equations, and to determine the parameter
values accurately at the phase boundaries.

We use the null-coordinate integration schemes introduced into spherically symmetric
gravitational collapse by Goldwirth and Piran <cit.>. 
The null techniques are particularly efficient because the coordinate grid flows
inwards with the collapsing matter. For example, the null methods can reproduce
the universal scaling phenomena in massless scalar collapse  <cit.>,
which otherwise is only possible with less efficient mesh refinement techniques 
<cit.>.


Throughout this work, we use units in which the reduced Planck constant 
ħ and velocity of light c are equal to unity. The reduced Planck mass 
M_p=(8 π G)^-1/2, where G is Newton's constant.



§ MODEL AND FIELD EQUATIONS


We take the generic axion potential V, which is typical of the potentials which represent axion dark 
matter <cit.>:

    V(ϕ) = m^2  f^2  1 - cosϕ/f.

The parameters m and f are related by (<ref>)
 if the cosmological dark matter density is in the form of coherent axion oscillations, but we will generally
take m and f as free parameters. The Lagrangian density  of the axion field is

    ℒ_ϕ = - g^μν/2 (_μϕ)  (_νϕ) - V ϕ,

where g_μν is the metric. 

The focus of this paper is on spherically-symmetric collapse. Following 
<cit.>, 
we use a very efficient integration scheme obtained by introducing the retarded time coordinate u and 
radial coordinate r, with metric

    s^2 = - g(u,r) g̅(u,r)  u^2 - 2  g(u,r)  u  r + r^2 Ω^2.

As usual, Ω^2 is the metric on the unit sphere, and we suppose that
g, g̅ are two smooth functions. 
Without loss of generality, up to a redefinition of u, we can impose boundary
conditions at the origin, g̅(u,0) = 1. Imposing that there is no conical singularity 
at r=0 then implies that g(u,0) = 1 <cit.>.

We follow the conventions of <cit.> and introduce the notation 
h̅ for the scalar field ϕ. Radial derivatives of h̅ are used to define an auxiliary field
h. One can show that the Einstein equations are fully equivalent to a system of first order equations:

    _u h - g̅/2_r h = 
    	h - h̅/2  r[  1 - 8 π G  r^2  V h̅ g - g̅]
    	- g/2 r  V' h̅, 
       _r ln (g) = 4 π G/r (h - h̅)^2, 
       _r  r g̅
    	=  1 - 8 π G  r^2  V h̅ g,
       _r  r h̅=h ,

The first of these equations is a form of the Klein-Gordon equation which can be integrated using 
the method of characteristics.
Starting from the initial data surface u=0, we label the ingoing radial null geodesics
by a coordinate v. The ingoing null geodesics for the metric Eq. (<ref>) satisfy
the characteristic equation for (<ref>),

    .∂_u r|_v=-g̅ 2.

Changing to null coordinates, so that h(u,r) becomes h(u,v), gives the evolution
along the characteristic surfaces of constant v,

    _u h = 
    	h - h̅/2  r[  1 - 8 π G  r^2  V h̅ g - g̅]
    	- g/2 r  V' h̅,

In order to solve these field equations, we have adapted the numerical procedure 
from Refs. <cit.>.   
Starting from given initial data for h̅ and r at u = 0, we first compute h(0,v), g(0,v), and g̅(0,v) 
using (<ref>-<ref>). We evolve r and h in the u direction using Equations (<ref>) and 
(<ref>), discarding points for which r becomes negative. At each step, the constraints 
(<ref>-<ref>) are solved by integrating in the v direction.

Evolution methods based on the 3+1 space and time coordinates solve their constraints
on the initial time hypersurface, usually as a boundary value problem, and can be subject to 
constraint violation at later times.
This is not an issue with the null coordinate formalism. 
The method only requires us to solve the ordinary differential equations, 
 (<ref>)  and (<ref>), with integrations over v at each time step. 
As a result, the method is remarkably accurate, fast and reliable.

When a black hole forms, it is possible to follow the evolution up to the null surface u=u_T 
which contains a marginally trapped surface at r = r_T. We define the final black hole mass M_H as the  
Bondi mass <cit.>,

    M_ H = lim_u→ u_Tlim_v →∞r/2  G 1 - g̅/g,

since this is appropriate for null coordinate systems. A Schwarzschild black hole metric with
mass M, for example, has g = 1, g̅=1-2GM/r and M_H=M.
At the marginally trapped surface g→∞, and the computational grid has to be compressed to 
counter the growth in the right hand side of (<ref>). In practice, the integration is stopped 
when g̅ / g reaches a predetermined value. The Bondi mass is calculated at the final
value of u and with v at the extreme edge of the coordinate grid.

Removing the limits from (<ref>) gives a local quantity M_B(u,v) which evolves according to 

    ∂_u M_ B = - 2 π r^2 2/g∂_uh̅^2 + g̅ V (h̅ ) .

When g, g̅, and V are positive, then ∂_u M_ B≤ 0.
We will use -∂_u M_ B as a measure of the energy flux from the collapsing star.
Any increase of M along in the ingoing null direction indicates (at least if V remains positive) an 
artefact from the numerical integration, and the corresponding runs are discarded.    We can also use (<ref>) to 
put bounds on the error in the black hole mass from truncating the integration before the trapped surface at u=u_T.
This gives better control of the black hole mass than we would have using the mass at the trapped
surface, r_T/2G, which was used in previous work.




§ NUMERICAL RESULTS


We preface the full analysis with some results on the collapse of a massive, real, scalar field 
without self-interaction. Depending on initial conditions, the system can collapse to a black hole 
or a stable oscillaton, i.e. an oscillating field configuration that maintains its radial profile 
<cit.>. The phase diagram for relativistic collapse
in terms of mass and radius was obtained semi-analytically 
in Ref. <cit.>. 
The fully relativistic collapse of a massive scalar field was studied in 
some detail using null coordinates in Ref. <cit.> and using a 3+1 approach in
Ref. <cit.>.  

We use the null coordinate approach to plot the phase diagram in terms of the axion mass m,
the initial radius R and Bondi mass M_ B.
The choice of initial density profile is somewhat arbitrary, but we choose to work with a Gaussian 
profile which has been used previously for Bose stars <cit.>. The scalar field
is oscillatory in time, and when projected on to the light-cone in flat space,

    h̅_i(r) = √(2  M/π^3/2 m^2  R^3)^-r^2 / (2  R^2)cos (m  r).

The pre-factor has been chosen so that the mass of the star is M in the non-relativistic limit Rm≫1.
The relationship between the radius and the ingoing null coordinate on the initial surface can be 
specified freely, but the uniform choice r=2v will be used for simplicity. Initial conditions on the
remaining fields are determined by the constraints (<ref>-<ref>), which ensure that we 
have a consistent set of initial conditions for the fully relativistic collapse. 


The metric and scalar fields are evolved using the method described above.
After dimensional rescaling, the solutions only depend on the initial parameters in the
combinations mR and mM_B/M_p^2. 
Figure <ref> shows the phase diagram for collapse in terms of these rescaled
parameters. When the mass and radius are appropriate for a stable axion star, the scalar 
field profile settles down to the axion star
field profile within a few oscillations. 
The phase transition boundary is traced out by dots, obtained by a bisection search technique. 
The condition for black hole collapse becomes nearly independent of radius when m  R ≳ 2, with critical 
Bondi mass M_B≈ 15.22 M_p^2m^-1 and mass parameter M≈ 15.55 M_p^2m^-1. 
The Bondi mass is in very good agreement with the results of the 3+1 approach, where the critical ADM mass is 
M_ ADM≈ 15.2 M_p^2m^-1 (after accounting for the switch from Planck mass to 
reduced Planck mass) <cit.>.

Inclusion of self-interaction leads to a third possible outcome of gravitational collapse,
a Bosenova, where the collapsing field loses mass in 
pulses of axion radiation. The possibility of collapse due to self-interaction was first noticed in
the non-relativistic limit <cit.>, but the radiation 
is highly relativistic, as pointed out in Ref. <cit.>. A fully relativistic
treatment of axion collapse with general relativity was given in Ref.  <cit.>. The phase plane of mass 
and axion scale was also discussed in <cit.>, 
where evidence was given for the existence of a tri-critical point between the three outcomes. 
We make use of our rapid integration scheme to give a clearer picture of the phase plane.

The Gaussian profile is used for the initial data as before, but with fixed initial radius. 
It is desirable to have stable axion stars settle down as quickly as possible to their
final form in order to keep down the integration time. In order to achieve this, we choose the initial radius
using the radius for non-relativistic Bose stars <cit.>,
R≈ 2(24π^3)^1/2M_p^2/Mm^2.
Figure <ref> shows that this radius lies in the region of the phase diagram
where the dependence on radius is weak.


The Bosenova is characterised by collapses of the stellar core followed by bursts
of axion radiation. The collapse and burst pattern is repeated until a significant portion
of the initial mass has been radiated away. An example is illustrated in figure <ref>, which shows
the central density and the radiation escaping at the edge of the integration volume 
as functions of retarded time. The radiation escapes a short retarded time after each collapse, 
indicating that the radiation is highly relativistic. Burst can be highly irregular, both in their timing
and amplitude. When a black hole forms instead of a Bosenova, this tends to happen 
at the same retarded time as the first spike in the central density.



The phase diagram is shown in figure <ref>. There is a tricritical point
in agreement with Ref. <cit.>, but the mass and axion scale at the tricritical point are
substantially different from the eralier results. The difference is too large to be explained by the
difference between using initial conditions on a null surface instead of
a timelike surface. Our results are broadly consistent with the non-relativistic limit though, where
there is a critical mass for the collapse of a Bose star <cit.>,
M≈ 50.77  f M_p /m, shown on figure <ref>. 


The phase boundary between axion stars and black holes is sharply
defined, and the mass of the black hole is discontinuous at the phase boundary.
However, the phase boundary between black holes and Bosenovas becomes diffuse
at small values of the axion scale parameter f, in the sense that there is a range of 
masses near the phase boundary where the outcome of gravitational collapse can go either way, 
as shown in figure <ref>. Some of the axion clumps with initial conditions near 
the phase boundary emit enough axion radiation to avoid forming a black hole.
This seems to happen erratically. A similar phenomenon was observed for 
the collapse of non-self interacting scalar fields in Ref.  <cit.>, 
where the effect was ascribed to gravitational cooling, using a term introduced for scalar field
emission by collapsing Boson stars in Ref. <cit.>. 








§ DISCUSSION


Gravitational collapse with nothing more than gravity and a scalar field is
a remarkably rich subject. It seems sensible to build up an understanding of it
in small steps, the simplest being spherically symmetric collapse.
We have considered three possible scenarios for axion collapse: axion stars, black holes and 
Bosenovas. The numerical results clearly point to a critical point with Bondi mass 
M_B≈ 10.6 M_p^2m^-1 and axion decay constant f≈ 0.25 M_p
when the initial conditions are presented on a phase diagram.
Mostly, the distinction between the different phases is clear, but in some parts of the 
phase diagram, there is no clean line between initial conditions which collapse to a black hole 
and those which remain non-singular. There may also exist special
final states like the self-similar solutions for massless scalar collapse <cit.>
which we have not considered.

The fate of a Bosenova is to eject mass until it eventually settles down into a
stable axion star. In terms of the eventual outcome, the stars and Bosenova's
are similar. However, the difference has a large physical significance for
dark energy scenarios, since the axion radiation from many bursts and many
sources would combine into a background of incoherent ALP particles. 

Just how dependent the results are on the initial density profile and the use of null initial
data surface may be determined through further work. The Gaussian initial density profile is
known to work to within a few percent for results on non-relativistic Bose stars <cit.>,
but we have done some runs with radically different density profiles and find O(1) changes in the
masses at the phase boundaries.

The null coordinate approach is fast, accurate and reliable, but it is limited to spherically symmetric
collapse. It provides an important check on the results obtained from
more sophisticated 3+1 integration schemes, but ultimately the 3+1 schemes are necessary for
handling general collapse situations. 
We hope the clearer view of the simplest-case scenario considered here will help guiding future studies in this direction. 

We would like to acknowledge useful discussions with Ruth Gregory and Gerasimos Rigopoulos.
We are grateful for the hospitality of the Perimeter Institute, Waterloo, Canada where
this paper was written. IGM and FM are supported by the Leverhulme grant RPG-2016-233. IGM also 
acknowledges some support from the Science and Facilities Council of the 
United Kingdom grant number ST/P000371/1.






