The Evolution and Properties of Rotating Massive Star Populations

TL;DR

This study models the evolution of rotating massive stars and their impact on ionizing radiation, revealing the importance of fast rotators in sustaining ionizing photon production in young stellar populations, with implications for cosmic reionization.

Contribution

It introduces new stellar population synthesis models incorporating rotation effects and very massive stars, extending previous models to better match observations and high-redshift galaxy data.

Findings

Fast-rotating stars significantly influence ionizing photon output.

Models align with observations of young nuclear star clusters.

Rotation effects are crucial for understanding high-redshift galaxy emission lines.

Abstract

We investigate the integrated properties of massive (>10 Msun), rotating, single-star stellar populations for a variety of initial rotation rates (v/vcrit=0.0, 0.2, 0.4, 0.5, and 0.6). We couple the new MESA Isochrone and Stellar Tracks (MIST) models to the Flexible Stellar Population Synthesis (FSPS) package, extending the stellar population synthesis models to include the contributions from very massive stars (>100 Msun), which can be significant in the first ~4 Myr after a starburst. These models predict ionizing luminosities that are consistent with recent observations of young nuclear star clusters. We also construct composite stellar populations assuming a distribution of initial rotation rates. Even in low-metallicity environments where rotation has a significant effect on the evolution of massive stars, we find that stellar population models require a significant contribution from fast-rotating (v/vcrit>0.4) stars in order to sustain the production of ionizing photons beyond a few Myr following a starburst. These results have potentially important implications for cosmic reionization by massive stars and the interpretation of nebular emission lines in high-redshift star-forming galaxies.

Equations

$$v(r)=v_{\infty }{\left(1-\frac{R}{r}\right)}^2,$$

$$R_{\mathrm{eff},\mathrm{WR}}=R_{\ast }+\frac{3\kappa |\dot{M}|}{8\pi v_{\infty }}.$$

$$\frac{d{p}_{\mathrm{rad}}}{dt}=(1+{\tau }_{\mathrm{IR}})\frac{L_{\mathrm{bol}}(t)}{c},$$

$$\frac{d{p}_{\mathrm{wind}}}{dt}=\dot{M}{v}_{\mathrm{wind}},$$

$$\begin{array}{r}\log (v_{\mathrm{wind}})[\mathrm{km}{\mathrm{s}}^{-1}]=1.23-0.30\log (L/L_{\odot })\\ +0.55\log (M/M_{\odot })+0.64\log (T_{\mathrm{eff}}[K])\\ +0.13\log (Z/Z_{\odot }).\end{array}$$

$$p_{\mathrm{SNII}}=\frac{2E_{\mathrm{SNII}}}{v_{\mathrm{SNII}}},$$

$$\frac{d{p}_{\mathrm{SNII}}}{dt}=p_{\mathrm{SNII}}\frac{d{N}_{\mathrm{SNII}}}{dt},$$

Full text

1. Introduction

Massive stars, though rare in number, are energetically dominant across a wide range of environments spanning star cluster to extragalactic scales. Stellar feedback from young, massive stars in H ii regions and star-forming galaxies is a crucial yet poorly understood process. Broadly speaking, stellar feedback refers to the deposition of energy, momentum, mass, and nuclear-burning products via channels that include type I and type II supernovae (SNe), stellar radiation, and winds (e.g., Murray et al., 2010; Hopkins et al., 2011). These processes influence the state of the interstellar medium (e.g., McKee & Ostriker, 1977), thereby regulating star formation (e.g., Williams & McKee, 1997; Mac Low & Klessen, 2004; McKee & Ostriker, 2007) and driving both turbulence and galactic-scale outflows (e.g., Dekel & Silk, 1986; Martin, 1999; de Avillez & Breitschwerdt, 2004; Joung & Mac Low, 2006; Oppenheimer & Davé, 2006; Agertz et al., 2009; Tamburro et al., 2009; Hopkins et al., 2012; Creasey et al., 2013).

Observations of star-forming giant molecular clouds (GMCs) have suggested that early feedback processes disperse the dense gas well before the first SNe explode, which may increase the overall efficiency of feedback and reduce the star formation efficiency (e.g., Evans et al., 2009; Krumholz et al., 2012). Proposed mechanisms of this early feedback include the destruction of the dense molecular regions by expanding H ii bubbles (e.g., Whitworth, 1979; Matzner, 2002; Walch et al., 2012; Lopez et al., 2014) and radiation pressure (e.g., Murray et al., 2005; Krumholz & Matzner, 2009). The complex interplay between the properties of young stellar populations and the dissolution of their birth gas has significant implications for the number of photons that are able to leak out of the host galaxies (e.g., Dove & Shull, 1994; Gnedin, 2000; Ma et al., 2015) and drive cosmic reionization (e.g., Haardt & Madau, 1996).

The influence of stellar feedback in a galactic context has been explored by many groups, both analytically and in simulations (e.g., Haehnelt, 1995; Murray et al., 2005; Nath & Silk, 2009; Murray et al., 2011; Hopkins et al., 2011, 2012; Agertz et al., 2013; Ma et al., 2015; Muratov et al., 2015; Fierlinger et al., 2016). Prescriptions for the time-dependent momentum and energy deposition of stellar populations require stellar population synthesis (SPS) models such as Starburst99 (SB99; Leitherer et al. 1999; Vázquez & Leitherer 2005; Leitherer et al. 2010, 2014), which in turn rely on stellar evolution models, e.g., the Geneva models (Ekström et al., 2012; Georgy et al., 2013) and realistic stellar spectra, especially in the EUV. Stellar population models that include the effects of stellar rotation lead to an overall increase in feedback efficiency (Levesque et al., 2012; Leitherer et al., 2014) because rotation makes stars hotter, brighter, and longer-lived (e.g., Maeder & Meynet, 2000). In other words, rotating stellar populations have higher bolometric luminosities and sustain harder radiation fields (i.e., more ionizing photons) over a longer period of time compared to non-rotating stellar populations.

Binary interaction is another important and complex aspect of massive star evolution. A large fraction of O- and B-type stars are found in binary or higher multiplicity systems, and an estimated $\gtrsim 70\%$ of O-type stars are believed to undergo mass exchange, a third of which likely end up in a binary merger product (e.g., Chini et al. 2012; Sana et al. 2012, 2013; de Mink et al. 2014; but see also Kobulnicky et al. 2014, where the authors found that the multiplicity fraction is closer to $\sim 55\%$ for orbital periods less than 5000 days, and that it likely depends on the orbital period and separation). Broadly speaking, binarity and rotation have similar effects on stellar evolution. This is a fortunate aspect since large grids of stellar evolution models generally account for only one of these effects due to the sheer size of the parameter space that would otherwise need to be explored (see Song et al. 2016 for a recent example that accounts for tidal interactions in a binary system of rotating stars). We focus on single star models with rotation in this work and we compare with predictions from the Binary Population and Spectral Synthesis model (BPASS; Eldridge & Stanway 2009).

In this work, we investigate two features of single massive stars in the context of SPS modeling. First, we explore the contribution from very massive stars (VMSs) to the integrated stellar population properties by extending the initial mass function (IMF) upper mass limit from the fiducial value of $100~{}{M_{\odot}}$ to $300~{}{M_{\odot}}$. Recent observations of star clusters have suggested the need to include VMSs in models (e.g., Crowther et al., 2016; Smith et al., 2016), though a top-heavy IMF has also been proposed (e.g., Turner et al., 2015).

Second, we explore the effects of rotation over a range of rotation rates. We build on the previous work by Levesque et al. (2012), who investigated the effect of stellar rotation on the resulting stellar population properties. Their conclusions were based on a single rotation value (surface velocity set to $40\%$ of the critical, or break-up, velocity, $v/v_{\rm crit}=0.4$), though the authors also considered a composite population by weighting $30\%$ non-rotating and $70\%$ rotating stellar populations. We also investigate the effects of stellar populations harboring a range of rotations rates, which is supported by observations (Huang et al., 2010), and explore their significance in the context of cosmic reionization. One of the distributions considered in this work includes a tail of fast-rotators with $v/v_{\rm crit}>0.4$, which may be of binary origin (Ramírez-Agudelo et al., 2013).

Both VMSs and a distribution of rotation rates have interesting implications in light of recent observational and theoretical work. Smith et al. (2016) advocated for VMSs in lieu of rapid rotators to resolve the discrepancy between photometric and spectroscopic ages inferred for a young stellar population in NGC 5253. The authors reasoned that large numbers of rapid rotators, though able to explain the observed nitrogen enhancement and ionizing flux, are in tension with the moderate rotational velocities reported in a recent survey of single O-type stars in 30 Doradus (Ramírez-Agudelo et al., 2013). We explore whether a combination of moderate rotation rates and VMSs can provide a solution to this problem. Furthermore, Steidel et al. (2014) analyzed rest-frame optical spectra of $z\sim 2\textrm{--}3$ star-forming galaxies and found that their position in the “Baldwin, Phillips, and Terlevich” (BPT) diagram—a nebular emission line diagnostic first introduced in Baldwin et al. (1981)—is offset relative to their $z=0$ counterparts. The authors concluded that this shift could be explained by a harder stellar ionizing radiation field ($T_{\rm blackbody}=50000\textrm{--}60000$ K), higher ionization parameter, and nitrogen-enhanced nebular gas, suggesting rapid stellar rotation or binaries to explain these observations. In subsequent work, Steidel et al. (2016) found that the BPASS binary star models can self-consistently explain the observed line ratios. This, together with the finding that binary models can sufficiently boost the escape fraction of ionizing photons, $f_{\rm esc}$, to explain cosmic reionization (Ma et al., 2016b), are based on a unique feature of binary models: the ability to sustain sufficiently long-lived ($\gtrsim 5$ Myr) far-UV stellar radiation field. We investigate whether this can be explained by rotating stellar populations, particularly those with contributions from fast-rotators.

The paper is organized as follows. In Sections 2 and 3, we describe the MESA Isochrones and Stellar Tracks (MIST) and Flexible Stellar Population Synthesis (FSPS) models. Next, in Section 4, we present the results from stellar population modeling, including an investigation of the effects of VMSs, metallicity, and rotation distributions. In Section 5, we discuss the implications of our models in the context of the cosmic reionization and high-redshift star-forming galaxies. We conclude with a summary in Section 6.

2. MIST Models

In this section we provide a brief overview of the MIST models, focussing on the details that are most salient to the evolution of massive, rotating stars. For a comprehensive description and summary, the reader is referred to Section 3 and Table 1 in Choi et al. (2016).

2.1. Modules for Experiments in Stellar Astrophysics (MESA)

The stellar evolutionary tracks are computed using MESA,111http://mesa.sourceforge.net/ an open-source stellar evolution package (Paxton et al., 2011, 2013, 2015). MESAstar, its 1D stellar evolution module, solves the fully coupled Lagrangian equations of stellar structure and composition. Some of the key advantages of MESA include its robust numerical methods, its modular structure that easily enables a user to adapt it to a wide range of problems in stellar astrophysics, and parallelization via OpenMP. The reader is encouraged to consult the original MESA instrument papers and Section 2 in Choi et al. (2016) for more details.

2.2. Abundances

All of the models computed for this work are initialized with the solar-scaled abundances and $Z=Z_{\odot,\rm protosolar}=0.0142$ from Asplund et al. (2009). The helium mass fraction, $Y$, is computed assuming a linear enrichment law from the primordial helium abundance $Y_{\rm p}=0.249$ (Planck Collaboration et al., 2015) to the protosolar helium abundance $Y_{\odot,\rm protosolar}=0.2703$ (Asplund et al., 2009), where $\Delta Y/\Delta Z=1.5$.

Although $\rm[Fe/H]$ ranges from $-4.0$ to $0.5$ in the full set of grids published on the MIST website,222http://waps.cfa.harvard.edu/MIST/ we focus on the range $-2.0\leq\rm[Fe/H]\leq 0.0$ in this work.

2.3. Rotation

The effects of rotation appear in MESA stellar evolution calculations in three main ways. First, rotation decreases the gravitational acceleration $g$ via the centrifugal force, which in turn affects the stellar structure, making the star appear more extended and cooler near the equator. Second, rotation can promote extra mixing in the interior, providing a boost to the transport of chemicals and angular momentum. As a result, the helium fraction—and hence the mean molecular weight $\mu$—is increased in the surface layers and more fuel can be introduced to the convective core, resulting in a hotter, brighter, and more long-lived star. MESA adopts the common approach of treating the chemical and angular momentum transport in a diffusion approximation (Endal & Sofia 1978; Zahn 1983; Pinsonneault et al. 1989; Heger et al. 2000 but see also Maeder & Meynet 2000; Eggenberger et al. 2008; Potter et al. 2012). Third, rotation enhances mass loss (e.g., Heger et al., 2000; Brott et al., 2011; Potter et al., 2012). MESA adopts the formulation from Langer (1998), where the mass loss rate $\dot{M}$ is multiplied by a factor that increases dramatically as the surface angular velocity $\Omega$ approaches critical, or break-up, angular velocity, $\Omega_{\rm crit}=[(1-(L/L_{\rm Edd}))(GM/R^{3})]^{0.5}$.333Note that this expression assumes that the surface brightness is uniform over the stellar surface. Maeder & Meynet (2000) worked out an alternative expression that takes into account gravity darkening effects, but the difference is small unless $L/L_{\rm Edd}$ is much larger than $\sim 0.65$. Because $\Omega_{\rm crit}$ depends on the ratio of the bolometric luminosity to the Eddington luminosity, $L/L_{\rm Edd}$, massive stars near the Eddington limit can experience a sizable increase in their mass loss rates even with small amount of rotation, resulting in a prompt removal of their surface hydrogen layer (e.g., Maeder & Meynet, 2000; Choi et al., 2016).

