Axion star collisions with black holes and neutron stars in full 3D numerical relativity
Abstract
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 headon collisions of relativistic ASs with black holes (BHs) and neutron stars (NSs). In the case of BHAS 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 NSAS 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, NSAS collisions could give rise to electromagnetic observables in addition to their gravitational wave signatures.
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 EinsteinKleinGordon 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 axionlike particles (ALPs) being a wellmotivated 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 selfinteraction 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 ultralight ALPs).
The original QCD axion emerges as a consequence of the PecceiQuinn (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, highamplitude 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 socalled 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 highmass tail of the distribution cannot be ruled out at present. A population of relativistic ASs could also have been produced by nonstandard primordial perturbations with enhanced smallscale power Widdicombe et al. (2018).
Ultralight 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ödingerPoisson 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 corehalo 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 SanchisGual 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 multimessenger events, with, for example, phenomena such as Fast Radio Bursts (FRBs). In the case of a NSAS 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 BHAS 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); SanchisGual 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 BHAS from mergers of ASs with spinning and nonspinning BHs.
In Sec. IV we discuss NSAS collisions, and
investigate the remnant of the headon collision as the NSAS 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 ^{1}^{1}1Note 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 BHAS evolutions, and the BAM code Brügmann et al. (2008); Thierfelder et al. (2011); Dietrich et al. (2018a) for the NSAS cases, but undertake comparisons to ensure they give consistent results.
ii.1 BHAS collisions using GRChombo
The GRChombo code Clough et al. (2015) (www.grchombo.org) is based on the methodoflines with fourth order spatial finite differences and fourth order explicit RungeKutta 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 nontrivial “manyboxesinmanyboxes” mesh hierarchies using block structured Berger Rigoutsos grid generation Berger and Rigoutsos (1991). The grid is made out of a hierarchy of cellcentered 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:
(1) 
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 BergerOliger scheme is employed Berger and Oliger (1984) to coordinate time stepping on the grid hierarchy and we use a CourantFriedrichsLewy factor of on each level. The evolution of the scalar fields is based on the KleinGordon 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.
Name  

