2.1.1. bpass single star models

The stellar evolution used in bpass is a derivative of the Cambridge stars code that was first developed by Eggleton (Reference Eggleton1971). It has since been updated by many authors but the most recent thorough description for the variant we employ was given in Eldridge et al. (Reference Eldridge, Izzard and Tout2008). It is a Henyey, Forbes, & Gould (Reference Henyey, Forbes and Gould1964) type evolution code that solves for the detailed stellar structure using an adaptive numerical mesh and time step to track the stellar evolution. For our single star evolution code, the models have initial masses ranging from 0.1 to 300 M_⊙. We calculate every decimal mass between 0.1 and 10 M_⊙, every integer mass between 10 and 100 M_⊙ and beyond that increase the stellar mass by 25 M_⊙ between models. We attempt to compute all our stellar models from the zero-age main sequence to the end of evolution in a single run. We first calculate the models with a resolution of 499 meshpoints. For a small number of models, numerical difficulties are encountered and we recalculate these model at a resolution of 199 meshpoints to overcome such problems.

We assume an initial uniform composition of the stars such that the initial compositions of the hydrogen and helium abundances by mass fraction are given by X = 0.75 − 2.5Z and Y = 0.25 + 1.5Z, respectively. Here, Z is the initial metallicity mass fraction which scales all elements from their Solar abundances as given in Grevesse & Noels (Reference Grevesse, Noels, Prantzos, Vangioni-Flam and Casse1993). This is selected to match the composition of the opacity tables included in the models. The opacity tables are described in Eldridge & Tout (Reference Eldridge and Tout2004a) and are based on the OPAL (Iglesias & Rogers Reference Iglesias and Rogers1996) and Ferguson et al. (Reference Ferguson, Alexander, Allard, Barman, Hauschildt, Heffner-Wong and Tamanai2005) opacities.

There is no strong consensus in the literature regarding the definition of Solar metallicity. Villante et al. (Reference Villante, Serenelli, Delahaye and Pinsonneault2014), for example, suggest the metal fraction in the Sun is rather higher than the Z = 0.020 usually assumed, while some authors (Allende Prieto, Lambert, & Asplund Reference Allende Prieto, Lambert and Asplund2002; Asplund Reference Asplund2005) suggest that Solar metal abundances should be revised downwards to closer to Z = 0.014 (also appropriate for massive stars within 500 pc of the Sun, Nieva & Przybilla Reference Nieva and Przybilla2012). We retain Z_⊙ = 0.020 for consistency with our empirical mass-loss rates which were originally scaled from this value. We note that at the lowest metallicities of our models, the small uncertainty in where we scale the mass-loss rates will cause similarly small changes in the mass-loss rates due to stellar winds. At the lowest metallicities, mass-loss is primarily driven by binary interactions. The metallicities we calculate have Z = 10⁻⁵, 10⁻⁴, 0.001, 0.002, 0.003, 0.004, 0.006, 0.008, 0.010, 0.014, 0.020, 0.030, and 0.040 (i.e. 0.05% Solar to twice Solar).

Convective mixing is modelled by standard mixing-length theory and the Schwarzchild criterion. We include convective overshooting with δ_OV = 0.12. This value for the amount of overshooting was calibrated against observed double-lined eclipsing binaries. The calibrations involved assuming the area of both components was the same and overshooting was varied until the surface parameters of the stars were reproduced (Pols et al. Reference Pols, Tout, Eggleton and Han1995; Schroder, Pols, & Eggleton Reference Schroder, Pols and Eggleton1997; Pols et al. Reference Pols, Tout, Schroder, Eggleton and Manners1997; Stancliffe et al. Reference Stancliffe, Fossati, Passy and Schneider2015).

We do not follow rotational mixing in detail but rather adopt a simple parameterisation in which we assume a star is either not mixed or fully mixed depending on binary mass transfer (see next section). If the mass transfer from a companion exceeds 5% of a star’s initial mass, it is assumed to lead to either rejuvenation and/or quasi-chemically homogeneous evolution (QHE). Rejuvenation is the consequence of rapid mixing of new material into the stellar core, resetting its evolution to the zero-age main sequence, thereafter the star evolves normally. QHE follows if the rotational mixing continues after mass transfer ends, disrupting the formation of shell burning and other internal processes. We have calculated simple QHE models by assuming certain stars are fully mixed throughout their main-sequence lifetimes. They therefore evolve to higher temperatures during this time, going ‘the wrong way’ on the Hertzsprung–Russell diagram. We assume this evolution is only possible at Z ⩽ 0.004 and for masses above 20 M_⊙ (Yoon, Langer, & Norman Reference Yoon, Langer and Norman2006).

The final detail of our single star models is the mass-loss scheme incorporated for stellar winds. There are two aspects of this scheme: first, the mass-loss rate to apply to a star during a certain phase of its evolution, and second, how the mass-loss rates scale with metallicity. In all our models, we apply the mass-loss rates of de Jager, Nieuwenhuijzen, & van der Hucht (Reference de Jager, Nieuwenhuijzen and van der Hucht1988), unless the star is an OB star where we use the mass loss rates of Vink, de Koter, & Lamers (Reference Vink, de Koter and Lamers2001). These are theoretical rates that do not include clumping but do match observed mass-loss rates (e.g. Shenar et al. Reference Shenar2015) and similar calculations including clumping do show that for the luminous O star the rates are unchanged (Muijres et al. Reference Muijres, Vink, de Koter, Müller and Langer2012). For Wolf–Rayet stars, when the surface hydrogen abundance is less than 40% and the surface temperature is above 10⁴ K, we use the mass-loss rates of Nugis & Lamers (Reference Nugis and Lamers2000). At different metallicities, we scale these mass-loss rates by $\dot{M}(Z)=\dot{M}(\text{Z}_{\odot })(Z/\text{Z}_{\odot })^{\alpha }$ and typically use α = 0.5, except in the case of OB stars where α = 0.69 (Vink et al. Reference Vink, de Koter and Lamers2001). There is some evidence that the value should also be greater for Wolf–Rayet stars (Hainich et al. Reference Hainich, Pasemann, Todt, Shenar, Sander and Hamann2015) but this is currently not considered.

All our models evolve from the zero-age main sequence up to the end of core carbon burning, or neon ignition for our most massive star models. Less massive models either form carbon–oxygen or oxygen–neon cores and evolve up the Asymptotic Giant Branch (AGB). We do not model AGB thermal pulses in detail nor incorporate specific AGB mass-loss rates at the current time. To do so is difficult with the stars code and requires care to overcome numerical difficulties (Stancliffe, Tout, & Pols Reference Stancliffe, Tout and Pols2004). At the current time, we limit the spatial and temporal evolution of our AGB models and so we do not observe thermal pulses. We find that the cores then grow up to the Chandrasekhar mass and fail when carbon is ignited in the CO core. This is unphysical and leads to overly luminous AGB stars in our populations. In future, we will recalculate these models with realistic AGB mass-loss rates to terminate the evolution earlier or add a routine with in the population synthesis to use rapid synthetic model of AGB evolution (e.g. Izzard et al. Reference Izzard, Tout, Karakas and Pols2004). The latter option would allow us to investigate the importance of AGB stars in more detail.

Our lowest mass models (below approximately 0.8 M_⊙) never ignite helium and end their evolution as helium white dwarfs. Models above this mass and up to approximately 2 M_⊙ fail at the onset of core-helium burning due to the core being degenerate at this time. Due to the rapid increase of luminosity in degenerate material with the stars code, we are unable to evolve through this event and the models fail. Therefore, to model the further evolution of these stars, we follow the method of Stancliffe, Izzard, & Tout (Reference Stancliffe, Izzard and Tout2005) and create more massive models where helium ignites in only a slightly degenerate core, we then decrease the mass of the star and allow the helium core to grow out until our model matches the parameters of the last failed helium-flash model. We combine these tracks to make the evolution track complete the possible evolutionary paths of these intermediate mass stars. This is a reasonable approximation as recent more detailed calculations indicate that the helium flash does not explode the star nor cause significant structural changes (Mocák et al. Reference Mocák, Müller, Weiss and Kifonidis2008).

2.1.2. bpass binary star models

Our binary models are identical to our single star models in nearly every respect but we additionally allow for extra mass loss or gain via binary interactions. We assume the orbits of the binary are circular and thus described by Kepler’s third law such that (M ₁ + M ₂)/M_⊙ ∝ (a/215R_⊙)³/(P/yr)², where a is the orbital separation and P is the orbital period. We assume mass lost in stellar winds removes orbital angular momentum in a spherically symmetric shell around the star losing the mass. Thus, the orbits widen over the evolution of the star. We only follow one star with detailed calculations during the evolution. This avoids wasting computational effort on calculating the evolution of a 1 M_⊙ secondary star at the same time as a more rapidly evolving 10 M_⊙ primary. During the evolution of the primary star models, we use the single star rapid evolution equations of Hurley et al. (Reference Hurley, Tout and Pols2002) to approximate the secondary’s evolution. When the end state of the primary’s evolution has been established, we then substitute the secondary’s evolution with a detailed model, either using a single star model or calculating a new detailed binary model with a compact remnant, depending on whether the binary is bound or unbound.

The grid of masses we use for primary evolution also ranges from 0.1 to 300 M_⊙. We model every decimal mass from 0.1 to 2.1 M_⊙, then 2.3, 2.5, 2.7, 3, 3.2, 3.5, 3.7 M_⊙, every half M_⊙ from 4 to 10 M_⊙, then every integer mass to 25 M_⊙, then every 5 to 40 M_⊙ then 50, 60, 70, 80, 100, 120, 150, 200, and 300 M_⊙. The grid of mass ratios, q = M ₂/M ₁, ranges uniformly from 0.1 to 0.9 in steps of 0.1. We establish initial binary periods with log (P/d) from 0 to 4 in steps of 0.2.

For the secondary models which have a compact companion, the grid of models we calculate is determined from the results of the population synthesis after the first supernova. The grid uses the same array of possible periods as for the primaries, however, the secondary masses are reduced to M ₂ = 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.8 M_⊙, every integer mass from 1 to 25 M_⊙, 25, 30, 35, 40, 50, 60, 70, 80, 100, 120, 150, 200, 300, 400, 500 M_⊙. The compact remnant masses are arranged in a grid of log (M _{rem, 1}/M_⊙) from −1 to 2 in steps of 0.1.

Our numerical method to account for binary evolution is inspired by the method of Hurley et al. (Reference Hurley, Tout and Pols2002) and was first described in detail in Eldridge et al. (Reference Eldridge, Izzard and Tout2008). However, we have had to modify some of the aspects as they are difficult to implement into a detailed stellar evolution code. The key point in our models when binary evolution becomes very different to that of our single star models is when the radius of the primary star increases and so it fills its Roche Lobe. This is the equipotential surface beyond which any material is no longer most strongly gravitationally attracted to the primary star, and Roche lobe overflow (RLOF) can occur. We allow for this in our models by using the effective Roche lobe radius defined by Eggleton (Reference Eggleton1983), where a star will have an equivalent volume to that of its Roche lobe. At radii beyond this, it will begin to overflow mass towards its companion star. The Roche lobe radius is given by

