Modeling the High-energy Ionizing Output from Simple Stellar and X-Ray Binary Populations

We present a methodology for modeling the joint ionizing impact due to a “simple X-ray population” (SXP) and its corresponding simple stellar population (SSP), where “simple” refers to a single age and metallicity population. We construct composite spectral energy distributions (SEDs) including contributions from ultraluminous X-ray sources and stars, with physically meaningful and consistent consideration of the relative contributions of each component as a function of instantaneous burst age and stellar metallicity. These composite SEDs are used as input for photoionization modeling with Cloudy, from which we produce a grid for the time- and metallicity-dependent nebular emission from these composite populations. We make the results from the photoionization simulations publicly available. We find that the addition of the SXP prolongs the high-energy ionizing output from the population—and correspondingly increases the intensity of nebular lines such as He ii λ1640,4686, [Ne v] λ3426,14.3 μm, and [O iv] 25.9 μm by factors of at least two relative to models without an SXP spectral component. This effect is most pronounced for instantaneous bursts of star formation on timescales >10 Myr and at low metallicities (∼0.1 Z ⊙), due to the imposed time- and metallicity-dependent behavior of the SXP relative to the SSP. We propose nebular emission line diagnostics accessible with JWST suitable for inferring the presence of a composite SXP + SSP, and we discuss how the ionization signatures compare to models for sources such as intermediate-mass black holes.


INTRODUCTION
Over the past decade, spectroscopic observations of high redshift (z > 6) galaxies and their nearby analogs have revealed the presence of strong high-ionization nebular emission lines (e.g., He II, [Ar IV], [Ne V]) and lines with high equivalent widths, indicating the presence of relatively hard ionizing spectra, recent ( 30 Myr) star formation, and relatively low metallicities (∼ 0.1 Z ) (e.g., Jaskot & Oey 2013;Stark 2016;Berg et al. 2018Berg et al. , 2019;;Senchyna et al. 2017Senchyna et al. , 2019;;Olivier et al. 2022).Attempts to model the observed nebular features of such "extreme emission line galaxies" (EELGs) have revealed that the spectral energy distributions (SEDs) from young, low-metallicity stellar populations have difficulty reproducing the observed strengths of a number of high-ionization emission line species, a problem that has been referred to as the "high-energy ionizing photon production problem" (Berg et al. 2021).
Addressing this problem is key to understanding how the heating and reionization of the Universe proceeds at z > 6, a regime now becoming more accessible spectroscopically thanks to JWST.In the coming decade, second generation interferometers such as the Hydrogen Epoch of Reionization Array (HERA) will probe even earlier cosmic epochs (6 < z < 50) through the cosmic 21-cm signal.This high-redshift 21-cm signal is sensitive to the UV to soft X-ray (< 2 keV) radiation field produced by the first galaxies' dominant ionizing populations.Initial results from HERA already suggest that the emergent 0.5-2.0keV X-ray luminosity per galaxy star formation rate (L X (0.5-2 keV)/SFR) from galaxies at z > 6 is at least an order of magnitude higher than measured for local galaxies (HERA Collaboration et al. 2023).These results point to increased production efficiency of photons with energies 500 eV from galaxies at high redshift, possibly due to increased formation efficiency of accreting compact objects at the lower metallicities and/or younger stellar populations ages characteristic of early galaxies (e.g., Fragos et al. 2013a).Developing a comprehensive accounting of the major sources of high-energy ionizing photons in EELGs, both near and far, therefore requires a framework for modeling how additional ionizing sources-such as fast shocks, superbubbles, and X-ray binaries (XRBs)-evolve as a function of redshift-dependent properties such as metallicity and stellar population age (e.g., Allen et al. 2008;Oskinova & Schaerer 2022;Simmonds et al. 2021).
One potential alternative source of high-energy ionizing photons is radiative shocks, which may be produced as a consequence of massive star evolution through winds and supernovae explosions (e.g., Izotov et al. 2012Izotov et al. , 2021)).For fast shocks, there are publicly available grids of emission line ratios for a range of shock properties (Allen et al. 2008).Based on such models, shocks appear capable of sufficient high-energy ionizing photon production to reproduce select observed nebular line intensities for some shock velocities (e.g., Thuan & Izotov 2005); however, it is not yet clear whether fast shocks are capable of simultaneously reproducing suites of highionization nebular lines (i.e., from UV to optical), particularly if there is a metallicity-dependence to the observed line ratios (e.g., Senchyna et al. 2017).
Numerous works have also investigated the need for changes to prescriptions used in stellar spectral population synthesis and photoionization modeling, such as changes to assumed abundance patterns (e.g., Berg et al. 2021), the initial mass function (IMF) and associated products of binary evolution (e.g., Götberg et al. 2019;Senchyna et al. 2021), and/or assumptions about stellar rotation and line-driven winds for massive stars at very low metallicities (e.g., Telford et al. 2021).There are publicly available and well-vetted libraries for stellar atmosphere models and isochrones that allow for testing how these different prescriptions affect the ionizing output from stellar populations (see e.g., Conroy 2013, for a review of models).With stellar spectral popu-lation synthesis, these ingredients can be flexibly combined for a given IMF to produce the SEDs for simple stellar populations (SSPs)-populations at a single burst age and metallicity-for use in photoionization modeling and spectrophotometric fitting (e.g., Conroy et al. 2009;Conroy & Gunn 2010;Byler et al. 2017;Johnson et al. 2021).These tools enable investigation of how tweaks to stellar population models can address the high-energy ionizing photon production problem (e.g., Berg et al. 2019;Senchyna et al. 2021).However, presently available stellar population models, including those that model some of the products of binary evolution, still have difficulty reproducing the observed intensities of select high-ionization nebular emission lines (e.g., He II; Stanway & Eldridge 2019).This highlights the need for more comprehensive accounting of the products of binary evolution that may contribute to the ionizing photon budget.
Accreting stellar mass compact objects (i.e., XRBs) are one such product of binary evolution, and have long been considered as a potential source of high-energy ionizing photons (e.g., Garnett et al. 1991;Schaerer et al. 2019;Kovlakas et al. 2022).Recent photoionization modeling that combines XRBs with stellar populations hints that XRBs may be sufficient to reproduce observed high-ionization emission line strengths in some cases (e.g., Senchyna et al. 2020;Simmonds et al. 2021).However, drawing a cohesive picture of the efficacy of XRBs to addressing the high-energy ionizing photon production problem remains challenging due to the following issues: (1) differing prescriptions for intrinsic XRB SEDs and; (2) minimal consideration towards how to scale XRB power output relative to the stellar population as a function of star formation history (SFH) and metallicity.Critically, for XRBs we currently lack many of the ingredients necessary for truly flexible spectral population synthesis, and therefore a means to understand the XRB contribution to high-energy ionizing photon production.In this work, we take steps to address this via the development of a flexible framework for including XRBs alongside SSPs for use in photoionization modeling.
XRBs, and specifically those sources at the extreme bright-end of the luminosity function known as ultraluminous X-ray sources (ULXs), are a compelling option for producing high-energy ionizing photons in high redshift galaxies and their local analogs for a few key reasons.Specifically, such sources (1) form swiftly ( 3 Myr) in multiple generations following a burst of star formation, and are therefore capable of producing ionizing photons on longer timescales than single massive stars (e.g., Linden et al. 2010); (2) dominate the X-ray power output from normal star-forming galaxies (i.e., those without an actively accreting supermassive black hole) at high specific SFRs (e.g., Lehmer et al. 2016); and (3) have integrated luminosities directly proportional to SFR that increase with decreasing gas-phase metallicity, resulting in increased contribution to the photon budget in low-metallicity, highly star-forming galaxies (e.g., Basu-Zych et al. 2013, 2016;Brorby & Kaaret 2017).For ULXs specifically, there is also increasing evidence for the presence of powerful outflows (e.g., Pinto et al. 2016), which may additionally produce shocks (López et al. 2019;Gúrpide et al. 2022).ULXs may therefore be important sources of both radiative and mechanical feedback at the low metallicities and high specific SFRs characteristic of high redshift galaxies and their local analogs.
A key consideration for modeling the ionizing contribution from accreting stellar mass compact objects such as ULXs is that they should form as a consequence of the evolution of the stellar population.To address this, we present a methodology for constructing a simple X-ray population (SXP), analogous to an SSP (i.e., single age and metallicity), and an approach to a physically consistent coupling of the SXP to the corresponding SSP for use in photoionization modeling.This approach is aimed at providing a more comprehensive framework for investigating the XRB contribution to the high-energy ionizing photon production problem.
This paper is organized as follows: in Section 2 we present the main elements used in constructing the SXP and in Section 3 we outline how the SXP is scaled and added to the corresponding SSP to create the SED for a composite population.In Section 4 we describe our gridbased photoionization modeling set-up using the composite SXP + SSP SEDs.In Section 5 we present key results from the photoionization modeling for the nebular emission from the composite populations, and describe how to access our simulation results.In Section 6 we discuss the SXP model and photoionization simulation results in the context of X-ray detectability, recent literature results for the contribution of X-ray sources to the production of high-ionization emission lines, and alternative ionizing sources.Finally, in Section 7, we summarize our key findings.

CONSTRUCTING A SIMPLE X-RAY POPULATION
To construct the inputs for the photoionization modeling that follows, we apply the formalism of SSPs to accreting stellar mass compact objects.In what follows, we refer to models for single age and metallicity accreting compact object "populations" as "SXPs."In constructing these SXPs we consider selection of the fol-lowing three key ingredients: (1) power output for the population as a function of stellar population age and metallicity; (2) the dominant source type in the population; and (3) physically motivated intrinsic source spectra.In this section, we describe the SXP used in this work based on selections for these three ingredients.

Normalization as a Function of Metallicity and Instantaneous Burst Age
Constructing a true SXP requires a prescription for the normalization of X-ray emitting population as a function of burst age and metallicity.To describe this behavior, we use the theoretical predictions for the X-ray power output per stellar mass from XRBs as a function of age and stellar metallicity from the binary population synthesis models presented in Fragos et al. (2013b).We opt to use these theoretical predictions, as opposed to empirical measurements, for two key reasons.First, the evolution of XRB power output-particularly on finely resolved timescales 100 Myr that are important for luminous XRB formation-is not yet well-constrained observationally (e.g., Lehmer et al. 2017).In addition, empirical constraints on XRB radiative power as a function of time or metallicity are typically reported in terms of L X /SFR, where the SFR is a quantity averaged over some extended timescale (i.e., 10 or 100 Myr).In this way, available empirical measurements are a timeaveraged L X for the observed XRB population.Such measurements are therefore not well-suited to simulating the ionizing output from the intrinsic XRB population formed from an instantaneous burst of star formation, as we aim to do here.
A critical component of the selected theoretical population synthesis models from Fragos et al. (2013b) is therefore that the predictions are provided in terms of radiative power (L X ) from the intrinsic XRB population per stellar mass formed in a burst of star formation (M 0 ).Nearly all other theoretical binary population synthesis models currently in the literature that focus on XRBs provide predictions in terms of number counts of sources per stellar mass as function of time and metallicity (e.g., Linden et al. 2010;Wiktorowicz et al. 2019).For the purposes of the photoionization modeling that follows, L X /M 0 is the preferred quantity since transforming number counts to total radiative power is nontrivial, and requires additional assumptions about how the luminosities of different sources should be calculated from theoretical mass supply rates.The models from Fragos et al. (2013b) track the mass supply rate to the compact object population-both black holes (BHs) and neutron stars (NS)-as a function of time.Following prescriptions in Fragos et al. (2008Fragos et al. ( , 2009)), the bolomet- .The adopted scaling of LX/M0 for the SXP as a function of instantaneous burst age and metallicity (Z : dark blue dashed line, 0.1 Z : purple solid line, and 0.01 Z : orange dash-dot line) from theoretical binary population synthesis models (Fragos et al. 2013b).The instantaneous burst ages at which we model the SXP contribution in the photoionization simulations are marked with dark blue squares, purple circles, and orange diamonds for the Z , 0.1 Z , and 0.01 Z cases, respectively.The SXP component in our models turns on for instantaneous bursts > 3 Myr.Empirical results for LX/M0 from star-forming galaxy stacks in the Chandra Deep Fields, for which the median metallicity (∼ 0.6 Z ) is intermediate to the theoretical scalings adopted here, are shown as a green dotted line for reference.
ric L X is calculated directly for persistent XRBs from the mass supply rate along with the mass and radius of the accretor, assuming different conversion efficiencies for BHs and NS.For transient XRBs, a modified version of this formalism is applied, considering typical duty cycles and outburst luminosities.
The Fragos et al. (2013b) binary population synthesis models used here provide published values of L X /M 0 for XRBs at three different stellar metallicities1 : Z = 0.002 (0.1 Z ), Z = 0.02 (Z ), and Z = 0.03 (1.5 Z ), and burst ages from 0-10 Gyr.To calculate the normalization of the SXP as a function of age and metallicity, we select a set of seven discrete burst ages and three metallicities from the reference model in Fragos et al. 2013b (i.e., their Figure 2, data provided via private communication).
The seven discrete time steps (t burst ∼ {1, 3, 5, 8, 10, 15, 20} Myr) are selected to maximize L X /M 0 (i.e., capture the bounding case for XRB ionization), and to cover relevant timescales for important stages of massive star evolution.On these short burst timescales, the XRB population power output is dominated by sources descended from the most massive stars, and such sources are therefore considered high-mass XRBs (HMXBs).As such, we calculate the normalization of the SXP to reproduce the theoretical L X (0.5-8 keV)/M 0 from the HMXB component of the reference model in Fragos et al. (2013b).We additionally elect to turn on the SXP with a delay time that accounts for the minimum time for the first BHs to form and begin accreting (e.g., Belczynski & Taam 2008;Linden et al. 2010).The first time step in the simulations with SXP contribution is 4 Myr.
Two of the three metallicities selected for our models include those available from the published theoretical population synthesis results, namely Z = 0.002 (0.1 Z ) and Z = 0.02 (Z ), where the 0.1 Z scaling for L X /M 0 is relevant for nearby EELGs and sources at z 3 (e.g., Shirazi & Brinchmann 2012;Madau & Dickinson 2014).We additionally add a third metallicity set at ∼ 0.01 Z , though this metallicity is not originally included in the theoretical population synthesis models.To create normalizations for the SXP for this extremely metal-poor case, we assume the same dependence of L X /M 0 on post-starburst timescale as the 0.1 Z model, but allow the absolute normalization of L X /M 0 to be a factor of 3× higher than for the 0.1 Z at each time step.This scaling corresponds roughly to the observed L X /SFR from galaxies with 12 + log(O/H) ∼ 7.0, the lowest metallicity for which the observed L X -SFR-Z scaling relation has been calibrated (Lehmer et al. 2021).Though speculative, we include this case for reference, noting that it will likely need to be revised either theoretically or observationally going forward.
In Figure 1, we show the adopted theoretical L X /M 0 as a function of instantaneous burst age for the three metallicities employed in the photoionization simulations that follow, with the burst ages for which the SXP contribution is modeled denoted with symbols.For reference, we also plot in this figure empirical results from Gilbertson et al. (2022) for L X /M 0 for a sample of star-forming galaxy stacks from Chandra Deep Fields.The galaxy stacks from the Chandra Deep Fields have median metallicity ∼ 0.6 Z , intermediate to our selected metallicity range, illustrating that the empirical results are in qualitative agreement with the behavior of the adopted theoretical relations on the timescales of interest for these simulations.The adopted theoretical values for L X /M 0 from Figure 1 can be considered as good approximations for the average radiative output from a population.For the purposes of our simulations, this is appropriate as we seek to simulate a well-sampled population; however, for observed samples the X-ray luminosity function (XLF) may not be as well-populated.Stochastic sampling of the XLF is particularly important to consider at very low SFRs ( 0.1 M yr −1 ), where observed values for L X /M or L X /SFR can be subject to scatter of order ∼ 1 dex, especially for sources at the bright end of the XLF (Lehmer et al. 2021).We return to this point in Section 6.1 in discussing how results from these simulations should be compared against observations.

