11institutetext: Universitäts-Sternwarte, Fakultät für Physik, Ludwig-
Maximilians-Universität München, Scheinerstr.1, 81679 München, Germany
22institutetext: INAF, Osservatorio Astronomico di Trieste, via Tiepolo 11,
I-34131, Trieste, Italy 33institutetext: Max-Planck-Institut für Astrophysik,
Karl-Schwarzschild-Straße 1, 85748 Garching bei München, Germany

# Decomposition of galactic X-ray emission with Phox

Contributions from hot gas and XRBs

S. Vladutescu-Zopp , 11<EMAIL_ADDRESS>V. Biffi 22 K. Dolag, 1133

(Received 08/2022; accepted xx/xx/xxxx)

###### Abstract

Context. X-ray observation of galaxies with high spatial resolution
instruments such as Chandra revealed that major contributions to their diffuse
emission originate from X-ray bright point sources in the galactic stellar
field. It is established that these point sources are accreting compact
objects with stellar donors in a binary configuration called X-ray binaries.
They are classified according to the predominant accretion process: wind-fed
in case of high-mass donors, Roche-lobe mass transfer in case of low-mass
donors. Observationally, it is challenging to disentangle these two
populations reliably because of their similar spectra.

Aims. We provide a numerical framework with which spatially and spectrally
accurate representations of X-ray binary populations can be studied from
hydrodynamical cosmological simulations. We construct average spectra
accounting for a hot gas component and verify the emergence of observed
scaling relations between galaxy wide X-ray luminosity ($L_{X}$) and stellar
mass ($M_{*}$) as well as star-formation rate (SFR).

Methods. Using simulated galaxy halos extracted from the
$(48\,h^{-1}\mathrm{cMpc})^{3}$ volume of the Magneticum Pathfinder
cosmological simulations at $z=0.07$ we generate mock spectra with the X-ray
photon-simulator Phox. We extend the Phox code to account for the stellar
component in the simulation and study the resulting contribution in composite
galactic spectra.

Results. Well known X-ray binary scaling relations with galactic SFR and
$M_{*}$ emerge self-consistently, verifying our numerical approach. Average
X-ray luminosity functions are perfectly reproduced up to the one-photon
luminosity limit. Comparing our resulting $L_{X}-\mathrm{SFR}-M_{*}$ relation
for X-ray binaries with recent observations of field galaxies in the Virgo
galaxy cluster we find significant overlap. Invoking a metallicity dependent
model for high-mass X-ray binaries yields an anti-correlation between mass-
weighted stellar metallicity and SFR normalized luminosity. The spatial
distribution of high-mass X-ray binaries coincides with star-formation regions
of simulated galaxies while low-mass X-ray binaries follow the stellar mass
surface density. X-ray binary emission is the dominant contribution in the
hard X-ray band (2-10 keV) in the absence of an actively accreting central
super-massive black hole with $\sim 50\%$ contribution in the soft X-ray band
(0.5-2 keV) rivaling the hot gas component.

Conclusions. We conclude that our modelling remains consistent with
observations despite uncertainties connected to our approach. The predictive
power and easily extendable framework hold great value for future
investigations of galactic X-ray spectra.

###### Key Words.:

X-ray: galaxies - X-ray: binaries - X-ray: ISM - methods: numerical

##  1 Introduction

X-ray properties of local normal galaxies have been extensively studied in the
past where normal describes galaxies without a luminous active galactic
nucleus (AGN). One of the most challenging tasks has been the identification
of physical processes producing observed X-ray emission from those galaxies.
It is a well established fact that most of their X-ray power-output originates
from diffuse hot gas in the inter-stellar medium (ISM) and accretion powered
point-sources called X-ray binaries (XRB) (see reviews Fabbiano, 1989, 2006,
2019) excluding the emission from actively accreting super massive black holes
(SMBH).

Numerous studies of the hot diffuse ISM in star-forming normal galaxies
revealed a tight linear relation between total X-ray luminosity $L_{X}$ and
star-formation rate (SFR) (Strickland & Stevens, 2000; Ranalli et al., 2003;
Tyler et al., 2004; Mineo et al., 2012b; Li & Wang, 2013; Kouroumpatzakis et
al., 2020) due to feedback from young stellar objects and energy injection
from super novae (SN) which heat up the ISM. Recent studies suggest that the
ISM emission normalization in X-rays declines with increasing total gas-phase
metallicities possibly due to the effect of metal absorption and changes in
stellar feedback within the ISM (Garofali et al., 2020; Lehmer et al., 2022).

The ISM of normal elliptical galaxies (EGs) depends on their dynamical state
(Kim & Fabbiano, 2015; Babyk et al., 2018) and strongly correlates with the
ISM temperature (Boroson et al., 2011; Kim & Fabbiano, 2015) and dynamical
mass (Kim & Fabbiano, 2013; Forbes et al., 2017; Babyk et al., 2018). The ISM
of core EGs, where core denotes slow rotating systems with cored central
surface brightness and overall old stellar population, appear to follow gas
X-ray luminosity scaling relations consistent with virialized systems similar
to brightest cluster galaxies (BCG) and groups. Other types of gas poor EGs
behave similar to disk galaxies in these relations: they do not show
significant correlations between hot gas temperature and X-ray luminosity of
hot gas suggesting secondary effects like rotation, flattening and SN feedback
dissipating the ISM (see e.g. Fabbiano, 2019).

Separate scaling relations are found for collective XRB X-ray luminosity
depending on the type of XRB: Low-Mass X-ray binaries (LMXB or LXB) are
accreting compact objects (CO) where the low-mass companion undergoes Roche-
Lobe-Overflow at the end of its life. Thus, LMXBs are generally correlated
with the integrated stellar light ($L_{k}$) and by extension with the total
stellar mass of the host galaxies because of the long evolutionary timescales
of the low-mass ($<1\,M_{\odot}$) donor star in the binary system (Gilfanov,
2004; Boroson et al., 2011; Zhang et al., 2012; Lehmer et al., 2019). This
relation seems to be enhanced by the globular cluster (GC) specific frequency,
$S_{N}=N_{\mathrm{GC}}10^{0.4(M_{V}^{T}+15)}$, where gravitationally
dissolving GCs seed LMXBs into the galactic field in addition to in situ
formed LMXBs Irwin (2005); Zhang et al. (2011); Boroson et al. (2011); Zhang
et al. (2012); Lehmer et al. (2020). In High-mass X-ray binaries (HMXB or HXB)
the donor star is a massive O/B star $(>8\,M_{\odot})$ fueling accretion onto
the CO with intense stellar winds. Due to the short lifetimes of such massive
stars $(<100\mathrm{Myr})$, HMXBs provide an independent tracer for the
galactic star-formation history (SFH). It was found that the combined X-ray
luminosity of HMXBs relates linearly to SFR (Grimm et al., 2003; Shtykovskiy &
Gilfanov, 2005b; Mineo et al., 2012a; Lehmer et al., 2010, 2016, 2019).
Broadband X-ray luminosity functions of XRBs follow distinct power-law (PL) or
broken power-law (BPL) distributions in the local universe (see e.g. Grimm et
al., 2003; Gilfanov, 2004). The exact shape of XRB luminosity functions in
normal galaxies is a function of stellar age and evolves on timescales
consistent with stellar evolution timescales of the XRB donor, transitioning
between HMXB dominated to LMXB dominated at a stellar population age of
$\gtrsim 100\,\mathrm{Myr}$ (Lehmer et al., 2017; Gilbertson et al., 2022).

By combining the scaling relations of both XRB types, we expect a distinct
relation between SFR normalized $L_{X}$ and specific SFR (sSFR = SFR/$M_{*}$)
of the form $L_{X}/\mathrm{SFR}=\alpha\,\mathrm{sSFR}^{-1}+\beta$ where
$\alpha$ describes the contribution from LMXBs and $\beta$ the contribution
from HMXBs. This relation is especially useful in galaxies where XRB
candidates can not be resolved: it shows which of the two XRB types is the
dominant contributor to the galactic $L_{X}$. In previous works it has been
shown that this relation generally holds for local galaxies (Lehmer et al.,
2010, 2016, 2019) with modifications for high redshift galaxies (Lehmer et
al., 2016). The turnover from LMXB to HMXB dominated occurs at
$\log\mathrm{sSFR}\sim-10.5$. Recently, Soria et al. (2022) performed Chandra
observations of normal galaxies in the Virgo cluster which were marginally
compatible with scaling relations for local galaxies indicating non-negligible
environmental effects on this relation. A similar dependency was found by
Inoue et al. (2021) in the form of a fundamental plane
$L_{X}\sim\mathrm{SFR}+\alpha\,M_{*}$.

Numerous studies on generalizing locally derived XRB scaling relations for the
high redshift universe use computationally expensive population synthesis
codes for XRB evolution (Fragos et al., 2013a, b; Madau & Fragos, 2017;
Wiktorowicz et al., 2017). Taking into account the metallicity evolution and
star-formation history (SFH) of the universe constrained with the help of
cosmological simulations (Millenium II in the case of Fragos et al. (2013b))
suggest a more complex evolution of local XRB scaling relations with redshift
and HMXB emissivity dominating above $z\sim 2$ where the cosmic SFR density
peaks. This is supported by constraints from the VANDELS survey on HMXB
emission in high redshift galaxies (Saxena et al., 2021). Additional evidence
for the metallicity dependence of HMXB emissivity comes from recent studies of
Lehmer et al. (2021) where the gas-phase oxygen fraction of their galactic
sample was used to further refine the $L_{X}$-SFR relation. They demonstrate
that galaxies with lower metal content consistently show an increased number
of luminous HMXBs and total luminosity. Similarly, earlier studies from Brorby
et al. (2016) and Fornasini et al. (2020) found anti-correlated metallicity
enhancements in the $L_{X}$-SFR relation consistent with findings by Lehmer et
al. (2021). Additionally, data from the eROSITA Final Equatorial Depth survey
(eFEDS) showed elevated total X-ray luminosity for low-metallicity dwarf
galaxies with high sSFR (Vulic et al., 2021). This paper introduces a
numerical modelling to study the composition of galactic X-ray emission from a
theoretical point of view, using state of the art hydrodynamical cosmological
simulations from the Magneticum set. Instead of applying computationally
expensive population synthesis codes, we associate XRB populations to the
stellar component within the cosmological simulation by making use of
observationally derived XRB luminosity functions. We use observationally
derived local scaling relations to constrain seeding parameters and describe
the numerical setup used to extend the capabilities of the virtual X-ray
photon simulator Phox. Phox has been successfully used in conjunction with
Magneticum to study galaxy cluster $L_{X}$-temperature relations as well as
AGN luminosity functions. Additionally, it was used to predict the
contamination of cluster X-ray emission by AGN for the eROSITA mission (see
Biffi et al., 2012, 2013, 2018). The inclusion of an XRB component enables the
study of galactic X-ray spectra from a theoretical standpoint. We validate our
approach by retrieving well known XRB and gas scaling relations of galaxies
while maintaining a low computational cost.

The paper is structured as follows: In section 2 we describe the state of the
art hydrodynamical cosmological simulation suite Magneticum on which we base
our analysis. In section 3 we describe the methodology we used to constrain
the XRB population within the stellar component of cosmological simulations
and a novel approach to estimate linear SFRs for simple stellar populations.
In chapter 4 we describe the process to generate photons from X-ray emitting
sources in the simulations using the virtual photon simulator Phox. In section
6 we show our results regarding the reconstruction of known XRB and hot gas
scaling relations as well as proper spatial distribution of the two XRB types
in the galactic field. In section 7 we compare average count ratios of
different contributors to galactic X-ray emission and their relative
contribution to the total X-ray emission. In section 8 we put our findings
into context with our methodology and constructed dataset by highlighting
caveats and uncertainties. Finally we summarize our findings in section 9.

##  2 Cosmological hydrodynamical simulation