where q ₁ = M ₁/M ₂. When the star fills its Roche lobe, we determine a mass-loss rate due to overflow of

where

The value here was chosen by Hurley et al. (Reference Hurley, Tout and Pols2002) to ensure mass transfer is stable within their code. We also found this to be the case in applying their rate. There has been significant amount of work on the exact details of RLOF (e.g. Ritter Reference Ritter1988; Kolb & Ritter Reference Kolb and Ritter1990; Sepinsky, Willems, & Kalogera Reference Sepinsky, Willems and Kalogera2007) and our model is rather simple. However, upon the onset of RLOF, the stellar model will continue to expand until either CEE occurs or the mass transfer becomes stable. In future, we plan to calibrate the RLOF against observed systems undergoing mass transfer.

This mass is assumed to be transferred to the secondary star or lost from the system. We limit the accretion rate onto the secondary star by assuming it can only accept mass at a rate determined by its thermal timescale, such that $\dot{M_2} \le M_2/\tau _{\rm KH}$ . The rest of the material is assumed to be lost from the system, taking with it angular momentum from its orbit. In the case where the companion is a compact remnant with an initial mass less than 3 M_⊙, the accretion rate is limited instead by the Eddington luminosity. Above this mass, we assume that super-Eddington accretion can occur due to the compact remnant being a black hole, where super-Eddington luminosities have been observed. In the Eddington-limited case, any excess material and its angular momentum is lost from the system. These limits are also approximate and a more detailed accretion model will be required in future.

We do follow the rotation rate of the stars in our code but this has no direct physical consequence in the code beyond QHE due to spin-up in mass transfer in binaries (see Section 2.2.2). We follow the rotation rate to keep track of the angular momentum held by the stars, since as stars approach CEE, the Darwin mechanism can cause CEE to occur (Darwin Reference Darwin1880). While we attempted to implement tidal models such as those in Hurley et al. (Reference Hurley, Tout and Pols2002) based on Hut (Reference Hut1981), we found that the most stable method was to only assume tidal synchronisation once a star fills its Roche lobe. Assuming that the stars synchronise their rotation with the orbit, transferring angular momentum from the stars to the orbit or vice-versa. In future, we will add in a more general model of tidal forces and also allow the rotation to induce extra mixing in the stellar interior.

If RLOF does not stop the increase of the donor star’s radius, or if too much angular momentum is required to spin-up the donor star, then CEE occurs. In our code, CEE is taken to occur in our models when the radius of the primary star is the same as the binary separation. At this point, we switch the orbital evolution to a formalism based on the typical population synthesis model where the result of CEE is determined from comparing the binding energy of the envelope with the orbital energy (Ivanova et al. Reference Ivanova2013). We continue to determine the mass-loss rate from Equation (2) but set an upper mass-loss rate limit of 0.1 M_⊙ yr⁻¹ to avoid numerical problems of having such a high mass-loss rate. This unfortunately means we overestimate the CEE timescale in our models. To compare to the typical CEE mechanism, it is assumed that the change in the orbital energy of the binary is equal to the binding energy of the envelope of the star that is lost in the interaction:

This can be expressed as

where m ₁ and R ₁ are the mass and radius of the primary, a _i and a _f are the initial and final orbital separation, α_CE is a const representing the efficiency of the conversion from binding energy to orbital energy and λ is a constant representing how the structure of the star affects its binding energy (de Kool Reference de Kool1990; Dewi & Tauris Reference Dewi and Tauris2000).

In a detailed stellar model, it is difficult to remove the entire envelope in a single timestep. Therefore, we instead allow the mass-loss rate to increase as high as possible and relate the binding energy of the material lost to the change in the orbital energy. We do this by equating the binding energy of the lost material to the change in the orbital energy such that

which can be rearranged to give

where δa is the change in orbital separation and δM ₁ is the mass loss of the surface. The advantage of this method is that the structure of the star is naturally taken account of as we integrate over the removal of the envelope. The disadvantage is that because our time for CEE is a few orders of magnitude too large, other sources of energy may limit how efficient CEE is. However, the efficiency of CEE is uncertain.

Our CEE mechanism is meant to reproduce the envelope ejection and the in-spiral simultaneously. Eventually, we find that CEE either results in removal of the stellar envelope or a merger. If the primary star eventually shrinks within the Roche lobe, the CEE is taken to have ended. Mergers are assumed to occur when the companion star begins to fill its own Roche lobe. At this point, the total mass of the secondary is added to the primary star. We typically find the stars that experience CEE on, or just after, the main sequence tend to merge, while post-main-sequence stars tend to only remove the hydrogen envelope. Also, even in systems that do merge, mass loss prior to the event of the merger means that the final stellar mass is less than the total pre-CEE mass of the binary.

In Figure 1, we show for our fiducial population (as described in Table 1) the fraction of stars that interact, splitting this into those that experience RLOF or RLOF and CEE. We see that at masses of about 5 < M < 40 M_⊙, 80% of binaries interact with around 90% of those interactions leading to CEE. We note that for smaller mass ratios, 100% of interactions lead to CEE, while for our largest mass ratios, only 80% lead to CEE. Outside this, mass range RLOF dominates the binary interactions. At high masses, the stellar wind mass loss is already significant and so once RLOF starts, only a little extra mass loss is required to avoid CEE, especially at higher metallicity where stellar winds reduce the number of systems that interact to 30%. Below 5 M_⊙ down to approximately 1M_⊙, the fraction of interacting systems decreases. This is due to the maximum radius of the stars decreasing so fewer stars fill their Roche lobe and interact. At the lowest masses, below 0.7 M_⊙, the interaction fraction decreases to zero due to the stars being smaller than the orbital separation for their evolution within the age of the Universe. This is likely due to our minimum period for our binary models of 1 d, as well simple binary model that does not include processes that would allow for periods to shorten and drive interactions. As a result, with the exception of those in the closest binaries, such sources do not overfill their Roche lobes and so never interact. Stars below ~1 M_⊙ never interact in our default model set as they remain compact during their main sequence lifetime, which is long compared to the age of the Universe.

Figure 1. The fraction of primary stars in our fiducial population that experience a binary interaction at two metallicities within the age of the Universe versus the initial mass of the primary star. The solid line represents the total fraction that experience Roche lobe overflow (RLOF). While the dashed line represents the number when RLOF progresses to common envelope evolution (CEE).

We show the mean duration over which RLOF and CEE episodes occur in our binary models in Figure 2. More massive stars interact on shorter timescales than low mass stars, with RLOF ranging from 10⁵ up to 10⁷ yrs, while CEE has much shorter timescales of 100 to 10⁴ yrs. These times for our CEE events are much longer compared to the very short timescale of days or years that are predicted by dynamical simulations (see Ivanova et al. Reference Ivanova2013) as well as those found from observed events. However, they are significantly shorter than the thermal and nuclear timescales of the stars so they should have little impact on the eventual predictions of binary evolution.

Figure 2. The mean times that binary interactions last for in our fiducial simulation versus the initial mass of the primary star. The solid lines are for RLOF and the dashed lines for CEE. The thick lines are for the mean, the thin lines are at ±1σ. As discussed below, the CEE time are significantly overestimated due to our method of including CEE in our detailed evolution models.

2.1.3. Supernovae and remnants

At the end of the evolution of our models, we check to see what remnant may be produced. Currently, we assume that either a white dwarf, neutron star, black hole, or no remnant will be formed. We assume that a core-collapse supernova will occur if central carbon burning has taken place and the CO core mass is greater than 1.38 M_⊙ and the total stellar mass is greater than 1.5 M_⊙. If this condition is not met, then the remnant left will be a white dwarf with the mass of the helium core at the end of the stellar model, the secondary star will continue to evolve with a white dwarf in orbit around it.

In the case that a supernova occurs the remnant mass is determined by calculating how much material can be ejected from the star given an energy input of 10⁵¹ ergs, a typical supernova energy. This method was first outlined in Eldridge & Tout (Reference Eldridge and Tout2004b) which provides more details. The mass that is not unbound goes into the remnant, thus providing a way to estimate the mass of any resultant black hole. This method is similar to others (e.g. Heger et al. Reference Heger, Fryer, Woosley, Langer and Hartmann2003; Ugliano et al. Reference Ugliano, Janka, Marek and Arcones2012; Sukhbold et al. Reference Sukhbold, Ertl, Woosley, Brown and Janka2015; Spera, Mapelli, & Bressan Reference Spera, Mapelli and Bressan2015). For simplicity, neutron star remnants are assumed to have a constant mass of 1.4 M_⊙, while remnants more massive than 3 M_⊙ are considered to form black holes. Finally, if the helium core mass at the end of evolution is between 64 and 133 M_⊙, a pair-instability supernova (PISN) is assumed to occur which completely disrupts the star and leaves no remnant (Heger & Woosley Reference Heger and Woosley2002).

For any system with a core-collapse supernova, we estimate the supernova category as type IIP, II-other, Ib, and Ic. We identify the supernova types as described in Eldridge et al. (Reference Eldridge, Langer and Tout2011) and Eldridge et al. (Reference Eldridge, Fraser, Smartt, Maund and Crockett2013). If the mass of hydrogen in the star is greater than 10⁻³ M_⊙, then we have a type II supernova. If the total mass of hydrogen is greater than 1.5 M_⊙ and the ratio of hydrogen mass to helium mass in the progenitor is greater than 1.05, we have a type IIP. When the hydrogen mass is less than 10⁻³ M_⊙, a type Ib/c supernova occurs. If more than 10% of the ejecta mass is helium, then the supernova is assumed to be a type Ib, otherwise it is Ic. This leads to roughly twice as many Ib as Ic supernovae which is in agreement with the most recent attempt to accurately determine this ratio (Shivvers et al. Reference Shivvers2017).

We also determine whether an event is likely to generate a Gamma-ray Burst (GRB) or PISN. For a long-GRB, the star must have experienced QHE and have a remnant mass greater than 3 M_⊙. Note that this is only one possible pathway for GRBs and limits their presence in our rate estimates to lower metallicities. Some type Ic supernovae arising from non-QHE systems also likely generate GRBs but we do not account for these. For PISN, we apply mass limits from Heger & Woosley (Reference Heger and Woosley2002) as detailed above. Finally, any progenitor with a final mass between 1.5 and 2 M_⊙ is identified as a low-mass progenitor with an uncertain outcome. These events may be faint and rapidly evolving, as described by Moriya & Eldridge (Reference Moriya and Eldridge2016).

We also estimate the type Ia thermonuclear supernova rate by two channels, first selecting out white dwarfs that begin with a mass below the Chandrasekhar limit but then reach or exceed this mass due to binary mass transfer. Second, we count double white-dwarf binaries that merge due to gravitational radiation with the merger time calculated assuming circular orbits and the analytic expressions of Peters (Reference Peters1964). Our resultant rates are highly approximate and we do not consider them rigorous or precise but include them for completeness. We caution that these should be used with care, and in future we may refine our type Ia models further for comparison with other predictions and observational data.

