As mentioned in the note above, the original document is titled SolarAxionConversionPoint.org. From the Org file a PDF and HTML version is produced.

These can be found on

http://phd.vindaar.de/docs/SolarAxionConversionPoint/

That document contains an attempt to compute the same thing analytically, which as I later realized produces wrong results of course (at least if done the way I did it there). See the document.

In order to calculate the conversion point, we need:

• random sampling logic
• sampling from exponential distribution depending on energy
• the axion flux, telescope effective area and window absorption

Let's start by importing the modules we need:

import helpers / sampling_helper # sampling distributions
import unchained                 # sane units
import ggplotnim                 # see something!
import xrayAttenuation           # window efficiencies
import math, sequtils

where the sampling_helpers is a small module to sample from a procedure or a sequence.

In addition let's define some helpers:

from os import `/`, expandTilde
const ResourcePath = "~/org/resources".expandTilde
const OutputPath = "~/phd/Figs/axion_conversion_point_sampling/".expandTilde

proc thm(): Theme =
  ## A shorthand to define a `ggplotnim` theme that looks nice 
  ## in the thesis
  result = themeLatex(fWidth = 0.9, width = 600, baseTheme = singlePlot)

Now let's read the LLNL telescope efficiency as well as the axion flux model. Note that we may wish to calculate the absorption points not only for a specific axion flux model, but potentially any other kind of signal. We'll build in functionality to disable different contributions.

let flux = "solar_axion_flux_differential_g_ae_1e-13_g_ag_1e-12_g_aN_1e-15.csv"
let dfAx = readCsv(ResourcePath / flux)
  .filter(f{`type` == "Total flux"})
let llnl = "llnl_xray_telescope_cast_effective_area_parallel_light_DTU_thesis.csv"
let dfLLNL = readCsv(ResourcePath / llnl)
  .mutate(f{"Efficiency" ~ idx("EffectiveArea[cm²]") / (PI * 2.15 * 2.15)})

Note: to get the differential axion flux use readOpacityFile from https://github.com/jovoy/AxionElectronLimit. It generates the CSV file.

Next up we need to define the material properties of the detector window in order to compute its transmission.

let Si₃N₄ = compound((Si, 3), (N, 4)) # actual window
const ρSiN = 3.44.g•cm⁻³
const lSiN = 300.nm                  # window thickness
let Al = Aluminium.init()            # aluminium coating
const ρAl = 2.7.g•cm⁻³
const lAl = 20.nm                    # coating thickness

With these numbers we can compute the transmission at an arbitrary energy. In order to compute the correct inputs for the calculation we now have everything. We wish to compute the following, the intensity \(I(E)\) is the flux that enters the detector

\[ I(E) = f(E) · ε_{\text{LLNL}} · ε_{\ce{Si3.N4}} · ε_{\ce{Al}} \]

where \(f(E)\) is the solar axion flux and the \(ε_i\) are the efficiencies associated with the telescope and transmission of the window. The idea is to sample from this intensity distribution to get a realistic set of X-rays as they would be experienced in the experiment. One technical aspect still to be done is an interpolation of the axion flux and LLNL telescope efficiency to evaluate the data at an arbitrary energy as to define a function that yields \(I(E)\).

Important note: We fully neglect here the conversion probability and area of the magnet bore. These (as well as a potential time component) are purely constants and do not affect the shape of the distribution \(I(E)\). We want to sample from it to get the correct weighting of the different energies, but do not care about absolute numbers. So differential fluxes are fine.

The idea is to define the interpolators and then create a procedure that captures the previously defined properties and interpolators.

from numericalnim import newLinear1D, eval
let axInterp = newLinear1D(dfAx["Energy", float].toSeq1D,
                           dfAx["diffFlux", float].toSeq1D)
let llnlInterp = newLinear1D(dfLLNL["Energy[keV]", float].toSeq1D,
                             dfLLNL["Efficiency", float].toSeq1D)

