# Clues on the presence and segregation of very massive stars in the Sunburst Lyman-continuum cluster at $z=2.37$

U. Meštric, E. Vanzella, A. Upadhyaya, F. Martins, R. Marques-Chaves, D.

Schaerer, J. Guibert, A. Zanella, C. Grillo, P. Rosati, et al.

► To cite this version:

U. Meštric, E. Vanzella, A. Upadhyaya, F. Martins, R. Marques-Chaves, et al. Clues on the presence and segregation of very massive stars in the Sunburst Lyman-continuum cluster at $z=2.37$. Astronomy and Astrophysics - A&A, 2023, 673, pp.A50. 10.1051/0004-6361/202345895. hal-03949079

HAL Id: hal-03949079

https://hal.science/hal-03949079

Submitted on 10 Jun 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Clues on the presence and segregation of very massive stars in the Sunburst Lyman-continuum cluster at $z = 2.37^*$

U. Meštrić<fnref n="1" />, E. Vanzella<fnref n="1" />, A. Upadhyaya<fnref n="2" />, F. Martins<fnref n="3" />, R. Marques-Chaves<fnref n="2" />, D. Schaerer<fnref n="2" />,<fnref n="4" />, J. Guibert<fnref n="2" />, A. Zanella<fnref n="5" />, C. Grillo<fnref n="6" />,<fnref n="7" />, P. Rosati<fnref n="8" />, F. Calura<fnref n="1" />, G.B. Caminha<fnref n="9" />,<fnref n="10" />, A. Bolamperti<fnref n="5" />,<fnref n="11" />,<fnref n="12" />, M. Meneghetti<fnref n="1" />, P. Bergamini<fnref n="1" />,<fnref n="6" />, A. Mercurio<fnref n="13" />,<fnref n="14" />, M. Nonino<fnref n="15" />, and R. Pascale<fnref n="1" />

Received 11 January 2023 / Accepted 9 March 2023

## ABSTRACT

We report the identification of very massive stars (VMS; mass >100 $M_\odot$) that may be segregated in the center of the young massive star cluster at $z = 2.37$ hosted in the lensed galaxy called Sunburst galaxy. This result is based on two pieces of evidence: (1) VLT/MUSE spectra of several multiple images of the same star cluster show key spectral signatures of VMS, such as the HeIIλ1640 broad emission, NIVλ1486 emission, and an NIVλ1720 P Cygni profile. In particular, HeIIλ1640 is broad (~1610 ± 300 km s⁻<fnref n="1" />), with an equivalent width of 3 Å, and asymmetric profile. These features require an extremely young (~2.5 Myr) stellar population component in which the masses of the stars exceed 100 $M_\odot$. When a Salpeter initial mass function and BPASS models for normal massive stars are assumed, the observed spectral features require ~400 VMS. (2) The same star cluster is detected at a signal-to-noise ratio of ~100 in the Lyman continuum domain (λ < 900 Å). The Lyman continuum emission emerges from a region with a radius that is at least twice smaller than what is observed at 1700 Å (independently of magnification) and is located in the center of the cluster. After delensing, the effective radii in absolute scales are $R_{\text{eff}}[\text{LyC}] \sim 4.7 \pm 1.5$ pc and $R_{\text{eff}}[1700] = 7.8 \pm 1.4$ pc. The Lyman continuum radiation is mainly produced by hot and massive stars, which implies that their spatial distribution (including that of VMS) is preferentially more confined in the central parts of the cluster. Approximately 400 VMS hosted by a cluster of ~10⁷ $M_\odot$ produce ~15% of the escaping Lyman continuum photons, and the remaining photons are produced by other massive early-type stars.

**Key words.** stars: massive – Galaxy: stellar content

## 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_\odot$ (Vink et al. 2015), and predominantly populate the central regions of young massive star clusters (within the core

radius $r_c \sim 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_\odot$ 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_\odot$ of the current population synthesis models to investigate and understand young massive star clusters and

<footnote n="1">¹ INAF – OAS, Osservatorio di Astrofisica e Scienza dello Spazio di Bologna, via Gobetti 93/3, 40129 Bologna, Italy
e-mail: uros.mestric@inaf.it</footnote>

<footnote n="2">² Geneva Observatory, Department of Astronomy, University of Geneva, Chemin Pegasi 51, 1290 Versoix, Switzerland</footnote>

