# Measuring the Hubble Constant with a sample of kilonovae

Michael W. Coughlin, Sarah Antier, Tim Dietrich, Ryan J. Foley, Jack

Heinzel, Mattia Bulla, Nelson Christensen, David A. Coulter, Lina Issa, Nandita Khetan

► To cite this version:

Michael W. Coughlin, Sarah Antier, Tim Dietrich, Ryan J. Foley, Jack Heinzel, et al. Measuring the Hubble Constant with a sample of kilonovae. Nature Commun., 2020, 11 (1), pp.4129. 10.1038/s41467-020-17998-5. hal-02934021

HAL Id: hal-02934021

https://hal.science/hal-02934021

Submitted on 13 Dec 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
ARTICLE

https://doi.org/10.1038/s41467-020-17998-5

OPEN

<figure><img src="image_2.png" /></figure>

Check for updates

# Measuring the Hubble constant with a sample of kilonovae

Michael W. Coughlin<fnref n="1" />,²✉, Sarah Antier³, Tim Dietrich⁴,⁵, Ryan J. Foley<fnref n="6" />, Jack Heinzel⁷,⁸, Mattia Bulla<fnref n="1" />✉,⁹, Nelson Christensen⁷,⁸, David A. Coulter<fnref n="6" />, Lina Issa⁹,¹⁰ & Nandita Khetan<fnref n="1" />✉,¹¹

Kilonovae produced by the coalescence of compact binaries with at least one neutron star are promising standard sirens for an independent measurement of the Hubble constant ($H_0$). Through their detection via follow-up of gravitational-wave (GW), short gamma-ray bursts (sGRBs) or optical surveys, a large sample of kilonovae (even without GW data) can be used for $H_0$ constraints. Here, we show measurement of $H_0$ using light curves associated with four sGRBs, assuming these are attributable to kilonovae, combined with GW170817. Including a systematic uncertainty on the models that is as large as the statistical ones, we find $H_0 = 73.8_{-5.8}^{+6.3} \text{ km s}^{-1} \text{ Mpc}^{-1}$ and $H_0 = 71.2_{-3.1}^{+3.2} \text{ km s}^{-1} \text{ Mpc}^{-1}$ for two different kilonova models that are consistent with the local and inverse-distance ladder measurements. For a given model, this measurement is about a factor of 2-3 more precise than the standard-siren measurement for GW170817 using only GWs.

<footnote n="1">¹School of Physics and Astronomy, University of Minnesota, Minneapolis, MN 55455, USA. ²Division of Physics, Math, and Astronomy, California Institute of Technology, Pasadena, CA 91125, USA. ³APC, UMR 7164, 10 rue Alice Domon et Léonie Duquet, 75205 Paris, France. ⁴Institut für Physik und Astronomie, Universität Potsdam, Haus 28, Karl-Liebknecht-Str. 24/25, 14476 Potsdam, Germany. ⁵Nikhef, Science Park 105, 1098 XG Amsterdam, The Netherlands.</footnote>

<footnote n="6">⁶Department of Astronomy and Astrophysics, University of California, Santa Cruz, CA 95064, USA. ⁷Artemis, Université Côte d'Azur, Observatoire Côte d'Azur, CNRS, CS 34229, F-06304 Nice Cedex 4, France. ⁸Physics and Astronomy, Carleton College, Northfield, MN 55057, USA. ⁹Nordita, KTH Royal Institute of Technology and Stockholm University, Roslagstullsbacken 23, SE-106 91 Stockholm, Sweden. ¹⁰Département de Phyisque, Université Paris-Saclay, ENS Paris-Saclay, 91190 Gif-sur-Yvette, France. ¹¹Gran Sasso Science Institute (GSSI), I-67100 L'Aquila, Italy. ✉email: cough052@umn.edu</footnote>
Since the discovery of the accelerating expansion rate of the universe¹,², cosmology surveys have tried to measure the properties of dark energy. One of the most common metrics, type Ia supernovae (SNe), which are standardizable, have been an important tool in this endeavour, with the particular benefit of being detectable throughout a large portion of cosmic time. It has been previously found that the cosmic microwave background (CMB) is consistent with $\Lambda_{\text{CDM}}$ cosmology, but predicts a value for $H_0$ in direct tension with other measurements³. The redshifts of type Ia SNe in hosts with distances, already determined according to Cepheid variables⁴, were used in combination with Hubble Space Telescope imaging⁵ to obtain a value $4.4\sigma$ distinct from the Planck Collaboration measurement. It is not yet clear whether this tension is due to the experimental procedures themselves—perhaps rooted in some hidden systematic error—or if it indicates a more exotic physics; additional independent measurements are necessary to assess the true source of the tension.

One of the possible independent measurement methods for $H_0$ connects to the multi-messenger observation of compact binary mergers in which at least one neutron star is present. This approach has been vitalized by the recent combined detection of the neutron star merger (BNS) GW170817⁶, GRB 170817A⁷,⁸, and the optical transient AT2017gfo⁹, found in the galaxy NGC 4993 12 h after the GWs and GRB. In addition to the resulting insight into the equation of state (EOS) of neutron stars¹⁰ and the formation of heavy elements¹¹, one of the most exciting results was that of the $H_0$ measurement¹². This measurement is particularly powerful because GWs are standard sirens¹³, which do not rely on a cosmic distance ladder and do not assume any cosmological model as a prior (outside of assuming general relativity is correct). The combination of the distance measurement by the GWs and redshift from the electromagnetic counterpart makes constraints on $H_0$ possible¹⁴. The distance ladder independent measurement using GWs and the host redshift was $H_0 = 68_{-8}^{+18}$ km/s/Mpc (68.3% highest density posterior interval with a flat-in-log prior)¹²; inclusion of all O2 events reduced this uncertainty to $H_0 = 68_{-7}^{+14}$ km/s/Mpc¹⁵. Improvements on this measurement using more electromagnetic information, such as high angular resolution imaging of the radio counterparts¹⁶ or information about the internal composition of the NSs¹⁷, are also possible.

Here, we show that the electromagnetic evolution of kilonovae—particularly their decay rate and color evolution—can be compared to theoretical models to determine their intrinsic luminosity, making kilonovae standardizable candles¹⁸,¹⁹. Along with the measured brightnesses, kilonovae can be used to measure cosmological distances. We apply two kilonovae models to sGRB light curves to constrain $H_0$ to $H_0 = 73.8_{-5.8}^{+6.3}$ km s⁻¹ Mpc⁻¹ and $H_0 = 71.2_{-3.1}^{+3.2}$ km s⁻¹ Mpc⁻¹ for the two models, improving on what has been achieved so far with GWs alone by about a factor of 2-3.

