Nebular emission from young stellar populations including binary stars

We investigate the nebular emission produced by young stellar populations using the new GALSEVN model based on the combination of the SEVN population-synthesis code including binary-star processes and the GALAXEV code for the spectral evolution of stellar populations. Photoionization calculations performed with the CLOUDY code confirm that accounting for binary-star processes strongly influences the predicted emission-line properties of young galaxies. In particular, we find that our model naturally reproduces the strong HeII/Hb ratios commonly observed at high Hb equivalent widths in metal-poor, actively star-forming galaxies, which have proven challenging to reproduce using previous models. Including bursty star formation histories broadens the agreement with observations, while the most extreme HeII equivalent widths can be reproduced by models dominated by massive stars. GALSEVN also enables us to compute, for the first time in a way physically consistent with stellar emission, the emission from accretion discs of X-ray binaries (XRBs) and radiative shocks driven by stellar winds and supernova explosions. We find that these contributions are unlikely to prominently affect the predicted HeII/Hb ratio, and that previous claims of a significant contribution by XRBs to the luminosities of high-ionization lines are based on models predicting improbably high ratios of X-ray luminosity to star formation rate, inconsistent with the observed average luminosity function of XRBs in nearby galaxies. The results presented here provide a solid basis for a more comprehensive investigation of the physical properties of observed galaxies with GALSEVN using Bayesian inference.


INTRODUCTION
The James Webb Space Telescope (JWST) has opened a new window on the rest-frame ultraviolet and optical emission of young galaxies at the epoch of reionization.The emission-line signatures of such galaxies should contain valuable information on the sources that potentially contributed to reionize the Universe.Over the past several years, numerous studies have been conducted at various redshifts to ★ E-mail: lecroq@iap.frcharacterize and analyse the ultraviolet and optical spectra of metalpoor, actively star-forming galaxies approaching the properties of these reionization-era galaxies (e.g., Stark et al. 2014;Steidel et al. 2016;Amorín et al. 2017;Senchyna et al. 2017;Nakajima et al. 2018;Berg et al. 2022).These observations have put a heavy strain on models designed to interpret the nebular emission from young galaxies, as some emission lines requiring highly energetic photons, such as He ii-recombination lines, occasionally exhibit intensities so high that they cannot be reproduced by standard models fitting lowerionization lines (e.g., Shirazi & Brinchmann 2012;Nanayakkara et al. 2019;Stanway & Eldridge 2019;Plat et al. 2019;Schaerer et al. 2019;Olivier et al. 2022).
Several production mechanisms have been suggested for these missing energetic photons.A most natural one is the hard radiation from hot, nearly pure-He stars produced through processes induced by binary interactions (such as envelope stripping, quasihomogeneous evolution, and common-envelope ejection), often neglected in stellar population synthesis models (e.g., Eldridge & Stanway 2012;Eldridge et al. 2017;Götberg et al. 2020).This possibility is all the more compelling in that ∼ 70 per cent of massive stars in nearby stellar populations are thought to be part of binary systems (e.g.Sana et al. 2012;Moe & Di Stefano 2017).Yet, the BPASS models of Stanway & Eldridge (2018), which incorporate the above processes, do not appear to produce enough energetic radiation to account for the observations (Stanway & Eldridge 2019).Other proposed origins of the surprisingly strong high-ionization lines in some metal-poor star-forming galaxies include fast radiative shocks, accretion on to compact objects and active galactic nuclei (AGNs; e.g., Izotov et al. 2012;Nakajima et al. 2018;Plat et al. 2019;Schaerer et al. 2019;Umeda et al. 2022;Katz et al. 2023).Such components are typically modelled independently of the stellar population and then scaled to reproduce the data under consideration.While insightful, this approach suffers from the weakness of not guaranteeing physical consistency between the different spectral components involved.
In this paper, we explore the emission-line properties of interstellar gas photoionized by young binary-star populations computed using a new spectral-synthesis model, built by combining the SEVN (Iorio et al. 2023) and GALAXEV (Bruzual & Charlot 2003) populationsynthesis codes, coupled with the CLOUDY photoionization code (Ferland et al. 2017).This spectral-synthesis model, called GAL-SEVN, enables us to compute, for the first time in a physically consistent way, the emission from stars, accretion discs of X-ray binaries (XRBs) and radiative shocks driven by stellar winds and supernova (SN) explosions.We do not consider here the emission from AGNs, which cannot be tightly related to the star formation properties of a galaxy, and whose effects have been studies in detail by Plat et al. (2019).We show how the GALSEVN model can better account for observations of ultraviolet and optical emission lines in metal-poor star-forming galaxies -in particular strong He ii emission -than previous models of both single and binary stars.When combined with a Bayesian-inference tool such as BEAGLE (Chevallard & Charlot 2016), the model presented here should provide valuable constraints on physical parameters of distant galaxies from observations of nebular emission.
The paper layout is as follows: in Section 2, we present the new GALSEVN model of spectral evolution of binary-star populations and our procedure to compute the associated nebular emission with CLOUDY.In Section 3, we describe a reference observational sample and examine the ability of models with different stellar and nebular parameters to reproduce these observations.Then, in Section 4, we consider additional effects, such as emission from XRB accretion discs and radiative shocks, potentially helpful to reproduce the most extreme emission-line properties of galaxies in the sample.We place our findings in the context of other recent studies in Section 5. Section 6 summarizes our results.

THE MODELS
In this section, we describe the models used to compute the nebular emission from young galaxies.We follow the approach of Charlot & Longhetti (2001, see also Gutkin et al. 2016) and express the luminosity per unit frequency  emitted at time  by a star-forming galaxy as where ( − ′ ) is the star formation rate at time  − ′ ,   [ ′ ,  ( − ′ )] the luminosity produced per unit frequency per unit mass by a single generation of stars of age  ′ and metallicity  ( − ′ ), and   (,  ′ ) the transmission function of the interstellar medium (ISM).We compute   for populations of single and binary stars using a combination of the SEVN (Iorio et al. 2023) and GALAXEV (Bruzual & Charlot 2003) population-synthesis codes (Section 2.1) and   using the CLOUDY photoionization code (Ferland et al. 2017, Section 2.2).

Stellar population modelling
To describe the emission from populations of single and binary stars, we appeal to a combination of the SEVN population-synthesis code (Spera et al. 2015;Spera & Mapelli 2017;Spera et al. 2019;Mapelli et al. 2020;Iorio et al. 2023) with the GALAXEV code for the spectral evolution of stellar populations (Bruzual & Charlot 2003).We briefly recall the main features of these codes in the following paragraph, referring the reader to the original studies for more details.SEVN (Stellar EVolution for -body) is an open-source binary population-synthesis code,1 based on the interpolation of stellar properties (e.g., mass, radius, luminosity, core mass, core radius) from interchangeable libraries of pre-computed evolutionary tracks to describe the temporal evolution of stellar populations including binary-evolution processes.In this work, we adopt the PARSEC (PAdova and TRieste Stellar Evolution Code) library of evolutionary tracks (including pre-main sequence evolution) for stars with initial masses between 2 and 600 M ⊙ and metallicities in the range 10 −11 ⩽  ⩽ 0.04 (Bressan et al. 2012;Chen et al. 2015;Costa et al. 2019Costa et al. , 2023;;Nguyen et al. 2022;Santoliquido et al. 2023).
For stars ending their lives as supernovae (SNe), SEVN includes several prescriptions to compute the compact-remnant (neutron-star or black-hole) mass and natal kick depending on SN type, i.e., electron-capture (Giacobbo & Mapelli 2019), core-collapse (Fryer et al. 2012) or pair-instability (Mapelli et al. 2020) SNe.We adopt here the delayed SN model, which predicts a smooth transition between maximum neutron-star mass and minimum black-hole mass.We generate natal kicks as described by Giacobbo & Mapelli (2020), in agreement with the proper-motion distribution of young Galactic pulsars (Hobbs et al. 2005), and with reduced kick magnitude for stripped and ultra-stripped SNe (Tauris et al. 2017).
The binary-evolution processes incorporated in SEVN include, most notably, mass transfer driven by winds and Roche-lobe overflow, common envelope evolution, removal of stellar angular-momentum by magnetic braking, the effect of stellar tides on orbital motions, orbital decay through gravitational-wave emission and stellar mergers.These are described in depth in section 2.3 of Iorio et al. (2023).For convenience, we provide in Appendix A a summary describing how SEVN handles binary-evolution products using the PARSEC tracks adopted in this study.In the simulations presented in this paper, mass transfer is assumed to be always stable for donor stars on the main sequence and in the Hertzsprung gap.For other stars, mass transfer stability depends on the binary mass ratio and the evolutionary stage of the donor star (following option 'QCRS' in table 3 of Iorio et al. 2023).Also, we set here the adjustable fraction of orbital energy converted into kinetic energy during common-envelope evolution to  CE = 3 (section 2.3.3 of Iorio et al. 2023, see also Hurley et al. 2002).
Finally, in this work, we adopt the formalism of quasihomogeneous evolution introduced by Eldridge et al. (2011, see section 2.3.2 of Iorio et al. 2023).According to this formalism, a low-metallicity, main-sequence star spun up by accretion of substantial material via stable Roche-lobe overflow mass transfer sees its core replenished with fresh hydrogen through rotational mixing and remains fully mixed until it burns all its hydrogen into helium (ending as a pure-He star), at nearly constant radius.Following Eldridge et al. (2011, see also Eldridge & Stanway 2012), we include this evolution for stars with metallicity  ⩽ 0.004 (which have weak-enough winds to prevent strong loss of angular momentum) accreting at least 5 per cent of their mass and with post-accretion masses greater than 10 M ⊙ .Stars meeting these criteria brighten and see their temperature increase as they burn hydrogen at nearly constant radius during quasi-homogeneous evolution (e.g., Eldridge & Stanway 2012).Their inclusion augments the population of compact, hot luminous stars, which has important implications for the results presented in this paper.
We set the initial conditions as follows.For a given metallicity , we produce with SEVN a stochastic population of 10 6 evolving binary pairs,2 with initial primary-star masses in the range 2 ⩽  ⩽ 300 M ⊙ drawn from a power-law IMF () ∝  −1.3 (where () is the number of stars created with masses between  and  + ).For each pair, we draw the ratio of initial secondary-star mass to initial primary-star mass (), the orbital period () and the eccentricity () from the corresponding probability density functions (PDFs) adopted in Sana et al. (2012): PDF() ∝  −0.1 with  ∈ [0.1, 1.0]; PDF(P) ∝ P −0.55 with P = log(/day) ∈ [0.15, 5.5]; and PDF() ∝  −0.42 with  ∈ [0, 0.9]. 3 The adopted primary-star IMF produces a larger proportion of massive stars than the standard Galactic-disc IMF, which has () ∝  −2.3 for  ⩾ 1 M ⊙ (e.g., Kroupa 2001;Chabrier 2003).We make this choice to appropriately sample the high-mass end of the IMF among this base collection of one million pairs, from which we can then extract populations with different IMFs (including Chabrier 2003).
To compute the spectral evolution of stellar populations obtained this way, we adopt an approach similar to that in the GALAXEV code (Bruzual & Charlot 2003).Specifically, we assign a spectrum to each SEVN star by selecting among the different stellar libraries listed in appendix A of Sánchez et al. (2022), supplemented by WM-BASIC model atmospheres computed by Chen et al. (2015) for stars hotter than 50,000 K.For a given metallicity, these models are interpolated at the mass-loss rate, effective temperature (log  eff ) and gravity (log ) of the considered star.The spectra of stars hotter than the hottest Chen et al. (2015) model (i.e., with  eff > 10 5 K) are approximated by black bodies.We identify Wolf-Rayet (WR) stars of different types (WNL, WNE, WC and WO) in the populations generated with the SEVN code following the procedure outlined in section 3.2 of Chen et al. (2015), based on the H, C, N and O surface abundances of stars hotter than 25,000 K.We appeal to the highresolution version of the Potsdam-WR (PoWR) models to describe the spectra of these stars (Hamann & Gräfener 2004;Sander et al. 2012;Todt et al. 2015;Hainich et al. 2019).For each WR star, we select the PoWR model with closest parameters in a way similar to that outlined in appendix A of Plat et al. (2019), with the refinement that we also match the surface hydrogen fraction of WNL-type stars.For core-He burning stars less massive than 8 M ⊙ and hotter than 20,000 K, we adopt CMFGEN spectral models of stripped-envelope stars by Götberg et al. (2018).All stellar spectra are re-sampled on a common wavelength scale ranging from 5.6 Å to 360 m in 16,902 steps using the SpectRes tool (Carnall 2017).
In the course of this procedure, we record cases of neutron stars and black holes accreting mass from a companion, where the accretion disc can produce X-ray emission (Section 4.3).Also, we record rates of pair-instability supernovae (PISNe), exploding (i.e., non-failed; Spera et al. 2015) core-collapse SNe and type-Ia SNe, which can lead to shock-driven nebular emission (Section 4.4).
In the remainder of this paper, we refer to the combined SEVN and GALAXEV spectral evolution model as simply 'GALSEVN'.More details on the coupling between SEVN and GALAXEV will be provided elsewhere (Bruzual et al., in preparation), as well as predictions for the spectral evolution   (Charlot et al., in preparation).