<footnote n="3">³ LUPM, Université de Montpellier, CNRS, Place Eugène Bataillon, 34095 Montpellier, France</footnote>

<footnote n="4">⁴ CNRS, IRAP, 14 avenue E. Belin, 31400 Toulouse, France</footnote>

<footnote n="5">⁵ INAF – Osservatorio Astronomico di Padova, Vicolo dell'Osservatorio 5, 35122 Padova, Italy</footnote>

<footnote n="6">⁶ Dipartimento di Fisica, Università degli Studi di Milano, via Celoria 16, 20133 Milano, Italy</footnote>

<footnote n="7">⁷ INAF – IASF Milano, via A. Corti 12, 20133 Milano, Italy</footnote>

<footnote n="8">⁸ Dipartimento di Fisica e Scienze della Terra, Università degli Studi di Ferrara, via Saragat 1, 44122 Ferrara, Italy</footnote>

<footnote n="9">⁹ Technical University of Munich, TUM School of Natural Sciences, Department of Physics, James-Franck-Str 1, 85748 Garching, Germany</footnote>

<footnote n="10">¹⁰ Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Str. 1, 85748 Garching, Germany</footnote>

<footnote n="11">¹¹ Dipartimento di Fisica e Astronomia, Università degli Studi di Padova, Vicolo dell'Osservatorio 3, 35122 Padova, Italy</footnote>

<footnote n="12">¹² European Southern Observatory, Karl-Schwarzschild-Strasse 2, 85748 Garching bei München, Germany</footnote>

<footnote n="13">¹³ Dipartimento di Fisica “E.R. Caianiello”, Università Degli Studi di Salerno, Via Giovanni Paolo II, 84084 Fisciano (SA), Italy</footnote>

<footnote n="14">¹⁴ INAF – INAF – Osservatorio Astronomico di Capodimonte, Via Moiariello 16, 80131 Napoli, Italy</footnote>

<footnote n="15">¹⁵ INAF – Osservatorio Astronomico di Trieste, via G. B. Tiepolo 11, 34143 Trieste, Italy</footnote>

<footnote>* Based on observations collected at the European Southern Observatory for Astronomical research in the Southern Hemisphere under ESO programmes DDT MUSE program ID 107.22SK.001 (PI E. Vanzella), X-shooter program ID 0103.A-0688 (PI E. Vanzella) and DDT MUSE program ID 297.A-5012(A) (PI Aghanim).</footnote>
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_\odot$ (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 HeIIλ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_\odot$, 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 Nivλ1486, broad HeIIλ1640 (and narrow; see Gräfener & Vink 2015) emission, and the Nivλ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, 2022b; Johnson et al. 2017; Rigby et al. 2017, 2018a,b; 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 $\Omega_M = 0.3$, $\Omega_\Lambda = 0.7$ and $H_0 = 70 \text{ km s}^{-1} \text{ Mpc}^{-1}$. Within this model, one arcsec at $z =$

<figure><img src="image_2.png" /></figure>

**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.

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 $\text{Ly}\alpha$ 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_\odot$ (Chisholm et al. 2019), with a physical size of $\approx 10$ pc and a stellar (and dynamical) mass value of $\approx 10^7 M_\odot$ (Vanzella et al. 2022a).

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 \sim 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, 2022a). 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
<figure><img src="image_3.png" /><figcaption>**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.11 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_\lambda \sim \lambda^\beta$, with $\beta = -1.71$. The shaded gray strip indicates the (absence of) $O\text{V}\lambda1371$ line, which usually is an indicator of VMS too. Prominent P Cygni of $\text{Nv}\lambda1240$ and the strong emission part of the $\text{Civ}\lambda1550$ indicate the young age of the stellar cluster (bold black markers). In the bottom panel, the $1\sigma$ errors of both spectra are shown. Other detected interstellar features and stellar features are marked with dashed red and green lines, respectively.</figcaption></figure>

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_{\text{UV}} = -18.6$ (1700 Å magnitude and ultraviolet slope $\beta = -1.71 \pm 0.01$, $F_\lambda \sim \lambda^\beta$), massive ($M \sim 10^7 M_\odot$) star cluster 5.1, which has high magnification values ($\mu \sim 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 ($\mu = 30$).

**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\text{III}]\lambda1661, 1666$, $N\text{III}]\lambda1750$, $[S\text{III}]\lambda1883, 1892$, and $C\text{III}]\lambda1907, 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 $\text{Nv}\lambda1486$ emission, the evident P Cygni profile of the $\text{Civ}\lambda1550$, the prominent broad and asymmetric $\text{HeII}\lambda1640$ line profile, and the P Cygni signature of $\text{Niv}\lambda1720$ clearly stand out. All these lines are associated with young, hot, and (very) massive stars.

We report the detection of $\text{HeII}\lambda1640$ emission with a measured rest-frame EW = $3.0 \pm 0.3$ Å and FWHM = $8.8 \pm 1.7$ Å (~1610 ± 300 km s⁻¹). The broad shape of the $\text{HeII}\lambda1640$ 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 $\text{HeII}\lambda1640$ line is consistent with the profile predicted by models and synthetic spectra (see Martins & Palacios 2022; Gräfener & Vink 2015). Broad $\text{HeII}\lambda1640$ 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 $\text{HeII}\lambda1640$ 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_\odot$ are crucial for producing the $\text{HeII}\lambda1640$ emission line (Gräfener & Vink 2015; Crowther et al. 2016). On the other
hand, Martins & Palacios (2022) have shown that HeIIλ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 NIVλ1486 and NIVλ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 NIVλ1486 and EW = 0.15 Å and FWHM ~ 2 Å for NIVλ1720. The helium enrichment also contributes to the strength of HeIIλ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 HeIIλ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<sub>Hα</sub> < 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 (Hα/(1 - f<sub>esc</sub><sup>abs</sup>)), the resulting EW increases to EW<sub>Hα</sub> ~ 1231 Å. We adopted f<sub>esc</sub><sup>abs</sup> values from Rivera-Thorsen et al. (2019), where the corresponding absolute escape fraction of LyC photons along the line of sight is f<sub>esc</sub><sup>abs</sup> = 32<sup>+2</sup><sub>-4</sub>% . 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 NIVλ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 NIVλ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 NIVλ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 HeIIλ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 Ovλ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 Ovλ1371 disappears as VMS evolve to a lower effective temperature. In their Fig. 4, no sign of Ovλ1371 exists after ~1 Myr. Together with NIVλ1486, this places a rather tight constraint on the cluster age.

## 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<sub>☉</sub>) with different VMS masses (150 M<sub>☉</sub> and 200 M<sub>☉</sub>). 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<sub>eff</sub>) and less on an exact wind prescription (e.g., Martins & Palacios 2022). Furthermore, because winds can strongly affect stellar evolution and thus T<sub>eff</sub>, 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<sub>☉</sub> 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<sup>6</sup> M<sub>☉</sub>. 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<sub>☉</sub> 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<sub>☉</sub> upper-mass limit within a few mass bins given by Eq. (1) from the BPASS manual<fnref n="1" />. Equation (1) gives the number of massive stars in the mass range [M<sub>a</sub>; M<sub>b</sub>],

$$N(M_a; M_b) = C \times M_1^{\alpha_1} \int_{M_a}^{M_b} M^{\alpha_2} dM. \quad (1)$$

Here, C is a constant and has the value of 1.23 × 10<sup>5</sup> for an arbitrary burst mass of 10<sup>6</sup> M<sub>☉</sub>. Moreover, M<sub>1</sub> = 0.5 M<sub>☉</sub> α<sub>1</sub> = -1.3, and α<sub>2</sub> = -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<sub>☉</sub> and 200 M<sub>☉</sub>. In this manner, we computed SEDs including VMS with IMFs extending up to 175 and 225 M<sub>☉</sub>, respectively, following Martins & Palacios (2022). In the mass bin from 100 M<sub>☉</sub> to 175 M<sub>☉</sub>, the extrapolation of the Salpeter IMF using Eq. (1) gives ~237 VMS, and for the mass bin from 175 M<sub>☉</sub> to 225 M<sub>☉</sub>, it provides us with ~60 VMS. We added SEDs of 150 M<sub>☉</sub> VMS as an average of the VMS population in the mass bin 100 M<sub>☉</sub> ≤ M ≤ 175 M<sub>☉</sub>, and we did the same by adding SEDs of 200 M<sub>☉</sub> VMS for the mass bin 175 M<sub>☉</sub> ≤ M ≤ 225 M<sub>☉</sub>. Finally, these numbers have to be rescaled to the appropriate cluster mass, which is ~10<sup>7</sup> M<sub>☉</sub> (Vanzella et al. 2022a).

The YMC shows significant NIVλ1486 emission, as shown in Fig. 2. From the VMS models and synthetic spectra, NIVλ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