With the interpolators defined let's write the implementation for \(I(E)\):

proc I(E: keV): float =
  ## Compute the intensity of the axion flux after telescope & window eff.
  ##
  ## Axion flux and LLNL efficiency can be disabled by compiling with
  ## `-d:noAxionFlux` and `-d:noLLNL`, respectively.
  result = transmission(Si₃N₄, ρSiN, lSiN, E) * transmission(Al, ρAl, lAl, E)
  when not defined(noAxionFlux):
    result *= axInterp.eval(E.float)
  when not defined(noLLNL):
    result *= llnlInterp.eval(E.float)
  

Let's test it and see what we get for e.g. \(\SI{1}{keV}\):

echo I(1.keV)

yields \(1.249e20\). Not the most insightful, but it seems to work. Let's plot it:

let energies = linspace(0.01, 10.0, 1000).mapIt(it.keV)
let Is = energies.mapIt(I(it))
block PlotI:
  let df = toDf({ "E [keV]" : energies.mapIt(it.float),
                  "I" : Is })
  ggplot(df, aes("E [keV]", "I")) +
    geom_line() +
    ggtitle("Intensity entering the detector gas") +
    margin(left = 3.0) + thm() + 
    ggsave(OutputPath / "intensity_axion_conversion_point_simulation.pdf")

shown in fig. 1. It looks exactly as we would expect.

Now we define the sampler for the intensity distribution \(I(E)\), which returns an energy weighted by \(I(E)\):

let Isampler = sampler(
  (proc(x: float): float = I(x.keV)), # wrap `I(E)` to take `float`
  0.01, 10.0, num = 1000 # use 1000 points for EDF & sample in 0.01 to 10 keV
)

and define a random number generator:

import random
var rnd = initRand(0x42)

First we will sample 100,000 energies from the distribution to see if we recover the intensity plot from before.

block ISampled:
  const nmc = 100_000
  let df = toDf( {"E [keV]" : toSeq(0 ..< nmc).mapIt(rnd.sample(Isampler)) })
  ggplot(df, aes("E [keV]")) +
    geom_histogram(bins = 200, hdKind = hdOutline) +
    ggtitle("Energies sampled from I(E)") +
    thm() + 
    ggsave(OutputPath / "energies_intensity_sampled.pdf")

This yields fig. 2, which clearly shows the sampling works as intended.

The final piece now is to use the same sampling logic to generate energies according to \(I(E)\), which correspond to X-rays of said energy entering the detector. For each of these energies then sample from the Beer-Lambert law (see sec. 6.1.1)

\[ I(z) = I_0 \exp\left[ - \frac{z}{l_{\text{abs}} } \right], \] where \(I_0\) is some initial intensity and \(l_\text{abs}\) the absorption length. The absorption length is computed from the gas mixture properties of the gas used at CAST, namely argon/isobutane 97.7/2.3 at \(\SI{1050}{mbar}\). It is the inverse of the attenuation coefficient \(μ_M\)

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

where the attenuation coefficient is computed via

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

with \(N_A\) Avogadro's constant, \(M\) the molar mass of the compound and \(σ_A\) the atomic absorption cross section. The latter again is defined by

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

with \(r_e\) the classical electron radius, \(λ\) the wavelength of the X-ray and \(f₂\) the second scattering factor. Scattering factors are tabulated for different elements, for example by NIST and Henke. For a further discussion of this see the README and implementation of xrayAttenuation (Schmidt and SciNim contributors 2022).

We will now go ahead and define the CAST gas mixture:

proc initCASTGasMixture(): GasMixture =
  ## Returns the absorption length for the given energy in keV for CAST
  ## gas conditions:
  ## - Argon / Isobutane 97.7 / 2.3 %
  ## - 20°C ( for this difference in temperature barely matters)
  let arC = compound((Ar, 1)) # need Argon gas as a Compound
  let isobutane = compound((C, 4), (H, 10))
  # define the gas mixture
  result = initGasMixture(293.K, 1050.mbar, [(arC, 0.977), (isobutane, 0.023)])
