© The Authors 2025

1. Introduction

Stars with masses of M_⋆ ∼ 10 − 100 M_⊙ (generally referred to as massive stars) are ubiquitous in the nearby and distant Universe. They power H II regions and star-forming galaxies and their respective radiative output. Together with supernovae, massive stars dominate the mechanical feedback and chemical evolution in star-forming galaxies during the first gigayear of their evolution (Conti et al. 2008). In recent years, the existence of stars with initial masses higher than 100 M_⊙ – so-called very massive stars (VMS, Vink et al. 2015) – has been established or suggested, both in the Local Group and at cosmological distances, but their general impact on star-forming galaxies remains largely to be explored. This is one of main aims of the present work.

The best evidence for individual VMS comes from the well-studied region R136 at the core of the giant H II region 30 Dor in the Large Magellanic Cloud, which is resolved into hundreds of massive OB stars and several peculiar stars (with spectral types WNh), whose masses have now been established between 100–200 M_⊙ (e.g., Crowther et al. 2010; Bestenlehner et al. 2020; Kalari et al. 2022). In R136, a very young region (age ∼ 1 ± 1 Myr), the seven VMS produce a strong He IIλ1640 emission line (with EW ∼ 4 Å), and they are found to contribute 32% of the far-UV continuum (18% of the entire Tarantula nebula) and ∼25% of the ionizing photon production (Crowther et al. 2016; Bestenlehner et al. 2020; Crowther & Castro 2024).

The presence of VMS is also fairly well established or suspected in several nearby star-forming regions, young clusters, or galaxies, as is testified by their strong He IIλ1640 emission and other indirect evidence (Leitherer et al. 2018; Smith et al. 2016, 2023; Senchyna et al. 2021; Wofford et al. 2023). A recent search and critical analysis of VMS in nearby objects and low-z galaxies has been presented by Martins et al. (2023).

At cosmological distances, signatures of VMS have possibly been found in a gravitationally lensed star cluster in the Sunburst arc at z = 2.37 (Meštrić et al. 2023; Rivera-Thorsen et al. 2024), and Upadhyaya et al. (2024) have found evidence for VMS in very UV-bright (M₁₅₀₀ ∼ −23 to −24.5) compact star-forming galaxies discovered in the Sloan Digital Sky Survey by Marques-Chaves et al. (2020, 2021, 2022). Interestingly, Upadhyaya et al. (2024) find that galaxies showing signatures of VMS are frequent (found in ∼69%) of UV-bright galaxies, whereas the presence of VMS is quite rarely seen otherwise (e.g., in ≲1% of LBGs at z ∼ 2 − 5). This suggests that the presence of VMS, and an hence an extension of the IMF beyond the classical upper mass limit of M_up ∼ 100 − 120 M_⊙, might be more common or even frequent in high redshift star-forming galaxies.

In parallel, early observations with the JWST have revealed an unexpectedly large number of UV-luminous galaxies with blue spectral slopes in the early Universe (z ≳ 5 − 7, see e.g., Topping et al. 2024; Cullen et al. 2024). To explain this tension, several physical scenarios have been proposed, including dust clearing by radiation-driven outflows (Ferrara et al. 2023), a higher star formation efficiency, possibly due to reduced feedback (Dekel et al. 2023), stochastic star formation histories that can lead to biases boosting the UV luminosity (Mason et al. 2023; Shen et al. 2023), and different IMFs favoring massive stars (Finkelstein et al. 2023; Harikane et al. 2023; Trinca et al. 2024; Woodrum et al. 2023). For all these reasons, examining the role and impact of VMS and the effect of different IMFs on the observable properties of star-forming galaxies is of interest in a variety of contexts. This is the main goal of the this paper.

To model the spectral evolution of stellar populations including VMS, it is important to consistently compute stellar evolution and atmosphere models with appropriate prescriptions for stellar mass loss and wind properties, and, ideally, to test and/or “calibrate” the models against observations. Several synthesis models already allow for extensions of the IMF beyond 100 M_⊙, into the regime of VMS. For example, both the BPASS and the updated Charlot & Bruzual models contain predictions for IMFs up to 100 and 300 M_⊙ (see Eldridge et al. 2017; Bruzual & Charlot 2003; Plat et al. 2019; Mayya et al. 2023), and references therein). However, both model sets show difficulties in reproducing several UV observations and fail to reproduce in particular the most important spectral diagnostic of VMS, the He IIλ1640 emission, as has been shown by Senchyna et al. (2021), Martins et al. (2023), and Smith et al. (2023). The main reason for this failure is most likely the lack of appropriate atmosphere models for VMS, since the models cited before treat these objects with classical WM-BASIC models developed for O stars. Combining stellar evolution and non-LTE atmosphere models including strong stellar winds, Martins & Palacios (2022) have been able to reproduce the spectra of individual VMS in the LMC, the integrated spectrum of R136, and they also reproduce UV spectra of other clusters with strong He II features. These VMS models, combined with BPASS for stars ≤100 M_⊙, have also been successfully compared to the observations of the Sunburst arc cluster at z = 2.37 and to a sample of UV-bright galaxies at z ∼ 2 − 4, which show UV spectra indicative of the presence of VMS (Meštrić et al. 2023; Upadhyaya et al. 2024).

In this work, we use these VMS models to explore in more detail the UV, ionizing, and related observable properties of stellar populations including stars in the range of ∼100 − 400 M_⊙. We compute in particular the UV and Lyman continuum emission of such populations, including the associated nebular emission (primarily nebular continuum and H recombination lines), for a range of IMFs, covering different upper and lower mass limits and two different IMF slopes (Salpeter and a flatter IMF). We then show how the UV luminosity, the ionizing photon production and efficiency, line equivalent widths, and other properties are “boosted” by the presence of VMS and how they depend on the IMF. We also examine the shape of the ionizing spectra of populations including VMS, showing that they have smaller Lyman breaks and softer spectra than normal populations. To understand what maximum predictions are possible, e.g. for the UV luminosity per stellar mass or star formation rate, for the ionizing photon efficienciency, ξ_ion, and others, we also consider integrated populations completely devoid of normal stars (i.e., populations of pure VMS), and we show that most observational properties are naturally bounded by some maximum values, which we derive. We also briefly discuss cases dominated by pure nebular continuum emission. The model predictions presented here are made available publicly and should hopefully contribute to a better understanding of the importance of VMS and effects of unusual IMFs in star-forming galaxies.

In Sect. 2, we briefly describe the inputs of our model calculations. The predictions of our model grids are shown and discussed in Sect. 3. In Sect. 4, we compare our work with other studies, discuss the maximum UV and ionizing photon efficiencies of populations including VMS, present observational tests for extreme IMFs, and compare our results to recent observation from JWST. Our main results are summarized in Sect. 5.

2. Models

To study the impact of VMS on the integrated properties of star-forming regions and galaxies we use the recent combined stellar evolution and atmosphere models of Martins & Palacios (2022) at 1/2.5 solar metallicity. The models include the specific mass loss rates of VMS, that is mass loss rates that are boosted compared to normal O-type stars (Gräfener & Hamann 2008; Vink et al. 2011; Bestenlehner 2020) and that are the main ingredient leading to the special spectroscopic appearance of VMS. The mass loss recipe for VMS was adopted from Gräfener (2021). It relies on a calibration anchored on the observed wind properties of the VMS in R136 at the core of the Tarantula nebula in the LMC. Stellar tracks were computed for masses between 150 and 400 M_⊙ with the code STAREVOL (Amard et al. 2016). Atmosphere models and synthetic spectra were subsequently calculated along each track, at positions corresponding to ages between 0 and 2.5 Myr (with steps of 0.5 Myr). The atmosphere models were computed with the code CMFGEN (Hillier & Miller 1998) that include the physics relevant for massive stars (winds, spherical geometry, line-blanketing, non-local thermodynamical equilibrium). The evolution and atmosphere models are consistent in terms of stellar properties (effective temperature, mass, radius, luminosity, surface chemical composition): the output of the evolutionary calculations are taken for input in the atmosphere calculations. The final synthetic spectra sample the main sequence of VMS. The optical spectra resemble those of O supergiants with strong winds (OIf stars), transition objects (OIf/WN), or hydrogen-rich WN (WNh) stars. In the UV, the spectra show the classical P-Cygni emission lines (including C IVλ1550) and several features unique to VMS, in particular strong broad He IIλ1640 emission (Martins & Palacios 2022). The late stages of stellar evolution, corresponding to about 10% of the star’s lifetime, were not considered for the reasons explained in Martins & Palacios (2022).

The VMS models are added to “normal” stars, defined here as stars with initial masses ≤100 M_⊙, making use of the recent BPASS evolutionary synthesis models of Eldridge et al. (2017), following the simple procedure adopted by Martins & Palacios (2022) and updates described in Upadhyaya et al. (2024). In practice, we use the BPASS models v.2.2.1, including binary stars for a metallicity Z = 0.006, the closest available to Z = 0.00536 = 1/2.5 Z_⊙ for Z_⊙ = 0.0134 adopted by Martins & Palacios (2022). We explore different stellar initial mass functions, including the “classical” IMF, described by a Salpeter slope of α₂ = −2.35 for masses M_⋆ ≥ 0.5 M_⊙ and α₁ = −1.3 at low masses (from the lower mass limit of M_low = 0.1 − 0.5 M_⊙) and a flatter IMF (α₂ = −2.0) giving more weight to massive stars. The maximum stellar mass considered varies between M_up = 100 and 400 M_⊙, and we also consider extreme cases with populations formed exclusively with very massive stars; that is, M_low = 100 M_⊙ (VMS-only in the following). The different models used in this study and their parameters are summarized in Table 1.