## Results

**Measuring $H_0$ using kilonovae.** Here, we focus on the kilonova observation happening in coincidence with sGRBs. This type of analysis is particularly prescient given the difficulty of searches for GW counterparts during Advanced LIGO and Advanced Virgo's third observing run (O3)²⁰. AT2017gfo, synthesized by the radioactive decay of r-process elements in neutron-rich matter ejected during the merger²¹,²², is certainly the best sampled kilonova observation to date. Significant theoretical modeling prior to and after GW170817 has made it possible to study AT2017gfo in great detail, including measurements of the masses, velocities, and compositions of the different ejecta types. These measurements rely on models employing both simplified

semi-analytical descriptions of the observational signatures²³ and modelling using full-radiative transfer simulations²⁴,²⁵.

In addition to the observation of AT2017gfo, GW170817 was associated with GRB 170817A, which proved that at least some of the observed sGRBs are produced during the merger of compact binaries. This multi-messenger observation revealed the possible connection between kilonovae and sGRBs. For both cases, the GRB is then followed by an afterglow visible in X-rays, optical, and radio for days to months after the initial prompt γ-ray emission derived from the shock of the jet with the external medium. Our sGRB/kilonova sample follows ref. ²⁶, which combined state-of-the-art afterglow and kilonova models, jointly fitting the observational data to determine whether there was any excess light from a kilonova. The analysis showed light curves consistent with kilonovae in the cases of GRB 150101B²⁷, GRB 050709²⁸, GRB 160821B²⁹, and GRB 060614³⁰. Naturally, the error bars on the kilonova parameters are larger for these objects than for GW170817, which have light curves with potentially significant contamination from the afterglow. We refer the reader to ref. ²⁶ for extensive discussions of the photometric data quality and light-curve parameters and modeling. On top of these GRB observations, we will also include measurements from GW170817 (GRB 170817A)¹⁸. While no spectra of the kilonova excesses exist and X-ray excesses may point to shock heating driving these near-infrared emission²⁹, kilonovae are one possible (if not even the most likely) interpretation of the excesses. We point out that for the purpose of this article, we assume that the light curves are solely caused by a kilonova emission and neglect the possible contamination due to the sGRB afterglow. While this assumption leads to possible biases if strong sGRB afterglows would contaminate the observed data, adding an sGRB afterglow model on top of the kilonova light curves would increase the dimensionality of the analysis significantly and thus no (or only limited) constraints could be obtained.

The idea is to use techniques borrowed from the type-Ia SNe community to measure distance moduli based on kilonova light curves. We use the light-curve flux and color evolution, which do not depend on the overall luminosity, compared to kilonova models, to predict the luminosity; when combined with the measured brightness, the distance is constrained (see Methods). Here, we develop a model for the intrinsic luminosity of kilonovae based on observables, such that the luminosity can be standardized. Given the potential of multiple components and the change in color depending on the lanthanide fraction, it is useful to use kilonova models to perform the standardization. While it may be possible to standardize the kilonova luminosities based on measured properties, as is done for SN Ia cosmology measurements, in this analysis, we assume that we can use quantities inferred from the light-curve models³¹; this assumption will be testable when a sufficiently large sample of high-quality kilonovae observations are available. In this analysis, we use models from Kasen et al.²⁴ and Bulla²⁵. We note that this analysis uses some of the sampling techniques in the kilonova hypothesis testing and parameter estimation as demonstrated in refs. ¹⁰,²⁶, but this is a fundamentally orthogonal exercise to the use of observations in one particular band to standardize the light curves based on measurements from theoretical models. The use of one of the kilonovae such as GW170817 to inform standardization could be used however with unknown systematic errors. While other models for kilonovae exist at this point, e.g., ref. ³², they have been shown to give similar light-curve fits to GW170817, and so the expectation is they would give similar results to those here.

We measure the distances to both GW170817 and the sGRBs in our sample using the posteriors for model parameters and the distributions for the measured parameters from the fit. For the sGRBs, we use two distance estimates based on GW170817 to
inform the standardization; we use GW170817's distance combined with the difference between the computed distance moduli to extract the distance moduli for the sGRBs (see Methods). One powerful aspect of this is that for GW170817 in particular, surface brightness fluctuations (SBF) of the host galaxy NGC 4993 (blue)³³ pin the distance to within 1 Mpc. This requirement is similar to SN Ia measurements, where local distance ladders are required to calibrate the measurement. We also perform a comparison where we use the GW-derived posteriors to anchor the distance distribution, which results in broader but consistent posteriors (see Methods). From there, the distance modulus for each sGRB is solved for, resulting in a distribution of distances.

**H₀ constraints.** In addition to the study of GW170817/GRB 170817A/AT2017gfo¹⁸, we compute the corresponding values of H₀ for the sGRBs²⁶. As described above, we use the Kasen et al.²⁴ and Bulla²⁵ models, and we assume systematic error bars with 0.1 mag and 0.25 mag errors for comparison; these are chosen to be similar to photometric errors (0.1 mag) and twice as large (0.2 mag) to establish robustness. This broadens the posterior distributions on the ejecta parameters and these systematic errors also reweight the dependence of the eventual H₀ measurement on individual objects. Results for all kilonovae and the combined analysis are listed in Table 1 and Fig. 1. To determine the

posterior distribution for H₀, we perform a simultaneous fit for two cosmology models. The first is an empirical model taken to be the following⁴

$$D = \frac{z \times c}{H_0} \left( 1 + \frac{1}{2}(1 - q_0)z - \frac{1}{6}(1 - q_0 - 3q_0^2 + j_0)z^2 + O(z^3) \right). \quad (1)$$

The second is Λ<sub>CDM</sub>, which depends on H₀, Ω<sub>m</sub>, and Ω<sub>Λ</sub>. We checked that both analyses give similar constraints on H₀, but do

