Calculating the coherence length and rotation measure for Gaussian turbulent field
This tutorial demonstrates how to the coherence length and rotation measure of a magnetic field with Gaussian turbulence on the example of the Perseus cluster with the central radio galaxy NGC 1275. The assumed B-field environments are the same as in Ajello et al. (2016).
If you haven’t installed gammaALPs
already, you can do so using pip
. Just uncomment the line below:
[1]:
#!pip install gammaALPs
First some imports:
[2]:
from gammaALPs.core import Source, ALP, ModuleList
import numpy as np
import matplotlib.pyplot as plt
import time
from scipy.integrate import simps
from scipy.stats import norm
Next we set the source properties (redshift and sky coordinates), the ALP, and the energy range. We only need those to init the ModuleList
later – we won’t actually run the photon-ALP conversion probability calculation here.
[3]:
ngc1275 = Source(z=0.017559, ra='03h19m48.1s', dec='+41d30m42s')
alp = ALP(m=1., g=1.)
[4]:
ml = ModuleList(alp, ngc1275)
We add the propagation within the intra-cluster Gaussian turbulent field with the same values as in Ajello et al. (2016).
[5]:
ml.add_propagation("ICMGaussTurb",
0, # position of module counted from the source.
nsim=10, # number of random B-field realizations
B0=10., # rms of B field in muG
n0=3.9e-2, # normalization of electron density in cm-3
n2=4.05e-3, # second normalization of electron density, see Churazov et al. 2003, Eq. 4 on cm-3
r_abell=500., # extension of the cluster in kpc
r_core=80., # electron density parameter, see Churazov et al. 2003, Eq. 4 in kpc
r_core2=280., # electron density parameter, see Churazov et al. 2003, Eq. 4 in kpc
beta=1.2, # electron density parameter, see Churazov et al. 2003, Eq. 4
beta2=0.58, # electron density parameter, see Churazov et al. 2003, Eq. 4
eta=0.5, # scaling of B-field with electron denstiy
kL=0.18, # maximum turbulence scale in kpc^-1, taken from A2199 cool-core cluster, see Vacca et al. 2012
kH=9., # minimum turbulence scale, taken from A2199 cool-core cluster, see Vacca et al. 2012
q=-2.80, # turbulence spectral index, taken from A2199 cool-core cluster, see Vacca et al. 2012
seed=0 # random seed for reproducability, set to None for random seed.
)
environs.py: 431 --- INFO: Using inputted chi
We peek at the transversal magnetic field and the electron density:
[6]:
plt.plot(ml.modules["ICMGaussTurb"].r,
ml.modules["ICMGaussTurb"].B * np.sin(ml.modules["ICMGaussTurb"].psi),
lw=1)
plt.plot(ml.modules["ICMGaussTurb"].r,
ml.modules["ICMGaussTurb"].B * np.cos(ml.modules["ICMGaussTurb"].psi),
lw=1, ls = '--')
plt.ylabel('$B$ field ($\mu$G)')
plt.xlabel('$r$ (kpc)')
[6]:
Text(0.5, 0, '$r$ (kpc)')
[7]:
plt.loglog(ml.modules["ICMGaussTurb"].r,ml.modules["ICMGaussTurb"].nel * 1e-3)
plt.ylabel('$n_\mathrm{el}$ (cm$^{-3}$)')
plt.xlabel('$r$ (kpc)')
[7]:
Text(0.5, 0, '$r$ (kpc)')
Spatial correlation and coherence legnth
The gammaALPs.bfields.Bgaussian
class has methods to calculate the spatial correlation of the magnetic field and the rotation measure. We can access these methods through the the magnetic field model mehtod, which we can access through ml.modules['ICMGaussTurb'].Bfield_model
.
The spatial correlation \(C(x_3) = \langle B_\perp(\vec{x}) B_\perp(\vec{x} + x_3 \vec{e}_3)\rangle\) of the transversal magnetic field along the line of sight \(z\) is computed like this:
[8]:
z = np.linspace(0.,50.,1000) # distance in kpc from cluster center
c = ml.modules["ICMGaussTurb"].Bfield_model.spatial_correlation(z)
plt.plot(z, c / c[0])
plt.xlabel("$z$ (kpc)")
plt.ylabel("$C(z) / C(0)$")
plt.grid(True)
This is turn can be used to calculate the coherence length of the field,
[9]:
z = np.linspace(0.,1e3,1000) # distance in kpc from cluster center
c = ml.modules["ICMGaussTurb"].Bfield_model.spatial_correlation(z)
Lambda_c = simps(c, z) / c[0]
print ("Coherence length of the field is Lambda_C = {0:.3e} kpc".format(Lambda_c))
Coherence length of the field is Lambda_C = 1.492e+00 kpc
Calculate the rotation measure of the field
The rotation measure describes the rotation of the polarization vector due to the propagation of an electromagnetic wave through a medium with magnetic field and an electron plasma.
The change in polarization angle \(\Delta\phi\) is proportianol to the wavewavelength and the rotation measure \(\mathrm{RM}\) which is given by the integral over the magnetic field parallel to the propagation direction multiplied with the electron density:
We can calcluate the \(\mathrm{RM}\) through the ml.modules['ICMGaussTurb'].Bfield_model.rotation_measure
function which computes \(B_{||}\) for a requested number nsim
of random realizations. In addition to the line of sight we’re interested in, we also have to provide the electron density along \(z\) and the scaling of the magnetic field with the electron density, i.e. the factor \((n_\mathrm{el}(z) / n_\mathrm{el}(z=0)^\eta\). The latter can be calculated through the
ml.modules["ICMGaussTurb"].nel_model.Bscale
function.
We calculate the electron density and the scaling:
[10]:
n_el = ml.modules["ICMGaussTurb"].nel
Bscale = ml.modules["ICMGaussTurb"].nel_model.Bscale(ml.modules["ICMGaussTurb"].r)
And provide this to the rotation_measure
function. The loop over the realizations is rather slow, so we calculate the \(\mathrm{RM}\) for only 500 random realizations:
[11]:
ml.modules["ICMGaussTurb"].Bfield_model.seed = 0
t1 = time.time()
nsim=500
rm = ml.modules["ICMGaussTurb"].Bfield_model.rotation_measure(ml.modules["ICMGaussTurb"].r,
n_el=n_el,
Bscale=Bscale,
nsim=nsim)
t2 = time.time()
print("Calculating RM for {0:d} realizations took {1:.2f} seconds".format(nsim, t2 - t1))
Calculating RM for 500 realizations took 17.96 seconds
Finally, we plot the histogram of the \(\mathrm{RM}\) values. Since we’re assuming a purely turbulent field, we expect a distribution that peaks close to zero and we can compare the spread to the \(\mathrm{RM}\) values reported by Taylor et al. (2006) who found RM values between 6500 and 7500 \(\mathrm{rad}\,\mathrm{m}^{-2}\).
[12]:
n, bins, _ = plt.hist(np.sort((rm)), bins=30, density=True, label="Simulated RM")
mean = np.mean(rm) # mean
var = np.var(rm) # variance
print("RM mean +/- sqrt(var) in rad m^-2: {0:.2f} +/- {1:.2f}".format(mean, np.sqrt(var)))
plt.plot(bins, norm.pdf(bins, loc=mean, scale=np.sqrt(var)),
lw=2,
label="Gaussian Fit\n$\mu = {0:.2f}$\n$\sigma={1:.2f}$".format(mean, np.sqrt(var)))
plt.legend()
plt.gca().tick_params(labelleft=False, left=False, right=False, top=False)
plt.xlabel("Rotation Measure (rad m${}^{-2}$)")
plt.ylabel("Density")
RM mean +/- sqrt(var) in rad m^-2: 49.62 +/- 2707.52
[12]:
Text(0, 0.5, 'Density')
With our chosen \(B\) field, \(\sigma(\mathrm{RM}) \sim 2700 \,\mathrm{m}^{-2}\).
[ ]: