# Evidence for very massive stars in extremely UV-bright star-forming galaxies at $z \sim 2.2-3.6$

A Upadhyaya, R Marques-Chaves, D Schaerer, Fabrice Martins, I

Pérez-Fournon, Ana Palacios, E.R Stanway

► To cite this version:

A Upadhyaya, R Marques-Chaves, D Schaerer, Fabrice Martins, I Pérez-Fournon, et al. Evidence for very massive stars in extremely UV-bright star-forming galaxies at $z \sim 2.2-3.6$. Astronomy & Astrophysics - A&A, 2024, 686, pp.A185. 10.1051/0004-6361/202449184. hal-04659639

HAL Id: hal-04659639

https://hal.science/hal-04659639v1

Submitted on 24 Jul 2024

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
# Evidence for very massive stars in extremely UV-bright star-forming galaxies at $z \sim 2.2-3.6^*$

A. Upadhyaya<fnref n="1" />,<fnref n="2" />, R. Marques-Chaves<fnref n="1" />, D. Schaerer<fnref n="1" />,<fnref n="3" />, F. Martins<fnref n="4" />, I. Pérez-Fournon<fnref n="5" />,<fnref n="6" />, A. Palacios<fnref n="4" />, and E. R. Stanway<fnref n="2" />

Received 8 January 2024 / Accepted 21 March 2024

## ABSTRACT

We present a comprehensive analysis of the presence of very massive stars (VMS $> 100 M_\odot$) in the integrated spectra of 13 UV-bright star-forming galaxies at $2.2 \le z \le 3.6$ taken with the Gran Telescopio Canarias (GTC). These galaxies have very high UV absolute magnitudes ($M_{UV} \approx -24$), intense star formation (star formation rate $\approx 100-1000 M_\odot \text{ yr}^{-1}$), and metallicities in the range of $12 + \log(O/H) \approx 8.10-8.50$ inferred from strong rest-optical lines. The GTC rest-UV spectra reveal spectral features indicative of very young stellar populations with VMS, such as strong P-Cygni line profiles in the wind lines N V $\lambda$1240 and C IV $\lambda$1550 along with intense and broad He II $\lambda$1640 emission with equivalent width ($EW_0 \approx 1.40-4.60 \AA$, and full width half maximum (FWHM) $\approx 1150-3170 \text{ km s}^{-1}$). A Comparison with known VMS-dominated sources and typical galaxies without VMS reveals that some UV-bright galaxies closely resemble VMS-dominated clusters (e.g., R136 cluster). The presence of VMS is further supported by a quantitative comparison of the observed strength of the He II emission with population synthesis models with and without VMS, where models with VMS are clearly preferred. Employing an empirical threshold for $EW_0$ (He II) $\ge 3.0 \AA$, along with the detection of other VMS-related spectral profiles (N IV $\lambda$1486, 1719), we classify nine out of 13 UV-bright galaxies as VMS-dominated sources. This high incidence of VMS-dominated sources in the UV-bright galaxy population ($\approx 70\%$) contrasts significantly with the negligible presence of VMS in typical $L_{UV}^*$ LBGs at similar redshifts ($<1\%$). Our results thus indicate that VMS are common in UV-bright galaxies, suggesting a different initial mass function (IMF) with upper mass limits between 175 $M_\odot$ and 475 $M_\odot$.

**Key words.** stars: massive – galaxies: high-redshift – galaxies: starburst – ultraviolet: galaxies

## 1. Introduction

The epoch of cosmic noon at a redshift around two to three plays an important role in the evolution of our Universe. It is the time when galaxies underwent intense star formation with a peak in the star formation rate density (Lilly et al. 1996; Feulner et al. 2005; Madau & Dickinson 2014; López Fernández et al. 2018; Sánchez et al. 2019; Koushan et al. 2021).

Typical star-forming galaxies at cosmic noon show strong Ly$\alpha$ $\lambda$1216 emission along with stellar wind features such as N V $\lambda$1240 and C IV $\lambda$1550 P-Cygni line profiles in their rest-frame UV spectra (e.g., Shapley et al. 2003; Berry et al. 2012; Steidel et al. 2016; Le Fèvre et al. 2019; Marques-Chaves et al. 2020a). Due to their faintness, with UV absolute magnitudes of $M_{UV} \ge -21$, the detailed characterization of their rest-frame UV has been achieved using stacking techniques of hundreds to thousands of individual spectra (e.g., Shapley et al. 2003; Le Fèvre et al. 2015; Steidel et al. 2016), very deep spectroscopy (e.g., Pentericci et al. 2018; Garilli et al. 2021),

or with the help of gravitational lensing (e.g., Pettini et al. 2000; Quider et al. 2009; Dessauges-Zavadsky et al. 2010; Marques-Chaves et al. 2017, 2020a). Even with many extensive surveys before the launch of JWST, only a small number of star-forming galaxies brighter than $M_{UV} \le -22$ were discovered at this cosmic epoch (e.g., Bian et al. 2012).