Photoionization modelling
The transmission function   (,  ′ ) of the ISM in equation ( 1) incorporates absorption of stellar radiation by gas and dust in different phases of the ISM (ionized interiors and neutral envelopes of stellar-birth clouds, diffuse inter-cloud medium), as well as the production of nebular line-plus-continuum emission and dust emission (e.g., Charlot & Longhetti 2001;da Cunha et al. 2008;Vidal-García et al. 2017).Here, we are primarily interested in young star-forming galaxies with age close to the typical timescale for dissipation of giant molecular clouds (∼ 10 Myr, ; e.g., Murray et al. 2010;Murray 2011).In this case, we can write   (,  ′ ) ≡  +  ( ′ ), where  +  ( ′ ) is the transmission function of the ionized gas (Gutkin et al. 2016).
We compute  +  ( ′ ) using version c17.00 of the CLOUDY photoionization code (Ferland et al. 2017). Following Charlot & Longhetti (2001), we assume for simplicity that galaxies are ionization-bounded and that the transfer of radiation through ionized gas can be described galaxy-wide by means of 'effective' parameters, the main ones being (Gutkin et al. 2016, see also section 2.1 of Plat et al. 2019): • The hydrogen gas density,  H .
• The gas-phase metallicity, which unless otherwise specified is taken to be equal to that of the ionizing stars,  ISM = .
• The carbon-to-oxygen abundance ratio, C/O.
• The dust-to-metal mass ratio reflecting the depletion of heavy elements on to dust grains,  d .
• The zero-age, volume-averaged ionization parameter, ⟨⟩ ≡ ⟨⟩( ′ = 0).In spherical geometry, ⟨⟩( ′ ) can be expressed in terms of the rate of ionizing photons produced by the evolving stellar population,  H ( ′ ), the volume-filling factor of the gas (i.e., the ratio of the volume-averaged hydrogen density to  H ), , and the case-B hydrogen recombination coefficient,  B , as (e.g., Panuzzo et al. 2003): We stop the photoionization calculations at the outer edge of the H ii region, when the electron density falls below 1 per cent of  H (we adopt a default inner radius of 0.1 pc).We adopt the prescription of Gutkin et al. (2016) for the abundances and depletions of interstellar elements, as well as for scaling the total (primary+secondary) nitrogen abundance with that of oxygen.In this framework, the present-day solar (photospheric) metallicity is  ⊙ = 0.01524, corresponding to a proto-solar metallicity  0 ⊙ = 0.01774 (Bressan et al. 2012).It is of interest to examine how  H ( ′ ) depends on the fraction of binary stars in a stellar population according to the models presented in Section 2.1 (see also section 2.4.1 of Eldridge et al. 2017).This is shown in the top panel of Fig. 1 for an instantaneous-burst, 'simple stellar population' [SSP; corresponding to ( ′ ) = ( ′ ) in equation 1], for the typical metallicity  = 0.001 (i.e., about 6 per cent of solar) of galaxies in the sample of Section 3.1.The middle and lower panels show the equivalent plots for He i-and He ii-ionizing photons.In each panel, the different curves correspond to different zero-age binary fractions,  bin ≡  bin ( ′ = 0), defined as  bin = B/(S + B), where B is the number of binary pairs and S the number of single stars at  ′ = 0. To isolate the effects of binary interactions, the IMF is taken to be exactly the same for binary and single stars.This is achieved by drawing single stars as binary systems whose member stars never interact, i.e., assuming a Chabrier (2003) IMF for primary stars and a Sana et al. (2012) distribution of secondary-to-primary mass ratios (see Section 2.1), with orbital parameters (large semi-major axes and circular orbits) such that primary and secondary stars never interact.For reference, the total IMF (including primary and secondary stars) drawn in this way has an effective near-Salpeter (1955) slope, i.e., roughly () ∝  −2.35 , at  > 10 M ⊙ .All models are normalized to a total initial mass in stars of 1 M ⊙ between 0.1 and 300 M ⊙ . 4 Fig. 1 shows that, at ages younger than 1 Myr, models with  bin > 0 produce slightly larger  H ,  He I and  He II than single-star models, due to the presence of merged stars with masses up to ∼ 600 M ⊙ on the upper main sequence.Then, when the most massive stars leave the main sequence, all three quantities start to drop.The rate of He ii-ionizing photons rises sharply again after 2 Myr in models with  bin > 0. This is caused by the appearance of the hot (> 10 5 K), pure-He (WNE-like) products of the first massive stars to lose their H-rich envelope, either by dumping mass on their companion through Roche-lobe overflow, or by being driven to quasi-homogeneous evolution and burning all of their H to He after receiving mass from their companion, or through common-envelope ejection (see, e.g., Spera et al. 2019;Iorio et al. 2023).These stars have a comparatively weak influence on the rates of He i-and H-ionizing photons, which remain dominated by the much more numerous H-burning stars on the upper main sequence, at least until ages of several million years.Then, as the He-burning main sequence continues to develop, stellar mergers conspire to maintain the tip of the H-burning main sequence at a brighter and hotter point in binary-star models than in single-star models, leading to significantly larger  H ,  He I and  He II .After a few tens of million years, the difference between models with  bin > 0 and  bin = 0 declines, especially for  He II after the disappearance of the brightest and hottest WNE-like stars (around 40 Myr).The faintest of these stars, together with stellar mergers and stars undergoing quasi-homogeneous evolution, contribute to maintain  H and  He I relatively strong in binary-star models relative to single-star models.The rise of  He I and  He II at ages around 30 Myr in singlestar models is due to the development of the white-dwarf cooling sequence.From about 50 to 90 Myr, the hot ionizing radiation from the most massive of these stars leads to slightly higher  He II than in binary-star models, where their population is reduced through binary-star interactions (their progenitors avoiding white-dwarf evolution through mass gain, or being driven to fainter evolution through mass loss).
A notable feature of Fig. 1 is that even small fractions of binary stars can have a major impact on the predicted rates of highly energetic photons.At an age of 3 Myr, for example, the ratio of He ii-to H-ionizing photons,  He II /  H , which is 2700 times larger in models with  bin = 1 than in those with  bin = 0, is already 54, 480 and 1300 times larger for  bin = 0.01, 0.1 and 0.3, respectively.It is 2200 times larger for  bin = 0.7, often adopted as typical of massive-star populations (e.g.Sana et al. 2012).
In Fig. 2, we compare in more detail the ionizing spectra   ( ′ ) of SSPs with pure binary stars (  bin = 1, solid lines) and pure single stars (  bin = 0, dotted lines) at ages between 1 and 10 Myr, for 4 This normalization is performed by drawing primary-star masses from a Chabrier (2003) IMF (and the associated secondary stars from a Sana et al. 2012 distribution) over the full range 0.1-300 M ⊙ , even though only pairs with both primary-and secondary-star masses greater than 2 M ⊙ are retained (the pairs drawn with both primary-and secondary-star masses greater than 2 M ⊙ are matched to those in the base collection of Section 2.1).This has no influence on the predicted  H ,  He I and  He II , since main-sequence stars less massive than 2 M ⊙ are cooler than ∼ 14, 000 K at any metallicity and do not produce significant ionizing radiation.Also, rejecting binaries with primary mass above 2 M ⊙ and secondary mass below this limit has a negligible impact on ionizing fluxes: artificially including such draws by assigning them the light of pairs with same primary mass and a secondary mass of 2 M ⊙ would increase  H ,  He I and  He II by less than one per cent at ages below 40 Myr in Fig. 1.  1) (in units of the luminosity per unit frequency at the Lyman limit) for the SSP models of Fig. 1 with pure binary stars (  bin = 1, solid lines) and pure single stars (  bin = 0, dotted lines) at the ages  ′ = 1, 3, 5, 8 and 10 Myr (colour-coded as indicated on the right-hand side).Also shown in grey is the area sampled by ionizing spectra of AGNs with power-law indices between  AGN = −2.0(bottom edge) and −1.2 (top edge; Feltre et al. 2016).The ionization energies of some common species are indicated for reference.a metallicity  = 0.001.As expected from Fig. 1, at an age of 1 Myr, the spectra of both models nearly overlap, while at older ages, the model with binary stars exhibits an excess of energetic photons relative to that with single stars.Populations of binary stars will therefore boost the production of highly-ionized species in the gas they irradiate (we indicate for reference the ionization energies of some common species in Fig. 2).In fact, Fig. 2 shows that at ages  ′ ⩾ 3 Myr in the binary-star model, the shape of   ( ′ ) at ionizing energies up to ≳ 50 eV is similar to that of hard ionizing spectra of AGNs,   ∝   AGN with −2.0 ⩽  AGN ⩽ −1.2 (e.g.Feltre et al. 2016, grey shaded area).

COMPARISON WITH OBSERVATIONS
In this section, we compare the predictions of the models presented in Section 2 with observations of ultraviolet and optical emission lines in metal-poor star-forming galaxies.We start by describing the observational sample (Section 3.1), after which we examine the predictions of models with standard star and gas parameters (Section 3.2) and explore the influence of some key adjustable parameters on the results (Section 3.3).
In this paper, we compare the predictions of our model with the observed intensity ratios and equivalent widths (EWs) of selected ultraviolet (C iv 1549, He ii 1640, O iii] 1664 and C iii] 1908) and optical ([O iii] 5007, [O ii]3727, [O i] 6300, He ii 4686, H, H and [N ii] 6584) emission lines commonly observed in metal-poor star-forming galaxies.We will pay particular attention to the He ii lines, which require ionization energies above 54.4 eV and represent a challenge for current models.