2.2.1. Binary distribution parameters

In the population synthesis, the individual stellar models are combined together to make a synthetic stellar population which can be observed and compared to observations of the real Universe. Together with the models described above, we need to know the distribution of initial parameters required to create the population. The three main parameters are the initial mass function (IMF), the initial period distribution, and the initial mass-ratio distribution. Each population is generated at a single initial metallicity. For our fiducial models, we base our IMF on Kroupa, Tout, & Gilmore (Reference Kroupa, Tout and Gilmore1993) with a power-law slope from 0.1 to 0.5 M_⊙ of −1.30 which increases to −2.35 above this. The power law slope extends to a maximum initial stellar mass of 300 M_⊙, although we note that due to mergers stars, more massive than this can form in our binary populations. Our most massive stars have masses up to 570 M_⊙, although these are very rare systems. In addition to this standard model, we also calculate models with two different upper IMF slopes of −2.00 and −2.70, as well as varying the upper maximum initial stellar mass to 100 M_⊙. Finally, we calculate a model set with a constant IMF slope of −2.35 from 0.1 to 100 M_⊙ [the standard (Salpeter Reference Salpeter1955) IMF], for comparison with previous work.

For the binary populations, we assume a flat distribution in initial mass ratio and log period in all our models. The result of this is that while all of our stars are technically part of a binary population, in many cases the stars evolve as isolated individuals and are never close enough to interact. The actual interacting binary fraction will depend on the physical size of stars at different ages relative to their Roche lobes, and thus on mass and metallicity. Our upper period limit of 10 000 days has been chosen so that about 80% of massive stars (M ≳ 5 M_⊙) interact, as shown in Figure 1. From 5 to 0.7 M_⊙, this fraction decreases to about 50%. At lower masses, the separation of the orbits with a 1 d minimum is too wide for stars to interact within the age of the Universe. At Solar metallicity, many of the most massive (M ~ 100 M_⊙) stars also avoid interactions due to strong mass loss in stellar winds.

Observational constraints on these numbers are uncertain but our estimates of massive star interaction are comparable to the 70% estimated for O stars by Sana et al. (Reference Sana2012). On the other hand, the flat distributions are probably over-simplistic. As discussed by De Marco & Izzard (Reference De Marco and Izzard2017), binary parameters (including the binary fraction) vary with initial mass. While nearly every massive star is in a binary, the binary fraction drops to 20% at lower initial masses below a Solar mass, but the precise binary fraction as a function of stellar mass remains poorly constrained in the literature. In the Milky Way, Sana et al. (Reference Sana2012) found that the observed period distribution is slightly steeper than flat, with a bias towards more close binary systems for O stars. However, Kiminki & Kobulnicky (Reference Kiminki and Kobulnicky2012) found a flatter period distribution in the Cygnus OB2 association that is consistent with Öpik (Reference Öpik1924)’s law, although their results did also suggest a slight preference for short period systems as well. The most recent study by Moe & Di Stefano (Reference Moe and Di Stefano2017) has quantified the period and mass ratio distribution in an extensive compilation of previous works, and finds a slightly different period and mass ratio distribution to that we use. We intend to modify our input distributions accordingly in a future version of BPASS.

The uncertainty in assumed period distribution is degenerate with uncertainties in the assumed model to handle RLOF, common envelope evolution, tides, and other binary specific processes. Hence, we adopt this simple approach based on the well-studied massive star population (Sana et al. Reference Sana2012). If required, it is possible to simulate the effect of more complex populations and binary distributions by varying the mix between our single star and binary populations with age (since stars of different initial masses dominate in different time bins).

We note that, because we assume orbits are circular, our closest observational comparison will be with the semi-latus rectum distribution not the period distribution. As shown by Hurley et al. (Reference Hurley, Tout and Pols2002), the outcome of the interactions of systems with the same semi-latus rectum is almost independent of eccentricity. This is equivalent to assuming that systems are circularised by tidal forces before interactions occur. We only include tides in our evolution models when a star fills its Roche lobe as this can move the system into CEE if there is not enough angular momentum in the orbit to spin up the primary. We assume tidal forces are strong and the star’s rotation quickly synchronises with the orbit. This is of course an approximation and there are recent studies have begun to explore how mass transfer may be different in eccentric systems (e.g. Sepinsky et al. Reference Sepinsky, Willems and Kalogera2007; Bobrick, Davies, & Church Reference Bobrick, Davies and Church2017; Dosopoulou & Kalogera Reference Dosopoulou and Kalogera2016a, Reference Dosopoulou and Kalogera2016b). However, we note that the bpass stellar models have been successfully tested to see if they can reproduce an observed binary system with a slight eccentricity even after mass transfer (Eldridge Reference Eldridge2009).

2.2.2. Binary evolution treatment

We combine all the primary star models together according to these distributions. We then must consider how to account for the secondary stars. We pick out the stars that explode in supernovae by and estimate the remnant mass as described in Section 2.1.3. For remnants less than 3 M_⊙, we assume that the remnant receives a kick taken at random from a Maxwell–Boltzmann distribution with σ = 265 km s⁻¹ (see Hobbs et al. Reference Hobbs, Lorimer, Lyne and Kramer2005). For remnants more massive than 3 M_⊙, we assume the star has a kick that is reduced by a factor of the remnant mass divided by 1.4 M_⊙, assuming that the kicks track a momentum distribution rather than a velocity distribution. There is some observational evidence that black-hole kicks should be weaker than for neutron stars (Mandel Reference Mandel2016). We note that we have also investigated other models of neutron-star kicks but have not yet included these in the released bpass models (see Bray & Eldridge Reference Bray and Eldridge2016).

If the binary is disrupted, then the companion star is modelled as a single star in its further evolution. Conversely, if it remains bound, then the post-supernova orbit is calculated, circularised and a secondary model of a star orbiting a compact remnant is used to represent its future evolution. Eventually, if a second supernova occurs, the same process is repeated so that the parameters of double-compact remnant binaries can be estimated. Note that the secondary may still undergo supernova, even if the primary does not: this possibility is followed and these later events are also considered.

We also check in our models for rejuvenation of the secondary stars. If a secondary star accretes more than 5% of its initial mass and it is more massive than 2 M_⊙, then the star is assumed to be spun-up to critical rotation and is rejuvenated due to strong rotational mixing. That is, it evolves from the time of mass transfer as a more massive, zero-age main-sequence star (this is similar to the method used by Vanbeveren, De Loore, & Van Rensbergen Reference Vanbeveren, De Loore and Van Rensbergen1998). At high metallicities, we assume that the star quickly spins down by losing angular momentum in its wind so that there are no further consequences for its evolution. However, with weaker winds at lower metallicity, Z ⩽ 0.004, we assume that the spin down occurs less rapidly and so the star remains fully mixed throughout the main sequence and burns all its hydrogen to helium. We assume the star must have an effective initial mass after accretion >20 M_⊙ for this evolution for occur, based on the work of Yoon et al. (Reference Yoon, Langer and Norman2006). We note that we have discussed the importance of these quasi-homogeneously evolving (QHE) stars in greater detail in Eldridge et al. (Reference Eldridge, Langer and Tout2011), Eldridge & Stanway (Reference Eldridge and Stanway2012), and Stanway et al. (Reference Stanway, Eldridge and Becker2016) showing there is strong observational evidence that they exist. We only consider QHE that arises from mass transfer, however, and do not yet include it in tidal interactions as done by Mandel & de Mink (Reference Mandel and de Mink2016) and Marchant et al. (Reference Marchant, Langer, Podsiadlowski, Tauris and Moriya2016). In such a scenario, stars can be spun up by tidal forces acting between two stars in orbital periods of the order of a day, leading to both stars experiencing QHE, this is not yet implemented in bpass.

2.2.3. Binary mass function as a function of age

We present some fundamental results from our populations in Figures 3 and 4 which concern the mass in the stars at different different stellar population ages. Figure 3 shows that from an initial 10⁶ M_⊙, the total mass of stars in a co-eval binary population slowly decreases with time as stars die and lose mass in stellar winds. We see that the single star populations tend to decrease faster than the binary populations initially, while the binary populations retain more of their mass until the first SN occurs. This leads to an offset between the binary and single star populations which persists until late times. Higher metallicity populations retain more mass at late times. This is because the higher metallicity stars have longer lifetimes due to being less compact with lower central temperatures so burning through their nuclear fuel takes a longer time. Although at early times, the far stronger winds of high metallicity stars means more mass is lost at early times for the higher metallicity population. After this, the stars will evolve along the lower mass pathways at an appropriate rate. The net effect of this is to lead to longer total stellar lifetimes and so a delay in the death of the stars and formation of remnants. By the age of 1 Gyr, a Solar metallicity starburst will have lost 37% (39%) of its initial mass in the binary (single) star evolution case. This is towards the upper end of the ‘return fractions’ of material to the ISM of R ~ 0.4 for a Chabrier (Reference Chabrier2003) IMF and R ~ 0.3 for a Salpeter (Reference Salpeter1955) IMF found by Madau & Dickinson (Reference Madau and Dickinson2014). At 20% of Solar metallicity, both binary and single star return fractions increase by ~2%.

Figure 3. The evolution of the fiducial stellar population mass (the mass contained in stars only) with time. We assume formation of an initial population with total mass 10⁶ M_⊙ within the first Myr. Solid lines are for binary populations and dashed lines are for the single-star populations, and results are shown at two metallicities.

Figure 4. Contour plots showing how the stellar mass function evolves with age. The greyscale contours represent a binary population with each contour being a change in number by an order of magnitude. The lines represent the maximum mass at any given time, the solid line for the initial mass of the star, and the dashed line the final mass of the star with that lifetime.

Figure 4 illustrates this effect in greater detail, showing how the stellar mass functions vary with age at the same two metallicities. In both cases, due to rapid mergers of the most massive members of a stellar population, binary populations can have more massive stars than were initially assumed and these can retain mass for longer. At late times, the most massive star is always beyond that expected from single-star populations: the phenomenon that leads to the observation of ‘blue stragglers’ in globular cluster populations. The mass of the most massive star is greater at late times in high metallicity populations, even though at early times the mass of individual stars rapidly decreases due to the stronger mass loss.

We note here that beyond approximately 3 Gyrs, our current models predict a population of more massive stars than would be expected. These are more than twice the maximum mass for a single star population at this age. These are systems where the primary evolved to a white dwarf after rejuvenating its companion star. Our simple method of estimating the rejuvenation age currently uses the age at the end of the model rather than the time of the mass transfer. While this is adequate if the primary explodes in a supernova, it is not for the case where our models include the time on the white-dwarf cooling track.

2.3.1. Existing model grids

Our base atmosphere models are those of the BaSeL v3.1 library (Westera et al. Reference Westera, Lejeune, Buser, Cuisinier and Bruzual2002). This is a set of atmosphere models that extend over a large range of log g, log T _eff, and metallicity. We interpolate between models to predict the spectrum of each star generated by our evolution code. This version only extends over metallicities from [Fe/H] = −2.0 to 0.5. We therefore at our lowest metallcities use the slightly older v2.2 grid of Lejeune, Cuisinier, & Buser (Reference Lejeune, Cuisinier and Buser1998) that extends down to [Fe/H] = −3.5.

