© The Authors 2023

1. Introduction

For many years, the existence and occurrence of very massive stars (VMS) was mostly associated with the early Universe and metal-free environments in the context of the so-called Population III stars (e.g., Abel et al. 2002). VMS are short-lived stars that are ∼2–3 Myr old (e.g., Yusof et al. 2013), have a mass M > 100 M_⊙ (Vink et al. 2015), and predominantly populate the central regions of young massive star clusters (within the core radius r_c ∼ 0.1 − 0.2 pc; Portegies Zwart et al. 2010). Studies of VMS in Milky Way star clusters are limited to only few targets because their lifetimes are so short. Example are the Arches cluster (Martins et al. 2008) or NGC 3603 (Crowther et al. 2010). Individual VMS have been investigated in the local Universe with high spatial resolution with the Hubble Space Telescope (HST; Cignoni et al. 2015; Crowther et al. 2016; Calzetti et al. 2015; Smith et al. 2016, 2020; Brands et al. 2022). Very massive stars with masses higher than 100 M_⊙ are recognized as objects that strongly affect the evolution of early galaxies in their chemical enrichment and star formation through feedback (e.g., Goswami et al. 2021). It is therefore essential to extend the upper masses beyond 100 M_⊙ of the current population synthesis models to investigate and understand young massive star clusters and VMS at different redshifts (Smith et al. 2016; Crowther et al. 2016).

Despite some progress, the maximum stellar mass that can be attained and the conditions determining the presence of VMS largely remain unknown. Recent observations of local star clusters reported initial stellar masses up to ∼270 M_⊙ (Brands et al. 2022) in the star cluster R136, with a cluster age of ∼1.5 Myr (Crowther et al. 2016). Furthermore, observations of young stellar clusters have revealed peculiar spectroscopic features such as unusually strong broad He IIλ1640 emission (with full width at half maximum (FWHM) > 1000 km s⁻¹), which suggests the presence of VMS in these objects (e.g., Wofford et al. 2014; Crowther et al. 2016; Senchyna et al. 2021).

Alongside with the observations, different models are trying to predict and trace the evolution of the various spectroscopic features characteristic of VMS through different ages and masses (e.g., Köhler et al. 2015; Gräfener 2021). For example, Martins & Palacios (2022) have generated new evolutionary models and synthetic spectra of stars with initial masses in the range 150–400 M_⊙, taking the specific winds of VMS into account that are stronger from a simple extrapolation of OB-type winds at high luminosity. The resulting models predict specific features in the UV and optical part of the spectra that are characteristic signatures of VMS. The most robust ultraviolet spectral features associated with VMS are N IVλ1486, broad He IIλ1640 (and narrow; see Gräfener & Vink 2015) emission, and the N IVλ1720 P Cygni profile (Martins & Palacios 2022). Equivalent widths of these lines are expected to span the interval of 0.1 − 7 Å in the rest frame. Spectra with a high signal-to-noise ratio (S/N) with well-detected continuum are therefore required to identify them, as was shown by Crowther et al. (2016) in the R136 stellar cluster in the local Universe, for example. At cosmological distance, strong gravitational lensing is necessary to detect these faint spectral features, which allows us to further gain in spatial resolution at a scale of tens of parsec and depth (see also Vanzella et al. 2016, 2017, 2021; Vanzella et al. 2022b; Johnson et al. 2017; Rigby et al. 2017, 2018a, 2018b; Meštrić et al. 2022).

We present convincing spectroscopic evidence for VMS in a stellar cluster at cosmological distance for the first time (z = 2.37; Vanzella et al. 2020a, 2022a). The host galaxy is called Sunburst (Rivera-Thorsen et al. 2019, 2017), and Lyman continuum (LyC) radiation is detected from the same clumpy regions in which the VMS were detected. These massive stars in the center of the stellar cluster are significant producers of LyC radiation, through which a porous interstellar medium (ISM) is created that enables LyC escape. For the purpose of this work, we performed a comprehensive analysis of deep Very Large Telescope (VLT)/Multi Unit Spectroscopic Explorer (MUSE), X-shooter, and synthetic spectra with the aim to confirm the presence of the VMS. Additionally, we investigated whether VMS are segregated by modeling the morphology of the young massive star cluster (YMC) that is hosted in the Sunburst galaxy.

The paper is organized as follows. In Sect. 2 we briefly describe the Sunburst galaxy and the available observational data. In Sect. 3 we analyze the spectral signatures of very massive stars using MUSE/IFU and X-shooter observations in combination with the latest evolutionary models and synthetic spectra. In Sect. 4 we discuss the morphological properties of the YMC (called source 5.1) and the possible segregation of the (very) massive stars in its central parts. We present our conclusion in Sect. 5.

We assume a flat cosmology with Ω_M = 0.3, Ω_Λ = 0.7 and H₀ = 70 km s⁻¹ Mpc⁻¹. Within this model, one arcsec at z = 2.37 corresponds to a projected physical scale of 8200 parsec. All magnitudes are given in the AB system.

2. Lensed galaxy Sunburst

