<figure><img src="image_1.png" /></figure>

HAL
open science

# Nebular emission from young stellar populations including binary stars

M Lecroq, S Charlot, A Bressan, G Bruzual, G Costa, G Iorio, M Spera, M

Mapelli, Y Chen, J Chevallard, et al.

► To cite this version:

M Lecroq, S Charlot, A Bressan, G Bruzual, G Costa, et al. Nebular emission from young stellar populations including binary stars. Mon.Not.Roy.Astron.Soc., 2024, 527 (3), pp.9480-9504. 10.1093/mn-ras/stad3838. hal-04352423

HAL Id: hal-04352423

https://hal.science/hal-04352423v1

Submitted on 28 Dec 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_3.png" /></figure>

# Nebular emission from young stellar populations including binary stars

Marie Lecroq<sup>1,1*</sup> Stéphane Charlot<sup>1,1*</sup> Alessandro Bressan,<sup>2,3</sup> Gustavo Bruzual,<sup>4</sup>

Guglielmo Costa<sup>1</sup>,<sup>3,5,6,7</sup> Giuliano Iorio<sup>1</sup>,<sup>3,5,6</sup> Mario Spera<sup>1,2,8,9</sup> Michela Mapelli<sup>1</sup>,<sup>3,5,6,10</sup> Yang Chen,<sup>11,12</sup>

Jacopo Chevallard<sup>10,13</sup> and Marco Dall'Amico<sup>5,6</sup>

<sup>1</sup>Sorbonne Université, CNRS, UMR 7095, Institut d'Astrophysique de Paris, 98 bis bd Arago, F-75014 Paris, France

<sup>2</sup>SISSA, via Bonomea 265, I-34136 Trieste, Italy

<sup>3</sup>INAF, Osservatorio Astronomico di Padova, Vicolo dell'Osservatorio 5, I-35122 Padova, Italy

<sup>4</sup>Instituto de Radioastronomía y Astrofísica, UNAM, Campus Morelia, Michoacán, México, C.P. 58089, México

<sup>5</sup>Dipartimento di Fisica e Astronomia Galileo Galilei, Università di Padova, Vicolo dell'Osservatorio 3, I-35122 Padova, Italy

<sup>6</sup>INFN-Padova, Via Marzolo 8, I-35131 Padova, Italy

<sup>7</sup>Univ Lyon, Univ Lyon1, ENS de Lyon, CNRS, Centre de Recherche Astrophysique de Lyon UMR5574, F-69230 Saint-Genis-Laval, France

<sup>8</sup>INAF, Osservatorio Astronomico di Roma, Via Frascati 33, I-00040 Monteporzio Catone, Italy

<sup>9</sup>INFN-Trieste, Via Valerio 2, I-34127 Trieste, Italy

<sup>10</sup>Universität Heidelberg, Zentrum für Astronomie, Institut für Theoretische Astrophysik, Albert-Ueberle-Str 2, D-69120 Heidelberg, Germany

<sup>11</sup>School of Physics and Optoelectronic Engineering, Anhui University, Hefei 230601, China

<sup>12</sup>National Astronomical Observatories, Chinese Academy of Sciences, Beijing 100101, China

<sup>13</sup>Department of Physics, University of Oxford, Denys Wilkinson Building, Keble Road, Oxford OX1 3RH, UK

Accepted 2023 December 7. Received 2023 November 29; in original form 2023 September 21

## ABSTRACT

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 He II λ4686/H β ratios commonly observed at high H β 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 He II λ1640 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 He II λ4686/H β 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.

**Key words:** binaries: general – galaxies: general – galaxies: high-redshift – galaxies: ISM – X-rays: binaries.

## 1 INTRODUCTION

The 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 characterize and analyse the ultraviolet and optical spectra of metal-poor, 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 lower ionization lines (e.g. Shirazi & Brinchmann 2012; Nanayakkara et al. 2019; Plat et al. 2019; Schaerer, Fragos & Izotov 2019; Stanway & Eldridge 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, quasi-homogeneous evolution, and common-envelope ejection), often

© The Author(s) 2023.

Published by Oxford University Press on behalf of Royal Astronomical Society. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.

<footnote>*E-mail: lecroq@iap.fr (ML); charlot@iap.fr (SC)</footnote>
neglected in stellar population synthesis models (e.g. Eldridge & Stanway 2012; Eldridge et al. 2017; Göttberg 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, Thuan & Privon 2012; Nakajima et al. 2018; Plat et al. 2019; Schäerer, Fragos & Izotov 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) population-synthesis codes, coupled with the CLOUDY photoionization code (Ferland et al. 2017). This spectral-synthesis model, called GALSEVN, 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 studied 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.

## 2 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, Charlot & Bruzual 2016) and express the luminosity per unit frequency $v$ emitted at time $t$ by a star-forming galaxy as

$$L_v(t) = \int_0^t dt' \psi(t-t') S_v[t', Z(t-t')] T_v(t, t'), \quad (1)$$

where $\psi(t-t')$ is the star formation rate at time $t-t'$, $S_v[t', Z(t-t')]$ the luminosity produced per unit frequency per unit mass by a single generation of stars of age $t'$ and metallicity $Z(t-t')$, and