Standard models
We start by exploring how models with 'standard' parameters compare with the observations described above.As in Section 2.2, we consider stellar populations with metallicity  = 0.001, typical of the observational sample.For the photoionized gas, we adopt the same metallicity ( ISM = 0.001) and a density ( H = 100 cm −3 ), ionization parameter (log ⟨⟩ = −2.0)and dust-to-metal mass ratio ( d = 0.3) typical of metal-poor star-forming galaxies (e.g., Berg et al. 2016;Gutkin et al. 2016;Mingozzi et al. 2022) 2016) corresponds to a typical carbonto-oxygen abundance ratio C/O ≈ 0.17, i.e., about 40 per cent of the solar ratio (C/O) ⊙ = 0.44 (Gutkin et al. 2016).This is close to the typical value C/O ≈ 0.18 reported by Izotov et al. (2023) in lowredshift LyC-leaking galaxies with 12 + log(O/H) ≲ 8.1.Yet, given that galaxies in the observational sample of Section 3.1 exhibit C/O values up to 120 per cent of solar, we expect some of these galaxies to show stronger carbon lines than predicted by our standard models.We report the above standard nebular parameters in Table 1.We perform the photoionization calculations at stellar-population ages up to 10 Myr, i.e., roughly the timescale of dissipation of giant molecular clouds in star-forming galaxies (Murray 2011;Ma et al. 2015).
In Fig. 3, we show the intensity ratios and equivalent widths of emission lines predicted by the SSP models of Fig. 1 with pure binary stars (  bin = 1, labelled 'GALSEVN') and pure single stars (  bin = 0, 'GALSEVN single') at ages up to 10 Myr (colour-coded as indicated on the right).It is instructive to start by examining the behaviour of these models in the He ii 4686/H-versus-EW(H) diagram (Fig. 3b), where observations are traditionally challenging to reproduce with photoionization models powered by young stellar populations (e.g., Plat et al. 2019;Schaerer et al. 2019).All models start at high EW(H) on the right-hand side of the diagram.As massive stars leave the main sequence and the first red super-giant stars appear in the single-star GALSEVN model,  H and  He II decrease while the continuum at H increases, causing both EW(H) and He ii 4686/H to drop without ever reaching the area of the diagram occupied by the observations.
Remarkably, in the binary-star GALSEVN model, the rise in  He II triggered by the appearance of hot WNE-like stars with blue spectra around 3 Myr (Fig. 1) makes both He ii 4686/H and EW(H) rise in such a way that the models spend all the time thereafter up to 10 Myr in the same location as the observations.It is worth recalling that, as shown by Fig. 1, similar results would be obtained with models containing only a minor fraction of binary stars (  bin ≳ 0.3).
Also shown for comparison in Fig. 3(b) are photoionization calculations performed using the BPASS v2.2.1 models of Stanway & Eldridge (2018) for single and binary stars, as well as the C&B single-star model used by Plat et al. (2019, based on an updated version of the Bruzual & Charlot 2003 GALAXEV code), for a metal-licity  = 0.001 and the same nebular parameters as in Table 1.The CLOUDY calculations for these models were conducted similarly to those for the GALSEVN models as part of the present study.We note the striking difference between the C&B and GALSEVN single-star models at ages around 2 Myr, when the brief (∼ 0.2 Myr) WR phase that makes He ii 4686/H rise up to the observed range in the C&B model is absent from the GALSEVN model.This can also be seen in the evolution of  H and  He II shown in Fig. 4.Both models being based on the same library of stellar spectra (described in Section 2.1), this difference originates from recent updates (in particular in the opacities) of the PARSEC-track library of Bressan et al. (2012) and Chen et al. (2015), leading to reduced mass loss, and hence fewer WR stars at low metallicity (Costa et al. 2021;Nguyen et al. 2022).
The brief WR phase at early ages is also present (although to a lesser extent) in the BPASS models for both single and binary stars, which are based on the Cambridge STARS code (Eggleton 1971;Eldridge et al. 2008, see also Fig. 4).After this phase and the disappearance of the most massive stars, EW(H) and He ii 4686/H both drop rapidly in the single-star model, as in the GALSEVN and C&B single-star models.Instead, in the BPASS binary-star model, He ii 4686/H is maintained at a fairly high level by the appearance of envelope-stripped and spun-up stars, still without reaching the region of the diagram probed by the observations.Given that in GALSEVN,  He II is dominated by WNE stars after a few Myr, the discrepancy between the GALSEVN and BPASS binary-star models in Fig. 3(b) can likely be attributed to a distinct population of these stars in BPASS, particularly as both models use similar spectra for them (from the PoWR library).
The behaviour of the models in the other diagrams of Fig. 3 follows from the same arguments.In the classical Baldwin et al. (1981, hereafter BTP) diagram to distinguish AGNs from star-forming galaxies (Fig. 3a), at early ages, all models lie at the high-[O iii] 5007/H, low-[N ii] 6584/H end characteristic of low metallicity and high ionization parameter (see, e.g., fig. 2 of Gutkin et al. 2016).The GALSEVN and BPASS binary-star models stay in this area at ages up to 10 Myr, while the [O iii] 5007/H ratio of single-star models rapidly drops after the disappearance of photons more energetic than ∼ 35.1 eV capable of ionizing O 2+ .The implied drop in [O iii] 5007 luminosity is also the reason for the difference between binary-and single-star models in the 3c).In Fig. 3(d), none of the models (plotted for clarity at slightly different abscissae even though all have 12 + log(O/H) = 7.53) reaches the highest observed EW(He ii 1640) ∼ 2-8 Å.In this diagram, only the GALSEVN binary-star models maintains EW(He ii 1640) > 0.5 Å at ages  ′ ≳ 3 Myr, while the BPASS binary-star model declines to EW(He ii 1640) ≈ 0.03 Å and all single-star models to EW(He ii 1640) ≈ 0.
The difference between the GALSEVN binary-star model and all other models is more striking in Fig. 3(e), where GALSEVN is the only model able to account for the properties of star-forming galaxies similar to those of AGNs at low (C iv 1549 +C iii] 1908)/He ii 1640 and high C iv 1549/C iii] 1908.Again, albeit without certainty, the difference between the GALSEVN and BPASS binary-star models likely arises from distinct populations of WNE stars, which dominate the production of He ii-ionizing photons.Similarly, in Fig. 3(f), this model at ages  ′ ≳ 2 Myr is closest to the observational sequence at He ii 1640/O iii] 1664 > 0.3, although it does not appear to reach sufficiently high EW(C iii] 1908) at a given He ii 1640/O iii] 1664.
On the whole, Fig. 3 shows the need to include binary stars and their inherent physical processes in modelling the emission-line proper-  2001) and Kauffmann et al. (2003), respectively, to separate AGN-dominated from star-forming galaxies, while in (e), the dashed line shows the equivalent criterion proposed by Nakajima et al. (2018).The curves interspersed with filled squares, triangles, pentagons, hexagons and crosses show the predictions of the GALSEVN binary-and single-star models, the BPASS binary-and single-star models and the C&B single-star model, respectively, colour-coded with age as indicated on the right (the BPASS models, available only at ages above 1 Myr, do not sample the youngest ages in the colour bars).In (d), the models are shown for clarity at slightly different abscissae, although all have 12 + log(O/H) = 7.53.In all models, the ionizing stellar population is an SSP with metallicity  = 0.001 and the same zero-age Chabrier (2003) IMF as in Fig. 1, while the parameters of the photoionized gas are those listed in Table 1.ties of metal-poor, star-forming galaxies.It also validates reasonable agreement between GALSEVN binary-star models and (non-extreme) observations of such galaxies in our sample, although discrepancies between some observed and modelled spectral features remain to be explored.

Influence of adjustable model parameters
So far, we have considered only models with metallicity  = 0.001 and nebular parameters fixed at the standard values of Table 1.Given the diversity of objects in the observational sample of Section 3.1, it is relevant to explore the effect of changing these parameters.Here, we focus in turn on the influence of metallicity, volume-averaged ionization parameter and carbon-to-oxygen abundance ratio.We do not explore changes in  H and  d , which have been shown by Plat et al. (2019) to have a minimal impact on the emission-line ratios presented in Fig. 3.We adopt the same Chabrier (2003) IMF as in the standard models.

Metallicity
Fig. 5 shows the evolution of GALSEVN binary-star models with metallicities  =  ISM = 0.0005, 0.001, 0.002 and 0.004 in the same diagrams as in Fig. 3, with all other parameters fixed at their values in Table 1.An increase in metallicity makes cooling through collisionally-excited metal transitions more efficient, which reduces the electronic temperature  e .This causes the luminosity ratios of metal lines to H and He lines to rise and the ratios of high-to low-ionization lines to drop (Figs 5a, c, e and f).The increase of [N ii] 6584/H with  in Fig. 5(a) is further enhanced by the inclusion of secondary nitrogen enrichment in the nebularemission models (Section 2.2).At  = 0.004, we note the emergence of classical WR stars at early ages following the increased mass loss of evolved massive stars with metal-rich envelopes, which causes a brief excursion of the model at high He ii 4686/H and He ii 1640/O iii] 1664 in Figs.5(b) and (f).
As Fig. 5(d) shows, while metallicities in the range 0.0005 ⩽  ⩽ 0.004 can fully describe the range of 12 + log(O/H) probed by the observational sample, metallicity has little influence on the maximum EW(He ii 1640) achievable with the GALSEVN binary-star model.

Ionization parameter and C/O abundance ratio
In Fig. 6, we compare the standard GALSEVN binary-star model of Fig. 3 with two GALSEVN binary-star models: one with ionization parameter log ⟨⟩ = −1 instead of −2; and one with C/O = (C/O) ⊙ = 0.44 instead of 0.17.All models have metallicity  = 0.001.
The effect of increasing log ⟨⟩ is to increase the probability of multiple ionizations, and hence to strengthen high-ionization lines (see equivalent widths of He ii 1640 and C iii] 1908 in Figs 6d and  f) and increase ratios of high-to-low ionization lines (e.g., C iv/C iii] and O iii]/[O ii]; Figs 6c and e).In Fig. 6(a), the drop in [N ii]/H as log ⟨⟩ increases results from the conversion of N + into N 2+ .In Fig. 6(b), the drop in EW(H) is caused by the increase in H-column density, and hence in the absorption of H-ionizing photons by dust at fixed  and  d , as the Strömgren radius increases (e.g., Plat et al. 2019).Overall, Fig. 6 shows that enhancing the ionization parameter does not significantly improve the agreement between GALSEVN binary-star models and the observations considered here.We note that, while ionization parameters as high log ⟨⟩ = −1 have been favoured for only a few objects in this observational sample (see table 1 of Plat et al. 2019), such parameters may be more common in young star-forming galaxies near the epoch of reionization (e.g., Cameron et al. 2023a).
Increasing the C/O abundance ratio at fixed metallicity means raising the abundance of carbon while reducing that of all other metals (e.g., Gutkin et al. 2016).This results in a strengthening of C lines and a slight weakening of O lines, with negligible effect on H and He lines, as Fig. 6 shows.Interestingly, the adoption of solar C/O in place of the standard value in Table 1 (favoured by the trend with 12 + log(O/H) in fig.6 of Berg et al. 2016) appears to provide better agreement with the data in Fig. 6(f).Such high C/O ratio at low metallicity is in agreement with the measured properties of some metal-poor galaxies in our sample, particularly those in the subsamples from Berg et al. (2016Berg et al. ( , 2019)), Amorín et al. (2017) and Stark et al. (2014).Recently, D'Eugenio et al. ( 2023) also reported the existence of a metal-poor star-forming galaxy with strong C iii] emission corresponding to slightly super-solar C/O ratio at redshift  ∼ 12.5.

