A Hubble Space Telescope Search for r-Process Nucleosynthesis in Gamma-Ray Burst Supernovae

The existence of a secondary (in addition to compact object mergers) source of heavy element (r-process) nucleosynthesis, the core-collapse of rapidly rotating and highly magnetized massive stars, has been suggested by both simulations and indirect observational evidence. Here, we probe a predicted signature of r-process enrichment, a late-time (≳40 days post-burst) distinct red color, in observations of gamma-ray burst supernovae (GRB-SNe), which are linked to these massive star progenitors. We present optical to near-IR color measurements of four GRB-SNe at z ≲ 0.4, extending out to >500 days post-burst, obtained with the Hubble Space Telescope and large-aperture ground-based telescopes. Comparison of our observations to models indicates that GRBs 030329, 100316D, and 130427A are consistent with both no enrichment and producing 0.01–0.15 M ⊙ of r-process material if there is a low amount of mixing between the inner r-process ejecta and outer supernova (SN) layers. GRB 190829A is not consistent with any models with r-process enrichment ≥0.01 M ⊙. Taken together the sample of GRB-SNe indicates color diversity at late times. Our derived yields from GRB-SNe may be underestimated due to r-process material hidden in the SN ejecta (potentially due to low mixing fractions) or the limits of current models in measuring r-process mass. We conclude with recommendations for future search strategies to observe and probe the full distribution of r-process produced by GRB-SNe.


INTRODUCTION
For decades, the formation sites for many of the Universe's heaviest elements remained unknown.These elements, including many of those that enable life on Earth such as Thorium and Iodine, are widely believed to form via rapid neutron capture ("r-process") nucleosynthesis (Burbidge et al. 1957;Cameron 1957).This process requires origin sites capable of producing high neutron abundance fractions that are common enough to create the observed yields within a Hubble time.The traditional core-collapse of massive stars, in which material close to the remnant is neutronized (e.g., Cowan et al. 1991;Woosley et al. 1994;Hoffman et al. 1997;Thielemann et al. 2011;Cowan et al. 2021), was initially proposed as an origin site but recently has fallen out of favor (as discussed further below).Binary neutron star (BNS) or neutron star-black hole mergers, in which a high neutron fraction is naturally achieved through tidal stripping of the outer layers of the neutron star in the final orbits before merging (e.g., Lattimer & Schramm 1974;Rosswog et al. 1999) or through outflows from the post-merger accretion disk (e.g., Metzger et al. 2009), became the favored channel to produce r-process elements.In 2017, BNS mergers were confirmed to produce at least a fraction of the Universe's r-process elements with the first BNS merger detected through gravitational waves (GW170817; Abbott et al. 2017a), coincident short gamma-ray burst (GRB 170817A; Abbott et al. 2017b;Goldstein et al. 2017;Savchenko et al. 2017) and thermal kilonova AT 2017gfo (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).
Despite the spectacular confirmation of r-process in BNS mergers, it remains an open question as to whether other heavy element formation channels are needed to explain the mass and abundance pattern of r-process elements observed in our Solar System.Indirect evidence in the form of early heavy element enrichment in our Galaxy provides a convincing argument for a site tied to the collapse of massive stars (e.g., Côté et al. 2017;Hotokezaka et al. 2018;Naidu et al. 2022).Observations of r-process-enhanced metal-poor stars in the Milky Way, some ultra-faint dwarf galaxies (e.g., Ji et al. 2016;Hansen et al. 2017;Frebel 2018), and globular clusters (e.g., Zevin et al. 2019;Kirby et al. 2020Kirby et al. , 2023) ) suggest the existence of a heavy element formation channel with a short delay from star formation (Kirby et al. 2023 find a delay ≲ 0.8 Myr).Current estimates of the minimum delay time of BNS mergers, inferred from short GRB host galaxy offsets, stellar population ages, and star formation histories (e.g., Fong et al. 2022;Nugent et al. 2022;O'Connor et al. 2022;Zevin et al. 2022;Nugent et al. 2023) and predicted by simulations (e.g., Belczynski et al. 2002;Dominik et al. 2012;Tauris et al. 2017;Mandhai et al. 2022) remain uncertain, but it is unlikely that these events can provide a ubiquitous rprocess enrichment source as quickly as massive star channels.
Core-collapse SNe (CCSNe), associated with the deaths of massive stars and the formation of compact objects, provide a natural source with a short delay from star formation (on a stellar evolutionary timescale).In the 1990s, it was hypothesized that ordinary, jet-less CCSNe could produce neutron-rich material due to the high-entropy neutrino wind formed around the protoneutron star remnant (e.g., Woosley et al. 1994).However, high event rates and observations of ordinary CC-SNe have led to a general consensus that these events do not significantly contribute to the Universe's r-process budget (Macias & Ramirez-Ruiz 2018;Wallner et al. 2015Wallner et al. , 2021;;although see Tsujimoto & Shigeyama 2001 who find that SN 1987A's Ba/Sr ratio is more consistent with r-process than slow neutron capture, s-process).In addition, simulations have struggled to create or eject any neutron-rich material before accretion onto the remnant object (e.g., Arcones et al. 2007;Fischer et al. 2010;Martínez-Pinedo et al. 2012, 2014).
One suggested CCSNe source that could potentially decrease the event rates and overcome the ejection problem is magneto-rotational SNe (MR SNe; also termed jet-driven SNe or magneto-rotational hypernovae).In these MR SNe, rapid rotation of the iron core amplifies the magnetic field (e.g., Cameron 2003;Mösta et al. 2015) launching jets and subsequent magnetocentrifugal winds (e.g., Thompson et al. 2004;Metzger et al. 2007) that provide a mechanism to eject neutronrich material (see however Mösta et al. 2014).More recent simulations of MR SNe have found that strong (B ≈ 10 13 G) and efficient magnetic fields are critical for this pathway to the r-process, though it is unknown how common magnetic fields of this strength are in nature (e.g., Halevi & Mösta 2018;Mösta et al. 2018;Thompson & ud-Doula 2018).At present, no models for observational signatures of r-process produced through the MR SN production channel have been published.
Recently, Siegel et al. (2019) demonstrated that accretion disks following rapidly-rotating massive stars undergoing core-collapse ("collapsars") may also create and successfully eject r-process elements.Using magnetohydrodynamic (MHD) simulations they determined that the post-collapse disk favors weak interactions that produce neutron-rich material capable of synthesizing the heaviest elements.Their simulations also demonstrate that heating due to the disk's magnetic turbulence is sufficient to unbind the neutron-rich material, ejecting material in winds that will undergo the r-process, mix with the associated SN's outer layers and produce a red color signature analogous to the reddening of a kilonova (though on longer timescales; Siegel et al. 2019;Zenati et al. 2020).Additional simulations have further explored the dependency of r-process yields on neutrino treatment and accretion rate, and have disfavored collapsars as sources of heavy (e.g., lanthanide-rich) rprocess material (Miller et al. 2020;Fujibayashi et al. 2022;Just et al. 2022).However, given the uncertainties on these parameters (e.g., assumed range of accretion rates, treatment of disk viscosity, assumed stellar structure and rotation profiles of progenitor models), rprocess enrichment from collapsars remains plausible.While the MR SN and collapsar mechanisms represent distinct stages in the process of stellar core collapse, they are not mutually exclusive for a given event and in fact are associated with similar progenitor stars.Barnes & Metzger (2022) develop a semi-analytic light curve model for collapsar SNe enriched with r-process, which predicts their photometric color evolution for sufficiently large r-process enrichment levels.Their models produce a distinguishable red excess that emerges several weeks to months following the initial explosion.The model suite spans a range of r-process enrichment masses and degrees of mixing between the inner disk ejecta and SN outer layers, which imprint themselves on the light curves in the form of distinct optical to near-IR (NIR) colors.
Despite extensive work on the theoretical end, few observational searches have been performed for rprocess enhancement in SNe associated with collapsars.Searches for radio flares following long-duration GRBs (typically the product of collapsars) have been performed, a possible signature of interaction of the collapsar's wind ejecta with the surrounding medium (Lee et al. 2022).While no late-time radio flares were uncovered, there are several potential alternate sources that could explain radio emission in the event of such a discovery.Additionally, the high and sustained photospheric velocities inferred from observations of SN 2020bvc have been interpreted as power by a heavy element mixing source such as the r-process (Li et al. 2023).However, other interpretations such as interaction with circumstellar material (CSM) or shock cooling remain plausible explanations for this event (e.g., Izzo et al. 2020;Ho et al. 2020;Jin et al. 2021) and there are no signs of reddening in the SN light curve.Anand et al. (2024) perform a comprehensive search for r-process signatures in collapsar SNe using contemporaneous optical-NIR color measurements for a sample of 25 nearby broad-lined, stripped-envelope SNe (SNe Ic-BL) mostly discovered by the Zwicky Transient Facility (including one associated with a GRB, GRB 190829A).SNe Ic-BL are often associated with collapsars (MacFadyen & Woosley 1999), but may be explained with other mechanisms (e.g., Kashiyama et al. 2016).Their sample of SNe Ic-BL light curves was best fit by r-process-free models, favoring no or low r-process yields from nearby SNe Ic-BL detected mostly without GRBs.
The SNe Ic-BL associated with long GRBs (GRB-SNe) are generally considered the best targets for observable r-process enrichment in collapsars, in part due to the high angular momentum required to produce large accretion disks capable of launching the GRB jet (e.g., MacFadyen & Woosley 1999;Siegel et al. 2019;Barnes & Duffell 2023).In addition, Barnes & Duffell (2023) find that GRB jets are likely to increase mixing between the inner r-process ejecta and the outer layers, thus producing a more prominent red color.As this effect is likely enhanced closer to the jet axis, GRB-SNe, which have relatively pole-on orientations, are strong candidates for observing r-process signatures.Blanchard et al. (2023) do not find signs of r-process enrichment from a late-time NIR spectrum of the SN counterpart to GRB 221009A, though this observation was complicated by a bright afterglow.
The absence of a previous search for photometric signatures of r-process enrichment in a sample of GRB-SNe is in part due to the low rates of events within the requisite volume to detect the faint signatures (e.g., z ≲ 0.4) .For many past low-redshift events, no published late-time NIR data exists.The paucity of these measurements reflects the NIR sensitivity often required to study even low-redshift SNe on the timescales of these signatures (30 ≲ δt ≲ 300 days, where δt is the time since the GRB trigger).
Here, we provide late-time Hubble Space Telescope (HST) and large-aperture ground-based color measurements for a sample of four nearby GRB-SNe extending to ≳ 500 days following the GRB trigger.In Section 2 we describe our sample selection and detail the observations.In Section 3 we describe our process of ascertaining the intrinsic SNe colors.In Section 4 we compare our observations to the r-process enriched SNe models of Barnes & Metzger (2022).In Section 5 we review the implications of our work and discuss future observing strategies.Throughout, 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.