let gm = initCASTGasMixture()  

To sample from the Beer-Lambert law with a given absorption length we also define a helper that returns a sampler for the target energy using the definition of a normalized exponential distribution

\[ f_e(x, λ) = \frac{1}{λ} \exp \left[ -\frac{x}{λ} \right]. \]

Note: The sampling of the conversion point is the crucial aspect of this. Naively we might want to sample between the detector volume from 0 to \(\SI{3}{cm}\). However, this skews our result. Our calculation depends on the energy distribution of the incoming X-rays. If the absorption length is long enough the probability of reaching the readout plane and thus not being detected is significant. Restricting the sampler to \(\SI{3}{cm}\) would pretend that independent of absorption length we would always convert within the volume, giving too large a weight to the energies that should sometimes not be detected!

Let's define the sampler now. It takes the gas mixture and the target energy. A constant SampleTo is defined to adjust the position to which we sample at compile time (to play around with different numbers).

proc generateSampler(gm: GasMixture, targetEnergy: keV): Sampler =
  ## Generate the exponential distribution to sample from based on the
  ## given absorption length
  # `xrayAttenuation` `absorptionLength` returns number in meter!
  let λ = absorptionLength(gm, targetEnergy).to(cm)
  let fnSample = (proc(x: float): float =
                    result = expFn(x, λ.float) # expFn = 1/λ · exp(-x/λ)
  )
  const SampleTo {.intdefine.} = 20 ## `SampleTo`, set via `-d:SampleTo=<int>`
  let num = (SampleTo.float / 3.0 * 1000).round.int # # of points to sample at
  result = sampler(fnSample, 0.0, SampleTo, num = num)

Note that this is inefficient, because we generate a new sampler from which we only sample a single point, namely the conversion point of that X-ray. If one intended to perform a more complex calculation or wanted to sample orders of magnitude more X-rays, one should either restructure the code (i.e. sample from known energies and then reorder based on the weight defined by \(I(E)\)) or cache the samplers and pre-bin the energies.

For reference let's compute the absorption length as a function of energy for the CAST gas mixture:

block GasAbs:
  let Es = linspace(0.03, 10.0, 1000)
  let lAbs = Es.mapIt(absorptionLength(gm, it.keV).m.to(cm).float)
  let df = toDf({ "E [keV]" : Es,
                  "l_abs [cm]" : lAbs })
  ggplot(df, aes("E [keV]", "l_abs [cm]")) +
    geom_line() +
    ggtitle(r"Absorption length of X-rays in CAST gas mixture: \\" & $gm) +
    margin(top = 1.5) +
    thm() + 
    ggsave(OutputPath / "cast_gas_absorption_length.pdf")

which yields fig. 3

So, finally: let's write the MC sampling!

const nmc = 500_000 # start with 100k samples
var Es = newSeqOfCap[keV](nmc)
var zs = newSeqOfCap[cm](nmc)
while zs.len < nmc:
  # 1. sample an energy according to `I(E)`
  let E = rnd.sample(Isampler).keV
  # 2. get the sampler for this energy
  let distSampler = generateSampler(gm, E) 
  # 3. sample from it
  var z = Inf.cm
  when defined(Equiv3cmSampling): ## To get the same result as directly sampling
                                  ## only up to 3 cm use the following code
    while z > 3.0.cm:
      z = rnd.sample(distSampler).cm 
  elif defined(UnboundedVolume): ## This branch pretends the detection volume
                                 ## is unbounded if we sample within 20cm
    z = rnd.sample(distSampler).cm 
  else: ## This branch is the physically correct one. If an X-ray reaches the
        ## readout plane it is _not_ recorded, but it was still part of the
        ## incoming flux!
    z = rnd.sample(distSampler).cm
    if z > 3.0.cm: continue # just drop this X-ray
  zs.add z
  Es.add E

