We present uniform modeling of eight kilonovae five following short gammaray bursts GRBs including GRB170817A and three following long GRBs. We model their broadband afterglows to determine the relative contributions of afterglow and kilonova emission. We fit the kilonovae using a threecomponent model in MOSFiT that accounts for ejecta geometry and find population median ejecta masses for the total blue kappaB  0.5 cm2  g purple kappaP  3 cm2  g and red kappaR  10 cm2  g components of Mej tot  0.0850.0400.110 Modot Mej B  0.0060.0040.015 Modot Mej P  0.0200.0100.034 Modot and Mej R  0.0510.0450.100 Modot 68 confidence. The kilonova of GW170817 is near the median of the sample in most derived properties while the sample indicates great diversity. We investigate trends between the ejecta masses and the isotropicequivalent and beamingcorrected gammaray energies Egamma iso Egamma as well as restframe durations T90 rest. We find long GRB kilonovae have higher median red ejecta masses Mej R gt 0.05 Modot compared to onaxis short GRB kilonovae Mej R lt 0.02 Modot. We also observe a weak scaling between the total and red ejecta masses with Egamma iso and Egamma though a larger sample is needed to establish a significant correlation. These findings imply a connection between mergerdriven long GRBs and larger tidal dynamical ejecta masses which may indicate that their progenitors are asymmetric compact object binaries. We produce representative kilonova light curves and find that the planned depths and cadences of the Rubin and Roman Observatory surveys will be sufficient for orderofmagnitude constraints on Mej B and for Roman Mej P and Mej R of future kilonovae at z lt 0.1.
['Astrophysics - High Energy Astrophysical Phenomena']
Uniform Modeling of Observed Kilonovae Implications for Diversity and the Progenitors of MergerDriven Long GammaRay Bursts
{'Uniform Modeling of Observed Kilonovae: Implications for Diversity and the Progenitors of Merger-Driven Long Gamma-Ray Bursts': "J. C. Rastinejad, 1 W. Fong, 1 C. D. Kilpatrick, 1 M. Nicholl, 2 and B. D. Metzger 3, 4  \n1 Center for Interdisciplinary Exploration and Research in Astrophysics (CIERA) and Department of Physics and Astronomy, Northwestern University, Evanston, IL 60208, USA  \n- 2 Astrophysics Research Centre, School of Mathematics and Physics, Queen's University Belfast, BT7 1NN, UK  \n3 Department of Physics and Columbia Astrophysics Laboratory, Columbia University, Pupin Hall, New York, NY 10027, USA 4 Center for Computational Astrophysics, Flatiron Institute, 162 5th Ave, New York, NY 10010, USA", 'ABSTRACT': 'We present uniform modeling of eight kilonovae, five following short gamma-ray bursts (GRBs; including GRB170817A) and three following long GRBs. We model their broadband afterglows to determine the relative contributions of afterglow and kilonova emission. We fit the kilonovae using a three-component model in MOSFiT , and report population median ejecta masses for the total, blue ( κ B = 0 . 5 cm 2 g -1 ), purple ( κ P = 3 cm 2 g -1 ), and red ( κ R = 10 cm 2 g -1 ) components. The kilonova of GW170817 is near the sample median in most derived properties. We investigate trends between the ejecta masses and the isotropic-equivalent and beaming-corrected γ -ray energies ( E γ, iso , E γ ), as well as rest-frame durations ( T 90 , rest ). We find long GRB kilonovae have higher median red ejecta masses ( M ej , R ≳ 0 . 05 M ⊙ ) compared to on-axis short GRB kilonovae ( M ej , R ≲ 0 . 02 M ⊙ ). We also observe a weak scaling between the total and red ejecta masses with E γ, iso and E γ , though a larger sample is needed to establish a significant correlation. These findings imply a connection between merger-driven long GRBs and larger tidal dynamical ejecta masses, which may indicate that their progenitors are asymmetric compact object binaries. We produce representative kilonova light curves and find that the planned depths and cadences of the Rubin and Roman Observatory surveys will be sufficient for order-of-magnitude constraints on M ej , B (and, for Roman, M ej , P and M ej , R ) of future kilonovae at z ≲ 0 . 1.  \nKeywords: kilonovae, r -process, gamma-ray bursts', '1. INTRODUCTION': "On August 17, 2017, the nearly coincident detection of a binary neutron star (BNS) merger through gravitational waves (GWs, GW170817; Abbott et al. 2017a,b) and a short gamma-ray burst (GRB 170817A; Goldstein et al. 2017; Savchenko et al. 2017) confirmed the long-theorized connection between these 'multimessenger' signals. Adding to this long-awaited breakthrough, an optical-near-IR counterpart ('AT 2017gfo') to the GW and GRB signals was observed (Arcavi et al. 2017; Coulter et al. 2017; Lipunov et al. 2017; Tanvir et al. 2017; Soares-Santos et al. 2017; Valenti et al. 2017) and strongly resembled theoretical predic-  \ntions for a kilonova, a thermal transient powered by the radioactive decay of elements beyond the iron-peak formed through rapid neutron capture nucleosynthesis (' r -process'; e.g., Burbidge et al. 1957; Cameron 1957; Lattimer & Schramm 1974; Li & Paczy'nski 1998; Metzger et al. 2010). Indeed, the counterpart's rapid temporal evolution pointed to a relatively low ejecta mass ( ≲ 0 . 1 M ⊙ ), expected for neutron star mergers, while its reddening on ∼ day timescales indicated a high (compared to that of iron-group elements) ejecta opacity produced by newly-created heavy elements (e.g., Barnes & Kasen 2013; Tanaka & Hotokezaka 2013; Kasen et al. 2017). Modeling of AT 2017gfo demonstrated that the light curve was well-explained with a two-component model, in which each component is parameterized by a different ejecta mass, velocity, and opacity (e.g., Cowperthwaite et al. 2017; Drout et al. 2017; Kasen et al.  \n2017; Kilpatrick et al. 2017; Tanvir et al. 2017; Troja et al. 2017; Villar et al. 2017b), providing evidence for multiple emission sources from the merger.  \nThe first definitive kilonova opened a host of new questions and predictions that may only be studied with a wider population of observed events. Rapid progress on the theoretical end of kilonova studies has resulted in models that incorporate additional physics, including updated neutrino schemes and viewing angle dependencies (e.g., Bulla 2019; Nicholl et al. 2021; Wollaeger et al. 2021; Chase et al. 2022; Curtis et al. 2023), and a mapping of progenitors, remnants, and novel emission mechanisms to predicted light curve properties (e.g., Arcavi 2018; Gottlieb et al. 2018; Metzger et al. 2018; Kawaguchi et al. 2020a; Gompertz et al. 2023). In parallel, GW observations of diverse BNS and NSBH progenitors (e.g., Abbott et al. 2020a,b; The LIGO Scientific Collaboration et al. 2024) and recent kilonovae discoveries following LGRBs (Rastinejad et al. 2022a; Troja et al. 2022; Yang et al. 2022; Gillanders et al. 2023; Levan et al. 2024; Yang et al. 2024) further motivate expectations for observed kilonova diversity. Looking forward, constraints on light curve diversity are critical for kilonova search strategies with next-generation facilities such as the Rubin Observatory and the Roman Space Telescope (Roman; Margutti et al. 2018; Andreoni et al. 2024).  \nDespite progress on the theoretical end, constraints on kilonova ejecta parameters from observations other than AT2017gfo remain limited due to the low rates of detected neutron star mergers (e.g., Mandel & Broekgaarden 2022; The LIGO Scientific Collaboration et al. 2021). Kilonovae are also relatively faint and fleeting transients, rendering them difficult to observe at distances beyond z ≈ 0 . 1 without rapid response and highly-sensitive telescopes. Several kilonovae following GRBs have been discovered contemporaneously and in archival data (e.g., Berger et al. 2013; Tanvir et al. 2013; Yang et al. 2015; Gompertz et al. 2018) and, when combined with deep upper limits, reveal a span of ≈ 100 in optical luminosity (e.g., Gompertz et al. 2018; Rossi et al. 2020; Rastinejad et al. 2021). However, the majority of previous fits to observed light curves have been with by multiple codes with varying assumptions (e.g., Lamb et al. 2019; Troja et al. 2019; O'Connor et al. 2021). Thus, direct comparisons of derived parameters (e.g., mass, velocity, composition) between individual events are inadvisable.  \nPrevious uniform modeling of GRB kilonovae used a single component model (Ascenzi et al. 2019), prohibiting a search for diversity among emission components. Further, a uniform analysis has not been performed  \nsince 2019 (Ascenzi et al. 2019), thus excluding three recent events (GRBs 200522A, 211211A and 230307A; Fong et al. 2021; O'Connor et al. 2021; Rastinejad et al. 2022a; Troja et al. 2022; Levan et al. 2024; Yang et al. 2024). Surprisingly, two of these kilonovae followed long-duration GRBs. The progenitors and/or emission mechanisms driving their longer γ -ray signals remain an unknown, with several observationally untested theories (e.g., Yang et al. 2022; Gottlieb et al. 2023). Beyond significantly expanding the sample size, these new events motivate an updated, multi-component modeling endeavor to explore kilonova diversity and compare γ -ray and kilonova properties across the short and long GRB populations.  \nHere, we perform uniform, multi-component modeling of a sample of eight kilonovae discovered over 2005-2023. We provide compiled multi-wavelength light curves for future fitting efforts with additional models. In Section 2 we describe our sample selection and provide details for each event. In Section 3 we detail our modeling of the synchrotron afterglow component to extricate the kilonova emission. In Section 4 we describe our kilonova modeling procedure and report our results. In Section 5 we discuss the implications of our results for future work. We assume a cosmology of H 0 = 69.6 km s -1 Mpc -1 , Ω M = 0.286, Ω vac = 0.714 (Bennett et al. 2014) and report magnitudes in the AB system throughout this work.", '2.1. Sample Selection': "We begin with the list of claimed kilonovae compiled in previous literature and covering the full sample of short GRBs discovered in the postSwift era (Gompertz et al. 2018; Rastinejad et al. 2021). These include GRBs050709 (Jin et al. 2016), 050724A (Gao et al. 2017), 060614 (Yang et al. 2015), 070714B (Gao et al. 2017), 070809 (Jin et al. 2020), 130603B (Berger et al. 2013; Tanvir et al. 2013), 150101B (Troja et al. 2018), 160821B (Lamb et al. 2019; Troja et al. 2019), 170817A (GW170817) and 200522A (Fong et al. 2021; O'Connor et al. 2021). Though GRBs080503 and 150424A have been claimed as a kilonova candidate (Perley et al. 2009; Bucciantini et al. 2012; Rossi et al. 2022), we do not include these events as their redshifts, and thus their kilonova luminosities, remain uncertain (Perley et al. 2009; Fong et al. 2022). In addition, we also consider the more recent kilonova candidates GRB 211211A (Rastinejad et al. 2022a; Troja et al. 2022), and GRB 230307A  \n(Levan et al. 2024; Yang et al. 2024), both of which were identified following recent merger-driven long GRBs 1 .  \nWe limit our sample to events with sufficient observations (typically, multiple X-ray detections, several optical-near-IR detections past 2 days and at least one radio observation) for robust afterglow and multicomponent kilonova modeling. To evaluate the availability of X-ray and radio observations, we use the catalog of Fong et al. (2015), supplemented with data from the Gamma-ray Circular Notices (GCNs) and the literature (Fox et al. 2005; Mangano et al. 2007; Londish et al. 2006; Berger et al. 2013; Tanvir et al. 2013; Fong et al. 2014; Lamb et al. 2019; Troja et al. 2019; Mei et al. 2022; Rastinejad et al. 2022a; Troja et al. 2022; Levan et al. 2024; Yang et al. 2024). We restrict our sample to events with two or more optical-near-IR observations past 2 days for two reasons. First, this is the timescale on which we roughly expect the kilonova to significantly contribute to the observed flux. Second, multiple observations allow us to assess fading, which in turn constrains the ejecta mass and velocity, and/or color, giving insight to the relative contribution from different components (e.g., Barnes & Kasen 2013). To evaluate available optical and near-IR data for each burst we consider data collected and published in short GRB kilonova compilations (Gompertz et al. 2018; Rastinejad et al. 2021) for bursts prior to 2021 and the literature for those from 2021-2023. We do not include GRBs 050724A, 070714B and 150101B in our sample as they each have only a single optical detection past 2 days. We further remove GRB070809 from our sample as it was not observed in the radio and no optical data is available past ≈ 1 . 5 days.  \nWe are left with eight events that meet our criteria: GRBs 050709, 060614, 130603B, 160821B, 170817A (GW170817), 200522A, 211211A, and 230307A. The redshift range of these events is relatively low for short GRBs ( z = 0 . 008 -0 . 554; Fong et al. 2022), but greater than the expected horizon for GW-detected compact object mergers in O4 and O5 (Abbott et al. 2016). We describe each burst and the dataset used in our analysis in the following sections. We list the main properties of each burst in Table 1, including the γ -ray durations, GRB fluence ( f γ ), isotropic-equivalent γ -ray energy ( E γ, iso ; calculated using the method of Fong et al. 2015), and galaxy type (Nugent et al. 2022; Levan et al. 2024). For Swift GRBs we collect values for f γ (15-350 keV) and the durations over which 50% and 90% of the  \ngamma-ray fluence was detected (T 50 and T 90 ) from the Swift -BAT catalog 2 (Lien et al. 2016). For GRBs detected with other satellites we collect values from the GCNs and literature (Fox et al. 2005; Stanbro & Meegan 2016; Goldstein et al. 2017; Veres et al. 2023; Dalessi et al. 2023). In the following sub-sections, we summarize the burst discovery for each event in our sample and the source of data used in our analysis.", '2.2. GRB050709': 'GRB050709 was detected by the High Energy Transient Explorer II (HETE-II; Ricker et al. 2003) at 22:36:37 UT on 2005 July 9, with a T 90 = 0 . 07 s (Hjorth et al. 2005). Follow-up observations with the Chandra X-ray Observatory ( Chandra ) revealed a new Xray source within the HETE-II localization (Fox et al. 2005). Follow-up observations by the Swift X-Ray Telescope (XRT; Burrows et al. 2005) and Chandra confirmed this counterpart as fading and thus, likely related to GRB050709 (Fox et al. 2005). At the position of the X-ray source, an optical counterpart was detected, embedded in a star-forming galaxy at z = 0 . 161 (Fox et al. 2005; Hjorth et al. 2005). The Very Large Array (VLA) observed the position of the X-ray counterpart over four epochs, but did not detect any significant emission (Fox et al. 2005).  \nImaging with the Hubble Space Telescope ( HST ) revealed long-lived emission in the F814W band, which has been attributed to a kilonova by numerous groups in the literature (Jin et al. 2016; Gompertz et al. 2018; Ascenzi et al. 2019). We combine the X-ray, radio (Fox et al. 2005) and optical-near-IR (Fox et al. 2005; Hjorth et al. 2005; Covino et al. 2006) datasets for our analysis. We present all observations used in our analysis in Table 5.', '2.3. GRB060614': 'GRB060614 was discovered by the Swift Burst Alert Telescope (BAT; Barthelmy et al. 2005) on 2006 June 14 at 12:43:48.5 UT with a T 90 = 109 ± 3 s (Barthelmy et al. 2006). Prompt follow-up by Swift -XRT and the Swift Ultra-Violet Optical Telescope (UVOT; Roming et al. 2005) revealed bright counterparts to the burst (Mangano et al. 2007). The early ultraviolet through optical afterglow spectral energy distribution (SED) provides evidence for low line-of-sight extinction and a z < 1 . 3 origin (Gehrels et al. 2006). Subsequent followup by ground-based observatories revealed an optical  \nTable 1. GRB Sample & Properties  \nNote -1 Milky Way extinction values taken from Schlafly & Finkbeiner (2011).  \ncounterpart with a spectroscopic redshift of z = 0 . 125 (Price et al. 2006; Fugazza et al. 2006) on the outskirts of a star-forming galaxy at the same distance. The optical counterpart was monitored to late times with HST (Yang et al. 2015).  \nNumerous deep imaging observations and spectra, extending out to ≈ 65 days, failed to reveal the expected counterpart to a long-duration GRB, a SN Ic-BL (Della Valle et al. 2006; Fynbo et al. 2006; Gal-Yam et al. 2006). This fact, combined with its intriguing gamma-ray properties, motivated a later analysis showing that the optical light curve reddened at later times, suggesting instead the presence of a kilonova (Yang et al. 2015).  \nFor our analysis, we use the XRT and UVOT light curves (though we do not use UVOT detections past 10 days as the exposures are on ≳ day timescales; Mangano et al. 2007), ATCA radio upper limits (Londish et al. 2006) and combine the optical datasets in the literature (Cobb 2006; Della Valle et al. 2006; Fynbo et al. 2006; Schmidt et al. 2006; Xu et al. 2009; Yang et al. 2015).', '2.4. GRB130603B': 'GRB130603B was detected by Swift -BAT and the KonusWind Observatory (Melandri et al. 2013; Golenetskii et al. 2013) on 2013 June 13 at 15:49:14 UT with T 90 = 0 . 18 ± 0 . 02 s (Berger et al. 2013). Prompt Swift -XRT observations revealed an X-ray counterpart (Melandri et al. 2013). Rapid optical follow-up revealed the optical afterglow within the X-ray localiza-  \nion (Levan et al. 2013). Spectroscopy of the afterglow identified a GRB origin of z = 0 . 356 (Foley et al. 2013; Thone et al. 2013). Additional multi-wavelength followup detected a radio counterpart and a well-sampled optical and X-ray afterglow (e.g., Fong et al. 2014). Later observations with HST revealed that the optical counterpart had significantly reddened, resulting in a bright near-IR detection and deep limits in the optical bands. This represented the first bona-fide claim of an r -process-enriched kilonova (Berger et al. 2013; Tanvir et al. 2013).  \nThe kilonova detection of GRB130603B has been modeled by numerous groups in the literature, though with only a single detection, precise ejecta parameters remain uncertain (e.g., Berger et al. 2013; Tanvir et al. 2013; Hotokezaka et al. 2013; Barnes et al. 2016). We employ the multi-wavelength dataset compiled in Fong et al. (2014) and combine the optical through near-IR datasets published in the literature (Berger et al. 2013; Cucchiara et al. 2013; de Ugarte Postigo et al. 2014; Tanvir et al. 2013; Fong et al. 2014).', '2.5. GRB160821B': 'GRB160821B was detected by the Fermi Space Telescope Gamma-ray Burst Monitor (GBM; Meegan et al. 2009) and the Swift -BAT on 2016 August 21 at 22:29:13 UT. It was classified as a short GRB with T 90 = 0 . 048 ± 0 . 07 s ( Swift ; Lamb et al. 2019) and T 90 ≈ 1 s ( Fermi ; Stanbro & Meegan 2016). Swift -XRT quickly localized  \na counterpart (Evans et al. 2016). UVOT observed the location but did not detect a counterpart (Breeveld & Siegel 2016). Prompt follow-up identified optical and radio counterparts on the outskirts of a bright galaxy at z = 0 . 1616 (Fong et al. 2016; Levan et al. 2016; Xu et al. 2016), motivating further multi-color follow-up by HST and ground observatories (e.g., Kasliwal et al. 2017; Lamb et al. 2019; Troja et al. 2019).  \nThe afterglow and kilonova of GRB160821B have been analyzed by multiple groups in the literature (e.g., Kasliwal et al. 2017; Lamb et al. 2019; Troja et al. 2019), and show evidence for an early reverse shock (e.g., Lamb et al. 2019) and a bluer kilonova relative to AT 2017gfo (e.g., Rastinejad et al. 2021). For our analysis, we collect Swift and XMM-Newton (Troja et al. 2019) and VLA radio observations (Lamb et al. 2019; Fong et al. 2021) from the literature. We combine the optical-NIR datasets (Kasliwal et al. 2017; Lamb et al. 2019; Troja et al. 2019) and report where we draw specific measurements from in Table 5.', '2.6. GRB170817A/GW170817': "GRB170817A was detected by the Fermi -GBM and INTernational Gamma-ray Astrophysics Laboratory (INTEGRAL; Winkler et al. 2003) on 2017 August 17 at 12:41:06 UT, 1.7 s after the LVK-detected BNS merger GW170817 (Abbott et al. 2017b,a; Goldstein et al. 2017; Savchenko et al. 2017). The gamma-ray duration was T 90 = 2 . 0 ± 0 . 5 s (50-300 keV; Fermi ; Goldstein et al. 2017). An optical counterpart (AT 2017gfo) was first found δt = 0 . 452 days following the event discovery (Arcavi et al. 2017; Coulter et al. 2017; Lipunov et al. 2017; Tanvir et al. 2017; Soares-Santos et al. 2017; Valenti et al. 2017). The community observed the transient in exquisite temporal and color detail, revealing a fast-fading, reddening source with spectroscopic evidence of the radioactive decay of lanthanide elements (e.g., Chornock et al. 2017; Gillanders et al. 2022; Kasen et al. 2017; Kilpatrick et al. 2017; Nicholl et al. 2017; Pian et al. 2017; Smartt et al. 2017; Watson et al. 2019; Hotokezaka et al. 2023).  \nWe utilize the optical-NIR dataset compiled in Villar et al. 2017b with observations from the literature (Andreoni et al. 2017; Arcavi et al. 2017; Coulter et al. 2017; Cowperthwaite et al. 2017; D'ıaz et al. 2017; Drout et al. 2017; Evans et al. 2017; Hu et al. 2017; Kasliwal et al. 2017; Lipunov et al. 2017; Pian et al. 2017; Pozanenko et al. 2018; Shappee et al. 2017; Smartt et al. 2017; Tanvir et al. 2017; Troja et al. 2017; Utsumi et al. 2017; Valenti et al. 2017) extending out to ≈ 25 days. Specifically, we employ the set of observations used in their analysis, as there are known inconsistencies amongst  \nthe full dataset (Villar et al. 2017b). We do not collect multi-wavelength data to model this component as it is well-established that the off-axis afterglow did not contaminate observations on the timescales of the kilonova ( δt ≲ 30 days; e.g., Lyman et al. 2018; Fong et al. 2019; Margutti & Chornock 2020; Kilpatrick et al. 2021).", '2.7. GRB200522A': "GRB200522A was discovered on 2020 May 22 at 11:41:34 UT by Swift -BAT with T 90 = 0 . 62 ± 0 . 08 s (Evans et al. 2020). A prompt XRT counterpart location was reported (Evans et al. 2020), within which a catalogued galaxy with photometric redshift z phot = 0 . 4 ± 0 . 1 was noted (Alam et al. 2015; Fong et al. 2020). Subsequent observations identified a radio counterpart (Schroeder et al. 2020; Fong et al. 2021) and secured the host redshift to z = 0 . 554 (Fong et al. 2021; O'Connor et al. 2021). Follow-up by HST and ground observatories uncovered a fading optical-near-IR counterpart embedded in the host galaxy (Strausbaugh & Cucchiara 2020; Fong et al. 2021; O'Connor et al. 2021).  \nWe draw the Swift -XRT observations from the United Kingdom Swift Science Data Centre (UKSSDC; Evans et al. 2007, 2009) and incorporate Chandra and VLA observations (Fong et al. 2021). We combine the optical-near-IR data in the literature (Fong et al. 2021; O'Connor et al. 2021).", '2.8. GRB211211A': 'GRB211211A was detected on 11 December 2021 at 13:09:59 UT by Swift -BAT, Fermi -GBM, INTEGRAL, and the CALET Gamma-ray Burst Monitor (Mangan et al. 2021; Stamatikos et al. 2021; Tamura et al. 2021; Minaev et al. 2021). It was reported as a bright, longduration GRB with T 90 = 50 . 7 ± 0 . 9 s ( Swift ) and T 90 = 34 . 3 ± 0 . 6 s ( Fermi ; Veres et al. 2023). Swift -XRT and UVOT promptly identified counterparts to the burst (Osborne et al. 2021; Belles et al. 2021). Notably, the detection of the afterglow in the UVW 2 filter limits the event origin to z < 1 . 4 (Rastinejad et al. 2022a). An early optical counterpart was discovered proximate to galaxy SDSS J140910.47+275320.8 (Zheng et al. 2021). Later spectroscopy revealed a featureless afterglow (Rastinejad et al. 2022a) and a galaxy redshift of z = 0 . 0763, rendering it one of the most nearby GRBs observed across all durations (Rastinejad et al. 2022a; Troja et al. 2022).  \nMotivated by the low redshift of the putative host galaxy, the counterpart was followed in the optical and near-IR. These observations revealed a fast-fading, red transient with similar luminosities and behavior to AT2017gfo (Rastinejad et al. 2022a). Later ( δt ≈ 2-3  \nweeks), deep optical upper limits revealed no sign of a supernova counterpart to a luminosity lower than that of any known GRB-SN (Rastinejad et al. 2022a; Troja et al. 2022), further motivating an interpretation of the red excess as a kilonova. The kilonova of GRB 211211A has been modeled in the literature by numerous groups (e.g., Rastinejad et al. 2022a; Troja et al. 2022; Mei et al. 2022; Yang et al. 2022; Kunert et al. 2024). For our analysis, we use Swift -XRT observations from UKSSDC, XMM-Newton upper limits (Mei et al. 2022) and a VLA radio observation (Rastinejad et al. 2022a). We use the Swift -UVOT and optical-near-IR datasets from Rastinejad et al. 2022a, and add early optical data from Troja et al. 2022.', '2.9. GRB230307A': "GRB230307A was detected on 7 March 2023 15:44:06.67 UT by Fermi -GBMwith a duration of ≈ 35 s (Fermi GBM Team 2023; Dalessi et al. 2023). The burst was also observed by GECAM, the InterPlanetary Network (IPN), AGILE, AstroSAT, and GRBAlpha (Xiong et al. 2023; Kozyrev et al. 2023; Casentini et al. 2023; Navaneeth et al. 2023; Dafcikova et al. 2023), and was quickly noted as the second-brightest GRB seen by Fermi to date (Burns et al. 2023). The ULTRACAM instrument mounted on the 3.5 m New Technology Telescope (NTT) and Swift -XRT undertook wide-field searches for an optical and X-ray counterpart, respectively, discovering coincident candidates at δt = 1 . 4 -1 . 7 days (Burrows et al. 2023; Levan et al. 2024). The counterpart was offset 30.2 '' (38.9 kpc) from a bright spiral galaxy confirmed at z = 0 . 0646 (Levan et al. 2024), and an extensive multi-wavelength followup campaign was initiated.  \nDespite GRB230307A's nominal long duration, its high-energy properties, including its spectral lag and Xray flux decay, provided evidence for a compact object merger origin similar to GRB 211211A. Near-IR followup from ground-based observatories and JWST revealed a late-time red excess, light curve shape, and spectral features expected for a kilonova (e.g., Gillanders et al. 2023; Gillanders & Smartt 2024; Levan et al. 2024; Yang et al. 2024). We combine the multi-wavelength datasets presented in the literature, which includes Chandra , Swift , and XMM X-ray observations, broad coverage by ground observatories, late-time HST and JWST observations, and ATCA and AMI-LA radio observations (Levan et al. 2024; Yang et al. 2024) for our analysis.", '3. AFTERGLOW MODELING': "In addition to a radioactive decay-powered kilonova, BNS mergers are expected to launch a relativistic jet  \nwhose interaction with the surrounding medium produces broadband synchroton emission, or the 'afterglow' (e.g., Sari et al. 1998). The afterglow flux is a potential contaminant for modeling kilonovae, especially for on-axis events in which the afterglow often dominates the total luminosity at early times, motivating us to model the afterglows of each event in our sample, with the exception of GW 170817 (see Section 2.6). Our afterglow model uses the formulae of Granot & Sari (2002) and methods described in Laskar et al. (2014) to describe synchrotron emission from a forward shock (FS), produced by the interaction of the GRB's collimated jet and the surrounding medium, incorporating the effects of Inverse Compton cooling (Sari & Esin 2001; Laskar et al. 2015).  \nThe parameters fit in this afterglow model are the jet isotropic-equivalent energy ( E K , iso ), the circumburst density of the surrounding medium ( n 0 ), the input electron distribution power law index ( p ), the fraction of energy deposited into non-thermal relativistic electrons ( ϵ e ), the line-of-sight extinction in B -band ( A B ; fixed to 0 in our fits except for GRB 130603B, see below), and the time of the jet break ( t jet ), which is directly related to the jet opening angle (Rhoads 1999; Sari et al. 1999). For each event we assume a radially-homogeneous ISMlike environment ( k = 0), as this is expected in the local environments for most NS mergers. We fix the energy deposited in magnetic fields ( ϵ B ) to 0.01, the median for short GRBs (Fong et al. 2015; Schroeder et al. 2024), as the model fails to converge to a reasonable solution when leaving this parameter free. We fit the afterglow model using the Markov Chain Monte Carlo (MCMC) emcee package (Foreman-Mackey et al. 2013), enforcing a minimum 10% uncertainty on all detections to capture realistic measurement errors. We run each fit using 128 walkers for 5000 iterations and discard the first 10% of steps as burn-in. For each event, we employ the redshifts and extinction values listed in Table 1 and use values from the literature as starting parameters.  \nAs the kilonova and afterglow both contribute flux in the optical-to near-IR wavelengths, disentangling emission between the two can be difficult, especially at early times. In particular, observations of AT 2017gfo demonstrated that kilonovae may have early ( δt ≲ 1 day) blue emission due to large quantities of fast-moving lanthanide-poor ejecta or, potentially, additional energy sources such as free neutron decay (e.g., Metzger et al. 2015) or shock-heating (e.g., Kasliwal et al. 2017; Villar et al. 2017b; Arcavi 2018; Nicholl et al. 2021). Thus, we attempt to remain agnostic to the precise afterglow contribution in the optical and mask data in the range 10 13 -10 16 Hz (effectively, fitting only the radio and X-  \nFigure 1. The X-ray and radio detections (circles) and 3 σ upper limits (triangles) of the seven events in our sample with on-axis afterglows. In each panel, we show the model light curves' median and 68% confidence range (shaded region) along with the X-ray (grey) and radio (colored and labeled) observations. Open symbols denote data that was masked in the afterglow fit (Section 3). We show two models for GRB 230307A. The first ('Model 1') is fit to the X-ray and radio data only and provides a poor fit to the early TESS observations. The second ('Model 2') is fit with the TESS, radio and X-ray observations. It provides a better constraint on the optical contribution but cannot fully account for the late-time X-ray and radio observations.  \n<!-- image -->  \nobservations) in our fits to GRBs 160821B, 200522A, 211211A and 230307A. For GRB230307A, we find this method significantly underestimates the flux of TESS optical observations at 0 . 01 ≲ δt ≲ 0 . 2 days (Figure 1; 'Model 1') while a fit including the TESS observations provides a worse fit to the X-ray and radio data ('Model 2'). We perform our analysis on the results from both models due to the uncertainty in the emission source. For GRB160821B, we exclude the radio detection at δt = 0 . 17 days as it is likely the result of an early reverse shock (Lamb et al. 2019), and incompatible with the standard FS model. We do not expect the reverse shock to significantly contaminate the optical flux on the timescales of observations we use to model the kilonova ( δt > 0 . 95 days).  \nFor the remaining three events we include some early optical data in our fit either due to sparse X-ray and radio detections (GRBs 050709, 060614) or high line-of-  \nsight extinction that will significantly affect our estimation of the afterglow flux contribution (GRB 130603B; Cucchiara et al. 2013; Fong et al. 2014, 2015). While it is possible that some kilonova flux may be contributing in these optical detections, previous works have shown that these points can be explained with an afterglow model only (e.g., Fox et al. 2005; Fong et al. 2014; Xu et al. 2009). Specifically, for GRB 050709, we include two early ( δt < 2 . 4 days) R -band optical observations. For GRB 060614, we include early RIJK -band data ( δt = 0 . 7 -2 days) in our afterglow fit. We exclude X-ray data of GRB060614 prior to δt = 0 . 5 days as numerous analyses favor an energy injection scenario to explain the afterglow plateau observed at δt ≲ 0 . 5 day that is incompatible with the standard afterglow model (Figure 1; e.g., Xu et al. 2009; Mangano et al. 2007). We do not expect the exclusion of the energy injection episode in GRB 060614 to affect our kilonova modeling  \nas the FS is expected to dominate on the timescales and in the filters we use to model the kilonovae. Finally, for GRB130603B we include early ( δt ≲ day) opticalnear-IR detections to measure the line-of-sight dust, A B , which we propagate to the output models we use for subtraction. We also exclude the final XMM-Newton observation of GRB 130603B as it is known to be contaminated by an unrelated X-ray source (Rouco Escorial et al. 2023). We denote which observations were used in the afterglow fitting in Table 5.  \nFor all events, we visually inspect the optical-near-IR afterglow models to ensure they are not brighter than measured values beyond the uncertainties. In general, we find that our derived afterglow physical parameters are consistent with those in the literature, though in some cases we find inconsistent values (typically, in the degenerate parameters E K , iso and n 0 ). Variations are likely the result of discrepancies between modeling codes or the fact that we primarily use X-ray and radio data only; however the inferred afterglow parameters are not used in any subsequent analysis 3 , so the specific values are not important for this work. We create model light curves at the same rest-frame wavelengths as the opticalnear-IR observations for 1000 random draws from the full posterior of each event. From these 1000 draws we calculate the median and 68% credible flux range, using the 68% flux range as our uncertainty on the afterglow model. We show our fits to the multi-wavelength data, along with their uncertainties, in Figure 1. Following interpolation of the afterglow light curves to the δt of each observation, we subtract the median afterglow model flux from each observed flux, producing an 'afterglowsubtracted' light curve. We combine the data uncertainty and model uncertainty at the time of the observation in quadrature. We present our afterglow-subtracted values and errors in Table 4.", '4.1. Description of Kilonova Model': "We employ the Python-based Modular Open Source Fitter for Transients ( MOSFiT ) code (Guillochon et al. 2017) to fit kilonova models to the afterglow-subtracted data, for which the kilonova contribution has nominally been isolated. From this modeling, we derive physical parameters describing the observed kilonova emission, which we parameterize in terms of ejecta mass ( M ej ), velocity ( V ej ) and temperature cooling floor ( T floor ). The  \nTable 2. Kilonova Model Parameters & Priors  \nlatter is the temperature below which the photosphere recedes into the ejecta. Hence, T floor represents an effective emission temperature for the optically-thin nebular phase to follow. We elect to use MOSFiT as its modular design affords us the flexibility to build new modules, freeze or adjust parameters and their priors, add constraints, and test several samplers without a high computational cost. In addition, MOSFiT is a well-tested method that has been used to fit or determine upper ejecta mass limits of several past kilonovae (e.g., Villar et al. 2017b; Nicholl et al. 2021; Rastinejad et al. 2022a; Coulter et al. 2024). All of the results presented in this section were performed with nested sampling, as implemented in the dynesty fitting routine (Speagle 2020).  \nWe begin with MOSFiT 's existing three-component (blue, purple and red; see below) kilonova model (Villar et al. 2017a; Metzger 2019), which assumes analytic forms for the radioactive heating rate (Korobkin et al. 2012; equation 1 of Villar et al. 2017b) and the thermalization efficiency (Barnes et al. 2016; equation 2 of Villar et al. 2017b). The model calculates the bolometric luminosity assuming a central energy source and following an updated Arnett (1982) formalism (Chatzopoulos et al. 2012; equation 3 of Villar et al. 2017b). Here, each of the three components is parameterized by a constant 'grey' ejecta opacity ( κ ), the value of which correlates with lanthanide or electron fractions in portions of the total ejecta. Kilonova models comprised of two or three components were found to provide a better fit to the well-sampled AT 2017gfo compared to single component models (Cowperthwaite et al. 2017; Drout et al. 2017; Kasen et al. 2017; Kilpatrick et al. 2017;  \nTanaka et al. 2017; Villar et al. 2017b). They are also physically motivated by simulations that predict multiple ejecta mechanisms with distinct elemental compositions prior to and following the NS merger (e.g., Metzger & Fern'andez 2014; Kasen et al. 2017).  \nFor our three-component model we employ κ R = 10 cm 2 g -1 , κ P = 3 cm 2 g -1 , and κ B = 0 . 5 cm 2 g -1 for the 'red', 'purple' and 'blue' components, respectively (Tanaka et al. 2017). We expect these components to roughly map to the red (lanthanide-rich or Y e ≲ 0 . 2) tidal dynamical ejecta (e.g., Rosswog et al. 1999), the purple (moderately lanthanide-rich or 0 . 2 ≲ Y e ≲ 0 . 3) disk wind ejecta (e.g., Metzger & Fern'andez 2014; Just et al. 2015; Fern'andez & Metzger 2016; Lippuner et al. 2017), and the blue (lanthanide-poor or Y e ≳ 0 . 3) dynamical ejecta shocked at the NS contact interface and ejected near the poles (e.g., Sekiguchi et al. 2015) or ejected in a magnetized wind from the neutron star remnant prior to black hole formation (e.g., Metzger et al. 2018; Fong et al. 2021; Combi & Siegel 2023; Curtis et al. 2024). For each component in our model, we measure M ej , V ej and T floor (Table 2). We pursue a three-component model (as opposed to a two-component model) as this provides a better mapping of ejecta mechanism to opacity, enabling a search for trends tied to physical properties.  \nThe geometry of the mechanism producing each ejecta component and the viewing angle of the observer are known to significantly impact the observed light curve (e.g., Darbha & Kasen 2020; Chase et al. 2022), and thus any parameter inference. Of particular relevance to this work, the assumption of an isotropic kilonova will likely introduce a bias in estimating the mass of material ejected along the line-of-sight. Thus, under this assumption, GRB events observed along the jet axis likely have blue and red ejecta components that are overestimated and underestimated, respectively. Instead, here we account for the geometry of the ejecta by modifying the original spherical kilonova model (Villar et al. 2017b) to an aspherical model, wherein a half-opening angle ( θ open ) defines a conical boundary between red ejecta, confined to the equatorial region, and the blue and purple ejecta, modeled in the direction of the poles. We use the half-opening angle prescription of Darbha & Kasen (2020), as implemented in MOSFiT by Nicholl et al. (2021). The viewing angle ( θ obs ; defined relative to the axis of the GRB jet) of each event in our sample is wellestablished, either as relatively pole-on due to the detection of a cosmological GRB or measured through highprecision astrometry in the case of AT 2017gfo (Mooley et al. 2022). Thus, we fix θ obs = 22 · (the central value in the range given) for AT 2017gfo (Mooley et al. 2022) and  \nθ obs = 0 · for all other events in our sample. We allow the kilonova ejecta half-opening angle to be a free parameter and include an additional parameter ( σ ) to account for white noise in the likelihood function. These, in addition to M ej , V ej and T floor for each of the three components, comprise the 11 parameters measured for each kilonova.  \nWe acknowledge that the number of free parameters exceeds the number of data points for some events in our sample. We perform several test fits fixing each component's T floor and θ open , finding consistent masses compared to runs where these parameters are left free. In the end, we opt to keep these parameters free in order to marginalize over the uncertainties on other parameters when measuring masses. Uniform use of the three-component model is critical to comparisons between kilonovae. In these cases we find very broad posteriors on the masses (e.g., Figure 4 and Table 3), but are still able to constrain them compared to the broad uniform prior.  \nWe list the parameters and priors used in our kilonova model in Table 2. For M ej and V ej we elect to use the widest priors possible that correspond to physical values based on simulations (e.g., Radice et al. 2018; Nedora et al. 2021). For T floor we use the range 1000-4000 K following the reasoning of Nicholl et al. (2021) based on previous fits to AT 2017gfo with MOSFiT (Cowperthwaite et al. 2017; Villar et al. 2017b). For GRB 230307A, JWST detections at δt ≈ 29 and 61 days likely occur at epochs when nebular emission is either significantly contributing or dominating the observed flux. As kilonova nebular emission remains a challenge to properly model, we do not include the observations at δt ≈ 61 days in our kilonova fit. To accommodate the observations at δt ≈ 29 days and the mid-IR coverage of the JWST detections, we allow the T floor , P and T floor , R prior range to extend down to 800 and 440 K, respectively. We base our choice of T floor , R = 440 K on the blackbody fit to synthetic nebular lanthanide-rich kilonova spectra and Spitzer detections of AT 2017gfo (Kasliwal et al. 2022; Hotokezaka et al. 2021; Barnes & Metzger 2022). We choose a median value of T floor , P = 800 K as we expect the blackbody temperature of moderately lanthaniderich ejecta to fall in between that of lanthanide-poor and lanthanide-rich ejecta.  \nWe fit the eight afterglow-subtracted optical-near-IR light curves at δt > 0 . 5 days (Table 4) with the threecomponent kilonova model presented above. We do not include data prior to 0.5 days in our kilonova fits as these timescales may be affected by emission mechanisms be-  \nFigure 2. Afterglow-subtracted observations of the eight kilonovae in our sample along with best-fit (median and 68% confidence) model light curves. 3 σ upper limits are marked with triangles, while detections are marked with circles. We show light curves in the filters where observations (including upper limits) are available for each event. Two models are shown for GRB230307A based on the two afterglow models described in Section 3 and shown in Figure 1.  \n<!-- image -->  \nond radioactive decay, including energy injection in the afterglow (GRB 060614; Mangano et al. 2007), central engine activity (e.g., GRB 211211A; Hamidani et al. 2024a), or shock cooling (e.g., Rastinejad et al. 2022a). For our fits, we include only data points whose combined statistical and systematic (from afterglow fitting; Section 4) errors are < 0 . 5 mag and treat all observations with combined errors > 0 . 5 mag as upper limits.  \nIn general, this results in detections prior to ≈ 1 day being treated as upper limits (e.g., Figure 2 and Table 4). For GRB050709, we employ a threshold of < 0 . 6 mag, as we find the inclusion of the two additional points provides a significantly better fit. We present the median and 68% confidence range for each parameter measured in our kilonova modeling in Table 3.  \nTable 3. Kilonova Model Posteriors  \nNote -'B', 'P', and 'R' subscripts refer to the blue, purple and red ejecta components, as described in Section 4.", '4.2. Results and Observed M ej Diversity': "In Figure 2 we plot median and 68% confidence range model light curves, constructed from 900 random draws of the full posterior of each GRB (with the exception of GRB170817A for which we use 50 random draws due to computational constraints). Overall, we find that the models provide reasonable fits to the data and follow the predicted behavior of kilonovae: rapid decay, especially in bluer bands, and reddening over time. As expected, events with better observational coverage correspond to tighter constraints on the light curves. For AT 2017gfo, we find that our model provides a good fit to the optical data and the majority of the near-IR observations, but overpredicts the late-time ( δt ≳ 20 day) K -band observations, a feature also noted in previous MOSFiT fits to this event (Villar et al. 2017b). We discuss comparisons to previous fits further in Section 4.4.  \nFor GRB230307A we find that our choice of afterglow model significantly affects the shape of the model optical light curve at later times, reflected in the disparate best-fit values found for M ej , R and T floor , R . The fit to data subtracted with afterglow Model 1 (Figures 1 and 2; ''KN Fit 1') follows a typical fading behavior at δt ≳ 10 days. In contrast, the fit to data subtracted with afterglow Model 2 ('KN Fit 2'; in particular the JWST F070W detection at ≈ 29 days which is an upper limit in KN Fit 1) produces a flattening in the optical decay past δt ≳ 10 days. To produce this shape, KN Fit 2 requires higher values of T floor , B and T floor , R than KN Fit 1 to explain the emission at later times (Table 3). We  \nfavor KN Fit 1 in our later analysis (though we show both for completeness), for several reasons. First, we prefer the afterglow model that provides a better fit to the X-ray and radio data than the early TESS observations (afterglow Model 1; Figure 1) as there are a number of proposed emission mechanisms, including a reverse shock, shock cooling and free-neutron decay (Metzger et al. 2015; Gottlieb et al. 2018), that may explain the early TESS excess but could not account for a late-time radio and X-ray excess. Second, the true TESS bandpass is wider than the nominal I c -band reported in the literature (e.g., Yang et al. 2024), and may also explain excess emission relative to the model. Third, KN Fit 1 results in more physically realistic values for T floor , R , as they are similar to the blackbody temperature that approximates kilonova nebular emission at the wavelengths where the nebular observations occur (e.g., 440 K; Barnes & Metzger 2022). Finally, fits to both afterglow-subtracted datasets without JWST observations are more consistent with KN Model 1.  \nIn Figure 3 we show the median and 68% confidence range for M ej and V ej for each component and the total ejecta. In keeping with previous works (e.g., Rastinejad et al. 2021), our fits place the tightest constraints on the blue and purple component parameters and the coarsest constraints on the red component. This is in large part due to the traditional use of more sensitive optical telescopes for afterglow searches, especially for the kilonova candidates detected prior to GW 170817. Notably, our fits to events with just two or three detec-  \nFigure 3. Median M ej and V ej (68% confidence) for the blue (top left), purple (top right), and red (bottom left) components as well as the total ejecta (bottom right). We mark long GRBs with diamonds, GW 170817/AT 2017gfo with a star, and the remaining short GRBs with circles. Our analysis highlights that GW 170817/AT 2017gfo is a 'typical' kilonova.  \n<!-- image -->  \ntions and deep upper limits (GRBs 050709 and 200522A) place order-of-magnitude or tighter constraints on M ej , B and M ej , P . In Section 5.1 we analyze and compare the kilonova properties of long and short GRBs.  \nFocusing on the blue component, our analysis finds that all except one event prefer M ej , B ≤ 0 . 01 M ⊙ (Table 3). We also observe a general trend between increasing M ej , B and V ej , B (Figure 3), though the large error bars on several events prohibit a firm conclusion. Taking 1000 draws from the posterior of each event (and discarding GRB230307A KN Fit 2), we calculate a median of M ej , B = 0 . 006 +0 . 015 -0 . 004 M ⊙ (68% confidence; Figure 4).  \nThe blue kilonova emissions may be attributed to either dynamical ejecta heated at the contact surface between the NSs and ejected along the the axis of the jet (e.g., Oechslin et al. 2007; Sekiguchi et al. 2015) or post-merger disk ejecta experiencing neutrino irradiation from a NS remnant, which lowers the lanthaniderichness of any ejecta (e.g., Metzger & Fern'andez 2014; Miller et al. 2019). Notably, GRB200522A is a significant outlier in M ej , B (Figure 4) with M ej , B = 0 . 046 ± 0 . 009 M ⊙ , consistent with past findings (Fong et al. 2021). We slightly favor a disk-wind source (rather than a shock-heated dynamical source) to explain the majority of GRB 200522A's larger M ej , B as BNS merger sim-  \nions that measure shock-heated dynamical ejecta do not produce M ej , B ≳ 0 . 01 M ⊙ , even when spanning a range in NS masses, mass ratios and two NS equations of state (Sekiguchi et al. 2015). In contrast, simulations measuring a disk-wind mass in the case of a long-lived NS remnant produce M ej , B ≈ 0 . 03 M ⊙ (e.g., Lippuner et al. 2017). The extreme luminosity of GRB 200522A's kilonova has previously been explained with the creation of a magnetar remnant, which may provide an additional blue emission source (Fong et al. 2021).  \nTurning to the purple component, we find a population median of M ej , P = 0 . 020 +0 . 034 -0 . 010 M ⊙ (68% confidence; Figure 4). This range is broadly consistent with expectations of disk component masses (e.g., Lippuner et al. 2017). We find that several events (GRBs 050709, 130603B and 230307A) prefer lower ejecta velocities compared to those found for the blue or red components. This trend supports the purple component's source as a disk wind which is likely ejected with slower speeds compared to the dynamical red and blue components ( V ej ≈ 0 . 01 -0 . 1 c ; e.g., Metzger & Fern'andez 2014; Fern'andez et al. 2015). We find that GRB211211A is an outlier on the higher end in V ej , P . As our fit to this event finds high velocities for all three components, we posit this is likely a reflection of the fast-decaying UV-  \nFigure 4. The median and 68% confidence range of the M ej found for each GRB in the blue (top left), purple (top right), and red (bottom left) components as well as the total ejecta mass (bottom left). Long GRBs, short GRBs and AT 2017gfo are marked with diamonds, circles, and a star, respectively. We show the median and 68% confidence range across all events in the grey vertical line and shaded region, and mark the log-uniform prior for each component in the grey horizontal line. We also mark upper limits on the component ejecta masses from observations of short GRB afterglows with left-facing light purple triangles (Section 4.3). We show our two fits for GRB 230307A, marking KN Fit 2 with an open symbol as we disfavor this fit in our analysis (Section 4). GW 170817/AT 2017gfo falls comfortable within the 68% confidence range for each ejecta mass.  \n<!-- image -->  \nl emission observed at δt ≈ 0 . 5 -2 days, which previous fits have explained with a shock cooling model (e.g., Rastinejad et al. 2022a).  \nFinally, for the red component, we find a population median of M ej , R = 0 . 051 +0 . 100 -0 . 045 M ⊙ (68% confidence; Figure 4). Amongst the three components, this is highest median, and the widest range in ejecta masses, though we note there are poor constraints on this value for GRBs050709, 130603B and 200522A. The red kilonova component can be ascribed to the neutron-rich ejecta tidally stripped from the NS surfaces as the compact objects slowly inspiral (e.g., Lattimer & Schramm 1974; Li & Paczy'nski 1998; Metzger et al. 2010; Tanaka & Hotokezaka 2013). Generally, it is expected that red dynamical ejecta mass will increase with larger asymmetry  \nin the progenitor masses (particularly for NSBH events; e.g., Foucart et al. 2014) or higher spins (e.g., Kyutoku et al. 2015; Shibata & Hotokezaka 2019). Notably, the three highest M ej , R median values are found for the three LGRB events, GRBs060614, 211211A and 230307A. We further discuss the implications of this finding and other trends with GRB properties in Section 5.1.  \nAcross the sample, we derive a median total ejecta mass of M ej , tot = 0 . 085 +0 . 110 -0 . 040 M ⊙ (68% confidence). In every component and the total, the ejecta mass of GW170817/AT2017gfo falls comfortably within the 68% credible range found for all events. This indicates that AT2017gfo may be considered a 'representative' kilonova (Figure 4) compared to the seven kilonovae analyzed here. In contrast, the kilonova of GRB 160821B  \nfalls below the median ejecta masses found for all GRBs, rendering it a critical point in probing kilonova diversity.", '4.3. Constraints on M ej from Additional Short GRB Observations': "Next, we briefly explore the possibility that our kilonova sample is observationally biased towards more luminous events. As luminosity roughly scales with a fractional power of M ej (e.g., Metzger et al. 2018), a missing population of low-luminosity kilonovae would translate to a population with lower M ej than those reported here.To evaluate this, we employ observations of seven short GRBs with upper limits and afterglow detections that are less luminous compared to AT 2017gfo when matched in rest-frame time and band at their known redshifts (see Rastinejad et al. 2021, for details). These bursts 4 are GRBs050509B (Bloom et al. 2006; Cenko et al. 2005; Hjorth et al. 2005; Castro-Tirado et al. 2005), 080905A (Rowlinson et al. 2010a; Nicuesa Guelbenzu et al. 2012), 090515 (Rowlinson et al. 2010b), 100206 (Nicuesa Guelbenzu et al. 2012; Perley et al. 2012), 130822A Cenko et al. (2013); Rastinejad et al. (2021), 150120A (Rastinejad et al. 2021) and 160624A (O'Connor et al. 2021; Rastinejad et al. 2021). For this analysis, we assume all detections are dominated by afterglow flux and treat them as upper limits.  \nIn the observed frame for each event (determined with the redshift catalog of Fong et al. 2022), we generate three sets of one-component kilonova light curves, characterized by either the blue, purple, and red fixed opacity value mentioned in Section 4.1. We fix T floor = 1000 K and use the corresponding component's median velocity found in Section 4.2. For each component, we produce kilonova models log-spaced in M ej over the range M ej = 0 . 001 -0 . 5 M ⊙ . We compare our GRB observations to the set of kilonova models, and record the highest M ej in each component allowed by the upper limits for each event. In Figure 4 we plot constraining ( < 0 . 5 M ⊙ ) upper limits on the respective component ejecta masses. We note that our use of one-component models translates to a conservative upper limit as flux from other components is not taken into account in our procedure.  \nSimilar to previous findings, short GRB observations are most constraining of M ej , B (Figure 4; purple triangles). Past, deep rest-frame near-IR coverage is sparse. Thus, historical short GRB upper limits do not place meaningful constraints on M ej , P and M ej , R (Figure 4). In the blue component, we find one event less massive  \nthan AT2017gfo (GRB080905A; M ej , B < 0 . 002 M ⊙ ) and an additional two events with upper limits below the median value found in Section 4.2 (GRBs 050509B and 130822A; M ej , B < 0 . 005 M ⊙ ). Overall, the existing upper limits span the range of blue ejecta masses. Thus, it is difficult to conclude at present if the population that do not have detected kilonovae also have lower ejecta masses.", '4.4. Comparison to Previous Kilonova Fits': 'Comparing our results in Table 3 to those from previous fits, we find our results are generally consistent with those in the literature. We find variation in absolute differences, an expected outcome given the range of kilonova modeling codes used, which we further discuss here.  \nFor AT2017gfo, our fit produces a larger M ej , R and a smaller M ej , B and M ej , P compared to a previous MOSFiT run (Villar et al. 2017b), though a similar value for the total ejecta mass is reached ( M tot ≈ 0 . 07 M ⊙ ). We ascribe this difference in the relative component masses to our addition of the geometry prescription (Section 4.1; Darbha & Kasen 2020), which is expected to increase the amount of red mass relative to the blue for viewing angles less than the θ open (and is observed in Villar et al. 2017b though a larger θ obs and different asymmetry prescription was used). Nicholl et al. (2021) also model AT2017gfo with MOSFiT , incorporating constraints from GWobservations, and find comparable M ej , B and M ej , P . Their analysis finds a smaller M ej , R ( ≈ 0 . 001 M ⊙ ), which we attribute to their constraints on the tidal dynamical ejecta based on the BNS mass ratio and chirp mass (Nicholl et al. 2021). Compared to two-component (red and blue components only) fits to AT 2017gfo, we obtain M ej , B , R values that are within the range but on the upper end of those measured in the literature (e.g., Arcavi et al. 2017; Chornock et al. 2017; Kasen et al. 2017; Kasliwal et al. 2017; McCully et al. 2017; Pian et al. 2017; Smartt et al. 2017; Troja et al. 2017; Anand et al. 2023).  \nOf the remaining kilonovae in our sample, the majority have measured M ej and V ej . However, previous fits to these events use a variety of models and methods to measure the kilonova parameters and exact comparisons are not advisable. For GRB 060614, a previous estimate found M ej , tot ≈ 0 . 1 M ⊙ , on the lower end of our estimate (Yang et al. 2015). For GRB130603B, previous fits find a wide span in M ej , tot = 0 . 01 -0 . 1 M ⊙ (Berger et al. 2013; Tanvir et al. 2013; Barnes et al. 2016), broadly consistent with our results. For GRB 160821B, Troja et al. (2019) constrain a lanthanide-rich mass to < 0 . 006 M ⊙ and a lanthanide-poor mass to ≈ 0 . 01 M ⊙ .  \nFor the same burst, Lamb et al. (2019) find a dynamical ejecta mass of (1 . 0 ± 0 . 6) × 10 -3 M ⊙ and a post-merger ejecta mass of (1 . 0 ± 0 . 6) × 10 -2 M ⊙ . Within errors, our results are consistent with both findings but are on the upper end of the given ranges.  \nGRB211211A was fit with a three-component model in MOSFiT that included a shock heating prescription, finding M ej , B = 0 . 01 ± 0 . 001 M ⊙ , M ej , P = 0 . 01 ± 0 . 02 M ⊙ and M ej , R = 0 . 02 +0 . 02 -0 . 01 M ⊙ (Rastinejad et al. 2022a). They find consistent results when performing joint afterglow and kilonova modeling (Rastinejad et al. 2022a). Within errors, the blue and purple component masses are consistent with our results, but we find a larger median M ej , R by a factor of two. Additional M ej , tot estimates of this kilonova span 0 . 02 -0 . 1 M ⊙ (Troja et al. 2022; Yang et al. 2022; Kunert et al. 2024). For GRB230307A, M ej , tot of 0 . 05 +0 . 15 -0 . 05 M ⊙ (Levan et al. 2024) and ≈ 0 . 08 M ⊙ (Yang et al. 2024) are found in the literature. Taken together, our results are generally consistent with those found in the literature. Our median values are on the upper end of previous estimates, consistent with comparisons to AT 2017gfo modeling.', '4.5. Caveats to Kilonova Model': "We find that in several cases the error bars on M ej (particularly for the blue component) are unrealistically small. While a powerful tool to infer physical properties, MOSFiT makes a series of simplying assumptions that likely result in an underestimation of the true errors, similar to the conclusion made for MOSFiT modeling of tidal disruption events (Mockler et al. 2019).  \nIn particular, the assumption of a constant grey opacity may significantly affect the model posteriors, as M ej , V ej and κ are degenerate parameters in predicting the shape of the light curve. Notably, κ may vary up to an order of magnitude on the timescales of our observations (0 . 1 ≲ δt ≲ 30 days) as the ejecta temperature and density directly impact the elements' ionization states (e.g., Tanaka et al. 2020; Banerjee et al. 2024). We quantify the minimum systematic error introduced by our assumed opacities by running fits of AT2017gfo in which each component's opacity is a free parameter whilst holding the other two components' κ values constant and keeping the same prior ranges listed in Table 2. We then determine the difference in derived component M ej relative to the masses inferred in Section 4.1 with fiducial, fixed opacities. Specifically, we explore the effects of the range κ R = 5 -30 cm 2 g -1 (compared to κ R = 10 cm 2 g -1 in our Section 4.1 fits), κ P = 1 -5 cm 2 g -1 ( κ P = 3 cm 2 g -1 ), and κ B = 0 . 2 -2 . 5 cm 2 g -1 ( κ B = 0 . 5 cm 2 g -1 ). We use a minimum of κ B = 0 . 2 cm 2 g -1 as it corresponds to  \nthe minimum value dictated by Thompson scattering of ionized elements (Paczynski 1983). We base the remaining ranges on calculations for kilonova opacities at δt ≈ 1 day (Banerjee et al. 2024; Tanaka et al. 2020).  \nBased on the three fits with free κ values, we find the median errors on the component masses are ∆ M ej , R = ( -0 . 04 , +0 . 0008) M ⊙ , ∆ M ej , P = ( -0 . 004 , +0 . 02) M ⊙ and ∆ M ej , B = ( -6 . 5 × 10 -5 , 0 . 0006) M ⊙ . For the total mass, fitting with κ as a free parameter results in only modestly lower total ejecta masses of which the median value is δM ej , tot = -0 . 001 M ⊙ . We find that our fits varying κ P result in the most significant source of uncertainty. This result can naturally be explained as values of κ P more similar to κ B or κ R will divert a portion of the luminosity typically explained by M ej , B or M ej , R to M ej , P . This results in a different ratio between the two component masses but an overall similar M ej , tot . We conclude that our choice of opacity (within the range of values explored here) does not significantly impact our ejecta mass results. We do not include these uncertainties where we compare between our uniform fits.  \nWe caution that, in addition to the effect explored above, the kilonova models used in this analysis do not account for the effect of jet-ejecta interaction, shock cooling, central engine activity, or magnetic fields, all of which may play a large role in determining ejecta masses (e.g., Radice et al. 2018; Ciolfi & Kalinani 2020; Nativi et al. 2021; Shrestha et al. 2023; Curtis et al. 2024; Hamidani et al. 2024b). In addition, several pieces of kilonova physics are still not well understood even in state-of-the-art simulations, such as the uncertainty in wavelength-dependent opacities, nuclear heating rate, and thermalization efficiencies (Barnes et al. 2021; Bulla 2023; Brethauer et al. 2024; Sarin & Rosswog 2024). Further, our two-stage analysis of the afterglow and kilonova may result in additional uncertainty and bias in the derived kilonova parameters compared to a joint models (e.g., Wallace & Sarin 2024). To account for this point, we exclude data at < 0 . 5 days in our kilonova fits, the timescale on which the afterglow is most likely to dominate, and propagate the uncertainties in our afterglow model to the data passed to the kilonova model (Sections 3 and 4.1).  \nIn light of these uncertainties, we emphasize that the aim of this work is to perform uniform modeling on a sample of kilonovae, allowing for an exploration of diversity and correlations with γ -ray and environment properties. For all objects in our sample besides AT 2017gfo, the datasets are relatively sparse (e.g., Section 2 and Figure 2), limiting our ability to investigate each uncertainty listed above. We provide all observations (Ta-  \nle 5), including afterglow-subtracted photometry (Table 4), for the community to model with other existing, or future, codes.", '5.1. Comparison to γ -ray Properties and Implications for Long/Short Progenitors': 'Motivated by the unknown progenitor properties and/or mechanism driving merger-origin long GRBs, we next examine any trends between the kilonova ejecta and γ -ray properties. Historically, one factor in explaining the divide in γ -ray duration between BNS and stellar progenitors is their order-of-magnitude difference in the mass reservoir surrounding the central compact object. The significantly higher masses and larger physical size of long GRB progenitor stars lead to a longer timescale for accretion, translating to a longer-lived jet (e.g., Tchekhovskoy et al. 2011). To produce a long-lived GRB jet from a neutron star merger, previous works have posited the progenitors are a white dwarf-neutron star binary (Yang et al. 2022; Wang et al. 2024) or an NSBH (e.g., Rastinejad et al. 2022a; Gompertz et al. 2022), or that the merger remnant is a magnetar (e.g., Rastinejad et al. 2022a; Gompertz et al. 2022).  \nHere, we explore the theory that a long-lived, massive ( ≈ 0 . 2 M ⊙ ) accretion disk, a product of an asymmetric binary merger, is capable of powering longer-lived and more energetic GRBs (Gottlieb et al. 2023). The merger of an asymmetric binary, whether it is a BNS or an NSBH with a favorable mass ratio ( Q ≈ 3 -5; e.g., Kawaguchi et al. 2016), will produce a greater amount of lanthanide-rich tidal dynamical ejecta compared to a symmetric binary (e.g., Hotokezaka et al. 2013; Kyutoku et al. 2018; Kawaguchi et al. 2020b). Within our modeling framework, this translates to an expected trend between M ej , R and the duration and/or energy of the GRB.  \nTo investigate these possible trends, we compare the kilonova ejecta masses (Table 3) with the values of E γ, iso (described in Section 2.1), beaming-corrected E γ , and T 90 , rest (converted to the rest-frame using their respective redshifts; Table 1). To calculate E γ , we gather jet opening angle measurements ( θ j ) from the literature (Xu et al. 2009; Rastinejad et al. 2022a; Rouco Escorial et al. 2023). All events in our sample have measured θ j values from X-ray observations, with the exception of GRB050709 for which a lower limit is reported (Rouco Escorial et al. 2023). For GRB 230307A we use 5000 random draws from our afterglow model 1 (Section 3), as it is fit to the combined X-ray light curves from the literature, therein providing a tighter constraint on the jet break than previous analyses (Levan et al. 2024; Yang  \net al. 2024). Our analysis finds θ j = 3 . 95 +1 . 72 -0 . 65 deg. We calculate the beaming-corrected energy, given by,  \nE γ = [1 -cos( θ j )] × E γ, iso . (1)  \nIn Figure 5, we show the results of this analysis. We observe that M ej , B (with the exception of the outlier GRB200522A; Section 4.2), M ej , R and M ej , tot generally increase with higher values of E γ, iso and E γ . We further observe that M ej , R increases with longer T 90 . The trends with E γ, iso , E γ , and T 90 are most apparent with M ej , R (Figure 5, third column), though the error bars for several short GRB masses and θ j preclude a firm conclusion. We do not observe any apparent trends between M ej , P (Figure 5, second column) and γ -ray properties.  \nWe briefly test the statistical significance of any trends between the ejecta masses and the γ -ray properties. We apply the Pearson correlation coefficient (r-score) test to the datasets in each panel of Figure 5. We caution that this test is agnostic to a model, and thus does not probe all underlying physical motivations. We randomly draw values from the ejecta mass posteriors and calculate r- and p-scores with the E γ, iso , E γ , and T 90 values for 1000 iterations, producing distributions of 1000 rand p-values. We then determine the fraction of random draws that imply a significant correlation between the ejecta mass and γ -ray properties using a threshold of p < 0 . 05. We do not find that a significant fraction of the p-scores favor a correlation between the ejecta masses and any γ -ray properties. The strongest correlation is between M ej , R and T 90 , where 32% of p-scores indicate a significant correlation.  \nThough we do not find any statistically significant correlations, we observe that GRBs with T 90 , rest ≳ 2 s have higher median red ejecta masses ( M ej , R ≳ 0 . 05 M ⊙ ) compared to typical short GRBs ( M ej , R ≲ 0 . 02 M ⊙ ) potentially hinting at a bimodality. We observe a similar pattern with M ej , B and E γ (with the exception of GRB200522A), though this component may also originate in the post-merger disk wind (Section 4.2). The implication that M ej , B seems to trend more strongly with γ -ray properties compared to M ej , P may also suggest that long GRB mergers produce relatively blue disk winds, perhaps due to energy injection from the GRB. At present, our results could indicate asymmetric binaries as the progenitors of merger-driven long GRBs. However, to establish any statistically significant correlations, a larger population of joint GRB-kilonova detections with well-constrained ejecta masses is necessary.  \nFinally, we note that while asymmetric binaries are uncommon in the observed population of BNS systems (e.g., Tauris et al. 2017), LVK observations have revealed an asymmetric BNS (GW190425; Abbott et al.  \nFigure 5. The blue (first column), purple (second), red (third) and total (fourth) ejecta masses for each on-axis GRB event plotted against the respective γ -ray properties: E γ, iso (top row), beaming-corrected E γ (middle row), and T 90 (bottom row). Long GRBs are marked with diamonds, while short GRBs are marked with circles. As there is only a lower limit on the jet half-opening angle of GRB 050709 (Section 5.1) we mark the resulting upper limit on E γ with a triangle. We observe potential trends between E γ, iso , E γ and the red and total ejecta mass.  \n<!-- image -->  \n2020a) and an NSBH (GW230519; The LIGO Scientific Collaboration et al. 2024) merger. These detections may point to these asymmetric binaries being common in the Universe, with potential rates able to explain the increasing (but highly uncertain) rates of long GRBs from mergers.', '5.2. Kilonova Contribution to Universal r -Process': "At present, kilonovae are the only observationallyconfirmed source of r -process production in the Uni-  \nverse. Indirect observational evidence may favor the existence of a second, 'faster' heavy element nucleosynthesis channel. Specifically, to explain observations of r -process-enhanced metal poor stars in the Milky Way and dwarf galaxies (e.g., Ji et al. 2016; Hansen et al. 2017; Frebel 2018) a significant fraction of NS mergers are required to have short delay times from enrichment to star formation (e.g., Zevin et al. 2022). Simulations have demonstrated that collapsar and magnetohydrodynamical (jet-driven) supernovae could be this  \nsecond channel (e.g., Siegel et al. 2019; Mosta et al. 2018; Halevi & Mosta 2018). At present, observations of these candidates are limited and those that exist do not support r -process enrichment (Blanchard et al. 2023; Anand et al. 2024; Rastinejad et al. 2024). Here, we explore if the average r -process yield from kilonovae, calculated using the median ejecta masses calculated in Section 4 along with the current NS merger rates, is capable of producing the estimated r -process abundance in the Milky Way. In this analysis, we assume the kilonovae in our sample are created by BNS mergers only, though we note in Section 5.1 the kilonovae following long GRBs are favored to be from asymmetric mergers that may be NSBH events. We make this assumption for simplicity, as only a small fraction of NSBH events are expected to produce kilonovae and the rate of NSBH events is subdominant compared to the rates of BNS mergers (e.g., Abbott et al. 2023; Mandel & Broekgaarden 2022).  \nTo calculate the r -process enrichment in the Milky Way from NS mergers, we employ the equation from Rosswog et al. (2018):  \nM r ∼ 17 000 M ⊙ [ R BNS 500Gpc -3 yr -1 ][ ¯ m ej 0 . 03 M ⊙ ][ τ gal 1 . 3 × 10 10 yr ] (2)  \nwhere R BNS is the rate of NS mergers, ¯ m ej is the average kilonova ejecta mass, and τ gal is the age of the Milky Way, which we fix to 1 . 3 × 10 10 yr. We employ the range in BNS merger rate calculated from the Gravitational Wave Transient Catalog 3 (GWTC-3) of R BNS = 10 -1700 (Abbott et al. 2023). For ¯ m ej we employ our median for the eight kilonovae in our sample, M ej , tot = 0 . 085 +0 . 108 -0 . 044 M ⊙ , including the uncertainties due to the opacity (Section 4.5). Taking these values together, we find a wide range of M r ≈ 500 - 400,000 M ⊙ for BNS mergers, mostly driven by the large uncertainty in BNS merger rate.  \nThis range for M r encompasses the estimate of total r -process mass in the Milky Way, M r, MW ≈ 23 , 000 M ⊙ (Venn et al. 2004; Battistini & Bensby 2016; Hotokezaka et al. 2018) and, on the lower end, leaves room for the existence of a second r -process source. As discussed in Sections 4.3-4.5, our ¯ m ej value is likely to be an overestimate of the true value for two reasons. First, less luminous, and thus less massive, kilonovae are likely to have been missed due to observational biases (Section 4.3). Second, in keeping with fits to AT 2017gfo, our MOSFiT model derives ejecta mass values on the upper end of ranges from previous fits (Section 4.4). However, we anticipate that the uncertainty in ¯ m ej in subdominant compared to the uncertainty in BNS merger rate, as our values are generally comparable with literature values,  \nwhere they exist, (Section 4.4) and likely do not vary beyond an order of magnitude.  \nThe specific star formation rate of the host galaxy is a dominant factor in governing what fraction of ejected r -process elements enrich later generations of stars (Nugent et al., in prep.). Notably, in comparison to the quiescent host galaxy of GW 170817/AT 2017gfo (Blanchard et al. 2017; Levan et al. 2017), the hosts of the kilonovae in our on-axis GRB sample are all star-forming (Nugent et al. 2022; Levan et al. 2024; Table 1). Though losses may still be significant for these latter events, it is more likely that they enriched their galaxies with heavy elements compared to AT 2017gfo.", '5.3. Future Kilonova Observations': "We next consider the implications of our sample and model light curves for future wide-field kilonova searches, either triggered with a GRB or GW event or untriggered (e.g, Smartt et al. 2017; Doctor et al. 2017; Kasliwal et al. 2017; Andreoni et al. 2020). For these searches, an understanding of the span in kilonova light curve diversity is critical in the vetting process, which has yet to be enabled with observations beyond AT2017gfo.  \nFor each of the seven on-axis GRB events in our sample, we create 900 light curves in the observed-frame rizJHK -bands. Each light curve is based on a random draw from the posterior (we use GRB 230307A Model 1 in our fit given the reasoning in Section 4.2) and matched to the approximate rest-frame band in Figure 6. We create 50 light curves for GRB 170817A/AT 2017gfo due to computational constraints and iterate over these 18 times. Due to the tight posteriors on this event, we do not expect this process to affect our conclusions. Across the random draws, we calculate the median, 68% and 90% credible region in luminosity space for the sample of kilonovae. In Figure 6 we show the median and kilonova luminosity range in each filter. Taken together, our results indicate that future kilonovae will span ≳ one order of magnitude in luminosity.  \nAs we are motivated by both targeted and untargeted kilonova searches in large surveys, we also plot the single image 5 σ depths for the Rubin Observatory (Bianco et al. 2022), the WINTER J -band limiting magnitude (Frostig et al. 2022) and the 'Wide Tier' ( RZ -band) and 'Deep Tier' ( JH -band) limiting magnitudes for Roman High Latitude Time Domain Survey (HLTDS; Rose et al. 2021), shifted to z = 0 . 1. We caution that our results are based mostly on on-axis events, and the bolometric luminosity of kilonovae may vary up to factor of ∼ 10 with viewing angle (Darbha & Kasen 2020). In the riz -bands, Rubin observatory will be sensitive  \nFigure 6. The median (solid line), 68% credible range (dotted lines) and 90% credible range (shaded region) in luminosity space shown in the rizJHK -bands and calculated using random draws (light grey lines) from each event. We also show the expected depths of several current and upcoming wide-field observatories that plan to engage in kilonova searches (Margutti et al. 2018; Rose et al. 2021; Bianco et al. 2022; Frostig et al. 2022; Andreoni et al. 2024).  \n<!-- image -->  \nto the full and upper end of the 68% credible range of kilonovae at z = 0 . 1 out to δt ≈ 3 and ≈ 7 days, respectively. As shown with our analysis GRB 200522A, in Section 4.2, a single epoch of simultaneous rest-frame optical observations on these timescales is sufficient for order-of-magnitude constraints on M ej , B . We therefore find that order-of-magnitude constraints on M ej , B are possible for kilonovae observed by Rubin at one epoch. We acknowledge that large outstanding challenges exist to distinguish these rare events from other transients, especially in just one or two epochs of observations, and encourage further development of automatic vetting tools.  \nRoman's HLTDS will be a powerful tool in measuring the properties of low-redshift kilonovae if events can be identified. Indeed, Roman is capable of detecting z = 0 . 1 kilonovae in the optical and near-IR out to δt ≈ 2 and 3 weeks, respectively. If the Roman HLTDS observations occur at a five-day cadence, Roman is poised to observe z = 0 . 1 kilonovae over ≈ 2-3  \nepochs in the optical and ≈ 3-5 epochs in the near-IR. Two to three epochs is sufficient for obtaining better than order-of-magnitude constraints on M ej -V ej , as we have demonstrated in Section 4. Finally, WINTER will be sensitive to the most luminous kilonovae, but not the 68% credible range, at z = 0 . 1.  \nIn addition to wide-field searches in the nearby Universe, well-localized Swift GRBs continue to be a promising method to detect kilonovae. These events come with their own set of challenges, including higher redshifts, bright afterglows and low rates. Further, the more coarsely localized detections of GRBs by Fermi , Space Variable Objects Monitor (SVOM; Atteia et al. 2022), the InterPlanetary Network (IPN; e.g., Atteia et al. 1987; Svinkin et al. 2022) and other γ -ray telescopes offer a second route to finding kilonovae using 'targeted' wide-field surveys. However, for all GRB kilonovae, the higher average redshifts render follow-up with large-aperture ground-based telescopes and spacebased observatories critical. Indeed, seven of the eight  \nkilonovae in our sample had key observations with HST that were critical to detections on ≳ week timescales. Looking to the future, JWST can obtain similar observations for events out to z ≈ 1 (e.g., Rastinejad et al. 2021). For lower-redshift events JWST has the power to capture the kilonova spectral energy distribution (SED) in the nebular phase, an important element in refining the mapping from observations to M ej .", '6. CONCLUSIONS': "We have compiled and collated the multi-wavelength light curves of eight kilonovae from the GCNs and literature. Five of these events follow short GRBs, while three events follow long GRBs, allowing us to explore trends between γ -ray and kilonova ejecta properties. We uniformly model the afterglows of seven events with on-axis GRBs, producing 'afterglow-subtracted' light curves. We fit the afterglow-subtracted light curves with a threecomponent kilonova model in MOSFiT that accounts for geometric viewing effects. Our fits provide reasonable fits to the data, and we compare our posteriors to those in the literature. Our major conclusions are as follows:  \n- · Our fits unveil a wide span in derived kilonova properties, namely M ej and V ej , implying that the progenitors and/or remnants of these mergers are also diverse. We determine that the luminous kilonova of GRB 200522A has a significantly more massive M ej , B compared to the sample of events (or the luminosity may be boosted by nonradioactive heating like a magnetar; Fong et al. 2021), while the kilonova of GRB 160821B is the least massive of the total sample.\n- · While well-sampled events provide the tightest constraints, we also find value in kilonovae with a single color measurement, particularly if their colors are unique (e.g., GRBs 130603B and 200522A).\n- · We discuss the main uncertainties in our modeling (Section 4.5) and compare our results to previous fits with other modeling codes. We emphasize that all observations used in this work, including our 'afterglow-subtracted' light curves, are provided for future modeling endeavors in Tables 4 and 5.\n- · We demonstrate that GW170817/AT2017gfo is a 'representative' kilonova in each components' M ej and V ej (Table 3), when compared to the seven kilonovae in our sample and deep upper limits from the literature. Given this, our estimate of the total Milky Way r -process mass produced by kilonovae does not change significantly when using the median ejecta mass of our sample compared to previous estimates made for AT 2017gfo.\n- · We explore trends between our derived ejecta masses and E γ, iso , beaming-corrected E γ and T 90 . Overall, we do not find any statistically significant correlations but observe that long GRB kilonovae have larger median M ej , R compared to short GRB kilonovae. We hypothesize that this is indicative of an asymmetric binary merger origin for longerlived GRBs. A larger sample of well-studied kilonovae following short and long GRBs will be critical to confirming this hypothesis.\n- · We produce median, 68% and 90% confidence range light curves in a variety of bands based on the posteriors of the eight events in our sample. Comparing these light curves to the expected depths of upcoming surveys, we anticipate that Rubin and Roman will be sensitive to the majority of the kilonova luminosity range for z ≲ 0 . 1 and are capable of order-of-magnitude mass constraints.  \nHere, we have shown that the existing sample of kilonovae, the majority of which are detected at a fixed viewing angle, demonstrates diversity and trends with γ -ray properties. Widening the sample of these events requires dedicated strategies for observational pointings (e.g., Margutti et al. 2018; Coulter et al. 2024) and kilonova candidate vetting (e.g., Rastinejad et al. 2022b) that take into account the full diversity of compact binary merger EM counterparts and environmental properties. The advent of next-generation GW detectors, deep wide-field surveys, and new γ -ray instruments, when combined with dedicated search strategies, opens the doors for unprecedented exploration into the physics of compact binary mergers, jets, and kilonovae.", '7. ACKNOWLEDGEMENTS': 'The authors thank Ore Gottlieb, Genevieve Schroeder, Anya Nugent, Tanmoy Laskar, Igor Andreoni, Sylvia Biscoveneau, Nick Kaaz, Ben Margalit and Michael Fausnaugh for helpful conversations regarding this manuscript.  \nJ.C.R. acknowledges support from the Northwestern Presidential Fellowship. The Fong Group at Northwestern acknowledges support by the National Science Foundation under grant Nos. AST-2206494, AST-2308182, and CAREER grant No. AST-2047919. W.F. gratefully acknowledges support by the David and Lucile Packard Foundation, the Alfred P. Sloan Foundation, and the Research Corporation for Science Advancement through Cottrell Scholar Award 28284.  \nTable 4 . Afterglow-subtracted Observations  \nTable 4 continued  \nTable 4 (continued)  \nTable 4 continued  \nTable 4 (continued)  \nTable 4 continued  \nTable 4 (continued)  \nTable 4 continued  \nTable 4 (continued)  \nTable 4 continued  \nTable 4 (continued)  \nTable 4 continued  \nTable 4 (continued)  \nTable 4 continued', 'Rastinejad et al.': "Table 5 (continued)  \nNote - Observations are not corrected for Galactic nor local extinction.  \n† Starred observations were employed in our afterglow analysis.  \nReferences : (1) Fox et al. 2005, (2) Hjorth et al. 2005, (3) Covino et al. 2006, (4) UKSSDC (Evans et al. 2007, 2009), (5) Mangano et al. 2007, (6) Schmidt et al. 2006, (7) Xu et al. 2009, (8) Cobb 2006, (9) Fynbo et al. 2006, (10) Della Valle et al. 2006, (11) Yang et al. 2015, (12) Londish et al. 2006, (13) Fong et al. 2014, (14) de Ugarte Postigo et al. 2014, (15) Cucchiara et al. 2013, (16) Berger et al. 2013, (17) Tanvir et al. 2013, (18) Troja et al. 2019, (19) Lamb et al. 2019, (20) Kasliwal et al. 2017, (21) Fong et al. 2021, (22) O'Connor et al. 2021, (23) Mei et al. 2022, (24) Ito et al. 2021, (25) Troja et al. 2022, (26) Xiao et al. 2024, (27) Kumar et al. 2021, (28) Rastinejad et al. 2022a, (29) Strausbaugh & Cucchiara 2021, (30) Mao et al. 2021, (31) Moskvitin et al. 2021, (32) Gupta et al. 2021, (33) Yang et al. 2024, (34) Levan et al. 2024.", 'REFERENCES': "Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016, Living Reviews in Relativity, 19, doi: 10.1007/lrr-2016-1 -. 2017a, ApJL, 848, L12, doi: 10.3847/2041-8213/aa91c9 -. 2017b, PhRvL, 119, 161101, doi: 10.1103/PhysRevLett.119.161101 -. 2020a, ApJL, 892, L3, doi: 10.3847/2041-8213/ab75f5 Abbott, R., Abbott, T. D., Abraham, S., et al. 2020b, ApJL, 896, L44, doi: 10.3847/2041-8213/ab960f Abbott, R., Abbott, T. D., Acernese, F., et al. 2023, Physical Review X, 13, 011048, doi: 10.1103/PhysRevX.13.011048 Alam, S., Albareti, F. D., Allende Prieto, C., et al. 2015, ApJS, 219, 12, doi: 10.1088/0067-0049/219/1/12 Anand, S., Pang, P. T. H., Bulla, M., et al. 2023, arXiv e-prints, arXiv:2307.11080, doi: 10.48550/arXiv.2307.11080 Anand, S., Barnes, J., Yang, S., et al. 2024, ApJ, 962, 68, doi: 10.3847/1538-4357/ad11df Andreoni, I., Ackley, K., Cooke, J., et al. 2017, PASA, 34, e069, doi: 10.1017/pasa.2017.65 Andreoni, I., Kool, E. C., Sagu'es Carracedo, A., et al. 2020, ApJ, 904, 155, doi: 10.3847/1538-4357/abbf4c Andreoni, I., Coughlin, M. W., Criswell, A. W., et al. 2024, Astroparticle Physics, 155, 102904, doi: 10.1016/j.astropartphys.2023.102904 Arcavi, I. 2018, ApJL, 855, L23, doi: 10.3847/2041-8213/aab267 Arcavi, I., Hosseinzadeh, G., Howell, D. A., et al. 2017, Nature, 551, 64, doi: 10.1038/nature24291 Arnett, W. D. 1982, ApJ, 253, 785, doi: 10.1086/159681  \nAscenzi, S., Coughlin, M. W., Dietrich, T., et al. 2019, MNRAS, 486, 672, doi: 10.1093/mnras/stz891  \nAtteia, J. L., Cordier, B., & Wei, J. 2022, International Journal of Modern Physics D, 31, 2230008,  \ndoi: 10.1142/S0218271822300087  \nAtteia, J. L., Barat, C., Hurley, K., et al. 1987, ApJS, 64, 305, doi: 10.1086/191198  \nBanerjee, S., Tanaka, M., Kato, D., & Gaigalas, G. 2024, ApJ, 968, 64, doi: 10.3847/1538-4357/ad4029  \nBarnes, J., & Kasen, D. 2013, ApJ, 775, 18,  \ndoi: 10.1088/0004-637X/775/1/18  \nBarnes, J., Kasen, D., Wu, M.-R., & Mart'ınez-Pinedo, G.  \n2016, ApJ, 829, 110, doi: 10.3847/0004-637X/829/2/110  \nBarnes, J., & Metzger, B. D. 2022, ApJL, 939, L29,  \ndoi: 10.3847/2041-8213/ac9b41  \n- -. 2023, ApJ, 947, 55, doi: 10.3847/1538-4357/acc384\n- Barnes, J., Zhu, Y. L., Lund, K. A., et al. 2021, ApJ, 918, 44, doi: 10.3847/1538-4357/ac0aec\n- Barthelmy, S., Barbier, L., Cummings, J., et al. 2006, GRB Coordinates Network, 5256, 1  \nBarthelmy, S. D., Barbier, L. M., Cummings, J. R., et al.  \n2005, SSRv, 120, 143, doi: 10.1007/s11214-005-5096-3  \nBattistini, C., & Bensby, T. 2016, A&A, 586, A49,  \ndoi: 10.1051/0004-6361/201527385  \n- Belles, A., D'Ai, A., & Swift/UVOT Team. 2021, GRB Coordinates Network, 31222, 1  \nBennett, C. L., Larson, D., Weiland, J. L., & Hinshaw, G. 2014, ApJ, 794, 135, doi: 10.1088/0004-637X/794/2/135 Berger, E., Fong, W., & Chornock, R. 2013, ApJL, 774, L23, doi: 10.1088/2041-8205/774/2/L23  \ndoi: 10.1093/mnras/stv009  \nRicker, G. R., Atteia, J. L., Crew, G. B., et al. 2003, in American Institute of Physics Conference Series, Vol. 662, Gamma-Ray Burst and Afterglow Astronomy 2001: A Workshop Celebrating the First Year of the HETE Mission, ed. G. R. Ricker & R. K. Vanderspek (AIP), 3-16, doi: 10.1063/1.1579291  \nRoming, P. W. A., Kennedy, T. E., Mason, K. O., et al. 2005, SSRv, 120, 95, doi: 10.1007/s11214-005-5095-4 Rose, B. M., Baltay, C., Hounsell, R., et al. 2021, arXiv e-prints, arXiv:2111.03081,  \ndoi: 10.48550/arXiv.2111.03081  \n- Rossi, A., Stratta, G., Maiorano, E., et al. 2020, MNRAS, 493, 3379, doi: 10.1093/mnras/staa479\n- Rossi, A., Rothberg, B., Palazzi, E., et al. 2022, ApJ, 932, 1, doi: 10.3847/1538-4357/ac60a2\n- Rosswog, S., Liebendorfer, M., Thielemann, F. K., et al. 1999, A&A, 341, 499,  \ndoi: 10.48550/arXiv.astro-ph/9811367  \nRosswog, S., Sollerman, J., Feindt, U., et al. 2018, A&A,  \n615, A132, doi: 10.1051/0004-6361/201732117  \n- Rouco Escorial, A., Fong, W., Berger, E., et al. 2023, ApJ, 959, 13, doi: 10.3847/1538-4357/acf830  \nRowlinson, A., Wiersema, K., Levan, A. J., et al. 2010a,  \nMNRAS, 408, 383, doi: 10.1111/j.1365-2966.2010.17115.x  \nRowlinson, A., O'Brien, P. T., Tanvir, N. R., et al. 2010b,  \nMNRAS, 409, 531, doi: 10.1111/j.1365-2966.2010.17354.x  \nSari, R., & Esin, A. A. 2001, ApJ, 548, 787,  \ndoi: 10.1086/319003  \n- Sari, R., Piran, T., & Halpern, J. P. 1999, ApJL, 519, L17, doi: 10.1086/312109\n- Sari, R., Piran, T., & Narayan, R. 1998, ApJL, 497, L17, doi: 10.1086/311269\n- Sarin, N., & Rosswog, S. 2024, arXiv e-prints, arXiv:2404.07271, doi: 10.48550/arXiv.2404.07271\n- Savchenko, V., Ferrigno, C., Kuulkers, E., et al. 2017, ApJL, 848, L15, doi: 10.3847/2041-8213/aa8f94\n- Schlafly, E. F., & Finkbeiner, D. P. 2011, ApJ, 737, 103, doi: 10.1088/0004-637X/737/2/103  \nSchmidt, B., Peterson, B., & Lewis, K. 2006, GRB  \nCoordinates Network, 5258, 1  \n- Schroeder, G., Fong, W., & Alexander, K. D. 2020, GRB Coordinates Network, 27786, 1  \nSchroeder, G., Fong, W.-f., Kilpatrick, C. D., et al. 2024,  \narXiv e-prints, arXiv:2407.13822,  \ndoi: 10.48550/arXiv.2407.13822  \nSekiguchi, Y., Kiuchi, K., Kyutoku, K., & Shibata, M. 2015, PhRvD, 91, 064059, doi: 10.1103/PhysRevD.91.064059 Shappee, B. J., Simon, J. D., Drout, M. R., et al. 2017, Science, 358, 1574, doi: 10.1126/science.aaq0186  \nZevin, M., Nugent, A. E., Adhikari, S., et al. 2022, ApJL, 940, L18, doi: 10.3847/2041-8213/ac91cd  \nZheng, W., Filippenko, A. V., & KAIT GRB Team. 2021, GRB Coordinates Network, 31203, 1", '8. APPENDIX': 'Table 5 . Multi-wavelength Observations of GRB Afterglows and Kilonovae  \nTable 5 continued  \nTable 5 (continued)  \nTable 5 continued  \nTable 5 (continued)  \nTable 5 continued  \nTable 5 (continued)  \nTable 5 continued  \nTable 5 (continued)  \nTable 5 continued  \nTable 5 (continued)  \nTable 5 continued  \nTable 5 (continued)  \nTable 5 continued  \nTable 5 (continued)  \nTable 5 continued  \nTable 5 (continued)  \nTable 5 continued  \nTable 5 (continued)  \nTable 5 continued  \nTable 5 (continued)  \nTable 5 continued  \nTable 5 (continued)  \nTable 5 continued  \nTable 5 (continued)  \nTable 5 continued  \nTable 5 (continued)  \nTable 5 continued  \nTable 5 (continued)  \nTable 5 continued'}