Magneticum Pathfinder Simulations111Project webpage: www.magneticum.org are a
series of state-of-the-art hydrodynamical cosmological simulations exploring
varying ranges in particle number, volume and resolution. They are based on an
improved version of the N-body code Gadget 3, which is an updated version of
the code Gadget 2 (Springel (2005)) including a Lagrangian method for solving
smoothed particle hydrodynamics (SPH). The code introduces several
improvements regarding the SPH implementation by including a treatment of
viscosity and artificial conduction (Dolag et al., 2005; Beck et al., 2016).
Additional physical processes describing the evolution of the baryonic
component have been implemented. They encompass radiative gas cooling as
described in Wiersma et al. (2009), heating from a uniform time-dependent
ultraviolet background (Haardt & Madau, 2001) and a sub-resolution model for
star-formation with mass-loading rate proportional to SFR and resulting wind-
velocities of $v_{w}=350\,km\,s^{-1}$(Springel & Hernquist, 2003). A treatment
for chemical enrichment of the gaseous component through stellar evolution has
been implemented following the prescription of Tornatore et al. (2004, 2007).
Furthermore, a prescription of SMBH growth and gas accretion, powering energy
feedback for AGN was implemented following Springel (2005) and Di Matteo et
al. (2005) with modifications following Fabjan et al. (2010). Previous studies
using the Magneticum simulations are consistent with observed kinematical and
morphological properties of galaxies (Teklu et al., 2015, 2017; Remus et al.,
2017; Schulze et al., 2018, 2020), chemical properties of galaxies and
clusters (Dolag et al., 2017) as well as statistical properties of AGN
(Hirschmann et al., 2014; Steinborn et al., 2016; Biffi et al., 2018). In this
paper we will make use of a high resolution run of Magneticum called Box4/uhr
which is a $(48\,h^{-1}\,\mathrm{cMpc})^{3}$ comoving volume with a mass
resolution of $m_{\mathrm{DM}}=3.6\times 10^{7}\,M_{\odot}$ and
$m_{\mathrm{gas}}=7.3\times 10^{6}\,M_{\odot}$ for dark matter and gas
respectively, and a total of $576^{3}$ particles. Initial conditions for the
simulations are generated using standard $\Lambda$CDM cosmology using results
from the Wilkinson Microwave Anisotropy Probe (WMAP7) (Komatsu et al., 2011),
with Hubble parameter $h=0.704$, matter density $\Omega_{M}=0.272$, dark
energy density $\Omega_{\Lambda}=0.728$, baryon density $\Omega_{b}=0.0451$
and normalization of the fluctuation amplitude at 8 Mpc $\sigma_{8}=0.809$.

##  3 Modelling of XRB populations within SSPs

In this section we describe the process of modelling the population size of
LMXBs and HMXBs, given a stellar resolution element in the cosmological
simulation. Stellar elements describe a collection of stars as a simple
stellar population (SSP) sharing the same time of birth and initial
metallicity, distributed according to a common initial mass function (IMF).
The main difficulty lies in connecting the X-ray luminosity functions (XLFs),
which are derived from observed galactic properties, with SSP properties
traced by the simulation. In order to preserve the generality of our approach,
we do not ’a priory’ assume a galactic environment of the SSP in the
simulation but rather investigate scaling relations and arising XLFs by
processing each SSP individually.

###  3.1 Theoretical estimation of the LMXB population size

The modelling of LMXB is based on the fact that they can not be associated
with previous star formation events in a galaxy which is due to the long
lifetime of their low mass stellar companions. Thus, LMXBs are more typically
found in old galactic stellar regions which connects the LMXB population to
the total stellar mass of their host galaxy (Gilfanov, 2004; Zhang et al.,
2012; Lehmer et al., 2017, 2019). By requiring that a SSP must be older than 1
Gyr for it to be eligible to host a LMXB population, we can simply obtain the
LMXB population size by renormalizing the LMXB XLF to the mass of the SSP. In
particular, we will use the broken power-law XLF from Zhang et al. (2012)
throughout this paper which has the form:

$\dfrac{dN_{\mathrm{LXB}}}{dL}=\frac{A_{\mathrm{LXB}}}{10^{11}M_{\odot}}\cases{L}^{-\alpha_{1}},&L<L_{b,1}\\\
L_{b,1}^{(\alpha_{2}-\alpha_{1})}L^{-\alpha_{2}},L_{b,1}<L<L_{b,2}\\\
L_{b,2}^{(\alpha_{3}-\alpha_{2})}L^{-\alpha_{3}},L_{b,2}<L<L_{cut}\\\
0\,,L\geq L_{cut}{}\,,$ (1)

where $L_{b,1}=0.546$ and $L_{b,2}=5.99$ are break luminosities, $L_{cut}=500$
is the cut-off luminosity, $A_{\mathrm{LXB}}=54.48$ is the normalization in
units of $10^{11}M_{\odot}$, and $\alpha_{1}=1.02$, $\alpha_{2}=2.06$,
$\alpha_{3}=3.63$ are the respective power-law slopes. Luminosities are in
units of $10^{38}\,\mathrm{erg\,s^{-1}}$. In fact, LMXBs appear to not only be
dependent on the stellar mass of the host galaxy, but are also connected to
the GC specific frequency ($S_{N}$) which enhances the XLF normalization by
accounting for dynamically formed LMXBs in GCs (Irwin, 2005; Boroson et al.,
2011; Zhang et al., 2012; Lehmer et al., 2020). In Lehmer et al. (2020) they
were successful in estimating the field LMXB contribution from different
formation channels while earlier studies like Zhang et al. (2012) combine
these channels into a single XLF. Since it is not feasible to account for GCs
in the context of a cosmological simulation, because of mass resolution
limits, integration of eq. (1) is sufficient to describe the LMXB population
in a SSP as

$\tilde{N}_{\mathrm{LXB}}=M_{\mathrm{SSP}}\times\int_{L_{\mathrm{min}}}^{L_{cut}}\dfrac{dN_{\mathrm{LXB}}}{dL}dL\,,$
(2)

with $M_{\mathrm{SSP}}$ being the current mass of the SSP in solar masses.

###  3.2 Theoretical estimation of the HMXB population size

Due to the short lifetime of their massive companions, HMXBs are typically
associated with young stellar regions of a galaxy which implies a connection
to the SFR. It was shown that HMXBs are a robust independent tracer for recent
star formation activity in their host galaxy (Grimm et al., 2003; Mineo et
al., 2012a). In particular, the XLF normalization of HMXBs depends linearly on
the SFR. A standard HMXB XLF which we will use throughout the paper is given
in Mineo et al. (2012a) as

$\dfrac{dN_{\mathrm{HXB}}}{dL}=\dfrac{A_{\mathrm{HXB}}}{M_{\odot}\,\mathrm{yr}^{-1}}\,\cases{L}^{-\gamma}\,,&L<L_{cut}\\\
0\,,L\geq L_{cut}\\\ {}\,,$ (3)

where $A_{\mathrm{HXB}}$ is the normalization constant in units of
$M_{\odot}\,\mathrm{yr}^{-1}$, $\gamma=1.59$ is the power-law slope, and
$L_{cut}=10^{3}$ is the cut-off luminosity. The luminosity $L$ is given in
units of $10^{38}\,\mathrm{erg\,s^{-1}}$. One can obtain the number of
individual HMXBs in a galaxy with

$N_{\mathrm{HXB}}=\mathrm{SFR}\times\int_{L_{\mathrm{min}}}^{L_{cut}}\dfrac{dN_{\mathrm{HXB}}}{dL}dL\,,$
(4)

The main difficulty in transferring the SFR scaling property to SSPs in a
cosmological simulation lies in determining the SFH of these SSPs. Since SSPs
are created in an instantaneous star-formation event in the simulation
(Springel et al., 2005), their SFH resembles a $\delta$-function at their time
of creation (Tornatore et al., 2007). In order to circumvent this issue, we
can adapt the procedure in Mineo et al. (2012a) by obtaining the fraction
$f_{X}$ of compact objects becoming a XRB, which is derived from the birth
rate of massive stars ($>8\,M_{\odot}$) in a star-formation event (SFE).
Instead of using the birthrate $\dot{N}_{\mathrm{CO}}$ of compact objects
following a SFE, which is directly dependent on the IMF and SFR, we can
equally use the supernova type II (SNII) rate $R_{\mathrm{SNII}}$, which hides
the SFR dependence in the stellar lifetime function. Following Tornatore et
al. (2007), we can express the SNII rate as

$R_{\mathrm{SNII}}(\tau)=\Phi(m(\tau))\times\dfrac{dm(\tau)}{d\tau}\,,$ (5)

where $\Phi(m)$ is a mass normalized power-law IMF, which gives the number of
stars given the mass $m$, and $m(\tau)$ is the inverse stellar lifetime
function of Padovani & Matteucci (1993) (PM93), which gives the mass of stars
dying at an age of $\tau$. The PM93 lifetime function has the form

$\tau(m)=\cases{1}0^{\left(1.338-\sqrt{1.79-0.2232\cdot(7.764-\log
m})\right)/0.1116-9}\,,&m\leq 6.6M_{\odot}\\\ 1.2\cdot
m^{-1.85}+0.003\,,m>6.6{}.$ (6)

The units of $R_{\mathrm{SNII}}$ are then given in $\mathrm{Gyr}^{-1}$ for a
SSP of $1\,M_{\odot}$. We will employ the modified Chabrier IMF (Chabrier,
2003) used in the simulation which has the form

$\Phi(m)\propto\cases{m}^{-2.3}\,,&1.0\,M_{\odot}<m<100\,M_{\odot}\\\
m^{-1.8}\,,0.5\,M_{\odot}<m<1.0\,M_{\odot}\\\
m^{-1.2}\,,0.1\,M_{\odot}<m<0.5\,M_{\odot}\\\ {}\,,$ (7)

Using the SNII rate, we follow the calculations in Mineo et al. (2012a) to
obtain an estimate of the HMXB population size.

First we determine the product of X-ray bright fraction $f_{X}$ and average
bright-phase duration $\bar{\tau}_{X}$ assuming a XLF

$f_{X}\bar{\tau}_{X}\sim\dfrac{N_{\mathrm{HXB}}(>\\!L_{\mathrm{min}})}{\dot{N}(>\\!8\,M_{\odot})}\,,$
(8)

where $N_{\mathrm{HXB}}$ is proportional to the SFR according to eq. (4),
$\bar{\tau}_{X}\sim 0.1\,\mathrm{Myr}$ is the average duration of the X-ray
bright phase in HMXBs, and $\dot{N}(>\\!8\,M_{\odot})$ is the birthrate of
massive stars which can be obtained from the IMF to first order using

$\dot{N}(>\\!8\,M_{\odot})=\dfrac{\int_{8}^{M_{\mathrm{u}}}\Phi(m)\,dm}{\int_{M_{\mathrm{l}}}^{M_{\mathrm{u}}}\Phi(m)\,m\,dm}\times\mathrm{SFR}\,,$
(9)

with integration limits defined by the IMF. The SFR dependence then disappears
for $f_{X}\bar{\tau}_{X}$. Having computed $f_{X}\bar{\tau}_{X}$ from eq. (8),
we can replace the birthrate $\dot{N}$ with the SNII rate $R_{\mathrm{SNII}}$,
given the age of the SSP, and multiply with the current SSP mass to get the
expected HMXB population size

$\tilde{N}_{\mathrm{HXB}}=\cases{M}_{\mathrm{SSP}}\times
R_{\mathrm{SNII}}(t_{\mathrm{SSP}})f_{X}\bar{\tau}_{X}\,,&t_{\mathrm{SSP}}\leq\tau(8\,M_{\odot})\\\
0\,,t_{\mathrm{SSP}}>\tau(8\,M_{\odot}){}.$ (10)

Additionally we require $t_{\mathrm{SSP}}$ to be smaller than the lifetime of
stars with mass $8\,M_{\odot}$ which according to eq. (6) gives $\approx
30\,\mathrm{Myr}$.

The same procedure can be repeated for different IMFs and XLF models.

In order to study the metallicity dependence of HMXBs, we also include the
model proposed by Lehmer et al. (2021) (L21)

$\dfrac{dN_{\mathrm{HZB}}}{dL}=A_{\mathrm{HZB}}\,\mathrm{SFR}\,\exp[-L/L_{c}(Z)]\cases{L}^{-\gamma_{1}}\,,&L<L_{b}\\\
L_{b}^{\gamma_{2}(z)-\gamma_{1}}L^{\gamma_{2}(Z)}\,,L\geq L_{b}\\\ {}\,,$ (11)

with

