# Modeling the High-Energy Ionizing Output from Simple Stellar and X-ray
Binary Populations

Kristen Garofali NASA Goddard Space Flight Center, Code 662, Greenbelt, MD
20771, USA  Antara R. Basu-Zych NASA Goddard Space Flight Center, Code 662,
Greenbelt, MD 20771, USA  Center for Space Science and Technology, University
of Maryland, Baltimore County, Baltimore, MD 21250, USA  Center for Research
and Exploration in Space Science and Technology, NASA Goddard Space Flight
Center, Greenbelt, MD 20771, USA  Benjamin D. Johnson  Center for Astrophysics
$|$ Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA
Panayiotis Tzanavaris  Center for Space Science and Technology, University of
Maryland, Baltimore County, Baltimore, MD 21250, USA  Center for Research and
Exploration in Space Science and Technology, NASA Goddard Space Flight Center,
Greenbelt, MD 20771, USA  American Physical Society, Hauppauge, New York, NY
11788, USA  Anne Jaskot Department of Astronomy, Williams College,
Williamstown, MA 01267, USA  Chris Richardson Elon University, Elon, NC 27278,
USA  Bret D. Lehmer Department of Physics, University of Arkansas,
Fayetteville, AR 72701, USA  Mihoko Yukita William H. Miller III Department of
Physics and Astronomy, Johns Hopkins University, Baltimore, MD 21218, USA
Edmund Hodges-Kluck NASA Goddard Space Flight Center, Code 662, Greenbelt, MD
20771, USA  Ann Hornschemeier NASA Goddard Space Flight Center, Code 662,
Greenbelt, MD 20771, USA  Andrew Ptak NASA Goddard Space Flight Center, Code
662, Greenbelt, MD 20771, USA  Neven Vulic Eureka Scientific, Inc., Oakland,
CA 94602-3017, USA

###### Abstract

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 ultra-luminous X-ray sources (ULXs) 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 $\lambda$1640,4686, [Ne V]
$\lambda$3426,14.3$\mu$m, and [O IV] 25.9$\mu$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 ($\sim$ 0.1 $Z_{\odot}$), 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 discuss how the ionization
signatures compare to models for sources such as intermediate mass black
holes.

††software:  CLOUDY (Ferland et al., 2017), cloudyFSPS (Byler, 2018), FSPS
(Conroy et al., 2009; Conroy & Gunn, 2010), python-fsps (Johnson et al.,
2022), XSPEC (Arnaud, 1996), HEASoft (Nasa High Energy Astrophysics Science
Archive Research Center (2014), Heasarc), astropy (Astropy Collaboration et
al., 2013, 2018, 2022), scipy (Virtanen et al., 2020), matplotlib (Hunter,
2007), astroquery (Ginsburg et al., 2019), numpy (Harris et al., 2020)

