© The Authors 2023

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 the physical processes that produce the 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 interstellar medium (ISM) and accretion-powered point sources called X-ray binaries (XRBs; see the reviews Fabbiano 1989, 2006, 2019), excluding the emission from actively accreting super-massive black holes (SMBHs).

Numerous studies of the hot diffuse ISM in star-forming normal galaxies have revealed a tight linear relation between the total X-ray luminosity, L_X, and the star-formation rate (SFR; Strickland & Stevens 2000; Ranalli et al. 2003; Tyler et al. 2004; Mineo et al. 2012a, hereafter M12a, Li & Wang 2013; Kouroumpatzakis et al. 2020) due to feedback from young stellar objects and energy injection from supernovae (SNe), 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 an overall old stellar population, appears to follow gas X-ray luminosity scaling relations consistent with virialized systems, similar to brightest cluster galaxies (BCGs) and groups. Other types of gas-poor EGs behave similarly to disk galaxies in these relations: they do not show significant correlations between hot gas temperature and the X-ray luminosity of hot gas, suggesting secondary effects such as 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 (LMXBs) are accreting compact objects (COs) 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_⊙) donor star in the binary system (Gilfanov 2004; Boroson et al. 2011; Zhang et al. 2012, hereafter Z12; Lehmer et al. 2019, hereafter L19). This relation seems to be enhanced by the globular cluster (GC) specific frequency, $S_N=N_{\text{GC}}{10}^{0.4(M_V^T+15)}$, where gravitationally dissolving GCs seed LMXBs into the galactic field, which join the LMXBs formed in situ (Irwin 2005; Zhang et al. 2011, 2012; Boroson et al. 2011; Lehmer et al. 2020). In high-mass X-ray binaries (HMXBs) the donor star is a massive O/B star ( > 8 M_⊙) that fuels accretion onto the CO with intense stellar winds. Due to the short lifetimes of such massive stars ( < 100 Myr), HMXBs provide an independent tracer of the galactic star-formation history (SFH). It was found that the combined X-ray luminosity of HMXBs relates linearly to the SFR (Grimm et al. 2003; Shtykovskiy & Gilfanov 2005b; Mineo et al. 2012b, hereafter M12b, Lehmer et al. 2010, 2016, 2019). Broadband X-ray luminosity functions (XLFs) 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 the stellar evolution timescales of the XRB donor, transitioning from HMXB dominated to LMXB dominated at a stellar population age of ≳100 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 the specific star-formation rate (sSFR = SFR/M_*) of the form L_X/SFR = α sSFR⁻¹ + β, where α describes the contribution from LMXBs and β the contribution from HMXBs. This relation is especially useful in galaxies where XRB candidates cannot 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, hereafter L16). The turnover from LMXB to HMXB domination occurs at logsSFR ∼ −10.5. Recently, Soria et al. (2022, hereafter S22) performed Chandra observations of normal galaxies in the Virgo cluster that were marginally compatible with scaling relations for local galaxies, indicating non-negligible environmental effects on this relation. A similar dependence was found by Inoue et al. (2021) in the form of a fundamental plane: L_X ∼ SFR + α M_*.

Numerous studies on generalizing locally derived XRB scaling relations for the high redshift Universe have used 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 SFH of the universe constrained with the help of cosmological simulations (MILLENIUM II in the case of Fragos et al. 2013b) suggests a more complex evolution of local XRB scaling relations, with redshift and HMXB emissivity dominating above z ∼ 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). The VANDELS project is an ESO funded spectroscopic survey within the Chandra Deep Field South targeting high-redshift galaxies. Additional evidence for the metallicity dependence of HMXB emissivity comes from a recent study by Lehmer et al. (2021, hereafter L21), who used the gas-phase oxygen fraction of their galactic sample 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 a higher total luminosity. Similarly, earlier studies from Brorby et al. (2016) and Fornasini et al. (2020) find anticorrelated metallicity enhancements in the L_X-SFR relation, which is consistent with findings by L21. Additionally, data from the eROSITA Final Equatorial Depth Survey (eFEDS) show an elevated total X-ray luminosity for low-metallicity dwarf galaxies with high sSFRs (Vulic et al. 2022). This paper introduces a numerical modeling 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 associated XRB populations with 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 Sect. 2 we describe the state-of-the-art hydrodynamical cosmological simulation suite Magneticum, on which we base our analysis. In Sect. 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 (SSPs). In Sect. 4 we describe the process for generating photons from X-ray-emitting sources in the simulations using the virtual photon simulator PHOX. In Sect. 6 we show our results regarding the reconstruction of known XRB and hot gas scaling relations as well as a proper spatial distribution of the two XRB types in the galactic field. In Sect. 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 Sect. 8 we put our findings into context with our methodology and constructed data set by highlighting caveats and uncertainties. Finally, we summarize our findings in Sect. 9.

2. Cosmological hydrodynamical simulation

