Hubble Space Telescope Observations of GW170817: Complete Light Curves and the Properties of the Galaxy Merger of NGC 4993

We present the complete set of Hubble Space Telescope imaging of the binary neutron star merger GW170817 and its optical counterpart AT 2017gfo. Including deep template imaging in F814W, F110W, F140W, and F160W at 3.4 yr post-merger, we reanalyze the full light curve of AT 2017gfo across 12 bands from 5 to 1273 rest-frame days after merger. We obtain four new detections of the short γ-ray burst 170817A afterglow from 109 to 170 rest-frame days post-merger. These detections are consistent with the previously observed β = −0.6 spectral index in the afterglow light curve with no evidence for spectral evolution. We also analyze our limits in the context of kilonova afterglow or IR dust echo emission but find that our limits are not constraining for these models. We use the new data to construct deep optical and IR stacks, reaching limits of M = −6.3 to −4.6 mag, to analyze the local environment around AT 2017gfo and low surface brightness features in its host galaxy NGC 4993. We rule out the presence of any globular cluster at the position of AT 2017gfo to 2.3 × 104 L ⊙, including those with the reddest V − H colors. Finally, we analyze the substructure of NGC 4993 in deep residual imaging and find shell features that extend up to 71.″8 (14.2 kpc) from NGC 4993. The shells have a cumulative stellar mass of 6.3 × 108 M ⊙, roughly 2% of NGC 4993, and mass-weighted ages of >3 Gyr. We conclude that it was unlikely that the GW170817 progenitor system formed in the galaxy merger.


Introduction
The discovery and localization of the binary neutron star (NS) merger GW170817 (Abbott et al. 2017a(Abbott et al. , 2017b) from a gravitational wave (GW) signal and its optical counterpart enabled the first detailed study of this rare phenomenon. The electromagnetic counterpart to this event was initially identified from a short γ-ray burst (GRB; Goldstein et al. 2017) called GRB 170817A, confirming the hypothesis that NS mergers are sources of these high-energy astrophysical phenomena and launch relativistic jets (Lattimer & Schramm 1976;Li & Paczyński 1998;Metzger et al. 2010). GRB 170817A was followed 11 hr later by the discovery of a kilonova, or an optical transient powered by the radioactive decay of r-process elements synthesized in the merger's neutron-rich ejecta (called AT 2017gfo; Chornock et al. 2017;Coulter et al. 2017;Cowperthwaite et al. 2017;Drout et al. 2017;Kilpatrick et al. 2017;Tanvir et al. 2017;Troja et al. 2017). When the counterpart became visible 111 days after merger and after Sun constraint, its optical and IR emission appeared dominated by a broadband, synchrotron-powered afterglow from a relativistic and structured jet Lazzati et al. 2018;Lyman et al. 2018;Margutti et al. 2018;Mooley et al. 2018;Troja et al. 2018;Fong et al. 2019), also observable at early times in both the radio and X-ray bands (Alexander et al. 2017;Haggard et al. 2017;Hallinan et al. 2017;Margutti et al. 2017;Troja et al. 2017).
Hubble Space Telescope (HST) imaging and spectroscopy of GW170817 contributed significantly to analysis of the kilonova (Cowperthwaite et al. 2017;Kasliwal et al. 2017;Tanvir et al. 2017;Troja et al. 2017), GRB afterglow (Lyman et al. 2018;Fong et al. 2019;Lamb et al. 2019), and its S0 host galaxy NGC 4993 (Blanchard et al. 2017;Coulter et al. 2017;Kilpatrick et al. 2017;Levan et al. 2017;Palmese et al. 2017;Pan et al. 2017). Despite being the closest known NS merger to date at ≈40 Mpc, optical and near-IR emission from GW170817 faded below the detection thresholds of the largest ground-based telescopes owing to a combination of poor observability, intrinsic faintness (>26 mag), and contaminating light from its bright host galaxy. Thus, all optical detections and the most constraining upper limits on AT 2017gfo at >100 days have been enabled by HST (Figure 1). This large, homogeneous set of high-resolution imaging, spanning serendipitous archival imaging from months before the merger to exhaustive follow-up campaigns years after detection of the GW170817, uniquely probes its optical light curve, local environment, and faint features in the host galaxy.
HST optical light curves enabled constraints on the bulk energetics, ejecta velocity, and opacity of kilonova ejecta, which was observed in distinct "blue" and "red" components at early and late times, respectively (Cowperthwaite et al. 2017;Tanvir et al. 2017;Troja et al. 2017). The latter component indicates that kilonovae produce a significant mass of lanthanides in their ejecta and may account for the bulk of heavy-element production in the "third peak" of r-process production Metzger et al. 2017).
Combined with contemporaneous data obtained in the radio and X-ray bands, the HST data of AT 2017gfo beyond ≈100 days primarily probe the jetted, relativistic outflow from the NS merger. These light curves exhibited a constant spectral index for the first ≈900 days with a peak at 160 days followed by a relatively rapid decline (Margutti et al. 2017;Lyman et al. 2018;Lamb et al. 2019;Fong et al. 2019). Modeling of these data is best described by a structured jet with a relatively narrow, collimated core (3°-5°) and a wider-angle outflow moving at slower velocity (Alexander et al. 2017;Margutti et al. 2017;Wu & MacFadyen 2018;Hajela et al. 2019;Wu & MacFadyen 2019;Margutti & Chornock 2020). However, variations in the spectral index at 1234 days post-merger as seen in recent X-ray detections suggest that the structured jet is evolving or some new emission component, such as a relativistic shock from the slower-moving kilonova ejecta, is beginning to dominate the afterglow light curve (Hajela et al. 2020;Balasubramanian et al. 2021;Hajela et al. 2021;Troja et al. 2022). With the exception of HST/F606W data (Fong et al. 2019), optical and near-IR measurements of the faint afterglow luminosities were performed without deep template images, introducing uncertain contributions from the host galaxy.
The joint set of HST observations also enabled the highestresolution analysis of the local environment around AT 2017gfo, globular clusters (GCs) proximate to its merger site, and the global properties of NGC 4993. Some of the most intriguing features of this galaxy are the shells of gas and stars extending ≈13 kpc from the center of the galaxy (Blanchard et al. 2017;Palmese et al. 2017;Ebrová et al. 2020), whereas the projected separation of AT 2017gfo is only 2.2 kpc. These shells likely indicate a merger between NGC 4993 and a less massive galaxy within the past few hundred million years (Palmese et al. 2017;Ebrová et al. 2020). This inference is based on the relatively limited amount of HST imaging that Figure 1. Three-color HST WFC3/IR imaging of a 120.9 × 116.8 arcsec 2 region covering site of GW170817 in F160W (red), F140W (green), and F110W (blue), composed of all data taken to date in these bands. The black box represents the region highlighted by the subpanels centered on AT2017gfo, a 9.1 × 7.9 arcsec 2 region in F814W, F110W, F140W, and F160W from the latest epochs in each band (obtained on 2021 February 22, February 7, January 4, and January 6, respectively, as described in Table 1).
was available within 2 yr from the event. Detailed analysis of all HST imaging obtained since then could reveal additional, fainter shells interior or exterior to the primary structure. There is currently over 140 ks of wide-band HST imaging from the ultraviolet (UV) to near-IR (Table 1), and a more detailed analysis of the shell structures' overall morphology, luminosity, and colors can constrain the time since they were accreted, their mass, and the age of the associated stellar population.
Here we present all HST observations of AT 2017gfo and its host galaxy obtained to date, including new data: 14 orbits of Wide Field Camera 3 (WFC3) F814W, F110W, F140W, and F160W imaging obtained from 2021 January to February. Throughout this paper, we assume a line-of-sight Milky Way reddening to AT 2017gfo and NGC 4993 of E(B − V ) = 0.109 mag (Schlafly & Finkbeiner 2011). All magnitudes given throughout this paper are on the AB magnitude system (Oke & Gunn 1983). All dates and times are given in Coordinated Universal Time (UTC).