<table><caption>Table 1 Summary of H₀ results.</caption><thead><tr><th>Kilonova</th><th>Kasen-<br>0.25</th><th>Kasen-<br>0.1</th><th>Bulla-<br>0.25</th><th>Bulla-0.1</th></tr></thead><tbody><tr><td>GW170817</td><td>75<sup>+9</sup><sub>-8</sub></td><td>78<sup>+9</sup><sub>-8</sub></td><td>75<sup>+7</sup><sub>-7</sub></td><td>75<sup>+6</sup><sub>-6</sub></td></tr><tr><td>GRB 060614</td><td>68<sup>+22</sup><sub>-16</sub></td><td>66<sup>+13</sup><sub>-12</sub></td><td>71<sup>+33</sup><sub>-20</sub></td><td>66<sup>+16</sup><sub>-11</sub></td></tr><tr><td>GRB 150101B</td><td>66<sup>+40</sup><sub>-20</sub></td><td>76<sup>+39</sup><sub>-20</sub></td><td>89<sup>+54</sup><sub>-39</sub></td><td>96<sup>+51</sup><sub>-38</sub></td></tr><tr><td>GRB 050709</td><td>78<sup>+29</sup><sub>-19</sub></td><td>63<sup>+14</sup><sub>-11</sub></td><td>62<sup>+93</sup><sub>-7</sub></td><td>61<sup>+8</sup><sub>-6</sub></td></tr><tr><td>GRB 160821B</td><td>64<sup>+11</sup><sub>-10</sub></td><td>63<sup>+12</sup><sub>-9</sub></td><td>62<sup>+8</sup><sub>-6</sub></td><td>62<sup>+8</sup><sub>-6</sub></td></tr><tr><td>Combined</td><td>71.9<sup>+6.4</sup><sub>-5.5</sub></td><td>73.8<sup>+6.3</sup><sub>-5.8</sub></td><td>69.9<sup>+3.6</sup><sub>-3.7</sub></td><td>71.2<sup>+3.2</sup><sub>-3.1</sub></td></tr></tbody></table>

We use units of km s⁻¹ Mpc⁻¹. Individual rows refer to the GRB/GW observations and individual columns to the Kasen et al. and Bulla et al. model assuming a 0.25 and 0.1 mag systematic uncertainty. We note that the GRB individual H₀ measurements use the SBF of the host galaxy NGC 4993³³ to pin the distance.

<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

**Fig. 1 Posterior distributions for H₀ for individual events.** We show GW170817, GRB 060614, GRB 150101B, GRB 160821B, GRB 050709, and their combined posteriors. Fig. a is the Kasen²⁴ model with 0.1 mag errors, b is the Bulla²⁵ model with 0.1 mag errors, c is the Kasen model with 0.25 mag errors, and d is the Bulla model with 0.25 mag errors. The 1- and 2-σ regions determined by the superluminal motion measurement from the radio counterpart (blue)¹⁶, Planck CMB (TT,TE,EE+lowP+lensing) (green)³⁴ and SHoES Cepheid-SN distance ladder surveys (orange)⁵ are also depicted as vertical bands.
<figure><img src="image_7.png" /><figcaption>**Fig. 2 Summary of posterior distributions for $H_0$.** In a, we show the GW-only analysis for GW170817, in addition to the Kasen²⁴ and Bulla²⁵ model analyses with 0.1 mag and 0.25 mag errors from this letter. In addition, the 1- and 2-$\sigma$ regions determined by the superluminal motion measurement from the radio counterpart (blue)¹⁶, Planck CMB (TT,TE,EE+lowP+lensing) (green)³⁴, and SHoES Cepheid-SN distance ladder surveys (orange)⁵ are also depicted as vertical bands. In b, we use GRB 060614 and the Kasen²⁴ model with 0.1 mag errors. In addition to the standard SBF measurement in the letter, we show an analysis where we systematically add 0.1 $M_\odot$ to $M_{ej}$ and 0.5 decades to $X_{lan}$. Finally, we show a distribution where we replace the standard SBF measurement with the distance measurement from the GW170817 high spin posteriors.</figcaption></figure>

not significantly constrain the other model parameters. The systematics from using a particular kilonova model will remain, but the idea is that a sufficient sample of kilonovae will average out variations in the kilonovae. However, the relative consistency between the results of the two models²⁴,²⁵ yields some confidence that we are still statistics dominated. We caution that the systematic uncertainty could still be significantly larger than what is assumed here, but the expectation is that the difference between the models should be resolved with future observations. The final, combined $H_0$ measurement of our analysis is $H_0 = 73.8_{-5.8}^{+6.3}$ km s⁻¹ Mpc⁻¹ for the Kasen model and $H_0 = 71.2_{-3.1}^{+3.2}$ km s⁻¹ Mpc⁻¹ for the Bulla model. These improve on what has been achieved so far with GWs alone by about a factor of 2–3 (see left panel of Fig. 2); the results are consistent with both Planck CMB³⁴ and SHoES Cepheid-SN distance ladder surveys analyses⁵. We also performed the same analysis without GW170817, due to possible systematic uncertainties from the peculiar velocity of the host; this analysis resulted in both larger error bars than the analysis with GW170817, while still being consistent with it (see Methods).

We can use our results to construct a so-called Hubble–Lemai^tre Diagram (see Fig. 3), where we plot distance modulus vs. redshift for the observed kilonovae. The main benefit of plotting distance modulus instead of apparent magnitude is that it is independent of the source. For comparison, the green dashed line shows $\Lambda_{CDM}$, showing the consistency in the results. We include a Hubble–Lemai^tre residual panel to show the error bars. The error bars are, of course, large relative to SNe samples, where $\sigma \sim 0.1$ mag.

Discussion

We perform a few different tests to assess the systematics of the analysis (see right panel of Fig. 2). The first is where we systematically change the estimated values for the ejecta mass and lanthanide fractions from the light-curve analysis to assess the dependence on those values. Based on ref. ³⁵, we take values of 0.1 $M_\odot$ to add to the ejecta mass and 0.5 decades to add to $X_{lan}$, corresponding to the approximate size of the error bars on those parameters. We show these curves along with the original analysis of GRB 060614 on the right of Fig. 2. While the shift in the derived distribution is clear, in particular for $X_{lan}$, the distributions

<figure><img src="image_8.png" /><figcaption>**Fig. 3 Hubble–Lemai^tre Diagram for the kilonova analysis.** We plot distance modulus vs. redshift. The error bars indicate the extent of the posterior samples. The green dashed line shows $\Lambda_{CDM}$. We include a Hubble–Lemai^tre residual panel to show the residuals.</figcaption></figure>

still remain consistent with one another; we find distributions of $67_{-11}^{+13}$ for SBF, $67_{-15}^{+19}$ for SBF—$M_{ej}$, and $62_{-11}^{+13}$ for SBF—$X_{lan}$. This indicates that we are still dominated by statistical errors. As a further test, we show a distribution where we replace the standard SBF measurement with the distance measurement from the high spin posteriors presented in ref. ³⁶. Given the wide posterior distributions, this leads to a broadening of the $H_0$ measurement, but again, the distributions are consistent with one another. In the following, we will use the SBF measurement as the distance anchor for the analysis, but in the future, the GW based distributions may be appropriate to minimize systematics.

