Lambert-Beer's law (Bouguer 1729; Lambert 1760; Beer 1852)

gives the intensity of radiation \(I(z)\) after traversing through a medium with constant attenuation \(μ\) of length \(z\), given a starting intensity of \(I_0\). Directly related is of course the absorption length \(l_{\text{abs}} = 1/μ\) (or mean free path), which is a useful property when considering typical absorption depths.

This law is of vital importance for the behavior of X-rays traversing through matter, which is needed to compute the efficiency of a gaseous detector with an entrance window. In addition, it is also related to the mean free path of X-rays in a gas, which is an important parameter in gaseous detectors to understand the absorption efficiency of X-rays at different energies and the resulting expected diffusion.

In the context of X-rays the factor \(μ\) is typically rewritten via the 'mass attenuation coefficient' \(μ_m = μ · ρ\) with \(ρ\) the density of the material, commonly in \(\si{g.cm^{-3}}\). \(μ_m\) is then defined by

\[ μ_m = \frac{N_A}{M} σ_A, \]

where \(N_A\) is Avogadro's number, \(M\) the molar mass of the medium in units of \(\si{g\per\mol}\) and \(σ_A\) is the photoabsorption cross section in units of \(\si{cm^2}\). Thus, the mass attenuation coefficient is usually given in \(\si{cm^2.g^{-1}}\) such that \(μ = μ_m · ρ\) is of inverse length as expected. This directly yields the definition of the absorption length,

\[ l_{\text{abs}} = \frac{1}{μ}. \]

Further, the photoabsorption cross section can be described via the atomic scattering factor \(f₂\)

\[ σ_A = 2 r_e λ f₂, \]

where \(r_e\) is the classical electron radius and \(λ\) the wavelength of the X-ray. \(f₂\) is the imaginary part of the forward scattering factor \(f\)

\[ f = f₁ - i f₂ \]

which itself is the simplification of the general atomic scattering factor that describes the atom specific part of the scattering cross section.

This way of expressing it has the nice property of relying on a well tabulated parameter \(f₂\). Together with \(f₁\) these tabulated values can be used to compute everything from the refractive index at a specific X-ray energy of a compound to the attenuation coefficient and even reflectivity of a multi-layer substrate. It generalizes from single element to compounds easily by

\[ μ_m = \frac{N_A}{M_c} \sum_i n_i σ_{A,i}, \]

with \(M_c\) the molar weight of the compound and \(n_i\) the number of atoms of kind \(i\). X-ray absorption and transmission properties can be calculated from this only requiring the atomic scattering factors, which can be found tabulated for different elements, for example by NIST and Henke.

There is an online calculator for calculations of X-ray transmission found under ² (Henke, Gullikson, and Davis 1993), as well as a library implementation developed during the course of this thesis under ³ (Schmidt and SciNim contributors 2022) for this purpose.

Fig. 1 shows an example of X-ray transmission through a \(\SI{300}{nm}\) thick layer of \ccsini as well as transmission through \(\SI{3}{cm}\) of argon at normal temperature and pressure (NTP), \(\SI{1}{atm}\), \(\SI{20}{°C}\). All information about the absorption lines and general transmission is encoded in \(f₂\).

The same atomic scattering factors \(f₁\) and \(f₂\) introduced in section 6.1.1 for the attenuation can also be used to compute the reflectivity of X-rays under shallow angles. A great overview of the relevant physics for X-ray reflectivity is found in (Windt 1998), which introduces the IMD program for the simulation of multilayer coatings for X-rays.

By defining the combined scattering factor

\[ f(E) = f₁(E) + i f₂(E) \]

at energy \(E\), the refractive index \(n\) of a medium can be computed using

\[ n(E) = 1 - r_e \frac{λ²}{2π} \sum_i n_{ai} f_i(E) \]

where \(n_{ai}\) is the number density of the \(i\text{-th}\) compound of the medium. By expressing the refractive index for X-rays in this fashion, reflectivity can be expressed using the Fresnel equations, just like for visible light. The reflectivity for s-polarization is calculated by

while for p-polarization it is done via

Here \(θ_i\), \(θ_k\) are the incident and refracted angles and \(n_i\), \(n_k\) the refractive indices on the incident and outgoing side \(i\) and \(k\), respectively.

The total reflected energy, the reflectance \(R\), is expressed as

\[ R = \frac{1}{2}\left( \left|r^s\right|² + \left|r^p\right|²\right) \]

for unpolarized light. This can be generalized to multiple layers of material on a substrate and including a surface roughness. Combined, these provide the essence for a realistic computation of the efficiency of an X-ray telescope mirror shell.

This is also implemented in ³ (Schmidt and SciNim contributors 2022) and ² (Henke, Gullikson, and Davis 1993) also provides an online calculator for such reflectivities.

Depth graded multilayers

One particular kind of surface, the depth graded multilayer, is used in certain kinds of modern X-ray telescopes, for example the LLNL telescope at CAST following the NuSTAR design. In such a multilayer, repeating layers of a low atomic number \(Z\) material and a high \(Z\) material are stacked at decreasing thicknesses.

A depth-graded multilayer is described by the equation: \[ d_i = \frac{a}{(b + i)^c} \] where \(d_i\) is the depth of layer \(i\) (out of \(N\) layers), \[ a = d_{\text{min}} (b + N)^c \] and \[ b = \frac{1 - N k}{k - 1} \] with \[ k = \left(\frac{d_{\text{min}}}{d_{\text{max}}}\right)^{\frac{1}{c}} \] where \(d_{\text{min}}\) and \(d_{\text{max}}\) are the thickness of the bottom and top most layers, respectively.

For example for the the LLNL telescope a Pt/C depth graded multilayer is used, in which a top layer of carbon is stacked on top of a platinum layer. Between 2 to 5 repetitions of decreasing thickness are stacked with a ratio of \(\SIrange{40}{45}{\%}\) carbon to platinum in thickness. More details on this will be discussed in appendix 37, as it is of vital importance to calculate the axion image required for the limit calculation correctly.

The reflectivity for a depth-graded multilayer is computed recursively from the bottom of the stack to the top layer using

\[ r_i = \frac{r_{ik} + r_k \exp(2 i β_i)}{1 + r_{ik} r_k \exp(2 i β_i)} \]

with

\[ β_i = 2π · d_i · \frac{\cos(θ_i)}{λ}, \]

where \(θ_i\) as seen from the normal axis and with the wavelength of the incoming X-rays \(λ\). The \(r_{ik}\) values are computed following equations \eqref{eq:theory:fresnell_reflectance_s} and \eqref{eq:theory:fresnell_reflectance_p}. Such multilayers work by summing the reflecting contributions from the different layer transitions. Different thicknesses of the different multilayers mean X-rays at different energies and angles are best reflected from different layers. Thus, a much improved overall reflectivity over a wider energy and angle range can be achieved compared to a normal single layer on a substrate (e.g. a gold coating as used for the XMM-Newton or ABRIXAS optics).

