# Copyright (c) 2017 Weitian LI # MIT license """ Helper functions References ---------- .. [arnaud2005] Arnaud, Pointecouteau & Pratt 2005, A&A, 441, 893; http://adsabs.harvard.edu/abs/2005A%26A...441..893 .. [cassano2005] Cassano & Brunetti 2005, MNRAS, 357, 1313 http://adsabs.harvard.edu/abs/2005MNRAS.357.1313C .. [cassano2007] Cassano et al. 2007, MNRAS, 378, 1565; http://adsabs.harvard.edu/abs/2007MNRAS.378.1565C .. [cassano2012] Cassano et al. 2012, A&A, 548, A100 http://adsabs.harvard.edu/abs/2012A%26A...548A.100C .. [zandanel2014] Zandanel, Pfrommer & Prada 2014, MNRAS, 438, 124 http://adsabs.harvard.edu/abs/2014MNRAS.438..124Z """ import numpy as np from scipy import integrate from ...configs import CONFIGS from ...utils import COSMO from ...utils.units import (Units as AU, Constants as AC, UnitConversions as AUC) def radius_virial(mass, z=0.0): """ Calculate the virial radius of a cluster at a given redshift. Parameters ---------- mass : float Total (virial) mass of the cluster Unit: [Msun] z : float, optional Redshift Default: 0.0 (i.e., present day) Returns ------- R_vir : float Virial radius of the cluster Unit: [kpc] """ Dc = COSMO.overdensity_virial(z) rho = COSMO.rho_crit(z) # [g/cm^3] R_vir = (3*mass*AUC.Msun2g / (4*np.pi * Dc * rho))**(1/3) # [cm] R_vir *= AUC.cm2kpc # [kpc] return R_vir def radius_halo(mass, z=0.0): """ Calculate the radius of (giant) radio halo for a cluster. The halo radius is assumed to linearly scale with the virial radius, and is estimated by: R_halo = R_vir / 4 * halo radius is ~3-6 times smaller than the virial radius; Ref.[cassano2007],Sec.(1) * halo half radius is ~R500/4, therefore, R_halo ~ R_vir/4; Ref.[zandanel2014],Sec.(6.2) Parameters ---------- mass : float Total (virial) mass of the cluster Unit: [Msun] z : float, optional Redshift Default: 0.0 (i.e., present day) Returns ------- R_halo : float Radius of the (expected) giant radio halo Unit: [kpc] """ R_vir = radius_virial(mass=mass, z=z) # [kpc] R_halo = R_vir / 4.0 # [kpc] return R_halo def mass_to_kT(mass, z=0.0): """ Calculate the cluster ICM temperature from its mass using the mass-temperature scaling relation (its inversion used here) derived from observations. The following M-T scaling relation from Ref.[arnaud2005],Tab.2: M200 * E(z) = A200 * (kT / 5 keV)^α , where: A200 = (5.34 +/- 0.22) [1e14 Msun] α = (1.72 +/- 0.10) Its inversion: kT = (5 keV) * [M200 * E(z) / A200]^(1/α). NOTE: M200 (i.e., Δ=200) is used to approximate the virial mass. Parameters ---------- mass : float Total (virial) mass of the cluster. Unit: [Msun] z : float, optional Redshift of the cluster Returns ------- kT : float The ICM mean temperature. Unit: [keV] """ # A = (5.34 + np.random.normal(scale=0.22)) * 1e14 # [Msun] A = 5.34 * 1e14 # [Msun] # alpha = 1.72 + np.random.normal(scale=0.10) alpha = 1.72 Ez = COSMO.E(z) kT = 5.0 * (mass * Ez / A) ** (1/alpha) return kT def density_number_thermal(mass, z=0.0): """ Calculate the number density of the ICM thermal plasma. NOTE ---- This number density is independent of cluster (virial) mass, but (mostly) increases with redshifts. Parameters ---------- mass : float Mass of the cluster Unit: [Msun] z : float, optional Redshift Returns ------- n_th : float Number density of the ICM thermal plasma Unit: [cm^-3] """ N = mass * AUC.Msun2g * COSMO.baryon_fraction / (AC.mu * AC.u) R_vir = radius_virial(mass, z) * AUC.kpc2cm # [cm] volume = (4*np.pi / 3) * R_vir**3 # [cm^3] n_th = N / volume # [cm^-3] return n_th def density_energy_thermal(mass, z=0.0): """ Calculate the thermal energy density of the ICM. Returns ------- e_th : float Energy density of the ICM Unit: [erg cm^-3] """ n_th = density_number_thermal(mass, z) # [cm^-3] kT = mass_to_kT(mass, z) * AUC.keV2erg # [erg] e_th = (3.0/2) * kT * n_th return e_th def density_energy_electron(spectrum, gamma): """ Calculate the energy density of relativistic electrons. Parameters ---------- spectrum : 1D float `~numpy.ndarray` The number density of the electrons w.r.t. Lorentz factors Unit: [cm^-3] gamma : 1D float `~numpy.ndarray` The Lorentz factors of electrons Returns ------- e_re : float The energy density of the relativistic electrons. Unit: [erg cm^-3] """ e_re = integrate.trapz(spectrum*gamma*AU.mec2, gamma) return e_re def velocity_impact(M_main, M_sub, z=0.0): """ Estimate the relative impact velocity between the two merging clusters when they are at a distance of the virial radius. Parameters ---------- M_main, M_sub : float Total (virial) masses of the main and sub clusters Unit: [Msun] z : float, optional Redshift Returns ------- vi : float Relative impact velocity Unit: [km/s] References ---------- Ref.[cassano2005],Eq.(9) """ eta_v = 4 * (1 + M_main/M_sub) ** (1/3) R_vir = radius_virial(M_main, z) * AUC.kpc2cm # [cm] vi = np.sqrt(2*AC.G * (1-1/eta_v) * (M_main+M_sub)*AUC.Msun2g / R_vir) # [cm/s] vi /= AUC.km2cm # [km/s] return vi def time_crossing(M_main, M_sub, z=0.0): """ Estimate the crossing time of the sub cluster during a merger. NOTE: The crossing time is estimated to be τ ~ R_vir / v_impact. Parameters ---------- M_main, M_sub : float Total (virial) masses of the main and sub clusters Unit: [Msun] z : float, optional Redshift Returns ------- time : float Crossing time Unit: [Gyr] References ---------- Ref.[cassano2005],Sec.(4.1) """ R_vir = radius_virial(M_main, z) # [kpc] vi = velocity_impact(M_main, M_sub, z) # [km/s] # Unit conversion coefficient: [s kpc/km] => [Gyr] uconv = AUC.kpc2km * AUC.s2Gyr time = uconv * R_vir / vi # [Gyr] return time def magnetic_field(mass): """ Calculate the mean magnetic field strength according to the scaling relation between magnetic field and cluster mass. Parameters ---------- mass : float Cluster mass Unit: [Msun] Returns ------- B : float The mean magnetic field strength Unit: [uG] References ---------- Ref.[cassano2012],Eq.(1) """ comp = "extragalactic/clusters" b_mean = CONFIGS.getn(comp+"/b_mean") b_index = CONFIGS.getn(comp+"/b_index") M_mean = 1.6e15 # [Msun] B = b_mean * (mass/M_mean) ** b_index return B