Standardizing kilonovae and their use as standard candles to measure the Hubble constant

Michael W. Coughlin ,1,2 Tim Dietrich,3,4 Jack Heinzel,5,6 Nandita Khetan ,7 Sarah Antier ,8 Mattia Bulla,9 Nelson Christensen,5,6 David A. Coulter,10 and Ryan J. Foley10 1School of Physics and Astronomy, University of Minnesota, Minneapolis, Minnesota 55455, USA 2California Institute of Technology, 1200 East California Boulevard, MC 249-17, Pasadena, California 91125, USA 3Institut für Physik und Astronomie, Universität Potsdam, Haus 28, Karl-Liebknecht-Str. 24/25, 14476, Potsdam, Germany 4Nikhef, Science Park, NL-1098 XG Amsterdam, The Netherlands 5Artemis, Université Côte d’Azur, Observatoire Côte d’Azur, CNRS, CS 34229, F-06304 Nice Cedex 4, France 6Carleton College, Northfield, Minnesota 55057, USA 7Gran Sasso Science Institute (GSSI), I-67100 L’Aquila, Italy 8APC, UMR 7164, 10 rue Alice Domon et Lonie Duquet, F-75205 Paris, France 9Nordita, KTH Royal Institute of Technology and Stockholm University, Roslagstullsbacken 23, SE-106 91 Stockholm, Sweden 10Department of Astronomy and Astrophysics, University of California, Santa Cruz, California 95064, USA

