Hubble Space Telescope Observations of Tadpole Galaxies Kiso 3867, SBS0, SBS1, and UM461

Gas accretion onto galaxies drives their average star formation rates (SFRs; Finlator & Davé, 2008; Brooks et al., 2009; Combes, 2014), with excess local accretion possibly connected to major bursts (Sánchez Almeida et al., 2014a, b; Ceverino et al., 2016). Tadpoles (van den Bergh et al., 1996; Elmegreen et al., 2005), or cometary galaxies (Markarian, 1969; Loose & Thuan, 1985; Cairós et al., 2001), are gas-rich dwarfs with one giant star-forming region that looks like the head of a tadpole when the galaxy is viewed at some inclination; the rest of the disk comprises the tail. This morphology is common at high redshift, where it is found in $\sim 10$% of resolved galaxies disk-l(Elmegreen et al., 2005; van den Bergh et al., 1996; Straughn et al., 2006; Windhorst et al., 2006), but it is rare locally at less than 0.2% (Elmegreen et al., 2012). Based on simulations by Ceverino et al. (2016) of low metallicity gas infall onto disk galaxies, which typically triggers a burst of star formation near one end of the disk, we expect that tadpole galaxies would appear disk-like if viewed face-on, with at least one prominent star-forming region. Their baryon mass is dominated by gas (e.g., Filho et al. (2013)), and their stellar mass distribution tends to be triaxial with axial ratios 1:0.7:0.4 (Putko et al., 2019).

Tadpole galaxies in the local Universe (out to $\sim 100$ Mpc) tend to be extremely metal poor (XMP) (Morales-Luis et al., 2011). Spectroscopic observations of 22 tadpoles indicate their heads typically have a drop in metallicity by a factor of 3 to 10 compared with their tails, consistent with accretion of metal-poor gas (Sánchez Almeida et al., 2013, 2014a, 2014b, 2015; Olmo-Garcia et al., 2017; Lagos et al., 2018), a process that may be generally important in XMP systems (e.g., Ekta & Chengalur, 2010). This relation between excess SFR and drop in metallicity is observed in most low mass galaxies (e.g., Hwang et al., 2019; Sánchez-Menguiano et al., 2019).

One of the prototype tadpole galaxies is Kiso 5639. Its bright head is a large, over-pressure, short-lived star-formation region with a rich population of star clusters (Elmegreen et al., 2016) and a massive molecular cloud (Elmegreen et al., 2018). In our previous HST imaging of it (Elmegreen et al., 2016) we found the properties of the head to be consistent with a model by Ceverino et al. (2016) where infalling gas forms a low-metallicity starburst clump in the disk.

In this paper we present an HST imaging study in U, V, I, and H$\alpha$ of four additional tadpole galaxies, SBS0 (KUG0940+544), SBS1 (SBS1129+576), Kiso 3867 (KUG0937+298) and UM461. These were selected from previous tadpole spectroscopic studies (Sánchez Almeida et al., 2015, 2013) for showing clear metallicity drops in the star-forming heads, as discussed further below. Galaxy coordinates, distances, masses, and rotational speeds are listed in Table 1.

Our objective is to study the properties of the dominant XMP star forming complexes and how they correlate with the characteristics of host galaxies. We are also interested in the populations of young star clusters formed in these bursts, and in their mass and age distributions. For example, the localized starbursts in these tadpoles could have been triggered by infalling gas streams, and in other galaxies where this has happened, such as NGC 1569 (Johnson et al., 2012) and NGC 5253 (Turner et al., 2015), there are superstar clusters near the points of impact. We would like to know if there are outlier superstar clusters in these tadpoles or only a smooth extrapolation of the usual power-law cluster mass function (Elmegreen & Efremov, 1997).

All four galaxies have lower metallicities in the bright star-forming regions than elsewhere. The metallicity drops in SBS0 (J094416.6+541134.4), SBS1 (J113202.4+572245.2), and UM461 (J115133.3-022221.9) were shown by Sánchez Almeida et al. (2015) using the long-slit spectrograph OSIRIS at the 10.4 m Gran Telescopio Canarias¹¹1http://www.gtc.iac.es/GTChome.php. The metallicity drop in Kiso 3867 was shown by Sánchez Almeida et al. (2013) using the spectrograph IDS of the 2.5 m Isaac Newton Telescope at the Observatorio del Roque de Los Muchachos.