Since nebular emission is significant for populations including young and massive stars we include nebular continuum emission (free-bound, free-free, and two-photon) from H and He and added to the stellar spectra. We discuss our assumptions for the nebular continuum below. Finally, spectral energy distributions (SEDs) are computed for two limiting star formation histories, namely instantaneous bursts and constant star formation rates (SFR). Our models are normalized to a total stellar mass of 10⁶ M_⊙ for bursts or to SFR = 1 M_⊙ yr⁻¹. A selection of the main quantities predicted by our models are tabulated in Table 2 below. Detailed model predictions, including SEDs, are available on request from the first author and made available electronically.

To compute nebular continuum emission (free-bound, free-free, and two-photon) from H and He we use pyneb (Morisset et al. 2020), assuming assume electron temperatures and densities (T_e = 15 000 K, n_e = 100 cm⁻³) typical of high redshift moderately low-metallicity sources (e.g., Sanders et al. 2020; Curti et al. 2023). For significantly lower metallicities (≲0.03 − 0.1 solar) than adopted here, the nebular continuum can in principle be enhanced compared to case B assumption by a factor proportional to the mean energy of the ionizing photons, due to enhanced photoionization rates as shown by Raiter et al. (2010) and discussed recently by Katz et al. (2024). Finally, the effect of varying electron temperatures are illustrated in the Appendix (see Fig. A.1). Significant increases in the electron density beyond n_e ≳ 10^3 − 4 cm⁻³ decrease the two-photon emission, in other words nebular emission in the rest-UV, but have little effect at longer wavelengths. For a more detailed discussion of these effects, which may thus affect individual objects especially at very low metallicity we refer, for example, to Bottorff et al. (2006), Raiter et al. (2010), Katz et al. (2024). However, since only few galaxies observed with JWST show metallicities ≲0.03 − 0.1 solar, we do not expect significantly enhanced two-photon emission in the majority of them, and hence consider our assumptions should indeed be quite representative of “typical” objects with the moderately low-metallicity, adopted here.

One of the limitations of this work originates from the availability of detailed stellar evolution and atmosphere models for VMS, which are so far essentially only available for a single metallicity. In this sense our work should be considered as a pilot study, before a broader metallicity range can be explored. We acknowledge that our models are tailored for the LMC metallicity. VMS are mostly characterized by their strong stellar winds that are only calibrated at that metallicity. There is so far no direct information about the behavior of VMS winds or calibration of VMS mass loss rates at lower metallicity. Some empirical studies argue that VMS mass loss rates do not change with Z (Smith et al. 2023), while some theoretical prediction of H-free WR stars call for complex metallicity dependence (Sander & Vink 2020). The spectroscopic evolution of VMS at low Z will be the topic of a subsequent publication (Martins et al., in prep.). We briefly discuss this issue further in Sect. 4.4.

3. Predicted UV and ionizing properties of populations including VMS

To illustrate the impact of VMS on the UV–optical SED of integrated stellar populations, we show some examples from our models at an age of 1 Myr in Fig. 1. In the left panel, which compares the SED of a populations with a Salpeter IMF extending up to M_up = 100 and 400 M_⊙, one can clearly see the main changes due to VMS: First, VMS boost the UV luminosity of stellar populations, both in the Lyman continuum (LyC, at λ < 912 Å) and in the non-ionizing UV (λ > 912 Å). Second, the Lyman break of population including VMS is smaller, in other words the ratio of LyC to non-ionizing flux is increased. Third, the ionizing spectrum of populations with VMS decreases more rapidly between 350 and 228 Å, which means that their ionizing spectra are softer between ∼35 and 54 eV. These properties will be discussed further and quantified in what follows.

Fig. 1. Comparison of predicted SEDs for a 1 Myr old population with and without VMS. Left panel: The standard BPASS model with Mup = 100 M⊙ is shown in gray (stellar emission only) and in black (including nebular continuum emission). Blue and red lines show the model extending up to 400 M⊙, the former including only stellar emission, the latter the total continuum (stellar plus nebular continuum). Right panel: zoom on the UV domain, comparison the standard BPASS model with a population including only VMS (from 100 to 400 M⊙). The dotted line shows the nebular continuum emission of the VMS-only-400 population.

The right panel of Fig. 1 shows the SED of the VMS-only model, populated exclusively with masses between 100 and 400 M_⊙, compared to the model with a classical IMF with M_up = 100 M_⊙, where in both cases the total stellar mass is identical (10⁶ M_⊙, by construction). In this case the UV luminosity boost is very high – more than a factor 10 at most wavelengths – and the Lyman break of the VMS-only-400 model is essentially absent. In the right panel we also show separately the nebular continuum, which obviously provides a significant contribution longward of Lyα (at λ > 1216 Å), as well known for populations with strong Lyman continuum emission (see e.g., Schaerer 2002, 2003; Schaerer & de Barros 2012).

3.1. Ionizing photon flux and recombination line strengths

The Lyman continuum photon production, Q_H, which measures the number of H-ionizing photons emitted per unit time, of the integrated populations including VMS is shown in Fig. 2 as a function of time and for the different maximum masses M_up considered. Overall we see that an extension of the IMF beyond 100 M_⊙ leads, as expected, to a stronger LyC flux per unit stellar mass or SFR. The increase in Q_H is up to ∼0.4 dex (a factor ∼2.5) at young ages when stars up to 400 M_⊙ are included. At equilibrium, in other words at ages ≳5 Myr, the increase is of ∼0.2 dex (factor ∼1.6). We note also that the increase predicted by the calculations including the VMS models of Martins & Palacios (2022) is somewhat larger than that predicted by the BPASS models which extend to 300 M_⊙. This is due to differences in the stellar evolution and atmosphere models used by these works.

Fig. 2. Predicted temporal evolution of the ionizing photon production, QH, of different models, for instantaneous bursts and constant SFR. Burst models are normalized to 106 M⊙, models for constant SFR to SFR = 1 M⊙ yr−1.

One of the immediate implications of the increased ionizing photon emission Q_H is that SFR indicators using H recombination lines (e.g., Hα) must be revised if VMS are present. Since SFR(Hα) = c_Hα × L(Hα) and L(Hα)∝Q_H (see also Sect. A), the conversion factor c_Hα is simply increased by the amounts discussed above.

The predicted strength of the Hβ recombination line, as measured by its equivalent width, EW(Hβ), is shown in Fig. 3. It shows a well-known monotonic decrease with age due to the combined decrease in the ionizing photon flux and the increasing flux in the optical domain. A maximum EW(Hβ) ≈ 700 Å is reached by models including VMS up to 400 M_⊙, a factor 1.3 higher than for M_up = 100 M_⊙ in the models accounting for nebular continuum emission. This “boost” is less than the increase in Q_H, since EW ∝ Q_H/L₄₈₆₁ and since both Q_H and the continuum luminosity L₄₈₆₁ at Hβ wavelength increase in the presence of VMS.

Fig. 3. Predicted temporal evolution of the Hβ equivalent width of different models, for instantaneous bursts and constant SFR. Solid (dashed) lines show models without (with) nebular continuum emission.

Indeed, as well known the nebular continuum contributes a significant fraction of the flux in the optical domain for young populations (see Fig. 1). If this contribution is artificially left out, the Hβ equivalent width increases by ∼60 − 70%, as shown by the models with no nebular continuum in Fig. 3.

3.2. Properties of the SEDs in the Lyman continuum

Both the strength of the Lyman break and the detailed shape of the SEDs in the Lyman continuum are altered by the presence of VMS as illustrated in Fig. 1. To quantify the strength of the Lyman break we show the ratio of the ionizing to non-ionizing continuum flux at 1500 Å, L₉₀₀/L₁₅₀₀ measured in F_ν units, in Fig. 4. In models including VMS the LyC flux increases, in other words the break decreases, by factors up to ∼1.4 at the youngest ages and ∼1.5 for constant SFR. Accounting for nebular continuum emission, these effects lead to a global increase in L₉₀₀/L₁₅₀₀ by no more than 30% when VMS are included. Earlier predictions of the break (L900/l1500), derived from BPASS or other models, are shown, for example, in Schaerer (2003), Siana et al. (2007), Steidel et al. (2018) and they overall agree with the standard BPASS models excluding VMS shown here.

Fig. 4. Predicted temporal evolution of the L900/L1500 ratio of different models, for instantaneous bursts and constant SFR. Solid (dashed) lines show models without (with) nebular continuum emission.

To illustrate the shape of the ionizing radiation fields predicted by the different models we show their ionizing photon fluxes above a certain energy, Q(> E), as a function of energy in Fig. 5. All models are normalized at the Lyman limit (13.6 eV) for comparison, and the ionization potential of several important ions are indicated in the figure. Clearly the SED shapes of all models are essentially undistinguishable between 13.6 and ∼35 eV and they start to diverge at higher energies. The surprising result, a priori, is that the ionizing spectrum of stellar populations including VMS is softer and not harder as could naively be expected when more massive hot stars are added. This behavior is due to physical processes in the atmospheres of these stars, where the presence of strong (dense) stellar winds causes a different ionization structure resulting in softer spectra. In particular He²⁺ recombines in the outer parts of these atmospheres, similarly to the case of WR stars discussed already by Schmutz et al. (1992), leading, for example, to strong and broad He IIλ1640 emission, which is in some way the main observational characteristics of VMS. Therefore, most ionizing photons above 54 eV are absorbed in the atmosphere of these stars, which explains the strong decrease of the ionizing flux at these energies.

Fig. 5. Comparison of the ionizing spectra of different models for age 1 Myr. Shown is the cumulative number of ionizing photons emitted above energy E, Q(> E), as a function of the energy (in eV). All ionizing spectra are normalized to unity at the Lyman limit (13.6 eV), to allow for easy comparison of their relative hardness. Vertical lines mark the ionization potentials of several important ions.

