Ultraviolet spectra of extreme nearby star-forming regions: evidence for an overabundance of very massive stars

By virtue of their selection from the SDSS spectroscopic survey, high-quality imaging and optical spectra are already available for all seven new targets. The broadband photometric magnitudes of the target regions provide important constraints on the total stellar mass and SFR enclosed in the spectroscopic apertures. We measure aperture photometry centered on each target region from SDSS sky-subtracted and calibrated $ugriz$ frames covering each target, adopting a 3^′′-diameter aperture corresponding to the size of the SDSS spectroscopic aperture.

The public SDSS spectra for these systems provides $R\simeq 1500$–$2500$ coverage over 3800–9200 Å. As our targets reside at very low-redshifts $z<0.03$, the [O ii] $\lambda\lambda 3727,3729$ doublet (hereafter [O ii] $\lambda 3727$ as we observe it as an unresolved single line) remains outside the SDSS spectroscopic range for all but the highest-redshift, SB 61.²²2While [O ii] $\lambda 3727$ is indeed barely shifted into the wavelength range of the SDSS spectrum for the second highest-redshift system SB 125 at $z=0.018$, the bluest part of the line profile is cut-off at the edge of the SDSS spectrum. The 3727 doublet is an important component of gas-phase oxygen abundance and ionization parameter determination. Therefore, in addition to the SDSS spectra, we obtained supplementary blue spectra with the Blue Channel spectrograph mounted on the 6.5m MMT telescope on Mount Hopkins. Spectra for four of the targets (SB 119, 125, 126, and 153) were obtained on the night of April 9, 2019 with the 800gpm grating and $1.5^{\prime\prime}\times 180^{\prime\prime}$ slit, yielding spectra spanning 3400–5200 Å. Seeing as-measured from the wavefront sensor ranged from 1.1^′′–1.5^′′ over the course of the night. In addition, SB 49 was observed on the night of October 8, 2019 with the 300gpm grating and $1.0^{\prime\prime}\times 180^{\prime\prime}$ slit, resulting in coverage of 3500–8500 Å. All observations were reduced using standard longslit data reduction techniques, including wavelength calibration using HeAr/Ne lamp exposures and flux calibration with standard star observations conducted at similar airmass (Feige 34 and BD+33 2642 on April 9, and HZ 14 on October 8). To correct for potential differences in aperture and flux calibration between the SDSS and MMT spectra, we correct the flux of [O ii] $\lambda 3727$ measured in the MMT data by the median ratio between other strong lines measured in both spectra. We find a median scatter in this ratio measured between different lines of 8%, which we add in-quadrature to the corrected flux uncertainty.

Nebular line flux and equivalent width measurements were performed using the technique described in more detail in S17. In brief, a base model consisting of a Gaussian emission line (or two in the case of doublets or fitting for multiple components to a single line, as specified) and a linear continuum are fitted to the spectroscopic data in the region of each emission line of interest using the Markov chain Monte Carlo sampling framework emcee (Foreman-Mackey et al., 2013). The fits are examined by-eye to ensure that the fit wavelength range is sufficient to accurately capture the continuum level and ensure that the correct features are identified. The resulting posterior distributions after removing a nominal burn-in period straightforwardly provide confidence intervals for the parameters of interest, including fluxes, equivalent widths, and line widths.

With this photometry and nebular line information in-hand, important constraints can be placed on the bulk star formation history and ionized gas properties of these star forming regions. As a first step towards constraining the timescale of recent star formation in these systems, we fit the broadband photometric measurements and H$\beta$ equivalent widths measured from the SDSS data. In particular, we use the BayEsian Analysis of GaLaxy sEds (beagle, v0.23.0; Chevallard & Charlot, 2016) tool to model each system. The underlying stellar models and population synthesis methodology will be discussed in more detail in Section 5 and have been previously described by Gutkin et al. (2016). We assume a constant star formation history (CSFH) with an SMC extinction curve (Gordon et al., 2003) and allow the age of this star formation period ($0<\log_{10}[t/\mathrm{yr}]<9$), metallicity ($-2.2<\log_{10}[Z/Z_{\odot}]<0$), stellar mass ($0<\log_{10}[M/M_{\odot}]<10$), nebular ionization parameter ($-4<\log U<-1$) and $V$-band attenuation ($0<\hat{\tau}_{V}<2$) to vary over the indicated uniform prior ranges. We inflate the photometric flux uncertainties by a conservative 5% added in-quadrature to both more heavily-weight the spectroscopic H$\beta$ equivalent width and account for possible differences in spectroscopic versus photometric aperture, and find good agreement with all fit fluxes in the posterior PDFs.