These results indicate that continued searches for and the analyses of sGRB afterglows have significant science benefits. In particular, our analysis shows that $H_0$ measurements may be improved given more detections of sGRB afterglows and their peak luminosities, requiring detections as early as possible to be most usable. These objects also have observational challenges, at least those identified by the Fermi Gamma-ray Burst Monitor³⁷, given the large sky areas that require coverage and the candidate vetting that follows³⁸. Observations of more kilonovae will
significantly constrain the kilonova peak luminosity distribution$^{26}$, which is clearly a driving force of this analysis. This will likely depend on both whether the original system producing a kilonova is a binary neutron star or neutron star-black hole, in addition to the inclination angle influence on the kilonova characteristics. Kilonovae, in particular, are less bright than SNe Ia, and could be more useful for constraining other distance indicators rather than directly as cosmological probes.

## Methods

**Kilonova analysis.** In addition to the photometric error bars arising from the measured signal to noise, the modeling also has associated errors, which we will add in quadrature. Building upon the first analysis$^{18}$, where we employed a Gaussian Process Regression (GPR) based interpolation$^{39}$ to create a surrogate model of the model of Kasen et al.$^{24}$ for arbitrary ejecta properties$^{10,35}$, we also use the model of Bulla$^{25}$ for comparison. The idea is that we can use these models to derive constraints on possible kilonova light curves. We chose to have large prior boundaries that allows to describe both, kilonova produced by BNSs and by BHNS systems. For the Kasen et al. model, each light curve depends on the ejecta mass $M_{ej}$, the mass fraction of lanthanides $X_{lan}$, and the ejecta velocity $v_{ej}$. We use flat priors for each parameter covering: $-3 \le \log_{10}(M_{ej}/M_{\odot}) \le -1$, $0 \le v_{ej} \le 0.3$ c, and $-9 \le \log_{10}(X_{lan}) \le -1$. For the 2D$^{25}$ model, each light curve depends on the ejecta mass $M_{ej}$, the half-opening angle of the lanthanide-rich component $\Phi$ (with $\Phi = 0$ and $\Phi = 90^\circ$ corresponding to one-component lanthanide-free and lanthanide-rich models, respectively) and the observer viewing angle $\theta_{obs}$ (with $\cos \theta_{obs} = 0$ and $\cos \theta_{obs} = 1$ corresponding to a system viewed edge-on and face-on, respectively). We again use flat priors for each parameter covering: $-3 \le \log_{10}(M_{ej}/M_{\odot}) \le -1$, $15^\circ \le \Phi \le 30^\circ$, and $0 \le \theta_{obs} \le 15$. We restrict $0^\circ \le \theta_{obs} \le 15^\circ$ for the sGRB analysis because the viewing angle is much closer to the polar axis than for GW170817$^{40}$. By observing the sGRB, we can assume that we are near or within the opening angle of the sGRB jet, which is taken to be less than $15^{\circ}$41–44. Because the distribution of the viewing angles of kilonovae from sGRBs are likely quite anisotropic, we would expect this to create an appearance of changing lanthanide fractions as the viewing angle changed for spherical geometries, such as in the model of Kasen et al.$^{24}$; this could cause a bias in the Hubble Constant measurements using spherical models. Asymmetric models such as that of Bulla$^{25}$ overcome this potential issue.

Compared to the models presented in ref. $^{25}$, those used here adopt thermalization efficiencies from ref. $^{45}$ and estimate the temperature at each time from the mean intensity of the radiation field in each region of the ejecta. In

addition, we assume that no mass is located below $v_{min} = 0.025$ c (where c is the speed of light) following ref. $^{46}$. Studying the predictions of two independent models with physical assumptions allows us to estimate the systematic uncertainties of our analysis. More specifically, some of these physical assumptions (e.g., spherical or axial symmetry) may give the false impression of special kilonovae properties. Doing this study with multiple different models is therefore critical to reveal such systematics.

In addition to the attempt of providing a measure of the systematic uncertainty, our parameter ranges for both models are very agnostic in the sense that we do not restrict us to particular parameter ranges that are predicted by numerical-relativity simulations$^{35,47}$, in fact, the proposed method works for an even larger parameter space, thus, it seems to be very general approach that can be employed for a variety of future events.

**Light curves.** All data presented here were compiled from public sources and collated as presented in ref. 10 for GW170817 and ref. 26 for the remaining SGRBs. For GRB 150101B, data can be found in refs. $^{27,40}$, for GRB 050709, data can be found in refs. $^{28,48–50}$, for GRB 160821B, data can be found in ref. $^{29,51}$, and in GRB 060614, data can be found in refs. $^{30,52,53}$. We remind the reader that for the SGRB analyses, we take the peak r-band observation for comparison, as opposed to the K-band as used in the GW170817 analysis due to the sparsity of available light curves in that band. We perform the parameter fits using the entire light curves, assuming the light curves are the result of kilonovae. While GW170817 remains the only completely unambiguous kilonova detection, ref. 26 and other analyses, such as refs. $^{54,55}$ for GRB 160821B, have indicated that these SGRBs are consistent with having kilonovae present. On the other hand, GRB 060614 in particular has a number of possible interpretations, including a tidal disruption event$^{56}$, a WD-NS system$^{57}$ or a long GRB$^{58}$. We include it for completeness, but further analysis may require its removal in future analysis.

**H₀ analysis.** The underpinning assumption for this paper is that kilonovae can be standardized using models for their luminosity and color evolution, which was first explored in ref. $^{18}$. While in principle this method could be used for other transient types which are more numerous, including core collapse or superluminous supernovae, models for their emission properties are not nearly developed as those for kilonovae. This requires, of course, believable models that we expect encode the evolution of the kilonovae we measure. Metzger$^{59}$ showed that the semi-analytic methods of Arnett$^{60}$ already reproduce much of the expected physics for these systems, and these models are sufficient for predicting the GW170817 lightcurves (e.g., ref. $^{61}$). For radioactivity models$^{59}$, it assumes a power term taken to be $P \propto t^\beta$, and is also dependent on the energy in the system (related to the velocity $v$), the ejecta mass $M$, and the opacity $\kappa$. Under these assumptions, the luminosity as a

