1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
|
# 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 ...utils import cosmo
from ...utils.units import (Units as AU,
UnitConversions as AUC,
Constants as AC)
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;
4. Determine the initial electron density
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.
Parameters
----------
M0 : float
Cluster virial mass at redshift z0
Unit: [Msun]
z0 : float
Redshift from where to simulate former merging history.
"""
def __init__(self, M0, z0):
self.M0 = M0
self.z0 = z0
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.zmax``.
zend : float, optional
The redshift where to stop solving the Fokker-Planck equation.
Default: ``self.z0``.
n0_e : 1D `~numpy.ndarray`, optional
The initial electron number distribution.
Unit: [cm^-3].
Default: accumulated constant-injected electrons until zbegin.
Returns
-------
gamma : `~numpy.ndarray`
The Lorentz factor grid adopted for solving the equation.
n_e : `~numpy.ndarray`
The solved electron spectrum at ``zend``.
Unit: [cm^-3]
"""
if zbegin is None:
tstart = cosmo.age(self.zmax)
else:
tstart = cosmo.age(zbegin)
if zend is None:
tstop = cosmo.age(self.z0)
else:
tstop = cosmo.age(zend)
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,
)
gamma = fpsolver.x
if n0_e is None:
# Accumulated constant-injected electrons until ``tstart``.
n_inj = np.array([self.fp_injection(gm) for gm in gamma])
n0_e = n_inj * tstart
n_e = fpsolver.solve(u0=n0_e, tstart=tstart, tstop=tstop)
return (gamma, n_e)
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 = cosmo.age(z_begin) # [Gyr]
t_end = t_begin + time
if t_end >= cosmo.age(0):
z_end = 0.0
else:
z_end = cosmo.redshift(t_end)
return z_end
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]
age = cosmo.age(self.z0)
term1 = (self.injection_index-2) * self.eta_e
term2 = e_th / (self.pmin * self.mec * AC.c) # [cm^-3]
term3 = 1.0 / (age * 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 = 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, 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
Advection coefficient, describing the energy loss/gain rate.
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 _loss_ion(self, gamma, t):
"""
Energy loss through ionization and Coulomb collisions.
Parameters
----------
gamma : float
The Lorentz factor of electrons
t : float
The cosmic time/age
Unit: [Gyr]
Returns
-------
loss : float
The energy loss rate
Unit: [Gyr^-1]
References
----------
Ref.[sarazin1999],Eq.(9)
"""
z = cosmo.redshift(t)
mass = self._mass(t)
n_th = helper.density_number_thermal(mass, z)
coef = -1.20e-12 * AUC.Gyr2s # [Gyr^-1]
loss = coef * 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)
coef = -1.37e-20 * AUC.Gyr2s # [Gyr^-1]
loss = (coef * gamma**2 *
((self.magnetic_field/3.25)**2 + (1+z)**4))
return loss
|