Another relevant aspect for gaseous detectors is the energy deposition of charged particles. In particular for experiments that sit near Earth's surface, a major source of background is due to cosmic radiation, with cosmic muons making up more than \(\SI{95}{\%}\) (Zyla and others 2020) of radiation (aside from neutrinos) at the surface, see sec. 6.2.

The energy loss of such muons can be calculated with the Bethe-Bloch equation, which describes the average energy loss per distance for a charged particle with charge \(z\) in a homogeneous medium with charge carriers \(Z\). (Zyla and others 2020) ⁴

where the different variables are as follows:

• \(K = 4π N_A r_e² m_e c² = \SI{0.307075}{MeV.mol^{-1}.cm^2}\)
• \(W_{\text{max}}\): maximum possible energy transfer to an electron in a single interaction
• \(I\): mean excitation energy of the absorber material in \si{\eV}
• \(δ(βγ)\): density-effect correction to energy loss
• and \(r_e = \frac{e²}{4π ε_0 m_e c²}\) the classical electron radius, \(N_A\) Avogadro's number, \(m_e\) electron mass, \(c\) speed of light in vacuum, \(z\) charge number of incident particle, \(Z\) atomic number of absorber material, \(A\) atomic mass of absorber material, \(β = v/c\) speed of incident particle, \(γ\) Lorentz factor

This interaction behavior of muons leads to a specific, expected energy loss per distance. Commonly they are 'minimum ionizing particles' (MIPs), as their energies lies between \(\SIrange{0.1}{100}{GeV}\), the large 'valley' in which the Bethe equation is applicable ⁵. For argon gas at normal conditions this is shown in fig. 4.

As the Bethe formula was derived from quantum mechanical perturbation theory, higher order corrections can be computed. For our purposes here the leading order is enough. It is important to keep in mind that the Bethe-Bloch equation gives the mean energy per distance. When considering short distances as typically encountered in particle detectors, this mean is skewed by rare interactions that deposit large amounts of energy (towards \(W_{\text{max}}\)). The energy deposition along short distances is typically described by a Landau-Vavilov distribution (similar, but different from a normal Landau distribution) (Zyla and others 2020; Bichsel 2006). The most probable energy loss is often a more appropriate number to look at. It can be expressed as

where \(ξ\) is

\[ ξ = \frac{1}{2} K z² \left⟨ \frac{Z}{A} \right⟩ \frac{x}{β²} \, \si{MeV}, \]

for a detector in which the material column the particle travels through is expressed as \(x = d · ρ\) of a distance \(d\) in \(\si{g.cm^{-2}}\). \(j = \num{0.200}\) is an empirical constant (Zyla and others 2020; Bichsel 1988). Further, \(⟨Z / A⟩\) is simply the average \(Z/A\) for a material compound \(⟨Z/A⟩ = \sum_i w_i Z_i / A_i\).

The large difference typically encountered between the most probable and the mean value for the energy loss in particle detectors, makes studying the expected signals a complicated topic. For a detailed description relevant for thin gaseous detectors, see especially (Bichsel 2006).

Fig. 4 shows the comparison of the most probable energy loss via equation \eqref{eq:theory:most_probable_loss} and the mean energy loss via the Bethe-Bloch equation \eqref{eq:theory:bethe_bloch_eq} for muons of different energies traversing \(\SI{3}{cm}\) of argon gas.

Cosmic muons in their interactions with matter can ionize atoms, leading to the possible emission of X-rays if the removed electron is part of an inner shell, mostly K (and some L) shell electrons. This leads to a form of background based on real X-rays and thus represents a kind of background that is impossible to distinguish from any kind of axion signal unless external scintillator based vetoes are used. Of course, to be relevant as a form of detector background the material must be close to the detector, as the X-rays will otherwise be absorbed. This makes the detector material, the gas itself and all material in the direction of the detectors' sensitivity a candidate for X-ray fluorescence background.

Table 1 lists the X-ray fluorescence lines for elements one commonly encounters in the context of gaseous detectors for a helioscope experiment. Table 2 lists the atomic binding energies of the same elements. This is intended as a useful reference to understand possible lines in background, potential targets for an X-ray tube for reference X-ray data and understand the absorption edges for materials.

The difference between the binding energy and the energies of the most likely fluorescence lines is the explanation for why commonly a material is more transparent for a fluorescence X-ray than energies slightly above it. This (among other effects) explains the 'escape peak' seen in many types of gaseous detectors. For example see fig. 1 for the argon transmission. The argon \(K 1s\) binding energy of \(\SI{3.2}{keV}\) is visible in the form of the sudden drop in transmission (which is of course directly proportional to the absorption length!). At the same time an X-ray produced by an ionized argon atom from the \(K 1s\) electron via the Kα line yields an X-ray of only \(\SI{2.95}{keV}\) and thus argon gas is extremely transparent for such X-rays (this is the cause of so called 'escape photons', see sec. 6.3.8).

The code here is directly based on code written in my notes ./../org/Doc/StatusAndProgress.html being tangled into a file /tmp/muon_flux.nim.

Aside from the figures included in the main thesis, this also produces fig. 6.

import math, macros, unchained
import seqmath, ggplotnim, sequtils, strformat

let K = 4 * π * N_A * r_e^2 * m_e * c^2 # usually in: [MeV mol⁻¹ cm²]

defUnit(cm³•g⁻¹)
defUnit(J•m⁻¹)
defUnit(cm⁻³)
defUnit(g•mol⁻¹)
defUnit(MeV•g⁻¹•cm²)
defUnit(mol⁻¹)
defUnit(keV•cm⁻¹)
defUnit(g•cm⁻³)
defUnit(g•cm⁻²)

proc electronDensity(ρ: g•cm⁻³, Z, A: UnitLess): cm⁻³ =
  result = N_A * Z * ρ / (A * M_u.to(g•mol⁻¹))

proc I[T](z: float): T =
  ## use Bloch approximation for all but Argon (better use tabulated values!)
  # 188.0 eV from NIST table 
  result = if z == 18.0: 188.0.eV.to(T) 
           else: (10.eV * z).to(T)

proc calcβ(γ: UnitLess): UnitLess =
  result = sqrt(1.0 - 1.0 / (γ^2))

proc betheBloch(z, Z: UnitLess,
                   A: g•mol⁻¹,
                   γ: UnitLess,
                   M: kg): MeV•g⁻¹•cm² =
  ## result in MeV cm² g⁻¹ (normalized by density)
  ## z: charge of particle
  ## Z: charge of particles making up medium
  ## A: atomic mass of particles making up medium
  ## γ: Lorentz factor of particle
  ## M: mass of particle in MeV (or same mass as `m_e` defined as)
  let β = calcβ(γ)
  let W_max = 2 * m_e * c^2 * β^2 * γ^2 /
    (1 + 2 * γ * m_e / M + (m_e / M)^2)
  let lnArg = 2 * m_e * c^2 * β^2 * γ^2 * W_max / (I[Joule](Z)^2)
  result = (K * z^2 * Z / A * 1.0 / (β^2) * (
    0.5 * ln(lnArg) - β^2
  )).to(MeV•g⁻¹•cm²)

proc mostProbableLoss(z, Z: UnitLess, A: g•mol⁻¹, γ: UnitLess,
                      x: g•cm⁻²): keV =
  ## Computes the most probable value, corresponding to the peak of the Landau
  ## distribution, that gives rise to the Bethe-Bloch formula.
  ##
  ## Taken from PDG chapter 'Passage of particles through matter' equation
  ## `34.12` in 'Fluctuations in energy loss', version 2020).
  ##
  ## `x` is the "thickness". Density times length, `x = ρ * d`. The other parameters
  ## are as in `betheBloch` above.
  let β = calcβ(γ)
  let ξ = K / 2.0 * Z / A * z*z * (x / (β*β))
  const j = 0.200
  let I = I[Joule](Z)
  result = (ξ * ( ln((2 * m_e * c^2 * β^2 * γ^2).to(Joule) / I) + ln(ξ.to(Joule) / I) + j - β^2)).to(keV) # - δ*(β*γ)

proc density(p: mbar, M: g•mol⁻¹, temp: Kelvin): g•cm⁻³ =
  ## returns the density of the gas for the given pressure.
  ## The pressure is assumed in `mbar` and the temperature (in `K`).
  ## The default temperature corresponds to BabyIAXO aim.
  ## Returns the density in `g / cm^3`
  let gasConstant = 8.314.J•K⁻¹•mol⁻¹ # joule K^-1 mol^-1
  let pressure = p.to(Pa) # pressure in Pa
  # factor 1000 for conversion of M in g / mol to kg / mol
  result = (pressure * M / (gasConstant * temp)).to(g•cm⁻³)

proc E_to_γ(E: GeV): UnitLess =
  result = E.to(Joule) / (m_μ * c^2) + 1

proc γ_to_E(γ: UnitLess): GeV =
  result = ((γ - 1) * m_μ * c^2).to(GeV)

type
  Element = object
    Z: UnitLess
    M: g•mol⁻¹
    A: UnitLess # numerically same as `M`
    ρ: g•cm⁻³
proc initElement(Z: UnitLess, M: g•mol⁻¹, ρ: g•cm⁻³): Element =
  Element(Z: Z, M: M, A: M.UnitLess, ρ: ρ)

# molar mass. Numerically same as relative atomic mass
let M_Ar = 39.95.g•mol⁻¹
let ρAr = density(1050.mbar, M_Ar, temp = 293.15.K)
let Argon = initElement(18.0.UnitLess, 39.95.g•mol⁻¹, ρAr)

proc intBethe(e: Element, d_total: cm, E0: eV, dx = 1.μm): eV =
  ## integrated energy loss of bethe formula after `d` cm of matter
  ## and returns the energy remaining
  var γ: UnitLess = E_to_γ(E0.to(GeV))
  var d: cm
  result = E0
  var totalLoss = 0.eV
  while d < d_total and result > 0.eV:
    let E_loss: MeV = betheBloch(-1, e.Z, e.M, γ, m_μ) * e.ρ * dx
    result = result - E_loss.to(eV)
    γ = E_to_γ(result.to(GeV))
    d = d + dx.to(cm)
    totalLoss = totalLoss + E_loss.to(eV)
  result = max(0.float, result.float).eV

proc plotDetectorAbsorption() =
  let E_float = logspace(-2, 2, 1000)
  let energies = E_float.mapIt(it.GeV)
  let E_loss = energies.mapIt(
    (it.to(eV) - intBethe(Argon, 3.cm, it.to(eV))).to(keV).float
  )
  let df = toDf(E_float, E_loss)
  ggplot(df, aes("E_float", "E_loss")) +
    geom_line() +
    xlab("μ Energy [GeV]") + ylab("ΔE [keV]") +
    scale_x_log10() + scale_y_log10() +
    ggtitle("Energy loss of Muons in 3 cm Ar at CAST conditions") +
    ggsave("/home/basti/phd/Figs/muons/ar_energy_loss_cast.pdf")
plotDetectorAbsorption()

let Atmosphere = @[(0.78084, initElement(7.0.UnitLess, 14.006.g•mol⁻¹, 1.2506.g•dm⁻³.to(g•cm⁻³))), # N2
                   (0.20964, initElement(8.0.UnitLess, 15.999.g•mol⁻¹, 1.429.g•dm⁻³.to(g•cm⁻³))),  # O2
                   (0.00934, initElement(18.0.UnitLess, 39.95.g•mol⁻¹, 1.784.g•dm⁻³.to(g•cm⁻³)))]  # Ar

proc plotMuonBethe() =
  let E_float = logspace(-2, 2, 1000)
  let energies = E_float.mapIt(it.GeV)
  var dEdxs = newSeq[float]()
  for e in energies:
    var dEdx = 0.0.MeV•g⁻¹•cm²
    for elTup in Atmosphere:
      let (w, element) = elTup
      let γ = E_to_γ(e)
      dEdx += w * betheBloch(-1, element.Z, element.M, γ, m_μ)
    dEdxs.add dEdx.float
  let df = toDf(E_float, dEdxs)
  ggplot(df, aes("E_float", "dEdxs")) +
    geom_line() +
    xlab("μ Energy [GeV]") + ylab("dE/dx [MeV•g⁻¹•cm²]") +
    scale_x_log10() + scale_y_log10() +
    ggtitle("Energy loss of Muons in atmosphere") +
    ggsave("/home/basti/phd/Figs/muons/energy_loss_muons_atmosphere.pdf")  
plotMuonBethe()
#if true: quit()
import math, unchained, ggplotnim, sequtils

const R_Earth = 6371.km
func distanceAtmosphere(θ: Radian, d: KiloMeter = 36.6149.km): UnitLess =
  ## NOTE: The default value for `d` is not to be understood as a proper height. It.s an
  ## approximation based on a fit to get `R_Earth / d = 174`!
  result = sqrt((R_Earth / d * cos(θ))^2 + 2 * R_Earth / d + 1) - R_Earth / d * cos(θ)

defUnit(cm⁻²•s⁻¹•sr⁻¹)  
defUnit(m⁻²•s⁻¹•sr⁻¹)
proc muonFlux(E: GeV, θ: Radian, E₀, E_c: GeV,
              I₀: m⁻²•s⁻¹•sr⁻¹,
              ε: GeV): m⁻²•s⁻¹•sr⁻¹ =
  const n = 3.0
  let N = (n - 1) * pow((E₀ + E_c).float, n - 1)
  result = I₀ * N * pow((E₀ + E).float, -n) *
    #pow((1 + E / ε).float, -1) *
    ( ( 1.0 / (1 + 1.1 * E * cos(θ) / 115.GeV).float) + (0.054 / (1 + 1.1 * E * cos(θ) / 850.GeV).float) ) * 
    pow(distanceAtmosphere(θ), -(n - 1))