The Sunburst galaxy is a galaxy at z = 2.37. It is strongly lensed by the Planck cluster PSZ1 G311.65-18.48 at z = 0.44, as initially reported by Dahle et al. (2016). The strong gravitational lensing effect deflects the light from the background high-z Sunburst galaxy into four bright arcs. These bright arcs harbor at least 13 star-forming knots that likely are stellar clusters. There are more than 50 multiple images of this system (Pignataro et al. 2021), whose physical properties are studied in detail in Vanzella et al. (2022a). One of the 13 young stellar clusters has been identified 12 times (called source 5.1) and is the subject of this work (see Fig. 1). The multi-peaked Lyα emission of source 5.1 is consistent with an optically thin medium and LyC leakage along the line of sight (Rivera-Thorsen et al. 2017). Furthermore, the detection of LyC radiation emerging from the 12 detected multiple images of the young massive star cluster in source 5.1 was confirmed by HST multiband observations (Rivera-Thorsen et al. 2019). Additional analyses of the 12 LyC multiple images of source 5.1 have revealed that the star cluster is younger than 3 Myr and has a stellar metallicity of 0.5 Z_⊙ (Chisholm et al. 2019), with a physical size of ≃10 pc and a stellar (and dynamical) mass value of ≃10⁷ M_⊙ (Vanzella et al. 2022a).

Fig. 1. Sunburst bright arcs with a visible young massive stellar cluster. Left: HST F555W band image, showing the six aperture positions from which a MUSE 1D spectrum is extracted (red contours). White arrows point to the multiple images of the young stellar cluster. Right: MUSE IFU image at ∼ 1800 Å of the same region as on the left. The same apertures are plotted in red.

The Sunburst was observed with the HST, providing multiband photometry in the F275W, F410M, F555W, F606W, F814W, F098M, F105W, F140W, and F160W filters under the programs 15101 (PI Dahle), 15949 (PI Gladders), and 15377 (PI Bayliss). Sunburst has also been targeted with ground-based high-resolution (R ∼ 5000 − 9000) VLT/X-shooter spectroscopy covering the spectral range 3000–22 000 Å in three main arms, UVB, VIS, and NIR. The observational strategy and the data reduction procedures applied to HST imaging and VLT/X-shooter spectroscopy have been presented in Vanzella et al. (2020b,a). VLT/MUSE integral field spectroscopy at a resolution R = 3000 and covering the spectral range 4800–9400 Å was obtained during 2016 (1h integration, DDT; PI: Aghanim) and 2021 (1h integration; PI: Vanzella) in the wide-field mode configuration. The final datacube that combines the two hours and the data reduction is described in Vanzella et al. (2022a). We also presented a first version of the lens model in Pignataro et al. (2021) based on the 62 spectroscopically confirmed multiple images in the redshift range 1 < z < 3.5 (see also Sharon et al. 2022; Diego et al. 2022). A revised lens model will be computed when the new VLT/MUSE observations (7h integration) that are planned for 2023 will be performed (prog. 110.249D.001; PI: Vanzella).

Here we focus on the ≃3 Myr old, UV-bright, and LyC source with M_UV = −18.6 (1700 Å magnitude and ultraviolet slope β = −1.71 ± 0.01, F_λ ∼ λ^β), massive (M ∼ 10⁷ M_⊙) star cluster 5.1, which has high magnification values (μ ∼ 10 − 70 over 12 multiple images; Pignataro et al. 2021; Vanzella et al. 2022a). In the following, we perform a new analysis focusing on the nature of the ionizing source (Sect. 3) and its morphology (Sect. 4).

3. Spectral signatures of very massive stars in the Sunburst star cluster at z = 2.37

We investigate the UV and optical spectroscopic properties of the young stellar cluster 5.1. The VLT/MUSE one-dimensional spectra were extracted from six apertures enclosing nine multiple images of source 5.1 (shown in Fig. 1) and were subsequently combined to produce a continuum-detected high S/N (> 60) weighted-average spectrum (Fig. 2). The stacked spectrum shown in Fig. 2 is equivalent to an integration time of (2 × 9)×30² > 16 000 hours without lensing amplification when the minimum amplification in the 9 multiple images is adopted (μ = 30).

Fig. 2. MUSE IFU spectrum of the young massive star cluster 5.1 extracted from six apertures (thin black line) and the X-shooter long-slit spectrum of 5.1l knot (bold blue line). The key confirmed features indicating the presence of VMS in the stellar cluster are marked with shaded light red strips, and the dark orange line shows the best-fit Fλ ∼ λβ, with β = −1.71. The shaded gray strip indicates the (absence of) O Vλ1371 line, which usually is an indicator of VMS too. Prominent P Cygni of N Vλ1240 and the strong emission part of the C IVλ1550 indicate the young age of the stellar cluster (bold black markers). In the bottom panel, the 1σ errors of both spectra are shown. Other detected interstellar features and stellar features are marked with dashed red and green lines, respectively.

3.1. Observed VMS features with VLT MUSE and X-shooter