Observations
We retrieved all observations of GW170817/AT2017gfo from the MAST archive across all filters, spanning from around 5 to 1285 days after merger (observer frame; note that in Table 1 and for figures showing the light curve we provide the epoch in restframe days), as well as the pre-merger HST/ACS image from 2017 April 28. These observations used 12 different filters from the UV to the IR bands, and they include data from the Advanced Camera for Surveys (ACS) Wide Field Camera (WFC) and the Wide Field Camera 3 (WFC3) UVIS and IR detectors. In particular, we present new WFC3 observations from 2021 January 4 to February 22 under Program 15886 (PI: Fong) in the F814W, F110W, F140W, and F160W filters (see Table 1). The primary purpose of these late-time observations is to serve as "template images" for the image subtraction procedure (discussed in Section 2.1) in these bands. The new templates comprise three to four orbits in each band, for a cumulative on-source exposure time of 7.8-10.4 ks in each visit.
Starting with the flc (for WFC3/UVIS and ACS/WFC) and flt (for WFC3/IR) files, we reduced all observations using the hst123 v1.0.0 analysis and reduction code 10 as described in Kilpatrick et al. (2021a) and Kilpatrick (2021). We first combine images from every unique visit and band as listed in Table 1. We used astrodrizzle to optimally stack and regrid each WFC3/UVIS observation to a pixel scale of 0 05 pixel −1 and each WFC3/IR observation to 0 064 pixel −1 with driz_sep_pixfrac=0.8 and final_pixfrac=0.8. Frame-to-frame alignment and alignment of the final drizzled images are performed in drizzlepac.tweakreg, with 100-700 sources on average per frame depending on depth. The final rms alignment precision in the drizzled frames is 0 002-0 017. We also used drizzlepac.photeq to ensure a uniform photometric zero-point across both WFC3/ UVIS chips before image combination. In this way, photometric precision is preserved in the combined images, which are known to exhibit 0.4% rms photometric variation compared with the original flux-calibrated HST images (see discussion in Fruchter et al. 1997;Fruchter & Hook 2002;McMaster et al. 2008).
We show an IR three-color image (F110W, F140W, F160W) composed of the template images and single-filter zoomed-in versions centered on the location of GW170817/AT2017gfo in the F814W, F110W, F140W, and F160W filters in Figure 1. The lack of any apparent source at the transient location means that they can be adequately used for templates against which we can subtract any earlier imaging (Section 2.1).
In addition to creating stacks for each unique visit, we create "deep stacks" by combining every available image observed in the F606W, F814W, F110W, F140W, and F160W filters (bottom five rows of Table 1). We use these "deep stacks" of NGC 4993 for our analysis of the galaxy structure and environment (Section 5).

Photometry and Image Subtractions for AT 2017gfo
We performed photometry in every flc/flt frame (for ACS/ WFC, WFC3/UVIS, and WFC3/IR where appropriate) using dolphot v2.0 (Dolphin 2016). Our reductions followed standard recommendations for each imager as described in dolphot (Dolphin 2016) and hst123 ). We use this photometry for imaging without late-time templates or images where AT 2017gfo is detected at >20σ, which comprise every detection between 2017 August 22 and 29. These data are reported in Table 1.
For the remaining data in which AT 2017gfo is detected at <20σ, comprising observations obtained from 2017 December 6 through 2018 August 14, we report photometry derived from subtracted imaging, which comes primarily from HST/WFC3. Our F606W photometry, including all imaging of AT 2017gfo observed by ACS during this time period, is taken from Fong et al. (2019), which follows a similar procedure to image subtraction described below.
Subtracting early-and late-time imaging stacks requires an understanding of the point-spread function (PSF) shape in both epochs. WFC3 has a relatively stable PSF with little change in shape observed over 11 yr of operation. The primary PSF variations are "breathing" modes due to thermal expansion over the orbital period of the spacecraft. The effect of these changes in the FWHM of the PSF is largest at bluer wavelengths, with at most 0.3% variation redward of 8000 Å (Kimble et al. 2008). Otherwise, the WFC3 PSF in this wavelength range is well approximated by a Gaussian profile (before pixelation). Therefore, the difference between two WFC3 frames observed in the same filter and instrumental configuration is dominated by the current position angle (PA) of the spacecraft on the sky, which in turn affects the PA of any nonaxially symmetric components of the PSF, such as the location of diffraction spikes (as in Figure 1, where the PA was ≈110°). Thus, in order to match PSF shape with imaging observed at different epochs, we use a relatively simple, Gaussian convolution kernel.
For all F814W, F110W, F140W, and F160W imaging, which comprises our late-time imaging where AT 2017gfo is detected at <20σ, we follow a similar procedure to Fong et al. (2019) using hotpants v5.1.11  to subtract our template images from the combined imaging in each visit. The specific parameters used in each subtraction were varied in order to reduce the rms residuals from stars observed close to AT 2017gfo, but in general our default parameters are bgo=0.1, ko=0.05, and nsx=nsy=5. In all cases, we convolve and normalize the input image to the template image (i.e., c=t, n=t). These parameters are similar across both WFC3/UVIS and WFC3/IR imaging. We show example subtractions from all four bands in Figure 2.
After image subtraction, we derive photometry for AT 2017gfo in all four bands by empirically reconstructing the instrumental PSF from isolated stars in the original science image using the python-based tool photutils v0.7.2 (Bradley et al. 2020). We then perform forced PSF photometry at the position of AT 2017gfo, which is derived by aligning each frame to early-time imaging of the kilonova. Photometric uncertainties were calculated from the χ 2 of the profile fit following methods described in Stetson (1987). All detections of the afterglow from subtracted imaging are reported in Table 1. 10 https://github.com/charliekilpatrick/hst123