$\gamma_{2}(Z)=\gamma_{2,\odot}+\dfrac{d\gamma_{2}}{d\log
Z}[12+\log(O/H)-8.69]\,,$ (12)

and

$\log L_{c}(Z)=\log L_{c,\odot}+\dfrac{d\log L_{c}}{d\log
Z}[12+\log(O/H)-8.69]\,.$ (13)

Both $\gamma_{2,\odot}=1.16$ and $\log L_{c,\odot}=1.98$ are reference values
at solar metallicity. The other parameter values are the normalization
$A_{\mathrm{HZB}}=1.29$, break luminosity $\log L_{b}=0.54$, slope
$\gamma_{1}=1.74$ and first order metallicity corrections
$\frac{d\gamma_{2}}{d\log Z}=1.34$ and $\frac{d\log L_{c}}{d\log Z}=0.6$.
Luminosities are given in units of $10^{38}\,\mathrm{erg\,s^{-1}}$ in eq.
(11). We labeled the L21 model as HZB in order to indicate the metallicity
dependence $Z$.

##  4 Synthetic X-ray emission

In this section we outline the procedure to obtain synthetic X-ray emission of
the hot ISM, AGNs and XRBs from the simulations, where we focus on the XRB
emission model.

###  4.1 PHOX X-ray photon simulator

In this work we generate synthetic X-ray observations of the hot gas, actively
accreting SMBHs and the stellar component in the simulations, using the Phox
code (Biffi et al., 2012, 2013), for further details). In summary, the X-ray
photon simulator Phox consists of three individual modules:

  * •

unit 1: for every emitting source in the simulation, an idealized emission is
computed from a model spectrum, which is statistically sampled by a discrete
number of photons. The generated photon data is stored such that it reflects
the position and velocity information of the emitting sources in the
simulation.

  * •

unit 2: the photon data is projected along a direction through the simulation
and energies are Doppler shifted according to the line-of-sight velocity of
the corresponding emitting source. Additionally, a spatial selection can be
considered.

  * •

unit 3: the produced photon list is convolved with the specific response of a
chosen X-ray telescope, with realistic observing time and detector area.

The Phox code has been applied to simulations of galaxy clusters to study the
properties of the hot diffuse ICM (Biffi et al., 2012, 2013, 2014; Biffi &
Valdarnini, 2015; Cui et al., 2016) and its contamination from AGN emission
(Biffi et al. (2018)) by simulating the X-ray emission of the hot gas (thermal
Bremsstrahlung with metal emission lines) and actively accreting SMBHs
(intrinsically absorbed power-law) in the simulated clusters. Additionally,
the general approach of the code allows for the treatment of different
potential X-ray sources within the simulation. The corresponding emission
model, however, must be constrained from source properties within the
simulation and included in unit 1.

###  4.2 Hot gas emission model

The synthetic X-ray emission for each hot-phase gas element in the simulations
is modelled based on their intrinsic thermal and chemical properties (density,
temperature, metallicity). A single-temperature VAPEC thermal emission model
with heavy element emission lines (Smith et al., 2001) is then assumed for
every gas element making use of the implementation in XSPEC (Arnaud, 1996).
For further details on the implementation, we refer the reader to the
description provided in Biffi et al. (2012, 2013).

###  4.3 AGN emission model

Similarly, the synthetic X-ray emission for accreting BHs in the simulations
is modelled by determining their accretion rate w.r.t. their Eddington
accretion rate. In order to represent the effect of dense torus absorption,
the X-ray spectrum is modelled using an intrinsically absorbed power-law with
column densities that are stochastically selected from the distribution given
in Buchner et al. (2014). For further details on the implementation, we refer
the reader to the description provided in Biffi et al. (2018). The code has
since been updated to model the AGN intrinsic absorption using the TBABS
photo-ionization cross-sections (Wilms et al., 2000) instead of the WABS
cross-sections (Morrison & McCammon, 1983), both implemented in XSPEC.

###  4.4 XRB emission model

Following the modelling of AGN and hot gas emission in the Phox code, we
similarly determine the X-ray emission coming from a stellar resolution
element of the simulation. We make use of the SSP properties determined in
section 3 to calculate an expected number of point sources $\tilde{N}$ per
stellar element. We then sample the respective XLF based on the following
recipe. We select the correct XLF from a stellar age criterion: if the age of
the SSP is $\tau<30$ Myr, we classify it as HMXB hosting. If $\tau>1$ Gyr we
classify the SSP as LMXB hosting. From the XLF we construct a pseudo
cumulative density function (pCDF). This requires setting a luminosity
interval $[L_{\mathrm{min}},L_{\mathrm{max}}]$ at which the pCDF $C(L)$ will
be defined

$C(L)=1-\dfrac{N(>L)}{N(>L_{\mathrm{min}})}\,,$ (14)

with $N(>\\!L)$ being the integrated XLF from step 1, and $C(L)\in[0,1]$.
Typically, $L_{\mathrm{max}}$ is chosen to be the cut-off luminosity of the
respective XLF. Note that $C(L)$ is sensitive to the choice of
$L_{\mathrm{min}}$ due to the power-law dependency of the XLF. We then
calculate $\tilde{N}$ from eq. (2),(10) and choose $\lfloor\tilde{N}\rfloor$
uniformly distributed random numbers $p_{i}\in[0,1]$. If a separate uniform
random number $\tilde{p}$ satisfies
$\tilde{p}<\tilde{N}-\lfloor\tilde{N}\rfloor$, we draw an additional random
number $p_{i}$. For each $p_{i}$ we determine a corresponding luminosity
$L_{i}$, at which $C(L_{i})=p_{i}$. From this we can calculate the total
luminosity of the SSP as

$L_{\mathrm{SSP}}=\sum_{i}L_{i}\,.$ (15)

This recipe yields a SSP with a single X-ray luminosity in the same energy
range the XLFs have been defined in. As such, it will be treated as a single
XRB-like point source. Specifically, we describe XRB emission spectra as
redshift dependent, absorbed power-laws of the form

$A(E)=w(E(1+z))\times K\left[E\left(1+z\right)\right]^{-\Gamma}\,,$ (16)

where $K$ is the spectrum normalization at $1\,\mathrm{keV}$ and $\Gamma$ is
the photon index as given in the XSPEC manual (Arnaud (1996)). Absorption will
be modelled after the Wilms et al. (2000) photo-ionization cross-sections
including gas, grain and molecule components in the ISM which is expressed in
eq. (16) with $w(E)$. The redshift dependence of the spectrum allows for
consistent modelling in a cosmological context. The total observed luminosity
of the SSP (eq. (15)) in a certain energy band $E\in[E_{1},E_{2}]$ is then
related to eq. (16) by:

$L_{\mathrm{SSP}}=F\int^{E_{2}}_{E_{1}}A(E)EdE\,,$ (17)

with $F=4\pi\mathcal{D}_{L}(z)^{2}\times 1.602\cdot
10^{-9}\mathrm{erg\,keV^{-1}}$ being the rescaling factor between flux and
luminosity and $\mathcal{D}_{L}(z)$ the luminosity distance inferred from the
underlying cosmology. The energy range $[E_{1},E_{2}]$ has to be chosen in
correspondence to the range at which the XLFs were defined which is commonly
adopted as $E_{1}=0.5\,\mathrm{keV}$ and $E_{2}=8\,\mathrm{keV}$ in the
observed frame. We are then able to constrain the spectrum normalization
$K_{\mathrm{SSP}}$ using eq. (15) and (17) where we get

$K_{\mathrm{SSP}}=\dfrac{\sum_{i}L_{i}}{F\int_{E_{1}}^{E_{2}}w(E)[E(1+z)]^{-\Gamma+1}dE}\,.$
(18)

A fully detailed description of XRB spectra is not possible in our modelling.
For instance, the cosmological context of our approach makes it impossible to
capture XRB variability on small time-scales such that changes in spectral
hardness have to be accounted for statistically. Therefore, we opt for an
average X-ray emission spectrum for the two XRB types, in line with typical
values found in the literature. For our analysis, we chose the power-law slope
$\Gamma_{\mathrm{LXB}}=1.7$ in the LMXB case and $\Gamma_{\mathrm{HXB}}=2$ in
the HMXB case. Both slopes are motivated from observational data where the
$\Gamma_{\mathrm{HXB}}$ follows assumptions of Mineo et al. (2012a) and
$\Gamma_{\mathrm{LXB}}$ follows the median slope of all high-confidence XRBs
in Lehmer et al. (2019). For the absorption component we assume the median
column density $N^{\mathrm{xrb}}_{H}=2\cdot 10^{21}\,\mathrm{cm}^{-3}$ from
the same sample of XRBs in Lehmer et al. (2019). This choice of parameters
results in spectra resembling typical XRBs in their low/hard-state (see
Remillard & McClintock, 2006; Done et al., 2007; Sazonov & Khabibullin,
2017b). Additionally, it was found that $\Gamma_{\mathrm{LXB}}\sim 1.4-1.8$ is
sufficient to describe the hot atmosphere of elliptical galaxies (Boroson et
al., 2011; Wong et al., 2014; Babyk et al., 2018). We do not include any other
form of stellar X-ray sources such as AB stars or cataclysmic variables since
their cumulative luminosity is at least an order of magnitude below the total
XRB luminosity (see e.g. Boroson et al., 2011; Babyk et al., 2018).

Having determined $K$, $\Gamma$ and $N^{\mathrm{xrb}}_{H}$ of a single SSP, we
can compute the resulting XRB spectrum by adopting the ztbabs and zpowerlw
model embedded in XSPEC, from which we calculate an ideal photon list
associated to each SSP.

##  5 Simulated data set

For our analysis we extracted galaxy sized halos at redshift $z=0.07$ from the
Box4/uhr cosmological volume of the Magneticum Pathfinder Simulation set and
artificially placed them at $z=0.01$. Assuming fiducial values for effective
area $A_{\mathrm{eff}}=1000\,\mathrm{cm}^{2}$ and exposure time
$T_{\mathrm{exp}}=10^{5}\,s$, we select all halos generated by the subfind
algorithm (Springel et al., 2001; Dolag et al., 2009) which have a total
stellar mass $M_{*}>10^{9.7}\,M_{\odot}$ within a sphere of $R_{2500}$ around
the halo center, which also includes few massive group-like objects with
stellar masses of $M_{*}\sim 10^{12}\,M_{\odot}$. Halos with stellar masses
below a threshold of $M_{*}<10^{10}\,M_{\odot}$ are resolution limited: the
SPH implementation fails to correctly reproduce the required star-formation
for low-mass halos. We illustrate the $M_{\mathrm{gas}}-M_{*}$ relation in
figure 1. Halos which have high stellar masses and low gas masses were
subjected to feedback processes and were depleted of their gas content which
is connected to resolved physical processes within the simulations.

Figure 1: $M_{\mathrm{gas}}-M_{*}$ relation of our simulated galaxy sample.
The dashed diagonal line indicates a one-to-one ratio.

The selected set consists of $1480$ objects of which $335$ are considered
actively star-forming, based on the limit
$\mathrm{SFR}>0.4\,M_{\odot}\,\mathrm{yr}^{-1}$, with SFR derived from the
mass of stars born in the past 100 Myrs within the simulation. X-ray
luminosities were calculated from the photon list obtained from a cylinder
around each galaxy with radius $R_{2500}$ and projected along the l.o.s. for a
length of $2R_{2500}$. The l.o.s. direction coincides with the $z$-axis of the
simulation box. We chose $R_{2500}$ as the smallest available scale invariant
size from the halo catalogue. We assume a weak galactic foreground absorption
of $N_{H}=10^{20}\,\mathrm{cm}^{-3}$ for the photon generation. Because Phox
was originally conceived as a tool to study the X-ray emission of galaxy
clusters, the calculated emission for the hot gas component in galaxies does
not account for self-absorption from the ISM. We attempt to model self-
absorption following observational derivations of intrinsic ISM emission by
employing an additional TBABS model with $N^{gas}_{H}=5\cdot
10^{21}\,\mathrm{cm}$ at the source redshift (see e.g. Mineo et al., 2012b;
Gilbertson et al., 2022; Lehmer et al., 2022). In figure 2 we show our
galactic sample on the SFR-$M_{*}$ plane color coded by the $b$-value which is
a measure of galaxy morphology from the intercept of $M_{*}-j_{*}$ relation
(Romanowsky & Fall, 2012; Teklu et al., 2015). We adopt the classification
scheme by Schulze et al. (2020) with $b>-4.35$ for disk galaxies (blue),
$-4.73<b<-4.35$ for intermediates (green) and $b<-4.73$ for spheroids (red).
In order to increase the sample size of low-metallicity star-forming galaxies
we follow the approach taken by Weinmann et al. (2010) where we assign a
uniformly distributed $\log\,\mathrm{sSFR}$ in the range of -12.4 and -11.6
for halos with $\mathrm{SFR}=0$. This can be seen in figure 2 as an apparent
strip of mostly spheroidal galaxies located around the
$\log\,\mathrm{sSFR}=-12$ line. We emphasize that galaxies distributed
according to this criterion will not be used to derive hot gas properties in
star-forming halos since the assigned SFR is not connected to gas properties
within the simulations.