Models are initialized following solid-body rotation near the zero-age main sequence (ZAMS), similar to the approach adopted by many other stellar evolution codes (Pinsonneault et al., 1989; Heger et al., 2000; Eggenberger et al., 2008). We compute models with $v/v_{\rm crit}=0.0$, 0.2, 0.4, 0.5, and 0.6 at ZAMS, where $v/v_{\rm crit}$ corresponds to the ratio of the stellar rotation velocity to the critical rotation velocity at the surface and $v_{\rm crit}=[(1-(L/L_{\rm Edd}))(GM/R)]^{0.5}$.

2.4. Mixing

Convective mixing of elements is treated as a time-dependent diffusive process with the diffusion coefficient provided by the mixing length formalism of Henyey et al. (1965) and assuming $\alpha_{\rm MLT}=1.82$ obtained from the solar calibration (Choi et al., 2016). Convectively unstable zones are identified according to the Ledoux criterion, which is similar to the Schwarzschild criterion but also includes composition effects. The effects of convective overshoot mixing beyond the fiducial Ledoux boundary are taken into account following Freytag et al. (1996) and Herwig (2000), which assumes that the diffusion coefficient decays exponentially with distance from the convective boundary. A free parameter that controls the efficiency of this extra mixing in the core, $f_{\rm ov}=0.016$, is empirically calibrated to reproduce the main-sequence turnoff (MSTO) morphology in M67 (Choi et al. 2016, but see also Magic et al. 2010). Semiconvection and thermohaline mixing are also included, where we adopt the efficiency parameters $\alpha_{\rm sc}=0.1$ and $\alpha_{\rm th}=666$. Finally, we account for five rotationally induced instabilities that lead to chemical and angular momentum transport, namely the dynamical shear instability, secular shear instability, Solberg-Høiland instability, Eddington-Sweet circulation, and Goldreich-Schubert-Fricke instability (Heger et al., 2000; Paxton et al., 2013). We do not include mixing due to internal magnetic fields generated by a Tayler-Spruit dynamo (Spruit, 2002).

2.5. Mass Loss

In hot stars, mass loss is driven by the momentum transfer from ultraviolet photons to metal ions in the atmosphere (Lucy & Solomon, 1970; Castor et al., 1975) and is therefore metallicity-dependent. Mass loss in massive stars is separated into three categories in the MIST model calculations. For stars with $T_{\rm eff}>1.1\times 10^{4}$ K and surface hydrogen mass fraction $X_{\rm surf}>0.4$, we adopt the mass loss prescription from Vink et al. (2000, 2001). If the star loses a considerable amount of its outer hydrogen layer ($X_{\rm surf}<0.4$) and becomes a Wolf-Rayet (WR) star, we use the wind prescription of Nugis & Lamers (2000) instead. Lastly, we use the de Jager et al. (1988) prescription for all stars that have effective temperatures below $10^{4}$ K, including red supergiant (RSG) stars. As discussed in 2.3, mass loss is enhanced in rotating models by a factor that depends on the rotation rate and the Eddington ratio.

2.6. Tracks and Isochrones

We compute MIST isochrones for five rotation rates—$v/v_{\rm crit}=0.0,0.2,0.4,0.5$, and $0.6$—for each metallicity. The stellar evolutionary tracks ranging from 0.1 to $300~{}{M_{\odot}}$ are computed by MESA and processed into isochrones following the procedure outlined in Dotter (2016). We modified the mass sampling in the isochrones in this work compared to that in the standard MIST grids to ensure that fast evolutionary phases at early times are particularly well-represented. Figures 1 and 2 show 1, 3, 5, and 10 Myr isochrones at $\rm[Fe/H]=-1.0$ and $0.0$, respectively. We omit $v/v_{\rm crit}=0.5$ from these figures for display purposes. There are several notable features. First, the rotating models tend to be hotter and more luminous overall. Second, WR stars from very massive progenitors appear between 1 Myr and 5 Myr, particularly for $v/v_{\rm crit}>0.2$. Third, some isochrones, e.g., $v/v_{\rm crit}=0.6$ at $\rm[Fe/H]=-1.0$, show signatures of quasi-chemically homogeneous evolution (QCHE; Maeder, 1987; Langer, 1992). QCHE occurs in massive, fast-rotating stars where the mixing timescale becomes shorter than or comparable to the nuclear burning timescale. In this scenario, nuclear burning products—mostly helium—from the core are mixed into the outer layers and a large fraction of the star undergoes nuclear fusion as fresh fuel is channeled into the core. Overall, the star becomes hotter due to the reduced mean opacity, and brighter and more long-lived due to enhanced mixing in the core (see Figures 5, 6, and 7 in Brott et al. 2011 for examples). QCHE occurs more readily at low metallicity due to the reduced angular momentum loss via stellar winds (Yoon & Langer, 2005; Woosley & Heger, 2006) and the more compact stellar structure resulting in a decreased mixing timescale ($\tau\sim R^{2}/D$, where $R$ is the radius of the star and $D$ is the diffusion coefficient for mixing). To verify that both processes indeed contribute, we performed a simple test where we ran a total of eight $60~{}{M_{\odot}}$ models, with and without rotation, with and without mass loss, at solar and one-tenth-solar metallicities. We found that even in the absence of mass loss, i.e., no angular momentum loss, the MS lifetime enhancement due to rotation was still larger for the low metallicity model ($\sim 20\%$) compared to the solar metallicity model ($\sim 10\%$), which suggests that the compactness of the star at low metallicity is indeed an important effect. In the rotating case, there was a small ($\lesssim 5\%$) enhancement in the MS lifetime for the models without any mass loss compared to those with mass loss. An extensive study examining the relative importance of the two effects as a function of metallicity, rotation rate, mass, etc. is beyond the scope of this paper, but would be worth exploring further.

In Figures 3 and 4 we compare the $v/v_{\rm crit}=0.0$ and 0.4 MIST isochrones with the non-rotating and rotating Geneva isochrones, zooming in around the MSTO for clarity. In detail, there are clear morphological differences between the MIST and Geneva isochrones, but both their evolution with time and their modifications due to rotation are qualitatively similar. A notable difference is the magnitude of the effect of rotation on the MS lifetimes and the trajectory in the HR diagram at $t\gtrsim 5$ Myr; rotating Geneva models are much brighter than the non-rotating models compared to their MIST counterparts.