<footnote n="1">¹ https://flexiblelearning.auckland.ac.nz/bpass/9.html</footnote>
<figure><img src="image_4.png" /><figcaption>**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 (Nivλ1486, Heiiλ1640, and Nivλ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 Nvλ1240 P Cygni line, which is characteristic of very young stellar populations.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**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⊙.</figcaption></figure>

Figs. 3, 4, and A.1). However, the Heiiλ1640 and Nivλ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 Heiiλ1640 and Nivλ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 include 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_\odot$, but with an IMF slope steeper than for a Salpeter IMF ($\alpha_2 < -2.35$) for $M > 100 M_\odot$.

### 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.1I YMC

We now address the morphological properties of image 5.11, which is the most strongly magnified of the multiple images of the star cluster ($\mu_{tot} \approx 76$; Pignataro et al. 2021). Image 5.11 is the brightest image and was detected with a high S/N in the F275W ($S/N \sim 90$) and F555W ($S/N \gg 100$), allowing us to investigate and compare the morphology in these two spectral regions: the emitting LyC ($\lambda < 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.11, 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 ($\mu_{tang} \approx 57$). We perform a relative comparison between images to ensure that the results do not depend on the magnification values.

As a first step, we computed a realistic model of image 5.11 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.11 into the F275W image in five different

<figure><img src="image_6.png" /></figure>

**Fig. 5.** Results from Galfit modeling after using a two-component fit. From the left, the first panel shows source 5.11 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.

<figure><img src="image_7.png" /></figure>

**Fig. 6.** HST F275W image showing source 5.11 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.

positions around source 5.11 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.11 (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.11. 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.11) appears nucleated, suggesting that the emitting LyC region is smaller than the region at 1700 Å.

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).

The same procedure was applied to the LyC-emitting source 5.11, 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.11 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.11, 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.11 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 \ll 0.05$. This means that the null hypothesis is not satisfied with the LyC profile of source 5.11, and it deviates from the CoG of a point-like source. Furthermore, we also find that the half-light size of 5.11 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.11 at 1700 Å of $\approx 7.8 \pm 1.4$ pc after delensing<fnref n="2" />. This agrees with Vanzella et al. (2020a), while the LyC image (5.11) has a smaller radius, $R_{\text{eff}} \approx 4.7 \pm 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

<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Three curves of growth normalized to 1. The cyan CoG corresponds to source 5.11 leaking LyC that is located in the Sunburst arc (F275W band). The violet CoG is the PSF-convolved best-fit model of source 5.11 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\sigma$. Vertical lines in the bottom left part of the panel mark the pixel radii within which 50% of the light is enclosed. The colors correspond to the CoG colors.</figcaption></figure>

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_{\text{relax}}$. When $t_{\text{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_{\text{relax}} \approx (0.065 \times \sigma^3)/(G^2 m \rho \ln \Lambda)$ (Spitzer 1987), where $m$ is the average stellar mass in the cluster (the mass of the cluster divided by the total number of stars), $\sigma$ is the central three-dimensional velocity dispersion, $\rho$ is the central mass density, $G$ is the gravitational constant, and $\ln \Lambda$ is the Coulomb logarithm, a number typically between 6 and 15. To compute the central relaxation time of source 5.11 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 $\sigma = \sqrt{GM/(6r_{0.5})}$. We adopted a mass of $M = 10^7 M_{\odot}$, 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 $\approx 30 \text{ km s}^{-1}$. 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.11 YMC is $N \sim 9 \times 10^6$. Combined with the total cluster mass, this gives an average stellar mass of $m \sim 1.1 M_{\odot}$. The central density in a Plummer model is $\rho = 3M/(4\pi r_{0.5}^3)$, and we adopted $\ln \Lambda = 15$. The resulting $t_{\text{relax}} \sim 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 $\sigma^3$ and $\rho$ is lowest. The same measure performed at any larger distance from the cluster center would have given an even longer $t_{\text{relax}}$.

<footnote n="2">² Adopting the pixel scale of 0.03"/pixel, 8200 pc per arcsecond at $z = 2.37$ and $\mu_{\text{tang}} \approx 57$, $R_{\text{eff}} = 0.03 \times 8200 \times ((2.6^2 - 1.9^2)^{0.5})/57$, adopting the same uncertainty on $\mu_{\text{tang}}$ as reported by Vanzella et al. (2022a).</footnote>
**5. Summary and conclusions**

We have presented a detailed spectroscopic and morphologi-
cal 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 mod-
els including VMS (Martins & Palacios 2022) to conduct exten-
sive 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) HeIIλ1640 emission with an
EW ≈ 3 Å in the rest frame and a line width of 1610 km s⁻¹,
and NIVλ1486 with an EW ≈ 0.2 Å in emission. In addition,
the P Cygni profile of NIVλ1720 (along with NV and CIV) is
also observed. All these features suggest the presence of very
massive (>100 M⊙) stars. The absence of Ovλ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) indi-
cates 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 HeIIλ1640 emission line is entirely due
to VMS. 2) A comparison of the observations with the mod-
els 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 emerg-
ing 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 anal-
ysis of the most magnified image of the star cluster shows that
the region emitting LyC is not point-like, with a light profile dif-
ferent 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 proba-
bility *p* ≪ 0.05) than the nonionizing region (1700 Å), suggest-
ing 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 angu-
lar half-light radii, the LyC region appears barely resolved with
R<sub>eff</sub> ≈ 4.7 ± 1.5 pc, while at 1700 Å, it is R<sub>eff</sub> ≈ 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 an ele-
ment that allowed carving the ionizing channel and the devel-
opment 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 sys-
tems for investigations of the most exotic phenomena in the
Universe.

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**