Differences in the shape/hardness of the ionizing radiation field have implications for the relative strengths of nebular emission lines. Since the SEDs are essentially identical for energies below 35 eV the relative intensities of main strong emission lines in the optical range (e.g., [O II] λ3727, [O III] λλ4959,5007, [N II] λλ6548,6584, [S II] λλ6717,6731), which have ionization potentials between 13.6 and 35 eV, are expected to be basically invariant with respect to these differences. Only for lines of higher ionization stages, such as C IV, N IV, O IV, and He II one expects changes, namely lower intensities compared, for example, to hydrogen lines or to other lines from lower ionization stages. We have indeed confirmed this qualitative behavior with calculations of photoionization models using CLOUDY (Ferland et al. 2017), which show that the strong optical line ratios are basically unchanged by the inclusion of the VMS models. Higher ionization lines are predicted to be weak in all models, with intensities most likely below the detection threshold of most observations. Therefore, the predicted changes due to VMS are probably difficult to detect, and we defer a more thorough analysis of the impact of VMS on nebular emission line properties to later.

3.3. Ionizing photon efficiency

An important quantity measuring the ionizing photon emission of star-forming galaxies is the so-called ionizing photon efficiency, ξ_ion = Q_H/L_ν, 1500, which measures the ionizing photon flux per unit monochromatic UV luminosity at 1500 Å, L_ν, 1500, in units of erg s⁻¹ Hz⁻¹. The predictions for the different models considered here are shown in Fig. 6, both including and excluding nebular continuum emission. At young ages, ξ_ion is very high and the presence of VMS boosts the ionizing photon efficiency by up to a factor 1.4–1.5 (0.15–0.17 dex, for M_up = 400 M_⊙). For constant SFR at equilibrium (ages ≳ 300 Myr) ξ_ion increases by up to ∼0.2 dex when VMS are included, reaching log(ξ_ion) ≈ 25.4 erg⁻¹ Hz, which is a factor 60% higher than the “canonical” value used by Robertson et al. (2013), derived from synthesis models from Bruzual & Charlot (2003).

Fig. 6. Ionizing photon efficiency ξion as a function of age (left panel) and Hα equivalent width (right panel) of different models and star formation histories. Solid (dashed) lines show models without (with) nebular continuum emission. In the right panel we also show for comparison predictions for metal-free (Pop III) populations extending up to 500 M⊙ (blue dotted line). Dash-dotted lines show the “canonical” value of ξion used by Robertson et al. (2013) (left panel), and the mean relation between ξion and EW(Hα) derived by Tang et al. (2019) and Izotov et al. (2021) from observations.

From previous studies it is well know that the ionizing photon efficiency of an integrated stellar population not only depends on age, star formation history and the IMF, but also on metallicity (see e.g., Schaerer 2003; Raiter et al. 2010; Eldridge & Stanway 2022). In general, ξ_ion is maximum at young ages (≲2 Myr), increases for IMFs favoring massive stars (e.g., flatter slopes α₂ and or larger maximum stellar masses M_up), and increases toward low metallicities. The highest ionizing photon efficiencies reported from model predictions and “normal” IMFs (Salpeter slope and M_up ∼ 100 − 120 M_⊙) are log(ξ_ion) = 25.55 erg⁻¹ Hz for constant SFR and 26.05 at young ages for metal-free stellar populations (Raiter et al. 2010). For a flatter IMF (α₂ = −2.0), Eldridge & Stanway (2022) report a maximum of log(ξ_ion) = 25.55 erg⁻¹ Hz for constant SFR. For more extreme IMF including only massive stars, the ionizing efficiency can reach up to log(ξ_ion) = 26.1 erg⁻¹ Hz for metal-free stars, as shown by Raiter et al. (2010), and illustrated on the right panel of Fig. 6.

It should be noted that there is a natural upper limit to ξ_ion from stars/stellar populations, since the emission in the LyC and the non-ionizing UV are related (e.g., described by a Planck function or a more realistic SED) and since increasing LyC emission (Q_H) will also lead to increased nebular continuum emission (L₁₅₀₀ in particular). As Fig. 6 shows, the models including VMS up to 300–400 M_⊙ show quite similar ionizing efficiencies at young age, with values log(ξ_ion) ≈ 25.8 erg⁻¹ Hz when nebular continuum emission is accounted for, close to the maximum value of log(ξ_ion) ≈ 25.9 erg⁻¹ Hz predicted by the most extreme models (VMS-only). In short, the comparison of our models with these extreme values clearly show that stellar populations at moderately low metallicity can already reach very high ionizing photon efficiencies during a short burst phase. In fact, if stellar populations including only VMS – that is devoid of stars below ∼100 M_⊙ – exist, their efficiency will be log(ξ_ion) ≈ 25.8 − 25.9 erg⁻¹ Hz at LMC metallicity, within a factor of ∼2 of the maximum possible efficiency of PopIII stars. The upper limit on ξ_ion and other properties are further discussed in Sect. 4.

In the right panel of Fig. 6 we show ξ_ion as a function of the predicted Hα equivalent width, together with two average relations derived from observations of star-forming galaxies at low-redshift and at z ∼ 1.3 − 2.4 (Tang et al. 2019; Izotov et al. 2021). Although these relations are plotted here beyond the range over which they were established, one clearly sees that our models are shifted toward too high EW(Hα), even with the inclusion of nebular continuum emission. Most likely this shift is due to the presence of some previous/older stellar populations in the observed galaxies, which reduces the observed EW(Hα). A contribution of the order of ≳50% to the optical continuum is already sufficient to reconcile the two. Alternatively, or in addition, the observationally derived ξ_ion values could also be overestimated if the UV attenuation of these galaxies is systematically underestimated, although there is currently no indication for this (see e.g. Tang et al. 2019).

3.4. The UV luminosity of populations with VMS

One of the most direct and obvious implications of the presence of VMS is to “boost” the UV luminosity of stellar populations. For example, in the LMC cluster R136 the seven most massive stars – with masses between 110–325 M_⊙ – contribute 32% of the UV continuum flux at 1500 Å, as shown by Crowther et al. (2016). In the left panel of Fig. 7 we show the temporal evolution of the monochromatic luminosity L₁₅₀₀ for all the models from Table 1, both for instantaneous bursts and constant star formation. For comparison we also show predictions for a “top-heavy”, metal-free population, using the PopIII models from Raiter et al. (2010) with an IMF containing only stars between 50–500 M_⊙.

Fig. 7. Temporal evolution of the monochromatic UV luminosity L1500 (left) and the UV star formation rate conversion factor κUV (right panel, for SFR = const) predicted for models with a Salpeter IMF slope (α2 = −2.35) and a flatter IMF (α2 = −2.). L1500 is normalized to 106 M⊙ for bursts and to 1 M⊙ yr−1 for SFR = const. The models with a Salpeter slope are shown using the same colors as in previous figures. Models with the flat IMF, including also nebular emission), are shown by yellow and green lines. VMS-only models and PopIII models with a top-heavy IMF are shown by dotted lines. All models shown here include nebular continuum emission. Dash-dotted horizontal lines in the right panel show the value of the canonical SFR(UV) conversion factor from Madau & Dickinson (2014) for a Salpeter and a Kroupa IMF.

As Fig. 7 shows, for bursts the UV luminosity peaks around ∼2 Myr before it decreases, typically by one order of magnitude over ∼10 Myr. For constant SFR the UV luminosity progressively builds up with time, reaching an equilibrium after ≳50 − 100 Myr. For populations with “extreme” IMFs, containing only VMS, this equilibrium is reached already at ∼3 Myr, and further evolution is not shown here. At ages ≲3 Myr, an increase in the maximum stellar mass, M_up leads to a higher UV luminosity, with boosts up to a factor ∼2 when VMS up to 400 M_⊙ are included. At later times the boost factor of the UV luminosity is smaller, only of the order of 1.14 at ∼100 Myr for constant SFR. For bursts there is no difference after ∼3 Myr, the lifetime of VMS.

The slope of the IMF has a strong impact on the resulting UV luminosity of integrated stellar populations, as shown in Fig. 7. Indeed, a small flattening of the IMF with a change of slope from Salpeter (α₂ = −2.35) to α₂ = −2.0 already increases the UV luminosity significantly. At young ages, for example, models with VMS up to 400 M_⊙ and α₂ = −2.0 are brighter by a factor 5–6 than “classical” models (M_up = 100 M_⊙, α₂ = −2.35). For models with constant SFR the boost in UV luminosity is a factor ∼2.2 at ages ≳10⁸ yr. Models including only VMS show, quite naturally, even higher UV luminosities per unit stellar mass or SFR, as also shown Fig. 7. For constant SFR, L₁₅₀₀/SFR then surpasses the models with a fully populated IMF with α₂ = −2.0 even further, by an additional factor ∼2. Overall, between the UV luminosity predicted by standard BPASS models with M_up = 100 M_⊙ and the maximum for a “VMS-only” population, there is a factor ∼17 increase in non-ionizing UV output. It is also interesting to note that metal-free populations PopIII do not have a higher UV luminosity than a VMS-only population at Z = 0.006 at these wavelengths, since their average temperature is hotter and hence a larger fraction of the flux is emitted at higher energies (in the Lyman continuum). The maximum values predicted by the models shown here are summarized in Table 2.

The resulting conversion factor between SFR and UV luminosity, κ_UV, defined as

$$\begin{array}{c}\hfill \mathrm{SFR}(\mathrm{UV})={\kappa }_{\mathrm{UV}}\times L_{\nu ,1500},\end{array}$$(1)