The very high S/N MUSE spectrum (Fig. 2) allows us to identify several emission and absorption lines. They include the nebular emission lines associated with the interstellar medium of the galaxy, such as O III]λ1661,1666, N III]λ1750, [Si IIIλ]1883,1892, and C III]λλ1907,1909. The well-detected continuum allows us to investigate faint line emissions (of a fraction of an Å rest-frame equivalent width), and to sample the details of the line profiles that are otherwise not accessible without lensing amplification. In particular, faint N IVλ1486 emission, the evident P Cygni profile of the C IVλ1550, the prominent broad and asymmetric He IIλ1640 line profile, and the P Cygni signature of N IVλ1720 clearly stand out. All these lines are associated with young, hot, and (very) massive stars.

We report the detection of He IIλ1640 emission with a measured rest-frame EW = 3.0 ± 0.3 Å and FWHM = 8.8 ± 1.7 Å (∼1610 ± 300 km s⁻¹). The broad shape of the He IIλ1640 emission line observed in the Sunburst cluster is asymmetric and resembles a typical P Cygni profile. The blue end of the emission line drops steeply, while the red end drops more gradually. The P Cygni profile of He IIλ1640 line is consistent with the profile predicted by models and synthetic spectra (see Martins & Palacios 2022; Gräfener & Vink 2015). Broad He IIλ1640 emission observed in galaxies is usually related to a non-nebular origin and is commonly associated with Wolf-Rayet (WR) stars (e.g., Schaerer & Stasińska 1999; Brinchmann et al. 2008; Gräfener & Vink 2015; Leitherer et al. 2018; Senchyna et al. 2021), although the failure of the synthesis models to reproduce some of the strong He IIλ1640 lines might be related to missing ingredients in stellar evolution models (see, e.g., Leitherer et al. 2018). However, far-UV spectroscopic investigation of ∼57 individual stars located within the R136 star cluster revealed that massive stars with M > 100 M_⊙ are crucial for producing the He IIλ1640 emission line (Gräfener & Vink 2015; Crowther et al. 2016). On the other hand, Martins & Palacios (2022) have shown that He IIλ1640 can be produced in significant amount only when stellar winds are stronger than in normal O stars (see also Gräfener & Vink 2015). VMS develop these strong winds because they lie close to the Eddington limit (Vink et al. 2011; Bestenlehner 2020; Gräfener 2021). At the same time, these winds peel off the external layers of the stars and expose the products of hydrogen burning to the surface through the CNO cycle. This results in a strong nitrogen (and helium) enrichment that boosts the strength of N IVλ1486 and N IVλ1720. This typically occurs after ∼1.5 Myr. The two emission lines are both detected in the spectrum of the Sunburst cluster at S/N > 15 (Fig. 2), with EW = 0.2 Å and FWHM = 2.9 Å for N IVλ1486 and EW = 0.15 Å and FWHM ∼ 2 Å for N IVλ1720. The helium enrichment also contributes to the strength of He IIλ1640. The same effects (strong winds combined with surface chemical enrichment) occur in normal evolved massive stars when they are seen as WR stars. The key difference to VMS is that helium enrichment takes place only after the main sequence (≳4 Myr), while the same process takes place at younger ages in VMS. Furthermore, VMS are more luminous than normal WR stars, and they accordingly contribute more strongly to integrated light.

The nebular Hα equivalent width provides constraints on the cluster age and hence on whether the He IIλ1640 line is primarily due to WR stars or to VMS. According to the BPASS models and results from Eldridge & Stanway (2012, their Fig. 3), normal and WR stars produce EW_Hα < 1 Å for ≲3 Myr. The X-shooter spectrum reveals a prominent Hα line and no continuum detection, which very conservatively implies an equivalent width larger than 200 Å rest frame at 1σ. However, when we assume for Hα the same continuum level as observed at λ ∼ 5000 Å rest frame in the photometric spectral energy distribution (SED) by Vanzella et al. (2022a), this limit increases to ∼840 Å. This value would be still a lower limit, even in the case of leaking ionizing photons. After correcting the Hα flux for the fraction of escaping LyC photons ()), the resulting EW increases to EWHα ∼ 1231 Å. We adopted values from Rivera-Thorsen et al. (2019), where the corresponding absolute escape fraction of LyC photons along the line of sight is %. This large Hα equivalent width is consistent with a star-forming burst younger than ∼3 Myr (e.g., Leitherer et al. 2014). Furthermore, as discussed in Chisholm et al. (2019), the N Vλ1240 stellar wind profile predominantly depends on the stellar age, while variations due to different metallicity are negligible, and it is related to the young stellar populations (< 5 Myr). From the comparison of the observed N Vλ1240 with the models Figs. 3 and 4, we additionally demonstrate that the age of the cluster is < 3 Myr. From our age analysis, we can conclude that the properties of both N Vλ1240 and Hα fit a < 3 Myr age of the stellar cluster well, which requires other sources than WR stars to explain the observed strong He IIλ1640 EW = 3.0 ± 0.3 Å. Therefore, these results strongly suggest that VMS produce the spectral ultraviolet features we observe in this young massive star cluster. Moreover, Wofford et al. (2014) and Smith et al. (2016) have argued that O Vλ1371 in the integrated light of the clusters is a key feature of VMS as well. This line is not seen in the Sunburst cluster; see Fig. 2. As demonstrated by Martins & Palacios (2022), this is not incompatible with the presence of VMS because O Vλ1371 disappears as VMS evolve to a lower effective temperature. In their Fig. 4, no sign of O Vλ1371 exists after ∼1 Myr. Together with N IVλ1486, this places a rather tight constraint on the cluster age.

