chap4_1.qmd

# Gamma-ray Astronomy

As we saw, the origin of CR remains a mystery due to the deflection of CRs in their travel. Only CR astronomy at $10^{19}$ eV will be possible if the composition of these UHECR are protons, but even in this case the GZK limits the horizon. 

Gamma-ray astronomy provides an fundamental tool to observe the Universe using neutral deep penetrating gamma-ray particles. Originally the distinction between gamma-rays and X-rays was due to the nuclear production mechanism. X-rays are the product of transition of electrons in the atomic shell, while gamma-rays are produced in the atomic nucleus. This distinction also leads to a classification in energy. X-rays are typically below 100 keV. Photons with E > 100 keV are called gamma-rays. 

![Source: Wikipedia](images/EM_spectrum.svg)

## Spectral Energy Distribution

In gamma-ray astronomy is usual to study an object emission by its spectral energy distribution (SED).  The SED on an object is its energy emitted plotted against some measure of the photon - frequency or wavelength. The reason astronomers do this is to see how much energy is produced by the object as a function of frequency or wavelength. The SED is typically characterized by $\nu F_\nu$ and it is measured in units of ergs $\mathrm{cm}^{-2}\; \mathrm{s}\;^{-1}$, ie. it indicates the rate of energy emitted per surface. The function $F_\nu$ is the flux density which indicates the rate of energy emitted per surface and also per frequency, so it is expressed as ergs $\mathrm{cm^{-2}\; s^{-1}\; Hz^{-1}}$. All of this doesn't need to be confused with the spectrum! which is what we are being using to characterize CRs. In CR physics spectrum in the rate of particles per unit energy and surface and it is measured on $\mathrm{GeV^{-1}\; cm^{-2}\; s^{-1}}$. The spectrum multiply by energy, gives you the rate of particles per unit surface. Multiplying again by energy give use the energy rate per unit surface. Therefore we have the relation:

$$\nu F_\nu = E^2 \frac{\mathrm{d}N}{\mathrm{d}E}$$

## Gamma-ray production mechanism
Sources of CR can also produce gamma-rays by different mechanisms roughly divided in two main categories:

### Leptonic models: 
In this models only leptons (mostly electrons) will produce gamma-ray emission. The mechanismare mostly:
   * **Synchrotron radiation.** 
   * **Inverse Compton Scattering.**

### Hadronic model:
In hadronic models, gamma-rays are produced as a result of acceleration to high energies of protons or other hadrons. The mechanism of gamma-ray production in hadronic models are mainly:
   * **$\pi^0$ decay.** 
   * **Proton synchrotron radiation.** 

Another categories for gamma-ray production can some exotic models like dark matter annihilaiton, or matter-antimatter annihilation and nuclear transformation.

## Leptonic models

### Synchrotron radiation

Synchrotron radiation is extremely important for astrophysics as it was realized by Shklovskii in 1957 when studing the non-thermal emission of the Crab remnant. In order to understand better the synchrotron radiation we will have to dig in a bit in electromagnetism. We already saw that a charge particle of charge $q$, for example an electron, moving with velocity $v$ in a magnetic field $\vec{B}$ feels an an external force:

$$\vec{F} = \frac{q}{c}(\vec{v}\times \vec{B}) $$

Now, because of the force on the particle is perpendicular to the motion, the magnetic field cannot do work on the particle, and so its speed does not change, i.e. $|v|$ = constant, but there is an acceleration since the direction will change. On the other hand an accelerated electrical charge radiates electromagnetic waves! which... will slow down the particle... 

The reason of this aparent inconsistency comes from the fact that we treat the electric field lines as the purely Coulombic action-at-a-distance. In reality we have to take into account that as a particle moves, the electric field lines need to *re-arrange* and this re-arrangement cannot happen at a speed faster than the speed of light. 

Here are going to derive the radiation emission for a particle with an acceleration and we are going to do it by messing around with the Coulomb fields. For a more formal argument derived from Maxwell equations you can see the book of Longair.

Let's assume that the particle is at rest at the moment $t=0$. The electric field lines clearly point away from the origin. At that moment the particle accelerates which brings the velocity of the particle to $\Delta v$ in a time $\Delta t$, after that the particle continues with uniform velocity. 

![](images/bremsstrahlung.svg)

After a certain time, the particle will be in the position $t\Delta v$. In a sphere located far way (with a radius larger than $c t$), the electric field lines are still those of the stacionary particle, since they field lines cannot "know" yet that the particle has moved, so the point radially towards the origin at $t=0$. Inside a sphere of radius $c(t - \Delta t$), the electric field lines are already those from the electron that moves at a constant velocity.

The perturbation of the electric field lines needs to propagate radially, this *kink* is nothing more than a radiation! From simple geometry relations we can get that:

$$\frac{E_\theta}{E_{r}} = \frac{\Delta v  t \sin\theta}{c\Delta t} = \frac{r \sin \theta}{c^2}\frac{\Delta v}{\Delta t}$$

where $r = c t$, and $E_r$ is simply the Coulomb field $E_r = q/4\pi\epsilon_0 r^2$ or in [gaussian units](https://en.wikipedia.org/wiki/Gaussian_units) since the kink is really small and the eletric field along the kink remains the same $E_r = q/r^2$ and therefore the transverse component becomes:

$$E_\theta = \frac{q a \sin\theta}{c^2 r}$$

where we used the fact that $\Delta v/\Delta t$ is just the acceleration $a$. An interesting fact appears here, $E_\theta$ depends on $r^{-1}$ and  not $r^{-2}$ so for larger $r$, $E_\theta$ is going to dominate over $E_r$.  Accompanying this transverse electric field there will be a magnetic field, which is a property of an electromagnetic wave. In other words, an electromagnetic pulse is generated by the accelerated charge particle. Since this is an electromagnetic radiation there is an energy flow per unit area, per second and the direction is given by the Poynting vector (with $|E|=|B|$ as in electromagnetic wave):


$$\vec{S}  = \frac{c}{4\pi}(\vec{E}\times\vec{B})$$

Which points in the radial direction. In this case can be reduced to:

$$|S| \equiv \frac{ \mathrm{d}E}{\mathrm{d}t \mathrm{d}A} = \frac{c}{4\pi} E^2_\theta \vec{n}$$

Which is the energy flow per unit area per second. The unit area $\mathrm{d}A$ can be rewritten in terms of the solid angle as $\mathrm{d} A= r^2\mathrm{d}\Omega$, and so the rate of energy loss through the area subtended by the solid angle $\mathrm{d} \Omega$ at distance $r$ is given by:

$$\left(\frac{\mathrm{d}E}{\mathrm{d}t \mathrm{d} \Omega}\right)_{rad}  = \frac{q^2a^2 \sin^2 \theta}{4\pi c^3}$$

Notice that this energy loss rate follows a *dipole* pattern $\mathrm{d}P/\mathrm{d}\Omega \propto \sin^2 \theta$ and that $\theta$ is defined along the acceleration line. If we now integrate over all solid angles we obtain that the emitted power is given by:


$$P \equiv -\left(\frac{\mathrm{d}E}{\mathrm{d}t}\right)_{rad} = \frac{2q^2 a^2 }{3c^3}$$

Which is call the *Larmor's formula*. Note that this is valid for any form of acceleration $a$.


In the relativistic case this formula can be rewritten as:

$$P  = \frac{2q^2}{3c^3} \gamma^4 [ \gamma^2 a_{\parallel}^2 + a^2_{\perp}] $$ 



### Single-electron in uniform magnetic field

In an uniform magnetic field, a high energy charged particle, for example an electron, moves in spiral path at a constant *pitch angle*, $\alpha$.

![](images/spiral.svg){fig-align="center"}

Its velocity along the field lines is constant $v_{\parallel} = v \cos\alpha = \mathrm{const.}$, but its circular component $v_\perp = v\sin\alpha$. Let's first attack the non-relativistic case.

### Non relativistic case: Cyclotron radiation

From Newton's law and the Lorentz force we have that:

$$ m_e a_\perp = m_e \frac{v_\perp^2}{r_{gyr}} = \frac{e v_\perp B}{c}$$ 

where $r_{gyr}$ is the gyroradius which can be written as:

$$r_{gyr} = \frac{v_\perp}{\omega_{gyr}} = \frac{v\sin\alpha}{\omega_{gyr}}$$

where we can define the gyrofrecuency as:

$$\omega_{gyr} \equiv \frac{eB}{m_e c}$$

Therefore according to the Larmor's equation of power emitted we can write:

$$P = \frac{2 e^2}{3c^3} \omega_{gyr}^2 v^2\sin^2\alpha$$

Note the Larmor formula does not tell us frequency spectrum, but if a particle is moving in a circular motion, then from an observer far way, the "aparent" motion will be sinosoidal as illustrated in the figure below. 

![](images/cyclotron.svg)

In particular since the power $P$ depends on  $|v^2|$ for a distance in the $x-y$ plane, power will varie as $\sin^2(\omega_{gyr}t)$ and so the electric field line will change as $E(t)\propto \sin(\omega_{gyr}t)$. The radiated emission will appear as *monochromatic* with an angular frequency given precisely by the particle circular frequency $\omega_{gyr}$. This is the known as the *cyclotron* radiation. In this case the dipole emission pattern is turning around with the particle and an observer in the laboratory frame will see the same frequency. As a consequence also the radiation is independent of the viewing angle. 

### Relativistic case: Synchrotron radiation.{#synchrotron}

Things are slighly different when we consider the relativistic case where electrons move at $\beta \rightarrow 1$. Using again Newton's law with the Lorentz force we have:

$$\gamma m_e \frac{\mathrm{d} v_\perp}{\mathrm{d} t} = \frac{e v_\perp B}{c}$$ 

where we can obtain the perpendicular acceleration as:

$$a_\perp = \frac{e v\sin\alpha}{\gamma m_e c}$$


where now the electrons's relatisvistic angular gyroradius frequency is given by:


$$\omega^{rel}_{gyr} \equiv \frac{eB}{\gamma m_e c}$$

which is exactly related to the classical angular frequency as $\omega_{gyr}/\gamma$. According to the relativistic Larmor's formula above we have that:

$$P = \frac{2 e^2}{3c^3} \gamma^4 \left(\frac{e v\sin\alpha}{\gamma m_e c}\right)^2 $$

We can re-arrange this formula by assuming ultra-relativistic electrons $\beta \rightarrow 1$ and that we have an average number of electrons with different pitch angles and so $\langle \sin^2 \alpha \rangle = 2/3$. We are going also to use the definition of the Thomson cross-section:

$$\sigma_T \equiv \frac{8\pi}{3}r^2_0$$ 

where $r_0$ is the electron classic radius defined (in gaussian units) as:

$$r_0 = \frac{e^2}{m_e c^2} $$

with that we can write the power loss as:

$$P = \frac{4}{3} \sigma_T c U_B \gamma^2$$

where $U_B = B^2/8\pi$ is the magnetic field energy density. The power depends on the electron mass via the Thomson cross-section and the $\gamma$. In total it gives $P\propto m_e^{-4}$ for the electrons. That's why for protons the syncrotron radiation is not as important as it is for electrons. The mean free-path for electrons given their losses sue to synchrotron is (for typical astrophysical values):

$$l_{syn} = \left(\frac{1}{E}\frac{\mathrm{d}E}{\mathrm{d}t}\right)^{-1}  \sim 100 \mathrm{\;pc}$$

That's the reason why electrons are not the main component of CR since the cannot travel long distances. 
If we want now to derive the power emitted per frequency, things get also a bit more complicated. One first thing to take into account is that in the relativistic regime the power radiated will be distorted. From Lorentz transformations we can derive that an angle $\Psi$ in the particle rest frame will be seen in laboratory frame as:

$$\sin \Psi = \frac{\sin\Psi^\prime}{\gamma(1+\beta \cos\Psi^\prime)}$$ 

(we saw this in the Exercises 1). For small angles we have that 

$$\sin \Psi \sim \Psi \sim \frac{1}{\gamma}$$

i.e, while the power is radiated nearly isotropically (dipole) in the particle's instantaneus rest frame, most of it will be *beamed* into a narrow cone of angle $\sim 1/\gamma$ in the laboratory frame. Note that I'm referring to an *instantaneous* rest frame, what does it even mean? Well it means that I in the particle rest frame, where the velocity of the particle is 0, the acceleration on the other hand is not 0, so eventually the particle will move. Ie, there is no particle rest frame, but there is an instantaneous rest frame.

This beaming in the lab reference frame has an impact on how the viewer sees the electric field variation. It is no longer a sinusoidal, instead the emission is pulsed every time the cone sweeps around the line of sight of an observer as illustrated below:


![](images/synchrotron_0.svg){fig-align="center"}


Since each cone beaming has a half amplidude of $1/\gamma$ the radiation from an observer is visible during a $\Delta \theta \sim 2/\gamma$. During this angular distance the particle travels from $A$ to $B$ at an angular frequency given by $\omega_{gyr}^{rel}$ which means that $\Delta t_{AB}$ is given by:

$$\Delta t_{AB} = \frac{\Delta \theta}{\omega_{gyr}^{rel}} = \frac{2}{\gamma \omega_{gyr}^{rel}}$$

Given the relativistic speed of the the electron, the frequency of the radiation is no longer a sinusoidal. When a photon is emitted at the point $B$ the photon emitted in $A$ has already travelled a distance $c\Delta t_{AB}$. Therefore the time interval between photons is given by:


$$\Delta t_{sync} = \frac{c\Delta t_{AB} - v_\perp \Delta t_{AB}}{c} \approx \Delta t_{AB}( 1 - \beta) \approx  \Delta t_{AB} \frac{1 - \beta^2}{1 + \beta} $$

$$ \Delta t_{sync} \approx \frac{ \Delta t_{AB} }{2\gamma^2} = \frac{1}{\gamma^3 \omega_{gyr}^{rel}} = \frac{1}{\gamma^2 \omega_{gyr}} $$

So synchrotron radiation is a very spiky series of widely spaced narrow pulses of time $2\Delta t_{sync}$, and adjacent spikes are separated in frequency $\Delta t_{gyr} = 1/\nu^{rel}_{gyr} = \gamma/\nu_{gyr}$. The emission from an observer is illustrated below:

![](images/synchrotron_1.svg)

The power spectrum (ie power as function of frequency) illustrated in the plot on the right is no long a monochromatic sinusoidal signal, but the Fourier transform of this time series of pulses. Taking of all the Fourier components into account results in the following expression for one electron (the formal derivation can be found for example in *Classical Electrodynamics* by J. D. Jackson):

$$P(\nu)\mathrm{d}\nu = \frac{\sqrt{3}e^3B\sin\alpha}{m_ec^2}\frac{\nu}{\nu_c}\int_{\nu/\nu_c}^\infty K_{5/3}(\psi)\mathrm{d}\psi$$

where $K_{5/2}$ is the Bessel function of order 5/2, and we expressed the power spectrum in terms of the frequency $\nu = \nu_{gyr} =\omega_{gyr}/2\pi$ and where we defined the critical frequency as:

$$\nu_c \equiv \frac{3}{2}\gamma^2\nu_{gyr} \sin \alpha$$

which the frequency at which the power spectrum will peak. In general terms we can describe the shape of the synchrotron power spectrum of a single electron has a logarithmic slope at low frequencies as $P \propto \nu^{1/3}$, a broad peak near the critical frequency $\nu_c$, and a sharp fall off at higher frequencies.

#### Spectrum for several electrons for optically thin sources

If a synchrotron source containing any arbitrary distribution of electron energies is optically thin ($\tau \ll 1$), then its spectrum is the superposition of the spectra from individual electrons and its flux density cannot rise faster than $\nu^{1/3}$ at any frequency $\nu$. The energy distribution of cosmic-ray electrons in most synchrotron sources is roughly a power law:

$$n(E) \mathrm{d}E \propto E^{-\delta}\mathrm{d}E$$

We make the very simple and crude approximation that each electron radiates all of its power at a single critical frequency:

$$\nu_c \simeq \gamma^2\nu_{gyr}$$

We can assume that the luminosity or total power emitted per unit volume per unit frequency, $L(\nu)\mathrm{d}\nu$ is given by the power emitted for each individual electron times the number of electrons:

$$ L (\nu) \mathrm{d}\nu = P n(E)\mathrm{d}E = -\left(\frac{\mathrm{d}E}{\mathrm{d}t}\right) n(E)\mathrm{d}E $$

where 
$$E = \gamma m_e c^2 \simeq \left(\frac{\nu}{\nu_{gyr}}\right)^{1/2}m_e c^2$$

and so

$$\frac{\mathrm{d}E}{\mathrm{d}\nu} \simeq \frac{m_e c^2 \nu^{-1/2}}{2\nu_{gyr}^{1/2}}$$

putting all together gives:

$$L_\nu \propto B^{(\delta+1)/2}\nu^{(1-\delta)/2}$$

That is, the synchrotron spectrum of a power-law energy distribution is itself a power law. This idea is represented in the figure below:

![](images/synchrotron_spectrum.svg){fig-align="center"}


For example in our Galaxy we expect the electron population to have an spectral index of $\delta \sim 2.4$, and so the synchrotron radiation should have an spectral index of $\sim 0.7$. 

#### Power spectrum for several electrons for optically thick sources

The argument of a power law synchrotron spectrum as a power-low with specltral index $\alpha = (1 - \delta)/2$ it's true as long as no absorption of photons by the emitting region happens. At low frequencies however synchrotron suffers from self-absorption, in which a photon interacts with a charge in a magnetic field. This also happens in *optically thick* sources. Therefore at low frequencies its emission is absorbed and re-emitted as blackbody radiation.
The derivation of the synchrotron self-absorption is complicated so we are going to give only an heuristic approach. The problem with the re-emission of low frequencies as a black-body radiation is that a power-law distribution does not have a characteristic temperature. We can however assume that each electrons emits energy at a given frequency (given that the synchrotron emission of individual electrons is highly peak at a given frequency):


$$\nu_c \equiv \frac{3}{2}\gamma^2\nu_{gyr} \sin \alpha$$


in order to calculate an *effective temperature* we can use that for an ultrarelativistic gas we have the relation:

$$E_e = 3kT_e(\nu)$$

where $E_e = \gamma m_e c^2$, since $\gamma \propto \nu^{1/2}$ we have that:

$$T_e(\nu) \propto \nu^{1/2}$$

The black-body radiation or Rayleigh-Jeans is proportional to $\nu^2$ but given the extra dependency of $T_e(\nu)$ it is changed as:

$$L(\nu) \propto T(\nu) \nu^2 B^{-1/2} \sim \nu^{1/2}\nu^{2}B^{-1/2} \sim \nu^{5/2} B^{-1/2}$$

A thing to note is that this part of the spectrum is independent of the orginal electron power-law spectral index.

:::{.callout title="Tutorial I: Working with SED using `naima`"}

We are going to plot some SED using a python module called [```naima```](https://github.com/zblz/naima). This package allows for the calculation of non-thermal radiation from relativistic particle populations.

```{python}
%matplotlib inline
import numpy as np
import matplotlib.pylab as plt
plt.rcParams['font.family'] = "STIXGeneral"
plt.rcParams.update({'axes.labelsize': 20})
plt.rcParams.update({'legend.fontsize': 20})
plt.rcParams.update({'figure.figsize': [8, 6]})
plt.rcParams['xtick.labelsize'] = 18
plt.rcParams['ytick.labelsize'] = 18
plt.rcParams['xtick.major.width'] = 1.5
plt.rcParams['xtick.minor.width'] = 1
plt.rcParams['xtick.major.pad'] = 8
plt.rcParams['xtick.direction'] = 'in'
plt.rcParams['ytick.major.size'] = 10
plt.rcParams['ytick.minor.size'] = 5
plt.rcParams['ytick.major.width'] = 1.5
plt.rcParams['ytick.minor.width'] = 1
plt.rcParams['ytick.major.pad'] = 8
plt.rcParams['ytick.direction'] = 'in'
plt.rcParams['legend.frameon'] = False
plt.rcParams['lines.linewidth'] = 1.5
plt.rcParams['axes.linewidth'] = 1.5



%config InlineBackend.figure_format = 'svg'
import naima
import astropy.units as u

# First we are going to define a population of electrons with a PowerLaw model                                    
e_pl = naima.models.PowerLaw(1e35*u.Unit('1/eV'), 1*u.TeV, 2.4)

#Now the equivlent synchrotron radiation for a magnetic field of 100 uG
syn = naima.models.Synchrotron(e_pl, B=100*u.uG)


spectrum_energy = np.logspace(-1,16,1000)*u.eV
sed_syn = syn.sed(spectrum_energy, distance=1.5*u.kpc)

fig = plt.figure(figsize=(10,4))
ax = plt.subplot(111)
ax.loglog(spectrum_energy,sed_syn,lw=3,label='Syncrotron')

ax.set_ylim(1e-9, 1e-5)
ax.set_xlabel('Photon energy [{0}]'.format(
        spectrum_energy.unit.to_string('latex_inline')))
ax.set_ylabel('$E^2 dN/dE$ [{0}]'.format(
        sed_syn.unit.to_string('latex_inline')))
ax.grid()
ax.legend(loc='upper right')
plt.show()
```
:::

### Inverse Compton scattering.{#inverse-compton}

Compton discovered that photons can transfer part of their energy to electrons in a collision. In astrophysics the inverse Compton effect electrons accelerated to high energy collide with photons from the blackbody radiation (thermal photons) or starlight photons and transfer energy to them. 

In astrophysics the inverse compton scattering is very important as a fast electrons that produce the synchrotron component can hit a low-energy photon and transfer a large fraction of its energy to the photon. For low energy photons, the scattering is elastic, and this regime is called the classical *Thomson regime*. However, for high energy photons where $E_\gamma \gg m_e c^2$ the scatting is not elastic and the energy of the scattered photons changes. This regime is called the *Klein-Nishina*. The cross-section for this *inelastic* scattering is given by:

$$\frac{\mathrm{ d} \sigma_{KN}}{\mathrm{ d}\Omega} = \frac{r_e^2}{2} \frac{E^2_{\gamma,out}}{E^2_{\gamma,in}}\left(\frac{E_{\gamma,out}}{E_{\gamma,in}} + \frac{E_{\gamma,in}}{E_{\gamma,out}} - \sin^2\theta\right)$$

In the low energy regime $E_{\gamma,in} = E_{\gamma,out}$ and the cross-section approaches that of the Thomson scattering.

#### Maximum energy of inverse Compton

Let's assume the lab reference system, $O$, and the electrons rest frame system $O^\prime$. The energy of the photon is then given by:

$$E_{\gamma,in}^\prime = \gamma E_{\gamma,in} (1 + \beta \cos\theta)$$

where $\gamma$ and $\beta$ refer obviously to the electron. Since in the electron's rest frame the scattering happens at low-energy we can assume that the scattering is described by the elastic Thomson regime, ie, the photon is scattered with the same energy in the electron's rest frame. Going back to the lab frame we have then:

$$E_{\gamma,out} = \gamma E_{\gamma,out}^\prime ( 1 + \beta\cos\theta^\prime)$$

but $E_{\gamma,in}^\prime = E_{\gamma,out}^\prime$. The maximum energy transfer will be that of a head-on collision in which the electron $\cos\theta = \cos \theta^\prime = 1$ is bounced backwards we have that:

$$\nu_{max} = \nu_0 \gamma^2 ( 1 + \beta)^2 \simeq 4 \nu_0 \gamma^2$$

As electrons are moving in an isotropic photon field the will see mostly head on collitions. So the energy spectrum of of scattered photons peak close to the maximum frecuency. This relation is very important for astrophysics, because it tell us a relation between the energy of electrons and the spectrum of photons. We know there are electrons with $\gamma \sim 100 - 1000$ and therefore the resulting photon emission is of very high energy. For example assuming electrons with $\gamma \approx 1000$ we have:


| Waveband | Frequency (Hz) $\nu_0$  | Scattered Frequency (Hz) |
|------|------|------ | 
|   Radio  | $10^9$ |   $10^{15}$ = UV    |
|   Far-infrared  | $3\times10^{12}$ |   $3\times10^{18}$ = X-rays    |
|   Optical  | $4\times10^{14}$ |   $4\times10^{21}$ = $\gamma$-rays    |


The hand-waving argument $\nu \sim \nu_0 \gamma^2$ has also implication on the spectrum of photons. If electrons follow a power-law in the form of:

$$\mathrm{d} n(E) \propto E^{-p} \mathrm{d} E$$

we can derive the spectrum of scattered photons which will follow a power-law as:

$$ I(\nu) \propto \nu^{-\frac{p-1}{2}} $$

because electrons losse energy by a factor $\gamma^2$ and the frequency of photons is $\nu \approx \nu_0 \gamma^2$

::: {.callout-note}
In some cases the synchrotron spectrum can interact with the same electron population that generated them via IC,  scenario is called *Synchrotron Self Compton* (SSC). When external photons are present IC does not need the synchrotron field. Such process is referred as *External Inverse Compton*
:::

:::{.callout title="Tutorial II: Adding the IC "}

We are going now to calculate the IC spectrum using the CMB as the radiation field:

```{python}
from astropy.constants import c

IC = naima.models.InverseCompton(e_pl, seed_photon_fields=['CMB'])

sed_ic_cbm = IC.sed(spectrum_energy, distance=1.5*u.kpc, seed="CMB")

fig2 = plt.figure(figsize=(10,4))
ax = plt.subplot(111)
ax.loglog(spectrum_energy,sed_syn,lw=3,label='Syncrotron',c=naima.plot.color_cycle[1])

ax.loglog(spectrum_energy,sed_ic_cbm,lw=3,
        label='Inverse Compton CBM', c=naima.plot.color_cycle[2])

ax.set_ylim(1e-12, 1e-5)
ax.set_xlabel('Photon energy [{0}]'.format(
        spectrum_energy.unit.to_string('latex_inline')))
ax.set_ylabel('$E^2 dN/dE$ [{0}]'.format(
        sed_syn.unit.to_string('latex_inline')))
ax.grid()
ax.legend(loc='upper right')
plt.show()
```

#### Case of study: Active Galactic Nuclei

As we saw, blazars are a particular classification of AGNs with their relativistic plasma jet oriented close to the line of sight. Because blazars have their jets pointing towards, the full emission comes from a relatively small region (the width of the jets) and so these objects can exhibit rapid variability (due to the smaller causally connected regions). The typical SED for AGNS is characterized by a synchrotron component extending from radio to X-ray frequencies, and a second component peaking at gamma-ray frequencies due to either inverse-Compton radiation (or from hadronic processes). The plot below shows the spectral energy distribution of nearby (z = 0.044) blazar 1ES 2344+514 during a flaring and a quiet state.
:::

#### Lightcurves and flares

In addition to the SED one could measure the light-curve or number of photons (or energy density) as a function of time. This is usually done to illustrate the variability of the source. For AGNs it is possible to make the lightcurve for different energy ranges, and therefore addressing different emission mechanism. 

Sometimes it is possible to observe  a *flare* in high energy photon range but with no activity detected in soft/hard X-ray. These are called **orphan flare**. This is interesting since in the SSC scenarios the X-rays (due to synchrotron) and TeV gamma-rays (due IC scattering) must be correlated. An orphan flare is therefore an indication that SSC scenario might be excluded, while the EC is still possible. As we will see also an hadronic scenario in which TeV photons come from $\pi^0$ is still possible.

## Hadronic models
### Pion decay

Accelerated protons in the source can produce charged and neutral pions via proton-proton interaction. In this context a CR proton with energy $E_p$ and momentum $p_p$ interacts with the ISM at rest in the process:

$$ p + p \rightarrow p + p + \pi^0$$ 


There is a threshold for this production given by the invariant:

$$s = E_{CoM}^2 = (E_p + m_p)^2 - p_p^2 = (2m_p^2 + m_{\pi^0}^2)^2 $$

As we saw in lesson 2, the available energy to produce particles can be derived as:

$$\epsilon = E_{CoM} - 2m_p $$

The energy threshod to produce a pion at rest is then:

$$E_{p,th} = \frac{m_{\pi^0}^2 + 4m_p m_{\pi^0} + 2m_p^2}{ 2 m_p} = \frac{m_{\pi^0}^2}{2m_p} + 2 m_{\pi^0} + m_p$$

```{python}
mp = 938.28 #MeV
mpion = 134.976 #MeV
print (f"Kinetic energy threshold for pion production {(mpion**2/(2*mp) + 2*mpion):.1f} MeV")
```

The pions are produced with the same power law as the parent proton. Neutral pions decay rapidly ($\tau = 8.4 \times 10^{-17}$ ns, compared to 26 ns of charged pions) in to two photons. In the CoM each photon from the pion decay has an energy of $E^{*}_\gamma = m_{\pi^0}/2 \approx 70$ MeV as we saw for the 2-body decay. In order to calculate the energy in the lab system we need to make a lorentz transformation:

$$E_\gamma = \gamma(E^{*}_\gamma + \beta_\pi p^{*}_\gamma \cos\theta^{*})$$

where ($^*$) denotes the CoM.  Since pion is an isoscalar the decay is isotropic the limits are determined by $\cos\theta^{*} = \pm 1$. These limits are:

$$\frac{m_{\pi^0}}{2}\sqrt{\frac{1-\beta}{1+\beta}} \leq E_\gamma \leq \frac{m_{\pi^0}}{2}\sqrt{\frac{1+\beta}{1-\beta}}$$

where $\beta$ is the velocity of the parent pion and we used $p^* = E^*$ for photons. 

#### Pion bump

But what is the distribution of photons? Since the decay is isotropic in the CoM:

$$\mathrm{d} N  = \frac{1}{4\pi} d\Omega = \frac{1}{2} \mathrm{d}\cos\theta^*$$

Using the lorentz transformation from CoM system and the lab system we can express: $\mathrm{d}\cos\theta^* = \mathrm{d}E_\gamma/(\gamma \beta_\pi p^*_\gamma) = 2\mathrm{d} E/ (\gamma \beta_\pi m_{\pi^0})$ and thus:

$$\mathrm{d} N = \frac{\mathrm{d} E}{\gamma \beta_\pi m_{\pi^0}}$$ 

Then the distribution of photons $\mathrm{d}N/\mathrm{d}E = \mathrm{\; const.}$ ie is constant on a box from $E_{min}$ to $E_{max}$ when plotted as $\log(E_\gamma)$. If many pions are decaying, the distribution of photons will be superpossition of boxes around $m_{\pi^0}$. This is the so-called *pion bump*.

![](images/pion_bump.png)

#### Pion bump in SNR with Molecular Cloud

In 2013 Fermi (a gamma-ray satellite) confirmed the pion-bump in two old SNRs with Molecular Cloud. 

![](https://inspirehep.net/files/6ff98851515f2a9d31b6faabd53dc1b7)

So did we find the sources of Galactic Cosmic Rays? No. The steep gamma-ray spectrum at high energies suggests that the acceleration is not very active any more (as expected from old SNRs) and hence, even though there might be some particle acceleration these are not the sources of Galactic Cosmic Rays. 

## Galactic gamma-ray diffuse emission

Diffuse gamma-ray emission is that not associated with a particular source. 
Fermi-LAT observed gamma-ray counts in the energy range from 200 MeV to 100 GeV. The signal is dominated by the diffuse Galactic emission, which is strongest in the plane of our Galaxy and toward the Galactic center but present all over the sky. The following shows the diffuse gamma-ray emission from 200 MeV to 100 GeV.

The gamma-ray diffuse emission are produced primarly the $\pi^0$ produced by interactions of cosmic rays protons with the ISM. The Inverse Compton scattering on star light is less important as well as bremsstrahlung (braking radiation). In the plot below the models are split into the three basic emission components:

* $\pi^0$-decay (red, long-dashed) 
* IC (green, dashed),
* and bremsstrahlung (cyan, dash-dotted). 

![Source: Taken from [arXiv:1202.4039](http://arxiv.org/abs/1202.4039)](images/GalacticCenter_Fermi.jpg)

## Extragalactic Background Light

As we saw, VHE gamma-rays ($E_\gamma =30$ GeV) have a limited horizon due to their interaction with the CMB and the [Extragalactic Background Light (EBL)](http://en.wikipedia.org/wiki/Extragalactic_background_light). The EBL is difficult to measure directly due to strong foregrounds from our solar system and the Galaxy. The TeV signal of distant AGNs are (partially) absorved at the highest energies by the ELB. The absortion is energy dependent and increases with distance. Observations of very far distant AGNs can be used to put constrains on the amount of EBL. However it is difficult to distinguish between an intrinsic softening of blazar spectra and a softening caused by the interaction with low energy EBL photons. The plot below shows the gamma-ray horizon and measurements from some Imaging Atmospheric Cherenkov Telescopes (IACTs) up to z=0.536.

![Source: Taken from [arXiv:0904.0774v2](http://arxiv.org/abs/0904.0774)](images/ELB.png)

## Gamma-ray Detection
### The atmosphere


The Earth's atmosphere is opaque to photons with energy above 10 eV, meaning that to observe gamma-rays directly we need to place the detector above the atmosphere. A turning point in gamma-ray astronomy was the launch of the first satellite-borne telescope, SAS-2, in 1972.

![Source: scienceofdoom.wordpress.com](http://scienceofdoom.files.wordpress.com/2009/11/800px-atmospheric_electromagnetic_opacity-svg.png)

The mayor limitation of satellite experiments is their low area which limits their use to $\leq 100$ GeV. At 100 GeV is when electromagnetic showers from the initial photon can be detected in ground-based telescopes.

### Cherenkov emission
When relativistic particles traverse a medium at a speed greater than the speed of light in that medium it can induce Cherenkov radiation. Cherenkov light is emitted in the UV and blue region in a narrow cone with angle:

![](images/Cherenkov.svg){fig-align="center"}

$$\cos \theta = \frac{ct/n}{\beta c t} = \frac{1}{\beta n}$$

In the relativistic limit $\beta \sim 1$ we can write:

$$\sin \theta = \sqrt{1 - \frac{1}{n^2}} = \sqrt{\frac{n^2 - 1}{n^2}} = \sqrt{\frac{(n-1)(n+1)}{n^2}} \approx \sqrt{2(n-1)}$$

In the atmosphere the spectral index depends on the altitude (same as density and pressure):

$$n = 1 + \epsilon_0 e^{-h/h_0} $$

where $\epsilon_0 \simeq 2.8 \times 10^{-4}$ at sea level. 

:::{.callout title="Tutorial III: Plot of the Cherenkov angle as funcion of the spectral index"}

```{python}
diffn = np.logspace(-5, 0, 100) #n - 1 form 0.0001 to 1
fig, ax = plt.subplots(figsize=(5,4))

#lets remove the values that will give errors in the arcsin:
diffn = diffn[np.where(np.sqrt(2*diffn) <= 1.)]

#In the relativistic limit 
ax.plot(diffn, np.degrees(np.arcsin(np.sqrt(2*diffn))))
ax.set_yscale("log")
ax.set_xscale("log")
ax.set_xlabel("n - 1 ")
ax.set_ylabel(r"$\theta_C$ ($^\circ$)")
ax.grid()
plt.show()
thetaw = np.degrees(np.arcsin(np.sqrt(2*(1.33 -1))))
print(r"Cherenkov angle in water: $\theta_C = %.2f^{\circ}$"%thetaw)
#For the atmosphere n = 1.0003
thetaatm = np.degrees(np.arcsin(np.sqrt(2*(1.0003 -1))))
print(r"Cherenkov angle in the atmosphere: $\theta_C = %.2f^{\circ}$"%thetaatm)


```
:::

#### Number of photons

The number of photons from Cherenkov radiation follows the *Frank-Tamn* formula:

$$\frac{\mathrm{d}N}{\mathrm{d}x} = 2 \pi \alpha z^2\int_{\lambda_1}^{\lambda_2} \sin^2\theta \frac{\mathrm{d}\lambda}{\lambda^2}$$

For relativistic particles, $\beta \sim 1$ in water $n = 1.33$ the Cherenkov spectrum is: 

```{python}
import scipy.constants as cte
alpha = cte.alpha #fine structure constant

fig, ax = plt.subplots(figsize=(5,4))
wl = np.arange(200, 500,1)
n = 1.33
ax.plot(wl, 2 * np.pi * alpha * ( 1- 1/n**2)/wl**2)
ax.set_xlabel(r"$\lambda$ [nm]")
ax.set_ylabel(r"$\frac{\mathrm{d}N_\gamma}{\mathrm{d}x\mathrm{d}\lambda}$ (photons nm$^{-2}$)")
ax.grid()
plt.show()
```

#### Atmospheric extintion.

Once in the atmosphere, the Cherenkov light can suffer absortion and scattering due to several processes:

* Rayleigh scattering:

$$ I = I_0 \frac{8\pi^4 \alpha^2}{\lambda^4 R^2}(1 - \cos^2\theta)$$

where $R$ is the distance to the particle and $\theta$ the scattering angle.

* Absortion by ozone, above 20 km.
* Aerosol, dust. Independent of $\lambda$.

![](images/Cherenkov_NSB_spectra.png)

This processes change the Cherenkov spectrum dependending on the altitude. In ice on the other hand, there is no UV absorption and the Cherenkov spectrum keeps the $1/\lambda^2$ tendency.