Source Type: Ultra-luminous X-ray Sources
With the adopted scaling for L X /M 0 , we next consider the source type that dominates the power output from the population.In star-forming galaxies without a central accreting supermassive black hole, XRBs dominate the galaxy-integrated X-ray point source emission.For predominantly young stellar populations ( 100 Myr) HMXBs are the dominant accretor population, while in galaxies with older stellar populations low-mass XRBs (i.e., sources with low-mass donor stars) dominate the radiative output (e.g., Garofali et al. 2018;Lehmer et al. 2019Lehmer et al. , 2021)).
Regardless of stellar population age, the total radiative power output from a population will be dominated by the sources that populate the bright end of the XLF.Such sources are ULXs, off-nuclear point sources with L X > 10 39 erg s −1 .These luminosities, under the assumption of isotropic emission, exceed the Eddington limit for a 10 M BH or 1.4 M NS: where M is the mass of the compact object (such that m ≡ M/M ), X is the hydrogen mass fraction (assumed to be X = 0.73, corresponding to solar composition), G is the gravitational constant, c is the speed of light, and κ T = 0.2(1 + X) is the Thomson scattering opacity.Models to explain the extreme apparent luminosities from ULXs fall broadly into two categories: (1) isotropic emission from sub-Eddington accretion onto intermediate mass black holes (IMBHs; e.g., Colbert & Mushotzky 1999) or (2) super-Eddington mass transfer via modified accretion disks onto either stellar mass BHs or strongly magnetized NS (e.g., King et al. 2001;Begelman et al. 2006;Poutanen et al. 2007;King 2009).Recent results from X-ray timing studies favor the second scenario, as the detection of pulsations in an increasing number of ULXs indicates the presence of a NS for at least some fraction of the ULX population (e.g., Bachetti et al. 2014;Israel et al. 2017;Brightman et al. 2018).Broadband and high-resolution X-ray spectroscopy of ULXs also support the picture of stellar mass accretors with high mass transfer rates.In particular, the observed spectral turnover for ULXs at energies 3-8 keV is incompatible with models for intermediate mass or supermassive BHs accreting at sub-Eddington rates (Gladstone et al. 2009a), and the detection of absorption lines for some ULXs provide evidence for the presence of fast outflows (Section 2.3.2;Pinto et al. 2016;Kosec et al. 2018aKosec et al. ,b, 2021)), and potentially strongly magnetized NS accretors (Walton et al. 2018;Brightman et al. 2022).
These findings suggest ULXs (L X ∼ 10 39 − 10 42 erg s −1 , i.e., below the hyper-luminous regime) are a distinct class of stellar mass accretors with super-Eddington mass supply rates2 .This in turn implies that ULXs are indeed the high luminosity extension of the XRB luminosity function.In this work, we treat ULXs as the bright extension of the HMXB luminosity function specifically, as we consider timescales relevant to formation and evolution of the most massive stars (≤ 20 Myr, Section 2.1).
Given that they occupy the bright end of the XLF, ULXs in star-forming galaxies can surpass by at least an order of magnitude the combined output of various HMXB sub-populations filling out the lower luminosity-and more highly populated-portion of the luminosity function (e.g., Lehmer et al. 2021).In this work, we therefore choose to model the SXP as a population of ULXs.Assuming a well-sampled XLF (i.e., high SFR) as we model here, a ULX-based SXP is a reasonable approximation for the total XRB population.Indeed, many highly star-forming galaxies (SFRs ≥ 3 M yr −1 ) show ULX-like spectra globally (Garofali et al. 2020).Because various HMXB sub-populations likely have different accretion state distributions and scaling relations with metallicity and stellar population age, approximating the SXP using ULXs considerably simplifies the modeling that follows.Hereafter, we refer to the SXP as "SXP ULX ", to indicate the source type modeled.
This selection captures a bounding case for HMXB3 contribution to production of high-energy ionizing photons in terms absolute normalization as emphasized above, as well as hardness of spectral shape (Section 2.3).

Model for Ultra-luminous X-ray Source SED
The efficacy of the SXP ULX in the production of high-energy ionizing photons is highly dependent on the shape of the intrinsic SED and normalization relative to the corresponding SSP (e.g., Simmonds et al. 2021).Unfortunately, this intrinsic shape for ULXs is difficult to measure observationally and therefore remains highly uncertain.At the typical extragalactic distances ( 3 Mpc) to ULXs, fluxes are < 10 −12 erg s −1 cm −2 , such that high-resolution X-ray spectroscopy is currently feasible for a limited number (∼ 10) of the brightest sources (Kosec et al. 2018a).At CCD-resolution, many more sources are available for spectral fitting, but these fits are often performed using phenomenological models with no obvious physical basis, which cannot be reliably extrapolated into wavelength regimes outside the observed bandpass.The lack of comprehensive multiwavelength catalogs for ULXs (from IR-to-X-ray) is a further limitation to measuring their intrinsic SEDs, resulting in large uncertainties especially in the extreme ultraviolet (EUV) to very soft X-ray regime (∼ 54-200 eV).Critically, this is the wavelength regime responsible for setting the intensity of high-ionization nebular emission lines.
We therefore opt to model our SED (hereafter SED ULX ) using analytic prescriptions from theoretical models for stellar mass compact objects with super-Eddington mass supply rates.In the sections that follow, we outline the key components of the SED ULX , and motivate our selections.The SED ULX employed throughout this work is qualitatively similar to empirical models derived from analyses of ULXs detected within their own optical nebulae, where high-ionization emission lines are also detected (Berghea et al. 2010;Berghea & Dudik 2012;Kaaret & Corbel 2009;Lebouteiller et al. 2017;Simmonds et al. 2021).In such analyses, both the broad-band X-ray data and the observed emission line strengths provide constraints on the form of the ULX SED, where the line emission is particularly useful for constraining the shape of the unseen EUV portion.We note that such a multi-wavelength approach has promise for constraining the shape of the intrinsic SED; however, as pointed out in Simmonds et al. (2021), there is as yet no comprehensive approach in this respect for ULXs given the lack of collated multiwavelength data sets.

Accretor Type
As a starting point to constructing the SED ULX , we must select the mass of the accretor as well as the type (BH or NS), as the SED shape will depend not only on accretor mass and mass supply rate, but also on whether the accretion flow encounters a hard surface and/or a strong magnetic field (see Section 6.3.1).Because the selected L X /M 0 normalizations from Section 2.1 span burst ages ≤ 20 Myr, we choose a stellar mass BH (< 100 M ) as the accretor for the SED ULX .Theoretical binary population synthesis models show that the fraction of NS accretors contributing to the observed ULX population is still low on timescales < 20 Myr post-starburst, and that accreting BHs may have ages ∼ 4-40 Myr during their ULX phase (Middleton & King 2017;Wiktorowicz et al. 2017Wiktorowicz et al. , 2019)).
Given that we likewise elect to model the population at three discrete stellar metallicities (Z , 0.1Z , and ∼ 0.01Z ), we construct the SED ULX model as a function of approximate BH masses at each of these metallicities: 10 M , 20 M , and 40 M from solar to subsolar metallicities, respectively (e.g., Belczynski et al. 2010;Fryer et al. 2012;Eldridge & Stanway 2016).The selected masses are likely on the high-end of the mass function for BHs in ULXs, again in order to capture the bounding case for SXP ULX ionization.At present it is not clear theoretically or observationally what the typical BH mass is in a ULX, particularly as a function of burst age and metallicity.Thus, we consider using approximate BH masses as a function of metallicity only to be a reasonable simplifying choice given presently available results from binary population synthesis models for ULXs.It is important to note that while these selected BH masses factor into the parametrization of the SED ULX shape (Section 2.3.3),they do not strongly affect the normalization of the SED ULX , which is key to setting the simulated nebular line intensities.Instead, the primary factor affecting SED ULX normalization, and therefore nebular line intensities, is the adopted scaling relation for L X /M 0 as a function of age and metallicity, as presented in Section 2.1.

Accretion Flow Geometry & Components
Under the selection of a stellar mass accretor with super-Eddington mass supply rate for the SED ULX , theory and observation alike suggest the presence of a strong outflow, which is a key component in defining the SED shape.In the supercritical regime the accretion flow consists of distinct zones delineated by the wind/outflow component (e.g., Shakura & Sunyaev 1973;Lipunova 1999;Poutanen et al. 2007;Abolmasov et al. 2009).We illustrate these zones in Figure 2, as follows: • The hot inner region of the accretion flow before the wind launching zone (R < R ph, in ).
• The region between the wind launching zone and so-called the spherization radius (R ph, in < R < R sph ), where the wind is optically thick.The spherization radius (R sph ) can be thought of as the region where the disk is unstable and super-Eddington (i.e., the actual accretion rate is Eddington within this radius, and limited through the expulsion of matter from the disk as a wind).
• The region beyond the spherization radius where the wind becomes optically thin (R > R sph ).The outer extent of the wind is given by the effective photospheric radius of the outflow, R ph, out .
Under such a model, the launching of the disk wind leads to an evacuated "wind funnel," wherein the hard inner disk emission is geometrically collimated and scattered by the wind.At much larger disk radii (R > R sph ), the wind becomes optically thin and can be approximated as a pseudo-photosphere (i.e., a quasi-spherical component with a blackbody spectrum).This picture of the accretion flow is broadly consistent with CCD-resolution spectra of ULXs, which typically require multi-component models to fit their observed spectra (e.g., Gladstone et al. 2009b;Sutton et al. 2013;Middleton et al. 2015).
In defining our SED ULX model, we therefore include a component representing the hot inner disk and wind funnel, as well as a component for the outflow.In addition, we consider the effects of Compton scattering and irradiation on the resultant SED shape.Irradiation-where one component of the accretion flow absorbs and reemits radiation from another component-may be critically important in the context of SXP ULX contribution to high-energy ionizing photon production.This is because reprocessed emission due to irradiation can produce a spectrum with a more substantial EUV component, where photon energies exceed ionization potentials for lines such as He II λ4686 and [Ne V] λ3426.Given the geometry of the supercritical accretion flow described above, it is unlikely that there is direct irradiation of the outer regions of the accretion disk by hot inner disk photons (e.g., Gierliński et al. 2009).In the case of a ULX with an outflow, the wind likely blocks and reprocesses the hot inner disk emission.Irradiation may therefore occur via the wind irradiating the outer disk (Vinokurov et al. 2013;Yao & Feng 2019), or via "self-irradiation" of inner disk and funnel wall emission within the wind funnel itself (Abolmasov et al. 2009).
For the SED ULX in this work, we choose to model the irradiated component as being due to self-irradiation within the wind funnel.This is because there are uncertainties in how effectively soft photons from the wind thermalize in the outer disk (Kaaret & Corbel 2009;Grisé et al. 2012a), and additional complexities in handling the geometry and radiative transfer for wind irradiation.To represent the supercritical accretion flow with self-irradiation, we opt to use the sirf, or "selfirradiated funnel", template as implemented in XSPEC (Abolmasov et al. 2009).This particular template provides a physical model coupling all components of the supercritical accretion flow, including the hot inner disk, outflow, and the effects of irradiation consistent with the previously described flow geometry.To additionally account for the effects of Compton scattering, we include the simpl(SIMple Power Law) convolution component, which Comptonizes flux output from the sirf model.Thus, our SED ULX is implemented as simpl(sirf) in In this respect, the sirf template is ideal, as it allows for changing the inclination angle to simulate changes to the apparent emergent spectrum for an observer, a point we return to in more detail in Section 6.1.