We supplement these with Wolf–Rayet (WR) stellar atmosphere models from the Potsdam PoWR groupFootnote ³ (Hamann & Gräfener Reference Hamann and Gräfener2003; Sander et al. Reference Sander, Shenar, Hainich, Gímenez-García, Todt and Hamann2015). We only use WR spectra once the surface hydrogen abundance drops below 40% and the surface temperature is above log (T/K) = 4.45. A key change implemented in bpassv2.0 and v2.1 is that there are now different metallicity nitrogen-dominated (WN) WR atmospheres available (Todt et al. Reference Todt, Sander, Hainich, Hamann, Quade and Shenar2015). In our previous work (such as Eldridge & Stanway Reference Eldridge and Stanway2009), only Solar metallicity were available and we artificially reduced the line luminosities of these models at lower metallicites. We now use the closest metallicity WN models to the stellar models being considered. When Z ⩾ 0.0095, we use the Solar PoWR models, below this limit we use their LMC models and when Z < 0.0035, we use their SMC models. For WC (i.e. carbon-dominated) WR stars, there are still only Solar metallicity models available, but we have not altered these spectra in any way in our current models, using the Solar models at all metallicities. This may lead to slight overestimates in the strength of carbon WC features in low metallicity spectra but the abundances in the stellar atmospheres during the WR phase are relatively insensitive to the initial metallicity (which has a larger effect on mass-loss rate). This will be modified as further PoWR models become available. We switch to the WC models when there is no surface hydrogen and the helium mass fraction drops below 70%.

As in Eldridge & Stanway (Reference Eldridge and Stanway2009), we add in an excess He ii flux for massive of stars as suggested by Brinchmann, Pettini, & Charlot (Reference Brinchmann, Pettini and Charlot2008), since these are under-represented in atmosphere models. For non-WR stars with a surface temperature above 33 000 K and when log g ⩽ 3.676 log (T/K) − 13.253, we increase the He ii stellar wind lines at 1 640 and 4 686 Å by 520 and 33.8 L_⊙, respectively. Due to the rarity of these cases, we find this has a minimal effect on the predicted population line fluxes. Interestingly, Crowther et al. (Reference Crowther, Barnard, Carpano, Clark, Dhillon and Pollock2010) have studied He ii emission from massive stars in the Tarantula Nebula and report a very high equivalent width for a young stellar population. We note that with our current prescription as described above, our v2.1 models produce a comparably high equivalent width at slightly older (but similar ages), see Section 6.1 and Figure 28 below.

2.3.2. New O star models

In previous distributions of bpass models, we supplemented the above with a set of O star models generated using the wm basic code by Smith, Norris, & Crowther (Reference Smith, Norris and Crowther2002). We have now recalculated this grid, and extended it, as summarised in Table 2 Footnote ⁴ .

For bpass v2.1 models, we have created new wm-basic (Pauldrach et al. Reference Pauldrach, Lennon, Hoffmann, Sellmaier, Kudritzki, Puls and Howarth1998) models based upon the grid determined by Smith et al. (Reference Smith, Norris and Crowther2002). The key differences are that we have extended the metallicity coverage and also added a new sequence of models with a higher surface gravity to better match the evolution code predictions. We have kept the range of temperatures the same (25–50 kK) but now added a grid of high gravity stars with log g = 4.5 (c.f. 4.0 in the earlier grid). This is because we noticed that on the zero-age main sequence, many of our O stars had gravities higher than 4.0 which has a significant effect on the strength of some stellar lines. At low metallicities, our QHE models also increase their surface gravity as they evolve with the stars shrinking as helium is mixed to the surface and the mean-molecular weight increases. Therefore, these new models are important for predicting the correct spectrum.

We also calculate models at every metallicity we consider in bpass rather than simply interpolating between fixed points as before. The key benefit is that we now have spectra at the lowest metallicities when Z = 10⁻⁵ and 10⁻⁴, where the influence of these massive stars on the synthetic spectrum is strongest. These are beyond the range for which observational validation data is available for stellar atmospheres and so constitute a theoretical extrapolation from better constrained atmospheres at higher metallicity. However, incorporating these gives greater predicting power for stellar populations at these extreme metallicities.

We have made these models available on the bpass website and they have already been used in other published works (e.g. Byler et al. Reference Byler, Dalcanton, Conroy and Johnson2017). We note that only the inclusion of the higher gravity models leads to any significant differences in model outcomes. We find that the ionising flux predictions at low metallicites decrease by about 0.1 dex, but that the far-ultraviolet luminosity of a population increase by a similar amount (see Figure 35 later). This is because the QHE evolution models have more correct atmospheres. These stars evolve to hotter temperatures as they burn hydrogen to helium. These models also allow us to predict O V 1 039 Å lines at very young ages that are comparable to those found in observations (A. Runnholm, private communication).

2.4.1. Spectral energy distributions

The basic output of bpass is a set of moderately high resolution model spectra, gridded from 1 to 100 000 Å in 1 Å bins and measuring total flux in each wavelength bin in units of L_⊙ (i.e. producing a spectrum of F _λ in L_⊙/Å spanning from the extreme-ultraviolet to far-infrared wavelengths). These are generated for stellar population ages of 1 Myr to 100 GyrFootnote ⁵ in increments of Δ log (age/yr) = 0.1 where this age is taken relative to the onset of star formation, assuming a co-eval stellar population of total mass 10⁶ M_⊙. Examples are shown in Figures 5 and 6, illustrating the optical-infrared spectral energy distributions as a function of age and the extreme ultraviolet as a function of metallicity and binary evolution, respectively.

Figure 5. The synthetic spectra produced for a co-eval population (i.e. instantaneous starburst) at times of 1, 3, 10, 30, 100, 300, 1 000, and 3 000 Myr after star formation. Spectra are shown for binary populations (bold, coloured lines) and single stars (pale, greyscale line), and at metallicities of Z = 0.002 (top) and Z = 0.020 (centre). In the bottom panel, we compare bpass binary models at the two metallicities directly, at ages of 3, 30, and 300 Myr.

Figure 6. The extreme ultraviolet (200–1 600 Å) spectral region in synthetic spectra produced for a co-eval population (i.e. instantaneous starburst) at a time 30 Myr after star formation. The two panels show the spectra as a function of metallicity (0.05, 0.5% Z_⊙ in the left-hand panel, 10 and 100% Z_⊙ in the right-hand panel) and bpass single star vs. binary evolution for each metallicity. The effect of binary evolution is to increase the hardness of the spectrum, and hence the ionising photon output (total flux emerging short of 912 Å), particularly at very low metallicities. Synthetic spectra have been scaled to a common luminosity at 1 500 Å.

In Figure 7, we show a direct comparison between the simple stellar population models (i.e. co-eval starbursts) produced by different stellar population synthesis codes at Solar metallicity and allowed to evolve over time. In three panels, bpass single star models are compared to the closest available matching models in metallicity, age, and IMF from publicly available synthesis codes. Perhaps, the most direct comparison is between bpass and the starburst99 models (Leitherer et al. Reference Leitherer1999; Leitherer et al. Reference Leitherer, Ekström, Meynet, Schaerer, Agienko and Levesque2014, hereafter S99), both of which were originally designed and optimised for young starbursts. The S99 models shown here were generated using the non-rotating Geneva 2012/2013 stellar model set and a broken power law IMF with slopes matching our default IMF but with a maximum stellar mass of 100 M_⊙.

Figure 7. A comparison between the simple stellar population models (i.e. instantaneous starbursts) produced by different stellar population synthesis codes at Solar metallicity. In three panels, bpass single star models are compared to the closest available match in metallicity, age, and initial mass function from publicly available synthesis codes. These are (top right) the galaxev models of Bruzual & Charlot (Reference Bruzual and Charlot2003), specifically the galaxy templates used to classify galaxies in the Sloan Digital Sky Survey by Tremonti et al. (Reference Tremonti2004), which use a Chabrier IMF; (bottom right) the Maraston (Reference Maraston2005) models, shown here using the Kroupa IMF; (top left) starburst99 models (Leitherer et al), generated with the non-rotating Geneva model set and an IMF matching our power law slopes, and (bottom left) bpass binary star models. Synthetic spectra have been scaled to a common luminosity at 5 000 Å.

We also compare against two SPS codes designed primarily for modelling galaxies. The galaxev models of Bruzual & Charlot (Reference Bruzual and Charlot2003, hereafter BC03) are widely used for modelling mature galaxies and specifically we compare against the galaxy templates used to classify galaxies in the Sloan Digital Sky Survey by Tremonti et al. (Reference Tremonti2004). These use a Chabrier (Reference Chabrier2003) IMF and are built on the Geneva and Padova single star stellar evolution models, combined with a large library of semi-empirical template spectra and additional atmosphere models from the BaSeL atmospheres (Westera et al. Reference Westera, Lejeune, Buser, Cuisinier and Bruzual2002). The Maraston (Reference Maraston2005, hereafter M05) models were designed to confront the thermally pulsating AGB stage in more detail than previous models and use a prescription for fuel consumption to determine the contribution of post-main-sequence stars and their role in determining the near-infrared colours of galaxies. They employ the Cassisi stellar evolution tracks (Cassisi & Salaris Reference Cassisi and Salaris1997; Cassisi et al. Reference Cassisi, Castellani, Ciarcelluti, Piotto and Zoccali2000), BaSeL atmospheres (Lejeune et al. Reference Lejeune, Cuisinier and Buser1998), and a Kroupa (Reference Kroupa2001) broken power law IMF similar to ours (with an upper power law slope of −2.30 where we use −2.35, and a maximum mass of 100 M_⊙).

As the figure demonstrates, the bpass single star models are very similar to all three other synthesis models at ages of <1 Gyr and in the ultraviolet, where the models are dominated by massive stars. The BC03 models are slight outliers, appearing redder at very young ages (<10 Myr) than any of the other model sets. Our single star models remain similar to the output of other codes at a matching age through the blue optical bands, with subtle differences in the IMF and evolutionary timescales assumed for the most massive stars by different evolution codes leading to our single star models being slightly bluer than S99 and M05 models but redder than BC03 models at ages of ~10–100 Myr.

bpass models begin to diverge significantly from the other model sets in the red-optical and longwards at stellar population ages of more than a few times 100 Myr. At these ages, the light from our single star populations are dominated by evolved massive stars, primarily AGB stars, whose molecular line features are clearly visible in the composite SED. Perhaps unsurprisingly, our models are most similar in the infrared to the M05 models, which are also heavily influenced by the AGB phase, but we see a stronger effect still. Realistically, this comparison suggests that AGB stars are over-represented in our co-eval single star populations (either in terms of number or luminosity) at ages >1 Gyr. In our favoured binary models (shown in the fourth panel of Figure 7), the number of AGB stars in a population is suppressed by the effect of binary interactions, but these interactions also keep the spectra rather bluer than those seen in other stellar population synthesis models at late times. We note that extreme caution should be applied when using this version of bpass in the red at ages >1 Gyr as a result, we return to this point in Section 7.10.

