Free Access
Issue
A&A
Volume 659, March 2022
Article Number A2
Number of page(s) 17
Section Extragalactic astronomy
DOI https://doi.org/10.1051/0004-6361/202141590
Published online 25 February 2022

© ESO 2022

1. Introduction

A common morphological property of high redshift star-forming galaxies is the presence of clumps (Elmegreen et al. 2007; Förster Schreiber et al. 2011; Wuyts et al. 2013; Bournaud 2016; Guo et al. 2018; Zanella et al. 2015, 2019; Vanzella et al. 2021). Strong gravitational lensing easily reveals such clumps down to a ∼100 pc scale (Livermore et al. 2015; Rigby et al. 2017; Cava et al. 2018; Dessauges-Zavadsky et al. 2019, 2017) and the sizes of massive stellar clusters (≲30 pc) in high magnification regimes, μ > 20 (Vanzella et al. 2019, 2020a; Johnson et al. 2017). Probing the presence of clumps at the smallest physical scales remains crucial to understand the evolution of high-redshift galaxies and it is only recently that observations have been capable of achieving the resolution needed to detect the predicted substructures and thus understand their formation (Faure et al. 2021; Elmegreen et al. 2021, and the references therein).

Searching for the presence of single star clusters within high redshift clumps represents the next frontier and, certainly, a key science goal of future extreme adaptive optics facilities (XAO, eXtreme Adaptive Optics, e.g., Vanzella et al. 2021). Nowadays, gravitational lensing gives us the opportunity to access such small scales (< 100 − 200 pc). The truth is, however, that even after relatively high amplification, μ < 10 − 20, it is still a challenge to probe scales smaller than a few tens of pc that are key for the identification of single star clusters. Therefore, higher magnification regimes are required, which, in turn, are more rare and difficult to handle as they are generally subject to high modeling uncertainties. Now, new-generation lensing models, based on an unprecedented number of constraints, are providing a significant advances in this direction (e.g., Bergamini et al. 2021), especially in the vicinity of (or on the) critical lines (Vanzella et al. 2020b).

The system described in this work (dubbed Sunburst in the following) belongs to the category of the super-lensed arcs, subject to amplification values spanning the range of 15 < μ < 100. Sunburst is a strongly magnified galaxy at redshift 2.37 discovered by Dahle et al. (2016), gravitationally lensed by the Planck cluster PSZ1 G311.65–18.48 (at z = 0.44) and split over four very bright arcs (with the brightest arc having magnitude R ∼ 18). Several multiple images (more than 50 in total) of the same galaxy and internal star-forming knots populate the giant arcs. Among them, a compact star-forming region, identified 12 times and showing a structured Lyα profile with multi-peaks and with non-zero emission at the systemic redshift, was identified as a candidate region with Lyman continuum leakage (optically thin to ionizing photons, LyC, λ < 912 Å, Rivera-Thorsen et al. 2017). Subsequent observations with Hubble confirmed the LyC emerging from the same 12 multiple regions, in agreement with the exceptionally high image multiplicity (Rivera-Thorsen et al. 2019). Chisholm et al. (2019) accurately estimated an age of 3 Myr, sub-solar metallicity of 0.5 Z and E(B − V) = 0.15, while Vanzella et al. (2020a) derived the dynamical age from the limits on the stellar mass and size of the system (as a function of magnification), recognizing it as a plausible massive young star cluster with a stellar mass of ∼107 M, an effective radius smaller than 25 pc, and showing evident spectral signatures of massive O-type and Wolf-Rayet stars from the VLT/MUSE spectrum. The resulting initial constraints on the star formation rate and stellar mass surface densities of the LyC emitting knot also resemble those of local gravitationally bound star clusters (Vanzella et al. 2019, see their Sect. 4.1.1). Sources like Sunburst are unique laboratories for improving our understanding of the physical mechanisms related to escaping radiation at cosmological distances. This makes them very useful in the search for ionizing sources at earlier epochs during cosmic reionization (z > 6), when the LyC radiation cannot be observed directly due to high cosmic opacity rapidly increasing in the first ≃1.5 Gyr (Worseck et al. 2014; Romano et al. 2019). They are also key for making comparisons with LyC emitters identified in the Local Universe (including candidates), considered analogs of those with high redshift (Izotov et al. 2021a; Schaerer et al. 2016; Jaskot et al. 2019; Wang et al. 2021).

An accurate model of the lens is required to unveil the intrinsic properties of such sources; this type of model is described in a companion paper (Pignataro et al. 2021). Here, we present refined estimates of the physical quantities of the LyC star cluster, including the other knots identified in the surrounding, and we derive the global properties of the host galaxy.

This work is organized as follows. Section 2 describes the data used and Sect. 3 presents the Sunburst lensed system, the lens model, the emergence of stellar clusters, and an estimate of the dynamical mass for the stellar cluster showing Lyman continuum emission. Section 4 describes the galaxy hosting the aforementioned star clusters, along with its intrinsic appearance. In Sect. 5, we compute the cluster formation efficiency and in Sect. 6, we compute its global (along line of sight) escape fraction of ionizing photons. Section 7 summarizes our results and describes the implications drawn from possible analogies with sources at z > 6.5 that are likely responsible for reionization (Sect. 7.1). We also present the current limitations and future prospects in that section (Sect. 7.2).

We assume a flat cosmology with ΩM = 0.3, ΩΛ = 0.7 and H0 = 70 km s−1 Mpc−1. With this assumption, one arcsec at z = 2.37 corresponds to a projected physical scale of 8200 parsec. Unless specified, all magnitudes are given in the AB system.

2. Data

The galaxy cluster PSZ1 G311.65–18.48 was observed in several bands with the Hubble Space Telescope between 2018–2020 under the programs 15101 (PI Dahle), 15949 (PI Gladders), and 15377 (PI Bayliss). We retrieved the WFC3/F390W (total integration time 5.8 ks), WFC3/F410M (13 ks), WFC3/F555W (5.6 ks), WFC3/F606W (5.9 ks), ACS/F814W (5.3 ks), WFC3/F098M (1.4 ks), WFC3/F105W (1.3 ks), WFC3/F140W (2.7 ks), and WFC3/F160W (1.3 ks) exposures from the MAST archive. We processed them with the GRIZLI software1 following the procedure described by Williams et al. (2021). Each of the groups of exposures in a given filter and epoch (a “visit” in the Hubble observation nomenclature) was aligned to the Gaia DR2 absolute astrometric frame (Gaia Collaboration 2016, 2018) and accounting for proper motions. The PSZ1 G311.65–18.48 field has a high density of Gaia sources and the absolute astrometric precision of all image mosaics is ∼10 milliarcsec (1σ). The morphological analyses below are performed on mosaics in each filter drizzled to a common 30 milliarcsec pixel grid with the DRIZZLEPAC software (Gonzaga et al. 2012).

Additional ground-based observations obtained with the VLT/X-shooter (PI Vanzella) are used in this work and presented in Appendices A and B. A comprehensive analysis of the full spectrum of Sunburst, from Lyα to Hα wavelengths, will be presented in a forthcoming work. VLT/X-shooter data acquisition and reduction on Sunburst have already been presented in Vanzella et al. (2020c) and described in Appendix A, where we focus on key observables useful for the analysis described in this work.

MUSE observations in the narrow field mode configuration of a portion of the lensed system are also presented in this work (PI Vanzella, see Appendix C), by highlighting the exceptional performances of the extreme adaptive optics, as well as some limitations in cases such as the lensed systems described here.

3. Anatomy of the Sunburst: Detecting star clusters at z = 2.37

3.1. Lensing model: A very amplified system

The work by Pignataro et al. (2021) presents the details of the lens model, along with the predicted positions of all 62 multiple images of various galaxies and the associated families currently confirmed in the redshift range of 1 < z < 3.5 (out of which 54 belongs to the Sunburst galaxy at z = 2.37). The high level of accuracy achieved for this new model is underscored by its ability to predict the observed positions of 62 multiple images with a rms error of only 0.14″.

We took advantage of this lens model to recognize each mirrored and multiply imaged star-forming knot of Sunburst. The geometry of the lens shows a rather complicated behavior of the critical lines at the redshift of Sunburst due to several perturbers placed along the line of sight, which add to the underlying mass distribution of the galaxy cluster. Such a configuration boosts the number of multiple images and total magnification values (μtot), with μtot spanning the range between 15 and 100. Figure 1 shows the magnification map along with the position of the relevant arcs discussed in this work. For example, one of the knots (labeled 5.1 in Figs. 1 and 2) is detected 12 times due to a combination of the galaxy cluster and perturbers. The set of multiple images of a single knot is called a “family” and each member of the family is denoted with a different letter. For example, each of the 12 images of 5.1, 5.1(a–n) is subjected to a different magnification. Along the giant arcs, the total magnification (denoted as μtot) is dominated by tangential stretch, such that the ratio between the total magnification μtot and the tangential one, μtang, is nearly constant and close to 1.3 − 1.4. The ratio μtot/μtang is referred to the radial magnification, which is therefore very low on the arcs. Figure 2 shows the panoramic view of the four arcs: I, II, III, and IV that contain several portions of the same galaxy. Thirteen individual SF knots replicate several times with different amplification values, according to the complex geometry of the lens. The list of knots is reported in Table 1, presenting those highlighted in green labels in Fig. 2.

thumbnail Fig. 1.

Panoramic view of the main arcs discussed in this work. Left panel (in F555W-band image): Arcs are marked with red contours to guide the eye. The arrows indicate 11 multiple images of the star-forming knots labeled 5.1a–m (see text and Fig. 2 for additional details). The magnification map (total magnification, μtot) is reported on the right panel (from Pignataro et al. 2021) coded according with the color bar in square root scale: it saturates at μtot ≃ 10 (black) and 100 (white). The black contours enclose the region where μtot > 500, marking the locus of the critical line. The arcs contain multiple images of the same portions of the Sunburst galaxy and are subjected to amplification spanning the range 15 − 100, with the critical line crossing the arcs several times.

thumbnail Fig. 2.

Identification of the multiple images of compact star-forming knots is shown. In the bottom-left of the figure, the negative gray-scaled HST ACS/F555W image shows the panoramic view of the four arcs, labeled as I, II, III, and IV, with the little red square in the middle indicating the source region subject to large amplification. The SF knots are shown in the blue border insets labeled as I, II, III, and IV, following the arcs notation in the bottom left inset and with increasing magnification from bottom to top (long arrow on the right). Each inset includes the negative gray scale F555W image and a color image obtained by combining the WFC3/UVIS F275W (blue), WFC3/UVIS F606W (cyan), ACS/WFC F814W (yellow), and WFC3/IR F160W (red) filters (ESA/Hubble, NASA, Rivera-Thorsen et al., CC BY 4.0). Star-forming regions are indicated following the nomenclature of Pignataro et al. (2021). White dashed lines mark the position of the critical lines. The main triplet, 5.1, 5.2, and 5.3 are indicated with blue, red, and magenta open circles, respectively. The objects labeled with green IDs are the knots reported in Table 1. Those marked in yellow show all the identified multiple images of the knots. The inset in the bottom right (red contour) shows the reconstructed image of the Sunburst galaxy on the source plane, where the triplet 5.1, 5.2, and 5.3 would merge into a single HST pixel (∼30 mas separation). The observed angular separation between the knots 5.1 and 5.3 is reported on the right side of each inset, and the typical physical scales subtended by the FWHM (0.1″) along tangential direction is shown over the long black arrow to the right (arc II includes star-forming region smaller 10 pc). The transient (“Tr”) discussed in Vanzella et al. (2020c) and corresponding to the knot 5.10 is indicated with a black arrow in the inset II.

Table 1.

Summary of the identified SF knots observed on Sunburst.