The results of these SED fits are displayed in Table 3. By selection, these systems display extremely high equivalent width optical nebular line emission. Our measurements of their SDSS H$\beta$ equivalent widths range from 125 to 285 Å, yielding correspondingly strong broadband nebular line contamination. In particular, the target $g$-band photometric measurements (contaminated by H$\beta$ and [O iii] $\lambda\lambda 4959,5007$) are brighter than measurements in the adjacent $u$-band (close to the continuum, contaminated only by weaker [O ii] $\lambda 3727$) by 0.3–1.0 (median 0.5) magnitudes. The results of our SED fitting reflect this, indicating that the photometry of all seven objects is consistent with entirely young stellar populations. The median inferred age of our assumed constant star formation history model range from a very young $3.7$ Myr in the most extreme case (SB 153) to $\sim 20$Myr for systems with lower equivalent-width optical emission (SB 125, 126), consistent with the presence of WR wind signatures in the optical spectra (Figure 1). These ages represent an upper limit for the young stellar component of the target systems, as adopting a single-age simple stellar population would require ages uniformly $<10$ Myr to reproduce the observed optical line equivalent widths. The stellar masses inferred from this fitting are uniformly low, spanning $10^{4.9}$–$10^{7.0}$ $M_{\odot}$, as expected for very young stellar populations with low mass-to-light ratios. This fitting strongly supports the idea that the spectroscopic apertures centered on these systems are entirely dominated by continuum and nebular gas emission from massive stars formed in a very recent star formation episode.

Measurements of the prominent collisionally-excited nebular lines provides direct constraints on the ionization state and chemical abundances of this ionized gas. As in S17, we use the software pyneb (Luridiana et al., 2015) to derive gas-phase oxygen abundances. We use the temperature-sensitive ratios of [O ii] $\lambda 7325$ / $\lambda 3727$ and [O iii] $\lambda 4363$ / $\lambda 5007$ alongside the density-sensitive [S ii] $\lambda 6731$ / $\lambda 6716$ ratio to derive electron temperatures and densities appropriate for $\mathrm{O^{+}}$ and $\mathrm{O^{2+}}$, respectively. ³³3We assume atomic transition probabilities for O ii from Wiese et al. (1996); Froese Fischer & Tachiev (2004) and effective collision strengths from Tayal (2007). For O iii we adopt the atomic data of Storey & Zeippen (2000); Froese Fischer & Tachiev (2004) and collision strengths of Storey et al. (2014). And for S ii, we utilize the transition probabilities presented by Podobedova et al. (2009); Tayal & Zatsarinny (2010) and collision strengths from Tayal & Zatsarinny (2010). In the one case for which we lack a measurement of [O ii] $\lambda 3727$ (SB 9), we instead estimate $\mathrm{T_{e}}$(O ii) from that measured for O iii using the empirical relationship derived by Izotov et al. (2006, equation 14, for the $12+\log\mathrm{O/H}>8$ subsample). Attenuation corrections are derived from the Balmer decrement assuming an SMC internal extinction curve (Gordon et al., 2003) and an intrinsic value of H$\alpha$/H$\beta$ computed with pyneb using our derived $\mathrm{T_{e},n_{e}}$ after correcting first for Galactic extinction towards each object using the Schlafly & Finkbeiner (2011) maps and the $R_{V}=3.1$ extinction curve of Fitzpatrick (1999).

The resulting metallicities and other optical nebular line properties are presented in Table 4. In addition to high equivalent-width emission indicative of an extremely young effective age, these systems exhibit signs of very highly-ionized gas, with dust-corrected $\mathrm{O}_{32}$$=$ [O iii] $\lambda 4959$ + $\lambda 5007$ / [O ii] $\lambda 3727$ values of 5.3–10.7. Their gas-phase metallicities span the regime $8.0<12+\log_{10}\mathrm{O/H}<8.3$, which assuming a solar abundance of $12+\log(\mathrm{O/H})_{\odot}=8.69$ corresponds to a relative metallicity range of $0.2<Z/Z_{\odot}<0.4$.

We secured moderate-resolution UV spectra with HST/COS for all seven target regions in a program approved in HST Cycle 25 (GO:15185, PI:Stark). Our observations were conducted in the NUV with the G185M grating and in the FUV with G160M, with a 2-orbit visit devoted to each object. Each target was first centered with a 7–81 second ACQ/IMAGE exposure with either Mirror A or B, with the optical element and exposure time chosen with the COS ETC based upon an assumed compact source with flat SED in $F_{\nu}$ normalized to the SDSS $u$-band fiber magnitude. The resulting target acquisition images are shown in Fig. 2, and reveal a range of morphologies at 10–100 pc scales (though note that the marked Mirror B images are complicated by the fainter, displaced secondary image this optical element generates). The remainder of the first orbit of each visit was devoted to G185M observation, with the second orbit occupied by G160M. The central wavelength of each observation was selected with the SDSS redshift in mind to ensure that the full stellar C iv $\lambda\lambda 1548,1550$ profile, He ii $\lambda 1640$ and O iii] $\lambda\lambda 1661,1666$ were located in G160M segments A and B, and that the [C iii] $\lambda 1907$ + C iii] $\lambda 1909$ doublet (hereafter C iii] $\lambda 1908$) was centered in the middle segment NUVB.