Recently, Marques-Chaves et al. (2020b, 2021, 2022) found star-forming galaxies with $M_{UV} \sim -24$. Their brightness and very small number density suggest an extreme star-formation rate (SFR) which consumes the gas rapidly. Additionally, these UV bright galaxies show intense and broad He II $\lambda$1640 emission with rest-frame equivalent widths ($EW_0 \ge 3 \AA$ and full width half maximum (FWHM) of $\approx 1000-3000 \text{ km s}^{-1}$, which is not common for star-forming galaxies at this redshift. Last, but not least, some of these galaxies have been found to be strong emitters of Lyman continuum radiation, posing new questions on the main contributors of cosmic reionization (see, Marques-Chaves et al. 2021).

Determining the stellar content of these UV-bright star-forming galaxies is crucial to understanding the origin and evolution of these sources. The broad He II $\lambda$1640 emission provides some clues on the stellar populations that dominate the rest-frame UV spectra of these galaxies. The R136

<footnote n="1">¹ Geneva Observatory, Department of Astronomy, University of Geneva, Chemin Pegasi 51, 1290 Versoix, Switzerland
e-mail: ankurupadhyaya7@gmail.com</footnote>

<footnote n="2">² Department of Physics, University of Warwick, Gibbet Hill Road, Coventry CV4 7AL, UK</footnote>

<footnote n="3">³ CNRS, IRAP, 14 Avenue E. Belin, 31400 Toulouse, France</footnote>

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

<footnote n="5">⁵ Instituto de Astrofísica de Canarias, C/Vía Láctea, s/n, 38205 San Cristóbal de La Laguna, Tenerife, Spain</footnote>

<footnote n="6">⁶ Universidad de La Laguna, Dpto. Astrofísica, 38206 San Cristóbal de La Laguna, Tenerife, Spain</footnote>

<footnote>* A copy of the reduced spectra is available at the CDS
via anonymous ftp to cdsarc.cds.unistra.fr (130.79.128.5)
or via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/
686/A185</footnote>
cluster in the Large Magellanic Cloud (LMC) shows a sim-
ilar broad He II λ1640 profile (EW₀ ~ 4.5 Å, FWHM ≈
1800 km s⁻<fnref n="1" />; Crowther et al. 2010, 2016). The R136 cluster has
been confirmed to host very massive stars (VMS > 100 M⊙)
(Crowther et al. 2010, 2016; Bestenlehner et al. 2011, 2020;
Hainich et al. 2014; Brands et al. 2022) with the most recent
studies showing its most massive star to have a mass of at least
200 M⊙ (Kalari et al. 2022; Shenar et al. 2023).

There have been other studies in the local Universe where the presence or potential presence of VMS has been detected (e.g., Massey & Hunter 1998; Crowther & Dessart 1998; Bruhweiler et al. 2003; Martins et al. 2008). Using *Hubble* Space Telescope (HST) UV spectroscopy, Wofford et al. (2014) studied the super star cluster A1 in NGC 3125 and inferred the presence of VMS; a further study strengthened this interpretation (Wofford et al. 2023). Smith et al. (2016) inferred the presence of VMS in one of the nuclear star clusters in the compact H II region of the blue dwarf galaxy NGC 5253 using HST UV and Very Large Telescope (VLT) optical spectroscopy. Senchyna et al. (2021) analyzed HST UV spectra of a few nearby star-forming regions showing Wolf-Rayet (WR) features in the Sloan Digital Sky Survey (SDSS) spectrum and argued that VMS may be present in some of them. Smith et al. (2023) found evidence for VMS in a super star cluster in the metal-poor galaxy Mrk 71. Recently, Meštrić et al. (2023) noted the presence of VMS in a gravitationally lensed star cluster, the Sunburst cluster at z = 2.37. The recent study of Martins et al. (2023) not only infers the presence of VMS in a few nearby star-forming regions but also provides guidance on how to separate VMS from WR sources using both UV and optical spectroscopy. All of these studies suggest that sources hosting VMS have intense and broad He II λ1640 emission with a relatively high EW₀, typically EW (He II) ≥ 2–3 Å and FWHM ≥ 1000 km s⁻<fnref n="1" />.

So far, little is known about the presence and occurrence of VMS in high-redshift galaxies. While broad He II λ1640 emission is a fairly common feature of distant star-forming galaxies, it is relatively weak (EW₀ ~ 1–2 Å, FWHM ≥ 1000 km s⁻<fnref n="1" />) and attributed to WR stars (see e.g., Shapley et al. 2003). Stronger and broad He II emission seems rare and has been noted for example in studies by Cassata et al. (2013), Nanayakkara et al. (2019), Saxena et al. (2020), and Wofford et al. (2023) in Lyman break galaxies (LBGs).

Here, we report the detection of VMS in a large fraction of extremely UV-bright star-forming galaxies at $z \sim 2.2-3.6$. Using observations taken with the Optical System for Imaging and low-Intermediate-Resolution Integrated Spectroscopy (OSIRIS) spectrograph at the Gran Telescopio Canarias (GTC), and the latest VMS models from Martins & Palacios (2022), we present a detailed analysis of the rest-UV spectra and the main spectral features of 13 UV-bright galaxies discovered by Marques-Chaves et al. (2020a, 2021, 2022). These objects consistently show very strong He II λ1640 emission and other stellar features, providing unique insight on their stellar populations and the high mass end of the initial mass function.

The paper is structured as follows. In Sect. 2, we describe our observations and how the data have been reduced. In Sect. 3, we describe the properties of the UV bright galaxies and present their He II λ1640 emission spectral profiles. In Sect. 4, we provide empirical evidence for the presence of VMS in our sources. In Sect. 5, we present the models and their analysis to complement the evidence for the presence of VMS in our sources. In Sect. 6, we discuss various implications of these results. Finally, Sect. 7 provides us with conclusions and summarizes the results of our work. Throughout this work, we

assume a concordance cosmology with Ω_m = 0.274, Ω_Λ = 0.726, and H_0 = 70 km s⁻<fnref n="1" /> Mpc⁻<fnref n="1" />. All magnitudes are given in the AB system.

## 2. Observations

### 2.1. Sample selection

The galaxies studied in this work are part of a large sample of ~70 UV-luminous star-forming galaxies at $z \simeq 2.0-3.6$ selected from the ~9000 deg<fnref n="2" />-wide extended Baryon Oscillation Spectroscopic Survey (eBOSS; Abolfathi et al. 2018) of the SDSS (Eisenstein et al. 2011). The sample, selection techniques, and overall properties will be presented in a separate work (Marques-Chaves, in prep.). Here, we explore the properties of 13 of these UV-bright sources that have deep follow-up spectroscopy. The properties of three of these sources were already analyzed in detail in Marques-Chaves et al. (2020b, 2021, 2022) and Álvarez-Márquez et al. (2021), while the remaining are new discoveries.

Table 1 shows several properties of these 13 sources, including coordinates, redshifts, and magnitudes. Overall, these sources are optically bright, $R \simeq 20.8-21.9$, corresponding to rest-frame UV absolute magnitudes from $M_{UV} = -23.3$ to $-24.6$. Given their brightness, they are ideal targets for deep, high signal-to-noise (S/N) spectroscopy.

### 2.2. Optical spectroscopy

Optical spectra were obtained with the OSIRIS<fnref n="1" /> instrument on the GTC between 2017 and 2021 (Table 1) under the GTC programs IDs: GTC67-17A, GTCMULTIPLE2F-18A, GTC50-18A, GTCMULTIPLE2E-18B, GTCMULTIPLE2E-19A, GTC21-20A, and GTC29-21A (PI: R. Marques-Chaves). Observations were performed under good seeing conditions (≤1.2", FWHM) using the low-resolution grism R1000B, which provides a spectral resolution $R \sim 700$ and coverage of 3600–7600 Å. Long-slits with 1.2"-width were centered on each object and oriented with the parallactic angle. Total on-source exposure times vary from 3000 to 9000 s, depending on the source (Table 1), and were split into at least four individual exposures for cosmic rays removal.

Data were reduced following standard reduction procedures using IRAF, starting from the subtraction of the bias and flat-field correction. The wavelength calibration is done using HgAr+Ne+Xe arc-lamp data. 2D spectra are background subtracted using sky regions on both sides of the trace of each source. Individual 1D spectra are extracted, stacked, and corrected for the instrumental response using observations of standard stars observed each night. We use the extinction curve of Cardelli et al. (1989) and the extinction map of Schlafly & Finkbeiner (2011) to correct for the reddening effect in the Galaxy. Finally, we corrected for telluric absorption using the IRAF telluric routine. Individual 1D spectra are shown in Fig. A.1.

### 2.3. Near-IR spectroscopy

We also obtained near-IR spectra for nine sources with the Espectrógrafo Multiobjeto Infra-Rojo (EMIR)<fnref n="2" /> on the GTC between 2018 and 2019. For each source, we use the *HK* grism

<footnote n="1">¹ http://www.gtc.iac.es/instruments/osiris/</footnote>

<footnote n="2">² http://www.gtc.iac.es/instruments/emir/</footnote>
<table><caption>**Table 1.** Summary of the GTC observations.</caption><thead><tr><th rowspan="2">Name</th><th>RA</th><th>Dec</th><th>z</th><th>R-band</th><th colspan="2">Optical spectroscopy</th><th colspan="2">Near-IR spectroscopy</th></tr><tr><th>(J2000)</th><th>(J2000)</th><th></th><th>(AB)</th><th>Date<br>(dd/mm/yyyy)</th><th>Exp. time<br>(s)</th><th>Date<br>(dd/mm/yyyy)</th><th>Exp. time<br>(s)</th></tr></thead><tbody><tr><td>J0006+2452</td><td>00:06:44.73</td><td>+24:52:53.19</td><td>2.379</td><td>20.98</td><td>16/07/2020</td><td>3600</td><td>12/07/2019</td><td>2560</td></tr><tr><td>J0031+3545</td><td>00:31:12.43</td><td>+35:45:56.12</td><td>2.816</td><td>20.84</td><td>16/07/2020</td><td>3600</td><td>—</td><td>—</td></tr><tr><td>J0036+2725</td><td>00:36:06.10</td><td>+27:25:39.27</td><td>2.171</td><td>21.07</td><td>19/08/2020</td><td>3600</td><td>—</td><td>—</td></tr><tr><td>J0110–0501</td><td>01:10:45.58</td><td>−05:01:39.27</td><td>2.368</td><td>21.80</td><td>30/08/2019</td><td>3720</td><td>20/09/2018</td><td>2560</td></tr><tr><td>J0115+1837</td><td>01:15:21.95</td><td>+18:37:44.47</td><td>2.322</td><td>21.59</td><td>23/09/2019</td><td>3720</td><td>20/09/2018</td><td>2560</td></tr><tr><td>J0121+0025 <sup>(a)</sup></td><td>01:21:56.09</td><td>+00:25:20.30</td><td>3.246</td><td>21.60</td><td>18/08/2020</td><td>7200</td><td>—</td><td>—</td></tr><tr><td>J0146–0220</td><td>01:46:37.02</td><td>−02:20:55.86</td><td>2.160</td><td>21.28</td><td>10/08/2018</td><td>4500</td><td>19/09/2018</td><td>2560</td></tr><tr><td>J0850+1549</td><td>08:50:38.86</td><td>+15:49:17.88</td><td>2.424</td><td>21.43</td><td>12/03/2018</td><td>3600</td><td>25/04/2018</td><td>2560</td></tr><tr><td>J1013+4650</td><td>10:13:28.76</td><td>+46:50:43.47</td><td>2.286</td><td>21.40</td><td>07/11/2018</td><td>4500</td><td>26/12/2018</td><td>2560</td></tr><tr><td>J1157+0113</td><td>11:57:34.12</td><td>+01:13:08.21</td><td>2.545</td><td>21.88</td><td>01/05/2019</td><td>7200</td><td>—</td><td>—</td></tr><tr><td>J1220+0842 <sup>(b)</sup></td><td>12:20:40.72</td><td>+08:42:38.14</td><td>2.470</td><td>20.86</td><td>19/06/2017</td><td>3000</td><td>23/03/2018</td><td>2560</td></tr><tr><td>J1316+2614 <sup>(c)</sup></td><td>13:16:29.61</td><td>+26:14:07.05</td><td>3.612</td><td>21.26</td><td>20/04/2021</td><td>9000</td><td>08/05/2018</td><td>2560</td></tr><tr><td>J1335+4330</td><td>13:35:15.79</td><td>+43:30:35.07</td><td>2.170</td><td>21.34</td><td>11/04/2018</td><td>3600</td><td>25/04/2018</td><td>2560</td></tr></tbody></table>

Notes. – data already presented in (a) Marques-Chaves et al. (2021), (b) Marques-Chaves et al. (2020b), and (c) Marques-Chaves et al. (2022).

with a 0.8"-width, providing a spectral resolution $R \sim 700$ and coverage of 1.45–2.41 µm. The slits were centered using a bright reference star. Observations were taken with on-source exposure times of 16 × 160 s. with a standard 10" ABBA dither. Reduction of near-IR spectra was performed using the official EMIR pipeline<fnref n="3" />.

### 3. Properties of UV-bright galaxies

From the results of Marques-Chaves et al. (in prep.) the UV-bright galaxies show very steep UV slopes ($f_λ \propto λ^β$) between $β \approx -1.7$ and $β \approx -2.8$ with a mean and scatter of $\langle β \rangle = -2.33 \pm 0.30$. This translates into little dust obscuration, with a mean value of $E(B-V) = 0.05 \pm 0.04$ for this sample assuming the Calzetti et al. (2000) extinction law with an intrinsic $β_0 = -2.44$. SFRs are derived from the UV-luminosity, using the specific conversion factor $κ_{UV} = 1.3 \times 10^{-28} M_\odot \text{yr}^{-1}/(\text{erg s}^{-1} \text{Hz}^{-1})$ derived using Binary Population and Spectral Synthesis (BPASS) binary models (Eldridge et al. 2017; Stanway & Eldridge 2018; Byrne et al. 2022) that assumes the Chabrier (2003) initial mass function (IMF) with an upper mass cutoff of 100 $M_\odot$, the LMC metallicity, and a continuous star-formation over 10 Myr (the typical age observed for these sources)<fnref n="4" />. We correct the derived UV-SFRs for the dust attenuation using the observed UV slopes and assuming the Calzetti et al. (2000) attenuation law. The dust-corrected SFRs of our sources range over $\approx 100–900 M_\odot \text{yr}^{-1}$ and are listed in Table 2. The stellar masses of the young stellar population are derived assuming a continuous SFR over 10 Myr, and are also listed in Table 2. We refer to Marques-Chaves et al. (in prep.) for the description of the general properties of this sample. The properties of a few of these sources, namely J1220+0842, J0121+0025, and J1316+2614, were already investigated in detail in Marques-Chaves et al. (2020b, 2021, 2022), for which they find very young stellar populations ($\approx 10$ Myr) with SFR $\approx 400–600 M_\odot \text{yr}^{-1}$ and stellar masses of $\log(M_\star/M_\odot) \approx 9.6–9.8$. J0121+0025 and J1316+2614 are also found to be strong Lyman continuum (LyC) emitters, with LyC

escape fractions of $\approx 40\%–90\%$ (Marques-Chaves et al. 2021, 2022).

#### 3.1. Strong HeII λ1640 and other wind lines

The GTC spectra present a high S/N in the continuum, revealing a wealth of spectral features arising from different galaxy components. These include ISM absorption lines (e.g., Si II λ1260, C II λ1334), photospheric absorption lines (e.g., S V λ1501), or nebular gas for some sources (O III] λ1666, C III] λ1908).

In particular, and most relevant for the present work, is the detection of intense emission in He II λ1640 for several sources in our sample. Figure 1 shows the He II profiles observed in the spectra of our sources. We measure the strength of He II emission in our sources by measuring the rest-frame equivalent width ($EW_0$) using an integration region from 1630 Å to 1655 Å (rest; blue in Fig. 1), and two pseudo-continuum regions on both sides of He II for the continuum estimation (1610–1629 Å and 1675–1690 Å, respectively; red in Fig. 1). The continuum regions were selected to avoid the contribution of the low-ionization ISM lines Fe II λ1608 and Al II λ1670 and the nebular emission from [O III] λ1661, 1666. We measure $EW_0$ (He II) between $\approx 1.4$ Å and $\approx 4.6$ Å for our sources (see Table 2). In particular, J0006+2452, J0110–0501, and J1316+2614 are the most extreme He II emitters in our sample showing $EW_0$ (He II) $\approx 4.5$ Å, which is similar to that measured in local star-forming regions with VMS (e.g., R136 star cluster with $EW_0$ (He II) $\approx 4.5$ Å; Crowther et al. 2016).

We also fit Gaussian profiles to the continuum-normalized He II emission. The best-fits are shown in green in Fig. 1. The He II emission appears broad for all sources with $FWHM \approx 1000–3000$ km s$^{-1}$. The He II line appears symmetric for some sources (e.g., J0006+2452 or J0115+1837), while for others the line appears asymmetric (J1220+0842 or J1335+4330) or with multiple emission/absorption peaks (e.g., J0031+3545 or J0121+0025).

The GTC spectra also reveal strong P-Cygni line profiles in the wind lines N V λ1240 and C IV λ1550, but also N IV λ1486, and N IV λ1719 for some sources (see Appendix A). These profiles are the result of strong stellar winds from O-type stars and indicate very young ages ($\lesssim 10$ Myr) of the stellar population

<footnote n="3">³ https://pyemir.readthedocs.io/en/latest/index.html</footnote>