Fig. 3. MUSE ultraviolet spectra of the young star cluster (blue, bottom panel) with the X-shooter spectrum (green). The shaded gray regions show specific UV features that are closely associated with the VMS (N IVλ1486, He IIλ1640, and N IVλ1720), and some of them show a P Cygni profile, as is characteristic of young and massive stars. Furthermore, the black line shows the single BPASS model including only normal stars at 2.5 Myr, and the red line shows the BPASS single-star model augmented by VMS with masses up to 150 M⊙. The upper panels show the zoom in VMS characteristic features compared with the models and the N Vλ1240 P Cygni line, which is characteristic of very young stellar populations.

Fig. 4. All symbols are the same as in Fig. 3, except for the red lines, which show the BPASS single-star model augmented by VMS with masses up to 200 M⊙.

3.2. Comparing observations with models

To investigate the rest-frame UV spectrum of the cluster, we created an integrated VMS model following Martins & Palacios (2022) that includes normal-mass stars (0.1–100 M_⊙) with different VMS masses (150 M_⊙ and 200 M_⊙). It is important to note that our models rely on the predictions of Martins & Palacios (2022), which are based on the empirical mass-loss prescriptions of Gräfener (2021). Different mass-loss prescriptions can affect the evolution of VMS, as shown by Sabhahit et al. (2022), for instance. The HI ionizing budget mostly depends on the effective temperature (T_eff) and less on an exact wind prescription (e.g., Martins & Palacios 2022). Furthermore, because winds can strongly affect stellar evolution and thus T_eff, the estimate of the contribution of VMS to the ionizing budget we presented in Sect. 3.3 is model dependent. However, we stress that our models are able to reproduce the integrated spectroscopic properties of R136 (Martins et al., in prep.) and are thus appropriate.

We used the SEDs of the BPASS (Eldridge et al. 2017; Stanway & Eldridge 2018) v2.2.1 single-star population synthesis model. The model has an upper mass limit of 100 M_⊙ with a Salpeter IMF, a metallicity of Z = 0.006 (where 0.02 corresponds to solar metallicity), and an instantaneous star formation history, with a burst of mass 10⁶ M_⊙. The adopted metallicity of the model is closest to our measured value based on the N2 index (Marino et al. 2013), which is ≃0.4 Z_⊙ and is consistent with the estimate provided by Mainali et al. (2022; see also Chisholm et al. 2019).

We extrapolated the Salpeter IMF to a 225 M_⊙ upper-mass limit within a few mass bins given by Eq. (1) from the BPASS manual¹. Equation (1) gives the number of massive stars in the mass range [M_a; M_b],

(1)

Here, C is a constant and has the value of 1.23 × 10⁵ for an arbitrary burst mass of 10⁶ M_⊙. Moreover, M₁ = 0.5 M_⊙α₁ = −1.3, and α₂ = −2.35. The mass bins were selected in a way to add the SEDs of appropriate numbers of single VMS stars, which are available for discrete sets of VMS with masses including 150 M_⊙ and 200 M_⊙. In this manner, we computed SEDs including VMS with IMFs extending up to 175 and 225 M_⊙, respectively, following Martins & Palacios (2022). In the mass bin from 100 M_⊙ to 175 M_⊙, the extrapolation of the Salpeter IMF using Eq. (1) gives ∼237 VMS, and for the mass bin from 175 M_⊙ to 225 M_⊙, it provides us with ∼60 VMS. We added SEDs of 150 M_⊙ VMS as an average of the VMS population in the mass bin 100 M_⊙ ≤ M ≤ 175 M_⊙, and we did the same by adding SEDs of 200 M_⊙ VMS for the mass bin 175 M_⊙ ≤ M ≤ 225 M_⊙. Finally, these numbers have to be rescaled to the appropriate cluster mass, which is ∼10⁷ M_⊙ (Vanzella et al. 2022a).

The YMC shows significant N IVλ1486 emission, as shown in Fig. 2. From the VMS models and synthetic spectra, N IVλ1486 emission only appears after 1.5 Myr of VMS evolution (see Martins & Palacios 2022), and VMS last approximately until 2.5 Myr. Based on this, we created the SEDs of integrated VMS models at 1.5 Myr, 2 Myr, and 2.5 Myr. We normalized the spectrum of the cluster and the models by fitting a UV power law by using the spectral windows provided by Rix et al. (2004).