Parametrization of SED Shape
We aim to construct the SED ULX model with as few selected parameter values as possible.Therefore, to the extent possible, we calculate parameter values for the sirf template analytically from a few key selected physical properties (see Shakura & Sunyaev 1973;Poutanen et al. 2007;Abolmasov et al. 2009, for prescriptions).The primary physical properties for which values are selected are the funnel opening angle (θ) and compact object mass (m).For the SED ULX used in this work, we select θ = 45 • and m = 10, 20, or 40 M for Z , 0.1 Z , and ∼ 0.01 Z , respectively.As shown by the black solid line in Figure 3, the choice of θ = 45 • results in a relatively flat SED in the EUV regime even without invoking high mass supply rates (i.e., this selection assumes the SXP ULX is not universally strongly beamed; Equations 5-6).This choice for θ is therefore important for capturing a bounding case for production of high-energy ionizing photons by the SXP ULX , while still producing good qualitative agreement with other ULX SED models across the electromagnetic spectrum, as discussed in Section 6.2.Below we outline how the remaining parameters defining the shape of the SED ULX are determined analytically.
The mass transfer rate is a parameter critical to defining the spectral shape.By definition, a source in an ultra-luminous accretion state will have a mass supply rate exceeding the Eddington rate: where m is the mass of the accretor in units of M .It is often convenient to define a dimensionless mass transfer rate ṁ, as follows: A source in the ultra-luminous state would therefore be expected to have ṁ 1.The intrinsic luminosity of the disk then depends logarithmically on the dimensionless mass transfer rate (Poutanen et al. 2007): Following King (2009), the observed, or apparent, luminosity may differ from the intrinsic luminosity due to geometrical beaming, such that L obs = L int bol /b.Under this formalism, the beaming factor b is related to the funnel opening angle (θ) as b = Ω/4π, where Ω = 2 × 2π(1 − cosθ) is the combined solid angle from both cones of the wind funnel.In this way, the beaming factor (b) also sets the probability that the source is observed: For large funnel opening angles (θ ∼ 90 • ), the beaming factor approaches unity, and there is a high probability the source is observed.For very small opening angles, by contrast, the beaming factor is small and the observed luminosity of the central disk is boosted, but the probability that the source is observed at X-ray wavelengths is low.King (2009) suggest that the beaming factor itself may depend on the mass transfer rate, and we adopt this view for our model to reduce the number of predeter-mined parameter values.Specifically, we parameterize b in terms of the mass transfer rate as follows: By this definition, high mass transfer rates ( ṁ > 200) lead to small funnel opening angles, correspondingly small beaming factors, and therefore low probabilities for observing hard emission from the source.The model for the SED ULX used in this work has θ = 45 • , and therefore b = P obs (θ) ∼ 0.3 and ṁ ∼ 16, using Equations 5-6.In Figure 3 we show the SED ULX in black for various viewing angles (i = 0 • for the intrinsic model) alongside SEDs under different assumptions for the funnel opening angle and mass transfer rate and the relation between these parameters; these are discussed in more detail in Section 6.1.With θ, ṁ, and m specified as above, it is relatively straightforward to calculate a number of the remaining critical parameters setting the shape of the sirf component of the SED ULX model (see Table 1).The spherization radius where the outflow is launched is calculated as R sph = 18 cosθ R S ṁ, where R S = 2 GM c 2 is the Schwarzschild radius4 .The minimum inner radius for the photosphere (r ph,in ) in units of this spherization radius (R sph ) can be calculated as as the innermost stable circular orbit (R isco = 6 GM c 2 ) divided by θ: The corresponding temperature at the inner edge of the wind funnel (T funnel, in ), and the outer radius of the photopshere (r ph,out ), also in units of the spherization radius, can be calculated analytically following Abolmasov et al. (2009): where ṁ3 = Ṁ in units of 10 3 ṀEdd , Ω = −8π ln sinθ, M 10 is the mass of the BH in units of 10 M , and α is the exponent describing the velocity law for the wind.

Wavelength [ Å]
Figure 3.The broad-band SED for the SXPULX (i.e., simpl(sirf) in XSPEC) under differing assumptions about the funnel opening angle (θ), the mass transfer rate onto the BH ( ṁ), and the viewing angle (i = 0 • for a face-on observer; i = 90 • for an edge-on observer).All SEDs correspond to a 20 M BH, and are normalized to the same intrinsic bolometric luminosity.The SEDULX (intrinsic model for the photoionization simulations) is shown as the solid black line, while SEDs with different viewing angles for the same θ and ṁ are shown as dashed and dotted gray lines.In the case of an edge-on observer (i = 90) the source emits mostly in the EUV.The solid green and orange lines show the effect of changing the funnel opening angle (θ), which changes the corresponding mass transfer rate, while the solid blue SED demonstrates the effect of uncoupling θ and ṁ and allowing for a very high mass transfer rate with a moderate funnel opening angle.et al. (2009), where we have assumed T funnel, in is equal to the temperature at the bottom of the funnel (T bot ) for the sirf template.
The remainder of the parameters for the sirf component of the SED ULX model are selected, as summarized in Table 1, where notation corresponds to symbols used in XSPEC.In particular, the SED ULX model employed in the photoionization simulations has i = 0 • (such that the cloud in the photoionization modeling described in Section 4 sees all components of the inflow/outflow), α = 0 (no acceleration for the wind), γ = 4/3 (default value, the adiabatic index for the inner parts of the accretion flow), and irrad.= 4 (number of iterations for irradiation).The normalization for the SED ULX model is calculated following the theoretical scaling of L X /M 0 with starburst age from Fragos et al. (2013b) as described in Section 2.1.
Finally, we convolve the sirf component of the SED ULX model with a Comptonization component (simpl) in order to reproduce the typical high-energy (> 2 keV) shape of observed ULX spectra.All parameters values of the simpl component are selected, namely the photon power-law index Γ = 2.5, and the scattering fraction (10%).For observed ULXs fit with the simpl convolution component, the recovered Comptonization parameters are broadly consistent with those listed in Table 1; however, such parameters are typically not tightly constrained from observations (e.g.Walton et al. 2013Walton et al. , 2014)).As such, we consider the SED ULX shape at energies 2 keV, which is set by the simpl component, critical for comparisons to apparent observed L X , but not particularly physically informative in terms of shape.
We make this SED ULX model available as part of this work (Section 5.4), noting that it represents only the accretion flow, and therefore does not account for any emission related to the putative donor star.We have implemented the SED ULX without selecting a unique donor star in part because the distribution of such sources for ULXs, though the subject of intense study, is still not well-constrained theoretically or observationally (e.g., Liu et al. 2004;Copperwheat et al. 2005;Tao et al. 2011;Motch et al. 2011;Soria et al. 2012;Grisé et al. 2012b;Gladstone et al. 2013;Heida et al. 2014Heida et al. , 2016Heida et al. , 2019a,b;,b;Lau et al. 2019;Wiktorowicz et al. 2017Wiktorowicz et al. , 2021)).Identifying unique optical, UV, or IR counterparts at the extragalactic distances to most ULXs remains challenging.Constraining donor stars SEDs through photometry or spectroscopy is likewise difficult given uncertainties such as the contribution of emission from the accretion flow itself, and effects of irradiation of the donor star.In the photoionization simulations that follow, the SED ULX is always coupled to the SED for a stellar population, thereby making a range of potential donor stars available to the accretor.We consider this a sensible choice for modeling the spectra or nebular emission from simple populations; however, comparison of this model to broad-band observations for an individual ULX would warrant matching the SED ULX with the spectrum for an appropriate donor star, and therefore further customization of the SED ULX model presented here.

CONSTRUCTING A PHYSICALLY CONSISTENT COMPOSITE POPULATION
In constructing a composite population with contributions from both stars and ULXs, we start from the assumption that the SXP ULX can be specified using the same parameters governing the evolution of the stellar population.Ideally, the SXP ULX contribution would be determined directly from stellar population synthesis models themselves; however, there is currently no publicly available set of models for which properties of both the stellar and accreting compact object population are easily recoverable.Absent such models, our method aims to stitch together the SXP ULX model with publicly available models for SSPs in a physically consistent way.
To specify the stellar component (i.e., SSPs) of the composite population, we opt to use Flexible Stellar Population Synthesis (FSPS) and the associated python-fsps bindings (Conroy et al. 2009;Conroy & Gunn 2010;Johnson et al. 2022).For the SSPs, we use the Binary Population and Spectral Synthesis (BPASSv2.2;Eldridge et al. 2017) models available in FSPS, which set the prescriptions for the IMF, as well as metallicity-and time-dependent properties of the spectra for the stellar population.We use the spectra from BPASS v2.2 that include products of binary evolution, noting that this version does not explicitly include accreting compact objects as part of the population.As a consequence of using the default BPASS SSPs included in FSPS, the IMF has an upper mass cutoff of 100 M , a low-mass (0.1-0.5 M ) IMF slope of −1.30, and a highmass (0.5-100 M ) IMF slope of −2.35 (matching the IMF used in the theoretical models from Fragos et al. 2013b).For the BPASS models Z = 0.02 corresponds to Z .By convention in FSPS, each SSP is normalized by stellar mass.
We construct a grid of BPASS SSPs in age and metallicity corresponding to the selected burst ages and metallicities5 at which we model the SXP ULX (t burst ∼ {1, 3, 5, 8, 10, 15, 20} Myr, and Z = {0.0001,0.002, 0.02}).Because the SXP ULX (t burst , Z) models we produce are already in terms of 1 M stellar mass, we simply add them to the corresponding SSP to produce the composite [SXP ULX + SSP](t burst , Z).This ensures some measure of physical consistency between the SSP and SXP, for example, accounting for the delay between formation of high mass stars and the appearance of the first X-ray emitting accreting BHs.This procedure also incorporates metallicity dependent effects on the SSP and SXP evolution.For example, the increase in radiative output from the SXP ULX with decreasing metallicity shown in Figure 1 is likely the consequence of weaker line-driven winds for massive stars at low metallicities, which result in less mass and angular momentum loss from the binary.This in turn results in XRBs with more massive compact objects and tighter orbits, which facilitate more efficient mass transfer, thereby producing more luminous sources (Linden et al. 2010;Mapelli et al. 2010).
Pairing a low metallicity SXP ULX (high L X /M 0 ) with a high metallicity SSP would therefore be inconsistent, i.e., would assume that massive stars with strong linedriven winds typically produce XRBs with tight orbits and massive accreting compact objects.
In Figure 4 we show the [SXP ULX + SSP](t burst , Z) models normalized to 1 M (solid lines), where each panel corresponds to the annotated instantaneous burst age.The SXP ULX and SSP components are shown as dash-dot and dashed lines, respectively.We do not display the 3 Myr burst, as it is very similar to the 1 Myr panel, prior to SXP ULX formation.In all panels, the 0.5-12 keV bandpass is highlighted via the dark blue shaded region, while the EUV range is shown in light blue.These [SXP ULX + SSP](t burst , Z) models, which span wavelengths 1-10 8 Å, are used as the spectral input for photoionization modeling with CLOUDY.

PHOTOIONIZATION SIMULATIONS WITH CLOUDY
We perform photoionization simulations using the [SXP ULX + SSP](t burst , Z) models as input to CLOUDY v17.02 (Ferland et al. 2017).To construct the CLOUDY-specific input files for our photoionization simulation grid and organize the CLOUDY output, we employ a modified version of the cloudyFSPS code (Byler 2018).Below we describe the CLOUDY simulation set-up, relevant grid parameters, and saved output from the simulations.
For the cloud geometry in all simulations, we assume a closed spherical shell with a fixed inner radius R = 10 19 cm and a constant density n H = 100 cm −3 .The cloud is ionized by a central source, which we set as the [SXP ULX + SSP](t burst , Z) SED, implying coincident mixing between the SXP ULX and SSP components.This cloud set-up is most appropriate for radiation-bounded regions, while the assumed density is appropriate to Hii regions.The chosen value for the inner radius (R) corresponds to ∼ 0.3-3 R S over the entire simulation grid, where R S is the Strömgren 6 radius.
We set i = 0 for the SED ULX component of the [SXP ULX + SSP](t burst , Z), which is equivalent to assuming an isotropic distribution of photons from the composite spectrum is incident on the cloud.This is a simplifying assumption, but reasonable for the production of high-energy ionizing photons as considered here.The outflow component of the SED ULX should be quasi-spherical, and it is this component in particular that emits strongly in the EUV (e.g., Figure 3).We consider more complex geometrical effects beyond the scope of the present work, but note that although the cloud in all simulations is subject to an isotropic distribution of photons (effective i = 0), simulation results can be post-processed with i > 0 to simulate different viewing angles for an observer.
In setting the chemical composition of the cloud, we follow the prescriptions outlined in Byler et al. (2017).In brief, we assume the gas-phase metallicity is tied to the stellar metallicity, with the gas-phase abundances set following Dopita et al. (2000).The metal abundances, except for nitrogen, are scaled to solar following the solar abundances in Anders & Grevesse (1989), while , where Q H is the hydrogen ionizing photon rate, α(Te) is the recombination rate, ne is the electron density, and n H is the hydrogen density.To obtain the range of R S from our simulations, we calculate Q H directly from the input SEDs, use α(Te) = 2.6 × 10 −13 cm 3 s −1 (corresponding to Te = 10 4 K, as in Jaskot & Ravindranath 2016), and assume ne = n H = 100 cm −3 .
the nitrogen scaling follows that of Dopita et al. (2000).We include depletion onto dust grains again using prescriptions from Dopita et al. (2000), where depletion factors do not depend on metallicity and are applied regardless of whether grains are included in the cloud.For these simulations, we do not include dust grains.
The [SXP ULX + SSP](t burst , Z) models specify the SED of the central ionizing source normalized by stellar mass.In order to set the intensity of the composite population for the CLOUDY simulations, we normalize each model by the dimensionless ionization parameter U: where fν hν dν is the rate of emitted photons capable of ionizing hydrogen (i.e., hν 0 = 13.6 eV), R is the inner cloud radius (hydrogen ionized region) in cm, n H is the hydrogen number density in cm −3 , and c is the speed of light.By this definition, the dimensionless ionization parameter can be thought of as the density of photons relative to the density of atoms.For the simulations presented here we select seven values of log U = [−4.0,−3.5, −3.0, −2.5, −2.0, −1.5, −1.0], a range that is typical for modeling Hii regions or starbursts (e.g., Dopita et al. 2000;Kewley et al. 2013).Each model is therefore specified as [SXP ULX + SSP](t burst , Z, log U) before being run through CLOUDY.
This intensity scaling is required to simulate a reasonable range of observed nebular continuum and emission line intensities, as the SED for a composite population per 1 M burst mass would produce relatively few ionizing photons.Because R and n H remain fixed7 in our models, changing log U amounts to changing Q H and therefore the effective stellar mass, as more or fewer sources are needed for older and younger bursts, respectively, to produce sufficient ionizing photons to achieve the intensity specified by log U. Following the same convention as Byler et al. (2017), we define the rate of ionizing photons per solar mass as QH ≡ Q H / M , where Q H is the rate of ionizing photons corresponding to the log U that sets the intensity for a given model in the CLOUDY simulations.With this convention, output quantities can be recovered for a composite population corresponding to 1 M stellar mass through multiplication by QH /Q H .This normalization of input or output intensities using QH /Q H is specified throughout, where appropriate.With these parameters (t burst , Z, and log U), the overall grid input to CLOUDY consists of 147 separate models, corresponding to different combinations of the seven values for t burst , three values for Z, and seven values for log U.In Figure 5, we show the [SXP ULX + SSP](t burst , Z, log U) models for select burst ages and the full range for log U for the 0.1 Z case.The models with log U = −1 have the highest intensities, or correspondingly, the largest stellar mass formed in a burst (Q H / QH ∼ 10 4 − 10 6 M for log U = −1 and t burst = 1-20 Myr).
We again employ a modified version of cloudyFSPS to create the relevant CLOUDY input files for each model in the grid, execute CLOUDY simulations, and organize the CLOUDY output files.For each grid point, we allow the CLOUDY simulations to iterate to convergence, up to a maximum number of five iterations, and set the stopping criteria for the calculations to when the cloud temperature falls below 100 K or the ionized fraction falls to 1%.Building on the nebular emission line lists presented in Byler et al. (2017) and Byler et al. (2018), we record intensities from 407 emission lines, spanning the far IR to the near UV.The full line list is included in the Ap-pendix in Table A.1, including vacuum wavelength ( Å), line name, and CLOUDY specific line ID.