##  1 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 ($\lesssim$ 30 Myr) star formation, and relatively
low metallicities ($\sim$ 0.1 $Z_{\odot}$) (e.g., Jaskot & Oey, 2013; Stark,
2016; Berg et al., 2018, 2019; Senchyna et al., 2017, 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.0 keV X-ray luminosity per galaxy star formation rate ($L_{\rm
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 $\gtrsim$ 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., 2012, 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 high-ionization 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 population 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 ultra-luminous 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 ($\gtrsim$ 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 grid-based 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.

##  2 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 following 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.

###  2.1 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 $\lesssim$ 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_{\rm 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 time-averaged $L_{\rm 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_{\rm 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_{\rm X}$/$M_{0}$ is the preferred quantity since
transforming number counts to total radiative power is non-trivial, 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. (2008, 2009), the bolometric
$L_{\rm 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.

Figure 1: The adopted scaling of $L_{\rm X}$/$M_{0}$ for the SXP as a function
of instantaneous burst age and metallicity ($Z_{\odot}$: dark blue dashed
line, 0.1 $Z_{\odot}$: purple solid line, and 0.01 $Z_{\odot}$: 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_{\odot}$, 0.1 $Z_{\odot}$, and 0.01
$Z_{\odot}$ cases, respectively. The SXP component in our models turns on for
instantaneous bursts $>$ 3 Myr. Empirical results for $L_{\rm X}$/$M_{0}$ from
star-forming galaxy stacks in the Chandra Deep Fields, for which the median
metallicity ($\sim$ 0.6 $Z_{\odot}$) is intermediate to the theoretical
scalings adopted here, are shown as a green dotted line for reference.

The Fragos et al. (2013b) binary population synthesis models used here provide
published values of $L_{\rm X}$/$M_{0}$ for XRBs at three different stellar
metallicities111Throughout, we assume $Z$ = 0.02 corresponds to $Z_{\odot}$
for stellar metallicities. For abundances from Anders & Grevesse (1989), this
corresponds to a gas-phase metallicity 12 + log(O/H) = 8.93.: $Z$ = 0.002 (0.1
$Z_{\odot}$), $Z$ = 0.02 ($Z_{\odot}$), and $Z$ = 0.03 (1.5 $Z_{\odot}$), 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_{\rm burst}\sim$ {1, 3, 5, 8, 10, 15, 20}
Myr) are selected to maximize $L_{\rm 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_{\rm 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_{\odot}$) and $Z=0.02$ ($Z_{\odot}$), where the 0.1 $Z_{\odot}$
scaling for $L_{\rm X}$/$M_{0}$ is relevant for nearby EELGs and sources at
$z\gtrsim 3$ (e.g., Shirazi & Brinchmann, 2012; Madau & Dickinson, 2014). We
additionally add a third metallicity set at $\sim$ 0.01 $Z_{\odot}$, 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_{\rm X}$/$M_{0}$ on post-
starburst timescale as the 0.1 $Z_{\odot}$ model, but allow the absolute
normalization of $L_{\rm X}$/$M_{0}$ to be a factor of 3$\times$ higher than
for the 0.1 $Z_{\odot}$ at each time step. This scaling corresponds roughly to
the observed $\langle$$L_{\rm X}$/SFR$\rangle$ from galaxies with 12 +
log(O/H) $\sim$ 7.0, the lowest metallicity for which the observed $L_{\rm
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_{\rm 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_{\rm
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
$\sim$ 0.6 $Z_{\odot}$, 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_{\rm 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
($\lesssim$ 0.1 $M_{\odot}$ yr-1), where observed values for $L_{\rm
X}$/$M_{\star}$ or $L_{\rm X}$/SFR can be subject to scatter of order $\sim$ 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.

###  2.2 Source Type: Ultra-luminous X-ray Sources

With the adopted scaling for $L_{\rm 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 ($\lesssim$ 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., 2019, 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_{\rm X}$ $>$
1039 erg s-1. These luminosities, under the assumption of isotropic emission,
exceed the Eddington limit for a 10 $M_{\odot}$ BH or 1.4 $M_{\odot}$ NS:

$L_{\rm Edd}=\frac{4\pi GMc}{\kappa_{T}}=1.5\times
10^{38}~{}m~{}\frac{1.7}{1+X}~{}{\rm ergs^{-1}}$ (1)

where $M$ is the mass of the compact object (such that $m\equiv M/M_{\odot}$),
$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 $\kappa_{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). Broad-band 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 $\gtrsim$ 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., 2018a, b, 2021), and potentially strongly magnetized
NS accretors (Walton et al., 2018; Brightman et al., 2022).

These findings suggest ULXs ($L_{\rm X}$ $\sim$ 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 rates222As discussed in King et al.
(2023), the phrase “super-Eddington accretion” can lead to ambiguity, as
stellar mass BHs and NS in ULXs do not actually accrete matter at highly
super-Eddington rates; instead, these authors suggest that for ULXs it would
be more appropriate to refer to the mass supply or transfer rate as the super-
Eddington quantity in question. We adhere to this convention throughout.. 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 ($\leq$ 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 $\geq$ 3 $M_{\odot}$ 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 “SXPULX ”, to indicate the source type modeled. This selection captures
a bounding case for HMXB333There is no explicit ULX population included in the
theoretical models from Fragos et al. (2013b), and we therefore use the HMXB
component, corresponding to the bright end of the XLF, for normalizing the
SXP; however, the theoretical models do allow for XRBs undergoing outbursts to
exceed Eddington. We also note that empirical values for $L_{\rm X}$/$M_{0}$
from galaxy-integrated studies, where the populations can be assumed to be
dominated by sources at the brightest end of the XLF, i.e., ULXs, appear
relatively consistent with the theoretical scaling relations ostensibly for
HMXBs employed here, as shown in Figure 1. 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).

###  2.3 Model for Ultra-luminous X-ray Source SED

The efficacy of the SXPULX 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 ($\gtrsim$ 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 ($\sim$ 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 ($\sim$ 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 SEDULX) 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 SEDULX, and motivate our selections. The SEDULX 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.

####  2.3.1 Accretor Type

As a starting point to constructing the SEDULX, 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_{\rm X}$/$M_{0}$ normalizations from Section 2.1 span
burst ages $\leq$ 20 Myr, we choose a stellar mass BH ($<$ 100 $M_{\odot}$) as
the accretor for the SEDULX. 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 $\sim$ 4-40 Myr during their ULX phase (Middleton
& King, 2017; Wiktorowicz et al., 2017, 2019).

Given that we likewise elect to model the population at three discrete stellar
metallicities ($Z_{\odot}$, 0.1$Z_{\odot}$, and $\sim$ 0.01$Z_{\odot}$), we
construct the SEDULX model as a function of approximate BH masses at each of
these metallicities: 10 $M_{\odot}$, 20 $M_{\odot}$, and 40 $M_{\odot}$ from
solar to sub-solar 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 SXPULX 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 SEDULX shape
(Section 2.3.3), they do not strongly affect the normalization of the SEDULX,
which is key to setting the simulated nebular line intensities. Instead, the
primary factor affecting SEDULX normalization, and therefore nebular line
intensities, is the adopted scaling relation for $L_{\rm X}$/$M_{0}$ as a
function of age and metallicity, as presented in Section 2.1.

####  2.3.2 Accretion Flow Geometry & Components

Under the selection of a stellar mass accretor with super-Eddington mass
supply rate for the SEDULX, 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_{\rm ph,~{}in}$).

  * •

The region between the wind launching zone and so-called the spherization
radius ($R_{\rm ph,~{}in}<R<R_{\rm sph}$), where the wind is optically thick.
The spherization radius ($R_{\rm 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_{\rm sph}$). The outer extent of the wind is given by the effective
photospheric radius of the outflow, $R_{\rm 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_{\rm 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).

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
($\theta$), the inner and outer photospheric radii ($R_{\rm ph,in}$, $R_{\rm
ph,out}$, respectively), the spherization radius ($R_{\rm 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.

In defining our SEDULX 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 re-emits radiation from another component—may be critically
important in the context of SXPULX 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 $\lambda$4686
and [Ne V] $\lambda$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 SEDULX 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 “self-irradiated 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 SEDULX is implemented as simpl(sirf) in XSPEC,
with relevant parameters listed in Table 1. Under this accretion flow model,
the mass supply rate relative to Eddington sets the physical parameters (i.e.,
characteristic radii and temperatures) that define the shape of the intrinsic
spectrum, while the viewing angle ($i$ in Figure 2) determines the shape of
the observed spectrum. 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.

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.  Component |  | Name | Value | Units
---|---|---|---|---
(1) |  | (2) | (3) | (4)
|  | $T_{\rm funnel,~{}in}$ | Eq. 8 | keV
|  | $r_{\rm ph,~{}in}$ | Eq. 7 | $R_{\rm sph}$
|  | $r_{\rm ph,~{}out}$ | Eq. 9 | $R_{\rm sph}$
|  | $\theta$ | 45 | deg
sirf |  | $i$ | 0 | deg
|  | $\alpha$ | 0 | …
|  | $\gamma$ | 4/3 | …
|  | $\dot{m}$ | Eq 5–6 | …
|  | irrad. | 4 | …
|  | norm. | Fragos et al. (2013b) | …
|  | $\Gamma$ | 2.5 | …
simpl |  | Frac. Sctr | 0.1 | …
|  | Up Sctr. Only | 0 | …

####  2.3.3 Parametrization of SED Shape

We aim to construct the SEDULX 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 ($\theta$) and compact object mass
($m$). For the SEDULX used in this work, we select $\theta=45^{\circ}$ and $m$
= 10, 20, or 40 $M_{\odot}$ for $Z_{\odot}$, 0.1 $Z_{\odot}$, and $\sim$ 0.01
$Z_{\odot}$, respectively. As shown by the black solid line in Figure 3, the
choice of $\theta=45^{\circ}$ results in a relatively flat SED in the EUV
regime even without invoking high mass supply rates (i.e., this selection
assumes the SXPULX is not universally strongly beamed; Equations 5–6). This
choice for $\theta$ is therefore important for capturing a bounding case for
production of high-energy ionizing photons by the SXPULX, 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 SEDULX 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:

$\begin{split}\dot{M}_{\rm Edd}=\frac{48\pi GM}{c\kappa_{T}}\simeq\\\ &3\times
10^{-8}~{}m~{}~{}M_{\odot}~{}{\rm yr}^{-1}\end{split}$ (2)

where $m$ is the mass of the accretor in units of $M_{\odot}$. It is often
convenient to define a dimensionless mass transfer rate $\dot{m}$, as follows:

$\dot{m}=\dot{M}/\dot{M}_{\rm Edd}$ (3)

A source in the ultra-luminous state would therefore be expected to have
$\dot{m}$ $\gg$ 1\. The intrinsic luminosity of the disk then depends
logarithmically on the dimensionless mass transfer rate (Poutanen et al.,
2007):

$L^{\rm int}_{\rm bol}=L_{\rm Edd}(1+0.6{\rm ln}(\dot{m}))$ (4)

Following King (2009), the observed, or apparent, luminosity may differ from
the intrinsic luminosity due to geometrical beaming, such that $L_{\rm obs}$ =
$L^{\rm int}_{\rm bol}/b$. Under this formalism, the beaming factor $b$ is
related to the funnel opening angle ($\theta$) as b = $\Omega/4\pi$, where
$\Omega=2\times 2\pi(1-{\rm cos}\theta)$ 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:

$b=P_{\rm obs}(\theta)=1-{\rm cos}\theta$ (5)

For large funnel opening angles ($\theta$ $\sim$ 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
predetermined parameter values. Specifically, we parameterize $b$ in terms of
the mass transfer rate as follows:

$b=\begin{cases}1&\dot{m}\leq 8.5\\\ (8.5/\dot{m})^{2}&\dot{m}>8.5\end{cases}$
(6)

By this definition, high mass transfer rates ($\dot{m}>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 SEDULX used in this work has $\theta=45^{\circ}$, and therefore $b=P_{\rm
obs}(\theta)\sim 0.3$ and $\dot{m}\sim 16$, using Equations 5–6. In Figure 3
we show the SEDULX 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.

Figure 3: The broad-band SED for the SXPULX (i.e., simpl(sirf) in XSPEC) under
differing assumptions about the funnel opening angle ($\theta$), the mass
transfer rate onto the BH ($\dot{m}$), 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_{\odot}$ 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 $\theta$ and $\dot{m}$ 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 ($\theta$), which changes the corresponding mass transfer rate,
while the solid blue SED demonstrates the effect of uncoupling $\theta$ and
$\dot{m}$ and allowing for a very high mass transfer rate with a moderate
funnel opening angle.

With $\theta$, $\dot{m}$, 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 SEDULX model (see Table 1). The
spherization radius where the outflow is launched is calculated as $R_{\rm
sph}$ = $\frac{18}{cos\theta}R_{\rm S}\dot{m}$, where $R_{\rm
S}=2\frac{GM}{c^{2}}$ is the Schwarzschild radius444This is proportional to
the dimensionless spherization radius from Poutanen et al. (2007), given as
$r_{\rm sph}=5/3\times\dot{m}$, but critically also depends on the geometry of
the funnel.. The minimum inner radius for the photosphere ($r_{\rm ph,in}$) in
units of this spherization radius ($R_{\rm sph}$) can be calculated as as the
innermost stable circular orbit ($R_{\rm isco}=6\frac{GM}{c^{2}}$) divided by
$\theta$:

$r_{\rm ph,in}=\frac{1}{6\dot{m}}$ (7)

The corresponding temperature at the inner edge of the wind funnel ($T_{\rm
funnel,~{}in}$), and the outer radius of the photopshere ($r_{\rm ph,out}$),
also in units of the spherization radius, can be calculated analytically
following Abolmasov et al. (2009):

$\begin{split}T_{\rm funnel,~{}in}&=0.038(\Omega^{\prime}{\rm
tan}^{2}\theta~{}r_{\rm ph,in})^{-1/4}\\\ &\times~{}(2.8+0.6{\rm
ln}~{}\dot{m_{3}})^{1/4}\dot{m_{3}}^{-1/2}M_{10}^{1/4}~{}keV\end{split}$ (8)

$r_{\rm
ph,out}\simeq\left[\frac{2\Omega^{\prime}\sqrt{\dot{m}}}{(1+\alpha)\sqrt{cos\theta}}\right]^{1/(1+\alpha)}$
(9)

where $\dot{m_{3}}=\dot{M}$ in units of 103 $\dot{M}_{\rm Edd}$,
$\Omega^{\prime}=-8\pi~{}\rm{ln}~{}\rm{sin}\theta$, $M_{10}$ is the mass of
the BH in units of 10 $M_{\odot}$, and $\alpha$ is the exponent describing the
velocity law for the wind. Our Equation 8 corresponds to Equation 35 from
Abolmasov et al. (2009), where we have assumed $T_{\rm funnel,~{}in}$ is equal
to the temperature at the bottom of the funnel ($T_{\rm bot}$) for the sirf
template.

The remainder of the parameters for the sirf component of the SEDULX model are
selected, as summarized in Table 1, where notation corresponds to symbols used
in XSPEC. In particular, the SEDULX 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), $\alpha=0$
(no acceleration for the wind), $\gamma=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 SEDULX model is
calculated following the theoretical scaling of $L_{\rm 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 SEDULX 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 $\Gamma=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., 2013,
2014). As such, we consider the SEDULX shape at energies $\gtrsim$ 2 keV,
which is set by the simpl component, critical for comparisons to apparent
observed $L_{\rm X}$, but not particularly physically informative in terms of
shape.

We make this SEDULX 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
SEDULX 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., 2014, 2016,
2019a, 2019b; Lau et al., 2019; Wiktorowicz et al., 2017, 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 SEDULX 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 SEDULX with the
spectrum for an appropriate donor star, and therefore further customization of
the SEDULX model presented here.

##  3 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 SXPULX can be specified using the
same parameters governing the evolution of the stellar population. Ideally,
the SXPULX 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 SXPULX 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_{\odot}$, a low-mass (0.1–0.5
$M_{\odot}$) IMF slope of $-1.30$, and a high-mass (0.5–100 $M_{\odot}$) 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_{\odot}$. By
convention in FSPS, each SSP is normalized by stellar mass.

Figure 4: The composite SXPULX \+ SSP SEDs (solid lines) for the sampled
metallicities (dark blue: $Z_{\odot}$, purple: 0.1 $Z_{\odot}$, and orange:
0.01 $Z_{\odot}$) normalized to 1 $M_{\odot}$ 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 $\gtrsim$ 54
eV (log(Wavelength Å) = 2.36), particularly for the low-metallicity models
with $t_{\rm burst}>$ 10 Myr.

We construct a grid of BPASS SSPs in age and metallicity corresponding to the
selected burst ages and metallicities555We consider the extremely metal poor
case speculative in terms of SXPULX normalization, and therefore refer to this
case as 0.01 $Z_{\odot}$ throughout for simplicity, noting the true stellar
metallicity from BPASS for this set of models is slightly lower (0.005
$Z_{\odot}$). at which we model the SXPULX ($t_{\rm burst}\sim$ {1, 3, 5, 8,
10, 15, 20} Myr, and $Z$ = {0.0001, 0.002, 0.02}). Because the SXPULX($t_{\rm
burst}$, $Z$) models we produce are already in terms of 1 $M_{\odot}$ stellar
mass, we simply add them to the corresponding SSP to produce the composite
[SXPULX \+ SSP]($t_{\rm 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 SXPULX 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 SXPULX (high $L_{\rm X}$/$M_{0}$) with a high
metallicity SSP would therefore be inconsistent, i.e., would assume that
massive stars with strong line-driven winds typically produce XRBs with tight
orbits and massive accreting compact objects.

In Figure 4 we show the [SXPULX \+ SSP]($t_{\rm burst}$, $Z$) models
normalized to 1 $M_{\odot}$ (solid lines), where each panel corresponds to the
annotated instantaneous burst age. The SXPULX 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 SXPULX 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 [SXPULX \+
SSP]($t_{\rm burst}$, $Z$) models, which span wavelengths 1–108 Å, are used as
the spectral input for photoionization modeling with CLOUDY.

##  4 Photoionization Simulations with CLOUDY

We perform photoionization simulations using the [SXPULX \+ SSP]($t_{\rm
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$ = 1019 cm and a constant density $n_{\rm H}=100$
cm-3. The cloud is ionized by a central source, which we set as the [SXPULX \+
SSP]($t_{\rm burst}$, $Z$) SED, implying coincident mixing between the SXPULX
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 $\sim$ 0.3–3 $R_{\rm
S}$ over the entire simulation grid, where $R_{\rm S}$ is the
Strömgren666$R_{\rm S}=(\frac{3\mathcal{Q}_{\rm H}}{4\pi~{}\alpha(T_{\rm
e})n_{\rm e}n_{\rm H}})^{1/3}$, where $\mathcal{Q}_{\rm H}$ is the hydrogen
ionizing photon rate, $\alpha(T_{\rm e})$ is the recombination rate, $n_{\rm
e}$ is the electron density, and $n_{\rm H}$ is the hydrogen density. To
obtain the range of $R_{\rm S}$ from our simulations, we calculate
$\mathcal{Q}_{\rm H}$ directly from the input SEDs, use $\alpha(T_{\rm e})$ =
2.6 $\times$ 10-13 cm3 s-1 (corresponding to $T_{\rm e}$ = 104 K, as in Jaskot
& Ravindranath 2016), and assume $n_{\rm e}$ = $n_{\rm H}$ = 100 cm-3. radius.

We set $i$ = 0 for the SEDULX component of the [SXPULX \+ SSP]($t_{\rm
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 SEDULX
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 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 [SXPULX \+ SSP]($t_{\rm 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 $\mathcal{U}$:

$\mathcal{U}\equiv\frac{\mathcal{Q}_{\rm H}}{4\pi R^{2}n_{\rm H}c}$ (10)

where $\mathcal{Q}_{\rm
H}\equiv\int_{\nu_{0}}^{\infty}\frac{f_{\nu}}{h\nu}d\nu$ is the rate of
emitted photons capable of ionizing hydrogen (i.e., $h\nu_{0}=13.6$ eV), $R$
is the inner cloud radius (hydrogen ionized region) in cm, $n_{\rm 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~{}\mathcal{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 [SXPULX \+ SSP]($t_{\rm burst}$, $Z$,
$\log~{}\mathcal{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_{\odot}$ burst mass would produce relatively few ionizing
photons. Because $R$ and $n_{\rm H}$ remain fixed777We note that fixing $R$
and varying only ionizing flux means certain geometric effects are not
accounted for in these models; such geometrical effects can substantially
change both high-ionization and neutral line emission (e.g., Ramambason et
al., 2022). in our models, changing $\log~{}\mathcal{U}$ amounts to changing
$\mathcal{Q}_{\rm 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~{}\mathcal{U}$. Following the same convention as Byler et al. (2017), we
define the rate of ionizing photons per solar mass as $\hat{\mathcal{Q}}_{\rm
H}\equiv\mathcal{Q}_{\rm H}$ / $M_{\odot}$, where $\mathcal{Q}_{\rm H}$ is the
rate of ionizing photons corresponding to the $\log~{}\mathcal{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_{\odot}$ stellar mass through multiplication by
$\hat{\mathcal{Q}}_{\rm H}$/$\mathcal{Q}_{\rm H}$. This normalization of input
or output intensities using $\hat{\mathcal{Q}}_{\rm H}$/$\mathcal{Q}_{\rm H}$
is specified throughout, where appropriate.

Figure 5: The composite SXPULX \+ SSP SEDs at 0.1 $Z_{\odot}$ for a selection
of instantaneous burst ages, normalized by a range of $\log~{}\mathcal{U}$
values, as noted in the colorbar. Models with high $\log~{}\mathcal{U}$ have
the highest intensities, and correspond to the largest stellar mass formed in
a burst. As in Figure 4, 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.

With these parameters ($t_{\rm burst}$, $Z$, and $\log~{}\mathcal{U}$), the
overall grid input to CLOUDY consists of 147 separate models, corresponding to
different combinations of the seven values for $t_{\rm burst}$, three values
for $Z$, and seven values for $\log~{}\mathcal{U}$. In Figure 5, we show the
[SXPULX \+ SSP]($t_{\rm burst}$, $Z$, $\log~{}\mathcal{U}$) models for select
burst ages and the full range for $\log~{}\mathcal{U}$ for the 0.1 $Z_{\odot}$
case. The models with $\log~{}\mathcal{U}$ = $-$1 have the highest
intensities, or correspondingly, the largest stellar mass formed in a burst
($\mathcal{Q}_{\rm H}$/$\hat{\mathcal{Q}}_{\rm H}$ $\sim 10^{4}-10^{6}$
$M_{\odot}$ for $\log~{}\mathcal{U}$ = $-$1 and $t_{\rm 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
Appendix in Table A.1, including vacuum wavelength (Å), line name, and CLOUDY
specific line ID.

##  5 Results

Figure 6: 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_{\odot}$ and
$\log~{}\mathcal{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.

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

###  5.1 Characteristics of Nebular Emission Due to SXPULX and Potential
SXPULX Diagnostics

A key result from the photoionization simulations is that the inclusion of the
SXPULX 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 SXPULX. 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 SXPULX in the bottom panel. The
nebular contribution to the total flux is high ($>$ 50%) at all wavelengths on
short burst timescales ($\leq$ 3 Myr, before the SXPULX 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
SXPULX is typically very small ($<$ 5%). In the broad-band sense, the major
effect of the SXPULX 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 SEDULX as extended through the EUV, a key outcome from
the addition of the SXPULX is a change to the intensity of lines with
excitation potentials $\geq$ 54 eV. The CLOUDY line list for our simulations
(Table A.1) includes 63 lines redward of Ly$\alpha$ with ionization potential
$\geq$ 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 high-excitation lines may be used as diagnostics for
SXPULX ionization, we consider lines redward of Ly$\alpha$ that satisfy the
following criteria: (1) are relatively strong in the models with SXPULX
contribution, as measured in terms of median line flux relative to median flux
in H$\beta$ or Pa$\beta$888We use H$\beta$ as the reference line for the high-
ionization lines, which are mostly in the UV and optical, and Pa$\beta$ for
coronal lines in the IR. (i.e., log($\widetilde{f}_{\rm
line}/\widetilde{f}_{{\rm H}\beta}$) $\geq-3$ or log($\widetilde{f}_{\rm
line}/\widetilde{f}_{{\rm Pa}\beta}$) $\geq-3$); and (2) have a median
enhancement in line flux by a factor of 2$\times$ relative to the SSP-only
models (log($\widetilde{f}_{\rm line,SXP+SSP}/\widetilde{f}_{\rm line,SSP}$)
$\geq$ 0.3). There are only a select few emission lines satisfying both of
these criteria: He II $\lambda$1640,4686, and [O IV] 25.8832$\mu$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 SXPULX 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 cross-calibration. In all line diagnostic diagrams that
follow, the [SXPULX \+ SSP]($t_{\rm burst}$, $Z$, $\log~{}\mathcal{U}$) models
(i.e., models with SXPULX contribution) are shown as the solid lines, where
orange to dark blue lines represent metallicities from $Z$ = 0.01–1
$Z_{\odot}$ and light blue to dark blue lines represent ionization parameters
from $\log~{}\mathcal{U}$ = $-4$ to $-1$. For all grids in the diagnostic
plots, the star denotes the model with lowest metallicity (0.01 $Z_{\odot}$)
and highest ionization parameter ($\log~{}\mathcal{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 SXPULX 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] $\lambda$1661,6)/(C III]
$\lambda$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
$\lambda$1640/C III] $\lambda$1907,9 ratio in this same diagram is sensitive
to the hardness of the ionizing spectrum. For burst timescales where the
SXPULX contributes ($\geq$ 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 SXPULX, resulting in stronger intensity of
He II $\lambda$1640 relative to C III] $\lambda$1907,9 at the same ionization
parameter.

Figure 7: Potential UV diagnostic diagram for SXPULX ionization, as a function
of instantaneous burst age. The [SXPULX \+ SSP]($t_{\rm burst}$, $Z$,
$\log~{}\mathcal{U}$) grid is shown as solid lines, where orange to dark blue
lines represent metallicities 0.01 $Z_{\odot}$–$Z_{\odot}$ and light blue to
dark blue lines represent $\log~{}\mathcal{U}$ = $-4$ to $-1$, as illustrated
by the colorbars. In all panels, the star represents the 0.01 $Z_{\odot}$ and
$\log~{}\mathcal{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~{}\mathcal{U}$, thereby increasing the
intensity of He II $\lambda$1640 relative to C III] $\lambda$1907,9. This
effect is more pronounced at solar metallicity for $t_{\rm burst}$ = 5 Myr,
following depletion of WR stars, and low metallicities with $t_{\rm burst}$
$>$ 10 Myr (bottom panels), following depletion of the most massive stars.

For this set of line diagnostics, we compare the [SXPULX \+ SSP]($t_{\rm
burst}$, $Z$, $\log~{}\mathcal{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 [SXPULX \+ SSP]($t_{\rm
burst}$, $Z$, $\log~{}\mathcal{U}$) grid for ionization parameters
$\log~{}\mathcal{U}$ $\geq$ $-3.0$, the inclusion of the SXPULX slightly
broadens the range of metallicities consistent with the data points; whereas
the BCD sample is consistent with the SSP-only grid for $Z$ $\leq$ 0.1
$Z_{\odot}$, the composite grid extends the consistency to 0.1 $Z_{\odot}$
$\leq$ $Z$ $<$ $Z_{\odot}$ for older burst ages. This BCD sample has a median
gas-phase metallicity of $\sim$ 0.15 $Z_{\odot}$, suggesting that if SXPULX
ionization is important in these galaxies the population may be dominated by
star formation bursts $\geq$ 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 SXPULX grid is
preferred for such galaxies. Nonetheless, the comparison serves to illustrate
how the addition of the SXPULX 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 SXPULX
hardens the ionizing spectrum at a given $\log~{}\mathcal{U}$, which in effect
mimics an SSP with slightly lower metallicity.

Figure 8: Same as Figure 7, but for a potential optical diagnostic diagram
accounting for the inclusion of SXPULX ionization. The data points are a
sample of star-forming galaxies from Shirazi & Brinchmann (2012) with strong
nebular He II $\lambda$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 from Shirazi & 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_{\rm burst}$ $>$ 10 Myr
and $Z$ $\leq$ 0.1 $Z_{\odot}$; however, neither the [SXPULX \+ SSP]($t_{\rm
burst}$, $Z$, $\log~{}\mathcal{U}$) nor the SSP-only grids are capable of
reproducing some of the more extreme galaxies, particularly those with
detected WR features and higher N II/H$\alpha$ near the empirical maximum
starburst line.

For an optical emission line diagnostic we present Figure 8, where the [N II]
$\lambda$6584/H$\alpha$ ratio is sensitive to metallicity and
$\log~{}\mathcal{U}$ (among other parameters, e.g., Kewley et al. 2013), and
the He II $\lambda$4686/H$\beta$ 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 $\lambda$4686/H$\beta$, as the dashed
grids in Figure 8 illustrate. The [SXPULX \+ SSP]($t_{\rm burst}$, $Z$,
$\log~{}\mathcal{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_{\rm burst}>$ 10 Myr and lower metallicity ($Z$
$\sim$ 0.1 $Z_{\odot}$). In this case, the ability of the [SXPULX \+
SSP]($t_{\rm burst}$, $Z$, $\log~{}\mathcal{U}$) grid to explain a subset of
the galaxies with strong He II/H$\beta$ where SSP-only models cannot is again
a consequence of a harder radiation field due to the addition of the SXPULX.
However, the magnitude of the increase to the ionizing photon rate due to the
SXPULX 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 SXPULX ionizing contribution or SSP-only
photoionization. This is primarily because such galaxies have large ratios of
both He II $\lambda$4686/H$\beta$, implying hard ionizing spectra, and [N II]
$\lambda$6584/H$\alpha$, 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 Figure 9 for lines
accessible with JWST (e.g., Weaver et al., 2010). Here the [O IV]
25.9$\mu$m/[Ne III] 15.66$\mu$m ratio is sensitive to hardness of the ionizing
spectrum, while [Ne III] 15.66$\mu$m/[Ne II] 12.81$\mu$m is primarily
sensitive to the ionization parameter. For this diagnostic, we include the
classification regions from Richardson et al. (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 103 $M_{\odot}$
IMBH with fractional contribution 4–16% relative to a 20 Myr BPASS SSP
(described in more detail in Section 6.3.2). The [SXPULX \+ SSP]($t_{\rm
burst}$, $Z$, $\log~{}\mathcal{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 SXPULX hardens the ionizing spectrum.
The models including SXPULX contribution veer into the composite region for
some of the solar metallicity $t_{\rm burst}$ = 5 Myr models, and into the AGN
ionization region only for lowest metallicities ($\lesssim$ 0.1 $Z_{\odot}$)
and highest ionization parameters ($\log~{}\mathcal{U}$ $>$ -3) for bursts
with $t_{\rm burst}$ $>$ 10 Myr. This indicates, not surprisingly, that the
[SXPULX \+ SSP]($t_{\rm burst}$, $Z$) photoionization signature looks less
like pure SF when the fractional contribution of the SXPULX 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
($\sim$ 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
SXPULX 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.

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 103
$M_{\odot}$ 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_{\rm burst}$, $Z$, $\log~{}\mathcal{U}$) grid for
values $\log~{}\mathcal{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.

Finally, we comment on high-ionization lines that do not satisfy both criteria
in terms of strength relative to H$\beta$ or Pa$\beta$ and enhancement
relative to the SSP-only case. These lines fall broadly into three categories:
(1) lines that are strong relative to H$\beta$ or Pa$\beta$, but not
significantly enhanced relative to the SSP-only models (flagged with † in
Table A.1); (2) lines that are weak relative to H$\beta$ or Pa$\beta$, 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($\widetilde{f}_{\rm line}/\widetilde{f}_{{\rm H}\beta}$) $\sim-3$ and
median enhancement of 1.004$\times$ in flux relative to the SSP-only models.
By contrast, lines flagged ‡ in Table A.1 have log($\widetilde{f}_{\rm
line}/\widetilde{f}_{{\rm Pa}\beta}$) $\sim-9$, but median enhancement of
47000$\times$ 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 SXPULX component
is substantial. As such, there is a higher diagnostic potential for discerning
the SXPULX ionization contribution using coronal lines, many of which are in
the IR, but only if the fractional contribution of the SXPULX 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 SXPULX and weak lines heavily enhanced due to the addition of the
SXPULX—via additional emission line diagnostics in Figures 10–12. In Figure
10, we show the [SXPULX \+ SSP]($t_{\rm burst}$, $Z$, $\log~{}\mathcal{U}$)
grid relative to the same sample of local BCD galaxies from as in Figure 7.
Here, the grid with SXPULX contribution is only slightly offset to higher
values of C IV $\lambda$1548,51/O III] $\lambda$ 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 SXPULX 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 $\sim$ 126 eV and 64 eV, respectively. In Figure
11, the SSP-only models rarely reach log([Ne V] $\lambda$3426/[Ne III]
$\lambda$3870) $>-5$, while the grid with SXPULX 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
[SXPULX \+ SSP]($t_{\rm burst}$, $Z$, $\log~{}\mathcal{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 SXPULX contribution to the ionizing photon budget comes
at the expense of line strength. The ratio of [Ne V] 14.3$\mu$m/[Ne III]
15.56$\mu$m traced out by the grid with SXPULX 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
SXPULX ionizing contribution for a narrow range in $\log~{}\mathcal{U}$,
$t_{\rm 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 SXPULX 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).

Figure 10: 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: Same as Figure 7, but for a potential UV/optical diagnostic diagram
for SXPULX contribution. For these high-ionization 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~{}\mathcal{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 from Cleri 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_{\rm
burst}$, $Z$, $\log~{}\mathcal{U}$) models in a very narrow range of grid
space. Figure 12: Same as Figure 11, but for a potential IR diagnostic, albeit
using a weak coronal line that is highly sensitive to hardness of the ratio
field. As in Figure 9, the IMBH models from Richardson et al. (2022) for a
selection of IMBH fractional contributions are shown for reference. Figure 13:
Nebular emission line luminosities normalized by stellar mass for He II
$\lambda$4686, H$\beta$, [Ne V] $\lambda$3426, and [Ne III] $\lambda$3870 as a
function of burst age. As before, in each panel dark blue corresponds to
$Z_{\odot}$, purple corresponds to 0.1 $Z_{\odot}$, and orange corresponds to
0.01 $Z_{\odot}$ models. Increasing line opacity denotes increasing
$\log~{}\mathcal{U}$. The ratio of a high-ionization line (e.g., He II, [Ne
V]) with significant SXPULX contribution relative to a lower ionization line
(e.g., H$\beta$, [Ne III]) 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.

###  5.2 Consequences of Time Dependence of the SXPULX Relative to the SSP

The key parameters affecting the time-dependence of the nebular emission from
the [SXPULX \+ SSP]($t_{\rm burst}$, $Z$, $\log~{}\mathcal{U}$) grid are (1)
the selection to impose a delay time for the SXPULX to turn on and; (2) the
combination thereafter with SSPs of corresponding $t_{\rm burst}$. Imposing a
delay time creates a scenario in which the SXPULX 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 SXPULX component is indeed
descended from the stellar component, allowing time for compact object
formation and accretion to begin. In this way, the SXPULX, 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 SXPULX
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_{\rm
X}$/$M_{0}$ scaling for the SXPULX is largest for younger instantaneous bursts
at all metallicities (Figure 1), the SXPULX generally hardens the incident
ionizing spectrum more profoundly relative to the BPASS SSPs at older burst
ages. That is, the fractional contribution of the SXPULX 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 SXPULX component remains appreciable.
Given the imposed delay time and coupling with BPASS SSPs, the SXPULX
contribution to the total rate of ionizing photons is $\lesssim$ 2% for burst
ages between 5–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 SXPULX contribution
for a narrow range of the grid corresponding to $\sim$ 0.1 $Z_{\odot}$,
$\log~{}\mathcal{U}$ $>-2$, and $t_{\rm burst}$ $>$ 10 Myr. This suggests
SXPULX 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~{}\mathcal{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 $\lesssim$ 10 Myr, and may therefore require
a source other than an SXPULX 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 SXPULX for $t_{\rm burst}>$ 10 Myr and low
metallicities. For SFHs that include a young ($<$ 10 Myr) burst component, the
SXPULX 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_{\rm burst}<$ 10 Myr,
resulting in strong emission from lines such as H$\beta$. This strong stellar
ionizing flux effectively dilutes the SXPULX 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
$\lambda$4686, H$\beta$, [Ne V] $\lambda$3426, and [Ne III] $\lambda$3870 as a
function of burst age. In comparing the evolution of the intensities of He II
and H$\beta$, it is clear both lines reach their absolute maximum in intensity
for $t_{\rm burst}\sim$ 3–5 Myr; however, while the H$\beta$ 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 SXPULX
contribution, resulting in a higher ratio of He II $\lambda$/H$\beta$ at later
times. The same is true for [Ne V] $\lambda$3426 relative to [Ne III]
$\lambda$3870\. For this ratio, it is also particularly evident when the
SXPULX contribution turns for $t_{\rm burst}>$ 3 Myr, as there is no
appreciable [Ne V] emission prior to this.

To further illustrate the timescale dependent nature of the relative SXPULX
ionizing contribution, we reran the photoionization simulations without the
strict delay-time dependence for the SXPULX component. For these “no-delay
time” simulations, we extrapolated $L_{\rm X}$/$M_{0}$ from Figure 1 to burst
ages 1–3 Myr, allowing the SXPULX to contribute to the ionizing budget at
every time step. As Figure A.1 shows, the immediate formation of the SXPULX
alongside a 1 Myr or 3 Myr burst only weakly increases He II/H$\beta$, 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 SXPULX may contribute significantly to high-
energy ionizing photon production only in some environments, at least for the
SEDULX considered in this work. Due to the time-dependent nature of the SXPULX
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 SXPULX ionizing contribution.

###  5.3 Consequences of Metallicity-Dependence of the SXPULX and Adopted
Abundance Patterns

The photoionization simulation results are based on the assumption that the
metallicity of the SXPULX corresponds to the stellar metallicity of the SSP,
as would be expected if the SXPULX 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
SXPULX with higher formation efficiency ($L_{\rm X}$/$M_{0}$) than the SXPULX
combined with a solar-metallicity SSP.

The imprint of the metallicity-dependent normalization of the SXPULX is
particularly evident in emission line diagnostics that include metallicity-
sensitive line ratios (e.g., [N II]/H$\alpha$). To illustrate the magnitude of
this effect, we run simulations where we remove the correlation between the
SXPULX normalization and stellar metallicity. In these “no stellar metallicity
dependence” simulations, we allow an SSP at a given $t_{\rm burst}$ to be
coupled with an SXPULX normalized using $L_{\rm 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
$\lambda$4686/H$\beta$ are produced at a given [N II]/H$\alpha$, particularly
for the solar metallicity case. This results in much more overlap between the
grid with the SXPULX contribution and the observed data points, for a broader
range of burst ages; however, we note that such an elevated $L_{\rm
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_{\rm X}$ $>$ 1039 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 SXPULX 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 simulations are also of critical importance for
select emission line diagnostics. For example, though $L_{\rm 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 [SXPULX \+ SSP]($t_{\rm burst}$, $Z$,
$\log~{}\mathcal{U}$) grid in Figure 12, where the lowest metallicity models
including SXPULX contribution do not uniformly have the strongest [Ne V]
14.3$\mu$m/[Ne III] 15.56$\mu$m ratios, as would be expected if SXPULX
production efficiency (i.e., $L_{\rm 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, star-forming 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 SXPULX 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.

###  5.4 Availability of Models

We make our models and simulation results available in two separate formats,
as outlined below.

The models are available directly from FSPS
999https://github.com/cconroy20/fsps; This version currently only supports
adding SXPULX contribution alongside BPASS SSPs.. From FSPS, one can build
SEDs for a population with or without SXPULX contribution, including
consideration of nebular line and continuum emission. In this way, users can
simulate spectro-photometric data for populations including SXPULX
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
repository101010https://github.com/kgarofali/sxp-cloudy. In this repository we
provide code examples in python for reading and plotting simulation results,
allowing users to quickly plot the [SXPULX \+ SSP]($t_{\rm burst}$, $Z$,
$\log~{}\mathcal{U}$) grid for different emission line diagnostics. Users can
therefore compare data or other model grids to these simulations including
SXPULX 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.

##  6 Discussion

We have already presented potential UV, optical and IR emission line
diagnostics for SXPULX ionization in Section 5. Here we present best practices
for diagnosing SXPULX 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 SXPULX model presented here.

###  6.1 Caveats For Applying SXPULX Formalism to X-ray Observations

The emergent $L_{\rm X}$ (0.5–8 keV) from the SXPULX models is very similar to
the input $L_{\rm X}$ (i.e., Figure 1), as the hard photons are attenuated
very little by the cloud. However, comparison of the simulation results with
X-ray observations requires a more careful transformation of the intrinsic
$L_{\rm 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_{\rm X}$ values from the
photoionization simulations is that the SEDULX 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 inner 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_{\rm X}$ (obs) is likely dependent on inclination angle, and possibly also
on the mass supply rate itself, which in the SEDULX 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 SEDULX in black with different linestyles,
ranging from a face-on viewing angle ($i$ = 0∘) to the edge-on case ($i$ =
90∘), wherein the SEDULX is primarily an EUV source. However, these are all
for a fixed funnel opening angle ($\theta$), 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 $\theta$
(or, correspondingly, $\dot{m}$): 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 self-irradiation 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 SEDULX 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 ultra-luminous 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 edge-on, 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 SEDULX 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 $\theta$
and $\dot{m}$ 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 age-dependence
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 SXPULX $L_{\rm X}$ via a
comparison of simulated and observed He II $\lambda$4686/H$\beta$ as a
function of $L_{\rm 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 star-forming 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_{\rm X}$ corresponds roughly to
the $L_{\rm 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_{\rm 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_{\rm X}$ (or upper limits) are
determined from the same spectroscopic aperture, or within 1.4$\arcsec$ 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 SXPULX $L_{\rm X}$ (both in terms
of stellar mass) assuming a constant SFR of 1 $M_{\odot}$ yr-1. The results
are shown in Figure 14, with the continuous SFH models as circles and the
burst models as shaded regions, where the extent of the region encompasses the
range of simulated $\log~{}\mathcal{U}$ and $t_{\rm burst}$. The most extreme
line ratios are achieved for high $\log~{}\mathcal{U}$ and older $t_{\rm
burst}$. The simulations can reproduce much of the observed range of He II
$\lambda$/H$\beta$ (green diamonds) for reasonable $L_{\rm 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_{\rm burst}>$ 10 Myr.
For the continuous SFH case, the models produce much more modest He II
$\lambda$/H$\beta$, 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$\alpha$ or
H$\beta$ equivalent width, are therefore integral to determining the
importance of SXPULX high-energy ionizing contribution.

Figure 14: The emergent $L_{\rm X}$ (0.5–8 keV)/SFR relative to the strength
of He II $\lambda$4686/H$\beta$ 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 from Senchyna 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 $L_{\rm X}$/SFR; however, this requires instantaneous
bursts where $t_{\rm burst}>$ 10 Myr. Models with continuous SFHs produce
$L_{\rm X}$/SFR in the range of the observed galaxies and upper limits, but
much lower ratios of He II $\lambda$4686/H$\beta$.

As some of the observed galaxies in Figure 14 have values of $L_{\rm 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_{\rm 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 ($\gtrsim$ 3 $M_{\odot}$ yr-1), where the effects
of stochastic sampling of the XLF are on the order of, or smaller than, the
model uncertainties on $L_{\rm X}$, and where the strict linear scaling of
$L_{\rm 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_{\rm X}$ (obs) can be subject to large ($\sim$ 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_{\rm X}$/SFR ($\gtrsim$ 1041 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_{\odot}$ 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 $\sim$ 2$\sigma$ 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_{\rm 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_{\rm X}$ from a given SXPULX 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 SXPULX 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_{\rm X}$. For highly star-forming galaxies, the
population average $L_{\rm 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.

###  6.2 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
$\lambda$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
$\lambda$4686\. We focus on the former here for parity with our own
implementation.

In Figure 15 we show the SEDULX for a 20 Myr instantaneous burst of star
formation at 0.1 $Z_{\odot}$ in yellow. Here the SEDULX has been normalized to
1 $M_{\odot}$ stellar mass formed (Section 3), corresponding to $L_{\rm X}$
(0.5–8 keV) $\sim$ 3 $\times$ 1032 erg s-1. The corresponding SSP ($t_{\rm
burst}$ = 20 Myr, 0.1 $Z_{\odot}$) 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 SXPULX framework.

The photoionization simulations performed in Senchyna et al. (2020) explored
the effect of HMXBs on production of He II $\lambda$4686 and [Ne V]
$\lambda$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
($\sim$ 10–100 $M_{\odot}$) 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_{\rm
X}$/SFR = 10${}^{40}-10^{44}$ erg s-1 ($M_{\odot}$ 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~{}\mathcal{U}$. With these
simulation inputs, the authors conclude that the observed range of He II
$\lambda$4686/H$\beta$ can only be reproduced for extremely high X-ray
production efficiencies ($L_{\rm X}$ (0.5–8 keV)/SFR $>$ 1042 erg s-1
($M_{\odot}$ 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 $\lambda$4686/H$\beta$. As
already discussed in Section 6.1, our simulations are able to reproduce more
of the observed range of He II $\lambda$4686/H$\beta$ within the observed
range of $L_{\rm X}$/SFR, but primarily for instantaneous bursts with $t_{\rm
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 $\lambda$4686/H$\beta$ at the simulated range of
$L_{\rm 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 $\lambda$4686/H$\beta$ at a given $L_{\rm 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 SEDULX is much flatter through
the EUV than the multi-color disk model employed, and therefore produces more
ionizing EUV photons per $L_{\rm X}$.

Figure 15: The intrinsic SED for the composite model (SXPULX \+ SSP, solid
black line) with $t_{\rm burst}$ = 20 Myr and $Z$ = 0.1 $Z_{\odot}$,
normalized to 1 $M_{\odot}$ 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 $L_{\rm X}$ (0.5–8 keV) (or, in the
case of the “BB”, $L_{\rm EUV}$) as the SEDULX. The grey vertical lines mark
ionization potentials for select lines, as labeled.

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.005111111These 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 $\sim$ 0.3
$Z_{\odot}$, which we assume corresponds to the BPASS $Z$ = 0.005 model.)
again for a similar range in $\log~{}\mathcal{U}$ as employed here. In their
simulations, the ULX contribution is scaled relative to the BPASS SSP for a
range of $L_{\rm X}$/SFR, motivated by empirical scalings and scatter therein.
These authors find they can reproduce some of the observed range of He II
$\lambda$4686/H$\beta$ for reasonable X-ray production efficiencies ($L_{\rm
X}$/SFR $\simeq$ 10${}^{40}-10^{41}$ erg s-1 ($M_{\odot}$ yr-1)-1), roughly in
agreement with our findings for instantaneous bursts and the SEDULX model.

In general, the simulations results presented in Simmonds et al. (2021) also
achieve higher values of He II $\lambda$4686/H$\beta$ than our grids with
SXPULX 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 ULX-type 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 SEDULX, leading to a different shape when
extended into the EUV, particularly around 54 eV where He II $\lambda$4686 has
its ionization potential. While our SXPULX simulations are run coupled to a
range of SSPs with $L_{\rm 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, $\sim$ 0.3
$Z_{\odot}$) using a range of $L_{\rm X}$/SFR. This amounts to relaxing the
assumptions used here for scaling of the SXPULX 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 SXPULX simulations: using a broad range of
$L_{\rm X}$/SFR for a given $Z$ and $t_{\rm 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 SEDULX component,
despite very different approaches in (re)constructing the intrinsic ionizing
spectrum.121212The ionizing spectra appear more discrepant in the UV–to–IR
range, however, this regime was not the focus of the modeling in Umeda et al.
(2022), and therefore had the smallest contribution (in terms of data sets
used) to constraining the best-performing model. In general, using their full
galaxy sample, Umeda et al. (2022) show that the intensity of the observed
high-ionization 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 [SXPULX \+ SSP]($t_{\rm 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 high-ionization emission lines in a set of
nearby EELGs with very young ($\lesssim$ 10 Myr) and metal-poor ($\sim$ 0.1
$Z_{\odot}$) 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 dash-dot line
(“BB”). As already discussed in Section 6.1, the SEDULX 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 SEDULX 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 SEDULX 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 best-performing models depend critically on the suite of
lines available for fitting as well as the set of parametric models
considered. Whether using parametric or physically-motivated 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.

###  6.3 Alternative Ionizing Sources

For the SXPULX model used in this work, we have assumed the accretor is a
stellar mass BH ($\sim$ 10–100 $M_{\odot}$) and have constructed the SEDULX
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.

####  6.3.1 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 SXPULX 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_{\rm burst}$ $\leq$ 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 SXPULX would therefore further prolong the ionizing
output of a population, beyond that of the most massive stars and binaries,
and the BH SXPULX 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_{\rm X}$/$M_{0}$ will be
both SFH- and metallicity-dependent, and such a population break-down 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 X-ray
source ($L_{\rm X}$ $>$ 1041 erg s-1), which implies the presence of an
accreting NS with very strong magnetic field (Brightman et al., 2018, 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_{\rm M}$) compares to the spherization
radius ($R_{\rm 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 SXPULX a promising
avenue for future study, but outside the scope of the present work.

####  6.3.2 Intermediate Mass Black Holes

Moving up in the mass scale, IMBHs (100 $M_{\odot}$ $\leq$ $m$ $\leq$ 106
$M_{\odot}$) 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 SEDULX with IMBH SEDs for a 103 $M_{\odot}$ BH and 105 $M_{\odot}$
BH based on the qsosed model (Kubota & Done, 2018), all normalized to the same
$L_{\rm 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 broad-band 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).

Figure 16: Similar to Figure 15, but for a comparison with 103 $M_{\odot}$ and
105 $M_{\odot}$ 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 $L_{\rm X}$ (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–to–soft 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.

The full set of IMBH photoionization simulations from Richardson et al. (2022)
is expansive, encompassing different 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~{}\mathcal{U}$. We