2.4.2. Released data products

The spectral energy distributions described above for co-eval simple stellar populations form the core of the bpass v2.1 data release. Composite (i.e. non-co-eval) stellar populations are not provided as a part of our standard release but can be constructed by assuming a star formation history (see Section 6.2). The simplest case is a population forming stars continuously at a constant rate. Here, the only potential difficulty lies in dealing with logarithmically spaced time intervals in the models since the total number of stars to be added to the composite spectrum depends on the width of the interval. Stars contributing to the first time bin (at log(age/yrs) = 6.0) include all members of the population up to age 10^6.05, the second bin stars in the age range 10^6.05 to 10^6.15 yrs and so forth.

We calculate magnitudes for individual stellar models and for our synthetic populations in the standard UBVRIJHK filters, as well as ugriz SDSS filters, FUV(1 566Å), NUV(2 267Å) fluxes, and several HST optical filters at each step. Additionally, we calculate the production rate of ionising photons based on the extreme ultraviolet flux since this is strongly affected by both binary evolution and metallicity [see Figure 6 and Stanway et al. (Reference Stanway, Eldridge and Becker2016)].

Given that our population synthesis tracks all stages of stellar evolution, we are also able to produce corresponding rates of supernovae (separated by type) and information on the compact remnant population and stellar class distribution at each time step.

1. 1. Stellar model outputs:

   1. a. Binary stellar models with photometric colours.

   2. b. New OB stars atmosphere models.

2. 2. Stellar population outputs (all versus age):

   1. a. Massive star type number counts.

   2. b. Core collapse supernova rates.

   3. c. Yields, energy output from winds and supernovae, and ejected yields of X, Y, and Z.

   4. d. Stellar population mass remaining.

   5. e. HR diagram (isochronal contours).

3. 3. Spectral synthesis outputs (all versus age):

   1. a. Spectral energy distributions.

   2. b. Ionising flux predictions.

   3. c. Broadband colours.

   4. d. Colour–magnitude Diagram (CMD) making code.

4. 4. Available on request due to unverified status:

   1. a. Approximate accretion luminosities from X-ray binaries.

   2. b. A limited set of nebula emission models.

All outputs of the current bpass v2.1 data release can be found at http://bpass.auckland.ac.nz and are mirrored at http://warwick.ac.uk/bpass. The results can also be found at the PASA datastore.

3.2.1. Red and blue supergiants

In Figure 11, we compare population density contour plots and 1 000-d period stellar tracks to observed RSGs in the SMC, LMC, M31, and M33 taken from Levesque et al. (Reference Levesque, Massey, Olsen, Plez, Josselin, Maeder and Meynet2005), Massey et al. (Reference Massey, Silva, Levesque, Plez, Olsen, Clayton, Meynet and Maeder2009), Neugent et al. (Reference Neugent, Massey, Skiff and Meynet2012), Drout et al. (Reference Drout, Massey and Meynet2012), and Massey & Evans (Reference Massey and Evans2016). The fit is qualitatively good. Based on the comparison between data and model contours, we infer that RSGs typically have masses ⩽25 M_⊙. At higher masses, the number of observed supergiants decreases, especially in the SMC and LMC, which are at lower typical metallicities than the Milky Way. At all metallicities, the presence of relatively long period (~1 000 d) binaries appears to prevent RSGs from getting too cool. The details of the RLOF and CEE experienced by these stars is described in Section 2.1.2. Binary interactions tend to only occur for the higher mass tracks at this orbital period, with RSGs below 12 M_⊙ largely unaffected by being in such a wide binary. The reason being that while the RSGs have similar temperatures, the more massive and luminous RSGs are larger in radius and thus more likely to experience a binary interaction.

Figure 11. HR diagrams focused on the red supergiant branch at three metallicities for single and binary stars. The observational data is from the samples taken from Levesque et al. (Reference Levesque, Massey, Olsen, Plez, Josselin, Maeder and Meynet2005), Massey et al. (Reference Massey, Silva, Levesque, Plez, Olsen, Clayton, Meynet and Maeder2009), Neugent et al. (Reference Neugent, Massey, Skiff and Meynet2012), Drout, Massey, & Meynet (Reference Drout, Massey and Meynet2012), and Massey & Evans (Reference Massey and Evans2016). For Z = 0.004, we use RSGs from the SMC, for Z = 0.008, we use those in the LMC and M33, for Z = 0.020, we use those from the Galaxy and M31. The tracks and observed WR stars are colour coded to show the surface hydrogen mass fraction. The tracks have masses of 8, 10, 12, 15, 20, 30, 40, and 60 M_⊙, the binary star tracks shown have initial periods of 1 000 d and an initial mass ratio of 0.5, the section outlined in white indicates when a binary interaction is taking place. Each greyscale contour represents an order of magnitude in number density of objects.

The next stars to consider are BSGs. These are post-main-sequence objects that lie just off the main sequence at high luminosities as the stars travel through the Hertzsprung gap. Meynet et al. (Reference Meynet, Kudritzki and Georgy2015) investigated these stars using rotating models; here, we investigate similar evolutionary tracks but for binary models. We demonstrate in Figure 12 that we reproduce the trend in bolometric luminosity versus flux weighted gravity for these sources. The BSG models contributing to contours are selected from stellar models that are initially more massive than 8 M_⊙, have M _bol ⩽ −4, 4.4 ⩾ log (T _eff/K) ⩾ 3.9, X _surf ⩾ 0.4, log g ⩽ 3.5, and have completed core hydrogen burning. The observed BSG sample (Kudritzki et al. Reference Kudritzki, Urbaneja, Bresolin, Przybilla, Gieren and Pietrzyński2008; Urbaneja et al. Reference Urbaneja, Kudritzki, Bresolin, Przybilla, Gieren and Pietrzyński2008; U et al. Reference U, Urbaneja, Kudritzki, Jacobs, Bresolin and Przybilla2009; Kudritzki et al. Reference Kudritzki, Urbaneja, Gazak, Bresolin, Przybilla, Gieren and Pietrzyński2012, Reference Kudritzki, Urbaneja, Gazak, Macri, Hosek, Bresolin and Przybilla2013; Hosek Jr. et al. Reference Hosek2014; Kudritzki et al. Reference Kudritzki, Urbaneja, Bresolin, Hosek and Przybilla2014) is split between the three plots at log (Z/Z_⊙) ⩽ −0.55 for Z = 0.004 and log (Z/Z_⊙) ⩾ −0.2 for Z = 0.020 and in between these limits for Z = 0.008. In general, the observed BSGs lie over the contours of greatest probability in our models. At the highest bolometric luminosities and at Z = 0.002, the data deviates to slightly higher flux-weighted gravity than our models predict. At Z = 0.020, this trend reverses, with the data suggesting slightly lower flux-weighted gravities. This may result from a small discrepancy in mass, or perhaps metallicity between models and data.

Figure 12. Flux-weighted gravity versus bolometric luminosity for blue supergiants. Observational data taken from the compilation of Meynet, Kudritzki, & Georgy (Reference Meynet, Kudritzki and Georgy2015). Left-hand panel is for Z = 0.004, central panel for Z = 0.008, and right-hand panel for Z = 0.020. The tracks have masses of 3, 5, 8, 10, 12, 15, 20, 30, 40, and 60 M_⊙, and an initial binary period of 100 d.

3.2.2. Wolf–Rayet stars

Further comparisons can be made to WR stars, especially those that are hydrogen-free. We show the HR diagrams for these stars in the panels of Figure 13 at three different metallicities and for single stars and binary populations. Here the contours represent where we expect to observe completely hydrogen-free stars which corresponds to when the tracks are red. The binary tracks in the right-hand panels now have an initial period of 100 d with a mass ratio of 0.5. These are compared against observed WR stars (Crowther et al. Reference Crowther, Dessart, Hillier, Abbott and Fullerton2002; Hamann, Gräfener, & Liermann Reference Hamann, Gräfener and Liermann2006; Sander et al. Reference Sander, Hamann and Todt2012; Tramper et al. Reference Tramper2013; Hainich et al. Reference Hainich2014; Hainich et al. Reference Hainich, Pasemann, Todt, Shenar, Sander and Hamann2015; Shenar et al. Reference Shenar2016). Figure 13 illustrates a known outstanding problem in that many WC stars and WN stars have cooler surface temperatures than predicted by stellar evolution models (Hamann & Gräfener Reference Hamann and Gräfener2003; Sander et al. Reference Sander, Hamann and Todt2012). One possible explanation is that clumping in the outer convective zone of the star may inflate the envelope (Gräfener, Owocki, & Vink Reference Gräfener, Owocki and Vink2012; McClelland & Eldridge Reference McClelland and Eldridge2016). When we look at the comparison between our stellar models and observed WR stars, especially for hydrogen-free objects, in Figure 13, we see the luminosity distribution is reproduced, but we see the expected temperature offset at LMC and Galactic metallicities, especially for the hydrogen-free WC stars (shown as triangles). Interestingly, we do not see so large an offset in the SMC where the abundances of the observed sources are similar to those of nearby stellar tracks.

Figure 13. HR diagrams with contours plots showing the expected location of hydrogen-free Wolf–Rayet stars in the SMC (top), LMC (middle), and Galaxy (bottom) panels. The left-hand panels are for the single star populations and the right-hand panels for binary populations. The tracks and observed WR stars are colour coded to show the surface hydrogen mass fraction. The observations are taken from Hamann & Gräfener (Reference Hamann and Gräfener2003), Sander, Hamann, & Todt (Reference Sander, Hamann and Todt2012), Hainich et al. (Reference Hainich, Pasemann, Todt, Shenar, Sander and Hamann2015), and Shenar et al. (Reference Shenar2016). Triangles indicate WC/WO stars, crossed WN stars, diamonds indicate the WN3/O3 stars from Neugent et al. (Reference Neugent, Massey, Hillier and Morrell2017). The tracks have masses of 15, 20, 25, 30, 40, 60, 100, and 300 M_⊙, and an initial binary period of 100 d.

McClelland & Eldridge (Reference McClelland and Eldridge2016) found that inflation can decrease the surface temperature of a WR stars by about 0.2 dex. As can be seen, this is roughly the temperature decrease required to secure better agreement between the models and the observed WR stars. We do not include this effect in our current bpass models as the physics is uncertain and we are still unclear whether the fitting method of the atmosphere models is entirely correct. If inflation effects are metallicity dependent, this may explain why the offset is less pronounced at LMC metallicities. The relatively minor difference in temperature has little effect on the integrated spectrum of a stellar population, but we note that caution is needed when evaluating the properties of an individual WR star through model fitting.