from numericalnim/integrate import simpson
proc plotE_vs_flux(θ: Radian, E₀, E_c: GeV, I₀: m⁻²•s⁻¹•sr⁻¹, ε: GeV,
                   suffix = "") =
  let energies = linspace(E_c.float, 100.0, 1000)
  let E = energies.mapIt(it.GeV)
  let flux = E.mapIt(muonFlux(it, θ, E₀, E_c, I₀, ε).float) # .to(cm⁻²•s⁻¹•sr⁻¹)
  let df = toDf(energies, flux)

  echo "Integrated flux: ", simpson(flux, energies)
  
  ggplot(df, aes("energies", "flux")) +
    geom_line() +
    xlab("Energy [GeV]") + ylab("Flux [m⁻²•s⁻¹•sr⁻¹]") +
    scale_x_log10() + scale_y_log10() +
    ggtitle(&"Flux dependency on the energy of muons at θ = {θ.to(°)}{suffix}") +
    ggsave(&"/home/basti/phd/Figs/muons/energy_vs_flux_cosmic_muons{suffix}.pdf")
plotE_vs_flux(0.Radian,
              2.5.GeV, #4.29.GeV,
              0.5.GeV, 70.7.m⁻²•s⁻¹•sr⁻¹, 854.GeV)


let E₀ = 25.0.GeV
let I₀ = 90.0.m⁻²•s⁻¹•sr⁻¹
let E_c = 1.GeV
let ε = 2000.GeV

proc plotFlux_at_CAST() =
  let energies = linspace(0.5, 100.0, 1000)
  let E = energies.mapIt(it.GeV)
  let flux = E.mapIt(muonFlux(it, 88.0.degToRad.Radian, E₀, E_c, I₀, ε).float)
  let df = toDf(energies, flux)
  ggplot(df, aes("energies", "flux")) +
    geom_line() +
    xlab("Energy [GeV]") + ylab("Flux [m⁻²•s⁻¹•sr⁻¹]") +
    scale_x_log10() + scale_y_log10() +
    ggtitle("Flux dependency on the energy at θ = 88° at CAST altitude") +
    ggsave("/home/basti/phd/Figs/muons/flux_at_cast_88_deg.pdf")
plotFlux_at_CAST()

proc computeMeanEnergyLoss() =
  let energies = linspace(0.5, 100.0, 1000)
  let E = energies.mapIt(it.GeV)
  let flux = E.mapIt(muonFlux(
    it, 88.0.degToRad.Radian, E₀, E_c, I₀, ε).float
  )
  let E_loss = E.mapIt(
    (it.to(eV) - intBethe(Argon, 3.cm, it.to(eV))).to(keV).float
  )
  let fluxSum = flux.sum
  let df = toDf(energies, E_loss, flux)
      .mutate(f{"flux" ~ `flux` / fluxSum},
              f{"AdjFlux" ~ `E_loss` * `flux`})
  echo "Mean energy loss: ", df["AdjFlux", float].sum
computeMeanEnergyLoss()

proc computeHeight(S: Meter, θ: Radian): KiloMeter =
  ## For given remaining distance distance along the path of a muon
  ## `S` (see fig. 1 in 1606.06907) computes the remaining height above
  ## ground. Formula is the result of inverting eq. 7 to `d` using quadratic
  ## formula. Positive result, because negative is negative.
  result = (-1.0 * R_Earth + sqrt(R_Earth^2 + S^2 + 2 * S * R_Earth * cos(θ)).m).to(km)

import algorithm
defUnit(K•m⁻¹)
proc barometricFormula(h: KiloMeter): g•cm⁻³ =
  let hs = @[0.0.km, 11.0.km]
  let ρs = @[1.225.kg•m⁻³, 0.36391.kg•m⁻³]
  let Ts = @[288.15.K, 216.65.K]
  let Ls = @[-1.0 * 0.0065.K•m⁻¹, 0.0.K•m⁻¹]
  let M_air = 0.0289644.kg•mol⁻¹
  let R = 8.3144598.N•m•mol⁻¹•K⁻¹
  let g_0 = 9.80665.m•s⁻²
  let idx = hs.mapIt(it.float).lowerBound(h.float) - 1
  case idx
  of 0:
    # in Troposphere, using regular barometric formula for denities
    let expArg = g_0 * M_air / (R * Ls[idx])
    result = (ρs[idx] * pow(Ts[idx] / (Ts[idx] + Ls[idx] * (h - hs[idx])), expArg)).to(g•cm⁻³)
  of 1:
    # in Tropopause, use equation valid for L_b = 0
    result = (ρs[idx] * exp(-1.0 * g_0 * M_air * (h - hs[idx]) / (R * Ts[idx]))).to(g•cm⁻³)
  else: doAssert false, "Invalid height! Outside of range!"

import random
randomize(430)
proc intBetheAtmosphere(E: GeV, θ: Radian, dx = 10.cm): eV =
  ## integrated energy loss using Bethe formula for muons generated at
  ## `15.km` under an angle of `θ` to the observer for a muon of energy
  ## `E`.
  # Main contributions in Earth's atmosphere
  const τ = 2.19618.μs # muon half life
  let elements = Atmosphere
  var γ: UnitLess = E_to_γ(E.to(GeV))
  result = E.to(eV)
  var totalLoss = 0.eV
  let h_muon = 15.km # assume creation happens in `15.km`
  let S = h_muon.to(m) * distanceAtmosphere(θ.rad, d = h_muon)
  var S_prime = S
  while S_prime > 0.m and result > 0.eV:
    let h = computeHeight(S_prime, θ)
    let ρ_at_h = barometricFormula(h)
    var E_loss = 0.0.MeV
    for eTup in elements: # compute the weighted contribution of the element fraction
      let (w, e) = eTup
      E_loss += w * betheBloch(-1, e.Z, e.M, γ, m_μ) * ρ_at_h * dx

    ## Add step for radioactive decay of muon.
    ## - given `dx` compute likelihood of decay
    ## - eigen time of muon: dx / v = dt. dτ = dt / γ
    ## - muon decay is λ = 1 / 2.2e-6s
    let β = calcβ(γ)
    # compute effective time in lab frame
    let δt = dx / (β * c)
    # compute eigen time
    let δτ = δt / γ
    # probability of a decay in this time frame
    let p = pow(1 / math.E, δτ / τ)
    # decay with likelihood `p`
    #echo "γ = ", γ, " yields ", p, " in δτ ", δτ, " for energy ", E
    if rand(1.0) < (1.0 - p):
      echo "Particle decayed after: ", S_prime
      return 0.eV
          
    result = result - E_loss.to(eV)
    S_prime = S_prime - dx
    γ = E_to_γ(result.to(GeV))
    totalLoss = totalLoss + E_loss.to(eV)
  echo "total Loss ", totalLoss.to(GeV)
  result = max(0.float, result.float).eV