where L_ν, 1500 is in conventional units of erg s⁻¹ Hz⁻¹, is shown in the right panel of Fig. 7. The canonical values of κ_UV = 1.15 × 10⁻²⁸ M_⊙ yr⁻¹/(erg s⁻¹ Hz⁻¹) taken from Madau & Dickinson (2014) for a Salpeter IMF from 0.1–100 M_⊙ and the one adjusted to the Kroupa IMF used here (κ_UV = 0.77 × 10⁻²⁸ in the same units, after scaling with a factor 0.67), are also shown in the figure for reference. As this figure shows, all of our models with a standard IMF slope (α₂ = −2.35) yield SFR conversion factors which are very close to the canonical values from the literature. To first order the effect of an extension of the IMF to VMS has therefore a small effect on the SFR(UV) calibration factor κ_UV. However, as discussed before, the IMF slope has a strong impact on κ_UV, which is reduced by a factor ∼2.4 to κ_UV = 0.32 × 10⁻²⁸ M_⊙ yr⁻¹/(erg s⁻¹ Hz⁻¹) for a slope of α₂ = −2.0. Furthermore, the timescale over which this equilibrium value is reached is also strongly reduced to ages of ≳20 − 30 Myr and the canonical value of κ_UV is already reached after ∼4 − 10 Myr, as shown in Fig. 7. Finally, models including only VMS (from 100–400 M_⊙) have the highest UV luminosity per unit SFR, yielding thus the lowest values found here with κ_UV ≈ 0.2 × 10⁻²⁸ M_⊙ yr⁻¹/(erg s⁻¹ Hz⁻¹).

As for ξ_ion, it is of interest to examine if there is a lower limit to the SFR(UV) conversion factor, or equivalently if there is an upper limit to the UV luminosity per unit SFR. The answer is indeed yes, and this limit is essentially given by the maximum UV luminosity L₁₅₀₀/M_⋆ per unit mass of massive stars. Now to first order L/M_⋆ is constant for stars above M_⋆ ≳ 100 − 150 M_⊙ (see e.g., Martins et al. 2020), and their effective temperature at the zero-age main sequence also, which implies that L₁₅₀₀/M_⋆ ≈ constant with stellar mass. Therefore, increasing M_up and/or flattening the IMF slope will not further increase the UV luminosity of a stellar population of fixed stellar mass or SFR beyond a certain limit. The VMS-only models shown here illustrate this limit for κ_UV and ξ_ion also, for moderately low metallicities.

We note that the value of κ_UV depends on metallicity but in a non-monotonic fashion, as has already been discussed by Raiter et al. (2010). Indeed, while the UV luminosity per unit stellar mass increases first from solar to low metallicities due to the increase of the average effective temperature of stars on the main sequence, L₁₅₀₀ decreases below ≲1/50 solar since the bulk of the flux of very hot stars is then emitted at wavelengths ≪1500 Å. For this reason the minimum value of κ_UV = 0.18 × 10⁻²⁸ M_⊙ yr⁻¹/(erg s⁻¹ Hz⁻¹) obtained here for Z = 0.006 for a very “top-heavy” IMF is even lower than for top-heavy IMFs at zero metallicity¹.

3.5. The UV slope

The predicted UV slopes of our models, β₁₅₅₀ measured here between 1300–1800 Å following Raiter et al. (2010), are plotted in Fig. 8, for models with and without nebular continuum emission, and for the two star formation histories. The effect of VMS can be appreciated at the youngest ages (≲3 Myr) and in models without nebular emission, where we see that the UV slope becomes somewhat redder (i.e., β increases) when VMS are included and with increasing M_up. This is due to the average effective temperature of the VMS, which is somewhat lower than that of stars with ∼100 M_⊙, as shown in Martins & Palacios (2022). In other words with increasing stellar mass the zero-age main sequence does not become hotter in the VMS regime. This behavior of the VMS leads to a non-monotonous evolution of β₁₅₅₀ at young ages (≲5 Myr). The resulting changes in β are of the order of Δβ ∼ 0.2 − 0.3, which is relatively small, and comparable to the uncertainties of many measurements of the UV slope at high-z.

Fig. 8. UV slope β1550 as a function of age for the main models presented in this work. Solid lines show the slopes derived for pure stellar models, neglecting nebular continuum emission. Dashed and dotted lines show β1550 including nebular continuum emission.

As is well known, the nebular continuum changes significantly the UV slope, making it redder by Δβ ∼ 0.3 − 0.6 at young ages (e.g., Schaerer 2003; Raiter et al. 2010). For models with constant SFR, the effect of VMS changes the UV slope by up to ∼0.1 after ∼10 Myr, and by smaller amounts over longer timescales. Again, the effect of nebular continuum emission, in other words including nebular emission or not, is more important than an extension of the IMF into the VMS regime.

Interestingly, for the reasons already mentioned above, the UV spectral slope of populations including or even dominated by VMS are not extremely blue, as one could naively have expected. For example, the stellar spectrum of populations hosting only VMS have −2.6 ≲ β₁₅₅₀ ∼ −2.0 (without nebular continuum). On the contrary, when nebular emission is added, VMS-only spectra become significantly redder, with β₁₅₅₀ ∼ −1.8. In practice, this means that UV slope should not be a good diagnostic for the presence/absence of VMS and other features should be examined for this purpose.

3.6. Direct signatures of VMS

Several observational features can potentially be used to detect the presence of VMS from integrated spectra (see e.g., discussions in Hadfield & Crowther 2006; Wofford et al. 2014; Smith et al. 2016, 2023; Martins et al. 2023; Crowther & Castro 2024). From our current knowledge of individual VMS studied in the Galaxy and the LMC, the clearest and strongest feature of VMS in the UV is their strong and broad He IIλ1640 emission line, which can significantly outshine the emission from Wolf-Rayet stars in this line (Martins et al. 2023; Upadhyaya et al. 2024). On the other hand, Crowther & Castro (2024) showed empirically that the contribution from classical Wolf-Rayet stars in mixed-aged populations (NGC 2070 including R136) can be significant (∼56%). There is, however, a concensus that existing stellar population models essentially all fail to reproduce the main features of observed UV-spectra of low-metallicity objects with strong He IIλ1640 emission without inclusion of VMS (see Senchyna et al. 2021; Smith et al. 2023; Martins et al. 2023; Crowther & Castro 2024).

The predicted strength of the He II line from our models is shown in Fig. 9, where the equivalent widths have been derived from the BPASS+VMS high resolution SEDs using the same spectral windows as in the study of Upadhyaya et al. (2024). Our results agree with theirs, but we present additional models accounting for nebular continuum emission and exploring other IMFs. A maximum EW(He II) ∼ 7 Å is predicted for models with a classical Salpeter IMF and M_up = 400 M_⊙ and the EW decreases with decreasing maximum stellar mass M_up. Steepening the IMF slope increases the maximum EW at young ages to ∼8 Å and leads to EW ∼ 4 − 5 Å for constant SFR after ≳10 − 100 Myr. Finally, the most “extreme” models considered here, where only VMS are present, predict He II equivalent widths between ∼8 − 14 Å depending on age and star formation history (bursts or constant SFR), as also shown in Fig. 9.

Fig. 9. Predicted equivalent width of the stellar He II emission as a function of age for the different models.

4. Discussion

4.1. Comparison with other studies

A detailed comparison between different models that allow for extensions of the IMF beyond ≳100 M_⊙ is beyond the scope of this paper. We only briefly mention several studies which discuss UV and ionizing properties of such populations and quickly comment on some aspects of interest.

The work of Raiter et al. (2010) presents extensive predictions for 8 different IMFs including in some cases VMS up to 500 M_⊙, covering a wide range of metallicities from PopIII (Z = 0) to solar metallicity. However, for metallicities Z ≥ 1/50 Z_⊙ the upper mass limit is 120 M_⊙. Stanway & Eldridge (2019) present BPASS models for a wide range of metallicities, M_up = 100 and 300 M_⊙, and two IMF slopes (identical as those used here). Their paper mostly focuses on the hardness of the ionizing spectra and discusses how the hardness depends on IMF variations. Our study differs from this work in that it uses tailored atmosphere models for VMS, which leads to softer spectra for VMS compared to those adopted in their models, for the reasons discussed in Sect. 3.2. Our study is not optimized to examine in detail the hardness of the ionizing spectra, since we neglect post main-sequence phases which can contribute to emission at higher energies, as included in the work of Stanway & Eldridge (2019). Overall, their results include those presented as our BPASSS-100 and BPASS-300 models.