Upper Limits on Emission near AT 2017gfo
For several epochs, we do not detect any significant (>3σ) emission when performing forced photometry at the location of AT 2017gfo. In these cases, we place upper limits on the presence of any optical or near-IR counterpart with the FakeStars procedure in dolphot. Following procedures described in Kilpatrick et al. (2017), we injected 1000 sources into magnitude bins of 0.1 mag from 20 to 30 mag (i.e., 100,000 sources in total) around the location of AT 2017gfo in that frame after our full image subtraction procedure. We varied the position of each source randomly by drawing a Gaussian random variable centered at (x,y) = (0,0) and with a Gaussian FWHM PSF theoretical S N = for that source. We then determined the magnitude threshold at which >99.7% of sources were recovered at >3σ, which we consider to be the limiting magnitude for that visit as reported in Table 1. We acknowledge that since the FakeStars images have not been through the same subtraction procedure as the F814W, F110W, F140W, and F160W frames used above to identify AT 2017gfo at >110 days, our limits may be slightly deeper than those achievable by the full subtraction and photometry procedure.
We repeated this procedure independently for both the template images and deep stacks. In the latter case, instead of injecting sources at the location of AT 2017gfo, we instead injected them randomly within a 2″ radius of that position. In this way, we avoid biasing our magnitude limit with residual flux from AT 2017gfo in each of the frames we stacked. Thus, we consider the limiting magnitude for the deep stacks to be the threshold at which we can detect point-like sources near AT 2017gfo, rather than AT 2017gfo itself. All of these limits are reported in Table 1.

The Complete HST Light Curve of the GW170817 Counterpart
We confirm previous HST detections of AT 2017gfo observed after 2017 December 6 (e.g., in Lyman et al. 2018;Margutti et al. 2018;Troja et al. 2018;Lamb et al. 2019;Ryan et al. 2020) in our difference imaging. In addition, the use of deep template imaging in our subtractions uncovers four new detections ( Figure 2 and Table 1) in F814W (2017 December 6 and 2018 February 6), F110W (2017 December 9), and F160W (2017 December 7). Below we compare the HST light curve of AT 2017gfo to models of kilonova and GRB afterglow emission. We also consider emission components not previously observed in the optical or near-IR but that can be constrained by our new limits, including late-time changes in the afterglow spectral index (as in Balasubramanian et al. 2021;Hajela et al. 2021;Troja et al. 2022) and an IR dust echo (Lu et al. 2021). All model magnitudes, wavelengths, and dates discussed are given in the rest frame accounting for the redshift and luminosity distance to NGC 4993 (z = 0.00980 and 40.7 Mpc as in Cantiello et al. 2018) and assuming a merger time of 2017 August 17 12:41:04 (Abbott et al. 2017b).