Figure 2: $\mathrm{SFR}-M_{*}$ relation of our sample. The dashed diagonal
lines correspond to constant $\log\mathrm{sSFR}$ of -9 (top) to -12 (bottom)
in increments of 1. Each galaxy was color-coded by its respective $b$-value
(Teklu et al., 2015) with classification scheme according to Schulze et al.
(2020) (red: elliptical, green: intermediate, blue: disk). Galaxies clustered
at $\log\mathrm{sSFR}=-12$ were given a SFR estimate based on the Weinmann et
al. (2010) approach (see text).

Subsequent analysis of the generated ideal photons will concentrate on the
unfolded photon lists without taking into account any instrumental response.

##  6 Results

###  6.1 Spatial distribution of XRB emission

As a first step we show that the XRB seeding follows the correct spatial
distribution of XRB emission. We select a poster-child late-type galaxy from
our sample with $M_{*}=10^{11.3}\,M_{\odot}$,
$\mathrm{SFR}=8.4\,M_{\odot}\,\mathrm{yr}^{-1}$ and
$R_{2500}=91\,\mathrm{kpc}$ which appears face on to the line-of-sight. We
focus our showcase to the inner 60 kpc in order to better capture the stellar
body of the galaxy. In figure 3 we show the percentile contours of photons
emitted by HXBs (black) and LXBs (red) on top of the line-of-sight projected
stellar mass map of the galaxy. As expected, the LXB emission spans the whole
field of the galaxy with its center over the galaxy bulge and decreasing with
distance from the center reflecting the stellar surface density
($\Sigma_{*}$). In contrast, HXB emission is confined to the stellar mass
overdensities surrounding the galaxy center. The reason for this can be seen
in figure 4 where we show the same X-ray contours on top of the mass-weighted
stellar age map of the galaxy. It shows that HXB emission is bound to much
younger stellar fields of the galaxy. In particular, the stellar age map shows
the complex intersection of the two XRB components in the galactic field (e.g.
Lehmer et al., 2017; Gilbertson et al., 2022). While this general behaviour
for the two XRB types is expected, the ability to clearly distinguish the
origin of X-ray emission of a galaxy is a major advantage and enables
predictions on theoretical aspects of X-ray spectra of unresolved galaxies in
high redshift observations. We note, that the spatial distribution of emitted
photons was smeared out according to a Gaussian positioned at the projected
particle location with standard deviation equal to the simulation intrinsic
smoothing length of the SSP. Because of that, XRBs are not visible as point
sources.

Figure 3: Projected stellar surface mass density ($\Sigma_{*})$ of poster-
child disk galaxy in Box4 of Magneticum with halo ID 13633. Contours indicate
the 68, 90 and 99 percentiles of X-ray photon position for HXBs (black) and
LXBs (red).  Figure 4: Mass weighted stellar age ($\tau_{*}$) of poster-child
disk galaxy in Box4 of Magneticum with halo ID 13633. Contours indicate the
68, 90 and 99 percentiles of X-ray photon position for HXBs (black) and LXBs
(red).

###  6.2 XLFs

For the XLF reconstruction we chose a lower luminosity limit of
$L_{\mathrm{min}}=10^{35}\,\mathrm{erg\,s^{-1}}$, since all XLFs are expected
to have a turn-over at luminosities $\lesssim 10^{35}\,\mathrm{erg\,s^{-1}}$.
This is due to neutron star XRBs generally experiencing lower Eddington
accretion rates than their BH counterparts, as well as a strong magnetic field
influencing the accretion flow (Shtykovskiy & Gilfanov, 2005a, b). In
addition, the integrated luminosity from the XLF will be dominated by the high
luminosity end such that lower value choices of $L_{\mathrm{min}}$ won’t
change the integrated luminosity significantly. In figure 5 we compare the
integrated XLFs (colored), reconstructed from the approach outlined in section
4.4, with their original model (black). Corresponding XLFs are indicated with
the same linestyle. In detail, each coloured line is derived from the total
photon box produced for each SSP in the simulation. The apparent turn-over of
the reconstructed lines is unrelated to any physical effects associated with
NS XRBs but can be explained by the following reasons: 1\. low luminosity XRBs
are unlikely to produce observable photons given the distance of the
simulation box (low flux); 2. in cases where the XRB population size is
greater than one per SSP, low luminosity XRBs are removed from the XLF through
the summation process of eq. (15). The corresponding luminosity for a XRB
emitting a single photon in our setup, assuming an absorbed power law
spectrum, is marked by the vertical grey band at $L_{X}\sim
10^{37}\,\mathrm{erg\,s^{-1}}$ which roughly coincides with the turnover given
the energy range of $0.5-8$ keV. The solid vertical grey line is the
luminosity of a photon with mean energy $\sim 3\,\mathrm{keV}$ according to
eq. (16). Thus, SSPs with $L_{\mathrm{SSP}}<10^{37.4}\,\mathrm{erg\,s^{-1}}$
will produce at most one single photon given the discrete photon sampling in
Phox. Note, that our setup is not dependent on background flux calibration,
since our goal in Unit 1 of Phox is to create an ideal photon list. We find
excellent agreement between the reconstructed XLFs (colored lines) and their
original model (black, same linestyle) above the single photon luminosity
limit. The average reconstructed XLF for the L21 model (orange) is consistent
with solar values for $12+\log[\mathrm{O/H}]$ which is almost
indistinguishable from the Mineo et al. (2012a) XLF. The dotted and dash-dot-
dotted lines indicate the L21 model at $12+\log[\mathrm{O/H}]=9.2$ and 7.2
respectively which are not a good description of the average HZB XLF.

Figure 5: Integrated XLFs from eqs. (1),(3),(11) normalized by their
respective scaling (SFR and $M_{*}$). The dash-dot-dotted line is the XLFs of
the L21 model, for $\log[\mathrm{O/H}]+12=7.2$ and the dotted line for
$\log[\mathrm{O/H}]+12=9.2$. Colored lines are average XLFs obtained from the
luminosity sampling in Phox: Mineo et al. (2012a) (M12) dashed blue, Lehmer et
al. (2021) (L21) dash-dotted orange, Zhang et al. (2012) (Z12) solid green.
Corresponding lines share the same linestyle. The grey vertical line at $\log
L\sim 37$ marks the maximum luminosity of a source emitting one single photon
with average energy 0.5-8 keV.

###  6.3 Scaling Relations

Since our approach to assign luminosities to SSPs is agnostic to its
environment, i.e if the SSP is part of a star-forming halo within the
simulation, we will first verify the self-consistent emergence of well known
galactic X-ray scaling relations.

Studies of the diffuse galactic X-ray emission revealed linear relationship
with SFR in the soft X-ray (SXR) band which is associated to ISM heating
through SN remnants and stellar wind from heavy stars (Ranalli et al., 2003;
Gilfanov et al., 2004a; Mineo et al., 2012b; Bogdán et al., 2013). However, it
is observationally challenging to separate pure ISM emission from unresolved
emission from point sources. In fact, Mineo et al. (2012b) and Anderson et al.
(2013) show that 2/3 of the galactic extended emission in the soft-band can be
attributed to XRB emission which they estimated from XRB scaling relations.
Since Phox was designed to reproduce X-ray emission from galaxy clusters, we
have to verify that galactic ISM properties are appropriately reproduced. In
fig. 6 we show the relationship between the hot gas X-ray luminosity and SFR
in the soft X-ray band (SXB) (0.5-2 keV). The blue contours enclose the 68, 90
and 99 percentiles of galaxies in our simulated sample also indicated by
underlying green dots. We only show galaxies with $\log\mathrm{sSFR}>-11.6$
which removes all galaxies falling in the Weinmann et al. (2010) criterion.
This is done to reduce the influence of unresolved galaxies on the relation as
well as removing weakly star-forming galaxies. The solid black line indicates
the best-fit linear relation from Lehmer et al. (2022) with
$\log(L_{X}^{0.5-2\,\mathrm{keV}}/\mathrm{SFR})=39.58^{+0.17}_{-0.28}$. We
include absorbed ISM luminosities from Mineo et al. (2012b) (red squares).
Their reported values in the SXB are smaller compared to the L22 relation.
They argue however, that their normalization might be higher by at least a
factor of 2, due to uncertainties in observed column densities and depending
on their fitting procedure. Additional bolometric gas luminosities are
obtained from Bogdán et al. (2013) (B13) (black pentagons) in the 0.5-2 keV
energy range and a radial range of 0.05-0.15 $R_{vir}$ for late-type galaxies.
Although B13 aimed to study extended X-ray coronae around large spiral
galaxies, the measured luminosities are consistent with a linear
$L_{X}-\mathrm{SFR}$ relation prompting a connection between SFR and galactic
outflows. From Strickland et al. (2004) (S04) we obtained data points for
absorption corrected total gas luminosities in LTGs in the 0.3-2 keV energy
range (black stars) which are also showing a linear SFR relation. Assuming the
same energy range of 0.5-2 keV and no absorption correction for S04 values,
resulting luminosities would be slightly lower. Within the 90 percentile
contours, our galaxy sample is consistent with observational data of the
extended galactic emission. We notice however, that while a linear trend
exists in our data, the relation is less tight and contaminated by low-
luminosity galaxies at $\mathrm{SFR\sim 1}$. Low-luminosity galaxies may have
been the result of feedback processes removing large fractions of gas from the
central halo, which we discuss in more detail in section 8. The intrinsic
scatter in our data of $>1\,\mathrm{dex}$ is concerning, but might be a result
of an uncurated sample. An important caveat is the limited resolution of the
underlying simulation. Halos with stellar mass below the milky way stellar
mass suffer greatly from the resolution restriction. In those cases, the total
gas and stellar halo are comprised by only a few hundred resolution elements
giving rise to a large intrinsic scatter due to the SPH implementation of the
simulation as well as an eventual cut-off where a halo was not able to form
stars (see figure 1).

Figure 6: Relation between hot gas X-ray luminosity (0.5-2 keV) and SFR, for
star forming galaxies in the simulations. We only show galaxies with
$\log(\mathrm{sSFR})>-11.6$ which removes all galaxies falling in the Weinmann
et al. (2010) criterion. This is done to reduce the influence of unresolved
galaxies on the relation as well as removing weakly star-forming galaxies. The
solid black line indicates the linear relation from Lehmer et al. (2022). The
gray shaded area is the corresponding 1-$\sigma$ standard deviation. Data
points are bolometric luminosities of the diffuse emission taken from
Strickland et al. (2004) (SO4) in the 0.3-2 keV band, Mineo et al. (2012b)
(M12b) in the 0.5-2 keV band and Bogdán et al. (2013) (B13) in the 0.5-2 keV
band.

The primary goal of this paper is to model the stellar X-ray emission in the
form of XRBs using data from cosmological simulations. In figure 7 and figure
8 we show the principal scaling relations for HXBs and LXBs respectively. The
contours have the same meaning as in figure 6. The solid black line indicates
the expected scaling relation from integrating the luminosity functions in eq.
(3) and (1). In the HXB case in figure 7 the reconstructed scaling relation is
consistent with observations from Mineo et al. (2012a) but shows significant
deviation from the expected linear SFR dependence at low SFR values. This
apparent deviation was noticed in observations of star-forming galaxies as
well (Grimm et al., 2003; Gilfanov et al., 2004a) and explained in Gilfanov et
al. (2004b) where the transition between linear and power-law regime of the
SFR relation is caused by low-number sampling of a single sloped XLF. The
approximate relation calculated from Gilfanov et al. (2004b) for the XLF in
eq. (3) is shown as the dashed black line and has a slope
$\frac{1}{\gamma-1}$. For the LXB case in figure 8, the linear stellar mass
dependence is well reproduced and does not suffer from low-number sampling,
since the corresponding XLF from eq. (1) is not single sloped with suppressed
contribution from low luminosity LXBs. Additionally, our stellar mass limit
($M_{*}=10^{9.7}\,M_{\odot}$) is relatively high compared to the XLF
normalization (eq. (1)) which allows for sufficient sampling points. In the
stellar mass regime of dwarf galaxies ($M_{*}\lesssim 10^{9}$) we expect
breaks in the $L_{X}-M_{*}$ relation as well (Gilfanov et al., 2004b).

Figure 7: Same as 6, for HXB luminosity in the 0.5-8 keV band. The filled blue
contours represent the same halos with HXBs sampled from the L21 model (eq.
(11)). The solid black line is the expected relation from the employed HXB
model in eq. (3), the dashed red line is the predicted relation for discretely
sampling data points from a power-law distribution Gilfanov et al. (2004b).
The solid blue line is expected relation from the employed L21 model. Data
points are taken from Mineo et al. (2012a).  Figure 8: Same as 7, for LXB
emission in the 0.5-8 keV band. The solid black line is the expected 1:1
relation from the employed HXB model in eq. (1). Data points (cyan triangles)
are taken from Lehmer et al. (2020).

When combining the SFR normalized emission from HXBs and LXBs we expect a sSFR
dependent relation of the form

$L_{X}/\mathrm{SFR}=\alpha\,\mathrm{sSFR}^{-1}+\beta$ (19)

where $\alpha$ contains the linear mass dependence of the total LXB luminosity
and $\beta$ contains the linear SFR dependence of the total HXB luminosity
(Lehmer et al., 2010, 2016, 2019). In figure 9 we show the sSFR dependence of
the SFR normalized combined XRB luminosity of our sample with contours as in
figure 6. Additionally, we show data points from Mineo et al. (2014) (M14)
from nearby galaxies and a study performed on star-forming galaxies within the
Virgo cluster by Soria et al. (2022) (S22) with SFR values derived from the
$12\,\mu m$ WISE W3 measurement. We include an observationally derived
relation (black, dash-dotted) for local galaxies from Lehmer et al. (2019)
(L19). The red dashed diagonal line represents the theoretical contribution
from LXBs given our model choice (fig. 8) and the blue dashed line corresponds
to contributions from HXBs (fig. 7). Because of the HXB scaling relation being
more similar to a broken power law, the combined luminosity of LXBs and HXBs,
normalized by SFR, takes a different shape (dashed purple) compared to the L19
relation: it reaches a minimum value at $\mathrm{sSFR}\sim-10$ after which it
begins to approach the linear HXB regime. Our sample is in agreement with the
purple relation and is broadly consistent with the relations from L19 for
$\log\mathrm{sSFR}\lesssim-10.5$. However, due to the lack of a high sSFR
halos in the simulated sample, we are not able to reliably demonstrate the
expected turnover predicted by our model in purple. Interestingly, we find
remarkable agreement between the simulated data and S22 data obtained from
star-forming galaxies in the Virgo cluster. Furthermore, data from M14
connects seamlessly to both the S22 sample and our sample. From a visual
inspection alone, the data from S22 appears to be more consistent with a
broken SFR-scaling relation. In their study S22 show that most of their data
points with $\log\mathrm{sSFR}$ lie below the median tracks of Monte-Carlo
simulations of L19 for various populations of galaxies. They argue that
inclination effects and differences in the used galaxy population are causing
this discrepancy.

Figure 9: Relation between SFR normalized total X-ray luminosity (both XRB
types) and specific SFR, for star-forming galaxies in Box4 of Magneticum
(contours: percentiles). The solid black line is the relation of observed
normal galaxies at $z\sim 0$ from Lehmer et al. (2016). The dash-dotted line
represents the relation for a set of local galaxies from Lehmer et al. (2019).
The dashed red line depicts the expected scaling relation for a pure LXB
contribution given our XLF choice. The dashed blue line shows the expected
broken SFR scaling relation of a pure HXB contribution. Data points are taken
from Soria et al. (2022) (S22) and Mineo et al. (2014) (M14).

###  6.4 Metallicity relation of HXB