BPASS predictions for different IMFs are also show in Eldridge & Stanway (2022). The ξ_ion values reported there (for the metallicity Z = 0.006 used here) are in good agreement with our values, when the same assumptions are made. Independent predictions of the ionizing photon production Q and κ_UV for different IMFs and metallicities have been presented in Lapi et al. (2024), who use these data to constrain the IMF during the epoch of reionization from comparisons to UV-luminosity functions and other observations. Overall their results appear to differ from ours. While their canonical value for Q (for the Chabrier IMF) agrees with our predictions, our top-heavy (VMS-only) IMFs predict significantly higher values of Q. Their maximum Q values are comparable to ours for α₂ = −2.35 and M_up = 400 M_⊙, which would indicate that despite the vast parameter space they explore, their model does not predict large enough variations in the ionizing photon production. Also, their standard value of κ_UV (which is defined as 1/κ_UV) seems incompatible with ours. Finally, their predictions for very top-heavy IMF (with characteristic masses m_c ≳ 100 M_⊙ seem incoherent, as they indicate lower UV and ionizing photons emission than for normal IMFs. In short, we find it difficult to understand their predictions.

4.2. Maximum UV and ionizing photon production of stellar populations and their implications

The maximum UV and ionizing photon emission from stellar populations is of interest to understand, e.g., the UV-bright high-z galaxies observed with the JWST, and in other contexts. The predictions from our models, exploring a range of IMFs, and including also for comparison the most extreme metal-free stellar population (IMF C from Raiter et al. 2010, a top-heavy IMF with stars between 50–500 M_⊙), are listed in Table 2.

Overall the table shows the following: First, for a young simple stellar population (burst) with a Salpeter IMF with different upper mass limits (up to 400 M_⊙), the ionizing photon production (Q) increases by up to a factor 2 and the rest-UV luminosity (here at 1500 Å) by up to a factor 1.6, whereas the maximum ionizing photon efficiency remains essentially constant, for reasons discussed earlier. Significantly stronger boosts are obtained with a flatter IMF, where a small change of the IMF slope (for example, changing α₂ = −2.35 to α₂ = −2.0) can lead to enhancements in the maximum ionizing photon emission by a factor 4–8 and a boost of the UV luminosity by a factor ∼4.8 (1.7 magnitude). Considering stellar populations exclusively made of VMS (M_⋆ ≳ 150 M_⊙), one obtains maximal values which are basically independent of the exact IMF (slope and upper mass limit), and which are ∼30 times higher for the ionizing photon flux and ∼21 times higher (3.3 magnitude brighter) for the UV luminosity, both per unit stellar mass. These are “absolute” maximal values (or a minimum UV mass-to-light-ratio) limited by fairly basic principles of stellar physics, which can therefore not significantly be surpassed by other stellar populations, irrespective of their metallicity and IMF. This also includes metal-free stellar populations (Pop III), as comparison with the model from Raiter et al. (2010) with the most extreme IMF shows. Only the maximum of ξ_ion and the hardness of the ionizing spectrum, measured for example by ratio of the ionizing photons emitted above 54 eV (the ionization potential of He⁺) compared to those above 13.6 eV, are higher for PopIII (for the reasons discussed earlier), exceeding the maximum values “VMS-only” populations.

Second, for constant star formation rates and assuming that one reaches an equilibrium after ≳10 − 100 Myr depending on the quantity considered, the ionizing photon production (Q) and ξ_ion increase by up to a factor 1.7, whereas the rest-UV luminosity is boosted by only ≲10 for a fully populated Salpeter IMF (α₂ = −2.35), extending down to the low-mass regime. Larger boosts, of a factor ∼2, are obtained with a flatter IMF (α₂ = −2.0). The maximum ionizing photon emission per unit SFR, obtained for populations consisting exclusively of VMS, is a factor ∼20 higher than with a classical IMF. For this case, the UV luminosity is boosted by a factor ∼4, implying a fairly fundamental value of min(κ_UV) ≈ 1.8 × 10⁻²⁹ M_⊙ yr⁻¹/(erg s⁻¹ Hz⁻¹) for the mimimum of the SFR(UV) conversion factor.

Finally, we also include for completeness the predictions neglecting nebular continuum emission in Table 2. In this case the UV emission at 1500 Å is decreased, the corresponding quantities are lower, and both ξ_ion and EW(Hβ) higher. In general and for the interpretation of integrated observations of star-forming galaxies in particular, nebular emission should be present, and hence models including it should preferentially be used. Values intermediate between the predictions with and without nebular emission are expected, for example if a fraction of the ionizing photons escape or are destroyed by dust inside the H II regions.

The “boost” in UV luminosity and ionizing production per unit stellar or star formation rate, has obviously potentially important implications for our understanding of star-forming regions and galaxies, and it could contribute to understanding or solving some of the puzzles mentioned in the introduction. As our models show, only small differences in the IMF slope have significant effects on the UV luminosity and ionizing photon production. And the presence of VMS further boosts these effects. Could the IMF be steeper and/or extend beyond the classical mass limit of ∼100 M_⊙? There is various evidence supporting both possibilities, at least in some environments/conditions. Individual stars with masses in the range of ∼100 − 300 M_⊙ have been found in the compact young star cluster R136 in the 30 Dor region (Crowther et al. 2010; Bestenlehner et al. 2020; Kalari et al. 2022). The presence of VMS is also suspected or established in several low-z star-forming regions (e.g., Leitherer et al. 2018; Smith et al. 2016, 2023; Senchyna et al. 2021; Wofford et al. 2023), in a lensed star cluster at z = 2.37 Meštrić et al. (2023), and in UV-bright galaxies at z ∼ 2 − 4 Upadhyaya et al. (2024).

The IMF slope of 30 Doradus, which includes R136, the prototypical nearby region with VMS, has extensively been studied. Surrounding the central region with R136, where crowding limited their analysis, Schneider et al. (2018) found indications for a flatter-than-Salpeter IMF, with a slope $\alpha =-1.{90}_{-0.26}^{+0.37}$ (in our units, where the Salpeter slope is α = −2.35), and an observed number of VMS in R136 compatible with this IMF. And many other studies indicate that the IMF could be flatter (at least for massive stars) in dense gas environments and/or at low metallicity (e.g., Dabringhausen et al. 2009; Marks et al. 2012; Dib 2023) Similarly at integrated galaxy scales, several studies have proposed a flatter or “top-heavy” IMF in galaxies with high star formation rates (e.g., Weidner et al. 2011) and in sub-mm galaxies (cf. Baugh et al. 2005; Zhang et al. 2018), although this result is disputed (see e.g., Safarzadeh et al. 2017). Other works have combined a variety of different constraints including supernova rates, the extragalactic background and others, concluding, however, that it is difficult to distinguish a Salpeter IMF from steeper ones (cf. Ziegler et al. 2022). Finally, further examples are works who examine possible shifts of the characteristic mass of the IMF with redshift, using detailed measurements of galaxy SEDs from large surveys. Sneppen et al. (2022) and Steinhardt et al. (2022), for example, suggest that the IMF shifts to somewhat higher characteristic masses with increasing redshift, becoming thus more top-heavy (or equivalently bottom-light), albeit with the same Salpeter slope for massive stars. Another recent approach to constrain the IMF from UV luminosity functions and observations probing the cosmic reionization history considers a more flexible description of the IMF, finding that IMF slopes flatter than Salpeter are excluded and favoring also a higher characteristic mass (see Lapi et al. 2024). Numerous other studies have also speculated about a non-universal IMF at high redshift (Trinca et al. 2024; Cueto et al. 2024; Woodrum et al. 2023), although no consensus exists on this question and clear observational tests are difficult to find.

4.3. Observational tests for extreme IMFs

4.3.1. Possible observational tests

Our spectral modeling of populations including VMS and considering non-standard IMFs can offer some hints on how to search for “extreme”, very top-heavy IMFs. The simplest signatures offering direct tests would be: (1) strong stellar He IIλ1640 emission, (2) very high ionizing photon efficiencies ξ_ion, (3) high L(Hα)/L(UV) ratios, (4) and/or H recombination lines with very high equivalent widths, as we shall now discuss and compare with recent JWST results.

The most top-heavy IMFs considered here, those including only VMS stars, can be tested quite easily and already excluded for many objects. Indeed, in such a case all galaxies should show very strong He IIλ1640 emission (with EW ≳ 8 Å), if the model predictions for the metallicity considered here are applicable, which we discuss below. Furthermore the ionizing photon efficiency and other quantities such as the Hα or Hβ EWs would be limited to very high values (log(ξ_ion) ≳ 25.7 erg⁻¹ Hz), since such populations would contain only short-lived and VMS. The fact that ξ_ion typically covers a relatively wide range including median values significantly lower than the above value excludes already such extreme IMFs for the majority of star-forming galaxies from of high to low-z (see e.g., Bouwens et al. 2016; Izotov et al. 2021; Prieto-Lyon et al. 2023, for measurements of ξ_ion from z ∼ 0 − 7).

Many other observations, including the color distribution of galaxies, the presence of Balmer break and others, testify to the presence of lower mass, longer or long-lived stars in most galaxies, which therefore exclude very bottom-light IMFs as a frequent or universal phenomenon. This is also the case for certain galaxies at high redshift, where Balmer breaks have been established with JWST both from photometry and spectroscopy (see e.g., Carnall et al. 2023, 2024; Looser et al. 2024). As is well known, the Balmer break, for example, probes A-type stars with temperatures T_eff ∼ 10 kK, which corresponds to main-sequence turn-off masses M_⋆ ∼ 2 − 3 M_⊙ (with lifetimes τ_H ∼ 300 − 700 Myr for metallicities (0.1 − 1) × solar; Ekström et al. 2012). If the presence of stars down to 2–3 M_⊙ is established it implies that ∼40% of the total mass of the population is probed, assuming a Kroupa/Chabrier IMF, which means that the UV luminosity per unit mass cannot be boosted by more than a factor ∼2.5.

In short, very bottom-light IMFs or extremely top-heavy IMF can be probed observationally, and IMFs with purely VMS are already excluded for bulk of high-z galaxies. Now, if we assume that the IMF is populated with stars over a wide mass range (e.g., M_low ≪ 2 − 3 M_⊙, and M_up ∼ 100 M_⊙ or higher) constraining the IMF, its exact slope and upper mass limit, is more challenging, as we shall now examine.

4.3.2. Comparisons with JWST observations

More specifically, we shall now discuss if there is some direct evidence for unusual IMFs and/or the presence of VMS in the high-z galaxies observed with JWST.

From stacked JWST spectra of galaxies from z ∼ 5 − 11, Roberts-Borsani et al. (2024) find that C III] λ1909 is the strongest line with EW ∼ 5 − 14 Å (increasing with z). This excludes already that these galaxies have on average a VMS-only IMF, since He IIλ1640 is weaker (or non-detected) than C III] λ1909, whereas VMS-only models predict EW(He II) ∼ 8 Å. On the other hand, the data may just be compatible with a flatter IMF and M_up ∼ 400 M_⊙ (our Flat-IMF+VMS-400), on average.