block MuonLimits:
  let τ_μ = 2.1969811.μs
  # naively this means given some distance `s` the muon can
  # traverse `s = c • τ_μ` (approximating its speed by `c`) before
  # it has decayed with a 1/e chance
  # due to special relativity this is extended by γ
  let s = c * τ_μ
  echo s
  # given production in 15 km, means
  let h = 15.km
  echo h / s
  # so a reduction of (1/e)^22. So 0.
  # now it's not 15 km but under an angle `θ = 88°`.
  let R_over_d = 174.UnitLess
  let n = 3.0
  let E₀ = 25.0.GeV
  let I₀ = 90.0.m⁻²•s⁻¹•sr⁻¹
  let E_c = 1.GeV
  let ε = 2000.GeV

  # distance atmospher gives S / d, where `d` corresponds to our `h` up there
  let S = h * distanceAtmosphere(88.0.degToRad.rad)
  # so about 203 km
  # so let's say 5 * mean distance is ok, means we ned
  let S_max = S / 5.0
  # so need a `γ` such that `s` is stretched to `S_max`
  let γ = S_max / s
  echo γ
  # ouch. Something has to be wrong. γ of 61?

  # corresponds to an energy loss of what?
  let Nitrogen = initElement(7.0.UnitLess, 14.006.g•mol⁻¹, 1.2506.g•dm⁻³.to(g•cm⁻³))
  echo "================================================================================"
  echo "Energy left: ", intBethe(Nitrogen, S.to(cm), 6.GeV.to(eV), dx = 1.m.to(μm)).to(GeV)
  proc print(E: GeV, θ: Radian) =
    let left = intBetheAtmosphere(E, θ = θ).to(GeV)
    echo "E = ", E, ", θ = ", θ, ", Bethe = ", E - left
  print(6.GeV, 0.Radian)
  #print(200.GeV, 0.Radian)  
  #print(200.GeV, 88.°.to(Radian))
  #print(200.GeV, 75.°.to(Radian))

  let E_loss75 = 100.GeV - intBetheAtmosphere(100.GeV, 75.°.to(Radian)).to(GeV)
  plotE_vs_flux(75.°.to(Radian),
                E_loss75, #23.78.GeV, #25.GeV, #E_loss75,
                1.0.GeV,
                90.m⁻²•s⁻¹•sr⁻¹, #65.2.m⁻²•s⁻¹•sr⁻¹,
                2000.GeV, # 854.GeV,
                "_at_75deg")

  
  echo "S@75° = ", h * distanceAtmosphere(75.0.degToRad.rad, d = 15.0.km)
  echo "================================================================================"  
echo E_to_γ(4.GeV)
echo E_to_γ(0.GeV)

proc plotE_vs_flux_and_angles(E_c: GeV, I₀: m⁻²•s⁻¹•sr⁻¹, ε: GeV,
                              suffix = "") =
  ## Generates a plot of the muon flux vs energy for a fixed set of different
  ## angles.
  ##
  ## The energy loss is computed using a fixed 
  let energies = logspace(log10(E_c.float), 2.float, 1000)
  let angles = linspace(0.0, 80.0, 9)
  block CalcLossEachMuon:
    var df = newDataFrame()
    for angle in angles:
      let E = energies.mapIt(it.GeV)
      let θ = angle.°.to(Radian)
      var flux = newSeq[float]()
      var E_initials = newSeq[float]()
      var E_lefts = newSeq[float]()
      var lastDropped = 0.GeV
      for e in E:
        let E_left = intBetheAtmosphere(e, θ).to(GeV)
        if E_left <= 0.0.GeV:
          echo "Skipping energy : ", e, " as muon was lost in atmosphere"
          continue
        elif E_left <= E_c:
          echo "Skipping energy : ", e, " as muon has less than E_c = ", E_c, " energy left"
          lastDropped = e
          continue
        let E₀ = e - E_left
        flux.add muonFlux(e, θ, E₀, E_c, I₀, ε).float
        E_initials.add e.float        
        E_lefts.add E_left.float
      let dfLoc = toDf({E_initials, E_lefts, flux, "angle [°]" : angle})
      #  .filter(f{`E_initials` >= lastDropped.float})
      df.add dfLoc

    var theme = themeLatex(fWidth = 0.5, width = 600, baseTheme = sideBySide)
    theme.tickWidth = some(theme.tickWidth.get / 2.0)
    theme.tickLength = some(theme.tickLength.get / 2.0)    
    theme.gridLineWidth = some(theme.gridLineWidth.get / 2.0)    
    ggplot(df, aes("E_initials", "flux", color = factor("angle [°]"))) +
      geom_line() +
      xlab(r"Initial energy [\si{GeV}]") + ylab(r"Flux [\si{m^{-2}.s^{-1}.sr^{-1}}]", margin = 2.5) +
      scale_x_log10() + scale_y_log10() +
      margin(right = 3.0) + 
      ggtitle(&"Muon flux at different angles{suffix}") +
      theme + 
      ggsave(&"/home/basti/phd/Figs/muons/initial_energy_vs_flux_and_angle_cosmic_muons{suffix}.pdf",
             useTeX = true, standalone = true, width = 600, height = 450)
  
    ggplot(df, aes("E_lefts", "flux", color = factor("angle [°]"))) +
      geom_line() +
      xlab(r"Energy at surface [\si{GeV}]") + ylab(r"Flux [\si{m^{-2}.s^{-1}.sr^{-1}}]", margin = 2.5) +
      scale_x_log10() + scale_y_log10() +
      margin(right = 3.0) +       
      theme + 
      ggtitle(&"Muon flux at different angles{suffix}") +
      ggsave(&"/home/basti/phd/Figs/muons/final_energy_vs_flux_and_angle_cosmic_muons{suffix}.pdf",
             useTeX = true, standalone = true, width = 600, height = 450)              
  block StaticLoss:
    var df = newDataFrame()
    for angle in angles:
      let E = energies.mapIt(it.GeV)
      let θ = angle.°.to(Radian)
      let E₀ = 100.GeV - intBetheAtmosphere(100.GeV, 0.0.Radian).to(GeV)    
      let flux = E.mapIt(muonFlux(it, θ, E₀, E_c, I₀, ε).float)
      let dfLoc = toDf({energies, flux, "angle [°]" : angle})
      df.add dfLoc
    ggplot(df, aes("energies", "flux", color = factor("angle [°]"))) +
      geom_line() +
      xlab("Energy [GeV]") + ylab("Flux [m⁻²•s⁻¹•sr⁻¹]") +
      scale_x_log10() + scale_y_log10() +
      ggtitle(&"Differential muon flux dependency at different angles{suffix}") +
      ggsave(&"/home/basti/phd/Figs/muons/energy_vs_flux_and_angle_cosmic_muons{suffix}.pdf")