$T_v(t, t')$ the transmission function of the interstellar medium (ISM). We compute $S_v$ 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 $T_v$ using the CLOUDY photoionization code (Ferland et al. 2017, Section 2.2).

### 2.1 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, Mapelli & Bressan 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 N-body) is an open-source binary population-synthesis code,<fnref n="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 (PaAdova and TRIeste Stellar Evolution Code) library of evolutionary tracks (including pre-main sequence evolution) for stars with initial masses between 2 and 600 M$_{\odot}$ and metallicities in the range $10^{-11} \le Z \le 0.04$ (Bressan et al. 2012; Chen et al. 2015; Costa et al. 2019, 2023; Nguyen et al. 2022; Santoliquido et al. 2023).

For stars ending their lives as 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 $\alpha_{CE} = 3$ (section 2.3.3 of Iorio et al. 2023, see also Hurley, Tout & Pols 2002).

Finally, in this work, we adopt the formalism of quasi-homogeneous evolution introduced by Eldridge, Langer & Tout

<footnote n="1">¹https://gitlab.com/sevncodes/sevn. The SEVN version used in this work is the release Iorio22 (https://gitlab.com/sevncodes/sevn/-/releases/iorio22).</footnote>
(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, Langer & Tout (2011, see also Eldridge & Stanway 2012), we include this evolution for stars with metallicity $Z \le 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_\odot$. 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 $Z$, we produce with SEVN a stochastic population of $10^6$ evolving binary pairs,$^2$ with initial primary-star masses in the range $2 \le m \le 300 M_\odot$ drawn from a power-law initial mass function (IMF) $\phi(m) \propto m^{-1.3}$ (where $\phi(m)dm$ is the number of stars created with masses between $m$ and $m+dm$). For each pair, we draw the ratio of initial secondary-star mass to initial primary-star mass ($q$), the orbital period ($P$), and the eccentricity ($e$) from the corresponding probability density functions (PDFs) adopted in Sana et al. (2012): PDF($q) \propto q^{-0.1}$ with $q \in [0.1, 1.0]$; PDF($P) \propto P^{-0.55}$ with $P = \log(P/\text{day}) \in [0.15, 5.5]$; and PDF($e) \propto e^{-0.42}$ with $e \in [0, 0.9]$.³ The adopted primary-star IMF produces a larger proportion of massive stars than the standard Galactic-disc IMF, which has $\phi(m) \propto m^{-2.3}$ for $m \ge 1 M_\odot$ (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 $T_{\text{eff}}$), and gravity (log $g$) of the considered star. The spectra of stars hotter than the hottest Chen et al. (2015) model (i.e. with $T_{\text{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 high-resolution version of the Potsdam-WR (PoWR) models to describe the spectra of these stars (Hamann & Gräfener 2004; Sander, Hamann & Todt 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_\odot$ and hotter than 20 000 K, we adopt CMFGEN spectral models of stripped-envelope stars by Göteborg 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, Mapelli & Bressan 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 $S_v$ (Charlot et al., in preparation).

## 2.2 Photoionization modelling

The transmission function $T_v(t, t')$ 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, Charlot & Elbaz 2008; Vidal-García et al. 2017). Here, we are primarily interested in young star-forming galaxies with age close to the typical time-scale for dissipation of giant molecular clouds (~10 Myr.; e.g. Murray, Quataert & Thompson 2010; Murray 2011). In this case, we can write $T_v(t, t') \equiv T_v^+(t')$, where $T_v^+(t')$ is the transmission function of the ionized gas (Gutkin, Charlot & Bruzual 2016).

We compute $T_v^+(t')$ 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, Charlot & Bruzual 2016, see also section 2.1 of Plat et al. 2019):

(i) The hydrogen gas density, $n_H$.

(ii) The gas-phase metallicity, which unless otherwise specified is taken to be equal to that of the ionizing stars, $Z_{\text{ISM}} = Z$.

(iii) The carbon-to-oxygen abundance ratio, C/O.

(iv) The dust-to-metal mass ratio reflecting the depletion of heavy elements on to dust grains, $\xi_d$.

(v) The zero-age, volume-averaged ionization parameter, $\langle U \rangle \equiv \langle U \rangle(t' = 0)$. In spherical geometry, $\langle U \rangle(t')$ can be expressed in terms of the rate of ionizing photons produced by the evolving stellar population, $\dot{N}_H(t')$, the volume-filling factor of the gas (i.e. the ratio of the volume-averaged hydrogen density to $n_H$), $\epsilon$, and the case-B hydrogen recombination coefficient, $\alpha_B$, as (e.g. Panuzzo et al. 2003):

$$ \langle U \rangle(t') = \frac{3\alpha_B^{2/3}}{4c} \left[ \frac{3\dot{N}_H(t')\epsilon^2 n_H}{4\pi} \right]^{1/3} . \quad (2) $$

We stop the photoionization calculations at the outer edge of the H II region, when the electron density falls below 1 per cent of $n_H$ (we adopt a default inner radius of 0.1 pc). We adopt the prescription of Gutkin, Charlot & Bruzual (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 $Z_\odot = 0.01524$,

<footnote>$^2$While we experimented with samples of up to $3 \times 10^6$ binary pairs, populations of $10^6$ pairs were found sufficient to provide stable evolution of the nebular properties of young galaxies.</footnote>

<footnote>$^3$Since the tables of PARSEC tracks implemented in SEVN for this work do not include the evolution of stars less massive than $2 M_\odot$, we redraw a new pair if a secondary star is found below this limit.</footnote>
<figure><img src="image_4.png" /><figcaption>**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, $f_{\text{bin}}$ (coloured as indicated), for a metallicity $Z = 0.001$. The models assume a Chabrier (2003) IMF for primary stars and a Sana et al. (2012) distribution of secondary-to-primary mass ratios at age zero (single stars being modelled as non-interacting binaries) and are normalized to a total initial stellar mass of $1\,\text{M}_\odot$ integrated over 0.1–300 $\text{M}_\odot$ (see the text for details).</figcaption></figure>

corresponding to a proto-solar metallicity $Z_\odot^0 = 0.01774$ (Bressan et al. 2012).

It is of interest to examine how $\dot{N}_H(t')$ 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 $\psi(t') = \delta(t')$ in equation 1), for the typical metallicity $Z = 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, $f_{\text{bin}} = f_{\text{bin}}(t' = 0)$, defined as $f_{\text{bin}} = B/(S+B)$, where B is the number of binary pairs and S the number of single stars at $t' = 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 semimajor 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 $\phi(m) \propto m^{-2.35}$, at $m > 10\,\text{M}_\odot$. All models are normalized to a total initial mass in stars of $1\,\text{M}_\odot$ between 0.1 and $300\,\text{M}_\odot$.<fnref n="4" />

Fig. 1 shows that, at ages younger than 1 Myr, models with $f_{\text{bin}} > 0$ produce slightly larger $\dot{N}_H$, $\dot{N}_{\text{HeI}}$, and $\dot{N}_{\text{HeII}}$ than single-star models, due to the presence of merged stars with masses up to ~600 $\text{M}_\odot$ 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 $f_{\text{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 He-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 $\dot{N}_H$, $\dot{N}_{\text{HeI}}$, and $\dot{N}_{\text{HeII}}$. After a few tens of million years, the difference between models with $f_{\text{bin}} > 0$ and $f_{\text{bin}} = 0$ declines, especially for $\dot{N}_{\text{HeII}}$ 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 $\dot{N}_H$ and $\dot{N}_{\text{HeI}}$ relatively strong in binary-star models relative to single-star models. The rise of $\dot{N}_{\text{HeI}}$ and $\dot{N}_{\text{HeII}}$ at ages around 30 Myr in single-star 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 $\dot{N}_{\text{HeII}}$ 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, $\dot{N}_{\text{HeII}}/\dot{N}_H$, which is 2700 times larger in models with $f_{\text{bin}} = 1$ than in those with $f_{\text{bin}} = 0$, is already 54, 480, and 1300 times larger for $f_{\text{bin}} = 0.01, 0.1$, and 0.3, respectively. It is 2200 times larger for $f_{\text{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 $S_v(t')$ of SSPs with pure binary stars ($f_{\text{bin}} = 1$, solid lines) and pure single stars ($f_{\text{bin}} = 0$, dotted lines) at ages between 1 and 10 Myr, for a metallicity $Z = 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

2012 distribution) over the full range 0.1–300 $\text{M}_\odot$, even though only pairs with both primary- and secondary-star masses greater than 2 $\text{M}_\odot$ are retained (the pairs drawn with both primary- and secondary-star masses greater than 2 $\text{M}_\odot$ are matched to those in the base collection of Section 2.1). This has no influence on the predicted $\dot{N}_H$, $\dot{N}_{\text{HeI}}$, and $\dot{N}_{\text{HeII}}$, since main-sequence stars less massive than 2 $\text{M}_\odot$ are cooler than ~14 000 K at any metallicity and do not produce significant ionizing radiation. Also, rejecting binaries with primary mass above 2 $\text{M}_\odot$ 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 $\text{M}_\odot$ would increase $\dot{N}_H$, $\dot{N}_{\text{HeI}}$, and $\dot{N}_{\text{HeII}}$ by less than 1 per cent at ages below 40 Myr in Fig. 1.

<footnote n="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.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Spectral energy distribution $S_v$ entering equation (1) (in units of the luminosity per unit frequency at the Lyman limit) for the SSP models of Fig. 1 with pure binary stars ($f_{\text{bin}} = 1$, solid lines) and pure single stars ($f_{\text{bin}} = 0$, dotted lines) at the ages $t' \sim 3$ 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 $\alpha_{\text{AGN}} = -2.0$ (bottom edge) and $-1.2$ (top edge; Feltre, Charlot & Gutkin 2016). The ionization energies of some common species are indicated for reference.</figcaption></figure>

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 $t' \ge 3$ Myr in the binary-star model, the shape of $S_v(t')$ at ionizing energies up to $\gtrsim 50$ eV is similar to that of hard ionizing spectra of AGNs, $S_v \propto \nu^{\alpha_{\text{AGN}}}$ with $-2.0 \le \alpha_{\text{AGN}} \le -1.2$ (e.g. Feltre, Charlot & Gutkin 2016, grey shaded area).

# 3 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).

## 3.1 Observational sample

We appeal to the compilation of emission-line measurements of star-forming galaxies by Plat et al. (2019, based on the samples of Fosbury et al. 2003, Erb et al. 2010, Jaskot & Oey 2013, Stark et al. 2014, Berg et al. 2016, Nakajima et al. 2016, Vanzella et al. 2016, Izotov et al. 2016a, b, Amorín et al. 2017, Chisholm et al. 2017, Izotov, Thuan & Guseva 2017, Schmidt et al. 2017, Senchyna et al. 2017, Vanzella et al. 2017, Berg et al. 2018, Izotov et al. 2018a, b, Berg et al. 2019, Nanayakkara et al. 2019, and Senchyna et al. 2019). From this data set, we select all metal-poor galaxies with confident oxygen abundances (i.e. determined using the direct-$T_e$ or the strong-lines methods) less than $12 + \log(O/H) \approx 8.2$, corresponding to $Z \lesssim 0.004$ (see e.g. table 2 of Gutkin, Charlot & Bruzual 2016). With this selection, we retain 102 galaxies, of which 20 are classified

as Lyman-continuum (LyC) leakers, with fractions of escaping LyC photons up to ~0.6. The galaxies in this sample span wide ranges in redshift, $0 \lesssim z \lesssim 7$, ionization parameter, $-3.6 \lesssim \log(U) \lesssim -1.0$, C/O abundance ratio, $0.1 \lesssim (\text{C/O})/(\text{C/O})_{\odot} \lesssim 1.2$, and specific star formation rate, $1 \lesssim \psi/M \lesssim 1,000 \text{ Gyr}^{-1}$ (where M denotes the galaxy mass). For reference, we also consider the emission-line properties of 84 AGNs and candidate AGNs from Diaz, Prieto & Wamsteker (1988), Kraemer et al. (1994), Dors et al. (2014), Stark et al. (2015), and Nakajima et al. (2018). We refer the reader to section 3 and tables 1 and 2 of Plat et al. (2019) for more details about this sample.

We complement this data set with recent emission-line measurements for a few more local analogues to primeval galaxies, including the extremely metal-poor ($Z < 0.1 \text{ Z}_\odot$) starburst-dwarf galaxies SBS 0335-052E (Wofford et al. 2021), J1631+4426, J104457, and I Zw 18 NW (Umeda et al. 2022), seven extremely metal-poor galaxy candidates from Kojima et al. (2021), and six C IV emitters from Senchyna et al. (2022).

In this paper, we compare the predictions of our model with the observed intensity ratios and equivalent widths (EWs) of selected ultraviolet (C IV $\lambda$1549, He II $\lambda$1640, O III] $\lambda$1664, and C III] $\lambda$1908) and optical ([O III] $\lambda$5007, [O II] $\lambda$3727, [O I] $\lambda$6300, He II $\lambda$4686, H $\beta$, H $\alpha$, and [N II] $\lambda$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.

## 3.2 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 $Z = 0.001$, typical of the observational sample. For the photoionized gas, we adopt the same
<table><caption>**Table 1.** Nebular parameters of the 'standard' models of Section 3.2.</caption><thead><tr><th>Parameter</th><th>Value</th><th>Description</th></tr></thead><tbody><tr><td>Z<sub>ISM</sub></td><td>0.001</td><td>Interstellar metallicity</td></tr><tr><td>n<sub>H</sub></td><td>100 cm<sup>-3</sup></td><td>Hydrogen density</td></tr><tr><td>C/O</td><td>0.17</td><td>Carbon-to-oxygen abundance ratio</td></tr><tr><td>log (U)</td><td>-2.0</td><td>Volume-averaged ionization parameter at age t' = 0</td></tr><tr><td>ξ<sub>d</sub></td><td>0.3</td><td>Dust-to-metal mass ratio</td></tr></tbody></table>

metallicity ($Z_{ISM} = 0.001$), density ($n_H = 100 \text{ cm}^{-3}$), ionization parameter ($\log(U) = -2.0$), and dust-to-metal mass ratio ($\xi_d = 0.3$) typical of metal-poor star-forming galaxies (e.g. Berg et al. 2016; Gutkin, Charlot & Bruzual 2016; Mingozzi et al. 2022). These choices of $Z_{ISM}$ and $\xi_d$ imply a gas-phase oxygen abundance $12 + \log(O/H) \approx 7.53$ (table 2 of Gutkin, Charlot & Bruzual 2016), which according to fig. 6 of Berg et al. (2016) corresponds to a typical carbon-to-oxygen abundance ratio $C/O \approx 0.17$, i.e. about 40 per cent of the solar ratio ($C/O)_\odot = 0.44$ (Gutkin, Charlot & Bruzual 2016). This is close to the typical value $C/O \approx 0.18$ reported by Izotov et al. (2023) in low-redshift LyC-leaking galaxies with $12 + \log(O/H) \lesssim 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 time-scale 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 EWs of emission lines predicted by the SSP models of Fig. 1 with pure binary stars ($f_{bin} = 1$, labelled 'GALSEVN') and pure single stars ($f_{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 $\lambda4686/H\beta$-versus-EW($H\beta$) 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, Fragos & Izotov 2019). All models start at high EW($H\beta$) on the right-hand side of the diagram. As massive stars leave the main sequence and the first red supergiant stars appear in the single-star GALSEVN model, $\dot{N}_H$ and $\dot{N}_{He II}$ decrease while the continuum at $H\beta$ increases, causing both EW($H\beta$) and He II $\lambda4686/H\beta$ to drop without ever reaching the area of the diagram occupied by the observations.

Remarkably, in the binary-star GALSEVN model, the rise in $\dot{N}_{He II}$ triggered by the appearance of hot WNE-like stars with blue spectra around 3 Myr (Fig. 1) makes both He II $\lambda4686/H\beta$ and EW($H\beta$) 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 ($f_{bin} \gtrsim 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 metallicity $Z = 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 this study. We note the striking difference between the C&B and GALSEVN single-star models at ages around 2 Myr, when the brief ($\sim 0.2$ Myr) WR phase that makes

He II $\lambda4686/H\beta$ 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 $\dot{N}_H$ and $\dot{N}_{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, Izzard & Tout 2008, see also Fig. 4). After this phase and the disappearance of the most massive stars, EW($H\beta$) and He II $\lambda4686/H\beta$ 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 $\lambda4686/H\beta$ 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, $\dot{N}_{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, Phillips & Terlevich (1981, hereafter BTP) diagram to distinguish AGNs from star-forming galaxies (Fig. 3a), at early ages, all models lie at the high-[O III] $\lambda5007/H\beta$, low-[N II] $\lambda6584/H\alpha$ end characteristic of low metallicity and high ionization parameter (see e.g. fig. 2 of Gutkin, Charlot & Bruzual 2016). The GALSEVN and BPASS binary-star models stay in this area at ages up to 10 Myr, while the [O III] $\lambda5007/H\beta$ ratio of single-star models rapidly drops after the disappearance of photons more energetic than $\sim 35.1$ eV capable of ionizing $O^{2+}$. The implied drop in [O III] $\lambda5007$ luminosity is also the reason for the difference between binary- and single-star models in the [O III] $\lambda5007/[O II]\lambda3727$-versus-[O I] $\lambda6300/[O III] \lambda5007$ diagram (Fig. 3c). In Fig. 3(d), none of the models (plotted for clarity at slightly different abscissæ even though all have $12 + \log(O/H) = 7.53$) reaches the highest observed EW(He II $\lambda1640$) $\sim 2-8$ Å. In this diagram, only the GALSEVN binary-star models maintain EW(He II $\lambda1640$) $> 0.5$ Å at ages $t' \gtrsim 3$ Myr, while the BPASS binary-star model declines to EW(He II $\lambda1640$) $\approx 0.03$ Å and all single-star models to EW(He II $\lambda1640$) $\approx 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 $\lambda1549 + C III] \lambda1908$)/He II $\lambda1640$ and high C IV $\lambda1549/C III] \lambda1908$. 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 $t' \gtrsim 2$ Myr is closest to the observational sequence at He II $\lambda1640/O III] \lambda1664 > 0.3$, although it does not appear to reach sufficiently high EW(C III] $\lambda1908$) at a given He II $\lambda1640/O III] \lambda1664$.

On the whole, Fig. 3 shows the need to include binary stars and their inherent physical processes in modelling the emission-line properties 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.
<figure><img src="image_6.png" /><figcaption>**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 EWs 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) 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 abscissæ, although all have 12 + log(O/H) = 7.53. In all models, the ionizing stellar population is an SSP with metallicity Z = 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.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Rates of H- (top) and He II- (bottom) ionizing photons as a function of age for SSPs of metallicity $Z = 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_{\odot}$ integrated over 0.1–300 $M_{\odot}$ as in Fig. 1.</figcaption></figure>

## 3.3 Influence of adjustable model parameters

So far, we have considered only models with metallicity $Z = 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 $n_{\text{H}}$ and $\xi_{\text{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.

### 3.3.1 Metallicity

Fig. 5 shows the evolution of GALSEVN binary-star models with metallicities $Z = Z_{\text{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 $T_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 $[\text{N II}] \lambda 6584/\text{H \alpha}$ with $Z$ in Fig. 5(a) is further enhanced by the inclusion of secondary nitrogen enrichment in the nebular-emission models (Section 2.2). At $Z = 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 $\lambda 4686/\text{H \beta}$ and He II $\lambda 1640/\text{O III}] \lambda 1664$ in Figs 5(b) and (f).

As Fig. 5(d) shows, while metallicities in the range $0.0005 \le Z \le 0.004$ can fully describe the range of $12 + \log(\text{O/H})$ probed by the observational sample, metallicity has little influence on the maximum EW(He II $\lambda 1640$) achievable with the GALSEVN binary-star model.

### 3.3.2 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 \langle U \rangle = -1$ instead of $-2$; and one with C/O $= (\text{C/O})_{\odot} = 0.44$ instead of 0.17. All models have metallicity $Z = 0.001$.

The effect of increasing $\log \langle U \rangle$ is to increase the probability of multiple ionizations, and hence to strengthen high-ionization lines (see EWs of He II $\lambda 1640$ and C III $\lambda 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 $\alpha$ as $\log \langle U \rangle$ increases results from the conversion of $N^+$ into $N^{2+}$. In Fig. 6(b), the drop in EW(H $\beta$) is caused by the increase in H-column density, and hence in the absorption of H-ionizing photons by dust at fixed Z and $\xi_{\text{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 \langle U \rangle = -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, Charlot & Bruzual 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(\text{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. (2016, 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 supersolar C/O ratio at redshift $z \sim 12.5$.

### 3.3.3 O/Fe abundance ratio

Several recent studies have found evidence for supersolar ratios of $\alpha$ 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 supersolar 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 Fig. 3 (with $Z = Z_{\text{ISM}} = 0.001$) with two new GALSEVN binary-star models: one exploring the influence of supersolar O/Fe by setting the stellar metallicity to 10 per cent of the ISM metallicity (i.e. $Z = 0.0001$ versus $Z_{\text{ISM}} = 0.001$) and the other investigating the impact of a global reduction of metallicity, with $Z = Z_{\text{ISM}} = 0.0001$.

At fixed $Z_{\text{ISM}} = 0.001$, lowering the stellar metallicity from $Z = 0.001$ to 0.0001 implies bluer stellar-population spectra and harder
<figure><img src="image_8.png" /><figcaption>**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 $Z = 0.0005, 0.001, 0.002$, and $0.004$, respectively.</figcaption></figure>

ionizing radiation, since stars evolve at higher effective temperature (e.g. Bressan et al. 2012; Plat et al. 2019). The H β EW increases significantly in Fig. 7(b) (due to both higher $N_{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 $Z = Z_{ISM} = 0.0001$ in Fig. 7 produces a luminosity ratio and EWs of H β and He II λ1640 comparable to those

in the model with the same $Z = 0.0001$ but higher $Z_{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 $Z_{ISM}$ is reduced from 0.0001 to 0.001 (Figs 7a, e, and f).

On the whole, Fig. 7 shows that adopting supersolar O/Fe significantly affects the predicted emission-line properties of metal-poor star-forming galaxies, although we note that the impact on the EW of He II λ1640 is relatively modest in Fig. 7(d).
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Same as Fig. 3, for three GALSEVN binary-star models with metallicity $Z = 0.001$: a model with the standard parameters of Table 1 ($\log\langle U\rangle = -2$, C/O = 0.17; filled squares), one with enhanced ionization parameter ($\log\langle U\rangle = -1$, C/O = 0.17; filled triangles), and one with enhanced C/O abundance ratio ($\log\langle U\rangle = -2$, C/O = (C/O)$_{\odot} = 0.44$; nearly overlapping with the standard model in (b) and (c); filled crosses), as indicated on the right.</figcaption></figure>

# 4 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 EWs of lines with high ionization energies, such as He II $\lambda$1640 and C III $\lambda$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 $Z = 0.001$ and the standard nebular parameters listed in Table 1.
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Same as Fig. 3, for three GALSEVN binary-star models: a model with the standard parameters of Table 1 (with $Z = Z_{ISM} = 0.001$; filled squares), one exploring the influence of supersolar O/Fe (with $Z = 0.0001$ and $Z_{ISM} = 0.001$; filled triangles), and one investigating the impact of a global reduction of metallicity (with $Z = Z_{ISM} = 0.0001$; filled crosses). See the main text for details. In (d), the model with $Z_{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.</figcaption></figure>

## 4.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, Hill & Tosi 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.
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Same as Fig. 3, for three GALSEVN binary-star models with metallicity $Z = 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 the text for details).</figcaption></figure>

We can explore more realistic histories of star formation for metal-poor 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-
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Same as Fig. 3, but for 800 GALSEVN models computed using the star formation histories of 100 SPHINX galaxies from Rosdahl et al. (2022), combined alternatively with four metallicities ($Z = Z_{ISM} = 0.0005, 0.001, 0.002, \text{and } 0.004$) and two C/O abundance ratios ($C/O = 0.17 \text{ 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 two-dimensional histograms of models are colour-coded according to the density scale shown on the right.</figcaption></figure>

line properties of 800 models computed using the star formation histories or 100 SPHINX galaxies, combined alternatively with four metallicities ($Z = Z_{ISM} = 0.0005, 0.001, 0.002, \text{and } 0.004$) and two C/O abundance ratios ($C/O = 0.17 \text{ 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 EWs 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).

## 4.2 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 generalized Rosin-Rammler distribution function (e.g. Chabrier 2003; Wise et al. 2012; Goswami et al. 2022),

$$ \phi(m) \propto m^{-2.3} \exp \left[ - \left( \frac{m_{\text{char}}}{m} \right)^{1.6} \right], \quad (3) $$

which approaches the standard Chabrier (2003) IMF at large m but is exponentially cut off at $m < m_{\text{char}}$. Wise et al. (2012) used this function with $m_{\text{char}} = 40 M_{\odot}$ to explore the formation of Population III stars in early dwarf galaxies, while Goswami et al. (2022) adopted 200 M$_{\odot} \le m_{\text{char}} \le 300 M_{\odot}$ to explore chemical enrichment by PISNe in extremely metal-poor galaxies. Here, we consider values of $m_{\text{char}}$ in the range 50–100 M$_{\odot}$ to explore the emission-line signatures of stellar populations strongly dominated by massive stars.

The upper panel of Fig. 10 shows two versions of the top-heavy IMF of equation (3), for $m_{\text{char}} = 50 M_{\odot}$ and 100 M$_{\odot}$, compared to the standard Chabrier (2003) IMF. In the lower panel, we show the evolution of the absolute AB magnitude at 1500 Å, $M_{\text{UV}}$, of GALSEVN binary-star models computed with these three IMFs, assuming a total initial mass in stars of $10^6 M_{\odot}$, metallicity $Z = 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 $M_{\text{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 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 $M_{\text{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 supergiant stars, at ages around 2 Myr, the implied strong drop in ionizing-photon rate and rise in 4861-Å continuum emission make the H β EW decline more abruptly than in the

<figure><img src="image_13.png" /><figcaption>**Figure 10.** Top: top-heavy IMFs computed using equation (3) for two values of $m_{\text{char}}$ (red: 50 M$_{\odot}$; salmon: 100 M$_{\odot}$), compared to the standard Chabrier (2003) IMF (in black). All IMFs are normalized to a total mass in stars of 1 M$_{\odot}$. Bottom: absolute ultraviolet AB magnitude at 1500 Å, $M_{\text{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_{\odot}$, metallicity $Z = 0.001$ and the standard nebular parameters of Table 1. For reference, the dotted horizontal line marks the magnitude $M_{\text{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).</figcaption></figure>

standard model (Fig. 11b). Meanwhile, after reaching a minimum, the He II λ4686/H β ratio (Fig. 11b) and He II λ1640 EW (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 EW 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 supergiant 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 $Z = 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 EWs 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.

## 4.3 Emission from accretion discs of XRBs

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 XRBs, where a compact object (neutron star or stellar-mass black hole) accretes material from a companion (e.g. Schaerer, Fragos & Izotov
<figure><img src="image_14.png" /><figcaption>**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 $m_{\text{char}} = 50$ and $100 M_{\odot}$, respectively. The cyan stars mark the points beyond which the models with top-heavy IMFs become fainter than $M_{\text{UV}} = -14$ mag and are not shown (see the text and Fig. 10 for details).</figcaption></figure>

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
<figure><img src="image_15.png" /></figure>

**Figure 12.** X-ray luminosity function (in the 0.5–8 keV energy band) of XRBs for a GALSEVN binary-star population of metallicity $Z = 0.008$ with binary fraction $f_{\text{bin}} = 0.7$, assuming 100 Myr of constant star formation at the rate $\psi = 1\text{M}_{\odot}\text{yr}^{-1}$. Accretion is assumed to proceed according to the two-phase model described by equation (B5), with $f_{\text{A}} = 100$ and $\epsilon_{\text{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).

disc photons by energetic coronal electrons (Steiner et al. 2009; Mirocha 2014). This model, which also accounts for X-ray variability (e.g. Tanaka & Shibazaki 1996; Chen, Shrader & Livio 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 $\psi = 1\text{M}_{\odot}\text{yr}^{-1}$. The model has a binary-star fraction $f_{\text{bin}} = 0.7$ (typical of massive-star populations; e.g. Sana et al. 2012) and the metallicity $Z = 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\beta$ and HeII $\lambda$4686 luminosities predicted by such a model for a GALSEVN binary-star population of standard metallicity $Z = 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 per cent 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 HeII emission from metal-poor star-forming galaxies in the sample of Section 3.1.

<figure><img src="image_16.png" /></figure>

**Figure 13.** Contribution by XRB accretion discs (in blue), compared to the stellar contribution (in green), to the total luminosities (in black) of $H\beta$ (top panel) and HeII $\lambda$4686 (bottom panel) plotted against age, for a GALSEVN binary-star SSP with metallicity $Z = 0.001$. The model has the same zero-age Chabrier (2003) IMF normalized to a total initial stellar mass of $1\text{M}_{\odot}$ integrated over 0.1–300 $\text{M}_{\odot}$ as in Fig. 1.

## 4.4 Emission from radiative shocks

Radiative shocks, such as those expected from massive-star winds and SN blast waves, have spectral signatures characterized by strong high-ionization lines (e.g. Dopita & Sutherland 1996; Alarie & Morisset 2019), similar to those observed in some metal-poor star-forming galaxies (Izotov, Thuan & Privon 2012; Plat et al. 2019). As noted by Plat et al. (2019), the high gas densities ($n_{\text{H}} \gtrsim 10^4 \text{cm}^{-3}$) measured from the [C III]$\lambda$1907 + C III]$\lambda$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 HII 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 ($V_s$) and pre-shock density ($n_{\text{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, Mapelli & Bressan 2015; Spera & Mapelli 2017). For reference, the rates of energy injection computed in this
<figure><img src="image_17.png" /><figcaption>**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 Z = 0.001, plotted against age (see the text Section 4.4 for details). The contributions are shown for the range of shock velocities 100 km s⁻¹ ≤ V_s ≤ 1000 km s⁻¹ and the density n_H = 100 cm⁻³. 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.</figcaption></figure>

way are roughly consistent with those reported by Leitherer, Robert & Drissen (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 Z = 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 ≤ V_s ≤ 1000 km s⁻¹, and n_H = 100 cm⁻³ (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).

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⁵¹ erg for the most massive PISNe (footnote 5) would increase the magnitude of the first peak, but only for a brief spike at the earliest ages.

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⁵¹ erg for the most massive PISNe (footnote 5) 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 compact-dwarf galaxies (e.g. Izotov, Thuan & Privon 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 this 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.

## 5 DISCUSSION

In Sections 3 and 4, 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 β EWs 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 XRBs to the ionizing radiation of a star-forming galaxy (Section 4.3). Interestingly, we find that reproducing the observed luminosity function of XRBs 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. 2017, 2020; Saxena et al. 2020), this finding contrasts with some recent claims (Schaerer, Fragos & Izotov 2019; Umeda et al. 2022; Katz et al. 2023), which we now examine in more detail.

Schaerer, Fragos & Izotov (2019) used an observed anticorrelation between X-ray luminosity per unit star formation rate ($L_X^{\text{tot}}$/SFR) and metallicity (Douna et al. 2015; Brorby et al. 2016), combined with an estimate of the rate of He II- ionizing photons ($\dot{N}_{\text{HeII}}$) per unit $L_X^{\text{tot}}$ in the metal-poor star-forming galaxy IZw 18, to conclude that XRBs are likely to be the main source of He II emission in metal-poor star-forming galaxies. The founding assumption of this reasoning, that He II emission in IZw 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 II-
emission peak. Also, the prediction by Schaerer, Fragos & Izotov (2019) that the He II λ4686/H β line-luminosity ratio should correlate with $L_X^{tot}/\text{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 $L_X^{tot}/\text{SFR}$ in excess of $10^{42} \text{ erg s}^{-1}/\text{M}_\odot \text{ 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 $L_X^{tot}/\text{SFR}$ ratios of galaxies with and without He II emission at redshift $z \sim 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 $L_X^{tot}/\text{SFR}$ ratio by arbitrarily scaling an ARES multicolour-disc model, for a fixed black-hole mass of 25 M$_\odot$, added to the BPASS ionizing spectrum of a 5-Myr-old binary-star SSP.<fnref n="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 JWSST/NIRSpec spectrum of the metal-poor star-forming galaxy S04590 at redshift $z = 8.5$ requires $L_{X,2-10\text{keV}}^{tot}/\text{SFR} > 10^{41} \text{ erg s}^{-1}/\text{M}_\odot \text{ yr}^{-1}$, which translates into $L_X^{tot}/\text{SFR} > 3 \times 10^{41} \text{ erg s}^{-1}/\text{M}_\odot \text{ yr}^{-1}$ when adopting the energy band 0.5–8 keV used by Senchyna et al. (2020) and in Fig. B3 of this paper. Again, this is much larger than observed in extremely metal-poor star-forming galaxies (Senchyna et al. 2017, 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 $\langle L_X^{tot}/\text{SFR} \rangle \approx 2.1 \times 10^{40} \text{ erg s}^{-1}/\text{M}_\odot \text{ yr}^{-1}$ produced by XRBs in the five metal-poor 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 \times 10^{40} \text{ erg s}^{-1}/\text{M}_\odot \text{ yr}^{-1}$ implies that this ratio can vary from galaxy to galaxy. The 10 star-forming galaxies with 12 + log(O/H) < 8.0 in the sample of Douna et al. (2015, see Fig. B3) exhibit a very similar mean $\langle L_X^{tot}/\text{SFR} \rangle \approx 2.7 \times 10^{40} \text{ erg s}^{-1}/\text{M}_\odot \text{ yr}^{-1}$ with a standard deviation of $2.4 \times 10^{40} \text{ erg s}^{-1}/\text{M}_\odot \text{ 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 IZw 18 NW. Umeda et al. (2022) then show that, for each galaxy, the

<figure><img src="image_18.png" /><figcaption>**Figure 15.** Relative offset between line luminosities of the best-fitting (i.e. minimum-$\chi^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 IZw 18 NW (from top to bottom) studied by Umeda et al. (2022). All luminosities are normalized to that of H β. The $1\sigma$ (3 σ) observational errors are shown in red (gold). For J1631+4426, the [O III] λ4363 point falls outside the plot, at $(L_{\text{mod}} - L_{\text{obs}})/L_{\text{obs}} \approx -0.50$, as indicated by the down-pointing black triangle.</figcaption></figure>

best-fitting combination of blackbody + power-law ionizing spectra can be approached by a combination of the accretion-disc model spectrum of Gierliński, Done & Page (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 $L_X^{tot}/\text{SFR}$ ratios of these models are $5.0 \times 10^{40} \text{ erg s}^{-1}/\text{M}_\odot \text{ yr}^{-1}$ (J1631+4426), $1.0 \times 10^{41} \text{ erg s}^{-1}/\text{M}_\odot \text{ yr}^{-1}$ (J104457), and $3 \times 10^{41} \text{ erg s}^{-1}/\text{M}_\odot \text{ yr}^{-1}$ (I Zw 18 NW), reaching again up to $11 \sigma$ above the mean observed value in the Douna et al. (2015) sample of metal-poor star-forming galaxies.<fnref n="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 self-consistent scaling of accretion-disc with stellar emission, compare with those of J1631+4426, J104457, and IZw 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 $n_H = 100 \text{ cm}^{-3}$, $\log(U) = -2$, and $\xi_d = 0.3$, only two C/O abundance ratios (C/O = 0.17 and 0.44), and four metallicities ($Z = Z_{\text{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(U) = -1.7$ and $-1.4$, C/O = 0.09, and $\xi_d = 0.1$. We present the results of this exercise in Fig. 15, which shows the

<footnote n="5">⁵To compute $L_X^{tot}/\text{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).</footnote>

<footnote n="6">⁶To compute $L_X^{tot}/\text{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).</footnote>
relative offset between luminosity of the best-fitting (i.e. minimum-$\chi^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 $\lambda$4686 luminosities of these galaxies, corresponding to He II $\lambda$4686/H $\beta$ 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] $\lambda$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] $\lambda$4363/[O III] $\lambda$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] $\lambda$4363 emission, such as exceptional heating perhaps related to temperature inhomogeneities (e.g. Cameron, Katz & Rey 2023b), is beyond the scope of the simple exercise considered here.

For reference, the $L_X$/SFR ratios of the best-fitting models in Fig. 15, computed by estimating the SFR from the 1500-Å luminosity by analogy with Umeda et al. (2022, see footnote 6), span the range 3.6-7.2 × 10$^{40}$ erg s$^{-1}$/M$_{\odot}$ yr$^{-1}$, compatible with the observations of Lehmer et al. (2019) and Douna et al. (2015). The best-fitting metallicities and volume-averaged ionization parameters are (0.0005, -2.0), (0.001, -1.4), and (0.0005, -1.7) for J1631+4426, J104457, and I Zw 18 NW, respectively, close to the values (0.0004, -1.6), (0.001, -2.0), and (0.0005, -1.7) derived by Umeda et al. (2022) using blackbody + power-law ionizing spectra. For reference, the best-fitting C/O and $\xi_d$ values are (0.09,0.1), (0.44,0.3), and (0.09,0.1), respectively.

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 star-forming galaxies, appears robust when confronted with observations of both the emission-line and X-ray properties of such galaxies.

# 6 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 XRBs, and SN-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.

(i) 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 star-star 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-, and He-

ionizing photons changing only little for stellar populations with 0.3 $\lesssim$ $f_{\text{bin}} \le 1.0$ (Fig. 1).

(ii) We find that the He II $\lambda$4686/H $\beta$ ratios and H $\beta$ EWs of GALSEVN SSP models with metallicity $Z \sim 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).

(iii) 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).

(iv) We find that the most extreme He II $\lambda$1640 EWs (and other high-ionization signatures) of galaxies in our sample can be reproduced by models with a top-heavy IMF with characteristic mass $\sim 50$ M$_{\odot}$, which could represent a transient phase in the evolution of a galaxy (Fig. 11).

(v) Reproducing the average luminosity function of XRBs 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. 2017, 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, Fragos & Izotov 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).

(vi) 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.

# ACKNOWLEDGEMENTS

We are grateful to M. Celeste Artale, Jeremy Blaizot, Andrea Lapi, Jean-Pierre Lasota, Jordan Mirocha, Christophe Morisset, and Adèle Plat for helpful discussions. AB acknowledges financial support from PRIN-MIUR 2017. GB acknowledges financial support from the National Autonomous University of México (UNAM) through grants DGAPA/PAPIIT IG100319 and BG100622. MM, GC, and GI acknowledge financial support from the European Research Council (ERC) through the ERC Consolidator Grant DEMOBLACK, under contract No. 770017. MM acknowledges financial support from the German Excellence Strategy via the Heidelberg Cluster of Excellence (EXC 2181–390900948) ‘STRUCTURES’. GC acknowledges sup-
port from the Agence Nationale de la Recherche grant POPSYCLE No. ANR-19-CE31-0022. YC acknowledges support from the National Natural Science Foundation of China (NSFC) No. 12003001. JC acknowledges financial support from the ERC through the ERC Advanced Grant ‘FirstGalaxies’, under contract No. 789056. MD acknowledges financial support from the Cariparo Foundation under grant 55440.

DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the corresponding author. The code SEVN is publicly available via the following link: https://gitlab.com/sevncodes/sevn.git.

REFERENCES

Alarie A., Morisset C., 2019, *Rev. Mex. Astron. Astrofis.*, 55, 377

Allen M. G., Groves B. A., Dopita M. A., Sutherland R. S., Kewley L. J., 2008, *ApJS*, 178, 20

Amorín R. et al., 2017, Nat. Astron., 1, 0052

Arnett D., 1996, Princeton: Princeton University Press, Supernovae and Nucleosynthesis: An Investigation of the History of Matter from the Big Bang to the Present

Bachetti M. et al., 2014, *Nature*, 514, 202

Baldwin J. A., Phillips M. M., Terlevich R., 1981, *PASP*, 93, 5 (BTP)

Begelman M. C., 2002, *ApJ*, 568, L97

Begelman M. C., King A. R., Pringle J. E., 2006, *MNRAS*, 370, 399

Berg D. A., Skillman E. D., Henry R. B. C., Erb D. K., Carigi L., 2016, *ApJ*, 827, 126

Berg D. A., Erb D. K., Auger M. W., Pettini M., Brammer G. B., 2018, *ApJ*, 859, 164

Berg D. A., Erb D. K., Henry R. B. C., Skillman E. D., McQuinn K. B. W., 2019, *ApJ*, 874, 93

Berg D. A. et al., 2022, *ApJS*, 261, 31

Bressan A., Marigo P., Girardi L., Salasnich B., Dal Cero C., Rubele S., Nanni A., 2012, *MNRAS*, 427, 127

Brorby M., Kaaret P., Prestwich A., 2014, *MNRAS*, 441, 2346

Brorby M., Kaaret P., Prestwich A., Mirabel I. F., 2016, *MNRAS*, 457, 4081

Bruzual G., Charlot S., 2003, *MNRAS*, 344, 1000

Cameron A. G. W., Mock M., 1967, *Nature*, 215, 464

Cameron A. J. et al., 2023a, A&A, 677, A115

Cameron A. J., Katz H., Rey M. P., 2023b, *MNRAS*, 522, L89

Carnall A. C., 2017, arXiv e-print, (p. arXiv:1705.05165)

Chabrier G., 2003, *PASP*, 115, 763

Charlot S., Longhetti M., 2001, *MNRAS*, 323, 887

Chen W., Shrader C. R., Livio M., 1997, *ApJ*, 491, 312

Chen Y., Bressan A., Girardi L., Marigo P., Kong X., Lanza A., 2015, *MNRAS*, 452, 1068

Chevallard J., Charlot S., 2016, *MNRAS*, 462, 1415

Chisholm J., Orlitová I., Schaerer D., Verhamme A., Worseck G., Izotov Y. I., Thuan T. X., Guseva N. G., 2017, *A&A*, 605, A67

Costa G., Girardi L., Bressan A., Marigo P., Rodrigues T. S., Chen Y., Lanza A., Goudfrooij P., 2019, *MNRAS*, 485, 4641

Costa G., Bressan A., Mapelli M., Marigo P., Iorio G., Spera M., 2021, *MNRAS*, 501, 4514

Costa G., Mapelli M., Iorio G., Santoliquido F., Escobar G. J., Klessen R. S., Bressan A., 2023, *MNRAS*, 525, 2891

Cullen F. et al., 2021, *MNRAS*, 505, 903

da Cunha E., Charlot S., Elbaz D., 2008, *MNRAS*, 388, 1595

D'Eugenio F. et al., 2023, preprint (arXiv:2311.09908)

Diaz A. I., Prieto M. A., Wamsteker W., 1988, A&A, 195, 53

Dopita M. A., Sutherland R. S., 1996, *ApJS*, 102, 161

Dors O. L., Cardaci M. V., Hägele G. F., Krabbe Â. C., 2014, *MNRAS*, 443, 1291

Douna V. M., Pellizza L. J., Mirabel I. F., Pedrosa S. E., 2015, *A&A*, 579, A44

Draine B. T., McKee C. F., 1993, *ARA&A*, 31, 373

Dray L. M., 2006, *MNRAS*, 370, 2079

Eggleton P. P., 1971, *MNRAS*, 151, 351

Eldridge J. J., Stanway E. R., 2012, *MNRAS*, 419, 479

Eldridge J. J., Izzard R. G., Tout C. A., 2008, *MNRAS*, 384, 1109

Eldridge J. J., Langer N., Tout C. A., 2011, *MNRAS*, 414, 3501

Eldridge J. J., Stanway E. R., Xiao L., McClelland L. A. S., Taylor G., Ng M., Greis S. M. L., Bray J. C., 2017, *Publ. Astron. Soc. Aust.*, 34, e058

Erb D. K., Pettini M., Shapley A. E., Steidel C. C., Law D. R., Reddy N. A., 2010, *ApJ*, 719, 1168

Faulkner J., Lin D. N. C., Papalozou J., 1983, *MNRAS*, 205, 359

Feltre A., Charlot S., Gutkin J., 2016, *MNRAS*, 456, 3354

Ferland G. J. et al., 2017, *Rev. Mex. Astron. Astrofis.*, 53, 385

Fosbury R. A. E. et al., 2003, *ApJ*, 596, 797

Fragos T. et al., 2013, *ApJ*, 764, 41

Fryer C. L., Belczynski K., Wiktorowicz G., Dominik M., Kalogera V., Holz D. E., 2012, *ApJ*, 749, 91

Fürst F. et al., 2010, *A&A*, 519, A37

García M., Herrero A., Najarro F., Lennon D. J., Alejandro Urbaneja M., 2014, *ApJ*, 788, 64

Giacobbo N., Mapelli M., 2019, *MNRAS*, 482, 2234

Giacobbo N., Mapelli M., 2020, *ApJ*, 891, 141

Gierliński M., Done C., Page K., 2009, *MNRAS*, 392, 1106

Goswami S. et al., 2022, *A&A*, 663, A1

Göteborg Y., de Mink S. E., Groh J. H., Kupfer T., Crowther P. A., Zapartas E., Renzo M., 2018, *A&A*, 615, A78

Göteborg Y., de Mink S. E., McQuinn M., Zapartas E., Groh J. H., Norman C., 2020, *A&A*, 634, A134

Gräfener G., Owocki S. P., Grassitelli L., Langer N., 2017, *A&A*, 608, A34

Gutkin J., Charlot S., Bruzual G., 2016, *MNRAS*, 462, 1757

Hainich R., Ramachandran V., Shenar T., Sander A. A. C., Todt H., Gruner D., Oskinova L. M., Hamann W. R., 2019, *A&A*, 621, A85

Hamann W. R., Gräfener G., 2004, *A&A*, 427, 697

Hameury J. M., 2020, *Adv. Space Res.*, 66, 1004

Hameury J. M., Lasota J. P., 2020, *A&A*, 643, A171

Heckman T. M., Thompson T. A., 2017, in *Alsabti A. W., Murdin P. eds, Cham: Springer International Publishing, Handbook of Supernovae*. p. 2431

Heger A., Woosley S. E., 2002, *ApJ*, 567, 532

Hertz P., Vaughan B., Wood K. S., Norris J. P., Mitsuda K., Michelson P. F., Dotani T., 1992, *ApJ*, 396, 201

Hobbs G., Lorimer D. R., Lyne A. G., Kramer M., 2005, *MNRAS*, 360, 974

Hurley J. R., Tout C. A., Pols O. R., 2002, *MNRAS*, 329, 897

Hynes R. I. et al., 2003, *MNRAS*, 345, 292

Iorio G. et al., 2023, *MNRAS*, 524, 426

Izotov Y. I., Thuan T. X., 2004, *ApJ*, 616, 768

Izotov Y. I., Thuan T. X., Privon G., 2012, *MNRAS*, 427, 1229

Izotov Y. I., Schaerer D., Thuan T. X., Worseck G., Guseva N. G., Orlitová I., Verhamme A., 2016a, *MNRAS*, 461, 3683

Izotov Y. I., Orlitová I., Schaerer D., Thuan T. X., Verhamme A., Guseva N. G., Worseck G., 2016b, *Nature*, 529, 178

Izotov Y. I., Thuan T. X., Guseva N. G., 2017, *MNRAS*, 471, 548

Izotov Y. I., Schaerer D., Worseck G., Guseva N. G., Thuan T. X., Verhamme A., Orlitová I., Fricke K. J., 2018a, *MNRAS*, 474, 4514

Izotov Y. I., Worseck G., Schaerer D., Guseva N. G., Thuan T. X., Fricke Verhamme A., Orlitová I., 2018b, *MNRAS*, 478, 4851

Izotov Y. I., Schaerer D., Worseck G., Berg D., Chisholm J., Ravindranath S., Thuan T. X., 2023, *MNRAS*, 522, 1228

Jaskot A. E., Oey M. S., 2013, *ApJ*, 766, 91

Kalogera V., Webbink R. F., 1998, *ApJ*, 493, 351

Katz H. et al., 2023, *MNRAS*, 518, 592

Kauffmann G. et al., 2003, *MNRAS*, 346, 1055

Kehrig C., Guerrero M. A., Vílchez J. M., Ramos-Larios G., 2021, *ApJ*, 908, L54

Kennicutt Robert C. J., 1998, *ARA&A*, 36, 189

Kewley L. J., Dopita M. A., Sutherland R. S., Heisler C. A., Trevena J., 2001, *ApJ*, 556, 121
Kojima T. et al., 2021, ApJ, 913, 22

Kraemer S. B., Wu C.-C., Crenshaw D. M., Harrington J. P., 1994, ApJ, 435, 171

Kroupa P., 2001, MNRAS, 322, 231

Lapi A., Raimundo S., Aversa R., Cai Z. Y., Negrello M., Celotti A., De Zotti G., Danese L., 2014, ApJ, 782, 69

Lasota J.-P., 2001, New Astron. Rev., 45, 449

Lehmer B. D. et al., 2019, ApJS, 243, 3

Leitherer C., Heckman T. M., 1995, ApJS, 96, 9

Leitherer C., Robert C., Drissen L., 1992, ApJ, 401, 596

Li L.-X., 2012, MNRAS, 424, 1461

López-Navas E., Degenaar N., Parikh A. S., Hernández Santisteban J. V., van den Eijnden J., 2020, MNRAS, 493, 940

Ma X., Kasen D., Hopkins P. F., Faucher-Giguère C.-A., Quataert E., Kereš D., Murray N., 2015, MNRAS, 453, 960

Mapelli M., Zampieri L., Ripamonti E., Bressan A., 2013, MNRAS, 429, 2298

Mapelli M., Spera M., Montanari E., Limongi M., Chieffi A., Giacobbo N., Bressan A., Bouffanais Y., 2020, ApJ, 888, 76

Martin R. G., Nixon C., Armitage P. J., Lubow S. H., Price D. J., 2014, ApJ, 790, L34

Meyer F., Meyer-Hofmeister E., 1981, A&A, 104, L10

Mingozzi M. et al., 2022, ApJ, 939, 110

Mirocha J., 2014, MNRAS, 443, 1211

Mitsuda K. et al., 1984, PASJ, 36, 741

Moe M., Di Stefano R., 2017, ApJS, 230, 15

Murray N., 2011, ApJ, 729, 133

Murray N., Quataert E., Thompson T. A., 2010, ApJ, 709, 191

Nakajima K., Ellis R. S., Iwata I., Inoue A. K., Kusakabe H., Ouchi M., Robertson B. E., 2016, ApJ, 831, L9

Nakajima K. et al., 2018, A&A, 612, A94

Nanayakkara T. et al., 2019, A&A, 624, A89

Nguyen C. T. et al., 2022, A&A, 665, A126

Olivier G. M., Berg D. A., Chisholm J., Erb D. K., Pogge R. W., Skillman E. D., 2022, ApJ, 938, 16

Panuzzo P., Bressan A., Granato G. L., Silva L., Danese L., 2003, A&A, 409, 99

Pinto C., Middleton M. J., Fabian A. C., 2016, Nature, 533, 64

Plat A., Charlot S., Bruzual G., Feltre A., Vidal-García A., Morisset C., Chevallard J., Todt H., 2019, MNRAS, 490, 978

Pottschmidt K. et al., 2003, A&A, 407, 1039

Rees M. J., 1988, Nature, 333, 523

Remillard R. A., McClintock J. E., 2006, ARA&A, 44, 49

Rodríguez Castillo G. A. et al., 2020, ApJ, 895, 60

Rósrdahl J. et al., 2022, MNRAS, 515, 2386

Rosen A. L., 2022, ApJ, 941, 202

Salpeter E. E., 1955, ApJ, 121, 161

Sana H. et al., 2012, Science, 337, 444

Sánchez S. F. et al., 2022, ApJS, 262, 36

Sander A., Hamann W. R., Todt H., 2012, A&A, 540, A144

Santoliquido F., Mapelli M., Iorio G., Costa G., Glover S. C. O., Hartwig T., Klessen R. S., Merli L., 2023, MNRAS, 524, 307

Saxena A. et al., 2020, MNRAS, 496, 3796

Schaerer D., Fragos T., Izotov Y. I., 2019, A&A, 622, L10

Schmidt K. B. et al., 2017, ApJ, 839, 17

Senchyna P. et al., 2017, MNRAS, 472, 2608

Senchyna P., Stark D. P., Chevallard J., Charlot S., Jones T., Vidal-García A., 2019, MNRAS, 488, 3492

Senchyna P., Stark D. P., Mirocha J., Reines A. E., Charlot S., Jones T., Mulchaey J. S., 2020, MNRAS, 494, 941

Senchyna P. et al., 2022, ApJ, 930, 105

Shakura N. I., Sunyaev R. A., 1973, A&A, 24, 337

Shirazi M., Brinchmann J., 2012, MNRAS, 421, 1043

Shull J. M., McKee C. F., 1979, ApJ, 227, 131

Sonbas E., Dhuga K. S., Göğüş E., 2019, ApJ, 873, L12

Spera M., Mapelli M., 2017, MNRAS, 470, 4739

Spera M., Mapelli M., Bressan A., 2015, MNRAS, 451, 4086

Spera M., Mapelli M., Giacobbo N., Trani A. A., Bressan A., Costa G., 2019, MNRAS, 485, 889

Stanway E. R., Eldridge J. J., 2018, MNRAS, 479, 75

Stanway E. R., Eldridge J. J., 2019, A&A, 621, A105

Stark D. P. et al., 2014, MNRAS, 445, 3200

Stark D. P. et al., 2015, MNRAS, 454, 1393

Steidel C. C., Strom A. L., Pettini M., Rudie G. C., Reddy N. A., Trainor R. F., 2016, ApJS, 826, 159

Steiner J. F., Narayan R., McClintock J. E., Ebisawa K., 2009, PASP, 121, 1279

Strickland D. K., Heckman T. M., 2009, ApJ, 697, 2030

Strom A. L., Steidel C. C., Rudie G. C., Trainor R. F., Pettini M., 2018, ApJ, 868, 117

Sutton A. D., Roberts T. P., Middleton M. J., 2013, MNRAS, 435, 1758

Tanaka Y., Shibazaki N., 1996, ARA&A, 34, 607

Tauris T. M. et al., 2017, ApJ, 846, 170

Todt H., Sander A., Hainich R., Hamann W. R., Quade M., Shenar T., 2015, A&A, 579, A75

Tolstoy E., Hill V., Tosi M., 2009, ARA&A, 47, 371

Topping M. W., Shapley A. E., Reddy N. A., Sanders R. L., Coil A. L., Kriek M., Mobasher B., Siana B., 2020, MNRAS, 499, 1652

Umeda H., Ouchi M., Nakajima K., Isobe Y., Aoyama S., Harikane Y., Ono Y., Matsumoto A., 2022, ApJ, 930, 37

van den Eijnden J., Degenaar N., Russell T. D., Miller-Jones J. C. A., Rouco Escorial A., Wijnands R., Sivakoff G. R., Hernández Santisteban J. V., 2022, MNRAS, 516, 4844

van Paradijs J., 1996, ApJ, 464, L139

Vanzella E. et al., 2016, ApJ, 821, L27

Vanzella E. et al., 2017, ApJ, 842, 47

Vidal-García A., Charlot S., Bruzual G., Hubeny I., 2017, MNRAS, 470, 3532

Watarai K.-y., 2006, ApJ, 648, 523

Wise J. H., Abel T., Turk M. J., Norman M. L., Smith B. D., 2012, MNRAS, 427, 311

Wofford A., Vidal-García A., Feltre A., Chevallard J., Charlot S., Stark D. P., Herenz E. C., Hayes M., 2021, MNRAS, 500, 2908

Yan Z., Yu W., 2015, ApJ, 805, 87

Yang P., Zhang G., Russell D. M., Gelfand J. D., Méndez M., Wang J., Lyu M., 2022, MNRAS, 514, 234

## 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 PARSEC 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.

(i) 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, Tout & Pols 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.

(ii) 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.

(iii) 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 XRBs 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, $\dot{M}_{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, Tout & Pols 2002). For accretion of fully ionized H + He material, with hydrogen mass fraction X, on to a compact object of mass $M_c$, this limit can be expressed as

$$L_{Edd} \approx 1.26 \times 10^{38} \left(\frac{2}{1+X}\right) \left(\frac{M_c}{M_{\odot}}\right) \text{erg s}^{-1}. \quad (B1)$$

The Eddington luminosity is related to an Eddington accretion rate, $\dot{M}_{Edd}$, through the formula (e.g. Shakura & Sunyaev 1973)

$$L_{Edd} = \eta \dot{M}_{Edd} c^2, \quad (B2)$$

where $\eta$ is the radiative efficiency (i.e. the fraction of binding energy radiated away by the accretion flow) and c 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, Middleton & Fabian 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 $\dot{M}_{acc} > \dot{M}_{Edd}$, the total accretion luminosity, noted $L_{acc}$, can exceed the Eddington luminosity by a logarithmic factor $\sim \ln(\dot{M}_{acc}/\dot{M}_{Edd})$ (see Shakura & Sunyaev 1973; Begelman, King & Pringle 2006).

We adopt here the convenient approximation (Watarai 2006, see also Lapi et al. 2014)

$$L_{acc} \approx 2 \ln \left( 1 + \frac{\dot{M}_{acc}}{2 \dot{M}_{Edd}} \right) L_{Edd}, \quad (B3)$$

which, for sub-Eddington accretion ($\dot{M}_{acc} \lesssim \dot{M}_{Edd}$), reduces to

$$L_{acc} \approx \eta \dot{M}_{acc} c^2. \quad (B4)$$

For super-Eddington accretion, the weaker-than-linear increase of $L_{acc}$ with $\dot{M}_{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 $\eta \approx 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 $\dot{M}_{acc}^{max} = f_{Edd} \dot{M}_{Edd}$, where $f_{Edd}$ is an adjustable factor (Spera et al. 2019; Iorio et al. 2023).<fnref n="7" /> This 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 $\dot{M}_c/(1-\eta)$, where $\dot{M}_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 $\dot{M}_{acc} > \dot{M}_{acc}^{max}$) is assumed to be lost from the system (other options are described in Iorio et al. 2023).

## B2 Spectral modelling

Observations of XRBs 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, Shrader & Livio 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_{\odot}$ (e.g. Yan & Yu 2015), it is also observed in high-mass XRBs (e.g. Martin et al. 2014; van den Eijnden 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, Dhuga & Göğüş 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, Lin & Papaloizou 1983; van Paradijs 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 angular-momentum 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.

<footnote n="7">⁷The quantity $f_{Edd}$ is noted $\eta_{Edd}$ in Iorio et al. (2023).</footnote>
<figure><img src="image_19.png" /><figcaption>**Figure B1.** X-ray luminosity functions of XRBs for GALSEVN populations with binary fraction $f_{\text{bin}} = 0.7$, assuming 100 Myr of constant star formation at the rate $\psi = 1\,\text{M}_\odot\,\text{yr}^{-1}$. Accretion is taken to proceed according to the two-phase model of equation (B5), for different choices of the fraction $\epsilon_A$ of mass accreted, while the rate is boosted by a factor $f_A$, as indicated. All models have the metallicity $Z = 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 the text for details).</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Figure B2.** Same as Fig. B1, for GALSEVN binary-star populations of different metallicities in the range $0.0001 \le Z \le 0.02$, as indicated. In all models, accretion is taken to proceed according to the two-phase model described by equation (B5), with $f_A = 100$ and $\epsilon_A = 0.99$.</figcaption></figure>

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, $\dot{M}_c/(1-\eta)$, and an outburst phase during which it is boosted by a factor $f_A$. We assume that a fraction $\epsilon_A$ of the mass is accreted during the outburst phase. For the accreted mass $\delta M_c$ to be conserved over the time interval $\delta t'$, we thus write

$$ \dot{M}_{\text{acc}}(t') = \begin{cases} \frac{1}{1-\eta} \frac{\delta M_c}{\delta t'} & \text{during } (1-\epsilon_A)\delta t', \\ \frac{f_A}{1-\eta} \frac{\delta M_c}{\delta t'} & \text{during } \frac{\epsilon_A}{f_A} \delta t', \end{cases} \quad (B5) $$

where $\delta M_c$ is the mass growth of the compact object in a time interval $\delta t'$ of the evolution of an SSP (adopting the same notation as in equation 1 of Section 2). For reference, $\delta M_c$ represents at most a few per cent of $M_c$.

We describe the spectral distribution of the power $L_{\text{acc}}(t')$ produced by the hot accretion disc at time $t'$ (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

<figure><img src="image_21.png" /><figcaption>**Figure B3.** Total X-ray luminosity $L_X^{\text{tot}}$ (per unit star formation rate) obtained by integrating the luminosity functions of Fig. B2, for metallicities in the range $0.0005 \le Z \le 0.02$ (red crosses). The data are from Chandra observations of a sample of 10 nearby star-forming galaxies with $12 + \log(\text{O/H}) < 8.0$ compiled by Douna et al. (2015, green points; see also Brorby, Kaaret & Prestwich 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(\text{O/H}) \approx 8.4$ shows the mean $L_X^{\text{tot}}/\text{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. The dotted horizontal green and blue lines represent the mean $L_X^{\text{tot}}/\text{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(\text{O/H})$ at 0.14.</figcaption></figure>

the disc. To this end, we use the ARES<fnref n="8" /> 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, $M_{\text{ARES}}$, and the maximum radius of the accretion disc, $R_{\text{max}}$. Following Senchyna et al. (2020), we fix $R_{\text{max}}$ to $10^4$ gravitational radii (we find that decreasing or increasing $R_{\text{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 $M_c = M_{\text{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⁻¹, Sutton, Roberts & Middleton (2013) conclude that the spectra of sources with sub-Eddington accretion rates may be well represented by multicolour discs, while those of sources with super-Eddington accretion rates tend to exhibit ‘hard-luminous’ and ‘soft-luminous’ components (in roughly equal proportions among the sample), presumably depending on inclination. For each $M_{\text{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 $L_{\text{acc}}(t') \le L_{\text{Edd}}$, we adopt the pure multicolour-disc spectrum. Otherwise, we randomly draw between the pure multicolour-disc and the Comptonized spectra.<fnref n="9" /> We note that the results presented in this paper about the X-ray

<footnote n="8">⁸https://ares.readthedocs.io/en/latest/</footnote>

<footnote n="9">⁹Since the ARES spectral library pertains to objects accreting at the Eddington rate, in practice, for a given $L_{\text{acc}}$, we adopt the spectrum corresponding to the ARES model for a compact object accreting at the Eddington luminosity $L_{\text{Edd}} = L_{\text{acc}}$ given by equation (B1) for a pure-H gas, i.e. $M_{\text{ARES}} = L_{\text{acc}}/(1.26 \times 10^{38} \text{erg s}^{-1}) M_{\odot}$.</footnote>
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.

The above approach allows us to compute the spectrum of every hot accretion disc of XRB in a GALSEVN binary-star SSP of age t', from which we extract the disc luminosity in the 0.5–8 keV energy band, noted $L_X$. For reference, the quantity $L_X$ amounts to typically between 30 and 90 per cent of $L_{acc}$. Fig. B1 shows X-ray luminosity functions (XLFs) of XRBs obtained in this way for populations with binary-star fraction $f_{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 $f_A$ and $\epsilon_A$ in equation (B5). All models have the metallicity $Z = 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 $L_X = 10^{37}$ erg s⁻¹ to about 3 per cent at $L_X = 10^{40}$ erg s⁻¹ (for details, see Lehmer et al. 2019). We therefore ignore this component in our model. Furthermore, we adopt here $f_{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 $\dot{M}_{acc} > 100 \dot{M}_{Edd}$.

We find that the predicted XLF depends in a somewhat degenerate way on the parameters $f_A$ and $\epsilon_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 $\epsilon_A = 1.00, 0.99$ and $0.95$ of the mass at $f_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/50th 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 $f_A f_{Edd} = 10^4–10^5$ times the Eddington rate $\dot{M}_{Edd}$, corresponding to a maximum accretion luminosity $L_{acc} \approx 17–22 L_{Edd}$ (equation B3), compatible with the $10–100 L_{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 $f_A = 100$).

The model with a single phase of enhanced accretion is the simplest conceptually. The one with $\epsilon_A = 0.99$ and $f_A = 100$ provides formally the best fit (minimum $\chi^2$) to the data, while the one with quiescent 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 $\epsilon_A = 0.99$ and $f_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 $f_A = 100$ and $\epsilon_A = 0.99$, for GALSEVN binary-star populations of different metallicities in the range $0.0001 \le Z \le 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 low-mass companions: nearly 80 per cent of overall accretion occurs in systems composed of a compact object lighter than $20 M_\odot$ (with progenitor main-sequence mass typically smaller than $30 M_\odot$; e.g. Spera, Mapelli & Bressan 2015) and a companion lighter than $25 M_\odot$. Since the evolution of stars lighter than $30 M_\odot$ 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, $L_X^{tot}/SFR$, obtained by integrating the XLFs of the binary-star models with different metallicities (converted to gas-phase O abundances using table 2 of Gutkin, Charlot & Bruzual 2016 for $\xi_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, Kaaret & Prestwich 2014) and the sample of nearby analogues of distant Lyman-break galaxies compiled by Brorby et al. (2016). Also shown is the mean $L_X^{tot}/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.<fnref n="10" /> As expected from Fig. B1, the model with $Z = 0.008$ (corresponding to $12 + \log(O/H) \approx 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, metal-poor 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 XRBs consistently with the spectral modelling of stars.

## APPENDIX C: LINE EMISSION FROM SHOCKS DRIVEN BY STELLAR WINDS AND SN EXPLOSIONS IN STAR-FORMING GALAXIES

The energy deposited into the ISM by stellar winds and SN explosions (e.g. Leitherer, Robert & Drissen 1992; Leitherer & Heckman 1995) can drive fast shocks contributing to line emission in star-forming galaxies (e.g. Shull & McKee 1979; Izotov, Thuan & Privon 2012). We express the rate of total energy deposition into the ISM at SSP age $t'$ (equation 1) as

$$ \dot{E}_{\text{tot}}(t') = \epsilon_W \dot{E}_W(t') + \epsilon_{\text{SN}} r_{\text{SN}}(t') E_{\text{SN}}, \quad (C1) $$

<footnote n="10">¹⁰For each galaxy, $L_X^{tot}/SFR$ was obtained by integrating the observed individual luminosity function in fig. 3 of Lehmer et al. (2019) and then dividing by the star formation rate in their table 1.</footnote>
where $\dot{E}_W$ is the rate of energy input by stellar winds, $E_{SN} \approx 10^{51}$ erg is the typical energy released by an SN explosion (e.g. Arnett 1996),$^{11}$ $r_{SN}(t')$ is the rate of SN explosions at time $t'$, and $\epsilon_W$ and $\epsilon_{SN}$ are thermalization efficiency factors. This energy injection will generate shocks, whose radiative properties we wish to estimate.

Emission-line predictions for fast radiative shocks have been computed using the MAPPINGS photoionization code (e.g. Dopita & Sutherland 1996; Allen et al. 2008; Alarie & Morisset 2019). These usually pertain to high-velocity, steady-flow, fully radiative, plane-parallel magneto-hydrodynamic shocks and are traditionally expressed in flux per unit area (erg cm$^{-2}$ s$^{-1}$). We need to scale these predictions to the total energy deposition rate from equation (C1) to obtain estimates of emission-line luminosities (in erg s$^{-1}$). This can be achieved by considering the total radiative flux of a shock. Assuming that the entire flux of mechanical energy through the shock is radiated, Dopita & Sutherland (1996) express this quantity as (see their equation 3.3)

$$F_{\text{tot}}^{\text{shock}} = 2.28 \times 10^{-3} \left( \frac{V_s}{100 \text{ km s}^{-1}} \right)^{3.0} \left( \frac{n_H}{\text{cm}^{-3}} \right) \text{erg cm}^{-2} \text{s}^{-1}, \quad (C2)$$

where $V_s$ denotes the shock velocity and $n_H$ the hydrogen density. The fraction of this flux radiated in the H $\beta$ recombination line in the plasma column behind the shock can be written (equation 3.4 of Dopita & Sutherland 1996)

$$F_{\text{H}\beta}^{\text{shock}} = 7.44 \times 10^{-6} \left( \frac{V_s}{100 \text{ km s}^{-1}} \right)^{2.41} \left( \frac{n_H}{\text{cm}^{-3}} \right) \text{erg cm}^{-2} \text{s}^{-1}. \quad (C3)$$

The photoionized region ahead of the shock also contributes to the total H $\beta$ emission. The corresponding flux can be expressed as (equation 4.4 of Dopita & Sutherland 1996)

$$F_{\text{H}\beta}^{\text{prec}} = 9.85 \times 10^{-6} \left( \frac{V_s}{100 \text{ km s}^{-1}} \right)^{2.28} \left( \frac{n_H}{\text{cm}^{-3}} \right) \text{erg cm}^{-2} \text{s}^{-1}. \quad (C4)$$

This is an upper limit to the contribution by the shock precursor, as the ISM can run out before the end of the H-Strömgren column ahead of the shock (especially for the fastest shocks). Equations (C2)–(C4) allow one to express the ratio of total (shock + precursor) H $\beta$ flux to total radiative flux generated by the shock in terms of the quantities $V_s$ and $n_H$. This ratio can also be regarded as the ratio of total H $\beta$ luminosity to total luminosity emitted by the shock, i.e. we write

$$\frac{F_{\text{H}\beta}^{\text{shock}} + F_{\text{H}\beta}^{\text{prec}}}{F_{\text{tot}}^{\text{shock}}} = \frac{L_{\text{H}\beta}^{\text{shock}} + L_{\text{H}\beta}^{\text{prec}}}{L_{\text{tot}}^{\text{shock}}}. \quad (C5)$$

Since we have assumed that the entire flux of mechanical energy through the shock is radiated, we can replace $L_{\text{tot}}^{\text{shock}}$ by the power injected into the ISM (equation C1) and derive from this the total H $\beta$ luminosity of the shock.

Unfortunately, analogues to equations (C2)–(C4) (derived by Dopita & Sutherland 1996 using MAPPINGS II) are not available for the more recent calculations of Alarie & Morisset (2019, based on MAPPINGS V). Instead, Alarie & Morisset (2019) provide tables of $F_{\text{H}\beta}^{\text{shock}}$ and $F_{\text{H}\beta}^{\text{prec}}$ versus shock velocity $V_s$ and density $n_H$, while

$$L_{\text{H}\beta}^{\text{shock}} + L_{\text{H}\beta}^{\text{prec}} = 4.39 \times 10^2 \epsilon_W \dot{E}_W \left( \frac{F_{\text{H}\beta}^{\text{shock}} + F_{\text{H}\beta}^{\text{prec}}}{\text{erg cm}^{-2} \text{s}^{-1}} \right) \times \left( \frac{V_s}{100 \text{ km s}^{-1}} \right)^{-3.0} \left( \frac{n_H}{\text{cm}^{-3}} \right)^{-1} \text{erg s}^{-1}, \quad (C6)$$

$$L_{\text{H}\beta}^{\text{shock}} + L_{\text{H}\beta}^{\text{prec}} = 1.39 \times 10^{46} \epsilon_{\text{SN}} \left( \frac{r_{\text{SN}}}{\text{yr}^{-1}} \right) \left( \frac{F_{\text{H}\beta}^{\text{shock}} + F_{\text{H}\beta}^{\text{prec}}}{\text{erg cm}^{-2} \text{s}^{-1}} \right) \times \left( \frac{V_s}{100 \text{ km s}^{-1}} \right)^{-3.0} \left( \frac{n_H}{\text{cm}^{-3}} \right)^{-1} \text{erg s}^{-1}. \quad (C7)$$

$$L_X^{\text{shock}} + L_X^{\text{prec}} = L_{\text{H}\beta}^{\text{shock}} \frac{F_X^{\text{shock}}}{F_{\text{H}\beta}^{\text{shock}}} + L_{\text{H}\beta}^{\text{prec}} \frac{F_X^{\text{prec}}}{F_{\text{H}\beta}^{\text{prec}}}, \quad (C8)$$

<footnote>$^{11}$This is generally the case for Type-II, Type-Ia, and low-mass pair-instability SNe, but the most massive PISNe can reach ~100 times this value (see</footnote>

<footnote>Heger & Woosley 2002). This has little effect on our conclusions (see Section 4.4).</footnote>

<footnote>equation (C2), which states that $F_{\text{tot}}^{\text{shock}}$ is equal to the flux of mechanical energy passing through the shock, can be assumed not to depend on the photoionization code. We therefore combine the results of Alarie & Morisset (2019) with the expression of $F_{\text{tot}}$ as a function of $V_s$ and $n_H$ from equation (C2) to compute the ratio on the left-hand side of equation (C5). Then, we consider the luminosity $L_{\text{H}\beta}^{\text{shock}} + L_{\text{H}\beta}^{\text{prec}}$ separately for radiative shocks powered by stellar winds and SN explosions, as the ejection velocities associated with SN explosions ($\gtrsim 10, 000 \text{ km s}^{-1}$; e.g. Draine & McKee 1993) are higher than those associated with stellar winds ($\lesssim 3000 \text{ km s}^{-1}$; e.g. Garcia et al. 2014; Gräfener et al. 2017). The implied radiative-shock velocities are therefore likely to differ on average.</footnote>

<footnote>We compute the total H $\beta$ emission from radiative shocks driven by stellar winds as (dropping the dependence on $t'$ for simplicity)</footnote>

<footnote>where the total flux $F_{\text{H}\beta}^{\text{shock}} + F_{\text{H}\beta}^{\text{prec}}$ for shocks with velocity $V_s$ propagating in a medium with density $n_H$ is taken from Alarie & Morisset (2019). We obtain the rate $\dot{E}_W$ of energy input from winds (equation C1) by summing the individual contributions $\frac{1}{2}\dot{m}\dot{v}_\infty^2$ of all stars, where $\dot{m}$ is the stellar mass-loss rate and $v_\infty$ the terminal wind velocity.$^{12}$ For OB stars, we adopt the widely used relation $v_\infty = 2.65 v_{\text{esc}}$, where $v_{\text{esc}}$ is the stellar escape velocity (e.g. Garcia et al. 2014). For WR stars, we take $v_\infty$ to be the minimum of 1.30 $v_{\text{esc}}$ and 1800 km s$^{-1}$ (e.g. Gräfener et al. 2017). For simplicity, we assume here that all of the kinetic energy injected by winds into the ISM is thermalized, i.e. we take $\epsilon_W = 1.0$ (e.g. Rosen 2022).</footnote>

<footnote>Similarly, we express the total H $\beta$ emission from radiative shocks driven by SN explosions as</footnote>

<footnote>Constraints from the prototypical starburst galaxy M82 suggest $0.3 \lesssim \epsilon_{\text{SN}} \lesssim 1.0$ (Strickland & Heckman 2009; Heckman & Thompson 2017). Here, we adopt for simplicity $\epsilon_{\text{SN}} = 1.0$, which should provide an upper limit on the shock contribution to line emission.</footnote>

<footnote>Once the H $\beta$ luminosity is known, whether for wind- or SN-driven shocks, the luminosity of any other line X can be computed as</footnote>

<footnote>where the ratios $F_X^{\text{shock}}/F_{\text{H}\beta}^{\text{shock}}$ and $F_X^{\text{prec}}/F_{\text{H}\beta}^{\text{prec}}$ are available from the tables of Alarie & Morisset (2019).</footnote>

<footnote>This paper has been typeset from a TeX/LA TeX file prepared by the author.</footnote>

<footnote>$^{12}$For simplicity, we compute this estimate for a population of single stars and ignore the influence of binary-star processes on energy injection by stellar winds into the ISM.</footnote>