Numerous studies on the total HXB emission in individual galaxies associate an
increase of HXB emissivity per unit SFR with a lower metallicity environment
of the HXB. (Vulic et al. (2021); Saxena et al. (2021); Lehmer et al. (2021);
Garofali et al. (2020); Fragos et al. (2013b)). In this section we try to
verify if our modelling allows us to capture the metallicity dependence of the
HXB luminosity function proposed by Lehmer et al. (2021) (L21) (see eq. (11)).
In figure 10 we show the SFR normalized X-ray luminosity of our star-forming
galaxy sample against the metallicity derived from the mass weighted oxygen
fraction of the young stellar population ($<30$ Myr). The solid black line
indicates the global model from L21 in the limit of perfect LF sampling
together with its 16-84 percentiles in gray. We binned our sample according to
metallicity and SFR with horizontal bars representing the bin-width in
metallicity, vertical bars representing the 16-84 percentile in luminosity.
Black crosses are the median values for the whole star-forming sample while
colored symbols correspond to different SFR bins. The SFR bins are displaced
w.r.t the metallicity bin center for illustrative purposes. Our full sample
(black crosses with errorbars) is consistent with the theoretical model and
shows the expected increase in luminosity for lower metallicities. Similarly,
the high SFR sample ($\gtrsim 5\,M_{\odot}\,\mathrm{yr}^{-1}$, blueish colors)
is also consistent with the global model up to a oxygen fraction of 8-8.5. Our
sample did not include highly star-forming galaxies with lower metallicities.
For the sample at $\mathrm{SFR}\sim 1\,M_{\odot}\,\mathrm{yr}^{-1}$ (green) we
find lower median luminosities compared to the global model but are consistent
within the 25 percentile margin. From the MCMC simulation performed in L21, it
is however expected that the median luminosity decreases with lower SFR values
due to worse HXB count statistics. Despite this, the median luminosity for the
low SFR sample (orange) lies more than an order of magnitude above the
relation predicted in L21 for $\mathrm{SFR}\lesssim 0.1$. Given our approach
to retrieve HXB luminosities for SSPs in the simulation, we find, that with
the given mass resolution of the simulation, a single eligible SSP in a low
SFR galaxy might already overcount the expected number of HXBs. Additionally,
only luminous SSPs will be sampled by Phox (see fig. 5. Since we did not take
into account galaxies which do not have HXB emission despite being star-
forming, we are biased towards a low SFR sample with comparatively large
luminosities and too large HXB population.

Figure 10: Metallicity dependence on the linear SFR scaling relation of HXB
luminosity. Oxygen fractions are derived from mass weighted stellar
metallicities of SSPs younger than 30 Myr. The black solid line is the mean
global model calculated by L21 in the limit of perfect LF sampling. The shaded
grey area is its 16-84 percentile margin. Our data is binned by oxygen
fraction and SFR where SFR bins are indicated by color. Vertical errorbars
correspond to the 25-75 percentile margin of the data within each metallicity
bin, horizontal black bars indicate the width of the metallicity bin. Colored
symbols corresponding to the SFR bins are displaced horizontally w.r.t the
metallicity bin center for illustrative purposes. The purple dashed line is
the corresponding M12a model where no explicit metallicity dependence is
assumed.

##  7 Relative emission of hot gas and XRBs

In this section, we use the ideal photon list generated with Phox to obtain
relative count ratios between the gas and XRB component in our galaxy sample.
Additionally, we build average spectra of normal galaxies by first looking at
the whole sample and subsequently dividing the sample according to different
galactic properties. Average spectra are obtained by first normalizing each
galaxy spectrum with its total emission in the 0.5-10 keV band and afterwards
taking the average in each energy bin of the stacked spectra.

###  7.1 Relative count ratio

In order to investigate the relative contribution of the hot gas and XRB
component in our galaxy sample, we compare the respective photon counts coming
from each component. We split the contribution in different energy bands
(0.5-2 keV, 2-8 keV, 0.5-8 keV) and bin our galaxies according to the count
ratio $r$ which is defined as

$r=1+\dfrac{c_{XRB}}{c_{GAS}}\,,$ (20)

where $c_{i}$ stands for the number of photon counts from the respective
component. A value of $r>2$ indicates a dominant XRB component in our setup.
In figure 11 we show the differential and cumulative fraction of galaxies in
our sample w.r.t. the ratio $r$. We show the LXB contribution in lime green,
the HXB contribution in dark green and the summed contribution of LXBs and
HXBs in orange. Different line styles correspond to the different energy bands
and are indicated in the figure legend. Regarding the combined contribution of
both XRB types we notice that $\sim 90\%$ of our galaxies are dominated by XRB
emission in the hard X-ray band (2-8 keV). Interestingly, this high fraction
is almost exclusively caused by the LXB emission. One has to keep in mind,
however, that HXBs only occur in galaxies with sufficient star-formation.
Thus, the low HXB contribution is rather an expression of the rareness of
star-forming galaxies compared to non-star-forming ones in our sample. If we
include the soft-band (0.5-2 keV), where the hot gas emission is stronger,
only $\lesssim 40\%$ of our galaxies are still dominated by XRB emission.

Figure 11: X-ray photon count ratio between the XRB and hot gas component
including all galaxies in our sample. Colors represent the component while
line styles represent the energy range in the observed frame. The label XRB
shows the sum of HXB and LXB emission. Top: Differential fraction of galaxies
with a given count ratio $r$. Bottom: Cumulative fraction of galaxies. Here, a
ratio of $r\geq 2$ means that the XRB component is dominant. Ratios exceeding
the given range are collected in the last bin.

###  7.2 Average galaxy spectra

In addition to comparing the total counts received from each component, we
construct an average spectrum from star-forming galaxy sample and quantify the
contribution of each component in each energy bin. To do this we first
normalize each galaxy spectrum to its 0.5-10 keV emission and taking the
average in every energy bin. The results can be seen in figure 12 where the
solid black line is the total average spectrum, the orange line is the average
combined spectrum of both XRB types (green) and the blue line is the average
spectrum of hot gas. In the bottom panel of figure 12 we show the relative
contribution of each component towards the total spectrum as a function of
energy. Additionally, we added the average spectrum by Lehmer et al. (2022)
which was derived for low metallicity dwarf galaxies (dashed red) and values
for the relative XRB contribution in different energy bands at redshift $z=0$
from Lehmer et al. (2016) (L16) as black symbols and errorbars in the lower
panel of figure 12. From a visual inspection, the combined average spectrum is
similar to the L22 spectrum in the energy range 0.8-10 keV despite differences
in mass and metallicity regimes in their sample. Below 0.8 keV our average
spectrum falls off more steeply than the L22 spectrum. This is connected to
differences in assumptions made for the different component spectra in L22.
They use a two temperature APEC model with intrinsic absorption column
densities similar to our $N^{gas}_{H}=5\cdot 10^{22}\,\mathrm{cm^{-3}}$ but
use a weaker absorption coefficient $N^{\mathrm{xrb}}_{H}=6\cdot
10^{20}\,\mathrm{cm^{-3}}$ compared to our $2\cdot 10^{21}\,\mathrm{cm^{-3}}$.
Because of this, the XRB contribution is much higher for low energies. The
sensitive behaviour of the spectrum at low energies can also be seen when
comparing the relative contributions of XRB with values from L16 in the lower
panel. In L16 they build an average spectrum for massive star-forming galaxies
from observations by Wik et al. (2014); Lehmer et al. (2015) and Yukita et al.
(2016). While we predict strong contributions of XRBs in the 0.5-1 keV energy
range with ¿30%, L16 find much lower values of $\sim 10\%$. Our predictions
agree well with L16 for harder energy bands (¿1.5 keV). Interestingly, Wik et
al. (2014) also find high XRB contributions in the 0.5-1 keV band which are
similar to our values again highlighting the strong model dependency of the
XRB contribution in the SXR. Given our simple choice for the XRB emission
model in the form of a single absorbed power law, we fail to capture the
expected steepening of the XRB power law for $E>6\,\mathrm{keV}$ which is
predicted by Lehmer et al. (2016) for normal galaxies, by population synthesis
performed in Fragos et al. (2013b) and by e.g. Persic & Rephaeli (2002) for
star-burst galaxies. A more sophisticated model of the XRB spectrum would be
needed in order to fully reproduce observed spectral shapes.

Figure 12: SFR normalized average X-ray spectrum of our star-forming sample
($\mathrm{SFR}>0.4\,M_{\odot}\,\mathrm{yr}^{-1}$). The combined spectrum
(black) is made up of the HXB (dark green), LXB (lime green) and the hot gas
(blue) component. The combined XRB emission is given in orange. The red line
is the L22 average spectrum derived from observations of dwarf galaxies.

Further disentanglement of the sample into increasing SFR bins yields a more
refined view on the underlying galactic properties influencing the spectral
composition. In figure 13 we show the same decomposition of galactic spectra
as in figure 12 for four SFR bins. While the average XRB contribution stays
consistent with previously discussed values, we see slight variations between
each SFR bin. Specifically, in the lowest SFR bin (top left) we notice a very
strong XRB component which generally exceeds estimates from L16. By comparing
to figure 6, where we found galaxies with $\mathrm{SFR}\lesssim
1\,M_{\odot}\,\mathrm{yr}^{-1}$ having low $L_{X}$ within the 99 percentile
contour, we attribute the strong XRB component to under-luminous gas-poor
galaxies in the probed SFR regime and normalization issues when determining
$\tilde{N}_{\mathrm{HXB}}$ from eq. (10), which leads to more HXBs per unit
SFR. In the next higher SFR bin (top right) the normalization of the spectrum
drops slightly, mainly in the HXB component because the normalization issue no
longer applies here. Simultaneously, the hot gas component increases
normalization in the hard band due to more massive galaxies contributing to
the bin with respectively hotter gaseous atmospheres. The relative
contribution of the XRB component is consistent with L16 values in this case
in the hard band. Further increasing SFR in each following bin yields even
more prominent hot gas contributions in the hard band (see bottom left and
bottom right). While the SFR normalized total XRB luminosity stays almost
constant in the two last bins, the XRB contribution is diminished by strong
contributions from the hot gas atmospheres around massive, BCG-like galaxies
with their deep gravitational potential. In these cases, the L16 values become
less suitable to compare against since the considered samples have vastly
different properties. In the highest SFR bin (bottom right) the inclusion of
group-like mass regimes also disfavors the comparison against the L22 spectrum
which was derived using high redshift analogues for normal galaxies. Since the
two lowest SFR bins present the largest contribution in terms of sample
number, the total average spectrum from figure 12 will mostly resemble to the
upper panels of figure 13. An important question associated with galactic
X-ray spectra is their influence on the cosmic X-ray background. Upcoming
planned X-ray missions like Athena or LEM are sensitive enough to receive
significant contributions from galactic emission to their background
calibration. We plan to investigate this issue in an upcoming study.

Figure 13: Mean X-ray spectrum solid black of our star-forming galaxy sample
binned by SFR from lowest (upper left panel) to highest (bottom right) panel
in the energy range 0.5-10 keV. The mean spectrum was calculated by
normalizing each galaxy spectrum to its total 0.5-10 keV emission and taking
the average in each energy bin. Dotted lines represent the different
components contributing to the total spectrum: hot gas emission (blue), LXB
(lime green) and HXB (dark green). XRB (orange) stands for LXB+HXB emission.
The red line is the L22 average spectrum derived from observations of dwarf
galaxies. The lower part of each panel shows the fractional contribution of
each component to the mean spectrum together with the XRB count ratios from
L16. The grey box highlights the bin size and the number of galaxies $N$
within each bin.

##  8 Discussion

###  8.1 Gas luminosity scaling relation

Compared to observational data we found that the $L_{X}-\mathrm{SFR}$ relation
in the 0.5-2 keV band for the hot gas component has a large intrinsic scatter
of $\sim 1\,\mathrm{dex}$ at $\mathrm{SFR}\lesssim 1$ with observers typically
reporting lower luminosities of $L_{X}\sim 10^{39}\times\mathrm{SFR}$ with
significantly smaller scatter (Strickland et al., 2004; Mineo et al., 2012b;
Bogdán et al., 2013). We propose several explanations as to why we expected
the scatter to be larger for the hot gas component.

  * •

In its core, the Phox code was perceived as a tool to study X-ray emission
from galaxy clusters with relatively hot atmospheres of $T>10^{7}\,\mathrm{K}$
(Biffi et al., 2012). Thus, an intrinsic dusty absorption component was not
required to realistically model cluster emission. In order to emulate a dusty
absorber intrinsic to each selected galaxy of our data set, we applied a
constant absorption model at the redshift of the selected halo consistent with
ISM absorption column densities. While it improved the global relation
compared to when no absorption is applied, it is not clear whether a strong
absorption component is suited for halo outskirts beyond $0.5\,R_{2500}$.
Furthermore, numerous studies explicitly mention that intrinsic absorption is
highly dependent on the observed galaxy and may vary significantly (e.g. Mineo
et al., 2012b; Lehmer et al., 2016, 2019, 2022).

  * •

In several galaxies with elevated $L_{X}$ we found isolated pockets of dense
hot gas clumps, despite specifically filtering out overheated star-forming gas
elements from the simulation. This issue has to be investigated further but
may be connected to the chosen projected volumes which we will explain in one
of the following points.

  * •

Recent studies showed that the extended hot gas atmosphere can significantly
impact the total observed luminosity depending on the chosen projection radius
(Nica et al., 2021; Comparat et al., 2022). In Comparat et al. (2022) they
studied the X-ray emission of the circum-galactic medium (CGM) for halos of
stellar mass $2\cdot 10^{10}-10^{12}\,M_{\odot}$ for different projection
radii (80 kpc, 300 kpc). They further differentiated between star-forming and
quiescent galaxies. They found that the enclosed luminosity in the 0.5-2 keV
band increases by a factor of at least 2 for star-forming galaxies and $>10$
for quiescent galaxies if the enclosing projected radius is increased from 80
kpc to 300 kpc. In addition, they find that while quiescent galaxies clearly
have an extended hot halo, star-forming galaxies share radial profiles with an
expected diffuse unresolved component and weak AGN. This may indicate that in
order to more properly study gas scaling relations for the ISM, the projected
radius has to be chosen more carefully.

  * •

Regarding the influence of AGN towards ISM properties, the Magneticum
simulations include a thermal feedback model for SMBH described in Hirschmann
et al. (2014). While SMBH properties are tracked by the simulations we did not
specifically check for AGN activity in each of our selected halos. It is
however known, that AGN experience different behaviour compared to normal
galaxies with respect to their X-ray radial profiles and scaling relations
(Comparat et al. (2019); Vulic et al. (2021)). A proper disentanglement of AGN
contribution to these relations goes beyond the scope of this paper. We prefer
to more thoroughly analyse their effects in a future study.

  * •

Another influencing factor for gas properties in the ISM are ongoing mergers
between galaxies. In some instances we found that our selected halos had
equally massive companions within $R_{2500}$, possibly introducing an
additional heating source through merger triggered interaction. We did not
exclude merging halos from our sample since we selected for a good
representation of different classes of objects for our XRB model.

  * •

Probably, the most driving factor of scatter especially towards low $L_{X}$ is
the resolution limit of the total gas mass in the cosmological simulation. SPH
simulations usually derive hydrodynamical quantities by interpolating between
neighbouring SPH elements with a set number of neighbours. As soon as a halo
gas mass reaches below a lower mass limit, hydrodynamical effects like star-
formation can no-longer be consistently resolved since the number of SPH
elements in the halo becomes smaller than the needed neighbour number. The
scatter is then dominated by the stochastic process with which sub-grid star-
formation is handled in the implementation (see Springel & Hernquist, 2003).

We want to emphasize again that the hot gas component was not a major focus of
this paper. Nonetheless, our approach shows that the simulation can recover a
linear trend between gas X-ray luminosity and SFR if we select sufficiently
resolved galaxies. However, the large scatter of over $\sim 1\,\mathrm{dex}$
shows that a simple approach to model the ISM is not optimal and needs more
careful investigation.

###  8.2 XRB scaling relations

The main purpose of this paper was to introduce the viability of a semi-
analytic approach to model X-ray emission from XRBs using hydrodynamical
cosmological simulations. We connected locally derived luminosity functions of
XRBs and their intrinsic behaviour to hydrodynamical quantities from the
simulations and showed the self-consistent emergence of expected scaling
relations and spatial distribution. Here we discuss some interesting aspects
regarding our results.

The strong scatter in the $L_{X}-\mathrm{SFR}$ relation of HMXBs is a major
concern. While the expected Poissonian noise from discrete sampling of the LFs
was predicted by Gilfanov et al. (2004b) additional uncertainties arise from
our SFH proxy (see section 3). For the birthrate of massive stars (eq. (9)) we
assumed a linear dependency on SFR which is only true as a first order
approximation and is highly sensitive to the underlying IMF. One would need to
model the complete SFH of the star-forming region to account for time-
dependence. Since we are not able to recover this information from SSPs
cosmological simulations, this first order approximation is the best we can
achieve. For the lifetime function of stars we assumed the prescription of
Padovani & Matteucci (1993) (PM93), which is also implemented in the
simulations (Tornatore et al., 2007). The advantage of the PM93 model is its
simple functional form. There are more recently developed stellar lifetimes
using evolutionary tracks calculated by Parsec (Bressan et al., 2012). Such
lifetime functions are nonetheless derived for single stars while it is known
that evolution of massive binary systems heavily depends on orbital, stellar
wind and common envelope interactions (see e.g. Lewin & van der Klis, 2006, as
a review on XRB evolution) which might drive inconsistencies in stellar age
determination. For this reason, we decided to keep the PM93 model. Another
minor inconsistency in our model assumptions stems from the expected duration
of the X-ray bright phase of HMXBs. We assume a constant timescale of 0.1 Myr
which is consistent with the most common type of Be-HMXB but does not account
for HMXB types with shorter bright-phases. One of the most important drivers
of uncertainty in the $L_{X}-\mathrm{SFR}$ relation for HMXBs at low SFR
values, is the limited mass resolution of SSPs in the simulation. If the mass
resolution decreases the probability to have SSPs in the required age range of
$\tau_{\text{SSP}}\in[0,30]\,\text{Myr}$ becomes smaller. Because there will
be fewer SSPs per galaxy to resolve SFR it adds another layer of Poisson
noise. With smaller SSP masses we would be able to more accurately estimate
the HMXB population size for each galaxy as is already seen in the LMXB case
where the scatter in the $L_{X}-M_{*}$ relation is much smaller. In fact,
Kouroumpatzakis et al. (2020) showed that the $L_{X}-\mathrm{SFR}$ relation on
sub-galactic scales should be consistent with galaxy-wide scaling relations.
They argue that their shallower slopes in the sub-galactic relation is caused
by LMXB contamination and differences in local SFH. Nonetheless, the overall
trend in our recovered $L_{X}-\mathrm{SFR}$ relation is still consistent with
observations.

In our investigation of the metallicity dependence we concluded that there is
only a weak link between SFR and mass-weighted stellar metallicity in the
simulation since the HMXB model from Mineo et al. (2012a) shows no significant
increase in total luminosity with decreasing oxygen fraction. In fact, it
instead suggests that there is a positive correlation between oxygen content
and total luminosity as indicated by the dashed purple line in figure 10 in
contrast to observations by Lehmer et al. (2021). We attribute this behaviour
of the M12 model to our small low Z sample and to the fact that our approach
does not account explicitly for metallicity. Sampling from the metallicity
dependent L21 model performs as expected and correctly reproduces
observational results. Note, that our low SFR sample (yellow) for
$12+\log[\mathrm{O/H}]\leq 8.5$ has much higher luminosities than predicted by
the MCMC chains of L21 for $\mathrm{SFR}\sim 0.1$. Apart from the number
normalization issues discussed above for low resolution SSPs, we speculate
that this difference might either be caused by a different zero-point in
metallicity or differences between mass-weighted stellar metallicities and
gas-phase metallicities of which the latter is used in observations. We have
to compare to the mass-weighted stellar metallicity in the simulations since
we rely on the metallicity estimates of the SSPs. Considering the formation
process of SSPs in the Springel & Hernquist (2003) prescription and the
Tornatore et al. (2007) metal distribution, both gas-phase and stellar
metallicities should be comparable in the simulations (see e.g. Dolag et al.,
2017).

In figure 9 we compared the $L_{X}-\mathrm{SFR}-M_{*}$ relation of our sample
with recent observations of galaxies within the Virgo cluster from Soria et
al. (2022). We find excellent agreement with our data. Due to the break in our
$L_{X}-\mathrm{SFR}$ relation for HMXBs our expected transition of LMXB to
HMXB dominated galaxies, around $\log\mathrm{sSFR}\eqsim-10.5$), allows for
slightly lower $L_{X}$ compared to a linear scaling in $L_{X}-\mathrm{SFR}$.
The normalization of our $L_{X}-\mathrm{SFR}$ relation in the limit of high
sSFR appears to be higher compared to the selected studies. In fact our model
does not rule out the purely linear $L_{X}-\mathrm{SFR}-M_{*}$ relation and is
equally consistent to the model in eq. (19) if we consider a lower
normalization for the $L_{X}-\mathrm{SFR}$ relation.