Two visits of our program required repeat observation. On 29 April 2018, guide star tracking was lost during reacquisition at the beginning of the second orbit for SB 126. While the G185M integration was successfully completed before fine guidance loss occurred, this prevented the G160M integration (LDI203020) from proceeding. A single-orbit repeat was scheduled as-per HOPR 90495, and a G160M spectrum was successfully recorded for this object on 27 May 2018 (LDI208010). During the original scheduled visit targeting SB 61 on 20 June 2018, the FGS acquired the guide stars after target acquisition had been attempted. As a result, the aperture was not correctly centered on-target for the spectroscopic exposures, resulting in an unexpectedly low continuum level and signal-to-noise in both the G160M and G185M spectra (LDI207020 and LDI207010, respectively). These observations were repeated following HOPR 90852 with a successful target acquisition on 9 August 2018, yielding spectra that reach the expected S/N which we use in this paper (LDI209020 and LDI209010).

All spectroscopic observations were taken in time-tagged photon-address mode (TIME-TAG) with wavelength calibration lamp exposures (FLASH=YES) to allow for optimal data correction post-observation. We cycle through all four focal plane offset positions for each grating (FP-POS=ALL) to mitigate the effect of fixed pattern noise. The data were reduced with the default settings for calibration and extraction using CALCOS v3.3.5 and CRDS v7.3.0 (HSTDP v2019.2). The final one-dimensional spectra, with resolution elements (resels) of 73.4 mÅ in the FUV/G160M and 102 mÅ in the NUV/G185M, achieve a typical effective spectral resolution of 0.25 Å as-inferred from the width of Milky Way absorption lines (identical to that found in S17, ) and median signal-to-noise of 11 (4) per resel at 1450 Å (1700 Å) in G160M and 1.7 per resel at 1900 Å in G185M.

We measure emission lines in the HST/COS spectra in the same manner as for the optical emission lines. The high spectral resolution afforded by the G160M and G185M gratings allows us to easily identify and separate Milky Way (MW) and ISM absorption lines, alleviating blending concerns except in the most extreme cases. In addition to our measurement of O iii] and C iii] semi-forbidden nebular line emission, we define two integration regions designed to measure the strength of the two key stellar wind features accessed by the G160M spectra. In order to ensure that ISM and MW absorption lines do not contaminate these measurements, we first mask-out the spectra $\pm 1.5$ Å around all identified strong absorption lines and interpolate the spectrum over these windows. We estimate the equivalent width of the C iv P-Cygni wind absorption trough by integrating this cleaned spectrum between 1530 and 1550 Å relative to a linear continuum determined by evaluating the median flux in two windows chosen to avoid wind and nebular line features: [1500:1525] and [1565:1575]. Likewise, we measure the strength of emission in the He ii $\lambda 1640$ wind and nebular emission line in the window 1630–1650 Å relative to a continuum estimated from flux at [1620:1630] and [1650:1660] ⁴⁴4We use only the blueward continuum window to integrated He ii $\lambda 1640$ for SB 80 from S17 as strong Al ii $\lambda 1670$ absorption contaminates a significant portion of the redward window.. We present the resulting measurements in Table 5.

The measured properties of the UV spectra of our sample (Figure 3) further evince a dominant population of massive stars in these galaxies. We detect C iii] emission in all seven systems, at equivalent widths ranging from 2–9 Å and O iii] $\lambda 1666$ in all but one (SB 126, where it is fully absorbed by the MW Al ii $\lambda 1671$ line). Including the three additional S17 objects reveals a similar detection rate, extending the C iii] equivalent width range up to just over 11 Å (in SB 191) with only one system undetected in this doublet (SB 80). The essentially unanimous detection of this doubly-ionized carbon and oxygen emission in the UV confirms the presence of a strong ionizing continuum at $\sim$20–40 eV.

As mentioned, the other high-ionization line complexes in our G160M spectra provide a direct window onto the massive stars powering this nebular gas emission. The C iv $\lambda\lambda 1548,1550$ doublet complexes show no clear sign of narrow nebular emission, instead revealing a characteristic P-Cygni profile consisting of broad redshifted emission and blueshifted absorption extending down to $\sim 1530$–1535 Å ($\sim 2500$–3500 km/s), reaching equivalent widths in absorption of 2.3–6.5 Å. This feature is formed in the dense winds of massive O stars, and is observed to extend to comparable terminal velocities in LMC stars at $Z/Z_{\odot}\simeq 0.5$ (e.g. Massey et al., 2004; Crowther et al., 2016). Likewise, rather than the nebular gas emission typically observed in He ii $\lambda 1640$ at metallicities below $\sim 0.2Z_{\odot}$, the target systems all show broad emission with FWHM $\sim 1500$–2000 km/s. This emission is characteristic of massive WN and O If stars, and is formed in extremely dense and optically thick stellar winds driven by stars approaching the Eddington limit. The prominence of this emission is particularly extraordinary, with three systems exceeding 4 Å. We will discuss the context and implications of these detections in more detail in the following section.

