Axion star collisions with black holes and neutron stars in full 3D numerical relativity

Axion star collisions with black holes and neutron stars in full 3D numerical relativity

Katy Clough    Tim Dietrich    Jens C. Niemeyer Institut für Astrophysik, Georg-August Universität, Friedrich-Hund-Platz 1, D-37077 Göttingen, Germany Nikhef, Science Park, 1098 XG Amsterdam, The Netherlands

Axions are a potential dark matter candidate, which may condense and form self gravitating compact objects, called axion stars (ASs). In this work, we study for the first time head-on collisions of relativistic ASs with black holes (BHs) and neutron stars (NSs). In the case of BH-AS mergers we find that, in general, the largest scalar clouds are produced by mergers of low compactness ASs and spinning BHs. Although in most of the cases which we study the majority of the mass is absorbed by the BH within a short time after the merger, in favourable cases the remaining cloud surrounding the final BH remnant can be as large as of the initial axion star mass, with a bosonic cloud mass of and peak energy density comparable to that obtained in a superradiant build up. This provides a dynamical mechanism for the formation of long lived scalar hair, which could lead to observable signals in cases where the axion interacts with baryonic matter around the BH, or where it forms the seed of a future superradiant build up in highly spinning cases. Considering NS-AS collisions we find two possible final states (i) a BH surrounded by a (small) scalar cloud, or (ii) a stable NS enveloped in an axion cloud of roughly the same mass as the initial AS. Whilst for low mass ASs the NS is only mildly perturbed by the collision, a larger mass AS gives rise to a massive ejection of baryonic mass from the system, purely due to gravitational effects. Therefore, even in the absence of a direct axion coupling to baryonic matter, NS-AS collisions could give rise to electromagnetic observables in addition to their gravitational wave signatures.

neutron stars, axion stars, black holes, numerical relativity

I Introduction

In the wake of multiple LIGO detections Abbott et al. (2016a, b, 2017a, 2017b, 2017c), including GW170817 Abbott et al. (2017d), the first combined detection of GWs and electromagnetic signals from the same astrophysical source, there has been renewed interest in the simulation of exotic compact objects (ECOs) which could mimick BH or NS observations, or provide altogether new, and as yet undetected, observational signatures Giudice et al. (2016).

One of the simplest potential ECOs is a boson star (BS), which is a stable solitonic solution to the coupled Einstein-Klein-Gordon equations for a complex scalar field with gravity. The idea of a self gravitating field configuration dates back to proposals by Wheeler for “geons” Wheeler (1955), but was first shown to work for complex scalar fields in Kaup (1968); Ruffini and Bonazzola (1969). The ideas were extended to real massive scalar fields in Seidel and Suen (1991), with the (quasi) stable objects later dubbed oscillotons. A non trivial self interaction potential, motivated by, for example, low energy effective theories from string theory or other models, modifies the stability and profile of the solutions, giving rise to new classes of self interacting BSs Eby et al. (2016); Colpi et al. (1986); Krippendorf et al. (2018).

The standard model of particle physics does not contain a bosonic particle that would allow the formation of BSs. Dark matter, on the other hand, may well be composed of bosons, with axion-like particles (ALPs) being a well-motivated class of candidates (see Marsh (2016) for a thorough review). ALPs are very light, weakly coupled particles that are produced with practically vanishing momenta and extremely high occupation numbers. They are usually treated as classical, real scalar fields, subject to a cosine potential parametrised by the axion decay constant and the axion mass . The leading order self-interaction for axions is thus attractive, but whilst here we simulate the full cosine potential, we use a large value of such that the the axion is effectively a massive boson and self interactions are negligible. In future work we hope to expand the study to quantify the effect of increasing self interactions.

In this work, we focus on ASs with masses comparable to the objects with which they collide since such setups will give rise to strong GW signals. Consequently, the AS masses in the NSAS configurations are , whilst in the BHAS case the physical mass of the BH sets the scale of the simulation such that the interpretation as solar mass BHs corresponds to solar mass ASs, and for supermassive BHs similarly the ASs have masses up to . With the AS mass fixed, we choose values for the axion mass such that the star is comparatively compact - ie, such that the de Broglie wavelength of the axion is comparable to the radius of the BH or NS. In the solar mass case this corresponds to eV (the lower end of the QCD axion scale), and for supermassive BHs, as low as eV (corresponding to ultra-light ALPs).

The original QCD axion emerges as a consequence of the Peccei-Quinn (PQ) symmetry breaking mechanism to solve the strong CP problem Peccei and Quinn (1977); Weinberg (1978). Its mass is constrained to lie between eV and eV, with decay constant . If dark matter consists of QCD axions whose PQ symmetry was broken after inflation, high-amplitude density fluctuations on the scale of the cosmological horizon at the time of the QCD phase transition are predicted to exist. They collapse during the radiation era and form so-called axion miniclusters with typical masses of Hogan and Rees (1988); Kolb and Tkachev (1993). A fraction of these masses may have formed ASs either directly during the collapse Schive et al. (2014); Veltmaat et al. (2018) or by scalar wave condensation Levkov et al. (2018). The subsequent evolution of the AS mass function as a result of minicluster mergers or ongoing condensation is largely unknown; although their typical masses at formation are well below those we consider in this work, the existence of a high-mass tail of the distribution cannot be ruled out at present. A population of relativistic ASs could also have been produced by non-standard primordial perturbations with enhanced small-scale power Widdicombe et al. (2018).