O/Fe abundance ratio
Several recent studies have found evidence for super-solar ratios of  to iron-peak elements (e.g., O/Fe) in metal-poor star-forming galaxies at various redshifts, as expected for chemically young objects enriched by core-collapse SNe (e.g., Steidel et al. 2016;Strom et al. 2018;Topping et al. 2020;Cullen et al. 2021;Senchyna et al. 2022).Since Fe dominates the opacity of stellar ultraviolet spectra (which is relatively insensitive to O/H), while O, the most abundant heavy element in H ii regions, dominates gas-phase metallicity (interstellar iron being almost entirely depleted on to dust grains), line emission of gas with super-solar O/Fe photoionized by metal-poor stars may be explored using models with fixed O/Fe, but adopting a lower total metallicity for the stars than for the gas (e.g., Steidel et al. 2016).
In Fig. 7, we compare the standard GALSEVN binary-star model of  At fixed  ISM = 0.001, lowering the stellar metallicity from  = 0.001 to 0.0001 implies bluer stellar-population spectra and harder ionizing radiation, since stars evolve at higher effective temperature (e.g., Bressan et al. 2012;Plat et al. 2019).The H equivalent width increases significantly in Fig. 7(b) (due to both higher  H and fainter  4861-Å continuum), as does the He ii 4686/H ratio at ages until ∼ 3 Myr, when the ionizing radiation becomes dominated by pure-He, WNE-type stars at both metallicities.A similar effect is seen at early ages in the (C iv 1549 +C iii] 1908)/He ii 1640 (Fig. 7e) and He ii 1640/O iii] 1664 ratios (Fig. 7f).
The model with  =  ISM = 0.0001 in Fig. 7 produces a luminosity ratio and equivalent widths of H and He ii 1640 comparable to those in the model with the same  = 0.0001 but higher  ISM = 0.001 (Figs 7b and d).This is because these emission lines are primarily sensitive to the shape of the stellar-population spectrum.Instead, as expected from Fig. 5, the ratios of metal-to-H and He lines decrease significantly when  ISM is reduced from 0.0001 to 0.001 (Figs 7a, e and f).
On the whole, Fig. 7 shows that adopting super-solar O/Fe significantly affects the predicted emission-line properties of metal-poor star-forming galaxies, although we note that the impact on the equivalent width of He ii 1640 is relatively modest in Fig. 7(d).

INCLUSION OF ADDITIONAL EFFECTS
In the previous section, we have explored the effects of stellar and several ISM parameters (ionization parameter, C/O and O/Fe abundance ratios) on the predicted emission-line properties of binary-star SSPs computed with the GALSEVN code.We showed that, while these models produce significantly more energetic radiation than previous single-and binary-star models, they do not appear to reach the largest observed equivalent widths of lines with high ionization energies, such as He ii 1640 and C iii] 1908.In the following paragraphs, we consider additional effects associated with the evolution of stellar populations, which could potentially influence these observables.Unless otherwise specified, in all models in this section, we consider binary-star populations with  = 0.001 and the standard nebular parameters listed in Table 1.

Multiple bursts of star formation
So far, we have considered the emission-line properties of SSPs only, while even the most metal-poor nearby star-forming galaxies show evidence of several episodes of star formation (e.g., Izotov & Thuan 2004;Tolstoy et al. 2009;Senchyna et al. 2019).In Fig. 8, we show the effect of adding, at the age of 3 Myr, a second burst of star formation of the same amplitude as the first one to the evolution of an SSP.The sudden input of massive main-sequence stars boosts the rate of ionizing photons at that age, leading to a temporary rise in not only EW(H), but also C iv 1549/C iii] 1908 and EW(C iii] 1908) (see the difference between the orange and purple curves in Figs 8b, e and  f).As a result, the model probes a different region of the observational space than sampled by a pure SSP, the effect being less noticeable in the other diagrams of Fig. 8.
We can explore more realistic histories of star formation for metalpoor galaxies by appealing to the SPHINX suite of cosmological (radiation-hydrodynamical) simulations of young galaxies in the epoch of reionization (Rosdahl et al. 2022).The green curve in Fig. 8 shows the emission-line properties obtained for a typical galaxy experiencing several bursts of different amplitudes at ages up to 10 Myr in these simulations.The effects of three major bursts following the onset of star formation are clearly distinguishable through the associated boosts in EW(H), C iv 1549/C iii] 1908 and EW(C iii] 1908).This enables the model to probe new regions of the observational space not sampled by the SSP and two-burst models.
To illustrate the observational space sampled by models spanning a wider range of parameters, we show in Fig. 9 the emission-line properties of 800 models computed using the star formation histories or 100 SPHINX galaxies, combined alternatively with four metallicities ( =  ISM = 0.0005, 0.001, 0.002 and 0.004) and two C/O abundance ratios (C/O =0.17 and 0.44), at ages between 2 and 10 Myr.The aim here is not to fully nor uniformly explore the space of adjustable parameters, but rather to provide global insight into the influence of star formation history, metallicity and C/O ratio on the observables of interest to us.The two (yellow) peaks of the model density distributions appearing in almost all diagrams of Fig. 9 correspond to the two values of C/O sampled (see Fig. 6).Fig. 9 shows that the restricted set of models considered here samples quite widely the observational space.In particular, it is worth emphasizing that the area of the He ii 4686/H-versus-EW(H) diagram with the highest density of models coincides with that where the bulk of observations lie, which was out of reach of previous models (Fig. 3).Yet, some areas remained out of reach of these models, most notably those populated by galaxies with both largest He ii 4686/H and EW(H) (the most extreme of which being classified as LyC leakers; Fig. 9b), He ii 1640 equivalent widths in excess of ∼ 1.5 Å (Fig. 9d), largest C iv 1549/C iii] 1908 (Fig. 9e) and both largest EW(C iii] 1908) and He ii 1640/O iii] 1664 (Fig. 9f).

Top-heavy initial mass function
Since the ionizing radiation of a young galaxy is controlled by massive stars, the upper part of the IMF could have a significant impact on the predicted emission-line properties of the models.So far, we have adopted in all models the standard Chabrier (2003) IMF truncated at 0.1 and 300 M ⊙ , which can lead in binary-star populations to the presence of merged stars with masses up to ∼ 600 M ⊙ (Section 2.2).To explore the influence of a top-heavy IMF, we adopt the generalised Rosin-Rammler distribution function (e.g., Chabrier 2003;Wise et al. 2012;Goswami et al. 2022),  3), for  char = 50 M ⊙ and 100 M ⊙ , compared to the standard Chabrier (2003) IMF.In the lower panel, we show the evolution of the absolute AB magnitude at 1500 Å,  UV , of GAL-SEVN binary-star models computed with these three IMFs, assuming a total initial mass in stars of 10 6 M ⊙ , metallicity  = 0.001 and the standard nebular parameters of Table 1.Since most of the mass in the models with top-heavy IMFs is concentrated in massive stars with low mass-to-light ratios, at fixed total initial stellar mass, these models are about 2 mag brighter in  UV than that with a standard Chabrier (2003) IMF at early ages.The brightening in all models at ages before 2 Myr arises in part from stellar evolution on the main sequence, and in part from stellar mergers.Then, once the bulk of stars in the top-heavy IMF models leave the main sequence and complete their evolution, the population fades rapidly.Instead, the ultraviolet luminosity decreases more gradually in the model with standard    2003) IMF (in black).All IMFs are normalized to a total mass in stars of 1 M ⊙ .Bottom: absolute ultraviolet AB magnitude at 1500 Å,  UV , plotted against age for three GALSEVN binary-star models computed with the three IMFs shown in the top panel, assuming a total initial mass in stars of 10 6 M ⊙ , metallicity  = 0.001 and the standard nebular parameters of Table 1.For reference, the dotted horizontal line marks the magnitude  UV = −14 mag typical of nearby, extremely metal-poor star-forming galaxies (reached by the top-heavy-IMF models at the ages marked by the cyan stars).Chabrier (2003) IMF, as stars of lower and lower mass progressively leave the main sequence.
In Fig. 11, we present the emission-line properties of the three models shown in the lower panel of Fig. 10.To avoid cluttering the diagrams, we do not show the properties of models fainter than  UV = −14 mag (limit marked by a cyan star), corresponding roughly to the typical luminosity of nearby, extremely metal-poor star-forming galaxies (see, e.g., table 1 of Senchyna et al. 2019).At early ages, when massive stars on the upper main sequence dominate the emission, the three models overlap in all diagrams.Once the bulk of stars in the top-heavy IMF models leave the main sequence to become red super-giant stars, at ages around 2 Myr, the implied strong drop in ionizing-photon rate and rise in 4861-Å continuum emission make the H equivalent width decline more abruptly than in the standard model (Fig. 11b).Meanwhile, after reaching a minimum, the He ii 4686/H ratio (Fig. 11b) and He ii 1640 equivalent width (Fig. 11d) promptly rise again due to the emergence of hot, pure-He products of massive stars (Section 2.2).In Fig. 11(f), the C iii] 1908 equivalent width shows similar behaviour to EW(He ii 1640) relative to the standard model.EW(H) takes longer than EW(He ii 1640) and EW(C iii] 1908) to increase again (Fig. 11b), as red super-giant stars have stronger continuum emission at 4861 Å than at ultraviolet wavelengths.
Remarkably, Fig. 11 shows that, in diagrams involving ultraviolet lines, the two GALSEVN binary-star SSP models with top-heavy IMFs, fixed metallicity  = 0.001 and C/O = 0.17 investigated here reach regions of the observational space not sampled by the more complete grid of models with standard Chabrier (2003) IMF considered in Fig. 9 (up to He ii 1640 equivalent withs of ∼ 10 Å; Fig. 11b).Although these SSP models spend less than a few Myr with such extreme emission-line properties, they can potentially be combined with more complex star formation histories.In fact, they could represent a transient phase in the evolution of a galaxy, when a cluster of massive stars happens to form and dominates the light for a while.

Emission from accretion discs of X-Ray binaries
Several recent studies have suggested that the unusually strong He ii 1640 emission observed in the spectra of many metal-poor star-forming galaxies may originate from accretion discs of X-ray binaries, where a compact object (neutron star or stellar-mass black hole) accretes material from a companion (e.g., Schaerer et al. 2019;Umeda et al. 2022;Katz et al. 2023).Emission from the hot accretion discs of these systems peaks at X-ray energies (e.g., Mitsuda et al. 1984;Fragos et al. 2013;Mirocha 2014;Senchyna et al. 2019).In detail, the accretion rate depends on the evolutionary stage and chemical evolution of the donor (companion) star, the mass ratio between the donor and the compact object and their orbital separation.
In this context, the GALSEVN model, which follows the evolution of the physical properties and the emission from stars in binary systems, provides an ideal framework to compute in a self-consistent way the contribution by accretion discs of XRBs to the integrated light of stellar populations.In fact, in Appendix B, we present a simple model to compute the spectral energy distribution of XRB accretion discs in GALSEVN binary-star populations, based on a library of multicolour-disc spectra including potential Compton upward scattering of soft disc photons by energetic coronal electrons (Mirocha 2014;Steiner et al. 2009).This model, which also accounts for Xray variability (e.g., Tanaka & Shibazaki 1996;Chen et al. 1997), provides good agreement with the observed X-ray luminosities of various samples of nearby, metal-poor star-forming galaxies (Douna et al. 2015;Brorby et al. 2016;Lehmer et al. 2019), as well as with the average X-ray luminosity function of XRBs in the five most metal-poor star-forming galaxies observed by Lehmer et al. (2019).
As shown in Appendix B, such agreement with observations can be obtained with models in which most of the mass is accreted during episodes of X-ray outbursts, at a much higher rate than during the quiescent phase, the actual strength and duration of outbursts being degenerate.Fig. 12 shows the X-ray luminosity function of XRBs computed assuming that 99 per cent of the mass is accreted by the compact object during outbursts at a rate 100 greater than in the quiescent phase, for a GALSEVN binary-star population corresponding to 100 Myr of constant star formation at the rate  = 1M ⊙ yr −1 .The model has a binary-star fraction  bin = 0.7 (typical of massive-star populations; e.g., Sana et al. 2012) and the metallicity  = 0.008, close to the typical metallicity of the five metal-poor galaxies for which Lehmer et al. (2019) derived the average luminosity function shown by the black data points with error bars, which is well reproduced by the model.In Appendix B, we further show that metallicity has only a minor influence on the predicted X-ray luminosity function (Fig. B2).
In Fig. 13, we show the contribution by XRB accretion discs to the total H and He ii 4686 luminosities predicted by such a model for a GALSEVN binary-star population of standard metallicity  = 0.001.The population of XRBs grows with time, as black holes and neutron stars form.Yet, the contribution by ionizing radiation from XRB accretion discs to the total line luminosities remains negligible, reaching at most a few percent at ages greater than about 30 Myr.We   conclude that, based on the XRB luminosity function observed in nearby galaxies, XRB accretion discs are unlikely to contribute significantly to the strong He ii emission from metal-poor star-forming galaxies in the sample of Section 3.1.