Great, now we have sampled the conversion points according to the correct intensity. We can now ask for statistics or create different plots (e.g. conversion point by energies etc.).

import stats, seqmath # mean, variance and percentile
let zsF = zs.mapIt(it.float) # for math
echo "Mean conversion position = ", zsF.mean().cm
echo "Median conversion position = ", zsF.percentile(50).cm
echo "Variance of conversion position = ", zsF.variance().cm

This prints the following:

Mean conversion position = 0.556813 cm
Median conversion position = 0.292802 cm
Variance of conversion position = 0.424726 cm

We see the mean conversion position is at about \(\SI{0.56}{cm}\). If we consider the median it's only \(\SI{0.29}{cm}\) (the number we use in 37.4). This number provides the target for the raytracing of the axion image as an offset from the focal length in sec. 37.4.

Let's plot the conversion points of all sampled (and recorded!) X-rays as well as what their distribution against energy looks like.

let dfZ = toDf({ "E [keV]" : Es.mapIt(it.float),
                 "z [cm]"  : zs.mapIt(it.float) })
ggplot(dfZ, aes("z [cm]")) +
  geom_histogram(bins = 200, hdKind = hdOutline) +
  ggtitle("Conversion points of all sampled X-rays according to I(E)") +
  thm() + 
  ggsave(OutputPath / "sampled_axion_conversion_points.pdf")
ggplot(dfZ, aes("E [keV]", "z [cm]")) +
  geom_point(size = 0.5, alpha = 0.2) + 
  ggtitle("Conversion points of all sampled X-rays according to I(E) " &
    "against their energy") +
  thm() + 
  ggsave(OutputPath / "sampled_axion_conversion_points_vs_energy.pdf",
         dataAsBitmap = true)

The former is shown in fig. 4. The overlapping exponential distribution is obvious, as one would expect. The same data is shown in fig. 5, but in this case not as a histogram, but by their energy as a scatter plot. We can clearly see the impact of the absorption length on the conversion points for each energy!

The code above is written in literate programming style. To compile and run it we use ntangle to extract it from the Org file:

ntangle <file>

which generates ./code/sample_axion_xrays_conversion_points.nim.

Compiling and running it can be done via:

nim r -d:danger code/sample_axion_xrays_conversion_points.nim

which compiles and runs it as an optimized build.

We have the following compilation flags to compute different cases:

• -d:noLLNL: do not include the LLNL efficiency into the input intensity
• -d:noAxionFlux: do not include the axion flux into the input intensity
• -d:SampleTo=<int>: change to where we sample the position (only to 3cm for example)
• -d:UnboundedVolume: if used together with the default SampleTo (or any large value) will effectively compute the case of an unbounded detection volume (i.e. every X-ray recorded with 100% certainty).
• -d:Equiv3cmSampling: Running this with the default SampleTo (or any large value) will effectively change the sampling to a maximum \SI{3}{cm} sampling. This can be used as a good crossheck to verify the sampling behavior is independent of the sampling range.

Configurations of note:

nim r -d:danger -d:noAxionFlux code/sample_axion_xrays_conversion_points.nim

\(⇒\) realistic case for a flat input spectrum Yields:

Mean conversion position = 0.712102 cm
Median conversion position = 0.445233 cm
Variance of conversion position = 0.528094 cm

nim r -d:danger -d:noAxionFlux -d:UnboundedVolume code/sample_axion_xrays_conversion_points.nim

\(⇒\) the closest analogue to the analytical calculation from section [BROKEN LINK: sec:axion_conversion_point:analytical] (outside of including isobutane here) Yields:

Mean conversion position = 1.25789 cm
Median conversion position = 0.560379 cm
Variance of conversion position = 3.63818 cm

nim r -d:danger code/sample_axion_xrays_conversion_points.nim

\(⇒\) the case we most care about and of which the numbers are mentioned in the text above.