A Radio-selected Population of Dark, Long Gamma-ray Bursts: Comparison to the Long Gamma-ray Burst Population and Implications for Host Dust Distributions

We present cm-band and mm-band afterglow observations of five long-duration $\gamma$-ray bursts (GRBs; GRB 130131A, 130420B, 130609A, 131229A, 140713A) with dust-obscured optical afterglow emission, known as"dark"GRBs. We detect the radio afterglow of two of the dark GRBs (GRB 130131A and 140713A), along with a tentative detection of a third (GRB 131229A) with the Karl G. Jansky Very Large Array (VLA). Supplemented by three additional VLA-detected dark GRBs from the literature, we present uniform modeling of their broadband afterglows. We derive high line-of-sight dust extinctions of $A_{V, \rm GRB} \gtrsim 2.2 - 10.6~{\rm mag}$. Additionally, we model the host galaxies of the six bursts in our sample, and derive host galaxy dust extinctions of $A_{V, \rm Host} \approx 0.3-4.7~{\rm mag}$. Across all tested $\gamma$-ray (fluence and duration) and afterglow properties (energy scales, geometries and circumburst densities), we find dark GRBs to be representative of more typical unobscured long GRBs, except in fluence, for which observational biases and inconsistent classification may influence the dark GRB distribution. Additionally, we find that $A_{V, \rm GRB}$ is not related to a uniform distribution of dust throughout the host, nor to the extremely local environment of the burst, indicating that a larger scale patchy dust distribution is the cause of the high line-of-sight extinction. Since radio observations are invaluable to revealing heavily dust-obscured GRBs, we make predictions for the detection of radio emission from host star formation with the next generation VLA.


