A tale of two mergers: constraints on kilonova detection in two short GRBs at z$\sim$0.5

We present a detailed multi-wavelength analysis of two short Gamma-Ray Bursts (sGRBs) detected by the Neil Gehrels Swift Observatory: GRB 160624A at $z=0.483$ and GRB 200522A at $z=0.554$. These sGRBs demonstrate very different properties in their observed emission and environment. GRB 160624A is associated to a late-type galaxy with an old stellar population ($\approx$3 Gyr) and moderate on-going star formation ($\approx$1 $M_{\odot}$ yr$^{-1}$). Hubble and Gemini limits on optical/nIR emission from GRB 160624A are among the most stringent for sGRBs, leading to tight constraints on the allowed kilonova properties. In particular, we rule out any kilonova brighter than AT2017gfo, disfavoring large masses of wind ejecta ($\lesssim$0.03 $M_\odot$). In contrast, observations of GRB 200522A uncovered a luminous ($L_\textrm{F125W}\approx 10^{42}$ erg s$^{-1}$ at 2.3~d) and red ($r-H\approx 1.3$ mag) counterpart. The red color can be explained either by bright kilonova emission powered by the radioactive decay of a large amount of wind ejecta (0.03 $M_\odot$ $\lesssim$ $M$ $\lesssim$ 0.1 $M_\odot$) or moderate extinction, $E(B-V)\approx0.1-0.2$ mag, along the line of sight. The location of this sGRB in the inner regions of a young ($\approx$0.1 Gyr) star-forming ($\approx$2-4 $M_{\odot}$ yr$^{-1}$) galaxy and the limited sampling of its counterpart do not allow us to rule out dust effects as contributing, at least in part, to the red color.

Kilonova emission following a NS-NS merger originates from the radioactive decay of freshly synthesized r-process elements in neutron-rich matter surrounding the remnant compact object (Li & Paczyński 1998;Metzger et al. 2010; Barnes & Kasen 2013;Tanaka & Hotokezaka 2013;Grossman et al. 2014;Kasen et al. 2017). Kilonovae are hallmarked by "blue" thermal emission within a day of merger (e.g., AT2017gfo) which fades and gives way to the "red" and near-infrared (nIR) emission, persisting for roughly a week post-merger. Neutron-rich material (electron fraction < 0.25) composed of high-opacity lanthanides produces the red component, while the blue component results from ejecta material with higher electron fractions (Barnes & Kasen 2013;Kasen et al. 2015Kasen et al. , 2017Tanaka et al. 2017;Wollaeger et al. 2018Wollaeger et al. , 2019Fontes et al. 2020;Even et al. 2020;Korobkin et al. 2020). Dynamical ejecta, tidally stripped from the approaching neutron star(s), primarily contributes to the red component. In addition, a portion of the matter that congregates in an accretion disk surrounding the remnant compact object is released as wind ejecta (Metzger et al. 2008;Dessart et al. 2009;Lee et al. 2009;Fernández & Metzger 2013;Perego et al. 2014;Just et al. 2015;Miller et al. 2019). Ejecta in the disk supports a wide range of electron fractions, enhancing either a blue or red kilonova. The identity of the merger remnant influences the disk, with a longer-lived high mass neutron star remnant increasing the electron fraction of disk ejecta (Kasen et al. 2015). The range of electron fractions of ejecta predicted from models of disk winds varies with the implementation of the neutrino transport (Miller et al. 2019).
Kilonovae exhibit near-isotropic emission, with viewing-angle dependent variations based on ejecta morphology ) and lanthanide curtaining effects (Kasen et al. 2015). Observations of AT2017gfo were possible thanks to the particular geometry of GW170817, whose relativistic jet was misaligned with respect to our line of sight (Troja et al. 2017;Lazzati et al. 2017;Mooley et al. 2018;Ghirlanda et al. 2019;Lamb et al. 2019). As a consequence, its afterglow appeared at later times and remained relatively dim, allowing for a complete view of the kilonova. Had the same event been seen closer to its jet's axis (on-axis), the GRB afterglow would have outshined any kilonova emission.
The majority of sGRBs are discovered at much larger distances than GW170817 and are observed close to their jet's axis (Ryan et al. 2015;Beniamini & Nakar 2019;Wu & MacFadyen 2019). Their bright broadband afterglow is often the dominant emission component, which complicates the identification of any associated kilonova (see, e.g., Yang et al. 2015;Ascenzi et al. 2019). Whereas the range of luminosities and timescales of kilonova emission largely overlaps with standard GRB afterglows, the red color of a kilonova (Barnes & Kasen 2013;Tanaka & Hotokezaka 2013), much redder than any typical afterglow, is one of its distinctive features. Nonetheless, even the color information may be insufficient for an unambiguous identification. A counterpart with unusually red colors was found for the short GRB 070724A ( -≈4; Berger et al. 2009) and GRB 071227 ( -≈1.5; Eyles et al. 2019) and, in both cases, attributed to dust effects at the GRB site. The rapid timescales and high luminosity (≈ 10 43 erg s −1 ) of these two sources did not match the predictions of a radioactive-powered kilonova, although they could fit within the expected range for a magnetar-powered kilonova (Yu et al. 2013).
Densely sampled multi-color observations, extending to the nIR range, proved to be essential in the identification of the kilonova candidates GRB 130603B ) and GRB 160821B . The counterpart of GRB 130603B was identified within the spiral arm of its bright host galaxy. The source appeared unusually red, in part due to significant presence of dust along the sightline ( ∼1 mag; de Ugarte Postigo et al. 2014), and was seen to evolve over the course of time, from -≈1.7 ± 0.15 at about 14 hr to -> 2.5 at about 9 d. GRB 160821B was instead located in the outskirts of a nearby spiral galaxy, and its counterpart was also identified as unusually red ( -≈1.9; Kasliwal et al. 2017a). A detailed modeling of the X-ray and radio afterglow was able to disentangle the presence of an additional emission component in the optical and nIR data, slightly less luminous than AT2017gfo and with similar timescales and color evolution . For both GRB 130603B and GRB 160821B, a good spectral sampling over multiple epochs was a fundamental ingredient to distinguish the kilonova candidate from the underlying bright afterglow.
In addition to these candidate kilonovae, there are a a number of claimed kilonova detections based on an optical excess, e.g., GRB 050709 (Jin et al. 2016), GRB 060614 (Yang et al. 2015), GRB 070809 (Jin, et al. 2020), GRB 080503 (Perley et al. 2009), and GRB 150101B . The situation for these events is less clear due to the lack of deep nIR observations, critical to distinguish kilonova emission from standard afterglow.
The number of sGRBs with well characterized afterglows and sensitive nIR observations is still restricted to a handful of cases, and for the majority of sGRBs no meaningful limits on the ejecta mass and composition can be placed (see, e.g., Gompertz et al. 2018;Rossi et al. 2020). In this work, we continue filling this observational gap by presenting a detailed multi-wavelength study of two distant sGRBs: GRB 160624A at = 0.483 and GRB 200522A at = 0.554. We complement the early Swift data with deep Chandra observations in order to characterize the afterglow temporal evolution up to late times, and use deep Gemini and HST imaging to search for kilonova emission. The paper is organized as follows. In §2, we present the observations and data analysis for GRBs 160624A and 200522A. In §3, we describe the methods applied for our afterglow and kilonova modeling, as well as the galaxy spectral energy distribution (SED) fitting procedure. The results are presented in §4, and our conclusions in §5. For each GRB, we provide the time of observations relative to the BAT trigger time, 0 , in the observer frame. All magnitudes are presented in the AB system. We adopt the standard Λ-CDM cosmology with parameters 0 = 67.4, Ω = 0.315, and Ω Λ = 0.685 (Planck Collaboration et al. 2018). All confidence intervals are at the 1 level and upper limits at the 3 level, unless otherwise stated. Throughout the paper we adopt the convention ∝ − − .

