# Copyright (c) 2017 Weitian LI <weitian@aaronly.me>
# MIT license
"""
Simulate (giant) radio halo originating from the last/ most recent
cluster-cluster major merger event, following the "statistical
magneto-turbulent model" proposed by [cassano2005]_, but with many
modifications and simplifications.
References
----------
.. [brunetti2011]
Brunetti & Lazarian 2011, MNRAS, 410, 127
http://adsabs.harvard.edu/abs/2011MNRAS.410..127B
.. [cassano2005]
Cassano & Brunetti 2005, MNRAS, 357, 1313
http://adsabs.harvard.edu/abs/2005MNRAS.357.1313C
.. [cassano2006]
Cassano, Brunetti & Setti, 2006, MNRAS, 369, 1577
http://adsabs.harvard.edu/abs/2006MNRAS.369.1577C
.. [cassano2012]
Cassano et al. 2012, A&A, 548, A100
http://adsabs.harvard.edu/abs/2012A%26A...548A.100C
.. [donnert2013]
Donnert 2013, AN, 334, 615
http://adsabs.harvard.edu/abs/2013AN....334..515D
.. [donnert2014]
Donnert & Brunetti 2014, MNRAS, 443, 3564
http://adsabs.harvard.edu/abs/2014MNRAS.443.3564D
.. [sarazin1999]
Sarazin 1999, ApJ, 520, 529
http://adsabs.harvard.edu/abs/1999ApJ...520..529S
"""
import logging
import numpy as np
from . import helper
from .solver import FokkerPlanckSolver
from .emission import SynchrotronEmission
from ...share import CONFIGS, COSMO
from ...utils.units import (Units as AU,
UnitConversions as AUC)
logger = logging.getLogger(__name__)
class RadioHalo:
"""
Simulate the extended radio halo emission from galaxy cluster
experiencing on-going/recent merger.
Description
-----------
1. Calculate the merger crossing time (t_cross; ~1 Gyr);
2. Calculate the diffusion coefficient (Dpp) from the systematic
acceleration timescale (tau_acc; ~0.1 Gyr). The acceleration
diffusion is assumed to have an action time ~ t_cross (i.e.,
only during merger crossing), and then been disabled (i.e.,
only radiation and ionization losses later);
3. Assume the electrons are constantly injected and has a power-law
energy spectrum, determine the injection rate by further assuming
that the total injected electrons has energy of a fraction (eta_e)
of the ICM total thermal energy;
4. Set the initial electron density/spectrum be the total injected
electrons during t_merger time;
5. Calculate the magnetic field from the cluster total mass (which
is assumed to be growth linearly from M_main+M_sub to M_obs);
6. Calculate the energy losses for the coefficients of Fokker-Planck
equation;
7. Solve the Fokker-Planck equation to derive the relativistic
electron spectrum at t_obs (i.e., z_obs);
8. Calculate the synchrotron emissivity from the derived electron
spectrum.
Parameters
----------
M_obs : float
Cluster virial mass at the current observation (simulation end) time.
Unit: [Msun]
z_obs : float
Redshift of the current observation (simulation end) time.
M_main, M_sub : float
The main and sub cluster masses before the (major) merger.
Unit: [Msun]
z_merger : float
The redshift when the (major) merger begins.
Attributes
----------
fpsolver : `~FokkerPlanckSolver`
The solver instance to calculate the electron spectrum evolution.
radius : float
The halo radius (scales with the virial radius)
Unit: [kpc]
"""
def __init__(self, M_obs, z_obs, M_main, M_sub, z_merger,
configs=CONFIGS):
self.M_obs = M_obs
self.z_obs = z_obs
self.M_main = M_main
self.M_sub = M_sub
self.z_merger = z_merger
self.configs = configs
self._set_configs()
self._set_solver()
def _set_configs(self):
comp = "extragalactic/halos"
self.eta_turb = self.configs.getn(comp+"/eta_turb")
self.eta_e = self.configs.getn(comp+"/eta_e")
self.gamma_min = self.configs.getn(comp+"/gamma_min")
self.gamma_max = self.configs.getn(comp+"/gamma_max")
self.gamma_np = self.configs.getn(comp+"/gamma_np")
self.buffer_np = self.configs.getn(comp+"/buffer_np")
self.time_step = self.configs.getn(comp+"/time_step")
self.injection_index = self.configs.getn(comp+"/injection_index")
def _set_solver(self):
self.fpsolver = FokkerPlanckSolver(
xmin=self.gamma_min, xmax=self.gamma_max,
x_np=self.gamma_np,
tstep=self.time_step,
f_advection=self.fp_advection,
f_diffusion=self.fp_diffusion,
f_injection=self.fp_injection,
buffer_np=self.buffer_np,
)
@property
def gamma(self):
"""
The logarithmic grid adopted for solving the equation.
"""
return self.fpsolver.x
@property
def age_obs(self):
return COSMO.age(self.z_obs)
@property
def age_merger(self):
return COSMO.age(self.z_merger)
@property
def time_crossing(self):
"""
The time duration of the sub-cluster crossing the main cluster,
which is also used to approximate the merging time, during which
the turbulence acceleration is regarded as effective.
Unit: [Gyr]
"""
return helper.time_crossing(self.M_main, self.M_sub,
z=self.z_merger)
@property
def radius(self):
"""
The halo radius derived from the virial radius by a scaling
relation.
Unit: [kpc]
"""
mass = self.M_main + self.M_sub # [Msun]
r_halo = helper.radius_halo(mass, self.z_merger) # [kpc]
return r_halo
@property
def angular_radius(self):
"""
The angular radius of the radio halo.
Unit: [arcsec]
"""
DA = COSMO.DA(self.z_obs) * 1e3 # [Mpc] -> [kpc]
theta = self.radius / DA # [rad]
return theta * AUC.rad2arcsec
@property
def volume(self):
"""
The halo volume, calculated from the above radius.
Unit: [kpc^3]
"""
return (4*np.pi/3) * self.radius**3
@property
def magnetic_field(self):
"""
The magnetic field strength at the simulated observation
time (i.e., cluster mass of ``self.M_obs``), will be used
to calculate the synchrotron emissions.
Unit: [uG]
"""
return helper.magnetic_field(self.M_obs)
@property
def kT_merger(self):
"""
The cluster ICM mean temperature at z_merger when the merger
begins.
Unit: [keV]
"""
mass = self.M_main + self.M_sub
kT = helper.mass_to_kT(mass, z=self.z_merger)
return kT
@property
def injection_rate(self):
"""
The constant electron injection rate assumed.
Unit: [cm^-3 Gyr^-1]
The injection rate is parametrized by assuming that the total
energy injected in the relativistic electrons during the cluster
life (e.g., ``age_obs`` here) is a fraction (``self.eta_e``)
of the total thermal energy of the cluster.
The electrons are assumed to be injected throughout the cluster
ICM/volume, i.e., do not restricted inside the halo volume.
Qe(γ) = Ke * γ^(-s),
int[ Qe(γ) γ me c^2 ]dγ * t_cluster = eta_e * e_th
=>
Ke = [(s-2) * eta_e * e_th * γ_min^(s-2) / (me * c^2 * t_cluster)]
References
----------
Ref.[cassano2005],Eqs.(31,32,33)
"""
s = self.injection_index
e_thermal = helper.density_energy_thermal(self.M_obs, self.z_obs)
term1 = (s-2) * self.eta_e * e_thermal # [erg cm^-3]
term2 = self.gamma_min**(s-2)
term3 = AU.mec2 * self.age_obs # [erg Gyr]
Ke = term1 * term2 / term3 # [cm^-3 Gyr^-1]
return Ke
def calc_electron_spectrum(self, zbegin=None, zend=None, n0_e=None):
"""
Calculate the relativistic electron spectrum by solving the
Fokker-Planck equation.
Parameters
----------
zbegin : float, optional
The redshift from where to solve the Fokker-Planck equation.
Default: ``self.z_merger``.
zend : float, optional
The redshift where to stop solving the Fokker-Planck equation.
Default: ``self.z_obs``.
n0_e : 1D `~numpy.ndarray`, optional
The initial electron number distribution.
Unit: [cm^-3].
Default: accumulated constantly injected electrons until zbegin.
Returns
-------
electron_spec : float 1D `~numpy.ndarray`
The solved electron spectrum at ``zend``.
Unit: [cm^-3]
"""
if zbegin is None:
tstart = COSMO.age(self.z_merger)
else:
tstart = COSMO.age(zbegin)
if zend is None:
tstop = COSMO.age(self.z_obs)
else:
tstop = COSMO.age(zend)
if n0_e is None:
# Accumulated constantly injected electrons until ``tstart``.
n_inj = self.fp_injection(self.gamma)
n0_e = n_inj * tstart
self.electron_spec = self.fpsolver.solve(u0=n0_e, tstart=tstart,
tstop=tstop)
return self.electron_spec
def calc_emissivity(self, frequencies, n_e=None, gamma=None):
"""
Calculate the synchrotron emissivity for the derived electron
spectrum.
Parameters
----------
frequencies : float, or 1D `~numpy.ndarray`
The frequencies where to calculate the synchrotron emissivity.
Unit: [MHz]
n_e : 1D `~numpy.ndarray`, optional
The electron spectrum (w.r.t. Lorentz factors γ).
If not provided, then used the cached ``self.electron_spec``
solved above.
Unit: [cm^-3]
gamma : 1D `~numpy.ndarray`, optional
The Lorentz factors γ of the electron spectrum.
If not provided, then used ``self.gamma``.
Returns
-------
emissivity : float, or 1D `~numpy.ndarray`
The calculated synchrotron emissivity at each specified
frequency.
Unit: [erg/s/cm^3/Hz]
"""
if n_e is None:
n_e = self.electron_spec
if gamma is None:
gamma = self.gamma
syncem = SynchrotronEmission(gamma=gamma, n_e=n_e,
B=self.magnetic_field)
emissivity = syncem.emissivity(frequencies)
return emissivity
def calc_power(self, emissivity):
"""
Calculate the halo synchrotron power (i.e., power *emitted* per
unit frequency) from emissivity.
Parameters
----------
emissivity : float, or 1D `~numpy.ndarray`
The synchrotron emissivity at multiple frequencies.
Unit: [erg/s/cm^3/Hz]
Returns
-------
power : float, or 1D `~numpy.ndarray`
The calculated synchrotron power w.r.t. each input emissivity.
Unit: [W/Hz]
"""
return helper.calc_power(emissivity, volume=self.volume)
def calc_flux(self, emissivity):
"""
Calculate the synchrotron flux density (i.e., power *observed*
per unit frequency) from emissivity.
Parameters
----------
emissivity : float, or 1D `~numpy.ndarray`
The synchrotron emissivity at multiple frequencies.
Unit: [erg/s/cm^3/Hz]
Returns
-------
flux : float, or 1D `~numpy.ndarray`
The calculated synchrotron flux w.r.t. each input emissivity.
Unit: [Jy] = 1e-23 [erg/s/cm^2/Hz] = 1e-26 [W/m^2/Hz]
"""
power = self.calc_power(emissivity) # [W/Hz]
return helper.calc_flux(power, z=self.z_obs)
def calc_brightness_mean(self, emissivity, frequency, pixelsize=None):
"""
Calculate the mean surface brightness (power observed per unit
frequency and per unit solid angle) expressed in *brightness
temperature* at the specified frequencies from emissivity.
Parameters
----------
emissivity : float, or 1D `~numpy.ndarray`
The synchrotron emissivity at multiple frequencies.
Unit: [erg/s/cm^3/Hz]
frequency : float, or 1D `~numpy.ndarray`
The frequencies where the synchrotron emissivity is calculated.
Unit: [MHz]
pixelsize : float, optional
The pixel size of the output simulated sky image.
Unit: [arcsec]
Returns
-------
Tb : float, or 1D `~numpy.ndarray`
The mean surface brightness at each frequency.
Unit: [K] <-> [Jy/pixel]
"""
omega = np.pi * self.angular_radius**2 # [arcsec^2]
flux = self.calc_flux(emissivity)
return helper.calc_brightness_mean(flux, frequency=frequency,
omega=omega, pixelsize=pixelsize)
def fp_injection(self, gamma, t=None):
"""
Electron injection (rate) term for the Fokker-Planck equation.
NOTE
----
The injected electrons are assumed to have a power-law spectrum
and a constant injection rate.
Qe(γ) = Ke * γ^(-s),
Ke: constant injection rate
Parameters
----------
gamma : float, or float 1D `~numpy.ndarray`
Lorentz factors of electrons
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, or float 1D `~numpy.ndarray`
Current electron injection rate at specified energies (gamma).
Unit: [cm^-3 Gyr^-1]
References
----------
Ref.[cassano2005],Eqs.(31,32,33)
"""
Ke = self.injection_rate # [cm^-3 Gyr^-1]
Qe = Ke * gamma**(-self.injection_index)
return Qe
def fp_diffusion(self, gamma, t):
"""
Diffusion term/coefficient for the Fokker-Planck equation.
NOTE
----
The diffusion coefficients cannot be zero or negative, which
may cause unstable or wrong results. So constrain ``chi_acc``
be a small positive number (e.g., 0.01 [Gyr^-1]).
Parameters
----------
gamma : float, or float 1D `~numpy.ndarray`
The Lorentz factors of electrons
t : float
Current (cosmic) time when solving the equation
Unit: [Gyr]
Returns
-------
diffusion : float, or float 1D `~numpy.ndarray`
Diffusion coefficients
Unit: [Gyr^-1]
References
----------
Ref.[donnert2013],Eq.(15)
"""
chi_acc = self._chi_acceleration(t) # [Gyr^-1]
diffusion = gamma**2 * chi_acc / 4
return diffusion
def fp_advection(self, gamma, 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 [donnert2014], which includes all relevant energy loss
functions and the energy gain function due to turbulence.
Returns
-------
advection : float, or float 1D `~numpy.ndarray`
Advection coefficients, describing the energy loss/gain rates.
Unit: [Gyr^-1]
"""
advection = (abs(self._loss_ion(gamma, t)) +
abs(self._loss_rad(gamma, t)) -
(self.fp_diffusion(gamma, t) * 2 / gamma))
return advection
def _mass(self, t):
"""
Calculate the main cluster mass at the given (cosmic) time.
NOTE
----
We assume that the main cluster grows (i.e., gains mass) linearly
in time from (M_main, z_merge) to (M_obs, z_obs).
Parameters
----------
t : float
The (cosmic) time/age.
Unit: [Gyr]
Returns
-------
mass : float
The mass of the main cluster.
Unit: [Msun]
"""
t_merger = self.age_merger
rate = (self.M_obs - self.M_main) / (self.age_obs - t_merger)
mass = rate * (t - t_merger) + self.M_main
return mass
def _chi_acceleration(self, t=None):
"""
Calculate the electron acceleration coefficient due to turbulent
waves at the given (cosmic) time.
Considering that the turbulence acceleration lies a 2nd-order
Fermi process, it has only a effective acceleration time of
several 1e8 years. Therefore, the turbulence is assumed to
only accelerate the electrons during the merging period, i.e.,
this coefficient "chi" is zero after "t_merger + time_cross".
NOTE
----
A zero diffusion coefficient may lead to unstable/wrong results,
so constrain this acceleration coefficient to be a small positive
but non-zero number.
Parameters
----------
t : float, optional
The (cosmic) time/age.
If not given, then assumed to be within the merging process.
Unit: [Gyr]
Returns
-------
chi : float
The electron acceleration coefficient.
Unit: [Gyr^-1]
References
----------
Ref.[cassano2005],Eq.(40,B12)
"""
# The minimum acceleration coefficient to avoid unstable results
chi_min = 0.01 # [Gyr^-1]
if (t is None) or (t < self.age_merger + self.time_crossing):
mass = self.M_main + self.M_sub
term1 = (mass / 2e15) ** 1.5
term2 = (self.kT_merger / 7) ** (-0.5)
term3 = 500 / self.radius
chi = 6.3 * self.eta_turb * term1 * term2 * term3
else:
chi = chi_min
return chi
def _loss_ion(self, gamma, t):
"""
Energy loss through ionization and Coulomb collisions.
Parameters
----------
gamma : float, or float 1D `~numpy.ndarray`
The Lorentz factors of electrons
t : float
The cosmic time/age
Unit: [Gyr]
Returns
-------
loss : float, or float 1D `~numpy.ndarray`
The energy loss rates
Unit: [Gyr^-1]
References
----------
Ref.[sarazin1999],Eq.(9)
"""
z = COSMO.redshift(t)
mass = self._mass(t)
n_th = helper.density_number_thermal(mass, z) # [cm^-3]
loss = -3.79e4 * n_th * (1 + np.log(gamma/n_th) / 75)
return loss
def _loss_rad(self, gamma, t):
"""
Energy loss via synchrotron emission and inverse Compton
scattering off the CMB photons.
References
----------
Ref.[sarazin1999],Eq.(6,7)
"""
z = COSMO.redshift(t)
mass = self._mass(t) # [Msun]
B = helper.magnetic_field(mass) # [uG]
loss = -4.32e-4 * gamma**2 * ((B/3.25)**2 + (1+z)**4)
return loss