The Kilonova Light Curve before 2017 August 30
The kilonova emission of AT 2017gfo dominated the observed UV and optical light at rest-frame epochs of 14.5 days (2017 August 30), although ground-based near-IR observations occurred up to 25 days from merger (Andreoni et al. 2017;Kasliwal et al. 2017;Pian et al. 2017), after which the field went into solar conjunction. The kilonova of AT 2017gfo was inferred to have multiple components to its UV to near-IR light curve following a thermal, radioactively  Table 1). We performed forced PSF photometry at the location of AT 2017gfo (noted with orange lines in all frames) as described in Section 2.1 and obtain detections of AT 2017gfo in F110W (top left), F160W (top right), F814W on 2017 December 6 (bottom left), and F814W on 2018 February 5 (bottom middle). In the remaining two subtractions, we calculate upper limits at the site of AT 2017gfo as described in Section 2.2. heated prescription as in Arnett (1982). These two components are commonly referred to as "blue" and "red" (see, e.g., Metzger et al. 2010;Kasen et al. 2013;Barnes et al. 2016), based on a rapidly declining component with opacity κ ≈ 1 cm 2 g −1 observed that dominated the spectral energy distribution (SED) within 2 rest-frame days of merger and a longer-lived component with κ ≈ 5-10 cm 2 g −1 observed at later times (Arcavi et al. 2017;Andreoni et al. 2017;Cowperthwaite et al. 2017;Drout et al. 2017;Lipunov et al. 2017;Kasliwal et al. 2017;Kilpatrick et al. 2017;McCully et al. 2017;Metzger et al. 2017;Nicholl et al. 2017;Smartt et al. 2017;Soares-Santos et al. 2017;Tanvir et al. 2017;Utsumi et al. 2017;Valenti et al. 2017;Villar et al. 2017). The characterization of the kilonova at these very early epochs was led by ground-based campaigns and the Neil Gehrels Swift Observatory (Swift), but at UV wavelengths AT 2017gfo had faded below the threshold detectable by Swift once HST started to observe (Figure 3).
The earliest set of HST observations occurred at ≈4.7-5.3 days after merger, spanning the UV to near-IR bands. First, these observations enabled the only detection at UV wavelengths (F336W) at >3 days from merger (whereas most of the UV detections before this epoch came from Swift/UVOT as in Cowperthwaite et al. 2017;Drout et al. 2017;Evans et al. 2017), providing an important anchor on the blue component. Subsequent HST observations also provided the deepest and latest constraints on this blue component emission. Second, HST probed the more slowly evolving, red kilonova emission with multiband observations continuing to ≈12 days.
To characterize the red and blue kilonova components of AT 2017gfo in the early-time HST data before 2017 August 30, we adopt the kilonova SEDs of Kasen et al. (2017). As in Kilpatrick et al. (2017), we use two models to characterize both kilonova components, consisting of a blue component with an ejecta mass M ej = 0.025 M e of ejecta moving with an ejecta velocity v k = 0.25c. The ejecta composition is broadly characterized by a lanthanide fraction X lan defined as the ratio of the mass of heavy lanthanide species (Z = 57-71) to the total ejecta mass (see, e.g., Kasen et al. 2013;. For the blue kilonova model, this corresponds to a radial gradient in X lan with X log 6 lan ( ) =in the outer layers and increasing to a higher but still low X log 4 lan ( ) =in the center (see Kasen et al. 2017;Kilpatrick et al. 2017, for details). The red kilonova model is simply characterized by M ej = 0.035 M e , v k = 0.15c, and X log 2 lan ( ) = -. The parameters used here for the two kilonova components are broadly consistent with other studies, and the corresponding light curves are well matched to the HST data.
We also combined the time-varying SEDs from these models and smoothed them over a window of 100 Å while masking out spectral bins where the flux drops to zero. We then passed the predicted spectrum at each time through the filter transmission functions for each ACS/WFC, WFC3/UVIS, and WFC3/IR band. We overplot these models for the bands in which AT 2017gfo was observed on the left side of Figure 3. In addition, in Figure 4 we show example SEDs at times when there were more than two broadband HST constraints on AT 2017gfo over a span of 1 day.
As shown in previous work on AT 2017gfo (Arcavi et al. 2017;Cowperthwaite et al. 2017;Drout et al. 2017;Kilpatrick et al. 2017;Smartt et al. 2017;Villar et al. 2017), the bestfitting model suggests that the total luminosity transitioned from blue to red kilonova dominated at around 3 rest-frame days from merger. Thus, all HST observations, which start at 4.74 rest-frame days, occurred at a time when the vast majority of kilonova emission was well characterized by our red kilonova model.
The only exceptions are the optical and UV bands, where a rapidly declining blue tail of the kilonova SED was observed from around 5 to 11 rest-frame days from merger (Figures 3 and 4). While the HST observations probed timescales when the red kilonova component was dominant, they also provided the only constraints on the rapidly declining blue tail. These data can uniquely probe this component of the kilonova at later times and its physical origin, which is still largely unsolved but may be due to energy injection from a long-lived NS remnant that lowers the electron fraction in the bluer ejecta (Lippuner et al. 2017), or possibly from accretion outflows from a disk that forms around the merger (Miller et al. 2019).

The GRB Afterglow Light Curve after 2017 December 6
After the field once again became observable with HST at >100 rest-frame days from merger, the optical and near-IR emission from AT 2017gfo was dominated by GRB afterglow (Lyman et al. 2018; Mooley et al. 2018;Troja et al. 2018;Fong et al. 2019;Lamb et al. 2019). Novel to this work are the late-time templates described in Section 2, which enabled four new detections in F814W, F110W, and F160W. To compare our updated photometry and upper limits from AT 2017gfo at these epochs, we compare its HST light curve to the afterglow model based on an off-axis relativistic structured jet and presented in Hajela et al. (2019). We adopt the updated parameters of Hajela et al. (2021) for a relativistic structured jet viewed at an angle of θ obs = 23°and interstellar medium density n 0 = 0.01 cm −3 . We choose these models for comparison over other afterglow models (e.g., JetFit models in Wu & MacFadyen 2018, with θ obs ≈ 30°) because the predicted observation angle is consistent with independent constraints from superluminal motion in the relativistic jet (≈20°in Mooley et al. 2018).
The resulting optical and near-IR light curves are shown on the right side of Figure 3, with the corresponding SEDs shown in  2019), we find no evidence for a change in spectral shape across the optical and near-IR SED (Figure 4). Our best constraints come from the afterglow light curve at 109.6 and 170.5 rest-frame days from merger, with two and three detections over a span of ≈2 days, respectively. In both cases, the observations are consistent with a constant spectral index of f ν ∝ ν −0.6 , reinforcing the broader constraints from radio to X-ray Figure 4. The SED of the kilonova (rest-frame days 5.06-11.31) and GRB afterglow (110.38-170.50 days) components of AT 2017gfo as constrained by HST detections and upper limits (circles) and described in Sections 3.1 and 3.2. The horizontal error bars correspond to the equivalent rectangular width of the corresponding filter as described in Rodrigo et al. (2012). We overplot the average kilonova and GRB afterglow models for data obtained within ± 0.5 days of the average day given next to each model. For the first model at 5.06 days (violet), there are two kilonova models from Kasen et al. (2017) within this time range, which are plotted as a shaded region between the brighter (upper) and fainter (lower) model. observations (Hajela et al. 2019) and a lack of any synchrotron curvature at intermediate frequencies.

Constraints on a Kilonova Afterglow
Recent detections of AT 2017gfo by Chandra at 0.3-10 keV suggest that the spectral index of the afterglow is changing even as it fades (Hajela et al. 2020;Balasubramanian et al. 2021;Hajela et al. 2021;Troja et al. 2022). The cause of this change in spectral index has been attributed to changes in the density of the circumburst medium, shock velocity, or microphysical parameters where the afterglow originates (Granot et al. 2018), or a new emission source arising from interaction between the slower kilonova ejecta and circumburst medium (called a "kilonova afterglow"; e.g., Hajela et al. 2019).
Although we do not detect AT 2017gfo in the latest HST observations, in principle these limits can be used to constrain the presence and nature of a kilonova afterglow (Kilpatrick et al. 2021b). Following a similar analysis to Hajela et al. (2021), we use each of our HST optical limits to constrain the optical to X-ray spectral index β OX (where F =-+ erg cm −2 s −1 . Based on our optical limits, the strongest constraints on β OX come from the F110W limit, which implies that β OX  0.85. This is moderately more constraining than the spectral index implied by F140W presented in Hajela et al. (2021) but significantly steeper than the spectral index implied by an afterglow spectrum with no evolution (∝ν −0.6 ), let alone the shallower spectrum implied by radio measurements (Balasubramanian et al. 2021). Without >2 mag rebrightening in the overall light curve of AT 2017gfo, it is unlikely that new optical observations will be constraining for the future evolution of this specific emission model.

Constraints on an Infrared Dust Echo
Our new late-time constraints on AT 2017gfo can provide unique constraints on IR dust echoes from the GRB afterglow (similar to those observed for supernovae and tidal disruption events; e.g., Graham et al. 1983;Jiang et al. 2021). IR dust echoes may be observed from the interaction between UV emission from the GRB afterglow and a sufficiently dense shell of dust proximate to the NS merger. As demonstrated in Lu et al. (2021), a shell of dust surrounding the merger would be sublimated up to some radius (r sub ) at which point radiation from the afterglow can no longer heat dust grains above their sublimation temperature. For an off-axis afterglow with a viewing angle θ obs , the timescale on which an IR dust echo from the heated dust grains becomes visible is therefore t r c 1 cos sub obs ( ) q = -, which for θ obs = 20°and r sub = 6 pc could be as long as 1.2 yr.
We show our limits compared with expected light curves for IR dust echoes in F110W, F140W, and F160W in Figure 5. The comparison light curves are constructed from the model presented in Lu et al. (2021). 11 They are parameterized primarily by the density of gas in the circum-merger medium, which was measured to be n H ≈ 10 −5 -10 −2 cm −3 depending on the exact afterglow model and microphysical parameters assumed Wu & MacFadyen 2018;Hajela et al. 2019;Wu & MacFadyen 2019). We note that we have assumed that the afterglow has an initial UV luminosity of L UV = 3 × 10 47 erg s −1 that peaks on a timescale of t 300 s max = as in Lu et al. (2021). For this model, our limits do not extend below the predicted flux for IR dust echoes propagating through an environment with n H = 1 cm −3 , and since L ν ∝ n H , we do not constrain this emission mechanism for realistic conditions around GW170817. However, we note that the predicted light curve peaks on later timescales in redder bands, and so observations with the James Webb Space Telescope may be able to detect an IR dust echo for GW170817 in the future.