In order to investigate the effect of rotation on the MS lifetimes of massive stars in more detail, we plot the MS lifetime–initial mass relations in Figure 5. The black dotted and dot-dashed lines correspond to the non-rotating and rotating Geneva (Ekström et al., 2012; Georgy et al., 2013) evolutionary tracks for $Z=0.002$ ($\rm[Fe/H]=-0.86$ assuming their adopted $Z_{\odot}=0.014$) and $Z=0.014$. At a fixed initial mass, a higher rotation rate lengthens the MS lifetime due to more efficient rotational mixing. Note that these plots demonstrate that the MS lifetime boost at a fixed initial mass is larger in the Geneva models compared to that in the MIST models (see also Figure 20 in Choi et al. 2016), suggesting that rotational mixing may be more efficient in the former. The MS lifetime flattens out at the high mass end since $t\propto L/M$ and luminosity becomes a shallower function of mass: for $L\propto M^{\alpha}$, we find $\alpha\sim 4$ for $M\lesssim 10~{}{M_{\odot}}$ and $\alpha\sim 1.5$ for $M\gtrsim 100~{}{M_{\odot}}$. At solar metallicity, there is an interesting non-monotonic behavior around $80~{}{M_{\odot}}$. Rapidly rotating stars with masses greater than $80~{}{M_{\odot}}$ begin to evolve more vertically and blueward in the HR diagram during the MS.444Although quasi-chemically homogeneous evolution appears more readily in low-metallicity environments, they may also occur in metal-rich systems. See Martins et al. (2013) for more details. Due to their high luminosities and temperatures and therefore heavy mass loss, they quickly shed their H-rich envelopes and become WR stars. As they continue to lose significant mass due to heavy WR winds, they effectively “reset” as lower-mass stars, which prolongs their MS lifetimes. In other words, the MS lifetime–initial mass relation systematically shifts upward to another branch beyond $80~{}{M_{\odot}}$.

3. Stellar Population Synthesis Models

3.1. Flexible Stellar Population Synthesis

Models of simple stellar populations (SSPs) in this work are computed using the Flexible Stellar Population Synthesis package (FSPS; Conroy et al., 2009; Conroy & Gunn, 2010). The primary stellar spectral library consists of the MILES empirical library (Sánchez-Blázquez et al., 2006; Falcón-Barroso et al., 2011), which is supplemented with the CMFGEN WR spectra (available for both WN and WC subtypes) from Smith et al. (2002) and the WM-Basic (Pauldrach, 2012) hot star spectra (J.J. Eldridge, priv. comm.). The WR spectra are assigned to points in the isochrone that are first identified as WR then further categorized into WR subtypes, e.g., WN and WC, according to the surface composition: stars with surface C/O ratio $>1$ and $\leq 1$ are labeled as WC and WN, respectively. The WM-Basic models are applied for MS stars with $T_{\rm eff}>2.5\times 10^{4}~{}\rm K$.

For all subsequent analyses, we simulate a population of total mass $10^{6}~{}{M_{\odot}}$ following an instantaneous burst of star formation. All integrated quantities are computed following the Kroupa 2001 initial mass function (IMF), assuming that the IMF is fully sampled. The IMF lower and upper mass limits are set to $0.08~{}{M_{\odot}}$ and $300~{}{M_{\odot}}$, respectively (see Section 4.4 for a discussion). Stochastic IMF sampling effects are likely to have a bearing on real stellar populations, especially in low-mass systems (e.g., Cerviño et al., 2000, 2001; da Silva et al., 2012, 2014). We revisit this point briefly in Section 3.2.

As discussed in Maeder (1990), assigning a spectrum to a WR star in the isochrone is a non-trivial task due to the high optical depth of its wind. As a result, the effective photosphere—the $\tau=2/3$ surface—has a radius larger than the hydrostatic radius reported in the stellar evolutionary track, and the inferred temperature from observations, $T_{\rm eff,\;WR}$, is cooler than the temperature from stellar evolution calculations, $T_{*}$. We follow the standard approach adopted by Maeder (1990) and Smith et al. (2002), where $T_{\rm eff,\;WR}$ is estimated via a weighted sum of $T_{*}$ and $T_{\rm eff,\;wind}$. $T_{\rm eff,\;wind}$ is computed by assuming a velocity law:

[TABLE]

where $v_{\infty}$ is the terminal velocity, set to $3\times 10^{8}~{}\rm cm\;s^{-1}$.

Starting with the definition of optical depth $d\tau=-\kappa\rho\;dr$ and integrating inward to $\tau=2/3$, the effective photosphere is obtained

[TABLE]

Since $L$ is constant, it is trivial to compute $T_{\rm eff,\;WR}$ in terms of $R_{\rm eff,\;WR}$, $R_{*}$, and $T_{*}$. Then the temperature used to retrieve the spectrum for a WR star with $T_{*}$ is $0.6T_{*}+0.4T_{\rm eff,\;WR}$. Smith et al. (2002) chose this weighting scheme to reproduce the range of observed temperatures of the Galactic and LMC WR population. In principle, the weighting factors should be re-derived for this work given the differences in the underlying stellar evolutionary models adopted in Smith et al. (2002) (Meynet et al. 1994 tracks) and MIST. However, the objective here is a simple assessment of the uncertainties introduced due to the ill-defined WR temperatures. We revisit this point in Section 4.2.

3.2. IMF Stochasticity

We evaluate the assumption that the IMF is smoothly sampled for this work. A key concern regarding IMF stochasticity in the context of SPS modeling is that many quantities of interest, e.g., ionizing flux, are dominated by a very small number of the most massive stars, a subpopulation that is also the rarest and thus most prone to sampling effects. Thus two stellar populations that are otherwise completely identical in parameter space (e.g., metallicity, total mass, age) may appear quite different depending on the stellar mass distribution. Figure 6, which shows the fractional flux contribution from stars in different mass ranges, provides a sense of how much variation there may be in spectral energy distributions (SEDs) due to the IMF sampling effects.

There are alternative approaches to SPS modeling (e.g., Barbaro & Bertelli, 1977; Cerviño et al., 2000; Villaverde et al., 2010; Eldridge, 2012; da Silva et al., 2012, 2014) that simulate the IMF stochastic sampling effects. Cerviño et al. (2002) constructed a statistical formalism to quantify the uncertainties due to IMF sampling errors and applied their results to model a starburst population from birth to $20~{}\rm Myr$. Overall, they found that a minimum total stellar mass of $10^{5}~{}{M_{\odot}}$ is necessary to ensure a relative dispersion of less than 10% in hydrogen and helium ionizing flux, though during the “WR phase” (between 2 and 5 Myr), this threshold mass is increased to $10^{6}~{}{M_{\odot}}$.

Since we are operating under the assumption that the total stellar mass of the cluster modeled in our work is $10^{6}~{}{M_{\odot}}$, random IMF sampling effects are likely unimportant for our conclusions. However, the exact value of this critical cluster mass depends on the adopted IMF (e.g., functional form, lower and upper mass limits). The conclusions from Cerviño et al. (2002) were based on a Salpeter 1955 IMF over the mass range $2\textrm{--}120~{}{M_{\odot}}$. Even though the increased mass range ($0.08\textrm{--}300~{}{M_{\odot}}$) in our work would push the threshold total mass upward, an ensemble of such clusters in e.g., a starburst galaxy, should collectively smear out the sampling effects (see e.g., Cerviño et al., 2000; Eldridge, 2012).

4. Stellar Population Properties

4.1. Spectral Energy Distributions (SEDs)