(Received 2 August 2019; accepted 26 February 2020; published 9 April 2020) The detection of GW170817 is revolutionizing many areas of astrophysics with the joint observation of gravitational waves and electromagnetic emissions. These multimessenger events provide a new approach to determine the Hubble constant, thus, they are a promising candidate for mitigating the tension between measurements of type-Ia supernovae via the local distance ladder and the cosmic microwave background. In addition to the "standard siren" provided by the gravitational-wave measurement, the kilonova itself has characteristics that allow one to improve existing measurements or to perform yet another, independent measurement of the Hubble constant without gravitational-wave information. Here, we employ standardization techniques borrowed from the type-Ia community and apply them to kilonovae, not using any information from the gravitational-wave signal. We use two versions of this technique, one derived from direct observables measured from the light curve, and the other based on inferred ejecta parameters, e.g., mass, velocity, and composition, for two different models. These lead to A precise knowledge of the Hubble constant (H 0 ), to determine the expansion rate of the Universe, is one of the most important measurements driving the study of cosmology [1,2]. It has been known for a long time that the combined detection of gravitational waves (GWs) and their potential electromagnetic counterparts are useful for measuring the expansion rate of the Universe [3]. These measurements are interesting since the GW standard siren measurements of H 0 do not rely on a cosmic distance ladder and do not assume any cosmological model. This measurement has been made possible by the detection of GW170817 [4] and AT2017gfo, a "kilonova," which is thermal emission produced by the radioactive decay of neutron-rich matter synthesized from the ejecta of the compact binary coalescence at optical, near-infrared, and ultraviolet wavelengths [5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21]. The analysis of GW170817 and the km/s/Mpc (median and symmetric 68% credible interval), where degeneracy in the GW signal between the source distance and the weakly constrained angle of inclination between the total angular momentum of the binary and the line of sight dominated the H 0 measurement uncertainty [22]. It has been estimated that ∼50-100 GW events with identified optical counterparts would be required to have a H 0 precision measurement of ∼2% [23]. Of course, these searches are a significant observational challenge because one has to cover, over a short time interval, a large localization region, typically larger than the ∼20 square degrees for GW170817 [24,25]. The resulting H 0 measurements can be improved with, for example, high angular resolution imaging of the radio counterpart. Hotokezaka et al. [26] applied this technique for GW170817 and obtained H 0 = 68.9 +4.7 −4.6 km/s/Mpc. In this Rapid Communication, we employ techniques borrowed from the type-Ia supernova community to measure distance moduli based on kilonova light curves. When constraining distances with type-Ia supernovae, it is required to anchor them to some other "primary" distance indicators in order to calibrate their absolute magnitude. In this work, we are using modeled light curves, and here the absolute magnitude is given by these models, and so we do not require to know their distances a priori by some other technique. The method relies on differences in the modeled light curves due to ejecta parameters such as the ejecta mass (M ej ), ejecta velocity (v ej ), and lanthanide fraction (X lan ). Recently, it was also shown that a similar approach can be used for the bolometric luminosity to standardize kilonovae [28]. Qualitatively, the imprint of the ejecta properties can be modeled by the semianalytic methods of Arnett [29], and more broadly by these models' success at predicting the light curves for GW170817 (e.g., Ref. [19]). Arnett-based models for kilonovae require a power term, which has generally been taken to be P ∝ t β as appropriate for radioactivity models [30], in addition to having other parameters such as the ejecta mass, the energy (or equivalently the velocity), and the opacity. Most important to the overall luminosity is the heating rate per mass of the ejecta, determined by the product of intrinsic decay power and thermalization efficiency. Estimates of the thermalization efficiency based on simulations exist [31], although they still are the largest systematic error budget, as the mass scales roughly inversely with the powering level. From Ref. [29], the diffusion timescale is τ ∝ ( κM v ) 1/2 , and similarly, the density 3 ; in this way, all of these quantities affect the observables. For the simplest model where all of the energy is injected at t = 0, corresponding to the time of peak luminosity, then the luminosity as a function of time L(t ) ∝ L 0 e −t/τ , which implies that log[L(t )/L 0 ] ∝ −t/τ . This argues that the change in magnitude will be proportional to τ .
Based on this, we explore color-magnitude diagrams for kilonovae, with the idea that measuring time constants and colors may be useful for determining the underlying luminosities. In the left-hand panel of Fig. 1, we show these quantities for all of the spherical models made available in Ref. [27] plotting the i band minus K band [32]. A few trends stand out: As expected, the simulations with lower X lan have lower absolute K-band magnitudes than those with higher X lan , with 1-2 mag differences seen depending on the lanthanide fraction. The much larger effect is on the color. From the lowest to highest X lan , the i-band minus K-band color can vary by up to 5-6 mag. In addition, the velocities predominantly change the color, but at a much lower level, changing the color 1 mag. The trend with M ej is a clear increase in peak magnitude, which is true of all lanthanide fraction and ejecta velocity pairs. Moreover, the overall K-band peak luminosity increases as X lan increases (or as one looks to the right in the grid). There is a similar but smaller trend with velocity.
In the right-hand panel, we plot the change in luminosity, m 7 , between the peak and 7 days later in the K band. As X lan decreases, the effects on m 7 increase, with m 7 1 mag at low X lan and m 7 1 mag at higher X lan . We plot the peak K-band magnitudes versus ejecta mass for the available lanthanide fractions and ejecta velocities of the employed simulation set in the Supplemental Material [33] (this is essentially the same plot as Fig. 1, but with the points separated out by ejecta velocity and lanthanide fraction).
The clear linear structure in Fig. 1 motivates the potential for their standardization, similar to type-Ia supernovae (SNe Ia) light curves (see, e.g., Ref. [34]). In the case of type-Ia supernova (SN Ia) light curves, they typically reach a peak 17 days after explosion and then decay on a timescale of a few months. This motivates the use of the peak brightness, the time of the peak, and the "width" of the light curve as characteristic variables that can be compared, as realized early on [1,2]. Similar to the decline-rate parameter used in the SN Ia community (typically m 15 ), we will define a K-band decay parameter over 7 days as discussed above. This has the benefit of being measured from the observed light curve, with a downside of being tied to a particular filter and photometric system. It also requires the peak in this passband to be well measured, which is perhaps more straightforward in the near infrared where the peaks occur after a few days and therefore may be identified more easily. One downside might be that this band is less likely to be imaged in typical follow-up observations (i.e., before a kilonova has been confidently identified) because of the lack of infrared imagers on typical telescopes. Therefore, while our analysis makes one choice, there could be others more suitable for the particular observational situation.
For the moment, we will ignore the so-called K corrections that arise from the fact that the observed spectral energy distribution is redshifted by a factor (1 + z) and are effectively observing with filters that have been blueshifted in the rest frame of the kilonova. Given the local sensitivity volume for these kilonovae, this is reasonable, although their inclusion could be straightforward by adopting a spectral energy distribution. A similar concern is that we implicitly ignore dust reddening in the analysis. In general, the changes in color due to the evolution of the kilonova is much faster than the evolution of the dust reddening timescales. This color term can be disentangled by reasonably high-cadence, multicolor imaging to measure the dust component that follows a typical reddening law and is constant in time.
We adopt two models where we do regressions in three dimensions in order to fit a distance for the kilonovae. The first is based purely on observable quantities, while the second will be on quantities inferred from the model-based light-curve fitting. We note that these points in parameter space are given equal weights in the model, and therefore implicitly assume that all portions of the parameter space are equally likely. To do this, we use a Gaussian process regression (GPR)-based interpolation [35] to create a fit to the peak K-band magnitude for arbitrary parameters. This takes the form We employ a GPR-based interpolation instead of a linear fit due to the significant covariance between parameters. As discussed above, the quantities we fit are m 7 in the K band, m 7 in the i band, and the i minus K-band color at the peak in the K band. In order to compare with the apparent magnitudes, we therefore compute a distance modulus, where m app is the apparent magnitude and μ = 5 log 10 ( D 10 pc ) is the desired distance modulus (here, D is the distance). In principle, after the m 7 and color-dependent corrections have been applied, there will remain an intrinsic dispersion in the light curves, perhaps arising from location in the galaxy or perhaps dependent on the galaxy type.
We compare this "measured" fit to a fit based on "inferred" quantities of ejecta mass, ejecta velocity, and lanthanide fraction, We note that this is not specific to these observables; for example, Ref. [36] use ejecta mass M ej , the temperature at 1 day after the merger T 0 [37], the half-opening angle of the lanthanide-rich component (with = 0 and = 90 • corresponding to one-component lanthanide-free and lanthaniderich models, respectively), and the observer viewing angle θ obs (with cos θ obs = 0 and cos θ obs = 1 corresponding to a system viewed edge on and face on, respectively). We will compare these two models in the following. For consistency, we also employ a GPR-based interpolation here, although we  1) and (3) and to the color-magnitude diagram shown in Fig. 1 for the models in Ref. [27]. The "measured" values are derived from direct observables measured from the light curve, and the "inferred" values are from model-dependent ejecta parameters (mass, velocity, and lanthanide fraction). The bottom panel shows the difference between the computed values from the model and the fits for the simulations analyzed here.
have found that a linear fit based on these variables gives comparable results. When performing the analysis, we include an overall error during the fit that represents scatter from intrinsic variability in the kilonovae models. These errors are ∼0.7 mag for the "measured" case and ∼0.4 mag for the inferred case, which are derived from the median error reported by the GPR across the parameter space. Figure 2 shows the performance of the fits of Eqs. (1) and (3) compared to the models in Ref. [27]. In general, there is some systematic scatter in the "measured" case, indicating that the m 7 in the K band, m 7 in the i band, and the i minus K-band color at the peak in the K band are not sufficient to completely standardize across the parameter space. In all likelihood, the proper choice of variables for this purpose will become more apparent as detections are made. The fits broken up by lanthanide fraction and velocity for the inferred and measured cases are shown in the Supplemental Material [33].
We use the fit of Eqs. (1) and (3) and apply them to the posteriors on ejecta mass, velocity, and lanthanide fraction derived previously [39][40][41]. In these fits, we assume a 1.0 mag systematic uncertainty on the model; in this sense, we do not require that the models are "perfect" but instead encode some of the systematic uncertainties associated with them. We sample from the posteriors, in addition to the distributions for the measured parameters, to derive a constraint of D = 31 +17
Following the analysis of Ref. [22], we compute the corresponding values of the local Hubble constant for the kilonova  [27] and Bulla et al. [36] differ in several aspects, with the latter assuming parametrized opacities and temperatures, adopting different density profiles and ejecta geometries, and taking into account the interplay between the two ejecta components. While the use of two different kilonovae models is not a truly independent analysis, the fact that they are different but the corresponding distributions consistent with one another gives confidence that the analysis is not particularly model dependent.  [27] and Bulla [36] models kilonova-only analyses, and the Kasen et al. [27] combined GW-EM analysis are shown. The 1-and 2-σ regions determined by the "superluminal" motion measurement from the radio counterpart (blue) [26], Planck CMB (TT,TE,EE+lowP+lensing) (green) [45], and SHoES Cepheid-SN distance ladder surveys (orange) [46] are also depicted as vertical bands.
We also compute the corresponding values of the local Hubble constant for an analysis which combines the kilonova and GW-derived distances (using the high-spin posteriors, as done in Ref. [22]). Because the GW and kilonova data are independent, the posterior probabilities for the distances can simply be multiplied. For the Kasen et al. [27] inferred analysis, we find a combined measurement of H 0 = 78 +14 −9 km s −1 Mpc −1 , while for the measured analysis, the results are the same as for the GW analysis. Altogether, this proves that the kilonova measurement is competitive with the GW measurements to obtain an independent constraint on H 0 .
In this Rapid Communication, we have demonstrated how to use parameters derived from radiative transfer simulations to give distance measurements using kilonovae (see also Ref. [28] for an alternative approach). We have adopted the peak luminosity in the K band, the decay in the K band over 7 days, m 7 , and the i minus K-band colors. We have shown that these distance estimates are consistent with other measurements of GW170817's host galaxy directly and provide competitive measurements of H 0 to GW distance measurements alone.
These techniques will play an important role in an era of both large-scale optical surveys, and radio-to-x-ray followup efforts, in identifying the electromagnetic counterparts of compact binary coalescences. The Large Synoptic Survey Telescope (LSST) [47] and the Wide-Field InfraRed Survey Telescope (WFIRST) [48] in particular would be able to identify candidates with optical/NIR emission similar to GW170817 beyond 300 Mpc (i.e., the limit of current Advanced LIGO GW surveys for compact binary coalescences). Techniques such as those described here will play a key role in making these detections useful for studies of the Hubble constant and thus dark energy.
To explore this further, we examine the potential constraints over a sample of similar events. First of all, we assume that the 26% measurement here is representative of the typical width of the H 0 measurement from an individual event; this might be optimistic as AT2017gfo has a very well-sampled light curve. We also assume that the constraints converge as 1 √ N , where N is the number of kilonovae. Reference [23] notes that for the GW counterpart measurement, the typical tails in the distributions "average out" to make a smaller effective width than is usual for any given single event, so this may be a conservative assumption. Under these assumptions, the number of events required to make a 6% and 2% measurement is ∼20 and ∼170, respectively. For the sake of comparison, Ref. [23] found that 50-100 binary neutron stars with counterparts are required to make a 2% measurement. They also found that for binary neutron stars without counterparts, where a galaxy-catalog-based statistical method is required, ∼50 binary neutron stars are required to make a 6% measurement. In this sense, it is less powerful than binary neutron stars with counterparts but more so than binary neutron stars without counterparts. It is also likely that in practice, the error bars on a kilonova-only analysis will become smaller with time as the models improve. Further detections of kilonovae will also make it possible to perform fits for Eq. (1) directly using data rather than using models, which may also yield improved constraints.
There are other techniques that could be employed that may yield further constraints. Reference [28] uses a combination of semianalytic kilonovae models with ejecta mass and velocity fits to numerical relativity simulations to show correlations between the maximum bolometric luminosity and decline from peak luminosity for a simulated population of binary neutron star mergers, showing that it may be possible to incorporate either models or future observations into the standardization of bolometric light curves. In addition, Ref. [49] laid the groundwork for techniques for performing joint standard candle-standard siren measurements that bring together the analyses of GWs and their electromagnetic counterparts, demonstrating in particular their importance in cases where the selection effects from kilonova and GW observations can be significant.
In addition, these techniques can be used to augment GW detector calibration, which are currently limited to a few percent in amplitude and a few degrees in phase in the latest observing runs and are reaching their fundamental limits in performance [50,51]. The match between the wave form predicted by general relativity and the GW strain signal alone can calibrate a single detector's relative amplitude and phase responses to a GW as a function of frequency [52]. This GW strain signal unaccompanied by other messenger signals can also calibrate relative responses between two GW detectors. However, when augmented by independent measurements of the event's distance and inclination angle, the detectors' absolute amplitude and phase responses to GWs can be calculated. Therefore, our techniques could be a critical method in the effort to calibrate responses of GW detectors using astrophysical signals. Although there will need to be many more GW detections with electromagnetic counterparts to improve the precision of the astrophysical calibration measurements over the current existing in situ measurements, single events can also be used to improve and verify in situ measurements. Reference [52] showed that GW data alone can constrain the relative amplitude calibration uncertainty to less than 10% with 10-20 events (less than 1% with 1000-2000 events); when combining conservative constraints with broad distance and inclination constraints from electromagnetic data, it takes 400-600 events. Given the similar distance constraints assumed to those derived here, this is the right ballpark for this level of calibration uncertainty.
Further work is needed to understand how our restriction to a spherical geometry with a single-component would change for multiple components including possibly reprocessing between ejecta components (e.g., Ref. [53]). Discovering more kilonovae will improve the understanding of the effect of ejecta geometries, supporting the use of multiple components, and informing how they should be included. In addition, the standardization assumes that the radiative-transfer models are consistent to produce proper absolute magnitudes and colors (at least within the assumed error bars), which motivates continued work to improve the accuracy of the models and the grids that they are simulated on. The potential for more sources, at the price of higher systematics, with this method leads to a trade-off between the use of kilonovae or purely GW measurements; GW measurements will have relatively lower levels of systematics at the price of fewer objects to use. This method can be used in particular for any kilonovae detection with measurements of the host galaxy redshift, by way of fits of the ejecta parameters (e.g., Ref. [54]). It may also be possible to constrain the inclination using associated GRB detections or upper limits on potential subthreshold gamma-ray transients. As both of these analyses will be directly testable by future kilonovae observations, the utility of analyses of this type will be dependent on the coming comparisons.