Abel, T., Bryan, G. L., & Norman, M. L. 2002, *Science*, **295**, 93

Bastian, N., Saglia, R. P., Goudfrooij, P., et al. 2006, A&A, 448, 881

Bestenlehner, J. M. 2020, MNRAS, 493, 3938

Binney, J., & Tremaine, S. 1987, *Galactic Dynamics* (Princeton, N.J.: Princeton University Press)

Brands, S. A., de Koter, A., Bestenlehner, J. M., et al. 2022, A&A, 663, A36

Brinchmann, J., Kunth, D., & Durret, F. 2008, A&A, 485, 657

Calzetti, D., Johnson, K. E., Adamo, A., et al. 2015, ApJ, 811, 75

Chisholm, J., Rigby, J. R., Bayliss, M., et al. 2019, ApJ, 882, 182

Cignoni, M., Sabbi, E., van der Marel, R. P., et al. 2015, ApJ, 811, 76

Crowther, P. A., Schnurr, O., Hirschi, R., et al. 2010, MNRAS, 408, 731

Crowther, P. A., Caballero-Nieves, S. M., Bostroem, K. A., et al. 2016, MNRAS, 458, 624

Dahle, H., Aghanim, N., Guennou, L., et al. 2016, A&A, 590, L4

Dib, S., Shadmehri, M., Padoan, P., et al. 2010, MNRAS, 405, 401

Diego, J. M., Pascale, M., Kavanagh, B. J., et al. 2022, A&A, 665, A134

Doran, E. I., Crowther, P. A., de Koter, A., et al. 2013, A&A, 558, A134

Eldridge, J. J., & Stanway, E. R. 2012, MNRAS, 419, 479

Eldridge, J. J., Stanway, E. R., Xiao, L., et al. 2017, PASA, 34

Goswami, S., Slemer, A., Marigo, P., et al. 2021, A&A, 650, A203

Gräfener, G. 2021, A&A, 647, A13

Gräfener, G., & Vink, J. S. 2015, A&A, 578, L2

Habibi, M., Stolte, A., Brandner, W., Hußmann, B., & Motohara, K. 2013, A&A, 556, A26

Harayama, Y., Eisenhauer, F., & Martins, F. 2008, ApJ, 675, 1319

Johnson, T. L., Rigby, J. R., Sharon, K., et al. 2017, ApJ, 843, L21

Khorrami, Z., Langlois, M., Clark, P. C., et al. 2021, MNRAS, 503, 292

Köhler, K., Langer, N., de Koter, A., et al. 2015, A&A, 573, A71

Leitherer, C., Ekström, S., Meynet, G., et al. 2014, ApJS, 212, 14

Leitherer, C., Byler, N., Lee, J. C., & Levesque, E. M. 2018, ApJ, 865, 55

Mainali, R., Rigby, J. R., Chisholm, J., et al. 2022, ApJ, 940, 160

Marino, R. A., Rosales-Ortega, F. F., Sánchez, S. F., et al. 2013, A&A, 559, A114

Martins, F., & Palacios, A. 2022, A&A, 659, A163

