Relativistic collapse of axion stars

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.


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 [1], but the concept has since been extended to a class of axion-like particles (ALP's) with ultra-light masses [2]. In ALP scenarios, the dark matter forms gravitationally bound objects which may form into galaxy cores [3], or for larger masses into axion mini-clusters [4][5][6]. 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 [3,7].
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 [1,7] ALP ∼ 0.
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 [2].
We follow the recent trend of referring to stable axion structures as axion stars (though the term Bose star is also frequently * Corresponding author. E-mail addresses: florent.c.michel@durham.ac.uk (F. Michel), ian.moss@newcastle.ac.uk (I.G. Moss). used in this context). So far three distinct scenarios of gravitational collapse for APL's have been identified [8,9]: 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 [10]. 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. [9] 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 [11]. 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 [12], which otherwise is only possible with less efficient mesh refinement techniques [13].
Throughout this work, we use units in which the reduced

Model and field equations
We take the generic axion potential V , which is typical of the potentials which represent axion dark matter [14][15][16]: The parameters m and f are related by (1) 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 where g μν is the metric. The focus of this paper is on spherically-symmetric collapse. Following [11,12,17], we use a very efficient integration scheme obtained by introducing the retarded time coordinate u and radial coordinate r, with metric ds 2 = −g(u, r)ḡ(u, r)du 2 − 2g(u, r)dudr + r 2 d 2 .
(4) As usual, d 2 is the metric on the unit sphere, and we suppose that g, ḡ are two smooth functions. Without loss of generality, up to a redefinition of u, we can impose boundary conditions at the origin, ḡ(u, 0) = 1. Imposing that there is no conical singularity at r = 0 then implies that g(u, 0) = 1 [11].
We follow the conventions of [11,12] 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: The first of these equations is a form of the Klein-Gordon equation which can be integrated using the method of characteristics. This is the only true evolution equation in the system, the other three equations are geometrical constraints.
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. (5) satisfy the characteristic equation for (5), Changing to null coordinates, so that h(u, r) becomes h(u, v), gives the evolution along the characteristic surfaces of constant v, In order to solve these field equations, we have adapted the numerical procedure from Refs. [11,18,12]. Starting from given initial data for h and r at u = 0, we first compute h(0, v), g(0, v), and ḡ(0, v) using (6)- (8). We evolve r and h in the u direction using Equations (9) and (10), discarding points for which r becomes negative. At each step, the constraints (6)-(8) 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, (10) and (9), 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 [19,20], since this is appropriate for null coordinate systems. A Schwarzschild black hole metric with mass M, for example, has g = 1, ḡ = 1 − 2G M/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 (10). In practice, the integration is stopped when ḡ/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 (11) gives a local quantity M B (u, v) which evolves according to When 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 (12) 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 [6,5]. The phase diagram for relativistic collapse in terms of mass and radius was obtained semi-analytically in Ref. [21]. The fully relativistic collapse of a massive scalar field was studied in some detail using null coordinates in Ref. [17] and using a 3+1 approach in Ref. [22].
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 [23]. The scalar field is oscillatory in time, and when projected on to the light-cone in flat space, 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 (6)- (8), 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 2 p . Fig. 1 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.  [22].
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 scalar collapse under gravity with a quartic self-interaction potential was first discussed many years ago [24]. More recent work on this model using the non-relativistic limit can be found in Refs. [23,25]. In actual fact, axion radiation turns out to be highly relativistic, as pointed out in Ref. [8]. This latter work considered relativistic fields with Newtonian gravity, as did the semi-analytic discussions of axion stars in Refs. [26,27]. A fully relativistic treatment of axion col-lapse with general relativity was given in Ref. [9]. The phase plane of mass and axion scale was also discussed in [9], 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 [23,25], R ≈ 2(24π 3 ) 1/2 M 2 p /Mm 2 . Fig. 1 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 Fig. 2, 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 Fig. 3. There is a tricritical point in agreement with Ref. [9], but the mass and axion scale at the tricritical point are substantially different from the earlier 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 [28,25], M ≈ 50.77 f M p /m, shown on Fig. 3.
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 Fig. 4. 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. [22], where the effect was ascribed to gravitational cooling, using a term introduced for scalar field emission by collapsing Boson stars in Ref. [6]. We have checked that the results are not due to numerical noise.

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. 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 [13] 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 ax-ion 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 [25], 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. We would expect properties like phase boundaries to retain their qualitative nature close to spherical symmetry, but in reality this can only be addressed with non-spherically symmetric codes. The strength of out method is that it provides an important check on the results obtained from more sophisticated 3+1 integration schemes. We hope the clearer view of the simplest-case scenario considered here will help guiding future studies in this direction.