INTRODUCTION
Long duration γ-ray bursts (GRBs), the most luminous, energetic transients in the universe (e.g. Racusin et al. 2008), are associated with the death of massive stars (MacFadyen & Woosley 1999) and Type Ic supernovae (e.g. Hjorth et al. 2003;Woosley & Bloom 2006). The relativistic shocks produced by these bursts interact with the surrounding medium, leading to the pro- * NASA Einstein Fellow duction of broad-band synchrotron emission from radio to X-rays, known as the "afterglow".
A subset of long GRBs have suppressed optical emission, earning them the moniker of optically "dark" GRBs, with the first documented dark GRB being GRB 970828 (Groot et al. 1998). It is estimated that ∼ 10%-50% of the long GRB population are among the dark GRB class (Jakobsson et al. 2004;Cenko et al. 2009;Fynbo et al. 2009;Greiner et al. 2011;Melandri et al. 2012;. One of the prominent causes of the darkness is attributed to dust extinction arXiv:2205.01124v2 [astro-ph.HE] 6 May 2022 along the line-of-sight to the GRB, though other causes, such as intrinsically faint and rapidly fading bursts (e.g. Berger et al. 2002), or Lyα absorption from a high redshift (z 6) origin (Haislip et al. 2006;Salvaterra et al. 2009;Tanvir et al. 2009;Cucchiara et al. 2011) can also produce optically faint or undectectable afterglows. For the dust-extinguished population of dark GRBs, the amount of extinction (A V,GRB ) along the line-of-sight is of interest as it can provide insight on the amount of dust and distribution within their host galaxies Zafar et al. 2011;Krühler et al. 2011;Zauderer et al. 2013c;. Typically, A V,GRB is estimated using simple power law arguments extending from the X-rays to the optical bands (e.g. . However, the addition of a detected radio afterglow can allow for proper afterglow modeling and provide a robust measurement of A V,GRB . Since the first radio afterglow of a GRB was discovered in 1997 (GRB 970508, Frail et al. 1997), radio follow-up has been vital to our general understanding of GRB afterglow behavior, helping to constrain the energetics and environment of GRBs. A fairly comprehensive census of radio follow-up of 304 GRB afterglows found that only ∼ 30% of GRB afterglows have radio detections (Chandra & Frail 2012). However, it was postulated that the low fraction of radio detections is likely attributed to the limited sensitivity of prior generations of radio telescopes (Chandra & Frail 2012;Osborne et al. 2021), and therefore more sensitive radio facilities could increase the fraction of detected radio afterglows of GRBs. The upgrade to the Karl G. Jansky Very Large Array (VLA), which concluded in 2012, increased the sensitivity of the radio array by a factor of ∼ 10 (Perley et al. 2011), and has allowed for well-sampled, multi-band, follow-up of long GRB radio afterglows (including, but not limited to: Chandra & Frail 2012;Laskar et al. 2013c;Zauderer et al. 2013c;Laskar et al. 2014Laskar et al. , 2015Laskar et al. , 2016Alexander et al. 2017;Laskar et al. 2018a,b,c). Radio observations of dust obscured dark GRBs are critical for properly modeling the afterglow of these bursts, especially when the optical afterglow is not detected (e.g. Jakobsson et al. 2005;Castro-Tirado et al. 2007;Rol et al. 2007;Zauderer et al. 2013c;van der Horst et al. 2015;Higgins et al. 2019;Kangas & Fruchter 2021), as the lowfrequency afterglow can constrain burst energetics and environment density through determination of the synchrotron break frequencies (e.g. Sari et al. 1998;Granot & Sari 2002).
In addition to the properties of the afterglow of dark GRBs, there has been much interest in the host galaxies of dark GRBs and how they compare to the typical long GRB population (e.g. Krühler et al. 2011;. The Swift X-ray Telescope (Swift/XRT) typically provides afterglow positions of ∼ 2 (Evans et al. 2009), while the detection of a radio afterglow with the VLA can often provide unambiguous association to a host galaxy via sub-arcsecond precision, especially in the case of optically-faint or non-detected dark GRBs (e.g. Zauderer et al. 2013c). Previous studies based on small numbers have found that the hosts of dust-obscured bursts are overall more massive, more luminous, more star forming, and dustier than other long GRB hosts (i.e. Krühler et al. 2011;. Long GRBs in general are already associated with star forming galaxies (i.e. Djorgovski et al. 1998;Christensen et al. 2004;Japelj et al. 2016;Palmerio et al. 2019), and the obscuration of dark GRBs may point to obscured star formation in their hosts (Blain & Natarajan 2000;Ramirez-Ruiz et al. 2002).
Here, we present the multi-wavelength observations for five dark GRBs, including VLA observations for all five bursts. We present the discovery of the VLA radio afterglow for one burst (GRB 130131A), present a nominal VLA radio detection for another burst (GRB 131229A), present new radio and millimeter detections for an additional event (GRB 140713A), and new upper limits for two events. For our three VLA detected/nomiminally detected dark GRBs, we uniformly model the afterglows and host galaxies, and include uniform models of the afterglow and host galaxies of three VLA detected dark bursts from the literature. We proceed to compare the afterglow and host properties of dark GRBs to the larger population of typical long GRBs that have not been classified as dark. In Section 2 we describe our sample and our criteria for classifying GRBs as dark. In Section 3 we describe our methods for self-consistent modeling of all available broadband afterglow data and we apply our afterglow modeling to the relevant bursts in our sample to extract their local environment, burst energetics and microphysical parameters. In Section 4 we present host galaxy modeling for the bursts with robust host galaxy associations and spectroscopic or photometric data. In Section 5 we compare our dark GRB sample to the broader population of long GRBs in terms of their γ-ray, afterglow, and host properties, and consider the detectability of obscured star formation in long GRB host galaxies. We conclude in Section 6. We present the details of the multi-wavelenth observations and data reduction in Appendix A. In this paper, we employ the ΛCDM cosmological parameters of H 0 = 68 km s −1 Mpc −1 , Ω M = 0.31, Ω Λ = 0.69.
Our primary goal is to uniformly model the afterglows and host galaxies of dark GRBs. We focus on dark GRBs with VLA observations taken after the upgrade (Perley et al. 2011), specifically those taken with our Programs 13A-046,  1 . This sample includes five dark GRBs, two of which have unambiguous radio afterglow detections, and one of which has a tentative radio afterglow detection. We supplement this sample with three other dark GRBs with VLA detections, GRB 110709B ), 111215A van der Horst et al. 2015) 2 , and GRB 160509A (Laskar et al. 2016, classify-ing this latter burst as dark for the first time in Section A.6.3). As 14 dark GRBs have been observed by the VLA since its upgrade (Zauderer & Berger 2012;Zauderer et al. 2013c;Laskar et al. 2013a;Veres et al. 2015;Horesh et al. 2015;Laskar et al. 2016, the VLA Data Archive, This Work), our sample of 8 comprises over half of the known dark GRBs with upgraded VLA observations. The eight 8 dark GRBs featured in this paper are listed in Table 1.
Standard measures of classifying long GRBs as "dark" in the standard synchrotron model (Granot & Sari 2002 and Section 3), involve the optical-to-X-ray spectral index β OX and the X-ray spectral index β X . One such method of classification uses the expected spectral index between the optical and X-ray afterglow fluxes of −p/2 < β OX < (1 − p)/2, where p is the power-law index of the electron energy distribution, F ν ∝ ν β , and the spectral index is dependent on the location of the cooling frequency (ν c ) in relation the optical and X-ray bands. As p > 2 is typically expected (e.g. Sari et al. 1998), the shallowest expected β OX = −0.5, where p = 2. Therefore, if β OX > −0.5, then the shallow spectral slope indicates the optical afterglow has been suppressed and the burst can be considered dark (Jakobsson et al. 2004). An alternate definition of afterglow darkness is β OX > β X + 0.5, corresponding to an optical flux that is even lower than the shallowest possible extrapolation from the X-rays, ν c ≈ ν X , in the synchrotron framework (van der Horst et al. 2009). For the purposes of this paper, we will consider any burst "dark" if they meet the Jakobsson et al. (2004) criterion of β OX > −0.5, though 1 We exclude GRB 130606A as the darkness of this burst is likely attributed to its high redshift of z ≈ 5.91 (Castro-Tirado et al. 2013;Chornock et al. 2013a;Littlejohns et al. 2015) 2 While the upgrade to the VLA completed at the end of 2012, the C-band (4-8 GHz) and K-band (18-26.5 GHz) receivers were upgraded by 2011 (Perley et al. 2011). As the majority of the VLA observations of GRB 110709B and GRB 111215A were taken with these upgraded receivers, we include these bursts in our sample.
we will note whether the bursts in our sample meet the van der Horst et al. (2009) classification criterion as well.
To determine β OX , we interpolate the X-ray light curve to the times of the optical observations, using least squares fits to the X-ray light curves to calculate the X-ray temporal index α X (where F ν ∝ t α ), and calculate β OX , correcting the optical observations for Galactic extinction in the direction of the burst (Schlafly & Finkbeiner 2011). To determine β X , we create timesliced spectra from the Swift tool 3 , which calculates the photon index, Γ X (Evans et al. 2009), from which we derive β X ≡ 1 − Γ X . On a per-burst basis, we exclude any times over which the X-ray light curve exhibits flaring activity super-imposed on the power-law afterglow, as this emission is not likely to originate from the external shock (Burrows et al. 2007;Margutti et al. 2010).
The full details of the X-ray, optical, near-infrared (NIR) and radio observations of the GRB afterglows and host galaxies, as well as their classifications as dark GRBs are presented in Appendix A. We show the fields and afterglow localizations of six dark GRBs in our sample in Figure 1, and their afterglow detectabilities, as well as whether we model their afterglows and host galaxies, are summarized in Table 1.

AFTERGLOW MODELING
We will consider the X-ray, optical, and radio afterglow light curves and spectral energy distributions (SEDs) of the GRBs in our sample in the context of synchrotron emission from the acceleration of electrons from a relativistic blast wave (i.e. Sari et al. 1998;Chevalier & Li 2000;Panaitescu & Kumar 2000;Granot & Sari 2002). These electrons are accelerated to a non-thermal power law distribution, N (γ e ) ∝ γ −p e . The afterglow SEDs can be described by power law segments which connect at three break frequencies (the self absorption frequency, ν a , the characteristic frequency, ν m , and the cooling frequency, ν c ) and the characteristic flux, F ν,m (Granot & Sari 2002). The SED and light curve temporal evolution are dependent on the following parameters: p, the isotropic-equivalent kinetic energy of the burst, E K,iso , the density of the environment 4 , ρ = Ar −k , and the fractional energy density imparted on the electrons, e , and the magnetic field, B .
3 https://www.swift.ac.uk/user objects/docs.php#specform 4 In a constant density interstellar medium (ISM) environment, k = 0 and A = n 0 mp. In a wind environment, k = 2, A = 5 × 10 11 gm cm −1 A * , normalized to a progenitor mass-loss rate of A * = 10 −5 M yr −1 and wind velocity of 1000 km s −1 ; (Chevalier & Li 2000).  In each panel, the blue solid-line circle indicates the XRT afterglow location (90% confidence). When available, the red circle indicates the radio afterglow location (1σ), the green circle denotes the optical or NIR position (1σ), and the dashed blue circle indicates the CXO position. Images have been smoothed with a 3-pixel Gaussian, and the orientation is North up and East to the left. An angular scale is provided in all images, in addition to a physical scale when the redshift is known. For GRBs 130420B and 130609A, we do not identify host galaxies.
In addition to the standard synchrotron frame work laid out in Granot & Sari (2002), we also consider the effects of beaming on the GRB afterglow light curves. Collimation is expected to manifest itself as a jet break, which occurs when the angular size of the relativistic beam of the GRB jet approaches the value of the true opening angle of the jet (θ jet ). At the time that this jet break occurs, t jet , the observed light curve is predicted to steepen achromatically (Rhoads 1999;Sari et al. 1999). Determining t jet from a GRB afterglow can therefore provide us with the opening angle of the jet, which allows for the correction of the energetics for beaming and in turn, we can derive the true kinetic (E K ), γ-ray (E γ ), and total energies of the burst (E tot = E K + E γ ). Additionally, the evolution of the afterglow light curves at δt > t jet (where δt is the time after the Swift/Burst Alert Telescope trigger) is F ν ∝ t −p when ν obs ≥ ν m (Sari et al. 1999), where ν obs is the observing frequency. Thus, identifying a jet break in the afterglow and measuring the post-break light curve slope provides an additional constraint on the value of p.
To properly model the X-ray afterglow, we include the effects of inverse Compton (IC) cooling in our model, which can lower the location of ν c (in comparison to the value predicted by a spherical blast-wave without IC cooling, e.g., Granot & Sari 2002) by a factor of (1+Y ) 2 , where Y is the Compton Y -parameter (for the detailed explanation see: Sari & Esin 2001;Laskar et al. 2015).
Prior to performing a full fit to the broad-band data, we first examine the SEDs and light curves of our afterglow observations and fit them using either a single power law function or a smoothly broken power law, given by where F b is the normalized flux of the break, ν b and t b are the break frequency and time, respectively, s is Note-a Whether or not the afterglow was detected in the optical/NIR or radio. b Whether or not we model the afterglow and host galaxy of the burst in this paper c We find a tentative 6.0 GHz detection for GRB 131229A (see Section A.4.2) d Due to the limited radio data, we present a simple analytical afterglow model for GRB 131229A in Section 3.4, and analytically derive afterglow properties. We model the host galaxy of GRB 131229A in this work. e The X-ray to radio afterglow of GRB 160509A was modeled within our modeling framework in Laskar et al. (2016). We model the host galaxy of GRB 160509A in this work. the break smoothness 5 , and α 1 , α 2 and β 1 , β 2 are the temporal or spectral slopes of the fits, respectively. We use the convention F ν ∝ t α ν β throughout. Where data quality allows, we use Equation 1 for light curves and Equation 2 for SEDs. When necessary, we use temporal power law fits to interpolate the flux to common times for our SED fitting. We use these basic broken power law considerations to place initial constraints on p, t jet , the nature of the burst environment (i.e. ISM vs. wind), and the location of the break frequencies with respect to our observing bands, when possible. We then model the GRB with the Markov Chain Monte Carlo (MCMC) modeling framework laid out in Laskar et al. (2014), including IC effects. We also take into consideration the scattering effects of scintillation, which may cause variability on short timescales at GHz frequencies (Rickett 1990). In the situations for which the burst environment cannot be constrained with initial considerations alone, we model the GRB afterglow with both an ISM and wind environment and choose the model that provides a better statistical fit.
To account for any potential systematic uncertainties in flux calibration for data taken across different facilities, we include an uncertainty floor of 10% on individual measurements prior to modeling. The free parameters for our model are p, E K,iso , n 0 in an ISM environment or A * for a wind environment, e , B , and t jet . We fit each GRB with available broad-band afterglow data with the MCMC model, using 128 walkers for 10,000 steps, and discard the first ∼ 1% of steps as "burn in", where the average likelihood across the chains have yet to reach a stable value.

GRB 110709B
We compile all available radio data of GRB 110709B, along with optical upper limits and the Swift X-ray light curve (Section A.6.1), to model the afterglow. The data consists of a VLA 5.8 GHz light curve spanning δt ≈ 2.1 − 69.7 days, a VLA 21.8 GHz upper limit at δt ≈ 12.1 days, GROND optical/NIR (g r i z JHK s ) non-detections at δt ≈ 0.11 days, and the Swift X-ray afterglow light curve spanning δt ≈ 10 −3 − 10 2 days.
The afterglow of GRB 110709B was previously modeled by Zauderer et al. (2013c), and then by Kangas & Fruchter (2021). Zauderer et al. (2013c) did not include IC effects in their modeling, while Kangas & Fruchter (2021) did include IC effects, though they did not fit for A V,GRB , which is a key parameter of interest for our study. Here, we model GRB 110709B to ensure consistency across our sample and to determine A V,GRB .
We first investigate the nature of the steepning in the X-ray light curve of GRB 110709B. Such steepenings are often explained by either the passage of ν c through the band, or a jet break. For the passage of ν c at δt ≈ Note-The top row for each parameter corresponds to the best fit forward shock value from our MCMC modeling. The bottom row for each parameter corresponds to the summary statistics from the marginalized posterior density functions (medians and 68% credible intervals), except in the case of GRB 131229A (see a ) a Values derived from analytical arguments (see Section 3.4) b A V,GRB value from afterglow model using median values 0.6 days, the expected change in temporal index is ∆α = 0.25, which is too shallow to explain the observed ∆α ≈ 0.7 − 0.9. Thus, we attribute the steepening instead to a jet break. We now use the X-ray spectral index and pre-break light curve to determine where ν X lies in relation to ν c . If ν m < ν X < ν c , the X-ray spectral index implies p = 3.12 ± 0.07. However, we find that the measured value of α X,1 yields p = 1.49 ± 0.04 in a wind environment or p = 2.16 ± 0.04 in an ISM environment, neither of which are consistent with the value derived from β X . On the other hand, for ν m , ν c < ν X , we require p = 2.12 ± 0.07 to match the X-ray spectral index and p = 1.82 ± 0.04 (in both the wind and ISM environment) to match the light curve. The values of p are in agreement to within 3σ, and we therefore conclude that ν m , ν c < ν X and p ≈ 1.8 − 2.1. In this regime, the X-ray observations cannot be used to discriminate between an ISM and wind environment.
The expected temporal decay after a jet break is α = −p, and our post-break X-ray light curve slope of −1.7 α X,2 −1.6 is shallower than the expected −2.1 α −1.8. However, the break time and slope are degenerate with the smoothness of the break, and a later break time is consistent with a steeper post break decline. This suggests that the break in the X-ray light curve is due to a jet break, whose onset occurs at δt 0.6 days. Additionally, if t jet ≈ 0.6 days, then the 5.8 GHz light curve, which does not have any observations prior to t jet , cannot be used to distinguish between the ISM and wind environment, as the behavior of the synchrotron model is the same regardless of environment after t jet .
In conclusion, the X-ray afterglow of GRB 110709B is consistent with ν m , ν c < ν X , p ≈ 1.8 − 2.1, and t jet 0.6 days. Additionally, the X-ray and radio light curves cannot be used to distinguish between the ISM and wind environment. . X-ray, optical, and radio afterglow light curves of GRB 110709B, together with the best-fit forward shock model in an ISM environment (lines). Circles represent detections, and triangles represent 3σ upper limits . Open symbols indicate data that are not included in the fit, and shaded regions represent variability due to scintillation. The inset shows the model Ks-band and Hband light curves (solid lines) as well as the non-extinguished models (dashed lines), with AV,GRB 3.9 mag necessary to be consistent with the upper limits.

MCMC Modeling
As we can not distinguish between the wind and ISM environments using preliminary analytical arguments, we therefore fit the GRB 110709B afterglow data with both a wind and ISM environment, and choose the solution that provides a better statistical fit.
We find the ISM environment model is marginally preferred by the data, with χ 2 /d.o.f. ≈ 386/1289 and likelihood (L) of ≈ 2092, compared to the wind environment best fit model with χ 2 /d.o.f. ≈ 434/1289 and a lower L of ≈ 2063 (see Appendix B.1 for wind model, provided for completeness). We present the best-fit (highest likelihood) ISM model in Figure 2 and list the parameters as well as the summary statistics from the marginalized posterior density functions (medians and 68% credible intervals) in Table 2.
We find that p ≈ 2.02 and t jet ≈ 1.4 days (consistent with the arguments laid out in Section 3.1.1), resulting in θ jet ≈ 1.4 • . The SED remains in the slow cooling phase with a break frequency ordering of ν sa < ν m < ν c for the entirety of the afterglow observations considered. Our fit also confirms that ν c < ν X , as expected. For our best fit ISM model, the GROND optical/NIR limits imply A V,GRB 3.94 mag.
We comment here on a few notable discrepancies between the data and model light curves. The 5.8 GHz model light curve under-predicts the last detection at δt ≈ 40 days by ≈ 11.6σ. This discrepancy is too large to be reconciled with scintillation effects alone. Moreover, the X-ray model light curve under-predicts the data at δt 7.6 days. This result is not unexpected, as we measure a shallower temporal index in the post-break X-ray light curve than would be expected for post jet-break behavior (see Section 3.1.1), and therefore we find an excess within our model. We first consider whether this excess flux in the X-ray afterglow is due to Klein Nishina (KN) effects. KN effects become important when ν m ≤ ν obs (typically the ν obs is the X-ray frequency), where ν m = ν m ( γ m ), and γ m is the critical Lorentz factor (Nakar et al. 2009). Given our best fit parameters, we find that at δt ≈ 12.6 days, ν X ν m ≈ 1.6 × 10 25 Hz. Therefore, we conclude that KN effects are not causing the excess flux in the X-ray afterglow. We next consider whether IC emission is the cause of the observed X-ray excess. We calculate the flux of the IC spectra at δt ≈ 12.6 days and find that the IC flux at ν X is ≈ 3.8 × 10 −10 mJy, which is a factor ≈ 10 4 times lower than the X-ray flux at that time. Therefore, we also conclude that IC effects are not contributing to the excess flux.
Instead, the discrepancy between the X-ray and 5.8 GHz model light curves and observations could be reconciled with a slightly later jet break at δt ≈ 2.0 days. However, this would violate the 22 GHz upper limit at δt ≈ 12.1 days. In summary, there is not a natural explanation for the late-time excess emission in these bands, although fixing the time of the jet break to be later does not significantly affect the parameters of interest (E K , n 0 , and A V,GRB ), and they remain consistent with the fit presented above within errors.
Both Zauderer et al. (2013c) and Kangas & Fruchter (2021) found the wind environment to best fit the afterglow of GRB 110709B, whereas we found the ISM environment to best fit the afterglow of GRB 110709B. The discussion of our wind environment model fit and comparison to the fits of Zauderer et al. (2013c) and Kangas & Fruchter (2021) can be found in Appendix B.1.

GRB 111215A
We compiled the radio and millimeter light curves of the afterglow of GRB 111215A, along with NIR/optical upper limits and the Swift X-ray light curve for our afterglow modeling (Section A.6.2). The radio afterglow of GRB 111215A is one of the best sampled of any dark GRB, with 1.4 GHz non-detections spanning δt ≈ 10 − 237 days, 4.8, 4.9, and 6.7 GHz observations spanning δt ≈ 1.4 − 238 days (all together C-band), 8.5 GHz observations spanning δt ≈ 15 − 88 days (X-band), 15 GHz observations spanning δt ≈ 34 − 79 days, 19.1 and 24.5 GHz observations spanning δt ≈ 3 − 88 days (all together K-band), 86.7, 93, 93.7, and 94.5 GHz observations spanning δt ≈ 1 − 74 days (all together 3mm), 104.7 GHz detection at δt ≈ 2.7 days, and a 230 GHz upper limit at δt ≈ 17 days. The deepest optical and NIR upperlimits of GRB 111215A span δt ≈ 0.2 − 0.3 days. The Swift X-ray afterglow of GRB 111215A spans δt ≈ 4.9 × 10 −3 − 16.7 days. GRB 111215A has previously been modeled by Zauderer et al. (2013c) and van der Horst et al. (2015), though neither of them included IC effects. Additionally, Kangas & Fruchter (2021) modeled GRB 111215A with the inclusion of IC effects, but they chose not to fit for A V,GRB . Like for GRB 110709B, we choose to model GRB 111215A in our modeling framework both for consistency across our sample, and to determine A V,GRB .

Basic Considerations
We first consider the radio and millimeter light curves of GRB 111215A to determine the location of the radio observations in relation to ν sa and ν m . At δt ≈ 8.0 days and δt ≈ 10.4 days, the spectral index between the Cband and 3mm afterglow is β ≈ 0.37 − 0.39, suggesting the spectral ordering of ν sa < C-band < 3mm < ν m (expected β = 1/3). The 3mm light curve over δt ≈ 8.5 − 12.0 days is roughly constant, with F ν,3mm ≈ 1.4 mJy, implying the burst occurred in a wind environment (expected α = 0). We find the same conclusion from the C-band light curve, where from δt ≈ 6.1 − 10.1 days the afterglow flux matches F ν ∝ t 0 . Thus the radio and millimeter afterglow of GRB 111215A is consistent with a wind environment, and indicates ν sa < Cband < 3mm < ν m at least from δt ≈ 8.0 − 12 days.
We now turn our attention to an interesting spectral phenomena in our radio and millimeter afterglow observations. As discussed, the spectrum between Cband and 3mm at δt ≈ 8.0 and 10.4 days indicates that ν sa < C-band < 3mm < ν m in a wind environment. This conclusion appears to still hold true for the spectrum at δt ≈ 45.25 days (β ≈ 0.12), as the passage of ν m through the millimeter and radio would result in a negative spec-tral slope. However, between δt ≈ 11.8 − 26.5 days, an additional bump in the spectrum peaking between X-band and K-band is inconsistent with the expected β = 1/3 ( Figure 3). A possible explanation for this spectral feature is a reverse shock. However, the exploration of this possibility is outside the scope of this work.
We now consider the X-ray spectrum and light curve to determine the location of ν c and p. The X-ray afterglow of GRB 111215A exhibits flaring until δt ≈ 0.02 days ( Figure 3). Thus, we only consider the X-ray afterglow for δt ≈ 0.02 − 16.7 days, where we generate a PC-mode time-sliced spectra from the Swift online tool, which found Γ X = 2.11 ± 0.06 (1σ, Evans et al. 2009). The X-ray afterglow of GRB 111215A is characterized by a spectral index β X = −1.11 ± 0.06 and can be fit with a single temporal power law characterized by α X = −1.43 ± 0.02 (Figure 3).
In a wind environment, the scenario of ν m < ν X < ν c β X implies p = 3.22 +0.12 −0.11 , which is inconsistent with the α X derived value of p = 2.25±0.03, and we therefore rule out this spectral frequency ordering. However, in the case of ν m , ν c < ν X , we require p = 2.22 ± 0.12 to match the measured value of β X , which is consistent within 3σ to the α X derived value of p = 2.58±0.03. Therefore, we conclude that ν m , ν c < ν X and p ≈ 2.2 − 2.6. Moreover, the X-ray light curve does not exhibit any break to δt ≈ 17 days, placing a lower limit on the time of the jet break to t jet 17 days ( Figure 3).
In conclusion, the radio afterglow of GRB 111215A is consistent with a wind environment where ν sa < Cband < 3mm < ν m . Additionally, the X-ray afterglow is consistent with ν m , ν c < ν X . We find a preliminary estimate of p ≈ 2.2 − 2.6, and expect t jet 17 days.

MCMC Modeling
We fit the afterglow data of GRB 111215A with a wind environment. We present the best fit parameters wind model in Figure 3 and list the parameters as well as the summary statistics from the marginalized posterior density functions (medians and 68% credible intervals) in Table 2.
For our best fit wind model, the optical/NIR limits imply A V,GRB 10.6 mag. Our model parameters are p ≈ 2.88 and t jet ≈ 26.0 days, resulting in a θ jet ≈ 2.6 • . The value of p is higher than our initial prediction of p ≈ 2.2 − 2.6 (see Section 3.2.1). We can reconcile the higher value of p found in the full modeling compared to the simplistic calculations performed earlier by investigating the IC cooling effects for our best fit parameters. We find that the Compton Y -parameter decreases from Y ≈ 6.4 at the time of the fast-to-slow cooling transition  Figure 3. Left: X-ray, optical, millimeter and radio afterglow light curves of GRB 111215A, together with the best-fit forward shock model in a wind environment (lines). Circles represent literature data van der Horst et al. 2015), and triangles represent 3σ upper limits. Open symbols indicate data that are not included in the fit, and shaded regions represent variability due to scintillation. The inset shows the model z -, J-, and Ks-band light curves (solid lines) as well as the non-extinguished models (dashed lines), indicating AV,GRB 10.6 mag to explain the upper limits. Right: Radio and millimeter spectral energy distributions (SEDs) of the afterglow of GRB 111215A from δt ≈ 1.5-88.23 days, together with the best-fit forward shock model (lines). Overall the model provides a good match to the broad-band temporal and spectral evolution of the afterglow.
(δt ≈ 0.2 days) to Y ≈ 0.7 at the time of the last X-ray detection (δt ≈ 14.3 days). This decrease in Y results in a faster evolution of ν c , and results in a shallower light curve at higher p. Consistent with the likely existence of an additional component (Section 3.2.1), our FS model under-predicts the light curves of the X-and K-band afterglow at δt 45 days, as well as the 3mm light curve at δt 8.0 days.
Similar to the work presented here, Kangas & Fruchter (2021) modeled GRB 111215A in a wind environment with IC effects included, finding similar parameter values to us, though our model most disagrees with Kangas & Fruchter (2021) on the placement of t jet , where they found t jet ≈ 10.8 days, driven mainly by the fit to the early data.

GRB 130131A
We compiled the radio light curves (6.0 GHz, 19.2 GHz, and 24.5 GHz, spanning δt ≈ 0.7 − 13.9 days) and millimeter light curve (85.5 GHz, spanning δt ≈ 0.9 − 2.8 days) of the afterglow of GRB 130131A (Section A.1.2), along with the optical/NIR detections (RJHK, spanning δt ≈ 0.04 − 0.67 days) and the Swift X-ray light curve spanning δt ≈ 7.5 × 10 −4 − 4.7 days (Section A.1.1) for our afterglow modeling. . Left: X-ray, optical, NIR, millimeter and radio afterglow light curves of GRB 130131A, together with the best-fit forward shock model in a wind environment (lines). Squares represent data newly reported here, and triangles represent 3σ upper limits. Open symbols indicate data that are not included in the fit, and shaded regions represent variability due to scintillation. The inset shows the model J-and R-band light curves (solid lines) as well as the non-extinguished models (dashed lines), indicating AV,GRB ≈ 2.2 mag to explain the detections. Right: Radio and millimeter spectral energy distributions (SEDs) of the afterglow of GRB 130131A from δt ≈ 0.7-13.9 days, together with the best-fit forward shock model (lines). Overall the model provides a good match to the broad-band temporal and spectral evolution of the afterglow.

Basic Considerations
To determine p and the location of ν c , we turn our attention to the X-ray afterglow. We ignore the X-ray flaring activity of GRB 130131A (δt ≈ 0.08 − 2.28 × 10 −2 days, Figure 4), and only consider the X-ray light curve and spectra from δt ≈ 0.02 − 4.7 days. We generate a PC-mode time-sliced spectra from the Swift online tool, which found Γ X = 2.4 +0.3 −0.2 (1σ, Evans et al. 2009). The X-ray spectra for this time range can be characterized with β X = −1.4 +0.3 −0.2 and the X-ray light curve for this time range can be fit with a single power law with α X = −1.14 ± 0.04. We first consider the scenario in which ν m < ν X < ν c . In this regime, β X would indi-cate p = 3.8 ± 0.5. However, the temporal decline of the X-ray light curve would indicate p = 1.85 ± 0.06 in the wind environment and p = 2.52 ± 0.06 in an ISM environment. While the values of p found for the ISM environment are consistent within 3σ, the nominal values of p for the spectral and light curve analysis differ by ≈ 1.3. On the other hand, for the regime ν m , ν c < ν X , α X = −1.14 ± 0.04 implies p = 2.19 ± 0.06, and β X = −1.4 +0.3 −0.2 implies p = 2.8 ± 0.5. These values are more consistent (2σ), with a nominal value difference of only ≈ 0.6. Therefore, the X-ray afterglow of GRB 130131A indicates ν m , ν c < ν X and p ≈ 2.2 − 2.8. Additionally, the X-ray light curve does not exhibit a break to δt ≈ 2.2 days, placing a lower limit on any jet break to t jet 2.2 days. In this regime we cannot discriminate between the ISM and wind environments based on the X-ray afterglow alone.
We now consider whether the radio afterglow can provide additional constraints on the properties of GRB 130131A. The radio light curves of GRB 130131A are not well sampled, with only 4 epochs of observations (and 7 detections) between the 3 frequencies. Additionally, the 6.0 GHz light curve exhibits variability, with detections at δt ≈ 0.7 and ≈ 4.7 days, interspersed with a non-detection at δt ≈ 2.0 days. Variability on short time scales at ν obs 10 GHz is likely attributable to scintillation (Rickett 1990), and we therefore ignore the 6.0 GHz light curve in our analytical arguments (but include it, along with the anticipated scintillation effects, in our MCMC modeling). Turning our attention to the 19.2 GHz and 24.5 GHz observations, we note that the observations have a positive spectral slope of β ≈ 0.6 at δt ≈ 0.7 days, and both the 19.2 GHz and 24.5 GHz light curves fade significantly at δt > 2.0 days, requiring α −0.7 to be consistent with the later nondetections. These spectral and temporal indices can be consistent with either a wind environment pre-jet break if ν sa < 19.2 − 24.5 GHz < ν c < ν m (expected β = 1/3 and α = −2/3), or an ISM or wind environment if t jet ≈ 2.0 days and ν sa < 19.2 − 24.5 GHz < ν m (expected β = 1/3 and α = −1/3). Therefore, we are unable discriminate between the ISM and wind environment with the 19.2 GHz and 24.5 GHz afterglow observations.

MCMC Modeling
As we cannot distinguish between the wind and ISM environments from our afterglow observations of GRB 130131A, we fit the data with both an ISM and wind environment. While the host galaxy of GRB 130131A does not have a spectroscopically determined redshift (see Section A.1.3), we assume a photometric redshift of z = 1.55 for both fits, based on host galaxy SED fitting (see Section 4). We find our wind environment model is marginally preferred by the data, providing a χ 2 /d.o.f. ≈ 25/187 and L ≈ 137, and we use this model for broader analysis in the rest of the paper. Comparatively, our ISM environment best fit model produced a χ 2 /d.o.f. ≈ 32/187 and L ≈ 132 (we provide the ISM fit in Appendix B.2 for completeness). We present the best fit wind model in Figure 4 and list the parameters as well as the summary statistics from the marginalized posterior density functions (medians and 68% credible intervals) in Table 2.
In confirmation of the arguments laid out in Section 3.3.1, we find p ≈ 2.3, and t jet ≈ 3.5 days from our best fit wind model. The spectrum remains in the fast cooling phase until δt ≈ 2.0 days, and ν m , ν c < ν X for the entirety of the X-ray afterglow that we consider. We find the optical/NIR afterglow detections are best fit with a line-of-sight extinction value of A V,GRB ≈ 2.2 mag.
The model matches the first epoch of 19.2 GHz and 24.5 GHz observations, but under-predicts the light curves at δt 2.0 days by 2 − 4σ. This discrepancy can be reconciled somewhat by setting t jet 14 days (past the time of all of our afterglow observations), but the model parameters of interest (E K , n 0 , and A V,GRB ) are not significantly affected by this change. As suggested in Section 3.3.1, we expect strong scintillation at 6.0 GHz based on our model, consistent with the large variability observed in the 6.0 GHz light curve.
We first investigate the X-ray light curve of GRB 131229A to determine p, the location of ν c , and place limits on t jet . The combined windowed timing (WT) mode and PC-mode X-ray light curve can be characterized by a broken power law (Equation 1), with α 1 ≈ −1.0 and α 2 ≈ −1.4. As the break in the light curve occurs between the WT-mode and PC-mode (δt ≈ 0.39 − 4.8 × 10 −2 days), we investigate each mode separately. The WT-mode X-ray light curve of GRB 131229A is characterized by a single power law of α WT = −1.03 ± 0.04, where as the PC-mode X-ray light curve of GRB 131229A is characterized by a single power law of α PC = −1.36 ± 0.03. The change in temporal index of ∆α (see Section 3.1.1) is ≈ 0.34, indicating this temporal break may be the passage of ν c through the X-ray band (expected ∆α = 0.25).
If the break in the X-ray light curve is indeed the passage of ν c , the derived p from the temporal and spectral indices should be consistent before and after the break. We create a WT-mode time-sliced spectra from the Swift online tool, and find the WT-mode (δt ≈ (1.2 − 3.9) × 10 −3 days) photon index to be Γ WT = 1.84 +0.05 −0.05 (1σ, Evans et al. 2009), corresponding to a WT-mode spectral index of β WT = −0.84 +0.05 −0.05 . The hardness of this spectrum implies ν m < ν X < ν c , and we derive p = 2.68 +0.09 −0.09 , assuming β WT = (1 − p)/2. In a wind environment, the measured value of α WT yields p = 1.71 +0.05 −0.05 , which is inconsistent with our β WT de- Right: Optical to X-ray SED of GRB 131229A at 0.0368 days, with optical/NIR upperlimits (purple) and estimated Fν c,min at νc,max (green). Dashed grey line represents the extrapolation from Fν c,max assuming νm < NIR/optical < νc with no extinction (AV,GRB = 0 mag). Black solid line represents the extinction curve with the minimum extinction required to be consistent with the upper limits (AV,GRB 9.6 mag).
rived p, and we therefore rule out a wind environment.
On the other hand, in an ISM environment, the measured value of α WT yields p = 2.37 +0.05 −0.05 , which is consistent with the β WT derived p within 3σ. We find a mean weighted WT-mode p of p WT = 2.44 ± 0.04. Past the break (δt ≈ 0.05−1.32 days), we create a PC-mode timesliced spectra from the Swift online tool, which finds the PC-mode photon index to be Γ PC = 2.14 +0.11 −0.10 (1σ, Evans et al. 2009), corresponding to a PC-mode spectral index of β PC = −1.14 +0.11 −0.10 . Assuming ν m , ν c < ν X , we find β PC yields p = 2.28 +0.21 −0.20 , consistent with p WT . Furthermore, α PC yields p = 2.49 +0.04 −0.04 , once again consistent with p WT , and we find a mean weighted PCmode p of p PC = 2.48 ± 0.04. We therefore conclude that an ISM environment is preferred for GRB 131229A, and that ν c passes through the X-ray band between δt ≈ (0.39 − 4.8) × 10 −2 days, yielding an overall mean weighted p of p = 2.46 ± 0.03. Furthermore, the X-ray light curve does not exhibit any additional steepening, placing a limit of t jet 1.3 days.
We next examine the radio to X-ray SED at the time of the radio observations (δt ≈ 0.91 days) to place constraints on the location of ν m . The shallow radio to X-ray spectral index at δt ≈ 0.91 days (β radio−X ≈ −0.39, dotted line, Figure 5) is inconsistent with a single, optically thin power law spectrum with index β = (1 − p)/2 ≈ −0.73 extending from the radio to the X-rays, and instead requires a spectral peak in between.
In the regime of 6.0 GHz < 21.8 GHz < ν m < ν c , we expect a positive spectral index of β = 1/3 for the radio band. However, the observed limit of β −0.38 between the 6.0 GHz and 21.8 GHz observations is inconsistent with this expectation. The discrepancy between the tentative detection at 6.0 GHz and the upper limit at 21.8 GHz may be explained by a variety of factors, such as the 6.0 GHz counterpart being unrelated to the FS afterglow of GRB 131229A, or due to scintillation effects which can cause variability on short timescales at ν obs 10 GHz (Rickett 1990). We instead utilize the 21.8 GHz non-detection to place constraints on the FS emission.
We use the 21.8 GHz non-detection and the X-ray light curve to place constraints on the location of ν c and ν m . Based on our analytical arguments, the latest time ν c can reasonably pass through the X-ray band is at the start of the PC-mode X-ray light curve (δt ≈ 4.8×10 −2 days), and therefore we assume ν c ≈ ν X at this time. Scaling ν c to the time of the radio observations (ν c ∝ t −1/2 ), we find the maximum value of ν c to be ν c,max ≈ 5.5 × 10 16 Hz at δt ≈ 0.91 days (resulting in F ν,c,min ≈ 0.65 µJy). We fit a broken power law (Equation 2) with the 21.8 GHz non-detection and F νc,min , fixing β 1 = 1/3, β 2 = (1 − p)/2 ≈ −0.73, and s = 1.84 − 0.40p ≈ 0.856. We find ν m 10 13 Hz (solid line, Figure 5).
Our X-ray analysis provides measured values of the synchrotron flux above ν c , and our radio analysis provides a constraint on the synchrotron flux between ν sa and ν m ( Table 1 in Granot & Sari 2002). Combined with our constraint on ν c ( Table 2 in Granot & Sari 2002), we place limits on the energy, density, and microphysics of the system. Assuming a redshift of z = 1.04 (see Section 4), we find E K,iso 3.21 × 10 1 1/3 B 10 52 erg, n 0 1.72 × 10 −5 −5/3 B cm −3 , and e 1.19 × 10 −2 −1/3 B . Using these constraints, along with our constraint on t jet , we place constraints on θ jet and find θ jet 1.32 (Sari et al. 1999).
We next place constraints on A V,GRB of GRB 131229A. Assuming ν m < NIR/optical < ν c , and using our assumption of ν c ≈ ν X at δt ≈ 4.8×10 −2 days, we interpolate the SED between ν m and ν c at the time of the most constraining optical/NIR upper-limits (δt ≈ 0.0368 days). We find that the extinction necessary to be consistent with the optical/NIR upper limits is A V,GRB 9.6 mag (see Figure 5).
In conclusion, we find for GRB 131229A p = 2.46 ± 0.03, t jet 1.3 days, and A V,GRB 9.6 mag. We also place constraints on the isotropic energy and circumburst density of GRB 131229A of E K,iso /10 52 erg 3.21 × 10 1 1/3 B and n 0 1.72 × 10 −5 −5/3 B , leading to further constraints on θ jet . With limited radio observations, it is difficult to derive meaningful constraints on the afterglow properties of GRB 131229A within our MCMC modeling framework, and therefore we instead report our analytically derived values in Table 2.

GRB 140713A
We compiled the radio and millimeter light curves of the afterglow of GRB 140713A (Sec A.5.2), with WSRT and GMRT upper limits at 1.4 GHz spanning δt ≈ 11.1 − 25.1 days, WSRT observations at 4.8 GHz spanning as well as VLA observations at 4.9 and 7.0 GHz δt ≈ 3.3 − 81.0 days (all together, C-band), VLA observations at 13.4 and 15.9 GHz, as well as AMI observations at 15.7 GHz spanning δt ≈ 0.1 − 81.0 days (all together, K u -band), and CARMA observations at 85.5 GHz and PdBI observations at 86.7 GHz spanning δt ≈ 1.5−63.1 days (all together, 3mm). In combination with the radio and millimeter light curves, we include the NOT optical limits at δt ≈ 0.2 days and the Swift X-ray light curve spanning δt ≈ 1.5 × 10 −3 − 1.9 days in our afterglow modelling (Sec A.5.1).
The afterglow of GRB 140713A has been previously modeled by Higgins et al. (2019), though they did not include IC effects. We choose to model GRB 140713A in our modeling framework both for consistency across the bursts in our sample, and because we are introducing new VLA and 3mm afterglow observations that have not previously been modeled.

Basic Considerations
We observe a steep decline in the 3mm light curve, with α ≈ −2.2 at δt 12 days. We cannot reconcile this steep decline with the standard synchrotron model of a spherical blast-wave, as the light curve slopes predicted by the standard model are too shallow. We therefore conclude that this decline is caused by a jet break, and that t jet 12 days. Furthermore, this decline indicates that p ≈ 2.2, ν m ≈ 3mm, and F ν,m ≈ 1.7 mJy at δt ≈ 12 days (see Section 3).
Assuming ν m ≈ 3mm at δt ≈ 12 days, and evolving this break frequency and F ν,m forward in time (ν m ∝ t −2 , F ν,m ∝ t −1 post jet break, Sari et al. 1999) from 12 days, we find that ν m ≈ K u -band at δt ≈ 28 days, with a characteristic flux of F ν,m ≈ 0.76 mJy. Indeed, we can fit the 15.7 GHz light curve with a steep power law of α ≈ −2.0 at δt 27 days, and find the flux at 27 days to be ≈ 0.78 mJy. Therefore, the K u -band data corroborates that t jet 12 days, and we conclude that p ≈ 2.0 − 2.2. With t jet 12 days, the majority of our radio and millimeter afterglow observations are taken at δt t jet , and we are unable to discriminate between an ISM and wind environment using these observations.
We now determine the location of the X-rays in relation to ν c . We ignore the X-ray flare (δt ≈ (0.15 − 1.9) × 10 −2 days, see Figure 6) of GRB 140713A, and as such we only consider the X-ray light curve and spectra from δt ≈ 0.06 − 1.88 days in the synchrotron framework to determine the location of the X-rays in relation to ν c . We create a time-sliced PC-mode spectra from the Swift online tool, which finds the X-ray photon index to be Γ X = 2.0±0.2 (1σ, Evans et al. 2009). The X-ray spectra for this time range is characterized by β X = −1.0 ± 0.2 and the X-ray light curve for this time range can be fit with a single power law with α X = −0.9 ± 0.1.
For p ≈ 2.1, derived from the identification of a jet break in the radio afterglow, we would expect a spectral index of β The measured value of β X = −1.0±0.2 is more consistent with the latter case, within 1σ, where as the former case is only consistent within 3σ. In the regime of ν m , ν c < ν X , we would expect α X = (2 − 3p)/4 ≈ −1.1, consistent with our measured value of α X ≈ −0.9 within 2σ. Therefore, we conclude ν m , ν c < ν X , and note that in this regime we can not discriminate between the ISM and wind environment with the X-ray light curve.
In conclusion, the radio and X-ray afterglow of GRB 140713A is consistent with p ≈ 2.0 − 2.2, t jet 12 days, and ν m , ν c < ν X . Neither the radio nor X-ray observations allow us to analytically distinguish between the ISM and wind environment. . Left: X-ray, optical, millimeter and radio afterglow light curves of GRB 140713A, together with the best-fit forward shock model in a wind environment (lines). Squares represent data newly reported here, circles represent literature data (Higgins et al. (2019)), and triangles represent 3σ upper limits. Open symbols indicate data that are not included in the fit, and shaded regions represent variability due to scintillation. The inset shows the model r-, i-, and z-band light curves (solid lines) as well as the non-extinguished models (dashed lines), indicating AV,GRB 3.5 mag to explain the upper limits. Right: Radio and millimeter spectral energy distributions (SEDs) of the afterglow of GRB 140713A from δt ≈ 3.0-31.5 days, together with the best-fit forward shock model (lines). Overall the model provides a good match to the broad-band temporal and spectral evolution of the afterglow.

MCMC Modeling
As we cannot distinguish between the wind and ISM environments from our afterglow observations of GRB 140713A, we fit the data with both an ISM and wind environment. We find that our wind environment model better fits the data, with the best fit model having a reduced χ 2 ≈ 114/102 and L ≈ 79. Comparatively, our ISM environment best fit model produced a reduced χ 2 ≈ 272/102 and L ≈ 11 (for completeness, we present the best fit ISM model in in Appendix B.3). We present the best fit wind model in Figure 6 and list the parameters as well as the summary statistics from the marginalized posterior density functions (medians and 68% credible intervals) in Table 2.
In confirmation of the arguments laid out in Section 3.5.1, the parameters of our best-fit model are p ≈ 2.17 and t jet ≈ 3.6 days, resulting in a θ jet ≈ 21 • . Additionally, the SED remains in the fast cooling phase until δt ≈ 35 days, and the ordering of the break frequencies at δt ≈ 1.65 days is ν c < ν m < ν X . For our best fit wind model, the NOT optical limits imply A V,GRB 3.52 mag.
The X-ray model light curve under-predicts the data at δt 0.7 days. This is not unexpected, as we found in Section 3.5.1 that the X-ray temporal slope was shal- Note-Host galaxy properties derived by Prospector.
lower than the expected slope for p ≈ 2.1. We first consider whether this excess flux in the X-ray light curve is due to KN effects, which become important when ν X ≥ ν m . Given our best fit parameters, we find that at δt ≈ 1.65 days, ν m ≈ 1.7 × 10 21 Hz ν X . Therefore, we conclude that KN effects are not causing the excess flux in the X-ray light curve. We next consider whether IC effects are the cause of the excess flux in the X-ray light curve. We calculate the flux of the IC spectra at ν X at δt ≈ 1.65 days, and find that the IC flux at ν X is ≈ 4.7 × 10 −7 mJy, ≈ 80 times smaller than the X-ray flux at that time. Therefore, IC effects cannot account for the excess X-ray flux, although we note that such excess emission has been seen in other events Margutti et al. 2015;Laskar et al. 2018aLaskar et al. , 2019. Higgins et al. (2019) have previously modeled the afterglow of GRB 140713A. Their model allowed for p < 2, and they found a value of p ≈ 1.85, smaller, but not far off from our value of p ≈ 2.17. Their other parameters are also similar to ours, with the biggest difference in our model being that we identify a jet break in the 3mm light curve, and therefore find t jet ≈ 3.61 days, where as they predict t jet ≈ 25 − 30 days.

HOST GALAXY MODELING
To model the stellar population properties of the host galaxies, we use the stellar population inference code Prospector (Leja et al. 2017). Prospector determines properties such as total mass formed, age of the galaxy at the time of observation (t age ), optical depth of old and young stars, stellar metallicity (Z * ), the star formation history, and redshift using the available photometric and/or spectroscopic data for each host. We apply a nested sampling fitting routine with dynesty (Speagle 2020) to the observational data of each host to produce posterior distributions in each property. Model SEDs are built using Python-fsps (Flexible Stellar population synthesis; Conroy et al. 2009;Conroy & Gunn 2010). Unless the redshift of a host is known, we allow redshift to be a sampled parameter. For hosts with spectra, we fit their spectral continuum with a 10 th order Chebyshev polynomial and add a gas-phase metallicity and gas ionization parameter to accurately fit the nebular emission lines. We also assume a Chabrier initial mass function (IMF) (Chabrier 2003), Milky-Way Dust Extinction Law (Cardelli et al. 1989), and a parametric delayed-τ SFH (SFH ∝ t * e −t/τ ), where the τ is a sampled parameter in the Prospector fitting. Furthermore, we apply the Gallazzi 2005 Mass-Metallicity relation (Gallazzi et al. 2005) and a 2:1 ratio in the dust attenuation between old and young stars respectively, as stellar populations are noticed to follow this trend (Calzetti et al. 2000;Price et al. 2014). The total dust attenuation in optical depth is converted to a V -band magnitude, and hence forth referred to as A V,Host . We follow the methods in Nugent et al. (2020) to determine the mass-weighted age t m , stellar mass (M * ), and star formation rate (SFR).
We model the host galaxy of GRB 110709B with photometry from Zauderer et al. (2013c) and Selsing et al. (2019) with a fixed z = 2.109 (Perley et al. 2016a;Selsing et al. 2019). We find that the stellar population has dust extinction A V,Host = 0.61 +0.26 −0.25 mag (Table 3). Though the Prospector SED model fits the photometric data well, the model is based on only three detections and one limit (Figure 7).
We fit the host galaxy of GRB 111215A with photometry from van der Horst et al. (2015), corrected for Galactic extinction in the direction of the burst at a fixed the redshift of z = 2.012 (van der Horst et al. 2015;Chrimes et al. 2019). We find the host has dust attenuation A V,Host = 0.91 +0.06 −0.08 mag. Our host galaxy properties are similar to those found by van der Horst et al. (2015).
For the host of GRB 130131A, the redshift is unknown, although the spectrum indicates 1.3 < z < 4 (see Section A.1.3). Thus, we leave redshift as a free parameter with a flat prior of 1.3 < z < 4. We find a photometric redshift of z = 1.55 +0.01 −0.05 , and A V,Host = 1.14 +0.06 −0.04 mag. For the host of GRB 131229A, the redshift is unknown, although the deep optical limits ( 23 − 26 mag) indicates a redshift of z 1 − 1.5 (see Section A.4.3). Thus, we leave redshift as a free parameter with a flat prior of 0.1 < z < 3. We find a photometric redshift of z = 1.04 +0.28 −0.36 , and A V,Host = 1.17 +0.72 −0.66 mag. We jointly fit the photometric and spectroscopic data (Table A3) of the host of GRB 140713A. We find that the stellar population has a low A V,Host = 0.21 +0.04 −0.04 mag. We also find that the photometry and spectrum of the host are overall well-fit by the Prospector SED model, especially the [OII] (λ3727Å) and Hγ spectral line strengths (Figure 8). Our host galaxy properties are similar to those found by Higgins et al. (2019).
Finally, we fit the host galaxy of GRB 160509A, with a redshift of z = 1.17 (Laskar et al. 2016;Kangas et al. 2020). When we fit the full host galaxy SED of Keck/LRIS, HST/WFC3, and Spitzer photometry, our Prospector model over-predicts the Keck/LRIS photometry by an order of magnitude. This may indicate that the Keck/LRIS observations are dominated by afterglow contribution, in contradiction with Laskar et al. (2016). However, this Prospector model also finds a high A V,Host of ≈ 4.7 mag, much higher than what is expected for the normal dark GRB host population (i.e. . This high value of A V,Host is driven by the color between the HST/WFC3 and Spitzer bands, and thus we also fit the host of GRB 160509A with only the Keck/LRIS and Spitzer photometry. This method 10000 5000 30000 50000 Observed Wavelength (Å) results in a more expected A V,Host ≈ 1.6 mag, but the model under-predicts the HST/WFC3 photometry by a factor of ≈ 5. Without further observations of the host galaxy of GRB 160509A, including re-observing at later times in g -and r -band, we cannot conclusively determine which fit is correct, and as such we quote the ranges for both fits in Table 3 and present both SEDs in Figure 7. The full details of the best-fit properties for all of the host galaxies of our sample can be found in Table 3 and the model fits are displayed in Figure 7 and Figure 8.

DISCUSSION
We have presented multi-wavelength observations and modeling of the afterglows of five dark GRBs (GRBs 110709B, 111215A, 130131A, 131229A, 140713A), as well as SED modeling of the host galaxies of six dark GRBs with Prospector (GRBs 110709B, 111215A, 130131A, 131229A, 140713A, 160509A). We have classified two additional long GRBs as dark (GRBs 130420B, 130609A) and presented their radio and millimeter observations. However, for these two events, there is insufficient afterglow and host galaxy follow-up to allow for more in-depth modeling. For the purposes of this discussion, we group bursts that have been classified as dark (either explicitly, or satisfying β OX < 0.5, i.e. Cenko et al. 2009;Melandri et al. 2012;Krühler et al. 2012;Rossi et al. 2012;Littlejohns et al. 2015) and/or "dusty" (i.e. Krühler et al. 2015;Perley et al. 2016a) and refer to both groups as "dark". Equipped with this larger sample, including the new bursts presented here with particularly high extinction site-lines, we now investigate whether dark GRBs differ from the broader long GRB population in terms of their γ-ray, afterglow, and host properties. In our comparisons we will refer to any long GRB not classified as dark as a "typical" GRB, for brevity.

The γ-ray and afterglow properties of Dark GRBs
To determine whether dark GRBs differ from the typical GRB population, we first examine their γ-ray properties: fluence (f γ , 15-150 keV band) and duration (T 90 ), as even the dark GRBs with sparse afterglow data often have uniformly derived γ-ray properties from Swift/BAT. The exception to this is GRB 160509A, which was instead discovered by the Fermi Large Area Telescope (LAT; Longo et al. 2016). We plot these properties from the catalog in Lien et al. (2016) in Figure 9 for 103 dark GRBs (as classified by Jakobsson et al. 2004;Castro-Tirado et al. 2007;Cenko et al. 2009;van der Horst et al. 2009;Krühler et al. 2011Krühler et al. , 2012Zauderer et al. 2013c;Hunt et al. 2014;Chrimes et al. 2019, and This Work) as red points, and all other typical GRBs as blue points. We find that there is a notable lack of dark GRBs in the parameter space corresponding to low fluence (f γ 2 × 10 −7 erg/cm 2 ) and short duration (2 s < T 90 5 s). For each parameter, we test the null hypothesis that the the dark GRB population is drawn from the same distribution as the typical GRB population using a two-sample Kolmogorov-Smirnov (KS) test from the scipy.stats package, where a value of p < 0.05 rejects the null hypothesis. We obtain p = 0.0001 for the f γ distribution and p = 0.0503 for the T 90 distribution, implying that  Lien et al. (2016). Dark GRB classification is based on the either explicit classification, or satisfying βOX < 0.5 (Cenko et al. 2009;Krühler et al. 2011Krühler et al. , 2012Melandri et al. 2012;Rossi et al. 2012;Hunt et al. 2014;Schady et al. 2014;Jeong et al. 2014;Littlejohns et al. 2015;Perley et al. 2016a;Chrimes et al. 2019, This Work). The dashed horizontal line represents the theoretical minimum fγ,min ≈ 2.5 × 10 −7 erg/cm 2 , below which the X-ray (and therefore optical) afterglow is typically too faint to be classified as dark with ground based rapid optical follow-up that reaches limits of R > 24 mag. A simple 2 sample Kolmogorov-Smirnov (KS) test from the scipy.stats package which tests the null hypothesis that dark GRBs are drawn from the same population as the typical GRB population, finds p = 0.0001 for the fγ distribution and p = 0.0503 for the T90 distribution, where a value of p < 0.05 rejects the null hypothesis.
dark GRBs are not drawn from the typical GRB population in terms of their fluence, but may be in terms of their duration.
One possible explanation for the lack of dark GRBs in the low f γ and short T 90 parameter space may be an observational bias in missing dark GRBs with low fluence 6 , as these bursts have been shown to have systematically fainter X-ray afterglows, and in turn fainter optical afterglows than the rest of the population (Gehrels et al. 2008). Thus, such bursts may have afterglows which are fainter than the sensitivity threshold of optical afterglow searches, preventing accurate classification of these lower f γ events as dark. To quantify this effect, we use the Jakobsson et al. (2004) darkness classification of β OX < 0.5, and assume prompt optical observations (δt ≈ 0.1 days) of R > 24 mag, which represents normal GRB follow-up capabilities. The minimum X-ray flux at 0.1 days needed to accurately classify a burst as dark is F X,min ≈ 4.1 × 10 −2 µJy. We extrapolate F X,min to δt = 11 hr assuming F X ∝ t −1 (see Nousek et al. 2006;Zhang et al. 2006;Evans et al. 2009), and use the derived F X,11 hr −f γ relation of Gehrels et al. (2008) to calculate the minimum fluence of f γ,min ≈ 2.5 × 10 −7 erg cm −2 necessary to produce an X-ray afterglow bright enough to accurately classify a GRB as dark. This limit lies just below the majority of dark GRBs with the lowest f γ (Figure 9). Therefore, it is plausible that the lack of low f γ dark bursts is due to an observational bias, as opposed to an intrinsic effect. Additionally, as f γ and T 90 are correlated (Balázs et al. 2004), this also provides a natural explanation for the lack of observed dark bursts at 2 s < T 90 < 5 s. If we exclude bursts with f γ < 2.5 × 10 −7 erg cm −2 then the dark GRB population does become more statistically similar to the long GRB population (p = 0.0002 for the f γ distribution, p = 0.0734 for the T 90 distribution). Finally, we note that there is not a complete catalogue of all dark GRBs, and we may be missing a significant fraction of the dark GRB population. This is made apparent when one considers that the estimated fraction of dark GRBs with respect to all long GRBs is 10-50% (Jakobsson et al. 2004;Cenko et al. 2009;Fynbo et al. 2009;Greiner et al. 2011;Melandri et al. 2012;, whereas only ∼ 10% of the long GRBs in Figure 9 have been classified as dark. We next explore the dark GRB inferred burst explosion properties (e.g., kinetic energies and opening angles) to investigate whether dark GRBs differ from the typical population. As most dark GRBs do not have extensive broadband afterglow modeling, we focus on the six dark GRBs (GRBs 110709B, 111215A, 130131A, 131229A, 140713A, 160509A) in our sample which are uniformly modeled. We find that the dark GRB sample spans a wide range of beaming-corrected kinetic energies (E K ≈ 0.01 − 6.6 × 10 51 erg). Compared to the values for other long GRBs with afterglow modeling (Panaitescu & Kumar 2002;Price et al. 2002;Yost et al. 2003;Frail et al. 2005;Chandra et al. 2008;Cenko et al. 2010Cenko et al. , 2011Laskar et al. 2013cLaskar et al. , 2014Laskar et al. , 2015Alexander et al. 2017;Tanvir et al. 2018;Laskar et al. 2018aLaskar et al. ,b,c, 2019, we find that the kinetic energies of the dark GRB sample are consistent with those of the typical GRB sample, E K ≈ 0.01 − 36.0 × 10 51 erg. Additionally, our dark GRB sample spans a wide range of jet opening angles (θ jet ≈ 1.6 − 21 • ), with a distribution again consistent with that of the typical GRB population (θ jet ≈ 1.1 − 50 • ). There is no clear evidence that dark GRBs are distinct from the typical GRB population in afterglow properties. Our conclusions do not change if we broaden our sample to include the small sample of long GRBs with A V,GRB > 1 mag that have not been classified as dark (GRBs 980329 and 980703; Yost et al. 2003, and GRB 011121;Price et al. 2002).

The origin of the dust along the line-of-sight
We next consider whether our dark GRB sample has different environmental properties than typical GRBs. In practice, the higher line-of-sight extinctions derived for dark GRBs could be a result of local environment (probed by the afterglow), larger structures such as starforming regions, global host galaxy dust distributions, or a combination of all three. An investigation into the global and local environmental properties of long GRBs will help determine the cause of the extinction of dark GRBs.
First, we investigate the global host properties of dark GRBs, to examine how their hosts differ from the typical GRB host population. We gather a sample of 150 GRB hosts that have inferred stellar mass measurements  Figure 10. We find that the host galaxies of dark GRBs tend to be more massive (median log(M * /M ) ≈ 9.9) than typical GRB host galaxies (median log(M * /M ) ≈ 9.2), in alignment with previous results based on smaller samples . Moreover, we find that only ∼ 16% (12/82, excluding upper limits) of the hosts of the typical GRB population are more massive than the median host mass of the dark GRB population. A natural explanation could be that high-mass galaxies have a larger number of obscured sight lines as they overall have larger dust contents (e.g. Santini et al. 2014;Calura et al. 2017) which would result in high-A V , dark GRBs in these types of hosts .
To . Available upper limits are denoted by dashed, unfilled distribution with leftward arrow. The hosts of the dark GRB population are generally more massive than the typical GRB population (median log 10 (M * /M ) ≈ 9.9 for the hosts of dark GRB population, vs median log 10 (M * /M ) ≈ 9.2 for the typical GRB population) Bottom: Distribution of extinction, AV,Host, of 92 typical GRB hosts (blue) and 60 dark GRB hosts (red), gathered from the literature (Savaglio et al. 2009;Hunt et al. 2014;Krühler et al. 2015;Piranomonte et al. 2015;Japelj et al. 2016, This Work), where upper limits are denoted by dashed, unfilled distribution with leftward arrow. The hosts of the dark GRB population are generally more dusty than the typical GRB population (median AV,Host ≈ 1.1 mag for the hosts of dark GRB population, vs median AV,Host ≈ 0.5 mag for the typical GRB population) population in Figure 10. Indeed, we find that dark GRBs typically occur in dustier host galaxies (median A V,Host ≈ 1.1 mag) compared to typical GRBs (median A V,Host ≈ 0.5 mag) (see also ). This trend, combined with the tendency of dark GRBs to occur in more massive galaxies, implies that the mass of dark GRB host galaxies is linked to a dustier host galaxy overall. Indeed, there are known positive correlations between stellar mass and dust mass (e.g. Santini et al. 2014;Calura et al. 2017), that support this idea. The trend of dark GRBs originating in massive, dusty hosts indicates that the high line-of-sight extinction of dark GRBs is, at least in part, attributed to global host extinction.
We next investigate whether the dust extinction of dark GRBs as inferred from their afterglows is linked to the dust content of the host galaxy by determining whether A V,GRB is directly linked to the A V,Host .  suggested that the majority of long GRBs have A V,GRB values within a factor of ∼ 2 − 3 of their host A V,Host , following a near 1 : 1 relation between A V,GRB and A V,Host . This correlation between A V,GRB and A V,Host could occur, for instance, if host galaxies of GRBs have uniform dust distributions, resulting in a homogenous screen of dust in the hosts' diffuse ISM (see . To test this, we gather available A V,GRB measurements from the literature for typical (Kann et al. 2006(Kann et al. , 2010 Figure 11, and find that the entire long GRB population spans a wide range of A V,Host (∼ 0 − 4.7 mag) and A V,GRB (∼ 0 to 12.3 mag).
To test whether there is a relation between A V,Host and A V,GRB , we randomly sample A V,Host -A V,GRB pairs from the long GRB sample for which there are published A V,Host and A V,GRB values, taking into account error bars and upper/lower limits, 7 and fit a line to the randomly drawn sample using curve_fit from the scipy 7 For cases in which there are quoted error bars for an A V value (either A V,Host or A V,GRB ), we randomly sample using an asymmetric Gaussian, with the A V value as the mean, µ, and the upper and lower errors on A V as the standard deviation, σ, for either side of the Gaussian, respectively. For cases in which there is an upper limit quoted for an A V value, we randomly sample using a tophat function from 0 to the A V value. For GRB 160509A we assume a tophat function for A V,Host from 1.60-4.66 mag (See Section 4) Finally, for cases in which there is a lower limit quoted for an A V,GRB value, we randomly sample from a half Gaussian distribution with the A V,GRB lower limit as µ, and a 5σ maximum A V,GRB of 25 mag, such that σ = (25 − µ)/5. package. We repeat this process 1000 times and produce a distribution of line slopes that are fit by the A V,Host -A V,GRB pairs. We find that our distribution of A V,Host -A V,GRB slopes has a median value of ≈ 0.05, a nearly flat correlation instead of a 1 : 1 relation. Moreover, within the dark GRB sample alone, we also find only weak correlation between A V,Host -A V,GRB (median slope of ≈ 0.03). The weak correlation between A V,Host and A V,GRB could indicate that either that the high line-ofsight extinction is caused by dust extinction from the extremely local (∼parsec) environment of the dark GRB, or that the host galaxies of long GRBs have patchy, rather than a uniform, dust distributions. In this latter case, the dust extinction driving the high A V,GRB is a geometrical line-of-sight effect that is probabilistic in nature.
To investigate the contribution of the extremely local environment of the GRB (within the blast wave radius ∼ 0.2 − 50 parsec) to the high line-of-sight extinctions, we compare the circumburst densities as inferred from the afterglow. The naive expectation is that if the dust that is providing obscuration of the afterglow originates on ∼parsec scales, dark GRBs will trace environments with higher inferred densities. The sample of six dark GRBs we model in this paper spans a wide range of densities (n 0 ∼ 10 −3 − 10 1 cm −3 ; in the cases of wind environments we calculate n 0 at 10 17 cm). Compared to typical long GRBs (Panaitescu & Kumar 2002;Price et al. 2002;Yost et al. 2003;Frail et al. 2005;Chandra et al. 2008;Cenko et al. 2010Cenko et al. , 2011Laskar et al. 2013cLaskar et al. , 2014Laskar et al. , 2015Alexander et al. 2017;Tanvir et al. 2018;Laskar et al. 2018aLaskar et al. ,b,c, 2019, the densities of the dark GRBs in our sample fall well within the bounds of the typical GRB sample as a whole (n 0 ∼ 10 −5 − 10 3 cm −3 ). Additionally, we find that A V,GRB is not correlated with the circumburst density. Indeed, the long GRBs with the highest measured densities (GRBs 050904 and 120404A, n 0 ∼ 3.5 − 6.3 × 10 2 cm −3 ) have low inferred A V,GRB values ( 0.05 − 0.13 mag; Laskar et al. 2014Laskar et al. , 2015. This implies that the higher inferred line-of-site extinction for dark GRBs is not a result of the extremely local environment (∼ 0.2 − 50 pc).
In summary, we find that dark GRBs tend to occur in more massive, dustier host galaxies than typical GRBs. Additionally, the observed dust obscuration of dark GRBs favors a patchy dust distribution over a uniform one in host galaxies, as A V,GRB is only weakly correlated to A V,Host . Furthermore, the dust obscuration of dark GRBs is not purely a result of the extremely local (∼parsec) environment of the GRB, as A V,GRB is not correlated to the circumburst density. The combination of high dust content and a patchy dust distribution  . Our sample of six dark GRBs with AV,Host and AV,GRB measurements are represented by red squares with dark outlines. When necessary, we convert EB−V values to AV,Host using the convention R = EB−V /AV = 3.14 (Schultz & Wiemer 1975). To be consistent with Figure 15 from , for host galaxies with upper limits on their median host extinction, we plot the upper limit at the upper bound (+1σ) of the allowed AV,Host range. Additionally, we choose the AV,GRB derived from the SMC dust fit from Table 1 of Kann et al. (2006) and Table 3 of Kann et al. (2010), both to be consistent with Figure 15 from  and to be consistent with our sample, for which we also use an SMC dust extinction law. For clarity, for GRBs with AV,Host or AV,GRB ≈ 0.0 mag, we plot their upper bounds at 0.05 mag. The dotted line denotes a 1 : 1 relation between median host extinction and GRB line-of-sight extinction. The dashed line denotes the correlation found between AV,Host-AV,GRB pairs, with a slope of ≈ 0.05 (see Section 5.2). These lines are meant to guide the readers eye. results in a higher probability of any given line-of-sight to intersect a patch of dust, leading to a high A V,GRB and dark GRB. The combination of high line-of-sight extinctions, patchy dust distributions, and association of long GRBs with star-foming galaxies, make dark GRBs exciting probes of obscured star formation.

SFR and Radio Limits
We next place limits on the amount of obscured star formation (SF) occurring in the host galaxies of the typical GRBs and dark GRBs. Long GRBs are inherently linked to SF due to their massive star progenitors and therefore their association with star-forming galaxies (Djorgovski et al. 1998;Christensen et al. 2004;Japelj et al. 2016;Palmerio et al. 2019). Additionally, massive, dusty, galaxies such as those that host dark GRBs often have higher SFRs (Santini et al. 2014;Calura et al. 2017), which has been corroborated by dark GRB host studies (e.g. , and patchy dust distributions within dark GRB hosts may increase the probability of ongoing SF in the host to be obscured (Perley & Perley 2013). This obscured SF may result in radio bright host galaxies, as the dust becomes transparent at radio wavelengths and reveals the true SFR (SFR Radio ), whereas the SFR measured from UV/optical (e.g. emission line) diagnostics and from stellar population synthesis modeling of broad-band, galaxy-integrated photometry is only sensitive to the SFR unobscured by dust (Perley & Perley 2013). The search for radio emission from typical GRB host galaxies and dark GRB hosts is of interest as the results help answer the question of whether long GRBs are biased or unbiased tracers of SF across the universe (i.e. Perley & Perley 2013;Perley et al. 2015;Gatkine et al. 2020).
Several studies have been conducted to search for radio emission from long GRB host galaxies (Perley & Perley 2013;Stanway et al. 2014;Perley et al. 2015;Greiner et al. 2016;Perley et al. 2017;Gatkine et al. 2020;Eftekhari et al. 2021), with a handful (∼ 10) of successful detections. These studies are not necessarily representative of the overall long GRB host population, as some focused exclusively on dark GRBs, some focused on typical GRBs at high redshift (z > 2.5), and few took an unbiased approach to all long GRB hosts. As there is no large, unbiased, survey of radio observations of long GRB host galaxies, we must rely on the these smaller studies, despite their potential biases in sample selection. We collect the radio fluxes from these studies to provide a broad view of the currently available observations, as well as to make comparisons to our sample of six dark GRBs. We calculate the SFR Radio of these observations using the relation in Greiner et al. (2016) (assuming a spectral index of −0.75), and plot them against the SED-derived SFR (SFR SED , Perley et al. , 2015, supplementing with UV-derived SFR or SFR tracers such as Hα, Hβ, and OII from Hunt et al. 2014;Krühler et al. 2015), in Figure 12. For comparison, we also plot a sample of star-forming field galaxies . Radio SFR (SFR Radio ) versus SED derived SFR (SFRSED) for typical GRB host galaxies (blue) and dark GRB host galaxies (red). Triangles correspond to radio upper limits. SFR Radio derived from radio observations in the literature (Perley & Perley 2013;Stanway et al. 2014;Perley et al. 2015;Greiner et al. 2016;Perley et al. 2017;Gatkine et al. 2020;Eftekhari et al. 2021) using the SFR Radio relation in Greiner et al. (2016) and assuming a spectral index of −0.75. SFRSED is taken from the literature  and supplemented with UV derived SFR or SFR tracers such as Hα (Hunt et al. 2014;Krühler et al. 2015). Our sample of 6 dark GRBs are shown as filled in red triangles, with SFR Radio from the most constraining radio afterglow upper limits and SFRSED from our Prospector host galaxy modeling. Also shown are a sample of star forming galaxies for comparable redshifts z ≈ 0.03 − 3.4 from the VLA-COSMOS source catalog (grey, Smolčić et al. 2017). The dashed and dotted line indicate SFR Radio = SFRSED and SFR Radio = 10 × SFRSED, respectively. at comparable redshifts (z ≈ 0.03 − 3.4) from the VLA-COSMOS source catalog (Smolčić et al. 2017).
The VLA-COSMOS sources typically have SFR ratios, R SFR = SFR Radio /SFR SED , on the order of ≈ 1 − 10. We define galaxies with "significant" obscured SF as galaxies that satisfy R SFR 10. Of the ∼ 80 long GRB host galaxies that have been observed at radio wavelengths, only ∼ 10 have unambiguous host detections (Perley & Perley 2013;Stanway et al. 2014;Perley et al. 2015Perley et al. , 2017. Of the ∼ 10 detected long GRB host galaxies, only 3 display significant amounts of obscured SF, and the majority of radio detected long GRB hosts have R SFR within the range we would expect compared to other star forming galaxies (R SFR ≈ 1 − 10, i.e. the VLA-COSMOS sources).
To assess the detectability of such obscured SF in the hosts of the six dark GRBs in our sample, we now place upper limits on the SFR radio of the host galaxies. We consider the most constraining afterglow upper limits of our sample van der Horst et al. 2015;Chandra & Nayana 2014a;Laskar et al. 2016, This Work), and we place limits of R SFR 8 − 700 (also plotted in Figure 12). With the afterglow radio limits, we are unable to rule out significant amounts of obscured SF for all bursts in our sample except GRB 130131A. We also searched the VLASS (Lacy et al. 2020) for radio emission at 3 GHz at the positions of the bursts, but found the non-detections to be less constraining than the radio afterglow upper limits; this holds true even for the projected total sensitivity of VLASS of the combined 3 epochs (R SFR 30 − 2000).
With the next generation VLA (ngVLA), we would be sensitive to radio emission at the level of ∼ 0.75 µJy with one hour of observation at 2.4 GHz (Carilli et al. 2015) 8 . This radio flux would correspond to SFR Radio 0.02 − 16 M /yr for z = 0.1 − 2.0. With the ngVLA, we would be able to detect at least 5 of the 6 dark GRB hosts in our sample (with the exception of GRB 110709B, whose SFR SED is ≈ 7× lower than that of the ngVLA limit at the host redshift), assuming R SFR ≈ 1 (corresponding to unobscured SF). If we consider an unbiased sample of long GRBs, such as the Swift/BAT6 sample (Salvaterra et al. 2012), we can calculate the expected detection fraction of the long GRB host galaxies by the ngVLA. Of the 24 Swift/BAT6 host galaxies at 0.1 < z < 2 with SFR SED (or lower limits on SFR SED ) presented in Japelj et al. (2016) and Palmerio et al. (2019), we calculate a detection rate with the ngVLA of ≈ 50%, assuming R SFR ≈ 1. A detection rate much larger than these estimates would indicate that long GRB hosts have some amount of obscured star formation, and exact measurements of SFR Radio would determine the fraction of long GRB host galaxies that house significant obscured SF.

CONCLUSIONS
We have newly classified 2 long GRBs as dark (GRB 130420B and GRB 160509A) and presented VLA, CARMA, and PdBI afterglow observations of five dark GRBs (GRB 130131A, 130420B, 130609A, 131229A, and 140713A). We uniformly modeled the radio to X-ray afterglow of five dark GRBs with VLA detections (GRB 110709B, 111215A, 130131A, 131229A, and 140713A), using our afterglow modeling software that incorporates effects due to jet breaks, scintillation, and IC cooling, and include one dark GRB from the literature which was modeled using the same method and software. The radio detections allowed us to determine the environment and break frequencies of the synchrotron afterglow, in turn constraining the burst energetics, circumburst density, and geometries. Additionally, we fit the host galaxies of 6 dark GRBs (GRB 110709B, 111215A, 130131A, 131229A, 140713A and 160509A) using Prospector, and present photometric redshifts for 2 of the dark GRBs (GRB 130131A at z ≈ 1.55 and GRB 131229A at z ≈ 1.04). We come to the following conclusions: • Dark GRBs are not distinct from typical long GRBs in terms of duration, burst kinetic energy, jet opening angle, or circumburst density. However, they are statistically distinct from typical long GRBs in terms of fluence, though this distinction may be attributed to observational biases and inconsistent classification of dark GRBs.
• Our sample of six uniformly modeled dark GRBs with VLA detections have line-of-sight extinction values of A V,GRB ≈ 2.2− 10.6 mag, demonstrating the importance of radio observations in revealing GRBs with heavily dust-obscured sightlines.
These values are 0.7 to 12.8 times greater than their median host galaxy values A V,Host ≈ 0.3 − 4.7 mag.
• While dark GRBs do occur in dustier and more massive galaxies than typical long GRBs, the lineof-sight extinction is not strongly correlated to the median host extinction, nor to the circumburst density. This indicates that the origin of the dust along the line-of-sight is due to a clumpy, rather than uniform, dust distribution within the host galaxy. This also disfavors a dust origin from the extremely local (∼ parsec) environment of the burst.
• Targeted radio searches with ∼ µJy sensitivity (e.g. the ngVLA) should be capable of detecting ≈ 50% of long GRB host galaxies at 0.1 < z < 2, where a higher detection rate and exact flux measurements will determine the amount of obscured SF within long GRB hosts.
Our work demonstrates the unique power of rapidresponse radio observations with the VLA in uncovering the most obscured GRBs via their afterglows. This is especially important given that these events by definition have extinguished optical emission. Additionally, observations of dark GRB environments, from parsec to kiloparsec scales, lends insight on the distribution of dust and star formation in the galaxies which give rise to these relatively rare transients. Looking forward, next generation radio facilities, in conjunction with UVoptical observations, can be leveraged to determine the degree of obscured star formation for a large population of GRB environments across redshift. 7. ACKNOWLEDGEMENTS G.S. acknowledges for this work was provided by the NSF through Student Observing Support award SOSP20B-001 from the NRAO. The Fong Group at Northwestern acknowledges support by the National Science Foundation under grant Nos. AST-1814782, AST-1909358 and CAREER grant No. AST-2047919. W.F. gratefully acknowledges support by the David and Lucile Packard Foundation. AJCT acknowledges support from the Spanish Ministry Project PID2020-118491GB-I00, Junta de Andalucía Project P20 01068 and the "Center of Excellence Severo Ochoa" award for the Instituto de Astrofísica de Andalucía (SEV-2017-0709) The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. This work made use of data supplied by the UK Swift Science Data Centre at the University of Leicester. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. This work is based on observations carried out under project number S14DD004 with the IRAM NOEMA Interferometer. IRAM is supported by INSU/CNRS (France), MPG (Germany) and IGN (Spain). This research was supported in part through the computational resources and staff contributions provided for the Quest high performance computing facility at Northwestern University which is jointly supported by the Office of the Provost, the Office for Research, and Northwestern University Information Technology. W. M. Keck Observatory and MMT Observatory access was supported by Northwestern University and the Center for Interdisciplinary Exploration and Research in Astrophysics (CIERA). Some of the data presented herein were obtained at the W. M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California and the National Aeronautics and Space Administration. The Observatory was made possible by the generous financial support of the W. M. Keck Foundation. The authors wish to recognize and acknowledge the very significant cultural role and reverence that the summit of Maunakea has always had within the indigenous Hawaiian community. We are most fortunate to have the opportunity to conduct observations from this mountain. Some observations reported here were obtained at the MMT Observatory, a joint facility of the University of Arizona and the Smithsonian Institution. The United Kingdom Infrared Telescope (UKIRT) was supported by NASA and operated under an agreement among the University of Hawaii, the University of Arizona, and Lockheed Martin Advanced Technology Center; operations are enabled through the cooperation of the East Asian Observatory. We thank the Cambridge Astronomical Survey Unit (CASU) for processing the WFCAM data and the WFCAM Science Archive (WSA) for making the data available. This paper includes data gathered with the 6.5 meter Magellan Telescopes located at Las Campanas Observatory, Chile. The LBT is an international collaboration among institutions in the United States, Italy and Germany. The LBT Corporation partners are: The University of Arizona on behalf of the Arizona university system; Istituto Nazionale di Astrofisica, Italy; LBT Beteiligungsgesellschaft, Germany, representing the Max Planck Society, the Astrophysical Institute Potsdam, and Heidelberg University; The Ohio State University; The Research Corporation, on behalf of The University of Notre Dame, University of Minnesota and University of Virginia. IRAF is distributed by the National Optical As-tronomy Observatory, which is operated by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation. This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation.  Here we present the X-ray to radio afterglow observations of the 8 dark GRBs in our sample, as well as their host galaxy observations. For the five dark GRBs with no previously published VLA observations, we summarize the Neil Gehrels Swift Observatory (Swift) X-ray Telescope (XRT) properties in Table 1. Unless otherwise stated, all VLA data were manually reduced using standard procedures with the Common Astronomy Software Applications (CASA, McMullin et al. 2007), and all Combined Array for Research in Millimeter Astronomy (CARMA) data were manually reduced using standard procedures with the Miriad software package (Sault et al. 1995). For VLA and CARMA observations, we measure the flux density and position of the afterglow using the imtool program under the pwkit package, which fits the afterglow to a point source (Williams et al. 2017b). The radio observations, including the configuration, gain, bandpass, and flux, calibrators, are summarized in Table A2.

A.1. GRB 130131A
A.1.1. Swift and Optical Observations GRB 130131A was discovered by the Burst Alert Telescope (BAT) on-board Swift on 2013 January 31.58 (Grupe et al. 2013). The XRT started observations of GRB 130131A at δt = 58.5 s (where δt is the time after BAT trigger), finding an uncatalogued X-ray source within the BAT position (Grupe et al. 2013;Evans et al. 2013).
An uncatalogued, fading, optical/near-infrared (NIR) source was found within the XRT error circle at δt ≈ 0.02 − 0.04 days in R-, J-, and K-band, and was determined to be the optical/NIR afterglow after the source faded by 2 mag in K−band by δt ≈ 0.85 days (Volnova et al. 2013;Tanvir et al. 2013a,b). We initiated MMT SAO Widefield InfraRed Camera (SWIRC) observations at δt ≈ 0.67 days in J-and H-band and detected a source coincident with the optical/NIR afterglow (Volnova et al. 2013;Tanvir et al. 2013a,b) in both bands. For photometric calibration, we use sources in the field in common with the 2MASS catalog, and perform aperture photometry using IRAF. We find afterglow magnitudes of 21.9 ± 0.1 mag in J-band and 21.2 ± 0.1 mag in H-band at a position of R.A.=11 h 24 m 30.35 s and Dec=+48 • 04 33.08 with a positional uncertainty of 0.10 .
Interpolating the XRT light curve to the times of the R-band observations, we find β OX ≈ −0.2 at 0.016 days and 0.041 days, meeting the Jakobsson et al. (2004) criterion and confirming the darkness classification first asserted by Volnova et al. (2013). We find β X = −1.4 +0.3 −0.2 , and therefore, β OX − β X > 0.5 for the R-band detections, indicating that GRB 130131A also meets the van der Horst et al. (2009)    We initiated a series of four observations of GRB 130131A starting on 2013 February 1.28 UT with the VLA under Program 13A-046 (PI: E. Berger) at 6.0 GHz and 21.8 GHz (side-band frequencies of 19.2 GHz and 24.5 GHz). We detect a clear source within the XRT error circle with 3.6σ significance of ≈ 29.4 µJy at 6.0 GHz (the information presented here supersedes the information from a preliminary reduction reported in Laskar et al. 2013b). We took further VLA 6.0 GHz and 21.8 GHz observations until 2013 February 14.47 UT. The first 3 epochs of observations were taken with 8-bit sampling, while we used 3-bit sampling for the final 21.8 GHz observations to improve sensitivity. Due to the increased complexity of the 3-bit observations, we used the CASA VLA pipeline (version 2020.1.0.36) for the reduction (McMullin et al. 2007). The details of the VLA observations are listed in Table A2. Using the observation with the highest signal to noise (the second 21.8 GHz epoch), we derive a position for the radio afterglow of R.A.=11 h 24 m 30.401 s (±0.739") and Dec=+48 • 04 33.15 (±0.357 ).
We initiated a series of three observations of GRB 130131A starting on 2013 February 1.45 UT using CARMA under program number c0999 (PI: A. Za-uderer). We observed the field at a mean frequency of 85.5 GHz and do not detect a source to a 3σ limit of 230 µJy (the information presented here supercedes the marginal detection reported in Zauderer et al. 2013a). We took further observations until 2013 February 3.40 UT, with no significant detection. The details of the CARMA observations are listed in Table A2.

A.1.3. Host Galaxy Observations
The faint host galaxy of GRB 130131A was first reported in Chrimes et al. (2019), within the XRT position (Evans et al. 2009). We initiated griz−band observations of the host galaxy on 2015 December 3 UT with the Large Binocular Camera (LBC) on the Large Binocular Telescope (LBT/LBC) atop Mount Graham, Arizona. We reduced these images using standard routines in the IRAF/ccdproc package (Tody 1986(Tody , 1993. We applied bias and flat-field corrections and co-added the dithered images in each filter. We performed absolute astrometry with IRAF/ccmap and ccsetwcs using sources in common with the SDSS DR12 catalog (Alam et al. 2015). Using SExtractor (Bertin & Arnouts 1996), we measure a position of the host: R.A.=11 h 24 m 30.317 s and Dec =+48 • 04 33.32 . The radio and NIR afterglows are at projected offsets of 0.857 and 0.401 from the host, respectively (Figure 1).
We initiated NIR Y JHK-band observations of the host galaxy in 2016 March, 2020 December, and 2021 January using the Wide-field Camera (WFCAM; Casali et al. 2007) mounted on the 3.8-m United Kingdom Infrared Telescope (UKIRT) and the MMT and Magellan Infrared Spectrograph (MMIRS) on the MMT telescope. For UKIRT data, we obtained pre-processed images from the WFCAM Science Archive (Hamly et al. 2008), which are corrected for bias, flat-field, and dark current by the Cambridge Astronomical Survey Unit 9 . For each filter, we co-added the images and performed astrometry relative to 2MASS using a combination of tasks in Starlink 10 and IRAF. For MMIRS data, we used a custom pipeline, POTPyRI 11 , to perform bias, flat-fielding, dark corrections, sky subtraction, and coaddition.
The host galaxy is detected in all bands. We performed aperture photometry of the host using IRAF/phot with a source radius of 2.5 × θ FWHM and a background annulus immediately surrounding the host. For photometric calibration, we use sources in the field in common with the SDSS DR12 and 2MASS catalogs 9 http://casu.ast.cam.ac.uk/ 10 http://starlink.eao.hawaii.edu/starlink 11 https://github.com/CIERA-Transients/POTPyRI and employ the standard Vega to AB conversions as necessary. The host galaxy is faint, with r ≈ 23.6 mag, but with red colors (r − J ≈ 1.7 mag). Our full photometric results are listed in Table A3. Using the projected offset from the radio and optical afterglows and the rband magnitude, we calculate a probability of chance coincidence , P cc ≈ 0.003 − 0.007, establishing a robust association with GRB 130131A.
We obtained 2×1200 s of spectroscopy with the Multi-Object Double Spectrograph (MODS) mounted on the LBT on 2015 December 4 UT with the 200l grating, to cover a wavelength range of λ ≈ 3200 − 10000Å. The spectral continuum is weakly detected with an average S/N ≈ 1.2. The spectrum does not exhibit any obvious features to λ ≈ 10000Å, yielding a tentative lower limit on the redshift of z 1.3 (due to the absence of [OII]λ3727). This is consistent with the redshift upper limit of z 4 from the detection of the host in the HST F160W and F606W bands (Chrimes et al. 2019).
The host galaxy is well detected in Spitzer /IRAC 3.6 µm and 4.5 µm observations taken on 2016 July 17 (PI: Perley). We downloaded the pipeline processed post-basic calibrated data (pbcd) mosaics and performed photometry using a 3 aperture and 3 -7 background annulus, masking nearby bright sources from the background region. We applied standard aperture corrections 12 , and obtain host magnitudes of ≈ 21.2 mag and ≈ 21.0 mag at 3.6 µm and 4.5 µm, respectively.

A.2. GRB 130420B
A.2.1. Swift and Optical Observations GRB 130420B was discovered by BAT on 2013 April 20.54 (Oates et al. 2013). The XRT started observations of GRB 130420B at δt = 54.5 s, finding an uncatalogued X-ray source within the BAT position (Oates et al. 2013). Ground-based follow-up at optical wavelengths to search for the optical afterglow of GRB 130420B reached limiting magnitudes of R > 22.6 mag at a mean time of δt ≈ 0.02 days using the 2.4m Gao-Mei-Gu (GMG) telescope (Zhao et al. 2013). We interpolate the X-ray light curve to the time of the most constraining optical limit, and we calculate β OX −0.35, meeting the criterion of a dark burst as determined by Jakobsson et al. (2004). We find β X = −0.8 ± 0.2, leading to β OX −β X 0.45, which does not meet the darkness criteria of van der Horst et al. (2009), due to GRB 130420B's shallow β X . Note-All values are in AB magnitudes and are corrected for the Galactic extinction, A λ , in the direction of the burst (Schlafly & Finkbeiner 2011). a We use this z-band limit in our host galaxy modeling as it is the most constraining for this burst (see Section 4) † No host identified.

A.2.2. Radio and Millimeter Observations
We initiated observations of GRB 130420B on 2013 April 23.14 (δt ≈ 2.60 days) with CARMA under program number c0999 (PI: A. Zauderer) at a mean frequency of 85.5 GHz. We do not detect any source within the XRT error circle to a 3σ limit of < 436 µJy (the information presented here supercedes the information from a preliminary reduction reported in Zauderer & Berger 2013).
We initiated observations of GRB 130420B on 2013 April 23.32 (δt ≈ 2.78 days) with the VLA under program 13A-046 (PI: E. Berger) at 6.0 GHz and 21.8 GHz. We do not detect a source within the XRT error circle to a 3σ limit of < 22.8 µJy at 6.0 GHz and < 34.8 µJy at 21.8 GHz (the information presented here supercedes the information from a preliminary reduction reported in Zauderer & Berger 2013).

A.2.3. Host Galaxy Search
We obtained MMT/Binospec r-band observations of the field of GRB 130420B in Nov 2020. While there is a point source directly to the west of the XRT position (Evans et al. 2009), we do not detect any source within the XRT position (90% confidence) to a 3σ limit of r 24.2 mag. The Binospec image is shown in Figure 1.

A.3. GRB 130609A
A.3.1. Swift and Optical Observations GRB 130609A was discovered by BAT (Cummings et al. 2013). XRT started observations at δt = 66.8 s, finding an uncatalogued X-ray source within the BAT position (Cummings et al. 2013). Deep optical and NIR observations to search for the afterglow of GRB 130609A were taken at a mean time of δt ≈ 0.07 days, using the Reionization and Transients Infrared Camera (RATIR) in the r i ZY JH-bands, reaching limits of > 23.34 mag to > 21.06 mag in r − and H−band, respectively (Butler et al. 2013). We interpolate the X-ray light curve to the time of the most constraining optical limit (i > 23.27 mag at δt ≈ 0.07 days), and we calculate β OX 0.49, meeting the criterion of a dark burst as determined by Jakobsson et al. (2004), and confirming the classification suggested by Perley & Cenko (2013). We find β X = −1.5 ± 0.2, leading to β OX − β X 1.00, confirming GRB 130609A is also dark according to the criterion of van der Horst et al. (2009).

A.3.2. Radio and Millimeter Observations
We initiated observations of GRB 130609A on 2013 June 9.9 UT (δt ≈ 0.78 days) with the VLA under program 13A-046 (PI: E. Berger) at 6.0 GHz and 21.8 GHz. We do not detect a source within the XRT error circle to a 3σ limit of < 28.5 µJy at 6.0 GHz and < 38.4 µJy at 21.8 GHz .
We initiated observations of GRB 130609A on 2013 June 10.15 (δt ≈ 1.02 days) with CARMA under program number c0999 (PI: A. Zauderer) at a mean frequency of 85.5 GHz. We do not detect any source within the XRT error circle to a 3σ limit of < 341 µJy (the information presented here supercedes the information from a preliminary reduction reported in Zauderer et al. 2013b).

A.3.3. Host Galaxy Search
We obtained deep MMT/Binospec r-band and MMIRS J-band observations of the field of GRB 130609A in Nov and Dec 2020 to search for a host galaxy. The r-band imaging reveals a source within the XRT position (Evans et al. 2009) with r = 23.29 ± 0.11 mag; this source is also weakly detected in J-band (labeled "S1" in Figure 1). However, this source has a PSF consistent with being point-like, indicating that it is a foreground star. We identify a second, fainter source to the northwest of the XRT position (labeled "S2" in Figure 1) with r = 24.48 ± 0.14 mag that is more clearly extended and is a potential host of GRB 130609A. Thus, it is difficult to draw any strong conclusions regarding the origin or redshift of GRB 130609A.

A.4. GRB 131229A
A.4.1. X-ray and Optical Observations GRB 131229A was discovered by BAT on 2013 December 29.28 (Page et al. 2013). The XRT started observations of GRB 131229A at δt = 93.8 s, finding an uncatalogued X-ray source within the BAT position (Page et al. 2013). We initiated deep optical ground-based observations to search for the optical afterglow of GRB 131229A with the 6.5 m Magellan Clay telescope in r i z -bands at a mean time of δt ≈ 0.03 days, reaching z > 24.3 mag (Chornock et al. 2013b). We interpolate the X-ray light curve to the time of the most constraining optical limit, and we calculate β OX 0.24, meeting the criterion of a dark burst as determined by Jakobsson et al. (2004), and confirming the classification of Chrimes et al. (2019). We find β X = −1.2 ± 0.1, leading to β OX − β X 1.41, confirming GRB 131229A is also dark according to the criterion of van der Horst et al. (2009).

A.4.2. Radio and Millimeter Observations
We initiated observations of GRB 131229A on 2013 December 30.18 (δt ≈ 0.91 days) with the VLA under program 13A-541 (PI: E. Berger) at 6.0 GHz and 21.8 GHz. We do not detect a source within the XRT error circle at 21.8 GHz to a 3σ limit of < 45.6 µJy. We detect a clear source within the XRT error circle of 3.9σ significance of 74 µJy at 6.0 GHz at a position of R.A.=05 h 40 m 55.649 s (±0.093 ) and Dec=−04 • 23 47.098 (±0.115 ). We did not re-observe the field, and thus the variability of the source could not be determined, therefore we cannot definitively claim this source as the radio afterglow of GRB 131229A. To place limits on the presence of a background radio source, we searched the VLA Sky Survey (VLASS, Lacy et al. 2020), and found no source within the XRT localization to a limit of 420 µJy at 3.0 GHz (with observations taken δt ≈ 3.9 yr).
GRB 131229A was observed with CARMA at 93.0 GHz at a mean time of δt ≈ 1.00 days after the burst. No millimeter afterglow emission was found within the XRT error circle to a limit of 0.6 mJy (Perley 2013).

A.4.3. Host Galaxy Observations
The faint host galaxy of GRB 131229A was first reported in Chrimes et al. (2019) from HST/F160W imaging, and it is the only detected object within the 90% XRT position (Figure 1, Evans et al. 2009). We obtained observations of the field with LBT/LBC (griz-bands), MMT/MMIRS (Y -band), and Magellan/Fourstar (JKbands). We reduced and co-added the data in a similar manner as described before. The host galaxy is not detected in any of our optical imaging to deep limits of 23.2 − 24.9 mag, and is weakly detected in our Y , J and K-band imaging with K = 22.95 ± 0.29 mag. Overall, the host is red, with r − K 1.8 mag, is at a 0.41 offset from the VLA position and has a P cc = 1.5×10 −3 . Our full photometric results are listed in Table A3. From the K-band image, we measure a position for the host of R.A.=5 h 40 m 55.632 s and Dec=−4 • 23 46.77 . The deep non-detection of the host in the i-band, coupled with a brightness at Y -band that is 1.3 mag brighter, suggest a 4000Å break in this wavelength regime, with a redshift of z 1 − 1.5.
To place the Chandra afterglow on the host galaxy image, we perform relative astrometry using three common sources between Chandra and Magellan K-band. We obtained the Chandra observation from the archive (PI: Levan; ObsID 15195). We find a tie uncertainty of σ Magellan→Chandra = 0.13 . The corrected position is R.A. = 05 h 40 m 55.64 s and Dec. = −04 • 23 46.824 with a positional uncertainty of 0.63 (including the uncertainty in the astrometric tie, the positional uncertainty, and the absolute astrometric uncertainty of 0.6 ). Our Chandra afterglow position is consistent with the VLA C-band afterglow position and intersects the host galaxy of GRB 131229A (Figure 1).

A.5. GRB 140713A
A.5.1. Swift and Optical Observations GRB 140713A was discovered by BAT on 2014 July 13.78 (Mangano et al. 2014).The XRT began observations starting at δt = 72.8 s, detecting an uncatalogued X-ray source within the BAT position (Mangano et al. 2014;Beardmore et al. 2014). Deep optical observations to search for the optical afterglow using the AL-FOSC instrument on the 2.5-m Nordic Optical Telescope (NOT) were taken at δt ≈ 0.14 − 0.17 days, resulting in 3σ upper limits on the afterglow flux of r 24.30 mag, i 23.50 mag, and z 22.60 mag (Cano et al. 2014;Higgins et al. 2019). We interpolate the X-ray light curve to the time of the deepest optical limit (r > 24.30 mag at δt ≈ 0.15 days), and we calculate β OX −0.26, meeting the criterion of a dark burst as determined by Jakobsson et al. (2004), and confirming the classification of Higgins et al. (2019). We find β X = −1.0 ± 0.2, leading to β OX − β X 0.73, indicating GRB 140713A also meets the darkness criterion of van der Horst et al. (2009). at 6.0 GHz and F ν ≈ 0.82 mJy at 14.7 GHz. Using the observation with the highest signal to noise (the third 14.7 GHz epoch), we derive a position for the radio afterglow of R.A.=18 h 44 m 25.481 s (±0.149 ) and Dec=+59 • 38 00.69 (±0.052 ), an improvement on our CARMA position.
In addition, we initiated observations of GRB 1407134A with the Plateau deBure Interferometer (PdBI) at 86.7 GHz as part of a long-term ToO program (Program number S14DD004, PI: A. Castro-Tirado). The PdBI observed the source at six separate epochs across 2014 Jul 18-Sep 14 UT (δt ≈ 4.4 − 63.1 days). We reduced the data with the standard CLIC and MAPPING software distributed by the Grenoble GILDAS group 13 , and use the carbon star MWC349 as the flux calibrator. We detect a source in all epochs except the final epoch, at a position consistent with the millimeter and radio afterglows. The flux measurements of these observations are listed in Table A2.
To supplement our CARMA, VLA, and PdBI data, we include literature data in our subsequent analysis from the Arcminute Microkelvin Imager (AMI) Large Array (mean frequency of 15.7 GHz) (Anderson et al. 2014(Anderson et al. , 2018 Observations by the 10.4m Gran Telescopio CA-NARIAS (GTC) telescope at δt ≈ 3.1 days revealed a faint, r ≈ 24 mag source as the candidate host galaxy ). Identification and further analysis of this source was also provided in Higgins et al. (2019). We initiated griz-band observations of the host galaxy with the Large Binocular Camera (LBC) on the Large Binocular Telescope (LBT/LBC) on 2014 Oct 23 UT. We reduced and co-added the data in a similar manner as described before. We calibrated the absolute astrometry using sources in common with the Pan-STARRS1 catalog (Chambers et al. 2016) using IRAF/ccmap and ccsetwcs. We performed aperture photometry on these images using IRAF/phot using a source radius of 2.5 × θ FWHM and a background annulus 13 https://www.iram.fr/IRAMFR/GILDAS at the host position. For photometric calibration, we use standard stars from PAN-STARRS1 in the field (Chambers et al. 2016), and then convert to the SDSS system using standard transformations Tonry et al. (2012). We find a host magnitude of r = 23.8±0.2 mag, and the host galaxy photometry in all bands is listed in Table A3.
In addition, we initiated J and K s -band NIR observations of the host galaxy in Oct 2018 using the Multi-Object Spectrometer for Infra-Red Exploration (MOS-FIRE) instrument mounted on the Keck I telescope (PI: Fong), We used POTPyRI 14 to perform bias, flat-fielding, dark corrections, sky subtraction, and co-addition. We calibrated the absolute astrometry to 2MASS using IRAF/ccmap and ccsetwcs. The host galaxy is well detected in each of the NIR bands. We calibrated to 2MASS, and converted to the AB system (Blanton & Roweis 2007); The magnitude values are listed in Table A3.
We obtained 3 × 1800 s of host galaxy spectroscopy using the Low Resolution Imaging Spectrometer (LRIS) mounted on the 10m Keck I telescope on 2018 October 6 UT. We used the 400/3400 grism and 400/8500 grating in combination with the D560 dichroic for an effective wavelength range of ∼3200-10100Å. The raw frames were corrected for bias from the overscan region, flattened, and stitched together using custom methods implemented in pyraf 15 . We then extracted onedimensional spectra of the host galaxy and applied a dispersion correction derived from arc-lamp spectra. Finally, we applied an atmospheric absorption correction and flux calibration derived from a spectrum of the spectrophotometric standard GD71 obtained on the same night. We detect a faint continuum (S/N∼10) with several clear emission features. We identify [OII]λ3727, Hβλ4861, and [OIII]λ5007 at a common redshift of z = 0.935 ± 0.002. The spectrum is shown in Figure 8.
To supplement our optical and NIR host galaxy observations, we include 3.6 and 4.5µm observations from Spitzer/IRAC photometry, published in Higgins et al. (2019) (Table A3). Our grizJHK s data, as well as the available Spitzer photometry are subsequently used in our host galaxy modeling.
A.6.2. GRB 111215A We gather radio and optical afterglow observations of GRB 111215A from Zauderer et al. (2013c) and van der Horst et al. (2015), and X-ray afterglow observations from Swift. GRB 111215A has previously been determined to be a dark burst through the β OX > −0.5 criterion (Jakobsson et al. 2004) in Zauderer et al. (2013c). We gather host galaxy observations of GRB 111215A from van der Horst et al. The radio, optical, and X-ray afterglow of GRB 160509A has previously been modeled with our modeling framework in Laskar et al. (2016), and we utilize the afterglow parameters within for our discussion (Section 5). The optical afterglow of GRB 160509A was heavily obscured, with a line-of-sight extinction of A V,GRB ≈ 3.4 mag, indicating it is likely a dark burst. Using the optical afterglow detections in r -band at δt ≈ 0.25 days (Laskar et al. 2016), we calculate β OX ≈ 0.03, meeting the darkness criterion of Jakobsson et al. (2004), and classifying GRB 160509A as a dark burst. Additionally, we find β X = −1.0 ± 0.1, leading to β OX − β X ≈ 0.98. Thus, GRB 160509A is also classified as a dark burst according to the van der Horst et al. (2009) criterion.
The host galaxy of GRB 160509A was previously observed with HST/WFC3 in the F110W and F160W filters on 2017 July 5 (PI: Kangas) (Kangas et al. 2020). We retrieved, aligned, and drizzled the individual exposures for each band using the HST reduction pipeline hst123 (Kilpatrick et al. 2022;Kilpatrick 2021). We also included a sky subtraction step as part of astrodrizzle (Hack et al. 2021) to remove large-scale background emission near the host galaxy. Using an elliptical aperture and the tabulated HST zeropoints, we calculated the F110W and F160W brightness of the host galaxy (Table A3 and Figure 1).
The host galaxy is also clearly detected in Spitzer /IRAC 3.6 µm observations taken on 2017 Nov 05 (PI: Perley). We downloaded and reduced the individual basic calibrated data frames (i.e., cbcd) using photpipe (Rest et al. 2005;Kilpatrick et al. 2018), including alignment, flux calibration, and optimal stacking to a pixel scale of 0.3 /pix. We performed final PSF photometry in the stacked frames using a custom version of DoPhot (Schechter et al. 1993). The host galaxy appears point-like in the stacked IRAC frames and is clearly separated from a nearby galaxy of similar brightness seen in the HST frames ( Figure 1). Therefore, we use the DoPhot photometry calibrated using the appropriate warm Spitzer sensitivity function for IRAC and obtain a magnitude of ≈ 19.6 mag at 3.6 µm (Table A3).
We supplement these host galaxy observations with the Keck/LRIS g -and r -band measurements taken at δt ≈ 28.2 days, where afterglow modeling found the host galaxy to dominate the flux (Laskar et al. 2016; Table A3).

B. ALTERNATIVE MCMC AFTERGLOW MODELS
Here we present the alternative MCMC afterglow models for the bursts whose afterglow observations did not distinguish between a wind or ISM environment through preliminary analytical arguments (GRB 110709B, 130131A, and 140713A).

B.1. GRB 110709B
We present the best fit (highest likelihood) wind environment model for GRB 110709B in Figure B1 and list the model parameters as well as the summary statistics from the marginalized posterior density functions (medians and 68% credible intervals) in Table 2 in Table A4.
Our best fit wind model for GRB 110709B finds a later t jet ≈ 2.4 days compared to our ISM model (Section 3.1.2). While the wind model better fits the X-ray light curve at δt 5 days, the marginal statistical preference for the ISM model is a result of the ISM model better matching the X-ray light curve at δt 5 days, where the majority of the X-ray data exist. The best fit parameters for our wind model are the same as Kangas & Fruchter (2021) within 2σ, though our E K and A * are orders of magnitude different than those found by Zauderer et al. (2013c), likely due to the inclusion of IC effects in our model. Our best fit wind model requires A V,GRB 5.01 mag to match the GROND optical/NIR afterglow upper limits.

B.2. GRB 130131A
We present the best fit wind environment for GRB 130131A model in Figure B2 and list the model parameters as well as the summary statistics from the marginalized posterior density functions (medians and 68% credible intervals) in Table A4.
Our ISM model of GRB 130131A finds an earlier jet break time of t jet ≈ 0.8 days compared to our wind model. As a consequence, while the radio afterglow of GRB 130131A is better fit by the ISM model, the Xray light curve is under-predicted by the ISM model at

a
Note-The top row for each parameter corresponds to the best fit forward shock value from our MCMC modeling. The bottom row for each parameter corresponds to the summary statistics from the marginalized posterior density functions (medians and 68% credible intervals) a A V,GRB value from afterglow model using median values δt 1 day, resulting in the statistical preference for the wind model of GRB 130131A. Our ISM model finds that the extinction required to match the optical and NIR observations is A V,GRB ≈ 2.41 mag.

B.3. GRB 140713A
We present the best fit ISM environment model for GRB 140713A in Figure B3 and list the model parameters as well as the summary statistics from the marginalized posterior density functions (medians and 68% credible intervals) in Table 2 in Table A4.
Similar to the wind environment model (Section 3.5.2), the X-ray light curve for the ISM model of GRB 140713A is under-predicted at δt 0.7 days. While the K u -band observations are better fit by the ISM model, the significant statistical preference for the wind model is attributed to the better fit C-band and 3mm observations light curves. Our ISM model for GRB 140713A requires A V,GRB 3.12 mag necessary to match the optical limits. This limit is the same as that derived by Higgins et al. (2019) for an SMC-like galactic extinction model.  Figure B1. X-ray, optical, radio afterglow light curves of GRB 110709B, together with the best-fit forward shock model in an wind environment (lines). Circles represent literature data , and triangles represent 3σ upper limits. Open symbols indicate data that are not included in the fit, and shaded regions represent variability due to scintillation. The inset shows the model Ks-band and Hband light curves (solid lines) as well as the non-extinguished models (dashed lines), indicating AV,GRB 5.0 mag to explain the upper limits.  Figure B2. Left: X-ray, optical, millimeter and radio afterglow light curves of GRB 130131A, together with the best-fit forward shock model in an ISM environment (lines). Squares represent data newly reported here, and triangles represent 3σ upper limits.
Open symbols indicate data that are not included in the fit, and shaded regions represent variability due to scintillation. The inset shows the model J-and R-band light curves (solid lines) as well as the non-extinguished models (dashed lines), indicating AV,GRB ≈ 2.4 mag to explain the detections. Right: Radio and millimeter spectral energy distributions (SEDs) of the afterglow of GRB 130131A from δt ≈ 0.7-13.9 days, together with the best-fit forward shock model (lines).  Figure B3. Left: X-ray, optical, millimeter and radio afterglow light curves of GRB 140713A, together with the best-fit forward shock model in an ISM environment (lines). Squares represent data newly reported here, circles represent literature data (Higgins et al. (2019)), and triangles represent 3σ upper limits. Open symbols indicate data that are not included in the fit, and shaded regions represent variability due to scintillation. The inset shows the model r-, i-, and z-band light curves (solid lines) as well as the non-extinguished models (dashed lines), indicating AV,GRB 3.1 mag to explain the upper limits. Right: Radio and millimeter spectral energy distributions (SEDs) of the afterglow of GRB 140713A from δt ≈ 3.0-31.5 days, together with the best-fit forward shock model (lines).