chap3_0.qmd

# Interactions of CR particles in the atmosphere
## Cosmic Ray Air Shower

In this lesson we will the different components of a cosmic air shower when cosmic rays interact with Earth's atmosphere. As seen in the figure below, we have 3 main components, the **electromagnetic component**, the **hadronic component** and muons and neutrinos which can be seen as the **muonic component**. 

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

### The Atmosphere

Before studying the interactions of cosmic rays in the atmosphere we need to setup a model that will describe our atmosphere. To study the cosmic rays interactions in the atmosphere it is useful to define a parameter that we will call the **vertical atmospheric depth** (sometimes also called column density) defined as the integral in altitude of the atmospheric density $\rho$ above the observation level $h$:    

$$X(h) =\int_h^{\infty} \rho(h^\prime)\mathrm{ d}h^\prime$$


![](images/x.svg){fig-align="center" width="50%"}

#### The Isothermal Model of the Atmosphere

In an isothermal hydrostatic atmosphere a particular layer of gas at some altitude is static. That means that the downward (towards the planet) force of its weight, plus the downward force exerted by pressure in the layer above it, and the upward force exerted by pressure in the layer below, all sum to zero. Assuming a segment of area $A$ and heigth $\mathrm{ d}h$ we can write this equilibrium of forces as:

$$P\cdot A - (P + \mathrm{ d}P)\cdot A - (\rho A \mathrm{ d}h)g_0  =0$$

$$\mathrm{ d}P = -g_0\rho(h) \mathrm{ d}h$$

Using the ideal gas law:

$$ P = \frac{\rho R T}{M}$$

where $R$ is the ideal gas constant, $T$ is temperature, $M$ is average molecular weight, and $g_0$ is the gravitational acceleration at the planet's surface. We get

$$\frac{\mathrm{ d}P}{P} = -\frac{g_0 M}{R T} \mathrm{ d} h$$

assuming a constant and isothermal gas (const $T$) we can integrate a pressure decreases exponentially with increasing height as:

$$P = P_0 e^{-\frac{g_0 M}{R T} h} $$

where the we can define the **scale height** as:

$$ h_0 = \frac{R T}{M g_0} $$

Since the temperature is assumed to be constant it follows that $\rho$ also changes exponentially as $\rho = \rho_0 e^{-h/h_0}$ and therefore the column density  can be written as:

$$X = X_0 \mathrm{ e}^{-h/h_0}$$ 

where $X_0$ is 1030 g/cm$^{2}$ is the atmospheric depth at sea level, h = 0. In particular for the isothermal model we have that the relation between atmospheric depth (aka column density) and density is:

$$ \rho(X) = \frac{X}{h_0} $$


##### The Scale Height

Using typical values ($T = 273$ K and $M = 29$ g/mol) we get that $h_0 \sim 8$ km  which coincidentally is the approximate height of Mt. Everest.  

In reality the temperature changes and hence the scale height decreases with increasing altitude until the tropopause.