###  8.3 Average Galaxy spectra

When we investigated the different components contributing to galactic X-ray
spectra, we only used the ideal photon lists and did not perform a full mock
observation using instrumental response with subsequent unfolding of spectra.
While different instruments are more or less suited to detect X-ray point
sources like XRBs, we did not think that a full mock observation would prove
to be beneficial in order to present the XRB algorithm. Since real
observations on X-ray spectra of galaxies report unfolded spectra, we skipped
the instrumental response and reconstruction in favor of cleaner and better
data. Our analysis of the average galactic spectra showed that XRBs can have
significant contribution for energies $E\leq 2\,\mathrm{keV}$. We note however
that our approach suffers from the same spectral simplification as discussed
in Sazonov & Khabibullin (2017a). They argue that variability of XRB spectra
has to be accounted for individually including intrinsic absorption and
spectral hardness for each source. Otherwise, the luminosity estimates for
soft sources can be significantly different. They suggest taking into account
ISM maps of the host galaxies and individually resolving spectral shapes to
improve upon luminosity estimates for XRBs. Given the modular design of our
XRB module for the photon simulator Phox, spectral models for XRBs can easily
be updated to account for XRBs in different hardness states, e.g. using
spectral characteristics from Sazonov & Khabibullin (2017b). From the same
arguments, the estimated count ratios might also be slightly different if some
sources had higher flux in the soft band.

The mean X-ray spectrum of our star-forming sample is consistent with
estimates from Lehmer et al. (2022) although their derived SED is based on
observations of dwarf galaxies with sub-solar metallicity as high redshift
analogues to local galaxies. The good agreement with our data is caused by our
hot gas luminosity estimates being in part 0.5 dex above conventional values
for the ISM (see figure 6). An increase of gas luminosity for lower
metallicity environments in from L22 thus improves the visual overlap but does
not correspond to a lower metallicity in our sample. Another issue is the
strong suppression of metal lines in the tabulated SED in L22. We accounted
for that by convolving our derived spectrum with a fairly wide Gaussian kernel
which still does not suppress some very prominent line features at high
energies. We also find that the shape of our simulated SEDs are flatter at
$E\sim 1-2\,\mathrm{keV}$ compared to the L22 SED. Differences in spectral
shape, especially in the soft energy regime, are highly dependent on the
adopted absorption models as discussed for figure 12 which leads in our case
to overestimation of the XRB component when compared to values from Lehmer et
al. (2016). In the hard band the L16 estimates are consistent with our derived
values. Splitting the average spectrum into different SFR bins revealed
stronger contributions from the hot gas component in high SFR galaxies. This
is due to the influence from BCG-like galaxies in the sample with hot gaseous
atmospheres. In those cases, the L16 values are no longer suited to compare
our derived values against. In Gilbertson et al. (2022) they showed that the
expected influence of the hot gas component is also a function of stellar
population age with the hot gas contribution declining to less than 5% in
older populations which was not investigated in our presented sample. In
addition to a SFR normalized emission model which is consistent with the L22
spectrum, Gilbertson et al. (2022) also report stellar mass normalized average
spectra which we did not compare against as it goes beyond the scope of this
paper. We will conduct a future, more detailed investigation of the influence
of the hot gas component to the average spectrum, where we we will also drill
down on line properties within the ISM and CGM.

##  9 Summary and Conclusions

Based on SNII rate estimates using SSP ages and masses tracked by numerical
cosmological simulations we are able to self-consistently reconstruct a SFR
proxy suitable for Monte-Carlo-like sampling of XRB luminosity functions. We
presented a possible implementation within the Phox code and verified the
validity of our approach using simulated galaxies from the Magneticum
Pathfinder simulations. We extracted galaxy-sized halos from the simulation
based on a stellar mass criterion ($10^{9.5}\,M_{\odot}\leq M_{*}\leq
10^{12}\,M_{\odot}$) and produced ideal photon lists for the hot gas and XRB
component using the stochastic photon simulator Phox. XRB emission was assumed
to follow a static power-law with photon index $\Gamma=1.7-2$ while hot gas
was modelled assuming a single temperature APEC model with thermal line-
broadening and varying metal abundance for each gas element in the simulation.
In order to account for self-absorption of the ISM at the source location we
additionally require an intrinsic absorption component for both hot gas and
XRB component with column density $N^{\mathrm{gas}}_{H}=5\cdot
10^{21}\mathrm{cm}^{-2}$ and $N^{\mathrm{xrb}}_{H}=2\cdot
10^{21}\mathrm{cm}^{-2}$ respectively. For each halo in our sample we
constructed cylindrical volumes around the center of mass with $R=R_{2500}$
and $h=2\,R_{2500}$ and projected the enclosed photons accounting for redshift
and peculiar motion of the source w.r.t. the line of sight. Our findings using
the resulting galactic X-ray spectra can be summarized as follows:

  1. 1.

The global reconstructed XRB luminosity functions perfectly resemble the
underlying analytic shape. The observed flattening of the reconstructed XLFs
below $\log L_{XRB}\approx 37$ is consistent with the expected one-photon
luminosity. The reconstructed XLF derived from the metallicity dependent model
from Lehmer et al. (2021) most closely resembles the XLF shape for solar
metallicity which is almost indistinguishable from the Mineo et al. (2012a)
model. A more curated sample of galaxies would be needed in order to reliably
disentangle the effects of metallicity on the XLF shape which may be analyzed
in a separate study.

  2. 2.

The $L_{X}$-$M_{*}$ relation for LXBs is tightly constrained and is in
excellent agreement with the Zhang et al. (2012) relation. More recent
measurements of this relation by Lehmer et al. (2020) are also consistent with
our data.

  3. 3.

The $L_{X}$-$\mathrm{SFR}$ relation for HXBs suffers greatly from low-number
statistics for $\mathrm{SFR}\lesssim 1\,M_{\odot}\,\mathrm{yr}^{-1}$.
Additionally, the SFR proxy used to determine the number of HXBs per SSP has
large intrinsic uncertainties stemming from the limited mass resolution of the
simulation and assumptions for stellar lifetime and IMF. Despite these
drawbacks, our approach is still consistent with a broken powerlaw relation
which is expected to arise from incomplete sampling of the HXB XLF dominated
by Poissonian noise (Gilfanov et al. (2004b)).

  4. 4.

The $L_{X}$-$\mathrm{SFR}$-$Z$ relation as proposed by Lehmer et al. (2021)
does not arise when sampling from the M12a relation indicating a weak
correlation between SFR and Z within the simulation. When employing the L22
model directly, the correlation becomes more obvious but continues to suffer
from the same uncertainties inherent to the HXB sampling.

  5. 5.

Our recovered $L_{X}/\mathrm{SFR}$-$\mathrm{sSFR}$ relation studying the
combined effect of both XRB types is in remarkable agreement with recent
studies of field galaxies within the Virgo cluster from Soria et al. (2022).
The shape of our derived relation shows more complex behaviour at
$\log\mathrm{sSFR}\sim-10$ compared to the conventional relation found in the
literature (Lehmer et al., 2010, 2016, 2019) due to the observed break in the
HXB relation.

  6. 6.

As expected we recover that HXBs spatially coincide with the star-forming
regions of their host galaxies while the spatial LXB distribution follows the
stellar-mass surface density.

  7. 7.

We find that approximately 40-50% of our galaxy sample is dominated by the
combined XRB emission compared to emission from the hot gas component in the
0.5-8 keV energy band. In the soft band only 30-40% of galaxies are dominated
by gas emission, while more than $90\%$ are dominated by XRB emission in the
hard band. The surprisingly low contribution of HXBs to these percentages
reflects the relatively small portion of high SFR galaxies within our sample.

  8. 8.

We constructed average SFR normalized galactic X-ray spectra form our complete
star-forming ($\mathrm{SFR}>0.4$) sample and measured the relative
contribution of the two XRB types and the hot gas component towards the total
spectrum. The average spectrum of the total star-forming sample is consistent
with spectral fits of Lehmer et al. (2022) in terms of magnitude and shape
assuming an average SFR of 1. Contributions from XRBs to the soft part of the
spectrum are inconsistent with values from Lehmer et al. (2016) but consistent
with Wik et al. (2014) which reflects different model assumptions for
absorption and XBR spectra. Additionally, the normalization of the spectrum
only shows minimal variations across different SFR bins and is mostly
influenced by the presence of massive galaxies in the sample.

In this paper we presented an approach to model galaxy X-ray spectra from
cosmological simulations. In addition to the hot gas component in the
simulation we take properties of stellar elements into account and implemented
a fast and consistent algorithm to compute X-ray binary emission from these
elements while reliably separating the two different XRBs types. The presented
method proves its viability as a supporting tool for observers with its highly
modular design which also allows for alterations of employed model spectra for
individual components. Future studies building upon the suggested
implementation will improve sample selection as well as study the influence of
X-ray AGN emission on scaling relations, background and on the retrieved
spectra.

###### Acknowledgements.

SVZ and VB acknowledge support by the _Deutsche Forschungsgemeinschaft, DFG_
project nr. 415510302. KD acknowledges support through the COMPLEX project
from the European Research Council (ERC) under the European Union’s Horizon
2020 research and innovation program grant agreement ERC-2019-AdG 882679. This
research was supported by the Excellence Cluster ORIGINS which is funded by
the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under
Germany’s Excellence Strategy – EXC-2094 – 390783311. The calculations for the
hydrodynamical simulations were carried out at the Leibniz Supercomputer
Center (LRZ) under the project pr83li. We are especially grateful for the
support by M. Petkova through the Computational Center for Particle and
Astrophysics (C2PAP)