We directly compared the cluster spectrum with the two VMS models at three different ages. The comparison shows that VMS are clearly needed to reproduce the observations (see Figs. 3, 4, and A.1). However, the He IIλ1640 and N IVλ1720 lines in the models appear stronger than observed even at an age of 1.5 Myr and with a maximum mass of 175 M_⊙. To match the observed He IIλ1640 and N IVλ1720 profiles, we therefore reduced the VMS contribution by decreasing their numbers. When the VMS contribution is reduced by a factor of 6 in the VMS model, which includes only 150 M_⊙ VMS at 2.5 Myr (Fig. 3), the results agree well. Alternatively, a similar match is also found by reducing the VMS contribution by a factor of 8 in models that also incliude the 200 M_⊙ VMS (Fig. 4). For a cluster mass of 10⁷ M_⊙, this implies a total of ∼370 − 400 VMS for both cases. In short, the observations are compatible with an IMF extending up to ∼175 or 225 M_⊙, but with an IMF slope steeper than for a Salpeter IMF (α₂ < −2.35) for M > 100 M_⊙.

3.3. VMS contribution to the Lyman continuum budget of the young stellar cluster

The contribution of the VMS to the total LyC emission can be determined in a straightforward way from the theoretical SEDs by comparing models with and without VMS by integrating the theoretical SEDs in the defined LyC range. From the favored model including VMS, but with a steep IMF slope, and a model excluding VMS, we find that VMS contribute ∼15% of the LyC radiation that is emitted over the range ∼730–900 Å rest frame (covered by HST F275W filter). Similar percentages are also obtained by integrating over the total LyC. For comparison, we note that this percentage is somewhat lower than the fraction of LyC photons produced by VMS in R136, which is found to be 25–40% (Doran et al. 2013; Crowther et al. 2010). This is consistent with our inference of fewer VMS than in a standard IMF. Since LyC photons are mainly produced by O-type and more massive stars, ∼15% of the LyC production is generated from only ∼1% of the stars that are capable of producing LyC photons.

4. Spatial segregation of the Lyman continuum radiation

4.1. Morphological analysis of 5.1l YMC

We now address the morphological properties of image 5.1l, which is the most strongly magnified of the multiple images of the star cluster (μ_tot ≃ 76; Pignataro et al. 2021). Image 5.1l is the brightest image and was detected with a high S/N in the F275W (S/N ∼ 90) and F555W (S/N ≫ 100), allowing us to investigate and compare the morphology in these two spectral regions: the emitting LyC (λ < 900 Å in HST F275W band) and the nonionizing radiation at 1700 Å (HST F555W band). We followed two approaches: (1) We ran simulations injecting the sources in the F275W band, and (2) we analyzed the curve of growth of the resulting images.

Figure 5 shows the F555W image of 5.1l, in which the elongation is clearly visible in the direction of the tangential stretch produced by gravitational lensing. As discussed in Pignataro et al. (2021; see also Vanzella et al. 2022a), the tangential amplification largely dominates along the arc (μ_tang ≃ 57). We perform a relative comparison between images to ensure that the results do not depend on the magnification values.

Fig. 5. Results from Galfit modeling after using a two-component fit. From the left, the first panel shows source 5.1l in the F555W band (UV 1700 Å). The second panel shows the model from Galfit. The third panel shows the residual, and the fourth panel is the normalized residual produced after dividing the residual with the original image. The white contour encloses the region we used to check the quality of the produced model.

As a first step, we computed a realistic model of image 5.1l in the HST F555W image using Galfit (Peng et al. 2010). The point spread function (PSF) was extracted by combining unsaturated stars in the field of view. While the fit with a single component did not produce acceptable residuals (higher than 20%), we reproduce the light profile of the object quite well by combining two components: a core with a Gaussian light profile and an effective radius (R_eff) smaller than 0.5 pixels (in practice nearly unresolved), and an extended component with R_eff = 6 pixels and a Sérsic index n = 1 (similar results are also obtained with n = 0.5). The combination of the two components produces an optimal shape that leaves normalized residuals lower than 10% (see Fig. 5). It is worth investigating whether this resolved shape (sampled at 1700 Å) is recovered when it is placed in the F275W LyC image. For this check, we injected mock images of source 5.1l into the F275W image in five different positions around source 5.1l that were not contaminated by the flux from other sources. These images were produced from the aforementioned two-component Galfit model constructed at 1700 Å (F555W), but now accounting for the F275W PSF (i.e., convolved by the F275W PSF) to allow for a proper comparison with the LyC source 5.1l (observed in the HST F275W band). These images were added to F275W after each of them was rescaled to the observed peak value of the LyC source 5.1l. This step was performed with the IRAF (Tody et al. 1986) tasks IMARITH and IMCOPY. Figure 6 shows the results, in which all the injected images show a spatially resolved morphology along the tangential magnification. Conversely, the observed LyC image (of source 5.1l) appears nucleated, suggesting that the emitting LyC region is smaller than the region at 1700 Å.