Then, individual galaxies should be better tests of VMS and extreme IMF, since the signatures of VMS (which have short lifeftimes) are best identified in objects dominated by young populations. Castellano et al. (2024) recently discovered a very bright galaxy at z = 12.34 showing numerous strong emission lines, including He IIλ1640. The possible detection of this line (albeit of low significance) with EW = 4.9 ± 3.1 Å is in the range of our models with VMS and a flat IMF. However, the line could also trace nebular emission (fully or partly) and there is thus no clear evidence for VMS in this object. Wang et al. (2024) reported the possible detection of He IIλ1640 in a blue galaxy at z = 8.16, suggesting the possible presence of Pop III stars co-existing with metal-enriched stellar populations. According to these authors He II is strong, with an equivalent width EW(He IIλ1640) = 21 ± 4 Å, and dominated by nebular emission, although the width of the line shown in their Fig. 3 appears relatively broad and comparable to observation by Upadhyaya et al. (2024) which indicate the presence of VMS. However, the strength of He II clearly surpasses those of previous observations in other objects and also model predictions for normal and extreme stellar populations. Further studies of this peculiar object would be warranted. Beyond these detections, most spectroscopic studies of high-z galaxies with JWST do not detect He IIλ1640 and provide only relatively loose upper limits. For example, Carniani et al. (2024) have observed two z ∼ 14 galaxies, detecting C III] emission but no He II, with upper limits EW(He IIλ1640) <  7 (16) Å (3σ) for the two objects. Similarly, for z ∼ 10 galaxies, Curtis-Lake et al. (2023) report upper limits of EW(He IIλ1640) <  13.5 − 15.4 Å for three galaxies (one has < 6 Å all 2σ limits). Clearly these data do not exclude the presence of VMS, flatter IMFs, and even the VMS-only scenario. Better, high S/N, medium resolution spectra are needed to provide tighter constraints.

Two interesting objects with clear nebular continuum emission and strong lines have recently been reported by Cameron et al. (2024). The authors suggest that these galaxies are “completely dominated” by nebular emission and would require a very top-heavy IMF. We find it unlikely that the galaxy GS-NDG-9422 at z = 5.943 which they discuss in detail, is dominated by nebular continuum emission, since the observed equivalent widths of the H recombination lines are significantly lower than the predicted values assuming a pure nebular continuum. Indeed, in this case one has EW = ϵ_ij/γ_tot, where ϵ_ij is the line emissivity and γ_tot the continuous emission coefficient at the wavelength of the line. For Hβ, for example, one has EW(Hβ)∼1000 Å for typical electron temperatures and assuming Case B, as shown in Fig. A.2. This value exceeds the observed equivalent width of GS-NDG-9422 (EW(Hβ)≈400 Å) by a factor ∼2.5, showing that the nebular continuum contributes approximately ∼40% of the total in the rest-optical. Also a boost of the two-photon continuum emission is accompanied by enhanced Hα/Hβ ratios (Raiter et al. 2010), opposite to the low Balmer decrement observed in this galaxy (Cameron et al. 2024). Furthermore, the observed ionizing photon efficiency, log(ξ_ion) ≈ 25.7 erg⁻¹ Hz as obtained from the Hβ flux and the observed UV continuum flux, falls short of the theoretical efficiency log(ξ_ion) ≈ 26.3 − 26.4 erg⁻¹ Hz for pure nebular continuum emission (see Appendix). Other sources of emission must therefore contribute, as has already been pointed out by Tacchella et al. (2024) and Li et al. (2024), who show that the spectrum of GS-NDG-9422 can be explained with a combination of nebular, stellar, and AGN emission, without the need for unusual stellar populations (see also Heintz et al. 2025).

Several studies have measured the ionizing photon production of high-z galaxies with JWST from photometry, by determining Q from the Hα luminosity, which can be obtained from “flux” excesses detected in specific filters, or in a more indirect fashion. For example, Prieto-Lyon et al. (2023) and Simmonds et al. (2024) find median values of log(ξ_ion) = 25.31 ± 0.43 and log(ξ_ion) >  25.5 erg⁻¹ Hz at z > 6, respectively, with some objects reaching log(ξ_ion) ≈ 26.0 erg⁻¹ Hz, close to the maximum predicted by our models. The median values are clearly higher than those predicted for standard IMFs at equilibrium; that is for constant SFR over sufficiently long timescales. However, this does not imply that IMFs favoring massive, or VMS are needed, if most of objects are dominated by populations younger than 10–50 Myr. According to our models log(ξ_ion) >  25.7 erg⁻¹ Hz requires either VMS (masses > 100 M_⊙) and or an IMF flatter than Salpeter. A fraction of the galaxies observed by Prieto-Lyon et al. (2023) and Simmonds et al. (2024) show indeed such high ionizing efficiencies, indicating possible IMF differences. Alternate explanations are an overestimate of ξ_ion if dust attenuation was underestimated, or lower metallicities on average. However, since the maximum of ξ_ion only increases by ∼0.1 dex when metallicity decreases from solar to 1/100 solar (cf. Raiter et al. 2010) this is insufficient to explain galaxies with log(ξ_ion) ≳ 25.8 erg⁻¹ Hz with a standard IMF. A more in-depth analysis of the high-ξ_ion objects could be of interest to understand if this is related to a peculiar IMF, or to additional but not accounted for ionizing sources (e.g. AGN). However, we note that the objects with the highest ξ_ion values are generally not UV-bright, which suggests that the latter are not strongly boosted in the UV and ionizing emission beyond a standard IMF.

Since by definition ξ_ion = Q_H/L_ν, 1500 = (L(Hα)/c_Hα) (κ_UV/L_ν, 1500), where L(Hα) is the Hα luminosity and c_Hα = ϵ₃₂/α_B (cf. Sect. A), the ionizing photon efficiency is also proportional to L(Hα)/L(UV) and the two quantities depend in the same way on the IMF slope and upper mass limits. Indeed, measurements of L(Hα)/L(UV) have long been used to measure and constrain the upper end of the IMF using observations of star clusters and integrated galaxy spectra (see e.g., Meurer et al. 2009; Lee et al. 2009; Bastian et al. 2010; Ashworth et al. 2017, and references therein). Alternatively, L(Hα)/L(UV) or analogous measurements of SFR(Hα)/SFR(UV) are often used to examine the “burstiness” or stochasticity of star formation (e.g., Domínguez et al. 2015; Emami et al. 2019; Atek et al. 2022), since this quantity also depends on the star formation history. Since L(Hα)/L(UV) and ξ_ion are equivalent, searches for IMFs favoring massive or VMS would therefore also look for objects exceeding the maximum L(Hα)/L(UV) predicted by standard models, as was already discussed above.

In principle, if the range between the maximum value of L(Hα)/L(UV) (or ξ_ion) at young ages and the value at equilibrium can be measured this can also be used to constrain the IMF, as can be seen from Table 2. Whereas for IMFs with a Salpeter slope these quantities decrease by 0.48 − 0.54 dex from maximum to equilibrium, this range is reduced for flatter IMFs (reaching e.g., 0.32 dex for Flat-IMF-400), down to ∼0.1 dex for VMS-only populations. The latter is already excluded for the majority of objects, which show log(ξ_ion) <  25.7 erg⁻¹ Hz from the JWST measurements. If stacking of multiple objects provides an average corresponding to constant SFR, measurements from stacked spectra could be compared to the maximum value (determined e.g., from individual objects) for this test. If this is feasible in practice remains, however, to be seen. Three recent studies have constructed stacked spectra for galaxies between z ∼ 3 and 10 from JWST NIRSpec observations (Roberts-Borsani et al. 2024; Langeroodi & Hjorth 2024; Kumari et al. 2024). For all galaxies Roberts-Borsani et al. (2024) measure log(ξ_ion) = 25.30 ± 0.01 erg⁻¹ Hz, in good agreement with the median log(ξ_ion) = 25.33 ± 0.47 erg⁻¹ obtained by Prieto-Lyon et al. (2023), and with our equilibrium predictions for constant SFR with a standard IMF extending to ∼200 − 300 M_⊙. From the data of Prieto-Lyon et al. (2023) the observed difference between their maximum log(ξ_ion) ≈ 26.0 erg⁻¹ Hz and this value is ∼0.7 dex, larger than the value predicted by the models for all IMFs considered here. Possibly this is due to uncertainties in the observational determination of ξ_ion, which are of the order of 0.1–0.2 dex or larger, and which can be overestimated if the dust attenuation is underestimated. Again, more in-depth studies and possibly additional/more data would be needed to constrain the IMF.

Finally, the last method to test for IMFs favoring massive stars or VMS is to search for objects with high equivalent widths of hydrogen lines (e.g., Hα, Hβ), which, if exceeding some maximum value for standard IMFs (say EW(Hβ)∼550 Å or even ≳800 Å if nebular continuum emission can somehow be avoided), would require non-standard IMFs. In principle this test is straightforward, although the amount of time during which EWs exceed this limit are short, less than ≲2 − 3 Myr at best, and therefore only very few objects are expected in such a phase. Other possibilities to encounter such high EWs include observations where stellar emission is spatially offset from the nebular emission region, in which case maximum EWs can reach very high values, as illustrated by the pure nebular case in Fig. A.2.

To summarize, several features accessible with rest-UV-to-optical observations can in principle be used to constrain the IMF or more specifically test if extreme IMFs showing deviations from the Salpeter slope and/or extending to significantly higher masses than commonly adopted. However, since these features also depend on other parameters, most importantly star formation histories, age and metallicity, some of which are difficult to constrain, these tests are not straightforward in practice. Whereas top-heavy IMF are therefore difficult to constrain, extremely bottom-light IMFs – for example populations consisting only of (very) massive stars which would be very UV-luminous per stellar mass – are probably fairly simple to exclude.

4.4. Caveats

Our predictions build on the BPASS synthesis models of Eldridge et al. (2017) including binary stars for a metallicity Z = 0.006, which are combined with the state-of-the-art stellar evolution and atmosphere models for VMS from (Martins & Palacios 2022) which have successfully been tested against observations of individual VMS in the LMC (R136). Our combined models have also been compared against other observations of low-redshifts star-forming regions and galaxies with metallicities comparable to the LMC (Martins et al. 2023). Finally, they have also successfully been applied to the model the Sunburst arc spectrum, a strongly lensed star cluster at z = 2.37, and they reproduce the observed rest-UV spectra of very UV-bright galaxies at z ∼ 2 − 4 that show signatures pointing toward the presence of VMS (Meštrić et al. 2023; Upadhyaya et al. 2024), although Wolf-Rayet features have been found in the Sunburst arc by Rivera-Thorsen et al. (2024).