Ultra-light ALPs with masses in the range of eV, motivated from string theory compactifications Svrcek and Witten (2006); Arvanitaki et al. (2010), are candidates for “fuzzy dark matter” (FDM) with interesting new phenomenology on scales of their de Broglie wavelength in dark matter halos on scales of kpc Hu et al. (2000); Hui et al. (2017). In particular, cosmological simulations using the nonrelativistic Schrödinger-Poisson equations to describe FDM show the formation of ASs in the form of solitonic halo cores Schive et al. (2014); Veltmaat et al. (2018) with masses of the order of and above. The evolution of their mass function as a result of halo mergers was studied in Du et al. (2017) and shown to approach the core-halo mass relation found by Schive et al. (2014).

Thus, even if axions make up all the dark matter, relativistic ASs would constitute only a small () fraction of the total mass (also, note that observational constraints from microlensing on compact ASs would be similar to those on primordial BHs). Nevertheless, there could still be a sufficient number density for collisions with known objects such as BHs and NSs to occur, motivating an exploration of their observable signatures.

Simulations of ECO mergers more generally have to date focused mainly on complex scalar field BSs, e.g. Balakrishna (1999); Palenzuela et al. (2007); Choptuik and Pretorius (2010); Bezares et al. (2017); Lai (2004); Palenzuela et al. (2008); Cardoso et al. (2016); Dietrich et al. (2018a), but other classes such as real scalar field oscillotons Helfer et al. (2018), and just recently, massive, complex vector field Proca stars in Sanchis-Gual et al. (2018, 2017) have been explored.

There has, to the best of our knowledge, not yet been a study of the merger of ASs with BHs and NSs, even though these collisions are of interest since the ASs could act as a potential NS or BH mimickers in GW signals. In addition, since some types of axions are expected to couple weakly to baryonic matter, such mergers could result in distinct multi-messenger events, with, for example, phenomena such as Fast Radio Bursts (FRBs). In the case of a NS-AS collision the axions may interact with the neutron star matter either during or after the merger (see for example Hook and Huang (2018); Eby et al. (2017); Raby (2016); Iwazaki (2015); Barranco et al. (2013) and the discussion in Dietrich et al. (2018b)). Similarly, the case of a BH-AS merger may provide a mechanism by which one can dynamically form long lived (quasibound) scalar clouds, or “wigs”, around BHs Barranco et al. (2012, 2017); Sanchis-Gual et al. (2016); Herdeiro and Radu (2015); Hod (2012). Such clouds could then interact with any baryonic matter present in a potential accretion disc, or, in appropriate cases, provide the seeds for a superradiant instability (see e.g. Brito et al. (2015a, b)) to develop. Several other novel methods for detecting such scalar clouds have also been proposed, see e.g. Ferreira et al. (2017); Hannuksela et al. (2018); Degollado and Herdeiro (2014). Consequently, these events have great potential to further constrain the properties of the axion sector.

For our simulations, we assume that the axion field is only coupled to baryonic matter via the gravitational interaction, i.e., they interact due to their mutual impact on the metric, and no additional couplings are implemented. As discussed, we also focus on relatively high values for the axion decay constant of , meaning that self interactions are negligible. These choices are most relevant to the case of ALPs with a negligible coupling to standard model physics.

The paper is organised as follows. Sec. II summarises the numerical methods used for the work and presents code tests, as well as a preliminary investigation to determine the initial conditions for the simulations shown in the remainder of the paper. In Sec. III we consider the impact of varying the AS compactness on the remnant clouds left by BH-AS from mergers of ASs with spinning and non-spinning BHs.

In Sec. IV we discuss NS-AS collisions, and investigate the remnant of the head-on collision as the NS-AS mass ratio is varied. We summarize our results and future plans in Sec. V.

Throughout this paper we employ geometric units, with and a mass scale which in the case of NSs is set to , but in the case of BHs is a free mass scale by which the results may be scaled for varying BH masses. Consequently code units for lengths, times and masses are multiples of the mass scale 111Note that the scale which appears in the potential function is the quantity , with dimension such that a value of in code units corresponds to a particle mass of in the NS case (note that in geometric units), we refer to Dietrich et al. (2018a) for further discussion..

Ii Numerical methods

In this section we describe the set ups used for the respective simulations. We use the GRChombo code Clough et al. (2015) for the BH-AS evolutions, and the BAM code Brügmann et al. (2008); Thierfelder et al. (2011); Dietrich et al. (2018a) for the NS-AS cases, but undertake comparisons to ensure they give consistent results.

ii.1 BH-AS collisions using GRChombo