<figure><img src="image_9.png" /><figcaption>**Fig. 4** Fit of Eq. (3) for the models in ref. $^{25}$. We vary the opening and viewing angles of the employed simulation set in r-band. The black points are the fits (with the measured error bar from the Chi-squared) while the blue points are computed directly from the models in ref. $^{25}$.</figcaption></figure>
function of time evolves as log (L(t)/L₀) ∝ −t/τ, where τ is the diffusion timescale τ ∝ (kM/v) 1/2 18 . Numerical models build upon these semi-analytic models with line- based opacities, radiative transport, thermalization efficiency and other parameters to make these models further realistic, although the key point that luminosity depends on these measurable parameters remains.

We now explain explicitly how to derive the standardization. For each model in the simulation set, each of which corresponds to a specific set of intrinsic parameters, we compute the peak magnitude in each passband. We choose K- band for the GW170817 analysis when analyzed on its own, given that the transient was observable for the longest in the near-infrared; we chose r-band for the sGRB analyses, including in the standardization of GW170817 in the magnitude differences computed below, as it was the band most commonly imaged between the various sGRBs. In this way, we have for both the Kasen et al. 24 model and the Bulla 25 model, a grid of the intrinsic parameters and the peak magnitude associated with the simulation. To map the intrinsic parameters to a peak magnitude, we use a GPR based interpolation (similar to the light- curve interpolation described above) $^{39}$.

<figure><img src="image_10.png" /><figcaption>Fig. 5 Posterior distributions for distance to GW170817. We show the results of the GW-only analyses (high spin) $^{36}$, the Kasen et al. 24 and the Bulla 25 kilonova analysis for both systematic errors assumed. The 1- and 2- σ regions determined by the surface brightness fluctuations (SBF) of NGC 4993 (blue) 33 and the Fundamental Plane (FP) of E and SO galaxies (red) 62 are also depicted as vertical bands.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 6 Corner plots for the model parameters. We show them for both Kasen et al. 24 model (a) and the Bulla 25 model (b) for GRB 060614. We also show the inferred distance modulus in the last column.</figcaption></figure>

For the Kasen et al. 24 model, the equation takes the form

$$
M_{r=\tau_{\text{max}}} = f(\log_{10} M_{\text{ej}}, v_{\text{ej}}, \log_{10} X_{\text{lan}}) \quad (2)
$$

and for the Bulla 25 model,

$$
M_{r=\tau_{\text{max}}} = f(\log_{10} M_{\text{ej}}, \Phi, \theta_{\text{obs}}) \quad (3)
$$

where f is a GPR based interpolation and the parameters are inferred quantities based on the light-curve fits. The idea is that these magnitude fits are related to the (observed) apparent magnitudes by the distance modulus µ = 5log$_{10}$( D 10pc ).

$$
M = m - \mu. \tag{4}
$$

We can evaluate the performance of these fits by comparing the peak K-band magnitudes to those predicted as a function of ejecta mass. Figure 4 uses the simulation set for the Bulla 25 model, showing the performance as a function of viewing and opening angle. While the performance tends to be worst for the most extreme viewing and opening angles, the fact that the Gaussian Process estimates errors, which change across the parameter space help to sufficiently reproduce the behavior.

We use the fit of Eqs. (2) and (3) and apply it to the posteriors on the model parameters, as were derived previously for GW170817 10,35 and the sGRBs $^{26}$. Applying directly to the GW170817 posteriors, we show the estimated distances for both models in Fig. 5, consistent with other measurements of the host galaxy, e.g., refs. 33,62,63, and the GW posteriors $^{36}$.

To apply the model to the sGRB light curves, we combine Eqs. (2) and (4) as

$$
\mu = m - f(\log_{10} M_{\text{ej}}, v_{\text{ej}}, \log_{10} X_{\text{lan}}) \quad (5)
$$

and Eqs. (3) and (4) as

$$
\mu = m - f(\log_{10} M_{\text{ej}}, \Phi, \theta_{\text{obs}}). \quad (6)
$$

We then take the difference between two observations,

$$
\mu_1 - \mu_2 = m_1 - m_2 + f(\log_{10} M_{ej,1}, v_{ej,1}, \log_{10} X_{lan,1}) - f(\log_{10} M_{ej,2}, v_{ej,2}, \log_{10} X_{lan,2}). \tag{7}
$$

and

$$
\mu_1 - \mu_2 = m_1 - m_2 + f(\log_{10} M_{ej,1}, \Phi_1, \theta_{\text{obs},1}) - f(\log_{10} M_{ej,2}, \Phi_2, \theta_{\text{obs},2}). \quad (8)
$$

Combined with either a SBF-based measurement of host galaxy NGC 4993 33 or the GW-derived posteriors to anchor the distance distribution for GW170817, this equation is used to measure the sGRB distance moduli (see Fig. 6 for GRB 060614 as an example). We note that the overall luminosity of the light curve is allowed to vary arbitrarily, which adds a linear offset. But, this offset does not affect the color evolution, which is mostly determining the intrinsic parameter distributions. Therefore, we are able to extract both, the intrinsic parameters and the distance to the source, which is a crucial prerequisite for our study.

<figure><img src="image_12.png" /></figure>
<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /></figure>

**Fig. 7** Posterior distributions for $H_0$ for individual events without GW170817. We show them for GRB 060614, GRB 150101B, GRB 160821B, GRB 050709, and their combined posteriors are also shown. **Fig. a** is the Kasen$^{24}$ model with 0.1 mag errors, **b** is the Bulla$^{25}$ model with 0.1 mag errors, **c** is the Kasen model with 0.25 mag errors, and **d** is the Bulla model with 0.25 mag errors. The 1- and 2-$\sigma$ regions determined by the superluminal motion measurement from the radio counterpart (blue)$^{16}$, Planck CMB (TT,TE,EE+lowP+lensing) (green)$^{34}$ and SHoES Cepheid-SN distance ladder surveys (orange)$^{5}$ are also depicted as vertical bands.

**$H_0$ analysis without GW170817.** In the main text, we included GW170817 when combining the posterior distributions for the Hubble constant. Although the combined analysis is more constraining, the inclusion of GW170817 increases the systematic uncertainties as its $H_0$ measurement depends on the peculiar velocity of the host. This problem will remain for all close-by kilonovae. Due to their distance, the sGRB analysis is not affected nearly as much by the peculiar motion of local galaxies. Removing GW170817's Hubble Constant constraints yields measurements of $H_0 = 71.9^{+2.7}_{-2.7}$ km s$^{-1}$ Mpc$^{-1}$ for the Kasen model and $H_0 = 68.2^{+4.3}_{-4.6}$ km s$^{-1}$ Mpc$^{-1}$ for the Bulla model (see Fig. 7); these measurements both have larger error bars than the analysis with GW170817, while still being consistent with it.