Limits on Stars and Unresolved Stellar Clusters
We take advantage of the deep stacks that reach limits of 26.8-28.8 mag to place stringent limits on stellar sources and unresolved clusters. Previous F606W imaging ruled out bluer GCs with L  6.7 × 10 3 L e (Fong et al. 2019; see also Lamb et al. 2019), well below the peak of the F606W GC luminosity function (GCLF) determined for NGC 4993 (Lee et al. 2018). However, observations of nearby quiescent galaxies find that GC populations exhibit color bimodalities (Larsen et al. 2001), with red clusters corresponding to a metal-rich population (Brodie & Strader 2006). Indeed, V − H ≈ 0.3-2 mag was derived for the giant elliptical galaxy NGC 1399 (Blakeslee et al. 2012).
At the distance of NGC 4993, our deep limit (Section 2.2) of m F160W  26.8 mag corresponds to M F160W  −6.2 mag corrected for Milky Way extinction. To directly compare this to an H-band luminosity function and place constraints on an unresolved red GC, we create a representative GCLF by sampling the V F606 − H F160 GC color distribution of NGC 1399 (Blakeslee et al. 2012) and apply this color correction to the NGC 4993 GCLF in F606W (Lee et al. 2018). Our limits rule out redder clusters to L  2.3 × 10 4 L e or M  4.6 × 10 4 M e assuming a standard mass-to-light ratio of 2M e /L e ; only ≈0.05% of mass in red GCs is below this limit. Although this is nominally less constraining in luminosity than the previous F606W limit, we have placed an independent and deep constraint on reddened clusters.
To place these limits in the context of sources identified near GW170817, we obtain photometry of seven point-like sources within a 2″ radius (380 pc) of the site of AT 2017gfo as shown in Figure 6. These objects were selected from F160W to be point-like (dolphot object type 1 for "bright stars"; Dolphin 2016). Given these selection criteria, they tend to be redder than most cluster candidates in nearby galaxies and fainter than evolved massive stars (e.g., from Drout et al. 2009;Li et al. 2015). The brightest two of these sources are consistent with the low end of the GCLF, but otherwise they all appear too faint to be GCs.
We also show our limits placed at a distance of 200 and 330 Mpc, which are approximately the predicted detection ranges for binary NS mergers during Observing Run 4 and Observing Run 5 (LIGO Scientific Collaboration et al. 2015; Abbott et al. 2018). High-resolution optical mapping of future binary NS mergers using HST-like sensitivity can probe part of the GCLF and address their potential origin in GCs (Pooley et al. 2003;Rodriguez et al. 2016).

GALFIT Surface Brightness Modeling
Our deep stacks allow us to characterize the morphological substructure of NGC 4993 in great detail. To analyze this structure, we used the complete set of F606W, F814W, F110W, F140W, and F160W imaging, consisting of 63.9, 18.5, 20.0, 12.6, and 14.5 ks, respectively, at the site of GW170817. Based on the kilonova and afterglow light curves at the time each individual exposure was obtained, we expect AT 2017gfo to have an exposure-weighted magnitude of 23.6, 25.9, 22.3, 26.7, and 21.5 mag in our stacks. This is nominally above the deep stack detection thresholds from Section 2.2 in every band.
However, the combined kilonova and afterglow flux is much fainter than the galaxy as a whole, and we mask out emission from point sources, and so we do not expect this flux to significantly affect the quality of the model or analysis described below.
In order to isolate the shell structure, we use the GALFIT v3.0.5 software package (Peng et al. 2010) to fit two-dimensional surface brightness profiles to the smooth galaxy light from NGC 4993 with a Sérsic model. A Sérsic profile is parameterized by Sérsic index n, effective radius r e , PA, and ellipticity. When performing these fits, we mask out bad pixels and light from stars and galaxies in the field using the segmentation map of sources derived by combining the astrodrizzle image mask and running SExtractor v2.19.5 (Bertin & Arnouts 1996) to identify point-like sources. Due to the large spatial extent of NGC 4993 in these images, we use DETECT_THRESH=15 and BACK_SIZE=16 (size of the mesh in pixels over which the background is estimated) to prevent masking flux from the galaxy itself. We also perform a PSF deconvolution using the PSF model described in Section 2. First, we model the surface brightness profile of NGC 4993 in each of the F606W, F814W, F110W, F140W, and F160W deep stacks. We allow the centroid of the galaxy and aforementioned parameters to vary for a single Sérsic profile. Next, we undertake a two-component fit to better identify and characterize the large-scale substructure. For the primary component, we fix the center position from the previous profile and use the other fitted parameters as input values for a new fit. For the secondary component, we use a new Sérsic profile with the same fixed center position and all other parameters free to vary. The addition of a secondary component results in an improved residual map and χ 2 (by a factor of 2-8 in each band) compared to the single Sérsic case. F606W is the only band that is not well fit by a double Sérsic profile. We find that the inclusion of a Sérsic profile modified by a Fourier mode 2 (which distorts the shape of the 2D ellipsoid) for the secondary profile improves the fitting compared to standard Sérsic profiles in F606W. In the near-IR bands (F110W, F140W, and F160W), the fit is further improved by the addition of a third, PSF-like component at the center to model the presence of a weak active galactic nucleus (AGN).
Across all bands, we find that the primary component is characterized by n ≈ 2.9-4.4 (representative of a de Vaucouleurs profile for elliptical galaxies) and r e ≈ 14″-18″, resulting in r e ≈ 2.9-3.5 kpc at the distance of NGC 4993. The PA evolves from optical to near-IR, suggesting the presence of different superimposed stellar populations, as also found in Palmese et al. (2017) using an independent data set. The other primary component results are also consistent with fits previously performed in the literature of a single Sérsic profile (Blanchard et al. 2017;Levan et al. 2017;Palmese et al. 2017). The secondary component aims at modeling the core of the galaxy, which is not properly accounted for by the primary component. Inclusion of this secondary component enables a better fit to the galaxy profile and the detection of shell features closer (<20″) to the center of NGC 4993. Modeling of this core component is more challenging in the optical bands, where dust lane obscuration complicates the geometry of the galaxy more than in the redder bands. In the near-IR bands F110W, F140W, and F160W, the fitting prefers a cored, less cuspy component (i.e., a low Sérsic index n 2 ; see Table 2) close to a Gaussian or exponential disk. For comparison, we show the corresponding limit assuming a distance of 200 and 330 Mpc (dashed red lines). We place these limits in the context of photometry of GCs modeled from those in NGC 1399 (green circles; see description in Section 4 and Blakeslee et al. 2012), stellar clusters in the starburst galaxy M82 (blue circles; here we use F555W photometry as a proxy for V band from Li et al. 2015), and sources within 2″ (≈380 pc) of the site of AT 2017gfo (purple squares).
The residual images, as a result of subtracting the best-fitting light profile, are shown in Figure 7, and the GALFIT parameters for each frame are given in Table 2.