Fig. 6. HST F275W image showing source 5.1l leaking LyC in its center. The other two multiple images of the same source are labeled 5.1i and 5.1h. In the region not contaminated by other sources, five models are injected around source 5.11 (see text for more details). They are enclosed in the white squares. In the bottom right corner, we show the largest (diameter of 34 pix) aperture size used to construct CoGs.

To quantify this result, we calculated the curve of growth (CoG) of the images shown in Fig. 6. The flux was then measured in the F275W band in 34 circular apertures. The smallest aperture had a radius of 0.1 pixel. Intermediate apertures were drawn with increasing radii, with a step of 0.5 pix, up to the largest aperture, which has a radius of 34 pixels. As a reference point-like source, we constructed the mean CoG from a selected sample of 20 unsaturated and uncontaminated stars. The resulting CoG is shown in Fig. 7, where the y-axis reports the fraction of the flux that is enclosed in the corresponding radius in pixels (x-axis).

Fig. 7. Three curves of growth normalized to 1. The cyan CoG corresponds to source 5.1l leaking LyC that is located in the Sunburst arc (F275W band). The violet CoG is the PSF-convolved best-fit model of source 5.1l observed in F555W, constructed by averaging 5 CoGs (see text for more details). The orange growth curve is used for comparison and was constructed by averaging 20 single CoGs from randomly selected stars. In both cases (cyan and violet), the error bars are 1σ. Vertical lines in the bottom left part of the panel mark the pixel radii within which 50% of the light is enclosed. The colors corresponds to the CoG colors.

The same procedure was applied to the LyC-emitting source 5.1l, while another CoG was constructed by averaging the five CoG of the injected models that resemble the morphology at 1700 Å. Figure 7 compares all the CoGs after normalizing them to the saturation value at the largest radius. The first result from this test is that the CoG of the observed 5.1l LyC source clearly deviates from the behavior of a point-like source (stars). This was not explored before and suggests that in the most magnified image of the star cluster, the LyC appears to be spatially resolved. This is the first evidence of resolved stellar LyC emission at cosmological distance. Second, this barely resolved LyC emission appears to be more nucleated than the emission at 1700 Å. Consequently, the sources of ionizing radiation appear to be located in the central part of the cluster. These curves show that 50% of the stellar flux is enclosed within a radius of ∼1.9 pixel, while for source 5.1l, the flux lies within ∼2.2 pixels. Additionally, we performed a Kolmogorov-Smirnov two-sample test (KS-test) to determine whether the CoGs derived from the stars and from source 5.1l follow the same distribution (null hypothesis). For this purpose, we used the statistical function ks_2samp from scipy.stats. After comparing the average CoG of the stars with cyan and violet CoGs, the KS-test gives p ≪ 0.05. This means that the null hypothesis is not satisfied with the LyC profile of source 5.1l, and it deviates from the CoG of a point-like source. Furthermore, we also find that the half-light size of 5.1l at 1700 Å is larger than the ionizing region, ∼2.6 pixels compared to the ∼2.2 pixels. When we correct these radii for the instrumental resolution (given by the stars), we obtain an effective radius for source 5.1l at 1700 Å of ≃7.8 ± 1.4 pc after delensing². This agrees with Vanzella et al. (2020a), while the LyC image (5.1l) has a smaller radius, R_eff ≃ 4.7 ± 1.5 pc. The LyC emission is therefore more compact than the nonionizing UV continuum, which we interpret as a spatial segregation of the most massive stars.

4.2. Segregation of (very) massive stars?

The segregation of massive stars is observed in young stellar clusters (Harayama et al. 2008; Habibi et al. 2013), and its origin is still an open question. Two scenarios for mass segregation within young star clusters have been proposed. The first scenario suggests that mass segregation is primordial (e.g., Dib et al. 2010). The second scenario that can explain mass segregation requires dynamical evolution. In this case, massive stars migrate to the central parts of the cluster as a consequence of two-body encounters with less massive stars (e.g., Binney & Tremaine 1987).

Mass segregation has so far been investigated and reported in several stellar clusters located in the Milky Way (e.g., NGC 3603; Harayama et al. 2008) or in the Local Universe (e.g., R136; Khorrami et al. 2021). However, at cosmological distances as in the case of source 5.1 YMC presented in this work, the analysis described in the previous section suggests that the regions emitting ionizing and nonionizing radiation have different sizes and that the bulk of LyC radiation arises from the core of the star cluster. This implies that the only producers of LyC radiation, that is, hot massive stars, are segregated in the central parts of the cluster and are surrounded by a more diffuse distribution of less massive, non-LyC emitting stars.