The GRChombo code Clough et al. (2015) ( is based on the method-of-lines with fourth order spatial finite differences and fourth order explicit Runge-Kutta timestepping. It is built on top of the Chombo Adams et al. () framework for solving partial differential equations with adaptive mesh refinement (AMR). It supports non-trivial “many-boxes-in-many-boxes” mesh hierarchies using block structured Berger Rigoutsos grid generation Berger and Rigoutsos (1991). The grid is made out of a hierarchy of cell-centered Cartesian grids consisting of a maximum of refinement levels labeled with resolution increasing by a factor of on successive levels according to . Given that the mesh is adaptive, the configuration of grids changes dynamically during the simulation. Around the BH all refinement levels are utilised, but for the scalar matter the number required varies between 4 and 9. The regridding is determined by setting the thresholds , , and cut offs , in the regridding condition:


such that smaller values of the thresholds force regridding for less varying/oscillatory data. In principle this allows us to maintain a consistent level of refinement on the scalar matter around the BH, even as it decreases in amplitude by several orders of magnitude during the merger; see Adams et al. (); Clough et al. (2015) for further details on the regridding process. The Berger-Oliger scheme is employed Berger and Oliger (1984) to coordinate time stepping on the grid hierarchy and we use a Courant-Friedrichs-Lewy factor of on each level. The evolution of the scalar fields is based on the Klein-Gordon Equation and uses the methods outlined in Ref. Helfer et al. (2017). The exact grid configurations for the three resolutions are given explicitly in Tab. 1.

R1 8 128 (1.0,0.20,0.06) 0.48 3e-6 1e-2 8.00 0.031250
R2 8 192 (0.6,0.16,0.04) 0.32 3e-6 1e-2 5.33 0.020833
R3 8 256 (0.5,0.10,0.03) 0.24 3e-6 1e-2 4.00 0.015625
Table 1: Grid configurations for the GRChombo simulations for spinning BHs. The columns refer to: total number of levels, number of points along each dimension on the coarsest level, the regridding threshold on (for cases ), the regridding threshold on , the regridding cutoff on , the regridding cutoff on , coarsest grid spacing, and finest grid spacing.

ii.2 NS-AS collisions using BAM

The BAM code Brügmann et al. (2008); Thierfelder et al. (2011); Dietrich et al. (2015a); Bernuzzi and Dietrich (2016); Dietrich et al. (2018a) is based on the method-of-lines, Cartesian grids and finite differencing. The grid is made out of a hierarchy of cell-centered nested Cartesian boxes consisting of refinement levels labeled with increasing resolution. Each level’s resolution increases by a factor of two leading to a resolution of . The inner levels employ points per direction and move following the technique of ‘moving-boxes’ where the star’s center are estimated by the minimum of the conformal factor . The outer levels remain fixed and employ grid points. Because of the symmetry of the problem, we employ bitant-symmetry to reduce the number of grid points and the computational costs by a factor of two. For the time integration a fourth order explicit Runge-Kutta type integrator is used with a Courant-Friedrichs-Lewy factor of . For the time stepping of the refinement level the Berger-Oliger scheme is employed Berger and Oliger (1984). Metric spatial as well as scalar field derivatives are approximated by fourth order finite differences. The general relativistic hydrodynamic equations are solved with standard high-resolution-shock-capturing schemes based on primitive reconstruction and the local Lax-Friedrich central scheme for the numerical fluxes; see Thierfelder et al. (2011); Bernuzzi et al. (2012) for more details. The evolution of the scalar fields based on the Klein-Gordon Equation uses techniques outlined in Ref. Dietrich et al. (2018a). The exact grid configurations are given explicitly in Tab. 2.

For the construction of the initial NS-AS configurations we use BAM’s multigrid solver. The solver employed the Conformal Thin Sandwich (CTS) formalism to obtain initial configurations in agreement with Einstein’s field equations. However, we note that as outlined in Dietrich et al. (2018a) we do not resolve the equations for general relativistic hydrodynamics or the Klein-Gordon equation which would be required to obtain fully consistent initial configurations.

R1 7 160 4 96 16.00 0.25
R2 7 240 4 144 10.67 0.167
R3 7 320 4 192 8.00 0.125
Table 2: Grid configurations for the BAM simulations. The columns refer to: total number of levels, number of points along each dimension, number of moving box levels using points per direction, coarsest grid spacing, and finest grid spacing.

ii.3 Code testing and comparison

Figure 1: Single axion star comparison between GRChombo (blue) and BAM (red). We employ three different axion star potentials varying .

For all presented configurations we performed numerical simulations with different resolutions to provide an error estimate. In addition, we have validated our results by comparing directly the outcome of the GRChombo and BAM codes. For this purpose we evolve a set of single star AS configurations as well as one example BH-AS merger. We also refer to Ref. Dietrich et al. (2018a) for a comparison of a BS-BS head-on merger between BAM and GRChombo.

Starting with isolated AS simulations, in Fig. 1 we present the time evolution of the central value of the scalar field for a number of axion star potentials of the form


Here is the axion decay constant and is an inverse length scale which is a function of the axion mass .

For the comparison we picked an initial configuration with . For an axion decay constant of we find sinusoidal oscillations between (top panel). Once we decrease the axion decay constant to the AS undergoes additional modulations (middle panel). Finally for the pressure support is not strong enough to avoid gravitational collapse and a BH forms (bottom panel). For all three cases is almost identical for BAM and GRChombo and we find excellent agreement.

Figure 2: Spatial profile of the conformal factor for BAM (blue) and GRChombo (red) for different instances of time for a BH-AS collision.

More challenging than the single star comparison is the evolution of a BH-AS merger simulation. For this purpose we perform a head-on collision with an initial separation of , an axion potential with , and an AS with an initial central field value of (as for the single AS simulations presented above). For comparison we present, although gauge dependent, the spatial profile of the conformal factor for different times (Fig. 2). We find overall good agreement between BAM and GRChombo although details in the simulation scheme and methods are different as, e.g., the usage of the Z4c for BAM or the CCZ4 evolution scheme for GRChombo.

ii.4 Dependence on initial separation

In order to determine for our simulations the optimal initial separation of the stars and BHs, we tested the impact of initial separation on the bosonic cloud mass which formed in the case of a zero spin BH-AS merger. The results are shown in Fig 3.

Figure 3: Evolution of the bosonic cloud mass surrounding the final BH remnant, for varying initial separations , for an AS with , and a BH of mass and . We present results for resolution R2 with solid lines and corresponding results for resolution R1 are shown dashed. Between and there is no significant difference in the cloud mass which forms.

We find that the mass of the remnant cloud gets larger for an increasing initial separation. However, the absolute difference between large distances () is small and below the uncertainty introduced by the numerical discretization, cf. dashed and solid lines. We thus choose an initial separation of 100 , where denotes the total mass of the system, i.e., the sum of the masses of the individual stars in isolation.

Iii BH-AS mergers

In this section we study the head-on mergers of BHs and ASs. The BH can take any physical mass , from a small primordial BH to supermassive BH, with the results scaled accordingly. We take the mass of the AS to be of order the BH mass, and set in geometric units with , such that the relevant axion mass varies according to eV. This sets the radius of the axion star to be roughly of order the Schwarzschild radius of the BH, which is the regime which is most favourable for GW production. However, we mainly focus on ASs with a relatively low mass and compactness, since we expect (and indeed find) that these are the best candidates to form large scalar clouds around the final BH remnant. Unlike in the case of superradiance, there is no physical requirement that , which essentially sets the de Broglie wavelength of the axion to the BH radius. In practise, it is numerically challenging to simulate ASs with vastly different sizes and masses to the BH, which is another reason that we focus in this work on such setups. These cases are of particular interest, since their remnants may go on to form a superradiant instability, but we emphasise that what we describe can be made more general, and does not require an exact mass correspondence or extremal spins. Note that, as shown in Barranco et al. (2012), the dynamical resonances around the BH should be longer lived for cases in which , in which case the half life of the scalar field can have a cosmological timescale. Given our choice of parameters, we should expect to see shorter lived solutions, but to confirm this we would need to evolve for longer periods beyond the merger, which is numerically challenging.

iii.1 Configurations

For the study of BH-AS mergers, we perform simulations for 6 different binary configurations using three AS compactnesses and two BH spins, cf. Table 3.

For the construction of the initial BH-AS configurations, we use the initial data for the ASs as in Ref. Helfer et al. (2017), with and set to an initial profile which is a stable solution for the free field oscilloton. Since we use a large value of (above the “triple point” as described in Ref. Helfer et al. (2017)), the stars are stable and far from collapsing to BHs, and the modulation of the base frequency due to the self interaction is small.

The AS’s ADM (Arnowitt-Deser-Misner) mass is computed by fitting the Schwarzschild solution to the numerical solution which is obtained via a shooting method, as in Helfer et al. (2017). The bosonic mass of the AS is determined by integrating the density of the field over the spatial volume of the initial slice.

We superpose this data with Bowen-York data Bowen and York (1980) for the spinning BHs, using the perturbative analytic solution for given in Gleiser et al. (1998), and then apply a relaxation procedure (see Clough et al. (2015)) to reduce the Hamiltonian constraint violation. Note that we first convert the AS data into conformally flat coordinates, so that it may be superimposed onto the conformally flat BH data with less deformation. The spin axis (the z-axis) is perpendicular to the merger direction (along the x-axis), and the initial distance is adjusted such that at it is approximately ; cf. Sec. II.4 for a study of the effect of different initial separations.

1.00 0.20 0.20 0.014 0.0 0.02 120
1.00 0.36 0.35 0.048 0.0 0.07 136
1.00 0.52 0.52 0.130 0.0 0.20 152
1.00 0.20 0.20 0.014 0.5 0.02 120
1.00 0.36 0.35 0.048 0.5 0.07 136
1.00 0.52 0.50 0.013 0.5 0.20 152
Table 3: Simulated BHAS setups. The columns refer to: the name of the configuration, the ADM mass of the BH in isolation , the AS’s ADM mass in isolation , the AS’s bosonic mass in isolation , the AS’s compactness in isolation (where is the effective radius of the oscillaton which encompasses 95% of its mass) the spin parameter for the BH , the central value of the axion scalar field , and the initial separation in units of .

iii.2 Qualitative merger dynamics

We find that the BH-AS merger provides a dynamical mechanism for the formation of scalar clouds around the BH.

We will see in the following section that we can classify the merger dynamics of NS-AS systems into three qualitatively different categories depending on the compactness (and mass) of the AS with which they collide. In the BH-AS cases there is no such distinction, but rather a gradual change in the formation of the bosonic clouds with AS compactness.

To illustrate this, we show in Fig. 5 the bosonic energy density, see Ref. Dietrich et al. (2018a)222Energy densities expressed in geometric units correspond to , and the conformal factor for different instances of time for the three setups (left row), (middle row), and (right row). For the final time, the result for the spinning case is also shown in purple for comparison. In Fig. 4 we present a volume rendering for the case, to visualise the full 3D configuration of the clouds formed.

Figure 4: A volume rendering of the bosonic energy density (blue) and the BH conformal factor (grey) for the case for the initial, merger and final configurations. The two clouds which form are aligned with the x-axis along which the merger occurs.

Case I []: The metric, cf. the evolution of the conformal factor, changes only slightly once the BH and the AS merge (most of the change shown here is a gauge effect). We see that before the AS and the BH centres collide, the bosonic matter already starts to be sucked into the BH. At the axion cloud settles to a fairly incoherent, slowly decaying configuration, with a maximum energy density of the order of and a spatial extent of . A similar process occurs in the spinning case, but, most likely due to the greater dispersal of the cloud before merger due to frame dragging of the bosonic matter, the final cloud is larger than in the non spinning case, with a peak energy density of .

Case II []: For the larger AS mass, one can see in Fig. 5 that the AS is more compact with a larger central density. During the merger process, it retains more of its shape, and, since it is more compact, a greater amount is sucked into the BH once the centres meet. At the axion cloud is more coherent than in the case. It also has a maximum energy density of the order of , but with a spatial extent of , which explains the approximate order of magnitude drop in the final cloud mass shown in Fig. 6. Again the spinning case results in a larger final cloud, with a peak energy density of .

Case III []: For the largest AS mass considered, the AS is barely perturbed until the final moments of the merger and the majority of the bosonic matter falls rapidly into the BH. The remaining cloud is negligible, with a peak energy density of ( for the spinning case), and a spatial extent of .

Figure 5: Evolution of the axion energy density in the non spinning (blue) and spinning (purple) cases, We also show the conformal factor as a gray dashed line, cf.  right axis. Different rows show different instances of time as labelled in each panel. We show the quantities at , a later time corresponding approximately to the start of the merger, and the final timestep . The different columns refer to the configurations: (left), (middle), (right), i.e., becoming more massive, and more compact, from left to right. Only the final state is shown for the spinning case.

Overall, we see that larger clouds are formed initially for ASs with a lower mass and correspondingly lower compactness, and for higher spins. This is quantified further in the following section.

iii.3 Bosonic clouds from merger

In the case of BH-AS mergers, we are particularly interested in regions which generate large scalar field clouds, as this may lead to observable effects post merger, as discussed in the introduction. The cases we consider span different compactness of the AS (which also corresponds to adjusting the mass ratio with the BH), and cases with and without spin ( and ). Our results are summarised in Fig. 5, Fig 6, and Fig. 7.

In Fig 6, we show the total mass of axions over the course of the merger, obtained from integrating


over the numerical domain at each timestep. Whilst we do not explicitly excise the mass within the event horizon, the particular choice of gauge condition causes the matter inside the BH to fall into the “puncture” and to leave the numerical domain. Consequently, the majority relates to that outside the horizon. This can be clearly seen in Fig. 5. We find, which agrees to intuition, that the largest scalar clouds are generated for low compactness ASs. This appears to be mainly because the bosonic matter in the stars is more diffuse and less gravitationally bound, and thus becomes dispersed about the BH before the actual merger occurs. In the higher mass cases, the AS is more compact and retains its shape during the merger. Thus when it overlaps the BH, the majority of the mass falls immediately into the horizon.

In Fig. 7, we show the effect of the BH spin on the merger outcome. We find that, in general, more massive scalar clouds form from low compactness ASs and spinning BHs. We observe that particular combinations appear to be more efficient at creating massive clouds, which might be caused by particular quasi normal mode (QNM) or quasi-stable bound states of the BH being excited.

Although over 98% of the mass is absorbed by the BH within 1000 after the merger for most of the studied configurations, for some favourable cases the remaining cloud can be as large as 30% of the initial AS mass, with a mass of order , and energy densities of . This is comparable to the values obtained in superradiant build up around BHs (see Brito et al. (2015b)).

As mentioned above, the results presented for BHs have a free mass parameter which means that the mass of the AS/BH system can be scaled by the mass of the axion. For the QCD axion eV, our results would correspond to a black hole of order several solar masses, whereas for so called “fuzzy dark matter” (FDM) with eV they would correspond to supermassive black holes.

Figure 6: Evolution of the bosonic cloud mass surrounding the BH over the course of the merger. We present results for resolution R3 with solid lines and corresponding results for resolution R2 are shown dashed. We see that less compact ASs result in larger masses of clouds post merger, despite their initially lower mass.
Figure 7: Evolution of the bosonic cloud mass surrounding the BH showing the effect of spin. We present results for resolution R3 with solid lines representing the non spinning case, , and dashed lines showing the results for the corresponding case with . We see that there is a clear increase in the mass of the remaining cloud with spin. Note that the spin axis is perpendicular to the merger direction.

iii.4 Gravitational wave signals

Figure 8: GW metric multipoles for the three example cases of Fig. 5, plotted on the same axis for comparison. The least massive case gives a barely distinguishable signal in comparison to the most compact ASs.

In Fig. 8 we show the dominant (2,2)-mode of the GW signal for the cases in Fig. 5. We find that for the Case I setup, e.g. , the emitted GW signals emit GW energy in the (2,2)-mode of . For increasing AS masses the total released energy can be increased significantly. For the setup we obtain amplitudes about an order of magnitude larger than for , and the emitted GW energy in the dominant mode is . For the largest case considered, , the amplitude is again an order of magnitude larger and the energy in the dominant mode is .

Iv NS-AS mergers

Following the study of BH-AS collisions, we will focus in this section on head-on mergers of NS-AS systems. In contrast to BHs, NSs have an intrinsic mass scale set by the Equation of State (EOS). This sets the axion mass to eV for the NSs which we consider, so that the radius of the ASs is comparable to that of the NS (ie, since we set in geometric units).

iv.1 Configurations

For the study of NS-AS mergers, we perform simulations for seven different binary configurations listed in detail in Table 4. All setups are evolved with resolutions R1, R2, R3 as presented in Tab. 2. To compute the initial configurations we keep the bosonic and baryonic energy density fixed and solve the CTS equations to obtain configurations consistent with general relativity. Initially, we pick a NS mass in isolation of , which due to changes in the conformal factor for the binary can slightly increase during the iterative procedure to solve the Einstein’s constraint equations; see Dietrich et al. (2018a) for details. The axion star configurations are calculated as in the BH-AS case and employ central values of the dominant zeroth mode of the scalar field of (Table  4). The initial distance is again adjusted such that at it is approximately . We employ for the simulation of the baryonic matter the SLy EOS Douchin and Haensel (2001); Read et al. (2009); Dietrich et al. (2015b) which is in agreement with most of the current EOS constraints directly inferred from GW170817 Abbott et al. (2017d, 2018); Coughlin et al. (2018); Bauswein et al. (2017); Radice et al. (2018); Annala et al. (2018); Most et al. (2018).

1.376 1.527 0.201 0.200 0.014 0.02 155.1
1.380 1.532 0.279 0.277 0.028 0.04 162.9
1.383 1.536 0.334 0.330 0.042 0.06 168.4
1.384 1.537 0.357 0.352 0.048 0.07 172.8
1.385 1.538 0.378 0.371 0.055 0.08 172.8
1.387 1.540 0.413 0.404 0.068 0.10 176.3
1.391 1.546 0.525 0.502 0.130 0.20 187.4
Table 4: Simulated NSAS setups. The columns refer to: the name of the configuration, the gravitational mass of the NS in isolation , the baryonic mass of the NS , the AS’s ADM mass in isolation , the AS’s bosonic mass in isolation , the AS’s compactness in isolation the central value of the axion scalar field , and the initial separation .

iv.2 Qualitative merger dynamics

Figure 9: Evolution of the bosonic energy density (blue) incorporating the determinant of the 3-metric to allow the computation of the bosonic cloud (Fig. 11) and the baryonic density (red). We also present the conformal factor as a grey dashed line, cf.  right axis. We show the conformal factor for at level and for later times at refinement level . The different columns refer to the configurations: (left), (middle), (right). Different columns show different instances of time as labeled in each panel.

We can classify the merger dynamics of NS-AS systems in three qualitatively different categories:

  1. In cases for which the AS has a small mass (and correspondingly low compactness), the NS is only weakly perturbed during the collision. Thus, only a small amount of matter is ejected and GW luminosity is small. The NS stabilises and becomes surrounded by a bosonic cloud which appears to be long lived.

  2. Increasing the AS mass, we obtain a merger remnant which is strongly excited, thus emitting GWs and ejecting large amounts of baryonic matter.

  3. In cases where the axion star is even more massive, the final remnant is a black hole. For these cases the GW luminosity is larger, but baryonic and bosonic ejecta are suppressed.

In the following we discuss in more detail these possible merger outcomes. For this purpose we show in Fig. 9 the bosonic and baryonic energy density, see Ref. Dietrich et al. (2018a), and the conformal factor for different times for the three setups (left row), (middle row), and (right row). Additionally, we also present the time evolution of the conformal factor in Fig. 10. This can be compared to the same cases for BHAS mergers, in Fig. 5.

Case I []: From the upper left panel we can immediately see qualitatively differences between bosonic and baryonic stars. While ASs do not have a sharp surface and are characterized by an exponentially decaying scalar field, the density at the NSs surface is only continuous. Considering we find that initially the central energy density of the AS is about orders of magnitude below the central density of the NS. However, the peak energy density significantly increases during the merger process becoming compatible to the NS’s density. Nevertheless, the overall change in the metric, cf. the evolution of the conformal factor, changes only slightly once the NS and the AS merge. In fact at around the minimum of the conformal factor returns to the initial value before the collision; see Fig. 10. Finally at the NS settles to a setup for which it is surrounded by a bosonic cloud, which again has maximum energy densities of the order of .

Case II []: For cases where the AS’s mass is larger than for Case I setups, we find that at merger the central bosonic energy density can be about order of magnitude larger. The central value of the conformal factor decreases from to , but no BH forms during the merger process. At this stage the NS’s central density increases and the radius of the NS shrinks, see middle panel of Fig. 9. We expect that the transition between BH formation and no BH formation is characterized by a type-I critical phenomenon as in the NSNS merger case Kellermann et al. (2010). We postphone a careful investigation employing a larger number of configurations for future work. Even at the remnant has not yet settled to a stable configuration and undergoes continuous oscillations (see green line in Fig. 10). Compared to the Case I setup, we find larger differences between resolutions R2 and R3, however, the qualitative shape agrees well, emphasising the robustness of the results.

Case III []: As the AS mass is further increased we find that (i) ASs are more compact with energy densities compatible with the NS density. In contrast to Case II setups, the pressure of the bosonic and baryonic matter is not able to counteract the gravitational collapse, which is clearly visible by the conformal factor becoming zero. At we find that for such cases there is no noticeable baryonic disk surrounding the final remnant and also the bosonic cloud is relatively small.

Figure 10: Evolution of the minimum value of the conformal factor . We present results for resolution R3 with solid lines and corresponding results for resolution R2 are shown dashed. Overall we find good agreement between different resolutions showing the robustness of the numerical methods.

iv.3 Bosonic cloud and baryonic outflow

Figure 11: Evolution of the bosonic cloud mass surrounding the final BH remnant. We present results for resolution R3 with solid lines and corresponding results for resolution R2 are shown dashed. We evaluate the cloud mass on the refinement level .

As for the BHAS mergers, we are interested in the potential formation of a BH surrounded by a cloud of bosonic particles. We present in Fig. 11 the cloud mass for the cases which undergo BH formation. In agreement with our previous investigations we find that the cloud mass is larger for the case of less compact ASs. In particular for the configuration the bosonic cloud reaches up to . We assume that close to the threshold of BH formation there will be a maximum achievable BH axionic cloud mass. We expect that as presented for the study of BH-AS systems an additional intrinsic spin of the NS or also the consideration of systems with orbital angular momentum will increase the cloud and disk mass.

Considering the amount of baryonic mass surrounding the final remnant, we find that even shortly after BH formation the baryonic disk mass drops below . Such small disk masses would generally not give rise to short gamma-ray-bursts (GRBs). We assume that the reason for this small disk mass is the missing angular momentum support due to the fact that we restrict our analysis to head-on collisions. Further studies investigating orbiting NS-AS simulations will be needed.

In addition to the baryonic disk acting as the central engine for a possible sGRB, another electromagnetic counterpart might be triggered by the neutron-rich outflow of the baryonic matter, namely a kilonova Lattimer and Schramm (1974); Li and Paczynski (1998); Metzger et al. (2010); Roberts et al. (2011); Kasen et al. (2017a). The kilonovae properties depend on the ejecta mass, geometry, and composition. Generally, more massive ejecta are able to trigger brighter kilonovae. Therefore an estimate of the ejected baryonic mass is of crucial importance. We show the ejecta mass estimates for all studied configurations in Fig. 12. For systems close to the threshold of BH formation the amount of ejected matter easily exceeds . Consequently this material should be able to create electromagnetic signatures which we will describe in detail in Dietrich et al. (2018b).

Figure 12: Evolution of the baryonic matter ejected during and after the merger process of the NSAS systems. Results for resolution R3 are shown solid, corresponding results for resolution R2 are dashed. We evaluate the ejecta mass on the refinement level .

iv.4 BH formation threshold

Figure 13: BH mass as a function of the initial central value of the AS. For the inset we show the black hole mass for supercritical configurations assuming a critical BH formation threshold at . Error bars refer to the difference between the resolutions R3 and R2.

With reference to the threshold of BH formation for the configurations in this section, we present the BH mass as a function of the initial central value of the AS in Fig. 13. As can be concluded from Fig. 10 the critical value for BH formation lies in . Motivated by studies of the critical collapse, see Gundlach and Martin-Garcia (2007) and references therein, we assume that the BH mass is proportional to


We fit the supercritical configurations to this function and obtain an aproximate threshold of BH formation of , as shown in Fig. 13. However, we would require many more simulations closer to the critical point to obtain accurate values. For type II collapses, as the BH mass should approach zero, and thus there is the possibility that at the threshold of BH formation all the baryonic and bosonic matter gets ejected from the system. However, we would expect that the collapse is type I, as is the case in NS-NS mergers Kellermann et al. (2010).

iv.5 Gravitational wave signal

Figure 14: GW metric multipoles for the three example cases of Fig. 9. Results for resolution R3 are shown solid, corresponding results for resolution R2 are dashed.

We close our discussion about NS-AS mergers by considering the emitted GW signal. The dominant (2,2)-mode of the GW signal is shown in Fig. 14 for the cases in Fig. 9. We find that for the Case I setup, e.g. , the emitted GW signals have amplitudes of about . Over the entire simulation the emitted GW energy is . For increasing AS masses the total released energy can be increased by orders of magnitudes. For setup we obtain amplitudes more than an order of magnitude larger than for and the emitted GW energy is about . Furthermore, the remnant is still highly dynamical and the GW luminosity has not decreased noticeably by the end of the simulation. Due to the highly dynamical postmerger regime we find that overall agreement between different resolutions is worse than for , cf. dashed and solid lines. However, different resolutions (including resolution R1) give similar results and lead to the same estimate of the emitted GW energy. Finally, for cases which form BHs during the evolution the amount of GW energy can grow up to . As an example we show the case for .

V Discussion

In this article we have presented what is, to the best of our knowledge, the first study of axion star collisions with black holes and neutron stars using full 3+1D numerical relativity simulations. Such a study seems timely in the gravitational wave astronomy era, in which multiple detections of compact binaries are expected in coming years Abbott et al. (2016c).

With respect to our black hole-axion star merger simulations, we have investigated the impact of the axion star’s compactness, and the BH spin, on the mass of the remnant bosonic cloud surrounding the black hole. Although in most of the considered cases of the axion star’s mass is absorbed by the black hole shortly after the merger, in favourable cases the remaining cloud can be as large as of the initial axion star mass, with a bosonic cloud of mass of and peak energy density of , comparable to that obtained in a superradiant build up. We find that the largest scalar clouds are generated for low compactness ASs and spinning black holes. We note that there appear to be particular combinations which are overall more efficient at producing large axion clouds. We speculate that this might be caused by the excitation of particular quasi-bound states of the black hole.

The presented results are important since they show (i) that axion star-black holes mergers can provide a dynamical mechanism for the formation of scalar hair around black holes and (ii) that faster spinning (but not yet extremal) black holes allow for relatively large cloud masses. The spinning case is especially interesting as it may provide the seed for a superradiant build up, which could lead to additional observable gravitational wave signatures post merger. However, superradiance requires extremal spins and an appropriate matching of the axion and BH masses, whereas the effects we observe here are in principle more general. It would be worth extending this study to a larger range of mass ratios, spins and spin orientations, to confirm the approximate trends observed in this paper and identify whether the proposal that particular mass ratios and spin combinations are favourable for forming clouds is consistent with a wider set of results. It would also be interesting to consider the effect of larger self interactions of the axion field, other values of and, in the longer term, interactions with baryonic matter in an accretion disc.

For our study of neutron star-axion star collisions, we restricted our investigations to the merger of axion stars of various compactnesses with a “typical” neutron star having a gravitational mass of and the SLy equation of state. We found that for the setups studied, there exists a critical mass threshold for the axion star required to form a BH during the collision. In the considered cases, the black hole formation is triggered by the axion star being perturbed within the potential well of the neutron star. Its collapse leads to a black hole within the neutron star, rather than collapse of the neutron star itself.

For sub-threshold axion star masses the merger remnant is a perturbed neutron star enveloped in an axion cloud. For super-threshold axion star masses the final remnant is a black hole with a scalar cloud surrounding it. We suggest that the black hole formation threshold may correspond to a type I critical phase transition, as in binary neutron star mergers, and therefore universality and scaling relations could exist near to the critical point. We present a first (although very approximate) estimate of the critical threshold parameter , but further simulations are required for a more stringent constraints on the critical parameters.

Interestingly, we found that in the marginally sub critical cases, a large amount of baryonic mass was released from the merger remnant due to the formation of shocks in the NS. These ejecta can give rise to a kilonova-like counterpart, such as ATF201gfo, e.g. Abbott et al. (2017e); Coulter et al. (2017); Cowperthwaite et al. (2017); Smartt et al. (2017); Kasliwal et al. (2017); Kasen et al. (2017b). The potential new type of transient produced by such a near-critical neutron star-axion star collision is discussed in more detail in Dietrich et al. (2018b). In cases where a black hole forms after the merger, the ejection of matter as well as the formation of a baryonic accretion disk or bosonic cloud is suppressed. However, in the most extreme case the final black hole remnant can be embedded in a bosonic cloud of mass .

In future, we plan to perform further numerical simulations in which we add a direct interaction between the axions and the neutron star fluid. Such couplings, which are necessary to correctly model the QCD axion, could also give rise to observable effects in NS-NS collisions that occur in a background of axions. Such an approach would allow us to further constrain the properties of axionic dark matter using observations of the merger of binary neutron stars within dark matter halos.

Vi Acknowledgements

KC thanks her GRChombo collaborators (, and masters student Sladja Radnovic for her initial work on the simulations. TD acknowledges support by the European Unions Horizon 2020 research and innovation program under grant agreement No 749145, BNSmergers. Computations have been performed on the supercomputer SuperMUC at the LRZ (Munich) under the project number pr48pu, the compute cluster Minerva of the Max-Planck Institute for Gravitational Physics, BSC Marenostrum IV via PRACE grant Tier0 PPFPWG, La Palma Astrophysics Centre via BSC/RES grants AECT-2017-2-0011 and AECT-2017-3-0009 and on SurfSara Cartesius under Tier-1 PRACE grant DECI14 14DECI0017, and the GWDG cluster in Goettingen.


Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description