We will conduct this comparison experiment against the updated C&B single-star models using the beagle inference framework. As implied by our findings in Section 4, the improvements in the treatment of very massive star evolution adopted by these models has alleviated some of the tension previously found with the stellar He ii wind line; but these models still neglect the physics of binary evolution and rotation, which may have a significant impact even at very young stellar population ages. The version of the C&B models used in this work is described in more detail in Vidal-García et al. (2017, their Section 2.1 and Appendix A, updating ). The most important features of this code for the purposes of this work concern the revised predictions for the evolution and atmospheres of the most massive stars (extended now to 300 $M_{\odot}$). The latest predictions from the PAdova and TRieste Stellar Evolution Code (parsec) are now adopted to predict the evolution of individual massive stars (Bressan et al., 2012; Chen et al., 2015). Notably, these latest parsec results include a revised formalism for the prediction of stellar wind driven mass loss rates based upon the stellar Eddington factor (Vink et al., 2011). This is in contrast to most other stellar evolutionary codes which predict O and WR mass loss rates using older calibrations (Vink et al., 2001) derived under assumptions about wind driving which likely break down for stars at the highest luminosities and sub-solar metallicities (see e.g. Gräfener & Hamann, 2008; Müller & Vink, 2008; Sander et al., 2020). In addition, the C&B models leverage the large library of non-LTE, line-blanketed, spherically-expanding WR model atmospheres produced by the Potsdam Wolf Rayet group (PoWR) spanning metallicities from $Z=0.001$–$0.014$, effective temperatures from $10^{4.4}$–$10^{5.3}$ K, and a range of CNO mass fractions⁶⁶6see http://www.astro.physik.uni-potsdam.de/~wrh/PoWR/powrgrid1.php. (Gräfener et al., 2002; Hamann & Gräfener, 2003; Hainich et al., 2014; Hainich et al., 2015; Sander et al., 2015; Todt et al., 2015). These updated mass loss and atmosphere prescriptions for massive stars near the Eddington limit have significantly improved bulk agreement with observations for stellar lines like He ii $\lambda 1640$ produced in the optically-thick winds such Of and WR stars drive, as we discussed in Section 4 above.

The spectral inference framework beagle (Chevallard & Charlot, 2016, introduced in Section 3.1) is readily able to fit spectrophotometric observables using the C&B stellar population models and associated cloudy nebular line predictions, and provides the inference framework we require for this experiment. For this analysis, we choose to adopt the fiducial Chabrier (2003) initial mass function (IMF) with upper mass cutoff of 300 $M_{\odot}$ (as we will discuss below, the contribution of very massive stars $>100M_{\odot}$ to high-ionization stellar wind lines can be dominant; e.g. Crowther et al., 2016). For all fundamental variables describing the stellar population, we adopt a uniform prior over the adopted range for simplicity. Since the initial photometric fits including nebular emission were well-fit by a CSFH model (Section 3.1) and this likewise populates the region of stellar wind parameter space occupied by our observations (Figure 4), we parameterize our assumed SFH as such. However, in addition to the variable start time (which we allow to vary from $10^{6.5}$–$10^{9.0}$ years) and stellar mass ($M\leq 10^{10}M_{\odot}$) for this recent period of star formation, we allow an early truncation of up to 10 Myr before the present time to allow additional flexibility in the mix of very massive stars present. We allow the stellar metallicity to freely vary with $-2.2\leq\log_{10}(Z/Z_{\odot})\leq 0$. While this framework explicitly couples the interstellar metallicity to this stellar metallicity, depletion of refractory elements onto dust grains can change the gas-phase oxygen abundance at fixed $Z$. We allow the dust-to-metal mass ratio $\xi_{d}$ to vary from 0.1–0.5, which effectively corresponds to a factor of 1.6 range in the ratio of gas-phase oxygen to stellar iron (Gutkin et al., 2016). This ratio plays an important role in modulating the joint nebular and stellar spectra of star-forming systems (e.g. Steidel et al., 2016; Strom et al., 2018; Sanders et al., 2020), since oxygen is the primary species probed in the nebular spectrum while the winds and ionizing spectra of massive stars are shaped primarily by iron. Thus a model preference for low $\xi_{d}$ can be interpreted in part as a preference for enhanced $\alpha$/Fe in this context. Since we are studying relatively small star-forming regions, we assume an SMC extinction curve with variable $0\leq\hbox{$\hat{\tau}_{V}$}{}\leq 2$. Finally, we allow the nebular ionization parameter to vary over $-4\leq\log U\leq-1$.

