Introduction

For decades, the formation sites for many of the Universe's heaviest elements have 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. 1991Cowan et al. , 2021;;Woosley et al. 1994;Hoffman et al. 1997;Thielemann et al. 2011), 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 postmerger 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 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;Soares-Santos et al. 2017;Tanvir 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 rprocess-enhanced metal-poor stars in the Milky Way, some ultrafaint 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 gamma-ray burst (GRB) host galaxy offsets, stellar population ages, and star formation histories (e.g., Fong et al. 2022;Nugent et al. 2022Nugent et al. , 2024;;O 'Connor et al. 2022;Zevin et al. 2022) 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 r-process 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 proto-neutron star remnant (e.g., Woosley et al. 1994). However, high event rates and observations of ordinary CCSNe have led to a general consensus that these events do not significantly contribute to the Universe's r-process budget (Wallner et al. 2015(Wallner et al. , 2021;;Macias & Ramirez-Ruiz 2018; although see Tsujimoto & Shigeyama 2001 who find that SN 1987A'sBa/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 magnetorotational (MR) supernovae (SNe; also termed jetdriven SNe or MR 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 neutron-rich 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 &u d -D o u l a2018). At present, no models for observational signatures of r-process produced through the MR supernova (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 rprocess elements. Using magnetohydrodynamic 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) r-process material (Miller et al. 2020;Fujibayashi et al. 2023;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), r-process 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 r-process enhancement in SNe associated with collapsars. Searches for radio flares following long-duration GRBs (LGRBs; 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., Ho et al. 2020;Izzo et al. 2020;Jin et al. 2021), and there are no signs of reddening in the SN light curve. Anand et al. (2024) performed a comprehensive search for rprocess 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 (ZTF; 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 fitb yr-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. (2024) 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.

Observations

Sample Selection and Data Description

As a starting point, we utilize the comprehensive GRB-SNe compilation of Dainotti et al. (2022), which includes 58 LGRBs 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 r-process 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;Yangetal.2022Yangetal. , 2024;;Levan et al. 2024;Gillanders 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), and GRB 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),s p a r s el a t e -t i m e sampling of observations, and afterglow contamination in early HST epochs (Levan et al. 2023). 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, by its GRB name.

In total, we collected 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 Mikulski Archive for Space Telescopes (MAST) archivefoot_1 and the Hubble Legacy Archive (HLA). 29 The observations span 8 < δt < 969 days post-burst (including template observations) and seven HST filters (Table 1). We also collect and reduce Very Large Telescope (VLT)/X-shooter (acquisition camera), VLT/HAWK-I, and