We discussed our prior choices in Section 5. In particular, we choose large prior boundaries that allows to describe both, kilonova produced by BNSs and by BHNS systems, and in particular for ejecta mass, $-3 \le \log_{10}(M_{ej}/M_\odot) \le -1$. We can compare this choice to a few other distributions. For example, taking the fit of ref. $^{35}$ and applying them to the posteriors of GW170817$^{6}$ and GW190425$^{64}$, we see a bi-modal distribution, the former peaked near to $\log_{10}(M_{ej}/M_\odot) = -1.3$ and the latter near to $\log_{10}(M_{ej}/M_\odot) = -2.2$; cf. Fig. 8. We also include distributions flat in the component masses with 100 nonparametric EOSs consistent with GW170817 as provided in ref. $^{65}$. This analysis shows that our flat priors cover both of the confirmed events so far, and in this sense, applying flat priors has the significant advantage that a larger range of the parameter space is covered. Further, it reduces systematic uncertainties since no numerical-relativity inferred relations have been applied.

**Choice of ejecta mass priors.** When performing the fits, we rely on the assumption that the standardization is equally valid across the parameter space. Until there are further unambiguous kilonova detections, in particular those that also have GW counterparts, it will be difficult to create more strongly informed

<figure><img src="image_17.png" /></figure>

**Fig. 8** Potential prior distributions for $M_{ej}$. In addition to the flat prior used in this analysis, we show ejecta mass distributions for GW170817$^6$ and GW190425$^{64}$ based on the fit of ref. $^{35}$. We also include distributions flat in the component masses with 100 nonparametric EOSs consistent with GW170817 from ref. $^{65}$.
<figure><img src="image_18.png" /><figcaption>**Fig. 9** Fundamental plane plot. The top panel shows the fit of Eq. (3) for the models in Bulla²⁵ model. The Broad Priors values are derived from the entire available grid, covering a range $-6 \le \log_{10}(M_{ej}/M_\odot) \le -0.15^\circ \le \Phi \le 75^\circ$, and $0^\circ \le \theta_{\text{obs}} \le 90^\circ$. The Realistic Priors values are derived from a fraction of the grid drawn from the sGRB informed priors, covering a range $-3 \le \log_{10}(M_{ej}/M_\odot) \le -1.15^\circ \le \Phi \le 30^\circ$, and $0^\circ \le \theta_{\text{obs}} \le 15^\circ$. The bottom panel shows the difference between the computed values from the model and the fits for the simulations analyzed here. The $1-\sigma$ error bars correspond to those assigned by the Gaussian Process Regression.</figcaption></figure>

priors; for now, doing so would require assumptions on both the distribution of BNS and BHNS masses and an estimate of the NS EOS. It is not so obvious what unphysical regions there are yet, outside of likely very high ejecta masses, which are already excluded by our priors. To evaluate the effect this may have, Figure 9 shows the performance of the fit of Eqs. (3) compared to the Bulla²⁵ model for both Broad and Realistic priors, i.e., those used in this analysis. Here, the Broad prior values are derived from the entire available grid, covering a range $-6 \le \log_{10}(M_{ej}/M_\odot) \le -0.15^\circ \le \Phi \le 75^\circ$, and $0 \le \theta_{\text{obs}} \le 90^\circ$. The Realistic prior values are derived from the prior range used in the sampling, $-3 \le \log_{10}(M_{ej}/M_\odot) \le -1.15^\circ \le \Phi \le 30^\circ$, and $0^\circ \le \theta_{\text{obs}} \le 15^\circ$, which was tuned for the sGRBs. In general, the estimated values from both versions of the fit are consistent with one another and the measured values from the model, and therefore, the effect on the standardization from this perspective is minimal. In the future, if there are regions which are deemed to be disfavored, the priors should be updated to reflect this.

As for assumptions about the ejecta geometry, the current version is generic enough such that it is applicable to both BNS and BHNS cases. However, an improved geometry would likely need to be split into a BNS and BHNS case, which would mean that the standardization for BNS systems could be different from the one for BHNS systems. However, it is currently difficult to assess what regions of the parameter space are favoured for each system and whether these are distinct or overlap. The detection of more kilonovae in the future will help pin down the ejecta geometry and ejecta mass ratio for BNS and BHNS, allowing us to update our priors on $\phi$ and $M_{ej}$ and investigate the possibility of different standardizations.

## Data availability

Upon request, the first author will provide posterior samples from these analyses. All photometric data used in this analysis are publically available from a variety of sources, specified in Methods Section 6 and compiled in https://github.com/mcoughlin/gwemlightcurves/tree/master/lightcurves. Spectral energy distributions for the grid used here will be made available at https://github.com/mbulla/kilonova Models.

## Code availability

The light-curve fitting and Hubble Constant code is available at: https://github.com/ mcoughlin/gwemlightcurves.

Received: 31 January 2020; Accepted: 29 July 2020; Published online: 17 August 2020

## References

1. Riess, A. G. et al. Observational evidence from supernovae for an accelerating universe and a cosmological constant. *Astron. J.* **116**, 1009–1038 (1998).

2. Perlmutter, S. et al. Measurements of $\Omega$ and $\Lambda$ from 42 high-redshift Supernovae. *Astrophys. J.* **517**, 565–586 (1999).

3. Planck Collaboration et al. Planck 2015 results. XIII. Cosmological parameters. *Astron. Astrophys.* **594**, A13 (2016).

4. Riess, A. G. et al. A 2.4% determination of the local value of the Hubble constant. *Astrophys. J.* **826**, 56 (2016).

5. Riess, A. G., Casertano, S., Yuan, W., Macri, L. M. & Scolnic, D. Large magellanic cloud Cepheid standards provide a 1% foundation for the determination of the Hubble constant and stronger evidence for physics beyond $\Lambda$CDM. *Astrophys. J.* **876**, 85 (2019).

6. Abbott, B. P. et al. GW170817: observation of gravitational waves from a binary neutron star inspiral. *Phys. Rev. Lett.* **119**, 161101 (2017).

7. Abbott, B. P. et al. Gravitational waves and gamma-rays from a binary neutron star merger: GW170817 and GRB 170817A. *Astrophys. J. Lett.* **848**, L13 (2017).