X-ray Observations
The Swift X-ray Telescope (XRT; Burrows, et al. 2005) began observing at 0 + 59 s and localized the X-ray afterglow at a position 1 RA, DEC (J2000) = 22 ℎ 00 46 .21, +29°38 37.8 with an accuracy of 1.7 (90% confidence level (CL); Evans et al. 2007Evans et al. , 2009). Data were collected in Window Timing (WT) mode during the first 150 s and, as the source rapidly decreased in brightness, in Photon Counting (PC) mode. A deeper observation was carried out at 0 + 8.5 d (PI: Troja; ObsId 18021) by the Chandra X-ray Observatory (ACIS-S3), but no X-ray counterpart was detected. We describe the observed temporal decay (Figure 7) with a broken power-law, consisting of two segments with ∝ − . The initial decay is 1 = 0.6 ± 0.3 which steepens to 2 = 4.0 ± 0.3 after break ∼ 140 s.
The Chandra data were re-processed using the CIAO 4.12 data reduction package with CALDB Version 4.9.0, and filtered to the energy range 0.5 − 7 keV. We corrected the native Chandra astrometry by aligning the image with the Gaia Data Release 2 (Gaia Collaboration, et al. 2018). We utilized CIAO tools to extract a count-rate within the XRT error region (1.7 source aperture radius), utilizing nearby source-free regions to estimate the background. We detect zero counts in the source region with an estimated background of 0.6 counts yielding a 3 upper limit of 1.2 × 10 −4 cts s −1 (Kraft et al. 1991). We convert this rate to 1.9 × 10 −15 erg cm −2 s −1 in the 0.3-10 keV band using the best fit spectral parameters. The derived X-ray fluxes are reported in Table 1.