Table 1 GRB Sample and Properties GRB SN T 90,rest a zA V,MW E(B -V ) loc Ref. HST Observed Filters b Program ID(s) (s)( mag)( mag) 030329 2003dh 18 0.1685 0.069 0.04 1 (1) F606W, F814W, F110W, F160W 9405 100316D c 2010bh 261 0.0591 0.319 0.14 (2) F555W, F814W, F125W, F160W 11709, 12323 130427A 2013cq 182 0.3399 0.055 0.05 (3) F606W, F160W 13110, 13117, 13230,13951 190829A 2019oyw 49 0.0785 0.133 1.04 (4) F606W, F110W, F140W, F160W 15089, 15510, 16042, 16320 Notes. Milky Way extinction values are taken from Schlafly & Finkbeiner (2011). a T 90 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. b Filters selected based on data availability on timescales relevant for observing r-process-enriched components (contemporaneous colors observed at δt  30 days). c Low-luminosity GRB. References: (1) Matheson et al. (2003); (2) Bufano et al. (2012); (3) Levan et al. (2014a); (4) Chand et al. (2020). 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.I nS e c t i o n3.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 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 significantly contribute to any of our templates, with the potential exception of GRB 130427A in F160W (see Section 2.6). is to compare the optical-NIR color evolution of the GRB-SNe in our sample to relevant models, we collect only hostsubtracted photometry in the rRiIJHK bands of GRBs 030329 (Matheson et al. 2003), 100316D (Olivares et al. 2012), 130427A (Perley et al. 2014), and190829A (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 data set at early times.

GRB 030329 GRB 030329 was discovered by the High Energy Transient

Explorer II (HETE-II) at 11:14:14.67 UT on 2003 March 29, with a duration of 21 s (Vanderspek et al. 2003). A bright optical afterglow counterpart was quickly localized. The burst's redshift of z = 0.1685 was identified through afterglow spectroscopy (Greiner et al. 2003). Subsequent spectroscopic and photometric observations of the counterpart revealed evidence of an SN Ic-BL (Hjorth et al. 2003;Kawabata et al. 2003;Matheson et al. 2003;Mazzali et al. 2003;Stanek et al. 2003;Lipkin et al. 2004).

GRB 030329 was observed with ACS/WFC in the F606W and F814W filters at several nearly contemporaneous epochs (within ≈24 hr) over 17  δt  228 days, and at δt = 428 days in F606W only (Program 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 2003 November 12 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 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 zero-points 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 Handbook.foot_3

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 2010 March 16. 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 an 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 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 corrections. 31 We note that the SN in the F125W and F160W images on 2010 April 4 (δ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).A tδt  24 days, we find our values are fainter by ∼0.3-1.2 mag compared to their values. This difference can be ascribed to our process of image subtraction in comparison to their direct photometry on the host-embedded 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.42 UT. 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 multiwavelength 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;Levan et al. 2014a;Melandri et al. 2014).

GRB 130427A was observed with ACS/WFC and WFC3/ UVIS in the F606W filter and WFC3/IR in the F160W filter for seven 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 did not see a 3σ significant residual and discarded the image. For all other epochs, we performed 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 2019 August 29 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 r band 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 pipeline. 32 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 Sloan Digital Sky Survey r (Alam et al. 2015) or Two Micron All Sky Survey 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 star-forming regions (e.g., Blanchard & Berger 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 fita 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 and 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 late-time 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.

Determination of the Intrinsic SN Colors

Correcting for Line-of-sight Dust Extinction

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),a n dw h e r e 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;Xuetal.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 (Matheson et al. 2003). Literature values for the 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 month-long 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. 2023). In the majority of cases, we expect that follow-up on timescales beyond approximately 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 late-time 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.45 measured from the multiwavelength afterglow (Perley et al. 2014). We note that Maselli et al. (2014) observe a jet break in the X-ray and optical light curves at δt ≈ 10 hr 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 data set 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 yr). Finally, for GRB 190829A, we utilize the r-band temporal index of α = -1.45 from 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 r-band afterglow contributions in Figure 2 as gray 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 significant 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).

SN 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 CCSNe remain uncertain, though new observations with JWST are beginning to constrain the dust abundance in Type II SNe on approximately 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 preexisting 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 yr 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 the reprocessing of emission in an eject-CSM shock interaction (Figure 2).

Exploring r-Process Enrichment in Our Observational Sample

Direct Comparison of the Sample to Models

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 of ∼10 -3 -10 -2 M e s -1 ; 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 (non-r-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 ),massofr-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 r-process material being homogeneously mixed all the way to the ejecta surface), resulting in a more pronounced and prolonged red color (Figure 4).T h eR -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 that we are left with only M rp and ψ mix . To determine the values of M ej , β ej , and M 56Ni , which 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. 2017b;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. 2017b).

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 bestfit model is enriched with M rp = 0.01 M e and highly mixed (ψ mix = 0.9; dotted lines in Figure 3). The values of M ej , β ej , and M 56Ni 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 r-process enrichment. We use the models described above and shown in Figure 3 enriched with M rp = 0.03 M e , 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 3 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 UBVRIJHK 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).W e explore how the GRB-SN color evolution compares to models as a function of mixing fraction, and also compare to an rprocess 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 rest-frame V -R, where model color differences are small between the rprocess-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 consistent with both the 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 rest-frame ground-based B -R and B -H data are not appropriate to distinguish between strong mixing or r-process-free 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 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 M ej , β ej , and M 56Ni 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 M ej , β ej , and M 56Ni for each GRB-SNe throughout our analysis.

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.03 M e (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.0 3M e ). Future color measurements with the cadence and long baseline similar to GRB 190829A would allow for more detailed population studies.

Observed Color Diversity among GRB-SNe Sample

Already, we observe diversity in color evolution within our sample of four GRB-SNe. In Figure 5, we plot SN-dominated 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 rest-frame 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ʼsS N is more slowly reddening, and GRB 190829A's SN is becoming bluer (Figure 5). This 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 fitt o 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 e (Barnes & Metzger 2022).W ec o n s i d e ro n l y 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.sf function.