Dynamical relaxation occurs on timescales comparable with the relaxation time t_relax. When t_relax exceeds the age of the cluster, the most likely origin for segregation is therefore primordial. The relaxation time in the center of the cluster can be estimated as t_relax ⋍ (0.065 × σ³)/(G²mρlnΛ) (Spitzer 1987), where m is the average stellar mass in the cluster (the mass of the cluster divided by the total number of stars), σ is the central three-dimensional velocity dispersion, ρ is the central mass density, G is the gravitational constant, and lnΛ is the Coulomb logarithm, a number typically between 6 and 15. To compute the central relaxation time of source 5.1l YMC, we assumed that its stellar distribution is well approximated by that of an isotropic self-gravitating Plummer (1911) model. In a Plummer model, the central velocity dispersion is . We adopted a mass of M = 10⁷ M_⊙, as estimated by Vanzella et al. (2022a), and we assumed that in the cluster, light traces mass, so that its half-mass radius is equal to the half-light radius r_0.5 = 8 pc, as measured in F555W. The resulting central velocity dispersion is ≃30 km s⁻¹. This value agrees with other extremely massive star clusters reported in the literature (Bastian et al. 2006; Östlin et al. 2007). The number of stars was calculated from Eq. (2.4) of Stanway & Eldridge (2018) and was rescaled to the mass of source 5.11. The estimated number of stars in source 5.1l YMC is N ∼ 9 × 10⁶. Combined with the total cluster mass, this gives an average stellar mass of m ∼ 1.1 M_⊙. The central density in a Plummer model is , and we adopted lnΛ = 15. The resulting t_relax ∼ 1.1 Gyr clearly implies that stars did not have enough time to change their initial trajectory, which further indicates that the observed segregation is most likely of primordial nature. We also note that our estimate of the relaxation time is a lower limit: the cluster center in our model lies in the place in which the ratio of σ³ and ρ is lowest. The same measure performed at any larger distance from the cluster center would have given an even longer t_relax.

5. Summary and conclusions

We have presented a detailed spectroscopic and morphological analysis of the massive and young stellar cluster hosted in the lensed galaxy called Sunburst galaxy at z = 2.37, for which LyC emission was also confirmed in the literature. We used results from recent stellar evolutions and atmosphere models including VMS (Martins & Palacios 2022) to conduct extensive comparisons with high spectral resolution observations performed with VLT/MUSE and X-shooter. The main results of this work can be summarized as follows: 1) In the spectroscopic observations, the high S/N MUSE and X-shooter spectra reveal features of broad (and asymmetric) He IIλ1640 emission with an EW ≃ 3 Å in the rest frame and a line width of 1610 km s⁻¹, and N IVλ1486 with an EW ≃ 0.2 Å in emission. In addition, the P Cygni profile of N IVλ1720 (along with NV and CIV) is also observed. All these features suggest the presence of very massive (> 100 M_⊙) stars. The absence of O Vλ1371 provides a lower age limit of 1 Myr. On the other hand, the large Hα EW (> 1231 Å after correcting for the escaping LyC radiation) indicates an age younger than ∼3 Myr. These narrow age constraints strongly favor the existence of VMS over WR stars, implying that the strength of the He IIλ1640 emission line is entirely due to VMS. 2) A comparison of the observations with the models revealed that the most plausible age of the star cluster is 2.5 Myr, and the estimated number is ∼370 − 400 VMS for a cluster mass of 10⁷ M_⊙. The observations are compatible with an IMF extending up to ∼175 − 225 M_⊙, but with a steeper slope than the Salpeter IMF. 3) The fraction of LyC radiation emerging from the VMS component is not negligible. We estimate that in the 730 Å–900 Å range (probed by the HST/F275W band), about 370–400 VMS (or roughly 1% of the total population of O-type stars in the star cluster) account for 15% of the escaping LyC photons. The remaining photons are mostly produced by the other less massive O-type stars. 4) Detailed morphological analysis of the most magnified image of the star cluster shows that the region emitting LyC is not point-like, with a light profile different from the average profile of stars present in the same field of view. This is the first evidence of a resolved LyC emission at any redshift. Remarkably, the physical scale of the LyC emitting region also appears to be smaller (with a significant K-S probability p ≪ 0.05) than the nonionizing region (1700 Å), suggesting that massive O-type stars that cause the LyC radiation, and likely the VMS (significantly contributing to it), are segregated in the central part of the star cluster. After delensing the angular half-light radii, the LyC region appears barely resolved with R_eff ≃ 4.7 ± 1.5 pc, while at 1700 Å, it is R_eff ≃ 7.8 ± 1.4 pc. The packaging of this large number of massive O-type stars per parsec cube in the central region, ≃70 pc⁻³, is likely a element that allowed carving the ionizing channel and the development of a high-speed outflowing gas (Rivera-Thorsen et al. 2017; Vanzella et al. 2022a; Mainali et al. 2022).

As a final remark, it is important to stress that we investigated this remarkable object and its exotic properties based on gravitational lensing. Otherwise, this investigation would be impossible with current instrumentation in blank fields at cosmological distances. This once again demonstrates the potential and importance of gravitationally lensed systems for investigations of the most exotic phenomena in the Universe.

Acknowledgments

We acknowledge financial support through grants PRIN-MIUR 2017WSCC32, 2020SKSTHZ and the INAF GO Grant 2022 “The revolution is around the corner: JWST will probe globular cluster precursors and Population III stellar clusters at cosmic dawn” (PI Vanzella). FC and RP acknowledge funding from PRIN INAF 1.05.01.85.01.

References

Appendix A: Initial models