SBS0, SBS1 and Kiso 3867 are morphological tadpoles, which suggests they are inclined disks. UM461 is presumably more face-on (inclination $\sim$30$\pm$ 10^∘; van Zee et al., 1998). UM461 is also unusual in having 2 luminous, compact, young stellar regions. The Hi line map of UM461 by van Zee et al. (1998) shows an extended Hi disk, typical of blue compact dwarf (BCD) systems, with no indication of a tidal interaction with the nearby UM462 system. The outer Hi shows evidence for filamentary structures. Similar features in other BCDs are suggested as possible indications of ongoing gas accretion (e.g., IC10; Wilcots & Miller (1998); Namumba et al. (2019), but see also comments by Lelli et al. (2014)).

Ekta et al. (2006) mapped SBS1 in the Hi 21 cm line using the Giant Metrewave Radio Telescope (GMRT). SBS1 is interacting with its more massive companion galaxy SBS 1129+577 and an Hi bridge joins the two galaxies. The interaction, however, seems unlikely to be the source of low metallicity material in SBS1 since the more massive and optically brighter SBS 1129+577 galaxy is expected to have higher metallicity than SBS1. Kiso 3868 (AGC 194068) was observed by Pustilnik & Martin (2007) using the Nançay telescope with Arecibo in the ALFALFA Hi survey (Haynes et al., 2018) with log(M_HI/M${}_{\odot})=7.6$ and M_HI/L${}_{B}\approx$3 (Durbala et al., 2020). The SBS0 system has only a weak upper limit on the Hi mass that still allows it to be a gas-rich system (Filho et al., 2013). In summary, the tadpoles in our sample resemble typical low luminosity BCDs in terms of their Hi properties and include the SBS1 interacting system and 3 examples of apparently single dwarf tadpole galaxies.

Our observations and data reduction methods are presented in Section 2, the analysis of large star-forming clumps, which include clusters and other young stars, is discussed in Section 3, the star clusters themselves are considered in Section 4, cluster mass and age distribution functions are in Section 5, H$\alpha$ emission is in Section 6, and the conclusions are in Section 7.

The four tadpoles in our sample were imaged with the Hubble Space Telescope for Cycle 27 GO proposal 15860 using the WFC3/UVIS and ACS cameras between November 2019 and April 2021. The filters used were F390W (U), F547M (V), and F657N (H$\alpha$) with WFC3/UVIS, and F814W (I) with ACS, with 1 orbit for each filter per galaxy, for a total of 16 orbits.

The 5$\sigma$ detection limits in AB mag from the HST exposure time calculator (ETC) for point sources are 27.2, 26.6, and 27.1 for filters F390W, F547M, and F814W, respectively. The 5$\sigma$ limit for the H$\alpha$ filter F657N is $4.5\times 10^{-17}$ erg cm² s^-1.

The WFC3/UVIS images were custom calibrated using the routines developed in Rafelski et al. (2015) and Prichard et al. (2022). First, the images were corrected with a pixel based correction for charge transfer efficiency (CTE) degradation to account for the radiation damage over time (MacKenty & Smith, 2012). We utilized the new improved CTE algorithm (Anderson et al., 2021) with reduced noise amplification and improved corrections. We used the new co-temporal UVIS darks to apply a custom hot pixel mask using a variable threshold as a function of distance from the amplifier readout to detect a constant number of hot pixels across the CCDs²²2https://github.com/lprichard/HST_FLC_corrections. We flagged readout cosmic rays, which are cosmic rays that land on the CCD while being read out, and are otherwise over-corrected by the CTE code which would result in negative divots in the images. We then equalized all four amplifier regions to produce images with constant background levels³³3https://github.com/bsunnquist/uvis-skydarks. Finally, due to the small number of exposures per filter, we applied a custom cosmic ray rejection based on the Astropy implementation of the LACOSMIC procedure⁴⁴4https://github.com/mrevalski/hst_wfc3_lacosmic (van Dokkum, 2001; McCully & Tewes, 2019).

The resulting custom calibrated science exposures were aligned to GAIA DR2 (Gaia Collaboration et al., 2018) with TweakReg and combined with AstroDrizzle with a 0.8 pixel fraction at a final pixel scale of 0.03^′′. Root mean square (RMS) uncertainty maps were created from the inverse variance (IVM) images as 1/sqrt(IVM), and were corrected for correlated pixel noise (Fruchter & Hook, 2002).