Identification of Shell Features around NGC 4993 and Analysis of Their Stellar Populations
NGC 4993 was classified as an S0 galaxy and part of the group LGG 332 (Garcia 1993) with at least two visible shells (it is also called MC 1307−231; Malin & Carter 1983). As previously shown in both ground-based (Palmese et al. 2017) and HST observations (Blanchard et al. 2017;Levan et al. 2017;Ebrová et al. 2020), NGC 4993 has clear shell structure in all of our deep HST stacked imaging (Figure 7). The apparently random distribution of shells around NGC 4993 is a characteristic of so-called "Type II" shell galaxies likely indicating a deviation from a perfectly radial orbit for the merger (as opposed to Type I shell galaxies with axisymmetric shells distributed in a double cone pattern; e.g., Hernquist & Quinn 1988;Sanderson & Helmi 2013;Bílek et al. 2015). Ebrová et al. (2020) identify 10 distinct shell features along the major photometric axis of the optical emission from NGC 4993 using stacked F814W imaging.
Here we use our GALFIT-subtracted residual imaging to systematically identify shell features at all PAs around NGC 4993. However, the quality of our GALFIT model results in large residuals and radial derivatives in those residuals close (<10″) to the center of NGC 4993. We ignore these features in the analysis below, and our detection of shells at the smallest radial separations from NGC 4993 only extends to features observed at the approximate radius of AT 2017gfo in Figure 7 (10 6). The shell features are the most prominent in the near-IR, and so we use the WFC3/IR F160W stacked frame with the largest IR spatial footprint to systematically identify shells around NGC 4993.
From the residual F160W stacked image, we also mask out emission due to stars and background galaxies using our SExtractor segmentation map of point-like features. To validate that the final shell image was relatively free of emission other than the shells, we visually inspected the masked image to ensure that there were no clearly detectable point-like features. We then filled in all masked pixels flagged by the SExtractor segmentation image with the median value for all nonmasked pixels within a 2″ radius of each masked pixel. Finally, we rebinned each frame into a grid of 1″ × 1″ pixels representing the median pixel value in each cell of the grid rescaled by the average number of pixels per 1 arcsec 2 cell. We repeat this process in each band, binning by 244 pixels for the WFC3/IR frames with 0 064 pixel −1 and by 400 pixels for the ACS and WFC3/UVIS frames with 0 05 pixel −1 .
Next, we performed aperture photometry using photutils to estimate the surface brightness of the shells in each cell of our grid. For F160W, the shell surface brightnesses range from 21.0 to 30.0 mag arcsec −2 , with brighter shells tending to be at smaller projected separations from NGC 4993. We also incorporate the statistical uncertainty from the GALFIT model into our shell flux estimates, which dominates the shell flux uncertainty in each cell of our grid. However, while analysis of the effect of varying GALFIT parameters on our results is beyond the scope of this paper, we acknowledge that systematic uncertainties in the GALFIT model possibly dominate the total shell flux and colors discussed below.
Finally, we segmented the gridded map of shell surface brightnesses into 24 individual shell features using a method analogous to CLUMPFIND (Williams et al. 2011) by identifying groups of local maxima in 4″ × 4″ rectangular apertures. From the 24 local maxima we identified, we segmented the remaining pixels into features by identifying adjacent pixels whose surface brightnesses are within 0.3 mag of the local maximum. We iterated on the previous step until no remaining pixels could be found within 0.3 mag of the previous step.
We emphasize that these shell features are segmented arbitrarily and do not necessarily correspond to distinct "shells" as defined in shell-type galaxies or previous analysis of NGC 4993 (e.g., Malin & Carter 1983;Garcia 1993;Ebrová et al. 2020). Thus, we refer to these structures (which appear to form a double spiral structure in Figure 8) as "shell features" for convenience. We order each shell feature from largest to smallest projected separation in Table 3, with that order overlaid on the F160W image in Figure 8. These features do not have a one-to-one correspondence with the 10 shells identified in Ebrová et al. (2020), in part because we use F160W as opposed to F814W, we consider shell features at all PAs around NGC 4993 rather than those only along the major photometric axis, and our imaging is significantly deeper. However, we are confident that this census of shell emission is complete to the depth of our masked residual F160W image.
The location of each shell, its photometric properties, and stellar population properties of each shell derived from Prospector (see discussion in the Appendix) are given in Table 3. The total stellar mass in the shell features is M 6.3 10 1.5 2.6 8  -+ , or approximately 1.4%-2.5% of the total mass in stars depending on the stellar mass estimate in Note. Morphological parameter estimates for a double Sérsic profile in the different bands using GALFIT. The subscripts 1 and 2 refer to the first and second Sérsic component, respectively. In the IR bands, we include an additional PSF component to account for the presence of a weak AGN as described in Section 5.1. This component is characterized solely by the integrated magnitude (m PSF ) provided here in AB mag. The effective radius (r e ) for both components is provided in kpc assuming D = 40.7 Mpc to NGC 4993. The model ellipticity is parameterized by the ratio between the semiminor and semimajor axes of the model ellipsoid (b/a). The effective radii r e are provided in pixels, the position angles (PA) in degrees east of north.
NGC 4993 as a whole ((3.0-4.5) × 10 10 M e in Blanchard et al. 2017;Palmese et al. 2017;Ebrová et al. 2020). We consistently find a median mass-weighted age across the shell features of >3 Gyr, similar to the smooth galaxy light profile of ≈5-11 Gyr. We note that these statistics are limited by our constraints on the SED blueward of F606W. In particular, the lack of deep UV observations limits our constraints on recent star formation in the shells. The poorer constraints on metallicity and recent star formation history (SFH) are a known limitation in Prospector analyses derived from broadband photometry in the absence of UV photometry (see also Leja et al. 2017). The utility of additional UV imaging is, however, limited by the inability of standard stellar population synthesis models to produce the "UV upturn," a well-known observational feature in these systems (e.g., Yoon et al. 2008;Yi et al. 2011;Vazdekis et al. 2016) thought to be caused by binary interactions (Han et al. 2007).
We therefore caution that our analysis of the shell population SFH is subject to significant systematic uncertainties on timescales in the youngest age bins (<100 Myr; see the Appendix). The conclusions we draw below are robust to these uncertainties given the low total mass observed in these age bins. showing the galaxy's shell structure in F606W, F814W, F110W, F140W, and F160W. All images are centered on AT 2017gfo and represent the same 140″ × 147″ region. We indicate the location of AT 2017gfo in each frame with a blue star. Note that in all frames light from the diffraction pattern of a bright star (USNOA2 0600-15448796; m V = 12.6 mag) is visible on the right-hand side and is unassociated with any shell emission. We mask out these features in our analysis of the shell features as described in Section 5.2. Figure 8. GALFIT-subtracted F160W image of NGC 4993 with masking and binning into 1″ × 1″ cells as described in Section 5.2. We segment this image into 24 shell features ordered from smallest to largest projected separation from the center of NGC 4993 (red circle), corresponding to the features given in Table 3. The location of AT 2017gfo is indicated with a blue star. However, at most 1.3% of all the stars in the shell features are <262 Myr based on the nonparametric SFH described in the Appendix. The maximum star formation rate (SFR) averaged across all of the shell features is approximately 0.04 M e yr −1 in the oldest (>4.7 Gyr) age bin (Figure 9). This is on the low side for the star formation main sequence at z > 1 (e.g., Noeske et al. 2007) but still within ≈2σ of expectations. The SFR declines monotonically until it levels off in the youngest two age bins at the current day rate of ≈0.02 M e yr −1 , roughly half of its peak rate.