Statistical errors on magnification have been calculated by extracting the total(tangential) magnification μtot (μtang) at the model predicted positions (of the star-forming knots discussed below) over 500 realizations of the lens model by sampling the posterior probability distribution function. This has been performed by using a Bayesian Markov chain Monte Carlo (MCMC) technique (as described in Bergamini et al. 2021). Total (μtot) and tangential (μtang) magnifications and the associated errors are reported in Table 1. The physical quantities reported in Table 1 are calculated based on the best-fit lens model and the uncertainties over the grid of 500 magnifications. Overall, such large magnifications imply that a scale of a few tens of pc can be probed along the tangential stretch, allowing us to identify single star clusters (discussed in the next section).

3.2. Emergence of star-forming regions at parsec scale

To provide the most stringent constraint on the size of each knot, we consider the brightest one within the associated family, characterized by a high magnification. All the identified multiple images of the knots, except for knot 5.1, appear as unresolved objects in the F555W-band (1600 Å rest-frame), independently of their magnification value. The F555W is the bluest band probing the ultraviolet continuum (not including Lyα emission or IGM attenuation). Under these conditions, this band is also the one with the narrower PSF, thus, it offers the best probe (spatial contrast) of the knots and therefore enables optimal estimates of the sizes in the ultraviolet. The current HST data in the near infrared bands do not allow us to address the sizes in the optical rest-frame with the same level of detail. As discussed in Sect. 7.2, this will be performed with ease thanks to future facilities based on XAO.

We modeled the two-dimensional light distribution with Galfit. Galfit fitting has been performed by assuming the Sérsic index n = 0.5 (Gaussian), 1.0 (exponential) up to n = 4 (de Vaucouleurs) light profiles. Given the circular symmetric shapes of the knots, we fix the axis ratio (b/a) and position angle to 1 and 0, respectively. The inferred sizes always converge to compact solutions, independently of the adopted Sérsic index, formally with PSF-deconvolved radii Reff < 1 pixel for all of them (1 pix = 0.03″), but one, the brightest 5.1l which has an Reff = 2.0 ± 0.6 pix. The modeled images and residuals of the modeling (observed – model) are shown in Fig. 3. Table 1 summarizes the inferred intrinsic magnitudes resulting from the fit and the physical sizes (effective radii), along with the magnifications and associated uncertainties (statistical errors). The estimated radii span the range < 22 pc, which likely represents upper limits (given the sources are unresolved), approaching sizes of a few pc (< 10 pc) for the cases very close to the critical lines. It is worth noting that all knots show a point-like shape also in the F390W band, implying a half width at half maximum (HWHM) of 0.04″ (as inferred from stars in the field), which is consistent with the small upper limits on the effective radii derived above with Galfit in F555W-band.

thumbnail Fig. 3.

Galfit fitting of the star-forming knots reported in Table 1 zooming on various regions (groups A, D, G, H and L, see Fig. 2). Left column: original F555W image, in which only the knots listed in Table 1 are marked with their IDs. Middle column: Galfit models and right column: residual after subtracting the models, with the positions of the knots indicated as in the left column. Besides the IDs in the reference sample (marked in blue), additional multiply imaged knots have been modeled, such as the cloud of clusters of the group D (see the same group described in Fig. 2). All knots are not resolved with resulting effective radii ≲1 pixel, independently from the Sérsic index (n = 0.5 − 4).

3.3. Cocoon of gravitationally bound star clusters

The small sizes inferred from modeling (Reff ≲ 22 pc) might suggest that such star-forming knots are stellar clusters. In order to investigate such a possibility, we consider whether such compact star-forming knots are gravitationally bound systems by calculating the dynamical age, Π (equal to the ratio of the age and crossing time).

The dynamical age Π can be calculated as described in Gieles & Portegies Zwart (2011) (see also Adamo et al. 2020a; Vanzella et al. 2021). In particular, the crossing time expressed in Myr is defined as T CR = 10 × ( R eff 3 / G M ) 0.5 $ T_{\mathrm{CR}} = 10 \times (R_{\mathrm{eff}}^{3}/GM)^{0.5} $, where M and Reff are the stellar mass and the effective radius, respectively, with G ≈ 0.0045 pc−3 M 1 $ M_{\odot}^{-1} $ Myr−2 is the gravitational constant. Stellar systems evolved for more than a crossing time have Π > 1, suggestive of being bound (Gieles & Portegies Zwart 2011). This criterion has been used extensively for the identification of star clusters in the local Universe (e.g., Calzetti et al. 2015a; Adamo et al. 2017; Ryon et al. 2017). It is valid under the following assumptions: the system is in virial equilibrium, follows a Plummer density profile, and the light traces the underlying mass. Therefore the calculation of Π for each knot requires the knowledge of three main parameters: age, stellar mass, and effective radius.

Age. The estimation of recent star-forming events on short timescales would require age indicators such as the Balmer lines, which are highly uncertain – despite the high magnification – since a detailed, spatially resolved spectroscopy of each star-forming knot is not yet achievable. On the other hand, the SED-fitting weakly constrains young ages (being essentially dependent on the assumed star formation histories, e.g., Carnall et al. 2019) and is affected by blurring in the near infrared bands, where the PSF degrades and the emission from the knots and the host galaxy cannot be easily deblended. Nonetheless, an estimate of the average ages comes from VLT/X-shooter observations (Vanzella et al. 2020c). As mentioned above, individual knots cannot be probed due to the relatively coarse resolution. However, we estimated the Hα equivalent width of a group of 10 knots (group “d”, Fig. 2), obtaining spatially-averaged limits on the age, assuming that stars formed with an instantaneous burst (see Appendix A). For knot 5.1, the Hα emission is prominent with respect to the underlying continuum, implying an age younger than 5 Myr and consistent with the 3 Myr inferred by Chisholm et al. (2019). For the other knots of the group “d” instead, the Hα EW smaller than 100 Å rest-frame suggests an average age larger than ≃7 Myr. In the following, we adopt an age of 3 Myr for 5.1 (Chisholm et al. 2019) and assume the other knots are older than 7 Myr (see Fig. A.1 and Table 1). As discussed below, we note that even when relaxing the ages to smaller values, Π still remains close to 1 or higher. Moreover, it is worth noting that assuming different SF histories (continuous or declining) would result in older ages of the knots (Zanella et al. 2015), further increasing the dynamical age Π.

Stellar mass. For the same reasons described above, the estimation of the stellar mass from SED fitting suffers from the resolution issue on PSF-matched images. In addition, the apertures the photometry is calculated within, on the basis of multiple images of the same family, are (lens) model-dependent. This makes the choice of the aperture shape over multiple images tricky. Such an issue will need dedicated simulations with forward modeling techniques, which will be described in a forthcoming work. In this work, we rely on the ultraviolet F555W-band image, which has a narrow PSF (of 0.1″), such that the majority of knots are well recognized (Fig. 2).

We assume all the knots formed through instantaneous bursts (as is the case for 5.1). We derive the stellar masses by comparing the dust-corrected UV luminosity with Starburst99 models (Leitherer et al. 2014). We adopt a dust extinction of E(B − V) = 0.15 as derived by Chisholm et al. (2019) for knot 5.1 (corresponding to A1600 ≃ 1, Reddy et al. 2016), whereas for the remaining knots we adopt the average E(B − V) inferred from SED fitting (see below) of all multiple images of the same family. In fact, while the colors extracted from the PSF-matched images are preserved, the derivation of the stellar masses from SED fitting is strongly limited by the aforementioned issues (see also Cava et al. 2018). The stellar masses have therefore been inferred from the dust-corrected F555W luminosity, assuming the aforementioned ages (3 Myr for 5.1 and 7 Myr for the rest), and span the interval 105 − 7 M, weakly dependent on the assumed metallicity: they vary by ≃0.18 dex if the metallicity spans the interval of Z = 0.001 − 0.02 (values reported in Table 1 assume Z = 0.008, with Z = 0.02 indicating solar metallicity in the Sunburst99 nomenclature). It is worth noting that the UV-based stellar mass for 5.1 is consistent with the inferred dynamical mass discussed below (Sect. 3.4.1). Despite the aforementioned limitations related to crowding and aperture size definition, we check the consistency between the stellar mass inferred for knot 5.1 with both SED fitting and the single ultraviolet band F555W. In fact, the expected biases mentioned above (i.e., the crowding and lens-dependent aperture photometry among multiple images tends to be reduced) if the knot is bright and compact, therefore dominating the signal almost independently from the chosen aperture. The two methods have been applied to eight multiple images of the same knot 5.1 (a,b,c,h,i,l,m,n). Knots 5.1d,e,f,g have been excluded from the computations because affected by a foreground perturber which might bias the estimate of the magnification and affects the photometry of the images (Pignataro et al. 2021). SED fitting has been performed on 0.2″ diameter apertures, adopting a metallicity of 0.4 Z (the value in the BC03 library closer to the metallicity of the youngest star cluster 5.1, which has likely inherited the metallicity of the host galaxy), E(B − V) in the range 0 − 1.1 and age < 50 Myr. The extracted quantities (e.g., luminosity and stellar masses) are compatible within multiple images and comparable among the two methods within the relevant uncertainties (the median stellar mass inferred from single 1600 Å band and SED fitting are Mass(UV) = 1.1(±0.5)×107 M and Mass(SED) = 0.9(±0.6)×107 M, respectively.

Sizes. As described in Sect. 3.2, we verify the compactness of the knots by fitting the light profile with Galfit. We obtain upper limits on the radii smaller than the half width at half maximum. All of them have radii smaller than 22 pc, whereas 5.1 is marginally resolved with a radius of ≃8.5 pc.

All the quantities with uncertainties are reported in Table 1 along with the dynamical ages. The resulting dynamical age Π exceeds 1 (gravitationally bound) for the majority of the knots, suggesting they indeed are likely to be bound stellar clusters. It is worth calculating for each of them the minimum age above which Π exceeds 1. Since the UV-based stellar mass depends on the assumed age, we recalculate the mass at each adopted age, starting from 1 Myr and increasing it by 0.1 Myr till the condition Π = 1 is reached. The last column of Table 1 reports the age and stellar mass at Π = 1. Even in the case of younger ages (1−3 Myr) most of the knots are still bound. Therefore, the age limit (> 7 Myr) based on the Hα equivalent width discussed above indicates that all knots are likely bound star clusters, including the knot 5.1, which is the most extreme case of a 3-Myr-old YMC with detected LyC leakage.

3.4. LyC emitting cluster

3.4.1. Dynamical mass estimate

We estimated the dynamical mass of the knot labeled 5.1 by taking advantage of the observed [O III]λ5007 line and the inferred velocity dispersion from high spectral resolution X-shooter spectroscopy, R ≃ 5600 (Rhoads et al. 2014, and see also Sect. 4.4 of Vanzella et al. 2017). The line profile and double Gaussian component fitting is shown in Fig. B.1 and the resulting dynamical mass of 107 M is discussed in Appendix B. The resulting ratio between the stellar (photometry-based) and the dynamical masses is ≃1 and weakly depends on the magnification μ. This is due to the fact that the photometric and dynamical masses scale with total and tangential magnifications, μtot and μtang, respectively. In the case of Sunburst, the tangential magnification largely dominates over the radial one, μrad, such that μtot/μtang = μrad ≃ 1.3 − 1.4, nearly constant along the arcs. The mass ratio, therefore, only includes μrad, which is well constrained (with < 10% uncertainty, being far from the radial critical line). The resulting stellar and virial masses are 1.1 × 107 M and 107 M, respectively. This indicates that the total mass is dominated by the stars, as was also found in most local bound clusters (e.g., McLaughlin & van der Marel 2005).

The combination of the stellar mass (107 M), age (3 Myr) and size (Reff ≃ 8 pc) implies a star formation rate (ΣSFR) and stellar mass (ΣM) surface densities quite large for the knot 5.1, Log10SFR)≃3.7 M yr−1 kpc−2 and Log10M)≃4.1 M pc−2, consistent with the values observed in local YMCs (e.g., Bastian et al. 2006; Östlin et al. 2007). A constant star formation history would imply a SFR of ∼3.3 M yr−1 and the sSFR ≃ 330 Gyr−1. Interestingly, such large densities and sSFR are lower limits since the SF history was likely not constant with time. The presence of such YMC influences the overall observed appearance (morphology and photometry) of the host galaxy (see Sect. 4); in this case, the star cluster is well recognized only because of strong gravitational lensing that boosts the observed flux and significantly increases the spatial resolution.