Figures 7 and 8 show SEDs from 1 to 20 Myr for $\rm[Fe/H]=-1.0$ and 0.0, respectively. The colored curves represent four different initial rotation rates for FSPS+MIST, while the gray and black curves show the BPASS v2.0555IMF slopes of $-1.30$ for $0.1<M_{i}/{M_{\odot}}<0.5$ and $-2.35$ for $0.5<M_{i}/{M_{\odot}}<300$. (Stanway et al., 2016) and SB99+Geneva (Leitherer et al., 2014) predictions666IMF slopes of $-1.3$ for $0.1<M_{i}/{M_{\odot}}<0.5$ and $-2.3$ for $0.5<M_{i}/{M_{\odot}}<100$. for comparison where available. The $v/v_{\rm crit}=0.5$ FSPS+MIST model is excluded for clarity. The three blue vertical lines mark the wavelengths of photons capable of ionizing hydrogen, singly ionizing helium, and doubly ionizing helium (912 Å, 504 Å, and 228 Å). For clarity, we do not include the BPASS single star or SB99+Geneva non-rotating models, but they are in broad agreement with the FSPS+MIST non-rotating models.

The FSPS+MIST model predicts a harder spectrum compared to the SB99+Geneva model, partially due to differences in the underlying isochrones. Moreover, the FSPS+MIST model presented here is slightly more metal-poor than the SB99+Geneva model. The large discrepancy at 3 Myr is due to the lack of stars more massive than $300~{}{M_{\odot}}$ in the SB99+Geneva models.777The Geneva group has published stellar evolutionary tracks of VMSs (120 to $500~{}{M_{\odot}}$) at solar, LMC, and SMC metallicities (Yusof et al., 2013), but the corresponding isochrones are not incorporated into SB99 at the time of writing. With the exception of $v/v_{\rm crit}=0.6$, the FSPS+MIST models cease to produce an appreciable amount of photons blueward of 228 Å beyond $\sim 6$ Myr, while the BPASS binary models continue to output significant EUV flux even at very late times (20 Myr). This can be understood by recalling that single star models rely exclusively on the most massive stars to produce ionizing photons. As a result, the principle ionizing sources vanish once the most massive stars disappear after the first few Myrs. In contrast, binary interactions and mergers can persist even after the most massive stars disappear, so binary models are capable of generating ionizing photons at late times.

4.2. Ionizing Photons

Starting from the SEDs as a function of age, $F_{\lambda}(t)$, we can compute the time evolution of the ionizing photon luminosity, $Q(t)$. Figures 9 and 10 show $Q(t)$ from 1 to 30 Myr for a $10^{6}~{}{M_{\odot}}$ population at $\rm[Fe/H]=-1.0$ and $0.0$. The top, middle, and bottom panels show photons with wavelengths below 912 Å, 504 Å, and 228 Å, respectively. We also show the BPASS binary and SB99+Geneva rotating model predictions for comparison. Note that the metallicities are not perfectly matched: for BPASS and SB99+Geneva, $Z=0.002$ and 0.001 for Figure 9 and $Z=0.014$ and 0.02 for Figure 10.

All of the models produce comparable hydrogen-ionizing photon output rates until $\sim 3~{}\rm Myr$. Note that $Q$ is systematically offset in the SB99+Geneva model at very early times due to the lack of stars more massive than $300~{}{M_{\odot}}$. The discrepancies between different models become more pronounced at higher energies: for photons blueward of 228 Å, the differences are substantial as early as $\rm 2~{}Myr$. This is because $Q$ becomes much more sensitive to the shape of the SED at progressively shorter wavelengths, where the predictions are sensitive to fast evolving stars and details of hot star SEDs. Overall, the SB99+Geneva models produce the softest radiation field (see also Figures 7 and 8). The FSPS+MIST models produce copious amounts of helium-ionizing photons after 2 Myr when the WR stars begin to appear, which is several orders of magnitude larger than the SB99+Geneva models. Again, the binary effects are responsible for the prolonged production of ionizing photons at very late times.

As described in Section 3.1, assigning an appropriate spectrum to a WR star in the isochrone is a nontrivial task due to the high optical depth of the WR wind beyond the standard hydrostatic surface computed by stellar evolutionary codes. We tested the effect of modifying the WR temperature according to the Smith et al. (2002) weighting scheme, which is the default choice in SB99. Consistent with what Levesque et al. (2012) found, the difference is negligible for the hydrogen ionizing luminosity but becomes more pronounced for harder photons, which may explain the large discrepancy between the SB99+Geneva and FSPS+MIST predictions for helium-ionizing photons. Since this pertains to WR stars only, any variation due to the choice of WR temperatures disappears after $\sim 5\textrm{--}6$ Myr. We proceed with our default WR $T_{\rm eff}$ assignment, but emphasize that the choice of this weighting scheme (or the lack thereof) will introduce some variation to the predicted output of the most energetic photons.

4.3. Momentum Output

Here we consider the momentum output from massive stars. Figures 11 and 12 show the cumulative (top) and instantaneous (bottom) specific momentum injected in radiation, wind, and type II supernovae. Momentum injected via radiation is computed as follows:

[TABLE]

where $\tau_{\rm IR}$ is the infrared optical depth. In detail there is also a factor of $(1-\exp(-\tau_{\rm UV/optical}))$, but given the large optical depth of UV/optical photons, it reduces to a factor of order unity (e.g., Agertz et al., 2013). The $\tau_{\rm IR}$ term accounts for multiple absorption and re-radiation of the infrared photons in a very optically thick medium. Its preferred fiducial value is still under debate (e.g., Murray et al., 2010; Andrews & Thompson, 2011; Hopkins et al., 2011; Krumholz & Thompson, 2012). For simplicity, we assume $\tau_{\rm IR}=1$ in order to isolate the effect of SPS models alone, but empirical relations between $\tau_{\rm IR}$ and cluster mass (see e.g., Figure 3 in Agertz et al. 2013) suggest that $\tau_{\rm IR}$ can be as large as $\sim 100$ for a $10^{6}~{}{M_{\odot}}$ cluster.

The momentum injection rate from a stellar wind is obtained using the mass loss rate and the wind speed:

[TABLE]

The wind speed is estimated using the relation adopted in SB99 (Leitherer et al., 1992), which is comparable to the escape speed from the stellar photosphere:

[TABLE]

Finally, the momentum deposition per type II SN (SNII) event is

[TABLE]

where we assume that a typical SNII explosion outputs $10^{51}~{}\rm erg$ of kinetic energy (thermalization via shocks) and the ejecta are launched at $v_{\rm SNII}=3\times 10^{8}~{}\rm cm\;s^{-1}$ (Dale, 2015). Moreover, we assume that only stars with initial masses between $8~{}{M_{\odot}}$ and $40~{}{M_{\odot}}$ terminate their lives as typical SNII and stars more massive than $40~{}{M_{\odot}}$ directly collapse to a black hole (Fryer 1999, but see also e.g., Sukhbold & Woosley 2014). The corresponding momentum injection rate is

[TABLE]