As described in Section 3.2, we compared the X-shooter and MUSE spectroscopic observations with the BPASS models at different ages (1.5 Myr, 2 Myr, and 2.5 Myr). We narrowed our models to the mentioned age range because the predicted age of the cluster exceeds 1.5 Myr (inferred from the N IVλ1486 emission line) and the lifetime of the VMS is about 2.5 Myr (Martins & Palacios 2022). We started with the BPASS that has an upper mass limit of 100 M_⊙ and added 236.56 stars in the 100 - 175 M_⊙ mass range and 60.30 stars in the mass range 175 – 225 M_⊙. The resulting models (at different ages) are shown in Fig. A.1. All models produce significantly strong spectroscopic features that are characteristic of VMS. To better match models with observations, we decreased the numbers of the VMS, and the final results are presented in Section 3 (Figs. 3 and 4).

Fig. A.1. Comparison of the observations (MUSE spectrum, red line) with BPASS models (black line). The left column shows three BPASS models with an IMF up to 100 M⊙, ages of 1.5 Myr, 2 Myr, and 2.5 Myr and an added number of 236.56 150 M⊙ VMS. The right column shows a BPASS model with an IMF up to 100 M⊙ at same ages (as shown in previous column), but with an added number of 236.56 and 60.30 VMS of 150 M⊙ and 200 M⊙.

All Figures

	Fig. 1. Sunburst bright arcs with a visible young massive stellar cluster. Left: HST F555W band image, showing the six aperture positions from which a MUSE 1D spectrum is extracted (red contours). White arrows point to the multiple images of the young stellar cluster. Right: MUSE IFU image at ∼ 1800 Å of the same region as on the left. The same apertures are plotted in red.
In the text

Fig. 2. MUSE IFU spectrum of the young massive star cluster 5.1 extracted from six apertures (thin black line) and the X-shooter long-slit spectrum of 5.1l knot (bold blue line). The key confirmed features indicating the presence of VMS in the stellar cluster are marked with shaded light red strips, and the dark orange line shows the best-fit Fλ ∼ λβ, with β = −1.71. The shaded gray strip indicates the (absence of) O Vλ1371 line, which usually is an indicator of VMS too. Prominent P Cygni of N Vλ1240 and the strong emission part of the C IVλ1550 indicate the young age of the stellar cluster (bold black markers). In the bottom panel, the 1σ errors of both spectra are shown. Other detected interstellar features and stellar features are marked with dashed red and green lines, respectively.

In the text

Fig. 3. MUSE ultraviolet spectra of the young star cluster (blue, bottom panel) with the X-shooter spectrum (green). The shaded gray regions show specific UV features that are closely associated with the VMS (N IVλ1486, He IIλ1640, and N IVλ1720), and some of them show a P Cygni profile, as is characteristic of young and massive stars. Furthermore, the black line shows the single BPASS model including only normal stars at 2.5 Myr, and the red line shows the BPASS single-star model augmented by VMS with masses up to 150 M⊙. The upper panels show the zoom in VMS characteristic features compared with the models and the N Vλ1240 P Cygni line, which is characteristic of very young stellar populations.

In the text

	Fig. 4. All symbols are the same as in Fig. 3, except for the red lines, which show the BPASS single-star model augmented by VMS with masses up to 200 M⊙.
In the text

Fig. 5. Results from Galfit modeling after using a two-component fit. From the left, the first panel shows source 5.1l in the F555W band (UV 1700 Å). The second panel shows the model from Galfit. The third panel shows the residual, and the fourth panel is the normalized residual produced after dividing the residual with the original image. The white contour encloses the region we used to check the quality of the produced model.

In the text

Fig. 6. HST F275W image showing source 5.1l leaking LyC in its center. The other two multiple images of the same source are labeled 5.1i and 5.1h. In the region not contaminated by other sources, five models are injected around source 5.11 (see text for more details). They are enclosed in the white squares. In the bottom right corner, we show the largest (diameter of 34 pix) aperture size used to construct CoGs.

In the text

Fig. 7. Three curves of growth normalized to 1. The cyan CoG corresponds to source 5.1l leaking LyC that is located in the Sunburst arc (F275W band). The violet CoG is the PSF-convolved best-fit model of source 5.1l observed in F555W, constructed by averaging 5 CoGs (see text for more details). The orange growth curve is used for comparison and was constructed by averaging 20 single CoGs from randomly selected stars. In both cases (cyan and violet), the error bars are 1σ. Vertical lines in the bottom left part of the panel mark the pixel radii within which 50% of the light is enclosed. The colors corresponds to the CoG colors.

In the text

Fig. A.1. Comparison of the observations (MUSE spectrum, red line) with BPASS models (black line). The left column shows three BPASS models with an IMF up to 100 M⊙, ages of 1.5 Myr, 2 Myr, and 2.5 Myr and an added number of 236.56 150 M⊙ VMS. The right column shows a BPASS model with an IMF up to 100 M⊙ at same ages (as shown in previous column), but with an added number of 236.56 and 60.30 VMS of 150 M⊙ and 200 M⊙.

In the text