#proc plotE_vs_flux_and_angles(E_c: GeV, I₀: m⁻²•s⁻¹•sr⁻¹, ε: GeV,
#                              suffix = "") =
#  ## Generates a plot of the integrated muon flux vs angles for a fixed set of different
#  ## energies.
#  let angles = linspace(0.0, 90.0, 100)
#  var df = newDataFrame()
#  let energies = linspace(E_c.float, 100.0, 1000)
#  let E = energies.mapIt(it.GeV)
#  for angle in angles:
#    let θ = angle.°.to(Radian)
#    let E₀ = 100.GeV - intBetheAtmosphere(100.GeV, θ).to(GeV)
#    let flux = E.mapIt(muonFlux(it, θ, E₀, E_c, I₀, ε).float) 
#    let dfLoc = toDf({energies, flux, "angle [°]" : angle})
#    df.add dfLoc
#  ggplot(df, aes("energies", "flux", color = factor("angle [°]"))) +
#    geom_line() +
#    xlab("Energy [GeV]") + ylab("Flux [m⁻²•s⁻¹•sr⁻¹]") +
#    scale_x_log10() + scale_y_log10() +
#    ggtitle(&"Differential muon flux dependency at different angles{suffix}") +
#    ggsave(&"/home/basti/phd/Figs/muons/energy_vs_flux_and_angle_cosmic_muons{suffix}.pdf")

# different angles!      
block MuonBehavior:
  plotE_vs_flux_and_angles(0.3.GeV, 90.m⁻²•s⁻¹•sr⁻¹, 854.GeV)

proc unbinnedCdf(x: seq[float]): (seq[float], seq[float]) =
  ## Computes the CDF of unbinned data
  var cdf = newSeq[float](x.len)
  for i in 0 ..< x.len:
    cdf[i] = i.float / x.len.float
  result = (x.sorted, cdf)

import random, algorithm
proc sampleFlux(samples = 1_000_000): DataFrame =
  randomize(1337)
  let energies = linspace(0.1, 100.0, 100_000)
  #let energies = logspace(0, 2, 1000)
  let E = energies.mapIt(it.GeV)
  let flux = E.mapIt(muonFlux(it, 88.0.degToRad.Radian, E₀, E_c, I₀, ε).float)
  # given flux compute CDF
  let fluxCS = flux.cumSum()
  let fluxCS_sorted = flux.sorted.cumSum()
  let fluxCDF = fluxCS.mapIt(it / fluxCS[^1])
  let fluxCDF_sorted = fluxCS_sorted.mapIt(it / fluxCS_sorted[^1])

  let (data, cdf) = unbinnedCdf(flux)

  let dfX = toDf(energies, fluxCS, fluxCS_sorted, fluxCDF, fluxCDF_sorted)
  ggplot(dfX, aes("energies", "fluxCS")) +
    geom_line() +
    ggsave("/t/cumsum_test.pdf")
  ggplot(dfX, aes("energies", "fluxCDF")) +
    geom_line() +
    ggsave("/t/cdf_test.pdf")    
  ggplot(dfX, aes("energies", "fluxCS_sorted")) +
    geom_line() +
    ggsave("/t/cumsum_sorted_test.pdf")    
  ggplot(dfX, aes("energies", "fluxCDF_sorted")) +
    geom_line() +
    ggsave("/t/cdf_sorted_test.pdf")

  ggplot(toDf(data, cdf), aes("data", "cdf")) +
    geom_line() +
    ggsave("/t/unbinned_cdf.pdf")
  
  #if true: quit()
  var lossesBB = newSeq[float]()
  var lossesMP = newSeq[float]()
  var energySamples = newSeq[float]()

  let dedxmin = 1.519.MeV•cm²•g⁻¹
  echo "Loss = ", (dedxmin * Argon.ρ * 3.cm).to(keV)
  
  for i in 0 ..< samples:
    # given the fluxCDF sample different energies, which correspond to the
    # distribution expected at CAST
    let idx = fluxCdf.lowerBound(rand(1.0))
    let E_element = E[idx]
    # given this energy `E` compute the loss
    let lossBB = (E_element.to(eV) - intBethe(Argon, 3.cm, E_element.to(eV), dx = 50.μm)).to(keV).float
    lossesBB.add lossBB
    let lossMP = mostProbableLoss(-1, Argon.Z, Argon.M, E_Element.E_to_γ(), Argon.ρ * 3.cm)
    lossesMP.add lossMP.float
    #echo "Index ", i, " yields energy ", E_element, " and loss ", loss
    energySamples.add E_element.float
  let df = toDf(energySamples, lossesBB, lossesMP)
    .gather(["lossesBB", "lossesMP"], "Type", "Value")
  ggplot(df, aes("Value", fill = "Type")) +
    geom_histogram(bins = 300, hdKind = hdOutline, alpha = 0.5, position = "identity") +
    margin(top = 2) +
    xlim(5, 15) +
    ggtitle(&"Energy loss of muon flux at CAST based on MC sampling with {samples} samples") +
    ggsave("/home/basti/phd/Figs/muons/sampled_energy_loss.pdf")

  ggplot(df, aes("energySamples")) +
    geom_histogram(bins = 300) +
    margin(top = 2) +
    ggtitle(&"Sampled energies for energy loss of muon flux at CAST") +
    ggsave("/home/basti/phd/Figs/muons/sampled_energy_for_energy_loss.pdf")
  let (samples, bins) = histogram(energySamples, bins = 300)
  let dfH = toDf({"bins" : bins[0 ..< ^1], samples})
    .filter(f{`bins` > 0.0 and `samples`.float > 0.0})
  ggplot(dfH, aes("bins", "samples")) +
    geom_line() +
    scale_x_log10() + 
    margin(top = 2) +
    ggtitle(&"Sampled energies for energy loss of muon flux at CAST") +
    ggsave("/home/basti/phd/Figs/muons/sampled_energy_for_energy_loss_manual.pdf")

  ggplot(toDf(energies, flux), aes("energies", "flux")) +
    geom_line() +
    scale_x_log10() +
    ggsave("/tmp/starting_data_e_flux.pdf")

  ggplot(toDf(energies, flux), aes("energies", "flux")) +
    geom_line() +
    ggsave("/tmp/linear_starting_data_e_flux.pdf")
    

discard sampleFlux(samples = 1_000_000)

ntangle thesis.org && nim c -d:release code/muons
code/muons

Of common interest when dealing with gas mixtures is the notion of partial pressures. In ideal gas theory, a mixture of gases at a pressure \(P\) can be considered to be the sum of the 'partial pressures' of each gas species

\[ P = \sum_i p_i, \]

initially noted by John Dalton in 1802 (Dalton 1802).

The contribution of each gas only depends on the species' mole fraction. Typically when considering gas mixtures the fractions of each gas species is given as a percentage. The percentage already refers to the mole fraction of each species. As such, the partial pressure of a single species can be expressed as:

\[ p_i = n_i P \]

where \(n_i\) is the mole fraction of species \(i\).

This is an extremely valuable property when computing different interactions of particles with gas mixtures. For example when computing the absorption of X-rays after propagating through a certain distance of a gas mixture, as one can compute absorption for each partial pressure independently and combine the contributions after (whether they be added, multiplied or something else of course depends on the considered process).

In order to understand the signals detected by a gaseous detector, the average number of electrons liberated by an incident particle should be known. This can be calculated given the energy loss in a gas required to produce a single electron-ion pair, called the \(W\text{-value}\),

\[ W = \frac{I}{⟨N⟩}. \]