Recently, Neugent et al. (Reference Neugent, Massey, Hillier and Morrell2017) discovered a new class of WR stars, WN3/O3 stars. We have included them in the LMC panels in Figure 13 as diamond symbols. We see that the models agree well with the observed luminosity, temperature, and surface composition of both binary and single star stellar tracks. This suggests therefore that these are typical WR stars, however their reported mass-loss rates indicate that the rates used in our models may be too high for such stars. We note that it has been suggested by Smith, Gotberg, & de Mink (Reference Smith, Gotberg and de Mink2017) that these may be less massive helium stars rather than WR, which would also be consistent with our binary models. However, the projected distance from O stars they find of a few hundred parsecs could be understood if these objects are also runaway stars. The results of Eldridge et al. (Reference Eldridge, Langer and Tout2011) for type Ib/c SN progenitors suggest that 10 to 20% of WR stars could travel 100 parsecs, and further, away from their birth places (see their Figure A3).

3.2.3. Quasi-chemically homogenous evolution

As discussed in Section 2.1, a small subset of evolved stars are expected to undergo quasi-homogenous evolution and these models too require observational verification. QHE stars for the most part should look identical to more normal WR stars so will be difficult to separate out, and indeed there are no clearly unambiguous examples known. However, their chemical composition may aid identification. Most typical WR stars form helium cores and then lose their hydrogen envelopes. By comparison, QHE stars, when hydrogen rich, are still core-hydrogen burning so will exist for a longer period of time and also have different compositions at locations on the HR diagram close to the traditional hydrogen-burning main sequence.

In Figure 14, we compare the QHE tracks (shown with bold lines and evolving towards higher temperatures) to standard binary evolution tracks (thin lines, evolving first towards lower temperatures). Contours give the probability distribution of all models in our population synthesis which experience QHE. We see that in general the WR stars are in better agreement with the standard evolution tracks than with QHE models. As noted by Hainich et al. (Reference Hainich, Pasemann, Todt, Shenar, Sander and Hamann2015) and Shenar et al. (Reference Shenar2016), it appears that mass-loss rates may have been slightly overestimated in current models. We see however that the lowest luminosity WR star in this sample is much cooler and more hydrogen rich than the other stars. Its location matches quite closely that of our 40 M_⊙ QHE track. This star, AB2, is perhaps the best candidate for a QHE (by mass transfer) star in the nearby Universe and deserves further study (Martins et al. Reference Martins, Hillier, Bouret, Depagne, Foellmi, Marchenko and Moffat2009; Hainich et al. Reference Hainich, Pasemann, Todt, Shenar, Sander and Hamann2015).

Figure 14. HR diagram showing the evolution of QHE models at Z = 0.004 compared to the properties of WR stars in the SMC. Crosses indicate WR stars, asterisks are the O star companions from Shenar et al. (Reference Shenar2016). Colour coding is the surface hydrogen abundance. The thicker tracks are for QHE stars, while thinner tracks are standard binary models. The tracks have masses of 20, 25, 30, 40, 60, 100, 300 M_⊙ and an initial binary period of 100 d.

3.2.4. Stellar-type ratios

An important validation test to consider is not just whether we reproduce the location of these stars in the HR diagram but whether we reproduce the correct relative numbers. We show in Figure 15 the predicted numbers of massive stellar types (log(L/L_⊙) > 4.9), specifically O stars, BSG, YSG, RSG, WR (total), WN, and WC stars as a function of metallicity for constant star formation and also with age assuming an instantaneous starburst. We define the surface parameters of each stellar type in Table 3. O stars and BSGs are the most populous classes at all metallicities, as expected for the stars that exist when stars burn hydrogen to helium. The numbers of the remaining subtypes vary smoothly with metallicity. In general, the WR stars decrease in number at lower metallicities since the formation of these requires mass loss and the action of stellar winds is less efficient in the absence of metals. We note that binary populations, which provide another avenue for mass loss, exhibit a weaker metallicity trend than the single star population.

Figure 15. Number of massive stars with log (L/L_⊙) > 4.9 of different stellar types from our fiducial populations. The stellar types are defined in Table 3. Dashed lines are for single stars, solid binary populations. The dot-dashed line for the WR population includes lower luminosity, lower mass stars which otherwise satisfy the temperature and surface abundance criteria for WR stars. Top panel shows the metallicity variation of stellar type ratios assuming a population undergoing continuous star formation at a rate of 1 M_⊙ yr⁻¹. Lower panels show time evolution of an instantaneous burst (i.e. a single-aged stellar population) at two different metallicities. Where lines are interrupted, zero stars in that type category remain in the (finite-sized) population in a given time bin.

Considering the evolution of an instantaneous burst, we see that the different stellar types appear at different stellar population ages. Binary interactions extend the lifetime of the WR stars, especially WC stars, to later times. This highlights the difficulty of precisely aging a stellar population from the fact that WR stars exist alone. In our RSG numbers, we see a late time peak at around 100 Myrs in both metallicities shown in the figure. This is due to a small contribution of AGB stars at these ages which meet the temperature and luminosity constraints for selection as RSGs in our model. Observationally, they would also be difficult to distinguish.

Similar predictions to these have been used for many years to constrain stellar models by comparing the model to the number of different stellar types observed in an entire galaxy, ideally at known metallicity (e.g. Maeder & Meynet Reference Maeder and Meynet1994). Here, we perform a similar test, comparing bpass predictions of the WR to O, BSG to RSG, WR to RSG, and WC to observed WN ratios in nearby galaxies (see figure caption for data sources). One problem with such comparisons is uncertainty in how complete each observational sample is, especially for the WR to O star ratio where both stellar types are hot and difficult to find in optical surveys. The comparison will also be somewhat dependent on the star formation history and here we assume continuous star formation at a constant rate.

The panels in Figure 16 show the current bpass v2.1 predictions (black and blue lines given two IMF mass limits, M _max = 100 or 300 M_⊙) for these ratio in comparison to the observed ratios which exist and also previous bpass v1.0 predictions (red line). For all these number and stellar type number ratios, we only include stars with log (L/L_⊙) > 4.9, since Massey & Olsen (Reference Massey and Olsen2003) used this as a luminosity cut-off to make sure that AGB stars are not included in their sample (although we note that a small number of our AGB models do indeed satisfy this criterion). WR stars do not typically have luminosities below this value (McClelland & Eldridge Reference McClelland and Eldridge2016). It is consistent to also apply this to our O star number predictions for comparison to the data, but we note that this excludes the many O stars that lie just below this limit in our models. In the upper left panel of Figure 16, we include predictions for the WR/O ratio that include all O stars as well as only those above our luminosity limit (dotted and dot-dashed lines for binary and single stars, respectively).

Figure 16. Stellar-type number ratios for massive star populations. Solid lines show bpass binary models and dashed lines show single models. Points with error bars show observational constraints taken from the literature. All ratios include only stars with log (L/L_⊙) > 4.9, except the WR/O star ratio where the dotted (binary) and dash-dotted (single) lines include all O stars. The red lines show results for constant star formation from bpass v1.0, while the blue and black lines show similar results from BPASS v2.1 with two different IMFs (M _max = 300 and M _max = 100 M_⊙, respectively). WR/O ratio data is taken from Maeder & Meynet (Reference Maeder and Meynet1994), the BSG/RSG ratio data is from Massey & Olsen (Reference Massey and Olsen2003) and the WR/RSG data comes from P. Massey (private communication). The WC/WN ratio data plotted in blue and red and is taken from Rosslowe & Crowther (Reference Rosslowe and Crowther2015) with the blue points omitting the dusty WC stars. The WC/WN ratios plotted in black are from Meynet & Maeder (Reference Meynet and Maeder2005).

At Solar metallicity, the effect of binaries, and of omitting the lower luminosity O stars, is small. These effects grow at decreasing metallicity, since stars at low metallicity have hotter effective temperatures for the same mass, and thus more O stars are present in the population. Comparison to the observed ratios suggests that models including all main-sequence O stars provide a good agreement to the data, except in the highest metallicity cases. The reason for this is probably related to the fact that the observed sample is old and may suffer from significant incompleteness. Modern, more complete, surveys may find quite different results, but it is encouraging that bpass models reproduce the observed trend in the observed ratios.

The BSG/RSG ratio observational sample (Massey & Olsen Reference Massey and Olsen2003) only contains ratios for the SMC and LMC. Again binary evolution predicts values within a factor of a few of the observed ratios, but somewhat underpredicts the observed BSG fraction, likely due to the lack of rotational mixing in our models. More rotational mixing would lead to more BSGs (Castro et al. Reference Castro, Fossati, Langer, Simón-Díaz, Schneider and Izzard2014) and so we would expect their inclusion to boost our predicted ratio in future and consider the fraction in bpass v2.1 to be a lower limit. We note that the mist models (Choi et al. Reference Choi, Dotter, Conroy, Cantiello, Paxton and Johnson2016) predict a higher BSG/RSG ratio than bpass.

The WR/RSG ratio was presented in our first study with the bpass stellar models (Eldridge et al. Reference Eldridge, Izzard and Tout2008). In v1.0, it was our worst fit of the stellar population ratios since the observational data showed a pronounced trend in ratio with metallicity which the models did not reproduce. Here we compare our models to updated observational numbers (P. Massey, private communication). We see that our binary populations are far closer to the new observed ratios which show a significantly weaker trend. The observational data now shows very little variation in population ratio with metallicity above ~0.3 Z_⊙, consistent with the ratios in our models. In the bpass v2.1 model population, the ratio continues to remain flat at lower metallicities. In the observational data set, the SMC is an outlier. However, this ratio is based on only 13 WR stars, suggesting a considerable uncertainty in small number statistics, and potential incompleteness. Alternatively, as the observed data point lies between the single star and binary population predictions from bpass, it is possible that in the case of the SMC, the effective binary fraction could be smaller.

Finally, the WC/WN ratio is the focus of our final comparison between stellar populations and is primarily determined by the WR mass-loss rates. In this figure, we show the latest observations from Rosslowe & Crowther (Reference Rosslowe and Crowther2015) as well as the older sample from Meynet & Maeder (Reference Meynet and Maeder2005). In Eldridge et al. (Reference Eldridge, Izzard and Tout2008), we found that bpass binary WC/WN ratio predictions lay below what was then thought to be the observed ratio. Vanbeveren, Van Bever, & Belkus (Reference Vanbeveren, Van Bever and Belkus2007) discussed in detail how binary populations affect this ratio: binary interactions produce WR stars from low-mass stars which would not otherwise evolve into the WR phase in the single star evolution case. The products of such evolutionary pathways tend to be WN stars rather than WC stars and hence more interacting binaries in a model set causes a drop in the WC/WN ratio. Recent observational surveys presented here show that the ratio is indeed lower than previously thought, in closer agreement to the binary populations predictions. The excellent agreement shown in the bottom right panel of Figure 16 is strong circumstantial evidence that binary interactions are responsible for a good number of WR stars and that single-star stellar winds are not strong enough to create every WR star we see in the sky.

In summary, our population models are consistent with the observed location of post-main sequence massive stars on the HR diagram, as well as the relative numbers of stars we see in galactic stellar populations. This is important as RSG and WR stars have unique spectral signatures in the emergent spectrum from a stellar population. RSGs dominate the longer redder wavelengths, while WR stars are bright in the UV, produce hard ionising radiation, and also strong broad stellar emission lines. The ability of our models to reproduce the stars in resolved populations is good evidence that we likely predict these spectral features accurately.