3.4.2. Escape velocity

Young massive star clusters are powerful producers of ionizing radiation, with an enormous sSFR and forging of hot and massive stars that demonstrate they can provide sufficient stellar feedback in such an early phase to perforate the interstellar medium of the galaxy. Observations of local young massive clusters indicate that even massive systems (M ≃ 105 − 106 M) appear essentially gas-free even at very young ages (a few Myr, Cabrera-Ziri et al. 2015), indicating that strong feedback processes must already be at play in the earliest phases, in some cases before the onset of SN explosions. Again, the YMC hosted in the Sunburst galaxy and the low-column-density of neutral gas along the line of sight (NHI < 1017.2 cm−2 given by a large Lyman continuum escape fraction) may represent such a phase of intense stellar feedback caught in the act at 3 Myr since the initial burst of star formation.

We estimate the escape velocity vesc of the YMC 5.1 from the following equation (Cabrera-Ziri et al. 2016):

v esc = f c M cl R eff km s 1 , $$ \begin{aligned} { v}_{\rm esc} = f_{\rm c} \sqrt{\frac{M_{\rm cl}}{R_{\rm eff}}}\;\mathrm{km\,s^{-1}}, \end{aligned} $$(1)

where Mcl is the cluster mass, Reff the effective radius, and fc accounts for the dependence of the escape velocity on the density profile of the cluster. The value of fc ranges between the minimum and maximum, namely, 0.076−0.130 (from Table 2 of Georgiev et al. 2009). Adopting Mcl = 107 M and Reff = 8 pc as determined above and the highest value of fc = 0.130, we get an upper limit to the escape velocity of the cluster, vesc ≃ 145 km s−1. The clear asymmetric blue tail of the [O III]λ5007 line (shown in Fig. B.1) implies an outflow component along the line of sight. Following Perna et al. (2015), we derive (1) the maximum velocity defined as the velocity at 2% (v02) of the cumulative flux, F ( V ) = a b F v ( v ) d v $ F(\mathrm{V}) = \int_{a}^{b} F_{\mathit{v}}(\mathit{v}^{\prime})\mathrm{d}\mathit{v}^{\prime} $ (where FV is the line profile in the velocity space and the position of v = 0 of the cumulative flux is set at the systemic redshift given by the narrow component, λ = 16878.6 Å, or z = 2.3702, adopting 5008.24 Å [O III] vacuum wavelength) and (2) the line width w80, that is, the width comprising 80% of the flux defined as the difference between the velocity at 90% (v90) and 10% (v10) of the cumulative flux of the entire line profile. The maximum velocity and the w80 are 380 and 260 km s−1, which are also consistent with the FWHM of the broad component of the double Gaussian fit mentioned above, 320 km s−1. Overall, such outflow velocities are consistent with what is inferred by Rivera-Thorsen et al. (2017) on the same object based on absorption lines of silicon and eventually higher than the escape velocity of the cluster inferred above. This also agrees with the scenario in which the most plausible source for driving the escape of ionizing radiation is stellar wind and radiation from hot massive stars (Izotov et al. 2018; Heckman et al. 2011), which are certainly hosted in the cluster core (Vanzella et al. 2020a).

4. The Sunburst galaxy

In this section, we investigate the physical properties and the appearance of the Sunburst galaxy on the source plane in the context of the aforementioned presence of star clusters. In particular, we emphasize the significant contribution of the YMC 5.1 to the optical nebular emission and the ultraviolet appearance of the entire galaxy.

4.1. Physical properties

SED-fitting has been performed on the image “m” of the galaxy (dubbed HostM hereafter), which offers the closest view of the entire Sunburst galaxy (see Fig. 2, arc III, image “m”). This is indeed the least magnified image of the target and therefore the one that resembles more closely the intrinsic shape and morphology of the entire galaxy. We use the HST images F390W, F410W, F555W, F606W, F814W, F098W, F105W, F140W, and F160W available in the HST archive (see Sect. 2). In order to obtain unbiased photometric measurements, we first estimated the point spread function for all the images using bright, unsaturated stars and then we PSF-matched all the bands to the F160W image using relevant convolution kernels. Photometry of the HostM is extracted in the elliptical aperture shown in Fig. 4 using the software A-PHOT (Merlin et al. 2019). The median(mean) magnification μtot within such an aperture is nearly constant, 19.1(19.3), with a standard deviation of 2.1. The SED-fitting on the multi-band photometry has been computed with the zphot.exe code (Fontana et al. 2000), as described in Castellano et al. (2016) (see also Vanzella et al. 2019) and the results are shown in Fig. 4, while the delensed physical quantities are reported in the last row of Table 1. Bruzual & Charlot (2003) templates have been adopted with the following priors: exponentially declining star-formation histories with e-folding times of 0.1 ≤ τ ≤ 15.0 Gyr, a Salpeter (1955) initial mass function, and the extinction laws from both Calzetti et al. (2000) and Small Magellanic Cloud (Prevot et al. 1984). We considered the following range of physical parameters: 0.0 ≤ E(B − V)≤1.1, Age > 20 Myr (defined as the onset of the star-formation episode) and metallicity Z/Z = 0.02, 0.2, 1.0. The resulting delensed ultraviolet magnitude at 1600 Å, stellar mass and SFR are 24.8, ∼1.1 × 109 M, ∼10 M yr−1, respectively, with an estimated age of ≃130 Myr and E(B − V) in the range 0.03 − 0.1 (uncertainties are reported in the same figure). An additional clear photometric signature is the apparent boost of flux in the F160W-band, indicated in Fig. 4. The zooming provided by strong lensing on arcs I and II coupled with VLT/X-shooter observations allow us to confirm that the origin of the photometric break is mainly due to knot 5.1 (i.e., the 3 Myr old YMC discussed above). Indeed, VLT/X-shooter spectroscopy shows the presence of nebular lines, Hβ and [O III]λλ4959, 5007 emerging prominently from 5.1 (see Appendix A). In particular, as a test case, we compared the F160W photometric excess measured on the YMC (5.1b + 5.1c, Fig. 4) and the corresponding line fluxes secured with VLT/X-shooter on the same pair. If interpreted as a boost from the nebular line group Hβ + [O III]λλ4959, 5007 (magnitude 20.25 on the F160W band and 21.10 at continuum), such photometric jump corresponds to (2.00 ± 0.05)×10−18 erg s−1 cm−2 Å−12, which is fully consistent with what is inferred from the spectrum of the same pair of images, (1.97 ± 0.10)×10−18 erg s−1 cm−2 Å−1 (such values also correspond to a rest-frame equivalent width of ≃1100 Å for the group of lines). We then compared what contribution the YMC provides to the F160W excess of the whole galaxy HostM. To this aim, we used A-PHOT to extract photometry on image 5.1 m within a circular aperture of 0.4″ diameter (≃2 × FWHM) and from the full image HostM (elliptical aperture, Fig. 4), adopting the F140W band as a reference probe of the underlying continuum. It turns out that the YMC contributes to about 30% of the observed excess on HostM. It is worth noting that its stellar mass accounts for only less than 1% of the total stellar mass of the galaxy.

thumbnail Fig. 4.

Nebular SED fitting of the Sunburst galaxy, referred to as HostM in the text is shown in the main left panel. It is moderately magnified and represents the emission from the entire galaxy. In the inset, the F555W-band image shows HostM with the elliptical aperture (red line) used to compute PSF-matched photometry and the knots 5.1, 5.2 and 5.3. The shaded green ellipses mark the photometric excess in the F160W band, due to nebular emission of [O III]λ4959 mainly due to the YMC 5.1. High S/N nebular SED-fit of YMC 5.1 is shown in the upper right panel, including both 5.1b and 5.1c to highlight the clear photometric jump in the F160W band, as confirmed from VLT/X-shooter observations (Appendix A). The inset shows the images of 5.1b,c and the adopted aperture in the F555W band (blue). Same SED-fit of the complex 5.3 h is shown in the bottom right panel as an example in which no evident nebular contribution is present in the F160W band. The inset shows the image of 5.3 h and the adopted aperture in the F555W band (blue ellipse).

Other regions of the same galaxy, like 5.3, do not show such an excess (Fig. 4). The overall SED of HostM includes the contribution of different star forming regions, with 5.1 showing the signature of a recent burst, whereas the rest of the galaxy might have more evolved characteristic features (e.g., relatively evolved star-forming regions with ages older than 3 Myr).

4.2. Morphology

The delensed image of the galaxy HostM is shown Fig. 5. This reconstruction is obtained as follows. We start with one of the multiple images of the Sunburst arc, namely: the arc segment shown in inset III of Fig. 2 (also displayed in Fig. 5). We super-sampled the arc on a regular grid whose pixel scale is 0.01 arcsec and we use the ray-tracing algorithm implemented in the software SkyLens (Meneghetti et al. 2008, 2010, 2017; Plazas et al. 2019) to map the surface brightness in each pixel onto the source plane. We chose to super-sample the image at a resolution of 0.01 arcsec/pixel in order to increase the number of support points on the source plane where we mapped the image surface brightness. The resulting unstructured data points are then interpolated on another regular grid, the resolution of which is again 0.01 arcsec. We convolved the image with a model of the HST PSF in the F555W band, obtained using the software Tiny Tim (Krist et al. 2011). Finally, we re-sampled the image at the resolution of 0.03 arcsec per pixel3.

thumbnail Fig. 5.

Reconstruction of the HostM galaxy in the source plane (delensed) is shown in the left panel, with almost all the details merged into a compact, barely resolved source. Middle panels: galaxy HostM in the F555W-band and the segmentation image with the blue contour at 5σ (left and right, respectively). The three regions 5.1, 5.2, and 5.3 (white apertures with 0.2″ diameter) containing star clusters account for at least 40% of the ultraviolet light of the galaxy. Right panel: same region of the galaxy is outlined with the blue contour in the F275W band probing LyC. It is clear the conspicuous escape of ionizing radiation from the YMC 5.1.

As expected, without strong lensing (i.e., assuming μ = 1), the HostM would appear as a galaxy sampled by ≲4 HST resolution elements (FWHM = 0.1″, in the F555W band, Fig. 2). The same figure also shows that the YMC (5.1) and the two 5.2 and 5.3 clusters merge into a nucleated unresolved component, which would be recognized as a ∼1 kpc-size star-forming clump (adopting the 0.1″ PSF in the UV, or 250 pc per pixel 0.03″). Future facilities assisted by extreme adaptive optics working at a resolution of a few tens of milliarcsec will be capable of discerning these structures also at the magnification of HostM (but see Sect. 7.2). In the framework of the LyC galaxies, the observation in unlensed fields with medium-resolution spectroscopy (R > 4000) of sources leaking LyC might capture ionized channels possibly associated to underlying (and totally unresolved) star clusters. In fact, in the case of Sunburst the detection of a prominent Lyα emission at the systemic redshift for the YMC (Rivera-Thorsen et al. 2017), in addition to a more regular and typical broader emission generated by radiative transfer effects, is consistent with the scenario in which an ionized channel was carved along the line of sight (corroborated by the LyC detection presented in Rivera-Thorsen et al. 2019, and as predicted by Behrens et al. 2014). In unlensed fields, such Lyα emission at systemic redshift has been identified for other LyC galaxies showing nucleated ultraviolet morphology; this, if compared to Sunburst (Vanzella et al. 2020a), might imply young massive star clusters are acting as prodigious ionizing photon producers and drillers within the interstellar medium. It is still not clear if this LyC escaping radiation mode is typical at redshift 2 − 3 (Matthee et al. 2021, and references therein) or even during reionization.