Properties of the Galaxy Merger of NGC 4993
We analyzed the time since first impact of the merger for the shells identified around NGC 4993 in the context of the shell radius−age relations derived in Ebrová et al. (2020) based on the radial profile model from Palmese et al. (2017). This relation assumes that accreted galaxies plunge into their host galaxies on radial trajectories and stars stripped from the infalling galaxies then move along close-to-radial orbits. Thus, the shells correspond to overdensities of stars near the apocenters of their orbits (see also Quinn 1984;Dupraz & Combes 1986). Due to the energy gradient in the satellite galaxy, the kinetic energy of the stars forming the shells, and thus the apocenter radii move outward with time, allows us to place a lower limit on the time from accretion by calculating the time it would take stars to reach the outermost shell in the gravitational potential of NGC 4993.
From the largest projected radius at which we detect any shell emission, we can calculate the look-back time since the merger, where a look-back time of zero corresponds to the redshift of NGC 4993. The outermost shell radius of 71 8 (or 14.2 kpc) corresponds to a minimum look-back time of 220 Myr, comparable to the 200 Myr estimate in Ebrová et al. (2020). We do not detect any significant emission beyond this radius in our deep stacked F160W image, and so we are confident that there are no shell features to deep luminosity and mass limits (<10 6 M e ) based on our Prospector analysis.
The SFH derived in the Appendix for the shell features levels off after our 262-685 Myr bin, and 2.5 1.3 2.5 -+ % of the total stellar mass is formed at this time. We infer from this finding that the merger occurred at most 685 Myr ago, largely depleting the secondary galaxy of star formation material and resulting in a small fraction of younger stars in the shell structure. This finding and the small population of stars with the youngest ages are still consistent with a relatively "dry" merger as noted by Levan et al. (2017). The presence of some additional gas from the merger is also indicated by the weak AGN observed in Chandra and radio imaging (Blanchard et al. 2017) and in the WFC3/IR bands. Our limit on the time since the merger is corroborated by the upper limit derived from dynamical models performed in Ebrová et al. (2020), which suggest that the Table 3 Shell Locations, Photometry, and Stellar Population Properties Note. Properties of the 24 shell features we present in Section 5. Here shell # corresponds to the features shown in Figure 8. The projected separation (r sep ) and position angle (PA; east of north) are averaged over each 1″ × 1″ cell in the F160W map and weighted by the F160W flux. Both quantities are computed with respect to the center of NGC 4993 in our deep F160W image (red circle in Figure 8), which we measure to be R.A. = 13:09:47.70, decl. = −23:23:02.305 (J2000) in the Gaia DR2 astrometric frame (Lindegren et al. 2018). S F160W is the average surface brightness across each shell in F160W. We also indicate the average colors of each shell, weighted by the flux in F160W, in F606W, F814W, F110W, and F140W with respect to F160W. M å and T å represent the total stellar mass and median mass-weighted stellar age, respectively, as determined from our Prospector fits in the Appendix and averaged across each shell feature. number and separation between the inner shells are consistent with a merger timescale of <600 Myr. Although our limit on the age is less constraining, our independent method directly probes the stellar population that is likely formed in the merger, and so we consider the upper age limit to be <685 Myr below.
Based on the morphological substructure of NGC 4993 and the lower and upper limits on the time from merger we derived, we consider the true origin of the binary NS progenitor of GW170817. Using the total estimated stellar mass of the shell features as a proxy for the secondary galaxy in the merger (i.e., the smaller galaxy that was accreted by NGC 4993), we find a minimum mass ratio of ≈1:50, classifying the past event as a minor merger. We note that it is possible that we are not completely characterizing the stellar mass of the secondary component owing to geometric orientation of the shells and our view in projection, a small fraction of the secondary's stars not being located in the observed shells (i.e., not located close to the apocenter of their orbit), and the limits of our HST observations preclude detection of the lowest surface brightness features. Still, none of these reasons would account for a major increase in the mass ratio.
From simulations, the stellar mass ratios of shell-forming features are typically at least 1:10, but somewhat lower ratios are occasionally observed (see, e.g., Pop et al. 2018, Figure 8). The incidence of shells is also more common in isolated galaxies and known to decrease in groups or rich clusters (Colbert et al. 2001), and NGC 4993 is part of a group (Garcia 1993).
We conclude that based on the stellar mass in shells, it is highly improbable that the binary NS progenitor of GW170817 originated from the low-mass secondary galaxy and instead originated from the primary. However, based on the typical delay times of binary NS mergers, which extend to several Gyr, and the minimum time since merger of 220 Myr, it is fully plausible that the binary NS progenitor was formed long before merger. In this scenario, the galaxy-scale merger may have affected the ultimate trajectory, orbit, and merger timescale of the progenitor of GW170817.