RESULTS
We highlight select results illustrating the importance of the time-and metallicity-dependence of the SXP ULX implementation, and provide potential diagnostics for investigating a SXP ULX ionizing contribution.We further describe the availability of our models and simulation results.

Characteristics of Nebular Emission Due to SXP ULX and Potential SXP ULX Diagnostics
A key result from the photoionization simulations is that the inclusion of the SXP ULX component does not significantly alter the broad-band colors in the UV-to-FIR relative to the case of the BPASS SSPs alone.This can easily be seen by looking at the enhancement in total nebular (line and continuum) emission due to the addition of the SXP ULX .In Figure 6 we show the fractional contribution of the nebular emission to the total flux as a function of wavelength and burst age in the top panel, and the percentage enhancement in this nebular flux due to the addition of the SXP ULX in the bottom panel.The nebular contribution to the total flux is high (> 50%) at all wavelengths on short burst timescales (≤ 3 Myr, before the SXP ULX has turned on), and in the MIR for all burst timescales simulated here.However, as the bottom panel illustrates, the enhancement in nebular emission due to the inclusion of the SXP ULX is typically very small (< 5%).In the broad-band sense, the major effect of the SXP ULX is therefore to slightly prolong the ionizing output relative to the SSP alone, a consequence of the imposed age-dependence of our models.
Due to the shape of the SED ULX as extended through the EUV, a key outcome from the addition of the SXP ULX is a change to the intensity of lines with excitation potentials ≥ 54 eV.The CLOUDY line list for our simulations (Table A.1) includes 63 lines redward of Lyα with ionization potential ≥ 54 eV (Kramida et al. 2022), of which 28 have ionization potentials in excess of 90 eV, including a number of so-called "coronal lines" (i.e., forbidden transitions with very high ionization potentials).In order to establish which of these highexcitation lines may be used as diagnostics for SXP ULX ionization, we consider lines redward of Lyα that satisfy the following criteria: ( 1 The fractional contribution of the nebular line and continuum emission to the total flux (stellar and nebular) as a function of wavelength and instantaneous burst age for models with Z = 0.1 Z and log U = −1.Bottom: The percentage enhancement in nebular line and continuum flux due to the addition of the SXPULX.For all burst timescales simulated, the SXPULX-driven enhancement in the total nebular emission is typically small (< 5%).The black labels denote high-ionization species whose nebular line intensity is particularly enhanced due to the inclusion of the SXPULX component.
or Paβ8 (i.e., log( (2) have a median enhancement in line flux by a factor of 2× relative to the SSP-only models (log( f line,SXP+SSP / f line,SSP ) ≥ 0.3).There are only a select few emission lines satisfying both of these criteria: He II λ1640,4686, and [O IV] 25.8832µm.We annotate these lines in Figure 6, and note them in Table A.1 (* flag).In Figures 7-9 we show emission line diagnostic diagrams in the UV, optical, and IR, respectively, to highlight the potential for uncovering ionization due to an SXP ULX using these lines passing both the criteria as outlined above.In general, diagnostic diagrams such as these are constructed using lines relatively close in wavelength to reduce the effects of reddening and make the lines used in the diagnostic accessible with a single instrument, reducing uncertainties due to crosscalibration.In all line diagnostic diagrams that follow, the [SXP ULX + SSP](t burst , Z, log U) models (i.e., models with SXP ULX contribution) are shown as the solid lines, where orange to dark blue lines represent metallicities from Z = 0.01-1 Z and light blue to dark blue lines represent ionization parameters from log U = −4 to −1.For all grids in the diagnostic plots, the star denotes the model with lowest metallicity (0.01 Z ) and highest ionization parameter (log U = −1).The BPASS SSP- only models are shown as slightly transparent, dashed lines in the background with the same color scheme.In cases where the models with SXP ULX contribution are indistinguishable from SSP-only models (e.g., burst ages < 5 Myr), only solid lines are shown.
For the UV diagnostic in Figure 7, the (O III] λ1661,6)/(C III] λ1907,9) ratio is sensitive to the ionization parameter, C/O ratio (i.e., abundance pattern, see discussion in Section 5.3), and uses lines that are relatively strong in the UV (Byler et al. 2018).The He II λ1640/C III] λ1907,9 ratio in this same diagram is sensitive to the hardness of the ionizing spectrum.For burst timescales where the SXP ULX contributes (≥ 5 Myr), this diagnostic illustrates that the ionizing spectrum is harder at all metallicities relative to the SSP-only models due to the inclusion of the SXP ULX , resulting in stronger intensity of He II λ1640 relative to C III] λ1907,9 at the same ionization parameter.
For this set of line diagnostics, we compare the [SXP ULX + SSP](t burst , Z, log U) grid to a sample of local blue compact dwarf (BCD) galaxies from Senchyna et al. (2017) for which these same lines have been measured.While most of the data points from the BCDs (cyan diamonds) are consistent with both the SSP-only and [SXP ULX + SSP](t burst , Z, log U) grid for ionization parameters log U ≥ −3.0, the inclusion of the SXP ULX slightly broadens the range of metallicities consistent with the data points; whereas the BCD sample is consistent with the SSP-only grid for Z ≤ 0.1 Z , the composite grid extends the consistency to 0.1 Z ≤ Z < Z for older burst ages.This BCD sample has a median gas-phase metallicity of ∼ 0.15 Z , suggesting that if SXP ULX ionization is important in these galaxies the population may be dominated by star formation bursts ≥ 10 Myr; however, given the relative sparsity of the metallicity sampling in our grids and in the absence of a more detailed age analysis-which we consider beyond the scope of this work-it is not clear from this qualitative comparison whether the SSP-only or SXP ULX grid is preferred for such galaxies.Nonetheless, the comparison serves to illustrate how the addition of the SXP ULX alters the parameter space spanned by the grid, therefore changing the range of physical properties consistent with the observed population.In the case of Figure 7, the addition of the SXP ULX hardens the ionizing spectrum at a given log U, which in effect mimics an SSP with slightly lower metallicity.
For an optical emission line diagnostic we present Figure 8, where the [N II] λ6584/Hα ratio is sensitive to metallicity and log U (among other parameters, e.g., Kewley et al. 2013), and the He II λ4686/Hβ ratio again probes the hardness of the radiation field.Here we compare to a sample of star-forming galaxies with strong nebular He II emission from SDSS identified by Shirazi & Brinchmann (2012).A little under half of this sample lack Wolf-Rayet (WR) star features in their spectra ("Non-WR"; red triangles), while the remaining galaxies show evidence for broad emission line features due to the presence of WR stars ("WR"; blue circles).In general, SSP-only models have a difficult time reproducing the observed range of He II λ4686/Hβ, as the dashed grids in Figure 8 illustrate.The [SXP ULX + SSP](t burst , Z, log U) grid is capable of reproducing the observed intensities for a subset of these data points, notably those without clear WR features in their spectra, but only for t burst > 10 Myr and lower metallicity (Z ∼ 0.1 Z ).In this case, the ability of the [SXP ULX + SSP](t burst , Z, log U) grid to explain a subset of the galaxies with strong He II/Hβ where SSP-only models cannot is again a consequence of a harder radiation field due to the ad- The models with SXPULX contribution overlap with the observed intensities for some of the Non-WR feature galaxies, particularly for t burst > 10 Myr and Z ≤ 0.1 Z ; however, neither the [SXPULX + SSP](t burst , Z, log U) nor the SSPonly grids are capable of reproducing some of the more extreme galaxies, particularly those with detected WR features and higher N II/Hα near the empirical maximum starburst line.
dition of the SXP ULX .However, the magnitude of the increase to the ionizing photon rate due to the SXP ULX is strongly dependent on SFH, as we discuss in more detail in Sections 5.2 and 6.1.Interestingly, the observed WR feature galaxies are difficult to reproduce using either the models with SXP ULX ionizing contribution or SSP-only photoionization.This is primarily because such galaxies have large ratios of both He II λ4686/Hβ, implying hard ionizing spectra, and [N II] λ6584/Hα, implying higher metallicities.We explore metallicity in the context of this diagnostic in Section 5.3 and alternative ionizing sources in Section 6.3.2, and here comment instead on the effect of stellar wind contamination of nebular emission.Lines such as He II can be produced either via stellar winds (broad emission line component) or photoionization (narrow emission line).Moderately good spectral resolution is therefore needed to distinguish the broad stellar from narrow nebular components for such lines (Brinchmann et al. 2008).If not fully resolvable, residual emission from the broader stellar wind component could enhance the line intensity assumed to be purely nebular.It is likely that stellar wind contamination, if present, operates preferentially to inflate line intensities at higher metallicities, where line driven winds from massive stars are stronger (Vink et al. 2001).It is therefore important to consider the presence and magnitude of potential stellar wind contamination in observations of such nebular line species, particularly at higher metallicities, when comparing to the purely nebular line emission determined from photoionization simulations (e.g., Byler et al. 2018).
In the IR, we consider the diagnostic presented in  (2022) for ionization due to pure star-formation (SF), SF/active galactic nuclei (AGN), and pure AGN, as well as the results from photoionization simulations for a 10 3 M IMBH with fractional contribution 4-16% relative to a 20 Myr BPASS SSP (described in more detail in Section 6.3.2).The [SXP ULX + SSP](t burst , Z, log U) grid is largely consistent with the pure SF region in this diagnostic, though offset from the SSP-only grid, again demonstrating how the addition of the SXP ULX hardens the ionizing spectrum.The models including SXP ULX contribution veer into the composite region for some of the solar metallicity t burst = 5 Myr models, and into the AGN ionization region only for lowest metallicities ( 0.1 Z ) and highest ionization parameters (log U > -3) for bursts with t burst > 10 Myr.This indicates, not surprisingly, that the [SXP ULX + SSP](t burst , Z) photoionization signature looks less like pure SF when the fractional contribution of the SXP ULX component relative to the SSP is most pronounced.As in Figure 7, this occurs when the strongest ionizing components of the stellar population are depleted, namely WR stars (∼ 5 Myr at roughly solar metallicities) and the most massive stars that will explode as supernovae or collapse to a BH (> 10 Myr).We note that the SXP ULX grid does not significantly overlap with IMBH grids, suggesting some power in this emission line diagnostic for distinguishing between ionization due to BHs separated by orders of magnitude in mass; however, this too depends on how the accreting BH (at any mass) is scaled relative to the stellar population, a point that we return to in more detail in Section 6.3.2.Finally, we comment on high-ionization lines that do not satisfy both criteria in terms of strength relative to Hβ or Paβ and enhancement relative to the SSP-only case.These lines fall broadly into three categories: (1) lines that are strong relative to Hβ or Paβ, but not significantly enhanced relative to the SSP-only models (flagged with † in Table A.1); (2) lines that are weak relative to Hβ or Paβ, but are significantly enhanced relative to the SSP-only models (flagged with ‡ in Table A.1); or (3) lines that have zero flux for thsee models (flagged with || in Table A.1). Lines flagged with † in Table A.1 have a median log( f line / f Hβ ) ∼ −3 and median enhancement of 1.004× in flux relative to the SSPonly models.By contrast, lines flagged ‡ in Table A.1 have log( f line / f Paβ ) ∼ −9, but median enhancement of 47000× in line flux relative to the SSP-only models.The latter case encompasses many coronal lines with ionization potentials > 90 eV, i.e., the regime where the SSP provides relatively little flux, but the SXP ULX component is substantial.As such, there is a higher diagnostic potential for discerning the SXP ULX ionization contribution using coronal lines, many of which are in the IR, but only if the fractional contribution of the SXP ULX relative to the SSP is high enough to boost the line strengths into an easily detectable range.
We illustrate these two cases-moderately strong lines only marginally enhanced by the SXP ULX and weak lines heavily enhanced due to the addition of the SXP ULXvia additional emission line diagnostics in Figures 10-12.In Figure 10, we show the [SXP ULX + SSP](t burst , Z, log U) grid relative to the same sample of local BCD galaxies from as in Figure 7. Here, the grid with SXP ULX contribution is only slightly offset to higher values of C IV λ1548,51/O III] λ 1661,6 compared to the SSP-only grid.The lack of overlap between some of the galaxies and the grids in this diagnostic is likely a function of the adopted C/O ratio for the simulations in this work, which we discuss in more detail in Section 5.3.
In Figures 11-12, we show emission line diagnostics in the optical and IR that include the high-ionization line [Ne V].For these emission line diagnostics, the addition of the SXP ULX has a profound effect on the grid relative to the SSP-only case, given the ratio of [Ne V]/[Ne III] compares lines with ionization potentials of of ∼ 126 eV and 64 eV, respectively.In Figure 11, the SSP-only models rarely reach log([Ne V] λ3426/[Ne III] λ3870) > −5, while the grid with SXP ULX contribution begins to populate the composite region for low metallicities and high-ionization parameters (classification regions from Cleri et al. 2023).In Figure 12, the [SXP ULX + SSP](t burst , Z, log U) grid traces out a relatively narrow vertical track, distinct from the SSP-only and IMBH photoionization cases.However, for both of these diagnostics, the potential power for discerning SXP ULX contribution to the ionizing photon budget comes at the expense of line strength.The ratio of [Ne V] 14.3µm/[Ne III] 15.56µm traced out by the grid with SXP ULX contribution spans roughly 10 orders of magnitude.Such lines are therefore undetectable for pure star-formation ionization and only potentially detectable (e.g., with JWST) with an SXP ULX ionizing contribution for a narrow range in log U, t burst , and Z.This implies that very high-ionization species such as [Ne V] or [Ar V] may be more reliable diagnostics of AGN or IMBH photoionization rather than SXP ULX ionizing contribution; however, the potential selection biases inherent in using particular lines should be considered in assessing the reliability of any such classification diagnostic (e.g., Richardson et al. 2022).

Consequences of Time Dependence of the SXP ULX
Relative to the SSP The key parameters affecting the time-dependence of the nebular emission from the [SXP ULX + SSP](t burst , Z, log U) grid are (1) the selection to impose a delay time for the SXP ULX to turn on and; (2) the combination thereafter with SSPs of corresponding t burst .The lack of overlap between some of the galaxies and either grid for this diagnostic is likely a function of adopted C/O ratio for the photoionization modeling in this work.
Imposing a delay time creates a scenario in which the SXP ULX prolongs the ionizing output of the population beyond that of the products of binary evolution already included in BPASS, as would be expected if the SXP ULX component is indeed descended from the stellar component, allowing time for compact object formation and accretion to begin.In this way, the SXP ULX , or accreting compact objects in general, are another way of rejuvenating the ionizing power from a population.In spite of this, the magnitude of the additional contribution of the SXP ULX to the emergent nebular emission has time-dependent limits by virtue of being coupled to the time evolution of the SSPs.Although the absolute L X /M 0 scaling for the SXP ULX is largest for younger instantaneous bursts at all metallicities (Figure 1), the SXP ULX generally hardens the incident ionizing spectrum more profoundly relative to the BPASS SSPs at older burst ages.That is, the fractional contribution of the SXP ULX to the ionizing photon budget is typically highest for burst ages > 10 Myr in the models considerd here.This corresponds to the timescale on which the harder components of the stellar ionizing spectrum have been depleted (e.g., WR and other very massive stars), while the SXP ULX component remains appreciable.Given the imposed delay time and coupling with BPASS SSPs, the SXP ULX contribution to the total rate of ionizing photons is 2% for burst ages between 5- ) is largest at later times (> 10 Myr), particularly for the lower metallicity models.This corresponds to the timescale on which the ionizing contribution from the stellar population has decreased significantly, while the SXPULX ionizing contribution remains appreciable.
10 Myr, and increases only to 5-7% on timescales > 10 Myr.The consequences of this time dependence are especially evident in Figure 11, where the local EELGs are consistent with the models with SXP ULX contribution for a narrow range of the grid corresponding to ∼ 0.1 Z , log U > −2, and t burst > 10 Myr.This suggests SXP ULX photoionization could contribute to their observed line intensities under this specific set of circumstances; however, such a range of physical properties (e.g., high log U coupled with older burst ages) may not be applicable in general to EELGs.For example, two of the galaxies in Figure 11 from have inferred ages 10 Myr, and may therefore require a source other than an SXP ULX or SSP to explain the observed line intensities, as discussed in Olivier et al. (2022).
In general, our models suggest select high-ionization emission lines can be produced by ionization due to an SXP ULX for t burst > 10 Myr and low metallicities.For SFHs that include a young (< 10 Myr) burst component, the SXP ULX contribution to the ionizing photon budget is not significant enough to reproduce the observed line ratios that include high-ionization lines such as He II or [Ne V].This is because, for the models considered here, the pure stellar ionizing continuum is substantial for t burst < 10 Myr, resulting in strong emission from lines such as Hβ.This strong stellar ionizing flux effectively dilutes the SXP ULX contribution, which primarily affects the high-ionization line (e.g., He II) rather than hydrogen ionization.We illustrate this point in two ways.First, in Figure 13 we show nebular emission line luminosities normalized by stellar mass for He II λ4686, Hβ, [Ne V] λ3426, and [Ne III] λ3870 as a function of burst age.In comparing the evolution of the intensi-ties of He II and Hβ, it is clear both lines reach their absolute maximum in intensity for t burst ∼ 3-5 Myr; however, while the Hβ intensity drops by almost two orders of magnitude thereafter, the evolution of the He II intensity is much flatter over the same time span, owing to the SXP ULX contribution, resulting in a higher ratio of He II λ/Hβ at later times.The same is true for [Ne V] λ3426 relative to [Ne III] λ3870.For this ratio, it is also particularly evident when the SXP ULX contribution turns for t burst > 3 Myr, as there is no appreciable [Ne V] emission prior to this.
To further illustrate the timescale dependent nature of the relative SXP ULX ionizing contribution, we reran the photoionization simulations without the strict delaytime dependence for the SXP ULX component.For these "no-delay time" simulations, we extrapolated L X /M 0 from Figure 1 to burst ages 1-3 Myr, allowing the SXP ULX to contribute to the ionizing budget at every time step.As Figure A.1 shows, the immediate formation of the SXP ULX alongside a 1 Myr or 3 Myr burst only weakly increases He II/Hβ, given that the strength of the stellar ionizing continuum is still quite substantial on these timescales, as in Figure 13.
In this way, the addition of the SXP ULX may contribute significantly to high-energy ionizing photon production only in some environments, at least for the SED ULX considered in this work.Due to the timedependent nature of the SXP ULX ionizing contribution from the models presented here, measured SFHs or other proxies for age (e.g., strong line equivalent widths) should be considered alongside emission line diagnostics when comparing with these models to discern SXP ULX ionizing contribution.

Consequences of Metallicity-Dependence of the SXP ULX and Adopted Abundance Patterns
The photoionization simulation results are based on the assumption that the metallicity of the SXP ULX corresponds to the stellar metallicity of the SSP, as would be expected if the SXP ULX and SSP form and evolve from the same parent stellar population.Consequently, a low-metallicity SSP, which already has a somewhat harder ionizing spectrum, will always be combined with an SXP ULX with higher formation efficiency (L X /M 0 ) than the SXP ULX combined with a solar-metallicity SSP.
The imprint of the metallicity-dependent normalization of the SXP ULX is particularly evident in emission line diagnostics that include metallicity-sensitive line ratios (e.g., [N II]/Hα).To illustrate the magnitude of this effect, we run simulations where we remove the correlation between the SXP ULX normalization and stellar metallicity.In these "no stellar metallicity dependence" simulations, we allow an SSP at a given t burst to be coupled with an SXP ULX normalized using L X /M 0 from the orange curve (i.e., maximum) in Figure 1, regardless of the SSP metallicity.The results are shown in Figure A.2, where higher ratios of He II λ4686/Hβ are produced at a given [N II]/Hα, particularly for the solar metallicity case.This results in much more overlap between the grid with the SXP ULX contribution and the observed data points, for a broader range of burst ages; however, we note that such an elevated L X /M 0 for ULXs at solar metallicities is not supported by observations.Even allowing for variations due to stochastic sampling of the XLF, sources with L X > 10 39 erg s −1 are more rare at higher metallicities, due to steepening of the bright end of the XLF with increasing metallicity (Lehmer et al. 2021).We therefore include this "no stellar metallicity dependence" example simply to illustrate the effects of our model choices.
While the SXP ULX may not be an efficient additional source of ionizing photons at high metallicity given the relations in Figure 1, other sources of high-energy ionizing photons may operate preferentially in this regime.For example, as pointed out in Shirazi & Brinchmann (2012), the sub-sample of galaxies with WR features in Figure 8 generally has higher metallicity, which correlates with higher stellar mass and SFR, and therefore potentially stronger stellar feedback.For such galaxies, this stronger stellar feedback (i.e., winds and supernovae) may drive shocks, which could be an additional source of high-energy ionizing photons (e.g., Section 6.3.3).
The choice to couple stellar and gas-phase metallicities and the adopted abundance patterns for these sim-ulations are also of critical importance for select emission line diagnostics.For example, though L X /M 0 increases with decreasing metallicity thereby increasing the number of high-energy ionizing photons, at very low metallicities the absolute abundance of some elements will be low enough that this becomes the dominant factor in setting the line intensity (e.g., discussion in Senchyna et al. 2020).This can easily be seen in the shape of the [SXP ULX + SSP](t burst , Z, log U) grid in Figure 12, where the lowest metallicity models including SXP ULX contribution do not uniformly have the strongest [Ne V] 14.3µm/[Ne III] 15.56µm ratios, as would be expected if SXP ULX production efficiency (i.e., L X /M 0 ) were the dominant factor in determining line strength.
Likewise, the adopted C/O ratio in the simulations affects the extent of the grids as presented in Figures 7 and  10.Our abundance pattern corresponds to log(C/O) = −0.45,which is a factor of approximately two higher than the average log(C/O) found for metal-poor, starforming dwarf galaxies from Berg et al. (2019), though there is a large (0.17 dex) intrinsic dispersion in this value.For the emission line diagnostic in Figure 7, adopting a lower C/O ratio would shift the grids to the right (Jaskot & Ravindranath 2016), thereby encompassing all the observed BCDs.Likewise, a lower C/O ratio would shift the grids in Figure 10 up and to the left (Byler et al. 2018).We consider a more detailed treatment of the metallicity-dependence of the C/O and N/O ratios in the context of SXP ULX ionization out of the scope of the present work, and refer the reader to Appendix B in Byler et al. (2018) for an excellent demonstration and discussion of the effect of differing abundance patterns on emission line diagnostics.We note that care should be taken in comparing observed line ratios to these simulations if using lines that are particularly sensitive to the assumed abundance pattern, such as C and N.

Availability of Models
We make our models and simulation results available in two separate formats, as outlined below.
The models are available directly from FSPS9 .From FSPS, one can build SEDs for a population with or without SXP ULX contribution, including consideration of nebular line and continuum emission.In this way, users can simulate spectro-photometric data for populations including SXP ULX contribution for use in SED fitting codes (e.g., Johnson et al. 2021;Doore et al. 2023).
We also provide the results through a standalone github repository 10 .In this repository we provide code examples in python for reading and plotting simulation results, allowing users to quickly plot the [SXP ULX + SSP](t burst , Z, log U) grid for different emission line diagnostics.Users can therefore compare data or other model grids to these simulations including SXP ULX contribution.In this repository, we likewise provide Table A.1 in a machine-readable format, such that lines can be selected on the basis of relevant flags included in the table notes.

DISCUSSION
We have already presented potential UV, optical and IR emission line diagnostics for SXP ULX ionization in Section 5.Here we present best practices for diagnosing SXP ULX ionization when using X-ray observations, and summarize the results from this work relative to recent literature results for photoionization due to high-energy sources.We additionally present a brief discussion of alternative sources of hard ionizing photons and their connection to the SXP ULX model presented here.