Here, \(I\) is the mean ionization energy of the gas and \(⟨N⟩\) the average number of electron-ion pairs generated per interaction. \(⟨N⟩\) is usually smaller than one, \(\numrange{0.6}{0.7}\) for noble gases and even below \(\num{0.5}\) for molecular gases. Not all energy of the incoming particle is always deposited in the form of generation of, for example, a photoelectron. Other forms of energy loss are possible. In molecular gases vibrational and rotational modes offer even more possibilities, resulting in the smaller values.

The mean excitation energy \(I\) of an element is the weighted combination of the most likely energy levels the element is ionized from. The exact values are dependent on the specific element and tabulated values like from NIST (Hubbell and Seltzer 1996) exist. Above some \(Z\) (roughly argon, \(Z = 18\)) a rough approximation of \(I = \SI{10}{eV} · Z\) can be substituted (Zyla and others 2020), developed by Bloch (Bloch 1933).

The precise number for \(W\) strongly depends on the used gas mixture and typically requires experimental measurements to determine. Monte Carlo based simulations can be used as a rough guide, but care must be taken interpreting the results as significant uncertainty can remain. Tools for such Monte Carlo simulations include GEANT4 (Agostinelli and others 2003) ⁷ and MCNelectron (Poškus 2015) ⁸. These are based on atomic excitation cross sections, which are well tabulated in projects like ENDF (Brown et al. 2018) (citation for the latest data release) ⁹ and LXCat (Pancheshnyi et al. 2012; Pitchford et al. 2016; Carbone et al. 2021) ¹⁰.

While we have already discussed the mean free path for X-rays in sec. 6.1.1, we should revisit it quickly for electrons in a gas. It is a necessary concept when trying to understand other gaseous detector physics concepts, in particular for the gas amplification. ¹¹

The mean free path of an electron in a gas can be described by

where \(n = N/V\) is the number density (atoms or molecules per unit volume) of the gas particles and \(σ\) the cross section of electrons interacting with them. Such cross sections are tabulated for different elements, for example again see LXCat (Pancheshnyi et al. 2012; Pitchford et al. 2016; Carbone et al. 2021).

Based on the ideal gas law, we can express it as a function of pressure \(p\) and temperature \(T\) to

\[ p V = N k_B T ⇔ p = n k_B T, \]

with the Boltzmann constant \(k_B\). Inserting this via \(n\) into eq. \eqref{eq:theory:mean_free_path_e_def} yields

Diffusion in the context of gaseous detectors is the process of the random walk of electrons either in longitudinal or transverse (to the electric field) direction they exhibit, when drifting towards the readout. Depending on the specific detector technology, some transverse diffusion may be a desired property. For long drift distances, magnetic fields can be used to significantly reduce the transverse diffusion.

For an overview of the mathematics of diffusion as a random walk process, see (Berg 1993). In the context of gaseous detectors see for example any of (Sauli 2014; Kolanoski and Wermes 2020; Hilke and Riegler 2020) In general, for precise numbers measurements must be taken. Monte Carlo simulations using tools like Magboltz (Biagi 1995b) or PyBoltz (Atoum et al. 2020) can be used as a general reference if no measurements are available.

The mean distance \(σ_t\) from a starting point after a certain amount of time \(t\) is given by

\[ σ_t = \sqrt{ 2 N D t } \]

for a random walk in \(N\) dimensions, where \(D\) is a diffusion coefficient specifying the movement in distance squared per time. For example for a point like source of multiple electrons, after diffusion of some distance a 2 dimensional gaussian distribution will be expected with a standard deviation of \(σ_t\). By relating the time to the drift velocity \(v\) along an axis and introducing the diffusion constant \(D_t\)

\[ D_t = \sqrt{ \frac{2 N D}{v} } \]

we can further express \(σ_t\) as

which is often practically useful to estimate the expected diffusion after a certain drift length.

Note that the terminology of "diffusion constant", "diffusion coefficient" and similar is often used ambiguously as to whether they refer to \(D\) or \(D_t\) (or sometimes even \(σ_t\)). Nor is the considered dimensionality \(N\) always clearly indicated. Keeping \(N = 1\) and handling multiple dimensions as independent random walks is a practical approach to take (as long as it is valid in the application). ¹²

The drift velocity is the average speed at which electrons move towards the anode in a gaseous detector under the influence of an electric field. It is required to understand how recorded time information relates to longitudinal shape information of recorded data. Based on the so called 'friction force model' an analytical expression for the drift velocity in an electromagnetic field can be written to:

\[ \vec{v} = \frac{e}{m_e} \frac{τ}{1 + ω²τ²}\left( \vec{E} + \frac{ωτ}{B} (\vec{E} × \vec{B}) + \frac{ω²τ²}{B²}(\vec{E} · \vec{B}) \vec{B} \right) \]

with the electron charge \(e\) and mass \(m_e\), in an electric field \(\vec{E}\) and magnetic field \(\vec{B}\), given the Lamor frequency \(ω = eB / m_e\) and the mean collision time \(τ\). (Zyla and others 2020)

For detectors without a magnetic field \(B = 0, ω = 0\) and a constant, homogeneous electric field \(E\), this reduces to the Townsend expression:

\[ v = \frac{e E τ}{m_e}. \]

If measurements are not available, these can also be computed by Magboltz (Biagi 1995b) or PyBoltz (Atoum et al. 2020), which solve the underlying transport equation, the Boltzmann equation.

In order to turn the individual electrons into a measurable signal, gaseous detectors use some form of gas amplification. Details vary, but it is usually a region in the gas volume close to the readout with a very strong electric field (multiple \(\si{kV.cm^{-1}}\)) such that each electron causes many secondary ionizations, leading to an avalanche effect. In case of the detectors described in this thesis, amplification factors between \(\numrange{2000}{4500}\) are desired.

An electron in a strong electric field \(\vec{E}\) will gain enough energy to ionize an atom of the gas after a distance

\[ l \geq \frac{I}{|\vec{E}|} \]