8. Goldstein, A. et al. An ordinary short gamma-ray burst with extraordinary implications: Fermi-GBM detection of GRB 170817a. *Astrophys. J.* **848**, L14 (2017).

9. Coulter, D. A. et al. Swope Supernova Survey 2017a (SSS17a), the optical counterpart to a gravitational wave source. *Science* **358**, 1556–1558 (2017).

10. Coughlin, M. W. et al. Constraints on the neutron star equation of state from AT2017gfo using radiative transfer simulations. *Mon. Not. R. Astron. Soc.* **480**, 3871–3878 (2018).

11. Abbott, B. P. et al. Estimating the contribution of dynamical ejecta in the kilonova associated with GW170817. *Astrophys. J. Lett.* **850**, L39 (2017).

12. Abbott, B. P. et al. A gravitational-wave standard siren measurement of the Hubble constant. *Nature* **551**, 85–88 (2017).

13. Schutz, B. F. Determining the Hubble constant from gravitational wave observations. *Nature* **323**, 310 (1986).

14. Doctor, Z. Thunder and lightning: using neutron-star mergers as simultaneous standard candles and siren to measure cosmological parameters. *Astrophys. J.* **892**, L16 (2020).

15. Abbott, B. P. et al. A gravitational-wave measurement of the Hubble constant following the second observing run of Advanced LIGO and Virgo. Preprint at https://arxiv.org/abs/1908.06060 (2019).

16. Hotokezaka, K. et al. A hubble constant measurement from superluminal motion of the jet in GW170817. *Nat. Astron.* https://www.nature.com/articles/s41550-019-0820-1?error=cookies_not_supported&code=e98ffae6-00e6-4882-a9a6-3443f30d831b. (2019).

17. Dietrich, T. et al. New constraints on the supranuclear equation of state and the Hubble constant from nuclear physics—multi-messenger astronomy. Preprint at https://arxiv.org/abs/2002.11355 (2020).

18. Coughlin, M. W. et al. Standardizing kilonovae and their use as standard candles to measure the Hubble constant. *Phys. Rev. Res.* **2**, 022006 (2020).

19. Kashyap, R., Raman, G. & Ajith, P. Can kilonova light curves be standardized? *Astrophys. J.* **886**, L19 (2019).

20. Coughlin, M. W. et al. GROWTH on S190425z: searching thousands of square degrees to identify an optical or infrared counterpart to a binary neutron star merger with the zwicky transient facility and palomar gattini-IR. *Astrophys. J.* **885**, L19 (2019).

21. Lattimer, J. M. & Schramm, D. N. Black-hole-neutron-star collisions. *Astrophys. J. Lett.* **192**, L145–L147 (1974).

22. Li, L.-X. & Paczynski, B. Transient events from neutron star mergers. *Astrophys. J. Lett.* **507**, L59 (1998).

23. Metzger, B. D. Kilonovae. *Living Rev. Rel.* **20**, 3 (2017).

24. Kasen, D., Metzger, B., Barnes, J., Quataert, E. & Ramirez-Ruiz, E. Origin of the heavy elements in binary neutron-star mergers from a gravitational-wave event. *Nature* **551**, 80 EP – (2017).

25. Bulla, M. POSSIS: predicting spectra, light curves, and polarization for multidimensional models of supernovae and kilonovae. *Mon. Not. R. Astron. Soc.* **489**, 5037–5045 (2019).

26. Ascenzi, S. et al. A luminosity distribution for kilonovae based on short gamma-ray burst afterglows. *Mon. Not. R. Astron. Soc.* **486**, 672–690 (2019).

27. Fong, W. et al. The afterglow and early-type host galaxy of the short GRB 150101b at $z = 0.1343$. *Astrophys. J.* **833**, L15 (2016).

28. Fox, D. B. et al. The afterglow of GRB 050709 and the nature of the short-hard $\gamma$-ray bursts. *Nature* **437**, 845 EP – (2005).

29. Kasliwal, M. M., Korobkin, O., Lau, R. M., Wollaeger, R. & Fryer, C. L. Infrared emission from Kilonovae: the case of the nearby short hard burst GRB 160821B. *Astrophys. J. Lett.* **843**, L34 (2017).

30. Zhang, B. et al. Making a short gamma-ray burst from a long one: implications for the nature of GRB 060614. *Astrophys. J. Lett.* **655**, L25 (2007).

31. Coughlin, M. W. et al. Implications of the search for optical counterparts during the first six months of the Advanced LIGO's and Advanced Virgo's third observing run: possible limits on the ejecta mass and binary properties. *Mon. Not. R. Astron. Soc.* https://academic.oup.com/mnras/article-abstract/492/1/863/5672639?redirectedFrom=fulltext (2019).

32. Hotokezaka, K. & Nakar, E. Radioactive heating rate of r-process elements and macronova light curve. *Astrophys. J.* **891**, L12 (2020).
33. Cantiello, M. et al. A precise distance to the host galaxy of the binary neutron star merger GW170817 using surface brightness fluctuations. *Astrophys. J.* **854**, L31 (2018).

34. Ade, P. A. R. et al. Planck 2015 results. XIII. Cosmological parameters. *Astron. Astrophys.* **594**, A13 (2016).

35. Coughlin, M. W. et al. Multimessenger Bayesian parameter inference of a binary neutron star merger. *Mon. Not. R. Astron. Soc. Lett.* **489**, L91–L96 (2019).

36. Abbott, B. P. et al. Properties of the binary neutron star merger GW170817. *Phys. Rev. X* **9**, 011001 (2019).

37. Meegan, C. et al. The fermi gamma-ray burst monitor. *Astrophys. J.* **702**, 791–804 (2009).

38. Coughlin, M. W. et al. 2900 square degree search for the optical counterpart of short gamma-ray burst GRB 180523b with the Zwicky Transient Facility. *Publ. Astron. Soc. Pac.* **131**, 048001 (2019).

39. Doctor, Z., Farr, B., Holz, D. E. & Pürrer, M. Statistical gravitational waveform models: what to simulate next? *Phys. Rev. D* **96**, 123011 (2017).

40. Troja, E. et al. A luminous blue kilonova and an off-axis jet from a compact binary merger at z=0.1341. *Nat. Commun.* **9**, 4089 (2018).

41. Fong, W. et al. Short GRB 130603B: discovery of a jet break in the optical and radio afterglows, and a mysterious late-time X-ray excess. *Astrophys. J.* **780**, 118 (2014).