In Figure 6, we show the M rp -ψ mix parameter space. We color code and label each cell according to the p-value 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 .F o r 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 a r ec o n s i s t e n tw i t ho b s e r v a t i o n s .M rp = 0.15 M e 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.01 M e for GRB 100316D).W ed on o tfind 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 rerun 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 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 e r-process material (grayscaled; color gradient corresponding to ψ mix ) or rprocess-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 mag of diversity in the color evolution of the bursts in our sample between 20  δt  80 days.

provide any models for the grid of ψ mix -M rp combinations with p > 0.02.

Discussion

Connection between ψ mix , Gamma-Ray, and SN Properties

We consider any potential connections between the gammaray, 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 gamma-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 gamma-ray durations listed in Table 1 and literature values of the SN explosion energy (e.g., Mazzali et al. 2003;Cano et al. 2017aCano et al. , 2017b;;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 gamma-ray duration (a lower limit of 261 s as seen by Swift) and the lowest SN explosion energy (Cano et al. 2017b) 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 gamma-ray duration (rest frame 18 and 182 s as seen by HETE-2 and Swift) and higher estimated SN explosion energy (Cano et al. 2017b). 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 ultralong GRBs (ULGRBs; Levan et al. 2014b) may be the ideal candidates for the production of r-process elements. 33 We 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 rprocess 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 r-process 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.15 M e for low values of ψ mix . W h i l ear a n g eo fr-process yields is expected, the median is unknown at this time (Siegel et al. 2019;Barnes & Metzger 2022). Specifically, we employ a yield of M rp = 0.05 M e b a s e do nt h e average of consistent M rp values with p > 0.05 (Figure 6).W e note that, with the exception of GRB 100316D, all of our events were best fitw i t ha nr-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:

  t ~´--MM m M 17 000 500 Gpc yr 0.03 1.3 10 yr 1 r rp 31 ej gal 10



where rp  is the event rate, mej is the average r-process 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 LGRBs produce rprocess such that

79 33 57 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.05 M e, this results in a total contribution of

 ~-+ MM 4500 r 7700 2600

. This is significantly below the calculated total r-process of the Milky Way of M r, MW ≈ 23,000 M e (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 solar mass (e.g., Bauswein et al. 2014;Rosswog et al. 2018). However, taking M r,MW = 23,000 M e 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 among 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 ( . Combined with our yield estimate of M rp = 0.05 M e , this gives M r ∼ 170,000 M e, 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 is 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 2024), 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 an 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 the 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 r-process-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;Hotokezaka et al. 2018;Rosswog 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 finetuned 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,w ec o n s i d e rt h e observability of enriched models consistent with our analysis in Sections 4 and 2.6. We plot GRBs 030329 (M rp = 0.03 M e , ψ mix = 0.1), 100316D (M rp = 0.08 M e , ψ mix = 0.3),a n d 130427A (M rp = 0.01 M e , ψ 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 ≈ 1d a y ; 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., Laskar et al. 2013;Perley et al. 2014). 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 largeaperture 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 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 demonstrate 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 groundand 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 NIR colors, while JWST is poised to cover a larger dynamic range in optical-NIR 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., Levan et al. 2024;Gillanders et al. 2023;Hotokezaka et al. 2023). 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 high-resolution 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 by each element.

Conclusions

We have presented optical-NIR observations of four GRB-SNe extending out to δt = 588 days and searched for signs of rprocess enrichment in the context of the Barnes & Metzger (2022) models. Our data set 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: ). 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 chosen 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 gray the time spans 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. An SN Ic-BL like the optimistic GRB 130427A-like case would be distinguishable by large-aperture ground-based observatories.

1. Comparing our observed colors to those from the models of Barnes & Metzger (2022),wefind that GRB 190829A is consistent with no r-process 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 rprocess. 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. 2. 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). 3. 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 low-luminosity GRBs. 4. 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 are 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 2024). The rates of BNS mergers (e.g., The LIGO Scientific Collaboration et al. 2021;Mandel & Broekgaarden 2022;Rouco Escorial et al. 2023) 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 2024). 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;Kirby et al. 2023;Simon et al. 2023).

Moving forward, to identify or place deep constraints on the 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-processenriched 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 rprocess yields. Though spectroscopic observations, especially with JWST, are critical to definitively establishing these events as sites of r-process element production, 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.