5. High cluster formation efficiency in the Sunburst galaxy

As discussed in Sect. 3.3, we identified at least 13 likely star clusters (see Table 1). Taking advantage of the new lens model (Pignataro et al. 2021), we can compare the integrated ultraviolet light of all star clusters to the luminosity of the galaxy HostM (known as TL(UV) parameter). We follow two ways of calculating such a fraction of light: (a) by comparing the delensed ultraviolet luminosity of the clusters and the host galaxy (dependent on the lens model) and (b) by comparing the ultraviolet light self consistently within the observed image HostM, independently from the lens model (this is the conservative approach). In case (a) the direct sum of delensed fluxes of the star clusters reported in Table 1 accounts for 55% of the UV light of HostM (M1600 = −20.3). The uncertainly on TL(UV) derived in this way is at least 50% of its value and is dominated by errors on the magnification factors (arcs I and II). In case (b) the fraction of the UV light is measured from the image HostM, which is nearly uniformly and moderately amplified (median μtot = 19). In this case, a relative measure of the UV light can be computed independently from magnification.

The ultraviolet flux of HostM shown in Fig. 5 is compared to the combined flux emerging from 5.1 m, 5.2 m and 5.3 m, calculated within circular apertures of diameter 0.2″ (≃2 × FWHM of the F555W band). The resulting TL(UV) from 5.1, 5.2, and 5.3 is at least ∼40% of the total 1600 Å ultraviolet light inferred from the host galaxy, while the single YMC 5.1 accounts for 16%. These values should be regarded as lower limits, since the other star clusters are not included in the calculation. Therefore the two inferred TL(UV) from case (a) and (b) are consistent, suggesting that TL(UV) spans the range between 40−60%. The fact that star clusters significantly contribute to the ultraviolet light might also suggest that a significant fraction of the star formation in the galaxy efficiently occurred in such objects.

Indeed, the cluster formation efficiency (Γ) is effectively a measurement of the total stellar mass forming in bound star clusters with respect to the total stellar mass forming in the galaxy, referring to the same interval of time (Bastian 2008). Following Adamo et al. (2015)Γ is defined as the star cluster formation rate (CFR) over the star formation rate of the host galaxy, both referred to the same interval of time, dt. The total stellar mass formed in clusters during dt is the CFR. In the Local Universe (≲20 Mpc distance), the CFR is typically estimated by integrating the stellar mass in bound clusters younger than 10 Myr (dt) above an established mass limit. The undetected part of the mass distribution is included by extrapolating over an assumed (and properly normalized) star cluster mass function down to 100 M (Adamo et al. 2017). Giving the complexity introduced by the lens we do not correct for such a mass limit, therefore, the inferred Γ can be considered a lower limit (such a limit does not alter the conclusion of this section). The age range used to derive the CFR is limited by the time scales over which the available SFR tracer are sensitive to. Finally, the resulting Γ is simply the ratio between the CFR and the SFR of the galaxy, expressed on a comparable time scale.

At high redshift such calculation is extremely challenging for two reasons: (I) gravitationally bound star clusters need to be identified and (II) the formation time scale for both clusters and the hosting galaxy need to be understood. The first point is partially addressed thanks to lensing amplification (Sect. 3.3), the second is the most uncertain since star formation rate indicators on short timescales (≲10 Myr) would be needed for each knot (e.g., Hα, see Sect. 7.2). The integrated stellar mass of bound star clusters reported in Table 1 is ≳70 × 106 M, under the aforementioned assumptions of instantaneous burst and age of 3 Myr for the knot 5.1 (Chisholm et al. 2019) and 7 Myr for the others (Appendix A). The time interval within which such a fraction of stellar mass formed cannot be constrained from the available data. However, it is worth noting that if we relax the ages of the clusters to 50 Myr (fixing the age of the knot 5.1 to 3 Myr), the integrated mass would exceed the total stellar mass of the galaxy derived from the SED fitting (≃109 M). Although this argument might appear stuck in a vicious circle, it might nonetheless suggest the sample of 13 clusters discussed in this work was likely formed in the recent past of the galaxy (< 50 Myr). Therefore, adopting a formation timescale of the star clusters of 10(20)(50) Myr, we have Γ10, 20, 50 > 60(30)(12)%. The lower limit is due to the fact that forming star clusters fainter than the detection limit are not considered. As mentioned above, the complex geometry of the lens currently prevents us from correcting for such an incompleteness and will require a second version of the lens model (e.g., by adding more multiple images or constraints from galaxies at different redshifts), so we rely on what is detected and consider it a lower limit (a future analysis will take into account this lack of completeness). Adopting a fiducial time scale of 20 Myr for the CFR and SFR ≃ 10 M yr−1 for the HostM, then the cluster formation efficiency Γ would exceed 30%. Both Γ and TL(UV) suggest there was a vigorous formation of star clusters in the Sunburst galaxy, in which a significant fraction of the stellar mass was produced in star clusters.

It is worth noting that such a high value of Γ is consistent with what is observed in the Local Universe. The observed area of HostM shown in Fig. 5 corresponds to ≃3 sq. kpc on the source plane, which implies a relatively high Log 10 ( Σ SFR ) 0 . 5 0.2 + 0.3 M $ _{10}(\Sigma_{\mathrm{SFR}}) \simeq 0.5_{-0.2}^{+0.3}\,M_\odot $ yr−1 kpc−2 (independent from magnification and with the 68% confidence interval on the SFR from SED fitting). The Sunburst galaxy falls in the upper right part of the Γ − ΣSFR relation, mimicking the local extreme forming galaxies described by Adamo et al. (2020b) (see Fig. 6). Interestingly, this is a region where similar local LyC leakers have also been identified (Keenan et al. 2017; Östlin et al. 2021).

thumbnail Fig. 6.

Cluster formation efficiency, Γ, as a function of the star-formation rate surface density ΣSFR. The compilation of gray symbols represent estimates in the local Universe taken from Adamo et al. (2020a,b). The solid blue line reproduces the Kruijssen (2012) fiducial model. The estimated lower limit on Γ (calculated over 20 Myr time scale, see text for details) for the Sunburst galaxy is reported with the filled red circle at 30%, while the transparent circle is the estimate fraction of the ultraviolet light (TL(UV)) computed from the delensed UV light integrated over all the star clusters reported in Table 1. The inferred quantities are likely lower limits, as additional clusters might have remained undetected.

6. Escaping Lyman continuum from the galaxy

Rivera-Thorsen et al. (2019) estimated the escape fraction, fesc, of the 12 multiple images of 5.1(a–n), assuming an intergalactic attenuation and Starburst99 models (Leitherer et al. 2014). Following the standard formalism for the estimation of the relative escape fraction (e.g., Steidel et al. 2001), fesc,rel, we have for 5.1:

fesc , rel (5.1) = F ( F 814 W ) / F ( 275 W ) intr F ( F 814 W ) / F ( 275 W ) obs × 1 T ( IGM ) , $$ \begin{aligned} \ \mathtt {fesc,rel(5.1)} = \frac{F(F814W)/F(275W)_{\rm intr}}{F(F814W)/F(275W)_{\rm obs}} \times \frac{1}{T(\mathrm{IGM})}, \end{aligned} $$(2)

where T(IGM) is the intergalactic attenuation of the LyC probed with the F275W band (T(IGM) = 1 meaning no attenuation) and the flux density ratios, intrinsic (“intr”) and the observed (“obs”), at the given positions a − n. The estimated escape fraction quantity has been observed to vary among the multiple images of 5.1, as it is likely modulated by the differential amount of H I gas along the different paths. Overall, Rivera-Thorsen et al. (2019) reported a line-of-sight fesc,rel = 93 11 + 7 $ 93^{+7}_{-11} $% with 46% as a robust lower limit.

The next step is to estimate the (observed) fesc,rel associated to the entire galaxy. We can adopt the same quantities presented above and include in Eq. (2) the observed flux density ratio referred to the whole galaxy observed in the segment, m, which has an observed magnitude: F814W = 21.60 (which corresponds to 24.8 delensed, Fig. 4). The measured F275W and F814W magnitudes of the knot 5.1 m are 27.58 (at S/N = 6.2) and 23.85 (at S/N > 40) (Table 1 of Rivera-Thorsen et al. 2019). The magnitude contrast in the F814W-band between the YMC (5.1 m, with observed F814W = 23.85) and the host galaxy (HostM, with observed F814W = 21.60) implies a reduction of fesc,rel by a factor 8. Therefore fesc,rel of the entire galaxy is ∼12%, with a lower limit of 6%. Such values should further be considered lower limits since there might be missing LyC signal emerging from other regions of the galaxy which are fainter than the current F275W depth. It is worth noting that the YMC has an intrinsic LyC and 1600 Å magnitudes of 30.6 and 26.9, respectively, after correcting for magnification at the predicted position of 5.1 m (μtot ≃ 16.1). Such values confirm the difficulty in securing relatively bright (sub-L) high redshift LyC galaxies in unlensed fields down to fesc,rel ∼ 10%. Those currently confirmed at high redshift are relatively bright LyC leakers and might be the result of hidden and multiple (spatially indistinguishable) star clusters (e.g., Steidel et al. 2018; Pahl et al. 2021; Vanzella et al. 2016, 2018, 2020a; de Barros et al. 2016). In the case of Sunburst galaxy, a single YMC dominates the observed LyC leakage. However, it is worth noting that in general multiple and coeval YMCs containing massive O-type stars might concur to produce an integrated large ionizing photon production efficiency, observable through prodigiously large photometric excess due to nebular lines in the optical rest-frame bands, much stronger than what is reported here (see Sect. 7.1). In this regard, it is worth referring to the recently discovered super-bright unlensed LyC galaxy by Marques-Chaves et al. (2021) showing QSO-like luminosity, MUV = −24.11, large sSFR ≃ 100 Gyr−1, and Log10SFR)∼2.2, implying a large stellar cluster formation efficiency (Γ) according to the relation shown in Fig. 6.

7. Summary and conclusions

The Sunburst galaxy is an exceptional laboratory for investigating high-z star formation for two main reasons: (1) individual forming stellar clusters can be identified at cosmological distance and (2) the link between the LyC leakage of the galaxy and the presence of such young (and bursty) star clusters can be seen. The main results of this study are reported in the following:

  1. Star clusters at cosmological distance. The Sunburst arcs are the result of a large lensing tangential amplification of a dwarf galaxy, in which star-forming regions down to a few parsec scale are probed. More than 50 multiple images of at least 13 star-forming knots belonging to the same galaxy have been identified (a panoramic view is shown in Fig. 2). In particular, taking advantage of a new accurate lens model (Pignataro et al. 2021), we find that such knots are likely gravitationally bound stellar clusters with stellar masses and effective radii spanning the range ≃105 − 7 M and ≃1 − 20 pc. The resulting cluster formation efficiency of the galaxy is Γ > 30%, placing it in the regime of extreme physical conditions when compared to local galaxies showing similar large Γ (Adamo et al. 2020b).

  2. The region emitting LyC is a young massive star cluster. The knot 5.1 is the youngest (3 Myr, Chisholm et al. 2019), the most massive (107 M both from photometric and dynamical estimates) and the brightest (MUV = −18.6) among the identified star clusters. Such a YMC shows an effective radius smaller than 10 pc and evident presence of hot massive (including Wolf-Rayet) stars (based on evident detection of P Cygni profiles of N Vλ1240, C IVλ1550 and broad He IIλ1640, Chisholm et al. 2019; Vanzella et al. 2020a); LyC emerges from this cluster with large relative escape fraction of 43 − 93%, depending on which multiple image is used for the calculation (Rivera-Thorsen et al. 2019). The delensed LyC(1600 Å) magnitude is 30.6(26.9). The sSFR is very large > 330 Gyr−1, consistently with the prominent optical nebular emissions lines, e.g., Hβ + [O III]λλ4959, 5007 rest-frame equivalent width of ≳1000 Å (as inferred from VLT/X-shooter spectroscopy) and as expected in such bursty events exhibiting ionizing photon production efficiency that is greater than the canonical values (Chevallard et al. 2018). The O32 index ([O III]λλ4959, 5007/[O II]λ3727, 3729) estimated on the YMC is 17 ± 3, consistent with the necessary condition of having a LyC leakage (Barrow et al. 2020, see also Izotov et al. 2018).

  3. The hosting galaxy. Sunburst would appear as a relatively small (≃3 sq. kpc) and low mass (≃109 M) galaxy with a dense star formation rate surface density (Log10SFR)≃0.5) and large sSFR ≳ 10 Gyr−1. It would be marginally resolved in the UV, nucleated and dominated by the YMC 5.1 SF knot UV emission, together with 5.2 and 5.3, which would be recognized as a kpc scale SF clump without lensing. The fesc,rel of the galaxy is > 6 − 12%. The lower limit is due to the fact that additional – possibly missing – LyC radiation might escape along the l.o.s. from other regions of the galaxy, but not probed with current depth yet. The photometric jump observed in the F160W-band (relative to the continuum probed in the F140W-band) implies a rest-frame equivalent width of the line complex Hβ + [O III]λλ4959, 5007 of ≃450 Å, of which ≃30% of the signal is due to the single YMC 5.1. It is worth noting that the same YMC 5.1 represents ≲1% of the total stellar mass of the galaxy.