Sample Selection and Data Description
As a starting point, we utilize the comprehensive GRB-SNe compilation of Dainotti et al. (2022), which includes 58 long-duration GRBs with claimed SNe observed from 1990-2021.Following our motivation to identify NIR photometric excesses in the SN light curve due to the presence of r-process material, we narrow this sample to GRB-SNe events (i) for which late-time (δt ≳ 30 days), nearly simultaneous optical-NIR observations are available, (ii) identified after the discovery of SN 1998bw, as we do not expect previous GRB-SNe to have well-sampled light curves, (iii) at z < 0.4, the approximate distance out to which HST is capable of observing the predicted color evolution of rprocess enriched GRB-SNe, and (iv) not associated with a putative kilonova (e.g., GRBs 211211A and 230307A; Rastinejad et al. 2022;Troja et al. 2022;Yang et al. 2022;Gillanders et al. 2023;Levan et al. 2023a;Yang et al. 2023).Applying these criteria, our final sample comprises four GRB-SNe: GRB 030329 (SN 2003dh), GRB 100316D (SN 2010bh), GRB 130427A (SN 2013cq), andGRB 190829A (SN 2019oyw).We also considered but ultimately did not include GRB 221009A in our sample, due to a combination of high Galactic extinction (e.g., Williams et al. 2023), sparse late-time sampling of observations, and afterglow contamination in early HST epochs (Levan et al. 2023b).We list the basic properties of these bursts in Table 1 and describe each GRB and our data reduction further in Sections 2.2-2.5.Throughout this work, we refer to each event, including the SN, with its GRB name.
In total, we collect 79 HST observations obtained with the Advanced Camera for Surveys (ACS) Wide Field Channel (WFC), Wide Field Camera 3 (WFC3) Ultraviolet-Visible (UVIS) and Infrared (IR) channels, and the Near Infrared Camera and Multi-Object Spectrometer (NICMOS) Camera 2 (NIC2) instruments from the MAST archive 1 and the Hubble Legacy Archive (HLA)2 .The observations span 8 < δt < 969 days post-burst (including template observations) and seven HST filters (Table 1).We also collect and reduce VLT/X-shooter (acquisition camera), VLT/HAWK-I and MMT/Binospec observations of GRB 190829A (Section 2.5).With the exception of the NICMOS/NIC2 imaging, all reported photometry is performed on image subtractions with a late-time (δt ≳ 420 days) template (see Section 2.6 for further discussion on potential template contamination).In Figure 1 we show representative HST images where the SN is detected (left column) and the template image used for image subtraction (right column).We report all photometry in Table 2 and plot the observations in Figure 2. In Section 3.1 we discuss our corrections for Galactic and local dust extinction.
To complement the photometry analyzed in this work, we gather additional relevant data from the literature.As our goal is to compare the optical-NIR color evolution of the GRB-SNe in our sample to relevant models, we collect only host-subtracted photometry in the rRiIJHK-bands of GRBs 030329 (Matheson et al. 2003), 100316D (Olivares E. et al. 2012), 130427A (Perley et al. 2014) and 190829A (Hu et al. 2021).This results in an additional 519 observations from the literature.The vast majority of these measurements occur at 0.1 ≲ δt ≲ 40 days, extending our dataset at early times.
GRB 030329 was observed with ACS/WFC in the F606W and F814W filters at several nearly contemporaneous epochs (within ≈ 24 hours) over 17 ≲ δt ≲ 228 days, and at δt = 428 days in F606W only (Pro- Note- † T90 duration as seen by Swift, except for GRB 030329 which was observed by HETE-II, converted to the rest frame using the redshifts listed here.‡ Filters selected based on data availability on timescales relevant for observing r-process-enriched component (contemporaneous colors observed at δt ≳ 30 days).* Low-luminosity GRB.Milky Way extinction values are taken from Schlafly & Finkbeiner (2011).References: (1) Matheson et al. (2003) gram 9405; PI: Fruchter).We download the flat-fielded, dark-subtracted and CTE-corrected images and combine them using astrodrizzle (Gonzaga et al. 2012) with a pixel scale of 0.05 ′′ .We note that while the F814W image observed on 12 November 2003 was partially contaminated by internally scattered light from the WFPC2 internal lamp, the position of the SN is not affected by the uneven background.
We perform subpixel alignments between coadded images using tweakreg (Gonzaga et al. 2012) and align in image coordinates using standard IRAF tasks.We employ HOTPANTS (Becker 2015), which uses point-spread function (PSF) convolution, and IRAF/imarith for image subtraction.In general, the subtraction methods produce consistent photometry.We use HOTPANTS as our default image subtraction software throughout this work as it often produces higher signal-to-noise residuals than imarith, likely because of the treatment of the slight variations in PSF due to focus and orientation changes.In select cases we find that HOTPANTS returned a pattern that does not resemble a point-source residuals, motivating us to employ imarith.For each residual image, we detect a ≳ 3σ residual consistent with the position of the SN (Matheson et al. 2003) in the subtractions.We utilize the tabulated HST zeropoints to calibrate our images and perform aperture photometry with a 3-4 pixel aperture on the subtracted images (corresponding to ∼ 1× FWHM), accounting for the appropriate encircled energy corrections (Bohlin 2016).
GRB 030329 was also observed in the F110W and F160W filters with NICMOS/NIC2 at δt ≈ 17, 23, 44 and 228 days.We download the drizzled images from the HLA.We discard the F110W and F160W images observed at δt = 17 days due to saturation of the SN.As noted in previous works, these images suffer from known artifacts ( Östlin et al. 2008) and have a narrow field-of-view, preventing robust alignment and thus, reliable image subtraction.Thus, we perform relative photometry to obtain our measurements by subtracting the flux of the host galaxy in the epoch at δt = 228 days from those in each of the initial three epochs.We correct our photometry for the NICMOS/NIC2 encircled energy corrections listed in the NICMOS Handbook3 .

