diff options
| -rw-r--r-- | fg21sim/extragalactic/clusters/halo.py | 698 |
1 files changed, 698 insertions, 0 deletions
diff --git a/fg21sim/extragalactic/clusters/halo.py b/fg21sim/extragalactic/clusters/halo.py new file mode 100644 index 0000000..26ff6de --- /dev/null +++ b/fg21sim/extragalactic/clusters/halo.py @@ -0,0 +1,698 @@ +# Copyright (c) 2017 Weitian LI <liweitianux@live.com> +# MIT license + +""" +Simulate (giant) radio halos following the "statistical +magneto-turbulent model" proposed by Cassano & Brunetti (2005). + +References +---------- +[1] Cassano & Brunetti 2005, MNRAS, 357, 1313 + http://adsabs.harvard.edu/abs/2005MNRAS.357.1313C +""" + +import logging + +import numpy as np +import astropy.units as au +import astropy.constants as ac +import scipy.interpolate +import scipy.integrate +import scipy.optimize + +from .cosmology import Cosmology +from .solver import FokkerPlanckSolver + + +logger = logging.getLogger(__name__) + + +class HaloSingle: + """ + Simulate a single (giant) radio halos following the "statistical + magneto-turbulent model" proposed by Cassano & Brunetti (2005). + + First, simulate the cluster merging history from the extended + Press-Schecter formalism using the Monte Carlo method; then derive + the merger energy and turbulence energy as well as its spectrum; + after that, calculate the electron acceleration and time evolution + by solving the Fokker-Planck equation; and finally derive the radio + emission from the electron spectra. + + References + ---------- + [1] Cassano & Brunetti 2005, MNRAS, 357, 1313 + http://adsabs.harvard.edu/abs/2005MNRAS.357.1313C + + Parameters + ---------- + M0 : float + Present-day (z=0) mass (unit: Msun) of the cluster. + configs : `ConfigManager` + A `ConfigManager` instance containing default and user configurations. + For more details, see the example configuration specifications. + + Attributes + ---------- + mec : float + Unit for electron momentum (p): mec = m_e * c, p = gamma * mec, + therefore value of p is the Lorentz factor. + cosmo : `~Cosmology` + Adopted cosmological model with custom utility functions. + mtree : `~MergerTree` + Merging history of this cluster. + """ + # Merger tree (i.e., merging history) of this cluster + mtree = None + # Unit for electron momentum (p), thus its value is the Lorentz factor + mec = ac.m_e.cgs.value*ac.c.cgs.value # [g cm / s] + # Mean molecular weight + # Ref.: Ettori et al, 2013, Space Science Review, 177, 119-154, Eq.(6) + mu = 0.6 + # Atomic mass unit (i.e., a.m.u.) + m_atom = ac.u.cgs.value # [g] + # Common units conversion + # TODO: move these to a separate module/class + Msun2g = au.solMass.to(au.g) + kpc2cm = au.kpc.to(au.cm) + keV2erg = au.keV.to(au.erg) + Gyr2s = au.Gyr.to(au.s) + + def __init__(self, M0, configs): + self.M0 = M0 # [Msun] + self.configs = configs + self._set_configs() + + def _set_configs(self): + """ + Set up the necessary class attributes according to the configs. + """ + comp = "extragalactic/halos" + self.zmax = self.configs.getn(comp+"/zmax") + # Mass threshold of the sub-cluster for a significant merger + self.merger_mass_th = self.configs.getn(comp+"/merger_mass_th") + self.radius_halo = self.configs.getn(comp+"/radius_halo") + self.magnetic_field = self.configs.getn(comp+"/magnetic_field") + self.eta_t = self.configs.getn(comp+"/eta_t") + self.eta_e = self.configs.getn(comp+"/eta_e") + self.pmin = self.configs.getn(comp+"/pmin") + self.pmax = self.configs.getn(comp+"/pmax") + self.pgrid_num = self.configs.getn(comp+"/pgrid_num") + self.time_step = self.configs.getn(comp+"/time_step") + self.injection_index = self.configs.getn(comp+"/injection_index") + # Cosmology model + self.H0 = self.configs.getn("cosmology/H0") + self.OmegaM0 = self.configs.getn("cosmology/OmegaM0") + self.cosmo = Cosmology(H0=self.H0, Om0=self.OmegaM0) + logger.info("Loaded and set up configurations") + + def simulate_mergertree(self): + """ + Simulate the merging history of the cluster using the extended + Press-Schechter formalism. + """ + raise NotImplementedError + + def calc_electron_spectrum(self): + """ + Calculate the relativistic electron spectrum by solving the + Fokker-Planck equation. + """ + fpsolver = FokkerPlanckSolver( + xmin=self.pmin, xmax=self.pmax, + grid_num=self.pgrid_num, + tstep=self.time_step, + f_advection=self.fp_advection, + f_diffusion=self.fp_diffusion, + f_injection=self.fp_injection, + ) + p = fpsolver.x + # Assume NO initial electron distribution + n0_e = np.zeros(p.shape) + tstart = self.cosmo.age(self.zmax) + tstop = self.cosmo.age0 + n_e = fpsolver.solve(u0=n0_e, tstart=tstart, tstop=tstop) + return (p, n_e) + + def kT_mass(self, mass): + """ + Estimate the cluster ICM temperature from its mass by assuming + an (observed) temperature-mass relation. + + TODO: upgrade this M-T relation. + + Parameters + ---------- + mass : float + Mass (unit: Msun) of the cluster + + Returns + ------- + kT : float + Temperature of the ICM (unit: keV) + + References + ---------- + [1] Nevalainen et al. 2000, ApJ, 532, 694 + Ettori et al, 2013, Space Science Review, 177, 119-154 + NOTE: H0 = 50 * h50 [km/s/Mpc] + """ + kT = 10 * (mass/1.23e15) ** (1/1.79) # [keV] + return kT + + def _radius_virial(self, mass, z=0.0): + """ + Calculate the virial radius of a cluster. + + Parameters + ---------- + mass : float + Mass (unit: Msun) of the cluster + z : float + Redshift + + Returns + ------- + Rvir : float + Virial radius (unit: kpc) of the cluster at given redshift + """ + Dc = self.cosmo.overdensity_virial(z) + rho = self.cosmo.rho_crit(z) # [g/cm^3] + R_vir = (3*mass*self.Msun2g / (4*np.pi * Dc * rho))**(1/3) # [cm] + R_vir /= self.kpc2cm # [kpc] + return R_vir + + def _radius_stripping(self, mass, M_main, z): + """ + Calculate the stripping radius of the sub-cluster at which + equipartition between static and ram pressure is established, + and the stripping is efficient outside this stripping radius. + + Note that the value of the stripping radius obtained would + give the *mean value* of the actual stripping radius during + a merger. + + Parameters + ---------- + mass : float + The mass (unit: Msun) of the sub-cluster. + M_main : float + The mass (unit: Msun) of the main cluster. + z : float + Redshift + + Returns + ------- + rs : float + The stripping radius of the sub-cluster. + Unit: kpc + + References + ---------- + [1] Cassano & Brunetti 2005, MNRAS, 357, 1313 + http://adsabs.harvard.edu/abs/2005MNRAS.357.1313C + Eq.(11) + """ + vi = self._velocity_impact(M_main, mass, z) * 1e5 # [cm/s] + kT = self.kT_mass(mass) * self.keV2erg # [erg] + coef = kT / (self.mu * self.m_atom * vi**2) # dimensionless + rho_avg = self._density_average(M_main, z) # [g/cm^3] + + def equation(r): + return coef * self.density_profile(r, mass, z) / rho_avg - 1 + + r_vir = self._radius_virial(mass, z) # [kpc] + rs = scipy.optimize.brentq(equation, a=0, b=r_vir) # [kpc] + return rs + + def _density_average(self, mass, z=0.0): + """ + Average density of the cluster ICM. + + Returns + ------- + rho : float + Average ICM density (unit: g/cm^3) + + References + ---------- + [1] Cassano & Brunetti 2005, MNRAS, 357, 1313 + http://adsabs.harvard.edu/abs/2005MNRAS.357.1313C + Eq.(12) + """ + f_baryon = self.cosmo.Ob0 / self.cosmo.Om0 + Rv = self._radius_virial(mass, z) * self.kpc2cm # [cm] + V = (4*np.pi / 3) * Rv**3 # [cm^3] + rho = f_baryon * mass*self.Msun2g / V # [g/cm^3] + return rho + + def density_profile(self, r, mass, z): + """ + ICM (baryon) density profile, assuming the beta model. + + Parameters + ---------- + r : float + Radius (unit: kpc) where to calculate the density + mass : float + Cluster mass (unit: Msun) + z : float + Redshift + + Returns + ------- + rho_r : float + Density at the specified radius (unit: g/cm^3) + + References + ---------- + [1] Cassano & Brunetti 2005, MNRAS, 357, 1313 + http://adsabs.harvard.edu/abs/2005MNRAS.357.1313C + Eq.(13) + """ + f_baryon = self.cosmo.Ob0 / self.cosmo.Om0 + M_ICM = mass * f_baryon * self.Msun2g # [g] + r *= self.kpc2cm # [cm] + Rv = self._radius_virial(mass, z) * self.kpc2cm # [cm] + rc = self._beta_rc(Rv) + beta = self._beta_beta() + norm = self._beta_norm(M_ICM, beta, rc, Rv) # [g/cm^3] + rho_r = norm * (1 + (r/rc)**2) ** (-3*beta/2) # [g/cm^3] + return rho_r + + @staticmethod + def _beta_rc(r_vir): + """ + Core radius of the beta model for the ICM density profile. + + TODO: upgrade this! + """ + return 0.1*r_vir + + @staticmethod + def _beta_beta(): + """ + Beta value of the beta model for the ICM density profile. + + TODO: upgrade this! + """ + return 0.8 + + @staticmethod + def _beta_norm(mass, beta, rc, r_vir): + """ + Calculate the normalization of the beta model for the ICM + density profile. + + Parameters + ---------- + mass : float + The mass (unit: g) of ICM + beta : float + Beta value of the assumed beta profile + rc : float + Core radius (unit: cm) of the assumed beta profile + r_vir : float + The virial radius (unit: cm) of the cluster + + Returns + ------- + norm : float + Normalization of the beta model (unit: g/cm^3) + + References + ---------- + [1] Cassano & Brunetti 2005, MNRAS, 357, 1313 + http://adsabs.harvard.edu/abs/2005MNRAS.357.1313C + Eq.(14) + """ + integration = scipy.integrate.quad( + lambda r: r*r * (1+(r/rc)**2) ** (-3*beta/2), + 0, r_vir)[0] + norm = mass / (4*np.pi * integration) # [g/cm^3] + return norm + + def _velocity_impact(self, M_main, M_sub, z=0.0): + """ + Calculate the relative impact velocity between the two merging + clusters when they are at a distance of virial radius. + + Parameters + ---------- + M_main : float + Mass of the main cluster (unit: Msun) + M_sub : float + Mass of the sub cluster (unit: Msun) + z : float + Redshift + + Returns + ------- + vi : float + Relative impact velocity (unit: km/s) + + References + ---------- + [1] Cassano & Brunetti 2005, MNRAS, 357, 1313 + http://adsabs.harvard.edu/abs/2005MNRAS.357.1313C + Eq.(9) + """ + eta_v = 4 * (1 + M_main/M_sub) ** (1/3) + R_vir = self._radius_virial(M_main, z) * self.kpc2cm # [cm] + G = ac.G.cgs.value + vi = np.sqrt(2*G * (1-1/eta_v) * + (M_main+M_sub)*self.Msun2g / R_vir) # [cm/s] + vi /= 1e5 # [km/s] + return vi + + def _time_crossing(self, M_main, M_sub, z): + """ + Calculate the crossing time of the sub-cluster during a merger. + + Parameters + ---------- + M_main : float + Mass of the main cluster (unit: Msun) + M_sub : float + Mass of the sub cluster (unit: Msun) + z : float + Redshift where the merger occurs. + + Returns + ------- + time : float + Crossing time (unit: Gyr) + """ + R_vir = self._radius_virial(M_main, z) # [kpc] + vi = self._velocity_impact(M_main, M_sub, z) # [km/s] + # Unit conversion coefficient: [s kpc/km] => [Gyr] + # uconv = au.kpc.to(au.km) * au.s.to(au.Gyr) + uconv = 0.9777922216731284 + time = uconv * R_vir / vi # [Gyr] + return time + + def _z_end(self, z_begin, time): + """ + Calculate the ending redshift from ``z_begin`` after elapsing + ``time``. + + Parameters + ---------- + z_begin : float + Beginning redshift + time : float + Elapsing time (unit: Gyr) + """ + t_begin = self.cosmo.age(z_begin) # [Gyr] + t_end = t_begin + time + if t_end >= self.cosmo.age(0): + z_end = 0.0 + else: + z_end = self.cosmo.redshift(t_end) + return z_end + + @property + def merger_events(self): + """ + Trace only the main cluster, and filter out the significant + merger events. + + Returns + ------- + mevents : list[dict] + List of dictionaries that records all the merger events + of the main cluster. + NOTE: + The merger events are ordered by increasing redshifts. + """ + events = [] + tree = self.mtree + while tree: + if (tree.major and tree.minor and + tree.minor.node.mass >= self.merger_mass_th and + tree.major.node.z <= self.zmax): + events.append({ + "M_main": tree.major.node.mass, + "M_sub": tree.minor.node.mass, + "z": tree.major.node.z, + "age": tree.major.node.age + }) + tree = tree.major + return events + + def _coef_acceleration(self, z): + """ + Calculate the electron-acceleration coefficient at arbitrary + redshift, by interpolating the coefficients calculated at every + merger redshifts. + """ + if not hasattr(self, "_coef_acceleration_interp"): + # Order the merger events by decreasing redshifts + mevents = list(reversed(self.merger_events)) + redshifts = np.array([ev["z"] for ev in mevents]) + chis = np.array([self._chi_at_zidx(zidx, mevents) + for zidx in range(len(redshifts))]) + self._coef_acceleration_interp = scipy.interpolate.interp1d( + redshifts, chis, kind="linear", + bounds_error=False, fill_value=0.0) + logger.info("Interpolated acceleration coefficients w.r.t. z") + return self._coef_acceleration_interp(z) + + def _chi_at_zidx(self, zidx, mevents): + """ + Calculate electron-acceleration coefficient at the specified + merger event which is specified with a redshift index. + + Parameters + ---------- + zidx : int + Index of the redshift where to calculate the coefficient. + mevents : list[dict] + List of dictionaries that records all the merger events + of the main cluster. + NOTE: + The merger events should be ordered by increasing time + (or decreasing redshifts). + + Returns + ------- + chi : float + The calculated electron-acceleration coefficient. + (unit: Gyr^-1) + + References + ---------- + [1] Cassano & Brunetti 2005, MNRAS, 357, 1313 + http://adsabs.harvard.edu/abs/2005MNRAS.357.1313C + Eq.(40) + """ + redshifts = np.array([ev["z"] for ev in mevents]) + zbegin = mevents[zidx]["z"] + M_main = mevents[zidx]["M_main"] + M_sub = mevents[zidx]["M_sub"] + t_crossing = self._time_crossing(M_main, M_sub, zbegin) + zend = self._z_end(zbegin, t_crossing) + try: + zend_idx = np.where(redshifts < zend)[0][0] + except IndexError: + # Specified redshift already the last/smallest one + zend_idx = zidx + 1 + # + coef = 2.23e-16 * self.eta_t / (self.radius_halo/500)**3 # [s^-1] + coef *= self.Gyr2s # [Gyr^-1] + chi = 0.0 + for ev in mevents[zidx:zend_idx]: + M_main = ev["M_main"] + M_sub = ev["M_sub"] + z = ev["z"] + R_vir = self._radius_virial(M_main, z) + rs = self._radius_stripping(M_sub, M_main, z) + kT = self.kT_mass(M_main) + term1 = ((M_main+M_sub)/2e15 * (2.6e3/R_vir)) ** (3/2) + term2 = (rs/500)**2 / np.sqrt(kT/7) + if rs <= self.radius_halo: + term3 = 1.0 + else: + term3 = (self.radius_halo/rs) ** 2 + chi += coef * term1 * term2 * term3 + return chi + + def fp_injection(self, p, t=None): + """ + Electron injection term for the Fokker-Planck equation. + + The injected electrons are assumed to have a power-law spectrum + and a constant injection rate. + + Qe(p) = Ke * (p/pmin)**(-s) + Ke = ((s-2)*eta_e) * (e_th/(pmin*c)) / (t0*pmin) + + Parameters + ---------- + p : float + Electron momentum (unit: mec), i.e., Lorentz factor + t : None + Currently a constant injection rate is assumed, therefore + this parameter is not used. Keep it for the consistency + with other functions. + + Returns + ------- + Qe : float + Current electron injection rate at specified energy (p). + Unit: [cm^-3 Gyr^-1 mec^-1] + + References + ---------- + [1] Cassano & Brunetti 2005, MNRAS, 357, 1313 + http://adsabs.harvard.edu/abs/2005MNRAS.357.1313C + Eqs.(31-33) + """ + if not hasattr(self, "_electron_injection_rate"): + e_th = self.e_thermal # [erg/cm^3] + term1 = (self.injection_index-2) * self.eta_e + term2 = e_th / (self.pmin * self.mec * ac.c.cgs.value) # [cm^-3] + term3 = 1.0 / (self.cosmo.age0 * self.pmin) # [Gyr^-1 mec^-1] + Ke = term1 * term2 * term3 + self._electron_injection_rate = Ke + else: + Ke = self._electron_injection_rate + Qe = Ke * (p/self.pmin) ** (-self.injection_index) + return Qe + + def fp_diffusion(self, p, t): + """ + Diffusion term/coefficient for the Fokker-Planck equation. + + Parameters + ---------- + p : float + Electron momentum (unit: mec), i.e., Lorentz factor + t : float + Current time when solving the equation (unit: Gyr) + + Returns + ------- + Dpp : float + Diffusion coefficient + Unit: [mec^2/Gyr] + + References + ---------- + [1] Cassano & Brunetti 2005, MNRAS, 357, 1313 + http://adsabs.harvard.edu/abs/2005MNRAS.357.1313C + Eq.(36) + [2] Donnert 2013, AN, 334, 615 + http://adsabs.harvard.edu/abs/2013AN....334..515D + Eq.(15) + """ + z = self.cosmo.redshift(t) + chi = self._coef_acceleration(z) # [Gyr^-1] + # NOTE: Cassano & Brunetti's formula misses a factor of 2. + Dpp = chi * p**2 / 4 # [mec^2/Gyr] + return Dpp + + def fp_advection(self, p, t): + """ + Advection term/coefficient for the Fokker-Planck equation, + which describes a systematic tendency for upward or downard + drift of particles. + + This term is also called the "generalized cooling function" by + Donnert & Brunetti (2014), which includes all relevant energy + loss functions and the energy gain function due to turbulence. + + Returns + ------- + Hp : float + Advection coefficient + Unit: [mec/Gyr] + + References + ---------- + [1] Donnert & Brunetti 2014, MNRAS, 443, 3564 + http://adsabs.harvard.edu/abs/2014MNRAS.443.3564D + Eq.(15) + [2] Cassano & Brunetti 2005, MNRAS, 357, 1313 + http://adsabs.harvard.edu/abs/2005MNRAS.357.1313C + Eqs.(30,36,38,39) + """ + Hp = (abs(self._dpdt_ion(p, t)) + + abs(self._dpdt_rad(p, t)) - + (self.fp_diffusion(p, t) * 2 / p)) + return Hp + + def _dpdt_ion(self, p, t): + """ + Energy loss through ionization and Coulomb collisions. + + References + ---------- + [1] Cassano & Brunetti 2005, MNRAS, 357, 1313 + http://adsabs.harvard.edu/abs/2005MNRAS.357.1313C + Eq.(38) + """ + z = self.cosmo.redshift(t) + n_th = self._n_thermal(self.M0, z) + coef = -3.3e-29 * self.Gyr2s / self.mec # [mec/Gyr] + dpdt = coef * n_th * (1 + np.log(p/n_th) / 75) + return dpdt + + def _dpdt_rad(self, p, t): + """ + Energy loss via synchrotron emission and IC scattering off the CMB. + + References + ---------- + [1] Cassano & Brunetti 2005, MNRAS, 357, 1313 + http://adsabs.harvard.edu/abs/2005MNRAS.357.1313C + Eq.(39) + """ + z = self.cosmo.redshift(t) + coef = -4.8e-4 * self.Gyr2s / self.mec # [mec/Gyr] + dpdt = (coef * (p*self.mec)**2 * + ((self.magnetic_field/3.2)**2 + (1+z)**4)) + return dpdt + + @property + def e_thermal(self): + """ + Calculate the present-day thermal energy density of the ICM. + + Returns + ------- + e_th : float + Energy density of the ICM (unit: erg/cm^3) + """ + mass = self.M0 + f_baryon = self.cosmo.Ob0 / self.cosmo.Om0 + kT = self.kT_mass(mass) # [keV] + N = mass * self.Msun2g * f_baryon / (self.mu * self.m_atom) + E_th = kT*self.keV2erg * N # [erg] + Rv = self._radius_virial(mass) * self.kpc2cm # [cm] + V = (4*np.pi / 3) * Rv**3 # [cm^3] + e_th = E_th / V # [erg/cm^3] + return e_th + + def _n_thermal(self, mass, z=0.0): + """ + Calculate the present-day number density of the ICM thermal plasma. + + Parameters + ---------- + mass : float + Mass (unit: Msun) of the cluster + z : float + Redshift + + Returns + ------- + n_th : float + Number density of the ICM (unit: cm^-3) + """ + f_baryon = self.cosmo.Ob0 / self.cosmo.Om0 + N = mass * self.Msun2g * f_baryon / (self.mu * self.m_atom) + Rv = self._radius_virial(mass, z) * self.kpc2cm # [cm] + V = (4*np.pi / 3) * Rv**3 # [cm^3] + n_th = N / V # [cm^-3] + return n_th |