Optical/nIR Imaging
The Ultra-Violet Optical Telescope (UVOT; Roming et al. 2005) on-board Swift began observations at 0 + 77 s although no optical afterglow is identified within the XRT position to ℎ ≥ 21.8 AB mag (de Pasquale & D'Ai 2016). The field was imaged with the Gemini Multi-Object Spectrograph (GMOS; Hook et al. 2004) on the 8.1-meter Gemini North telescope (PI: Cucchiara) starting at 0 + 31 min. An initial 180 second r-band exposure led to the identification of a candidate host galaxy within the XRT error region (SDSS J220046.14+293839.3; Cucchiara & Levan 2016); for further discussion of the host association see §4.1.3. This was followed by deeper observations (900 s and 1440 s, respectively) at 0 + 1 hr and 0 + 1 d. Seeing during the observations was ∼ 0.5 with mean airmass 1.1 and 1.0, respectively. We retrieved the data from the Gemini archive. Data were analyzed following standard Optical/nIR fluxes were corrected for Galactic extinction due to interstellar reddening ( − ) = 0.06 mag (Schlafly & Finkbeiner 2011) using the extinction law by Fitzpatrick (1999); Indebetouw et al. (2005). X-ray fluxes were corrected for Galactic absorption = 9.14 × 10 20 cm −2 (Willingale et al. 2013), and converted into flux densities at 1 keV using photon index Γ = 1.76.
CCD reduction techniques and using the Gemini IRAF 2 reduction package. 2 IRAF is distributed by the National Optical Astronomy Observatory, which is operated by the Association of Universities for Research in Astronomy (AURA) under cooperative agreement with the National Science Foundation (NSF).
At later times, we performed two epochs of observations (PI: Troja; ObsId 14357) with the Hubble Space Telescope (HST) Advanced Camera for Surveys (ACS) Wide Field Camera (WFC) and Wide Field Camera 3 (WFC3) in Infrared (IR). See Table 1 for a log of observations. The HST data were processed using the  sndrizpipe 3 pipeline, which makes use of standard procedures within the DrizzlePac package, in order to align, drizzle, and combine exposures. The final image pixel scale was 0.09 /pix for WFC3 (i.e., 125 and 160 ) and 0.04 /pix for ACS (i.e., 606 ). We identify no candidate counterpart within the XRT localization in either Gemini or HST images. Since the XRT localization overlaps significantly with the candidate host galaxy (see Figure 2), we performed image subtraction between epochs using the High Order Transform of Psf ANd Template Subtraction code (HOTPANTS; Becker 2015) to search for transient sources embedded within the host galaxy's light. Due to the short time delay between Gemini epochs, our analysis may not reveal a slowly evolving transient. We therefore verified our results using the late ( 0 + 8.2 d) HST/ 606 3 https://github.com/srodney/sndrizpipe image as the template. Furthermore, as kilonovae can dominate at either early or late times, depending on the composition of the ejecta, we performed image subtraction between the HST epochs using each epoch as a template image. No significant residual source was uncovered in either the Gemini or HST difference images at any epoch, as shown in Figure 2.
In order to determine the upper limit on a transient source in these images, we injected artificial point-like sources within the XRT position and performed image subtraction to detect any residual signal. Gemini magnitudes were calibrated to nearby Sloan Digital Sky Survey Data Release 12 (SDSS; Fukugita et al. 1996) stars. HST magnitude zeropoints were determined with the photometry keywords obtained from the HST image headers, and were corrected with the STScI tabulated encircled energy fractions. The upper limits derived for the field are presented in Table 1. Upper limits within the galaxy's light are shallower by 0.3-0.5 mag.
Finally, we obtained imaging of the candidate host galaxy on May 20, 2018 ( 0 + 1059 d) with the Large Monolithic Imager (LMI) mounted on the 4.3-meter Lowell Discovery Telescope (LDT) in Happy Jack, AZ. Observations were taken in the griz filters, with seeing ∼ 1.65 at a mean airmass of ∼ 1.3. We applied standard procedures for reduction and calibration of these images. We obtain the galaxy apparent magnitude in each filter using the SExtractor MAG_AUTO parameter, which utilizes Kron elliptical apertures (Bertin & Arnouts 1996). Magnitudes were calibrated to nearby SDSS stars, and reported in Table 2. Near-infrared imaging in the Y bands was carried out on July 25, 2020 with the Near-Infrared Imager (NIRI; Hodapp et al. 2003) on the 8-m Gemini North telescope. Data were reduced using standard procedures within the DRAGONS package. The photometry was calibrated to nearby sources from the Panoramic Survey Telescope and Rapid Response System (Pan-STARRS; Chambers et al. 2016) and Two Micron All Sky Survey (2MASS; Skrutskie et al. 2006) for the Yand -band images, respectively. We used the offsets from Blanton & Roweis (2007) to convert 2MASS Vega magnitudes to the AB system.

Optical Spectroscopy
A spectrum of the candidate host galaxy was obtained using Gemini/GMOS (PI: Cucchiara) starting at 0 + 46 min. GMOS was configured with the R400 grating at a central wavelength of 600 nm. We reduced and analyzed the data using the Gemini IRAF package (v. 1.14). The resulting spectrum is shown in Figure 3. Emission features observed at obs ≈ 5533, 7211 and 7428 Å associated with the [OII] doublet, H , and [OIII] 5008 transitions, respectively, yield a redshift = 0.4833 ± 0.0004 in agreement with the preliminary estimate of Cucchiara & Levan (2016). At this redshift we also observe a low significance feature at the expected location of H . Line properties were derived by fitting the lines with Gaussian functions using the specutils package in Python.

Radio Observations
Radio observations were carried out with the Karl J. Jansky Very Large Array (JVLA) starting at ∼ 0 + 1 d (PI: Berger; project code: 15A-235) with the array in the B configuration. The observations were taken in the X-band, with a central frequency 10 GHz and a bandwidth of 2 GHz. The time on source was 47 minutes. Data were downloaded from the National Radio Astronomical Observatory (NRAO) online archive, and processed locally with the JVLA CASA pipeline v1.3.2 running in CASA v4.7.2. Galaxies 3C48 and J2203+3145 were used as primary and phase calibrators, respectively. We do not detect a radio transient coincident with the enhanced XRT position with a flux density upper limit < 18 Jy.

X-ray Observations
Swift/XRT observations were delayed due to the South Atlantic Anomaly, and began at 0 + 406 s . The X-ray counterpart was detected at RA, DEC (J2000) = 00 ℎ 22 43 .7, −0°16 59.4 with an accuracy of 2.2 (90% CL) 5 . XRT followup observations lasted 3 days for a total exposure of 17.5 ks in PC mode. We performed two ToO observations (PI: Troja; ObsIds 22456, 22457, and 23282) with Chandra/ACIS-S3 in order to track the late-time evolution of the X-ray lightcurve. During the first epoch ( 0 + 5.6 d), we detect the X-ray afterglow at RA, DEC (J2000) = 00 ℎ 22 43 .74, −0°16 57.53 with accuracy 0.5 , consistent with the XRT enhanced position. A second bright X-ray source lies ∼ 10.4 from the GRB position, and is coincident with the known quasar SDSS J002243.61-001707.8 at redshift = 1.44862± 0.00079 (Krawczyk et al. 2013). Due to their proximity, the two sources are not resolved in the XRT images and both contribute to the observed X-ray flux.
Analysis of the Swift and Chandra data was performed using Optical/nIR fluxes were corrected for Galactic extinction due to interstellar reddening ( − ) = 0.02 mag (Schlafly & Finkbeiner 2011) using the extinction law by Fitzpatrick (1999); Indebetouw et al. (2005). X-ray fluxes were corrected for Galactic absorption = 2.9 × 10 20 cm −2 (Willingale et al. 2013), and converted into flux densities at 1 keV using photon index Γ = 1.45. the methods described in §2.1.2. We use our Chandra observations to characterize the nearby quasar, and estimate its contribution to the measured Swift/XRT flux. Using XSPEC, we derive a photon index Γ = 1.53 ± 0.14 (C-stat 162 for 156 dof). This yields a flux = (6.6 ± 0.6) × 10 −14 erg cm −2 s −1 (0.3-10 keV). To constrain the impact of this second source on the XRT observations, we folded the quasar spectrum with the XRT response function to obtain the expected count rate with Swift/XRT, (1.5 ± 0.2) × 10 −3 cts s −1 (0.3-10 keV). We subtract this mean count rate from all XRT observations.

Optical/nIR Imaging
The Swift/UVOT began settled observations in the wh filter at 0 + 448 s (Kuin et al. 2020). Subsequent observations were performed in all optical and UV filters. There was no source detected within the enhanced XRT position. The observations were analyzed using HEASoft v6.27.2. The photometry was performed using circular apertures with a 3 radius and calibrated using the standard UVOT zeropoints (Breeveld et al. 2011); see Table 3.
We imaged the field of GRB 200522A with GMOS on the 8-m Gemini North telescope (PI: Troja). A first set of -band exposures was obtained at 0 +2.1 d under poor weather conditions (Dichiara et al. 2020a), and repeated at 0 + 3.1 d. A last observation at 0 + 9.1 d serves as a template for image subtraction. The identification of a counterpart is complicated by the presence of a bright galaxy. Image subtraction, using HOTPANTS, between the second ( 0 + 3.1 d) and third ( 0 + 9.1 d) epochs finds a weak (≈3 ) residual source within the Chandra localization. By performing aperture photometry on the difference image we estimate a magnitude of = 26.0 ± 0.4 AB, calibrated against nearby SDSS stars.
Near-IR imaging was carried out with the HST/WFC3 using the 125 and 160 filters at three epochs (PI: Berger; ObsId 15964): 0 + 3.5, 0 + 16.3, and 0 + 55.2 d. The data was processed using sndrizpipe to a final pixel scale of 0.06 /pix. Image subtraction, using HOTPANTS, between the first ( 0 + 3.5 d) and third ( 0 + 55.2 d) epoch uncovers a significant residual source in both filters, at a location consistent with the optical and X-ray positions (see Figure 5). The absolute position of the nIR transient is RA, DEC (J2000) = 00 ℎ 22 43 .737, −0°16 57.481 with a 1 un-   Table 3. There are no significant residuals detected at the afterglow location in the 125 difference image between the second ( 0 + 16.3 d) and third epochs. Following the procedure outlined in §2.1.3, we inject artificial point sources to determine an upper limit 125 > 27.2 AB mag at the afterglow position.
An independent analysis of the HST data was recently reported by Fong et al. (2020), confirming our detection of the optical/nIR counterpart. The analysis of Fong et al. (2020) reports a source brighter by ≈ 0.4 mag in the 125 filter, which we can reproduce by using a different template image derived by combining the two epochs at 0 + 16 and 55 days using the astrodrizzle package. However, some of our models (see §4.2) still predict a weak signal at 0 + 16 d, and small variations in the nearby galaxy's nucleus may also influence the photometry and the measured color. Therefore, in our work we adopt the results derived using the late time epoch at 0 + 55 d, available for both the 160 and 125 filters, and verify that all our conclusions also hold for the alternative result of a slightly brighter transient (see §4.2.1).
Late-time optical and nIR images were acquired to characterize the host galaxy's properties. Observations were carried out in the ugiz filters with the LDT/LMI on July 30, 2020, and in the YK filters with Gemini/NIRI on July 17, 2020. Data were reduced following the same procedures described in §2.1.3, and photometry was calibrated using nearby sources from SDSS and the United Kingdom Infrared Telescope Infrared Deep Sky Survey (UKIDSS; Lawrence et al. 2007). The results are listed in Table 4.

Optical Spectroscopy
We obtained a spectrum of the putative host galaxy using GMOS on the 8-m Gemini North telescope (PI: Troja) on June 27, 2020 ( 0 +36 d). We performed a set of 6×600 s exposures using the R400 grating with central wavelength ≈740 nm. The resulting combined spectrum is shown in Figure 6. We identify emission features at obs ≈ 5795, 6747, 7556, 7782, 7708 Å which are associated with the [OII] doublet, H , H , [OIII] 4959 , and [OIII] 5008 transitions, respectively. This yields a redshift = 0.5541 ± 0.0003, in agreement with our preliminary estimate (Dichiara et al. 2020b). Line properties were derived through the methods outlined in §2.1.4.

Afterglow Modeling
We model the observed afterglows within the standard fireball model (Mészáros & Rees 1997;Sari et al. 1998;Wĳers & Galama 1999;Granot & Sari 2002), described by a set of five parameters: the isotropic-equivalent kinetic energy 0 , the circumburst density 0 , the fraction of energy in magnetic fields and in electrons , and the slope of the electron energy distribution ( ) ∝ − . We assume that the environment surrounding the binary merger has a uniform density profile, consistent with the interstellar medium (ISM). Three more parameters account for the outflow's collimated geometry: the jet core width , the observer's viewing angle , and the jet's angular profile. We apply two angular profiles: (i) a uniform (tophat) jet profile with ≈ 0 and (ii) a Gaussian function in angle from the core described by ( ) = 0 exp(− 2 /2 2 ) for ≤ , where is the truncation angle of the Gaussian wings. The beaming corrected kinetic energy, is given by for a Gaussian angular profile. We include the effect of intrinsic dust extinction assuming a Fitzpatrick (1999) reddening law parametrized by = / ( − ) = 3.1. We utilize a Bayesian fitting method in conjunction with the afterglowpy software 6 , described in Ryan et al. (2020), to determine the GRB jet parameters. We apply the (version 2. We compare models by evaluating their predictive power with the Widely Applicable Information Criterion (WAIC, also known as the Watanabe-Akaike Information Criterion; Watanabe 2010). The WAIC score is an estimate of a model's expected log predictive density (elpd): roughly the probability that a new set of observations would be well described by the model's fit to the original data. A model with a high elpd (and WAIC score) produces accurate, constraining predictions. A model with a low elpd (and WAIC score) produces predictions which are inaccurate or not constraining. This naturally penalizes over-fitting: over-fit models typically show wide variability away from the observations which leads to poor predictive power and a low WAIC score. The WAIC score estimates the elpd by averaging the likelihood of each observation over the entire posterior probability distribution. As a model comparison tool, this incorporates the uncertainties in the model parameters (unlike the reduced 2 or the Akaike Information Criterion) and does not require the posterior probability distribution to be normal (unlike the the Deviance Information Criterion).
We compute the WAIC score following Gelman et al. (2013) using the recommended WAIC 2 estimate for the effective number of parameters. The WAIC score is computed at every datapoint and the total WAIC score is the sum of these contributions. In order to compare the WAIC score of two different models, we compute the standard error, , of their difference, ΔWAIC elpd (Vehtari et al. 2015). One model is favored over another if it has a higher overall WAIC score and the difference between the WAIC scores ΔWAIC elpd is significantly larger than its uncertainty . As can underestimate the true standard deviation of ΔWAIC elpd by up to a factor of two (Bengio & Grandvalet 2004), we report a conservative significance on the WAIC score difference using ΔWAIC ≈ 2 .

Empirical Constraints
Optical and infrared observations constrain the properties of possible kilonovae associated with each GRB. In the case of GRB 160624A, no optical or nIR counterpart was detected and kilonova models are directly constrained by our photometric upper limits. GRB 200522A presents instead a complex phenomenology, characterized by a bright X-ray afterglow and an optical/nIR counterpart with a red color ( − ≈1.3, Table 3). Such red color could be the result of dust along the sightline or the telltale signature of a kilonova. In order to constrain the contribution of the latter, we add to our afterglow fit an additional thermal component. Korobkin et al. (2020) demonstrated that a simple analytical fit to kilonova lightcurves can lead to order of magnitude uncertainties in the inferred ejected mass, depending on the unknown geometry of the system. We therefore use a different approach which combines empirical constraints and detailed numerical simulations. At the time of our optical/nIR observations (∼ 3.5 d), the kilonova component is roughly described by a simple blackbody. Therefore, we use a parameterized blackbody component, included in addition to the standard forward shock (FS) emission, to determine the possible thermal contribution from a kilonova. The range of fluxes allowed by this fit is then compared to an extensive suite of simulated kilonova lightcurves ( §3.2.2) in order to derive the ejecta properties. The blackbody component is described by two parameters, its temperature and emission radius , with uniform priors between [0-8,000 K] and [0-5×10 15 cm], respectively. The upper limit on radius is chosen to prevent a superluminal expansion velocity. Since our optical and nIR data are nearly coeval, we assume that there is negligible temporal and spectral evolution between observations.

Constraints on Kilonova Ejection Properties
We compare optical and infrared constraints to simulated kilonova lightcurves of varying input parameters, consistent with a wide range of plausible ejecta morphologies, compositions, masses, and velocities. For this study, we use a grid of two-component kilonova models from the LANL group (Wollaeger et al., in prep). This data set was previously used in Thakur et al. (2020) to constrain ejecta parameters for GW190814. Simulations include timedependent spectra, as early as three hours post-merger, which are subsequently converted to lightcurves for various filters. We simulate kilonovae with SuperNu (Wollaeger & van Rossum 2014), a multi-dimensional, multi-group Monte Carlo transport code, which has previously been used in a wide range of kilonova studies (Wollaeger et al. 2018(Wollaeger et al. , 2019Even et al. 2020;Korobkin et al. 2020;Thakur et al. 2020). Our simulations rely on the WinNet nucleosynthesis network ) to simulate heating from radioactive decay in addition to the latest LANL opacity database ). We consider a full set of lanthanide opacities, while uranium acts as a proxy for all actinide opacities.
We model kilonovae with two ejecta components: dynamical ejecta including heavy r-process elements and wind ejecta emanating from the resultant accretion disk surrounding the remnant compact object. We consider two disparate wind models, representing ejecta with either high-latitude ( = 0.37) or mid-latitude ( = 0.27) compositions, both having negligible lanthanide contributions. Wind ejecta assumes either a spherical or "peanut-shaped" morphology, while dynamical ejecta remains toroidal. These models correspond to the TS and TP morphologic profiles in Korobkin et al. (2020).
The grid of models includes a range of mass and velocity parameters, in addition to two morphologies and two wind compositions. Both the dynamical and wind ejecta span five possible masses: 0.001, 0.003, 0.01, 0.03, 0.1 . We ascribe one of three possible velocity distributions to both the dynamical and wind ejecta components, with median velocities of either 0.05 , 0.15 , or 0.3 , corresponding to maximum ejecta velocities of 0.1 , 0.3 , and 0.6 . The grid spans the anticipated range of ejecta properties expected from numerical simulations and observations of GW170817 (Korobkin et al. 2012;Kasen et al. 2017;Côté et al. 2018;Metzger 2019;Krüger & Foucart 2020). Each multi-dimensional simulation computes kilonova emission for 54 different polar viewing angles. Our axisymmetric simulations report spectra and lightcurves for separate viewing angles, distributed uniformly in sine of the polar angle from edge-on, on-axis viewing angles (0°) to edge-off (180°). Including all viewing angles and kilonova properties, we have 48,600 different sets of time-dependent spectra to compare to optical and infrared observations. We limit our simulation grid to kilonovae observed on-axis, when considering optical and infrared observations in conjunction with a GRB, observed with viewing angle ≈ 0. Our on-axis simulations consider viewing angles from 0°to 15.64°. We then compare the 900 on-axis kilonova simulations to optical and nearinfrared observations. We restrict plausible kilonova parameters to the range of simulated lightcurves consistent with observations, either lightcurves dimmer than the measured upper limits (see §4.1.2) or lightcurves residing in the range of inferred kilonova emission from analytic fits (see §4.2.5).