Color composite figures were made in the Image Reduction and Analysis Facility (IRAF) using RGB in DS9, with the F390W (for blue), F547M (for green), and F814W (for red) images. These are shown in Figures 1 and 2 for linear and logarithmic stretches, respectively. As traced by the relative blue colors, star formation is prominent in the heads of SBS0 and Kiso 3867. The head is not as prominent in SBS1, which also has some bright star-forming regions along its tail. UM461 has two large star-forming regions in the E and W portions in the northern half of the galaxy. Star clusters are evident in several regions throughout the galaxies, discussed below. The logarithmic stretches show extended emission well beyond the brightest parts of the disks. In UM461, the southwest part of the galaxy is extended away from the star-forming peaks and a ridge containing young stellar populations extends westwards from the north side of the system.

To quantify the continuum associated with the H$\alpha$ images, a linear interpolation of the F547M and F814W filters was performed. First, each broad-band image was normalized to the H$\alpha$ image by dividing it by the ratio of the corresponding PHOTFLAM values stored in the image header. This bring the three images into a common photometric system, where a source with constant flux per unit wavelength would have exactly the same number of counts. Aperture photometry was then carried out for those stellar objects in the field visible in the three filters, excluding those objects too faint in the H$\alpha$ image to give meaningful results, and those that appeared to be saturated in the F814W frame. There were only two or three useful stellar objects for each galaxy, so that a continuum fit on a galaxy by galaxy base was not feasible. Instead, the results for the four galaxies were merged together, for a total of 10 useful sources, and the best fit was computed to the relation:

where $f$ are the fluxes (on the normalized images) derived from the aperture photometry. The best-fit coefficients are: $a=0.524\pm 0.139$ and $b=0.456\pm 0.130$, with the uncertainties computed using the bootstrap method. Equation (1) assumes the calibration stars to have negligible H$\alpha$ emission and an average spectral energy distribution (SED) similar to the stellar parts of the galaxies.

As an example, Figure 3 shows the H$\alpha$ image of SBS1 before and after continuum subtraction. The H$\alpha$ image was calibrated considering the throughput of the whole optical system, including the F657N filter, at the observed H$\alpha$ wavelength.

Photometric measurements were made of the large-scale star-forming complexes, or clumps, larger than 20 pixels ($0.6^{\prime\prime}$) in diameter for each galaxy using the task imstat in IRAF. Boxes were defined around the clumps based on isophotal contours in the F390W image at $20\sigma_{sky}$, labeled in Figure 4. The background was subtracted from each clump measurement to derive the clump magnitudes. Photometric zeropoints calculated from the WFC3 website⁵⁵5http://www.stsci.edu/hst/wfc3/phot_zp_lbn and ACS website⁶⁶6https://www.stsci.edu/hst/instrumentation/acs/data-analysis/zeropoints were used to convert counts to AB mag.

Ages and extinctions for the tadpole clumps were determined by fitting stellar population models to the observed spectral energy distributions, using the (U-V) and (V-I) colors. The observed colors were corrected for background by subtracting in each filter the counts for nearby fields from the counts in the clumps. Uncertainties in the corrected counts were set equal to the square roots of the sums of the standard deviations in the clump and background counts, weighted by the number of pixels for each, and the relative uncertainties were taken to be these corrected uncertainties divided by the total background-subtracted counts in the clump. The relative uncertainties in the colors were taken to be the square roots of the sum of the squares of the relative uncertainties in each passband. On average, the uncertainties in the colors, converted to magnitudes, were 0.69 mag for (U-V) and 0.71 mag for (V-I).

The intensities were modelled with single stellar population models at 0.2 solar abundance and a Kroupa IMF in Bruzual & Charlot (2003) using the full model spectrum at each trial age multiplied by $\exp(-\tau)$ for wavelength-dependent dust opacity $\tau$ and by the wavelength-dependent filter throughput. We used extinction curves from Calzetti et al. (2000) and Leitherer et al. (2002) with $A_{V}$ ranging from 0 to 4 mags in steps of 0.1 mag. Throughputs for the UVIS filters are from the WFC3 website⁷⁷7http://www.stsci.edu/hst/wfc3/ins_performance/throughputs/Throughput_Tables. The two observables, (U-V) and (V-I), were thus converted to two physical quantities, age and extinction.