In this experiment, we are interested first in fitting the nebular line fluxes which constrain the stellar ionizing spectrum and gas-phase metallicity. In particular, we focus on the set of strong optical nebular lines necessary to measure the gas-phase oxygen abundance via the direct method: [O ii] $\lambda 3727$, [O iii] $\lambda 4363$, H$\beta$, [O iii] $\lambda\lambda 5007$, and H$\alpha$ Since the gas densities inferred from the [S ii] $\lambda\lambda 6716,6731$ doublet are all consistent with 100–200 $\,\mathrm{cm}^{-3}$, for simplicity we exclude these lines from the fits and instead fix the gas density to the nearest grid value of 100 $\,\mathrm{cm}^{-3}$. For additional leverage on the dust attenuation, we add H$\gamma$; and to better anchor the recent star formation history of the models, we also fit equivalent widths measured in SDSS for H$\gamma$, H$\beta$, [O iii] $\lambda 5007$, and H$\alpha$. In each case, the observed line flux or equivalent width is compared to the flux measured via simple integration using a linear continuum in the model spectra. We choose the integration region and continuum windows for each line manually to ensure that they capture the entire line flux and that the continuum is not biased by other features. Since our objective is to determine whether the UV stellar wind lines can be reproduced, we also perform a second fit for each object including the equivalent widths of these features measured as for the data (as in Section 3.2, though obviously without the need for ISM or MW absorption line masking in the model spectra). For each object, we extract constraints on the underlying model parameters from the posterior (displayed in Table 6), as well as the flux distributions for the fitted lines to examine fit quality (Figure 5). Finally, we also sample the posterior model distribution and extract predicted UV spectra at C iv $\lambda\lambda 1548,1550$ and He ii $\lambda 1640$. Analysis of the UV spectra in a continuum-normalized space minimizes issues with potential flux calibration mismatch between the optical and UV. We fit a cubic spline to the UV spectrum with the stellar wind lines masked both for the data and the model spectra, and normalize the spectra by this continuum before comparison in Figure 6.

Thus far, we have focused on fits to only the optical nebular lines and the equivalent widths of the UV stellar wind lines. This relatively straightforward approach allows us to ignore the complications that reddening and aperture matching uncertainties introduce when fitting both UV and optical line fluxes from HST/COS and SDSS. However, we are also interested in whether the UV O iii] and C iii] lines can be reproduced, and the ratio of these lines provides valuable constraints on the gas-phase C/O ratio. Thus, we perform an additional set of fits with the above-described optical line information along with the flux and equivalent width of both O iii] $\lambda 1661,1666$ and the total flux of the C iii] $\lambda 1908$ doublet (with upper limits employed where lines are undetected or impacted by MW absorption), and allow the C/O abundance to vary over the full range of the models (0.1–1.4). We verify that the strength of the O iii and C iii lines are reasonably well-reproduced by this fit within the measurement uncertainties, confirming that the prominent C iii] detections here can be reproduced by our population synthesis models and that the resulting C/O abundances are reasonable. We present these C/O constraints in Table 5. We also verify that the results of the experiment we present in the following section, comparing nebular line fits with or without including constraints on the UV stellar equivalent widths, produce qualitatively-similar results for the stellar wind lines when the UV nebular lines are included. For simplicity of interpretation, we focus on fits only to the optical nebular line and UV equivalent width fits for the remainder of this paper. However, we utilize the results of the joint optical and UV nebular line fits to first fix the C/O abundance assumed in the population synthesis fits presented and discussed here, adopting the median values presented in Table 5 for each.

The parameter estimates from the optical nebular line fits (Table 6) for the new systems in this paper are generally in reasonable agreement with the fits to the photometry and $W_{0}(\mathrm{H\beta})$ described in Section 3.1 (Table 3). The star formation rates are generally well-matched between the two methods within their uncertainty. The nebular line fits without UV information generally prefer a slightly earlier initiation of the current period of star formation than the photometry fits (median age of 20 Myr compared to 4 Myr), and the inferred stellar masses tend to be correspondingly higher by $\sim 0.3$–0.5 dex. ⁷⁷7With the exception of SB 153, where the nebular line fit prefers a star formation rate 0.4 dex lower and a lower mass by 0.3 dex. The nebular-only fits uniformly prefer a truncation of the current period of star formation within the last Myr, and suggest gas-phase oxygen abundances $12+\log\mathrm{O/H}\simeq 8.3$–$8.5$ modestly larger than those inferred from the direct-$T_{e}$ method (Table 4). With $\xi_{d}\simeq 0.1$–0.4, this places the estimates of the stellar metallicity $Z_{\star}$ at $25$–50% solar. Adding the UV stellar equivalent widths changes this picture, as will be discussed in detail in the following section and Appendix B.

In this subsection, we will analyze the results of the population synthesis fits described above and explore possible explanations for any discrepancies with the observations this uncovers. We will focus first on the ability of the C&B models to reproduce the nebular lines and massive stellar wind signatures in 5.2.1. In 5.2.2 we will compare the metallicities derived using the C&B photoionization models and those computed with the direct-$T_{e}$ method (Section 3.1). Finally, we will discuss the peculiar profile of broad He ii emission detected in two of the strongest emitters in our sample in 5.2.3.