GRB 100316D
The sub-energetic, low-luminosity GRB 100316D was discovered with the Neil Gehrels Swift Observatory (Swift; Gehrels et al. 2004) at 12:44:50 UT on 16 March 2010.The burst duration was 277 s and its spectrum was noticeably soft (Sakamoto et al. 2010).The Swift X-Ray Telescope (XRT) promptly identified and localized an X-ray counterpart.Shortly thereafter, multiple telescopes observed an optical counterpart embedded in a galaxy at z = 0.0591 (e.g., Chornock et al. 2010;Starling et al. 2011).Spectroscopic features of a SN Ic-BL were observed in the optical counterpart a few days following the initial trigger (Chornock et al. 2010;Starling et al. 2011;Bufano et al. 2012;Olivares E. et al. 2012).
GRB 100316D was observed with WFC3/UVIS in the F555W and F814W filters, and with WFC3/IR in the F110W and F160W filters (Programs 11709, 12323; PI: Bersier) over 8 ≲ δt ≲ 505 days.We reduce the HST data in the same manner as described in Section 2.2.We perform image subtraction for each epoch using HOTPANTS and photometry on the HOTPANTS residual images using a 3-4 pixel aperture (∼ 1 − 2× FWHM) and account for the appropriate encircled energy correc- We do not expect the SN or afterglow to be significantly contributing to any of our templates, with the potential exception of GRB 130427A in F160W (see Section 2.6).
tions4 .We note that the SN in the F125W and F160W images on 4 April 2010 (δt ≈ 16 days) is saturated and we do not include these images in our analysis.We find our early photometry is in reasonable agreement with the analysis of this data reported in Cano et al. (2011).At δt ≳ 24 days, we find our values are fainter by ∼ 0.3−1.2mag compared to their values.This difference can be ascribed to our process of image subtraction in comparison to their direct photometry on the hostembedded SN (Cano et al. 2011).

GRB 130427A
GRB 130427A was detected by Swift and the Fermi Space Telescope (Fermi ;Meegan et al. 2009) on 27 April 2013 at 07:47:06.42UT.Its gamma-ray properties were unprecedented: at the time, the GRB had the highest observed fluence and most energetic photon recorded to date (von Kienlin 2013;Ackermann et al. 2014; holding the fluence record until the recent discovery of GRB 221009A; e.g., Burns et al. 2023).The burst duration as seen by Swift was 244 s.Prompt optical afterglow spectroscopy identified the GRB at z = 0.3399 (Levan et al. 2013) and multi-wavelength follow-up revealed an extraordinarily luminous afterglow (e.g., Laskar et al. 2013;Maselli et al. 2014;Perley et al. 2014).Spectroscopy of the counterpart confirmed the presence of an associated SN Ic-BL (Xu et al. 2013;Melandri et al. 2014;Levan et al. 2014a).
GRB 130427A was observed with ACS/WFC and WFC3/UVIS in the F606W filter and WFC3/IR in the F160W filter for 7 contemporaneous epochs over 74 ≲ δt ≲ 969 days (Programs 13110, 13117, 13230, 13951;PIs: Fruchter, Levan).We perform image subtractions between the initial six observations and the template image observed at δt = 969 days, then perform photometry on the residual with a 3 pixel aperture.Due to striping on the F160W image observed on 2014 December 5, we do not see a ≳ 3σ significant residual and discard the image.For all other epochs, we perform photometry at the position of the residual with a 3 pixel aperture and correct for the encircled energy.

GRB 190829A
GRB 190829A was identified by the Swift, Fermi and Konus-Wind satellites on 29 August 2019 at 19:55:53 UT (time as discovered by Fermi) and promptly localized to a galaxy with a known redshift of z = 0.0785 (Dichiara et al. 2019;Fermi GBM Team 2019;Tsvetkova et al. 2019).The burst had a gamma-ray duration of 53 s.The early afterglow was detected across the electromagnetic spectrum (e.g., Rhodes et al. 2020;H. E. S. S. Collaboration et al. 2021;Dichiara et al. 2022), and later spectroscopic observations revealed features consistent with an SN Ic-BL (Hu et al. 2021;Anand et al. 2024).
GRB 190829A was observed with HST with WFC3/UVIS in F606W and WFC3/IR in F140W over six contemporaneous epochs spanning 87 ≲ δt ≲ 500 days (Programs 15510, 16042, 16320;PIs: Levan, Tanvir).We also incorporate early (29 < δt < 58 days) images observed in the F110W and F160W filters (Program 15089; PI: Troja).We download, combine and align the F606W and F140W images using the procedure described above, and perform image subtractions with HOTPANTS.A residual at the SN position is detected at ≳ 3σ significance in all subtracted images.We perform photometry on all residual images (and directly on the F110W and F160W images as host contamination is insignificant at these epochs) with a 3-4 pixel aperture and correct for the encircled energy using the tabulated values.
In addition, we present VLT and MMT imaging of GRB 190829A observed over 25 < δt < 141 days.VLT observations were obtained in the rband with the X-shooter acquisition camera and in the JHK s -band with the HAWK-I instrument (Programs 0103.D-0819(A), 2103.D-5067(A), 105.20N7.001,0104.D-0600(E), 103.202P.002;PIs: Levan, Tanvir).We retrieve the X-shooter images from the ESO archive facility and reduce r-band images using standard IRAF/ccdproc tasks.We use the fully reduced HAWK-I images from the ESO archive.We reduce the r-band MMT/BINOSPEC (Program 2019C-UAO-G199, UAO-G205-23B; PIs: Fong, Rastinejad) images using a custom Python pipeline5 .For all bands, we utilize a template image obtained with the same instrument at δt ≳ 500 days for HOTPANTS image subtractions.We calibrate our images using SDSS (r; Alam et al. 2015) or 2MASS (JHK; Skrutskie et al. 2006) stars in the field.Finally, where a residual at the position of the SN is detected at > 3σ level, we perform aperture photometry with IRAF.We report all our observations in Table 2.