where the supernovae rate is obtained by integrating the IMF weight over the most massive star still alive in the previous and current time steps divided by the time interval.

We compare FSPS+MIST to SB99+Geneva, which is a popular choice for stellar population model in many galaxy simulations that attempt to include the effects of stellar feedback (e.g., Hopkins et al., 2014; Agertz & Kravtsov, 2015). We use the SB99 outputs directly for the wind momentum and SNII rates, but we compute ourselves the radiation momentum using the bolometric luminosity reported by SB99 along with Equation 3, and convert the SNII rates to the SNII momentum injection rate using Equations 6 and 7. The first SN explosion occurs at an earlier time in the SB99+Geneva model because their rate calculation assumes that all stars above $8~{}{M_{\odot}}$ explode as SNII. For the FSPS+MIST model, the onset of SNII is delayed at higher rotation rates due to the increased lifetimes (see also Figure 5). Radiative momentum is generally dominant, though this may be an underestimate since we set the $\tau_{\rm IR}$ enhancement factor to unity. At a fixed rotation rate, the metal-poor population ($\rm[Fe/H]=-1.0$) outputs more radiative momentum compared to the metal-rich population, because metal-poor stars tend to be hotter and brighter.

4.4. The Effects of Very Massive Stars

The IMF dictates how frequently stars within a certain mass range occur in nature. Understanding the origin of and quantifying the high-mass IMF slope and the cutoff mass (and their dependence on environmental factors such as metallicity) are active areas of research (e.g., Krumholz et al., 2011; Kroupa et al., 2013; Hopkins, 2013; Andrews et al., 2013; Narayanan & Davé, 2013; Dib, 2014; Weisz et al., 2015). There are ongoing efforts to constrain these properties, but it is a difficult task due to their short lifetimes, rare numbers, and stochastic sampling effects.

Here we investigate the importance of VMSs ($M>100~{}{M_{\odot}}$) in SPS modeling by comparing two otherwise identical models with different cutoff masses. The default option in SB99, a widely-used SPS package that is commonly paired with the Geneva stellar evolutionary models (Ekström et al., 2012; Georgy et al., 2013; Yusof et al., 2013), accounts for stars with initial masses up to $100~{}{M_{\odot}}$, though the underlying Geneva models are available up to $500~{}{M_{\odot}}$. The BPASS models are available with the upper mass limit set to either $100~{}{M_{\odot}}$ or $300~{}{M_{\odot}}$. In a recent observational work, Smith et al. (2016) found that VMSs, rather than less massive but very fast-rotating stars, are necessary to simultaneously explain the large ionizing flux and spectral emission line features observed in the nuclear star clusters in a blue compact dwarf galaxy NGC 5253. They stressed the need for SPS models to include stars more massive than $100~{}{M_{\odot}}$ in order to correctly predict the properties of young massive star clusters, including those residing in high-redshift star-forming galaxies.

Since the MIST model grid ranges from $0.1$ to $300~{}{M_{\odot}}$, we are able to examine the contribution from VMSs to the SED, more specifically the ionizing luminosity over time. Recall from Figure 6 that the most massive stars contribute significantly to the flux blueward of 228 Å at 1 and 2 Myr, and completely dominate at 3 Myr. Beyond 3 Myr, they contribute zero flux because they have all but disappeared.

Figure 13 shows the ratio of ionizing luminosity emitted by stellar populations with IMF cutoffs at 300 and $100~{}{M_{\odot}}$. For simplicity, we show only the $v/v_{\rm crit}=0.0$ and 0.6 models. Similarly, Figure 14 shows the ratio of the total bolometric luminosity and momentum output. As expected, the ratio is nearly unity beyond $\sim 4~{}\rm Myr$, by which time stars more massive than $100~{}{M_{\odot}}$ have disappeared. In detail, it is slightly less than unity due to the IMF weights. Moreover, note that the ratio of the momentum injection from supernovae is unity at all times because we assume that only stars less massive than $40~{}{M_{\odot}}$ end their lives as SNII. These plots confirm what could be gleaned from Figure 6; the inclusion of VMSs can have an important effect on the integrated stellar population properties at very early times. Interestingly, this amounts to a factor of a few difference in some cases, e.g., hydrogen ionizing flux, bolometric luminosity, and radiative momentum output, whereas there can be an order of magnitude or more difference in the hardest EUV photon flux and wind momentum output. The dramatic difference in the helium ionizing photon luminosity can be understood by examining the contribution to the flux blueward of 228 Å in Figure 6; the large peak between 2 and 4 Myr is directly linked to the inclusion (or exclusion) of massive stars completely stripped of their hydrogen envelope.

These results are in broad agreement with those reported by Smith et al. (2016), who investigated the discrepancy in ages inferred from SED-fitting ($\sim 1$ Myr; Calzetti et al. 2015) and WR spectroscopic features ($3\textrm{--}5$ Myr; Turner et al. 2015) of two nuclear star clusters in NGC 5253. They concluded that VMSs can explain the observed WR-like spectroscopic features, bringing the two age estimates into agreement ($1\textrm{--}2$ Myr) without needing to invoke older ages or extreme rotators. We separately computed the ionizing photon luminosity of a $v/v_{\rm crit}=0.4$, $\rm[Fe/H]=-0.50$ (the metallicity of NGC 5253 is $35\%$ of the solar value assuming the Asplund et al. 2009 solar oxygen abundance; Monreal-Ibero et al. 2012) stellar population both with and without stars more massive than $100~{}{M_{\odot}}$. Following Smith et al. (2016), we assumed a total stellar mass of $3.3\times 10^{5}~{}{M_{\odot}}$ for the two clusters and compared the ionizing flux at 2 and 4 Myr to the observed value of $2.2\times 10^{52}~{}\rm s^{-1}$ for the central 5 pc. When we set the IMF upper limit at $100~{}{M_{\odot}}$, we obtain $1.7\times 10^{52}~{}\rm s^{-1}$ and $1.0\times 10^{52}~{}\rm s^{-1}$ at 2 and 4 Myr, respectively. Importantly, the maximum value of the ionizing flux attained by a population excluding VMSs is only $\sim 1.9\times 10^{52}~{}\rm s^{-1}$. We confirm the Smith et al. (2016) result that models that do not account for VMSs underpredict the ionizing flux. However, when we extend the upper mass limit to $300~{}{M_{\odot}}$, we obtain $4.0\times 10^{52}~{}\rm s^{-1}$ (maximum) and $0.9\times 10^{52}~{}\rm s^{-1}$ at 2 and 4 Myr. The conclusions do not change qualitatively when we instead consider a non-rotating stellar population. However, some rotational mixing resulting from a moderate rotation is necessary in order to match the observed factor of $\sim 3$ enhancement in the nebular nitrogen abundance reported by Smith et al. (2016). Both the excess and the enrichment timescale for nitrogen in the MIST rotating model are in agreement with the Köhler et al. (2015) evolutionary model, which was demonstrated in Smith et al. (2016) to be capable of producing the level of enrichment observed in the clusters.