This equations are valid for vertical particles, for zenith angles < 60$^\circ$ (for which we can ignore the Earth's curvature) the formula is scaled with $1/\cos{\theta}$ giving the *slant depth*

$$X_{slant-depth} = \frac{X}{\cos \theta} $$


### Energy Losses in the Atmosphere

Charge particles when entering in the atmosphere will suffer different process of energy losses. We are going to review some of them

#### Ionization Losses

The **ionization energy loss** of high energy charged particles with colission with atomic electrons is given by the Bethe-Block formula:

$$\left(\frac{\mathrm{ d} E}{\mathrm{ d} x}\right)_{ion}= -\left(\frac{4\pi N_0z^2e^4}{mv^2}\right)\left(\frac{Z}{A}\right)\left\{\log\left[\frac{2mv^2\gamma^2}{I}\right]-\beta^2\right\}$$

where $m$ is the mass of the electron, $v$ and $ze$ are the velocity and charge of the incoming particle, $N_0$ is the Avogadro's number, $Z$ and $A$ are the atomic and mass numbers of the atmos in the medium and $x$ the path travelled, and $I$ is the ionization potential of the medium is approximatively 10 Z eV.

* Since $Z/A \sim \frac{1}{2}$ in most materials it depends little on the medium.
* It varies as $1/v^2$ at low speed and independent of the incident particle mass.
* It reaches a minimum at about $3Mc^2$ and it increases logarithmically until it reaches a plateu value.

####  Radiation Losses
In addition to ionization losses, charge particles also undergo **bremsstrahlung** or braking radiation when travelling through a material given by:

$$\left(\frac{\mathrm{ d} E}{\mathrm{ d} x}\right)_{rad} = -\frac{E}{X_0}$$ 
where for electrons the **radiation length** is:

$$\frac{1}{X_0} = 4 \alpha\left(\frac{Z}{A}\right)(Z + 1)^2r_e^2N_0\log\left(\frac{183}{Z^{1/3}}\right)$$

where $r_e = e^2/4\pi m_e c^2$ is the classical electron radius and $\alpha = \frac{1}{137}$ is the fine structure constant. 

* Bremsstrahlung is proportional to $\frac{1}{X_0} \propto r^2_{e} \propto 1/m_e^2$. The radiation length of a muon will be $(m_\mu/m_e)^2$ times that for an electron.
* Bremsstrahlung is proportional to the energy.

The critical energy is the energy at $(\mathrm{ d}E/\mathrm{ d}x)_{ion} = (\mathrm{ d}E/\mathrm{ d}x)_{rad}$. Above this energy the radiation process dominates, below the ionization. For electrons this is roughly $\epsilon_c \sim 600/Z$ MeV, and for the atmosphere this is $\epsilon_e \sim 85$ MeV. 

####  Cherenkov Radiation

When relativistic particles traverse a medium at a speed greater than the speed of light in that medium it can induce Cherenkov radiation. 

![](images/Cherenkov.svg)


Cherenkov light is emitted in the UV and blue region in a narrow cone with angle:

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

so the threshold for production is $\beta > \frac{1}{n}$. Most of the components in the air shower will produced abundant Cherenkov light.

We will see more on Cherenkov radiation on the next lesson about $\gamma$-ray astronomy.


#### Pair Production

If a photon from bremsstrahlung has enougth energy $E_\gamma > 2 m_e$ it can produced a pair of electron positron. The cross-section rises quickly at the threshod of $2 m_e$ but in the high energy part it can be approximated to:

$$\sigma = \frac{7}{9} r_0^2 Z (Z+ 1)\log \left(\frac{183}{3\sqrt{Z}}\right)$$

The pair production cannot occure in vacum, a photon desintegrated to the pair $e^-e^+$ will have a null momentum in the CoM system, therefore a nucleus has to be present to absorb the momemtum. 
As can be seen the radiation length $X_0$ is very similar to the one from radiation losses. In fact we can write:

$$\frac{1}{X_{pair}} = \frac{7}{9}\frac{1}{X_0}$$

Which means that the **radiation lengths for braking radiation and pair production are comparable**.

## Electromagnetic Shower

When photons from radiation losses of electrons, have enough energy to produce pairs of positrons electrons, these hey also can produce photons which, in turn, can also produces pairs, etc, etc. This is called an **electromagnetic shower**.

### The Heitler Toy-Model 

The Heitler toy model explains very well the development of an electromagnetic shower.
As we saw, in the ultrarelativistic limit the radiation lenghs for pair production and bremsstrahung are comparable. We can defined a distance $d_{split} = X_0 \log2$ where an electron will looses, on average, half of its energy. 

An electron with inital energy $E_0$ in a medium will generate a photon in a $d_{split}$ length of energy $E_0/2$, in the next radiation length the photon can convert into $e^+e^-$ each with energies $E_0/4$. After $t$ steps the electrons, positrons will have energies of $E(t)=E_0/2^t$. This continues until the electrons, positrons fall below the critical energy of electrons, $\epsilon_e$, and ionization dominates. The process is illustrated in the figure below, where each step $n$ corresponds to one $d_{split}$ length.

![](images/em-shower.svg){fig-align="center" width="50%"}

The Heitler model has the following properties:

* The shower has maximum at:

$$t_{max} = \frac{\log(E_0/\epsilon_e)}{\log2}$$

* The maximum number of particles is:

$$N_{max} = 2^{t_{max}} = \frac{E_0}{\epsilon_e}$$


* The shower maximum will be at a depth $X_{max}$:

$$X_{max}=d_{split} \frac{\log(E_0/\epsilon_e)}{\log2} = X_0 \log(E_0/\epsilon_e)$$

For air $\epsilon_e = 85$ MeV and the **radiation length** $X_0 = 36.7 \mathrm{ g/cm}^2$. Actual showers also spread laterally mostly due to Coulomb scattering. The lateral spread is a few times the so-called **Moliere unit** equal to $21/\epsilon_e$ (MeV).

The $X_{max}$ prediction of the Heitler model is in good agreement with Monte Carlo simulations. However, the electron to photon ratio of 2 is not in agreement given that the model predicts only one photon emitted by bremsstrahlung. Simulations show a ratio of 1/6 since in reality several photons are emitted and electrons lose energy much faster than photons do. 

## Hadronic Showers

Before modeling the baryon-induced showers or hadronic showers, we need to estimate the *nuclean mean free path* in the atmosphere. We can write the **standard mean free path** for $\sigma^{air}_N$ nucleon-air cross section as:

$$l_N = \frac{1}{n \sigma_N^{air}}$$

where as we saw in the introduction $n$ is the **number density** of targets, in our case air nuclei. This number density can be expresed as $n = N_T/V $ where $N_T$is the total number of air nuclei in volume $V$.

The **density mean free path**, is defined as $\lambda_N = \rho_{air} l_N$, where $\rho_{air} = n \cdot m_{air}$ where $m_{air}$ is the mass of air nuclei, that can be written as $m_{air} \sim A m_p$. Puttin everything together we have that:

$$\lambda^{air}_N = \rho_{air} l_N = n A m_p\frac{1}{n\sigma_N^{air}} = \frac{A m_{p}}{\sigma_N^{air}}$$


For air $A$ is average mean the mass number of air nuclei components (mainly nitrogen, oxygen) and we can assume it to be $A \sim 14.5$ and $\sigma_N^{air} \approx$ 300 mb, which corresponds to $\lambda_N^{air} \approx 80$g/cm$^{2}$.


Note that this definition of the density mean free path is independent of the mass density of the medium, so if the density changes with altitude, like in the case of our atmosphere, the density mean free path is the same.


### The Heitler-Matthew Model

In the Heitler model can be adapted also for hadronic showers. This is what Matthew did. We can imagine a proton initiating the cascade instead of a photon/electron, in this case a hadronic air shower will develope:

![](images/proton_eas.png){width=50% fig-align="center"}

We can assume that the first interaction is defined by the proton mean free path $\lambda^{air}_N$. Defining the first interaction point where the proton will lose (on average) half of its energy this first interaction length is given as $\lambda^{air}_{N} \log 2$ where for protons $\lambda^{air}_N \approx 80$g/cm. The following general interaction is expected:

$$ p + p/n \rightarrow p + n/p + \pi^0 + \pi^{\pm} + K^0 + K^{\pm} + .. $$

We will focus only on pion production (same argument can be done for kaons). After the first interaction we can use the simplified assumption that the hadronic interaction produces only $3N_\pi$ pions. Of those $ 2 N_\pi$ will be charged pions and $N_\pi$ will be $\pi^0$. 

![](images/hadronic-shower.svg){fig-align="center" width=80%}


We also assume the energy is equally distributed among them, so $2/3E_0$ will go to charge pions and $1/3E_0$ will go to the neutral pions. The $\pi^0$ has a very short decay time, so it will decay and produce an electromagnetic shower. Charge pions will continue generating hadronic shower in each $d_{split} = \lambda_{\pi}^{air} \log 2$ with the mean free path of pions $\lambda_\pi^{air} \sim 120$g/cm$^2$ until they reach the critical energy where pions decay is more probable than interactions $\epsilon_\pi$. On each step we assume that energy is equally divided among the $3 N_\pi$ pions. Thefore at each step $t$ the energy of the pions is:

$$E_\pi = \frac{E_0}{(3N_\pi)^t} $$

The number of radiation lengths $t$ to reach the critical energy ie $E_\pi = \epsilon_\pi$, and is given (as in the case of EM showers):

$$t_{max} = \frac{\log(E_0/\epsilon_\pi)}{\log(3N_\pi)}$$ 


Assuming that after that energy all charged pions (ie $2N_\pi$) decay to muons, the number of muons is given by:

$$N_\mu =(2N_\pi)^{t_{max}}$$

introducing $\beta = \log(2N_\pi)/\log(3N_\pi)$ we have:

$$N_\mu =(E_0/\epsilon_\pi)^\beta $$

This is also called the *multiplicity* and corresponds to the muon bundles as we will see later.
For pions between  1 GeV and 10 TeV an appropiate number is $N_\pi = 5$ and in that case $\beta = 0.85$. Therefore the number of muons doesn't grow linearly with the intial energy but a slower rate.

The definition of $X_{max}$ is somehow less clear than in an EM shower. Hadronic showers are still dominated by electromagnetic processes, so we can assume that $X_{max}$ depends dominantly on the first generation of $\pi^0$ EM subshowers. For proton primaries, the first interaction will be given by the nucleon mean free path where in this first interaction the proton splits in $3N_\pi$ particles, so $pi^0$ with initial energy $E_0/3N_\pi$ will initiate an EM shower. The depth of maximum is then obtained as the sum of the first proton interaction length and the shower maximum of the first EM sub-shower: 

$$X_{max} = \lambda^{air}_N \log2 + X_0 \log\left(\frac{E_0}{3N_\pi\epsilon_e}\right) $$

where again $\epsilon_e$ is the critical energy of electrons. The expected values of this formula are low when compared to detailed simulations because it neglects the contributions of the next one or two generation of $\pi^0$ production.

:::{.callout-note title="Superposition model for heavy nuclei air showers"}

We can extend the discussion to heavy nuclei by adopting the *Superposition model* in which a nucleus of mass A and energy $E_0$ essentially generates $A$ subshowers of energy $E_0/A$. In that case the muon multiplicity will be:

$$N_\mu = A \left(\frac{E_0}{A\epsilon_\pi}\right)^\beta \propto E_0^\beta A^{1 - \beta} $$

therefore the muon multiplicity will depend on the CR composition. Likewise the shower maximum is given by:

$$X^{A}_{max}(E_0) = X_{max}^{p}(E_0) - X_0\log{A}$$

ie, for a given energy $E_0$ the shower max depends on the mass of the CR primary and it is typicall smaller than for protons (ie reach the maximum sooner in the atmosphere). For composition studies therefore it is necessary to measure $X_{max}$ and the energy of the shower, which can be estimated from fluorescence techniques.
:::