Evaluating the impact of binary parameter uncertainty on stellar population properties

Binary stars have been shown to have a substantial impact on the integrated light of stellar populations, particularly at low metallicity and early ages - conditions prevalent in the distant Universe. But the fraction of stars in stellar multiples as a function of mass, their likely initial periods and distribution of mass ratios are all known empirically from observations only in the local Universe. Each has associated uncertainties. We explore the impact of these uncertainties in binary parameters on the properties of integrated stellar populations, considering which properties and timescales are most susceptible to uncertainty introduced by binary fractions and whether observations of the integrated light might be sufficient to determine binary parameters. We conclude that the effects of uncertainty in the empirical binary parameter distributions are likely smaller than those introduced by metallicity and stellar population age uncertainties for observational data. We identify emission in the He II 1640 Angstrom emission line and continuum colour in the ultraviolet-optical as potential indicators of a high mass binary presence, although poorly constrained metallicity, dust extinction and degeneracies in plausible star formation history are likely to swamp any measurable signal.


INTRODUCTION
Stellar population and spectral synthesis (SPS) models provide the bridge between unresolved observations of galaxies or stellar clusters, and interpretation of their stellar content. Such models are routinely fit to both photometric and spectroscopic data, yielding stellar masses, population ages, star formation histories and sometimes metallicity or even redshift (see e.g. Tinsley 1975;Tinsley & Gunn 1976;Leitherer et al. 1995;Bruzual & Charlot 2003;Conroy et al. 2009;Conroy & Gunn 2010;Conroy 2013). Each SPS model is constructed from models for the time evolution of individual stars, tracking changes in both their physical properties and their observable characteristics. A simple stellar population (SSP) accounts for the fragmentation of an initial star-forming region into stars according to an initial mass function (IMF) and traces the evolution of the resulting single-aged stellar population. A composite stellar population (CSP) combines SSPs of different ages and masses to construct a population with a known star formation history.
The implicit assumption underlying this process is that the evolution and resultant atmospheres of individual stars ⋆ E-mail: e.r.stanway@warwick.ac.uk are well understood, and that their distribution in initial mass can be described by a simple function. While the initial mass function is largely a valid construct for the ensemble of typical stellar populations in the local Universe, there is increasing evidence for variation between galaxies (see , and low mass systems present a challenge since their initial mass function must be sampled stochastically rather than statistically (see Eldridge 2012). Similarly, there are regimes in which theoretical understanding of stellar evolution is still undergoing active and rapid development. In particular, the evolution of stars at very low metallicity, at very high mass, undergoing rotation or influenced by binary interactions, remain subject to significant uncertainties (e.g. Langer 2012;de Mink et al. 2013;Paxton et al. 2015;Vanbeveren & Mennekens 2017;Eldridge et al. 2017;Chrimes et al. 2020;Shenar et al. 2020, and references therein). The determination that the majority of massive stars interact with a binary companion during their evolutionary lifetime (Sana et al. 2012(Sana et al. , 2014, combined with advancements in extragalactic astronomy which push towards the young, metal-poor stellar populations at ever higher redshifts (e.g. Steidel et al. 2018;Strom et al. 2018;Reddy et al. 2018;Theios et al. 2019), now require these uncertainties to be addressed, or at least acknowledged, in SPS modelling and the fitting of galaxy data. In particular, interest has focused on the role of binary and multiple stars in shaping the evolution of a population. This has been shown to have a particularly significant effect on the ionizing photon production and its time evolution (e.g. Stanway et al. 2016;Götberg et al. 2019), and is evidently crucial in interpreting the rates of stellar compact mergers (e.g. de Mink & Belczynski 2015;Tang et al. 2020) The Binary Population and Spectral Synthesis (BPASS, Eldridge et al. 2017;Stanway & Eldridge 2018) project produces simple stellar population models from a custom grid of detailed stellar evolution tracks which incorporate not only the evolution of isolated single stars, but also binary interactions resulting from a range of binary initial mass ratio and initial separations. The use of detailed stellar models allows the response of individual stars to interactions, mass loss and certain rotational effects to be traced in detail, rather than approximated. This is important in regimes in which the response of a star does not vary smoothly with mass, mass ratio, etc. and in determining the mass of core-collapse remnants created at the end of stellar evolution. However, as a consequence, the detailed approach lacks the speed of synthetic model populations based on algorithmic interpolation between a sparser grid of models (e.g. Hurley et al. 2002;. A typical detailed binary population and spectral synthesis model requires about 8 hours of processor time to run (assuming the underlying stellar models already exist and have been pre-processed to the correct format and resolution). The current BPASS data release model set (as of v2.2) incorporates 9 representative initial mass functions (IMFs), with or without including binary evolution tracks, each calculated at 13 metallicities -a total of 208 models which must be run independently, with the results output at 51 time steps. Our recent studies on gravitational wave transient rates have also involved trialling two different supernova kick velocity models (doubling the number of models again), albeit for a single IMF (Tang et al. 2020). In order to make this manageable with current computing resources, each model is calculated with only one set of assumptions for the properties of the binary population as a function of mass. However as computing power, and more importantly, the number of cores and RAM available per core, continues to increase, more ambitious model runs become possible.
In Stanway & Eldridge (2019) we explored the impact of varying our initial mass function, for a grid of variant models which were truncated at 1 Gyr to reduce the run-time, but otherwise followed the BPASS v2.2 prescription described by Stanway & Eldridge (2018). In particular, we assumed that the binary population parameters as a function of mass were independent of metallicity, and remained fixed as the initial mass function varied. Here we reverse these assumptions: we fix the initial mass function and instead question the impacts of the binary population parameters and the uncertainties associated with them. Previous work in this field has been limited and has either considered simple uniform binary fractions with mass (Dorn-Wallenstein & Levesque 2018) or has made use of rapid population synthesis methods (e.g. de Mink & Belczynski 2015) and has considered primarily resolved stellar populations or stellar transients.
In this paper we quantify uncertainties on key aspects of a binary stellar population and spectral synthesis model grid which arise directly as a result of uncertainties on the input binary population parameters. In section 2 we vary the initial binary parameters used as an input by resampling the empirical constraints given their known uncertainty distribution. In section 3 we evaluate the impact of such variation on the output properties of the integrated stellar light, considering ionizing photon production, photometric colour and the ultraviolet-optical spectrum. In section 4 we set aside the observational constraints on binary fraction and vary this distribution through a larger range in order to search for signatures which may be observable, particularly if binary fraction evolves towards low metallicity. The stellar populations contributing most to the resulting variations and uncertainties are discussed in section 5. Finally we present our conclusions in section 6.

Binary Parameters
The models presented here are calculated as a full population and spectral synthesis as described in Eldridge et al. (2017) and Stanway & Eldridge (2018), although each model is truncated at an age of 1 Gyr. Where a star in a binary undergoes a supernova, we estimate the survival probability of the binary and weight models for the subsequent evolution of its stars assuming the Hobbs et al (2006) neutron star kick distribution, sampled 2000 times (sufficient to populate all but the rarest subsequent evolution pathways; a full BPASS data release model typically uses 100,000 kick iterations). All models presented here are calculated with the default initial mass function for single and primary stars "imf135 300"a broken power law model with a slope of -1.3 for primary stars between 0.1 and 0.5 M ⊙ and a slope of -2.35 between 0.5 M ⊙ and an upper mass cut-off of 300 M ⊙ . Models are evaluated at three selected metallicities: Z=0.002, Z=0.010 and Z=0.020 (where Z is the metallicity mass fraction, and Z=0.020 is conventionally Solar metallicity). For the sake of brevity, these are labelled as "z002", "z010" and "z020" respectively where appropriate. The same binary distribution parameters are taken to apply at all metallicities. The binary populations implemented by v2.2 of the BPASS models are defined through the application of empirical constraints on five aspects of the population statistics, and these in turn are constrained by 65 individual quantities, each derived from a meta-analysis of the observed local stellar binary population and associated with an uncertainty by Moe & Di Stefano (2017, hereafter MS17). These comprise the binary fraction as a function of primary mass (five empirical constraints), the companion probability as a function of period (three relevant period ranges, in five mass bins: 15 empirical constraints), and the dependence of companion mass ratio on mass and period, modelled as a broken power law (giving two power law slopes, each in three relevant period ranges, in five mass bins: 30 empirical constraints). A fourth period bin exists in the MS17 data, but stars in this regime would be considered single in BPASS, and this data is only utilised if a functional form is constructed as a function of period. We also utilise an additional 15 constraints on excess twin fraction (i.e. the number of systems with a mass ratio near unity, as a function of primary mass and period) although these are not considered further here, due to the limitations of BPASS stellar evolution models in precise handling of equal mass binaries.
The current BPASS v2.2 implementation uses the bestfit values from table 13 of MS17. Specifically: (i) We implement a binary fraction for the purposes of BPASS as the complement of the single star fraction (i.e. 1 − f sin ) at the median masses of the five stellar mass bins considered by MS17. Masses are linearly interpolated between these fixed mass points. The value given at 16 M ⊙ is loosely extrapolated to unity at 50 M ⊙ and applied at all primary masses above this, while that at 1 M ⊙ is extrapolated to lower masses (with the binary fraction constrained to lie between zero and unity in all cases).
(ii) We implement a period distribution in each MS17 mass bin given by step-wise interpolation between the frequencies given at log (P/days) = 1, 3 and 5. A second linear interpolation is then performed between adjacent MS17 mass bins at fixed period, in order to populate a BPASS mass-period grid.
(iii) We implement a mass ratio distribution by linearly interpolating the upper and lower mass-ratio power-law indices separately between fixed empirical data points. Again interpolation is performed initially in period at each MS17 mass bin and then between mass bin centres to populate the full BPASS grid in period and mass.
(iv) Similarly the fraction of twin systems is interpolated first as a function of period and then of mass, and an appropriate additional weighting is added to our highest massratio bin (q = 0.9).
Care is taken at each stage to ensure that the initial period and mass ratio distributions at each mass sum to unity, and that the weighting of each primary or single star is drawn from the underlying initial mass function, together with the binary fraction for its mass.

Uncertainties on inputs
The current population synthesis models do not take into account the uncertainties on these measurements, but simply apply the best fit values reported by MS17. As Figs 1 and 2 demonstrate, the uncertainties on the MS17 best-fit binary population parameters can be considerable. For example, the companion frequency for O-type stars with log(P/days)=3 is 0.32±0.11, while the power law slopes on the mass ratio distributions have typical uncertainties of ±0.4 or more. BPASS v2.2, as described above, is also restricted to linear interpolation between fixed mass and period points, rather than considering any functional relationship between parameters. Moe & Di Stefano (2017) suggested that the fraction of multiples at each log(initial primary mass) could be fit with a linear function, while the period distribution at each mass is well approximated by a normal distribution in log(Period). Fitting such functions and implementing them directly as a function of mass and period mitigates the scatter associated with each individual measurement and would remove some of the unphysically abrupt changes in binary parameters which result from our interpolation, but rely on these being the correct functional forms. Given the computational burden of each synthesis model, undertaking the tens of thousands of iterations likely to be required by a Bayesian MCMC approach to calculate the posterior probability distribution of all modelled values remains unfeasible. This is particularly true given how little is known about the prior probability distributions on each possible input parameter. Instead we apply a bootstrap resampling of the empirical binary distribution parameters in order to establish the likely impact of their uncertainties.

Calculating Perturbed Models
We evaluate an extensive set of variant stellar population and spectral synthesis models, each with a slightly different set of input binary population parameters. To define these models, perturbations are randomly drawn from a Gaussian distribution with a width scaled to match the uncertainty reported on the parameter by MS17. These values represent a combination of Poisson detection statistics, measurement uncertainties, fitting uncertainties and systematic uncertainties, reported by MS17 as a symmetric 1 σ error range on each parameter. The perturbations are applied to the best fit parameters and a new synthesis performed, either using the existing interpolation scheme as described above, or by fitting a functional form to the newly-perturbed parameters. We vary categories of parameters separately (leaving remaining parameters untouched at the BPASS v2.2 defaults unless otherwise stated), before performing a simultaneous perturbation to all the input parameters.
Variant models are calculated as follows: (i) We perturb the binary fraction as a function of initial mass only. Mass ratio and period distributions remain fixed. A linear fit is used to binary fraction versus log(primary mass). 120 iterations are calculated per metallicity. The resultant binary fraction distributions are shown in Fig. 1. (ii) We perturb the period distribution simultaneously in  all five mass bins. Binary fraction and mass ratio distributions are fixed.
Step-wise linear interpolation is used between mass and period bins. 120 iterations are performed per metallicity. The resultant binary period distributions at three masses are shown in Fig. 2 (upper panels).
(iii) We perturb the period distribution simultaneously in all five mass bins. Binary fraction and mass ratio distribution remain fixed. A Gaussian functional form is fit to perturbed values of log(P) at each mass, where these are linearly interpolated between mass bins as necessary. 120 iterations are calculated per metallicity. The resultant binary period distributions at three masses are shown in Fig. 2 (lower panels).
(iv) We perturbing the two mass ratio power law indices simultaneously in all five mass bins, three period bins and two mass ratio bins. Binary fraction and period distribution remain fixed, with 120 iterations per metallicity. The resultant binary mass ratio distributions at two combinations of mass and period are shown in Fig. 3.
(v) We simultaneously perturb all of the above. 300 random draws are taken from the uncertainty distributions on each empirical parameter, at each metallicity. Input distributions are calculated from the perturbed parameters using a linear fit to the binary fraction with log(primary initial mass), a Gaussian fit to the distribution in log(Period) at each mass, and the perturbed broken power law for the mass ratio distribution at each mass and period. Uncertainties are assumed to be uncorrelated.
As the figures demonstrate, the observational uncertainties in binary fraction as a function of mass are relatively small, allowing for little spread in the models except at the high and low mass limits of the observational data. By contrast, the uncertainties on the period distribution as a function of mass permit substantial variation, particularly when a step-wise linear interpolation method is employed. This relatively crude approach permits fractions of close, high mass stars ranging from near zero to twice the canonical fraction, with potentially large impacts on the output population. If a Gaussian fit is used instead, the more extreme models become disfavoured, and the evolution in period far smoother and likely more physical. Uncertainties on the mass ratio distribution are also large, particularly for mass ratios approaching zero -the distributions for mass ratios q > 0.3 are generally well constrained, with a slowly declining frequency as the mass ratio approaches unity. By contrast those at q < 0.3 show substantial variation permitted by current observational uncertainties, ranging from a decline with increasing mass ratio to a rapid growth.
We note that the uncertainties on binary parameters are likely correlated in reality (i.e. a lower binary fraction at a given mass may be degenerate with a period distribution biased towards larger separations). Thus the estimates on output scatter from variant sample (v) represent a worstcase scenario given current observational constraints.
Outputs of each variant model include stellar type number counts, an integrated stellar spectrum, and various parameters derived from these, including the ionizing photon flux and colours. Models output SSPs at intervals in log(age/years)=0.1 between 10 6 and 10 9 years, these can be combined to produce complex CSPs given a star formation history. Older populations are not considered here.

Ionizing photon output and UV continuum
A key output from an SPS model is the ionizing photon production rate, and its efficiency compared to continuum  Figure 3. Examples of interpolation of MS17 initial mass ratio distributions onto the BPASS mass grid. MS17 provide two power law slopes of the form p q ∝ q γ at each of five mass bins and four initial periods. The thick line is the result of the BPASS v2.2 interpolation algorithm. Faint lines result from perturbing each of the 20 initial power law slopes (simultaneously) before interpolation. In each case, the sum of the mass ratio frequencies is normalised to unity. emission longwards of the 912Å Lyman break. For any stellar population which incorporates young stars, the ionizing spectrum is likely to be reprocessed by nebular gas in the interstellar medium, and heavily modify the resulting spectrum. A full calculation of the observed spectrum requires the stellar spectra derived from SPS codes to be passed through an independent radiative transfer code. We do not undertake this here, but instead consider the ionizing photon production directly as a function of metallicity for a case in which a stellar population is undergoing star formation at a constant rate of 1 M ⊙ yr −1 . Since it takes time for the rate of stellar birth to balance that of stellar death, and so for the output spectrum to stabilise, we evaluate the photon production at a fixed epoch of 100 Myr after the onset of star formation.
The variation in the ionizing photon output in this continuous star formation case is illustrated in Fig. 4 and tabulated in Table 1 for variants (i) to (iv) above. In each case, the scatter introduced by variation in the binary inputs is well below the differences between the three metallicities considered here, and typically below 0.05 dex (except in the case of set (ii)). The effect of switching from step-wise interpolation of the binary fraction between mass bins, as implemented in BPASS v2.2, to a linear fit to the fraction as a function of mass is a small, metallicity-dependent systematic shift of [0.04, 0.02, 0.02] dex at Z=0.002, 0.010 and 0.020 respectively. This is larger than the random error introduced by perturbations on the binary fractions themselves, which are <0.01 dex at all metallicities. The uncertainties introduced by varying the binary period distribution (while holding other parameters fixed) depend on the method of interpolation. The step-wise linear interpolation approach currently adopted in the BPASS v2.2 model grids leads to both a systematic offset in the ionizing photon rate at each metallicity and an increased standard deviation (although this is driven primarily by large outliers which have either very large or very small numbers of interacting binaries at the highest masses), both systematic and random uncertainties are comparable in scale at ∼0.1 dex. By contrast, a log-gaussian interpolation approach leads to a much smaller scatter in ionizing photon rates with perturbation of the input parameters, resulting in only <0.015 dex variation, while also producing a more physically plausible input period distribution.
By contrast, perturbing the power law indices of the mass ratio distribution results in a broader distribution of ionizing photon rates than any of the other parameters, with no clear systematic offset but standard deviations of [0.041, 0.026, 0.025] dex at Z=0.002, 0.010 and 0.020 respectively. While not obviously related, the large scatter may result from the interaction of this parameter with the IMF sampling in a binary population synthesis (see Appendix A in the supplementary information for a more detailed discussion of this point). During the population synthesis, primary stars are weighted by mass according to the selected IMF. The weighting for each primary mass bin is then subdivided between single and binary models, and the binary weighting is further divided between models with a range of periods and mass ratios. The mass of secondary companions are therefore not included in the initial IMF weighting, but are included in the final normalisation of the weighting scheme to a total of 10 6 M ⊙ . Increasing the mean companion mass for low mass stars, or decreasing that for high mass stars, would have the effect of decreasing the total mass fraction occupied by stars above any given mass threshold, and thus decreasing the ionizing flux for the integrated population.
The time dependence of the uncertainty on the ionizing  photon production rate, the FUV continuum luminosity and the FUV-NUV spectral slope for each of variants are shown in Fig. 5, for the case on a single-aged, evolving starburst population (i.e. an SSP). Since the constant star formation case is effectively an integral of the emission from populations less than 100 Myr in age, this allows us to examine at which ages the populations are most sensitive to binary parameters, and hence later consider the populations involved.
The uncertainty on both the UV continuum spectral slope power law index β UV < 0.02 and that on the rest-UV luminosity L UV < 0.02 dex for a simple stellar population at all three metallicities and ages less than about 500 Myr (after which the UV flux is intrinsically very faint and the uncertainty on these parameters grows significantly). Before this point, neither L UV nor β UV show significant time or metallicity evolution in their uncertainties, except in the case of a step-wise interpolation of the period distribution. Additional uncertainty occurs at early times in this case, because the prescription allows for both zero close binary fractions and for very high numbers of close binaries -a range which is more restricted in the other variants -and hence there is   (v)). The panels show the variation in ionizing photon production rate (left) and the standard deviation in photon production, UV spectral slope and UV luminosity (right) as in Figs. 4 and 5. also large scatter in the number of very massive, luminous, blue stars which interact or merge early in their evolution.
Ionizing photon production rate shows more time dependence in its scatter between model variants. In all cases the lowest metallicity we consider, Z=0.002 (0.1 Z ⊙ ), demonstrates the largest scatter in ionizing photon rate. This output is sensitive both to the fraction of massive stars at early times (log(age/yrs)<6.5), and the fraction of interacting binaries or their products at later times.
As was the case for L UV and β UV , the ionizing photon rate of variant (ii) shows large early time scatter due to large uncertainty in the close binary population of massive stars, which may not be physical. However both period interpolation variants also show a relatively large variance at late log(age/yrs)∼8 at all three metallicities. This occurs because moving the stars from closer initial orbits to wider orbits, and vice versa, changes the fraction of relatively numerous intermediate mass stars which are likely to interact during their main sequence lifetimes or soon thereafter. A subset of these will either be stripped (primaries), or spun-up and rejuvenated (secondaries), and will become hotter and bluer, producing more ionizing photons. Interestingly, the uncertainty on the ultraviolet continuum does not track the ionizing photon production in this case, suggesting that most of the variation is occurring blueward of the Lyman limit, and thus an important role for stripped-envelope stars.
For variant (iv), the mass-ratio distribution, there is a short-lived rise in the uncertainty on ionizing photon production rate at the lowest metallicity and ages of log(age/yrs)∼7.
The same results when all parameter distributions are sampled simultaneously and independently (i.e. variant (v), the maximum possible scatter permitted by current observational constraints in the local Universe) are shown in Fig. 6. Unsurprisingly the uncertainties on the observable parameters are slightly larger than in the variants where all but one distribution is held fixed. The time evolution of those uncertainties also show contributions from the combination of the same effects shown above, including a short-lived rise in uncertainty on the ionizing flux around log(age/years)∼ 7 and a longer duration epoch of uncertainty at log(age/years)∼ 8. The uncertainty on the ionizing photon flux varies from 0.02 dex to almost 0.12 dex over the first Gyr for a simple stellar population, depending on age and metallicity, but this variation is reduced in a continuously star forming population, where uncertainties on ionizing photon production are 0.03-0.04 dex. The uncertainties on both L UV nor β UV remain relatively time and metallicity insensitive implying a typical uncertainty on these parameters of 0.02-0.03 (dex for L UV and linearly in β UV ) for binary SSPs.

UV and Optical Absorption Lines
The sensitivity of the synthetic population spectrum to variation in all binary parameters (i.e. variant (v)) is illustrated in Fig. 7, where the one sigma flux uncertainty at each metallicity is shown as a function of wavelength for both a simple stellar population at 10 Myr and for a population which has been constantly forming stars at 1 M ⊙ yr −1 for 100 Myr. The largest uncertainties at all metallicities arise in the rest-frame far-ultraviolet, shortwards of the Lyman limit at 912Å, and again are largest at low metallicity, reaching almost 30 per cent at Z=0.002 and typically 20 per cent at higher metallicities for a population with constant star formation. At wavelengths longwards of 1200Å, perturbations in the binary parameters tend to introduce a systematic offset in the normalisation of the spectrum at any wavelength, with the uncertainty in flux per angstrom varying only slowly and of order 10-20 per cent at all three metallicities.
As Fig. 8 demonstrates, absorption lines and continuum strength have comparable uncertainties. An exception to this is seen in the strength of certain absorption lines in the ultraviolet, including Si II 1261Å, C I * 1657Å and Al II 1671Å at the higher metallicities in a composite stellar population, and He II 1640Å and C II 1335Å at the lower metallicity for a young starburst. Each of these shows a slight excess uncertainty relative to the continuum, although this is unlikely to be useful as a diagnostic of binary parameters: the variation is still small (e.g. 12 rather than 10 per cent scatter relative to the fiducial population), these lines are intrinsically very weak, and they are subject themselves to additional uncertainty due to factors such as α-element enhancement and other abundance ratio or formation timescale effects. They are also likely to vary between different atmosphere models in use in the population synthesis (e.g. a different choice of input atmosphere models could plausibly change minor lines such as these) and so are unlikely to prove diagnostic given observational uncertainties. Small differences in the age or metal abundance of the stellar population will overwhelm any variation due to binary fraction. A simple stellar population (shown at 10 Myr in Fig. 7 and 8) exhibits larger typical uncertainties than a composite population, in which an effective time-average smooths over extremes in any given binary parameter set.
There is no single spectral feature longwards of 912Å which is clearly and strongly diagnostic of variations in the binary parameter distributions, given the current level of uncertainty on binary parameters in the local Universe.

Colours
Unsurprisingly, the small scatter in the continuum as shown in section 3.2 results in a similarly small scatter in the predicted rest-frame colours of the integrated stellar light. This is demonstrated in Fig. 9 which gives the distribution in V − I vs I − K colours for a sample of model variants, when all binary parameters are varied (i.e. variant (v)). The time evolution is given both for a simple stellar population as it ages from 1 Myr to 1 Gyr, and for a population with constant star formation as a function of time since onset of star formation. Note that the reddening effects of nebular continuum emission are not considered in this direct comparison. At most population ages, the uncertainty on rest-frame optical and near-infrared colour is <0.05 mag, significantly below the uncertainty that would be introduced by assumptions regarding stellar population age, star formation history and metallicity. However at stellar population ages between 10 and 100 Myr, binary products have a significant effect, causing the population to turn back towards bluer colours. The effect of varying the binary distribution parameters is to cause this blue loop to begin either slightly earlier (corresponding to a large interacting binary fraction at higher mass) or later (if interactions are less common at the highest masses). This leads to degeneracy between model age and the input binary parameter distributions, such that at Solar metallicity, this colour combination cannot distinguish the two. At significantly sub-Solar metallicities, the scatter in population colour is instead more pronounced at early ages, where the most massive stars dominate.
In the constant star formation case, binary parameter uncertainties result in a ∼0.015 mag variation in V-I colour at most ages, with the blue loop undertaken by the simple stellar populations leading to a larger colour range (∼0.2 mag of scatter) for high metallicity populations at ages between 10 and 100 Myr.

SIGNATURES OF THE BINARY FRACTION
In the above tests, we have assumed that the empirical constraints on the binary population parameters derived from the local Universe apply at all metallicities. There are rea-sons to question that assumption. Metallicity likely affects the cooling and fragmentation of star forming molecular clouds. While direct evidence for the resultant binary fractions is difficult to obtain, there is now growing evidence that the close binary fraction in Solar-mass stars is higher at low metallicities (Moe et al. 2019;El-Badry & Rix 2019;Price-Whelan et al. 2020), which in turn would suggest a flatter relation between binary fraction and primary mass (or alternatively a period distribution more skewed towards close binaries at low mass).
Given that the empirical uncertainty on the extant binary population parameters is unlikely to lead to a clear observational signature, it may be instructive to investigate whether such a deviation from the local parameters (for example in very low metallicity environments, or in the distant Universe) might be detectable. This question has been explored in past work by Dorn-Wallenstein & Levesque (2018) who used a grid of models in which the binary fraction var- ied between 0 and 1, but in which -crucially -the binary fraction was independent of mass (i.e. a weighted combination of single star models and models in which every star was in a binary). Their analysis focused on observables derived from resolved stellar populations, identifying the Red Supergiant (RSG) star population and its abundance relative to stripped stars as diagnostic of massive star multiplicity. Here we extend a similar analysis to consider a range of mass-dependent binary fractions which differ significantly from the empirical constraints considered up to this point, and focus instead on plausible signatures in the integrated light of unresolved stellar populations.

Extended Model Grid
The parameter distribution used in BPASS v2.2 and shown in Fig. 1 is characterised by a binary fraction that increases with the logarithm of stellar mass, from approximately 40 per cent at 1 M ⊙ to 100 per cent at approximately 100 M ⊙ . To evaluate the effects of a substantial change to this distribution (i.e. one far exceeding the current empirical constraints) we evaluate 20 additional models at each metallicity with artificially modified binary fractions, as shown in Fig. 10.
These models each assume a linear relation between multiple (in this case binary) fraction and log(primary mass) and fall into two categories: ten models in which the low stellar mass binary fraction rises (set 1) and ten models in which the high mass binary fraction falls (set 2).
We make no claims for the physical relevance of these models. While observational evidence may favour the set 1 scenario at low metallicities, as discussed above, the constraints are weak and it is difficult to be certain which regime Each line indicates a model binary fraction distribution which either raises the binary fraction at low stellar mass (set 1, dashed lines) or lowers it at high mass (set 2, solid lines). Data points are drawn from MS17 (as in Fig. 1) and the thick line indicates the fiducial model applied in BPASS v2.2. may dominate in the distant Universe. Other binary parameters (i.e. binary mass ratio and initial period distributions) are held fixed at the fiducial values.

Stellar Continuum
The ultraviolet and optical spectra resulting from these artificial variant models are shown in Fig. 11 Figure 11. The effect of changing the binary fraction distribution on the ultraviolet-optical spectrum at Z=0.002. Left hand panels show set 2 (decreasing high mass binary fraction) and right hand panels set 1 (increasing low mass binary fraction) as defined in Fig. 10. The top row shows results for a 10 Myr old SSP, while the bottom row shows the case for constant star formation. All models are normalised by the fiducial BPASS v2.2 model. Models are colour coded as in Fig. 10. shown at the lowest metallicity (Z=0.002) and each model is normalised by the flux as a function of wavelength in the fiducial BPASS v2.2.1 model set. The shift from a piecewise interpolation to a linear function for binary fraction as a function of mass is such that none of these models precisely reproduces the BPASS v2.2 models. Nonetheless, clear trends are visible with binary fraction. Models with small high mass binary fractions unsurprisingly produce lower ionizing fluxes than those with large high mass binary fractions, and also show a lower continuum normalisation in the nearultraviolet.
Interestingly, increasing the binary fraction of low mass stars has a similar (although not identical) effect. As before, this likely arises from the application of the initial mass function in BPASS. The broken power-law IMF used uniformly in all these models assigns a weighting to binary systems first by their primary mass, and then subdivides this weighting between variants with different mass ratios and period distributions. If the number of low mass binaries are increased, the weighting for each primary is unchanged, but the overall model normalisation must allocate mass to the (low mass) companions of low mass stars, which would otherwise have been allocated to high mass primaries and their (moderateto-high mass) companions. Since the number of hot, lumi-nous stars falls as a result, the overall model luminosity also drops.
In the optical, the behaviour diverges between the two sets of binary fraction variants. Models in which the low mass binary fraction is varied (set 1) shows the same trend in the optical as the ultraviolet. By contrast models in which the high mass binary fraction is decreased (set 2) show much redder spectra in the optical than those with a large high mass binary fraction. This results from a larger population of massive stars reaching the cool giant branch without being subject to stripping by a binary companion (i.e. a binary dependence in the Wolf-Rayet to Red Supergiant ratio), leading to a much redder integrated light spectrum. In each case, this behaviour occurs in both young simple stellar populations, and in composite stellar populations forming stars at a constant rate, although in the latter case the variations between models are typically smaller.
In addition to considering the continuum luminosity and shape, we can also measure the strengths of ultraviolet emission and absorption lines after normalisation by a pseudo-continuum (obtained by fitting a third-order polynomial to the wavelength range λ = 1000 − 2000Å). After this normalisation, the strengths of emission and absorption features relative to the continuum prove surprisingly insensitive to changes in the binary fraction, with modifications to line equivalent widths typically < 0.01Å (as Fig. 12 demonstrates). This is particularly true in the continuously star-forming composite stellar spectra in which any variation in line strength is averaged out over time. A notable exception occurs only at the lowest metallicity under consideration (Z=0.002) and in the case of the broad, stellarwind-driven He II 1640Å feature which arises predominantly from massive stars with stripped atmospheres. This shows a clear trend in strength with high mass binary fraction (see Fig. 12) but no clear dependence on low mass binary fraction variation.
The origin of this dependence is suggested by Fig. 13 which presents the time evolution of the line equivalent width, evaluated as a summation of the excess flux within ±1000 km s −1 of the line centre. It demonstrates that the dependence in the composite population arises primarily from simple stellar populations at ages between 3 and 20 Myr. Unsurprisingly this is the epoch at which the ultraviolet spectrum is dominated by massive stars (at around 10-30 M ⊙ ) which are approaching the end of their main sequence lifetimes. In the high metallicity cases, these stars lose their hydrogen envelopes primarily as the result of strong radiationdriven winds. By contrast, at low metallicities, the primary mechanism for stripping is through binary mass transfer onto a companion. This both delays the emission of He II line flux and causes it to be dependent on the binary fraction at low metallicity.
It should be noted however that the difference between the largest high mass binary fractions considered here (∼1) and the lowest (∼0.4) changes the line equivalent width from 0.6Å to only 0.8Å for a constantly star forming population, and from 0.9 to 1.6Å for a simple stellar population at peak emission, although the latter is only clearly diagnostic if the age of the population is constrained to around 0.1 dex or better. Degeneracies in interpretation of this line as a binary indicator are likely to be substantial due to star for-mation history and metallicity uncertainties. We also note that the atmosphere models in this very low metallicity, high mass, stripped star regime remain an active area of study (e.g. Hainich et al. 2019;Götberg et al. 2018Götberg et al. , 2019. As a result, while the stellar-wind-driven broad He II is sensitive to massive binaries, it is unlikely to be a clear and unambiguous indicator of their quantitative fractions in realistic observations.

Ionizing fluxes
If the stellar spectrum shows little dependence on binary fraction (over the range considered here) then can we instead identify strong changes in the ionizing spectrum (i.e. shortwards of the Lyman limit at 912Å)? Fig. 11 suggests that the ionizing spectrum varies significantly between models in both set 1 and set 2. However this is primarily in terms of normalisation, just as is the case with the ultraviolet-optical stellar continuum. This means that the overall number of ionizing photons emitted by a population of known age and star formation history will be higher when the high mass stellar binary fraction increases, but the observable UV continuum is also stronger. The ionizing photon production efficiency parameter ξ ion = N ion /L U V can be measured as a ratio between the observable 1500Å continuum and the strength of hydrogen recombination emission lines (as a proxy for ionizing flux). This parameter shows a maximum difference between models in set 2 of < 0.05 dex for a continuously star forming population at Z=0.002 and < 0.03 dex at Z=0.020, as shown in Fig. 14. While slightly larger variation can be seen in simple stellar populations, this is highly time dependent.
At any given age, populations with lower binary fractions amongst high mass stars show lower ionizing flux. The models in set 2 diverge significantly after log(age/years)∼6.8 in a simple stellar population model and thereafter show a ∼ 0.2 dex spread in ξ ion . Interestingly, Fig. 14 demonstrates that binary fraction dominates over metallicity in determining the ionizing photon production efficiency at ages above log(age/years)=7.4 (25 Myr). If simple stellar populations can be identified in this age regime, then ξ ion may prove a useful probe of binary fraction. However, as Fig. 14 also demonstrates, any more complex star formation history drowns out this signal (as would any significant variation in binary mass ratio and period distributions) leading to ambiguity between population age, metallicity and binary fraction.
The helium ionizing photon flux (i.e. integrated light emitted below 227Å) is also dependent on massive binary fraction. The time evolution of this flux, which will power a narrow nebular emission line component overlaying the stellar He II 1640Å line discussed above, is shown in Fig. 15. For a low metallicity composite stellar population with constant star formation, by the age of 100 Myr the ratio of helium ionizing flux to UV continuum, ξ He II , varies by 25 per cent, although the variation is substantially smaller at higher metallicities.

Prospects for Observation
While the He II 1640Å line (and also its counterpart in the optical at 4686Å) emerges as a potential binary indicator  Figure 14. The strength of ionizing photon flux shortwards of 912Å. Flux is shown relative to the 1500Å UV continuum luminosity, i.e. the widely-used ξ ion parameter. Models are colour coded as in Fig. 10. Set 1 (increasing low mass binary fraction) models are essentially indistinguishable from the highest set 2 (reducing high mass binary fraction) models.
both through its broad stellar wind-driven component and the narrow recombination component emitted by ionized nebular gas, the prospects for using this line as a binary fraction indicator remain poor, particularly for sources in the distant Universe, where the metallicities and population ages are expected to favour binary population effects.
The He II emission feature has been identified in distant galaxies ever since the first stacked spectra of galaxies at z ∼ 3 were constructed with sufficient signal to noise (Shapley et al. 2003). As this and studies since have demonstrated, typical ultraviolet-selected star-forming galaxies in the early Universe show evidence for weak, broad He II emission en masse (e.g. Shapley et al. 2003;Steidel et al. 2016;Cullen et al. 2019). A handful of objects in the distant Universe have also been identified as individual emitters of strong but usually narrow (nebular) He II emission (e.g. Erb et al. 2010;Sobral et al. 2019;Bayliss et al. 2014). Similar strong emission features are also seen in both the ul-traviolet and optical He II lines in extreme, low metallicity starbursts (e.g. Berg et al. 2019;Kehrig et al. 2018).
Such strong He II emission is by no means universal even in young starbursts: Saxena et al. (2019), working with extremely deep spectra from the VANDELS survey suggest that only about 20 galaxies show strong (W(He II)> 0.5Å) narrow emission, and 6 broad emission, amongst a much larger sample (899, < z >∼ 3.5) of galaxies with similar observed properties, masses and inferred star formation rates. Typical uncertainties on the line flux are at the 25 to 40 per cent level, even for these detections. Similarly, Nanayakkara et al. (2019) identify a small sample of 13 individually detected sources amongst almost eight hundred galaxies observed in deep (>9 hr) VLT/MUSE spectroscopic fields.
This suggests that such emission is highly stochastic, occurring in short-lived bursts, more consistent with the ∼10 Myr timescale for emission from simple stellar populations shown in figure 13 (left), rather than a slow but contin- . Flux is shown relative to the H I ionizing flux (below 912Å) and relative to the 1500Å UV continuum luminosity, i.e. the He II equivalent of the widely-used ξ ion parameter. Models are colour coded as in Fig. 10. Set 1 (increasing low mass binary fraction) models are essentially indistinguishable from the highest set 2 (reducing high mass binary fraction) models.
uous star formation episode. Such stochasticity in the star formation history will make the interpretation of any given line detection as a binary indicator challenging. The time variation in the SSP line strength on 10 Myr timescales exceeds the difference in line strength introduced by changing the binary fraction, even if the line is measured at sufficient signal-to-noise to distinguish between binary fraction models (i.e. uncertainties <<10 per cent, beyond most current observations in the distant Universe, and challenging given the faintness of the lines even at low redshift).
A further complication is hinted at by the strengths of the observed emission lines. In extreme cases, these reach equivalent widths of W(He II)> 10Å in the rest-frame. As we discussed in Stanway & Eldridge (2019), such lines are beyond the capacity of any current stellar population synthesis model, including BPASS, to reproduce and imply the presence of a source of ionizing photons other than stellar photospheres. This, together with the strong metallicity dependence observed in X-ray binary models (e.g. Fornasini et al. 2019), has motivated ongoing work by us and others to incorporate XRB populations in stellar population synthesis (Eldridge, Stanway et al, in prep). Since emission spectra of such sources are not implemented in the current version of BPASS (or any other code) we caution against overinterpretation of the He II feature.
A full analysis of nebular radiative transfer from these ionizing populations lies well beyond the scope of this paper. However, we note that the comparison of flux below 912Å (hydrogen ionizing) and that below 227Å (helium ionizing) may fail to detect changes in the ionizing spectrum shape between these limits. Other nebular emission lines may be sensitive to this behaviour, for example the ratios between the semi-forbidden C III] doublet, the He II line, Lyman-α and [O III] are diagnostic of the relative flux below the respective critical wavelengths for these features. All of these features save Lyman-α, itself a resonantly scattered line and difficult to interpret, are weak in the spectrum of typical high redshift stellar populations. Nonetheless spectroscopic surveys such as VANDELS (McLure et al. 2018) are now achieving signal to noise ratios ∼ 20 per resolution element. The prospect of NIRSPEC on the James Webb Space Telescope extending similar high precision measurements to larger samples and into the rest-frame optical, suggests that detailed reconstruction of the ionizing spectrum (and the role of binaries in shaping it) may be possible in the near future.
Potentially more promising as a high mass binary fraction diagnostic is the rest-frame ultraviolet through optical colour, identified in Fig. 11. In Fig. 16 we estimate photometric colours for the variant models of section 4. Given that the details of filter profile will depend on instrument and target redshift, we model the observations as simple top-hat filters with a width of 100Å, centred at rest-frame wavelengths of 1500, 4000 and 7000Å. At z = 5, such a combination might be approximated by combining NIRCAM filters F090W (or HST/WFC3 F606W), F162M and F300M, while at z = 3 the equivalent combination would be F090W, F250M and F430M (note: broader filters will typically reduce photometric spread between models).
The top left panel of Fig. 16 suggests that this colour combination (or a similar filter selection) is indeed a potential diagnostic of both metallicity (at the level of precision explored here) and high mass binary fraction for continuously star forming populations, if ∼0.01 mag photometric precision can be achieved. Distinguishing between low mass binary fractions requires significantly higher (millimag) precision. Unfortunately, as the remaining panels of Fig. 16 make plain, such a simple interpretation will be compromised by any young stellar populations or complex star formation history. The predicted colours of simple stellar populations at different metallicities overlap over much of their parameter space, introducing a stellar age-metallicity-binary fraction degeneracy.
We also note that the colours presented in Fig. 16 are for the stellar continuum only and do not include the reddening effects of the nebular gas continuum. The effects of nebular gas reprocessing are demonstrated (for continuous star formation and a single selected set of gas parametersn e = 10 2 cm −3 , spherical, complete covering, fixed geometry, dust grain depletion 1 ) in Fig. 17. This acts as a reddening term to reduce the spread in predicted photometric colour. The models still form a grid which is potentially diagnostic of metallicity and binary fraction given ∼0.01 mag photometric precision. However this is challenging to achieve in UV to optical colours due to systematic calibration uncertainties (more typically a few percent), corrections for dust extinction (with uncertainties often >0.1 mag in the ultraviolet), and the degeneracies introduced by star formation history, with additional possible uncertainties due to the gas parameters also introduced.
Upcoming observations, such as the planned NIRCam-NIRSpec Galaxy Assembly Survey and The Cosmic Evolution Early Release Science (CEERS) Survey, will plausibly be able to eliminate some of the more extreme potential variants in binary fraction at moderate redshift, through a combination of extremely high signal to noise and uniformly calibrated imaging, the use of multiple lines to determine dust extinction, and high precision spectroscopic data. Nonetheless the difficulty of distinguishing systems with high binary fractions from those with lower metallicities or slightly younger stellar populations will not be easily solved. It remains to be seen whether combining full spectral-fitting (including nebular gas radiative transfer) with photometric spectral energy distributions will be sufficient to break these degeneracies. The analysis presented here suggests that doing so may be challenging, but further analysis may point to more hopeful strategies in future.

THE EFFECTS OF UNCERTAINTIES ON STELLAR POPULATIONS
The time evolution of the variant models explored above has already provided hints regarding the components of the stellar populations most vulnerable to uncertainties in the input distributions. However, it is also informative to consider the variation in these populations directly. Fig. 18 shows the effect of varying binary population parameters on the numbers of different stellar types as a function of stellar population age for a simple stellar population of 10 6 M ⊙ at Z=0.002. Number counts are presented for the dominant hot star types responsible for ionizing radiation in star-forming populations: type O main sequence stars, partially-stripped WNH stars and stripped-envelope Wolf-Rayet (WR) stars in the WC and WN sub-classes. All stellar models in BPASS are classified by effective temperature and surface composition (see Eldridge et al. 2017). Behaviour at other metallicities is qualitatively similar, although evolution tends to be a little faster as metallicity increases (due to stronger winds) and the number of luminous WR stars produced is lower at the highest metallicities.
The top panels of Fig. 18 illustrate the ±1 σ range of number counts which are permitted by current observational constraints on binary parameters (i.e. the variant (v) models of section 2.3), while the lower panels show explicitly the dependence on binary fraction while keeping all other parameters fixed (the models of section 4). In each case, the total population by stellar type is shown as a function of age.
1 Modelled with CLOUDY c17.01 (Ferland et al. 2017) It should be noted that stripped and partially-stripped populations are almost entirely comprised of luminous, massive stars (log(L/L ⊙ )>4.9) at early ages and slower evolving, often post-interaction, lower mass stars at late times, resulting in the characteristic double-humped time distribution.
As the comparison of this figure with Fig. 15 makes clear, the origin of time (and binary fraction) dependence in He II line strength can be attributed almost entirely to the luminous (log(L/L ⊙ )>4.9) population of classic Wolf-Rayet stars. Indeed, the time evolution of the He II stellar line equivalent width hints that WN stars likely dominate over WC stars (which peak earlier). These are present in the population in a well-defined temporal window between log(age/years)∼6.3 and 7.3, and are strongly dependent on the binary fraction of massive stars at low metallicity, as is also seen in the equivalent width of the He II 1640Å emission line. At higher metallicities, where winds contribute significantly, the dependence on binary fraction is less marked. While other stellar populations may contribute weakly to this line at other times, the behaviour of the He II 1640Å line in these models is clearly dominated by this handful of very massive stars which have lost their envelope.
The hydrogen ionizing photon production efficiency at low metallicities may also receive a significant contribution from this population. As Fig 14 demonstrated, at Z=0.002 the early-time decline in ionizing flux is much gentler than that seen at Z=0.010 and Z=0.020, before abruptly dropping at log(age/years)∼7.3. This would be consistent with Fig 14 primarily tracking the O star population (which dominates in both total luminosity and number counts) but being boosted at low metallicities by the WR population, before these drop sharply away.
In the case of low high-mass binary fractions, the decline in the number of giant stars whose envelopes are stripped through Roche Lobe overflow is compensated for by an increase in the number of red supergiants (inferred types K and M with log(L/L ⊙ )> 4.9, labelled RSG on Fig. 18). These stars retain their envelopes through their expansion onto the giant branch and cool significantly, resulting in the redder colours seen in the optical in Figs. 11 and 16. In common with Dorn-Wallenstein & Levesque (2018), we find that the WR-to-RSG ratio is diagnostic of high mass binary fraction in resolved populations, and note that the relative fraction of type II core-collapse supernovae to stripped-envelope type Ib/c supernovae also tracks this ratio.
The role of the lower luminosity hot population which arises primarily from binary interactions is less clear cut. Comparison of the timescales in Fig. 18 with Fig. 5 suggests that the binary initial period distribution affects these populations most strongly. The ionizing photon production rate at early times is always dominated by the (relatively binaryinsensitive) short-lived massive stars, but at late times it is the low luminosity population that dominates. In particular, the large population of partially stripped WNH stars are sensitive to binary fraction and show a marked increase from low numbers at early times, to a peak in population numbers at log(age/years)∼ 8.5. While fully-stripped low luminosity helium stars also grow in number until their population is comparable to that of WNH stars, they are more peaked in time, falling away sharply after log(age/years)∼ 8.1. By contrast the uncertainty in ionizing photon rate from the populations shown in Fig. 5 varies more slowly, suggesting .020 (bottom right). Selected age steps are annotated with log(age/years). In all cases, the 4000-7000Å colour for set 2 is higher than that for set 1. No nebular gas reprocessing is considered that it is the partially-rather than fully-stripped stars which are dominating this uncertainty, and thus dominating the uncertainty in the ionizing photon rate from a continuously star forming population.

CONCLUSIONS
In this paper we have explored the impact of empirical uncertainties on the parameter distributions in binary fraction as a function of initial primary mass, initial orbital period and initial mass ratio. We have sampled the range of values permitted by current empirical constraints on these distributions and also explored the effect of modifying the binary fraction as a function of primary mass into regimes outside of the current constraints. We have focused on unresolved stellar populations, since the majority of circumstances in which the products of SPS models are used do not permit the identification of individual stars, but rather considered the integrated light of stellar populations or galaxies. Our principle conclusions can be summarised as follows: (i) The current constraints on binary distribution param-eters are relatively tight on the binary fraction as a function of mass and period, but allow more variation in mass ratio. Fitting a functional form to the empirical constraints as a function of mass or period prevents unphysical distributions and reduces the impact of uncertainties.
(ii) The inferred uncertainties on the ionizing photon flux from a young stellar population are < 0.1 dex, while those on the UV spectral slope and 1500Å continuum luminosity are <0.03 and <0.03 dex respectively. Low metallicity models show a slightly increased sensitivity to the binary parameters.
(iii) Uncertainty on the ultraviolet-optical continuum luminosity is typically 10-15 per cent, with the uncertainties acting primarily on the continuum normalisation rather than slope or emission/absorption lines.
(iv) Photometric colours are typically uncertain at the < 0.05 mag level, except in a narrow range of ages.
(v) The He II emission line at 1640Å shows increased sensitivity to the binary fraction relative to the continuum, both due to stellar emission and nebular emission from ionized gas. However this emission is highly time dependent  Figure 17. The synthetic photometric colours of populations for stellar populations with varying binary fraction, colour coded as in Fig. 10. Colours are given in AB magnitudes, and shown for all three metallicities given a continuously star forming population at 100 Myr. The reprocessing of the stellar spectrum by nebular gas is included, as described in section 4.4. and small uncertainties in population age or metallicity will compromise attempts to interpret measurements.
(vi) The 1500-4000Å colour relative to the 4000-7000Å colour of the stellar continuum is also mildly sensitive to extreme binary fraction distributions. Again, age and metallicity uncertainties will complicate interpretation, as will uncertainties on dust and nebular gas reprocessing of the stellar continuum.
(vii) Most of the variation seen with binary parameters arises from changes in the number of stripped-envelope stars (particularly luminous Wolf-Rayet stars) in the population, with contributions also arising from partially-stripped WNH stars.
While the absence of any clear and unambiguous binary parameter indicators in the integrated light of stellar populations is disappointing, the very small uncertainties permitted by current observational constraints on the nearby stellar population are reassuring: interpretation of these populations are unlikely to be compromised by binary parameter uncertainties, and is only mildly sensitive to massive star binary fractions ranging between 40 and 100 per cent. It is notable however that in every parameter, the lowest metallicity models show the largest sensitivity to the input binary parameter distributions. Since it is in low metallicity environments that binary fractions are also expected to vary, and for which the stellar atmospheres of hot stars are also most uncertain, these results that caution should be used in their interpretation, and are a reminder that efforts should be intensified to fully characterise very low metallicity populations wherever possible.  Figure 18. The distribution of stars between stellar types as a function of age, for SSP models at Z = 0.002. Top: mean ±σ uncertainty range in number counts permitted by current uncertainties on binary parameters; Bottom: varying binary fraction, colour coded as in Fig. 10. A1   APPENDIX A: THE INTERACTION OF  BINARY PARAMETERS AND IMF (ONLINE  ONLY) As discussed in detail by , the initial mass function of a stellar population presents significant difficulty, both conceptually and in terms of its observational determination. Increasing evidence suggests that there may be no such thing as a universal initial mass function, and those empirical functions determined from observations rely on stellar evolution and population synthesis models to infer the mass and age of individual stars and estimate the changes in the population which have occurred since the stars reached the zero-age main sequence. Such inference is typically made with single star population synthesis models, which may introduce systematic offsets, and suffer from additional uncertainties in mass and age determination from unresolved binaries, or high mass-ratio systems in which the secondary is not observed.

Uncertainties in Binary Pop Synth
The application of an initial mass function also presents certain complications in binary population synthesis. Since the companion probability and mass ratio distribution is strongly dependent on primary mass, altering these parameters will result in higher or lower numbers of comparatively low mass companions for each star, altering the overall mass function. Thus each perturbation in binary parameters is also an effective perturbation on the IMF power law slope.
Kroupa and coworkers (e.g. Kroupa et al. 1993;Kroupa 2001;Kroupa & Weidner 2003;Kroupa et al. 2013;Kroupa & Jerabkova 2018) have extensively discussed the effects of stellar multiplicity on initial mass function estimates. They advocate the random pairing of stars drawn from a canonical single star IMF, with a period assigned statistically to each binary, in a so called Birth Binary population. However such an approach does not permit the mass ratio distribution to be controlled or varied, and implicitly assumes that this distribution is determined only by the initial mass function (i.e. low mass companions are more common as a result of random pairings, simply because low mass stars are more common). By contrast, other binary population synthesis models including BSE , binary c  and StarTrack Banerjee et al. 2019) treat the mass function as applying only to single and primary stars, with secondary stars lying outside of the initial mass function weighting and their mass determined instead by the mass ratio distribution.
The implementation of initial mass function in BPASS follows the latter approach. The IMF is presumed to apply to the sum of primary stars and single stars (which together dominate the mass), but not to the secondary stars. As discussed in section 3.1, during the population synthesis, primary stars are weighted by mass according to the selected IMF. The weighting for each primary mass bin is then subdivided between single and binary models, and the binary weighting is further subdivided between models with a range of periods and mass ratios. The mass of secondary companions are therefore not included in that initial IMF weighting, but are included in the final normalisation of the weighting scheme to a total of 10 6 M ⊙ . The models are sampled statistically rather than stochastically, and as such may be inappropriate in situations where a low starburst mass results  Figure A1. The role of single, primary and secondary stars in building the effective IMF using the canonical prescription for BPASS v2.2. The nominal IMF (applied to single stars and primaries) has a slope of -2.350 between 0.5 and 300 M ⊙ . Secondary stars perturb this, such that a power law fit between these mass limits yields an effective slope of -2.357.  Figure A2. The effective IMF slopes (determined from the full stellar population between 0.5 and 300 M ⊙ ) for the variant models described in section 4. In each case, the nominal IMF has a slope of -2.350, and secondary stars act to reduce the IMF slope. The slope is shown against the binary fraction at 1 M ⊙ for set 1 and the binary fraction at 30 M ⊙ for set 2 and BPASS v2.2.
in a truncated or sparsely sampled upper mass limit to the initial mass function (although this only applies in stellar populations with total masses 10 3 M ⊙ , see Eldridge 2012, for discussion). Fig. A1 illustrates this process for the canonical "imf135 300" version of BPASS v2.2. The number (i.e. final model weighting) of stars in each mass bin is presented separately for single stars, primaries and secondaries, as well as for the total population. The resultant lines are not entirely smooth -the result of the limited grid of mass ratios  Figure A3. The effective IMF slopes (determined from the full stellar population between 0.5 and 300 M ⊙ ) for the variant (v) models described in section 2. In each case, the nominal IMF has a slope of -2.350, and secondary stars act to reduce the IMF. The slope is shown against the binary fraction at 1 M ⊙ and at 30 M ⊙ , with the pairs of values at these two masses linked by a line. The BPASS v2.2 canonical distribution is shown by the thick line.
permitted for each primary star mass, which result in some secondaries being shifted to the nearest available model bin. A power law function of the form N(M) ∝ M −α is then fit to the total population between 0.5 and 300 M ⊙ to produce an effective power law slope. For the canonical BPASS v2.2. model the nominal slope of -2.350 is altered to an effective power law slope of -2.357: an insignificant change compared to the large slope differences between IMFs being modelled in that data release (-2.7, -2.35, -2.0).
The small change is largely the result of the dependence of binary fraction on primary mass. The high mass stars which have many companions are relatively few in number to begin with, while the abundant low mass stars have relatively few companions. Thus the perturbation to the IMF is small. The same cannot be said for all the variant prescriptions explored in this study. The extreme variants explored in section 4 include those where the binary fraction is near constant with mass, such that the fraction of low mass stars with companions approaches unity. Fig. A2 demonstrates the effect of these models on the effective initial mass function (again with a nominal slope of -2.35, and assuming a fixed mass ratio distribution). As the figure shows, the effect of either decreasing the high mass binary fraction, or increasing the low mass binary fraction is to steepen the effective IMF slope in the models, with the extreme cases shown in the figure reaching a slope of -2.40.
We have also calculated the effective IMF slope measured in the variant (v) models of section 2 and show these in Fig. A3. Note that in this case there is no clear trend in IMF slope with overall binary fraction for primary stars at either 1 M ⊙ or 30 M ⊙ , since the binary fraction itself is varying over a relatively narrow range (as permitted by empirical constraints) and perturbations in the mass ratio distribution dominate the effective IMF slope determination.