6.6.1. Ionising photon production efficiency

A key parameter for constraining the reionisation history of the Universe, and one of the most significant to change behaviour between bpass v2.0 and v2.1, is the ionising photon production efficiency $\xi _{\mathrm{ion}} = \dot{N}_\mathrm{ion}/f_\mathrm{esc}\,L_\mathrm{UV}$ , i.e. the ratio between the ionising photon production rate and the rest-frame ultraviolet continuum luminosity for a given fraction of such photons which escape their host galaxy. While the luminosity can be directly observed in the distant Universe, the ionising photon production rate must be inferred from either models or indirect observations of the strong nebular emission lines generated from the galaxy in question, assuming that a fraction (1 − f _esc) of the ionising continuum is absorbed by the interstellar medium.

For simplicity, an escape fraction f _esc = 0 is often assumed, yielding ξ_0,ion, but in practice the paradigm that the intergalactic medium is ionised by early star formation, and that ionisation state is maintained over cosmic time, requires a much higher escape fraction. Observations aimed at constraining that escape fraction compare measured luminosities at wavelengths shortwards of the ionisation edge of hydrogen at λ(rest) = 912 Å with models of what flux the stellar population is expected to produce, given the observed 1 500 Å luminosity. ξ_0,ion is thus a crucial parameter both for estimation of escape fraction and for interpretation of the effects of star-forming galaxies on the reionisation epoch.

Between v2.0 and v2.1, we have made two important changes to the models that affect this parameter: we have extended our models to still lower metallicities, and we have improved the stellar atmosphere models in use for the highest gravity O stars (see Section 2.3.2). These models are particularly important at low metallicities (where mass loss is lowest and stars can thus remain more massive), and in the rotationally mixed stars which dominate the ionising photon production.

In Figure 35, we present a direct comparison between our previous estimates of ionising photon production rate, ultraviolet luminosity, and thus ξ_0,ion and new estimates from v2.1 for the case of continuous star formation observed at 100 Myr after the onset of the star-formation epoch. There are three clear differences in behaviour between our improved v2.1 models and the older v2.0 models. First, v2.1 does not show the dramatic increase in ionising photon production rate at metallicities below Z = 0.004, which was associated with the onset of quasi-homogenous chemical evolution—while this still leads to increased photon production, it is moderated by the new atmosphere models. Second, the v2.1 models are more luminous in the rest-frame ultraviolet than the v2.0 models and the difference between single and binary models has narrowed due to a more rapid build up of continuum flux than in our previous models. Third, as a result of these behaviours, the metallicity evolution of the ionising photon production efficiency has become less clear cut. In v2.0, the binary population always exceeded the single star population in ξ_0,ion at a given metallicity. In v2.1, we see lower photon production efficiencies in the binary population than in single stars at metallicities above ~0.3 Z_⊙ in our standard IMF, and near-equal efficiencies in IMFs excluding the most massive stars. This behaviour reverses at lower metallicities, where the binary population continues to dominate.

Figure 35. A comparison between the ionisation characteristics of bpass v2.0 and v2.1 models for a population which has been continuously forming stars at a rate of 1 M_⊙ yr⁻¹ over a 100 Myr period, shown as a function of metallicity. In the top panel, we show the small effect of our improved atmosphere models on the ionising photon production rate. By contrast, the centre panel shows a strong change in rest-ultraviolet luminosity between v2.0 and v2.1. As a result, the lower panel shows a change in the behaviour in ξ_0,ion. Single star models are shown with dashed lines, while binary models are shown with solid lines. Our standard IMF is shown (labelled imf300), with an identical IMF with M _max = 100 M_⊙ also given for comparison.

This behaviour is entirely driven by the behaviour of the youngest stars in a given population. As Figure 36 demonstrates, the ionising photon production output of a single-aged, instantaneous-starburst binary population exceeds that of the single star population at all metallicities for ages >3 Myr but the stars below this age produce photons with a far higher efficiency. In the low metallicity populations, the longer ionising photon production lifetime of the binary stars leads to the binary population overtaking the single stars at late ages for any ongoing or evolving starburst, while at near-Solar metallicity the rapid mass loss due to stellar winds dominates over binary effects.

Figure 36. The time evolution of ionising photon production efficiency as a function of star-formation history shown for metallicities of ~Z_⊙ and 0.1 Z_⊙. Single star models are shown in greyscale, with dashed lines, while binary models are shown in colour with solid lines. The star-formation histories shown are those defined in Figure 29.

The decrease in overall ionising photon efficiency in bpass v2.1 relative to v2.0 reflects an improvement in stellar atmosphere modelling, but also emphasises the significant uncertainties in the stellar models, particularly in the low metallicity regime, and thus the difficulty of interpreting distant galaxies. We note that an improved treatment of rotation (see Section 7.2) could reverse the shift seen in this version.

In Figure 37, we compare our current predictions for ionising photon efficiency with observational constraints. These are drawn from the literature and based on a standard analytic conversion from Hα emission line flux to ionising photon production rate (which itself is somewhat model dependent). Distant galaxies are shown as shaded regions indicating the mean and standard deviation of the (often small) measured sample in ξ_0,ion. The metallicity is poorly constrained in these observations and samples are shown offset to indicate a plausible range in metallicity mass fraction. We also show data for a sample of extreme local starburst galaxies identified as analoguous to the z ~ 5 galaxy population for which metallicity information can be inferred from the optical spectra and allows the data to be binned (Greis et al. Reference Greis, Stanway, Davies and Levan2016). In previous versions of this comparison, we have considered a model population continuously forming stars at 1 M_⊙ yr⁻¹ over 100 Myr; now we compare with models at 10 Myr after the onset of star formation although, as Figure 29 demonstrated, an alternative approach would be to assume a more bursty and less constant star formation history.

Figure 37. The ionising photon production efficiency as a function of metallicity and star-formation history for an ongoing star-formation event which has formed stars continuously over the previous 10 Myr. Grey regions indicate measured values of ξ_0,ion in the distant Universe together with their uncertainties. These samples are unconstrained in metallicity and are shown offset to indicate plausible metallicity ranges given BPASS population synthesis models. Data points show results for a sample of local extreme star-forming galaxies (Lyman Break analogues) from Greis et al. (Reference Greis, Stanway, Davies and Levan2016), binned by optical emission line metallicity. Star-formation histories are colour-coded as in Figure 36.

As can be seen, a single age starburst would struggle to reproduce the high ξ_0,ion values inferred from both distant and local extreme starburst observations, unless captured in the very first few million years. By contrast, an ongoing or rising starburst history provides a reasonable match to the observed data at ages of order ~10 Myr.

6.6.2. Optical recombination line ratios

In Stanway et al. (Reference Stanway, Eldridge, Greis, Davies, Wilkins and Bremer2014), we demonstrated that, subject to the usual uncertainties in nebular gas conditions, our v1.1 single age starburst models were capable of reproducing the optical strong line ratios seen in the z ~ 2 − 3 galaxy population. In Figure 38, we update this comparison incorporating our newer v2.1 models. For consistency, we use the same nebular emission model as in the above work, with a fixed electron density (100 cm⁻³) and cloud geometry (spherical, 10 pc inner radius), and allow stellar population age to vary along tracks. We confirm our earlier finding that young stellar populations incorporating binary evolution pathways can reproduce both the strong line ratios and the Hβ equivalent widths seen in the distant population, while single star models struggle to do so.

Figure 38. The evolution in Balmer emission line strength and [O iii]/Hβ line ratio with stellar population age. Age increases along the tracks from high equivalent width to low. Tracks for binary models are shown in colour and those for single star models in greyscale, with different line styles. Overplotted are the data for star-forming galaxies at z = 2 − 3 reported by Holden et al. (Reference Holden2016, triangles) and Schenker et al. (Reference Schenker, Ellis, Konidaris and Stark2013, squares). Greyscale indicates the measured properties of the local SDSS star-forming galaxy population.

We also show the measured properties of the SDSS star-forming galaxy population, as defined above. Interestingly, these show an offset from both the binary and single star model populations, lying between the two but more consistent with the binary tracks. We note that these have not been corrected for internal dust extinction (which may affect lines more strongly than continuum in some extinction models, and so would push the data towards the single star tracks in equivalent width) but are corrected for fibre losses. We consider only a simple star-formation history in this case, since the line emission spectrum is heavily dominated by the youngest starburst. A more complex star-formation history will tend to slightly reduce the model Hβ equivalent width, again pushing the data towards the single star models. However, as Figure 38 also demonstrates, the evolution of the single star models in line equivalent width is extremely rapid, and so the SDSS population would have to be caught in a very narrow time window some 10 Myr after the onset of star formation. By contrast, the binary models provide high Hβ equivalent width over a more extended time period increasing the probability of populating the parameter space occupied by the SDSS sources. If binary evolution models are preferred, there are two likely explanations for the offset of the SDSS galaxies from our models. The first may be use of an inappropriate gas density or cloud geometry for the local population—adjusting these will tend to reduce the model [O iii]/Hβ ratio at a given Hβ equivalent width, bringing the binary models better in line with the data (see Stanway et al. Reference Stanway, Eldridge, Greis, Davies, Wilkins and Bremer2014). Alternately, the binary distribution parameters assumed (see Table 1), which appear to reproduce those in the distant Universe and simple massive stellar populations, may overestimate the interacting binary fraction in high metallicity galaxies in the local Universe with more complex star-formation histories and old underlying stellar populations. The result of assuming a lower binary fraction in low mass stars would be to move the binary models towards the local galaxy data, as further discussed in Section 7.10.

In Figure 39, we show an alternate ionisation diagnostic—the Baldwin, Phillips, & Terlevich (Reference Baldwin, Phillips and Terlevich1981, BPT) diagram widely used to distinguish emission powered by star-forming regions from that powered by accretion onto an AGN. We show the locus occupied both by local star-forming galaxies (greyscale) and the distant galaxy population reported by Steidel et al. (Reference Steidel, Strom, Pettini, Rudie, Reddy and Trainor2016, data points). Overlaid are evolutionary tracks for a co-eval simple stellar population with age, at a range of metallicities, assuming the same simple nebular model as above. At ages of ~10 Myr, the tracks cross the Kauffmann et al. (Reference Kauffmann2003) boundary for separating local star-forming galaxies, and, as before, more complex star-formation histories are required to overlay the SDSS galaxy population, which typically lies somewhere between the ratios expected for a single age starburst and constant, ongoing star formation. A full exploration of the BPT parameter space occupied by bpass models with different assumed nebular gas conditions will be presented in Xiao et al. (Reference Xiao2017).

Figure 39. The evolution of strong optical emission line ratios with stellar population age. Age increases along the tracks from 1 Myr to 30 Myr for instantaneous burst models. Small circles are data for star-forming galaxies at z = 2 − 3 reported by (Steidel et al. Reference Steidel, Strom, Pettini, Rudie, Reddy and Trainor2016). An indicative typical error bar for the data is shown at the bottom right. Greyscale indicates the measured properties of the local SDSS star-forming galaxy population. The solid grey line shows the Kauffmann et al. (Reference Kauffmann2003) condition for distinguishing regions ionised by AGN from star forming galaxies. Triangles indicate the line ratios expected for a binary stellar population with constant star formation over the last 100 Myr.