Emission from radiative shocks
Radiative shocks, such as those expected from massive-star winds and SN blast waves, have spectral signatures characterized by strong highionization lines (e.g., Dopita & Sutherland 1996;Alarie & Morisset 2019), similar to those observed in some metal-poor star-forming galaxies (Izotov et al. 2012;Plat et al. 2019).As noted by Plat et al. (2019), the high gas densities ( H ≳ 10 4 cm −3 ) measured from the [C iii]1907+C iii]1909 doublet in some metal-poor star-forming galaxies could be suggestive of the presence of radiative shocks from SNe.So far, however, no predictive phenomenological model has been proposed to link shocks to other galaxy properties.
In Appendix C, we propose a simple prescription to compute the emission from radiative shocks generated by stellar winds and SN blast waves consistently with the emission from stars and photoionized H ii regions in young star-forming galaxies.This prescription allows line luminosities to be computed from knowledge of the rates of energy injection into the ISM by stellar winds and SN explosions, based on line fluxes tabulated as a function of shock velocity ( s ) and preshock density ( H ) in grids of radiative-shock models (e.g., Dopita & Sutherland 1996;Alarie & Morisset 2019).We estimate the rate at which stellar winds inject energy into the ISM by summing the individual contributions of all stars (ignoring for simplicity the influence of binary-star processes on this rate; see Appendix C).For energy injection by SNe, we consider the rate of SN explosions predicted by SEVN, including pair-instability, type-II and type-Ia SNe.We ignore non-exploding, 'failed' SNe leading to the formation of black holes (Spera et al. 2015;Spera & Mapelli 2017).For reference, the rates of energy injection computed in this way are roughly consistent with those reported by Leitherer et al. (1992, see their fig.6) under slightly different assumptions.
The top-two panels of Fig. 14 show the resulting contributions of radiative shocks driven by stellar winds and SN explosions to the total H and He ii 4686 luminosities of a GALSEVN binary-star population with standard metallicity  = 0.001.These contributions, obtained by inserting line fluxes from the recent radiative-shock calculations of Alarie & Morisset (2019, based on MAPPINGS V) into equations ( C6)-(C8), are shown for the full available range of shock velocities, 100 ⩽  s ⩽ 1000 km s −1 , and  H = 100 cm −3 (Table 1).In each panel, the black dotted line shows the total line emission obtained by adding the maximum possible shock contribution (independently for stellar winds and SNe) to the emission from stars and XRB accretion discs (green solid line).
Fig. 14 shows that radiative shocks in this model account at most for only a few per cent of the total H and He ii 4686 line emission at ages below 10 Myr, with no significant impact on the global He ii 4686/H ratio, except for an extremely brief period around 2 Myr, when shocks from stellar winds can produce around 40 per cent of the total He ii 1640 emission.At ages beyond 10 Myr, the shock contribution increases as the rates of ionizing photons produced by the stars begin to decrease, but it remains with no significant impact on the line luminosities and emission-line ratios.The two peaks in SN-shock contribution originate from PISNe at ages around 4 Myr, and type-II SNe at ages around 10 Myr, type-Ia SNe appearing only after about 30 Myr and in small number.We note that adopting energies greater than 10 51 erg for the most massive PISNe (footnote 11) would increase the magnitude of the first peak, but only for a brief spike at the earliest ages.
It is instructive to also examine the predictions of our model for the luminosity of the [Ne v] 3426 emission line, which has been suggested as a potential signature of radiative shocks in blue compactdwarf galaxies (e.g., Izotov et al. 2012).The bottom panel of Fig. 14 shows that, indeed, radiative shocks driven by stellar winds and SN explosions can significantly contribute to, and even dominate, [Ne v] 3426 emission in metal-poor star-forming galaxies.
Overall, for what concerns the direct focus of the present study, we conclude from Fig. 14 that radiative shocks driven by stellar winds and SN explosions are unlikely to contribute prominently to the strong He ii emission observed in metal-poor star-forming galaxies in the sample of Section 3.1.

DISCUSSION
In Sections 3 and 4 above, we have seen that the new GALSEVN spectral-evolution model of binary-star populations presented in Section 2 predicts significantly harder ionizing radiation than previous models of single-and binary-star populations, making it possible to reproduce, in particular, the high He ii 4686/H ratios and large H equivalent widths commonly observed in metal-poor star-forming galaxies.The most extreme He ii 1640 emitters can be accounted for by invoking populations entirely dominated by massive stars (Section 4.2).The model also allows us to compute, in a physically consistent way, the contribution by accretion discs of X-ray binaries to the ionizing radiation of a star-forming galaxy (Section 4.3).Interestingly, we find that reproducing the observed luminosity function of X-ray binaries in nearby, metal-poor star-forming galaxies implies that XRB accretion discs must have, on average, a negligible impact on the predicted luminosities of high-ionization lines.While consistent with the results from several previous studies (e.g., Jaskot & Oey 2013;Senchyna et al. 2017Senchyna et al. , 2020;;Saxena et al. 2020), this finding contrasts with some recent claims (Schaerer et al. 2019;Umeda et al. 2022;Katz et al. 2023), which we now examine in more detail.Schaerer et al. (2019) used an observed anti-correlation between Xray luminosity per unit star formation rate ( tot X /SFR) and metallicity (Douna et al. 2015;Brorby et al. 2016), combined with an estimate of the rate of He ii-ionizing photons (  He II ) per unit  tot X in the metal-poor star-forming galaxy I Zw 18, to conclude that XRBs are likely to be the main source of He ii emission in metal-poor starforming galaxies.The founding assumption of this reasoning, that He ii emission in I Zw 18 is produced by XRBs, has since been ruled out by Kehrig et al. (2021), who find that the high-mass binary dominating the X-ray emission lies about 200 pc away from the He iiemission peak.Also, the prediction by Schaerer et al. (2019) that the He ii 4686/H line-luminosity ratio should correlate with  tot X /SFR does not appear to be confirmed observationally by the sample of 11 metal-poor star-forming galaxies with high-quality constraints on both X-ray and He ii emission studied by Senchyna et al. (2020, see their fig.9).
In fact, by arbitrarily scaling ARES multicolour-disc models for a wide range of black-hole masses (Mirocha 2014) added to BPASS ionizing spectra of binary-star populations with constant star formation, Senchyna et al. (2020) find that  tot X /SFR in excess of 10 42 erg s −1 /M ⊙ yr −1 is required for XRBs to significantly boost He ii emission.This is 10-100 times larger than observed in the 11 metal-poor star-forming galaxies of their sample, suggesting that XRBs are not primarily responsible for He ii emission in these galaxies.Similarly, Saxena et al. (2020) find no significant difference in the  tot X /SFR ratios of galaxies with and without He ii emission at redshift  ∼ 3 in the Chandra Deep Field South.These results are consistent with our finding that accretion discs of XRBs negligibly affect He ii emission (Fig. 1).Katz et al. (2023) also explore the emission-line properties of gas photoionized by sources with different  tot X /SFR ratio by arbitrarily scaling an ARES multicolour-disc model, for a fixed black-hole mass of 25 M ⊙ , added to the BPASS ionizing spec-trum of a 5 Myr-old binary-star SSP. 5 They find that reproducing the surprisingly high [O iii] 4363/[O iii] 5007 ratio of 0.048 (0.055 when dust-corrected) in the JWST/NIRSpec spectrum of the metal-poor star-forming galaxy S04590 at redshift  = 8.5 requires  tot X,2−10keV /SFR > 10 41 erg s −1 /M ⊙ yr −1 , which translates into  tot X /SFR > 3 × 10 41 erg s −1 /M ⊙ yr −1 when adopting the energy band 0.5-8 keV used by Senchyna et al. (2020) and in Fig. B3 of the present paper.Again, this is much larger than observed in extremely metal-poor star-forming galaxies (Senchyna et al. 2017(Senchyna et al. , 2020)), and than predicted by our physically consistent model of the emission from stars and XRB accretion discs, anchored on the average observed XRB luminosity function of Fig. 12.
It is important to bear in mind that the simple prescription presented in Section 4.3 pertains to the average ⟨ tot X /SFR⟩ ≈ 2.1 × 10 40 erg s −1 /M ⊙ yr −1 produced by XRBs in the five metalpoor star-forming galaxies considered by Lehmer et al. (2019) to build the X-ray luminosity function of Fig. 12, and that the associated standard deviation of 1.4 × 10 40 erg s −1 /M ⊙ yr −1 implies that this ratio can vary from galaxy to galaxy.The 10 starforming galaxies with 12 + log(O/H) < 8.0 in the sample of Douna et al. (2015, see Fig. B3 below) exhibit a very similar mean ⟨ tot X /SFR⟩ ≈ 2.7 × 10 40 erg s −1 /M ⊙ yr −1 with a standard deviation of 2.4 × 10 40 erg s −1 /M ⊙ yr −1 .Yet, the X-ray luminosity per unit star formation rate required by Katz et al. (2023) to significantly impact the [O iii] 4363/[O iii] 5007 ratio of S04590 appears improbably high, as it is, respectively, 20  and 11  larger than the mean observed values for metal-poor galaxies in the samples of Lehmer et al. (2019) and Douna et al. (2015).
Using a slightly different approach, Umeda et al. ( 2022) appeal to Markov Chain Monte Carlo techniques to select the combination of a blackbody and power-law ionizing spectra and gas parameters providing the best fits to the emission-line properties of three extremely metal-poor galaxies with strong He ii 4686 emission (included in the sample of Section 3.1): J1631+4426, J104457 and I Zw 18 NW.Umeda et al. (2022) then show that, for each galaxy, the best-fitting combination of blackbody+power-law ionizing spectra can be approached by a combination of the accretion-disc model spectrum of Gierliński et al. (2009) and the spectrum of a BPASS single-star SSP with an age between about 2 and 8 Myr (we note that the emission-line properties actually produced by these spectra are not compared to the observations in their paper).The  tot X /SFR ratios of these models are 5.0 × 10 40 erg s −1 /M ⊙ yr −1 (J1631+4426), 1.0×10 41 erg s −1 /M ⊙ yr −1 (J104457) and 3×10 41 erg s −1 /M ⊙ yr −1 (I Zw 18 NW), reaching again up to 11  above the mean observed value in the Douna et al. (2015) sample of metal-poor star-forming galaxies. 6 It is of interest to check how the emission-line and X-ray properties of the collection of models shown in Fig. 9, which include selfconsistent scaling of accretion-disc with stellar emission, compare with those of J1631+4426, J104457 and I Zw 18 NW.These 800 models, built using the star formation histories of 100 SPHINX galaxies (Rosdahl et al. 2022), are far from sampling the whole available space of stellar and gas parameters, as they have fixed   =0.17 and 0.44) and four metallicities ( =  ISM = 0.0005, 0.001, 0.002 and 0.004), for ages between 2 and 10 Myr.Therefore, the use of sophisticated fitting algorithms (e.g., BEAGLE, Chevallard & Charlot 2016) is not appropriate in this case, and we simply explore the properties of the models that, from this limited library, are best able to approximate the observed emission-line properties of the three galaxies considered by Umeda et al. (2022).After some experimentation and to better account for the physical and chemical conditions in these galaxies, it appeared opportune to probe higher ionization parameters, carbon-to-oxygen and dust-to-gas ratios, so we extended the library to include a few more values: log ⟨⟩ = −1.7 and −1.4,C/O =0.09 and  d = 0.1.
We present the results of this exercise in Fig. 15, which shows the relative offset between luminosity of the best-fitting (i.e., minimum- 2 ) model and observed luminosity, for the lines in common between the two sets.Remarkably, among the restricted library of models at hand, those shown in Fig. 15 appear to provide reasonable fits to all but one emission line.In particular, the models reproduce the strong observed He ii 4686 luminosities of these galaxies, corresponding to He ii 4686/H ratios of 0.023 (J1631+4426), 0.018 (J104457) and 0.034 (I Zw 18 NW), i.e., totally out of reach of previous stellar population synthesis models (Fig. 3b).
The only line that the limited set of models considered here does not seem to naturally reproduce is [O iii] 4363, whose observed luminosities in two of the three galaxies are higher than predicted by the best-fitting models, especially in the case of J1631+4426.The observed [O iii] 4363/[O iii] 5007 ratios are 0.048 (J1631+4426), 0.032 (J104457) and 0.036 (I Zw 18 NW).Intriguingly, the value for J1631+4426 is the same as that found in the high-redshift galaxy S04590 studied by Katz et al. (2023), for which no simple explanation, other than an improbably high contribution by XRBs (see above) or a high cosmic-ray background (not applicable to nearby galaxies), could be found.Exploring other potential causes of strong [O iii] 4363 emission, such as exceptional heating perhaps related to temperature inhomogeneities (e.g., Cameron et al. 2023b), is beyond the scope of the simple exercise considered here.
We therefore conclude that the prediction from the GALSEVN model presented in this paper, that XRB accretion discs have little influence on the high-ionization line properties of metal-poor starforming galaxies, appears robust when confronted with observations of both the emission-line and X-ray properties of such galaxies.