Nebular spectroscopy at cosmological distances, soon to be extended to large samples at $z>6$ with JWST, represents an opportunity to study early star-forming systems in much greater physical detail than photometry alone can afford. However, the accuracy of physical inferences drawn from fitting this emission will be fundamentally limited by the uncertainties in stellar population synthesis models, which still vary substantially with different prescriptions for stellar evolution especially at ionizing energies (e.g. Stanway & Eldridge, 2019). Evidence from both photometric bands contaminated by rest-optical line emission and the first rest-UV spectra in the reionization era suggests that a significant fraction of $z>6$ galaxies are dominated by populations of massive stars formed in bursts initiated within $<20$ Myr of observation (e.g. Stark et al., 2015b; Endsley et al., 2020). Thus, our understanding of the most distant galaxies JWST will observe will depend sensitively on models for the evolution of populations of massive stars on much shorter timescales than typical star-forming galaxies even at $z\sim 2$ (e.g. Steidel et al., 2016).

Study of stars and clusters in the Local Group has provided important clues as to the physics shaping the evolution of very young stellar populations. In particular, the 30 Doradus star-forming region in the LMC and its central young cluster R136 ($\sim 1$–3 Myr, $Z/Z_{\odot}\sim 0.5$) has long represented a benchmark analogue for intense starbursts observed at high redshift. Intriguingly, the resolved stellar contents of this regions reveal evidence for a significant overabundance of massive stars $>30M_{\odot}$ relative to a Salpeter IMF (Schneider et al., 2018), including a sizable population of very massive stars with initial masses 100–320 $M_{\odot}$ in R136 (Crowther et al., 2016). Simulations suggest that such an overpopulation can be readily produced by mass transfer and mergers in binary systems, yielding numerous very massive stars even exceeding the upper limit of the initial mass function on timescales of Myrs (e.g. Schneider et al., 2014; de Mink et al., 2014). These very massive Of stars drive dense winds yielding prominent He ii and C iv (Crowther et al., 2016) and, if they acquire sufficient angular momentum during mass transfer, may undergo quasi-chemically homogeneous evolution and effectively become an extremely massive Wolf-Rayet star. The stars produced by this mechanism are very good candidates for long gamma-ray burst progenitors (e.g. Cantiello et al., 2007), and may have an outsized impact on the properties of recently-formed star clusters. However, without additional resolved star-forming region approaching the age and luminosity of 30 Doradus (e.g. Kennicutt, 1984), the prevalence of these extreme mass transfer products and their impact on the integrated spectra of young star-forming regions as a function of metallicity remains unclear.

While unresolved, young star-forming regions such as those presented in this work provide the next best laboratories in which the evolution of populations of massive stars can be constrained. The systems presented in this paper are by-selection dominated by very young stellar populations, with [O iii]+H$\beta$ equivalent widths and photometric sSFRs at or exceeding the highest values inferred at $z>6$ and suggesting these systems provide a clear window onto very recently-formed $\lesssim 10$ Myr stellar populations (Tables 3 and 4). And also by selection, all reside at sufficiently high metallicities ($Z/Z_{\odot}\sim$ 20–50%) for radiatively-driven stellar wind signatures to be detectable in the integrated spectra alongside strong ionized gas emission. The He ii $\lambda 1640$ stellar wind line is a particularly useful probe of both initially very massive and luminous massive stars rejuvenated by mass transfer (see also Stanway et al., 2020). Reproducing this stellar wind line along with the other available constraints is a crucial litmus test for stellar population synthesis, but matching it at high equivalent width in both local and $z\sim 2$ galaxies has long stymied stellar models (Section 4).

In Section 5, we described an experiment designed to leverage both the nebular and stellar wind lines to search for evidence of stars missing from current stellar population synthesis models. We focus on the C&B stellar population synthesis models, which adopt the latest prescriptions for single-star evolution but do not include binary mass transfer or rotation. While BPASS incorporates a model for the impact of mass transfer on very massive stars, the predicted impact on He ii on $<10$ Myr timescales in this framework is minimal, and insufficient to reach the equivalent widths observed in our local young star-forming regions (Section 4 and Eldridge & Stanway, 2012). Significant recent improvements in wind-driven mass loss predictions have allowed the C&B models to resolve much of the tension with the strength of He ii $\lambda 1640$ relative to older stellar population models (Section 4), so residual issues in matching this line may be a relatively clean probe of additional physics these models are missing at very young ages.

Our comparisons with the C&B models revealed significant tension with the observed star-forming regions. When the nebular lines alone were fit, we found that the models systematically underestimated the strength of both the UV stellar C iv P-Cygni and He ii emission lines (Section 5.2.1). Explicitly including the UV stellar equivalent widths in the fits failed to identify a consistent model able to match both the nebular line measurements and the UV–optical stellar wind strengths, suggesting that the models would require unreasonably-high metallicity stars to reach the observed wind strengths. We discussed other possible explanations for this offset, but find no satisfactory alternative explanation, suggesting that even with their improved treatment of massive stellar winds, the C&B models simply cannot self-consistently reproduce the nebular emission and observed wind line equivalent widths.