7.1. Plausible indirect evidence of a high cluster formation efficiency at high redshift (and during reionization)

Studies in the local Universe have already demonstrated the key role of super star clusters as sources of intense stellar feedback that significantly affects the interstellar medium of the host galaxies (ionization, cavities, or outflows, e.g., James et al. 2016; Adamo et al. 2020b), eventually carving ionized optically thin channels to LyC photons (e.g., Micheva et al. 2017; Bik et al. 2018). The aforementioned properties of the Sunburst galaxy resemble those of confirmed unlensed LyC leakers at z ≲ 4 (e.g., Vanzella et al. 2016, 2018; de Barros et al. 2016; Schaerer et al. 2016; Izotov et al. 2016, 2021b; Marques-Chaves et al. 2021) and those of z ∼ 2 − 3 analogs of z > 6.5 sources contributing to cosmic reionization (Du et al. 2020; Tang et al. 2021; Katz et al. 2020; Matthee et al. 2021). For example, the photometric discontinuity observed with Spitzer/IRAC at z ∼ 6.5 − 8 due to optical nebular emission lines is also a distinctive signature present in the Sunburst SED for the same rest-optical lines (e.g., Castellano et al. 2017; Smit et al. 2014, 2015, 2016; Roberts-Borsani et al. 2016; Laporte et al. 2014; Finkelstein et al. 2013; De Barros et al. 2019). In particular, Endsley et al. (2021) find a large median (Hβ + [O III]) EW of 760 ± 100 Å, with a significant fraction (23 ± 7%) showing extreme values larger than 1200 Å. Similarly, Castellano et al. (2017) derive EW of 1500 ± 500 Å at z ≃ 7 from Spitzer/IRAC color excess, as well as De Barros et al. (2019) at z ≃ 8. Such results imply a high specific star formation rate (sSFR > 10−35 Gyr−1) and suggest intense bursts are common in these sources during reionization. Conversely, the occurrence of strong (Hβ + [O III]) emission is only a few percent at lower redshift, z ≃ 2 − 3 (e.g., Fig. 9 of Du et al. 2020, and references therein), suggesting that the frequency of such phases of intense star formation evolves with cosmic time, which appears to cycle regularly in the reionization era (Endsley et al. 2021). A possible explanation is that the net increase of the SFR surface density with increasing redshift (up to z ∼ 8, Naidu et al. 2020) would imply a high cluster formation efficiency (as observed for dense star-forming galaxies in the local Universe, Adamo et al. 2020b), namely: a prodigious “burstiness” due to actively forming star clusters. An increase of Γ with redshift has also been calculated by Pfeffer et al. (2018) in cosmological simulations that follow the co-evolution of galaxies and their star cluster populations, assuming the observed local Γ − ΣSFR relation.

The Sunburst galaxy described in this work represents the first example to demonstrate that prominent optical-rest emission lines, a high cluster formation efficiency Γ > 30 − 50%, relatively large sSFR ∼ 10 Gyr−1, and Lyman continuum leakage can co-exist together in the same system. From the observational point of view, such objects as Sunburst represent the “Rosetta stone” of stellar ionization, which might be rare at z ∼ 2, but possibly more frequent at z > 6.5 (e.g., Matthee et al. 2021). The presence of older (than 3 Myr) massive star clusters in the galaxy suggests that similar (bursty) conditions where in place in the past, potentially favoring an intermittent behavior of the escaping LyC radiation and ionizing production (Trebitsch et al. 2017; Wise et al. 2014). Sunburst is not as extreme as other LyC sources observed at z ≃ 3 or in the local Universe, for which (Hβ + [O III]) EW exceeds 1000−2000 Å rest-frame (e.g., at z = 3.2 Vanzella et al. 2016; de Barros et al. 2016 or in the Local Universe, e.g., Izotov et al. 2021b).

It is worth noting that young star clusters are the sites where: (1) the majority (> 90%) of ionizing massive stars (O-type) are forged Vargas-Salazar et al. (2020) and (2) stellar feedback consequently originated (e.g., Bik et al. 2018). Therefore, if star-formation drives reionization, then such aggregates are likely the key ionizing agents whose occurrence might be higher at high redshift when nebular emission is observed to be most prominent (Endsley et al. 2021; Castellano et al. 2017; De Barros et al. 2019).

It is worth considering the possibility that the presence of multiple and almost coeval massive young star clusters might concur to elevate the ionizing photon production, especially if the truncation mass of the initial cluster mass function increases with redshift (e.g., Pfeffer et al. 2018). Such conditions might also increase the fesc,rel, as multiple ionized channels can be carved in the ISM. Local examples rich of young star clusters and high cluster formation efficiency are the dwarf starburst galaxies NGC 5253 (Calzetti et al. 2015b), NGC 3125 (Wofford et al. 2014), or NGC 1569 (Anders et al. 2004). Finally, it is worth noting that local LyC leakers show a rather nucleated morphology – implying high SFR densities (Izotov et al. 2018), which in turn, follows from the findings of Adamo et al. (2020a) as well as the details in Fig. 6; this would suggest young star clusters are present (though not spatially recognizable) and possibly playing a key role in carving ionizing channels, as in the case of Sunburst.

7.2. Current limitations and future prospects

The presence of a bright star (H < 14 in Vega mag) only 3.3″ away from arc IV of Sunburst allowed us to use the VLT/MUSE narrow field mode (NFM) with optimal AO correction (see Fig. 7). Such an XAO correction provides PSF in line with the expectations, even better than 60 mas. Details of these observations are presented in Appendix C. There are two key aspects regarding MUSE-NFM worth mentioning: the one hand, we found the achieved angular resolution was optimal (FWHM ≲ 60 mas at λ > 6000 Å) despite the relatively high airmass (∼1.7), showing that MUSE-NFM equals HST imaging at similar wavelengths (see Fig. C.1). A resolution of 60 mas on the counter-arc corresponds to ∼60 parsec on the source plane along the tangential direction.

thumbnail Fig. 7.

Main panel: HST/F555W-band image (30 mas/pixel) with superimposed the field of views (FoV) of future instruments, such as VLT/MAVIS (IFU and imaging modes), ELT/MAORY-MICADO (imaging), and HARMONI (IFU), along with the layout of the MUSE-NFM targeting arc IV (transparent green box in the bottom-left corner, presented in Appendix C). The four arcs of Sunburst are indicated with I, II, III, IV. MUSE-NFM requires the presence of a very bright star within 3.4″ from the target, whilst future AO facilities will definitely probe stellar clusters down to the parsec scale in the most magnified regions. Right panels: zoomed region indicated with the orange shaded rectangle located in the main panel over the objects 5.2 h and 5.3 h (left). In particular, the region enclosing the knots 5.2 h and 5.3 h will be covered with MAVIS and HARMONI IFUs which will probe 1−5 parsec per spaxel, while imaging with MAVIS and MAORY-MICADO (wide blue and magenta squares in the main panel) will probe the parsec scale in the optical and near-infrared, respectively.

On the other hand, the very limited sky coverage offered by MUSE-NFM (e.g., the lack of a very close bright star) prevents us from targeting very magnified regions of the Sunburst arcs, such arc II, in which a 60 mas resolution would have secured spectroscopy at 5−10 parsec scale. MUSE-WFM on that arc, assuming seeing 0.8″ and μtang = 50 − 100 probes 131 − 65 pc scale.

For this reason, an increased sky coverage in which XAO can be performed on, e.g. arc II, is key in such studies, allowing us to reach a spatial resolution of a few pc. In general, for super-lensed targets like Sunburst (e.g., Sharon et al. 2020; Rigby et al. 2018a,b), VLT/MAVIS and ELT MAROY/MICADO or HARMONI are complementary (in the wavelength domain) and will provide a census of the cluster formation efficiency across cosmic epochs, their role in shaping the galaxy growth and ionizing output. In this respect, ELT HARMONI will be ideal for probing the Hα of each SF knot in the near infrared wavelengths, while VLT/MAVIS will complement on similar angular resolution in the ultraviolet and optical bands. The JWST integral field spectroscopy will also map Hα line emission, which will be free from sky lines and background.


2

Calculated as (−2)×10−0.4(m + 48.59), where c is the speed of light, λ the effective wavelength of the F160W-band (15405.2 Å) and m the AB magnitude.

3

The observed image of the arc used to reconstruct the galaxy HostM is already convolved with the HST PSF. When re-convolving the unlensed image, our procedure does not account for this pre-existing convolution. By de-lensing the arc, its size, including that of the PSF, is reduced by a factor μtot ≃ 19. Thus, the resulting PSF size on the source plane is so small that we can safely neglect it.

Acknowledgments

We thank the anonymous referee for the careful reading and constructive comments. We thank A. Adamo for very stimulating discussion and providing the measures of Γ in the local Universe (shown in Fig. 6). E.V. thanks D. Calzetti for the precious clarifications on the formation channels of O-type stars. E.V. thanks A. Renzini for illuminating discussions about any aspect of star, star-cluster and galaxy formation, and J. Chisholm for the stimulating interactions about LyC galaxies. We thank A. Fontana for making available the zphot.exe code. This project is partially funded by PRIN-MIUR 2017WSCC32 “Zooming into dark matter and proto-galaxies with massive lensing clusters”. We acknowledge funding from the INAF for “interventi aggiuntivi a sostegno della ricerca di main-stream” (1.05.01.86.31). P.B. acknowledges financial support from ASI though the agreement ASI-INAF n. 2018-29-HH.0. M.M. acknowledges support from the Italian Space Agency (ASI) through contract “Euclid – Phase D” and from the grant MIUR PRIN 2015 “Cosmology and Fundamental Physics: illuminating the Dark Universe with Euclid”. F.C. acknowledges support from grant PRIN MIUR 2017 – 20173ML3WW_s. C.G. acknowledges support by VILLUM FONDEN Young Investigator Programme through grant no. 10123. K.C. acknowledges funding from the ERC through the award of the Consolidator Grant ID 681627-BUILDUP. G.B.C. acknowledges the Max Planck Society for financial support through the Max Planck Research Group for S.H. Suyu and the academic support from the German Centre for Cosmological Lensing. This research made use of the following open-source packages for Python and we are thankful to the developers of these: Matplotlib (Hunter 2007), MPDAF (Piqueras et al. 2019), PyMUSE (Pessa et al. 2020), Numpy (van der Walt et al. 2011).