SUMMARY
We have explored the emission-line properties of interstellar gas photoionized by young binary-star populations computed with the new GALSEVN model built by combining the SEVN (Iorio et al. 2023) and GALAXEV (Bruzual & Charlot 2003) population-synthesis codes.This model allows us to compute in a physically consistent way, for the first time, the emission from stars, accretion discs of Xray binaries and supernova-driven radiative shocks.We compared the predictions of this model with observed UV and optical emission-line properties of a sample of about 120 metal-poor star-forming galaxies in a wide redshift range.We can summarize our main results as follows.
• We confirm that, as well known from previous work (e.g., Eldridge & Stanway 2012;Eldridge et al. 2017;Götberg et al. 2020), the inclusion of binary-star processes (such as envelope stripping, quasi-homogeneous evolution, common-envelope ejection and starstar mergers) boosts and hardens the ionizing emission from young stellar populations.Such effect is achieved assuming even a relatively low fraction of binary stars, the rates of H-, He i-and He iiionizing photons changing only little for stellar populations with 0.3 ≲  bin ⩽ 1.0 (Fig. 1).
• We find that the He ii 4686/H ratios and H equivalent widths of GALSEVN SSP models with metallicity  ∼ 0.001 reproduce remarkably well the observations of galaxies in our sample at all ages between 3 and 10 Myr, providing better agreement with the data than achievable with previous single-star and binary-star models (Fig. 3b).
• The inclusion of bursty star formation histories, such as those predicted by the SPHINX simulations of young galaxies in the epoch of reionization (Rosdahl et al. 2022), allows broader sampling of the observational space than achievable with simpler models, and hence, should be encouraged when interpreting observations (Figs 8-9).
• We find that the most extreme He ii 1640 equivalent widths (and other high-ionization signatures) of galaxies in our sample can be reproduced by models with a top-heavy IMF with characteristic mass ∼ 50 M ⊙ , which could represent a transient phase in the evolution of a galaxy (Fig. 11).
• Reproducing the average luminosity function of X-ray binaries per unit star formation rate in the most metal-poor star-forming galaxies observed by Lehmer et al. (2019) requires that the bulk of mass in such systems be accreted during episodes of X-ray outbursts, at a much higher rate than during the quiescent phase, the actual strength and duration of these outbursts being degenerate (Appendix B).This implies that XRBs contribute negligibly to the H-and He ii-line emission from young star-forming galaxies, in agreement with expectations from several previous studies (e.g., Jaskot & Oey 2013;Senchyna et al. 2017Senchyna et al. , 2020;;Saxena et al. 2020, see Fig. 13).We find that the claims in some recent studies that XRB accretion discs might contribute significantly to the emission from high-ionization UV and optical lines (Schaerer et al. 2019;Umeda et al. 2022;Katz et al. 2023) are based on models predicting improbably high ratios of X-ray luminosity to star formation rate (Section 5).
• Finally, by combining the rates of energy injection into the ISM by stellar winds and (pair-instability, type-II and type-Ia) SN explosions predicted by the GALSEVN model with the radiative-shock models of Alarie & Morisset (2019, see also Dopita & Sutherland 1996), we find that shocks are unlikely to contribute significantly to the H-and He ii-line emission from young star-forming galaxies (Fig. 14).
The work presented in this paper provides only first examples of the properties of the new GALSEVN model of the emission from binary-star populations.The application of this model to interpret the spectral properties of young star-forming galaxies using Bayesian inference (e.g., with the BEAGLE tool) requires building a comprehensive library spanning full ranges of stellar and nebular parameters and the exploration of a comprehensive set of emission lines.This will be the subject of a forthcoming paper.

APPENDIX A: HANDLING OF BINARY-EVOLUTION PRODUCTS IN SEVN
The SEVN code evolves stars in a binary system by interpolating from a set of precomputed single-star tracks (in this work, we use the PAR-SEC track library; see Section 2.1).The evolution of a star involves interpolation from four tracks bracketing its mass and metallicity (see Iorio et al. 2023 for details).However, binary processes, such as wind mass transfer, Roche-Lobe overflow and stellar mergers, can significantly alter the properties of a star.In such cases, SEVN searches for new tracks to interpolate from to better match the current stellar properties.This is achieved through different strategies tailored to specific situations, as follows.
• Mass loss/accretion: for main-sequence stars, SEVN checks whether the net cumulative mass variations due to binary processes exceed 1 per cent of the current stellar mass.In such cases, SEVN searches for new interpolating tracks providing a better match to the total mass of the star at the same percentage of stellar life.The evolution of stars with decoupled He and CO cores is driven by core properties (Hurley et al. 2002).For this reason, we do not allow stars outside the main sequence to change stellar tracks unless the core mass has changed (e.g., during stellar mergers) or the star has completely lost its envelope.
• Envelope stripping: if a star completely loses its envelope, SEVN switches to pure-He stellar-evolution tables, selecting the pure-He interpolating tracks that best match the mass of the bare He core at the same percentage of stellar life.In this work, we use the PARSEC pure-He stellar-evolution tables described in Iorio et al. (2023).When a pure-He star loses its He envelope, SEVN models the bare CO core as a remnant: if the core is massive enough to produce a black hole or neutron star, its properties remain constant until remnant formation, otherwise it immediately becomes a white dwarf.
• Stellar merger: depending on the properties of the merging stars, SEVN employs different strategies to select the interpolating tracks for the merger product.If two main-sequence stars merge, SEVN changes stellar evolution tracks similarly to cases of mass loss or accretion.
Instead, if at least one of the merging stars is evolved, SEVN searches for new interpolating tracks by matching the core mass at the same percentage of stellar life.For mergers between a 'standard' star and a stripped, pure-He star, SEVN searches for new interpolating tracks in the H-rich tables matching the new core mass.We assume that mergers between stars and remnants (white dwarfs, neutron stars, black holes and bare-CO cores) completely destroy the star, leaving only the remnant, and no mass is accreted on to the remnant.
Additional details about the strategy for selecting interpolating tracks of binary-evolution products in SEVN can be found in Iorio et al. (2023).

APPENDIX B: EMISSION FROM XRB ACCRETION DISCS B1 Context
We identify X-ray binaries among a population of evolving binary pairs generated with the SEVN code (Section 2.1) by searching for black holes and neutron stars whose masses increase with time, i.e., with positive mass accretion rate across a dynamical time step,  acc > 0. For spherical infall, the maximum luminosity that can theoretically be produced by an accretion disc in hydrostatic equilibrium is the Eddington luminosity corresponding to the balance between radiation pressure and gravity (Cameron & Mock 1967, see also Hurley et al. 2002).For accretion of fully-ionized H+He material, with hydrogen mass fraction X, on to a compact object of mass  c , this limit can be expressed as The Eddington luminosity is related to an Eddington accretion rate,  Edd , through the formula (e.g., Shakura & Sunyaev 1973) where  is the radiative efficiency (i.e. the fraction of binding energy radiated away by the accretion flow) and  the speed of light.
Observationally, accretion on to compact objects is known to occasionally produce luminosities in excess of the theoretical Eddington limit, perhaps as a result of radiation-driven inhomogeneities (e.g., Begelman 2002).Examples include tidal disruption events (e.g., Rees 1988) and ultra-luminous X-ray sources (e.g., Bachetti et al. 2014;Pinto et al. 2016;Rodríguez Castillo et al. 2020), while some evidence suggests that most low-mass XRBs may have undergone a phase of super-Eddington accretion (Kalogera & Webbink 1998).For super-Eddington accretion, corresponding to  acc >  Edd , the total accretion luminosity, noted  acc , can exceed the Eddington luminosity by a logarithmic factor ∼ ln(  acc /  Edd ) (see Shakura & Sunyaev 1973;Begelman et al. 2006).We adopt here the convenient approximation (Watarai 2006, see also Lapi et al. 2014) which, for sub-Eddington accretion (  acc ≲  Edd ), reduces to For super-Eddington accretion, the weaker-than-linear increase of  acc with  acc in equation (B3) (illustrated by fig.7 of Watarai 2006) reflects how the accreting gas becomes optically too thick to radiate away all the dissipated energy, causing some radiation to be trapped and advected inward.For simplicity, we fix here the radiative efficiency to the standard value  ≈ 0.1 corresponding to the fraction of energy liberated by a gas particle falling from far away on to the innermost stable circular orbit of a stationary (non-spinning) black hole (e.g., Shakura & Sunyaev 1973).
In SEVN, the maximum rate at which a compact object can accrete mass is  max acc =  Edd  Edd , where  Edd is an adjustable factor (Spera  et al. 2019; Iorio et al. 2023). 7This factor can be set to a value greater than unity to account for super-Eddington accretion.In practice, we compute the mass accretion rate as  c /(1 − ), where  c is the mass-growth rate of the compact object recorded in SEVN.In the setup of the code used here, the mass lost by the donor that is not accreted by the compact object (if  acc >  max acc ) is assumed to be lost from the system (other options are described in Iorio et al. 2023).