The mismatch between the C&B predictions and our data in the UV stellar wind lines is consistent with a significant overabundance of very massive stars relative to the models. In R136, stars in excess of $100M_{\odot}$ completely dominate He ii and contribute strongly to C iv due to the particularly strong winds these luminous stars drive (Crowther et al., 2016). A top-heavy initial mass function could explain this, though the existing evidence for significant IMF variations at subsolar metallicities is not clear-cut (e.g. Gunawardhana et al., 2011; Marks et al., 2012; Hopkins, 2018); but as outlined above, a high incidence of mass transfer and mergers among initially very massive stars would also naturally yield such a distribution. If binary evolution has substantially affected the massive stars present, we would expect many to be rapidly rotating due to angular momentum transfer. The detection of peculiar double-peaked profiles which may be produced by rapidly-rotating stars in two of the most discrepant systems presented lends additional support to the notion that mass transfer occurring on very short timescales may be at work (Section 5.2.3).

Our results suggest that the stellar populations in very young systems with sSFR comparable to bursts observed at $z\sim 6$ may host populations of very massive stars significantly enhanced by mass transfer and mergers. Much of the work on the observational consequences of binary evolution for the SEDs of star-forming galaxies has focused on the impact of these processes on $\sim 10$–100 Myr timescales when numerous initially lower-mass stars ($\lesssim 20M_{\odot}$) become significantly hotter than canonically expected (e.g. Eldridge & Stanway, 2012; Ma et al., 2016; Xiao et al., 2018). However, our difficulty in matching the ultraviolet spectra of these young star-forming regions suggests that binary evolution may profoundly alter the properties and relative number of very massive stars on significantly shorter timescales as well. The fact that this results in detectable differences in the integrated spectra of young star-forming regions has dramatic consequences for attempts to model systems dominated by similar populations at high-redshift.

These data indicate that neither C&B nor BPASS can fully describe the properties of young stellar populations dominated by very massive stars. This suggests that the mapping from metallicity and age to stellar population composition and the emergent ionizing spectrum may be significantly in-error at ages $<10$ Myr. At yet lower metallicities ($<20\%\;Z_{\odot}$) where stellar winds weaken and become increasingly transparent to hard ionizing radiation, the same stellar products of binary evolution will likely have an even more dramatic impact on the integrated ionizing spectrum and thus observed nebular emission. It is crucial that these models be carefully tested at these young ages to mitigate the impact of these systematic uncertainties on our ability to understand the spectra of high-sSFR systems in the early Universe.

The nebular emission and rich continuum signatures accessible simultaneously in extreme local objects such as these likely represent one of our best opportunities to directly constrain models for the evolution of populations of very massive stars. Combined with observations of additional UV wind lines which can provide more detailed leverage on the detailed makeup and age of the massive star population (e.g. Smith et al., 2016; Chisholm et al., 2019), self-consistently reproducing UV–optical spectra for young star-forming regions such as those presented here will be an important benchmark for next-generation stellar population synthesis models incorporating the substantial improvements made both to prescriptions for single and binary star evolution over the past decade. Successfully matching the stellar and nebular features in these systems would be strong evidence that the mapping between stellar metallicity, age, and ionizing spectrum predicted by the models is accurate at very young ages, directly addressing one of the chief sources of systematic uncertainty in the modeling of high-redshift galaxies.

The growing body of rest-frame UV line detections at $z\gtrsim 6$ has ignited substantial interest in locating and understanding nearby galaxies that power similar emission. Particularly puzzling among the four C iii] detections at these redshifts are those in EGS-zs8-1 and z7_GND_42912, two very bright $L^{\star}_{UV}\simeq 3$, 2 reionization-era galaxies with C iii] measured at 22, 16 Å (respectively; Stark et al., 2015a; Hutchison et al., 2019). This is an order of magnitude higher than typical galaxies at similar luminosities at lower redshift (e.g. Shapley et al., 2003; Du et al., 2017). The emerging picture from extensive work at both $z\sim 2-4$ (e.g. Nakajima et al., 2018; Mainali et al., 2019; Du et al., 2020) and in the local $z\sim 0$ Universe (e.g. Rigby et al., 2015; Senchyna et al., 2017, 2019) suggests that such prominent C iii] emission requires a combination of both very high specific star formation rates $\gtrsim 10-100\mathrm{Gyr}^{-1}$ and very metal-poor ($Z/Z_{\odot}\lesssim 0.2$, $12+\log\mathrm{O/H}\lesssim 8.0$) stars and gas. If such low metallicities are indeed necessary, these prominent C iii] detections at $z>6$ imply that these $\gtrsim L^{\star}$ systems reside at sub-SMC metallicities, implying a large shift in the luminosity-metallicity relationship from $z\sim 2$–6.