Despite these successes and confrontations with observations our models have also some limitations. First, the VMS models used here follow only their main-sequence evolution, neglecting thus later evolutionary stages, primarily for computational reasons. In general post main-sequence evolution is typically ∼10% of the total stellar lifetime, which implies that the impact of these phases should overall be small or limited. The line intensity of He IIλ1640 depends mostly on luminosity, mass loss, helium mass fraction, and ionization (thus mostly T_eff). The most advanced models computed by Martins & Palacios (2022) already have more than 80% of mass in the form of helium. Their mass loss rates are high, of the order 10⁻⁴ M_⊙ yr⁻¹, which is the maximum observed for normal Wolf-Rayet stars in the LMC (Hainich et al. 2015; Aadland et al. 2022). The mass loss recipe of Gräfener (2021) do not predict significantly larger values for VMS. Finally, the luminosity of post main-sequence VMS in the LMC should be within 0.2 dex of that on the MS (e.g., Martinet et al. 2023). Consequently we do not expect a quantitative boost of the He IIλ1640 intensity in the latest phases of evolution, at least not strong enough to compensate for the reduced duration of these phases compared to the main-sequence. For comparison, the VMS models of Yusof et al. (2013) have a phase of very strong mass loss, which leads to the formation of WR stars after core H-burning accompanied by a rapid decrease in luminosity and an evolution to higher temperatures. We therefore primarily expect that such objects would increase the ionizing flux at high energies but not strongly alter the stellar He IIλ1640 emission which will be dominated by the more luminous VMS during H-burning. However, combined stellar evolution and atmosphere models including also these phases will be needed to properly quantify these statements. In any case, we note that from normal Wolf-Rayet stars we typically expect He II emission reaching EW(He IIλ1640) ≈ 2 − 3 Å during a short phase at ∼4 − 5 Myr, according to the models of Schaerer & Vacca (1998) at Z = 0.008, and that none of the existing stellar synthesis models with normal IMFs (i.e., without VMS) predict strong enough He IIλ1640 lines to explain the observations (see also Wofford et al. 2014; Smith et al. 2016, 2023; Martins et al. 2023; Crowther & Castro 2024).

Second, the present work discusses only models at a single metallicity Z = 0.006 ≈ 0.44 solar, since detailed, well-calibrated VMS models treating consistently both stellar evolution and atmosphere, are so far only available for this metallicity. This means that it is somewhat uncertain how our results can be transposed to other metallicities, in particular to lower metallicities of prime interest for high-z studies. Indeed, since several quantities, in particular the strength of the He IIλ1640 emission line, are very sensitive to the stellar wind properties of VMS, the main unknown is probably how mass loss depends on metallicity for VMS, a topic currently under investigation (Sander 2023; Sabhahit et al. 2023). Interestingly, from observations of Mrk 71, a star-burst region with ∼16% solar metallicity showing strong stellar He IIλ1640 emission, and comparisons with R136 and other higher-metallicity regions, Smith et al. (2023) find that the He IIλ1640 emission does not significantly change with metallicity and therefore suggest that the mass loss properties of VMS should be similar at low metallicity. First results from Martins et al. (in prep.) indeed suggest that this main VMS signature remains present and strong, at least down to 0.1 solar metallicity. Then, the metallicity-dependence of the properties of massive stars (i.e., UV luminosity, ionizing photon production etc.) is accounted for in standard synthesis models (BPASS and others), and, as discussed above, metallicity effects appear overall smaller than changes of the IMF, such as an extension toward VMS and/or changes of the IMF slope. We therefore expect that the relative effects of VMS on stellar populations are, to very first order at least, similar and independent of metallicity between ∼0.1 − 0.5 solar. Future work will be needed to examine this in detail and more quantitatively.

5. Conclusion

Combining the stellar evolution and atmosphere models of VMS from Martins & Palacios (2022) for a moderately-low metallicity (Z = 0.00536 = 1/2.5 Z_⊙) with BPASS evolutionary synthesis models from Eldridge et al. (2017) we have computed extensive sets of spectral energy distributions for integrated stellar populations hosting VMS with masses between 150 and 400 M_⊙. This work expands on the previous models and analysis of UV spectra of nearby star-forming regions, the Sunburst arc at z = 2.37, and UV-bright galaxies at z ∼ 2 − 4 presented earlier by Martins & Palacios (2022), Meštrić et al. (2023), Martins et al. (2023), and Upadhyaya et al. (2024).

We considered two limiting cases of star formation histories, instantaneous bursts and constant star formation and explored a wide range of IMFs, including different upper mass limits, and an IMF slope between Salpeter (α₂ = −2.35) and flatter (α₂ = −2). We also considered stellar populations with extreme IMFs, for example some consisting only of VMS stars with initial masses between 100–400 M_⊙. Accounting also for nebular continuum emission and H recombination lines, we predict integrated SEDs and derived quantities, such as the UV luminosity (L_UV), ionizing photon production (Q) and efficiency (ξ_ion), the strength of the Lyman break, H line luminosities and equivalent widths, the UV slope, and the strength of the He IIλ1640 line (which is dominated by VMS); the main quantities are presented in Table 2. The main predictions from our paper can be summarized as follows:

• VMS contribute significantly to the UV luminosity and Lyman continuum emission of young stellar populations: for example, for IMFs with a Salpeter slope the inclusions of stars up 400 M_⊙ boost the luminosity at 1500 Å, L_UV, by a factor of 1.63 at young ages and by ∼11% at equilibrium, for constant SFR over long timescales. Similarly, the ionizing photon production, Q, is boosted by factors of 2.45 (for young bursts) and 1.74 (constant SFR). Larger boosts are obtained with flatter IMFs, and only small changes in the IMF slope (adopting α₂ = −2 instead of the Salpeter slope) lead to increases by a factor of 4.8 (8.1) at young ages and 2 (5) for SFR = const for the non-ionizing and ionizing UV, L_UV (Q).

• The LyC emission of VMS leads to an increase in the ionizing photon efficiency, ξ_ion, by 0.14–0.2 dex for young bursts or constant SFR and a Salpeter slope. For α₂ = −2, the boost reaches up to 0.18–0.4 dex when the IMF is extended up to 400 M_⊙. Increased LyC emission also leads to stronger emission lines, with H recombination line luminosities boosted by the same factors as Q, and equivalent widths reaching EW(Hβ)∼700 Å at the youngest ages.

• The SEDs including VMS show reduced Lyman breaks. The shape of the ionizing spectra of populations with or without VMS is similar, between the LyC limit and up to ∼35 eV. At higher energies the spectra are softer with VMS, due to the increased opacity in the winds of VMS, caused by the boosted mass loss rate they experience compared to normal massive stars. The differences in ionizing spectra do not affect the classical, optical emission line (so-called BPT) diagrams.

• Populations with VMS show strong stellar He II emission, with EW(He IIλ1640) up to 4–8 Å at young ages or ∼2.5 − 4 Å for constant SFR, considering IMF slopes between α₂ = −2.35 and −2. (Fig. 9).

• Extreme IMFs populated only by stars in the VMS regime lead to “extreme” values for most of observational properties considered here. Per unit stellar mass (for young bursts) the UV luminosity (Q) can be boosted by a factor ∼20 (30) approximately. For constant SFR at equilibrium the maximum boost (compared to a classical Kroupa/Chabrier IMF extending to 100 M_⊙) is ∼4 (16) for L_UV (Q) by unit SFR. The ionizing photon efficiency of such populations is log(ξ_ion) ∼ 25.85 − 25.95 erg⁻¹ Hz, and even ∼0.1 dex higher if nebular continuum emission is absent. These values are essentially “fundamental” maximum limits, which are independent of the exact IMF slope and maximum stellar mass. This is the case since stellar models in the VMS domain follow L ∝ M and show similar zero-age main sequence properties. They are therefore to first order scale-free, and populations consisting purely of VMS are therefore indistinguishable, that is insensitive to the detailed shape and upper mass limit of the IMF, and only dependent on the total amount of stars formed as VMS. Compared to very top-heavy IMFs considered sometimes for metal-free populations, our moderately metal-poor VMS-only models emit ∼2 times more UV photons per stellar mass or SFR in the non-ionizing UV (1500 Å for example), since the average temperature is less extreme than those of Pop III, whose emission peaks well in the Lyman continuum. This comparison illustrates the relatively weak metallicity dependence of the predictions of such extreme populations (see also Raiter et al. 2010).

Regarding the need to explain UV-bright galaxies in the early Universe (see Sect. 1), we found that the addition of VMS in the range of 100–400 M_⊙ provides only relatively modest boosts of the UV luminosity, if the IMF slope remains Salpeter. However, small changes of the IMF slope can significantly increase the UV and ionizing flux output of galaxies.

Finally, we have presented several methods to observationally search for unusual IMF which would favor (very) massive stars. We showed how measurements of the He IIλ1640 line, ξ_ion, L(Hα)/L(UV), and H recombination lines can constrain the high mass tail of the IMF, and also how these quantities and others, such as the Balmer break, can be used to exclude “bottom-light” IMFs. Using published JWST observations we have examined if high redshift (z ≳ 5 − 6) galaxies show some evidence of the presence of VMS and/or signs of non-standard IMFs. Extreme IMFs containing only VMS can be excluded on average, but the current data does not exclude IMFs with Salpeter slope or somewhat steeper (e.g., α₂ = −2) including VMS for the majority of galaxies. Rest-UV spectra with higher spectral resolution and better S/N are needed to establish the existence of VMS in high-z galaxies. Furthermore, consistent stellar evolution and atmosphere models for VMS at lower metallicities than those currently available will be needed to examine if and how “universal” their spectral features are, and to quantify more broadly their impact.