with the ionization potential of the gas, \(I\). This needs to be put into relation with the mean free path, eq. \eqref{eq:theory:mean_free_path_e}. If \(l \ll λ\) no secondary ionization takes place and if \(l \gg λ\) every interaction leads to ionization, likely resulting in a breakdown causing an arc (see also Paschen's law (Paschen 1889), not covered here). In the intermediate range some interactions cause secondary ionization, some does not.

We can make the statistical argument that

\[ \mathcal{N} = e^{-l/λ} \]

is the relative number of collisions with \(l > λ\). This allows to define the probability of finding a number of ionizations per unit length to be

\[ P(l) = \frac{1}{λ} e^{-l / λ} = α, \]

where we introduce the definition of the 'first Townsend coefficient', \(α\). We can insert the definition of the mean free path, eq. \eqref{eq:theory:mean_free_path_e}, and \(l\) into this equation to obtain

With this we have an expression for the temperature and pressure dependency of the first Townsend coefficient. ¹³ This derivation followed (Scharenberg 2019) based on (Engel and Marton 1965), but see (Aoyama 1985) for a more general treatment. Similarly to diffusion and drift parameters, the first Townsend coefficient can be computed numerically using tools like Magboltz.

Note that (Sauli 2014) introduces the Townsend coefficient as \(α = \frac{1}{λ}\), introducing a specific \(λ\) and \(σ\) referring to the ionization mean free path and ionizing cross sections. This can be misleading as it makes it seem as if the Townsend coefficient is inversely proportional to the temperature. This is of course only the case in the regime where each interaction actually causes another ionization.

The first Townsend coefficient can be used to express the multiplication factor – or gas amplification – used in a detector,

\[ G = e^{α x} \]

as the number increases by \(α\) after each step \(\mathrm{d}x\).

The statistical distribution describing the number of electrons after gas amplification is the Pólya distribution

where \(N\) is a normalization constant, \(θ\) is another parameter performing shaping of the distribution and \(G\) is the effective gas gain. \(Γ\) refers to the gamma function. It is to note that the term "polya distribution" in this context is different from other mathematical definitions, in which polya distribution usually means a negative binomial distribution. The above definition goes back to Alkhazov (Alkhazov 1970) and in slight variations is commonly used. Due to the complexity of this definition, care needs to be taken when performing numerical fits to data with this function (using bounded parameters and preferring more general non-linear optimization routines instead of a normal Levenberg-Marquardt non-linear least squares approach).

Based on eq. \eqref{eq:theory:townsend_coefficient} the largest impacts on the expected gas amplification have the electric field, the choice of gas and the temperature of the gas. While the former two parameters are comparatively easy to control, the temperature in the amplification region may vary and is hard to measure. As such depending on the detector details and application, gas gain variations are expected and corrections based on a running gas gain value may be necessary.

As the large number of interactions in the amplification region can excite many atoms of the (typically) noble gas, UV photons can be produced. Their mean free path is comparatively long relative to the size of the amplification region. They can start further avalanches, potentially away from the location of the initial avalanche start, lowering spatial resolution and increasing apparent primary electron counts. Molecular gases are added – often only in small fractions – to the gas mixture to provide rotational or vibrational modes to dissipate energy without emitting UV photons. In this context the molecular additions are called 'quencher gases'.

Because of the statistical processes associated with the full detection method used in gaseous detectors, even a perfect delta like signal in energy, will lead to a natural spread of the measured signal. The convolution of different ionization yields, potential losses and gas gain variation all contribute to such a spread.

As such a typical measure of interest for gaseous detectors is the energy resolution, which is commonly defined by

\[ ΔE = \frac{σ}{E} \]

where \(σ\) is the standard deviation of the normal distribution associated with the resolved peak in the detectors data, assuming a – for practical purposes – delta peak as the input spectrum. Sometimes definitions using the full width at half maximum (FWHM) are also used in place of \(σ\). Typical values for the energy resolution defined like this are smaller than \(\SI{15}{\%}\).

If the absolute magnitude of the \(σ\) at a given energy is constant, which at least is partially reasonable as the width is not fully due to energy dependent effects, the energy resolution is proportional to \(1/E\).

The Fano factor (Fano 1963), defined by the variance over the mean of a distribution \(F = σ² / μ\) (typically within some time window), improves the ideal energy resolution. It arises in the associated statistical processes, because there are a finite number of possible interaction states. For X-rays an additional aspect is due to the maximum energy transferable into the gas due to the X-rays energy. In practice though the energy resolution of gaseous detectors is usually limited by other effects.

Finally, gaseous detectors need to be calibrated, monitored and tested. This is commonly done with a \cefe source. \cefe is a radioactive isotope of iron, which decays under inverse beta decay to \(\ce{^{55}Mn}\). Due to the required restructuring of the electronic shells, the manganese is in an excited state. While the emission of an Auger electron with \(\SI{5.19}{keV}\) dominates with a probability of \(\SI{67.2}{\%}\) (Krause 1979; Schoonjans et al. 2011; Brunetti et al. 2004), as an X-ray source the Kα₁ and Kα₂ lines with combined energies of about \(\SI{5.9}{keV}\) are of note.

When using such a \cefe source as a calibration source for an argon filled gaseous detector, the \(\SI{5.9}{keV}\) photons will produce a photoelectron in argon. Mostly an inner shell electron will be liberated, producing a photoelectron of around \(\SI{2.7}{keV}\). The excited argon atom can now emit an Auger electron in about \(\SI{88.2}{\%}\) of cases (Krause 1979; Schoonjans et al. 2011; Brunetti et al. 2004), which will fully deposit its remaining energy into the surrounding gas via further ionizations, resulting in the 'photopeak' at around \(\SI{5.9}{keV}\).

If however another photon is produced with an energy below the \(K 1s\) energy of argon (\(\SI{3.2}{keV}\)) – for example via Kα₁ or Kα₂ fluorescence of argon, both at about \(\SI{2.95}{keV}\) – such a photon has a very long absorption length in the gas volume, about \(l_{\text{abs}} = \SI{3.5}{cm}\) (cf. fig. 1). This can cause it to easily escape the active detection region, especially if the sensitive region of the detector is comparatively small. The result is a measured signal of \(E_i - E_k = \SI{5.9}{keV} - \SI{2.95}{keV} \approx \SI{2.9}{keV}\), called the 'escape peak'. (Sauli 2014; Kolanoski and Wermes 2020; Hubbell and Seltzer 1996)

Such an additional escape peak is useful as a calibration tool, as it gives two distinct peaks in a \cefe calibration spectrum, which can be utilized for an energy calibration. One important consideration of such an escape peak is however, that this escape peak of about \(\SI{3}{keV}\) is not equivalent to a real \(\SI{3}{keV}\) X-ray. As the original X-ray is still the \(\SI{5.9}{keV}\) photon with its distinct absorption length, the geometric properties of ensembles of escape peak events and real \(\SI{3}{keV}\) photons is different. This will be important later in sec. 12.4.7.

The fluorescence yield is the fraction with which an excited nucleus decays via a radiative emission of an X-ray instead of an Auger electron.

So the relevant terms are either 'fluorescene yield' or 'Auger yield' depending.

The main resource still often compared to is this paper by Krause (Krause 1979). The paper contains a very useful table of the fluorescence yields of the different elements.

There is a C library for such calculations, (Schoonjans et al. 2011; Brunetti et al. 2004). They have produced a PDF with similar tables for different elements, which can be found here:

http://ftp.esrf.fr/pub/scisoft/xraylib/xraylib_tables_v2.3.pdf

One might encounter the term 'jump factor' in this context: It refers to the ratio of how the absorption changes from directly below to directly above the absorption edge.

So for Argon see page 31 in this PDF (or page 9 in Krause's paper) to find \(ω_K = \num{0.118}\), where

\[ ω_K + a_K = 1 \]

the fluorescence yield \(ω_K\) + auger yield \(a_K\) needs to be one. NOTE: Technically \(ω_K\), \(a_K\) is only for the K shell, so just an approximation. You will find \(ω_1, ω_2\), … etc for L shells.

For Manganese (because Manganese is the excited nucleus after \cefe decay) you'll find it on page 47 (or page 9 in Krause's paper). They sum to \(ω_{K+L} = 0.345\).