However, the most extreme of the local galaxies with UV spectroscopy hint at a possible alternate solution. The detection of very strong $\sim 9$– $11$ Å C iii] emission alongside prominent stellar He ii in SB 179 and 191, very young local star-forming regions at $12+\log\mathrm{O/H}=8.3$ ($Z/Z_{\odot}\simeq 0.4$), suggested that very massive stars could potentially power stronger emission than expected at well above SMC metallicities (S17). The spectra presented in this work for seven additional systems selected in the same manner as SB 179 and 191 allow us to test this possibility. In addition to their WR wind signatures in the optical, our targets extend to very high sSFRs, and reach equivalent widths in [O iii] $\lambda 5007$ up to $1719\pm 133$ Å (in SB 153) higher than in any of the targets studied in S17. According to the trend towards more prominent C iii] emission at high sSFR or $W_{0}$([O iii]) displayed by that sample and observed by others (e.g. Mainali et al., 2019), such young effective stellar populations should be capable of powering equivalent widths in C iii] approaching that observed in the reionization era.

In Figure 9 we plot the observed equivalent width of C iii] emission in these systems and other UV observations of $z\lesssim 0.3$ star-forming galaxies. Despite their extraordinarily high equivalent width optical emission lines and correspondingly young effective stellar ages, none of these objects power C iii] emission above 12 Å. Even SB 153, with extremely prominent wind emission at $3.9\pm 0.1$ Å in He ii $\lambda 1640$ and an effective young stellar population age of $\lesssim 5$ Myr (assuming a recent CSFH), only reaches $8.7\pm 0.4$ Å in C iii]. This suggests that we are observing the upper envelope of the expected strength of C iii] at metallicities $12+\log\mathrm{O/H}>8.0$.

Besides metallicity and sSFR, several other factors play a role in regulating C iii] emission. The ISM abundance of carbon has a direct effect on the strength of C iii], and can suppress this emission significantly at very low C/O. However, our galaxies reside at moderate C/O abundances ($-0.7<\log_{10}(\mathrm{C/O})<-0.2$, with median $-0.31$; Table 5), consistent with other local systems at comparable oxygen abundances but significantly higher than the bulk of more metal-poor galaxies at $12+\log\mathrm{O/H}<8$ which extend down to $\log_{10}\mathrm{C/O}=-1.0$ (e.g. Esteban et al., 2014; Berg et al., 2019a). Finding star-forming galaxies at high redshift with substantially higher C/O than these systems is unlikely given this observed pseudo-secondary behavior and low-metallicity plateau in C/O in local systems. We do expect to encounter $\alpha$/Fe enhancement in high-redshift galaxies dominated by a rising star formation history, which can lead to a lower-metallicity stellar population and thus harder ionizing radiation field than expected from the gas-phase oxygen abundance (e.g. Steidel et al., 2016; Strom et al., 2018; Sanders et al., 2020). The presence of lower-metallicity stars may act to further enhance C iii] even at moderate gas-phase oxygen abundance, though the magnitude of the effect will depend on the interplay between both gas ionization shaped by these stars and cooling efficiency determined largely by O/H.

Even for systems dominated by very young stellar populations with prominent optical [O iii] emission, it appears unlikely that C iii] significantly in-excess of 10 Å can be powered for stars and gas at metallicities above $12+\log\mathrm{O/H}>8.0$ ($Z/Z_{\odot}\gtrsim 0.2$). The number of systems at very high sSFR studied at metallicities $12+\log\mathrm{O/H}>8.4$ remains relatively small (Fig. 9), but all are consistent with even weaker C iii] $\lesssim 5$ Å. This is broadly consistent with the picture painted by photoionization models (e.g. Jaskot & Ravindranath, 2016; Nakajima et al., 2018; Plat et al., 2019), which indicate that both softer ionization radiation fields and lower electron temperatures encountered at higher metallicities conspire to suppress collisionally-excited emission from the C iii] doublet. While C iii and O iii have similar ionization potentials and our systems power extremely prominent [O iii] emission, the C iii] doublet is far more sensitive to electron temperature than the optical [O iii] lines (e.g. Jaskot & Ravindranath, 2016; Tang et al., 2020). As a result, C iii] is suppressed more strongly by the drop in temperatures associated with the increased efficiency of metal-line cooling at these high gas-phase metal abundances, producing the relatively sharp evolution with metallicity we observe. The possibility remains that some of the strongest emission $\gtrsim 20$ Å encountered in this line at high-redshift may require additional photoionization and heating from AGN or fast radiative shocks (Nakajima et al., 2018; Le Fèvre et al., 2019; Plat et al., 2019). However, in galaxies with nebular line ratios consistent with star formation, detection of C iii] emission at equivalent widths in-excess of 15 Å can be considered strong positive evidence for metallicities below that of the SMC ($12+\log\mathrm{O/H}\lesssim 8.0$). Along with C iv and He ii emission which become prominent at yet lower metallicities ($12+\log\mathrm{O/H}<7.7$; e.g. Senchyna et al., 2019; Berg et al., 2019b), these strong UV nebular lines may be a powerful probe of early chemical evolution in the JWST era when detection of [O iii] $\lambda 4363$ is intractable due to faintness or when the optical lines are entirely inaccessible at the highest redshifts.