Conclusions
We present the full set of HST observations of GW170817/ AT 2017gfo obtained to date, including new template observations obtained from 2021 January 4 to February 22. The full HST data set, representing over 140 ks of wide-band observations obtained in a range of filters and instruments, is one of the deepest sets of observations of a 40 Mpc galaxy ever obtained with HST. These observations have enabled four new detections of the nonthermal afterglow from 2017 December 7 to 2018 February 6, as well as new limits on the presence of optical and near-IR sources around AT 2017gfo and extended shell emission in NGC 4993 as a whole.
The new detections of the nonthermal GRB afterglow, which span from 109 to 170 rest-frame days post-merger, remain consistent with an unchanging spectral index of β ≈ −0.6. However, similar constraints on the evolution of NS merger counterparts out to later times (>100 days) can yield insight into their origin and the processes driving their emission mechanisms, including kilonova afterglows (Hajela et al. 2021), magnetar-boosted events (Fong et al. 2021a), IR dust echoes (Lu et al. 2021), and variations in the circum-merger density (e.g., Ramirez-Ruiz et al. 2019).
The limits on sources near AT 2017gfo derived from the full HST data set are also significantly constraining for nearby GCs, which have been suggested as a potential origin for the binary NS progenitor system of the merger (Belczynski et al. 2018;Baillot d'Etivaux et al. 2019;Ye et al. 2020). We can detect all GCs, including those with the reddest colors, down to 4.6 × 10 4 M e , representing the large majority of the GCLF. Future binary NS mergers, especially those near the distance limit of the LIGO during Observing Runs 4 and 5, will likely not have such strong limits on the presence of coincident GCs. However, even at 200 Mpc, limits similar to those for AT 2017gfo can rule out 27% of the GCLF, precluding an origin from the most massive and luminous GCs.
The total HST data set also enables the most complete sample of low surface brightness features around NGC 4993, and in particular the shell structure first noted in Malin & Carter (1983). These data provide an observational blueprint for a Type II shell galaxy, which can be studied irrespective of its connection to binary NS mergers. Using F606W through F160W imaging, we are able to identify shell structure out to ≈71 8, the most complete census of such emission around NGC 4993 to date. Fitting this photometry with Prospector stellar population models, we constrain the total stellar mass in the shells to be 6.3 × 10 8 M e , approximately 2% of that in NGC 4993 as a whole, with a mass-weighted stellar age across all of the shells >3 Gyr. The geometry of the shells supports an age from 220 to 600 Myr (Ebrová et al. 2020), and the SFH of the shells supports a maximum time since merger of ∼685 Myr. Given the lack of evidence for a very young stellar component in the shells, we consider it unlikely that the progenitor of AT 2017gfo originated in this stellar population.
We thank W. Lu for providing the IR dust echo models presented in this paper and Jay Strader for helpful discussions. Support for program No. 15886 was provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. The Fong Group at Northwestern acknowledges support by the National Science Foundation under grant Nos. Figure 9. The SFH averaged across all shell features in NGC 4993. We include the median SFR in M e yr −1 (black line), as well as the limits indicated by the 16th and 84th percentile most likely SFR (gray region), roughly corresponding to our 1σ uncertainty. The SFR declines monotonically with time from its peak rate of ≈0.04 M e yr −1 and begins to level off in the two youngest bins with ages of <262 Myr. For reference, we draw a dashed line at the point where the predicted SFR is 1/2 the peak rate where this leveling-off begins. Shell Stellar Population Modeling We constrained the properties of the stellar population in the NGC 4993 shells using the stellar population inference code Prospector v1.0.0 Johnson et al. 2021). This analysis is based on the broadband photometry from each of our deep stacked and GALFIT-subtracted images in Section 5. Following the masking and gridding procedure described in Section 5.2, we gridded each subtracted image into 4″ × 4″ cells. We then performed photometry and identified 268 such cells in which there was flux detected at >3σ in the F606W, F814W, F110W, F140W, and F160W frames. This set of photometry formed the basis for 268 independent Prospector fits.
We fit this photometry using a nonparametric SFH with seven age bins, with the first two spaced from 0 to 30 Myr and from 30 to 100 Myr and the remaining five log-spaced in time with an upper limit at the age of the universe (13.63 Gyr) at the redshift of NGC 4993. We assumed an initial mass function from Chabrier (2003). In addition, we adopt a continuity prior as described in Fong et al. (2021b) and originally presented in Leja et al. (2019) such that the SFH does not sharply deviate from a flat distribution, but otherwise we do not place any constraints on SFH.
We used this SFH model to fit for a total stellar mass (  M M log  ) and metallicity (  Z Z log  ). We also include two dust components to model the attenuation of light due to dust in stellar birth clouds and affecting only young stars (dust1) and attenuation of light in the diffuse interstellar medium (dust2, parameterized as dust2/dust1). Both parameters can be interpreted as the additional optical depth due to each dust component (see, e.g., Conroy 2013;Kriek & Conroy 2013;Price et al. 2014). This model is a power-law perturbation from the dust attenuation curve in Calzetti et al. (2000), and so we fit a power-law index for the attenuation curve for dust1 and dust2, referred to as δ, for which we adopt a flat prior from −1.0 to 0.4. For all 268 Prospector fits, we assume the fixed distance and redshift given above, and our input photometry was corrected for Milky Way extinction. We performed the fit in each of the 268 cells by jointly fitting all five bands in that cell with the nested sampling routine dynesty v1.0.1 . The in-band magnitudes are inferred for each WFC3/ UVIS or WFC3/IR band using python-fsps v0.4.1rc1 Conroy & Gunn 2010). We adopt an error floor of 5% in the photometric uncertainties to avoid overfitting to the shell flux in each cell, some of which have photometric uncertainties <0.01 mag.
To combine the output quantities for each cell into individual shell features as indicated in Table 3 and Figure 8, we concatenated the sampled parameters for the Prospector fit in each cell. Each stellar mass in the concatenated samples was rescaled by the total stellar mass inferred across all of the cells we included in that shell feature.
Finally, we note that the 4″ × 4″ cells do not span the same solid angle as the 1″ × 1″ cells we used to segment the shell emission in Section 5.2. This is due to the fact that we require a >3σ detection in all five photometric bands, and the shells are detected at higher significance in the F160W band. Therefore, we rescaled the total stellar mass in each shell feature by the ratio between the F160W flux in the 1″ × 1″ map and the same flux in the 4″ × 4″ map. This ratio is 1.1-2.5 for each shell, with features at larger projected separations tending to have larger ratios. The final stellar mass inferred for each shell feature is given in Table 3, along with the median massweighted age.