Martins, F., Hillier, D. J., Paumard, T., et al. 2008, A&A, 478, 219

Meštrić, U., Vanzella, E., Zanella, A., et al. 2022, MNRAS, 516, 3532

Östlin, G., Cumming, R. J., & Bergvall, N. 2007, A&A, 461, 471

Peng, C. Y., Ho, L. C., Impey, C. D., & Rix, H.-W. 2010, AJ, 139, 2097

Pignataro, G. V., Bergamini, P., Meneghetti, M., et al. 2021, A&A, 655, A81

Plummer, H. C. 1911, MNRAS, 71, 460

Portegies Zwart, S. F., McMillan, S. L. W., & Gieles, M. 2010, ARA&A, 48, 431

Rigby, J. R., Johnson, T. L., Sharon, K., et al. 2017, ApJ, 843, 79

Rigby, J. R., Bayliss, M. B., Chisholm, J., et al. 2018a, ApJ, 853, 87

Rigby, J. R., Bayliss, M. B., Sharon, K., et al. 2018b, AJ, 155, 104

Rivera-Thorsen, T. E., Dahle, H., Gronke, M., et al. 2017, A&A, 608, L4

Rivera-Thorsen, T. E., Dahle, H., Chisholm, J., et al. 2019, Science, 366, 738

Rix, S. A., Pettini, M., Leitherer, C., et al. 2004, ApJ, 615, 98

Sabbahit, G. N., Vink, J. S., Higgins, E. R., & Sander, A. A. C. 2022, MNRAS, 514, 3736

Schaerer, D., & Stasińska, G. 1999, A&A, 345, L17

Senchyna, P., Stark, D. P., Charlot, S., et al. 2021, MNRAS, 503, 6112

Sharon, K., Mahler, G., Rivera-Thorsen, T. E., et al. 2022, ApJ, 941, 203

Smith, L. J., Crowther, P. A., Calzetti, D., & Sidoli, F. 2016, ApJ, 823, 38

Smith, L. J., Bajaj, V., Ryon, J., & Sabbi, E. 2020, ApJ, 896, 84

Spitzer, L. 1987, *Dynamical Evolution of Globular Clusters* (Princeton, N.J.: Princeton University Press)

Stanway, E. R., & Eldridge, J. J. 2018, MNRAS, 479, 75

Tody, D. 1986, in *Instrumentation in astronomy VI*, ed. D. L. Crawford, Proc. SPIE, 627, 733

Vanzella, E., De Barros, S., Cupani, G., et al. 2016, ApJ, 821, L27

Vanzella, E., Calura, F., Meneghetti, M., et al. 2017, MNRAS, 467, 4304

Vanzella, E., Caminha, G. B., Calura, F., et al. 2020a, MNRAS, 491, 1093

Vanzella, E., Meneghetti, M., Pastorello, A., et al. 2020b, MNRAS, 499, L67

Vanzella, E., Caminha, G. B., Rosati, P., et al. 2021, A&A, 646, A57

Vanzella, E., Castellano, M., Bergamini, P., et al. 2022a, ApJ, 940, L53

Vanzella, E., Castellano, M., Bergamini, P., et al. 2022b, A&A, 659, A2

Vink, J. S., Muijres, L. E., Anthonisse, B., et al. 2011, A&A, 531, A132

Vink, J. S., Heger, A., Krumholz, M. R., et al. 2015, Highlights Astron., 16, 51

Wofford, A., Leitherer, C., Chandar, R., & Bouret, J.-C. 2014, ApJ, 781, 122

Yusof, N., Hirschi, R., Meynet, G., et al. 2013, MNRAS, 433, 1114

*Acknowledgements.* We acknowledge financial support through grants PRIN-MIUR 2017WSCC32, 2020SKSTHZ and the INAF GO Grant 2022 "The
## 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 NIVL1486 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_\odot$ and added 236.56 stars in the 100 - 175 $M_\odot$ mass range and 60.30 stars in the mass range 175 - 225 $M_\odot$. 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).

<figure><img src="image_9.png" /><figcaption>**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_\odot$, ages of 1.5 Myr, 2 Myr, and 2.5 Myr and an added number of 236.56 $150 M_\odot$ VMS. The right column shows a BPASS model with an IMF up to 100 $M_\odot$ at same ages (as shown in previous column), but with an added number of 236.56 and 60.30 VMS of $150 M_\odot$ and $200 M_\odot$.</figcaption></figure>