Galaxy SED Modeling
We compare the observed photometry to a range of synthetic spectral energy distributions (SEDs) generated with the flexible stellar population synthesis (FSPS) code (Conroy et al 2009). We adopt the same models used in Mendel et al. (2014) to describe SDSS galaxies: a Chabrier (2003) initial mass function (IMF) with integration limits of 0.08 and 120 (imf_type = 1); an intrinsic dust attenuation using the extinction law of Calzetti et al. (2000, dust_type = 2); a smoothly declining star-formation history characterized by an e-folding timescale, . We apply a delayed-model (sfh=4) for the star-formation history. The contribution of nebular emission is computed using the photoionization code C (Ferland et al. 2013) and added to the spectrum. These choices result in 5 free parameters: the total stellar mass formed , the age age of the galaxy, the e-folding timescale , the intrinsic reddening ( − ), and the metallicity * . From these parameters, we derive the stellar mass, * = , where is the ratio of the surviving stellar mass to the formed mass, and the star-formation rate, SFR, computed as: where ( , ) is the lower incomplete gamma function. The massweighted stellar age, , is then derived as: We sampled the posterior probability density function of these parameters by using the affine-invariant ensemble MCMC sampler implemented in the package (Foreman-Mackey et al. 2013). We adopted uniform priors in log age , log , log , ( − ) over the same parameter range as Mendel et al. (2014, cf. their Table 2), and ran each fit with 128 walkers for 4096 steps, dropping the initial 100,000 steps as an initial burn-in phase and generating ≈400,000 posterior samples. The MCMC walkers were initialized near the maximum of the posterior, calculated through optimization of the likelihood function. Fits were performed with the code (Johnson et al. 2019), customized to use our chosen cosmology.

RESULTS
GRB 160624A and GRB 200522A are two short duration GRBs located at a similar distance ( = 0.483 and 0.554, respectively), which display very different properties in their observed emission. GRB 160624A is characterized by a bright, short-lived X-ray afterglow that is no longer detected after a few hundreds seconds. The faint afterglow and lack of any optical/nIR counterpart allow for stringent constraints on kilonova emission from the deep Gemini and HST observations (see §4.1). GRB 200522A displays instead a bright and long-lived counterpart. The red color of its optical/nIR emission (r-H ≈ 1.3 mag) represents tantalizing evidence for a kilonova, but its interpretation is complicated by the uncertain contribution of the standard afterglow. The burst location, close to the galaxy's center, and evidence for active star-formation suggest that extinction along the sightline could also contribute to the observed color (see §4.2).