Caveats For Applying SXP ULX Formalism to X-ray Observations
The emergent L X (0.5-8 keV) from the SXP ULX models is very similar to the input L X (i.e., Figure 1), as the hard photons are attenuated very little by the cloud.However, comparison of the simulation results with Xray observations requires a more careful transformation of the intrinsic L X from the simulations into a model for observed X-ray counts.We must consider of factors such as line of sight absorption (e.g., Wilms et al. 2000), viewing angle (e.g., Abolmasov et al. 2009;Kovlakas et al. 2022), source variability (e.g., Earnshaw et al. 2018), and instrumental effects.Such considerations likely require a full forward-modeling approach (e.g., Gilbertson et al. 2022), but could allow the use of relatively low count X-ray data in spectro-photometric fitting (e.g., Doore et al. 2023).Nevertheless, implementation will be non-trivial, and we defer full consideration to a future work.Here, we instead provide caveats relevant for interpreting the simulation results with respect to the observations, given factors affecting X-ray detectability.
A primary consideration in interpreting L X values from the photoionization simulations is that the SED ULX model assumes a specific geometry appropriate to a supercritical accretion flow.In such a model, the wind or outflow component creates an effective funnel geometry, which collimates and scatters the harder 10 https://github.com/kgarofali/sxp-cloudyinner disk emission for face-on viewers and may obscure this component for more edge-on viewing angles (e.g., Begelman et al. 2006;Poutanen et al. 2007;King et al. 2023).Under this geometry, L X (obs) is likely dependent on inclination angle, and possibly also on the mass supply rate itself, which in the SED ULX model is set by the fixed funnel opening angle.
The effects of changing inclination angle and mass supply rate are illustrated in Figure 3.We show the adopted SED ULX in black with different linestyles, ranging from a face-on viewing angle (i = 0 • ) to the edge-on case (i = 90 • ), wherein the SED ULX is primarily an EUV source.However, these are all for a fixed funnel opening angle (θ), tied to the mass transfer rate and BH mass.We therefore also show the effect of changing the funnel opening angle and corresponding mass supply rate in green for a small opening angle, and orange for a large funnel opening angle (both for i = 0 • , corresponding to a face-on viewing angle).These two extremes illustrate that the shape of the SED is also, not surprisingly, highly dependent on the assumed θ (or, correspondingly, ṁ): for a small funnel opening angle, where the mass transfer rate is very high, the peak of the SED shifts into the UV/optical owing to the dominance of the outflow component, while for a large funnel opening angle the peak shifts to the soft X-ray regime, and there is a loss of the very flat EUV component due to lack of relatively isotropic outflow component and associated selfirradiation within the wind funnel.Finally, we relax the assumption that the funnel opening angle is tied to the mass supply rate, allowing for a high mass transfer rate and modest funnel opening angle, shown by the blue line in Figure 3.This results in an intrinsic SED ULX with the strongest EUV component, as well as an appreciable optical/UV contribution, though it is unclear if such a case is physical.
The diversity of observed spectral shapes for ultraluminous sources may therefore be explained by differences in both funnel opening angle and viewing angle.For example, ultra-luminous supersoft sources (e.g., Soria & Kong 2016; Urquhart & Soria 2016) may represent accreting compact objects with modest mass supply rates and large funnel opening angles viewed edgeon, while ultra-luminous UV sources may be compact objects with higher mass transfer rates likewise viewed edge-on (e.g., Kaaret et al. 2010).
Unfortunately, the dependence of mass transfer rate (or lack thereof) on funnel opening angle, and corresponding distribution of funnel opening angles is not well-known for either the intrinsic or the observed ULX population.Some binary population synthesis models suggest that the dominant accretor population (i.e., BH versus NS) and the timescale on which they appear will affect the fraction of the population that is strongly beamed (e.g., Wiktorowicz et al. 2019).In this respect, the form of the SED ULX used in this work is not exhaustive.As Figure 3 illustrates, it is possible an SXP could be an even more extreme EUV source under certain assumptions about the mass supply rate and distribution of funnel opening angles for the population.Different assumptions about typical values for ULX parameters such as θ and ṁ will therefore propagate to the intensities of high-ionization nebular emission lines.In practice, this could work to either increase or decrease the line intensities, or perhaps even alter the strict agedependence for particular line ratios.A more detailed investigation of the effect of these parameters on the resultant nebular emission is therefore beyond the scope of the present work.
Despite these caveats, we demonstrate the applicability of our models to observations in terms of line intensities and SXP ULX L X via a comparison of simulated and observed He II λ4686/Hβ as a function of L X /SFR (a typical observable for HMXB and ULX populations; e.g., Mineo et al. 2012;Kovlakas et al. 2020).For the observational comparison set, we use a sample of starforming dwarf galaxies for which there are spectra from SDSS, HST /COS, and MMT, and X-ray coverage from Chandra (Senchyna et al. 2020).Given the aforementioned caveats to X-ray detectability, comparison of our simulations with this observed sample rests on the assumption that the observed L X corresponds roughly to the L X from the intrinsic (simulated) population.Additionally, given the discussion in Section 5.2, we must consider whether the SFH assumed for the simulated results is appropriate to the observed sample (i.e., line intensities, SFRs, and L X are time-averaged in a similar way).For the observed sample, line intensities are derived from the spectroscopic data, and the corresponding SFRs and values for L X (or upper limits) are determined from the same spectroscopic aperture, or within 1.4 of the aperture for the X-ray data.SFRs are measured using calibration constants derived from BPASS appropriate to the effective age and metallicity of the stellar population assuming continuous star formation.
To produce a simulated comparison set, we consider both burst and continuous SFH cases.The burst properties are easily recoverable from the simulations, as shown in figures throughout.To recover the continuous SFH case, we integrate the simulated line intensities and SXP ULX L X (both in terms of stellar mass) assuming a constant SFR of 1 M yr −1 .The results are shown in Figure 14, with the continuous SFH models as circles and the burst models as shaded regions, where  2020) are shown as green diamonds.The low metallicity simulations are consistent with all but the most extreme observed values and upper limits for the line ratios and LX/SFR; however, this requires instantaneous bursts where t burst > 10 Myr.Models with continuous SFHs produce LX/SFR in the range of the observed galaxies and upper limits, but much lower ratios of He II λ4686/Hβ.
the extent of the region encompasses the range of simulated log U and t burst .The most extreme line ratios are achieved for high log U and older t burst .The simulations can reproduce much of the observed range of He II λ/Hβ (green diamonds) for reasonable L X /SFR (i.e., overlap with detections, or are below upper limits); however, this is true only for instantaneous bursts of star formation where t burst > 10 Myr.For the continuous SFH case, the models produce much more modest He II λ/Hβ, as discussed in Section 5.2, and there is much less overlap with the observed sample.This demonstrates that the assumed SFH is of critical importance in such comparisons.SFHs recovered from SED fitting (e.g., Johnson et al. 2021), or other proxies for age such as Hα or Hβ equivalent width, are therefore integral to determining the importance of SXP ULX high-energy ionizing contribution.As some of the observed galaxies in Figure 14 have values of L X /SFR in excess of the simulated results, we briefly comment on the effects of sampling in making such comparisons.In general, testing for agreement between the simulated results for L X and an observed population is best achieved when the observed population has a well-sampled XLF.This typically holds in the high SFR regime ( 3 M yr −1 ), where the effects of stochastic sampling of the XLF are on the order of, or smaller than, the model uncertainties on L X , and where the strict linear scaling of L X with SFR holds (Grimm et al. 2003;Justham & Schawinski 2012;Lehmer et al. 2021).In the low SFR regime, on the other hand, integrated L X (obs) can be subject to large (∼ 0.7 dex) scatter due to stochastic sampling of an XLF that, for HMXBs and ULXs, flattens at the bright end (e.g., Lehmer et al. 2021).Similar but far less severe considerations about statistical sampling apply when determining and comparing SFRs.Empirical or model-derived calibration constants for different SFR indicators are based on the assumption of a fully populated IMF and SFR indicators which have equilibrated.In other words, typical calibration constants are appropriate for regions physically large enough for the IMF to be well-sampled and for which the SFH has been continuous long enough for a given SFR indicator to reach equilibrium (e.g., 5-100 Myr Kennicutt 1998;Kennicutt & Evans 2012).The galaxies in Figure 14 that are X-ray detected with the highest L X /SFR ( 10 41 erg s −1 ) are some of the most distant galaxies (D > 200 Mpc) in the sample, such that the spectroscopic apertures encompass large areas on the sky, making it unlikely that statistical sampling affects the measured SFRs.These galaxies have measured SFRs in the range 0.003-6 M yr −1 , and therefore may be subject to significant sampling effects at the low SFR end; however, stochastic sampling is likely not enough to explain their high X-ray production efficiencies, even at the very low SFR end, making them ∼ 2σ outliers relative to empirical scaling relations at low metallicities (Lehmer et al. 2021).These galaxies may therefore represent cases where alternative high-energy ionizing sources with L X /SFR well in excess of average XRB scaling relations contribute to the observed line ratios (e.g., Section 6.3.2).
Taken together, the discussion presented here implies that in order to determine whether L X from a given SXP ULX is consistent with observations requires modeling the effects of inclination, both on the observable spectral components and their apparent luminosities (e.g., Kovlakas et al. 2022), as well as consideration of the effects of sampling.The former is a critical consideration in the context of these simulations, as it implies an SXP ULX may be intrinsically very luminous, and therefore efficient at producing high-energy ionizing photons, while potentially being difficult to detect in X-rays.The latter is important because all line intensities (either for bursts or calculated for continuous SFH) are simulated assuing a population average L X .For highly starforming galaxies, the population average L X may be a good approximation; however, if these conditions do not hold for an observed sample, comparison with simulation results may result in erroneous inferences about population age, metallicity, and/or source(s) of ionization.

Comparison with Recent Literature
There have been a number of works in the literature exploring the efficacy of X-ray sources in the production of high-energy ionizing photons, with a particular recent emphasis on the X-ray contribution to production of He II λ4686 in high-redshift galaxies and their analogs (e.g., Garnett et al. 1991;Schaerer et al. 2019;Senchyna et al. 2020;Simmonds et al. 2021;Oskinova & Schaerer 2022;Kovlakas et al. 2022;Umeda et al. 2022;Ramambason et al. 2022).The most recent investigations fall broadly into two categories: (1) full photoionization simulations of the ionizing effect of an X-ray source coupled to a stellar population; and (2) analytic approximations for the X-ray source contribution to select line intensities, such as He II λ4686.We focus on the former here for parity with our own implementation.
In Figure 15 we show the SED ULX for a 20 Myr instantaneous burst of star formation at 0.1 Z in yellow.
Here the SED ULX has been normalized to 1 M stellar mass formed (Section 3), corresponding to L X (0.5-8 keV) ∼ 3 × 10 32 erg s −1 .The corresponding SSP (t burst = 20 Myr, 0.1 Z ) is shown in green, while the composite of these two components is shown in black.The additional labeled SEDs ("MCD", "DIS", "BB", and "I Zw18") are selections from the literature that have previously been used in photoionization simulations.Below, we briefly summarize salient points from recent investigations using these SEDs, and compare to our findings based on the SXP ULX framework.
The photoionization simulations performed in Senchyna et al. (2020) explored the effect of HMXBs on production of He II λ4686 and [Ne V] λ3426.The authors couple an HMXB model, represented via a multi-color disk SED ("MCD", green dotted line in Figure 15) at different BH masses (∼ 10-100 M ) to BPASS stellar populations (Z = 0.001-0.020)assuming continuous star formation.The HMXB component is scaled to the stellar population for a range of X-ray production efficiencies (L X /SFR = 10 40 − 10 44 erg s −1 (M yr −1 ) −1 ).Much like our approach, they used a modified version of cloudyFSPS to initialize the CLOUDY input grid, employing a similar range in log U.With these simulation inputs, the authors conclude that the observed range of He II λ4686/Hβ can only be reproduced for extremely high X-ray production efficiencies (L X (0.5-8 keV)/SFR > 10 42 erg s −1 (M yr −1 ) −1 ).Such extreme X-ray production efficiencies are not observed, even in nearby extremely metal-poor galaxies, and are likely excluded by reionization-era constraints (Lehmer et al. 2021;HERA Collaboration et al. 2023).In Figure 14 we compare our model line predictions with the observational sample from Senchyna et al. (2020) which they used, in part, to draw the conclusion that HMXB ionization is insufficient to explain the observed range of He II λ4686/Hβ.As already discussed in Section 6.1, our simulations are able to reproduce more of the observed range of He II λ4686/Hβ within the observed range of L X /SFR, but primarily for instantaneous bursts with t burst > 10 Myr at low metallicities.Under the assumption of a continuous SFH, as in Senchyna et al. (2020), our models are not able to reproduce the more extreme values of He II λ4686/Hβ at the simulated range of L X /SFR.This echoes the results presented in Senchyna et al. (2020), underscoring the importance of the assumed SFH in setting line ratios; however, even for similar assumptions about the SFH, our models produce stronger He II λ4686/Hβ at a given L X /SFR as compared with Senchyna et al. (2020), which we attribute to differences in the assumed SED for the X-ray photoionizing component.Our SED ULX is much flatter through the EUV than the multi-color disk model employed, and therefore produces more ionizing EUV photons per L X .
The importance of the form of the SED to the production of high-energy ionizing photons was highlighted in the CLOUDY photoionization simulations from Simmonds et al. (2021).These authors explored the photoionizing effect of ULXs using several empirical models for the ULX SED coupled to a BPASS SSP (1 Myr instantaneous burst with Z = 0.005 11 ) again for a similar range in log U as employed here.In their simulations, the ULX contribution is scaled relative to the BPASS SSP for a range of L X /SFR, motivated by empirical scalings and scatter therein.These authors find they can reproduce some of the observed range of He II λ4686/Hβ for reasonable X-ray production efficiencies (L X /SFR 10 40 − 10 41 erg s −1 (M yr −1 ) −1 ), roughly in agreement with our findings for instantaneous bursts and the SED ULX model.
11 These authors list the metallicity of the BPASS SSP as 12 + log(O/H) = 8.1.For their adopted abundances from Grevesse et al. (2010), this corresponds to ∼ 0.3 Z , which we assume corresponds to the BPASS Z = 0.005 model.

Wavelength [ Å]
Figure 15.The intrinsic SED for the composite model (SXPULX + SSP, solid black line) with t burst = 20 Myr and Z = 0.1 Z , normalized to 1 M stellar mass formed.The stellar population and SEDULX components of the composite are shown in green and yellow, respectively.Various SEDs from the literature that have previously been used in photoionization simulations are shown in different colors and linestyles for comparison, all normalized to the same LX (0.5-8 keV) (or, in the case of the "BB", LEUV) as the SEDULX.The grey vertical lines mark ionization potentials for select lines, as labeled.
In general, the simulations results presented in Simmonds et al. ( 2021) also achieve higher values of He II λ4686/Hβ than our grids with SXP ULX contribution, at least for their most optimistic SED model ("DIS" from Berghea & Dudik 2012, as shown via the orange dashed line in Figure 15).We attribute differences in high-energy ionizing photon production efficiency between simulations, ostensibly both for ULXtype sources, in part to differences in adopted SED shape, and in part to treatment of how the ULX component is scaled relative to the stellar population as a function of metallicity and burst age.As is evident from Figure 15, the "DIS" model has a stronger soft X-ray excess than our SED ULX , leading to a different shape when extended into the EUV, particularly around 54 eV where He II λ4686 has its ionization potential.While our SXP ULX simulations are run coupled to a range of SSPs with L X /M 0 scaled as a function of burst age and metallicity, the simulations in Simmonds et al. (2021) are performed by scaling the ULX component to a single BPASS SSP (1 Myr, ∼ 0.3 Z ) using a range of L X /SFR.This amounts to relaxing the assumptions used here for scaling of the SXP ULX formation efficiency with SSP metallicity and age, and the requirement for a delay time for ULX formation.In this way, the approach in Simmonds et al. (2021) complements the average efficiencies employed in the SXP ULX simulations: using a broad range of L X /SFR for a given Z and t burst is akin to simulating a large degree of scatter in the X-ray production efficiency for that SSP.
While the simulations in Senchyna et al. (2020), Simmonds et al. (2021), and this work explore HMXBs and ULXs specifically, there are other recent approaches to the high-energy ionizing photon production question that take a more agnostic approach to source type.In Umeda et al. (2022), the authors attempt to parametrically constrain the EUV shape of the intrinsic ionizing spectra for a sample of extremely metal-poor galaxies.To do so, they construct SEDs consisting of parametrized powerlaw and blackbody components which they run through CLOUDY.They then use MCMC to determine the parameters of these SED components that best describe a set of observed high-ionization emission lines from the sample of extremely metal-poor galaxies.Their best-performing models for each galaxy are then compared with more physical SEDs, in this case BPASS SSPs coupled to different ionizing sources such as ULX or AGN, to infer physical properties of the ionizing population.In Figure 15, we show the results of this approach for one of their galaxies, I Zw18 (blue dashed line), for which there exist excellent broad-band multi-wavelength coverage (i.e., X-ray, optical, and high-ionization nebular emission lines).We note that the results are remarkably consistent with our SED ULX component, despite very different approaches in (re)constructing the intrinsic ionizing spectrum.12In general, using their full galaxy sample, Umeda et al. (2022) show that the intensity of the observed highionization emission lines requires different combinations of BPASS SSPs (in terms of burst age) and ULX-like components with various scaling factors relative to the SSPs.Our [SXP ULX + SSP](t burst , Z) grid is constructed with these various scaling factors determined from theoretical binary population synthesis models, while in Umeda et al. (2022) the characteristic burst ages and necessary scaling factors for the high-energy ionizing photon production are inferred from the best parametric fits after comparison with physical models.The two approaches are therefore quite complementary, and highlight that determining the sources that contribute to high-energy ionizing photon production may require consideration of not only different sources of ionizing photons, but also their relative contributions as a function of stellar population age and metallicity.
As a final point of comparison, we highlight the results from Olivier et al. (2022), where the authors fit the highionization emission lines in a set of nearby EELGs with very young ( 10 Myr) and metal-poor (∼ 0.1 Z ) stellar populations.These authors find that BPASS SSPs alone cannot reproduce the observed strengths of some of the highest ionization potential species, such as He II and [O IV].To reproduce the strengths of such lines, they opt to add an 80,000 K blackbody to the BPASS SSPs scaled to a relatively large fractional contribution (45-55%) to the total luminosity.We show such a blackbody component in Figure 15 as the pink dashdot line ("BB").As already discussed in Section 6.1, the SED ULX includes a pseudo-blackbody component, albeit at slightly different temperature and normalization than the blackbody employed in Olivier et al. (2022).However, the temperature of the roughly blackbody component of the SED ULX and its relative contribution to the total composite SED depends on the specific assumptions about the accretion flow geometry, accretor type and mass, and mass supply rate, and our SED ULX model is certainly not exhaustive in these respects.Thus, an EUV blackbody source with no corresponding X-ray emission as employed in Olivier et al. (2022) is not necessarily incompatible with an SXP origin (e.g., Figure 3).
These results highlight a few crucial considerations for simulating photoionization due to an SXP-like source.Namely, the results depend heavily on the assumed SED shape for this additional ionizing component, and on how the component is scaled relative to the stellar population.Unfortunately, we still do not have tight empirical constraints on either of these ingredients for the case of ULXs specifically, or XRBs more broadly.However, as these photoionization simulation results suggest, select high-ionization emission lines may offer an additional lever arm for inferring the shape of unseen portions of the SED.The approaches presented in Umeda et al. (2022) and Olivier et al. (2022), for example, have the potential to be very powerful, but the recovered bestperforming models depend critically on the suite of lines available for fitting as well as the set of parametric models considered.Whether using parametric or physicallymotivated models as inputs to the photoionization simulations, larger grids are likely necessary going forward.We consider the recent approaches highlighted here important first steps in this process.

