xray.py 6.48 KB
Newer Older
Guang's avatar
Guang committed
1
2
3
4
5
6
7
8
9
10
11
# -*- coding: utf-8 -*-
# Copyright (C) 2013 Centre de données Astrophysiques de Marseille
# Copyright (C) 2013 Institute of Astronomy, University of Cambridge
# Copyright (C) 2014 University of Crete
# Licensed under the CeCILL-v2 licence - see Licence_CeCILL_V2-en.txt
# Author: Guang Yang

"""
X-ray module
=============================

12
This module implements the X-ray emission from the galaxy and AGN corona.
Guang's avatar
Guang committed
13
14
15
16
17
18
19
20
21
22
23
24
25

"""

from collections import OrderedDict

import numpy as np
import scipy.constants as cst

from . import SedModule

class Xray(SedModule):
    """X-ray emission

Guang's avatar
Guang committed
26
    This module computes the X-ray emission from the galaxy and AGN corona.
Guang's avatar
Guang committed
27
28
29
30
31
32
33

    """

    parameter_list = OrderedDict([
        ("gam", (
            "cigale_list()",
            "The photon index (Gamma) of intrinsic X-ray spectrum.",
34
            1.7
Guang's avatar
Guang committed
35
        )),
36
        ("max_dev_alpha_ox", (
Guang's avatar
Guang committed
37
            "cigale_list()",
38
39
40
41
42
43
            "Maximum deviation of alpha_ox from the empirical alpha_ox-Lnu_2500A relation (Just et al. 2007), "
            "i.e. |alpha_ox-alpha_ox(Lnu_2500A)| <= max_dev_alpha_ox. "
            "alpha_ox is the power-law slope connecting L_nu at rest-frame 2500 A and 2 keV, "
            "defined as alpha_ox = 0.3838*log(Lnu_2keV/Lnu_2500A), "
            "which is often modeled as a function of Lnu_2500A."
            "The alpha_ox-Lnu_2500A relation has a 1-sigma scatter of ~ 0.1.",
44
            0.2
Guang's avatar
Guang committed
45
46
47
48
49
50
51
        ))
    ])

    def _init_code(self):
        """Build the model for a given set of parameters."""

        self.gam = float(self.parameters["gam"])
52
        self.max_dev_alpha_ox = float(self.parameters["max_dev_alpha_ox"])
Guang's avatar
Guang committed
53
54

        # We define various constants necessary to compute the model. For
Guang's avatar
Guang committed
55
        # consistency, we define speed of light in units of nm s¯¹
Guang's avatar
Guang committed
56
        c = cst.c * 1e9
Guang's avatar
Guang committed
57
        self.c = c
58
        # Define wavelenght corresponding to some energy in units of nm.
Guang's avatar
Guang committed
59
        lam_0p5keV = c*cst.h / (5e2*cst.eV)
Guang's avatar
Guang committed
60
        lam_2keV = c*cst.h / (2e3*cst.eV)
61
        lam_300keV = c*cst.h / (3e5*cst.eV)
Guang's avatar
Guang committed
62
63
        self.lam_2keV = lam_2keV
        self.lam_10keV = c*cst.h / (1e4*cst.eV)
Guang's avatar
Guang committed
64
        # Define frequency corresponding to 2 keV in units of Hz.
Guang's avatar
Guang committed
65
        self.nu_2keV = c / lam_2keV
Guang's avatar
Guang committed
66

Guang's avatar
Guang committed
67
        # We define the wavelength range for the non thermal emission
68
        # corresponding to 0.4-1200 keV
69
        self.wave = np.logspace(-3, 0.5, 1000)
Guang's avatar
Guang committed
70
71
72
73

        # 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
Guang's avatar
Guang committed
74
        self.lumin_hotgas = self.wave**0 / (lam_0p5keV - lam_2keV)
Guang's avatar
Guang committed
75
76
77
78
        # 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)
79
        # Assuming gamma=2. for both components (Wu et al. 2008)
Guang's avatar
Guang committed
80
        # normalized such that L(2-10 keV) = 1.
81
        self.lumin_xrb = self.wave**-1 / np.log(5.)
Guang's avatar
Guang committed
82

Guang's avatar
Guang committed
83
        # We compute the unobscured AGN corona X-ray emission
84
        # The shape is power-law with high-E exp. cutoff
85
        # with cut-off E=300 keV (Ueda+2014; Aird+2015)
86
        self.lumin_corona = self.wave**(self.gam - 3.) * \
87
                                np.exp(-lam_300keV/self.wave)
88
89
        # Normaliz the SED at 2 keV
        self.lumin_corona /= lam_2keV**(self.gam - 3.) * \