Afterglow Properties
As shown in Figure 7, GRB 160624A displays a bright and rapidly fading X-ray afterglow. This feature is common among sGRBs and is often interpreted as long-lived activity of the central engine (e.g., Rowlinson et al. 2013). No evidence for a standard FS component is found by deep optical and X-ray follow-up observations. At early times (<2 hrs), this event is characterized by the deepest available optical limits for a sGRB, as demonstrated in Figure 8 (see also, e.g., Sakamoto, et al. 2013;Berger 2014;Fong et al. 2015). These observations would have detected nearly all the known sGRB optical afterglows, with the only exception of GRB 090515. Using afterglowpy, we explore the range of afterglow parameters allowed by the broadband upper limits. We fixed = 2.2, and left the other parameters ( 0 , 0 , , ) free to vary. Although low density (≈ 10 −3 cm −3 ) solutions are favored, ∼ 15% of the allowed models are consistent with 0 1 cm −3 . The faintness of the afterglow therefore does not necessarily imply a rarefied environment.
Our Chandra observation sets an upper limit to the X-ray luminosity 1.1 × 10 42 erg s −1 (0.3-10 keV) at 0 + 5.9 d (rest frame). This limit is compared in Figure 7 to the late-time X-ray excess detected for short GRBs 080503 (Perley et al. 2009) and 130603B (Fong et al. 2014), which is often attributed to a longlived and highly magnetized NS remnant (e.g., Gao et al. 2013;Metzger & Piro 2014). The interaction of the magnetar spin-down radiation with the merger ejecta could power a bright X-ray transient on timescales of a few days after the merger. The predicted peak luminosity is 10 43 − 10 45 erg s −1 with a decay following the temporal behavior of the spin-down emission, sd ∝ −2 . In order to be consistent with these models, our non-detection of X-rays favors a newborn NS with a large magnetic field 10 15 G for ejecta masses ej 10 −2 (Metzger & Piro 2014). Alternatively, the early steep decay of the X-ray afterglow may mark the NS collapse to a BH.

Kilonova Constraints
We constrain parameters of a potential kilonova associated with GRB 160624A by comparing optical and infrared upper limits to 900 on-axis kilonova simulations, introduced in §3.2.2. We only consider observations after 0 + 0.125 d (rest frame), as spectra are not computed prior to this time. Using spectra simulated at various times, we convert kilonova emission to observer frame lightcurves in Gemini/r-band, HST/ 606 , HST/ 125 , and HST/ 160 filters.
The four panels of Figure 9 show the range of simulated lightcurves in each filter. Colored regions indicate the range of lightcurves eliminated by upper limits, while gray regions indicate lightcurves consistent with observations. AT2017gfo lightcurves are included for comparison, and show that HST/ 125 and HST/ 160 observations may be sensitive to AT2017gfo-like kilonovae. The AT2017gfo photometry was compiled from Arcavi et al. Our HST/ 160 observations provide the most stringent constraints on the range of plausible lightcurves, disallowing 30% of on-axis lightcurves from the simulation grid. Individually, the Gemini/r-band upper limit at ∼ 1 d (observer frame) eliminated 12% of lightcurves, while the earlier HST/ 606 and HST/ 125 observations eliminate 8% and 23% of lightcurves, respectively. The HST/ 125 upper limit at ∼ 8 d post-merger (observer frame) and the HST/ 160 upper limit at ∼ 9 d post-merger (observer frame) only provide redundant constraints, eliminating lightcurves otherwise constrained by the four aforementioned upper limits. The HST/ 606 upper limit at ∼ 8 d post-merger (observer frame) places no constraint on the range of kilonova parameters. In total, 31% of simulated on-axis kilonovae are ruled out by the observational constraints on GRB 160624A. Figure 10 indicates the fraction of simulated lightcurves consistent with observations for all combinations of dynamical and wind ejecta mass. Constraints indicate that high ejecta masses are strongly disfavored, with 86% of simulations with 0.2 total ejecta mass (0.1 of dynamical ejecta + 0.1 wind ejecta) excluded by observational constraints. High wind ejecta masses (≥ 0.1 ) are disfavored with over 80% of models disallowed by the upper limits. Wind mass dictates lightcurve behavior at lower (optical) wavelengths, while dynamical ejecta mass dominates at higher (nIR) wavelengths. However, due to the cosmological distance of this GRB, our reddest filter, HST/ 160 , corresponds only to the rest frame y-band. As a result, our observations can only weakly constrain the dynamical ejecta with 53% of 0.1 models consistent with the data. Nearly all models with ejecta masses below 0.04 are consistent with the data.
For high ejecta masses, we are able to place strong constraints on the range of ejecta velocities. For example, lightcurves with wind ejecta masses of 0.1 and low velocities (≤ 0.15 ) are strongly disfavored. Similarly, models with wind ejecta 0.03 and low velocity (0.05 ) are predominately disfavored, while the majority of high mass models (with velocity 0.3 ) are consistent with the data. Velocity constraints are primarily due to the timing of our observations. For constant mass, higher velocities result in earlier and brighter peak emission (see, e.g., Thakur et al. 2020). As a result, many high velocity models have dimmed by the time of these observations and cannot be ruled out by the upper limits, while several low wind velocity models coincide with observations and are thus ruled out.  Figure 11. The best fit model spectrum (solid line) and photometry (squares) describing the host galaxy SED for GRB 160624A compared to the observed photometry (circles), which is corrected for Galactic extinction.

Environment
The best localization for GRB 160624A is its XRT position with an error radius of 1.7 (90% CL). This position intercepts a bright galaxy (Figure 2), which we identified as the GRB host galaxy. Using the XRT localization, the maximum projected physical offset from this host galaxy is 21.5 kpc (90% CL). Using the GALFIT package (Peng et al. 2002), we fit this galaxy with a Sersic profile of index =1 and derive an optical half-light radius ≈ 1.1 (6.8 kpc). Therefore, the maximum host-normalized offset is / 3.2. Following Bloom et al. (2002), we calculate the chance probability using the R-band number counts from Beckwith et al. (2006) for optical observations, and the H-band number counts from Metcalfe et al. (2006) for nIR observations. We derive =0.03 using the observed magnitude = 22.18 ± 0.02, and =0.02 using HST/ 160 = 20.566 ± 0.004 (see Table 2). We searched the field for other candidate host galaxies. We identify three bright SDSS galaxies at offsets of 21 , 22.4 , and 39 with ≈ 0.2, 0.5, and 0.9, respectively. There are a few dim extended objects, uncovered in the HST observations, at moderate offsets 4 with 0.5. Additionally, a faint source with 160 = 26.2 ± 0.3 mag is observed within the XRT position. Due to its faint nature, we cannot determine whether this is a star or a galaxy. In the latter case, the source's probability of chance coincidence is ≈ 0.7. Therefore, the association of GRB 160624A with the bright galaxy SDSS J220046.14+293839.3 remains the most likely.

