<jupyter_start><jupyter_text>Example to calculate photon-ALP oscillations from NGC 1275 in a cavity field This notebook demonstrates how to calculate the photon-ALP transition probability for NGC 1275, the central AGN of the Perseus cluster. The assumed B-field environments are the same as in https://arxiv.org/abs/1908.03084, and include the cluster field and the magnetic field of the Milky Way. Electron density model parameters taken from https://iopscience.iop.org/article/10.1086/374923/fulltext/57155.text.html.Note that in order to obtain the field used in 1908.03084, you have to set $\theta$ to 225° instead of 45° (adding 180°). This effectively inverts the sign of all field components. While this change of $\theta$ and the sign is irrelevant for the conversion probability, it changes the field along the line of sight and inverts the sign of the RM, as well.<jupyter_code>from gammaALPs.core import Source, ALP, ModuleList
from gammaALPs.base import environs, transfer
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline<jupyter_output><empty_output><jupyter_text>Set the ALP Initialize an ALP object, that stores the ALP mass $m$ (in neV) and the coupling $g$ (in $10^{-11}\mathrm{GeV}^{-1}$).<jupyter_code>m, g = 2, 100
alp = ALP(m,g)<jupyter_output><empty_output><jupyter_text>Set the source Set the source properties (redshift and sky coordinates) in the ```Source``` containier<jupyter_code>ngc1275 = Source(z = 0.017559, ra = '03h19m48.1s', dec = '+41d30m42s')<jupyter_output><empty_output><jupyter_text>Init the module list Initialize the list of transfer modules that will store the different magnetic field environments.  Energies are supplied in GeV as ```numpy.ndarray```<jupyter_code>logEMeV = np.linspace(2, 5, num=150)
EGeV = np.power(10, logEMeV - 3)
EMeV = np.power(10, logEMeV)<jupyter_output><empty_output><jupyter_text>Now initialize the initial photon polarization. Since we are dealing with a gamma-ray source, no ALPs are initially present in the beam (third diagonal element is zero). The polarization density matrix is normalized such that its trace is equal to one, $\mathrm{Tr}(\rho_\mathrm{in}) = 1$.<jupyter_code>pin = np.diag((1.,1.,0.)) * 0.5<jupyter_output><empty_output><jupyter_text>Add modules: Now we add propagation modules for the cluster, the EBL, and the Galactic magnetic field; run the module. Repeat for only the ICM field.<jupyter_code>ml = ModuleList(alp, ngc1275, pin = pin, EGeV = EGeV)
ml.add_propagation("ICMStructured", 
                   0, # position of module counted from the source. 
                   B0 = 8.3,  # rms of B field
                   R = 93,
                   theta = 225,
                   theta_rad=False,
                   pa=147,
                   pa_rad=False,
                   n0 = 3.9e-2,  # normalization of electron density
                   n2 = 4.05e-3, # second normalization of electron density, see Churazov et al. 2003, Eq. 4
                   r_abell = 500., # extension of the cluster
                   r_core = 80.,   # electron density parameter, see Churazov et al. 2003, Eq. 4
                   r_core2 = 280., # electron density parameter, see Churazov et al. 2003, Eq. 4
                   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, # scaling of B-field with electron denstiy
                   cell_num=1000
                   )
ml.add_propagation("EBL",1, eblmodel = 'dominguez') # EBL attenuation comes second, after beam has left cluster
ml.add_propagation("GMF",2, model = 'pshirkov', model_sym = 'BSS') # finally, the beam enters the Milky Way Field

px, py, pa = ml.run()


# again with only structured field

ml2 = ModuleList(alp, ngc1275, pin = pin, EGeV = EGeV)
ml2.add_propagation("ICMStructured", 
                   0, # position of module counted from the source. 
                   B0 = 8.3,  # rms of B field
                   R = 93,
                   theta = 225,
                   theta_rad=False,
                   pa=147,
                   pa_rad=False,
                   n0 = 3.9e-2,  # normalization of electron density
                   n2 = 4.05e-3, # second normalization of electron density, see Churazov et al. 2003, Eq. 4
                   r_abell = 500., # extension of the cluster
                   r_core = 80.,   # electron density parameter, see Churazov et al. 2003, Eq. 4
                   r_core2 = 280., # electron density parameter, see Churazov et al. 2003, Eq. 4
                   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 # scaling of B-field with electron denstiy
                   )