Assessing Transient Contamination in Template Images
Across all events, we consider the potential for afterglow or SN contamination in our HST template images, which may affect image subtraction residuals and bias our reported colors.Given that our sources are embedded in their host galaxies and likely to occur in starforming regions (e.g., Blanchard et al. 2016;Lyman et al. 2017), thus making simple visual inspection difficult, we apply an analytic model to determine the expected magnitude of the SN or afterglow at the time of the template image.
As we will show in Section 3.2, for GRBs 030329, 100316D and 190829A we expect the SN to dominate at late times, while for GRB 130427A we expect the afterglow to be the main source of emission at the time of the template image (Figure 2).To assess the SN contribution for the former GRB-SNe, we fit a simple 56 Ni and 56 Co decay model (e.g., Arnett 1982;Tinyanont et al. 2022;Kilpatrick et al. 2023) to the late-time (δt ≳ 110 days) optical (typically, the r, R, F606W and/or F555W bands) light curve.We then extrapolate the expected magnitude to the time of the template image.We find that the SN contribution in the template images ranges between m = 27.8 − 30.7 mag with the greatest potential contribution for GRB 100316D.However, as all of the GRB 100316D photometry is significantly brighter than this limit (m ≤ 24.6 AB mag) we consider any template contribution negligible.To assess the afterglow contribution for GRB 130427A we extrapolate the afterglow decay to late times using the parameters described in Section 3.2, and predict any SN or afterglow contamination to be m F606W ≳ 28.6 mag and m F160W ≳ 28.0 mag.We note that this may indicate transient contamination that would impact later (δt ≳ 300 days) F160W observations.However, as we justify further in Section 3.2, we do not use these latetime color measurements to probe SN r-process enrichment as they are likely dominated by afterglow emission.Further, as our earlier (δt ≈ 22, 74 days) SN-dominated observations of GRB 130427A are brighter than 23 mag, we do not expect contamination in the template to affect our SN colors.Dust extinction along the GRB's line of sight may significantly contaminate the observed optical-NIR color of the SNe in our sample though, notably, it does not impact the relative color evolution.However, it is important to take into account because particularly high extinction may result in a reddened light curve that will affect any identification or constraints on r-process material.In addition to corrections for Galactic extinction, we correct for extinction from the local environment of the GRB (Table 1), the assumptions of which we discuss below.
We utilize local or, when not available, host galaxy extinction values (E(B − V ) loc ) from the literature (Matheson et al. 2003;Cano et al. 2011;Bufano et al. 2012;Levan et al. 2014a;Chand et al. 2020;Huang et al. 2023) and apply a Small Magellanic Cloud (SMC) extinction law (Cardelli et al. 1989;Gordon et al. 2003) to determine corrections in each filter.We employ the SMC extinction law as it is a well-studied dust relation for a stellar population with a higher specific star formation rate compared to the Milky Way, appropriate for long GRB host galaxies (Kann et al. 2006;Schady et al. 2012).We note that the choice of extinction law will not significantly impact our optical-NIR corrections (as opposed to the ultraviolet regime).
For GRBs 030329, 100316D and 130427A, the derived values are low (E(B − V ) loc < 0.15 mag; Table 1) and where multiple values exist in the literature, they are consistent within ≈ 0.10 mag.For GRBs 100316D and 130427A, we take values of E(B − V ) loc measured from spectroscopy of the afterglows covering the Na I D doublet (λλ5890, 5896) and calculated using a Milky Way gas-to-dust ratio (Table 1; e.g., Bufano et al. 2012;Xu et al. 2013).For GRB 130427A this value is consistent with E(B − V ) loc calculated from the multiwavelength afterglow spectrum (Perley et al. 2014).For GRB 030329 there is no measurement from afterglow spectroscopy or fitting, so we utilize a value derived from the Hα/Hβ ratio in the host galaxy spectrum (Mathe-son et al. 2003).Literature values for local extinction of GRB 190829A are significantly higher.Using the XRT and Swift/Ultra-Violet Optical Telescope (UVOT) spectrum, Chand et al. (2020) derive E(B − V ) loc = 1.04.We note there is also a measurement E(B −V ) loc = 0.64 derived from the XRT light curve alone (Huang et al. 2023), but we prefer the value that incorporates UVOT data as this encompasses the regime where dust is most pronounced.

Characterizing the Afterglow Contribution
The GRB afterglow, typically modeled as synchrotron emission from the interaction of the jet with the surrounding medium, is expected to dominate over any SN emission in the few days following a burst (though this timescale may be highly variable).On week-to monthlong timescales, SN emission typically dominates in the optical band given the afterglow's steep power law decay (though there are notable exceptions in which the afterglow dominates for several months; e.g., GRB 221009A; Levan et al. 2023b).In the majority of cases, we expect that follow-up on timescales beyond a ∼month will probe the SN's color evolution.Therefore, the HST observations extending to very late times should be free from afterglow contamination.However, as using color evolution to determine r-process-enrichment in a GRB-SN is our primary objective, it is critical to determine which observations are clearly dominated by the SN emission.
Hence, we undertake a literature search to determine the temporal evolution (F ∝ t α , where F is the observed flux) for GRB afterglows to extrapolate as needed to the timescales of our observations.For GRB 030329 we utilize the post jet-break (observed at ≈ 0.5 day) afterglow decay of α = −2 (Lipkin et al. 2004;Moss et al. 2023).We note that at early times this afterglow is complex and contaminated by flares (e.g., Lipkin et al. 2004;Tiengo et al. 2004;Kamble et al. 2009), but we expect its latetime behavior to decline smoothly.For GRB 100316D, the optical afterglow was faint compared to the SN even at early times, and we use the late-time X-ray temporal slope of α = −0.87(Margutti et al. 2013).For GRB 130427A we use α = −1.45measured from the multi-wavelength afterglow (Perley et al. 2014).We note that Maselli et al. (2014) observe a jet break in the Xray and optical light curves at δt ≈ 10 hours post-burst, followed by a decline of α = −1.36.We prefer the values of Perley et al. (2014) which do not detect a jet break as their dataset includes greater temporal coverage (see also De Pasquale et al. 2016 who do not find evidence for a jet break in X-ray observations out to ∼ 2.5 years).Finally, for GRB 190829A we utilize the r-band temporal GRB 030329 t -2 t -0.87 t -1.45 t -1.45 index of α = −1.45from early GTC imaging (Hu et al. 2021).Apart from GRB 130427A, none of the GRBs in our sample have late-time jet breaks suggested in the literature.
We employ the above values for α and early observations to extrapolate the afterglow decay in each band and determine in which observations the predicted afterglow flux dominates our observations or is within 0.3 mag of the observed value.We plot the expected rband afterglow contributions in Figure 2 as grey dashed lines.In Figure 2 and Table 2 we denote observations that are likely afterglow-dominated, which we ignore in the analysis that follows.We do not anticipate any sig-nificant afterglow contamination in the light curves of GRBs 030329, 100316D and 190829A past δt ≈ 4 days.However, under our assumption of no observed jet break, the afterglow significantly contributes to or dominates the SN in the light curve of GRB 130427A at δt ≳ 200 days (Figure 2).