References

  1. Adamo, A., Kruijssen, J. M. D., Bastian, N., Silva-Villa, E., & Ryon, J. 2015, MNRAS, 452, 246 [NASA ADS] [CrossRef] [Google Scholar]
  2. Adamo, A., Ryon, J. E., Messa, M., et al. 2017, ApJ, 841, 131 [Google Scholar]
  3. Adamo, A., Zeidler, P., Kruijssen, J. M. D., et al. 2020a, Space Sci. Rev., 216, 69 [Google Scholar]
  4. Adamo, A., Hollyhead, K., Messa, M., et al. 2020b, MNRAS, 499, 3267 [Google Scholar]
  5. Anders, P., de Grijs, R., Fritze-v Alvensleben, U., & Bissantz, N. 2004, MNRAS, 347, 17 [NASA ADS] [CrossRef] [Google Scholar]
  6. Bacon, R., Accardo, M., Adjali, L., et al. 2012, The Messenger, 147, 4 [NASA ADS] [Google Scholar]
  7. Barrow, K. S. S., Robertson, B. E., Ellis, R. S., et al. 2020, ApJ, 902, L39 [NASA ADS] [CrossRef] [Google Scholar]
  8. Bastian, N. 2008, MNRAS, 390, 759 [Google Scholar]
  9. Bastian, N., Saglia, R. P., Goudfrooij, P., et al. 2006, A&A, 448, 881 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  10. Behrens, C., Dijkstra, M., & Niemeyer, J. C. 2014, A&A, 563, A77 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  11. Bergamini, P., Rosati, P., Vanzella, E., et al. 2021, A&A, 645, A140 [EDP Sciences] [Google Scholar]
  12. Bik, A., Östlin, G., Menacho, V., et al. 2018, A&A, 619, A131 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  13. Bournaud, F. 2016, in Bulge Growth Through Disc Instabilities in High-Redshift Galaxies, eds. E. Laurikainen, R. Peletier, & D. Gadotti, 418, 355 [NASA ADS] [Google Scholar]
  14. Bruzual, G., & Charlot, S. 2003, MNRAS, 344, 1000 [NASA ADS] [CrossRef] [Google Scholar]
  15. Cabrera-Ziri, I., Bastian, N., Longmore, S. N., et al. 2015, MNRAS, 448, 2224 [NASA ADS] [CrossRef] [Google Scholar]
  16. Cabrera-Ziri, I., Bastian, N., Hilker, M., et al. 2016, MNRAS, 457, 809 [NASA ADS] [CrossRef] [Google Scholar]
  17. Calzetti, D., Armus, L., Bohlin, R. C., et al. 2000, ApJ, 533, 682 [NASA ADS] [CrossRef] [Google Scholar]
  18. Calzetti, D., Lee, J. C., Sabbi, E., et al. 2015a, AJ, 149, 51 [Google Scholar]
  19. Calzetti, D., Johnson, K. E., Adamo, A., et al. 2015b, ApJ, 811, 75 [CrossRef] [Google Scholar]
  20. Caminha, G. B., Grillo, C., Rosati, P., et al. 2017, A&A, 600, A90 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  21. Carnall, A. C., Leja, J., Johnson, B. D., et al. 2019, ApJ, 873, 44 [Google Scholar]
  22. Castellano, M., Amorín, R., Merlin, E., et al. 2016, A&A, 590, A31 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  23. Castellano, M., Pentericci, L., Fontana, A., et al. 2017, ApJ, 839, 73 [NASA ADS] [CrossRef] [Google Scholar]
  24. Cava, A., Schaerer, D., Richard, J., et al. 2018, Nat. Astron., 2, 76 [Google Scholar]
  25. Chevallard, J., Charlot, S., Senchyna, P., et al. 2018, MNRAS, 479, 3264 [Google Scholar]
  26. Chisholm, J., Rigby, J. R., Bayliss, M., et al. 2019, ApJ, 882, 182 [Google Scholar]
  27. Dahle, H., Aghanim, N., Guennou, L., et al. 2016, A&A, 590, L4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  28. de Barros, S., Vanzella, E., Amorín, R., et al. 2016, A&A, 585, A51 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  29. De Barros, S., Oesch, P. A., Labbé, I., et al. 2019, MNRAS, 489, 2355 [NASA ADS] [CrossRef] [Google Scholar]
  30. Dessauges-Zavadsky, M., Schaerer, D., Cava, A., Mayer, L., & Tamburello, V. 2017, ApJ, 836, L22 [Google Scholar]
  31. Dessauges-Zavadsky, M., Richard, J., Combes, F., et al. 2019, Nat. Astron., 3, 1115 [NASA ADS] [CrossRef] [Google Scholar]
  32. Du, X., Shapley, A. E., Tang, M., et al. 2020, ApJ, 890, 65 [NASA ADS] [CrossRef] [Google Scholar]
  33. Elmegreen, D. M., Elmegreen, B. G., Ravindranath, S., & Coe, D. A. 2007, ApJ, 658, 763 [Google Scholar]
  34. Elmegreen, D. M., Elmegreen, B. G., Whitmore, B. C., et al. 2021, ApJ, 908, 121 [NASA ADS] [CrossRef] [Google Scholar]
  35. Endsley, R., Stark, D. P., Chevallard, J., & Charlot, S. 2021, MNRAS, 500, 5229 [Google Scholar]
  36. Faure, B., Bournaud, F., Fensch, J., et al. 2021, MNRAS, 502, 4641 [CrossRef] [Google Scholar]
  37. Finkelstein, S. L., Papovich, C., Dickinson, M., et al. 2013, Nature, 502, 524 [CrossRef] [Google Scholar]
  38. Fontana, A., D’Odorico, S., Poli, F., et al. 2000, AJ, 120, 2206 [NASA ADS] [CrossRef] [Google Scholar]
  39. Förster Schreiber, N. M., Shapley, A. E., Genzel, R., et al. 2011, ApJ, 739, 45 [Google Scholar]
  40. Gaia Collaboration (Prusti, T., et al.) 2016, A&A, 595, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  41. Gaia Collaboration (Brown, A. G. A., et al.) 2018, A&A, 616, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  42. Georgiev, I. Y., Hilker, M., Puzia, T. H., Goudfrooij, P., & Baumgardt, H. 2009, MNRAS, 396, 1075 [NASA ADS] [CrossRef] [Google Scholar]
  43. Gieles, M., & Portegies Zwart, S. F. 2011, MNRAS, 410, L6 [Google Scholar]
  44. Gonzaga, S., Hack, W., Fruchter, A., & Mack, J. 2012, The DrizzlePac Handbook (Baltimore: STScI) [Google Scholar]
  45. Guo, Y., Rafelski, M., Bell, E. F., et al. 2018, ApJ, 853, 108 [NASA ADS] [CrossRef] [Google Scholar]
  46. Heckman, T. M., Borthakur, S., Overzier, R., et al. 2011, ApJ, 730, 5 [Google Scholar]
  47. Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [NASA ADS] [CrossRef] [Google Scholar]
  48. Izotov, Y. I., Orlitová, I., Schaerer, D., et al. 2016, Nature, 529, 178 [Google Scholar]
  49. Izotov, Y. I., Worseck, G., Schaerer, D., et al. 2018, MNRAS, 478, 4851 [Google Scholar]
  50. Izotov, Y. I., Guseva, N. G., Fricke, K. J., et al. 2021a, A&A, 646, A138 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  51. Izotov, Y. I., Worseck, G., Schaerer, D., et al. 2021b, MNRAS, 503, 1734 [NASA ADS] [CrossRef] [Google Scholar]
  52. James, B. L., Auger, M., Aloisi, A., Calzetti, D., & Kewley, L. 2016, ApJ, 816, 40 [Google Scholar]
  53. Jaskot, A. E., Dowd, T., Oey, M. S., Scarlata, C., & McKinney, J. 2019, ApJ, 885, 96 [Google Scholar]
  54. Johnson, T. L., Rigby, J. R., Sharon, K., et al. 2017, ApJ, 843, L21 [Google Scholar]
  55. Katz, H., Ďurovčíková, D., Kimm, T., et al. 2020, MNRAS, 498, 164 [Google Scholar]
  56. Keenan, R. P., Oey, M. S., Jaskot, A. E., & James, B. L. 2017, ApJ, 848, 12 [NASA ADS] [CrossRef] [Google Scholar]
  57. Krist, J. E., Hook, R. N., & Stoehr, F. 2011, in Optical Modeling and Performance Predictions V, ed. M. A. Kahan, SPIE Conf. Ser., 8127, 81270J [Google Scholar]
  58. Kruijssen, J. M. D. 2012, MNRAS, 426, 3008 [Google Scholar]
  59. Laporte, N., Streblyanska, A., Clement, B., et al. 2014, A&A, 562, L8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  60. Leitherer, C., Ekström, S., Meynet, G., et al. 2014, ApJS, 212, 14 [Google Scholar]
  61. Livermore, R. C., Jones, T. A., Richard, J., et al. 2015, MNRAS, 450, 1812 [Google Scholar]
  62. Marques-Chaves, R., Schaerer, D., Álvarez-Márquez, J., et al. 2021, MNRAS, 507, 524 [NASA ADS] [CrossRef] [Google Scholar]
  63. Matthee, J., Sobral, D., Hayes, M., et al. 2021, MNRAS, 505, 1382 [NASA ADS] [CrossRef] [Google Scholar]
  64. McLaughlin, D. E., & van der Marel, R. P. 2005, ApJS, 161, 304 [NASA ADS] [CrossRef] [Google Scholar]
  65. Meneghetti, M., Melchior, P., Grazian, A., et al. 2008, A&A, 482, 403 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  66. Meneghetti, M., Rasia, E., Merten, J., et al. 2010, A&A, 514, A93 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  67. Meneghetti, M., Natarajan, P., Coe, D., et al. 2017, MNRAS, 472, 3177 [Google Scholar]
  68. Merlin, E., Pilo, S., Fontana, A., et al. 2019, A&A, 622, A169 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  69. Micheva, G., Oey, M. S., Jaskot, A. E., & James, B. L. 2017, ApJ, 845, 165 [NASA ADS] [CrossRef] [Google Scholar]
  70. Naidu, R. P., Tacchella, S., Mason, C. A., et al. 2020, ApJ, 892, 109 [NASA ADS] [CrossRef] [Google Scholar]
  71. Östlin, G., Cumming, R. J., & Bergvall, N. 2007, A&A, 461, 471 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  72. Östlin, G., Rivera-Thorsen, T. E., Menacho, V., et al. 2021, ApJ, 912, 155 [CrossRef] [Google Scholar]
  73. Pahl, A. J., Shapley, A., Steidel, C. C., Chen, Y., & Reddy, N. A. 2021, MNRAS, 505, 2447 [NASA ADS] [CrossRef] [Google Scholar]
  74. Perna, M., Brusa, M., Cresci, G., et al. 2015, A&A, 574, A82 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  75. Pessa, I., Tejos, N., & Moya, C. 2020, in Astronomical Data Analysis Software and Systems XXVII, eds. P. Ballester, J. Ibsen, M. Solar, & K. Shortridge, ASP Conf. Ser., 522, 61 [NASA ADS] [Google Scholar]
  76. Pfeffer, J., Kruijssen, J. M. D., Crain, R. A., & Bastian, N. 2018, MNRAS, 475, 4309 [Google Scholar]
  77. Pignataro, G. V., Bergamini, P., Meneghetti, M., et al. 2021, A&A, 655, A81 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  78. Piqueras, L., Conseil, S., Shepherd, M., et al. 2019, in Astronomical Data Analysis Software and Systems XXVI, eds. M. Molinaro, K. Shortridge, & F. Pasian, ASP Conf. Ser., 521, 545 [Google Scholar]
  79. Plazas, A. A., Meneghetti, M., Maturi, M., & Rhodes, J. 2019, MNRAS, 482, 2823 [Google Scholar]
  80. Prevot, M. L., Lequeux, J., Maurice, E., Prevot, L., & Rocca-Volmerange, B. 1984, A&A, 132, 389 [Google Scholar]
  81. Reddy, N. A., Steidel, C. C., Pettini, M., & Bogosavljević, M. 2016, ApJ, 828, 107 [NASA ADS] [CrossRef] [Google Scholar]
  82. Rhoads, J. E., Malhotra, S., Richardson, M. L. A., et al. 2014, ApJ, 780, 20 [Google Scholar]
  83. Rigby, J. R., Johnson, T. L., Sharon, K., et al. 2017, ApJ, 843, 79 [Google Scholar]
  84. Rigby, J. R., Bayliss, M. B., Chisholm, J., et al. 2018a, ApJ, 853, 87 [Google Scholar]
  85. Rigby, J. R., Bayliss, M. B., Sharon, K., et al. 2018b, AJ, 155, 104 [Google Scholar]
  86. Rivera-Thorsen, T. E., Dahle, H., Gronke, M., et al. 2017, A&A, 608, L4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  87. Rivera-Thorsen, T. E., Dahle, H., Chisholm, J., et al. 2019, Science, 366, 738 [Google Scholar]
  88. Roberts-Borsani, G. W., Bouwens, R. J., Oesch, P. A., et al. 2016, ApJ, 823, 143 [NASA ADS] [CrossRef] [Google Scholar]
  89. Romano, M., Grazian, A., Giallongo, E., et al. 2019, A&A, 632, A45 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  90. Ryon, J. E., Gallagher, J. S., Smith, L. J., et al. 2017, ApJ, 841, 92 [Google Scholar]
  91. Salpeter, E. E. 1955, ApJ, 121, 161 [Google Scholar]
  92. Schaerer, D., Izotov, Y. I., Verhamme, A., et al. 2016, A&A, 591, L8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  93. Sharon, K., Bayliss, M. B., Dahle, H., et al. 2020, ApJS, 247, 12 [CrossRef] [Google Scholar]
  94. Smit, R., Bouwens, R. J., Labbé, I., et al. 2014, ApJ, 784, 58 [NASA ADS] [CrossRef] [Google Scholar]
  95. Smit, R., Bouwens, R. J., Franx, M., et al. 2015, ApJ, 801, 122 [Google Scholar]
  96. Smit, R., Bouwens, R. J., Labbé, I., et al. 2016, ApJ, 833, 254 [NASA ADS] [CrossRef] [Google Scholar]
  97. Steidel, C. C., Pettini, M., & Adelberger, K. L. 2001, ApJ, 546, 665 [Google Scholar]
  98. Steidel, C. C., Bogosavljević, M., Shapley, A. E., et al. 2018, ApJ, 869, 123 [Google Scholar]
  99. Stuik, R., Bacon, R., Conzelmann, R., et al. 2006, New Astron. Rev., 49, 618 [CrossRef] [Google Scholar]
  100. Tang, M., Stark, D. P., Chevallard, J., et al. 2021, MNRAS, 503, 4105 [NASA ADS] [CrossRef] [Google Scholar]
  101. Trebitsch, M., Blaizot, J., Rosdahl, J., Devriendt, J., & Slyz, A. 2017, MNRAS, 470, 224 [Google Scholar]
  102. van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Comput. Sci. Eng., 13, 22 [Google Scholar]
  103. Vanzella, E., de Barros, S., Vasei, K., et al. 2016, ApJ, 825, 41 [NASA ADS] [CrossRef] [Google Scholar]
  104. Vanzella, E., Calura, F., Meneghetti, M., et al. 2017, MNRAS, 467, 4304 [Google Scholar]
  105. Vanzella, E., Nonino, M., Cupani, G., et al. 2018, MNRAS, 476, L15 [Google Scholar]
  106. Vanzella, E., Calura, F., Meneghetti, M., et al. 2019, MNRAS, 483, 3618 [Google Scholar]
  107. Vanzella, E., Caminha, G. B., Calura, F., et al. 2020a, MNRAS, 491, 1093 [Google Scholar]
  108. Vanzella, E., Meneghetti, M., Caminha, G. B., et al. 2020b, MNRAS, 494, L81 [Google Scholar]
  109. Vanzella, E., Meneghetti, M., Pastorello, A., et al. 2020c, MNRAS, 499, L67 [CrossRef] [Google Scholar]
  110. Vanzella, E., Caminha, G. B., Rosati, P., et al. 2021, A&A, 646, A57 [EDP Sciences] [Google Scholar]
  111. Vargas-Salazar, I., Oey, M. S., Barnes, J. R., et al. 2020, ApJ, 903, 42 [NASA ADS] [CrossRef] [Google Scholar]
  112. Wang, B., Heckman, T. M., Amorín, R., et al. 2021, ApJ, 916, 3 [NASA ADS] [CrossRef] [Google Scholar]
  113. Williams, C. C., Spilker, J. S., Whitaker, K. E., et al. 2021, ApJ, 908, 54 [NASA ADS] [CrossRef] [Google Scholar]
  114. Wise, J. H., Demchenko, V. G., Halicek, M. T., et al. 2014, MNRAS, 442, 2560 [NASA ADS] [CrossRef] [Google Scholar]
  115. Wofford, A., Leitherer, C., Chandar, R., & Bouret, J.-C. 2014, ApJ, 781, 122 [Google Scholar]
  116. Worseck, G., Prochaska, J. X., O’Meara, J. M., et al. 2014, MNRAS, 445, 1745 [NASA ADS] [CrossRef] [Google Scholar]
  117. Wuyts, S., Förster Schreiber, N. M., Nelson, E. J., et al. 2013, ApJ, 779, 135 [Google Scholar]
  118. Zanella, A., Daddi, E., Le Floc’h, E., et al. 2015, Nature, 521, 54 [Google Scholar]
  119. Zanella, A., Le Floc’h, E., Harrison, C. M., et al. 2019, MNRAS, 489, 2792 [Google Scholar]