The main conclusion here is that for most integrated quantities of interest, such as the hydrogen ionizing luminosity, increasing the IMF upper mass cutoff to $300~{}{M_{\odot}}$ results in a factor of a few increase at very early times ($t\lesssim 3$ Myr; see also Stanway et al. 2016). This difference may still have important observational consequences, as discussed in Smith et al. (2016). For other quantities, such as the production of the most extreme EUV photons, the difference can be as large as several orders of magnitude, though these quantities are more sensitive to the underlying models of VMSs, which still have large uncertainties.

4.5. The Effects of Metallicity

Although we focus on $\rm[Fe/H]=-1.0$ and 0.0 stellar populations in this work to enable a roughly equal comparison across different models,888As a reminder, we compare with SB99+Geneva models at $Z=0.002$ and $Z=0.014$ and BPASS models at $Z=0.001$, and $Z=0.020$. it is important to explore the effects of metallicity over a broader range. Metal-poor populations are particularly worth investigating since the effect of rotation on stellar population properties becomes more significant at lower metallicities. Moreover, as we will discuss in Section 5, high-redshift galaxies believed to be the principle sources of ionizing photons have stellar metallicities of $\log(Z/Z_{\odot})\lesssim-2$ at $z=6$ (e.g., Ma et al., 2016a).

Figure 15 illustrates the effect of metallicity on the ionizing photon production. In the top panel, we show the total number of hydrogen ionizing photons produced by 10 Myr. The different colors correspond to rotating stellar populations at different metallicities. As expected, the ionizing photon production is more efficient in low metallicity environments. For comparison, the equivalent point for the lowest metallicity BPASS binary model available ($Z=0.001$) is shown in black. We note that for the BPASS model, integrating out to 30 Myr instead of 10 Myr makes a $\sim 10\%$ difference in the total number of photons. The bottom panel shows the fraction of total hydrogen-ionizing photons emitted by 10 Myr for the same set of models. The decrease in ionizing photon production with time is more gradual at lower metallicities. We revisit these points in Section 5.

4.6. The Effects of a Rotation Distribution

Observations of massive stars in Galactic and local group clusters show a distribution of rotation rates (e.g., Penny, 1996; Huang et al., 2010; Ramírez-Agudelo et al., 2013), which may hold important clues for shedding light on their formation processes (e.g., Zinnecker & Yorke, 2007). A distribution of rotation rates is interesting in the context of this work for the simple reason that it would introduce a range of lifetimes, and therefore a range of MSTO masses at a fixed age. For example, by the time the dense GMC is sufficiently dispersed, allowing for ionizing photons to escape, the most massive slowly rotating stars will have already died, perhaps as SNe. However, the most massive rapidly rotating stars will still be present and producing copious quantities of ionizing photons.

For this work, recall that we computed MIST isochrones at five different rotation rates—$v/v_{\rm crit}=0.0,0.2,0.4,0.5$, and $0.6$—for each metallicity. From these, we create three composite populations by combining SSPs with different $v/v_{\rm crit}$ according to the probability distribution functions (PDFs) shown in Figure 16. “FLAT” represents a simple case of a flat distribution. “PL” is a simple power law that decreases with the rotation rate. “SGAUSS” approximates a skewed Gaussian distribution centered at $v/v_{\rm crit}=0.2$, which is supposed to mimic the observed distribution of O-type stars in the 30 Doradus region of the Large Magellanic Cloud (Ramírez-Agudelo et al., 2013). We simplify our approach by assuming that the rotation distribution is mass-independent over the mass range of interest ($M>20~{}{M_{\odot}}$). This is an acceptable assumption for two reasons. First, the rotation rate measurements of stars more massive than $M\sim 40~{}{M_{\odot}}$ are scarce—for instance, $<20\%$ of the Ramírez-Agudelo et al. 2013 sample was estimated to have masses above $40~{}{M_{\odot}}$—so their rotation rates, let alone their mass-dependence, are still relatively uncertain. Second, the observed rotation rates do not necessarily represent the rotation rates at birth: processes such as stellar winds, magnetic braking, and binary interactions modify the surface rotation rates, and the angular momentum is redistributed in the interior through mechanisms that are not fully understood in stars of all masses (e.g., Wolff et al., 2004; Cantiello & Langer, 2010; Gallet & Bouvier, 2013; Fuller et al., 2015). However, studying a select sample of young MS stars somewhat alleviates this issue (for example, Huang et al. 2010 uses a $\log g$ criterion). Further complications include macroturbulence (e.g., Penny & Gies, 2009) and $\sin i$ projection effects, which make the task of inferring the initial conditions and the subsequent evolution much more complex. We attempt to mitigate some of these problems by exploring a range of rotation PDFs. Nevertheless, given these complexities, we emphasize that this is an attempt to explore the broad effects of a rotation distribution on observed quantities.

In Figure 17, we plot the total number of hydrogen-ionizing photons emitted from time $t$ to 30 Myr for $\rm[Fe/H]=-1.0$ (left) and $0.0$ (right). The blue and orange curves correspond to $v/v_{\rm crit}=0.0$ and 0.6 single-rotation models and the gray dashed line represents the BPASS binary model, which assumes an IMF cutoff of $300~{}{M_{\odot}}$ for $Z=0.001$ (left) and $Z=0.014$ (right). In the left panel, we additionally show the $v/v_{\rm crit}=0.4$ model at $\rm[Fe/H]=-2.0$ for comparison. The gray vertical line marks roughly when GMCs begin to disperse, yielding order unity escape fraction of ionizing photons out of GMCs (e.g., Ma et al., 2015).

A composite population with a flat rotational distribution (“FLAT”)—which has the largest contribution from $v/v_{\rm crit}=0.6$—still has three times fewer photons in its budget compared to the fastest-rotating and binary scenarios by the time GMCs begin to disperse. The very low-metallicity case ($\rm[Fe/H]=-2.0$), where the rotation effects are expected to be large and therefore the fall-off in the production rate of ionizing photons is more gradual in time (Figure 15), has a comparable ionizing photon budget compared to the fast-rotating and binary scenarios by 3 Myr. However the number of available photons falls by an order of magnitude below that of the binary model by 8 Myr.

Interestingly, a very rapidly rotating population with $v/v_{\rm crit}=0.6$ is able to mimic the time-dependence of ionizing luminosity in a binary population. However, observations of single O-type stars in the LMC (Ramírez-Agudelo et al., 2013) effectively rule out populations of very fast-rotating massive stars, at least in the local Universe. It is also interesting to note that the ionizing luminosity in a moderately rotating stellar population at very low metallicity can rival that of the binary model and rapidly rotating model, though the production rate also falls off steeply beyond $\sim 5\textrm{--}6$ Myr.

5. Discussion

In this section we discuss the implications of the rotating, massive star models in two main contexts: cosmic reionization and the interpretation of high-redshift star-forming galaxies.