The resulting intensities for each passband were integrated over time for population ages back to some start time, which is identified with the age of the region. This integration assumes a constant star formation rate. The population models assume a Chabrier IMF with a metallicity of 0.2 solar. Trial models varied the population age from $10^{6}$ to $10^{10}$ years in 16 equal steps of $\log({\rm age})$. For each age and extinction, the model colors were set equal to $-2.5$ times the log of the ratio of intensities.

The best simultaneous fit for age and extinction was determined by minimizing the differences between the model and observed colors. One method weighted all of the age and extinction values by $\exp(-\chi^{2}/2)$ within a limited range of $\chi$ for $\chi^{2}$ equal to the sum of the squares of the differences between the observed and model colors divided by their respective uncertainties. We then used the average of these weighted values as the fitted solution. A second method took the $\chi^{2}$ weighted average of the age plus extinction values with the lowest root mean square differences between the model and observed colors, binning these differences in steps of 0.1 mag. The uncertainties for both methods were taken to be the rms uncertainties in these averages. Both methods gave about the same results, so here we report those from the latter method.

After the ages and extinctions were fit from the colors, the masses of the clumps were determined by setting the observed absolute intensity in the I band to the model absolute intensity, using the mass versus age tabulation in Bruzual & Charlot (2003).

Table 2 lists the masses, ages, and extinctions derived from model fits to the photometry for the large-scale clumps in the galaxies, as well as SFR and stellar surface density. The star formation rates listed are the ratios of the total stellar masses formed (not the present masses) divided by the ages.

The clump average log (Mass) ($M_{\odot}$) is 4.1$\pm 0.7$. The total log mass in clumps in each galaxy is 5.4, 5.3, 5.0, and 5.7, in listed galaxy order. The clump average log (Age) (yrs) is 6.7$\pm 0.4$, with a range from 6.1 to 7.4. Extinction A_V ranges from 0.3 to 1.1 mag for all of the clumps, with an average of 0.6$\pm 0.2$ mag. The log SFR average for the clumps is -2.5$\pm 0.8$. Histograms showing the distributions for all 4 galaxies for mass, age, and SFR are shown in Figure 5. The sizes of the clumps (from the square root of the area of the clumps) range from 16 to 162 pc, averaging 65$\pm 47$ pc. Stellar surface densities range from 1.1 to 26.5 $M_{\odot}$ pc^-2, with an average of 7.0$\pm 5.7$ $M_{\odot}$ pc^-2. For comparison, the background surface densities range from 22 to 97 $M_{\odot}$ pc^-2, with an average of 43.5$\pm 29.7$ $M_{\odot}$ pc^-2.

For comparison, Lagos et al. (2011) used Gemini/NIRI J, H, K images of UM461 to measure star complexes and modeled their masses and ages. Their clumps 2 and 5 in their figure 4 correspond to our biggest clumps 1 and 2, respectively. Their (log (Age), log (Mass)) for these are (6.3, 5.6) and (6.4, 4.6) compared with our values of (6.7, 5.6) and (6.5, 4.8), so the infrared and optical results are similar. Their clumps 1 and 9 are similar to our 4 and 3 also. In addition, they find several clumps with log mass ($M_{\odot}$) $\sim 4.6$ that do not stand out in our optical images. In particular, their clump 4 has no obvious optical counterpart, indicating obscuration south of our clump 1.

Our simultaneous fitting of age and extinction to only two broadband colors, (U-V) and (V-I), produces uncertainties in both quantities, especially since they compensate for each other with older ages giving redder colors like higher extinctions. To check the inaccuracies, we forced all the fits for star complexes to a constant extinction of $A_{V}=0.1$ mag, which is less than the fitted values. As expected, the ages increase and the masses decrease a little. For the four galaxies SBS0, SBS1, Kiso 3867 and UM461, respectively, the averages in the log (Mass) ($M_{\odot}$) decrease from 4.4 to 4.1, 4.4 to 4.2, 3.7 to 3.5, and 4.7 to 4.4. The averages in the log (Age) (yr) increase from 6.4 to 6.8, 6.5 to 6.9, 6.8 to 7.4, and 6.9 to 7.6. These represent a 5% decrease in log (Mass) and a 8% increase in log (Age). The star formation rates obtained by the ratios of mass to age decrease with these changes, by an average factor of 5.8 (-0.76 in the log) for all galaxies if $A_{V}$ is forced to equal 0.1 mag.