Supernova Dust Contribution
We briefly consider the possibility that dust produced by the SN may contribute to any observed reddening.The timescales for the production of dust in CC-SNe remain uncertain, though new observations with JWST are beginning to constrain the dust abundance in Type II SNe on ∼decade timescales (e.g., Hosseinzadeh et al. 2023;Shahbandeh et al. 2023).However, unlike these Type II SNe, our sample includes only relativistic, stripped-envelope GRB-SNe, as identified by the broad-lined features in their spectra.The mean absorption velocity measured for a large sample of GRB-SNe is ∼ 0.07c (Modjaz et al. 2016; see also Mazzali et al. 2021), indicating a high shock speed that would destroy any precursor dust grains before they may amass significantly.In addition, with the exception of a few rare cases, stripped-envelope SNe like GRB-SNe are rarely observed to have significant CSM, making any contribution from pre-existing dust unlikely (e.g., Prentice et al. 2019;Szalai et al. 2021).
In principle, newly-formed dust may contaminate our later observations.We do not find strong observational or theoretical evidence in the literature that significant dust is formed in GRB-SNe on ≲ 3-year timescales.There is evidence for reddening due to dust formation in other stripped envelope (e.g., SN 2013ge; Drout et al. 2014) or highly energetic Type I superluminous SNe (SN 2017ens; e.g., Sun et al. 2022).However, these events cannot be directly compared to GRB-SNe due to either lower ejecta velocities (in the case of SN 2013ge) or evidence for greater CSM interaction (Margutti et al. 2023).Notably, the early optical light curves of our sample of GRB-SNe resemble that of an afterglow with no excess luminosity due to reprocessing of emission in an ejecta-CSM shock interaction (Figure 2).We compare our extinction-corrected SN observations to the semi-analytic radiation transport models of Barnes & Metzger (2022) for a collapsar SN enriched with r-process material.In this collapsar r-process scenario, weak interactions within the dense and hot accretion disk feeding the black hole favor neutronization of the disk material (above a critical accretion rate De & Siegel 2021).A fraction of this neutron-rich midplane material is then ejected in disk winds, which undergoes r-process nucleosynthesis in the outflow on large scales.This process is similar to the disk outflows which contribute significantly to r-process production and kilonova emission in neutron star mergers (Siegel et al. 2019;Barnes & Metzger 2022).However, unlike in a merger, r-process disk wind ejecta collide and subsequently mix with the outer (nonr-process enriched) layers of the SN.The degree and radial extent of this mixing are uncertain and may be enhanced in the presence of a jet, resulting in a viewing angle dependence (Barnes & Duffell 2023).Barnes & Metzger (2022) predict that r-process material mixed with typical SN Ic-BL ejecta will produce a distinguishable red photometric color (distinct from the natural blue to red evolution of r-process-free models) that becomes pronounced on timescales of a few weeks to months.Naturally, larger dynamic ranges in filters hold larger discriminating power between models.The existing suite of models assumes a spherical distribution of ejecta parameterized by total SN ejecta mass (M tot ), 56 Ni ejecta mass (M 56Ni ), average ejecta expansion velocity as a fraction of the speed of light (β ej ), mass of r-process material (M rp ) and the mixing coordinate (ψ mix ).ψ mix describes the radial extent out through the SN ejecta to which the r-process enriched layers are mixed homogeneously (Barnes & Metzger 2022;Anand et al. 2024).The limit ψ mix = 1 corresponds to the rprocess material being homogeneously mixed all the way to the ejecta surface), resulting in a more pronounced and prolonged red color (Figure 4).The R − H color difference between the r-process-free and enriched cases may be up to ≈ 3 mag but is typically ≲ 0.1 mag prior to δt ∼ 30 days and in the cases ψ mix ≲ 0.2 (Barnes & Metzger 2022).
We attempt to constrain the number of free parameters so we are left with only M rp and ψ mix .To determine the values of M ej , β ej , and M 56Ni that produce model light curves most comparable to our observations, we convert our observations to the rest-frame, and compare observations at δt rest = 12 − 200 days to the grid of r-process-free (M rp = 0) models spaced in values of M ej , β ej , and M 56Ni (described in Barnes & Metzger 2022).We determine the best-fit parameters using χ 2 minimization.We show the best-fit models and their parameters along with relevant rest-frame observations in Figure 3.Our best-fit parameters are comparable to those found in the literature for GRBs 100316D, 130427A and 190829A (e.g., Cano et al. 2017a;Hu et al. 2021) though we find a lower M ej for GRB 030329 compared to previous analyses (Mazzali et al. 2003;Cano et al. 2017a).
To ensure that our initial choice of best-fit model based on a grid of r-process-free models does not bias our later conclusions about heavy element enrichment, we also run χ 2 minimization over the full grid of models, including both enriched and unenriched cases.For GRBs 030329, 130427A, and 190829A the best-fit model remains the r-process-free case shown in Figure 3.However, for GRB 100316D the best-fit model is enriched with M rp = 0.01M ⊙ and highly mixed (ψ mix = 0.9; dotted lines in Figure 3).The values of M ej , β ej , and M 56Ni Figure 3. Best-fit r-process-free (solid lines) models of four GRB-SNe in our sample and relevant rest-frame observations from our sample.We show the best-fit values of Mej, βej, and M56Ni for each GRB-SNe.We also show the best-fit model for GRB 100316D, the only event whose minimum χ 2 model was enriched with r-process, when compared to the total (including r-process enriched) grid of models.We fix these values of Mej, βej, and M56Ni for each GRB-SNe throughout our analysis.
for this model only vary slightly from those found for the r-process-free best-fit model, and we do not expect the colors to vary significantly based on these parameters.We thus conclude that our use of SN parameters determined from fits to r-process-free models in subsequent analysis will not significantly affect our conclusions about enrichment.
We next plot these models against our observations in color space to determine if the color evolution of any GRB-SNe in our sample resembles that modeled for rprocess enrichment.We use the models described above and shown in Figure 3 enriched with M rp = 0.03M ⊙ , a moderate value consistent with theoretical yields (e.g., Siegel et al. 2019).In Figure 4 we show the three best-sampled colors for each burst.We calculate colors for observations taken within three days of each other across ground-and space-based facilities.We combine color errors in quadrature, and incorporate an additional 0.1 mag error term to account for differences between the model (output in Johnson filters) and HST bandpasses.Models are available in the Johnson U BV RIJHK-bands.We compare our observations to the nearest approximate rest-frame model band and correct for time dilation using the redshift of the GRB (Table 1).We explore how the GRB-SN color evolution compares to models as a function of mixing fraction, and also compare to an r-process free model.We consider color measurements to be consistent with the model if they fall within 2σ errors.
For GRB 030329, the latest available colors are restframe V − R, where model color differences are small between the r-process-free and enriched cases.Still, we find that half the V − R measurements past 40 days are consistent with the highly-mixed case, while all are consistent with enrichment and ψ mix < 0.5.GRB 030329's optical-NIR (V − J and I − J) colors are only consistent with the r-process-free model or low mixing ψ mix ≲ 0.2 model.
For GRB 100316D, our rest-frame V −I HST measurements are also both consistent with r-process-enriched (ψ mix ≲ 0.5) and unenriched models.For this GRB, the discrepancy between F814W and ground-based i-band measurements can likely be ascribed to their different bandpass coverages and the steepening of the observed spectrum around 7500 − 8000 Å (Chornock et al. 2010).Optical-NIR observations of GRB 100316D do not appear to favor either enriched or r-process-free models, though, as we note above, results are inconsistent between telescopes.Later NIR observations of this burst would be necessary to distinguish between enriched and unenriched models.In GRB 130427A, the uncertainties and timing of restframe ground-based B−R and B−H data are not appropriate to distinguish between strong mixing or r-processfree models.The one optical-NIR color measurement at late times, (B − J at δt rest = 55 days) is consistent with only the enriched model with ψ mix = 0.1 − 0.2.Unfortunately, observations at ≳ 180 days are likely contaminated by the afterglow and are thus not appropriate for this analysis (Section 3.2; Figure 2).
Finally, the V − J colors of GRB 190829A are reasonably matched to both the r-process-free and low ψ mix models.On the other hand, the V − H and V − K colors are bluer than even the bluest r-process-free model.We note that this effect is also found for a number of SNe Ic-BL found without GRBs (Anand et al. 2024).Overall, we find that this GRB-SN provides a poor fit to the models, though the colors are more in line with the trends of the r-process-free case.
In summary, we find that, with a few exceptions, observations of GRBs 030329, 100316D and 130427A favor no enrichment or low values of ψ mix for M rp = 0.03M ⊙ (Figure 4), and inference between filters may vary.Most of these GRBs' color measurements are not on sufficient timescales to distinguish between r-process-free and low ψ mix values.On the other hand, the well-sampled color measurements of GRB 190829A provide a strong case for no r-process enrichment at the level of the models used in this section (M rp = 0.03M ⊙ ).Future color measurements with the cadence and long baseline similar to GRB 190829A would allow for more detailed population studies.