In our current understanding of the high-redshift universe, cosmic reionization was largely driven by high energy photons leaking out from star-forming galaxies (e.g., Madau et al., 1999; Haardt & Madau, 2012). The escape fraction of these ionizing photons, $f_{\rm esc}$, required to explain the observed ionization state of the $z=6$ intergalactic medium exceeds $10\%$ (e.g., Finkelstein et al., 2012; Robertson et al., 2013; Gnedin & Kaurov, 2014), well above the values inferred from their lower redshift analogs (e.g., Iwata et al., 2009; Boutsia et al., 2011; Siana et al., 2010; Leitet et al., 2013). Galaxy simulations have produced a wide range of predicted $f_{\rm esc}$ (e.g., Gnedin et al., 2008; Razoumov & Sommer-Larsen, 2010; Wise et al., 2014). In particular, a recent suite of simulations from the Feedback in Realistic Environment (FIRE) project produced a time-averaged $f_{\rm esc}$ of only $\lesssim 5\%$ (Ma et al., 2015). Since the production of ionizing photons is dominated by young massive stellar populations embedded in optically-thick natal environments, $f_{\rm esc}$ from these simulations depends sensitively on the complex connection between the properties of the underlying stellar population model999The FIRE simulations use SB99 with the “Padova+Asymptotic Giant Branch” option (Bressan et al., 1993; Fagotto et al., 1994a, b; Girardi et al., 2000). and the mechanisms that disrupt the GMCs (e.g., Murray et al., 2010). In their follow-up work, Ma et al. (2016b) found that replacing the underlying stellar population model with one that includes binary effects (BPASS) in their post-processing radiative transfer calculations can boost $f_{\rm esc}$ significantly due to its ability to sustain ionizing photon production at late times, well after the disruption of the GMCs. As emphasized by the authors, the key difference between the fiducial single star models and the binary models is the time-dependence of the ionizing photon production. Other mechanisms have also been proposed to increase the escape fraction of ionizing radiation, including runaway massive stars (Conroy & Kratter, 2012; Kimm & Cen, 2014).

In Section 4.6, we explored the same concept by considering low metallicities and velocity distributions as means to prolong the ionizing photon production. Together, Figures 15 and 17 demonstrate that stellar populations in sufficiently low-metallicity environments (e.g., high redshift) require only moderate rotation rates in order to produce significant ionizing photons. Moreover, at the time of order unity escape fraction ($\sim 3$ Myr), the ionizing photon budget can rival that of a binary population or a fast-rotating but slightly more metal-rich population. However, it is very challenging to p͡rolong the emission of ionizing photons by an amount required to reproduce the relatively large values of $f_{\rm esc}$ implied by the cosmic reionization models (Ma et al., 2016b), and we conclude that binary interactions may indeed play a critical role. Detailed numerical simulations using these low-metallicity rotating stellar population models are required to address this point more quantitatively.

Rotating massive stars have also been considered in the context of high-redshift star-forming galaxies. Recently, Steidel et al. (2014) have found that the locus of $z\sim 2\textrm{--}3$ star-forming galaxies is offset relative to their $z=0$ counterparts in the BPT diagram (Baldwin et al., 1981), an empirical diagnostic for probing the physical conditions and ionizing sources of nebular gas. The authors concluded that this shift could be explained by a harder stellar ionizing radiation field ($T_{\rm blackbody}=50000\textrm{--}60000$ K), higher ionization parameter (the ratio of the number densities of hydrogen ionizing photons to hydrogen atoms in a H ii region), and nitrogen-enhanced nebular gas, all of which could be explained by rapid rotation or binaries. In subsequent work, Steidel et al. (2016) found that the BPASS binary models are able to generate sufficiently hard ionizing radiation fields to reproduce the various observed line ratios, while the Geneva+SB99 cannot self-consistently explain the observations at any gas-phase metallicity or ionization parameter. Figure 7 shows that the FSPS+MIST models predict EUV flux that falls right in between the two models considered in Steidel et al. (2016). A more detailed comparison between the emission line predictions for the FSPS+MIST, Geneva+SB99, and BPASS models will be carried out in Byler et al. (2017, submitted).

6. Conclusions

In this work we explored the integrated properties of massive, rotating single-star stellar populations in both metal-rich and metal-poor environments. Our main conclusions are as follows:

We confirm that rotation leads to longer MS lifetimes and brighter and hotter stars. We also find that the magnitude of the effects is non-linear with the rotation rate. In particular, the $v/v_{\rm crit}=0.5$ and 0.6 models show a significant enhancement in the ionizing luminosity due to the blueward evolution of massive stars in the HR diagram at ages greater than 3 Myr. A comparison of the predicted MS lifetimes for the Geneva (Ekström et al., 2012; Georgy et al., 2013; Yusof et al., 2013) and MIST (Dotter, 2016; Choi et al., 2016) stellar evolutionary tracks demonstrates that the MS lifetime boost at a fixed mass is larger in the Geneva models compared to that in the MIST models, which suggests that rotational mixing may be more efficient in the former. 2. 2.

Rotation effects become more significant at lower metallicities as the star becomes more compact and angular momentum loss due to winds becomes less important. We found that the $\rm[Fe/H]=-2.0$ population produces twice as many photons as the $\rm[Fe/H]=0.0$ population by 10 Myr. From comparisons among models ranging from $\rm[Fe/H]=0.0$ to $-2.0$, we conclude that rotation leads to a more gradual decline in ionizing luminosity with time at low metallicities. As demonstrated in Ma et al. (2016b), this time-dependence has interesting ramifications for the escape fraction of ionizing photons, which in turn could impact high redshift, metal-poor galaxies at the time of cosmic reionization. 3. 3.

The contribution from VMSs ($\geq 100~{}{M_{\odot}}$) is significant for a short period of time ($t\lesssim 4$ Myr), increasing the ionizing luminosity by a factor of a few (H i) to a few orders of magnitude (He ii). Although the cumulative effect is small in most cases, taking VMSs into account in SPS modeling may be important for understanding a number of recent observations, including the ages and ionizing photon budget in nuclear star clusters (Smith et al., 2016). 4. 4.

We explored composite stellar populations with three different rotation PDFs to investigate whether a small number of fast rotators can lead to a significant boost in the ionizing photon luminosity at late times ($\sim 10$ Myr) compared to a non-rotating population. We find that the effect is small, amounting to a factor of two or smaller for a skewed distribution peaking at low $\rm v/v_{\rm crit}$ and a factor of $\sim 3$ for the extreme case of a flat PDF. Numerical simulations are required to assess whether or not this effect is important for the disruption of GMCs and/or the escape of ionizing radiation from galaxies.

We thank J.J. Eldridge, Mason Ng, and Georgie Taylor for sharing with us their unpublished WM-Basic models. We thank Matteo Cantiello, Aaron Dotter, and Daniel Weisz for useful comments and suggestions that have improved the manuscript. We also thank Bill Paxton and the rest of the MESA community who have made this project possible. C.C. acknowledges support from NASA grant NNX13AI46G, NSF grant AST- 1313280, and the Packard Foundation.