galaxy X-ray emission added

pcigale/sed_modules/xray.py

@@ -9,8 +9,7 @@
 X-ray module
 =============================
 
-This module implements the intrinsic (unobscured) X-ray emission from
-AGN corona.
+This module implements the X-ray emission from the galaxy and AGN corona.
 
 """
 
@@ -24,7 +23,7 @@ from . import SedModule
 class Xray(SedModule):
     """X-ray emission
 
-    This module computes the X-ray emission from AGN corona.
+    This module computes the X-ray emission from the galaxy and AGN corona.
 
     """
 
@@ -52,17 +51,32 @@ class Xray(SedModule):
         # consistency, we define speed of light in units of nm s¯¹
         c = cst.c * 1e9
         # Define wavelenght corresponding to some energy in units of nm.
+        lam_0p5keV = c*cst.h / (5e2*cst.eV)
         lam_2keV = c*cst.h / (2e3*cst.eV)
         self.lam_2keV = lam_2keV
         lam_10keV = c*cst.h / (1e4*cst.eV)
         self.lam_10keV = lam_10keV
         lam_100keV = c*cst.h / (1e5*cst.eV)
-        # Define wavelenght corresponding to 2 keV in units of Hz.
+        # Define frequency corresponding to 2 keV in units of Hz.
+        nu_0p5keV = c / lam_0p5keV
         nu_2keV = c / lam_2keV
 
         # We define the wavelength range for the non thermal emission
         # corresponding to 0.4-1200 keV
         self.wave = np.logspace(-3, 0.5, 1000.)
+
+        # X-ray emission from galaxies: 1.hot-gas & 2.X-ray binaries
+        # 1.Hot-gas, assuming power-law index gamma=3
+        # normalized such that L(0.5-2 keV) = 1
+        self.lumin_hotgas = 1./ (lam_0p5keV - lam_2keV)
+        # 2. X-ray binaries (XRB)
+        # also have two components:
+        #   2.1 high-mass X-ray binaries (HMXB)
+        #   2.2 low-mass X-ray binaries (LMXB)
+        # Assuming gamma=2 for both components (Lehmer et al. 2016)
+        # normalized such that L(2-10 keV) = 1.
+        self.lumin_xrb = self.wave**-1 / np.log(5)
+
         # We compute the unobscured AGN corona X-ray emission
         # The shape is power-law with high-E exp. cutoff
         self.lumin_corona = self.wave**(self.gam - 3.) * \
@@ -84,15 +98,36 @@ class Xray(SedModule):
         sed: pcigale.sed.SED object
 
         """
+        # Stellar info.
+        # log stellar age, units: Gyr
+        logT = np.log10( sed.info['sfh.age']/1e3 )
+        # star formation rate, units: M_sun/yr
+        sfr = sed.info['sfh.sfr100Myrs']
+        # metallicity, units: none
+        Z = sed.info['stellar.metallicity']
+        # stellar mass, units: 1e10 M_sun
+        mstar = sed.info['sfh.mass_total']/1e10
         # AGN 2500A intrinsic luminosity
         if 'agn.agn_intrin_Lnu_2500A' not in sed.info:
             sed.add_info('agn.agn_intrin_Lnu_2500A', 1., True)
         Lnu_2500A = sed.info['agn.agn_intrin_Lnu_2500A']
 
+        # Add the configuration for X-ray module
         sed.add_module(self.name, self.parameters)
         sed.add_info("xray.gam", self.gam)
         sed.add_info("xray.alpha_ox", self.alpha_ox)
 
+        # Calculate 0.5-2 keV hot-gas luminosities
+        # Mezcua et al. 2018, Eq. 5
+        l_hotgas_xray_0p5to2keV = 8.3e31 * sfr
+        # Calculate 2-10 keV HMXB luminosities
+        # Mezcua et al. 2018, Eq. 1
+        l_hmxb_xray_2to10keV = sfr * \
+            10**(40.28 - 62.16*Z + 569.44*Z**2 - 1833.80*Z**3 + 1968.23*Z**4)
+        # Calculate 2-10 keV LMXB luminosities
+        # Mezcua et al. 2018, Eq. 2
+        l_lmxb_xray_2to10keV = mstar * \
+            10**(40.276 - 1.503*logT - 0.423*logT**2 + 0.425*logT**3 + 0.136*logT**4)
         # Calculate total AGN corona X-ray luminosity
         l_agn_xray_total = np.trapz(Lnu_2500A * self.lumin_corona,
                              x=self.wave)
@@ -102,10 +137,17 @@ class Xray(SedModule):
         l_agn_xray_2to10keV = np.trapz(Lnu_2500A * self.lumin_corona[lam_idxs],
                               x=self.wave[lam_idxs])
 
+        # Save the results
+        sed.add_info("xray.hotgas_Lx_0p5to2keV", l_hotgas_xray_0p5to2keV, True)
+        sed.add_info("xray.hmxb_Lx_2to10keV", l_hmxb_xray_2to10keV, True)
+        sed.add_info("xray.lmxb_Lx_2to10keV", l_lmxb_xray_2to10keV, True)
         sed.add_info("xray.agn_Lx_total", l_agn_xray_total, True)
         sed.add_info("xray.agn_Lx_2to10keV", l_agn_xray_2to10keV, True)
-        sed.add_contribution('xray_corona', self.wave,
-                             self.lumin_corona * Lnu_2500A)
+        # Add the SED components
+        sed.add_contribution('xray_hotgas', self.wave, self.lumin_hotgas * l_hotgas_xray_0p5to2keV)
+        sed.add_contribution('xray_hmxb', self.wave, self.lumin_xrb * l_hmxb_xray_2to10keV)
+        sed.add_contribution('xray_lmxb', self.wave, self.lumin_xrb * l_lmxb_xray_2to10keV)
+        sed.add_contribution('xray_corona', self.wave, self.lumin_corona * Lnu_2500A)
 
 
 # SedModule to be returned by get_module