px2, py2, pa2 = ml2.run()<jupyter_output>[0;36menvirons.py:[0;35m 569[0;0m --- [1;36mINFO[1;0m: Using inputted chi
[0;36menvirons.py:[0;35m1175[0;0m --- [1;36mINFO[1;0m: Using inputted chi
[0;36m   core.py:[0;35m 652[0;0m --- [1;36mINFO[1;0m: Running Module 0: <class 'github.galps.gammaALPs.base.environs.MixICMStructured'>
[0;36m   core.py:[0;35m 652[0;0m --- [1;36mINFO[1;0m: Running Module 2: <class 'github.galps.gammaALPs.base.environs.MixGMF'>
[0;36menvirons.py:[0;35m 569[0;0m --- [1;36m[1;36mINFO[1;0m[1;0m: Using inputted chi
[0;36m   core.py:[0;35m 652[0;0m --- [1;36m[1;36mINFO[1;0m[1;0m: Running Module 0: <class 'github.galps.gammaALPs.base.environs.MixICMStructured'><jupyter_text>Plot the output<jupyter_code>plt.plot(EMeV, px[0] + py[0], label='ICM + EBL + GMF')
plt.plot(EMeV, px2[0] + py2[0], label='ICM')
plt.xscale('log')
plt.yscale('linear')
plt.ylim(0.4, 1.05)
#plt.xlim(1e0, 1e2)
plt.xlabel('Energy (MeV)')
plt.ylabel('Photon survival probability')
plt.legend(loc = 0, fontsize = 'medium')

plt.annotate(r'$m_a = {0:.1f}\,\mathrm{{neV}}, g_{{a\gamma}} = {1:.1f} \times 10^{{-11}}\,\mathrm{{GeV}}^{{-1}}$'.format(ml.alp.m,ml.alp.g),
             xy = (0.95,0.1), size = 'x-large', xycoords = 'axes fraction', ha = 'right')

plt.subplots_adjust(left = 0.2)
#plt.savefig("pgg.png", dpi = 150)<jupyter_output>[0;36mfont_manager.py:[0;35m1355[0;0m --- [1;31m[1;31mWARNING[1;0m[1;0m: findfont: Font family ['serif'] not found. Falling back to DejaVu Sans.<jupyter_text>Plot the B fieldsPlot r vs B field components of structured field, also plot electron density<jupyter_code>r = ml.modules[0].r
b_r = ml.modules[0]._Bfield_model.b_r
b_theta = ml.modules[0]._Bfield_model.b_theta
b_phi = ml.modules[0]._Bfield_model.b_phi
plt.plot(r, b_r, label='$B_r$')
plt.plot(r, b_theta, label=r'$B_\theta$')
plt.plot(r, b_phi, label=r'$B_\phi$')
plt.xlabel(r'$r [\SI{}{\kilo\parsec}]$')
plt.ylabel(r'$B [\SI{}{\micro\gauss}]$')
plt.legend()
nel = ml.modules[0]._nelicm(r)   # in cm^-3
plt.plot(r, nel)
plt.yscale('log')
plt.xscale('log')
plt.ylabel('$n_\mathrm{el}$ (cm$^{-3}$)')
plt.xlabel('$r$ (kpc)')<jupyter_output><empty_output><jupyter_text>Calculate rotation measureFaraday RM can be calculated via a method of the structured field model, argument is electron density evaluated at r = structured_field.r in cm^-3<jupyter_code>ml.modules[0]._Bfield_model.rotation_measure(nel)   # in rad * m^-2<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13214722259948392, 'fraction_numerical': 0.05717778079587125, 'mean_word_length': 3.261574074074074, 'pattern_counts': {'":': 0, '<': 38, '<?xml version=': 0, '>': 38, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '257702', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 2911, 'n_tokens_neox': 2513, 'n_words': 726, 'path': 'docs/tutorials/mixing_ICM_structured_field.ipynb', 'repo_name': 'specktakel/gammaALPs', 'size': '7363'}