An additional validation test for our models, and an important demonstration of the ambiguity that arises from nebular gas parameters, is their ability to reproduce the line ratios seen in GRB host galaxies. The presence of a GRB requires these galaxies to have formed massive stars recently, and thus (as with the WR-galaxies discussed above) their nebular emission tend to be dominated by a young stellar population. We consider here the TOUGH sample, observed spectroscopically by Krühler et al. (Reference Krühler2015) with X-SHOOTER on the ESO VLT, for which metallicity and extinction estimates are available. In Figure 40, we consider a third diagnostic line ratio, [OII] λ3726,3729 Å/[OIII]λ5007 Å. The evolution of this line ratio is not entirely smooth in either stellar population age or nebular electron density, with both playing an important role in exciting these forbidden lines. Nonetheless, considering a range of plausible ages (1 to 100 Myr) and electron densities (0.1 to 1000 cm⁻³) makes clear the locus of permitted line ratios relative to [OIII]/Hβ and their good agreement with the observed properties of GRB hosts (upper panel). Solid lines indicate a polynomial fit to the variation with age at each metallicity and demonstrates the relative insensitivity of this parameter space to metallicity in the range Z ~ 0.5–0.7 Z_⊙, but also the significant scatter around the mean locus given a range of assumptions regarding age and nebular gas conditions. This point is perhaps still clearer in Figure 40 (right-hand panel), which shows the evolution in [OII]/[OIII] as a function of metallicity at three stellar population ages and three gas densities. There is significant degeneracy between different plausible models. If interpreted at fixed age and electron densities, the distribution of [OII]/[OIII] line ratios seen in the GRB hosts of Krühler et al. (Reference Krühler2015) is consistent with a metallicity sequence. However, a variation in stellar population age, or electron density, or both, at fixed metallicity could also provide a good fit to the observations. Given that the metallicities shown for the data points were themselves derived from the optical recombination line spectrum, and thus ultimately derive from a photoionisation model and stellar population fit to a calibration data set, the correct interpretation is unclear, and likely involves some evolution in all three parameters.

Figure 40. The emission line properties of GRB host galaxies compared to BPASS v2.1 models. Both panels consider the [OII]/[OIII] excitation diagnostic. In the left-hand panel small coloured points indicate the predicted line fluxes for models with ages from 1 Myr to 100 Myr, and at electron densities ranging from 0.1 to 1000 cm⁻³. Solid lines are a polynomial fit to all models with n _e = 100 cm⁻³ as a function of metallicity and indicate the general trend in the permitted parameter space. The right-hand panel shows the [OII]/[OIII] ratio explicitly as a function of metallicity at three stellar population ages, and three electron densities, showing the ambiguity between an older population and a sparser ISM at any given metallicity. GRB host galaxy data are taken from Krühler et al. (Reference Krühler2015) and are shown in black for hosts at z > 1 and in grey for hosts at z < 1, metallicities for data points are those derived from optical strong line ratios by Krühler et al. (Reference Krühler2015). Only objects with measurements in all relevant quantities, and an estimated E(B − V) < 0.2 are shown. Line ratios have been corrected for internal extinction using the Calzetti dust law.

6.6.3. Ultraviolet spectral features

Ultraviolet emission lines, redshifted into the observed frame optical and near-infrared are crucial evidence for the stellar populations and properties of distant galaxies. A recent study by Bowler et al. (Reference Bowler, McLure, Dunlop, McLeod, Stanway, Eldridge and Jarvis2016) has already made use of beta versions of our lowest metallicity v2.1 models to explore the interpretation of the anomalous source CR7 (Sobral et al. Reference Sobral, Matthee, Darvish, Schaerer, Mobasher, Röttgering, Santos and Hemmati2015) which has been proposed as either a Population III starburst or, potentially, a direct collapse black hole. While none of the model sets explored in that work provide an unambiguous interpretation for CR7, bpass low metallicity, young starburst models are broadly consistent with the colours of the source in question (particularly if some degree of enhancement in α-elements is assumed) suggesting a stellar origin for the light. The extreme properties of CR7 (specifically its He ii line strength) have recently been challenged (Shibuya et al. Reference Shibuya2017). Nonetheless, this source and the similar high redshift emission line sources now being identified provide one of the very few constraints on our models at the lowest metallicities. Crucially, they are representative of a regime that may soon be more widely explored by JWST.

As is the case for optical lines, interpretation of the ultraviolet spectrum of a complex stellar population requires an assumption regarding the nebular gas conditions, which are largely unconstrained by observations. Steidel et al. (Reference Steidel, Strom, Pettini, Rudie, Reddy and Trainor2016) have suggested that the spectra of galaxies at z ~ 2 − 3 are best fit by combining very low metallicity stellar spectra (the observed properties of which are most strongly influenced by the iron abundance) with somewhat higher metallicity nebular gas (the properties of which are dominated by oxygen abundance). Again, this may be indicative of strong α-element enhancement, as would be expected in the young stellar populations at this redshift.

The difficulty in evaluating nebular effects is particularly true of resonantly scattered lines such as the Lyman-α feature, the ratio of which to the Balmer series is shown in Figure 41 for a continuously star-forming population, assuming our usual gas geometry and local electron densities ranging from 1 to 10⁴ cm⁻³. In Figure 42, we simply assume n _e = 10² cm⁻³ and Z(nebular) = Z(stellar), in order to demonstrate emission line strengths as a function of age and metallicity for our standard cloudy models, normalised relative to the Lyman-α line strength in the same galaxy. Perhaps, the most striking effect of decreasing metallicity is the increase in strength of rest-frame ultraviolet emission lines that are seldom seen in normal stellar populations in the local Universe, and which are usually associated with emission from AGN when observed. These lines are not significantly excited in nebular gas irradiated by our equivalent single star models.

Figure 41. The ratio of intrinsic Lyman-α emission to optical Balmer line ratio assuming a range of metallicities and gas parameters. While the Balmer line ratio is only weakly sensitive to the environment, Lyman-α is strongly affected by the metallicity and physical conditions of the nebular gas. Electron densities increase from 1 cm⁻³ to 10⁴ cm⁻³ along coloured (fixed metallicity) lines and are marked at 1 dex intervals. We assume constant star formation, observed at age 10⁸ years.

Figure 42. The evolution in rest-frame ultraviolet emission line flux ratios for a single-aged simple stellar population at four representative metallicities. Line strengths shown are normalised relative to the Hydrogen Lyman-α λ1216 Å feature, and are He ii λ1640 Å, C iii] λ1909Å, and O iii] λ1665 Å. For doublets, the total flux is given.

At the lowest metallicity in our models (Z = 1 × 10⁻⁵ = 0.05% Z_⊙), He ii λ1640 Å is the strongest emission line in the rest-UV longwards of the Lyman-α line. At the low to moderate metallicities, more typical of the distant Universe (0.5 to 10% Z_⊙), the semi-forbidden doublets of C iii] λ1909 Å and O iii] λ1665 Å increasingly exceed the helium line in line flux. All three species decline in strength towards Solar metallicity due to the decline in ionising photon strength emitted by stellar populations.

The [C iii] λ1907Å, Ciii] λ1909 Å doublet has attracted particular attention of late as the result of its detection in a handful of very distant (z > 6) star-forming galaxies (Stark et al. Reference Stark2015, Reference Stark2017), and its position in the rest-frame ultraviolet, which means it is observable in the optical in the distant Universe. This line may arise from the atmospheres and winds of massive stars (in which case it will typically be broad) and from ionised nebular gas in H ii regions (in which it is typically narrow). Both of these are seen locally in the environs of WR stars and planetary nebulae (e.g. Guerrero & De Marco Reference Guerrero and De Marco2013). Observed lines in distant sources tend to be narrow, suggesting a nebular origin, but the constraints on its velocity width are relatively poor in most cases. In Figure 43, we show the predicted total equivalent width of this doublet as a function of stellar population age and metallicity for a population continuously forming stars at 1 M_⊙ yr⁻¹, after passing through nebular gas. These are compared to the distribution of measured rest-frame equivalent widths in this doublet arising from star-forming galaxies at z = 0–8 (Stark et al. Reference Stark2014; Rigby et al. Reference Rigby, Bayliss, Gladders, Sharon, Wuyts, Dahle, Johnson and Peña-Guerrero2015; Vanzella et al. Reference Vanzella2016; Patrício et al. Reference Patrício2016; Du et al. Reference Du, Shapley, Martin and Coil2017, Maseda et al. Reference Maseda2017). We also indicate the rest-frame equivalent widths measured for four stacks of z ~ 3 LBGs, sorted into quartiles by Lyman-α equivalent width and each containing ~200 galaxies, by Shapley et al. (Reference Shapley, Steidel, Pettini and Adelberger2003). Unsurprisingly, the bulk of the observed measurements are consistent with the very low equivalent widths we expect for mature stellar populations at moderate metallicities. It is notable however that, at both z < 3 and z > 3, there exists a tail of sources extending to $W_{\rm {0,C {\rm iii]}}}>10$ Å. Comparison with our standard cloudy models suggest that these are most naturally explained by stellar populations with significantly sub-Solar metallicities, Z < 0.2 Z_⊙, or very young dominant stellar populations (<10 Myr). In the most extreme observed cases, both of these conditions may be necessary. At the very lowest metallicities in the bpass model set (Z < 0.05% Z_⊙), the hard stellar ionising spectrum excites very little C iii] emission simply because of the scarcity of carbon atoms in the nebular gas. As observations approach this regime in the very distant Universe, it is worth noting that, in the presence of a carbon-enriched ISM, the irradiating spectrum will likely excite strong line emission and that such enrichment may result from massive stellar explosions very rapidly after the onset of star formation. As discussed before, there is scope for further populating this parameter space by varying electron density, local ionisation parameter, and potentially the α-element enhancement of the nebular gas. However, the C iii] doublet is of particular interest since its line ratio is itself sensitive to electron density and thus potentially allows some such degeneracies to be broken. As samples of deep rest-frame ultraviolet spectra grow, it will likely be possible to break many such degeneracies with observations of both the doublet line ratio and further emission lines in future.

Figure 43. The rest-frame equivalent width of the (unresolved) [C iii] λ1907 Å, C iii λ1909 Å doublet, as a function of stellar population age and metallicity for a population continuously forming stars at the rate of 1 M_⊙ yr⁻¹. In the right-hand panel, we compare these predictions to the distribution of measured C iii 1909 lines in star-forming galaxies compiled from the literature (see Section 6.6.3). Measurements from sources at z > 3 are highlighted in red, while arrows indicate the measurements obtained from stacks of z ~ 3 Lyman break galaxies differing in Lyman-α line strength (Shapley et al. Reference Shapley, Steidel, Pettini and Adelberger2003).

A fuller exploration of the limited observational data on nebular line emission at high redshift is deferred to later work.