Magneticum Pathfinder simulations¹ are a series of state-of-the-art hydrodynamical cosmological simulations that explore 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 UV background (Haardt et al. 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⁻¹(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 kinematic 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 make use of a high resolution run of Magneticum called Box4/uhr, which is a (48 h⁻¹ cMpc)³ comoving volume with a mass resolution of m_DM = 3.6 × 10⁷ M_⊙ and m_gas = 7.3 × 10⁶ M_⊙ for dark matter and gas, respectively, and a total of 576³ particles. Initial conditions for the simulations are generated using standard Λcold dark matter cosmology using results from the Wilkinson Microwave Anisotropy Probe (WMAP7; Komatsu et al. 2011), with Hubble parameter h = 0.704, matter density Ω_M = 0.272, dark energy density Ω_Λ = 0.728, baryon density Ω_b = 0.0451 and normalization of the fluctuation amplitude at 8 Mpc σ₈ = 0.809.

3. Modeling of XRB populations within SSPs

In this section we describe the process of modeling 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 SSP that shares the same time of birth and initial metallicity, distributed according to a common initial mass function (IMF). The main difficulty lies in connecting the 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 priori 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 modeling of LMXB is based on the fact that they cannot 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; Z12; 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 re-normalizing the LMXB XLF to the mass of the SSP. In particular, we use the BPL XLF from Z12 throughout this paper, which has the form

$$\begin{array}{c}\hfill \frac{\mathrm{d}{N}_{\mathrm{LMXB}}}{\mathrm{d}L}=A_{\mathrm{LMXB}}\frac{M_{\ast }}{{10}^{11}{M}_{\odot }}\{\begin{array}{cc}{L}^{-{\alpha }_1},\hfill & L<L_{b,1}\hfill \\ L_{b,1}^{({\alpha }_2-{\alpha }_1)}{L}^{-{\alpha }_2},\hfill & L_{b,1}<L<L_{b,2}\hfill \\ L_{b,2}^{({\alpha }_3-{\alpha }_2)}{L}^{-{\alpha }_3},\hfill & L_{b,2}<L<L_{\mathrm{cut}}\hfill \\ 0,\hfill & L\ge L_{\mathrm{cut}}\hfill \end{array},\end{array}$$(1)

where L_b, 1 = 0.546 and L_b, 2 = 5.99 are break luminosities, L_cut = 500 is the cutoff luminosity, A_LMXB = 54.48 is the normalization, M_* is the total stellar mass of a galaxy, and α₁ = 1.02, α₂ = 2.06, α₃ = 3.63 are the respective PL slopes. Luminosities are given in units of 10³⁸ erg s⁻¹. In fact, LMXBs appear to not only be dependent on the stellar mass of the host galaxy, but also be 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; Z12; 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 Z12 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

$$\begin{array}{c}\hfill {\stackrel{\sim }{N}}_{\mathrm{LMXB}}={\displaystyle {\int }_{L_{\min }}^{L_{\mathrm{cut}}}}\frac{\mathrm{d}{N}_{\mathrm{LMXB}}}{\mathrm{d}L}(M_{\mathrm{SSP}})\mathrm{d}L,\end{array}$$(2)

where we use the current SSP mass M_SSP in solar masses instead of M_*. The total stellar mass of a galaxy would then be recovered by the sum over each single SSPs within its boundary.

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; M12b). In particular, the XLF normalization of HMXBs depends linearly on the SFR. A standard HMXB XLF that we use throughout the paper is given in M12b as

$$\begin{array}{c}\hfill \frac{\mathrm{d}{N}_{\mathrm{HMXB}}}{\mathrm{d}L}=A_{\mathrm{HMXB}}\frac{\mathrm{SFR}}{M_{\odot }{\mathrm{yr}}^{-1}}\{\begin{array}{cc}{L}^{-\gamma },\hfill & L<L_{\mathrm{cut}}\hfill \\ 0,\hfill & L\ge L_{\mathrm{cut}}\hfill \end{array},\end{array}$$(3)

where A_HMXB is the normalization constant, SFR of a galaxy is given in units of M_⊙ yr⁻¹, γ = 1.59 is the PL slope, and L_cut = 10³ is the cutoff luminosity. The luminosity L is given in units of 10³⁸ erg s⁻¹. One can obtain the number of individual HMXBs in a galaxy with

$$\begin{array}{c}\hfill N_{\mathrm{HMXB}}={\displaystyle {\int }_{L_{\min }}^{L_{\mathrm{cut}}}}\frac{\mathrm{d}{N}_{\mathrm{HMXB}}}{\mathrm{d}L}\mathrm{d}L.\end{array}$$(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 δ-function at their time of creation (Tornatore et al. 2007). To circumvent this issue, we can adapt the procedure in M12b by obtaining the fraction f_X of COs becoming a XRB, which is derived from the birth rate of massive stars (> 8 M_⊙) in a star-formation event. Instead of using the birthrate Ṅ_CO of COs following a star-formation event, which is directly dependent on the IMF and SFR, we can equally use the type II supernova (SNII) rate R_SNII, which hides the SFR dependence in the stellar lifetime function. Following Tornatore et al. (2007), we can express the SNII rate as

$$\begin{array}{c}\hfill R_{\mathrm{SNII}}(\tau )=\mathrm{\Phi }(m(\tau ))\times \frac{\mathrm{d}m(\tau )}{\mathrm{d}\tau },\end{array}$$(5)

where Φ(m) is a mass-normalized IMF, which gives the number of stars given the mass m, and m(τ) is the inverse stellar lifetime function of Padovani & Matteucci (1993; hereafter PM93), which gives the mass of stars dying at an age of τ. The SNII rate R_SNII is given in units of Gyr⁻¹ for a SSP of mass 1 M_⊙. The PM93 lifetime function has the form

$$\begin{array}{c}\hfill \tau (m)=\{\begin{array}{cc}{10}^{(1.338-\sqrt{1.79-0.2232·(7.764-logm}))/0.1116-9},\hfill & m\le 6.6M_{\odot }\hfill \\ 1.2·m^{-1.85}+0.003,\hfill & m>6.6M_{\odot }\hfill \end{array},\end{array}$$(6)

with mass m of the star given in M_⊙ and τ in units of Gyr. We employed the modified Chabrier IMF (Chabrier 2003) used in the simulation, which has the form

$$\begin{array}{c}\hfill \mathrm{\Phi }(m)\propto \{\begin{array}{cc}{m}^{-2.3},\hfill & 1.0M_{\odot }<m<100M_{\odot }\hfill \\ m^{-1.8},\hfill & 0.5M_{\odot }<m<1.0M_{\odot }\hfill \\ m^{-1.2},\hfill & 0.1M_{\odot }<m<0.5M_{\odot }\hfill \end{array}.\end{array}$$(7)

Using the SNII rate, we followed the calculations in M12b to obtain an estimate of the HMXB population size.

First we determined the product of the X-ray-bright fraction, f_X, and the average bright-phase duration, ${\overline{\tau }}_X$, assuming an XLF

$$\begin{array}{c}\hfill f_X{\overline{\tau }}_X\sim \frac{N_{\mathrm{HMXB}}(>L_{\min })}{\dot{N}(>8M_{\odot })},\end{array}$$(8)

where Ṅ(> 8 M_⊙) is the birthrate of massive stars. The birthrate can be obtained from the IMF to the first order using

$$\begin{array}{c}\hfill \dot{N}(>8M_{\odot })=\frac{{\int }_8^{M_{\mathrm{u}}}\mathrm{\Phi }(m)\mathrm{d}m}{{\int }_{M_{\mathrm{l}}}^{M_{\mathrm{u}}}\mathrm{\Phi }(m)m\mathrm{d}m}\times \mathrm{SFR},\end{array}$$(9)

with integration limits defined by the IMF. M12b determined instead the term in Eq. (8) as $f_X\sim 0.18\frac{0.1\mathrm{Myr}}{\overline{\tau }}$, based on binary evolution calculations for the most common type of Be HMXB and assuming ${\overline{\tau }}_X\sim 0.1\mathrm{Myr}$, therefore making the factor explicitly dependent on the bright phase duration. In our case, by computing directly the product of f_X and $\overline{\tau }$ as in Eq. (8), we can conveniently bypass uncertainties connected to assumptions for the bright phase duration. Given the proportionality to SFR in Eq. (9) and of N_HMXB (see Eq. (4)), then the dependence of $f_X{\overline{\tau }}_X$ on SFR disappears.

Then, from Eq. (8), we can replace the birthrate, Ṅ, with the SNII rate, R_SNII, given the age of the SSP, and multiply by the current SSP mass to get the expected HMXB population size:

$$\begin{array}{c}\hfill {\stackrel{\sim }{N}}_{\mathrm{HMXB}}=\{\begin{array}{cc}{M}_{\mathrm{SSP}}\times R_{\mathrm{SNII}}(t_{\mathrm{SSP}})f_X{\overline{\tau }}_X,\hfill & t_{\mathrm{SSP}}\le \tau (8M_{\odot })\hfill \\ 0,\hfill & t_{\mathrm{SSP}}>\tau (8M_{\odot })\hfill \end{array}.\end{array}$$(10)

Additionally, we require t_SSP to be smaller than the lifetime of stars with mass 8 M_⊙, in line with mass limits for HMXBs in the literature (see, e.g., Lewin & van der Klis 2006), which, according to Eq. (6), gives ≈30 Myr. Therefore, in the simulation, a stellar particle older than 30 Myr should no longer represent a SSP containing stars with masses greater than 8 M_⊙. The same procedure can be repeated for different IMFs and XLF models.

In order to study the metallicity dependence of HMXBs, we also included the model proposed by L21:

$$\begin{array}{c}\hfill \frac{\mathrm{d}{N}_{\mathrm{HZB}}}{\mathrm{d}L}=A_{\mathrm{HZB}}\mathrm{SFR}exp[-L/L_c(Z)]\{\begin{array}{cc}{L}^{-{\gamma }_1},\hfill & L<L_b\hfill \\ L_b^{{\gamma }_2(z)-{\gamma }_1}{L}^{{\gamma }_2(Z)},\hfill & L\ge L_b\hfill \end{array},\end{array}$$(11)

with

$$\begin{array}{c}\hfill {\gamma }_2(Z)={\gamma }_{2,\odot }+\frac{\mathrm{d}{\gamma }_2}{\mathrm{d}logZ}[12+log(\mathrm{O}/\mathrm{H})-8.69]\end{array}$$(12)

and

$$\begin{array}{c}\hfill log{L}_c(Z)=log{L}_{c,\odot }+\frac{\mathrm{d}log{L}_c}{\mathrm{d}logZ}[12+log(\mathrm{O}/\mathrm{H})-8.69].\end{array}$$(13)

Both γ_2, ⊙ = 1.16 and log L_c, ⊙ = 1.98 are reference values at solar metallicity. The other parameter values are the normalization A_HZB = 1.29, break luminosity log L_b = 0.54, slope γ₁ = 1.74 and first-order metallicity corrections $\frac{\mathrm{d}{\gamma }_2}{\mathrm{d}logZ}=1.34$ and $\frac{\mathrm{d}log{L}_c}{\mathrm{d}logZ}=0.6$. Luminosities are given in units of 10³⁸ erg s⁻¹ in Eq. (11). We labeled the L21 model “HZB” 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, AGN, and XRBs from the simulations, where we focus on the XRB emission model.

4.1. PHOX X-ray photon simulator

In this work we generated synthetic X-ray observations of the hot gas, actively accreting SMBHs, and the stellar component in the simulations using the PHOX code (see 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 are stored such that they reflect the position and velocity information of the emitting sources in the simulation.

Unit 2. The photon data are 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 intra-cluster medium (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 PL) 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 modeled based on their intrinsic thermal and chemical properties (density, temperature, metallicity). A single-temperature thermal emission model with heavy element emission lines based on the Astrophysical Plasma Emission Code (APEC; 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 black holes in the simulations is modeled by determining their accretion rate with respect to their Eddington accretion rate. In order to represent the effect of dense torus absorption, the X-ray spectrum is modeled using an intrinsically absorbed PL 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 modeling of AGN and hot gas emission in the PHOX code, we similarly determined the X-ray emission coming from a stellar resolution element of the simulation. We made use of the SSP properties determined in Sect. 3 to calculate an expected number of point sources $\stackrel{\sim }{N}$ per stellar element. We then sampled the respective XLF based on the following recipe. We selected the correct XLF from a stellar age criterion: if the age of the SSP is τ < 30 Myr, we classified it as HMXB hosting. If τ > 1 Gyr we classified the SSP as LMXB hosting. Both choices are based on expectations from binary evolution simulations (see, e.g., the review by Lewin & van der Klis 2006) and are well motivated from observations. In the HMXB case Shtykovskiy & Gilfanov (2007) found peak formation efficiency in the Small Magellanic Cloud at ≈50 Myr, Garofali et al. (2018) found peak formation efficiency in M33 at ≈40 Myr, and Antoniou et al. (2019) find peak formation efficiency again in the Small Magellanic Cloud at ≈30 − 40 Myr, all by correlating the specific SFH of observed stellar regions with HMXBs. Since these values are dependent on the stellar models assumed for the observed regions we settled for the value given by the lifetime function (see Sect. 3). In the LMXB case the 1 Gyr boundary yields a donor mass limit of 1.75 M_⊙, which reflects the Roche-lobe overflow scenario of accretion onto a CO.

From the XLF we constructed a pseudo cumulative density function (pCDF). This requires setting a luminosity interval [L_min, L_max] at which the pCDF C(L) will be defined:

$$\begin{array}{c}\hfill C(L)=1-\frac{N(>L)}{N(>L_{\min })},\end{array}$$(14)

with N(> L) being the integrated XLF from step 1, and C(L) ∈ [0, 1]. Typically, L_max is chosen to be the cutoff luminosity of the respective XLF. It should be noted that C(L) is sensitive to the choice of L_min due to the PL dependence of the XLF. We then calculated $\stackrel{\sim }{N}$ from Eqs. (2) and (10) and chose $\lfloor \stackrel{\sim }{N}\rfloor$ uniformly distributed random numbers p_i ∈ [0, 1]. If a separate uniform random number $\stackrel{\sim }{p}$ satisfies $\stackrel{\sim }{p}<\stackrel{\sim }{N}-\lfloor \stackrel{\sim }{N}\rfloor$, we draw an additional random number p_i. For each p_i we determined a corresponding luminosity L_i, at which C(L_i) = p_i. From this we can calculate the total luminosity of the SSP as

$$\begin{array}{c}\hfill L_{\mathrm{SSP}}={\displaystyle \sum _i}{L}_i.\end{array}$$(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 described XRB emission spectra as redshift-dependent, absorbed PLs of the form

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

where K is the spectrum normalization at 1 keV and Γ is the photon index as given in the XSPEC manual (Arnaud 1996). Absorption was modeled after the Wilms et al. (2000) photo-ionization cross sections, and the model includes 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 modeling in a cosmological context. The total observed luminosity of the SSP (Eq. (15)) in a certain energy band E ∈ [E₁, E₂] is then related to Eq. (16) by

$$\begin{array}{c}\hfill L_{\mathrm{SSP}}=F{\displaystyle {\int }_{E_1}^{E_2}}A(E)E\mathrm{d}E,\end{array}$$(17)

with F = 4π𝒟_L(z)² × 1.602 ⋅ 10⁻⁹erg keV⁻¹ being the rescaling factor between flux and luminosity and 𝒟_L(z) the luminosity distance inferred from the underlying cosmology. The energy range [E₁, E₂] has to be chosen in correspondence to the range in which the XLFs were defined, which is commonly adopted as E₁ = 0.5 keV and E₂ = 8 keV in the observed frame. We are then able to constrain the spectrum normalization K_SSP using Eqs. (15) and (17), obtaining

$$\begin{array}{c}\hfill K_{\mathrm{SSP}}=\frac{{\sum }_i{L}_i}{F{\int }_{E_1}^{E_2}w(E){[E(1+z)]}^{-\mathrm{\Gamma }+1}\mathrm{d}E}.\end{array}$$(18)

A fully detailed description of XRB spectra is not possible in our modeling. For instance, the cosmological context of our approach makes it impossible to capture XRB variability on small timescales, so 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 PL slope Γ_LMXB = 1.7 in the LMXB case and Γ_HMXB = 2 in the HMXB case. Both slopes are motivated from observational data where the Γ_HMXB follows assumptions of M12b and Γ_LMXB follows the median slope of all high-confidence XRBs in L19. For the absorption component we assume the median column density $N_H^{\text{xrb}}=2\times {10}^{21}{\text{cm}}^{-3}$ from the same sample of XRBs in L19. 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 Γ_LMXB ∼ 1.4 − 1.8 is sufficient to describe the hot atmosphere of EGs (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 active binaries 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). Typically applied corrections for unresolved X-ray emission of cataclysmic variables and active binaries are $\sim 8\cdot {10}^{27}\text{erg}{\text{s}}^{-1}{M}_{\odot }{}^{-1}$ according to Babyk et al. (2018). For young stellar objects, M12a estimate a collective unresolved emission of ∼2 × 10³⁸ erg s⁻¹ per unit SFR.

Having determined K, Γ and $N_H^{\text{xrb}}$ of a single SSP, we can compute the resulting XRB spectrum by adopting the ztbabs and zpowerlw model embedded in XSPEC, from which we calculated an ideal photon list associated with 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_eff = 1000 cm² and exposure time T_exp = 10⁵ s, we selected all halos generated by the SUBFIND algorithm (Springel et al. 2001; Dolag et al. 2009) that have a total stellar mass M_* > 10^9.7 M_⊙ within a sphere of R₂₅₀₀ around the halo center, which also includes a few massive group-like objects with stellar masses of M_* ∼ 10¹² M_⊙. Halos with stellar masses below a threshold of M_* < 10¹⁰ M_⊙ are resolution limited: the SPH implementation fails to correctly reproduce the required star-formation for low-mass halos. We illustrate the M_gas − M_* relation in Fig. 1. Halos that 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.

Fig. 1. Mgas − 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 SFR > 0.4 M_⊙ yr⁻¹, with SFR derived from the total mass of stars born in the past 100 Myrs within the projected volume of each galaxy. This SFR estimator thus probes similar timescales of the SFH as typical tracers used in the literature such as in M12b. X-ray luminosities were calculated from the photon list obtained from a cylindrical volume around each galaxy with radius R₂₅₀₀ and projected along the l.o.s. for a length of 2R₂₅₀₀. The l.o.s. direction coincides with the z-axis of the simulation box. We chose R₂₅₀₀ as the smallest available scale invariant size from the halo catalog and to focus our analysis on the inner regions of selected galaxies in our sample. Reported comparison values probe projected radii on ISM scales consistently smaller than R₂₅₀₀: Strickland et al. (2004, hereafter S04) and M12b within the D₂₅ isophotes, Bogdán et al. (2013, hereafter B13) within 0.05–0.15R₂₀₀ (20–60 kpc), L16 within ≳10 kpc (see chapter on stacking procedure) and Lehmer et al. (2022) within ∼50 kpc.

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_H^{\text{gas}}=5\times {10}^{21}\text{cm}$ at the source redshift (see, e.g., M12a; Gilbertson et al. 2022; Lehmer et al. 2022). In Fig. 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). To increase the sample size of low-metallicity star-forming galaxies, we followed the approach taken by Weinmann et al. (2010), assigning a uniformly distributed log sSFR in the range –12.4 to –11.6 for halos with SFR = 0. This can be seen in Fig. 2 as an apparent strip of mostly spheroidal galaxies located around the log 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.

Fig. 2. SFR − M* relation of our sample. The dashed diagonal lines correspond to constant logsSFR 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 the classification scheme from Schulze et al. (2020; red: elliptical, green: intermediate, blue: disk). Galaxies clustered at logsSFR = −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. Validation of the modeling

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 selected a poster-child late-type galaxy from our sample with M_* = 10^11.3 M_⊙, SFR = 8.4 M_⊙ yr⁻¹ and R₂₅₀₀ = 91 kpc that appears face on to the line of sight. We focus our showcase to the inner 60 kpc in order to better capture the main stellar body of the galaxy. In Fig. 3 we show the percentile contours of photons emitted by HMXBs (black) and LMXBs (red) on top of the line-of-sight projected stellar mass map of the galaxy. As expected, the LMXB 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 (Σ_*). In contrast, HMXB emission is confined to the stellar mass overdensities surrounding the galaxy center. The reason for this can be seen in Fig. 4 where we show the same X-ray contours on top of the mass-weighted stellar age map of the galaxy. It shows that HMXB 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 behavior 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.

Fig. 3. Projected stellar surface mass density (Σ*) of the poster-child disk galaxy in Box4 of Magneticum with halo ID 13633. Contours indicate the 68th, 90th, and 99th percentiles of X-ray photon positions for HMXBs (black) and LMXBs (red).

Fig. 4. Mass-weighted stellar age (τ*) of the poster-child disk galaxy in Box4 of Magneticum with halo ID 13633. Contours indicate the 68th, 90th, and 99th percentiles of X-ray photon positions for HMXBs (black) and LMXBs (red).

6.2. XLFs

For the XLF reconstruction we chose a lower luminosity limit of L_min = 10³⁵ erg s⁻¹, since all XLFs are expected to have a turn-over at luminosities ≲10³⁵ erg s⁻¹. This is due to neutron star XRBs generally experiencing lower Eddington accretion rates than their black hole 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, so lower value choices of L_min will not change the integrated luminosity significantly. In Fig. 5 we compare the integrated XLFs (colored), reconstructed from the approach outlined in Sect. 4.4, with their original model (black). Corresponding XLFs are indicated with the same linestyle. Luminosities of individual XRBs were calculated from the simulated photon lists associated with the SSP element hosting the XRB. In detail, each colored 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 neutron star XRBs but can be explained by the following: (1.) low luminosity XRBs are unlikely to produce observable photons given the distance of the simulation box (low flux); and (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 PL spectrum, is marked by the vertical gray band at L_X ∼ 10³⁷ erg s⁻¹, which roughly coincides with the turnover given the energy range of 0.5 − 8 keV. The solid vertical gray line is the luminosity of a photon with mean energy ∼3 keV according to Eq. (16). Thus, SSPs with L_SSP < 10^37.4 erg s⁻¹ will produce at most one single photon given the discrete photon sampling in PHOX. We 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 line style) above the single photon luminosity limit. The average reconstructed XLF for the L21 model (orange) is consistent with solar values for 12 + log[O/H], which is almost indistinguishable from the M12b XLF. Also shown is the L21 model at 12 + log[O/H] = 9.2 and 7.2, which are not good descriptions of the average HZB XLF.

Fig. 5. Integrated XLFs from Eqs. (1), (3), and (11) normalized by their respective scaling (SFR and M*). The dash-dot-dotted line is the XLFs of the L21 model, for log[O/H]+12 = 7.2, and the dotted line for log[O/H]+12 = 9.2. Colored lines are average XLFs obtained from the luminosity sampling in PHOX: from M12b in dashed blue, L21 in dash-dotted orange, and Z12 in solid green. Corresponding lines share the same line style. The gray vertical line at log L ∼ 37 marks the maximum luminosity of a source emitting one single photon with an average energy of 0.5–8 keV.

6.3. Scaling relations

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

Studies of the diffuse galactic X-ray emission revealed a linear relationship with SFR in the soft X-ray band, which is associated with ISM heating through SN remnants and stellar wind from heavy stars (Ranalli et al. 2003; Gilfanov et al. 2004b; M12a; B13). However, it is observationally challenging to separate pure ISM emission from unresolved emission from point sources. In fact, M12a and Anderson et al. (2013) show that two-thirds 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 (0.5– 2 keV). The blue contours enclose the 68th, 90th, and 99th percentiles of galaxies in our simulated sample also indicated by underlying green dots. We only show galaxies with logsSFR > −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. Also shown is best-fit linear relation from Lehmer et al. (2022) with $log(L_X^{0.5-2\mathrm{keV}}/\mathrm{SFR})=39.{58}_{-0.28}^{+0.17}$, using spectral energy distribution (SED) fitting to derive the SFH for each galaxy in their sample. We include absorbed ISM luminosities from (M12a; red squares) with SFR derived from a combination of UV and IR tracers. Their reported values in the soft X-ray band are smaller compared to the L22 relation, which may be caused by differences in sample properties like metallicity. 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 B13 (black pentagons) in the 0.5–2 keV energy range and a radial range of 0.05–0.15R_vir for late-type galaxies, with SFRs derived from the total IR luminosity. Although B13 aimed to study extended X-ray coronae around large spiral galaxies, the measured luminosities are consistent with a linear L_X − SFR relation prompting a connection between SFR and galactic outflows. From S04 we obtained data points for absorption-corrected total gas luminosities in late-type galaxies in the 0.3–2 keV energy range (black stars), which also show a linear SFR relation. They derived SFRs from the far IR band luminosity. 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 90th 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 SFR ∼ 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 Sect. 8. The intrinsic scatter in our data of > 1 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 composed of 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 cutoff where a halo was not able to form stars (see Fig. 1).

Fig. 6. Relation between hot gas X-ray luminosity (0.5–2 keV) and the SFR for star-forming galaxies in the simulations. We only show galaxies with log(sSFR) >  − 11.6, which removes all galaxies that fulfill the Weinmann et al. (2010) criterion. This is done to reduce the influence of unresolved galaxies on the relation as well as remove 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σ standard deviation. Data points are bolometric luminosities of the diffuse emission taken from S04 in the 0.3–2 keV band, M12a in the 0.5–2 keV band, and 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 Figs. 7 and 8 we show the principal scaling relations for HMXBs and LMXBs, respectively. The contours have the same meaning as in Fig. 6. The solid black line indicates the expected scaling relation from integrating the luminosity functions in Eqs. (3) and (1). In the HMXB case in Fig. 7 the reconstructed scaling relation is consistent with observations from M12b 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. 2004b) and explained in Gilfanov et al. (2004a) where the transition between linear and PL regime of the SFR relation is caused by low-number sampling of a single sloped XLF. The approximate relation calculated from Gilfanov et al. (2004a) for the XLF in Eq. (3) is shown as the dashed black line and has a slope $\frac{1}{\gamma -1}$. For the LMXB case in Fig. 8, the linear stellar mass dependence is well reproduced but our data systematically underestimates the expected relation by a factor of ≈0.15 dex. This is a manifestation of the one-photon limit associated with the XLF sampling. Also shown in Fig. 8 is the integrated luminosity function of Z12 (Eq. (1)) for L_min = 5 × 10³⁷erg s⁻¹, which is approximately the upper bound of the one-photon luminosity limit. Accounting for the under-sampling of the LMXB XLF caused by the one-photon limit thus improves the relation. A similar effect is not seen in the HMXB since the under-sampled XLF only contributes less than 5%. We do not observe a break in the LMXB-M_* relation since our stellar mass limit (M_* = 10^9.7 M_⊙) 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_* ≲ 10⁹) we expect breaks in the L_X − M_* relation as well (Gilfanov et al. 2004a).

Fig. 7. Same as Fig. 6, but for HMXB luminosity in the 0.5–8 keV band. The filled blue contours represent the same halos with HMXBs sampled from the L21 model (Eq. (11)). The solid black line is the expected relation from the employed HMXB model in Eq. (3), and the dashed red line is the predicted relation for discretely sampling data points from a PL distribution (Gilfanov et al. 2004a). The solid blue line is the expected relation from the employed L21 model. Data points are taken from M12b.

Fig. 8. Same as Fig. 7, but for LMXB emission in the 0.5–8 keV band. The solid black line is the expected relation from the employed LMXB model in Eq. (1). The dashed black line is the same relation for a luminosity limit Lmin of 5 × 1037 erg s−1. Data points (cyan triangles) are taken from Lehmer et al. (2020).

When combining the SFR-normalized emission from HMXBs and LMXBs, we expect a sSFR-dependent relation of the form

$$\begin{array}{c}\hfill L_X/\mathrm{SFR}=\alpha {\mathrm{sSFR}}^{-1}+\beta ,\end{array}$$(19)

where α contains the linear mass dependence of the total LMXB luminosity and β contains the linear SFR dependence of the total HMXB luminosity (Lehmer et al. 2010, 2016, 2019). In Fig. 9 we show the sSFR dependence of the SFR-normalized combined XRB luminosity of our sample with contours as in Fig. 6. Additionally, we show data points from Mineo et al. (2014; hereafter M14; SFRs as in M12b) from nearby galaxies and a study performed on star-forming galaxies within the Virgo cluster by S22 with SFR values derived from the 12 μm WISE W3 measurement. We include an observationally derived relation (black, dash-dotted) for local galaxies from L19 where SFRs were derived from a combination of far UV and 24 μm maps. The red dashed diagonal line represents the theoretical contribution from LMXBs given our model choice (Fig. 8) and the blue dashed line corresponds to contributions from HMXBs (Fig. 7). Because of the HMXB scaling relation being more similar to a BPL, the combined luminosity of LMXBs and HMXBs, normalized by SFR, takes a different shape (dashed purple) compared to the L19 relation: it reaches a minimum value at sSFR ∼ −10 after which it begins to approach the linear HMXB regime. Our sample is in agreement with the purple relation and is broadly consistent with the relations from L19 for logsSFR ≲ −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. We also compare model 5 of Aird et al. (2017; solid cyan) in which they quantify the total L_X of a star-forming galaxy as a function of M_*, SFR and redshift. They take into account their whole 0.1 < z < 2 sample of galaxies and convert the total count rate of each galaxy from the 0.5–2 keV band to the 2–10 keV band using a constant conversion factor based on an absorbed PL with Γ = 1.9. It takes the form

$$\begin{array}{c}\hfill L_X=\alpha {(1+z)}^{\gamma }{M}_{\ast }+\beta {(1+z)}^{\delta }{\mathrm{SFR}}^{\theta },\end{array}$$(20)

Fig. 9. Relation between SFR-normalized total X-ray luminosity (both XRB types) and sSFR, for star-forming galaxies in Box4of Magneticum (contours: percentiles). The solid black line is the relation of observed normal galaxies at z ∼ 0 from L16. The dash-dotted line represents the relation for a set of local galaxies from L19. The solid cyan line is the relation obtained from model 5 (m5) in Aird et al. (2017), where we fixed z = 0.01. The dashed red line depicts the expected scaling relation for a pure LMXB contribution given our XLF choice. The dashed blue line shows the expected broken SFR scaling relation of a pure HMXB contribution. Data points are taken from S22 and M14.

with log α[erg s⁻¹]=28.81, log β[erg s⁻¹]=39.5, γ = 3.9, δ = 0.67, θ = 0.86, M_* in M_⊙ and SFR in erg s⁻¹. We fixed the redshift to z = 0.01 in line with our sample. Due to the nonlinear SFR dependence, model 5 falls off for higher SFR values. From a visual inspection alone, both our broken SFR scaling relation and model 5 appear to adequately describe the SFR-normalized luminosities from S22 at high sSFR. However, model 5 underpredicts luminosities in the low sSFR regime compared to our simulated data. In their study S22 show that most of their data points with logsSFR 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.

6.4. Metallicity relation of HMXB

Numerous studies on the total HMXB emission in individual galaxies associate an increase in HMXB emissivity per unit SFR with a lower-metallicity environment of the HMXB. (Vulic et al. 2022; Saxena et al. 2021; L21; Garofali et al. 2020; Fragos et al. 2013b). In this section we try to verify if our modeling allows us to capture the metallicity dependence of the HMXB luminosity function proposed by L21 (see Eq. (11)). They obtained SFRs similar to L19 using far-UV and 24 μm maps and gas-phase metallicities from emission-weighted oxygen line ratios. In Fig. 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). We binned our sample according to metallicity and SFR with horizontal bars representing the bin width in metallicity and vertical bars representing the 25th–75th 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 with respect to the metallicity bin center for illustrative purposes. Our full sample (black crosses with error bars) is consistent with the theoretical model and shows the expected increase in luminosity for lower metallicities. Similarly, the high SFR sample (≳5 M_⊙ yr⁻¹, 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 SFR ∼ 1 M_⊙ yr⁻¹ (green) we find lower median luminosities compared to the global model but are consistent within the 25th percentile margin. From the Markov chain Monte Carlo simulation performed in L21, it is however expected that the median luminosity decreases with lower SFR values due to worse HMXB 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 SFR ≲ 0.1. Given our approach to retrieve HMXB 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 HMXBs. Additionally, only luminous SSPs will be sampled by PHOX (see Fig. 5). Since we did not take into account galaxies that do not have HMXB emission despite being star-forming, we are biased toward a low SFR sample with comparatively large luminosities and too large HMXB population.

Fig. 10. Metallicity dependence on the linear SFR scaling relation of HMXB 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 luminosity function sampling. The shaded gray area is its 16th-84th percentile margin. Our data are binned by oxygen fraction and SFR; SFR bins are indicated by color. Vertical error bars correspond to the 25th–75th percentile margin of the data within each metallicity bin, and horizontal black bars indicate the width of the metallicity bin. Colored symbols corresponding to the SFR bins are displaced horizontally with respect to the metallicity bin center for illustrative purposes. The dashed purple line is the corresponding M12b model where no explicit metallicity dependence is assumed.

7. Relative emission of hot gas and XRBs

In this section we used the ideal photon list generated with PHOX to obtain relative count ratios between the gas and XRB component in our galaxy sample. Additionally, we built 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 afterward 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 compared the respective photon counts coming from each component. We split the contribution into different energy bands (0.5–2 keV, 2–8 keV, and 0.5–8 keV) and binned our galaxies according to the count ratio, r, which is defined as

$$\begin{array}{c}\hfill r=1+\frac{c_{\mathrm{XRB}}}{c_{\mathrm{GAS}}},\end{array}$$(21)

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 Fig. 11 we show the differential and cumulative fraction of galaxies in our sample with respect to the ratio r. We show the LMXB contribution, the HMXB contribution, and the summed contribution of LMXBs and HMXBs. Different line styles correspond to the different energy bands. Regarding the combined contribution of both XRB types we notice that ∼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 LMXB emission. One has to keep in mind, however, that HMXBs only occur in galaxies with sufficient star-formation. Thus, the low HMXB 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 ≲40% of our galaxies are still dominated by XRB emission.

Fig. 11. X-ray photon count ratio between the XRB and the hot gas component including all galaxies in our sample. Colors represent the component and line styles the energy range in the observed frame. The label “XRB” shows the sum of HMXB and LMXB emission. Top: differential fraction of galaxies with a given count ratio, r. Bottom: cumulative fraction of galaxies. Here, a ratio of r ≥ 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 Fig. 12, which shows the total average spectrum, the average combined spectrum of both XRB types, and the average spectrum of hot gas. In the bottom panel of Fig. 12 we show the relative contribution of each component toward the total spectrum as a function of energy. Additionally, we added the average spectrum from 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 L16 in the lower panel of Fig. 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_H^{\text{gas}}=5\times {10}^{22}{\text{cm}}^{\text{-3}}$ but use a weaker absorption coefficient $N_H^{\text{xrb}}=6\times {10}^{20}{\text{cm}}^{\text{-3}}$ compared to our 2 × 10²¹ cm⁻³. Because of this, the XRB contribution is much higher for low energies. The sensitive behavior 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 ∼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 that are similar to our values, again highlighting the strong model dependence of the XRB contribution in the soft X-ray. Given our simple choice for the XRB emission model in the form of a single absorbed PL, we fail to capture the expected steepening of the XRB PL for E > 6 keV that is predicted by L16 for normal galaxies, by population synthesis performed in Fragos et al. (2013b), and by, for example, Persic & Rephaeli (2002) for star-burst galaxies. A more sophisticated model of the XRB spectrum would be needed to fully reproduce observed spectral shapes.

Fig. 12. SFR-normalized average X-ray spectrum of our star-forming sample (SFR > 0.4 M⊙ yr−1). The combined spectrum (black) is made up of the HMXB (dark green), the LMXB (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 Fig. 13 we show the same decomposition of galactic spectra as in Fig. 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 that generally exceeds estimates from L16. By comparing to Fig. 6, where we found galaxies with SFR ≲ 1 M_⊙ yr⁻¹ having low L_X within the 99th percentile contour, we attribute the strong XRB component to under-luminous gas-poor galaxies in the probed SFR regime and normalization issues when determining ${\stackrel{\sim }{N}}_{\mathrm{HMXB}}$ from Eq. (10), which leads to more HMXBs per unit SFR. In the next higher SFR bin (top right) the normalization of the spectrum drops slightly, mainly in the HMXB 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 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 that was derived using high redshift analogs for normal galaxies.

Fig. 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) 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 that contribute to the total spectrum: hot gas emission (blue), LMXB (lime green), and HMXB (dark green). “XRB” (orange) stands for LMXB+HMXB 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 gray box highlights the bin size and the number of galaxies, N, within each bin.

Since the two lowest SFR bins present the largest contribution in terms of sample number, the total average spectrum from Fig. 12 will mostly resemble to the upper panels of Fig. 13. An important question associated with galactic X-ray spectra is their influence on the cosmic X-ray background. Upcoming planned X-ray missions such as Athena or the Light Element Mapper are sensitive enough to receive significant contributions from galactic emission to their background calibration. We plan to investigate this issue in an upcoming study.

8. Discussion

8.1. Gas luminosity scaling relation

Compared to observational data we found that the L_X − SFR relation in the 0.5–2 keV band for the hot gas component has a large intrinsic scatter of ∼1 dex at SFR ≲ 1 with observers typically reporting lower luminosities of L_X ∼ 10³⁹ × SFR with significantly smaller scatter (S04, M12a, B13). We propose six explanations as to why we expected the scatter to be larger for the hot gas component.

First, 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⁷ 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₂₅₀₀. Furthermore, numerous studies explicitly mention that intrinsic absorption is highly dependent on the observed galaxy and may vary significantly (e.g., M12a; Lehmer et al. 2016, 2019, 2022).

Second, 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 explain in one of the following points.

Third, recent studies have shown that the extended hot gas atmosphere can significantly impact the total observed luminosity depending on the chosen projection radius (Nica et al. 2022; Comparat et al. 2022). In Comparat et al. (2022) they studied the X-ray emission of the circumgalactic medium for halos of stellar mass 2 × 10¹⁰ − 10¹² M_⊙ 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.

Fourth, regarding the influence of AGN toward 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 behavior compared to normal galaxies with respect to their X-ray radial profiles and scaling relations (Comparat et al. 2019; Vulic et al. 2022). A proper disentanglement of AGN contribution to these relations goes beyond the scope of this paper. We prefer to more thoroughly analyze their effects in a future study.

Fifth, 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₂₅₀₀, 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.

Finally, the most driving factor for uncertainty, especially toward low L_X, is probably the resolution limit of the total gas mass in the cosmological simulation. SPH simulations usually derive hydrodynamical quantities by interpolating between neighboring SPH elements with a set number of neighbors. 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 neighbor 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 ∼1 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 behavior 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 − SFR relation of HMXBs is a major concern. While the expected Poissonian noise from discrete sampling of the luminosity functions was predicted by Gilfanov et al. (2004a), additional uncertainties arise from our SFH proxy (see Sect. 3). For the birthrate of massive stars (Eq. (9)) we assumed a linear dependence on SFR that 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 in cosmological simulations, this first-order approximation is the best we can achieve. For the lifetime function of stars we assumed the prescription of 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 despite it being known that the 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.

One of the most important drivers of uncertainty in the L_X − 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 of having SSPs in the required age range of τ_SSP ∈ [0, 30] 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 − 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 − SFR relation is still consistent with observations.

In our investigation of the metallicity dependence we concluded that there is only a weak link between the SFR and the mass-weighted stellar metallicity in the simulation since the HMXB model from M12b 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 in Fig. 10, in contrast to observations by L21. We attribute this behavior of the M12b 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. We note that our low SFR sample (yellow) for 12 + log[O/H] ≤ 8.5 has much higher luminosities than predicted by the Monte Carlo Markov chains of L21 for SFR ∼ 0.1. Apart from the number normalization issues discussed above for low resolution SSPs, we speculate that this difference might be caused by either a different zero-point in metallicity or differences between mass-weighted stellar metallicities and gas-phase metallicities, the latter of which 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 Fig. 9 we compare the L_X − SFR − M_* relation of our sample with recent observations of galaxies within the Virgo cluster from S22. We find excellent agreement with our data. Due to the break in our L_X − SFR relation for HMXBs our expected transition of LMXB to HMXB dominated galaxies, around logsSFR ≂ −10.5), allows for slightly lower L_X compared to a linear scaling in L_X − SFR. The normalization of our L_X − 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 − SFR − M_* relation and is equally consistent to the model in Eq. (19) if we consider a lower normalization for the L_X − SFR relation. Interestingly, model 5 of Aird et al. (2017) shows a sublinear normalization in the SFR component that also leads to a much better representation of the S22 data compared to the L19 model. They argue that the nonlinear SFR dependence may arise from differences in the properties of the stellar population, such as age and metallicity (see, e.g., Fragos et al. 2013a), across a given mass range of galaxies. Combining galaxies of different redshifts thus neglects changes in stellar properties contributing to differences in the total X-ray emission per unit SFR. A similar line of argument is raised by S22 who claim that their observed incompatibility with the L19 models is caused either by the galaxy cluster environment or by unaccounted properties such as metallicity.

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 ≤ 2 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 ISM maps of the host galaxies into account 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, for example 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 analogs 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 Fig. 6). An increase of gas luminosity for lower-metallicity environments 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 ∼ 1 − 2 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 Fig. 12, which leads in our case to overestimation of the XRB component compared to values from L16. 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 as well as a higher normalization in general. This is due to the influence from massive 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. Gilbertson et al. (2022) show 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 that 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 circumgalactic medium.

9. Summary and conclusions

Based on SNII rate estimates made 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 have 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_⊙ ≤ M_* ≤ 10¹² M_⊙) and produced ideal photon lists for the hot gas and the XRB component using the stochastic photon simulator PHOX. We assumed XRB emission followed a static PL with photon index Γ = 1.7 − 2, while hot gas was modeled 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 required an intrinsic absorption component for both hot gas and the XRB component with column density $N_H^{\mathrm{gas}}=5\times {10}^{21}{\mathrm{cm}}^{-2}$ and $N_H^{\mathrm{xrb}}=2\times {10}^{21}{\mathrm{cm}}^{-2}$, respectively. For each halo in our sample, we constructed cylindrical volumes around the center of mass with R = R₂₅₀₀ and h = 2 R₂₅₀₀ and projected the enclosed photons, accounting for the redshift and peculiar motion of the source with respect to the line of sight. Our findings from the resulting galactic X-ray spectra can be summarized as follows:

1. The global reconstructed XRB luminosity functions perfectly resemble the underlying analytic shape. The observed flattening of the reconstructed XLFs below log L_XRB ≈ 37 is consistent with the expected one-photon luminosity. The reconstructed XLF derived from the metallicity-dependent model from L21 most closely resembles the XLF shape for solar metallicity, which is almost indistinguishable from the M12b model. A more curated sample of galaxies would be needed to reliably disentangle the effects of metallicity on the XLF shape, which may be analyzed in a separate study.

2. The L_X-M_* relation for LMXBs is tightly constrained and is in excellent agreement with the Z12 relation. More recent measurements of this relation by Lehmer et al. (2020) are also consistent with our data.

3. The L_X-SFR relation for HMXBs suffers greatly from low-number statistics for SFR ≲ 1 M_⊙ yr⁻¹. Additionally, the SFR proxy used to determine the number of HMXBs 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 BPL relation, which is expected to arise from incomplete sampling of the HMXB XLF dominated by Poissonian noise (Gilfanov et al. 2004a).

4. The L_X-SFR-Z relation as proposed by L21 does not arise when sampling from the M12b 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 HMXB sampling.

5. The L_X/SFR-sSFR relation we recovered from studying the combined effect of both XRB types is in remarkable agreement with recent studies of field galaxies within the Virgo cluster from S22. The shape of our derived relation shows more complex behavior at logsSFR ∼ −10 compared to the conventional relation found in the literature (Lehmer et al. 2010, 2016, 2019) due to the observed break in the HMXB relation.

6. As expected, we find that HMXBs spatially coincide with the star-forming regions of their host galaxies and that the spatial LMXB distribution follows the stellar-mass surface density.

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 HMXBs to these percentages reflects the relatively small portion of high SFR galaxies within our sample.

8. We constructed average SFR-normalized galactic X-ray spectra from our complete star-forming sample (SFR > 0.4) and measured the relative contribution of the two XRB types and the hot gas component toward the total spectrum. The average spectrum of the total star-forming sample is consistent with the spectral fits from 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 L16 but consistent with those from Wik et al. (2014), which reflects different model assumptions for absorption and XRB 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 have presented an approach to modeling galaxy X-ray spectra from cosmological simulations. In addition to the hot gas component in the simulation, we took the properties of stellar elements into account and implemented a fast and consistent algorithm to compute XRB emission from these elements while reliably separating the two different XRB types. The presented method proves its viability as a supporting tool for observers with its highly modular design, which also allows for alterations to employed model spectra for individual components. Future studies that build upon the suggested implementation will improve the sample selection and focus on the influence of X-ray AGN emission on scaling relations, on the background, and on the retrieved spectra.

Acknowledgments

We thank the anonymous referee for their helpful comments and feedback which helped improving this paper. S.V.Z. and V.B. acknowledge support by the Deutsche Forschungsgemeinschaft, DFG project nr. 415510302. K.D. 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

All Figures

	Fig. 1. Mgas − M* relation of our simulated galaxy sample. The dashed diagonal line indicates a one-to-one ratio.
In the text

Fig. 2. SFR − M* relation of our sample. The dashed diagonal lines correspond to constant logsSFR 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 the classification scheme from Schulze et al. (2020; red: elliptical, green: intermediate, blue: disk). Galaxies clustered at logsSFR = −12 were given a SFR estimate based on the Weinmann et al. (2010) approach (see text).

In the text

	Fig. 3. Projected stellar surface mass density (Σ*) of the poster-child disk galaxy in Box4 of Magneticum with halo ID 13633. Contours indicate the 68th, 90th, and 99th percentiles of X-ray photon positions for HMXBs (black) and LMXBs (red).
In the text

	Fig. 4. Mass-weighted stellar age (τ*) of the poster-child disk galaxy in Box4 of Magneticum with halo ID 13633. Contours indicate the 68th, 90th, and 99th percentiles of X-ray photon positions for HMXBs (black) and LMXBs (red).
In the text

Fig. 5. Integrated XLFs from Eqs. (1), (3), and (11) normalized by their respective scaling (SFR and M*). The dash-dot-dotted line is the XLFs of the L21 model, for log[O/H]+12 = 7.2, and the dotted line for log[O/H]+12 = 9.2. Colored lines are average XLFs obtained from the luminosity sampling in PHOX: from M12b in dashed blue, L21 in dash-dotted orange, and Z12 in solid green. Corresponding lines share the same line style. The gray vertical line at log L ∼ 37 marks the maximum luminosity of a source emitting one single photon with an average energy of 0.5–8 keV.

In the text

Fig. 6. Relation between hot gas X-ray luminosity (0.5–2 keV) and the SFR for star-forming galaxies in the simulations. We only show galaxies with log(sSFR) >  − 11.6, which removes all galaxies that fulfill the Weinmann et al. (2010) criterion. This is done to reduce the influence of unresolved galaxies on the relation as well as remove 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σ standard deviation. Data points are bolometric luminosities of the diffuse emission taken from S04 in the 0.3–2 keV band, M12a in the 0.5–2 keV band, and B13 in the 0.5–2 keV band.

In the text

Fig. 7. Same as Fig. 6, but for HMXB luminosity in the 0.5–8 keV band. The filled blue contours represent the same halos with HMXBs sampled from the L21 model (Eq. (11)). The solid black line is the expected relation from the employed HMXB model in Eq. (3), and the dashed red line is the predicted relation for discretely sampling data points from a PL distribution (Gilfanov et al. 2004a). The solid blue line is the expected relation from the employed L21 model. Data points are taken from M12b.

In the text

	Fig. 8. Same as Fig. 7, but for LMXB emission in the 0.5–8 keV band. The solid black line is the expected relation from the employed LMXB model in Eq. (1). The dashed black line is the same relation for a luminosity limit Lmin of 5 × 1037 erg s−1. Data points (cyan triangles) are taken from Lehmer et al. (2020).
In the text

Fig. 9. Relation between SFR-normalized total X-ray luminosity (both XRB types) and sSFR, for star-forming galaxies in Box4of Magneticum (contours: percentiles). The solid black line is the relation of observed normal galaxies at z ∼ 0 from L16. The dash-dotted line represents the relation for a set of local galaxies from L19. The solid cyan line is the relation obtained from model 5 (m5) in Aird et al. (2017), where we fixed z = 0.01. The dashed red line depicts the expected scaling relation for a pure LMXB contribution given our XLF choice. The dashed blue line shows the expected broken SFR scaling relation of a pure HMXB contribution. Data points are taken from S22 and M14.

In the text

Fig. 10. Metallicity dependence on the linear SFR scaling relation of HMXB 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 luminosity function sampling. The shaded gray area is its 16th-84th percentile margin. Our data are binned by oxygen fraction and SFR; SFR bins are indicated by color. Vertical error bars correspond to the 25th–75th percentile margin of the data within each metallicity bin, and horizontal black bars indicate the width of the metallicity bin. Colored symbols corresponding to the SFR bins are displaced horizontally with respect to the metallicity bin center for illustrative purposes. The dashed purple line is the corresponding M12b model where no explicit metallicity dependence is assumed.

In the text

Fig. 11. X-ray photon count ratio between the XRB and the hot gas component including all galaxies in our sample. Colors represent the component and line styles the energy range in the observed frame. The label “XRB” shows the sum of HMXB and LMXB emission. Top: differential fraction of galaxies with a given count ratio, r. Bottom: cumulative fraction of galaxies. Here, a ratio of r ≥ 2 means that the XRB component is dominant. Ratios exceeding the given range are collected in the last bin.

In the text

	Fig. 12. SFR-normalized average X-ray spectrum of our star-forming sample (SFR > 0.4 M⊙ yr−1). The combined spectrum (black) is made up of the HMXB (dark green), the LMXB (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.
In the text

Fig. 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) 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 that contribute to the total spectrum: hot gas emission (blue), LMXB (lime green), and HMXB (dark green). “XRB” (orange) stands for LMXB+HMXB 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 gray box highlights the bin size and the number of galaxies, N, within each bin.

In the text