It is our hope that the present study will help researchers to better understand the stellar content and IMF of the most distant galaxies and their impact on the ISM, cosmic reionization, and other related topics.

Data availability

The results from our models are made available in electronic format on Zenodo https://doi.org/10.5281/zenodo.14176049.

References

Appendix A: Nebular continuum and H recombination lines

We briefly illustrate the properties of the nebular continuum adopted in the present work and discuss several simple properties related to nebular continuum emission and H recombination lines. All predictions have been computed using pyneb (Morisset et al. 2020).

The monochromatic luminosity of the gas is given by

$$\begin{array}{c}\hfill L_{\lambda }^{\mathrm{neb}}=\frac{{\gamma }_{\mathrm{total}}}{{\alpha }_B}(1-f_{\mathrm{esc}})Q(\mathrm{H}),\end{array}$$(A.1)

where α_B is the case B recombination coefficient for hydrogen. The continuous emission coefficient γ_total, including free-free and free-bound emission by H and He, as well as the two-photon continuum of hydrogen is given by

$$\begin{array}{c}\hfill {\gamma }_{\mathrm{total}}={\gamma }_{\mathrm{HI}}+{\gamma }_{2\mathrm{q}}+{\gamma }_{\mathrm{HeI}}\frac{n({\mathrm{He}}^{+})}{n({\mathrm{H}}^{+})}+{\gamma }_{\mathrm{HeII}}\frac{n({\mathrm{He}}^{2+})}{n({\mathrm{H}}^{+})}.\end{array}$$(A.2)

The emission coefficients γ_i (in units of erg cm³ s⁻¹ Å⁻¹) are taken from pyneb. We assume a standard He abundance of 0.1 (in number) and that singly ionized He dominates (i.e., n(He⁺⁺)/n(H⁺) = 0.). The resulting continuum emission is illustrated in Fig. A.1 for different electron temperatures, T_e. The models discussed in this paper assume a typical value T_e = 15000 K. Higher (lower) electron temperatures would result in somewhat lower (higher) values for the predicted ξ_ion and equivalent widths of Hα and Hβ. The effect of doubly ionized He, if present, should be negligible in most cases.

Fig. A.1. Nebular continuum emission γtotal/αB of H plus He (with standard abundances) as a function of wavelength and for different electron temperatures. The default value adopted here is shown in black. The dotted line shows γtotal/αB for a highly ionized mixture with n(He+)/n(H+) = n(He2+)/n(H+) = 0.05.

At metallicities ≲0.03 − 0.1 solar, well below those adopted here, the nebular continuum can in principle be enhanced compared to the predictions shown here. Such a boost is due to collisional effects, which overpopulate the n = 2 level and enhance the photoionization rate, boosting thus the two-photon continuum emission by a factor which is proportional to the mean energy of the ionizing photons. At the same time the collisions also modify the Balmer decrement, increasing, for example, the Hα/Hβ ratio (see Raiter et al. 2010). At the metallicity adopted here these effect are, however, negligible.

Significant increases of the electron density beyond n_e ≳ 10^3 − 4 cm⁻³ decrease the two-photon emission, in other words nebular emission in the rest-UV, but has little effect at longer wavelengths. On the other hand, this collisional effect leads to a boost of Lyα emission. (For a more detailed discussion of these effects, which may affect the nebular continuum and recombination line properties of individual objects at very low metallicity see e.g., Bottorff et al. 2006; Raiter et al. 2010; Katz et al. 2024).

By definition, the ionizing photon efficiency is

$$\begin{array}{c}\hfill {\xi }_{\mathrm{ion}}=\frac{Q_H}{L_{1500}}=\frac{Q_H}{L_{1500}^{\star }+L_{1500}^{\mathrm{neb}}}\end{array}$$(A.3)

Hence, for a case where the stellar continuum can be neglected, one has

$$\begin{array}{c}\hfill {\xi }_{\mathrm{ion}}^{\mathrm{neb}}=\frac{{\alpha }_B}{{\gamma }_{\mathrm{total}}(1-f_{\mathrm{esc}})}.\end{array}$$(A.4)

For electron temperatures between 10-20 kK one obtains log(${\xi }_{\mathrm{ion}}^{\mathrm{neb}})=26.3-26.4$ erg Hz⁻¹.

For completeness, the fraction of energy contributed by nebular continuum emission at 1500 Å, f_neb, can easily be inferred from

$$\begin{array}{c}\hfill f_{\mathrm{neb}}=1-\frac{{\xi }_{\mathrm{ion}}}{{\xi }_{\mathrm{ion}}^{\star }},\end{array}$$(A.5)

where ξ_ion^⋆ = Q_H/L₁₅₀₀^⋆ is the ionizing photon efficiency neglecting nebular continuum emission, whose values are reported at the bottom of Table 2.

Similarly, for H recombination lines the luminosity is

$$\begin{array}{c}\hfill L^{\mathrm{ij}}=\frac{{ϵ}_{\mathit{ij}}}{{\alpha }_B}(1-f_{\mathrm{esc}})Q(\mathrm{H}),\end{array}$$(A.6)

where ϵ_ij is the line emissivity. Hence the equivalent width, neglecting the stellar continuum becomes

$$\begin{array}{c}\hfill E{W}^{\mathrm{neb}}(\mathrm{ij})=\frac{{ϵ}_{\mathit{ij}}}{{\gamma }_{\mathrm{total}}}.\end{array}$$(A.7)

These values and their dependence on electron temperature are shown in Fig. A.2 for selected Balmer lines. For Hβ, one has typically EW(Hβ)≈1000 Å. These purely nebular EWs are, for example, convenient to infer the contribution of the stellar (and/or other) continuum sources from the measured EW of a H recombination line.

Fig. A.2. Predicted equivalent width of the equivalent widths of Hα, Hβ, and Hγ for a nebular continuum as function of the electron temperature

All Tables

All Figures

Fig. 1. Comparison of predicted SEDs for a 1 Myr old population with and without VMS. Left panel: The standard BPASS model with Mup = 100 M⊙ is shown in gray (stellar emission only) and in black (including nebular continuum emission). Blue and red lines show the model extending up to 400 M⊙, the former including only stellar emission, the latter the total continuum (stellar plus nebular continuum). Right panel: zoom on the UV domain, comparison the standard BPASS model with a population including only VMS (from 100 to 400 M⊙). The dotted line shows the nebular continuum emission of the VMS-only-400 population.

In the text

	Fig. 2. Predicted temporal evolution of the ionizing photon production, QH, of different models, for instantaneous bursts and constant SFR. Burst models are normalized to 106 M⊙, models for constant SFR to SFR = 1 M⊙ yr−1.
In the text

	Fig. 3. Predicted temporal evolution of the Hβ equivalent width of different models, for instantaneous bursts and constant SFR. Solid (dashed) lines show models without (with) nebular continuum emission.
In the text

	Fig. 4. Predicted temporal evolution of the L900/L1500 ratio of different models, for instantaneous bursts and constant SFR. Solid (dashed) lines show models without (with) nebular continuum emission.
In the text

Fig. 5. Comparison of the ionizing spectra of different models for age 1 Myr. Shown is the cumulative number of ionizing photons emitted above energy E, Q(> E), as a function of the energy (in eV). All ionizing spectra are normalized to unity at the Lyman limit (13.6 eV), to allow for easy comparison of their relative hardness. Vertical lines mark the ionization potentials of several important ions.

In the text

Fig. 6. Ionizing photon efficiency ξion as a function of age (left panel) and Hα equivalent width (right panel) of different models and star formation histories. Solid (dashed) lines show models without (with) nebular continuum emission. In the right panel we also show for comparison predictions for metal-free (Pop III) populations extending up to 500 M⊙ (blue dotted line). Dash-dotted lines show the “canonical” value of ξion used by Robertson et al. (2013) (left panel), and the mean relation between ξion and EW(Hα) derived by Tang et al. (2019) and Izotov et al. (2021) from observations.

In the text

Fig. 7. Temporal evolution of the monochromatic UV luminosity L1500 (left) and the UV star formation rate conversion factor κUV (right panel, for SFR = const) predicted for models with a Salpeter IMF slope (α2 = −2.35) and a flatter IMF (α2 = −2.). L1500 is normalized to 106 M⊙ for bursts and to 1 M⊙ yr−1 for SFR = const. The models with a Salpeter slope are shown using the same colors as in previous figures. Models with the flat IMF, including also nebular emission), are shown by yellow and green lines. VMS-only models and PopIII models with a top-heavy IMF are shown by dotted lines. All models shown here include nebular continuum emission. Dash-dotted horizontal lines in the right panel show the value of the canonical SFR(UV) conversion factor from Madau & Dickinson (2014) for a Salpeter and a Kroupa IMF.

In the text

	Fig. 8. UV slope β1550 as a function of age for the main models presented in this work. Solid lines show the slopes derived for pure stellar models, neglecting nebular continuum emission. Dashed and dotted lines show β1550 including nebular continuum emission.
In the text

	Fig. 9. Predicted equivalent width of the stellar He II emission as a function of age for the different models.
In the text

	Fig. A.1. Nebular continuum emission γtotal/αB of H plus He (with standard abundances) as a function of wavelength and for different electron temperatures. The default value adopted here is shown in black. The dotted line shows γtotal/αB for a highly ionized mixture with n(He+)/n(H+) = n(He2+)/n(H+) = 0.05.
In the text

	Fig. A.2. Predicted equivalent width of the equivalent widths of Hα, Hβ, and Hγ for a nebular continuum as function of the electron temperature
In the text