B2 Spectral modelling
Observations of X-ray binaries show that these stars commonly undergo X-ray outbursts, during which their luminosity can exceed that in the quiescent state by up to several orders of magnitude (e.g., Tanaka & Shibazaki 1996;Chen et al. 1997).While this transience phenomenon is often discussed in the context of low-mass XRBs, where the donor star is typically less massive than a few M ⊙ (e.g., Yan & Yu 2015), it is also observed in high-mass XRBs (e.g., Martin et al. 2014;van den Eĳnden et al. 2022).Even high-mass XRBs classified as 'persistent' (i.e., non-transient sources emitting X-rays continuously over long periods) exhibit variability accompanied by (sometimes giant) X-ray flares (e.g., Hertz et al. 1992;Pottschmidt et al. 2003;Fürst et al. 2010).Such X-ray outbursts appear to be correlated with variability at all wavelengths from ultraviolet to infrared (e.g., Hynes et al. 2003;Sonbas et al. 2019;López-Navas et al. 2020;Yang et al. 2022).
The physical origin of X-ray variability is unclear.For several decades, the so-called disc-instability model has been put forward, according to which, as material falls on to the disc, a thermal-viscous instability develops, leading to recurrent burst-accretion events (e.g., Meyer & Meyer-Hofmeister 1981;Faulkner et al. 1983;van Paradĳs 1996;Hameury & Lasota 2020, see the reviews by Lasota 2001; Hameury 2020).While the realistic incorporation of thermodynamics and the ability to reproduce fairly well observations of X-ray outbursts have made this model popular, it does suffer from several weaknesses and shortcomings, notably with regard to angularmomentum transport and the inclusion of several poorly-constrained parameters (see Hameury 2020 for details).Hence, any attempt to implement the disc-instability model in SEVN would be highly uncertain.
In this context, we prefer to adopt here a more empirical approach to account for variability when modelling the spectral properties of XRB accretion discs.We assume that any accretion event spreads over two phases: a quiescent phase during which accretion proceeds at the rate provided by SEVN,  c /(1 − ), and an outburst phase during which it is boosted by a factor  A .We assume that a fraction  A of the mass is accreted during the outburst phase.For the accreted mass  c to be conserved over the time interval  ′ , we thus write where  c is the mass growth of the compact object in a time interval  ′ of the evolution of an SSP (adopting the same notation as in equation 1 of Section 2).For reference,  c represents at most a few percent of  c .We describe the spectral distribution of the power  acc ( ′ ) produced by the hot accretion disc at time  ′ (equation B3) by means of a multicolour-disc model (Mitsuda et al. 1984), i.e., a modified blackbody spectrum reflecting the gas-temperature distribution in the disc.To this end, we use the ARES8 code of Mirocha (2014) and compute a comprehensive library of multicolour disc models for compact objects accreting at the Eddington limit, parametrized in terms of the accreting-object mass,  ARES , and the maximum radius of the accretion disc,  max .Following Senchyna et al. (2020), we fix  max to 10 4 gravitational radii (we find that decreasing or increasing  max by an order of magnitude would, respectively, lower by about 20 per cent or increase negligibly the X-ray luminosities and rates of He ii-ionizing photons predicted by our models).By design, each spectrum is normalized to the Eddington luminosity given by equation (B1) for  c =  ARES , for a pure-H gas (X = 1).
XRB spectra also often present a high-energy tail component, whose strength can vary over time, and which is thought to arise from Compton upward scattering of soft disc photons by energetic coronal electrons (e.g., Remillard & McClintock 2006;Steiner et al. 2009).Based on the analysis of 20 XRBs with multi-epoch X-ray luminosities in the range ∼ 0.5-50 × 10 39 erg s −1 , Sutton et al. (2013) conclude that the spectra of sources with sub-Eddington accretion rates may be well represented by multi-colour discs, while those of sources with super-Eddington accretion rates tend to exhibit 'hardluminous' and 'soft-luminous' components (in roughly equal proportions among the sample), presumably depending on inclination.For each  ARES , we therefore compute both a spectrum including Comptonization of disc photons (implemented in ARES using the SIMPL model by Steiner et al. 2009) and a spectrum not including it.If accretion is sub-Eddington, i.e. for  acc ( ′ ) ⩽  Edd , we adopt the pure multi-colour-disc spectrum.Otherwise, we randomly draw between the pure multi-colour-disc and the Comptonized spectra. 9 We note that the results presented in this paper about the X-ray luminosities and He ii-ionizing photon rates of star-forming galaxies do not depend sensitively on this refinement of spectral modelling at the highest energies. 9Since the ARES spectral library pertains to objects accreting at the Eddington rate, in practice, for a given  acc , we adopt the spectrum corresponding to the ARES model for a compact object accreting at the Eddington luminosity  Edd =  acc given by equation (B1) for a pure-H gas, i.e.,  ARES =  acc /(1.26 × 10 38 erg s −1 ) M ⊙ .
The above approach allows us to compute the spectrum of every hot accretion disc of XRB in a GALSEVN binary-star SSP of age  ′ , from which we extract the disc luminosity in the 0.5-8 keV energy band, noted  X .For reference, the quantity  X amounts to typically between 30 and 90 per cent of  acc .Fig. B1 shows X-ray luminosity functions (XLFs) of XRBs obtained in this way for populations with binary-star fraction  bin = 0.7 (typical of massive-star populations; e.g., Sana et al. 2012), assuming 100 Myr of constant star formation, for different choices of the parameters  A and  A in equation (B5).All models have the metallicity  = 0.008, corresponding to the metallicity of the five most metal-poor star-forming galaxies observed with Chandra for which Lehmer et al. (2019) derived the average, completeness-corrected XLF of XRBs shown by the black data points with error bars (normalized per unit star formation rate in the same way as the models).The black dotted line shows an empirical model, incorporating also potential contamination by background point sources from the cosmic X-ray background, whose contribution to the XLF is however estimated to be minor, decreasing from about 20 per cent at  X = 10 37 erg s −1 to about 3 per cent at  X = 10 40 erg s −1 (for details, see Lehmer et al. 2019).We therefore ignore this component in our model.Furthermore, we adopt here  Edd = 100 to set the maximum allowed accretion rate in SEVN (Section B1), as significantly lower values do not produce enough high-luminosity sources, while adopting higher values negligibly affects the predicted XLFs, indicating that only a few systems can potentially have  acc > 100  Edd .
We find that the predicted XLF depends in a somewhat degenerate way on the parameters  A and  A of the two-phase model, the main requirement for reproducing the observations being that the bulk of mass be accreted during outbursts.This is illustrated by the three models shown in Fig. B1, accreting, during transient outbursts, fractions  A = 1.00, 0.99 and 0.95 of the mass at  A = 100, 100 and 1000 times the quiescent rate, respectively.In the first model, the entire mass is accreted during outbursts, while in the second model, outburst-and quiescent-accretion phases have similar durations, and in the third model, outburst accretion lasts less than 1/50 th of quiescent accretion.All three models provide roughly similar agreement with the observed XLF of XRBs compiled by Lehmer et al. (2019) in Fig. B1.The maximum effective accretion rate for these models is  A  Edd = 10 4 -10 5 times the Eddington rate  Edd , corresponding to a maximum accretion luminosity  acc ≈ 17-22 Edd (equation B3), compatible with the 10-100 Edd obtainable from radiation-driven inhomogeneities in accretion discs (e.g., Begelman 2002).We note that about a third of XRBs exhibit super-Eddington accretion rates during the phase of standard accretion, compared to three quarters during the phase of enhanced accretion (for  A = 100).
The model with a single phase of enhanced accretion is the simplest conceptually.The one with  A = 0.99 and  A = 100 provides formally the best fit (mininum  2 ) to the data, while the one with quisecent accretion ∼ 50 times longer than enhanced accretion is the closest conceptually to that proposed by, e.g., Li (2012) for the growth of supermassive black holes in quasars.Since all three models reproduce the observed XLF of XRBs measured by Lehmer et al. (2019) in nearby metal-poor galaxies, they predict similar contributions of XRB accretion discs to the ionizing-radiation budget of a star-forming galaxy modelled with the GALSEVN code.Here, we adopt the model with  A = 0.99 and  A = 100, which provides the best fit to the observations in Fig. B1.Any of the other two models would provide similar results on the emission from XRB accretion discs presented in this paper.
In Fig. B2, we show the XLFs obtained at fixed  A = 100 and  A = 0.99, for GALSEVN binary-star populations of different metal-  licities in the range 0.0001 ⩽  ⩽ 0.02.Remarkably, all XLFs remain compatible with the observational determination of Lehmer et al. (2019).This is because the XLF is largely dominated by systems composed of low-mass compact objects with relatively lowmass companions: nearly 80 per cent of overall accretion occurs in systems composed of a compact object lighter than 20 M ⊙ (with progenitor main-sequence mass typically smaller than 30 M ⊙ ; e.g., Spera et al. 2015) and a companion lighter than 25 M ⊙ .Since the evolution of stars lighter than 30 M ⊙ is not significantly affected by line-driven stellar winds, we do not expect a strong dependence of XRB properties on metallicity.The main metallicity-dependent effect comes from the reduced population of massive compact objects at high metallicity (e.g., Dray 2006;Mapelli et al. 2013), apparent in the slight decline in the XLF at the highest X-ray luminosities in Fig. B2, while agreement with the observed luminosity function is maintained.
In Fig. B3, we compare the total X-ray luminosities per unit star formation rate,  tot X /SFR, obtained by integrating the XLFs of the binary-star models with different metallicities (converted to gasphase O abundances using table 2 of Gutkin et al. 2016 for  d = 0.3) of Fig. B2 with those of the sample of 10 nearby star-forming galaxies with 12 + log(O/H) < 8.0 compiled by Douna et al. (2015, green points ; see also Brorby et al. 2014) and the sample of nearby analogues of distant Lyman-break galaxies compiled by Brorby et al. (2016).Also shown is the mean  tot X /SFR (and associated standard deviation) of the five metal-poor star-forming galaxies used by Lehmer et al. (2019) to derive the average XLF in Fig. B1. 10 As expected from Fig. B1, the model with  = 0.008 (corresponding to 12 + log(O/H) ≈ 8.4) lies near this mean value.More generally, Fig. B3 shows that the total X-ray luminosity per unit star formation rate predicted by our models is well within the ranges found by Douna et al. (2015, green points) and Lehmer et al. (2019, light blue points) and near the upper end of that reported by Brorby et al. (2016, dark blue points).
We consider the reasonable agreement between GALSEVN models of binary-star populations and X-ray observations of nearby, metalpoor star-forming galaxies and analogues of distant Lyman-break galaxies in Figs B2-B3 as support for our simple approach to model the emission from accretion discs of X-ray binaries consistently with the spectral modelling of stars.

Figure 1 .
Figure 1.Rates of H-(top), He i-(middle) and He ii-(bottom) ionizing photons as a function of age for 'simple' (i.e., instantaneous-burst) stellar populations (SSPs) with different fractions of binary stars,  bin (colour-coded as indicated), for a metallicity  = 0.001.The models assume a Chabrier (2003) IMF for primary stars and a Sana et al. (2012) distribution of secondary-toprimary mass ratios at age zero (single stars being modelled as non-interacting binaries) and are normalized to a total initial stellar mass of 1 M ⊙ integrated over 0.1-300 M ⊙ (see text for details).

Figure 2 .
Figure 2. Spectral energy distribution   entering equation (1) (in units of the luminosity per unit frequency at the Lyman limit) for the SSP models of Fig.1with pure binary stars (  bin = 1, solid lines) and pure single stars (  bin = 0, dotted lines) at the ages  ′ = 1, 3, 5, 8 and 10 Myr (colour-coded as indicated on the right-hand side).Also shown in grey is the area sampled by ionizing spectra of AGNs with power-law indices between  AGN = −2.0(bottom edge) and −1.2 (top edge;Feltre et al. 2016).The ionization energies of some common species are indicated for reference.