Alternative Ionizing Sources
For the SXP ULX model used in this work, we have assumed the accretor is a stellar mass BH (∼ 10-100 M ) and have constructed the SED ULX accordingly.In this section, we discuss the effects of changing these assumptions about the accretor (both mass and type), and additional processes capable of producing high-energy ionizing photons.

Neutron Star ULXs
While the simulations presented here consider only the case of stellar mass BH accretor for ULXs, there is now direct and growing observational evidence for NS accretors in ULXs, via the detection of pulsations (e.g.Bachetti et al. 2014;Fürst et al. 2016;Israel et al. 2017;Carpano et al. 2018).Swapping the BH for a NS accretor in the SXP ULX would result in a few fundamental differences, both in terms of longevity of the ionizing output relative to the stellar population and potentially the SED shape.
For the SSPs considered in this work (t burst ≤ 20 Myr) BHs should dominate the accretor demographics, while NS likely become the main accretor population at later times (e.g., Wiktorowicz et al. 2019).The addition of NS ULXs to the SXP ULX would therefore further prolong the ionizing output of a population, beyond that of the most massive stars and binaries, and the BH SXP ULX presented here.Furthermore, because the stellar ionizing contribution begins to wane considerably on timescales > 20 Myr after the massive star population is depleted, any additional source of high-energy photons present on these timescales will provide a larger fractional contribution to the overall ionizing photon budget.However, the contribution of different accretor demographics to the total L X /M 0 will be both SFHand metallicity-dependent, and such a population breakdown for the scaling relations is not yet well-constrained theoretically or observationally.
In addition to timescale differences, the spectral shape and upper limit for the luminosity may well be different for NS accretors relative to stellar mass BHs in the supercritical regime.This is primarily due to the fact, in the case of a NS, the Eddington limit and accretion flow geometry may be modified by the presence of a strong magnetic field.In calculating the Eddington limit, the electron scattering cross section is typically given as the Thomson scattering cross section; however, in the presence of a strong magnetic field, this cross section may be reduced, effectively changing the radiation pressure term in Equation 1, and increasing the Eddington limit for a given mass (Paczynski 1992).Indeed, there is recent evidence for cyclotron resonance scattering features in the spectra of both a ULX and a hyper-luminous Xray source (L X > 10 41 erg s −1 ), which implies the presence of an accreting NS with very strong magnetic field (Brightman et al. 2018(Brightman et al. , 2022)).It is important to note, however, that an increase in Eddington limit due to a NS with a strong magnetic field is not mutually exclusive with a strongly super-Eddington mass supply rate (e.g., Bachetti et al. 2022), which may also lead to driving an outflow (e.g., King & Lasota 2016;King et al. 2017).
A strong dipole field also funnels accreted material onto the poles of the NS (i.e., creates an accretion column) and possibly further modifies the picture of the supercritical accretion flow presented in Section 2.3.2,depending on how the extent of the Alfvén radius (R M ) compares to the spherization radius (R sph ).King et al. (2023) present a very detailed review and discussion of the current landscape of models for supercritically accreting NS, both pulsing and non-pulsing.Despite this attention to modeling, there exists no, to our knowledge, energetically self-consistent model in e.g., XSPEC for a NS ULX that allows for modeling the effects of a strong magnetic field as well as inclination-dependence on the resultant SED.We consider further investigation of the NS ULX contribution to the SXP ULX a promising avenue for future study, but outside the scope of the present work.

Intermediate Mass Black Holes
Moving up in the mass scale, IMBHs (100 M ≤ m ≤ 10 6 M ) have been proposed as an additional source of ionizing photons relative to stellar populations (e.g., Hatano et al. 2023).As relatively massive BHs with moderate to low mass transfer rates relative to Eddington, IMBHs may also have SEDs with a substantial EUV component.In Figure 16, we compare the SED ULX with IMBH SEDs for a 10 3 M BH and 10 5 M BH based on the qsosed model (Kubota & Done 2018), all normalized to the same L X (0.5-8 keV) for plotting purposes.As Figure 16 shows, it can be difficult to distinguish between a ULX and IMBH for sources of similar bolometric luminosity that are sub-dominant with respect to the stellar population, particularly in the absence of broadband X-ray coverage.However, differences in the SED shapes in the EUV imply there is diagnostic power with high-ionization lines for discerning ionization by accreting BHs of different masses.We provide a preliminary exploration of this via a comparison with select results from the IMBH photoionization simulations presented in Richardson et al. (2022).
The full set of IMBH photoionization simulations from Richardson et al. (2022) is expansive, encompassing dif-   15, but for a comparison with 10 3 M and 10 5 M IMBH SEDs (orange dotted and green dashed lines, respectively), as employed in the photoionization simulations from Richardson et al. (2022).All SEDs are normalized to the same LX (0.5-8 keV) for plotting purposes.In the absence of broad-band X-ray coverage, it can be difficult to distinguish between the SXPULX and IMBH SEDs, particularly the sources are similar in bolometric luminosity and sub-dominant relative to the stellar population.However, differences in the intrinsic SEDs in the EUV-tosoft X-ray regime may result in some diagnostic power for distinguishing between these different ionizing sources using high-ionization nebular emission lines.Ionization potentials for a selection of lines are labeled and marked by grey vertical lines.
ferent BH masses, prescriptions for the IMBH SED, mixing methodologies between the IMBH and stellar population, and a relatively fine gridding in stellar metallicity and log U. We compare only with the subset of their results most closely matching our simulation inputs and set-up.Namely, we show only their simulations for a 10 3 M IMBH (qsosed model) combined with a 20 Myr BPASS SSP assuming coincident mixing and a closed cloud geometry for grid points with -2 ≤ log(Z) ≤ 0 and -4 ≤ log U ≤ 1.In Figures 9 and 12 we plot their simulation results for IMBHs with fractional contributions ≤ 16% relative to the 20 Myr SSP (bottom right panels).
The comparisons in Figures 9 and 12 illustrate that there is strong diagnostic potential with very highionization line species for distinguishing between the ionizing impact of BHs across orders of magnitude in mass (e.g., Cann et al. 2018;Richardson et al. 2022).Offsets between the IMBH grids and grid with SXP ULX contribution in these emission line diagnostics are driven primarily by SED shape, suggesting these particular line ratios are sensitive to BH mass.However, this implies the diagnostic power hinges on assumed SED shape for a given BH mass, as well as typical fractional contribution of the BH ionizing component relative to the stellar population.
To more clearly show the differences due to changing IMBH ionizing fraction, in Figure 17 we instead plot different fractional contributions of the IMBH (from 4-100%) relative to the SSP in each panel, as compared only with the 20 Myr [SXP ULX + SSP](t burst , Z, log U) grid from this work.As before, the black solid line represents the empirically derived star formation/AGN demarcation, which the grid with SXP ULX contribution never exceeds, and the IMBH grids only exceed for ionizing contribution ≥ 16%.In addition to the fractional contribution relative to the stellar population, differences between the grids in this diagnostic space also depend critically on how the additional ionizing component is scaled with respect to the metallicity of the stellar population.For example, the IMBH simulations that we compare to do not include a metallicity-dependence for the scaling of the IMBH component relative to the SSP.In this way, the IMBH grids are produced in a fashion similar to our "no metallicity-dependence" model (e.g., Figure A.2), and can therefore achieve more extreme line ratios across a range in metallicities from solar to sub-solar.If instead the typical IMBH mass scales with stellar population metallicity, some regions of the IMBH grids as presented here would potentially be excluded.Harnessing the full power of such emission line diagnostics for investigating BH ionizing impact across the mass scale will therefore likely require careful consideration of typical host galaxy environments for ULXs and IMBHs, and how this affects their scaling and mixing with stellar populations (e.g., Polimera et al. 2022;Richardson et al. 2022).

Shock Ionization
While the grid with SXP ULX contribution presented here is capable of reproducing some of the range of observed line ratios for high-ionization emission lines, the models fail elsewhere.This implies the need either for tweaks to the SXP ULX implementation, or consideration of some alternative energetic process.
Here we explore the ionizing impact of fast radiative shocks, which have been considered by many others as a potential source of the requisite high-energy ionizing photons to explain the observed line intensities in select EELGs (e.g., Thuan & Izotov 2005;Izotov et al. 2012;Plat et al. 2019;Izotov et al. 2021).However, instead of suggesting fast shocks as an alternative source of ionizing photons relative to the SXP ULX , we consider them as an additional source related to the SXP ULX itself.
The accretion flow model prescribed for the SXP ULX in this work includes a quasi-spherical outflow component, which will likely interact with the surrounding medium.The interaction of ULX outflows with the interstellar medium is fairly well motivated by observations.A number of nearby ( 10 Mpc) ULXs have been observed to sit inside large (> 100 pc) bubbles.Such nebulae may be radiatively or mechanically inflated (e.g., Cseh et al. 2012;King et al. 2023).Several very bright ULXs show more direct evidence for the presence of fast outflows via the detection of blueshifted absorption lines in high-resolution X-ray spectra (Poutanen et al. 2007;Pinto et al. 2016;Kosec et al. 2018a,b).Very recently, a detailed optical spectroscopic study of a ULX with direct evidence for fast outflows, as well as a bubble detection, confirmed the ULX wind or jet likely drives bubble expansion (Gúrpide et al. 2022).And indeed, investigation of the emission line diagnostics for such ULX bubbles reveal signatures of shock and photoionization (e.g., Abolmasov et al. 2007;López et al. 2019;Gúrpide et al. 2022).Therefore, in at least some cases, ULX feedback is sufficient to inflate a bubble and drive both shock and photoionization.
The SED ULX model employed here, by virtue of including the outflow component, is therefore capable of driving shock ionization as well as photoionization, though we model only the latter in this work.Detailed shock and precursor ionization model grids are available from MAPPINGS (Allen et al. 2008), which may in principle be matched to a corresponding SXP ULX to simulate the combined signature of shock and photoionization.However, including the effects of shocks self-consistently alongside SXP ULX photoionization, assuming the shock is driven by the SXP ULX itself, will likely require more than the addition of separate grids.In all simulations presented here we assume the gas is uniformly distributed in the cloud, but the presence of shocks related to the SXP ULX (or indeed the stellar population) could preferentially displace gas from certain regions of the cloud.This may lead to shock ionization occurring in regions with different density than where the primary photoionization is occurring (e.g., Izotov et al. 2021).Such a non-uniform density may also alter the intensities of typical strong lines such as [O III], which is produced in the inner regions of the cloud.If gas is displaced from the central regions, this would amount to effectively lowering the ionization parameter for species such as [O III], relative to species produced elsewhere in the cloud (e.g., Kewley et al. 2013).It is clear a more careful treatment of joint consideration of SXP ULX shock and photoionization is warranted in terms of timescales, cloud properties, and mixing methodologies to understand the overall effects on typical line diagnostics.Given the modeling complexities, we defer consideration of joint SXP ULX shock and photoionization to a future related work.
As a final note, we remark that the presence of relatively ubiquitous ULX outflows implies the potential for feedback that can strongly shape the interstellar medium, akin to feedback from stellar winds and supernovae.For example, a 10 7 M instantaneous burst of star formation yields ∼ 10 40 erg s −1 mechanical luminosity over the course of 10 Myr as the most massive stars evolve and explode as supernovae (e.g., Leitherer et al. 1999;Prestwich et al. 2015).From Figure 1, a 0.1 Z SXP ULX provides ∼ 10 39 -10 40 erg s −1 mechanical luminosity for the same instantaneous burst, assuming ULX mechanical luminosity on par with radiative output.For at least some observed ULXs, it appears mechanical power from the outflow actually exceeds radiative power, suggesting the above estimate may even be a lower limit for SXP ULX feedback (e.g.Justham & Schawinski 2012;Pinto et al. 2016;Gúrpide et al. 2022).Feedback from an SXP ULX could therefore contribute to favorable conditions for ionizing photon escape (e.g., similar to the scenario for stellar populations from Jaskot & Oey 2013).Indeed, resolved studies of Lyman continuum leakers at X-ray wavelengths suggest ULXs could be an important source of mechanical feedback, facilitating the escape of ionizing photons by carving channels in the surrounding medium (Prestwich et al. 2015;Kaaret et al. 2017).An SXP ULX with strong outflow component could therefore be important for (1) producing high-energy ionizing photons via both shock and photoionziation; (2) facilitating ionizing photon escape; and (3) producing the soft (< 2 keV) photons that can more efficiently heat the intergalactic medium prior to the epoch of reionization (e.g., Das et al. 2017;HERA Collaboration et al. 2023).As such, continued attention to modeling SXP ULX -like sources is extremely relevant to high redshift studies.