<footnote n="4">⁴ Note that this value of $κ_{UV}$ implies SFR values higher by a factor $\sim 1.7$ than the classical assumption of SFR = const over $\sim 100$ Myr.</footnote>
<table><caption>**Table 2.** Global properties and measurements of the strengths of stellar He II λ1640 emission, N v λ1240, and C IV λ1550 absorption.</caption><thead><tr><th>Name</th><th><i>M</i><sub>UV</sub><br>(AB)</th><th>SFR<br>(<i>M</i><sub>☉</sub> yr<sup>-1</sup>)</th><th><i>M</i><sub>★</sub><sup>young</sup><br>(log [<i>M</i><sub>☉</sub>])</th><th>12 + log(O/H)</th><th>EW<sub>0</sub> (N v)<br>(Å)</th><th>EW<sub>0</sub> (C IV)<br>(Å)</th><th>EW<sub>0</sub> (He II)<br>(Å)</th><th>FWHM (He II)<br>(km s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>J0006+2452</td><td>-24.17</td><td>316 ± 31</td><td>9.50 ± 0.05</td><td>8.49 ± 0.10</td><td>-3.78 ± 0.40</td><td>-6.70 ± 0.19</td><td>4.54 ± 0.38</td><td>2070 ± 180</td></tr><tr><td>J0031+3545</td><td>-24.62</td><td>390 ± 60</td><td>9.59 ± 0.07</td><td>—</td><td>-2.11 ± 0.48</td><td>-5.40 ± 0.34</td><td>2.29 ± 0.58</td><td>1800 ± 430</td></tr><tr><td>J0036+2725</td><td>-23.91</td><td>612 ± 62</td><td>9.79 ± 0.05</td><td>—</td><td>-5.40 ± 0.34</td><td>-6.29 ± 0.19</td><td>2.69 ± 0.35</td><td>2810 ± 470</td></tr><tr><td>J0110-0501</td><td>-23.33</td><td>110 ± 15</td><td>9.04 ± 0.06</td><td>8.48 ± 0.15</td><td>-4.87 ± 0.64</td><td>-4.32 ± 0.35</td><td>4.57 ± 0.59</td><td>1990 ± 510</td></tr><tr><td>J0115+1837</td><td>-23.52</td><td>210 ± 35</td><td>9.33 ± 0.07</td><td>8.37 ± 0.12</td><td>-2.08 ± 0.76</td><td>-4.22 ± 0.29</td><td>3.80 ± 0.47</td><td>2870 ± 165</td></tr><tr><td>J0121+0025</td><td>-24.11</td><td>460 ± 90</td><td>9.65 ± 0.08</td><td>—</td><td>-3.54 ± 0.20</td><td>-5.51 ± 0.20</td><td>3.06 ± 0.33</td><td>3170 ± 580</td></tr><tr><td>J0146-0220</td><td>-23.68</td><td>340 ± 43</td><td>9.53 ± 0.06</td><td>8.47 ± 0.20</td><td>-3.38 ± 0.46</td><td>-3.87 ± 0.27</td><td>2.99 ± 0.42</td><td>2480 ± 135</td></tr><tr><td>J0850+1549</td><td>-23.76</td><td>190 ± 26</td><td>9.28 ± 0.06</td><td>8.27 ± 0.13</td><td>-1.18 ± 0.46</td><td>-2.27 ± 0.31</td><td>1.37 ± 0.44</td><td>2030 ± 1351</td></tr><tr><td>J1013+4650</td><td>-23.67</td><td>210 ± 27</td><td>9.32 ± 0.05</td><td>8.45 ± 0.10</td><td>-2.64 ± 0.53</td><td>-6.16 ± 0.27</td><td>3.45 ± 0.46</td><td>2520 ± 310</td></tr><tr><td>J1157+0113</td><td>-23.40</td><td>195 ± 68</td><td>9.29 ± 0.13</td><td>—</td><td>-2.32 ± 1.12</td><td>-2.41 ± 0.62</td><td>3.02 ± 0.92</td><td>2280 ± 260</td></tr><tr><td>J1220+0842</td><td>-24.36</td><td>360 ± 32</td><td>9.56 ± 0.04</td><td>8.13 ± 0.19</td><td>-3.16 ± 0.29</td><td>-2.00 ± 0.21</td><td>2.52 ± 0.33</td><td>1150 ± 20</td></tr><tr><td>J1316+2614</td><td>-24.65</td><td>415 ± 115</td><td>9.62 ± 0.11</td><td>8.45 ± 0.12</td><td>-1.87 ± 0.68</td><td>-3.82 ± 0.54</td><td>4.48 ± 0.14</td><td>2370 ± 85</td></tr><tr><td>J1335+4330</td><td>-23.64</td><td>920 ± 150</td><td>9.96 ± 0.08</td><td>8.33 ± 0.17</td><td>-2.70 ± 0.88</td><td>-3.64 ± 0.34</td><td>3.40 ± 0.51</td><td>2350 ± 200</td></tr></tbody></table>

**Notes.** The He II λ1640 emission EW<sub>0</sub> is measured by integrating normalized flux in the spectral regions from 1630 Å to 1655 Å. The N v λ1240 and C IV λ1550 absorption EW<sub>0</sub> are measured by integrating normalized flux in the spectral regions from 1230 Å to 1240 Å, and 1530 Å to 1543 Å respectively.

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

**Fig. 1.** Normalized GTC spectra of the 13 UV-bright galaxies studied in this work around the He II λ1640 line (black, and 1σ uncertainty in gray). Spectra are normalized using the continuum regions marked in red. Equivalent widths of the He II line are measured by integrating the spectral regions from 1630 Å to 1655 Å (blue regions). Green dashed lines show Gaussian fits of the He II λ1640 emission. The x-axes represent the rest-frame and observed wavelengths (bottom and top axis, respectively).

(e.g., Leitherer et al. 2011). Well-defined P-Cygni line profiles in Si IV λ1400, which are characteristic of O-type supergiants (e.g., Walborn et al. 1985; Garcia & Bianchi 2004), are also detected in some sources (e.g., J0006+2452, J0110-0501, J1316+2614).

## 3.2. Metallicity

An important property worth discussing is the metallicity of these sources, which governs the mass-loss rate of massive stars and thus the shape and strength of wind lines. Strong

P-Cygni profiles are detected in the GTC spectra and, in particular, in C IV which is known to be the most sensitive feature to the metallicity (Chisholm et al. 2019). This suggests that the metal content in these sources cannot be extremely low (hence $Z/Z_{\odot} \ge 0.1$). However, an accurate determination of the stellar metallicity using wind lines is rather difficult given the natural degeneracy between metallicity, age, and the shape of the IMF (slope and upper mass limit). Similarly, metallicity indicators using the strength of photospheric absorption features (e.g., Rix et al. 2004; Sommariva et al. 2012; Calabrò et al. 2021) may
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Relationship between the strength of He II λ1640 (EW in Å) and the oxygen abundance of our sources. The LMC metallicity is marked with the vertical dashed line.</figcaption></figure>

not work properly for these sources as well since they were calibrated using stellar models assuming standard IMF and ages of 100 Myr.

To overcome this, we perform a sanity check on the nebular metallicity using rest-frame optical lines. Specifically, we extract the flux of $H\alpha$ and [N II] λ6585 emission lines using Gaussian profiles and relate the observed $N_2 = \log([\text{N II}]/H\alpha)$ line ratio with metallicity using the strong-line calibrator of Marino et al. (2013). For the seven sources with observations of $H\alpha$ and [N II], we measure values of $N_2$ between -1.02 and -0.55, yielding $12 + \log(O/H) = 8.27-8.49$. We also used the metallicity measurements of J1220+0842 ($12 + \log(O/H) = 8.13 \pm 0.19$) and J1316+2614 ($12 + \log(O/H) = 8.45 \pm 0.12$) obtained in Marques-Chaves et al. (2020b, 2022), respectively, using the R23 metallicity calibrator. The derived metallicities for the nine sources have mean value and scatter of $12 + \log(O/H) = 8.36 \pm 0.11$, that is, $Z/Z_{\odot} \approx 0.5$ (taking solar value of $12 + \log(O/H) = 8.69$), and are listed in Table 2. Figure 2 shows the relationship between the strength of the He II λ1640 emission and $12 + \log(O/H)$. We find that sources with higher metallicities tend to have stronger He II λ1640 emission, suggesting a tentative correlation between (O/H) and $EW_0$ (He II).

### 3.3. No indications of AGN activity

Finally, we investigate the presence and possible contribution of AGN activity since these sources are extremely bright in the UV, presenting apparent magnitudes similar to QSOs ($M_{\text{UV}} \simeq -23.3$ to $\simeq -24.6$, see Table 1). In addition to the already mentioned P-Cygni wind lines, the GTC spectra of these sources also reveal photospheric absorption lines (see Appendix A), several of them detected with high significance. Examples of these are O IV λ1343 or S V λ1501. The identification of these inherently faint lines provides clear evidence that the UV luminosity is predominantly governed by starlight (e.g., González Delgado et al. 1998; Shapley et al. 2003) as they originate from the photospheres of hot and massive stars. It is worth noting that even a minor contribution from an AGN to the UV continuum, which is featureless in these spectral regions, would result in the attenuation

of these lines to a level imperceptible given the SNR of our spectra. As an additional test, we also look at the spectral profile of the Balmer lines in the EMIR rest-optical spectra. All sources show narrow profiles in Balmer lines ($H\beta$ or $H\alpha$) with intrinsic line widths between $FWHM \lesssim 400 \text{ km s}^{-1}$ (unresolved) and $\simeq 470 \text{ km s}^{-1}$. This contrasts with the much broader line profiles of the non-resonant He II λ1640 which show $FWHM \simeq 1150-2900 \text{ km s}^{-1}$ (Table 2), clearly indicating stellar origin (WR and/or VMS). We thus conclude that the AGN contribution to the UV luminosity and strength of He II in these sources is likely to be residual or null.

## 4. Empirical analysis on the presence of VMS

### 4.1. Evidence of VMS in UV-bright galaxies

Several empirical arguments suggest a significant contribution of VMS in our sources, or at least in some of them, that we now describe. The most important one refers to the very high $EW_0$ (He II) observed in our sources and the comparison with normal typical star-forming galaxies where WR stars are expected. Indeed, WR stars are formed in basically all conditions and, thus, are expected to be present in all star-forming galaxies (of course, knowing that their contribution to the He II line will be dependent on star-formation histories, age, metallicity, and other factors). In fact, broad He II emission is recurrently observed in the rest-frame UV spectra of normal star-forming galaxies of enough S/N (e.g., Shapley et al. 2003; Noll et al. 2004; Cabanac et al. 2008; Dessauges-Zavadsky et al. 2010; Jones et al. 2012; Marques-Chaves et al. 2020a), but its strength is significantly weaker than in our sources.

To investigate the contribution of VMS and/or WR stars in our UV-bright galaxies, we first create a composite spectrum to show the resulting average spectral shape of He II in our sources and compare it with the average rest-frame UV spectrum of normal LBGs from Shapley et al. (2003) where WR stars are expected. First, the spectra of the 13 sources were de-redshifted using the systemic redshifts shown in Table 1 and were resampled onto a common wavelength grid using linear interpolation. Next, we normalized the spectra using several spectral regions that are relatively free of emission or absorption features (i.e., excluding regions with ISM absorption, wind features, or nebular emission). Finally, we stacked all spectra by averaging the flux in each spectral bin.

Figure 3 shows the resulting stacked GTC spectrum of UV-bright galaxies and the comparison with that of typical ($L_{\text{UV}}^*$) star-forming galaxies (Shapley et al. 2003). As seen from Fig. 3, the composite of UV-bright galaxies shows much stronger He II emission than normal LBGs, by a factor of $\ge 2$. Using the same integration windows as used in Sect. 3.1 to derive the strength of He II λ1640 in our individual spectra, we measure $EW_0$ (He II) = $3.19 \pm 0.07$ Å for the stacked spectrum of UV-bright galaxies and $EW_0$ (He II) = $1.37$ Å for the Shapley et al. (2003) composite. Since WR stars produce the bulk of the He II emission observed in the Shapley et al. (2003) composite spectrum, as quantitatively demonstrated by Brinchmann et al. (2008) and Eldridge & Stanway (2012), we suggest that the excess of He II observed in UV-bright galaxies is due to the contribution of VMS. Furthermore, the C IV P-Cygni appears much stronger in the composite of UV-bright galaxies than in normal LBGs. As shown later in Sect. 5.3, VMS contribute significantly to the strength of C IV in the integrated spectrum.

The differences in the strengths of the He II and C IV features seen in the composites of UV-bright and normal galaxies may
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Comparison between our composite spectrum of UV-bright galaxies (black) and that of typical ($L_{UV}^{*}$) Lyman break galaxies at $z \sim 3$ from Shapley et al. (2003) (red). Measurements of the EW$_{0}$ (He II) for both composites are shown in the caption. The region near the He II λ1640 emission profile is normalized with a different normalization factor to match the level of the continuum. The measured EW$_{0}$ (He II) for the composite of UV bright galaxies in this work is $3.19 \pm 0.07$ Å while for the Shapley composite is $1.37$ Å.</figcaption></figure>

also arise from differences in metallicities between these two galaxy populations, recognizing the metallicity dependence on the strength of the He II and C IV wind lines. While precise measurements of the metallicity of the Shapley et al. (2003) composite are not available, Steidel et al. (2014) and Sanders et al. (2021) have obtained and studied the rest-frame optical spectra of similar galaxies (LBGs at $z \sim 2-3$). These works derive $12+\log(O/H)$ between 8.2 and 8.6 which is broadly consistent with that obtained for most of our sources (Fig. 2). So, overall we expect the O/H abundance of the Shapley et al. (2003) composite and our UV bright galaxies to be broadly similar, suggesting that the differences between the He II and C IV profiles in our sources and normal LBGs are mainly due to the presence/absence of VMS.

In addition, a significant fraction of our sources show other spectral features that are present in the spectra of individual VMS (e.g., R136-a1, -a2, -a3, or R146 see: Crowther et al. 2016; Brands et al. 2022; Martins & Palacios 2022). In particular, broad emission in N IV λ1486 and a significant P-Cygni line profile in N IV λ1719 are clearly detected in J0006+2452, J0036+2725, J0110-0501, J0146-0220, J1157+0113, J1220+0842, and partially/barely detected in J0115+1837, J0121+0025, and J1335+4330 (see Appendix A). The detection of these profiles is also reflected in the composite spectrum of our UV-bright galaxies shown in Fig. 3. These profiles are also seen in the integrated spectra of clusters or compact star-forming regions where VMS are suspected, including NGC 3125-A1 analyzed in Wofford et al. (2014, 2023), or J1129+2034, J1200+1343, and J1215+2038 analyzed by Senchyna et al. (2017, 2021). Furthermore, these features are also predicted by the VMS models of Martins & Palacios (2022). However, we note that some WR-dominated clusters also show these spectral features, and thus they are not restricted to VMS as pointed out by Martins et al. (2023). As such, the N IV λ1486 and N IV λ1719 features could be also expected in the composite spectra of $L_{UV}^{*}$ LBGs, but these are faint or not detected as shown in Fig. 3 possibly due to a dilution effect from the contribution of the more numerous OB-type stars.

Finally, the last argument is related to the star-formation histories. The He II line is boosted in bursty/instantaneous star-

formation histories when compared to continuous star-formation histories, regardless of whether He II originates from VMS or WR stars (see Sect. 5). The rest-frame UV spectra of our sources show signs of smooth/continuous star-formation histories, as they show strong P-Cygni line profiles in both N V λ1240 and Si IV λ11393, 1402 which are basically impossible to fit simultaneously with models of instantaneous bursts. In addition, SED analysis of the multiwavelength photometry and nebular emission performed by Marques-Chaves et al. (2020b, 2021, 2022) for a few of the sources studied here also suggest continuous star-formation histories, yet with very short ages ($\sim 10$ Myr).

## 4.2. Spectral comparison of UV-bright sources with other VMS-dominated clusters

Here we present a comparative analysis of the rest-frame UV spectra of some of our UV-bright sources with those of known VMS-dominated clusters. Figure 4 shows the spectra of J0006+2452 ($z = 2.377$, $M_{UV} = -24.17$, top), J0110-0501 ($z = 2.368$, $M_{UV} = -23.33$, middle), and J1220+0842 ($z = 2.469$, $M_{UV} = -24.40$, bottom). Overlaid are the spectra of three young star-clusters where the presence of VMS is confirmed or suspected. Specifically, we include the spectrum of the R136 cluster within 30 Dor/LMC (Crowther et al. 2016, blue), SB 179 (Senchyna et al. 2017, green), and the highly magnified Sunburst cluster at $z \approx 2.371$ (Meštrić et al. 2023, red). This comparative analysis serves primarily for illustrative purposes, aimed at highlighting the resemblances between some of our UV-bright galaxies and established examples of VMS-dominated clusters.

Overall, the spectra of our sources resemble those of VMS-dominated clusters. In particular, the spectral features typically associated with WN and WNh stars (including VMS), such as the intense and broad He II λ1640 emission and the N IV profiles at λ1486 and λ1719 (highlighted in gray regions), are present in both samples presenting also similarities in their spectral shapes. A closer examination also reveals differences. For example, the N IV λ1486 appears stronger in J0006+2452 and J0110-0501 than in R136 and SB 179, although effects of different resolution and SNR between spectra prevent a fair comparison. The blueshifted absorption component of the N IV λ1719
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Comparison of the normalized GTC rest-UV spectrum of some of our sources with star clusters containing VMS (R136 in blue: Crowther et al. 2016; SB 179 in green: Senchyna et al. 2017; and the Sunburst cluster in red: Meštrić et al. 2023). Gray-shaded regions mark the location of the most important spectral features associated with VMS: broad emission in N IV λ1486 and He II λ1640, and a P-Cygni line profile in N IV λ1719. Some physical properties of these systems, including M<sub>UV</sub>, stellar mass, and metallicity, are also shown (see text for details).</figcaption></figure>

P-Cygni line profile is also more pronounced in J0006+2452 and J0110–0501, showing substantially larger terminal velocities than what is seen in the spectra of R136 or SB 179. On the other hand, the He II profiles are almost identical between our sources and the comparison sources. Of course, variations in the spectral shapes and intensities of He II or N IV profiles or others are expected and can be attributed to factors such as age, metallicity, IMF, or other intrinsic properties of massive and very massive stars (e.g., mass-loss, terminal wind velocities, rotation, etc.). As already highlighted in Fig. 1, the spectra of our 13 sources exhibit diverse He II profiles, ranging from nearly Gaussian/symmetric profiles (e.g., J0006+2452 or J0115+1837, although with different line widths) to asymmetric profiles (J1220+0842 or J1335+4330), or more complex ones with multiple emission/absorption peaks (e.g., J0031+3545 or J0121+0025). Similarly, the spectra of nearby star-forming regions analyzed by Senchyna et al. (2021) or Martins et al. (2023) also reveal differences in the shape in the He line (and

N IV as well) from source to source (VMS- or WR-dominated sources).

### 4.3. Empirical classification of VMS- or WR-dominated sources

Martins et al. (2023) provide observational criteria to distinguish sources dominated by VMS or WR stars. In their scheme, VMS-dominated sources show intense and broad He II λ1640 emission, with EW<sub>0</sub> ≥ 3.0 Å at least. This limit is motivated by the empirical predictions of Schaerer & Vacca (1998) of the intensity of the He II λ1640 line from various types of WR stars in young stellar populations, which reach a maximum of EW<sub>0</sub> ≈ 3.0 Å for an instantaneous burst at Z/Z<sub>☉</sub> ≈ 0.5. Furthermore, the spectra of VMS-dominated sources also present broad optical He II λ4686, but on the other hand, the emission in N III+C III 4640 Å–4650 Å and C IV 5801 Å–5812 Å, should be weak or absent (Martins et al. 2023). Given that our rest-optical
spectra are too shallow to detect any of these optical lines or
they are not available at all, we restrict our classification to the
strength of the 1640 Å line.

Following the observed strengths in the He II λ1640 line in
our sources (see Table 2), 8 of them show EW<sub>0</sub> (He II) > 3.0 Å
and thus are likely to be dominated by VMS. An addi-
tional source, J1220+0842, represents a special case. Its
EW<sub>0</sub> (He II) = 2.52 ± 0.33 Å falls below our selection thresh-
old. However, the source shows strong and broad N IV λ1486
and a P-Cygni line in N IV λ1719, both of which are also VMS
signatures. This galaxy also has the lowest metallicity in our
sample (12 + log(O/H) = 8.13 ± 0.19, Fig. 2), which might be
expected to weaken the He II line (due to weaker stellar winds).
Given its spectral similarity to the Sunburst cluster and others
(Meštrić et al. 2023, see Fig. 4), we add this source to our sam-
ple of candidate VMS hosts.

In summary, we conclude that among the 13 UV-bright
galaxies analyzed in this work, nine of them are likely
VMS-dominated sources. We conservatively consider them as
candidates for the reasons described above, and in particu-
lar due to the lack of deep observations in the rest-optical
blue and red bumps, which according to Martins et al. (2023)
are necessary to unambiguously distinguish VMS- and WR-
dominated sources. The three strongest He II emitters in our
sample, J0006+2452 (EW<sub>0</sub> = 4.54 Å), J0110–0501 (EW<sub>0</sub> =
4.57 Å), and J1316+2614 (EW<sub>0</sub> = 4.48 Å) stand out to be
the best VMS candidates in our sample. The observed EW<sub>0</sub>
(He II) in these sources are similar to those in the spectra of
local star-clusters or star-forming regions dominated by VMS
(e.g., R136, NGC 3125-A1, J1129+2034/SB179, Crowther et al.
2016; Senchyna et al. 2021; Wofford et al. 2023), even noting
that these local clusters may have their He II emission enhanced
by bursty/single-age star-formation histories.

**5. Population synthesis models**

Population synthesis models can be used to infer stellar pop-
ulations from the comparison of synthetic and observed UV
spectra of star-forming galaxies. This approach was applied
in a few local star clusters like NGC3125-A1, NGC5253-
5, or II Zw 40, where intense He II λ1640 emission is detected
(EW₀ ≥ 4.0 Å; Wofford et al. 2023, 2014; Leitherer et al. 2018;
Smith et al. 2016). The strength of He II λ1640 observed in these
clusters is well above that any population synthesis models
without VMS can predict, leading these authors to suggest the
presence of VMS. For these reasons, it is important to discuss
first the validity and limitations of different synthesis models
with/without VMS in predicting the He II λ1640 emission. In this
Section, we investigate the ability of various population synthe-
sis models to account for the strong UV lines of our sources. We
consider models available in the literature and apply new models
of VMS from Martins & Palacios (2022).

5.1. BPASS models

We first retrieve the spectral energy distributions (SEDs)
produced by BPASS v2.2.1 models (Eldridge et al. 2017;
Stanway & Eldridge 2018; Byrne et al. 2022). We use BPASS
models with binary stellar populations and upper mass limits
of 100 M⊙ and 300 M⊙ with a Salpeter IMF slope of -2.35.
Furthermore, we select a metallicity of Z = 0.006 since
it is the closest to the average metallicity of these sources
(12 + log(O/H) = 8.36 ± 0.11, Sect. 3) and also the closest

to the metallicity of the Large Magellanic Cloud (LMC) for
which the new VMS models of Martins & Palacios (2022) were
developed. This is also consistent with the previous works of
Marques-Chaves et al. (2020b, 2021, 2022) on a subset of our
sources.

BPASS models have instantaneous star formation histories (ISFH) normalized to a burst mass of $10^6 M_\odot$. The detection of strong P-Cygni line profiles in both N V λ1240 and Si IV λ1400 in the spectra of these sources, which originate respectively from O-type main-sequence and supergiant phases, is difficult to be explained by single ISFH model (e.g., see Chisholm et al. 2019) and suggests a more smooth/continuous star-formation history. Hence we convert the SEDs of ISFH models to constant star formation history (CSFH) models, which we expect in our sources. We discuss how the models perform at reproducing the strength of the main UV lines in Sects. 5.3–5.5.

5.2. Models including VMS self-consistently

5.2.1. Stellar models

Spectra of four single-star VMS are presented in Martins & Palacios (2022) at masses of 150, 200, 300, and 400 $M_\odot$, and were computed using nonrotating evolutionary models using the code STAREVOL (Siess et al. 2000; Lagarde et al. 2012; Amard et al. 2019). VMS models use the mass loss prescription of Gräfener (2021) which are developed for a metallicity of 0.4 $Z_\odot$, considering $Z_\odot = 0.0134$. The model traces the evolution of each VMS from the zero-age main sequence to the end of the H-burning phase which lasts somewhere from 2 to 2.5 Myr depending on the birth mass of the VMS.

In these VMS models, the proximity of the star to the
Eddington limit makes the scaling of mass loss rates to the
Eddington factor steeper, which differs from the prescription
used in Vink et al. (2001). This leads to a higher mass loss rate
with optically thicker winds. It results in higher helium surface
abundance which, combined with higher wind density, produces
stronger He II λ1640 emission. The He II appears as a P-Cygni
line profile that strengthens with age (see, Martins & Palacios
2022). Apart from He II λ1640, the N V λ1240 and C IV λ1550
profiles also appear as P-Cygni line profile.

Other important features that show up in the theoretical
UV spectra of VMS are the N IV λ1719 and N IV λ1486 pro-
files. The N IV λ1719 appears weak at the ZAMS but evolves
into a stronger P-Cygni line profile with age. The N IV λ1486
profile appears as an emission only after 1.5 million years of
evolution of the VMS and becomes stronger with age. The
N IV λ1486 profile is not seen in the less massive O-type super-
giants even at solar metallicity (Walborn et al. 1985; Bouret et al.
2012) but it appears on the UV spectra of WN and WNH stars
(Hamann et al. 2006; Hainich et al. 2014). The stellar emission
in N IV λ1486 appears broad (FWHM > 1000 km s⁻¹) in VMS
models (Martins & Palacios 2022) and in VMS-dominated clus-
ters (e.g., R136; Crowther et al. 2016) and can be distinguish-
able from the nebular emission observed in some narrow-lined
UV-selected AGNs (e.g., Hainline et al. 2011) and star-forming
galaxies. This emission feature only rarely occurs in broad-lined
quasars in the SDSS (see, Bentz et al. 2004; Jiang et al. 2008),
and is also detected in a few radio galaxies (see, Vernet et al.
2001; Humphrey et al. 2008). The N IV λ1719 P-Cygni line pro-
file has been observed in some hot Of-type stars (e.g., Conti et al.
1996).

BPASS v2.2.1 models extending to 300 $M_\odot$ do incorporate
evolutionary tracks for both single and binary very massive
stars. However, these adopt the same wind mass loss rates in
the VMS regime as for massive O-Stars (Vink et al. 2001). The BPASS spectral synthesis also assigns these stars template spectra derived from the same stellar atmosphere template grid as other hot, massive stars (derived from the WMBasic atmosphere code), rather than producing custom spectra in this regime to reflect the different wind regime. Recent works on VMS have suggested that both the mass loss prescription and the resulting atmospheric radiative transfer of VMS might differ substantially from hot stars below 100 $M_\odot$ (see e.g., Gräfener 2021; Björklund et al. 2023).

### 5.2.2. Producing synthetic population spectra with VMS

To describe the integrated UV spectra of star-forming populations we combine the SEDs of VMS models from Martins & Palacios (2022) that are available for VMS of 150, 200, 300, and 400 $M_\odot$ to normal stellar populations from the BPASS v2.2.1 models that include stars up to 100 $M_\odot$; the method is similar to that adopted in Martins & Palacios (2022). We have defined mass bins in the range [100, 175], [175, 225], [225, 275], [275, 325], and [325, 475] and calculated the number of VMS required in these mass bins by extrapolating the IMFs to the upper bin mass. In this approach, we have used a corrected equation for a continuous IMF, which produces significantly different numbers than that has been used in Martins & Palacios (2022). The details of this process is described in Appendix B. In short, we take the BPASS models up to 100 $M_\odot$ and assume the same IMF (Salpeter) for stars with higher masses. Then, the VMS spectra are resampled to a common wavelength grid using the SPECTRES<fnref n="5" /> (Carnall 2017) python library. We consider five different maximum masses of VMS, between 175 and 475 $M_\odot$. This is done for all burst models with ages up to 2.5 Myr, after which no VMS is left in the population. The ISFH models including VMS are then also used to compute models for constant SFR. We note that given the assumptions of Martins & Palacios (2022), the post-main sequence evolution of VMS is neglected here, and we also neglect VMS in binary systems.

### 5.2.3. Validation of synthesis models with/without VMS

The VMS models used in this work were introduced by Martins & Palacios (2022). When coupled with synthetic population spectra from BPASS, Martins & Palacios (2022) demonstrated that their models successfully explain the strength of He II $\lambda$1640 observed in the R136 cluster ($EW_0 \approx 4.5$ Å), which is produced primarily from the few VMS located in the core of R136 (Crowther et al. 2016). Other synthetic population spectra including VMS like the new Bruzual & Charlot models (Plat et al. 2019), which incorporate VMS to 300 $M_\odot$, have been also used in studies to investigate the presence of VMS in local star-forming regions (e.g., Senchyna et al. 2021; Wofford et al. 2023; Smith et al. 2023). Even noting that Bruzual & Charlot and Martins & Palacios (2022) models have different VMS prescriptions (see references therein for details), they successfully reproduce several UV spectral features observed in VMS, and in particular, by boosting the He II $\lambda$1640 strength in the integrated synthetic population spectra. We thus conclude that the VMS models used in this work are valid to search for VMS in our sources.

A different situation may apply to standard synthetic spectra without VMS, where He II is expected to be produced primarily by classical WR stars. Figure 5 shows that the BPASS models with the upper mass cutoff of 100 $M_\odot$ (i.e., without VMS) struggle to exceed 1 Å and 2 Å in He II $\lambda$1640 for continuous and instantaneous star formation histories, respectively, at the adopted metallicity $Z = 0.006$. These are slightly lower than those found in Brinchmann et al. (2008) using Starburst99 continuous star formation models at similar metallicity, but including also the contribution of stars (reaching a maximum of $EW_0 \approx 1.4$ Å), or the burst models of Schäerer & Vacca (1998) that rely on empirical line luminosities of different types of WR stars and reach a maximum of $EW_0 \approx 3.0$ Å (at $Z/Z_\odot \approx 0.5$). Furthermore, the UV spectroscopic analysis by Martins et al. (2023) of seven WR-dominated local star-forming regions at $Z \approx 0.5 Z_\odot$ show somehow stronger He II $\lambda$1640 emission than the ones measuring here using BPASS models. They measure strengths in the 1640 Å line of $\approx 1.64-2.35$ Å for six of them, and an extreme value of $EW_0$ (He II) $\approx 3.89$ for Tol 89-A, a known local WR-dominated cluster (Sidoli et al. 2006). Adopting a metallicity of 0.008 instead of 0.006 for the $M < 100 M_\odot$ BPASS population only marginally increases the He II emission. From these results and the analyses of the optical blue and red bumps, Martins et al. (2023) suggested that the treatment of WR stars should be improved in current population synthesis models (see further details in Martins 2023).

We now examine the predicted UV spectra from synthesis models with and without VMS and compare them to our observations. To do so we use the models for constant star-formation, except if otherwise stated.

### 5.3. Strength of $N\nu\lambda$1240, and $C\nu\lambda$1550 P-Cygni line profiles

The rest-frame UV spectra of our sources (Appendix A) show indications of the presence of young and massive stars. In particular, strong P-Cygni line profiles are observed in $N\nu\lambda$1240 and $C\nu\lambda$1550 which are sensitive to age and metallicity in the UV (see e.g., Chisholm et al. 2019). Massive stars and VMS are very hot stars with surface temperatures ranging from 40 000 K to 60 000 K. As such, the dominant ionization state of stellar wind changes from $N^{3+}$ to $N^{4+}$, at a given metallicity (see, Kudritzki et al. 1998; Lamers & Cassinelli 1999; Leitherer et al. 2010). This results in a strong P-Cygni line profile with a blueshifted absorption and a redshifted emission. As stars evolve to relatively older ages the dominant ionization state changes to $N^{3+}$ resulting in a decrease in the strength of $N\nu\lambda$1240 P-Cygni line profile, that is, decreasing the strength of the corresponding absorption and emission. This makes the strength of the $N\nu\lambda$1240 P-Cygni line profile a potential tracer of the age of the stellar population in young star-forming galaxies.

We perform a quantitative study on the strength of $N\nu\lambda$1240 P-Cygni line profile by measuring the strength of the absorption component of the $N\nu\lambda$1240 P-Cygni profile, which is the most affected by the change in age both in our sources and models (Chisholm et al. 2019). We measure the $EW_0$ of the $N\nu$ absorption in the spectral window from 1230 Å to 1240 Å, avoiding as much as possible the contamination from the $Ly\alpha$ emission, the left-most part of this integration is chosen to be around ~3500 km s$^{-1}$ away from the $Ly\alpha$ rest-wavelength. However, due to nebular $Ly\alpha$ emission from the galaxy, which is not accounted for by our models and which can extend redward by several hundred km s$^{-1}$, we refrain from using the $N\nu\lambda$1240 absorption as

<footnote n="5">⁵ https://spectres.readthedocs.io/en/latest/</footnote>
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Variation of the equivalent width (EW<sub>0</sub>) of the absorption component of the N v λ1240 (middle), C IV λ1550 (bottom) P-Cygni line profiles and the stellar He II 1640 Å emission (top) as a function of age assuming a continuous star formation history (left) and instantaneous burst models (right). Different tracks represent models with different IMF upper mass limits (M<sub>up</sub>). Original BPASS models with M<sub>up</sub> = 100 M<sub>☉</sub> and 300 M<sub>☉</sub> (using the wind prescription from Vink et al. 2001) are shown in red and orange, respectively, while BPASS models with the new VMS template spectra of Martins & Palacios (2022) are shown in yellow, green, blue, pink, and brown for M<sub>up</sub> = 175 M<sub>☉</sub>, 225 M<sub>☉</sub>, 275 M<sub>☉</sub>, 325 M<sub>☉</sub>, and 475 M<sub>☉</sub> respectively.</figcaption></figure>

an absolute age indicator. For the estimation of the continuum level, we select the spectral region from 1269 Å to 1272.5 Å to avoid the contribution of the low-ionization Si II λ1260 absorption line.

The middle panels of Fig. 5 show how the N v λ1240 absorption strength varies as a function of age for models assuming a CSFH (left) and ISFH (right) with different upper mass limits.

At a very young ages (~1 Myr), N v λ1240 shows very strong absorption, between EW<sub>0</sub> ≈ −4.4 Å for a non-VMS IMF (i.e., M<sub>up</sub> = 100 M<sub>☉</sub>, red curve in Fig. 5) and EW<sub>0</sub> ≈ −5.6 Å when VMS are included. At later ages (≥2 Myr), the strength of the absorption component of N v decreases, reaching a plateau of EW<sub>0</sub> ≈ −1.5 Å at ≳20 Myr for CSFH models. This figure also shows that the
inclusion of VMS provides only a marginal contribution to the N v λ1240 absorption strength at any age, although this effect is more significant at ≤2 Myr. It is also worth noting that the original BPASS models with $M_{up} = 300 M_{\odot}$ (orange curve in Fig. 5) and the ones using VMS from Martins & Palacios (2022) predict similar strengths of N v, even noting that they assume different wind prescriptions: Vink et al. (2001) for the former, and Gräfener (2021) for the latter. Applying the same methodology to our sources, we measure $EW_0$ of the N v absorption between ≈−5.5 Å and ≈−1.1 Å (Table 2), indicating different ages according to Fig. 5. In the absence of any contamination by the Lyα emission, the sources J0006+2452, J0036+2725, and J0110−0501 have the strongest absorption in N v with $EW_0$ ≲ −3.5 Å and, therefore, stand out to be the youngest systems in our sample. The rest of the sources appear slightly older, in particular, those with $EW_0$ ≳ −3.0 Å.

We perform a parallel quantitative investigation on the C IV λ1550 P-Cygni profile, measuring its absorption strength. Given the presence of the interstellar medium (ISM) component in C IV, we measure the $EW_0$ of the stellar C IV λ1550 absorption within the spectral range spanning from 1530 Å to 1543 Å, that is, at least ~1000 km s⁻¹ away from the rest-velocity of C IV λ1548. Even with the selection of these spectral windows, we cannot rule out the ISM contamination in the stellar C IV λ1550 absorption, in particular in the presence of strong outflows. To establish the continuum level, we designate the spectral region from 1509 Å to 1519 Å on the left and a region from 1570 Å to 1580 Å on the right side of the C IV λ1550 profile. For J1316+2614, the continuum window on the right is exclusively selected from 1578 Å to 1580 Å due to atmospheric line contamination in the optical GTC spectra.

The bottom panels of Fig. 5 show how the C IV λ1550 absorption strength varies as a function of age for models assuming a CSFH (left) and ISFH (right) with different upper mass limits of the IMF ($M_{up} = 100–475 M_{\odot}$). Unlike N v λ1240 absorption profile, the C IV λ1550 absorption is susceptible to different mass loss prescriptions. The Gräfener (2021) mass loss prescriptions used in the Martins & Palacios (2022) models significantly boost C IV λ1550 absorption; for the CSFH model with IMF $M_{up} = 300 M_{\odot}$ at ~1 Myr, its strength is boosted by 2.4 Å in comparison to the BPASS model at the same age and IMF. At this age, the C IV shows very strong absorption between $EW_0$ ≈ −4.7 Å to −6.2 Å for models that include VMS with Gräfener (2021) mass loss prescription. The strength only decreases slightly till 2.5 Myr. While the BPASS model with IMF $M_{up} = 300 M_{\odot}$ that includes Vink et al. (2001) mass loss prescription, the strength starts with $EW_0$ ≈ −2.8 Å and reaches its strongest at $EW_0$ ≈ −4 Å at 2.5 Myr before reducing again.

We measure the C IV absorption strength for sources between $EW_0$ ≈ −6.7 Å to −2.0 Å (Table 2). We can not rule out the contamination from the ISM in the C IV for our sources, and the contamination results in the C IV absorption strength being overestimated in our sources. But, the degeneracy related to metallicity has also been affecting the strength of C IV absorption in our sources, as mentioned earlier, it is a metallicity-sensitive feature in the UV spectra. The two sources J1220+0842 and J0850+1549 on the right-most end of Fig. 6 with lowest C IV absorption strength are the lowest metallicity sources in our sample with $12 + \log(O/H) \approx 8.13$ and 8.27 respectively. While the sources with stronger C IV absorption strength have relatively higher metallicities.

### 5.4. Strength of He II λ1640 emission: VMS indicator

VMS have been already identified and characterized in detail throughout UV and optical spectroscopy, either individually (e.g., Massey & Hunter 1998; Bestenlehner et al. 2014; Crowther et al. 2016) or in integrated spectra of unresolved star-forming regions (e.g., Wofford et al. 2014, 2023; Smith et al. 2016, 2023; Meštrić et al. 2023). From these empirical results and models, broad and intense He II λ1640 emission appears to be ubiquitous in VMS, making it the best indicator of the presence of VMS in the rest-frame UV (see Martins & Palacios 2022). Since classical WR stars can also produce broad He II emission, the strength of He II depends thus on the relative contribution of VMS and WR stars in integrated light spectra.

The He II λ1640 emission equivalent width has been measured by using both Gaussian fitting and flux integration methods, as shown in Fig. 1. However, the He II profiles exhibit a non-Gaussian nature in the majority of our sources, while manifesting asymmetrical P-Cygni features in the models. To compute the He II λ1640 emission equivalent width for both the sources and the models, we have employed the flux integration method, as described in Sect. 3.1. The top panels of Fig. 5 show the variation of the equivalent width of the stellar He II emission as a function of age for CSFH (left) and ISFH (right) models assuming different IMF upper mass limits. As shown in this figure, BPASS models with $M_{up} = 100 M_{\odot}$ predict the weakest He II intensity within our different models, with a maximum of $EW_0$ (He II) ≈ 0.6 Å at ≥6–8 Myr for CSFH, that is, when the contribution of WR stars reaches its maximum. The original BPASS model with an upper mass limit of 300 $M_{\odot}$ (in orange), which has the wind prescription from Vink et al. (2001), can only predict a maximum of $EW_0$ (He II) ≈ 0.8 Å. For comparison, earlier models of Brinchmann et al. (2008) predict a maximum $EW_0$ (He II) ≈ 1.4 Å for CSFH at similar metallicity. On the other hand, BPASS coupled with VMS models of Martins & Palacios (2022) show much stronger He II emission, with $EW_0$ (He II) ranging from ≈1.8–7.6 Å, depending on the age and $M_{up}$. For these models, the He II strength peaks near 2.5 Myr, approximately the lifetime of the VMS. From Fig. 5, it is also evident that the impact of VMS in integrated spectra is much stronger in the He II emission (boosting it by a maximum factor of ≈7.5) than in the strength of N v P-Cygni line profile (left panel of Fig. 5).

### 5.5. He II λ1640 against N v λ1240 and C IV λ1550: Evidence of VMS

Putting together the results obtained in Sects. 5.3 and 5.4, we show in Fig. 6 the relationship between the strength of the absorption component of the N v and C IV P-Cygni line profiles and the stellar He II emission. Figure 6 also shows the measurements of $EW_0$ (N v) and $EW_0$ (He II) of our sources (listed in Table 2). Figure 6 provides strong evidence for the presence of VMS in most of our sources. The sources that lie on the left-hand side of this figure, J0006+2452, J0036+2725, and J0110−0501, show $EW_0$ (N v) ≲ −4.5 Å and are the youngest systems in our sample. These sources show $EW_0$ (He II) ≈ 2.7–4.6 Å, requiring a significant number and contribution of VMS, or IMF upper mass limits between $M_{up} = 225 M_{\odot}$ and $M_{up} = 275 M_{\odot}$. For the remaining sources showing $EW_0$ (N v) ≈ −4.5 Å to ≈−1.0 Å and $EW_0$ (He II) ≈ 2.0–4.5 Å, stellar population with VMS are also clearly preferred and possibly required, although the IMF upper mass limit is not well constrained given the large uncertainties in our measurements. An exception may be J1316+2614, showing
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Variation of the EW$_0$ of He II 1640 Å emission (y-axis) with the absorption component of the N v λ1240 P-Cygni line profile (x-axis) on the left and the same with the absorption component of C IV λ1550 P-Cygni line profile (x-axis) on the right as a function of age assuming a continuous star formation history. Measurements of our UV-bright galaxies are marked with solid circles. Different tracks represent models with different IMF upper mass limits ($M_{up}$). Original BPASS models with $M_{up} = 100 M_⊙$ and 300 $M_⊙$ (using the wind prescription from Vink et al. 2001) are shown in red and orange, respectively, while BPASS models with the new VMS template spectra of Martins & Palacios (2022) are shown in yellow, green, blue, pink, and brown for $M_{up} = 175 M_⊙, 225 M_⊙, 275 M_⊙, 325 M_⊙$, and $475 M_⊙$ respectively.</figcaption></figure>

EW$_0$ (N v) ≈ −1.9 Å and EW$_0$ (He II) ≈ 4.5 Å, which requires an IMF upper mass limit of $M_{up} \ge 475 M_⊙$.

Overall Fig. 6 shows that different BPASS+VMS models of Martins & Palacios (2022) can predict relatively well both the observed strength of He II λ1640 emission and N v λ1240 absorption of these sources. On the other hand, models without VMS are not able to reproduce the strength of He II λ1640 emission and C IV λ1550 absorption. In line with the models shown in Fig. 6, most of the sources in our sample require a significant contribution of VMS in integrated stellar populations to explain the observed strengths of wind lines, in particular He II.

## 5.6. Spectral comparison

In this section, we have compared the strengths of the observed He II stellar emission in our sources to those from different synthetic models. We now perform a qualitative comparison of the spectral profiles of the He II line between our sources and models.

The top panel of Fig. 7 shows the resulting stack spectrum of our 13 UV-bright galaxies. The most prominent VMS features are highlighted in red (He II λ1640 and the N IV profiles at 1486 Å and 1719 Å), as well as other stellar wind features produced by normal massive stars and VMS (N v λ1240, Si IV λ11393, 1402, and C IV λ1550). Using the same methodology described in Sects. 5.3 and 5.4, we measure the EW$_0$ of the absorption component of the N v P-Cygni line profile and the stellar emission of He II of the stacked spectrum, finding EW$_0$ (N v) = −3.0 ± 0.07 Å and EW$_0$ (He II) = 3.19 ± 0.07 Å, respectively. We also show in Fig. 7 the normalized spectrum of

a BPASS+VMS model (Martins & Palacios 2022) assuming an IMF with an upper mass cutoff of $M_{up} = 225 M_⊙$ with a continuous star-formation over 10 Myr (blue). The model spectrum was also convolved to match the spectral resolution of the GTC spectra (R ~ 700). Overall, the synthetic spectrum can reproduce several important features associated with massive stars (e.g., N v, Si IV, or C IV). The spectral features related to VMS are also relatively well-reproduced by the model spectra, but differences between their spectral shapes are also evident. Given the fact that the spectrum of the composite of UV-bright galaxies is well reproduced by a constant star formation history model with ≈10 Myr, it is more representative of a VMS+WR spectrum, where the contribution from the VMS in the rest-frame UV spectrum is significant.

The bottom right panel of Fig. 7 shows a close look at the He II line. As seen in the figure, BPASS models including VMS from Martins & Palacios (2022) do a better job in reproducing the strength of He II in our sources than models without VMS. However, the He II emission appears much narrower in the models than in our sources. The right panel of Fig. 7 shows the He II line widths (FWHM) obtained for our sources which range from ≈1000–3000 km s⁻¹ with a mean value of FWHM ≈ 2300 km s⁻¹. These are much larger on average than the line widths obtained from the models which range from ≈1000–1500 km s⁻¹. The shape of He II λ1640 profile of the sources is not fully reproduced by theoretical VMS models from Martins & Palacios (2022); while the shape in C IV λ1550 profile matches well. This may have to do with details of the wind properties such as clumping or the exact shape of the velocity structure. This issue should be further addressed in future studies of VMS and related populations.
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** A spectral and quantitative comparison of FWHM between the model and sources. Top panel shows the normalized composite spectrum using the 13 UV-bright sources used in this work (black). The synthetic spectrum of a BPASS+VMS model is also shown in blue and assumes an upper mass cutoff of the IMF of $M_{up} = 225 M_{\odot}$ with a continuous star-formation over 10 Myr (using the VMS models of Martins & Palacios 2022). Spectral features associated with only VMS are marked in red while the spectral features produced by VMS and normal massive stars are marked in green. The bottom left panel shows a zoom-in at the He II line. BPASS models with and without VMS are also shown (blue and red, respectively). The bottom right panel illustrates the differences between the He II line widths (FWHM) measured in our sources (black) and model spectra (blue region).</figcaption></figure>

# 6. Discussion

## 6.1. Incidence of VMS in UV-bright galaxies and in other galaxy populations

We now investigate the incidence of VMS in different galaxy populations. The left panel of Fig. 8 shows the intensity of the He II λ1640 line for different types of sources. These include young star clusters and H II regions (green) from the compilation by Martins et al. (2023), for which we only show VMS- and WR-dominated sources from their classification<fnref n="6" />, normal LBGs (red), and the UV-bright galaxies studied in this work (blue), as indicated by their respective UV absolute magnitudes.

As illustrated in the left panel of Fig. 8, a significant fraction of UV-bright galaxies exhibit prominent He II λ1640 line with $EW_0 \approx 3.0–5.0$ Å. Similar strengths are also observed in the spectra of local star clusters where VMS are suspected but are significantly higher than those found in WR-

dominated clusters. Furthermore, additional spectral features characteristic of VMS, such as the broad N IV λ1486 emission and a significant P-Cygni line profile in N IV λ1719 (e.g., Crowther et al. 2016; Martins & Palacios 2022), are observed in the spectra of several UV-bright galaxies. While acknowledging that follow-up observations of the rest-optical blue and red bumps are necessary to definitively establish their VMS nature, as suggested by Martins et al. (2023), we classify sources with $EW_0 \ge 3.0$ Å (including J1220+0842, as discussed in Sects. 4 and 5) as candidates to host a significant number of VMS.

It is important to emphasize that our selection of UV-bright sources for GTC follow-up was not based on the presence of strong He II emission lines. Indeed, He II λ1640 is not detected or barely detected only in their SDSS spectra. Furthermore, the average $EW_0$ (He II) = $3.24 \pm 0.91$ Å found for our GTC sample is also similar to that of the SDSS stack composed of ≈70 UV-bright galaxies, where our sources were initially selected ($EW_0$ (He II) ≈ $3.5$ Å, Marques-Chaves, in prep.). This means that our GTC subsample is representative of UV-bright galaxies in general, and suggests that the incidence of VMS-dominated sources in UV-bright galaxies is fairly high, around ~60–70% (right panel of Fig. 8).

<footnote n="6">⁶ Here we consider VMS candidates based on sources bearing the labels “VMS” or “VMS or WR” in the classification scheme of Martins et al. (2023), as detailed in their Table 4, thereby indicating at least sources where VMS are suspected.</footnote>
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** EW$_0$ (He II) and incidence rate as a function of $M_{UV}$ for UV bright galaxies, along with other sources. Left: intensity of the He II λ1640 line for different types of sources: young star clusters and H II regions (green), normal or typical ($L_{UV}^*$) Lyman break galaxies (red), and the UV-bright galaxies studied in this work (blue). The R136 cluster and other VMS-source candidates are marked in yellow. Right: incidence of VMS in different galaxy populations: normal LBGs in red and UV-bright galaxies in blue.</figcaption></figure>

A different situation is found in normal/typical ($L_{UV}^*$) star-forming galaxies. While He II λ1640 is recurrently observed in $L_{UV}^*$ galaxies, its strength is substantially weaker than the EW$_0$ = 3.0 Å threshold used to differentiate the VMS or WR contributions. As already discussed in previous sections, the composite spectrum of Shapley et al. (2003) from almost 1000 individual spectra of z ~ 3 LBGs shows EW$_0$ ≈ 1.4 Å, and has been suggested to be due to WR stars (Brinchmann et al. 2008; Eldridge & Stanway 2012). In addition, the composite spectrum of $L_{UV}^*$ LBGs does not show any hint of broad emission in N IV λ1486 nor N IV λ1719 P-Cygni line profiles. This suggests that the average contribution of VMS in $L_{UV}^*$ LBGs is likely marginal. These findings are corroborated by the works of Cassata et al. (2013), Nanayakkara et al. (2019), and Saxena et al. (2020), which provide measurements of EW$_0$ (He II) for individual z ~ 2–5 $L_{UV}^*$ LBGs using ultra-deep spectroscopy. For example, among nearly 950 sources observed within VANDELS, Saxena et al. (2020) find intense (EW$_0$ ≥ 3.0 Å) and broad (FWHM > 1000 km s⁻¹) emission in the 1640 Å line for only three of them (Fig. 8, left). This suggests that the relative number of VMS over OB stars in $L_{UV}^*$ star-forming galaxies is low or negligible, with a VMS incidence below <1% (right panel of Fig. 8).

Understanding the underlying factors that account for the prevalence of VMS in UV-bright galaxies but their rarity in $L_{UV}^*$ LBGs is certainly important, albeit challenging. Various factors, such as differences in metallicity, age, and star-formation histories between these two classes of sources, may potentially influence our results. We note, however, that metallicity estimations for $L_{UV}^*$ LBGs are generally similar to those inferred for our sources (e.g., Pettini et al. 2001; Steidel et al. 2014). Furthermore, LBGs are, by selection, actively star-forming galaxies, although they can be substantially older than our UV-bright galaxies (≈10 Myr). Nevertheless, even considering continuous star formation over 100 Myr, synthesis models with VMS predict EW$_0$ (He II) ≈ 2.0 Å, which still surpasses the EW$_0$ ≈ 1.4 Å observed in $L_{UV}^*$ LBGs. Therefore, a more natural explanation could be due to intrinsic differences in the IMF, which may enhance the formation of VMS in UV-bright galaxies.

Recent results do indeed suggest that the IMF can grow toward top-heavy with increasing density and decreasing metallicity (e.g., Marks et al. 2012; Haghi et al. 2020; Weatherford et al. 2021). While the metallicity argument may be not valid for our UV-bright galaxies, they appear very compact considering their high stellar masses and SFRs. At least two of these sources (J0121+0025 and J1220+0842) have ground-based optical imaging with very good seeing conditions, for which Marques-Chaves et al. (2020b, 2021) inferred characteristic sizes of $r_{eff}$ ≲ 1 kpc. This results in very high stellar mass and SFR surface densities of log(Σ$_M$/M$_⊙$ pc⁻²) > 2.8 and log(Σ$_{SFR}$/M$_⊙$ yr⁻¹ kpc⁻²) > 2.0, which are substantially higher than those found typically in $L_{UV}^*$ LBGs at similar redshifts (e.g., Shibuya et al. 2015).

## 6.2. Possible caveats

Since standard methods/calibrations used to derive physical parameters do not include the effects of VMS one may wonder how accurate, for example, the SFR values listed here are. Comparing our models with/without VMS, we find that SFR(UV) could be reduced by ~10–30% for ages between ~10–100 Myr and maximum stellar masses up to 400 M$_⊙$, since VMS boost the UV luminosity but not by very large factors in the case of extended (constant) star-formation. This uncertainty is less than the age-dependence of the SFR(UV) conversion factor $\kappa_{UV}$ over the same age interval (which is a factor ~1.7). A more detailed discussion of this and other effects of VMS will be presented elsewhere.

We discuss if chemical enrichment from VMS overestimates the metallicities of our sources. Higgins et al. (2023) and Vink (2023) have recently suggested that VMS could significantly pollute the ISM, also in nitrogen, which is used here to determine the metallicity from the N2 indicator. If correct, this would imply that the true stellar metallicity would be lower than inferred. In this case, the strong He II λ1640 emission would probably be even more difficult to explain with normal stellar populations since at least the emission from WR stars is known to diminish with decreasing metallicity. The inference about the exact VMS content would then depend on the metallicity-dependence
of VMS, both on their stellar evolution and atmosphere prop-
erties, which are essentially unknown, although Smith et al.
(2023) suggest that VMS signatures could be similar between
LMC metallicity and ~0.15 solar. If the wind density of VMS
decreases with metallicity their He II emission should be weaker,
and hence the amount of VMS required to explain the observed
emission line would be higher than in the present models. In
short, if the metallicity of our sources was significantly lower
than that of the LMC, we would probably underestimate the
amount of VMS.

Another limitation of the present work is that we have no good measure of the IMF, its slope, and maximum at the high mass end. For simplicity, we have assumed that the classical Salpeter slope extends to higher masses. Constraining independently the slope and upper mass limit appears difficult and is probably a degenerate problem. Also, an exact quantification of these properties will require better models of “normal” stellar populations, including a proper description of emission from WR stars, as discussed earlier. In any case, although the contribution of WR stars to He II λ1640 is uncertain, we think that the presence of VMS is quite securely established in the objects studied here.

Clearly, to progress on these issues and more firmly quantify
the contribution of VMS in our objects and in general, one needs
accurate metallicity measurements, ideally abundances of sev-
eral species to see if any peculiar abundance patterns are found,
deep optical spectra to search for the signatures of WR stars
which can be distinguished from VMS (cf. Martins et al. 2023),
improved spectral synthesis models of normal stellar popula-
tions, and VMS models at different metallicities, as also noted,
for example, by Senchyna et al. (2021) and Smith et al. (2023).

7. Conclusions

In this work, we have investigated the presence of very mas-
sive stars (VMS > 100 M⊙) in 13 UV-bright star-forming galax-
ies using deep GTC optical and near-IR spectroscopy. These
galaxies, with redshifts between 2.2 ≲ z ≲ 3.6, are among the
UV-brightest sources known at high redshift, with UV abso-
lute magnitudes ranging between −23.30 to −24.70. They also
present very large SFRs ≈ 100–1000 M⊙ yr⁻¹ and nebular metal-
licities of 12 + log(O/H) = 8.10–8.50, with a mean value of
≈8.36. We have analyzed the integrated rest-frame UV spectra
with GTC using empirical templates and population synthesis
models with and without VMS. From the analysis of these data,
we obtained the following results:

- The very high S/N rest-frame UV spectra reveal intense and broad He II λ1640 emission for all sources, with EW₀ (He II) between ≈1.40 Å and ≈4.60 Å and line widths of ≈1200–3200 km s⁻¹ (FWHM). These sources exhibit diverse He II spectral profiles, from nearly Gaussian/symmetric profiles to asymmetric profiles, or more complex ones with multiple emission/absorption peaks. We find a tentative correlation between (O/H) and EW₀ (He II) so that stronger He II λ1640 emission is found predominantly at higher metallicities. The rest-frame UV spectra also show other spectral features originating from very young stellar populations, such as strong P-Cygni line profiles in the wind lines N V λ1240, Si IV λλ1393, 1401 and C IV λλ1548, 1550, indicating very young ages on the order of ~10 Myr assuming continuous star-formation histories.

- We compare the GTC spectra and the He II λ1640 pro-
files of our UV-bright galaxies with those of known VMS-
dominated sources and other empirical spectra of typical

galaxies with normal stellar populations, that is, without
VMS. We find that the rest-UV spectra of some of the UV-
bright galaxies closely resemble those of VMS-dominated
clusters, like the local R136/LMC and SB179 clusters or
the Sunburst cluster at z = 2.37, including the strength and
spectral shape of the He II λ1640 line. On the other hand,
the spectra of UV-bright galaxies differ significantly from
those of typical (L<sup>*</sup><sub>UV</sub>) galaxies where only normal mas-
sive stars are expected (i.e., <100 M<sub>☉</sub> including WR stars).
Furthermore, the spectra of UV-bright galaxies also reveal
other VMS signatures, such as N IV λ1486 emission and
N IV λ1719 P-Cygni line profiles, providing additional evi-
dence for the presence of VMS in these sources.

- We also compare the strengths of the observed He II λ1640 emission and the absorption component of the N V λ1240 and C IV λλ1548,1550 of our sources with those of synthetic population spectra. For that, we use both standard BPASS models with an IMF upper mass cutoff of 100 M⊙ and updated models incorporating VMS self-consistently (Martins & Palacios 2022) with upper mass cutoffs up to 475 M⊙. We find that the majority of UV-bright galaxies require a contribution of VMS to explain the observed strengths of He II, N V, and C IV. Therefore, our results suggest that UV-bright galaxies have a different IMF with upper mass limits between M<sub>up</sub> = 175–475 M⊙, assuming a Salpeter slope.

- Using an empirical threshold of EW₀ (He II λ1640) = 3.0 Å to differentiate VMS or WR contributions, along with the detection of other VMS spectral profiles (N IV λ1486 and N IV λ1719), we classify nine out of 13 UV-bright galaxies as VMS-dominated sources. This suggests that the incidence of VMS-dominated sources in the UV-bright galaxy population is high, around ≈70%, and is much higher than in typical L<sup>*</sup><sub>UV</sub> Lyman break galaxies at similar redshifts where the incidence of VMS appears negligible (<1%).

Acknowledgements. We thank the anonymous referee for critical reviews of the manuscript. Based on observations made with the Gran Telescopio Canarias (GTC) installed in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias, in the island of La Palma. We thank the GTC staff for their help with the observations. We thank Dr. Eros Vanzella for providing us with the spectrum of the Sunburst cluster. A.U. is grateful for support from the Warwick Astronomy Prize Studentship Fund. E.R.S. is supported in part by UK STFC grants ST/X001121/1 and ST/T000406/1.

References

Abolfathi, B., Aguado, D. S., Aguilar, G., et al. 2018, ApJS, 235, 42

Álvarez-Márquez, J., Marques-Chaves, R., Colina, L., & Pérez-Fournon, I. 2021, A&A, 647, A133

Amdar, L., Palacios, A., Charbonnel, C., et al. 2019, A&A, 631, A77

Bentz, M. C., Hall, P. B., & Osmer, P. S. 2004, AJ, 128, 561

Berry, M., Gawiser, E., Guaita, L., et al. 2012, ApJ, 749, 4

Bestenlehner, J. M., Vink, J. S., Gräfener, G., et al. 2011, A&A, 530, L14

Bestenlehner, J. M., Gräfener, G., Vink, J. S., et al. 2014, A&A, 570, A38

Bestenlehner, J. M., Crowther, P. A., Caballero-Nieves, S. M., et al. 2020, MNRAS, 499, 1918

Bian, F., Fan, X., Jiang, L., et al. 2012, ApJ, 757, 139

Björklund, R., Sundqvist, J. O., Singh, S. M., Puls, J., & Najarro, F. 2023, A&A, 676, A109

Bouret, J. C., Hillier, D. J., Lanz, T., & Fullerton, A. W. 2012, A&A, 544, A67

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

Brinchmann, J., Pettini, M., & Charlot, S. 2008, MNRAS, 385, 769

Bruhweiler, F. C., Miskey, C. L., & Smith Neubig, M. 2003, AJ, 125, 3082

Byrne, C. M., Stanway, E. R., Eldridge, J. J., McSwiney, L., & Townsend, O. T. 2022, MNRAS, 512, 5329

Cabanac, R. A., Valls-Gabaud, D., & Lidman, C. 2008, MNRAS, 386, 2065

Calabrò, A., Castellano, M., Pentericci, L., et al. 2021, A&A, 646, A39

Calzetti, D., Armus, L., Bohlin, R. C., et al. 2000, ApJ, 533, 682
Cardelli, J. A., Clayton, G. C., & Mathis, J. S. 1989, ApJ, 345, 245

Carnall, A. C. 2017, arXiv e-prints [arXiv:1705.05165]

Cassata, P., Le Fèvre, O., Charlot, S., et al. 2013, A&A, 556, A68

Chabrier, G. 2003, ApJ, 586, L133

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

Conti, P. S., Leitherer, C., & Vacca, W. D. 1996, ApJ, 461, L87

Crowther, P. A., & Dessart, L. 1998, MNRAS, 296, 622

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

Dessauges-Zavadsky, M., D'Odorico, S., Schaerer, D., et al. 2010, A&A, 510, A26

Eisenstein, D. J., Weinberg, D. H., Agol, E., et al. 2011, AJ, 142, 72

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

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

Feulner, G., Gabasch, A., Salvato, M., et al. 2005, ApJ, 633, L9

Garcia, M., & Bianchi, L. 2004, ApJ, 606, 497

Garilli, B., McLure, R., Pentericci, L., et al. 2021, A&A, 647, A150

González Delgado, R. M., Heckman, T., Leitherer, C., et al. 1998, ApJ, 505, 174

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

Haghi, H., Safaei, G., Zonoozi, A. H., & Kroupa, P. 2020, ApJ, 904, 43

Hainich, R., Rühling, U., Todt, H., et al. 2014, A&A, 565, A27

Hainline, K. N., Shapley, A. E., Greene, J. E., & Steidel, C. C. 2011, ApJ, 733, 31

Hamann, W. R., Gräfener, G., & Liermann, A. 2006, A&A, 457, 1015

Higgins, E. R., Vink, J. S., Hirschi, R., Laird, A. M., & Sabhahit, G. N. 2023, MNRAS, 526, 534

Humphrey, A., Villar-Martín, M., Vernet, J., et al. 2008, MNRAS, 383, 11

Jiang, L., Fan, X., & Vestergaard, M. 2008, ApJ, 679, 962

Jones, T., Stark, D. P., & Ellis, R. S. 2012, ApJ, 751, 51

Kalari, V. M., Horch, E. P., Salinas, R., et al. 2022, ApJ, 935, 162

Koushan, S., Driver, S. P., Bellstedt, S., et al. 2021, MNRAS, 503, 2033

Kudritzki, R. P., Springmann, U., Puls, J., Pauldrach, A. W. A., & Lennon, M. 1998, in Properties of Hot Luminous Stars, ed. I. Howarth, ASP Conf. Ser., 131, 299

Lagarde, N., Decressin, T., Charbonnel, C., et al. 2012, A&A, 543, A108

Lamers, H. J. G. L. M., & Cassinelli, J. P. 1999, Introduction to Stellar Winds (Cambridge: Cambridge University Press)

Le Fèvre, O., Tasca, L. A. M., Cassata, P., et al. 2015, A&A, 576, A79

Le Fèvre, O., Lemaux, B. C., Nakajima, K., et al. 2019, A&A, 625, A51

Leitherer, C., Ortiz Otálvaro, P. A., Bresolin, F., et al. 2010, ApJS, 189, 309

Leitherer, C., Tremonti, C. A., Heckman, T. M., & Calzetti, D. 2011, AJ, 141, 37

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

Lilly, S. J., Le Fèvre, O., Hammer, F., & Crampton, D. 1996, ApJ, 460, L1

López Fernández, R., González Delgado, R. M., Pérez, E., et al. 2018, A&A, 615, A27

Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415

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

Marks, M., Kroupa, P., Dabringhausen, J., & Pawlowski, M. S. 2012, MNRAS, 422, 2246

Marques-Chaves, R., Pérez-Fournon, I., Shu, Y., et al. 2017, ApJ, 834, L18

Marques-Chaves, R., Pérez-Fournon, I., Shu, Y., et al. 2020a, MNRAS, 492, 1257

Marques-Chaves, R., Álvarez-Márquez, J., Colina, L., et al. 2020b, MNRAS, 499, L105

Marques-Chaves, R., Schaerer, D., Álvarez-Márquez, J., et al. 2021, MNRAS, 507, 524

Marques-Chaves, R., Schaerer, D., Álvarez-Márquez, J., et al. 2022, MNRAS, 517, 2972

Martins, F. 2023, A&A, 680, A22

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

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

Martins, F., Schaerer, D., Marques-Chaves, R., & Upadhyaya, A. 2023, A&A, 678, A159

Massey, P., & Hunter, D. A. 1998, ApJ, 493, 180

Meštrić, U., Vanzella, E., Upadhyaya, A., et al. 2023, A&A, 673, A50

Nanayakkara, T., Brinchmann, J., Boogaard, L., et al. 2019, A&A, 624, A89

Noll, S., Mehlert, D., Appenzeller, I., et al. 2004, A&A, 418, 885

Pentericci, L., McLure, R. J., Garilli, B., et al. 2018, A&A, 616, A174

Pettini, M., Steidel, C. C., Adelberger, K. L., Dickinson, M., & Giavalisco, M. 2000, ApJ, 528, 96

Pettini, M., Shapley, A. E., Steidel, C. C., et al. 2001, ApJ, 554, 981

Plat, A., Charlot, S., Bruzual, G., et al. 2019, MNRAS, 490, 978

Quider, A. M., Pettini, M., Shapley, A. E., & Steidel, C. C. 2009, MNRAS, 398, 1263

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

Sánchez, S. F., Avila-Reese, V., Rodríguez-Puebla, A., et al. 2019, MNRAS, 482, 1557

Sanders, R. L., Shapley, A. E., Jones, T., et al. 2021, ApJ, 914, 19

Saxena, A., Pentericci, L., Mirabelli, M., et al. 2020, A&A, 636, A47

Schaerer, D., & Vacca, W. D. 1998, ApJ, 497, 618

Schlafly, E. F., & Finkbeiner, D. P. 2011, ApJ, 737, 103

Senchyna, P., Stark, D. P., Vidal-García, A., et al. 2017, MNRAS, 472, 2608

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

Shapley, A. E., Steidel, C. C., Pettini, M., & Adelberger, K. L. 2003, ApJ, 588, 65

Shenar, T., Sana, H., Crowther, P. A., et al. 2023, A&A, 679, A36

Shibuya, T., Ouchi, M., & Harikane, Y. 2015, ApJS, 219, 15

Sidoli, F., Smith, L. J., & Crowther, P. A. 2006, MNRAS, 370, 799

Siess, L., Dufour, E., & Forestini, M. 2000, A&A, 358, 593

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

Smith, L. J., Oey, M. S., Hernandez, S., et al. 2023, ApJ, 958, 194

Sommariva, V., Mannucci, F., Cresci, G., et al. 2012, A&A, 539, A136

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

Steidel, C. C., Rudie, G. C., Strom, A. L., et al. 2014, ApJ, 795, 165

Steidel, C. C., Strom, A. L., Pettini, M., et al. 2016, ApJ, 826, 159

Vernet, J., Fosbury, R. A. E., Villar-Martín, M., et al. 2001, A&A, 366, 7

Vink, J. S. 2023, A&A, 679, L9

Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 2001, A&A, 369, 574

Walborn, N. R., Nichols-Bohlin, J., & Panek, R. J. 1985, NASA Ref. Publ., 1155

Weatherford, N. C., Fragione, G., Kremer, K., et al. 2021, ApJ, 907, L25

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

Wofford, A., Sixtos, A., Charlot, S., et al. 2023, MNRAS, 523, 3949
**Appendix A: GTC spectra of individual sources**

Figure A.1 and A.2 show the GTC optical spectra of the 13 sources located at redshift (z) ~ 2.2–3.6 providing the rest-frame UV coverage.

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

**Fig. A.1.** GTC spectra of the six out of 13 UV-bright galaxies studied in this work are in black with the $1\sigma$ uncertainty in gray. The calibrated spectra have a resolution of around, R ~ 700. The x-axes represent the rest frame and observed wavelengths (bottom and top axis, respectively). The three signature VMS profiles, He II λ1640 emission, N IV λ1486 emission, and N IV λ1719 P-Cygni line profile are highlighted in gray shaded regions.
<figure><img src="image_11.png" /><figcaption>**Fig. A.2.** GTC spectra of the seven out of 13 UV-bright galaxies studied in this work are in black with the $1\sigma$ uncertainty in gray. The calibrated spectra have a resolution of around, R ~ 700. The x-axes represent the rest frame and observed wavelengths (bottom and top axis, respectively). The three signature VMS profiles, He II λ1640 emission, N IV λ1486 emission, and N IV λ1719 P-Cygni line profiles are highlighted in gray shaded regions.</figcaption></figure>
## Appendix B: Population synthesis by extrapolating the IMF

Martins & Palacios (2022) uses an approach IMFs have been extrapolated to different upper mass limits by adding the SEDs of BPASS 100 $M_\odot$ stars with that of the VMS SEDs from Martins & Palacios (2022). In this approach, the number of VMS is calculated in the following mass bins: [100 $M_\odot$, 175 $M_\odot$], [175 $M_\odot$, 225 $M_\odot$], [225 $M_\odot$, 275 $M_\odot$], [275 $M_\odot$, 300 $M_\odot$]. These bins were chosen to add the SEDs of VMS which are only available for 150, 200, 250, and 300 $M_\odot$ stars. The following equation prescribed in Stanway & Eldridge (2018) has been used to calculate the number of stars between the mass bin 0.1 $M_\odot$ and a maximum mass $M_2$,

$$N(0.1, M_2) = C \times \left( \int_{0.1}^{M_1} M^{\alpha_1} dM + M_1^{\alpha_1} \int_{M_1}^{M_2} M^{\alpha_2} dM \right).$$

Upon closer inspection, we noticed that the equation had been calibrated only for a Chabrier IMF. In the case of Salpeter IMF, there is a discontinuity in the IMF slope transition at 0.5 $M_\odot$ which has been shown in Figure B.1. The discontinuity arises because of the change in slope of the IMF while transitioning from the lower mass range in 0.1 to 0.5 $M_\odot$ to a higher mass range in 0.5 $M_\odot$ to $M_2$ regime. To produce a continuous IMF in the Salpeter form, the equation should take the following form:

$$N(0.1, M_2) = C \times \left( \int_{0.1}^{M_1} M^{\alpha_1} dM + M_1^{\alpha_1-\alpha_2} \int_{M_1}^{M_2} M^{\alpha_2} dM \right).$$

The total mass in the mass bin (0.1, $M_2$) is given by

$$M(0.1, M_2) = C \times \left( \int_{0.1}^{M_1} M^{(1+\alpha_1)} dM + M_1^{(\alpha_1-\alpha_2)} \int_{M_1}^{M_2} M^{(1+\alpha_2)} dM \right)$$

$$M(0.1, M_2) = C \times \left( \frac{M_1^{(2+\alpha_2)}}{(2+\alpha_2)} - \frac{0.1^{(2+\alpha_2)}}{(2+\alpha_2)} + \frac{M_1^{\alpha_1-\alpha_2}}{(2+\alpha_2)} [M_2^{2+\alpha_2} - M_1^{2+\alpha_1}] \right).$$

We have used the new equation calibrated for Salpeter IMF to calculate the number of VMS in different mass bins. We calculate the normalization term C in each case when the IMF upper mass has been extended to different upper mass limits.

In this work, we have defined the following mass bins: [100 $M_\odot$, 175 $M_\odot$], [175 $M_\odot$, 225 $M_\odot$], [225 $M_\odot$, 275 $M_\odot$], [275 $M_\odot$, 325 $M_\odot$], and [325 $M_\odot$, 475 $M_\odot$] for which we calculate the normalization constant C each time. We then calculate the total mass in these mass bins [$M_a, M_b$] using the above equation calibrated for Salpeter IMF. The number of stars in the mass bin [$M_a, M_b$] is given by the following equation:

<figure><img src="image_12.png" /><figcaption>Fig. B.1. Visual illustration to show the discontinuity in the IMF when using the uncalibrated equation for the Salpeter IMF used in Martins & Palacios (2022) and taken from Stanway & Eldridge (2018) (in black). The discontinuity goes away when the equation is calibrated for the Salpeter IMF (in green) which is the form used in this work.</figcaption></figure>

$$N(M_a, M_b) = C \times \left( \frac{\frac{M_a^{(\alpha_1-\alpha_1)}}{2+\alpha_2} [M_b^{(2+\alpha_2)} - M_a^{(2+\alpha_2)}]}{\frac{M_a + M_b}{2}} \right).$$

This gives us 153.52 number of 150 $M_\odot$ VMS in the mass bin [100 $M_\odot$, 175 $M_\odot$], 44.41 number of 200 $M_\odot$ VMS in the mass bin [175 $M_\odot$, 225 $M_\odot$], 26.04 number of 250 $M_\odot$ VMS in the mass bin [225 $M_\odot$, 275 $M_\odot$], 16.85 number of 300 $M_\odot$ VMS in the mass bin [275 $M_\odot$, 325 $M_\odot$], and 34.45 number of 400 $M_\odot$ VMS in the mass bin [325 $M_\odot$, 475 $M_\odot$].

We note that, the SEDs of stars up to 100 $M_\odot$ from BPASS are calibrated for a total stellar mass of $10^6 M_\odot$. We add these SEDs to the SEDs of VMS from Martins & Palacios (2022) to get the new VMS models. This implies in our new VMS models, the luminosity contribution from stars up to 100 $M_\odot$ is slightly overestimated. In our new VMS models, the VMS contributes around 2.3% to 5.7% of the total mass of $10^6 M_\odot$. To incorporate this change in mass to SEDs, the SEDs of lower mass stars ($M < 100 M_\odot$) need to be calibrated by the same percentage of change in mass. As most of these masses will be distributed to the low mass end of the lower mass stars ($M < 100 M_\odot$) because of the shape of the IMF, we expect little change in overall SED contributions. In any case, the next-generation population synthesis models that will use the Gräfener (2021) mass loss recipe for the VMS to produce SEDs by incorporating all the mass ranges in their evolution, would be able to solve this discrepancy.