R1  8  128  (1.0,0.20,0.06)  0.48  3e6  1e2  8.00  0.031250 
R2  8  192  (0.6,0.16,0.04)  0.32  3e6  1e2  5.33  0.020833 
R3  8  256  (0.5,0.10,0.03)  0.24  3e6  1e2  4.00  0.015625 
ii.2 NSAS 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 methodoflines, Cartesian grids and finite differencing. The grid is made out of a hierarchy of cellcentered 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 ‘movingboxes’ 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 bitantsymmetry to reduce the number of grid points and the computational costs by a factor of two. For the time integration a fourth order explicit RungeKutta type integrator is used with a CourantFriedrichsLewy factor of . For the time stepping of the refinement level the BergerOliger 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 highresolutionshockcapturing schemes based on primitive reconstruction and the local LaxFriedrich 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 KleinGordon 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 NSAS 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 KleinGordon equation which would be required to obtain fully consistent initial configurations.
Name  

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 
ii.3 Code testing and comparison
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 BHAS merger. We also refer to
Ref. Dietrich et al. (2018a) for a comparison of a BSBS headon
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
(2) 
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.
More challenging than the single star comparison is the evolution of a BHAS merger simulation. For this purpose we perform a headon 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 BHAS merger. The results are shown in Fig 3.
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 BHAS mergers
In this section we study the headon 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 BHAS 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 BHAS 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 (ArnowittDeserMisner) 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 BowenYork 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 zaxis) is perpendicular to the merger direction (along the xaxis), 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.
Name  

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 
iii.2 Qualitative merger dynamics
We find that the BHAS 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 NSAS systems into three qualitatively different categories depending on the compactness (and mass) of the AS with which they collide. In the BHAS 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)^{2}^{2}2Energy 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.
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 .
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 BHAS 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
(3) 
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 quasistable 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.
iii.4 Gravitational wave signals
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 NSAS mergers
Following the study of BHAS collisions, we will focus in this section on headon mergers of NSAS 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 NSAS 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 BHAS 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).
Name  

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 
iv.2 Qualitative merger dynamics
We can classify the merger dynamics of NSAS systems in three qualitatively different categories:

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.

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

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 typeI 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.
iv.3 Bosonic cloud and baryonic outflow
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 BHAS 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 gammaraybursts (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 headon collisions. Further studies investigating orbiting NSAS 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 neutronrich 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).
iv.4 BH formation threshold
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 MartinGarcia (2007) and references therein, we assume that the BH mass is proportional to
(4) 
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 NSNS mergers Kellermann et al. (2010).
iv.5 Gravitational wave signal
We close our discussion about NSAS 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 holeaxion 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 quasibound states of the black hole.
The presented results are important
since they show (i) that axion starblack 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 staraxion 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 subthreshold axion star masses the merger remnant is a perturbed neutron star enveloped in an axion cloud. For superthreshold 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 kilonovalike 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 nearcritical neutron staraxion 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 NSNS 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 (www.grchombo.org/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 MaxPlanck Institute for Gravitational Physics, BSC Marenostrum IV via PRACE grant Tier0 PPFPWG, La Palma Astrophysics Centre via BSC/RES grants AECT201720011 and AECT201730009 and on SurfSara Cartesius under Tier1 PRACE grant DECI14 14DECI0017, and the GWDG cluster in Goettingen.
References
 Abbott et al. (2016a) B. P. Abbott et al. (Virgo, LIGO Scientific), Phys. Rev. Lett. 116, 061102 (2016a), arXiv:1602.03837 [grqc] .
 Abbott et al. (2016b) B. P. Abbott et al. (Virgo, LIGO Scientific), Phys. Rev. Lett. 116, 241103 (2016b), arXiv:1606.04855 [grqc] .
 Abbott et al. (2017a) B. P. Abbott et al. (VIRGO, LIGO Scientific), Phys. Rev. Lett. 118, 221101 (2017a), arXiv:1706.01812 [grqc] .
 Abbott et al. (2017b) B. P. Abbott et al. (Virgo, LIGO Scientific), Phys. Rev. Lett. 119, 141101 (2017b), arXiv:1709.09660 [grqc] .
 Abbott et al. (2017c) B. P. Abbott et al. (Virgo, LIGO Scientific), Astrophys. J. 851, L35 (2017c), arXiv:1711.05578 [astroph.HE] .
 Abbott et al. (2017d) B. P. Abbott et al. (Virgo, LIGO Scientific), Phys. Rev. Lett. 119, 161101 (2017d), arXiv:1710.05832 [grqc] .
 Giudice et al. (2016) G. F. Giudice, M. McCullough, and A. Urbano, JCAP 1610, 001 (2016), arXiv:1605.01209 [hepph] .
 Wheeler (1955) J. A. Wheeler, Phys. Rev. 97, 511 (1955).
 Kaup (1968) D. J. Kaup, Phys. Rev. 172, 1331 (1968).
 Ruffini and Bonazzola (1969) R. Ruffini and S. Bonazzola, Phys. Rev. 187, 1767 (1969).
 Seidel and Suen (1991) E. Seidel and W. M. Suen, Phys. Rev. Lett. 66, 1659 (1991).
 Eby et al. (2016) J. Eby, C. Kouvaris, N. G. Nielsen, and L. C. R. Wijewardhana, JHEP 02, 028 (2016), arXiv:1511.04474 [hepph] .
 Colpi et al. (1986) M. Colpi, S. L. Shapiro, and I. Wasserman, Phys. Rev. Lett. 57, 2485 (1986).
 Krippendorf et al. (2018) S. Krippendorf, F. Muia, and F. Quevedo, (2018), arXiv:1806.04690 [hepth] .
 Marsh (2016) D. J. E. Marsh, Phys. Rept. 643, 1 (2016), arXiv:1510.07633 [astroph.CO] .
 Peccei and Quinn (1977) R. D. Peccei and H. R. Quinn, Phys. Rev. Lett. 38, 1440 (1977), [,328(1977)].
 Weinberg (1978) S. Weinberg, Phys. Rev. Lett. 40, 223 (1978).
 Hogan and Rees (1988) C. J. Hogan and M. J. Rees, Phys. Lett. B205, 228 (1988).
 Kolb and Tkachev (1993) E. W. Kolb and I. I. Tkachev, Phys. Rev. Lett. 71, 3051 (1993), arXiv:hepph/9303313 [hepph] .
 Schive et al. (2014) H.Y. Schive, T. Chiueh, and T. Broadhurst, Nature Phys. 10, 496 (2014), arXiv:1406.6586 [astroph.GA] .
 Veltmaat et al. (2018) J. Veltmaat, J. C. Niemeyer, and B. Schwabe, (2018), arXiv:1804.09647 [astroph.CO] .
 Levkov et al. (2018) D. G. Levkov, A. G. Panin, and I. I. Tkachev, (2018), arXiv:1804.05857 [astroph.CO] .
 Widdicombe et al. (2018) J. Y. Widdicombe, T. Helfer, D. J. E. Marsh, and E. A. Lim, (2018), arXiv:1806.09367 [astroph.CO] .
 Svrcek and Witten (2006) P. Svrcek and E. Witten, JHEP 06, 051 (2006), arXiv:hepth/0605206 [hepth] .
 Arvanitaki et al. (2010) A. Arvanitaki, S. Dimopoulos, S. Dubovsky, N. Kaloper, and J. MarchRussell, Phys. Rev. D81, 123530 (2010), arXiv:0905.4720 [hepth] .
 Hu et al. (2000) W. Hu, R. Barkana, and A. Gruzinov, Phys. Rev. Lett. 85, 1158 (2000), arXiv:astroph/0003365 [astroph] .
 Hui et al. (2017) L. Hui, J. P. Ostriker, S. Tremaine, and E. Witten, Phys. Rev. D95, 043541 (2017), arXiv:1610.08297 [astroph.CO] .
 Du et al. (2017) X. Du, C. Behrens, J. C. Niemeyer, and B. Schwabe, Phys. Rev. D95, 043519 (2017), arXiv:1609.09414 [astroph.GA] .
 Balakrishna (1999) J. Balakrishna, A Numerical study of boson stars: Einstein equations with a matter source, Ph.D. thesis, Washington U., St. Louis (1999), arXiv:grqc/9906110 [grqc] .
 Palenzuela et al. (2007) C. Palenzuela, I. Olabarrieta, L. Lehner, and S. L. Liebling, Phys. Rev. D75, 064005 (2007), arXiv:grqc/0612067 [grqc] .
 Choptuik and Pretorius (2010) M. W. Choptuik and F. Pretorius, Phys. Rev. Lett. 104, 111101 (2010), arXiv:0908.1780 [grqc] .
 Bezares et al. (2017) M. Bezares, C. Palenzuela, and C. Bona, Phys. Rev. D95, 124005 (2017), arXiv:1705.01071 [grqc] .
 Lai (2004) C.W. Lai, A Numerical study of boson stars, Ph.D. thesis, British Columbia U. (2004), arXiv:grqc/0410040 [grqc] .
 Palenzuela et al. (2008) C. Palenzuela, L. Lehner, and S. L. Liebling, Phys. Rev. D77, 044036 (2008), arXiv:0706.2435 [grqc] .
 Cardoso et al. (2016) V. Cardoso, S. Hopper, C. F. B. Macedo, C. Palenzuela, and P. Pani, Phys. Rev. D94, 084031 (2016), arXiv:1608.08637 [grqc] .
 Dietrich et al. (2018a) T. Dietrich, S. Ossokine, and K. Clough, (2018a), arXiv:1807.06959 [grqc] .
 Helfer et al. (2018) T. Helfer, E. A. Lim, M. A. G. Garcia, and M. A. Amin, (2018), arXiv:1802.06733 [grqc] .
 SanchisGual et al. (2018) N. SanchisGual, C. Herdeiro, J. A. Font, and E. Radu, (2018), arXiv:1806.07779 [grqc] .
 SanchisGual et al. (2017) N. SanchisGual, C. Herdeiro, E. Radu, J. C. Degollado, and J. A. Font, Phys. Rev. D95, 104028 (2017), arXiv:1702.04532 [grqc] .
 Hook and Huang (2018) A. Hook and J. Huang, JHEP 06, 036 (2018), arXiv:1708.08464 [hepph] .
 Eby et al. (2017) J. Eby, M. Leembruggen, J. Leeney, P. Suranyi, and L. C. R. Wijewardhana, JHEP 04, 099 (2017), arXiv:1701.01476 [astroph.CO] .
 Raby (2016) S. Raby, Phys. Rev. D94, 103004 (2016), arXiv:1609.01694 [hepph] .
 Iwazaki (2015) A. Iwazaki, Phys. Rev. D91, 023008 (2015), arXiv:1410.4323 [hepph] .
 Barranco et al. (2013) J. Barranco, A. C. Monteverde, and D. Delepine, Phys. Rev. D87, 103011 (2013), arXiv:1212.2254 [astroph.CO] .
 Dietrich et al. (2018b) T. Dietrich, F. Day, K. Clough, M. Coughlin, and J. Niemeyer, in prep. (2018b).
 Barranco et al. (2012) J. Barranco, A. Bernal, J. C. Degollado, A. DiezTejedor, M. Megevand, M. Alcubierre, D. Nunez, and O. Sarbach, Phys. Rev. Lett. 109, 081102 (2012), arXiv:1207.2153 [grqc] .
 Barranco et al. (2017) J. Barranco, A. Bernal, J. C. Degollado, A. DiezTejedor, M. Megevand, D. Nunez, and O. Sarbach, Phys. Rev. D96, 024049 (2017), arXiv:1704.03450 [grqc] .
 SanchisGual et al. (2016) N. SanchisGual, J. C. Degollado, P. Izquierdo, J. A. Font, and P. J. Montero, Phys. Rev. D94, 043004 (2016), arXiv:1606.05146 [grqc] .
 Herdeiro and Radu (2015) C. A. R. Herdeiro and E. Radu, Proceedings, 7th Black Holes Workshop 2014: Aveiro, Portugal, December 1819, 2014, Int. J. Mod. Phys. D24, 1542014 (2015), arXiv:1504.08209 [grqc] .
 Hod (2012) S. Hod, Phys. Rev. D86, 104026 (2012), [Erratum: Phys. Rev.D86,129902(2012)], arXiv:1211.3202 [grqc] .
 Brito et al. (2015a) R. Brito, V. Cardoso, and P. Pani, Lect. Notes Phys. 906, pp.1 (2015a), arXiv:1501.06570 [grqc] .
 Brito et al. (2015b) R. Brito, V. Cardoso, and P. Pani, Class. Quant. Grav. 32, 134001 (2015b), arXiv:1411.0686 [grqc] .
 Ferreira et al. (2017) M. C. Ferreira, C. F. B. Macedo, and V. Cardoso, Phys. Rev. D96, 083017 (2017), arXiv:1710.00830 [grqc] .
 Hannuksela et al. (2018) O. A. Hannuksela, R. Brito, E. Berti, and T. G. F. Li, (2018), arXiv:1804.09659 [astroph.HE] .
 Degollado and Herdeiro (2014) J. C. Degollado and C. A. R. Herdeiro, Phys. Rev. D90, 065019 (2014), arXiv:1408.2589 [grqc] .
 Clough et al. (2015) K. Clough, P. Figueras, H. Finkel, M. Kunesch, E. A. Lim, and S. Tunyasuvunakool, Class. Quant. Grav. 32, 245011 (2015), [Class. Quant. Grav.32,24(2015)], arXiv:1503.03436 [grqc] .
 Brügmann et al. (2008) B. Brügmann, J. A. Gonzalez, M. Hannam, S. Husa, U. Sperhake, et al., Phys.Rev. D77, 024027 (2008), arXiv:grqc/0610128 [grqc] .
 Thierfelder et al. (2011) M. Thierfelder, S. Bernuzzi, and B. Brügmann, Phys.Rev. D84, 044012 (2011), arXiv:1104.4751 [grqc] .
 (59) M. Adams, P. Colella, D. T. Graves, J. N. Johnson, N. D. Keen, T. J. Ligocki, D. F. Martin, P. W. McCorquodale, D. Modiano, P. O. Schwartz, T. D. Sternberg, and B. Straalen, .
 Berger and Rigoutsos (1991) M. J. Berger and I. Rigoutsos, IEEE Trans. Sys. Man & Cyber. 21, 1278 (1991).
 Berger and Oliger (1984) M. J. Berger and J. Oliger, J.Comput.Phys. 53, 484 (1984).
 Helfer et al. (2017) T. Helfer, D. J. E. Marsh, K. Clough, M. Fairbairn, E. A. Lim, and R. Becerril, JCAP 1703, 055 (2017), arXiv:1609.04724 [astroph.CO] .
 Dietrich et al. (2015a) T. Dietrich, S. Bernuzzi, M. Ujevic, and B. Brügmann, Phys. Rev. D91, 124041 (2015a), arXiv:1504.01266 [grqc] .
 Bernuzzi and Dietrich (2016) S. Bernuzzi and T. Dietrich, Phys. Rev. D94, 064062 (2016), arXiv:1604.07999 [grqc] .
 Bernuzzi et al. (2012) S. Bernuzzi, A. Nagar, M. Thierfelder, and B. Brügmann, Phys.Rev. D86, 044030 (2012), arXiv:1205.3403 [grqc] .
 Bowen and York (1980) J. M. Bowen and J. W. York, Phys. Rev. D 21, 2047 (1980).
 Gleiser et al. (1998) R. J. Gleiser, C. O. Nicasio, R. H. Price, and J. Pullin, Phys. Rev. D57, 3401 (1998), arXiv:grqc/9710096 [grqc] .
 Douchin and Haensel (2001) F. Douchin and P. Haensel, Astron. Astrophys. 380, 151 (2001), astroph/0111092 .
 Read et al. (2009) J. S. Read, B. D. Lackey, B. J. Owen, and J. L. Friedman, Phys. Rev. D79, 124032 (2009), arXiv:0812.2163 [astroph] .
 Dietrich et al. (2015b) T. Dietrich, N. Moldenhauer, N. K. JohnsonMcDaniel, S. Bernuzzi, C. M. Markakis, B. Brügmann, and W. Tichy, Phys. Rev. D92, 124007 (2015b), arXiv:1507.07100 [grqc] .
 Abbott et al. (2018) B. P. Abbott et al. (Virgo, LIGO Scientific), (2018), arXiv:1805.11579 [grqc] .
 Coughlin et al. (2018) M. W. Coughlin et al., (2018), arXiv:1805.09371 [astroph.HE] .
 Bauswein et al. (2017) A. Bauswein, O. Just, H.T. Janka, and N. Stergioulas, Astrophys. J. 850, L34 (2017), arXiv:1710.06843 [astroph.HE] .
 Radice et al. (2018) D. Radice, A. Perego, F. Zappa, and S. Bernuzzi, Astrophys. J. 852, L29 (2018), arXiv:1711.03647 [astroph.HE] .
 Annala et al. (2018) E. Annala, T. Gorda, A. Kurkela, and A. Vuorinen, Phys. Rev. Lett. 120, 172703 (2018), arXiv:1711.02644 [astroph.HE] .
 Most et al. (2018) E. R. Most, L. R. Weih, L. Rezzolla, and J. SchaffnerBielich, (2018), arXiv:1803.00549 [grqc] .
 Kellermann et al. (2010) T. Kellermann, L. Rezzolla, and D. Radice, Class. Quant. Grav. 27, 235016 (2010), arXiv:1007.2797 [grqc] .
 Lattimer and Schramm (1974) J. M. Lattimer and D. N. Schramm, The Astrophysical Journal Letters 192, L145 (1974).
 Li and Paczynski (1998) L.X. Li and B. Paczynski, The Astrophysical Journal Letters 507, L59 (1998).
 Metzger et al. (2010) B. D. Metzger, G. MartínezPinedo, S. Darbha, E. Quataert, A. Arcones, D. Kasen, R. Thomas, P. Nugent, I. V. Panov, and N. T. Zinner, Monthly Notices of the Royal Astronomical Society 406, 2650 (2010), arXiv:1001.5029 [astroph.HE] .
 Roberts et al. (2011) L. F. Roberts, D. Kasen, W. H. Lee, and E. RamirezRuiz, The Astrophysical Journal Letters 736, L21 (2011).
 Kasen et al. (2017a) D. Kasen, B. Metzger, J. Barnes, E. Quataert, and E. RamirezRuiz, Nature 551, 80 EP (2017a).
 Gundlach and MartinGarcia (2007) C. Gundlach and J. M. MartinGarcia, Living Rev. Rel. 10, 5 (2007), arXiv:0711.4620 [grqc] .
 Abbott et al. (2016c) B. P. Abbott et al. (Virgo, LIGO Scientific), Astrophys. J. 833, L1 (2016c), arXiv:1602.03842 [astroph.HE] .
 Abbott et al. (2017e) B. P. Abbott et al. (Virgo, FermiGBM, INTEGRAL, LIGO Scientific), Astrophys. J. 848, L13 (2017e), arXiv:1710.05834 [astroph.HE] .
 Coulter et al. (2017) D. A. Coulter et al., Science (2017), 10.1126/science.aap9811, [Science358,1556(2017)], arXiv:1710.05452 [astroph.HE] .
 Cowperthwaite et al. (2017) P. S. Cowperthwaite et al., Astrophys. J. 848, L17 (2017), arXiv:1710.05840 [astroph.HE] .
 Smartt et al. (2017) S. J. Smartt et al., Nature 551, 75 (2017), arXiv:1710.05841 [astroph.HE] .
 Kasliwal et al. (2017) M. M. Kasliwal et al., Science 358, 1559 (2017), arXiv:1710.05436 [astroph.HE] .
 Kasen et al. (2017b) D. Kasen, B. Metzger, J. Barnes, E. Quataert, and E. RamirezRuiz, Nature (2017b), 10.1038/nature24453, [Nature551,80(2017)], arXiv:1710.05463 [astroph.HE] .