Afterglow Properties
Before introducing our broadband modeling, we start with some basic considerations on the afterglow properties. In X-rays, the observed spectral index =0.45±0.18 suggests that < X < , otherwise 1.5. Here, is the injection frequency of electrons, and is the cooling frequency. A consistent spectral index (at ∼ 2 ) is derived from the optical and nIR observations, =1.1±0.3, suggesting the optical/nIR and X-ray data could lie on the same spectral segment. We therefore fit the broadband SED (X-ray/optical/nIR) at 0 + 3.5 d with an absorbed powerlaw model redden * tbabs(zdust(powerlaw)) 7 within XSPEC. We fix the Galactic column density =2.9×10 20 cm −2 , and = / ( − )=3.1 (Rieke & Lebofsky 1985). As there is no X-ray observation at this time, we re-scale the flux of the earlytime XRT spectrum to the predicted value at 3.5 d. This fit yields OX =0.77±0.05 and ( − )=0.12 +0.12 −0.08 mag ( 2 =71 for 78 dof), consistent with < IR < X < and = 2.54±0.10. However, for this value the predicted temporal slope, =3 /2≈1.16, is steeper than the slope, =0.84±0.04, measured in X-rays. In GRB afterglows, it is not uncommon to observe a shallow temporal decay, not accounted for by the simple FS model. In general, this behavior is observed in the early (<1000 s) lightcurve, whereas later observations tend to be consistent with standard closure relations. Indeed, in this case too, we find that, by excluding the early X-ray data, the temporal slope steepens to =1.04±0.08 in agreement with the FS predictions. Therefore, in our modeling we only consider the late (>1000 s) X-ray data, as earlier epochs could be affected by complex factors (e.g., energy injection, jet structure) not included in the basic FS model.
In the radio band, the afterglow is detected at a nearly constant level in the two early epochs (at 0 + 0.23 and 0 + 2.2 d), and then seen to fade (see Fong et al. 2020). In the simple FS model, the radio emission is expected to either rise as 1/3 , if below the spectral peak ( < ), or decay as −( −1)/2 when above the peak ( > ). The observed flattening can therefore be explained only if the synchrotron peak crosses the radio band between 0.1 and 2 days. Alternatively, the flat radio lightcurve could reveal the presence of a reverse shock (RS) component contributing at early times (Fong et al. 2020), as observed in other bright short GRBs (e.g., Soderberg et al. 2006;Becerra et al. 2019). Both scenarios are considered in our modeling.
We include an additional systematic uncertainty in the radio data to account for the scattering and scintillation of radio wavelengths due to interstellar scintillation (ISS). We adopt the 'NE2001' model (Cordes & Lazio 2002), which yields a scattering measure SM = 1.9×10 −4 kpc m −20/3 and a transition frequency 0 ≈ 8 GHz for the direction of GRB 200522A. Radio observations at < 0 can be effected by strong scattering when the angular extent of the GRB jet is GRB < 0 ≈ 1 µas. In our afterglow modeling, we find that this condition, GRB < 0 , is satisfied at 2.5 d. Therefore, we include a 30% systematic uncertainty, added in quadrature with the statistical uncertainty, to account for the effects of ISS.
Using an MCMC Bayesian fitting approach, outlined in §3.1, we explore four different models to describe the broadband after-7 The command redden corrects optical/nIR wavelengths for extinction within the Milky Way (Cardelli et al. 1989), whereas zdust accounts for intrinsic dust extinction within the host galaxy (Pei 1992). We selected the method that uses Milky Way extinction curves. Table 5. The fit results of our afterglow modeling for GRB 200522A. The median and 68% confidence interval of the marginalized posterior probability for each parameter corresponding to each mode determined from our MCMC fitting for the FS+Ext, Gauss+FS+Ext, FS+BB, and FS+BB+Ext models (see §4.2.1). The FS+Ext model has only one solution mode, whereas the Gauss+FS+Ext model has two modes which are degenerate. The three modes for the FS+BB and FS+BB+Ext models are: a solution with a late radio peak at 10 d (Mode I), a radio peak between 2-5 d that fits the radio detection at 2 d (Mode II), and an early radio peak at 0.3-0.8 d which describes both VLA radio detections at 0.2 and 2 d (Mode III). Row 15 shows the minimum 2 value associated to each mode. Rows 16-17 display the WAIC score of the expected log predictive density (elpd), and the WAIC score difference, ΔWAIC elpd , compared to the FS+Ext model. The WAIC score is computed for the overall model fit (see Table A1 This is the only mode of the solution for this model. This mode is comprised of two solutions, but due to the high degree of degeneracy it is not possible to further sort these solutions. This mode of the solution appears only when the first radio detection at 0 + 0.2 d (Fong et al. 2020) is included within the fit.
glow of GRB 200522A from radio to X-rays. By assuming a top-hat jet, we consider (i) a forward shock with intrinsic extinction from the host galaxy (hereafter denoted by FS+Ext), (ii) a forward shock with an additional blackbody component (FS+BB, see §3.2.1), and (iii) a forward shock and simple blackbody with the addition of intrinsic extinction (FS+BB+Ext). For a structured jet, we consider (iv) a Gaussian profile, standard forward shock emission, and intrinsic extinction (Gauss+FS+Ext). The results are tabulated in Tables 5  and A1. We discuss the results of these fits for the models with only FS emission in §4.2.2 and models with an additional BB component in §4.2.3. A comparison of the WAIC score difference between these models is presented in §4.2.4. Lastly, in §4.2.5, we explore the consistency of the predicted flux from our BB models with detailed kilonova simulations ( §3.2.2).

Forward Shock Afterglow Models
Here, we discuss the inferred parameters for the model with standard FS emission, including intrinsic extinction from the host galaxy, for two jet configurations (see §3.1): a top-hat jet viewed on-axis ( ≈ 0) and a Gaussian angular profile viewed at an arbitrary angle . A comparison of these models to the broadband dataset of GRB 200522A is shown in Figure 12. In Figures A1 and A2 we present the marginalized posterior distributions for each parameter from our MCMC fit for the tophat and Gaussian jet, respectively. From this we identify that the Gaussian jet model has a bi-modal solution which we refer to as Mode I and Mode II; see Table 5. Mode I yields an extremely large beaming-corrected energy ≈ 2 × 10 52 erg and a very low circumburst density, 0 ≈ 3×10 −6 cm −3 , likely inconsistent with the observed offset . Mode II instead yields more typical values ≈ 2 × 10 49 erg, 0 ≈ 2 × 10 −3 cm −3 , consistent with those inferred for the top-hat jet model. Although both modes are statistically equivalent, we consider Mode I a less realistic description of the afterglow parameters. Figure 12 (top panel) shows that the standard FS model (FS+Ext) can describe the broadband dataset, except for the early radio and X-ray detections. An achromatic jet-break (Rhoads 1999;Sari et al. 1999) at ≈5 d is required by the data, which constrains the jet opening angle to ≈ 0.16 rad (9 • ). The Gaussian angular profile (bottom panel) leads to a slightly better description of the X-ray light curve. Due to its shallower temporal decline, it can more easily reproduce (within the 2 uncertainty) the first X-ray detection at <1000 s without requiring additional energy injection. However, unlike the top-hat jet model, the Gaussian angular profile does not provide a tight constraint on the jet's opening angle, . Although the ratio / is better constrained to / 2.3.
Both these jet models underestimate the radio detections at 0 + 0.23 d and at 0 + 2.2 d. This could be attributed to synchrotron self Compton (SSC) and Klein-Nishina effects (e.g., Jacovich et al. 2020), not included in our code, and/or to a bright RS component.

Afterglow Models including a Blackbody Component
In this section, we present the multi-modal solutions for our afterglow models that include a simple blackbody component, with (FS+BB+Ext) and without extinction (FS+BB). The fit to the broadband dataset for each model is shown in Figure 13, and the parameter values are presented in Table A1. Marginalized posterior distributions for each parameter are displayed in Figures A3 and A4 for the FS+BB and FS+BB+Ext models, respectively. Each model exhibits three modes, referred to as Mode I, Mode II, and Mode III, which are presented individually in Table 5. Mode I is characterized by a late ( >10 d) radio peak and no jet-break, Mode II shows an earlier peak (2 d< <5 d) and requires a jet-break, and Mode III 8 also 8 We note that Mode III appears with high significance only when the first radio detection is included in the fit. Figure 13. Same as in Figure 12, but for the forward shock model with a blackbody component (Top: model FS+BB) and with intrinsic extinction (Bottom: model FS+BB+Ext). The shaded regions display only the afterglow contribution to the flux. Multiple solutions are consistent with the data: a late radio peak at 10 d (Mode I), a radio peak between 2-5 d (Mode II), and an early radio peak at 0.3-0.8 d (Mode III; dotted line). The excess emission (HST/F125W) is compared to a simulated kilonova lightcurve (dashed maroon line) with properties: ej,d =0.001 , ej,d =0.3 , ej,w =0.1 , and ej,w =0.15 . finds a jet-break and describes the first and second radio detections without a RS component due to an early radio peak (0.3-0.5 d). Although the MCMC algorithm cannot distinguish one of these modes as providing a better description of the data, we disfavor Mode I based on the extremely low circumburst density ( ≈ 10 −5 cm −3 ) and the phenomenology of sGRB afterglows. Such late radio peaks have not been observed in sGRBs, and Mode I is likely an artifact of the poorly constrained late-time radio dataset. We therefore favor either Mode II or Mode III as a more likely description of the jet parameters. Both of them constrain the jet-opening angle to ≈ 0.10 rad (6 • ).
When including extinction, the only difference in the fit is in the temperature and radius of the blackbody component: the model without dust extinction requires a cooler thermal component (T≈4,000 K) to match the optical data. A higher temperature (T 6,500 K), as predicted, e.g., in the magnetar-boosted model (Fong et al. 2020), tends to overpredict the optical flux, unless dust extinction contributes to attenuate the observed emission.

Model Comparison
We perform a comparison of the models applied in this work using their WAIC scores, described in §3.1 (see also Troja et al. 2020;Cunningham et al. 2020). We note that the WAIC analysis is not applicable to individual modes within the models, and that we only compare the overall model fits presented in Table A1. For the four models considered, the difference between the WAIC scores is not significant enough to statistically favor any of them. In particular, the addition of a blackbody to the FS fit is not required by the data. We find that the WAIC score of the FS+BB+Ext model only marginally (at the 1.2 level) improves over the FS+Ext model.
The most significant difference in WAIC score (ΔWAIC elpd = −62 ± 40) is between the FS+Ext and Gauss+FS+Ext models. A larger WAIC score implies a better description of the data (see §3.1), and, therefore, the FS+Ext model is marginally preferred at the ≈ 1.4 level.
Our findings do not depend on the details of the data analysis, which yield slightly different magnitudes for the nIR counterpart (see §2.2.3). In particular, we verified that, by using the values presented by Fong et al. (2020), our results are unchanged. There is no significant difference between the posterior distributions of the fit parameters, and the WAIC score comparison continues to not favor any particular model fit.
We utilize our MCMC algorithm to determine a posterior distribution on the temperature and radius of the blackbody (Table  5), from which we compute a distribution for the emitted flux from the blackbody in each filter. The flux posterior distributions are then compared to the LANL suite of kilonova simulations (see §3.2.2) in order to estimate the ejecta masses and velocities required to reproduce the observations. Figure 14 presents five kilonova lightcurves (solid lines) consistent with the black body flux estimated in the FS+BB+Ext model. They lie within the inner 50% credible interval of both HST/ 125 and HST/ 160 posteriors. As the FS contribution likely dominates at optical wavelengths, we do not require consistent lightcurves to reside in the 50% credible interval of the Gemini/r-band constraint.
These results indicate that any thermal component, as constrained from our broadband modeling, agrees with a radioactively powered kilonova emission. The five consistent lightcurves span a wide range of dynamical ejecta masses (0.001 ej,d 0.1 ), but share many other properties, including a wind ejecta mass of ej,w =0.1 , wind ejecta velocity of ej,w =0.15 , dynamical ejecta velocity of ej,d =0.3 , and spherical wind morphology. These observations provide stronger constraints on the wind ejecta mass as the 125 and 160 filters probe rest frame optical wavelengths. We emphasize that these 5 lightcurves represent only a small subset of kilonova parameters capable of reproducing the flux posteriors.
This result differs from the findings of Fong et al. (2020), who argue for an additional power source (e.g., an enhanced radioactive heating rate or a magnetar) to boost the kilonova luminosity to values higher than AT2017gfo. This difference arises despite comparable predictions for the nIR thermal emission, with their predicted luminosities ( F125W ≈ (9.5 − 12.3) × 10 41 erg s −1 and F160W ≈ (8.9 − 11.4) × 10 41 erg s −1 ) fully encapsulated within our luminosity posterior distribution. Our multi-dimensional kilonova models can reproduce this range of values by incorporating the same physics adopted to model AT2017gfo (Troja et al. 2017;Evans et al. 2017;Tanvir et al. 2017;Wollaeger et al. 2018), including a thermalizable heating rate of ≈ 10 10 erg g −1 s −1 at = 1 day. In addition, they capture the multi-component character of kilonova ejecta, leading to an enhancement of the emission via photon reprocessing (Kawaguchi et al. 2020). This is illustrated in Figure 14, which compares our models (solid lines), produced by a spherical wind component girdled with a toroidal lanthanide-rich ejecta, with the emission produced by a single-component spherical morphology (dashed line), with properties ( ej =0.1 , ej,w =0.15 , =0.27) similar to the models used in Fong et al. (2020). The addition of the toroidal belt produces a ≈ 1 mag enhancement of the HST/F160W flux compared to the single-component model. This is attributed to the reprocessing of photons emitted from the spherical wind by the high-opacity ejecta. Photons absorbed by the toroidal component preferentially diffuse towards the polar regions, and are re-emitted at redder wavelengths. The flux enhancement at optical wavelengths is instead negligible. This effect is more prominent in events viewed along the polar axis, such as GRB 200522A.

Environment
GRB 200522A is located within a bright galaxy, SDSS J002243.71-001657.5, at = 0.554. The projected offset from the galaxy's nucleus is = 0.16±0.04 (1.07±0.27 kpc); in the bottom 15% of the sGRB offset distribution (Fong & Berger 2013). Using GALFIT, we derive a projected half-light radius = 0.60 ± 0.02 , corresponding to 4.0 ±0.1 kpc, and a normalized offset of / ≈ 0.27 ± 0.07. The chance alignment between the GRB and the galaxy is small: we determine =0.002 using the r-band magnitude, and =0.003 using the 160 filter. Both the GRB offset and the galaxy's angular size contribute to determine this value (Bloom et al. 2002). As shown in Figure 5, there are two nearby red galaxies seen at projected angular offsets of 2.5 and 3.9 with 160 ≈ 23.5 and 20.8 AB mag, respectively. The probability of chance coincidence is =0.11 and 0.05, respectively. Therefore, we consider SDSS J002243.71-001657.5 to be the likely host galaxy of GRB 200522A.
We determine the physical properties of the galaxy using the methods described in §3.3. The resulting best fit model is shown in Figure 15. We find an intrinsic extinction ( − ) = 0.02 ± 0.01 mag, a metallicity * / = 1.35 +0.16 −0.25 , an -folding time = 0.13 +0.07 −0.03 Gyr, a stellar mass log( * / ) = 9.44 ± 0.02, a young stellar population = 0.35 +0.08 −0.04 Gyr, and a star-formation rate SFR = 2.00 ± 0.12 yr −1 . These values are broadly consistent with those inferred by Fong et al. (2020), although they derive a slightly higher stellar mass. The stellar mass and age derived here are on the low end of the distributions for sGRB host galaxies, but are not unique to the population (Leibler & Berger 2010).
Moreover, we perform standard emission line diagnostics on the spectrum described in §2.2.4. The line flux ratio H /H = 0.25 ± 0.11 does not provide a strong constraint on the intrinsic extinction ( − ). Fong et al. (2020) infer negligible extinction, however, based on their reported line ratio H /H = 2.9 ± 0.9, we derive ( − ) 1.3 mag (3 ) which does not rule out the presence of moderate extinction along the GRB line of sight. The dominant emission lines imply significant on-going star formation. We derive a SFR(H ) ≈ 3.5 ± 0.9 yr −1 under the assumption H /H = 2.86, and SFR([OII])= 4.8 ± 1.3 yr −1 using the [OII] line luminosity (Kennicutt 1998).

CONCLUSIONS
We have presented an analysis of the multi-wavelength datasets for two sGRBs at ∼ 0.5: GRB 160624A at = 0.483 and GRB 200522A at = 0.554. These two events demonstrate the wide range of diversity displayed by sGRBs in terms of both their emission and environment. We utilize the broadband datasets for these two events to constrain the presence of kilonova emission arising after the compact binary merger.
Gemini and HST observations of GRB 160624A place some of the deepest limits on both optical and nIR emission from a short GRB. For this event, we find that the bright short-lived X-ray counterpart is likely related to long-lasting central engine activity, whereas emission from the forward shock is not detected. Thanks to the negligible afterglow contribution, we can robustly constrain emission from a possible kilonova. By comparing our limits to a large suite of detailed simulations, we derive a total ejecta mass 0.1 , favoring wind ejecta masses ej,w 0.03 . Any kilonova brighter than AT2017gfo (∼30% of the simulated sample) can be excluded by our observations. Very different is the phenomenology of GRB 200522A, which has a long-lived X-ray emission, largely consistent with standard forward shock models, and a bright nIR counterpart. Its unusual red color ( − ≈ 1.3) at a rest frame time ∼ 2.3 d is suggestive of a kilonova, although the inferred luminosity of F125W ≈ (7 − 19) × 10 41 erg s −1 is above the typical range observed in other candidate kilonovae Yang et al. 2015;Ascenzi et al. 2019). The identification of a kilonova is complicated by the bright, long-lived afterglow contribution, as well as by the limited color information available for this event.
Our thorough modeling of the multi-wavelength dataset finds that a standard forward shock model represents a good description of the X-ray, optical and nIR dataset provided that a modest amount of extinction, ( − ) ∼0.1-0.2 mag, is present along the GRB sightline. This value is consistent with the constraints from optical spectroscopy, ( − ) 1.3 mag. The location of GRB 200522A within its host (∼ 1.1 kpc from the center) and evidence for ongoing star-formation (∼2-4 yr −1 ) also support that dust effects might not be completely negligible, as instead assumed by Fong et al. 2020. We constrain the afterglow parameters to: ≈ 6 × 10 48 erg, 0 ≈ 2 × 10 −3 cm −3 , ≈ 2.3, and identify the presence of a jetbreak at ≈ 5 d. From this we derive an opening angle of ≈ 0.16 rad (9 • ), providing additional evidence that short GRB outflows are collimated into jets with a narrow core (e.g. Burrows et al. 2006;Troja et al. 2016). Since this GRB was likely observed close to its jet-axis ( / 2; for comparison, GW170817 had / ≈ 6; Ryan et al. 2020;Beniamini et al. 2020), no significant constraint can be placed on the jet structure and a simple top-hat jet profile seen on-axis ( ≈ 0) already provides an adequate description. Our best fit model cannot reproduce the early radio detection with a simple forward shock emission, suggesting the presence of an early reverse shock component (see also Jacovich et al. (2020) for a discussion of SSC effects).
A different interpretation is of course that the red color of the optical/nIR emission marks the onset of a luminous kilonova. Statistically, this scenario (models FS+BB and FS+BB+Ext in Table A1) is not preferred over a simple afterglow model (FS+Ext in Table A1) and therefore a kilonova component is not required by the data. Nevertheless, we explore the range of kilonova models consistent with our dataset. A radioactively-powered kilonova with wind ejecta mass ej,w = 0.1 , wind velocity ej,w = 0.15 , and dynamical ejecta velocity ej,d = 0.3 is capable of reproducing the observed nIR emission when considering a toroidal morphology for the lanthanide-rich ejecta ). This geometry can naturally produce brighter transients than spherically symmetric counterparts of equal mass, expansion velocity, and radioactive heating. The range of dynamical ejecta masses is however not well constrained, as the rest frame wavelengths probed by the observations lie in the optical band.
The ejecta mass implied by the nIR luminosity is slightly larger than the values derived for AT2017gfo ( ej,w ≈0.02-0.07 ), and pushes the boundaries of a standard NS merger model. Numerical simulations of accretion discs indicate that, following a NS merger, ≈10-40% of the disc can become unbound and form a massive outflow along the polar axis (see, e.g., Perego et al. 2017). Since these mergers can form discs up to ≈0.3 , a wind ejecta mass ej,w = 0.1 is still within the range of possible outcomes (e.g., Fujibayashi et al. 2020). Moreover, our derived ejecta mass is limited by the resolution of the LANL simulations which do not fully probe the range 0.03-0.1 . It is conceivable that a finer sampled grid of simulated light curves could find additional solutions in this mass range. Based on these considerations, we find no compelling evidence for a magnetar-powered kilonova discussed by Fong et al. 2020.
Our study of GRB 160624A and GRB 200522A demonstrates that deep HST observations can probe an interesting range of kilonova behaviors out to ∼ 0.5. However, whereas sensitive HST imaging can detect the bright kilonova emission, it is not sufficient to unambiguously identify a kilonova and disentangle it from the standard afterglow. Observations of GRB 200522A, based on a single multi-color epoch, can not break the degeneracy between the different models and, overall, the presence of a kilonova can not be confidently established. As shown by previous cases, such as GRB 130603B and GRB 160821B, multi-epoch multi-color observations are essential for the identification of a kilonova bump. Moreover, at these distances, we found that the component of lanthaniderich ejecta is only weakly constrained by the HST observations, with an allowed range of masses that spans two orders of magnitude (0.001-0.1 ). Future IR observations with JWST will be pivotal to constrain the properties of lanthanide-rich outflows from compact binary mergers. acknowledges financial support from the IDEAS Fellowship, a research traineeship program funded by the National Science Foundation under grant DGE-1450006. Partial support to this work was provided by the European Union Horizon 2020 Programme under the AHEAD2020 project (grant agreement number 871158). G.R. acknowledges support from the University of Maryland through the Joint Space Science Institute Prize Postdoctoral Fellowship. The work of E.A.C., C.L.F., R.T.W., C.J.F and O.K. was supported by the US Department of Energy through the Los Alamos National Laboratory, which is operated by Triad National Security, LLC, for the National Nuclear Security Administration of US Department of Energy (Contract No. 89233218CNA000001).
This work made use of data supplied by the UK Swift Science Data Centre at the University of Leicester. The scientific results reported in this article are based on observations made by the Chandra X-ray Observatory. This research has made use of software provided by the Chandra X-ray Center (CXC) in the application package CIAO. Based on observations obtained at the international Gemini Observatory, a program of NSF's OIR Lab, which is managed by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation on behalf of the Gemini Observatory partnership: the National Science Foundation (United States), National Research Council (Canada), Agencia Nacional de Investigación y Desarrollo (Chile), Ministerio de Ciencia, Tecnología e Innovación (Argentina), Ministério da Ciência, Tecnologia, Inovações e Comunicações (Brazil), and Korea Astronomy and Space Science Institute (Republic of Korea). The HST data used in this work was obtained from the Mikulski Archive for Space Telescopes (MAST). STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. Support for MAST for non-HST data is provided by the NASA Office of Space Science via grant NNX09AF08G and by other grants and contracts. These results also made use of Lowell Observatory's Lowell Discovery Telescope (LDT), formerly the Discovery Channel Telescope. Lowell operates the LDT in partnership with Boston University, Northern Arizona University, the University of Maryland, and the University of Toledo. Partial support of the LDT was provided by Discovery Communications. LMI was built by Lowell Observatory using funds from the National Science Foundation (AST-1005313). We additionally made use of Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration et al. 2018). The afterglow modeling was performed in part on the George Washington University (GWU) Pegasus computer cluster and on the YORP cluster administered by the Center for Theory and Computation, which is part of the Department of Astronomy at the University of Maryland (UMD).