SUMMARY
In this work, we have presented a framework for combining a "simple X-ray population" alongside the corresponding simple stellar population in a physically consistent and meaningful manner in order to simulate and explore their combined ionizing impact as a function of (1) instantaneous burst age and (2) metallicity.Using our combined SED as input to the photoionization code CLOUDY, we have produced a nebular line and continuum emission grid, which we make publicly available.Our principal findings can be summarized as follows: • The addition of the SXP ULX prolongs the ionizing output relative to single massive stars or the products of binary evolution already included in BPASS.This is a consequence of the physically robust combination of an SXP ULX to the corresponding BPASS SSP with relative scaling factors determined as a function of burst age and stellar metallicity.Namely, the SXP ULX is scaled with an increasing fractional contribution as a function of decreasing stellar metallicity, in line with empirical constraints, and combined only with an SSP of the corresponding age, following the expectation that the SXP ULX evolves from the parent stellar population.
• For the burst timescales modeled here (t burst ≤ 20 Myr) the ionizing photons due to the SXP ULX contribute no more ∼ 5% to the total nebular line and continuum emission.Broad-band FIRto-UV colors are therefore not strongly affected by the addition of this component, relative to the case of a stellar-only ionizing component.
• Given that the SXP ULX has an SED which is very flat through the EUV, the addition of this component increases the intensity lines with ionization potentials 54 eV such as He II λ1640,4686 and [O IV] 25.9µm by a factor of at least two relative to the BPASS SSPs alone.Emission line diagnostics including such lines can therefore be used to infer the presence of SXP ULX ionization.
• The intensity of very high-ionization lines (ionization potentials > 90 eV) such as [Ne V] λ3426, 14.3µm are strongly enhanced (> 10 5 ×) due to the addition of the SXP ULX ; however, these lines are typically very weak in the simulations presented here (e.g., log( f line / f Paβ ) ∼ -9 for coronal lines in the IR).This is due in large part to the fixed SXP ULX normalization relative to the SSP as a function of burst age and metallicity, where the maximum relative scaling is achieved for the lowest stellar metallicities.Because gas-phase and stellar metallicities are coupled in all simulations, the regime in which the SXP ULX reaches maximum ionizing contribution relative to the SSP is also the one where many elements have their lowest absolute abundances in the cloud, which can be the dominant factor in setting line intensity.
• SXP ULX photoionization is capable of reproducing the observed strengths of high-ionization emission line ratios (e.g., He II λ4686/Hβ), but primarily for instantaneous burst models with Z 0.1 Z and t burst > 10 Myr.This is again a consequence of the coupling of the SXP ULX to the SSP as a fixed function of stellar metallicity and burst age.On timescales < 10 Myr, the stellar ionizing continuum is still substantial, effectively diluting the SXP ULX ionizing output, particularly for highionization lines as measured relative to strong nebular lines for which there is little SXP ULX ionizing contribution.In this way, the assumed star formation history strongly affects the efficacy of the SXP ULX to high-energy ionizing photon production in terms of typical observed line ratios.
• Comparison of the simulation outputs to observations in terms of emission line ratios and/or typical X-ray observables (e.g., L X /SFR) should be performed with careful consideration of the underlying assumptions, including SXP ULX SED shape, cloud abundance patterns, assumed star formation history, and SXP ULX detectability as a function of viewing angle and variability.
• Though the X-ray detectability of the SXP ULX depends on viewing angle and duty cycle, the SXP ULX production of high-ionization nebular emission lines is ubiquitous for the SED ULX modeled here.This is due to the fact that the EUV and soft X-ray photons from the SXP ULX emanate from the outflow component, which is assumed to be relatively isotropic.In this way, high-ionization emission lines can be indirect tracers of a transient or unfavorably oriented SXP ULX .
• The outflow component of the SXP ULX may additionally drive fast shocks and therefore shock ionization.We discuss complexities in joint modeling of shock and photoionization due to an SXP ULX , and consider it a promising avenue for future study.This is particularly true given that SXP ULX mechanical feedback due to a disk wind or outflow could be similar in magntiude, if not timescale, to feedback from stellar winds and supernovae.As such, the SXP ULX could contribute to conditions in the interstellar medium conducive to ionizing photon escape, and additionally provide a substantial number of soft (0.5-2 keV) photons that could directly contribute to heating the intergalactic medium prior to the epoch of reionization.Note-Criteria for strength relative to Hβ or Paβ, or enhancement relative to the SSP-only models are described in Section 5.1.Only lines with ionization potentials ≥ 54 eV are marked in the table according to the above criteria.
Figure1.The adopted scaling of LX/M0 for the SXP as a function of instantaneous burst age and metallicity (Z : dark blue dashed line, 0.1 Z : purple solid line, and 0.01 Z : orange dash-dot line) from theoretical binary population synthesis models(Fragos et al. 2013b).The instantaneous burst ages at which we model the SXP contribution in the photoionization simulations are marked with dark blue squares, purple circles, and orange diamonds for the Z , 0.1 Z , and 0.01 Z cases, respectively.The SXP component in our models turns on for instantaneous bursts > 3 Myr.Empirical results for LX/M0 from star-forming galaxy stacks in the Chandra Deep Fields, for which the median metallicity (∼ 0.6 Z ) is intermediate to the theoretical scalings adopted here, are shown as a green dotted line for reference.

Figure 2 .
Figure 2. Cartoon depiction of the accretion flow geometry for the SEDULX model.Key components and parameters are labeled, including the accretor (BH) and accretion disk, evacuated wind funnel and associated funnel opening angle (θ), the inner and outer photospheric radii (R ph,in , R ph,out , respectively), the spherization radius (R sph ), and the inclination angle (i) with respect to an observer.The wind or outflow, which is launched from the accretion disk and assumed to be a quasi-spherical component, is shown by the blue arrows.Irradiating photons emanate from the hot inner accretion flow and funnel walls, as labeled.

Figure 4 .
Figure4.The composite SXPULX + SSP SEDs (solid lines) for the sampled metallicities (dark blue: Z , purple: 0.1 Z , and orange: 0.01 Z ) normalized to 1 M stellar mass formed in an instantaneous burst, where the burst age is annotated in each panel.In each panel, the dash-dot line shows the SXPULX component of the composite, while the dashed line denotes the SSP contribution.The dark blue shaded region indicates the 0.5-12 keV range, while the light blue shaded region shows the EUV regime.Grey dotted vertical lines mark the ionization potentials for select lines, from high to low energies (left to right):[Ne V], [Ne IV], [O IV], C IV, He II, [Ar II], and H.The addition of the SXPULX to the corresponding SSP substantially adds to the intensity at energies 54 eV (log(Wavelength Å) = 2.36), particularly for the low-metallicity models with t burst > 10 Myr.

Figure 5 .
Figure5.The composite SXPULX + SSP SEDs at 0.1 Z for a selection of instantaneous burst ages, normalized by a range of log U values, as noted in the colorbar.Models with high log U have the highest intensities, and correspond to the largest stellar mass formed in a burst.As in Figure4, the darker blue shaded region shows the X-ray bandpass, the lighter blue shaded region encompasses the EUV regime, and the grey vertical lines denote a selection of relevant ionization potentials (13.6-126 eV, from right to left).Each composite SED is used as input for the CLOUDY photoionization simulations.

%
Figure6.Top: The fractional contribution of the nebular line and continuum emission to the total flux (stellar and nebular) as a function of wavelength and instantaneous burst age for models with Z = 0.1 Z and log U = −1.Bottom: The percentage enhancement in nebular line and continuum flux due to the addition of the SXPULX.For all burst timescales simulated, the SXPULX-driven enhancement in the total nebular emission is typically small (< 5%).The black labels denote high-ionization species whose nebular line intensity is particularly enhanced due to the inclusion of the SXPULX component.

Figure 7 .
Figure 7. Potential UV diagnostic diagram for SXPULX ionization, as a function of instantaneous burst age.The [SXPULX + SSP](t burst , Z, log U) grid is shown as solid lines, where orange to dark blue lines represent metallicities 0.01 Z -Z and light blue to dark blue lines represent log U = −4 to −1, as illustrated by the colorbars.In all panels, the star represents the 0.01 Z and log U = −1 grid point.The dashed transparent lines in the background show the corresponding grid with SSP ionizing contribution only.A sample of local BCDs for which these UV lines have been measured are shown as cyan diamonds, where the uncertainties in line ratios are smaller than the plotted points(Senchyna et al. 2017).The addition of the SXPULX hardens the ionizing spectrum at a given log U, thereby increasing the intensity of He II λ1640 relative to C III] λ1907,9.This effect is more pronounced at solar metallicity for t burst = 5 Myr, following depletion of WR stars, and low metallicities with t burst > 10 Myr (bottom panels), following depletion of the most massive stars.

Figure 8 .
Figure8.Same as Figure7, but for a potential optical diagnostic diagram accounting for the inclusion of SXPULX ionization.The data points are a sample of star-forming galaxies fromShirazi & Brinchmann (2012) with strong nebular He II λ4686 emission (red triangles: galaxies without WR star features, blue circles: galaxies with WR star features).The black solid line is the empirically derived relation fromShirazi & Brinchmann (2012) to separate ionization due to star formation (below the line) versus ionization due to active galactic nuclei or composite populations (above the line).The models with SXPULX contribution overlap with the observed intensities for some of the Non-WR feature galaxies, particularly for t burst > 10 Myr and Z ≤ 0.1 Z ; however, neither the [SXPULX + SSP](t burst , Z, log U) nor the SSPonly grids are capable of reproducing some of the more extreme galaxies, particularly those with detected WR features and higher N II/Hα near the empirical maximum starburst line.
Figure 9 for lines accessible with JWST (e.g., Weaver et al. 2010).Here the [O IV] 25.9µm/[Ne III] 15.66µm ratio is sensitive to hardness of the ionizing spectrum, while [Ne III] 15.66µm/[Ne II] 12.81µm is primarily sensitive to the ionization parameter.For this diagnostic, we include the classification regions from Richardson et al.

Figure 9 .
Figure 9. Same as Figure 7, but for a potential IR diagnostic diagram for SXPULX ionization.The dot-dash grids are photoionization simulations of a 10 3 M IMBH at different fractional contributions relative to a 20 Myr SSP from Richardson et al. (2022), where the color scheme is the same as for the [SXPULX + SSP](t burst , Z, log U) grid for values log U and Z in common between the grids.The dashed lines show potential classification regions for pure star-formation ionization, composite star-formation and AGN ionization, and pure AGN ionization, also from Richardson et al. (2022), as annotated in the center top panel.

Figure 10 .
Figure10.A UV emission line diagnostic where the SXPULX contribution does not change the high-ionization line species appreciably relative to the SSP-only grid.The cyan diamonds are again a sample of local BCDs for which these UV lines have been measured(Senchyna et al. 2017).The lack of overlap between some of the galaxies and either grid for this diagnostic is likely a function of adopted C/O ratio for the photoionization modeling in this work.

Figure 11 .Figure 12 .Figure 13 .
Figure11.Same as Figure7, but for a potential UV/optical diagnostic diagram for SXPULX contribution.For these highionization lines, the SXPULX contribution is substantial relative to the SSP-only case, but the high-ionization line intensities are generally weak except for select grid points (low metallicities and high log U).The dashed lines illustrate classification regions for pure star-formation ionization, ionization from composite populations, ionization due to AGN, and potential ionization due to Pop III stars fromCleri et al. (2023).The cyan squares are the observed line ratios from a sample of local EELGs(Izotov et al. 2012;Berg et al. 2021;Olivier et al. 2022) that are consistent with select [SXPULX + SSP](t burst , Z, log U) models in a very narrow range of grid space.

Figure 14 .
Figure14.The emergent LX (0.5-8 keV)/SFR relative to the strength of He II λ4686/Hβ from photoionization simulations with SXPULX contribution for both continuous star formation (circles) and instantaneous bursts (shaded regions).As before, metallicities are color-coded and labeled.Observed values for a sample of star-forming dwarfs fromSenchyna et al. (2020) are shown as green diamonds.The low metallicity simulations are consistent with all but the most extreme observed values and upper limits for the line ratios and LX/SFR; however, this requires instantaneous bursts where t burst > 10 Myr.Models with continuous SFHs produce LX/SFR in the range of the observed galaxies and upper limits, but much lower ratios of He II λ4686/Hβ.

Figure 16 .
Figure 16.Similar to Figure15, but for a comparison with 10 3 M and 10 5 M IMBH SEDs (orange dotted and green dashed lines, respectively), as employed in the photoionization simulations fromRichardson et al. (2022).All SEDs are normalized to the same LX (0.5-8 keV) for plotting purposes.In the absence of broad-band X-ray coverage, it can be difficult to distinguish between the SXPULX and IMBH SEDs, particularly the sources are similar in bolometric luminosity and sub-dominant relative to the stellar population.However, differences in the intrinsic SEDs in the EUV-tosoft X-ray regime may result in some diagnostic power for distinguishing between these different ionizing sources using high-ionization nebular emission lines.Ionization potentials for a selection of lines are labeled and marked by grey vertical lines.

Figure 17 .
Figure 17.Similar to Figure 8, but with each panel showing only the 20 Myr [SXPULX + SSP](t burst , Z, log U) grid (solid lines), relative to grids covering similar log U and Z with different fractional IMBH contributions (dashed lines) from Richardson et al. (2022), where the IMBH contribution is annotated in the lower right of each panel.Emission line diagnostics including high-ionization species such as He II λ4686 can be useful for distinguishing between ionization due to SF versus accreting black holes across the mass scale (i.e., SXPULX, IMBH, AGN); however, classification regions depend on the fractional contribution of the additional ionizing source (e.g., SXPULX or IMBH) relative to the stellar population, and whether this contribution scales strongly as a function of metallicity and burst age.
609765[C I] 609.765µmC 1 609.590m* Line that is strong relative to Hβ or Paβ, and strongly enhanced by the addition of the SXPULX, relative to SSP-only models.† Line that is strong relative to Hβ or Paβ, but not strongly enhanced relative to SSP-only models.‡ Line that is weak relative to Hβ or Paβ, but strongly enhanced relative to SSP-only models.|| Line with zero flux.

Table 1 .
Summary of calculated and selected parameters for the SEDULX model, as implemented via the simpl(sirf) templates in XSPEC.Parameter names follow nomenclature used in XSPEC, as described in the text.