## References

  * Anderson et al. (2013) Anderson, M. E., Bregman, J. N., & Dai, X. 2013, ApJ, 762, 106
  * Arnaud (1996) Arnaud, K. A. 1996, in Astronomical Society of the Pacific Conference Series, Vol. 101, Astronomical Data Analysis Software and Systems V, ed. G. H. Jacoby & J. Barnes, 17
  * Babyk et al. (2018) Babyk, I. V., McNamara, B. R., Nulsen, P. E. J., et al. 2018, ApJ, 857, 32
  * Beck et al. (2016) Beck, A. M., Murante, G., Arth, A., et al. 2016, MNRAS, 455, 2110
  * Biffi et al. (2013) Biffi, V., Dolag, K., & Böhringer, H. 2013, MNRAS, 428, 1395
  * Biffi et al. (2012) Biffi, V., Dolag, K., Böhringer, H., & Lemson, G. 2012, MNRAS, 420, 3545
  * Biffi et al. (2018) Biffi, V., Dolag, K., & Merloni, A. 2018, MNRAS, 481, 2213
  * Biffi et al. (2014) Biffi, V., Sembolini, F., De Petris, M., et al. 2014, MNRAS, 439, 588
  * Biffi & Valdarnini (2015) Biffi, V. & Valdarnini, R. 2015, MNRAS, 446, 2802
  * Bogdán et al. (2013) Bogdán, Á., Forman, W. R., Vogelsberger, M., et al. 2013, ApJ, 772, 97
  * Boroson et al. (2011) Boroson, B., Kim, D.-W., & Fabbiano, G. 2011, ApJ, 729, 12
  * Bressan et al. (2012) Bressan, A., Marigo, P., Girardi, L., et al. 2012, MNRAS, 427, 127
  * Brorby et al. (2016) Brorby, M., Kaaret, P., Prestwich, A., & Mirabel, I. F. 2016, MNRAS, 457, 4081
  * Buchner et al. (2014) Buchner, J., Georgakakis, A., Nandra, K., et al. 2014, A&A, 564, A125
  * Chabrier (2003) Chabrier, G. 2003, PASP, 115, 763
  * Comparat et al. (2019) Comparat, J., Merloni, A., Salvato, M., et al. 2019, MNRAS, 487, 2005
  * Comparat et al. (2022) Comparat, J., Truong, N., Merloni, A., et al. 2022, arXiv e-prints, arXiv:2201.05169
  * Cui et al. (2016) Cui, W., Power, C., Biffi, V., et al. 2016, MNRAS, 456, 2566
  * Di Matteo et al. (2005) Di Matteo, T., Springel, V., & Hernquist, L. 2005, Nature, 433, 604
  * Dolag et al. (2009) Dolag, K., Borgani, S., Murante, G., & Springel, V. 2009, MNRAS, 399, 497
  * Dolag et al. (2017) Dolag, K., Mevius, E., & Remus, R.-S. 2017, Galaxies, 5, 35
  * Dolag et al. (2005) Dolag, K., Vazza, F., Brunetti, G., & Tormen, G. 2005, MNRAS, 364, 753
  * Done et al. (2007) Done, C., Gierliński, M., & Kubota, A. 2007, A&A Rev., 15, 1
  * Fabbiano (1989) Fabbiano, G. 1989, ARA&A, 27, 87
  * Fabbiano (2006) Fabbiano, G. 2006, ARA&A, 44, 323
  * Fabbiano (2019) Fabbiano, G. 2019, in The Chandra X-ray Observatory, ed. B. Wilkes & W. Tucker, 7–1
  * Fabjan et al. (2010) Fabjan, D., Borgani, S., Tornatore, L., et al. 2010, MNRAS, 401, 1670
  * Forbes et al. (2017) Forbes, D. A., Alabi, A., Romanowsky, A. J., et al. 2017, MNRAS, 464, L26
  * Fornasini et al. (2020) Fornasini, F. M., Civano, F., & Suh, H. 2020, MNRAS, 495, 771
  * Fragos et al. (2013a) Fragos, T., Lehmer, B., Tremmel, M., et al. 2013a, ApJ, 764, 41
  * Fragos et al. (2013b) Fragos, T., Lehmer, B. D., Naoz, S., Zezas, A., & Basu-Zych, A. 2013b, ApJ, 776, L31
  * Garofali et al. (2020) Garofali, K., Lehmer, B. D., Basu-Zych, A., et al. 2020, ApJ, 903, 79
  * Gilbertson et al. (2022) Gilbertson, W., Lehmer, B. D., Doore, K., et al. 2022, ApJ, 926, 28
  * Gilfanov (2004) Gilfanov, M. 2004, MNRAS, 349, 146
  * Gilfanov et al. (2004a) Gilfanov, M., Grimm, H. J., & Sunyaev, R. 2004a, MNRAS, 347, L57
  * Gilfanov et al. (2004b) Gilfanov, M., Grimm, H. J., & Sunyaev, R. 2004b, MNRAS, 351, 1365
  * Grimm et al. (2003) Grimm, H. J., Gilfanov, M., & Sunyaev, R. 2003, MNRAS, 339, 793
  * Haardt & Madau (2001) Haardt, F. & Madau, P. 2001, in Clusters of Galaxies and the High Redshift Universe Observed in X-rays, ed. D. M. Neumann & J. T. V. Tran, 64
  * Hirschmann et al. (2014) Hirschmann, M., Dolag, K., Saro, A., et al. 2014, MNRAS, 442, 2304
  * Inoue et al. (2021) Inoue, Y., Yabe, K., & Ueda, Y. 2021, PASJ[arXiv:2107.07775]
  * Irwin (2005) Irwin, J. A. 2005, ApJ, 631, 511
  * Kim & Fabbiano (2013) Kim, D.-W. & Fabbiano, G. 2013, ApJ, 776, 116
  * Kim & Fabbiano (2015) Kim, D.-W. & Fabbiano, G. 2015, ApJ, 812, 127
  * Komatsu et al. (2011) Komatsu, E., Smith, K. M., Dunkley, J., et al. 2011, ApJS, 192, 18
  * Kouroumpatzakis et al. (2020) Kouroumpatzakis, K., Zezas, A., Sell, P., et al. 2020, MNRAS, 494, 5967
  * Lehmer et al. (2010) Lehmer, B. D., Alexander, D. M., Bauer, F. E., et al. 2010, ApJ, 724, 559
  * Lehmer et al. (2016) Lehmer, B. D., Basu-Zych, A. R., Mineo, S., et al. 2016, ApJ, 825, 7
  * Lehmer et al. (2021) Lehmer, B. D., Eufrasio, R. T., Basu-Zych, A., et al. 2021, ApJ, 907, 17
  * Lehmer et al. (2022) Lehmer, B. D., Eufrasio, R. T., Basu-Zych, A., et al. 2022, arXiv e-prints, arXiv:2203.16566
  * Lehmer et al. (2017) Lehmer, B. D., Eufrasio, R. T., Markwardt, L., et al. 2017, ApJ, 851, 11
  * Lehmer et al. (2019) Lehmer, B. D., Eufrasio, R. T., Tzanavaris, P., et al. 2019, ApJS, 243, 3
  * Lehmer et al. (2020) Lehmer, B. D., Ferrell, A. P., Doore, K., et al. 2020, ApJS, 248, 31
  * Lehmer et al. (2015) Lehmer, B. D., Tyler, J. B., Hornschemeier, A. E., et al. 2015, ApJ, 806, 126
  * Lewin & van der Klis (2006) Lewin, W. H. G. & van der Klis, M. 2006, Compact Stellar X-ray Sources, Vol. 39
  * Li & Wang (2013) Li, J.-T. & Wang, Q. D. 2013, MNRAS, 435, 3071
  * Madau & Fragos (2017) Madau, P. & Fragos, T. 2017, ApJ, 840, 39
  * Mineo et al. (2014) Mineo, S., Gilfanov, M., Lehmer, B. D., Morrison, G. E., & Sunyaev, R. 2014, MNRAS, 437, 1698
  * Mineo et al. (2012a) Mineo, S., Gilfanov, M., & Sunyaev, R. 2012a, MNRAS, 419, 2095
  * Mineo et al. (2012b) Mineo, S., Gilfanov, M., & Sunyaev, R. 2012b, MNRAS, 426, 1870
  * Morrison & McCammon (1983) Morrison, R. & McCammon, D. 1983, ApJ, 270, 119
  * Nica et al. (2021) Nica, A., Oppenheimer, B. D., Crain, R. A., et al. 2021, arXiv e-prints, arXiv:2112.14919
  * Padovani & Matteucci (1993) Padovani, P. & Matteucci, F. 1993, ApJ, 416, 26
  * Persic & Rephaeli (2002) Persic, M. & Rephaeli, Y. 2002, A&A, 382, 843
  * Ranalli et al. (2003) Ranalli, P., Comastri, A., & Setti, G. 2003, A&A, 399, 39
  * Remillard & McClintock (2006) Remillard, R. A. & McClintock, J. E. 2006, ARA&A, 44, 49
  * Remus et al. (2017) Remus, R.-S., Dolag, K., Naab, T., et al. 2017, MNRAS, 464, 3742
  * Romanowsky & Fall (2012) Romanowsky, A. J. & Fall, S. M. 2012, ApJS, 203, 17
  * Saxena et al. (2021) Saxena, A., Ellis, R. S., Förster, P. U., et al. 2021, MNRAS, 505, 4798
  * Sazonov & Khabibullin (2017a) Sazonov, S. & Khabibullin, I. 2017a, MNRAS, 466, 1019
  * Sazonov & Khabibullin (2017b) Sazonov, S. & Khabibullin, I. 2017b, MNRAS, 468, 2249
  * Schulze et al. (2020) Schulze, F., Remus, R.-S., Dolag, K., et al. 2020, MNRAS, 493, 3778
  * Schulze et al. (2018) Schulze, F., Remus, R.-S., Dolag, K., et al. 2018, MNRAS, 480, 4636
  * Shtykovskiy & Gilfanov (2005a) Shtykovskiy, P. & Gilfanov, M. 2005a, A&A, 431, 597
  * Shtykovskiy & Gilfanov (2005b) Shtykovskiy, P. & Gilfanov, M. 2005b, MNRAS, 362, 879
  * Smith et al. (2001) Smith, R. K., Brickhouse, N. S., Liedahl, D. A., & Raymond, J. C. 2001, ApJ, 556, L91
  * Soria et al. (2022) Soria, R., Kolehmainen, M., Graham, A. W., et al. 2022, arXiv e-prints, arXiv:2201.07242
  * Springel (2005) Springel, V. 2005, MNRAS, 364, 1105
  * Springel et al. (2005) Springel, V., Di Matteo, T., & Hernquist, L. 2005, MNRAS, 361, 776
  * Springel & Hernquist (2003) Springel, V. & Hernquist, L. 2003, MNRAS, 339, 289
  * Springel et al. (2001) Springel, V., White, S. D. M., Tormen, G., & Kauffmann, G. 2001, MNRAS, 328, 726
  * Steinborn et al. (2016) Steinborn, L. K., Dolag, K., Comerford, J. M., et al. 2016, MNRAS, 458, 1013
  * Strickland et al. (2004) Strickland, D. K., Heckman, T. M., Colbert, E. J. M., Hoopes, C. G., & Weaver, K. A. 2004, ApJS, 151, 193
  * Strickland & Stevens (2000) Strickland, D. K. & Stevens, I. R. 2000, MNRAS, 314, 511
  * Teklu et al. (2015) Teklu, A. F., Remus, R.-S., Dolag, K., et al. 2015, ApJ, 812, 29
  * Teklu et al. (2017) Teklu, A. F., Remus, R.-S., Dolag, K., & Burkert, A. 2017, MNRAS, 472, 4769
  * Tornatore et al. (2007) Tornatore, L., Borgani, S., Dolag, K., & Matteucci, F. 2007, MNRAS, 382, 1050
  * Tornatore et al. (2004) Tornatore, L., Borgani, S., Matteucci, F., Recchi, S., & Tozzi, P. 2004, MNRAS, 349, L19
  * Tyler et al. (2004) Tyler, K., Quillen, A. C., LaPage, A., & Rieke, G. H. 2004, ApJ, 610, 213
  * Vulic et al. (2021) Vulic, N., Hornschemeier, A. E., Haberl, F., et al. 2021, arXiv e-prints, arXiv:2106.14526
  * Weinmann et al. (2010) Weinmann, S. M., Kauffmann, G., von der Linden, A., & De Lucia, G. 2010, MNRAS, 406, 2249
  * Wiersma et al. (2009) Wiersma, R. P. C., Schaye, J., & Smith, B. D. 2009, MNRAS, 393, 99
  * Wik et al. (2014) Wik, D. R., Lehmer, B. D., Hornschemeier, A. E., et al. 2014, ApJ, 797, 79
  * Wiktorowicz et al. (2017) Wiktorowicz, G., Sobolewska, M., Lasota, J.-P., & Belczynski, K. 2017, ApJ, 846, 17
  * Wilms et al. (2000) Wilms, J., Allen, A., & McCray, R. 2000, ApJ, 542, 914
  * Wong et al. (2014) Wong, K.-W., Irwin, J. A., Shcherbakov, R. V., et al. 2014, ApJ, 780, 9
  * Yukita et al. (2016) Yukita, M., Hornschemeier, A. E., Lehmer, B. D., et al. 2016, ApJ, 824, 107
  * Zhang et al. (2012) Zhang, Z., Gilfanov, M., & Bogdán, Á. 2012, A&A, 546, A36
  * Zhang et al. (2011) Zhang, Z., Gilfanov, M., Voss, R., et al. 2011, A&A, 533, A33