Appendix A: Hα and [O III]λ5007 emissions from star-forming knots

X-Shooter observations have been performed during April-August 2019 (Prog. 0103.A-0688, PI Vanzella). The two arcs I and II have been targeted with two slits oriented along the tangential stretch. The data reduction is described in Vanzella et al. (2020c). Briefly, no dithering has been used since the arcs develop along the spatial direction probed with the slits. Suitable free windows without bright objects have been included such that sky was modeled with polynomial fitting and subtracted. Here we focus on a specific region of arc I where two mirrored regions containing 10 SF knots each have been observed (see Figure A.1). In particular knot 5.1c, b and 5.1d, e, f, as well as the groups “c” and “d” of blue knots also shown in Figure 2 (5.4, 5.5, 5.6, 5.7, 5.3, 5.9, 5.2, and 5.11). Figure A.1 shows the spectrum at [O III]λ5007 and Hα wavelengths. Both lines emerge prominently from the YMC 5.1, while for the groups we can set only upper limits on their line fluxes and equivalent widths. As a test case, we verified the consistency between the [O III]λ5007 line flux inferred from VLT/X-Shooter and the photometric jump observed in the SED fitting performed on the same aperture probing 5.1b, c (as indicated in Figure 4). The 0.6 magnitude boost in the H band is fully consistent with the inferred [O III]λ5007 equivalent width, also including the [O III]λ4959 and Hβ. What is inferred from the spectrum is slightly fainter than what observed with Hubble photometry because of slit losses.

thumbnail Fig. A.1.

Main panel (thick line box): the [O III]λ5007 (left) and the Hα (right) two-dimensional zoomed spectra extracted from X-Shooter observations at spectral resolution R ≃5600. The slit orientation is indicated with the yellow box in the inset showing the HST color image. The magenta dotted lines mark the position of the knots 5.1b, c and 5.1d, e, f over the spatial direction along the slit. The blue dotted lines show the expected positions of the region collecting the group of fainter knots, named “group” in the figure (see also Figure 2). Hα and [O III]λ5007 clearly emerge from knots 5.1b, c, d, e, f, conversely the same lines are not detected for the two mirrored “groups” at the given depth. The inferred rest-frame equivalent widths are also quoted. The transparent vertical blue stripes mark the position of night sky emission lines. The right panel shows the Hα equivalent width as a function of time from a Starburst99 model of an instantaneous burst and Z = 0.008Z. Magenta and blue dotted arrows indicate the case of the knot 5.1 and “Group”, respectively. Knot 5.1 is younger than 5 Myr (in agreement with the 3 Myr old age derived by Chisholm et al. 2019), while the complex “group” is likely older than 7 Myr.

Similarly, the Hα emission is also prominent from knot 5.1. It is worth noting that the multiple images within the slit allow us to formally quintuple the integration time on 5.1, and duplicate it over the groups “c” and “d”. The resulting equivalent width of Hα is higher than 350Å rest-frame for knot 5.1 and lower than 100Å for the group of knots. Such values under the assumption of an instantaneous burst of star formation imply ages younger than 5 Myr for the YMC knot 5.1 and older than ≃7 Myr on average for the group of knots (see rightmost side of Figure A). Such a lower limit on the age for the group of knots is the value assumed in the text and reported in Table 1. The case of YMC 5.1 is fully consistent with the age inferred from Chisholm et al. (2019) of 2.9 Myr (from Starburst99 modeling).

Appendix B: Dynamical mass of the LyC star cluster

The [O III]λ5007 line is detected at S/N> 100, it is free from significant sky emission and shows an evident asymmetric profile with a blue tail (see Figure B.1). A fit with a double Gaussian component has been performed and reproduces well the observed profile, suggesting there is outflowing ionized gas along the line of sight (broad component) and the emission associated to the H II region close to the stellar component and tracing the systemic redshift (core of the line). Such an outflow is consistent with what was already observed by Rivera-Thorsen et al. (2017) based on absorption lines of Silicon IV. Here we focus on the core of the line. The resulting velocity dispersion of such an emission from the fit is σv = 37 ± 5 km s−1 (the resolution element in the NIR arm is  ≃ 53 km s−1, or σv ≃ 22 km s−1). Adopting the inferred σv and the known effective radius Reff ≃ 8 pc of the YMC (5.1), the dynamical mass can be calculated following Rhoads et al. (2014) as Mdyn = 4σ2Reff/G (see also Vanzella et al. 2017), which turns out to be (1.0 ± 0.2)×107 M. Therefore, both the photometric and dynamical masses are compatible, with a ratio close to 1. As discussed in the text, the presence of the magnification factors in the numerator and denominator implies the ratio between stellar and virial masses weakly depends on lensing amplification.

thumbnail Fig. B.1.

[O III]λ5007 line emerging from the LyC knot 5.1 (gray band) and the best-fit solution made of the sum of two Gaussian components (red line), a broad one tracing the outflowing ionized gas (blue dotted) and the other tracing the core of the line; the latter arises from the H II region surrounding the star cluster and tracing the redshift of the stellar component (green dotted). On top axis, the velocity is referred to the peak position of the line. The clear blue tail extends up to −400 km s−1.

Appendix C: MUSE Narrow Field Mode spectroscopic observations