Figure 3 .
Figure 3. Optical and ultraviolet emission-line properties of the reference observational sample of star-forming galaxies (filled grey circles: non LyC leakers; filled grey stars: LyC leakers) and AGNs (open grey diamonds) described in Section 3.1.The diagrams show different combinations of equivalent widths and ratios of the C iv 1549, He ii 1640, O iii] 1664 and C iii] 1908, [O iii] 5007, [O ii]3727, [O i] 6300, He ii 4686, H, H and [N ii] 6584 nebular emission lines, and the gas-phase oxygen abundance in (d).All line fluxes are corrected for attenuation by dust, as prescribed in the original studies (arrows show 1 upper limits).In (a), the dotted and dashed lines show the criteria of Kewley et al. (2001) andKauffmann et al. (2003), respectively, to separate AGN-dominated from star-forming galaxies, while in (e), the dashed line shows the equivalent criterion proposed byNakajima et al. (2018).The curves interspersed with filled squares, triangles, pentagons, hexagons and crosses show the predictions of the GALSEVN binary-and single-star models, the BPASS binary-and single-star models and the C&B single-star model, respectively, colour-coded with age as indicated on the right (the BPASS models, available only at ages above 1 Myr, do not sample the youngest ages in the colour bars).In (d), the models are shown for clarity at slightly different abscissae, although all have 12 + log(O/H) = 7.53.In all models, the ionizing stellar population is an SSP with metallicity  = 0.001 and the same zero-ageChabrier (2003) IMF as in Fig.1, while the parameters of the photoionized gas are those listed in Table1.

Figure 4 .
Figure 4. Rates of H-(top) and He ii-(bottom) ionizing photons as a function of age for SSPs of metallicity  = 0.001 computed with different models, as indicated (the BPASS models are available only at ages above 1 Myr).In each panel, the grey shaded area shows the area covered by the GALSEVN models with different fractions of binary stars from Fig. 1.All models have the same zero-age Chabrier (2003) IMF normalized to a total initial stellar mass of 1 M ⊙ integrated over 0.1-300 M ⊙ as in Fig. 1.

Figure 5 .
Figure 5. Same as Fig. 3, but for GALSEVN binary-star models with different metallicities, as indicated at the tip of the colour bars.The curves interspersed with filled squares, triangles, crosses and pentagons refer to the metallicities  = 0.0005, 0.001, 0.002 and 0.004, respectively.

Fig. 3 (
Fig. 3 (with  =  ISM = 0.001) with two new GALSEVN binary-star models: one exploring the influence of super-solar O/Fe by setting the stellar metallicity to 10 per cent of the ISM metallicity (i.e.,  = 0.0001 versus  ISM = 0.001) and the other investigating the impact of a global reduction of metallicity, with  =  ISM = 0.0001.

Figure 7 .
Figure 7. Same as Fig. 3, for three GALSEVN binary-star models: a model with the standard parameters of Table 1 (with  =  ISM = 0.001; filled squares), one exploring the influence of super-solar O/Fe (with  = 0.0001 and  ISM = 0.001 ; filled triangles), and one investigating the impact of a global reduction of metallicity (with  =  ISM = 0.0001 ; filled crosses).See main text for details.In (d), the model with  ISM = 0.0001, which corresponds to 12 + log(O/H) = 6.53, is shown at the plot limit of 12 + log(O/H) = 6.80 to retain the same scale as in Fig. 3.
which approaches the standardChabrier (2003) IMF at large  but is exponentially cut off at  <  char .Wise et al. (2012) used this function with  char = 40 M ⊙ to explore the formation of Population III stars in early dwarf galaxies, whileGoswami et al. (2022) adopted 200 ⩽  char ⩽ 300 M ⊙ to explore chemical enrichment by PISNe in extremely metal-poor galaxies.Here, we consider values of  char in the range 50-100 M ⊙ to explore the emission-line signatures of stellar populations strongly dominated by massive stars.The upper panel of Fig.10shows two versions of the top-heavy IMF of equation (

Figure 8 .
Figure 8. Same as Fig. 3, for three GALSEVN binary-star models with metallicity  = 0.001 and the standard nebular parameters of Table 1: a plain SSP model (in green, with filled squares); one with a second burst of star formation of same amplitude as the first one added at the age of 3 Myr (in orange, with filled triangles); and one with the star formation history of a typical galaxy experiencing several bursts of different amplitudes at ages up to 10 Myr in the SPHINX cosmological simulations of Rosdahl et al. (2022, in green, with filled crosses; see text for details).

Figure 9 .
Figure9.Same as Fig.3, but for 800 GALSEVN models computed using the star formation histories of 100 SPHINX galaxies fromRosdahl et al. (2022), combined alternatively with four metallicities ( =  ISM = 0.0005, 0.001, 0.002 and 0.004) and two C/O abundance ratios (C/O =0.17 and 0.44), at ages between 2 and 10 Myr (younger models are not shown to avoid the artificial peak in model properties at early ages).The 2-dimensional histograms of models are colour-coded according to the density scale shown on the right.

Figure 10 .
Figure 10.Top: top-heavy IMFs computed using equation (3) for two values of  char (red: 50 M ⊙ ; salmon: 100 M ⊙ ), compared to the standard Chabrier (2003) IMF (in black).All IMFs are normalized to a total mass in stars of 1 M ⊙ .Bottom: absolute ultraviolet AB magnitude at 1500 Å,  UV , plotted against age for three GALSEVN binary-star models computed with the three IMFs shown in the top panel, assuming a total initial mass in stars of 10 6 M ⊙ , metallicity  = 0.001 and the standard nebular parameters of Table1.For reference, the dotted horizontal line marks the magnitude  UV = −14 mag typical of nearby, extremely metal-poor star-forming galaxies (reached by the top-heavy-IMF models at the ages marked by the cyan stars).

Figure 11 .
Figure 11.Same as Fig. 3, for the three GALSEVN binary-star models with different IMFs shown in the lower panel of Fig. 10 (colour-coded as indicated on the right).The curves interspersed with filled squares, triangles and crosses refer to the models with a Chabrier (2003) IMF and two top-heavy IMFs (equation 3) with characteristic masses of  char = 50 and 100 M ⊙ , respectively.Cyan stars mark the points beyond which the models with top-heavy IMFs become fainter than  UV = −14 mag and are not shown (see text and Fig. 10 for details).

Figure 12 .
Figure 12.X-ray luminosity function (in the 0.5-8 keV energy band) of XRBs for a GALSEVN binary-star population of metallicity  = 0.008 with binary fraction  bin = 0.7, assuming 100 Myr of constant star formation at the rate  = 1M ⊙ yr −1 .Accretion is assumed to proceed according to the two-phase model described by equation (B5), with  A = 100 and  A = 0.99.Also shown for comparison is the average luminosity function derived from Chandra observations of the five most metal-poor star-forming galaxies in the sample studied by Lehmer et al. (2019, data points with error bars).The black dotted line shows an empirical model proposed by these authors (see Section 4.3 for details).

Figure 13 .
Figure 13.Contribution by XRB accretion discs (in blue), compared to the stellar contribution (in green), to the total luminosities (in black) of H (top panel) and He ii 4686 (bottom panel) plotted against age, for a GALSEVN binary-star SSP with metallicity  = 0.001.The model has the same zeroage Chabrier (2003) IMF normalized to a total initial stellar mass of 1 M ⊙ integrated over 0.1-300 M ⊙ as in Fig. 1.

Figure 14 .
Figure 14.Contribution from shocks driven by stellar winds (in red) and SNe (in blue) to the total H (top panel), He ii 4686 (middle panel) and [Ne v] 3426 (bottom panel) luminosities of a GALSEVN binary-star SSP with metallicity  = 0.001, plotted against age (see text Section 4.4 for details).The contributions are shown for the range of shock velocities 100 ⩽  s ⩽ 1000 km s −1 and the density  H = 100 cm −3 .In each panel, the black dotted line shows the total line emission obtained by adding the maximum possible shock contribution (independently for stellar winds and SNe) to the emission from stars and XRB accretion discs (green solid line).The model has the same zero-age Chabrier (2003) IMF normalized to a total initial stellar mass of 1 M ⊙ integrated over 0.1-300 M ⊙ as in Fig. 1.
H = 100 cm −3 , log ⟨⟩ = −2 and  d = 0.3, only two C/O abundance ratios (C/O 5 To compute  tot X /SFR, Katz et al. (2023) define an 'effective' star formation rate based on the rate of H-ionizing photons emitted by the 5 Myr-old BPASS SSP (H.Katz, private communication). 6To compute  tot X /SFR, Umeda et al. (2022) define an 'effective' star formation rate based on the 1500 Å luminosity produced by the BPASS single-star SSP of their best-fitting model (using the conversion in Kennicutt 1998).

Figure 15 .
Figure 15.Relative offset between line luminosities of the best-fitting (i.e., minimum- 2 ) GALSEVN model among the limited library described in Section 5, and observed line luminosities in the three metal-poor galaxies J1631+4426, J104457 and I Zw 18 NW (from top to bottom) studied by Umeda et al. (2022).All luminosities are normalized to that of H.The 1- (3-) observational errors are shown in red (gold).For J1631+4426, the [O iii] 4363 point falls outside the plot, at ( mod −  obs )/ obs ≈ −0.50, as indicated by the down-pointing black triangle.

Figure B1 .
Figure B1.X-ray luminosity functions of XRBs for GALSEVN populations with binary fraction  bin = 0.7, assuming 100 Myr of constant star formation at the rate  = 1M ⊙ yr −1 .Accretion is taken to proceed according to the two-phase model of equation (B5), for different choices of the fraction  A of mass accreted while the rate is boosted by a factor  A , as indicated.All models have the metallicity  = 0.008, corresponding to the metallicity of the five most metal-poor star-forming galaxies observed with Chandra for which Lehmer et al. (2019) derived the average, completeness-corrected luminosity function shown by the black data points with error bars.The black dotted line shows an empirical model proposed by these authors (see text for details).

Figure B2 .
Figure B2.Same as Fig.B1, for GALSEVN binary-star populations of different metallicities in the range 0.0001 ⩽  ⩽ 0.02, as indicated.In all models, accretion is taken to proceed according to the two-phase model described by equation (B5), with  A = 100 and  A = 0.99.

Figure B3 .
Figure B3.Total X-ray luminosity  tot X (per unit star formation rate) obtained by integrating the luminosity functions of Fig. B2, for metallicities in the range 0.0005 ⩽  ⩽ 0.02 (red crosses).The data are from Chandra observations of a sample of 10 nearby star-forming galaxies with 12 + log(O/H) < 8.0 compiled by Douna et al. (2015, green points ; see also Brorby et al. 2014) and a sample of nearby analogues of distant Lyman-break galaxies compiled by Brorby et al. (2016, dark blue points).The light-blue point (and error bars) at 12 + log(O/H) ≈ 8.4 shows the mean  tot X /SFR (and associated standard deviation) of the five metal-poor star-forming galaxies used by Lehmer et al. (2019) to derive the average XLF of Fig. B1.Dotted horizontal green and blue lines represent the mean  tot X /SFR of galaxies in the Douna et al. (2015) and Brorby et al. (2016) samples, respectively.Brorby et al. (2016) estimate the typical uncertainty on 12 + log(O/H) at 0.14.

Table 1 .
Nebular parameters of the 'standard' models of Section 3.2.