90
                                np.exp(-lam_300keV/lam_2keV)
Guang's avatar
Guang committed
91
92
93
94
95
96
        # Calculate total AGN corona X-ray luminosity
        self.l_agn_xray_total = np.trapz(self.lumin_corona, x=self.wave)
        # Calculate 2-10 keV AGN corona X-ray luminosity
        # 2-10 keV corredponds to self.wave[598:798]
        self.l_agn_xray_2to10keV = np.trapz(self.lumin_corona[598:798], x=self.wave[598:798])

Guang's avatar
Guang committed
97
98

    def process(self, sed):
Guang's avatar
Guang committed
99
        """Add the X-ray contribution.
Guang's avatar
Guang committed
100
101
102
103
104
105

        Parameters
        ----------
        sed: pcigale.sed.SED object

        """
Guang's avatar
Guang committed
106
107
        # Stellar info.
        # star formation rate, units: M_sun/yr
Guang's avatar
Guang committed
108
        sfr = sed.info['sfh.sfr']
Guang's avatar
Guang committed
109
110
111
112
        # stellar mass, units: 1e10 M_sun
        mstar = sed.info['stellar.m_star']/1e10
        # log stellar age, units: Gyr
        logT = np.log10( sed.info['stellar.age_m_star']/1e3 )
Guang's avatar
Guang committed
113
114
        # metallicity, units: none
        Z = sed.info['stellar.metallicity']
115
116
        # alpha_ox, added internally
        alpha_ox = self.parameters['alpha_ox']
117
        # AGN 2500A intrinsic luminosity
Guang's avatar
Guang committed
118
        if 'agn.agn_intrin_Lnu_2500A' not in sed.info:
119
            sed.add_info('agn.agn_intrin_Lnu_2500A', 1., True)
120
        Lnu_2500A = sed.info['agn.agn_intrin_Lnu_2500A']
121

Guang's avatar
Guang committed
122
        # Add the configuration for X-ray module
123
124
        sed.add_module(self.name, self.parameters)
        sed.add_info("xray.gam", self.gam)
125
126
        sed.add_info("xray.max_dev_alpha_ox", self.max_dev_alpha_ox)
        sed.add_info("xray.alpha_ox", self.parameters['alpha_ox'])
127

Guang's avatar
Guang committed
128
129
130
131
132
133
        # 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 * \
Guang's avatar
Guang committed
134
            10**(33.28 - 62.12*Z + 569.44*Z**2 - 1833.80*Z**3 + 1968.33*Z**4)
Guang's avatar
Guang committed
135
136
137
        # Calculate 2-10 keV LMXB luminosities
        # Mezcua et al. 2018, Eq. 2
        l_lmxb_xray_2to10keV = mstar * \
Guang's avatar
Guang committed
138
            10**(33.276 - 1.503*logT - 0.423*logT**2 + 0.425*logT**3 + 0.136*logT**4)
Guang's avatar
Guang committed
139
140

        # Calculate L_lam_2keV from Lnu_2500A
141
        Lnu_2keV = 10**(alpha_ox/0.3838) * Lnu_2500A
Guang's avatar
Guang committed
142
        L_lam_2keV = Lnu_2keV * self.nu_2keV**2/self.c
143
        # Calculate total AGN corona X-ray luminosity
Guang's avatar
Guang committed
144
        l_agn_xray_total = self.l_agn_xray_total * L_lam_2keV
145
        # Calculate 2-10 keV AGN corona X-ray luminosity
Guang's avatar
Guang committed
146
        l_agn_xray_2to10keV = self.l_agn_xray_2to10keV * L_lam_2keV
Guang's avatar
Guang committed
147

Guang's avatar
Guang committed
148
149
150
151
        # 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)
152
153
        sed.add_info("xray.agn_Lx_total", l_agn_xray_total, True)
        sed.add_info("xray.agn_Lx_2to10keV", l_agn_xray_2to10keV, True)
Guang's avatar
Guang committed
154
        sed.add_info("xray.agn_Lnu_2keV", Lnu_2keV, True)
Guang's avatar
Guang committed
155
        # Add the SED components
Guang's avatar
Guang committed
156
157
158
159
        sed.add_contribution('xray.galaxy', self.wave,
                self.lumin_hotgas*l_hotgas_xray_0p5to2keV + \
                self.lumin_xrb*l_hmxb_xray_2to10keV + \
                self.lumin_xrb*l_lmxb_xray_2to10keV)
Guang's avatar
Guang committed
160
        sed.add_contribution('xray.agn_corona', self.wave, self.lumin_corona * L_lam_2keV)
Guang's avatar
Guang committed
161
162

# SedModule to be returned by get_module
Guang's avatar
Guang committed
163
Module = Xray