42. D'Avanzo, P. Short gamma-ray bursts: a review. *J. High. Energy Phys.* **7**, 73–80 (2015).

43. Saleem, M. Prospects of joint detections of neutron star mergers and short GRBs with Gaussian structured jets. *Mon. Not. R. Astron. Soc.* **493**, 1633–1639 (2020).

44. Beniamini, P., Petropoulou, M., Barniol Duran, R. & Giannios, D. A lesson from GW170817: most neutron star mergers result in tightly collimated successful GRB jets. *Mon. Not. R. Astron. Soc.* **483**, 840–851 (2018).

45. Barnes, J., Kasen, D., Wu, M.-R. & Martínez-Pinedo, G. Radioactivity and thermalization in the ejecta of compact object mergers and their impact on kilonova light curves. *Astrophys. J.* **829**, 110 (2016).

46. Kawaguchi, K., Shibata, M. & Tanaka, M. Radiative transfer simulation for the optical and near-infrared electromagnetic counterparts to GW170817. *Astrophys. J. Lett.* **865**, L21 (2018).

47. Radice, D. et al. Binary neutron star mergers: mass ejection, electromagnetic counterparts and nucleosynthesis. *Astrophys. J.* **869**, 130 (2018).

48. Hjorth, J. et al. The optical afterglow of the short γ-ray burst GRB 050709. *Nature* **437**, 859–861 (2005).

49. Covino, S. et al. Optical emission from GRB 050709: a short/hard GRB in a star-forming galaxy. *Astron. Astrophys.* **447**, L5–L8 (2006).

50. Jin, Z.-P. et al. The Macronova in GRB 050709 and the GRB-macronova connection. *Nat. Commun.* **7**, 12898 (2016).

51. Jin, Z.-P. et al. Short GRBs: opening angles, local neutron star merger rate, and off-axis events for GRB/GW Association. *Astrophys. J.* **857**, 128 (2018).

52. Jin, Z.-P. et al. The light curve of the macronova associated with the long-short burst GRB 060614. *Astrophys. J. Lett.* **811**, L22 (2015).

53. Yang, B. et al. A possible macronova in the late afterglow of the long-short burst GRB 060614. *Nat. Commun.* **6**, 7323 (2015).

54. Lamb, G. P. et al. Short GRB 160821b: a reverse shock, a refreshed shock, and a well-sampled kilonova. *Astrophys. J.* **883**, 48 (2019).

55. Rossi, A. et al. A comparison between short GRB afterglows and kilonova AT2017gfo: shedding light on kilonovae properties. *Mon. Not. Roy. Astron. Soc.* **493**, 3379–3397 (2020).

56. Lu, Y., Huang, Y. F. & Zhang, S. N. A tidal disruption model for the gamma-ray burst of GRB 060614. *Astrophys. J.* **684**, 1330–1335 (2008).

57. Caito, L. et al. GRB060614: a “fake” short GRB from a merging binary system. *Astron. Astrophys.* **498**, 501–507 (2009).

58. Xu, D. et al. In search of progenitors for supernovaless gamma-ray bursts 060505 and 060614: re-examination of their afterglows. *Astrophys. J.* **696**, 971–979 (2009).

59. Metzger, B. D. Kilonovae. *Living Rev. Relativ.* **20**, 3 (2017).

60. Arnett, W. D. Type I supernovae. I—ānalytic solutions for the early part of the light curve. *Astrophys. J.* **253**, 785–797 (1982).

61. Smartt, S. J. et al. A kilonova as the electromagnetic counterpart to a gravitational-wave source. *Nature* **551**, 75 EP – (2017).

62. Hjorth, J. et al. The distance to NGC 4993: the host galaxy of the gravitational-wave event GW170817. *Astrophys. J.* **848**, L31 (2017).

63. Levan, A. J. et al. The environment of the binary neutron star merger GW170817. *Astrophys. J.* **848**, L28 (2017).

64. Abbott, B. P. et al. GW190425: observation of a compact binary coalescence with total mass ~ 3.4m☉. *arXiv* https://iopscience.iop.org/article/10.3847/2041-8213/ab75f5 (2020).

65. Landry, P. & Essick, R. Nonparametric inference of the neutron star equation of state from gravitational wave observations. *Phys. Rev. D* **99**, 084049 (2019).

## Acknowledgements

M.W.C. acknowledges support from the National Science Foundation with grant number PHY-2010970. T.D. acknowledges support by the European Union's Horizon 2020 research and innovation program under grant agreement No 749145, BNSmergers. N.C. and J.H. acknowledge support from the National Science Foundation with grant number PHY-1806990. S.A. is supported by the CNES Postdoctoral Fellowship at Laboratoire Astroparticle et Cosmologie. N.C., M.C., and J.H. gratefully acknowledge support from the Observatoire Côte d'Azur, including hospitality for M.C. and J.H. in Summer 2019. The UCSC team is supported in part by NASA grant NNG17PX03C, NSF grant AST-1911206, the Gordon & Betty Moore Foundation, the Heising-Simons Foundation, and by fellowships from the David and Lucile Packard Foundation to R.J.F. D.A.C. acknowledges support from the National Science Foundation Graduate Research Fellowship under Grant DGE1339067.

## Author contributions

M.W.C. conducted the light-curve analysis and was the primary author of the manuscript. M.W.C., S.A., T.D., R.J.F., J.H., M.B., N.C., D.C., L.I., and N.A. contributed to the data analysis procedures and edits to the manuscript.

## Competing interests

The authors declare no competing interests.

## Additional information

Supplementary information is available for this paper at https://doi.org/10.1038/s41467-020-17998-5.

Correspondence and requests for materials should be addressed to M.W.C.

Peer review information *Nature Communications* thanks the anonymous reviewers for their contribution to the peer review of this work. Peer reviewer reports are available.

Reprints and permission information is available at http://www.nature.com/reprints

Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_19.png" /></figure>

Open Access This article is licensed under a Creative Commons

Attribution 4.0 International License, which permits use, sharing,

adaptation, distribution and reproduction in any medium or format, as long as you give
appropriate credit to the original author(s) and the source, provide a link to the Creative
Commons license, and indicate if changes were made. The images or other third party
material in this article are included in the article's Creative Commons license, unless
indicated otherwise in a credit line to the material. If material is not included in the
article's Creative Commons license and your intended use is not permitted by statutory
regulation or exceeds the permitted use, you will need to obtain permission directly from
the copyright holder. To view a copy of this license, visit http://creativecommons.org/
licenses/by/4.0/.

© The Author(s) 2020