DATA AVAILABILITY
The data underlying this article will be shared on reasonable request to the corresponding author. Table A1. Fit results of our afterglow modeling for GRB 200522A. We consider four model scenarios to describe the broadband transient: (i) a tophat jet model with standard forward shock emission and intrinsic extinction from the host galaxy (FS+Ext), (ii) forward shock emission with the addition of a simple blackbody component (FS+BB), (iii) forward shock emission and a blackbody component with intrinsic extinction (FS+BB+Ext), and (iv) a Gaussian jet model with forward shock emission and intrinsic extinction. This table represents the results of the posterior distribution for each fit, and in Table 5 we demonstrate the posterior distributions corresponding to the multiple modes allowed by each fit (see §4.2.2 and 4.2.3). The median and 68% confidence interval of the marginalized posterior distribution for each parameter is shown in rows 3-14. Row 15 shows the minimum 2 value for each model. Rows 16-17 display the WAIC score of the expected log predictive density (elpd), and the WAIC score difference, ΔWAIC elpd , compared to the FS+Ext model.  The contours (off-diagonal) and dashed lines (diagonal) correspond to the 16%, 50%, and 84% quantiles of the marginalized posterior probability for each parameter.  Figure A1 for the Gauss+FS+Ext model corresponding to a Gaussian structured jet profile. The two modes (see Table 5) of the fit are shown in red and blue for Mode I and Mode II, respectively. Each mode is preferred with the same weight by the MCMC fit, i.e., ∼ 50% of walkers prefer each mode.  Figure A1 for the FS+BB+Ext model.