The extended diffuse components of the galaxies are expected to be older than the clumps, which photometric fits confirm. These outer diffuse regions are highlighted in the logarithmic stretch images in F814W shown in Figure 6. SBS0 shows extensions in the NE and SW, SBS1 shows extensions in the N and S, Kiso 3867 shows extensions in the W and particularly in the E, and UM461 shows a long extension in the SW. Rectangular boxes devoid of obvious clumps or star clusters were measured for background subtraction of the clumps, described earlier, and to determine surface brightnesses and ages. A histogram of their ages is shown in Figure 7; their average log (Age) (yrs) is 9.3$\pm 0.4$.

Cluster mass distribution functions are typically power laws with a slope of around $-1$ on a log-log plot, corresponding to a function $n(M)dM\propto M^{-2}dM$ for equal intervals of mass, M (Krumholz et al., 2019). Of interest for these galaxies is whether the mass functions have a statistically significant deviation from a power law at high mass, either a decrease as in a Schechter function, or an increase because of high-mass cluster outliers.

The cluster age distribution function is also interesting as this results from a combination of the cluster formation and destruction rates.

In the present study, there are too few clusters to determine a statistically significant probability distribution function for cluster mass in each galaxy, so a composite function was made for all the galaxies together by suitable normalization, assuming that this function is a power law. A composite age distribution function, or cluster decay rate, was determined as well. As before, we included only clusters with photometric uncertainties in the F390W band that are less than 0.2 mag, giving cluster numbers of 38 for SBS0, 96 for SBS1, 135 for Kiso 3867, and 267 for UM461.

Figure 10 shows the distribution of cluster mass versus age for each galaxy. The characteristic lower limit to the mass at each age is from the limiting apparent magnitude. Concentrations in age around 4 Myr are due to this age being the lowest available age in our model grid, whereas concentrations at 100 Myr and 1 Gyr are due to the variation of (V-I) color with population age. At $\sim 100$ Myr, the (V-I) color changes rapidly due to the appearance of the first red supergiants, followed by a roughly flat color evolution between $10^{8}$ and $10^{9}$ years. When the color changes rapidly, a large range of color due to measurement uncertainties corresponds to a small range of age, concentrating the clusters at that age. In practice, due to the limited SED coverage of our observations, the SED fitting likely favors an age at either the young or the old limits of this color-plateau.

To make a composite mass-age diagram, we multiply all the cluster masses by $(10/D)^{2}$ for galaxy distance $D$ in Mpc. This increases the masses for closer galaxies and decreases them for distant galaxies, making the detection limit the same for each, as if they were all at 10 Mpc.

After this, the lower limits to the mass-age distributions are all the same and they can be stacked into a composite diagram. This composite is shown in Figure 11. At each age, the vertical distribution of points is the cluster mass function. At each mass, the horizontal distribution of points is the age distribution, which for a constant formation rate is the survival distribution function, i.e., the number of clusters of a particular mass that survive as a function of age. Typically the vertical density of points tapers off like $\sim 1/M$, which means the mass function is $dn/dM\sim M^{-2}$, and the horizontal density of points is approximately constant or increasing slightly, meaning the cluster number decreases like $1/$age or shallower (Whitmore et al., 2007).

The cluster mass function for distance-normalized cluster masses was determined in two ways, first for a few narrow age intervals and then by stacking these intervals with appropriate vertical shifts. The age intervals were taken to be the three discrete ages at small age, namely, 4, 8, and 12 Myr, and the two other intervals with age concentrations, $\log({\rm age})=$ 7.6 to 8.3, and $\log({\rm age})=$8.8 to 9.2, for ages in years. These five intervals give five mass functions, which are each fairly noisy because of the small numbers of clusters in each. The left-hand panel of Figure 12 shows the resulting composite mass functions in the five age intervals, where each one combines all the galaxies by the distance normalization. There is a systematic shift toward higher masses for older clusters from the rising lower limit to the mass-age plot in Figure 10.

The age intervals can be combined by shifting the masses of the younger clusters upward by an appropriate amount. Stacking the result essentially assumes the mass function is a power law, but it does not assume any particular power. To do this stacking, we multiply all the cluster masses in each age interval by the ratio of the minimum cluster mass in the $\log({\rm age})=7.6$ to 8.3 interval to the minimum in that interval. This makes all the minimum masses have the same minimum as the 7.6 to 9.3 log-age interval. All five shifted mass functions are then stacked to make one composite mass function. The result is shown on the right in Figure 12.