Observed Color Diversity Amongst GRB-SNe Sample
Already, we observe diversity in color evolution within our sample of 4 GRB-SNe.In Figure 5 we plot SNdominated color measurements for each GRB that are most closely matched to the V − H filters (chosen due to late-time data availability).While for GRBs 030329, 100316D and 190829A the rest-frame filters are reasonably comparable, the higher redshift of GRB 130427A results in bluer rest-frame filters (≈ B − J) for the data plotted.Using the Barnes & Metzger (2022) models we predict the k-corrections of low-and high-redshift restframe bands produce color differences up to ≈ 0.6 mag (up to ≈ 0.3 mag for the three lower-redshift GRBs).
Figure 5 highlights the diversity within the GRB-SNe sample in terms of their color evolution.At 20 ≲ δt ≲ 60 days GRB 100316D's SN is reddening rapidly, GRB 030329's SN is more slowly reddening, and GRB 190829A's SN is becoming bluer (Figure 5).This The approximate rest-frame V − H color evolution of SN-dominated observations of GRB 030329 (yellow), GRB 100316D (blue), GRB 130427A (red), and GRB 190829A (green), corrected for Galactic and local extinction.Ground-based observations are shown with squares, while observations obtained with HST are represented with circles, and extend the majority of the color evolution curves significantly.Against these observations we plot color evolution models for a moderate-mass SN Ic-BL enriched with 0.03 M⊙ r-process material (greyscaled; color gradient corresponding to ψmix) or r-process-free (purple; Barnes & Metzger 2022).Broadly, the observations are consistent with both models free of r-process material and those with low mixing fractions.There is ≳ 1 magnitude of diversity in the color evolution of the bursts in our sample between 20 ≲ δt ≲ 80 days.
behavior may be explained by differences in M ej or V ej , or perhaps M rp or ψ mix between the GRBs.
In Figure 5 we also compare the GRB-SNe to the large SNe Ic-BL compilation of Anand et al. (2024).This sample includes color measurements for 25 nearby SNe Ic-BL, none of which show strong evidence for r-process enrichment when fit to models.Observations from this sample are corrected for Galactic extinction but not local extinction, though none of their SNe spectra indicate strong local dust.In general, the color evolution of the GRB-SNe in our sample is consistent with those of the Anand et al. (2024) sample, indicating there is no difference in V − H/R − H color evolution and observed diversity between SNe Ic-BL discovered with and without GRB counterparts.This observed trend contrasts with theoretical expectations that GRB-SNe are more likely to produce larger M rp and ψ mix (and thus more pronounced red colors; Siegel et al. 2019;Barnes & Duffell 2023) due to their high angular momentum and accretion disk sizes.Rather, we do not find evidence for differences in color evolution between GRB-SNe and SNe Ic-BL observed without jets.

Quantitative Constraints on M rp and ψ mix
We next consider a larger range of M rp values and constrain M rp and ψ mix for each GRB-SN.To do this, we employ large grids of the best-fit models that are parameterized by combinations of M rp and ψ mix (Barnes & Metzger 2022).We employ the same best-fit parameters for GRBs 030329, 100316D, 130427A, and 190829A as in Section 4 (e.g., Figure 3) and hold these constant across the grid.The grids are linearly spaced in ψ mix between 0.1 and 0.9, and at fixed values, M rp = 0.01, 0.03, 0.08 and 0.15 M ⊙ (Barnes & Metzger 2022).We consider only observations taken between 20 < δt < 200 × (1 + z GRB ) days as the color differences on early timescales are negligible between models.For each model, we determine a χ 2 value using the 1σ color error.As in Section 4, we account for bandpass differences between the models and observations with an additional 0.1 mag error.Finally, we convert the χ 2 value to a p-value using the scipy.chi2.sffunction.
In Figure 6 we show the M rp -ψ mix parameter space.We color code and label each cell according to the pvalue of the model parameterized by the corresponding M rp -ψ mix pair.Black space indicates that models for that M rp -ψ mix pair are not consistent with observed colors (p < 0.02), and are ruled out in our analysis.White space indicates that no models were created for those coordinates.
For GRB 030329, consistent models have ψ mix ≤ 0.3.Our analysis does not provide strong constraints on M rp .For GRB 100316D a wide portion of the parameter space is consistent with observations, although in general lower values of both ψ mix and M rp are favored.For GRB 130427A we favor lower values for both parameters, though multiple combinations are consistent with observations.M rp = 0.15M ⊙ is consistent only in the case where ψ mix = 0.1.Across the GRBs in our sample that are consistent with enriched models, observations disfavor high values of ψ mix (except in some cases of M rp = 0.01M ⊙ for GRB 100316D).We do not find strong constraints on M rp across ψ mix values from our observations.Finally, for GRB 190829A none of the parameter combinations are consistent with observations, in line with our finding in Section 4 that enriched models are, in general, bluer than the observed colors.As GRB 190829A has a high line-of-sight dust extinction that has not been measured with afterglow spectroscopy (Section 3.1), we re-run our analysis (including finding a new best-fit model) for our data corrected for a lower extinction value of E(B − V ) loc = 0.64.Though the lower extinction value produces redder overall colors, it still does not provide any models for the grid of ψ mix − M rp combinations with p > 0.02.

Connection Between ψ mix , γ-ray and SN Properties
We consider any potential connections between the γray, SN and r-process enrichment properties and compare them against theoretical predictions from the literature.From 2D hydrodynamical simulations, Barnes & Duffell (2023) propose that GRBs with the longest rest-frame γ-ray durations will be accompanied by SNe with higher ψ mix values, as both properties are correlated with a long-lived, more massive disk wind.In addition, they find that a higher initial SN explosion energy is associated with lower values of ψ mix .
In Section 4, we concluded that, within the sample, GRB 100316D observations are consistent with the highest mixing values (up to ψ mix ≈ 0.7) while the observed colors of GRBs 030329 and 130427A prefer ψ mix ≲ 0.3 (Figure 6).Comparing the results of these events to the T 90 rest-frame γ-ray durations listed in Table 1 and literature values of the SN explosion energy (e.g., Mazzali et al. 2003;Cano et al. 2017a,b;Hu et al. 2021), we explore trends in our sample.GRB 100316D (for which observations are consistent with ψ mix ≈ 0.4) has the longest rest-frame T 90 γ-ray duration (a lower limit of 261 s as seen by Swift) and the lowest SN explosion energy (Cano et al. 2017a) in our sample.We note that this GRB belongs to the low-luminosity class, which may be indicative of a less energetic central engine or shock breakout.GRBs 030329 and 130427A (for which we deduce ψ mix ≤ 0.3) have somewhat shorter γ-ray duration (rest-frame 18 and 182 s as seen by HETE-2 and Swift) and higher estimated SN explosion energy (Cano et al. 2017a).Though these are just three examples, they align with the trends predicted by Barnes & Duffell (2023).An expansion of M rp -ψ mix constraints for GRB-SNe is necessary to determine if these trends hold within a statistically significant sample.
The findings of Barnes & Duffell (2023) suggest that ultra-long GRBs (ULGRBs; Levan et al. 2014b) may be the ideal candidates for production of r-process el- ements. 6We consider but ultimately do not include several ULGRBs in our sample due to their high redshifts (z ≳ 0.6; e.g., Levan et al. 2014b;Greiner et al. 2015) or an absence of confirmed SN counterpart (e.g., GRB 130925A; Evans et al. 2014;Piro et al. 2014).However, these factors may also limit future rates of nearby ULGRBs suitable for deep NIR follow-up.

Universal r-Process Enrichment Implications
As an exercise, we use our constraints to infer an average r-process mass from GRB-SNe, quantify the contribution of these events to the Universe's r-process budget, and compare our estimate to the Milky Way rprocess enrichment.We caution that these results are highly model dependent and there are high uncertainties in these constraints.For this exercise, we separately consider the derived yields of the three GRBs in our sample consistent with M rp ≈ 0.01 − 0.15M ⊙ for low values of ψ mix .While a range of r-process yields are expected, the median is unknown at this time (Siegel et al. 2019;Barnes & Metzger 2022).Specifically, we employ a yield of M rp = 0.05M ⊙ based on the average of consis-tent M rp values with p > 0.05 (Figure 6).We note that, with the exception of GRB 100316D, all of our events were best-fit with an r-process-free model, though they return p > 0.05 values for some enriched models (e.g., Section 4).
We modify the equation of Rosswog et al. (2018) for events that produce r-process to determine the Milky Way contribution: where R rp is the event rate, mej is the average rprocess ejecta per event, and τ gal is the Milky Way age.We fix τ gal = 1.3 × 10 10 yr for all calculations.
We first consider the case in which only CCSNe associated with long-duration GRBs (LGRBs) produce r-process such that Gpc −3 yr −1 (Ghirlanda & Salvaterra 2022).This rate is modeled using the distributions of observed parameters such as fluence, T 90 , and jet opening angle from Fermi, the Compton Gamma Ray Observatory and Swift (Ghirlanda et al. 2007;Ghirlanda & Salvaterra 2022).
For a yield of M rp = 0.05M ⊙ this results in a total contribution of M r ∼ 4500 +2600 −7700 M ⊙ .This is significantly below the calculated total r-process of the Milky Way of M r,MW ≈ 23 000 M ⊙ (elements of nucleon number A ≥ 69; Hotokezaka et al. 2018

measured from
Europium abundances of local stars from Venn et al. 2004;Battistini & Bensby 2016).We note that other works find the total M r,MW may vary by several thousand M ⊙ (e.g., Bauswein et al. 2014;Rosswog et al. 2018).However, taking M r,MW = 23 000M ⊙ and considering the yields of both GRBs 100316D and 190829A, we determine that GRB-SNe produce 11 − 34% of the Milky Way's r-process abundance.If only GRB-SNe produce r-process amongst CCSNe, either their average M rp yields must be significantly higher than what we derive for GRB 100316D and GRB 190829A, the rates of GRB-SNe are higher than used in our estimate (indeed, the LGRB rate of Ghirlanda & Salvaterra 2022 is based on bright GRBs and the rate of low-luminosity GRBs such as GRB 100316D may be higher) or GRB-SNe are a subdominant r-process production channel compared to BNS mergers.
As our above analysis indicates that GRB-SNe alone cannot account for M r,MW , we consider the case in which collapsars, identified by SN Ic-BL alone, synthesize r-process elements.From the ZTF Bright Transient Survey, the total CCSNe rate is (10.1 +5.0 −3.5 ) × 10 4 Gpc −3 yr −1 of which SNe Ic-BL represent ≈ 4% (Perley et al. 2020).Combined, this results in a rate R rp = R Ic−BL (z = 0) = 3 030 Gpc −3 yr −1 .Combined with our yield estimate of M rp = 0.05M ⊙ , this gives M r ∼ 170 000 M ⊙ which significantly overpredicts Milky Way r-process enrichment.Thus, either only a fraction of SNe Ic-BL without GRBs produce r-process, or they are not a significant contributor, deduced from the significant mismatch in total r-process mass.This finding is consistent with the analysis of Anand et al. (2024).
Our above calculations are dependent on model assumptions of r-process observables in GRB-SNe.Notably, the abundance pattern of elements produced by CCSNe are highly uncertain which likely affects the observed spectral energy distribution (SED) and, thus, the resulting model colors.On the modeling side, abundances depend on observationally unconstrained parameters such as the black hole accretion rate in the collapsar scenario.Typical CCSNe are not expected to produce the heaviest r-process elements (but may produce up to the first peak; e.g., Wang & Burrows 2023) while collapsars may produce up to the third peak elements under high accretion rates or strong magnetic fields (but may still struggle to produce actinides; Siegel et al. 2019).We note that the models used in our analysis assume a r-process SED based on the kilonova AT 2017gfo, which likely produced up to and beyond the lanthanide elements (A ≥ 140; e.g., Kasen et al. 2017).Future models that account for mixing of only lighter r-process elements combined with JWST NIR spectroscopy (see Section 5.3) may help to identify and delineate the abundance pattern.Notably, unlike BNS mergers, long GRBs are localized to galaxies with young ages and high specific star formation rates.These locations lend themselves to enriching early generations of stars, especially compared to BNS mergers which more frequently occur on the outskirts of galaxies (Fong et al. 2022;Mandhai et al. 2022;O'Connor et al. 2022;van de Voort et al. 2022).
Our above derived rates are ballpark estimates and would benefit from additional modeling of r-process enriched collapsar light curves.Further modeling of rprocess-enriched MR SNe light curves could also provide benchmarks for assessing heavy element production through a second mechanism.Similar to BNS mergers and kilonovae (e.g., Bauswein et al. 2014;Shen et al. 2015;Rosswog et al. 2018;Hotokezaka et al. 2018 and references therein), our deduced M r values will be improved with more constrained rate estimates and characterization of the observed diversity of GRB-SNe rprocess yields.

Future Observations
The observational sample presented in this work was not fine-tuned for comparison to the recent models of Barnes & Metzger (2022).Here, we consider if and how future observational strategies of GRB-SNe can be designed to best observe and constrain M rp and ψ mix .In Figure 7 we consider the observability of enriched models consistent with our analysis in Section 4 and 2.6.We plot GRBs 030329 (M rp = 0.03M ⊙ , ψ mix = 0.1), 100316D (M rp = 0.08M ⊙ , ψ mix = 0.3) and 130427A (M rp = 0.01M ⊙ , ψ mix = 0.3) in both luminosity and color space, scaled to z = 0.2.We also plot our r-band afterglow extrapolations of GRBs 030329, 100316D and 190829A shifted to z = 0.2 (spanning m r ≈ 21 − 26 AB mag at δt ≈ 1 day; Section 3.2 and Figure 4) to understand how significantly this component will contaminate future low-redshift GRB-SNe at all times.We do not include GRB 130427A's afterglow in the range as it was known to be superlative in its brightness (e.g., Perley et al. 2014;Laskar et al. 2013).At peak brightness, many nearby GRB afterglows and SNe are observable with 1-2m-class telescopes.However, on the timescales that r-process color signatures may emerge (δt ≈ 40 − 200 days), the expected GRB-SNe brightness spans m = 22 − 27 AB mag, thus requiring large-aperture ground-based telescopes or sensitive space facilities.In considering the sensitivity of large-aperture ground-based telescopes, we account for the fact that most GRB-SNe are embedded in their host Figure 7. Top: expected R-and H-band light curves of a SN Ic-BL that is r-process-free (dotted lines) and enriched with rprocess material (solid lines) set at z = 0.2 (Barnes & Metzger 2022).We show the expected afterglow brightness at z = 0.2 and the depths of large ground and space-based facilities at which ≈ 5σ detections should be possible (y-coordinate chose arbitrarily; more details in Section 5.3).We show three models based on expected or potential outcomes of our analysis in Section 4: our favored parameters for GRB 190829A (left), our favored parameters for GRB 100316D (middle) and an observationally optimistic model that is consistent with observations of GRB 130427A (right).In the bottom panel, we plot R − H over time for corresponding models, shading in grey the timespans at which the afterglow may dominate and in purple the timescales at which HST is necessary.To distinguish and monitor the color evolution of GRB 190829A-like and GRB 100316D-like models from the r-process-free case, space-based observations on the timescales of ∼ 60 days are required.A SN Ic-BL like the optimistic GRB 130427A-like case would be distinguishable by large-aperture ground-based observatories.
galaxies, thus requiring image subtraction which reduces the source's signal-to-noise.
From Figure 7 we show that for a typical afterglow, we do not expect significant contamination on the timescales of the models (δt ≲ 300 days).Our models demonstrates that, to observe similar events at z = 0.2, the SN color evolution could be monitored at peak and distinguished using large-aperture ground-based telescopes, but HST sensitivity is necessary to capture the full color evolution and distinguish between mixing fractions.Overall, we find a wide diversity in expected luminosity and color evolution, but find that both ground-and space-based facilities can play an important role.Looking at upcoming facilities, serendipitous or target-of-opportunity observations with the Nancy Grace Roman Telescope can provide complementary near-IR colors, while JWST is poised to cover a larger dynamic range in optical-near-IR colors and better constrain models.As discussed in Section 4, an observational cadence similar to that of GRB 190829A is critical to constraining the r-process enrichment.
Beyond photometric searches, spectroscopy of GRB-SNe may definitively identify or place deep limits on observed emission lines from transitions of r-process elements.Spectroscopic observations of kilonovae have identified absorption and emission lines from individual elements, most notably Sr II (λ10500; e.g., Watson et al. 2019) and [Te III] (λ21500; e.g., Gillanders et al. 2023;Hotokezaka et al. 2023;Levan et al. 2023a).Caution should be taken when interpreting individual lanthanide line features in GRB-SNe, most notably the Sr II feature as it is coincident with known He lines (Tarumi et al. 2023).Further lines may be identified by future highresolution optical-NIR spectroscopy of kilonovae, most notably by JWST.Future late-time spectroscopy during the nebular phase of GRB-SNe can then be used to measure or place deep limits on the ejecta mass produced of each element.

CONCLUSIONS
We have presented optical-NIR observations of four GRB-SNe extending out to δt = 588 days and search for signs of r-process enrichment in the context of the Barnes & Metzger (2022) models.Our dataset primarily consists of newly-analyzed observations from HST, with additional data from the VLT and the MMT.This analysis provides a template for considerations (e.g., afterglow, local extinction contamination) that may affect observational designs aimed at obtaining future color measurements of GRB-SNe.Moreover, combined with literature data, these observations can be used in comparison to future, enhanced models of r-process-enriched GRB-SNe, including those that account for viewing angle.Our analysis presents the first GRB-SNe sample of color measurements extending to the NIR at late times and allows us to make constraints on M rp and ψ mix for individual objects.Based on our analysis we make the following conclusions: • Comparing our observed colors to those from the models of Barnes & Metzger (2022), we find that GRB 190829A is consistent with no rprocess enrichment, while GRBs 030329, 100316D and 130427A are consistent with models for some enrichment, especially at low mixing fractions (ψ mix ≲ 0.5).However, we note these latter three bursts are also consistent with models for no r-process.In nearly all cases, we find that models parameterized by high mixing fractions (ψ mix ≳ 0.3) are inconsistent with observations.We caution that our analysis is based on a fiducial set of semi-analytic models (Barnes & Metzger 2022) that is unable to account for ≳ 30% of our observations.Future modeling will improve on these conclusions.
• We do not observe differences in R − H GRB-SNe color measurements compared to those of SN Ic-BL without GRBs (Anand et al. 2024), contrasting with predictions that the presence of a jet and a pole-on viewing angle result in larger observed ψ mix (Barnes & Duffell 2023).
• We calculate the r-process enrichment from GRB-SNe and compare it to the observed Milky Way value from the literature.Considering our derived yields for GRBs 100315D and 190829A and the total Milky Way yield of Hotokezaka et al. (2018), GRB-SNe may produce 11 − 33% of the Milky Way's r-process abundance.However, this fraction is highly uncertain, and will be improved with further understanding of r-process yields from GRB-SNe and rate estimates that consider lowluminosity GRBs.
• Building a statistically significant sample of inferences on ejecta mass from GRB-SNe requires large-aperture ground-based and space-based telescopes to monitor events at δt ≳ 70 days.The cadence and long baseline of GRB 190829A observations is the archetype for future GRB-SN monitoring.
Despite the landmark discovery of r-process nucleosynthesis in a neutron star merger, it remains an open question as to whether other heavy element formation channels exist (e.g., Rosswog & Korobkin 2022).The rates of BNS mergers (e.g., The LIGO Scientific Collaboration et al. 2021;Mandel & Broekgaarden 2022;Rouco Escorial et al. 2022) as well as the mass and composition yields of kilonovae remain highly uncertain but likely vary widely (e.g., Gompertz et al. 2018;Metzger 2019;Kawaguchi et al. 2020;Rastinejad et al. 2021).The uncertainties in cosmological and Galactic BNS merger rates and their heavy element yields leave room for the existence of a second r-process formation channel (e.g., Holmbeck & Andrews 2023).At the same time, emerging observations of r-process-enriched metal-poor stellar populations indicate a source of heavy elements closely tracing star formation (e.g., Ji & Frebel 2018;Naidu et al. 2022;Ji et al. 2023;Simon et al. 2023;Kirby et al. 2023).
Moving forward, to identify or place deep constraints on r-process in rare classes of CCSNe requires significant effort and resources on both the theoretical and observational end.The development of observational predictions for r-process-enriched collapsar and MR SNe (or newly developed theories) is critical to establishing these sources as sites from photometric color measurements alone.Further late-time color and/or spectroscopic observations of GRB-SNe, SNe Ic-BL and, potentially, superluminous SNe (e.g., Reichert et al. 2023) will provide additional measurements or upper limits of their r-process yields.Though spectroscopic observations, especially with JWST, are critical to definitively establishing these events as sites of r-process element pro-duction, color measurements are possible for a greater volume of GRB-SNe, improving our understanding of the distribution of ejecta masses.Dedicated programs on space-based facilities such as HST and JWST are necessary for this late-time NIR follow-up.With unprecedented sensitivity in the NIR bands, the upcoming Nancy Grace Roman Space Telescope will be a critical facility for this field.Finally, these studies of some of the most promising candidates for r-process enrichment are not possible without the continued and ensured detection of well-localized GRBs by satellites such as Swift and its successors.

Figure 1 .
Figure1.Example HST images of the fields of the four GRBs in our sample in which the SN is detected (left) and the late-time templates (right).In each panel we show the position of the SN with pink crosshairs and note the time since GRB detection.We do not expect the SN or afterglow to be significantly contributing to any of our templates, with the potential exception of GRB 130427A in F160W (see Section 2.6).
Correcting for Line-of-Sight Dust Extinction

Figure 2 .
Figure 2. All data used in this work for GRBs 030329, 100316D, 130427A and 190829A in apparent magnitude (mAB) versus observed time after the GRB triggers (δt obs ), corrected for Milky Way and local dust extinction.Ground-based observations are shown with squares while HST observations are represented with circles.In each panel we show the extrapolated r-band (and H-band for GRB 130427A) afterglow decay with a dashed line (further described in Section 3.2) and mark observations that are dominated or significantly (> 0.3 mag difference) contaminated by afterglow flux with open symbols.For instance, we expect that our HST observations of GRB 130427A observed beyond 200 days are contaminated or dominated by afterglow flux.
4. EXPLORING R-PROCESS ENRICHMENT IN OUR OBSERVATIONAL SAMPLE 4.1.Direct Comparison of the Sample to Models

Figure 4 .
Figure 4. Assorted colors of the four GRB-SNe in our sample versus rest-frame time compared to best-fit models from Barnes & Metzger (2022) (Figure3and Section 4) with fixed r-process mass (Mrp = 0.03) and varying values of ψmix (lines).We also plot models without r-process enrichment in purple.Afterglow-dominated observations are shown as open symbols, and measurements dominated by SN emission are filled symbols.Though color measurements between filters vary, none of our late-time measurements favor high values of ψmix.Many measurements, particularly those that extend to late times, are unable to distinguish between models with ψmix ≲ 0.2 and those with no r-process.GRB 190829A is, in the best-sampled V − H filters, bluer than the unenriched models.
Figure 5.The approximate rest-frame V − H color evolution of SN-dominated observations of GRB 030329 (yellow), GRB 100316D (blue), GRB 130427A (red), and GRB 190829A (green), corrected for Galactic and local extinction.Ground-based observations are shown with squares, while observations obtained with HST are represented with circles, and extend the majority of the color evolution curves significantly.Against these observations we plot color evolution models for a moderate-mass SN Ic-BL enriched with 0.03 M⊙ r-process material (greyscaled; color gradient corresponding to ψmix) or r-process-free (purple;Barnes & Metzger 2022).Broadly, the observations are consistent with both models free of r-process material and those with low mixing fractions.There is ≳ 1 magnitude of diversity in the color evolution of the bursts in our sample between 20 ≲ δt ≲ 80 days.

Figure 6 .
Figure 6.Comparison of our datasets to large grids of models for a range of r-process masses and mixing fractions.Each cell in the grid represents a model from Barnes & Metzger (2022) parameterized by the corresponding Mrp and ψmix.Using the total color measurements per GRB meeting our criteria (Section 4.3), we evaluate a χ 2 value and report one-sided p-values for each model in the grid.The color shade and labeled number of each cell corresponds to the p-value of the model parameterized by that combination of Mrp and ψmix.Models whose p-values are less than 0.02 are shaded black.Together, our observations favor low ψmix in nearly all cases but are consistent with the full range of Mrp.We do not show a grid for GRB 190829A as no models have p-values > 0.02.

Table 2 .
Sample of Observations

Table 2 continued
Table 2 (continued) Indicates observation is dominated or significantly affected by afterglow emission.( †) Observations are not corrected for Galactic nor local extinction.