It is worth reporting here on the exceptional MUSE observation in the Narrow Field Mode (PI Vanzella) on arc IV, which provides spectroscopy and imaging at resolution smaller than 80 mas and better than HST blue bands. The proximity of a sufficiently bright star (H < 14 Vega, also dubbed as natural guide star, NGS) within 3.35″ angular separation from the targeted counterarc (arc IV) allowed for dedicated nearly diffraction-limited observations with MUSE (Bacon et al. 2012) Narrow Field Mode, offered by the GALACSI module implementing the laser tomography AO correction (Stuik et al. 2006). In this configuration, MUSE provides a field of view of 7.5″ × 7.5″ sampled at 0.025″/pixel scale. The relatively high (best) airmass of ∼1.7 for the target is typically not ideal for applying extreme AO corrections, however, we requested proper weather conditions as suggested in the manual for such airmass (seeing < 0.6″ at zenith). A total of 2.6h divided into two observing blocks (OB) of 3532s each on science (each OB split into 883 × 4) have been performed during April 29 2019 with the exceptional average seeing conditions of 0.4″ and producing a Strehl on the NGS of ≃20%. Data reduction have been performed following the prescription described in (Caminha et al. 2017; Vanzella et al. 2021).

Figure C.1 shows the FWHM of the Moffat profile fitted on the NGS (on-axis) star as a function of the wavelength. Remarkably, the FWHM decreases at values lower than 60 mas at λ > 6000Å, approaching 50 − 55 mas over the range 7000 − 9250Å. The same figure also shows the comparison on arc IV, between HST/ACS F606W and MUSE collapsed image in the range 6500 − 8000Å and 8000 − 9300Å. Even though the off-axis (∼3″ apart) FWHM is probably worse than the on-axis one on the NGS, MUSE with its FWHM < 100 mas appears to outperform HST F606W (∼100 mas) resolutions and allow us to probe smaller physical scales, at least in this specific arc.

thumbnail Fig. C.1.

Counter arc IV observed with HST ACS/F606W (left) and VLT/MUSE NFM after collapsing the datacube along the indicated wavelength intervals (middle and right images). The ground-based AO-assisted observations of arc IV outperform those from space (see narrow-band imaging after collapsing the cube in the wavelength ranges 6500-8000Å and 8000-9000Å), providing further evidences for the compactness of the star-forming region and in particular the LyC star cluster 5.1 (in this case, the image 5.1n has been targeted). The right most figure shows the FWHM of the OGS (On-axis Guide Source), measured within a 1.5″ region, as a function of the wavelength. The black line show the measurements in the final stacked data-cube, whereas the blue region indicates the scatter (68% confidence level) of the measurements in each single exposure. The FWHM is < 60 mas at λ > 7000Å and we assume a similar resolution is achieved on arc IV (middle).

Despite arc IV is not the most magnified, an upper limit on the size of the star-forming knot 5.1 of 40 pc (adopting FWHM < 60 mas) represents a remarkable technical achievement, considering the relatively high airmass (1.7). If, on the one hand, this MUSE-NFM shows its exceptional technical capabilities, which matches the HST imaging; on the other hand the poor flexibility (in terms of sky coverage) due to the need for a very close (within 3.35″) and bright (H < 14 Vega magnitude) natural guide star makes the observations of the most magnified arcs I or II inaccessible. Extreme AO-correction with much wider sky coverage is therefore essential for the future VLT instrumentation, such as MAVIS (7 mas/pixel, FWHM ∼25 mas in the optical bands), which will provide spatial details down to 1 − 3 pc/pix on arcs I or II (see Sect. 7.2).

All Tables

Table 1.

Summary of the identified SF knots observed on Sunburst.

All Figures

thumbnail Fig. 1.

Panoramic view of the main arcs discussed in this work. Left panel (in F555W-band image): Arcs are marked with red contours to guide the eye. The arrows indicate 11 multiple images of the star-forming knots labeled 5.1a–m (see text and Fig. 2 for additional details). The magnification map (total magnification, μtot) is reported on the right panel (from Pignataro et al. 2021) coded according with the color bar in square root scale: it saturates at μtot ≃ 10 (black) and 100 (white). The black contours enclose the region where μtot > 500, marking the locus of the critical line. The arcs contain multiple images of the same portions of the Sunburst galaxy and are subjected to amplification spanning the range 15 − 100, with the critical line crossing the arcs several times.

In the text
thumbnail Fig. 2.

Identification of the multiple images of compact star-forming knots is shown. In the bottom-left of the figure, the negative gray-scaled HST ACS/F555W image shows the panoramic view of the four arcs, labeled as I, II, III, and IV, with the little red square in the middle indicating the source region subject to large amplification. The SF knots are shown in the blue border insets labeled as I, II, III, and IV, following the arcs notation in the bottom left inset and with increasing magnification from bottom to top (long arrow on the right). Each inset includes the negative gray scale F555W image and a color image obtained by combining the WFC3/UVIS F275W (blue), WFC3/UVIS F606W (cyan), ACS/WFC F814W (yellow), and WFC3/IR F160W (red) filters (ESA/Hubble, NASA, Rivera-Thorsen et al., CC BY 4.0). Star-forming regions are indicated following the nomenclature of Pignataro et al. (2021). White dashed lines mark the position of the critical lines. The main triplet, 5.1, 5.2, and 5.3 are indicated with blue, red, and magenta open circles, respectively. The objects labeled with green IDs are the knots reported in Table 1. Those marked in yellow show all the identified multiple images of the knots. The inset in the bottom right (red contour) shows the reconstructed image of the Sunburst galaxy on the source plane, where the triplet 5.1, 5.2, and 5.3 would merge into a single HST pixel (∼30 mas separation). The observed angular separation between the knots 5.1 and 5.3 is reported on the right side of each inset, and the typical physical scales subtended by the FWHM (0.1″) along tangential direction is shown over the long black arrow to the right (arc II includes star-forming region smaller 10 pc). The transient (“Tr”) discussed in Vanzella et al. (2020c) and corresponding to the knot 5.10 is indicated with a black arrow in the inset II.

In the text
thumbnail Fig. 3.

Galfit fitting of the star-forming knots reported in Table 1 zooming on various regions (groups A, D, G, H and L, see Fig. 2). Left column: original F555W image, in which only the knots listed in Table 1 are marked with their IDs. Middle column: Galfit models and right column: residual after subtracting the models, with the positions of the knots indicated as in the left column. Besides the IDs in the reference sample (marked in blue), additional multiply imaged knots have been modeled, such as the cloud of clusters of the group D (see the same group described in Fig. 2). All knots are not resolved with resulting effective radii ≲1 pixel, independently from the Sérsic index (n = 0.5 − 4).

In the text
thumbnail Fig. 4.

Nebular SED fitting of the Sunburst galaxy, referred to as HostM in the text is shown in the main left panel. It is moderately magnified and represents the emission from the entire galaxy. In the inset, the F555W-band image shows HostM with the elliptical aperture (red line) used to compute PSF-matched photometry and the knots 5.1, 5.2 and 5.3. The shaded green ellipses mark the photometric excess in the F160W band, due to nebular emission of [O III]λ4959 mainly due to the YMC 5.1. High S/N nebular SED-fit of YMC 5.1 is shown in the upper right panel, including both 5.1b and 5.1c to highlight the clear photometric jump in the F160W band, as confirmed from VLT/X-shooter observations (Appendix A). The inset shows the images of 5.1b,c and the adopted aperture in the F555W band (blue). Same SED-fit of the complex 5.3 h is shown in the bottom right panel as an example in which no evident nebular contribution is present in the F160W band. The inset shows the image of 5.3 h and the adopted aperture in the F555W band (blue ellipse).

In the text
thumbnail Fig. 5.

Reconstruction of the HostM galaxy in the source plane (delensed) is shown in the left panel, with almost all the details merged into a compact, barely resolved source. Middle panels: galaxy HostM in the F555W-band and the segmentation image with the blue contour at 5σ (left and right, respectively). The three regions 5.1, 5.2, and 5.3 (white apertures with 0.2″ diameter) containing star clusters account for at least 40% of the ultraviolet light of the galaxy. Right panel: same region of the galaxy is outlined with the blue contour in the F275W band probing LyC. It is clear the conspicuous escape of ionizing radiation from the YMC 5.1.

In the text
thumbnail Fig. 6.

Cluster formation efficiency, Γ, as a function of the star-formation rate surface density ΣSFR. The compilation of gray symbols represent estimates in the local Universe taken from Adamo et al. (2020a,b). The solid blue line reproduces the Kruijssen (2012) fiducial model. The estimated lower limit on Γ (calculated over 20 Myr time scale, see text for details) for the Sunburst galaxy is reported with the filled red circle at 30%, while the transparent circle is the estimate fraction of the ultraviolet light (TL(UV)) computed from the delensed UV light integrated over all the star clusters reported in Table 1. The inferred quantities are likely lower limits, as additional clusters might have remained undetected.

In the text
thumbnail Fig. 7.

Main panel: HST/F555W-band image (30 mas/pixel) with superimposed the field of views (FoV) of future instruments, such as VLT/MAVIS (IFU and imaging modes), ELT/MAORY-MICADO (imaging), and HARMONI (IFU), along with the layout of the MUSE-NFM targeting arc IV (transparent green box in the bottom-left corner, presented in Appendix C). The four arcs of Sunburst are indicated with I, II, III, IV. MUSE-NFM requires the presence of a very bright star within 3.4″ from the target, whilst future AO facilities will definitely probe stellar clusters down to the parsec scale in the most magnified regions. Right panels: zoomed region indicated with the orange shaded rectangle located in the main panel over the objects 5.2 h and 5.3 h (left). In particular, the region enclosing the knots 5.2 h and 5.3 h will be covered with MAVIS and HARMONI IFUs which will probe 1−5 parsec per spaxel, while imaging with MAVIS and MAORY-MICADO (wide blue and magenta squares in the main panel) will probe the parsec scale in the optical and near-infrared, respectively.

In the text
thumbnail Fig. A.1.

Main panel (thick line box): the [O III]λ5007 (left) and the Hα (right) two-dimensional zoomed spectra extracted from X-Shooter observations at spectral resolution R ≃5600. The slit orientation is indicated with the yellow box in the inset showing the HST color image. The magenta dotted lines mark the position of the knots 5.1b, c and 5.1d, e, f over the spatial direction along the slit. The blue dotted lines show the expected positions of the region collecting the group of fainter knots, named “group” in the figure (see also Figure 2). Hα and [O III]λ5007 clearly emerge from knots 5.1b, c, d, e, f, conversely the same lines are not detected for the two mirrored “groups” at the given depth. The inferred rest-frame equivalent widths are also quoted. The transparent vertical blue stripes mark the position of night sky emission lines. The right panel shows the Hα equivalent width as a function of time from a Starburst99 model of an instantaneous burst and Z = 0.008Z. Magenta and blue dotted arrows indicate the case of the knot 5.1 and “Group”, respectively. Knot 5.1 is younger than 5 Myr (in agreement with the 3 Myr old age derived by Chisholm et al. 2019), while the complex “group” is likely older than 7 Myr.

In the text
thumbnail Fig. B.1.

[O III]λ5007 line emerging from the LyC knot 5.1 (gray band) and the best-fit solution made of the sum of two Gaussian components (red line), a broad one tracing the outflowing ionized gas (blue dotted) and the other tracing the core of the line; the latter arises from the H II region surrounding the star cluster and tracing the redshift of the stellar component (green dotted). On top axis, the velocity is referred to the peak position of the line. The clear blue tail extends up to −400 km s−1.

In the text
thumbnail Fig. C.1.

Counter arc IV observed with HST ACS/F606W (left) and VLT/MUSE NFM after collapsing the datacube along the indicated wavelength intervals (middle and right images). The ground-based AO-assisted observations of arc IV outperform those from space (see narrow-band imaging after collapsing the cube in the wavelength ranges 6500-8000Å and 8000-9000Å), providing further evidences for the compactness of the star-forming region and in particular the LyC star cluster 5.1 (in this case, the image 5.1n has been targeted). The right most figure shows the FWHM of the OGS (On-axis Guide Source), measured within a 1.5″ region, as a function of the wavelength. The black line show the measurements in the final stacked data-cube, whereas the blue region indicates the scatter (68% confidence level) of the measurements in each single exposure. The FWHM is < 60 mas at λ > 7000Å and we assume a similar resolution is achieved on arc IV (middle).

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.