Figure 12-right shows a smooth power law mass function with a slope of $-1.15\pm 0.17$ (90% confidence interval by the student-t distribution). This agrees with the idealized case resulting from wide-spread hierarchical structure in star formation, where the slope on a log-log plot would be -1 (Elmegreen & Efremov, 1997). There are no outlier masses beyond the power law at the high end, and no drop-off there either, as in a Schechter function (see also Mok et al., 2020).

The average age distribution can be found in a similar way, assuming the actual distribution is a power law in age, by taking the age distributions in narrow distance-normalized mass windows and shifting them all to have the same maximum value, again following the detection limit in the mass-age diagram. The distance-normalized mass intervals are given by log (Mass) $(M_{\odot})=2.5-3$, $3-3.5$, $3.5-4$, $4-4.5$ and $4.5-5$. The result is shown in Figure 13. The slope of the distribution is $0.24\pm 0.21$ (90% confidence interval by the student t-distribution). This slope implies that the number of clusters decays with time as $dn(t)/dt\propto t^{-0.76\pm 0.21}$.

Multiwavelength imaging of four tadpole galaxies with HST WFC3 and ACS reveal a pattern similar to that previously found for Kiso 5639. The three largest galaxies contain one or two giant stellar complexes with masses of $10^{5}\;M_{\odot}$ or more and ages of a few million years. The smallest galaxy, Kiso 3867, has a maximum complex mass of $\sim 5\times 10^{4}\;M_{\odot}$ and a similarly young age. There are typically smaller complexes nearby that are also young.

The complexes are accompanied by large populations of compact young star clusters. The mass distribution function for the combined clusters was determined by a double-normalization technique, which assumes only that the intrinsic function is a power law. In this technique, the cluster masses were first scaled as if all the galaxies were at the same distance. Then the cluster masses in five age bins were scaled a second time to make their distributions have the same lower mass cutoff from sensitivity limits. The composite mass function determined in this way was indeed a power law, confirming the initial assumption, and the slope on a log-log plot was $-1.15\pm 0.17$, similar to the expected value of $-1$ for hierarchical distributions. There were no significant deviations from this power law at high mass, i.e., no outliers with masses higher than the extrapolation of the function, and no suggestion of an upper mass cutoff.

The composite age distribution function was determined in a similar way, although there is no a priori reason for this to be a power law. Combining all masses between $10^{2.5}\;M_{\odot}$ and $10^{5}\;M_{\odot}$, the result showed a slight increase in cluster counts with log (Age) $t$, as $t^{0.24\pm 0.21}$, which corresponds to a loss of clusters over time as $t^{-0.76\pm 0.21}$.

The tadpole galaxies in our sample share the common feature of reduced metallicity in the most prominent star-forming complex as compared with its surroundings, as determined in previous papers. This suggests that recent accretion of low-metallicity gas triggered the main starburst (e.g., Sánchez Almeida et al., 2013). The gas could accrete in several ways. The off-center locations of the main complexes could be the result of tadpole motions relative to a surrounding lower-metallicity medium if the head is more compressed than the rest of the disk by ram pressure. Alternatively, the metal-poor gas could be accreted from a gaseous halo or cosmic filament with non-zero specific angular momentum. The quiescent kinematics of most tadpoles appears to rule out fueling from minor mergers. Aside from this accretion, these tadpole galaxies appear to be typical blue compact dwarfs.

The summed theoretical Lyman continuum (LyC) emission from the clusters is comparable to the LyC needed to excite the observed H$\alpha$, but the summed LyC emission including regions around the clusters is more than 20x larger. This excess LyC suggests a loss of LyC photons, which could be explained by dust absorption although these galaxies have low metallicities. Alternatively there could be additional H$\alpha$ emission in a halo that we do not see, or some escape of LyC photons from the galaxy.

Acknowledgments:

The imaging data presented in this paper were obtained from the Mikulski Archive for Space Telescopes (MAST) at the Space Telescope Science Institute. The specific observations analyzed can be accessed via https://doi.org/10.17909/kjp2-sh47 (catalog 10.17909/kjp2-sh47). Based on observations with the NASA/ESA Hubble Space Telescope obtained at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555. Support for HST-GO-15860 was provided through a grant from NASA/STScI. JSA, CMT, and NC acknowledge support from the Spanish Ministry of Science and Innovation, project PID2019-107408GB-C43 (ESTALLIDOS), and from Gobierno de Canarias through EU FEDER funding, project PID2020010050. This research has made use of the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration.