4. Theory of axions    Theory

A deep understanding of the axion requires a diverse set of knowledge of different aspects of the theoretical foundations of modern physics. The Standard Model of particle physics (SM), quantum field theory (QFT) and quantum chromodynamics (QCD). The \cpt invariance of the SM and related \cp violation 1 of the weak force. The common algebraic structure of the weak and strong force via \(\mathrm{SU}(2)\) and \(\mathrm{SU}(3)\). Anomalies in QFT, in particular the Adler-Bell-Jackiw anomaly. The structure of the QCD vacuum and the related \(\mathrm{U}(1)\) problem, which is solved by instantons and the concepts of Goldstone's theorem. To do justice to all these aspects in the context of a non-theory PhD thesis is not possible 2. As such, the theory part of this thesis will be kept short and instead a focus is placed on referencing useful material for an interested reader, in particular see sec. 4.1.

At low energies the Standard Model can be described by a combination of three different forces, the electromagnetic, the weak and the strong force. These can be represented mathematically by an internal group structure of \(\mathrm{U}(1) \times \mathrm{SU}(2) \times \mathrm{SU}(3)\), respectively. 3 The weak force, represented by \(\mathrm{SU}(2)\), has long since been known to exhibit a \cp violation 1. Due to the similar group structure between the weak and the strong force (\(\mathrm{SU}(2)\) vs. \(\mathrm{SU}(3)\)) many parallels exist between the mathematical descriptions of the two forces in the Standard Model. In particular, the following term, 4

\[ \mathcal{L}_θ = θ \frac{g_s²}{32π²} G^{μν}_a \tilde{G}_{aμν}, \]

is allowed under the requirements for a Standard Model conform term (i.e. gauge invariant, Lorentz invariant and so on). Here, \(θ\) is an angle, \(g²_s\) is the strong coupling constant and \(G^{μν}_α\) the gluon field strength tensor. As a matter of fact, this term is even required if the instanton solution to the \(\mathrm{U}(1)\) problem is considered (Hooft 1986; Peccei 2008b) and is a result of the complex vacuum structure of QCD (’T Hooft 1976). This term violates \(\mathrm{P}\) and \(\mathrm{T}\) transformations and as a result of \cpt symmetry also violates \cp. Peculiarly, any effect expected from this \cp violation has still not been observed. One such effect is an expected electric dipole moment of the neutron (Crewther et al. 1979, 1980; Baluni 1979). Such a dipole moment may naively be expected plainly from the fact that the constituent quarks of a neutron are charged after all. However, very stringent limits place an extremely low upper bound on it at (Baker et al. 2006; Pendlebury et al. 2015)

\[ d_{\text{NEDM}} \leq \SI{3e-26}{\elementarycharge \cm}, \]

where \(e\) is the electron charge. Nature's deviation from our expectation in this context is coined the strong \cp problem of particle physics. One possible solution to the strong \cp problem would be a massless up or down quark, in which case the QCD Lagrangian would feature a global \(\mathrm{U}(1)\) axial shift symmetry, which could shift \(θ ↦ 0\). Even in the late 1980's this was not entirely ruled out, despite the already understood mass ratio \(m_u / m_d = 5/9\) (Weinberg 1977), due to \(2^{\text{nd}}\) order chiral effects (Kaplan and Manohar 1986). However, nowadays based on lattice QCD calculations this has been ruled out (Cohen, Kaplan, and Nelson 1999; Dine, Draper, and Festuccia 2015; Group et al. 2022).

While it is possible that our universe is simply one in which the effect of the strong \cp violation is suppressed (or even exactly zero) "by chance", Helen Quinn and Roberto Peccei realized in 1977 (Peccei and Quinn 1977a, 1977b) that this behavior can be explained in the presence of an additional scalar field. Shortly after both Weinberg and Wilczek (Weinberg 1978; Wilczek 1978) realized the implication of such an additional field, namely a pseudo Nambu-Goldstone boson. Wilczek named it the axion, after a washing detergent, as it "washes the Standard Model clean" of the strong \cp problem. The most straightforward axion model based on the work by Wilczek and Weinberg yields a coupling of the axion to matter that is already excluded, because they associate the spontaneous symmetry breaking with the electroweak scale.

Models for an 'invisible axion' manage to unify the solution to the strong \cp problem with the current lack of experimental evidence for an axion-like particle. There are two main models for the invisible axion, the KSVZ (Kim-Shifman-Vainshtein-Zakharov) (Kim 1979; Shifman, Vainshtein, and Zakharov 1980) and the DFSZ (Dine-Fischler-Srednicki-Zhitnitskii) (Dine, Fischler, and Srednicki 1981; Zhitnitskii 1980) models.

For the most comprehensive modern overview of the theory of axions, different models, best bounds on different axion couplings and a general axion reference, make sure to look into the aptly named "landscape of QCD axion models" (Luzio, Giannotti, et al. 2020)!

From here we will first say a few more words about the two invisible axion models in sec. 4.3. Then we will look at the implications an effective axion Lagrangian has on axion interactions. We will sketch the derivation of the axion-photon conversion probability in sec. 4.4. This conversion probability leads to a discussion of the expected solar axion fluxes in sec. 4.5. At this point a very succinct interlude will introduce the chameleon particle, sec. 4.6. Finally, we will briefly go over the current relevant bounds on different axion couplings in sec. 4.7.

4.1. Useful reading material    optional

The following is a list of materials I would recommend, if someone wishes to understand the theory underlying axions better. The focus here is on things not already mentioned in the rest of the text. Your mileage may vary, of course.

First of all, at least have a short look at the original papers by Roberto Peccei and Helen Quinn (Peccei and Quinn 1977a, 1977b), as well as the 'responses' by Wilczek and Weinberg (Wilczek 1978; Weinberg 1978).

If your interest is in understanding the reasoning behind the origin of the strong CP problem via the QCD vacuum structure and the \(\mathrm{U}(1)\) problem, it may be worthwhile to look at (some of) the original papers. See Weinberg's paper introducing it (Weinberg 1975) and 't Hooft's paper on symmetry breaking via the Adler-Bell-Jackiw anomaly (’T Hooft 1976) as the next logical step. In (Hooft 1986) 't Hooft then shows how instantons solve the \(\mathrm{U}(1)\) problem. In 1999 't Hooft also wrote a 'historical review' (’T Hooft 1999) about renormalization of gauge theories, which touches on the \(\mathrm{U}(1)\) problem as well. It may be easier to follow (as these things so often are in hindsight).

In the vein of review papers, Roberto Peccei wrote a nice review paper in 2008 (Peccei 2008a) covering the connection between the \(\mathrm{U}(1)\) problem, the QCD vacuum structure, the strong CP problem and axions. Similarly, Jihn Kim (the K in the KSVZ model) wrote another very nice review in 2010 (Kim and Carosi 2010). I would probably recommend to read these two for a good overview. They give enough references to dive deeper, if desired.

On the side of axion detection, you should of course look into the two papers by Pierre Sikivie, introducing all our modern axion detection experiment techniques (Sikivie 1983, 1985). Sikivie wrote another review in 2020 (Sikivie 2021) covering search methods for invisible axions.

Another great overview, covering both theory and experimental parts of the axion, are Igor Irastorza's notes here: (Irastorza 2022). And there is also the text book 'Axions' (Kuster, Raffelt, and Beltrán 2007), which may be worth a look as a summary of different aspects of axion physics.

If you only want to look into a single document covering all aspects of axions, from theory to experimental searches, look no further than the 'landscape of QCD axion models' (Luzio, Giannotti, et al. 2020), also mentioned in the main text.

Finally, in my master thesis (Schmidt 2016) I attempted to introduce the relevant physics for the axion (and chameleon) to a level satisfactory to me at the time. Whether I succeeded or not, you can be the judge. This may (or may not) be a good introduction if you are a student, for whom courses in QFT are fresh in their mind.

4.2. Illustration of the \cpt symmetry [/]    extended

Fig. 1 is an illustration of the three discrete transformations part of the \cpt symmetry in a 1 dimensional spacetime. The three transformations are:

  • \(C\) - charge conjugation: replaces each particle by its anti particle and thus reversing the charge \(q ↦ -q\)
  • \(P\) - parity transformation: mirrors all positions at some origin, \(\vec{x} ↦ -\vec{x}\).
  • \(T\) - time reversal: reverses the arrow of time, \(t ↦ -t\)

The \cpt symmetry of the Standard Model states that a hypothetical mirror universe to ours - obtained by the three transformations together - follows the exact same physical laws and thus 'evolves' identically (due to the time reversal 'evolution' in quotation marks).

cpt_explanation_extended_linear.svg
Figure 1: Figure 1: Schematic showcasing the three transformations \(C, P\) and \(T\) in a 1 dimensional spacetime. Each of the three are discrete transformations essentially reversing the value along its axis. The combination of all three operations yields the schematic on the right.

4.3. Invisible axion models and axion couplings

Kim-Shifman-Vainshtein-Zakharov model
The so called KSVZ model (Kim 1979; Shifman, Vainshtein, and Zakharov 1980) is the simplest invisible axion model. It adds a scalar field \(σ\) and a superheavy quark \(Q\), which \(σ\) couples to via a Yukawa coupling. The main problem with the standard axion is that its energy scale is the electroweak scale \(v_F \approx \SI{250}{GeV}\) resulting in too strong interactions. The KSVZ model effectively achieves a symmetry breaking scale \(f_a \gg v_F\) and thus results in an 'invisible axion'. It contains a tree-level axion-photon coupling \(g_{aγ}\), but no axion-electron couplings. The latter can be found at one-loop level (Srednicki 1985).
Dine-Fischler-Srednicki-Zhitnitskii model
The DFZS axion model (Dine, Fischler, and Srednicki 1981; Zhitnitskii 1980) is another axion model, in which the scalar field \(σ\) couples to two Higgs doublet fields, \(H_u\) and \(H_d\). It does not require an extra superheavy quark, in this case the coupling of the scalar to the doublets achieves the decoupling of the axion symmetry breaking scale \(f_a\) from the electroweak scale \(v_F\). The end result is the same, an 'invisible axion' which is very light and has only light interactions with other matter. However, in contrast to the KSVZ model axion-lepton couplings appear at tree level!
Generalizations

Further generalizations of axion models over the KSVZ and DFSZ models are possible, resulting in more flexible couplings. From a practical standpoint of an axion experiment it is usually better to consider axion interactions via effective couplings, as the limiting factor is detecting something rather than determining its properties. In particular because the two models above yield very small expected coupling constants in regions of axion masses easily accessible via laboratory experiments. Therefore, an effective Lagrangian like the following

\begin{equation} \label{eq:theory:axion:general_axion_couplings} \mathcal{L}_{a,\text{eff}} = \frac{1}{2} \partial^{μ} a \partial_{μ} a - \frac{1}{2} m_a^2 a^2 - \frac{g_{aγ}}{4} \widetilde{F}^{μν} F_{μν} a - g_{ae} \frac{\partial_{μ} a}{2 m_e} \overline{ψ}_e γ^5 γ^{μ} ψ_e, \end{equation}

which contains an axion-photon coupling \(g_{aγ}\) and an axion-electron coupling \(g_{ae}\) is useful for experimental searches. Limits given on one of these parameters can always be converted to the specific couplings of one of the existing models if needed.

In principle other couplings exist, for example the axion-nucleon coupling \(g_N\). For brevity we ignore these as they are not relevant in the context of this thesis. Future helioscopes like IAXO may be able to be sensitive to at least \(g_N\) though (Luzio et al. 2022).

4.3.1. Notes on KSVZ, DFSZ    extended

KSVZ

The Lagrangian for this model can be written (Shifman, Vainshtein, and Zakharov 1980)

\begin{equation} \mathcal{L}_{\text{KSVZ}} = \overline{Q}\slashed{D}Q - h \left( σ \overline{Q}_R Q_L + σ^{\dag} \overline{Q}_L Q_R \right) + \partial^{μ} σ^{\dag} \partial_{μ} σ + m^2 σ^{\dag} σ - λ \left( σ^{\dag} σ \right)^2, \label{eq:theory:axion:KSZV_lagrangian} \end{equation}

with the dimensionless Yukawa coupling \(h\). The vacuum expectation value for the field \(σ\) calculates to

\begin{equation} f_a = ⟨σ⟩ \equiv σ_0 = \frac{m}{\sqrt{2 λ}}. \end{equation}

The mass of this axion turns out to be exactly as for the standard axion eq. \ref{eq:theory:axion:effective_axion_mass}, with the replacement \(\nu_{\text{EW}} \rightarrow ⟨σ⟩\). This means, the KSVZ model can be seen as the simplest extension of the standard axion, which allows for an arbitrary symmetry breaking scale \(f_a\).

DFSZ

While not needing an additional superheavy quark it relies on two scalar Higgs doublet fields. \(\Phi_u\) has hypercharge \(-1\) and couples to \(u\) type right-handed quarks, whereas \(\Phi_d\) has hypercharge \(+1\) and couples to \(d\) type right-handed quarks as well as leptons. The different Higgs field have different vacuum expectation values, with the requirement

\begin{equation} \nu^2_{\text{EW}} = \nu^2_u + \nu^2_d, \end{equation}

while the additional degree of freedom allows for \(\sqrt{2} ⟨σ⟩ \equiv f_{σ} \gg \nu_{\text{EW}}\). The mass term is again similar to the standard axion and KSVZ axion

\begin{equation} m_a = m_{a0} / N_g, \end{equation}

where \(N_g\) is the number of quark generations of the theory. Again \(\nu_{\text{EW}}\) is replaced by \(⟨σ⟩\). One interesting property of the DSFZ axion is its coupling directly to electrons (and other leptons (Kim and Carosi 2010, 2019)), given by (Redondo 2013; Peccei 2008b):

\begin{equation} \mathcal{L}_{al} = -i \frac{\nu_d^2}{\nu^2_{\text{EW}}} \frac{m_l}{f_a} a \overline{l} \gamma^5 l. \end{equation}

where \(l\) refers to lepton. This type of coupling may allow for easier detection of axions than models only including axion photon couplings at tree level.

4.4. Implications for axion interactions - conversion probability

Starting with the Lagrangian in eq. \eqref{eq:theory:axion:general_axion_couplings} and extending it by the Lagrangian for a free photon \(A_ν\),

\[ \mathcal{L} = \mathcal{L}_{a,\text{eff}} + \mathcal{L}_γ = \mathcal{L}_{a,\text{eff}} - \frac{1}{4} F_{μν} F^{μν}, \]

where \(F_{μν} = ∂_μ A_ν - ∂_ν A_μ\) is the electromagnetic field strength tensor. And noting that the axion interaction term of eq. \eqref{eq:theory:axion:general_axion_couplings} can be rewritten as,

\[ \mathcal{L}_{aγ} = \frac{1}{4}g_{aγ} \tilde{F}^{μν} F_{μν} a = -g_{aγ} a \vec{E}·\vec{B}, \]

we can apply the Euler Lagrange equations to both the axion \(a\) and photon \(A_ν\) to derive a modified Klein-Gordon equation for the axion,

\[ \left(\Box + m_a²\right) a = \frac{1}{4}g_{aγ} F_{μν} \tilde{F}^{μν}, \]

which has a photon source term. Similarly, for the photon equation of motion we derive the homogeneous Maxwell equations with an axion source term,

\[ ∂_μ F^{μν} = g_{aγ} (∂_μ a) \tilde{F}^{μν}. \]

Without going into too much detail, let's shortly sketch how one derives the axion-photon conversion probability from here.

By choosing a suitable gauge and specifying directions of electric and magnetic fields in a suitable coordinate system, we can then derive the mixing between photon and axion states. For example if the propagation of particles is along the \(z\) axis and we fix the two degrees of freedom of \(A_ν\) by the Lorenz gauge (\(∂_μ A^μ = 0\)) and Coulomb gauge (\(\vec{\nabla}·\vec{A} = 0\)) we can derive a single equation of motion for the axion and \(A_ν\) field by starting with a plane wave approach. We obtain

\[ \left[ (ω² + ∂²_z) \mathbf{1} - \mathbf{M} \right] \vektor{A_{\perp}(z) \\ A_{\parallel}(z) \\ 0} = 0, \]

with the parallel and orthogonal polarization of the photon \(A_{\parallel}\) and \(A_{\perp}\), respectively and matrix \(\mathbf{M}\):

\[ \mathbf{M} = \mtrix{ m²_γ & 0 & 0 \\ 0 & m²_γ & -ω g_{aγ} B_T \\ 0 & -ω g_{aγ} B_T & m²_a \\ } \text{ where } m²_γ = ω²_p. \]

Here effects of quantum electrodynamic (QED) vacuum polarization and other polarization effects are ignored. \(ω_p\) refers to the plasma frequency, \(ω_p = 4πα n_e / m_e\), where \(n_e\) the electron density and \(m_e\) the electron mass. The mass \(m_γ\) refers to an effective photon mass that can appear in media, \(ω\) is the frequency and \(B_T\) the transverse magnetic field. The constant magnetic field \(B_T\) appears, because we assume for our purpose the magnetic field is constant along \(z\), the propagation direction of the photon.

By recognizing that the orthogonal component \(A_{\perp}\) is decoupled from the other two, the problem reduces to a 2 dimensional equation. Note that a side effect of this decoupling is that photons produced from an incoming axion in a magnetic field are always linearly polarized in the direction of the external magnetic field! Further, this equation can be linearized in the ultra relativistic limit \(m_γ² \ll ω²\) to

\[ \left[ \left( ω + i∂_z \right) \mathbf{1} - \frac{\mathbf{M}}{2ω} \right] \vektor{ A_{\parallel}(z) \\ a(z) } = 0. \]

As the mass matrix \(\mathbf{M}\) is non diagonal, the fields \(A_{\parallel}\) and \(a\) are interaction eigenstates and not propagation eigenstates. If we wish to compute the axion to photon conversion probability, we need the propagation eigenstates however. Transforming from one to the other is done by a regular rotation matrix \(\mathbf{R}\), which diagonalizes \(\mathbf{M}/2ω\). In the basis of the propagation eigenstates the fields \(A'_{\parallel}\) and \(a'\) are then decoupled and can be easily solved by a plane wave solution. The fields we can measure in an experiment are those of the interaction eigenstates of course. 5

The interaction eigenstates after a distance \(z\) can therefore be expressed by

\[ \vektor{ A_{\parallel}(z) \\ a(z) } = \mathbf{R}^{-1} \mathbf{M_{\text{diag}}} \mathbf{R} \vektor{ A_{\parallel}(0) \\ a(0) }, \]

where \(\mathbf{M_{\text{diag}}}\) is the diagonalized mass matrix,

\[ \mtrix{ e^{-i λ_+ z} & 0 \\ 0 & e^{-λ_- z} }, \]

where \(λ_{+, -}\) are its eigenvalues and coefficients in the exponential of the plane wave solutions of the propagation eigenstate fields \(A'_{\parallel}\) and \(a'\) given by

\[ λ_{+,-} = \pm \frac{1}{4ω} \sqrt{ \left(ω²_P - m²_a\right)² + \left(2 ω g_{aγ} B_T\right)²}. \]

Then finally, one can compute the conversion probability by starting from initial conditions where no electromagnetic field is present, \(A_{\parallel}(0) = 0, a(0) = 1\). Computing the resulting \(A_{\parallel}(z)\) with these conditions yields the expression which needs to be squared for the probability to measure a photon at distance \(z\) when starting purely from axions in an external, transverse magnetic field \(B_T\)

\begin{equation} \label{eq:theory:axion_interaction:conversion_probability} P_{a↦γ}(z) = |{A_{\parallel}(z)}|² = \left( \frac{g_{aγ} B_T z}{2} \right)² \left(\frac{\sin\left(\frac{q z}{2}\right)}{\frac{q z}{2}}\right)², \end{equation}

with \(q = \frac{m²_γ - m²_a}{2ω}\) and we dropped additional terms \(∝ g_{aγ} B_T\) as arguments to \(\sinc = \sin(x)/x\), because they are extremely small compared to \(q z\) for reasonable axion masses, magnetic fields and coupling constants.

If further the coherence condition \(qL < π\) is such that \(qL \ll π\), the \(\sinc\) term approaches 1 and the relevant conversion probability is finally:

\begin{equation} \label{eq:theory:conversion_prob} P_{a↦γ, \text{vacuum}} = \left(\frac{g_{aγ} B L}{2} \right)^2, \end{equation}

where we dropped the \(T\) suffix for the transverse magnetic field and replaced \(z\) by the more apt \(L\) for the length of a magnet. This is the case for long magnets and/or low axion masses, which is generally applicable in this thesis. The point at which this condition does not strictly hold anymore is the axion mass at which helioscope experiments start to lose sensitivity.

Note that the above conversion probability is given in natural units, specifically Lorentz-Heaviside units, \(c = \hbar = ε_0 = 1\) (meaning \(α = e²/4π \approx 1/137\)). Arguments (\(B, L\)) need either be converted to natural units as well or the missing factors need to be added. The same equation in SI units is given by:

\[ P_{a↦γ, \text{vacuum}} = ε_0 \hbar c^3 \left( \frac{g_{aγ} B L}{2} \right)^2. \]

A detailed derivation for the above can be found in (Masaki, Aoki, and Soda 2017). 6 An initial derivation for the first axion helioscope prototype is found in (Van Bibber et al. 1989) based on (Raffelt and Stodolsky 1988). Sikivie gives expected rates in his groundbreaking papers about axion experiments, (Sikivie 1983, 1985) but is extremely short on details. Another source in the form of (Raffelt 1996) in which G. Raffelt covers a very large number of topics relevant to axion searches.

4.4.1. Effects of a buffer gas

As seen in the conversion probability above, there is a term for an effective photon mass \(m_γ\) as part of \(q\). And indeed, \(q\) becomes zero if \(m_γ = m_a\), which means the suppression effect of the \(\sinc\) term disappears. This is something that can be used to increase the conversion probability inside of a magnet, by filling it with a buffer gas (for example helium), as introduced in (Raffelt and Stodolsky 1988; van Bibber et al. 1989). However, one also needs to account for the attenuation effect of the gas on the produced X-rays. As such the derivation above needs to include this as part of the evolution of the field \(\vec{A}\) 7. By doing this and following the rest of the derivation, the conversion probability in the presence of a buffer gas comes out to:

\begin{equation} \label{eq:theory:full_conversion_prob} P_{a\rightarrow\gamma} = \left(\frac{g_{a\gamma} B}{2}\right)^2 \frac{1}{q^2 + \Gamma^2 / 4} \left[ 1 + e^{-\Gamma L} - 2e^{-\frac{\Gamma L}{2}} \cos(qL)\right], \end{equation}

where \(\Gamma\) is the inverse absorption length for photons (or attenuation length), \(B\) the transverse magnetic field, \(L\) the length of the magnetic field and \(q\) the axion-photon momentum transfer given by:

\[ q = \left|\frac{m_{\gamma}^2 - m_a^2}{2E_a}\right| \]

One can easily see that this reduces to the vacuum case we discussed before, where the attenuation length \(Γ\) and effective photon mass \(m_γ\) are zero.

Filling the magnet with a buffer gas to induce an effective photon mass and thus minimizing \(q\) has been done at CAST with \(\ce{He4}\) and \(\ce{He3}\) fillings, as we will mention in chapter 5.1.

Note: for a potential buffer gas run in BabyIAXO I did some calculations about the required gas pressure steps and the effect of different possible filling configurations. These are not really relevant for this thesis, but can be found in 8 both as a PDF as well as the original Org file plus the tangled source code.

4.4.1.1. Simplification    extended

The conversion probability simplifies to:

\begin{align} P_{a\rightarrow\gamma, \text{vacuum}} &= \left(\frac{g_{a\gamma} B}{2}\right)^2 \frac{1}{q^2} \left[ 1 + 1 - 2 \cos(qL) \right] \\ P_{a\rightarrow\gamma, \text{vacuum}} &= \left(\frac{g_{a\gamma} B L}{2} \right)^2 \left(\frac{\sin\left(\frac{qL}{2}\right)}{ \left( \frac{qL}{2} \right)}\right)^2 \end{align}

4.4.2. Dropping an additional \(Δ\) term    extended

Notation:

\[ q = \frac{m²_γ - m²_a}{2ω} \]

and

\[ Δ = \frac{-g_{aγ} B_T}{2} \]

In principle the conversion probability contains

\[ \sinc²(\frac{πz}{L_{\text{osc}}}) \]

with the osciallation length

\[ L_{\text{osc}} = \frac{π}{\sqrt{ \left( \frac{q}{2} \right)² + Δ² } } \]

I assume we can safely drop the term \(Δ²\), as it is likely orders of magnitude smaller than the first term. Let's check that:

import unchained

defUnit(GeV⁻¹)
proc q(m_a: meV, ω: keV): keV =
  result = (m_a*m_a / (2 * ω)).to(keV)

proc Δ(g_aγ: GeV⁻¹, B: Tesla): keV =
  result = (-g_aγ * B.toNaturalUnit() / 2.0).to(keV)

echo "q² = ", q(1.meV, 3.keV)^2
echo "Δ² = ", Δ(1e-12.GeV⁻¹, 9.T)^2

Yeah, as expected that's the reason.

4.4.3. Deriving the missing constants in the conversion probability [/]    extended

  • [ ] Move this into an appendix?

The conversion probability is given in natural units. In order to plug in SI units directly without the need for a conversion to natural units for the magnetic field and length, we need to reconstruct the missing constants.

The relevant constants in natural units are:

\begin{align*} ε_0 &= \SI{8.8541878128e-12}{A.s.V^{-1}.m^{-1}} \\ c &= \SI{299792458}{m.s^{-1}} \\ \hbar &= \frac{\SI{6.62607015e-34}{J.s}}{2π} \end{align*}

which are each set to 1.

If we plug in the definition of a volt we get for \(ε_0\) units of:

\[ \left[ ε_0 \right] = \frac{\si{A^2.s^4}}{\si{kg.m^3}} \]

The conversion probability naively in natural units has units of:

\[ \left[ P_{aγ, \text{natural}} \right] = \frac{\si{T^2.m^2}}{J^2} = \frac{1}{\si{A^2.m^2}} \]

where we use the fact that \(g_{aγ}\) has units of \(\si{GeV^{-1}}\) which is equivalent to units of \(\si{J^{-1}}\) (care has to be taken with the rest of the conversion factors of course!) and Tesla in SI units:

\[ \left[ B \right] = \si{T} = \frac{\si{kg}}{\si{s^2.A}} \]

From the appearance of \(\si{A^2}\) in the units of \(P_{aγ, \text{natural}}\) we know a factor of \(ε_0\) is missing. This leaves the question of the correct powers of \(\hbar\) and \(c\), which come out to:

\begin{align*} \left[ ε_0 \hbar c^3 \right] &= \frac{\si{A^2.s^4}}{\si{kg.m^3}} \frac{\si{kg.m^2}}{\si{s}} \frac{\si{m^3}}{\si{s^3}} \\ &= \si{A^2.m^2}. \end{align*}

So the correct expression in SI units is:

\[ P_{aγ} = ε_0 \hbar c^3 \left( \frac{g_{aγ} B L}{2} \right)^2 \]

where now only \(g_{aγ}\) needs to be expressed in units of \(\si{J^{-1}}\) for a correct result using tesla and meter.

4.4.4. Further notes on units of conversion probability    extended

The conversion probability follows the derivation by Biljana and Kreso: Axion-photon_conversion_report_biljana.svg Fortunately, they state explicitly, which units they use, to quote:

We stress that the equations above are written in terms of natural, rationalized electromagnetic units (natural Lorentz-Heaviside units) where hbar = c = 1 and the fine-structure constant is given as α = e² / 4π ≈ 1/137.

This is very important to know that we did not miss some 4π or √4π factor.

The fine structure constant is defined by

\[ α = \frac{e²}{4π ε_0 \hbar c} \]

in SI units. This means, \(ε_0 = 1\) as well (but explicitly not something like \(4π ε_0 = 1\), which one also finds, e.g. here http://ilan.schnell-web.net/physics/natural.pdf).

A PDF about natural unit systems I personally like a lot:

https://www.seas.upenn.edu/~amyers/NaturalUnits.pdf

It describes natural units in the same way as we use them here, i.e. 'Lorentz-Heaviside'. These are also implemented as the (current only) natural units in unchained.

In particular in these units, the conversion factors for tesla and meter are as follows:

import unchained
echo 1.T, " in natural units = ", 1.T.toNaturalUnit()
echo 1.m, " in natural units = ", 1.m.toNaturalUnit()

which are the relevant conversion factors, if one wishes to work with the natural unit version (and not relying on a library with an option to convert units for you).

See the table of the aforementioned NaturalUnits.pdf to deduce these conversion factors yourself.

Also, there exists a natural unit conversion map for GNU Units:

https://github.com/misho104/natural_units

which also explicitly mentions it uses Lorentz-Heaviside units. Therefore, feel free to use that to help with your conversions.

4.4.5. Full simplification of conversion probability in vacuum    extended

The full simplification for the vacuum case from the buffer gas conversion probability is as follows:

\begin{align*} P_{a\rightarrow\gamma} &= \left(\frac{g_{a\gamma} B}{2}\right)^2 \frac{1}{q^2 + \Gamma^2 / 4} \left[ 1 + e^{-\Gamma L} - 2e^{-\frac{\Gamma L}{2}} \cos(qL)\right] \\ \text{for vacuum } Γ &= 0, m_γ = 0 \text{ and thus} \\ P_{a\rightarrow\gamma, \text{vacuum}} &= \left(\frac{g_{a\gamma} B}{2}\right)^2 \frac{1}{q^2} \left[ 1 + 1 - 2 \cos(qL) \right] \\ P_{a\rightarrow\gamma, \text{vacuum}} &= \left(\frac{g_{a\gamma} B}{2}\right)^2 \frac{2}{q^2} \left[ 1 - \cos(qL) \right] \\ P_{a\rightarrow\gamma, \text{vacuum}} &= \left(\frac{g_{a\gamma} B}{2}\right)^2 \frac{2}{q^2} \left[ 2 \sin^2\left(\frac{qL}{2}\right) \right] \\ P_{a\rightarrow\gamma, \text{vacuum}} &= \left(g_{a\gamma} B\right)^2 \frac{1}{q^2} \sin^2\left(\frac{qL}{2}\right) \\ P_{a\rightarrow\gamma, \text{vacuum}} &= \left(\frac{g_{a\gamma} B L}{2} \right)^2 \left(\frac{\sin\left(\frac{qL}{2}\right)}{ \left( \frac{qL}{2} \right)}\right)^2 \\ P_{a\rightarrow\gamma, \text{vacuum}} &= \left(\frac{g_{a\gamma} B L}{2} \right)^2 \left(\frac{\sin\left(\delta\right)}{\delta}\right)^2 \\ \end{align*}

4.4.6. Details on the polarization of axion induced photons    extended

These notes refer to the notes by Bijlana and Kreso Axion-photon_conversion_report_biljana.svg

UPDATE: <2023-09-08 Fri 20:14> I just noticed the following sentence in Biljana's notes:

The Faraday rotation term \(∆_R\) , which depends on the energy and the longitudinal component of the external magnetic field, couples the photon polarization states A⊥ and Ak . While Faraday rotation is important when analyzing polarized sources of photons, it plays no role in the problem at hand.

which refers to off diagonal elements \(Δ_R\), which were dropped in their calculation! While we don't consider polarized sources, do we still need to consider this?

See eq. (3.4) (reproduced here)

\[ a(z, t) = a(z) e^{-i ω t} \:\: , \:\: \vec{A}(z, t) = \vektor{A_{\perp}(z) \\ A_{\parallel}(z) \\ 0} e^{-iωt} \]

for the definition of the axion and photon in the context. The equation of motion is eq. (3.8)

\[ \left[ (ω² + ∂²_z) \mathbf{1} - \mathbf{M} \right] = \vektor{A_{\perp}(z) \\ A_{\parallel}(z) \\ 0} = 0 \]

and with the definition of eq. (3.14)

\[ \mathbf{M} = \mtrix{ m²_γ & 0 & 0 \\ 0 & m²_γ & -ω g_{aγ} B_T \\ 0 & -ω g_{aγ} B_T & m²_a \\ } \text{ where } m²_γ = ω²_p \]

we can (as mentioned in the text) see that indeed the orthogonal \(\vec{A}\) component is independent of the axion field (the first row of \(\mathbf{M}\) only has an entry in the first column, i.e. the product only yields an equation for \(A_{\perp}\) alone without the axion field.

Due to the difference between the interaction (\(A_{\parallel}, a\)) and propagation eigenstates (\(A'_{\parallel}, a'\)) (connected via a rotation matrix) the initial separate eigenstates \(A_{\parallel}\) and \(a\) end up mixing in the propagation eigenstates (see page 9). Eq. 3.34

\[ i∂_z \vektor{A'_{\parallel}(z) \\ a'(z)} = \mathbf{H_D} \vektor{A'_{\parallel}(z) \\ a'(z)} \]

is then the equation of motion for the propagation eigenstates. The solutions are then eq. 3.38

\begin{align*} A_{\parallel}(z) &= A_{\parallel}(0) \left( \cos² θ e^{-iλ_+z} + \sin² θ e^{-iλ_-z}\right) + a(0) \frac{\sin 2θ}{2} \left( e^{-iλ_+z} - e^{-iλ_-z}\right) \\ a(z) &= A_{\parallel}(0) \frac{\sin 2θ}{2} \left( e^{-iλ_+z} - e^{-iλ_-z}\right) + a(0) \left( \sin² θ e^{-iλ_+z} + \cos² θ e^{-iλ_-z}\right) \\ \end{align*}

which - with our assumption in our experiment

\[ A_{\parallel}(0) = 0 \text{ and } a(0) = 1 \]

i.e. we start with purely axions and no photons before the magnet - can then be simplified to 3.39 and 3.40

\begin{align*} A_{\parallel}(z) &= a(0) \frac{\sin 2θ}{2} \left( e^{-iλ_+z} - e^{-iλ_-z}\right) \\ a(z) &= \sin² θ e^{-iλ_+z} + \cos² θ e^{-iλ_-z} \\ \end{align*}

What this implies is that the photon contribution after mixing that can end up as a detected physical photon is only of \(A_{\parallel}\) type, which (again if my rusty understanding is not failing me) implies that the produced photons all have the same polarization, the one parallel to the constant \(\vec{B}\) field (compare with fig. 1).

4.5. Solar axion flux

The effective Lagrangian as shown in eq. \eqref{eq:theory:axion:general_axion_couplings} allows for multiple different axion interactions, which result in the production of axions in the Sun. For KSVZ-like axion models (models with only \(g_{aγ}\)) the only interaction allowing for axion production in the Sun is the Primakoff effect 9 for axions. For DFSZ models with an axion-electron coupling \(g_{ae}\) multiple other production channels are viable.

One of the first papers to look at the implications of the axion in terms of astrophysical phenomena is (Mikaelian 1978) by K. Mikaelian. G. Raffelt expanded on this later in (Raffelt 1986) with calculations for the Compton and Bremsstrahlung production rates for DFSZ axion models. In (Raffelt 1988) he further calculates the production rate for the Primakoff effect (and later reviews the physics (Raffelt 1996)). J. Redondo combined all production processes in (Redondo 2013) to compute a full solar axion flux based on the axion-electron coupling \(g_{ae}\) using numerical calculations of the expected metallicity contents at different points in the Sun. This is done by making use of the opacities for different elements at different pressures and temperatures as tabulated by the 'Opacity Project' (Team 1995; Hummer and Mihalas 1988; Seaton 1987, 2005; Seaton et al. 1994; Badnell et al. 2005), based on the values provided by the numerical AGSS09 (Asplund et al. 2009) solar model. All these considered axion production channels are:

  • (\(P\)) Primakoff production via \(g_{aγ}\) in both KSVZ and DFSZ axion models \[ γ + γ ↦ a \]
  • (\(\text{ff}\)) electron ion bremsstrahlung (in radio astronomy low energetic cases with photons are also called free-free radiation) \[ e + Z \longrightarrow e + Z + a \]
  • (\(ee\)) electron electron bremsstrahlung (Raffelt 1986), \[ e + e \longrightarrow e + e + a \]
  • (\(\text{fb}\)) electron capture (also called recombination or free-bound electron transitions) \[ e + Z \longrightarrow a + Z^- \]
  • (\(C\)) Compton scattering (Raffelt 1986) \[ e + \gamma \longrightarrow e + a \]
  • (\(\text{bb}\)) and de-excitation (bound-bound electron transitions) via an axion \[ Z^* \longrightarrow Z + a \]

See fig. 2 for the corresponding Feynman diagrams.

axion_prod_channels_javi.svg
Figure 2: Figure 2: Feynman diagrams of all contributing axion production channels in the Sun for non-hadronic models. In hadronic models only the Primakoff interaction has meaningful contributions, because axion-electron couplings only arise at loop level. Taken from Redondo_2013.

With these production channels we can write down the expected axion flux on Earth, based on the production rate per volume in the Sun as the integral

\begin{equation} \frac{d Φ_a}{dω} = \frac{1}{4π R²_{\text{Earth}}} \int_{\text{Sun}} \mathrm{d}V\, \frac{4π ω²}{ (2π)³ } Γ_a(ω). \end{equation}

where

\begin{equation} Γ_a(\omega) = Γ^{\text{ff}}_a + Γ^{\text{fb}}_a + Γ^{\text{bb}}_a + Γ^C_a + Γ^{ee}_a + Γ^P_a \end{equation}

are all contributing axion production channels. The superscripts correspond to the bullet points above. Most important without going into the details of how each \(Γ\) can be expressed, is that they scale with the coupling constant squared. \(g²_{aγ}\) for \(Γ^P_{e}\) and \(g²_{ae}\) for the others. In (Redondo 2013) these production rates are expressed by relating them to the corresponding photon production rates for these processes, which are well known.

For a detailed look, see (Redondo 2013) and the master thesis of Johanna von Oy (Von Oy 2020) in which she – among other things – reproduced the calculations done by Redondo. Her work is used as part of this thesis to compute the axion production as required for the expected axion flux in the limit calculation (and provides the data for the plots in this section). The code responsible for computing the emission rates for different axion models is (Von Oy 2023), in particular the readOpacityFile program. It also uses the Opacity Project as the basis to compute the opacities for different elements.

Fig. 3(a) shows the differential axion flux arriving on Earth based on the different contributing interactions, in this case using \(g_{ae} = \num{1e-13}, g_{aγ} = \SI{1e-12}{GeV⁻¹}\). We see that the total flux (blue line) peaks at roughly \(\SI{1}{keV}\), with the small spikes from atomic interactions visible. The Primakoff flux has been amplified by a factor of \(100\) here to make it more visible, as for this choice of coupling constants its contribution is negligible.

As a result of the expected extremely low mass of the axion, the expected solar axion spectrum for Primakoff only models is essentially a blackbody spectrum corresponding to the temperatures near the core of the Sun, \(\mathcal{O}(\SI{15e6}{K})\), see fig 3(b) and compare it to the Primakoff flux on the left 10. Partially for this reason, in many cases analytical expressions are given to describe the Primakoff axion flux, which are solutions obtained for specific solar models. One such recent result is from (Dafni et al. 2019),

\begin{equation} \label{eq:theory:solar_axion_flux:primakoff} \frac{\dd Φ_a}{\dd E_a} = Φ_{P10} \left( \frac{g_{aγ}}{\SI{1e-10}{GeV⁻¹}} \right)² \frac{E^{2.481}_a}{e^{E_a / \SI{1.205}{keV}}} \end{equation}

where \(Φ_{P10} = \SI{6.02e10}{keV⁻¹.cm⁻².s⁻¹}\) and \(E_a\) the energy of the axion in \(\si{keV}\). (Dafni et al. 2019) also contains analytic expressions for the Compton and Bremsstrahlungs components. For the atomic processes this is not possible.

Figure 3(a): Differential axion flux
Figure 3(b): $\SI{15e6}{K}$ blackbody spectrum
Figure 3: 3(a)Differential solar axion flux based on the different interaction types using $g_{ae} = \num{1e-13}, g_{aγ} = \SI{1e-12}{GeV^{-1}}$. The Primakoff contribution was scaled up by a factor 100 to make it visible, as at these coupling constants the $g_{ae}$ contributions dominate. 3(b) shows a blackbody spectrum corresponding to $\SI{15e6}{K}$, roughly the temperature at the solar core. Up to a scaling factor it is essentially the Primakoff flux.

Fig. 4(a) shows how the solar axion flux (for DFSZ models) depends both on the energy and the relative radius in the Sun. We can see clearly that the major part of the axion flux comes from a region between \(\SIrange{7.5}{17.5}{\%}\) of the solar radius. The reason is the cubic scaling of the associated volumes per radius on the lower end and dropping temperatures and densities at the upper end. Interesting substructure due to the details of the axion-electron coupling is visible. The radial component alone comparing it to KSVZ models is seen in fig. 4(b), where we can see that the DFSZ flux drops off significantly at a specific radius resulting in the net flux from slightly smaller radii.

Figure 4(a): Flux
Figure 4(b): Radial emission
Figure 4: 4(a)Solar axion flux for DFSZ models showing the flux dependence on both the energy and relative solar radius. Dominant contribution below $0.3 R_{\odot}$. 4(b)Difference in the radial contributions of axion flux for KSVZ models against DFSZ models. DFSZ model production is constrained to slightly smaller radii.

The used solar model is an important part of the input for the calculation of the emission rate and thus differential flux. (Hoof, Jaeckel, and Thormaehlen 2021) conclude – based on a similar code (Hoof and Thormaehlen 2021) – that the statistical uncertainty of the solar models is \(\mathcal{O}(\SI{1}{\%})\), while the systematic uncertainty can reach up to \(\mathcal{O}(\SI{5}{\%})\). This is an important consideration for the systematic uncertainties later on.

4.5.1. Primakoff flux    extended

Including analytical equation for flux… :)

import unchained, ggplotnim, math, chroma, ginger
defUnit(keV⁻¹•m⁻²•yr⁻¹)
defUnit(keV⁻¹•cm⁻²•s⁻¹)
defUnit(GeV⁻¹)

proc axionFluxPrimakoff(E_a: keV, g_aγ: GeV⁻¹): keV⁻¹•cm⁻²•s⁻¹ =
  ## dΦ_a/dE taken from paper about first CAST results cite:PhysRevLett.94.121301
  let g₁₀ = g_aγ / 1e-10.GeV⁻¹ # * 10e10.GeV¹ #
  result = g₁₀^2 * 3.821e10.cm⁻²•s⁻¹•keV⁻¹ * (E_a / 1.keV)^3 / (exp(E_a / (1.103.keV)) - 1)

proc axFluxPerYear(E_a: keV, g_aγ: GeV⁻¹): keV⁻¹•m⁻²•yr⁻¹ =
  result = axionFluxPrimakoff(E_a, g_aγ).to(keV⁻¹•m⁻²•yr⁻¹)

proc axionFluxPrimakoffMasterThesis(ω: keV, g_ay: GeV⁻¹): keV⁻¹•m⁻²•yr⁻¹ =
  # axion flux produced by the Primakoff effect
  # in units of m^(-2) year^(-1) keV^(-1)
  # From the CAST 2013 paper on axion electron coupling, eq 3.1
  result = 2.0 * 1e18.keV⁻¹•m⁻²•yr⁻¹ * (g_ay / 1e-12.GeV⁻¹)^2 * pow(ω / 1.keV, 2.450) * exp(-0.829 * ω / 1.keV)

let E = linspace(1e-3, 14.0, 1000)
let df = seqsToDf(E)
  .mutate(f{float: "Flux" ~ axionFluxPrimakoff(`E`.keV, 1e-11.GeV⁻¹).float})
  .mutate(f{float: "FluxYr" ~ axFluxPerYear(`E`.keV, 1e-11.GeV⁻¹).float})
  .mutate(f{float: "FluxMSc" ~ axionFluxPrimakoffMasterThesis(`E`.keV, 1e-11.GeV⁻¹).float})    
ggplot(df, aes("E", "Flux")) +
  geom_line() +
  #geom_line(aes = aes(y = "FluxMSc"), color = some(parseHex("0000FF"))) + 
  ggtitle("Solar axion flux due to Primakoff production, g_aγ = 10⁻¹¹·GeV⁻¹") +
  xlab("Energy [keV]") +
  #ylab("Axion flux [keV⁻¹·cm⁻²·s⁻¹]") +
  ylab("Axion flux [keV⁻¹·m⁻²·yr⁻¹]") +  
  ggsave("/tmp/primakoff_axion_flux.pdf")

ggplot(df.mutate(f{"Flux" ~ `Flux` / 1e8}), aes("E", "Flux")) +
  geom_line() +
  #xlab("Energy [keV]", tickFont = font(12.0), margin = 1.5) +
  xlab(r"\fontfamily{lmss}\selectfont Energy [$\si{\keV}$]", margin = 2.0, font = font(16.0),
       tickFont = font(16.0)) +
  xlim(0, 14) + 
  #ylab("Axion flux [10¹⁰ keV⁻¹·cm⁻²·s⁻¹]", margin = 1.5) +
  ylab(r"\fontfamily{lmss}\selectfont Axion flux [\SI[print-unity-mantissa=false]{1e11}{\keV^{-1} \cm^{-2} \second^{-1}}]",
       margin = 2.0,
       font = font(16.0)) + 
  #     tickFont = font(12.0)) +
  #ggtitle(r"Expected solar axion flux, g_aγ = 10⁻¹⁰ GeV⁻¹", titleFont = font(12.0)) +
  annotate(r"\fontfamily{lmss}\selectfont Expected solar axion flux" &
    r"\\$g_{aγ} = \SI[print-unity-mantissa=false]{1e-11}{\GeV^{-1}}$", #10⁻¹⁰ GeV⁻¹",
           x = 6.2, y = 6.2, 
           font = font(16.0),
           backgroundColor = transparent) +
  #ggtitle(r"Expected solar axion flux, $g_{aγ} = \SI{1e-11}{\GeV^{-1}}$", titleFont = font(12.0)) + 
  #ggsave("/tmp/cristina_primakoff_axion_flux.pdf", width = 400, height = 300) #, useTeX = true, standalone = true)
  ggsave("/tmp/cristina_primakoff_axion_flux.pdf", useTeX = true, standalone = true)
  

defUnit(m⁻²•yr⁻¹)  
echo 1.cm⁻²•s⁻¹.to(m⁻²•yr⁻¹)

There are different analytical expressions for the solar axion flux for Primakoff production. These stem from the fact that a solar model is used to model the internal density, temperature, etc. in the Sun to compute the photon distribution (essentially the blackbody radiation) near the core. From it (after converting via the Primakoff effect) we get the axion flux.

Different solar models result in different expressions for the flux. The first one uses an older model, while the latter ones use newer models.

4.5.2. Generate figures for Primakoff and axion flux from readOpacityFile output [4/5]    extended

  • Use logic we use in TrAXer to compute differential fluxes as a function of radii.
  • [X] Do it for each type of flux independently
  • [X] Compute a smooth version of the radial distribution comparing axion-electron to Primakoff!
  • [X] Black body spectrum -> Produced below in a separate plot!
  • [X] Radial emission vs energy heatmap!!
  • [ ] ADD COMMAND TO PRODUCE CSV FILE!

First we produce the differential flux CSV file. For more details, see sec. 37.4.3.

./readOpacityFile \
    --suffix "_0.989AU" \
    --distanceSunEarth 0.9891144450781392.AU \
    --fluxKind fkAxionElectronPhoton \
    --plotPath ~/phd/Figs/readOpacityFile/ \
    --outpath ~/phd/resources/readOpacityFile/

Then the differential flux by type:

import ggplotnim
const fluxPath = "~/phd/resources/readOpacityFile/solar_axion_flux_differential_g_ae_1e-13_g_ag_1e-12_g_aN_1e-15_0.989AU.csv"
let fCol = "Flux / keV⁻¹ m⁻² yr⁻¹"
let df = readCsv(fluxPath)
  .filter(f{string -> bool: `type` notin ["LP Flux", "TP Flux", "57Fe Flux"]})
  .mutate(f{float: "diffFlux" ~ (if (idx("type", string) == "Primakoff Flux"): idx(fCol) * 100 else: idx(fCol))})
  .mutate(f{string: "type" ~ (if (`type` == "Primakoff Flux"): "Primakoff·100" else: `type`)})
ggplot(df, aes("Energy [keV]",f{`diffFlux` / 1e20}, color = "type")) +
  geom_line() +
  xlab(r"Energy [$\si{keV}$]") + ylab(r"Flux [$\SI{1e20}{keV^{-1}.m^{-2}.yr^{-1}}$]") +
  xlim(0, 15) + 
  ggtitle("Differential axion flux arriving on Earth") +
  themeLatex(fWidth = 0.5, width = 600, baseTheme = sideBySide) +  # width 360?
  # ggshow(800, 480)
  ggsave("~/phd/Figs/axions/differential_solar_axion_flux_by_type.pdf", useTeX = true, standalone = true, width = 600, height = 360)

Now the radial distribution of the flux by coupling constant.

import std / [sequtils, algorithm]
import ggplotnim, unchained
import ggplotnim/ggplot_sdl2
type
  FluxData* = object
    fRCdf*: seq[float]
    diffFluxR*: seq[seq[float]]
    fluxesDf*: DataFrame
    radii*: seq[float]
    energyMin*: float
    energyMax*: float

const
  alpha = 1.0 / 137.0
  g_ae = 1e-13 # Redondo 2013: 0.511e-10
  gagamma = 1e-12 #the latter for DFSZ  #1e-9 #5e-10 #
  ganuclei = 1e-15 #1.475e-8 * m_a #KSVZ model #no units  #1e-7
  m_a = 0.0853 #eV
  m_e_keV = 510.998 #keV
  e_charge = sqrt(4.0 * PI * alpha)#1.0
  kB = 1.380649e-23
  r_sun = 696_342_000_000.0 # .km.to(mm).float # SOHO mission 2003 & 2006
  hbar = 6.582119514e-25 # in GeV * s
  keV2cm = 1.97327e-8 # cm per keV^-1
  amu = 1.6605e-24 #grams
  r_sunearth = 150_000_000_000_000.0

const factor = pow(r_sun * 0.1 / (keV2cm), 3.0) /
               (pow(0.1 * r_sunearth, 2.0) * (1.0e6 * hbar)) /
               (3.1709791983765E-8 * 1.0e-4) # for units of 1/(keV y m²)

import ggplotnim
proc getFluxRadiusCDF*(path: string): FluxData =
  var emRatesDf = readCsv(path)
  # get all radii and energies from DF so that we don't need to compute them manually (risking to
  # messing something up!)
  # sort both just to make sure they really *are* in ascending order
  let radii = emRatesDf["Radius"]
    .unique()
    .toTensor(float)
    .toSeq1D
    .sorted(SortOrder.Ascending)
  let energies = emRatesDf["Energy [keV]"]
    .unique()
    .toTensor(float)
    .toSeq1D
    .sorted(SortOrder.Ascending)
  var emRates = newSeq[seq[float]]()
  ## group the "solar model" DF by the radius & append the emission rates for all energies
  ## to the `emRates`
  for tup, subDf in groups(emRatesDf.group_by("Radius")):
    let radius = tup[0][1].toFloat
    #doAssert subDf["Energy [keV]", float].toSeq1D.mapIt(it.keV) == energies
    #doAssert radius == radii[k], "Input DF not sorted correctly!"
    emRates.add subDf["emRates", float].toSeq1D
  var
    fluxRadiusCumSum: seq[float] = newSeq[float](radii.len)
    diffRadiusSum = 0.0

  template toCdf(x: untyped): untyped =
    let integral = x[^1]
    x.mapIt( it / integral )

  var fluxesDf = newDataFrame()
  var diffFluxR = newSeq[seq[float]](radii.len)
  var rLast = 0.0
  for iRad, radius in radii:
    # emRates is seq of radii of energies
    let emRate = emRates[iRad]
    var diffSum = 0.0
    var diffFlux = newSeq[float](energies.len)
    for iEnergy, energy in energies:
      let dFlux = emRate[iEnergy] * (energy.float*energy.float) * radius*radius * (radius - rLast) * factor
      diffFlux[iEnergy] = dFlux
      diffSum += dFlux
    fluxesDf.add toDf({"Energy" : energies.mapIt(it.float), "Flux" : diffFlux, "Radius" : radius})      
    diffRadiusSum += diffSum
    fluxRadiusCumSum[iRad] = diffRadiusSum
    diffFluxR[iRad] = diffFlux
    rLast = radius
  result = FluxData(fRCdf: fluxRadiusCumSum.toCdf(),
                    fluxesDf: fluxesDf,
                    diffFluxR: diffFluxR,
                    radii: radii,
                    energyMin: energies.min,
                    energyMax: energies.max)

proc fluxToDf(data: FluxData, typ: string): DataFrame =
  result = toDf({ "Radius" : data.radii,
                  "fluxPerRadius" : data.diffFluxR.mapIt(it.sum),
                  "Type" : typ })
    .mutate(f{"fluxPerRadius" ~ `fluxPerRadius` / col("fluxPerRadius").max})
proc main =
  const ksvzPath = "~/CastData/ExternCode/AxionElectronLimit/resources/solar_model_dataframe_fluxKind_fkAxionPhoton_0.989AU.csv"
  const dfszPath = "~/CastData/ExternCode/AxionElectronLimit/resources/solar_model_dataframe_fluxKind_fkAxionElectronPhoton_0.989AU.csv"
  let ksvz = getFluxRadiusCdf(ksvzPath)
  let dfsz = getFluxRadiusCdf(dfszPath)
  var df = newDataFrame()
  df.add ksvz.fluxToDf("KSVZ"); df.add dfsz.fluxToDf("DFSZ")
  ggplot(df, aes("Radius", "fluxPerRadius", color = "Type")) +
    geom_line() +
    xlab("Relative solar radius") + ylab("Relative emission") +
    ggtitle("Radial emission: KSVZ ($g_{aγ}$), DFSZ ($g_{aγ}, g_{ae}$)", titleFont = font(11.0)) +
    themeLatex(fWidth = 0.5, width = 600, baseTheme = sideBySide) + # golden ratio or height = 360, ?
    #ggshow() 
    ggsave("/home/basti/phd/Figs/axions/solar_axion_radial_emission.pdf", useTeX = true, standalone = true) # , width = 800, height = 480)

  echo dfsz.fluxesDf
  let dfFlux = dfsz.fluxesDf
    .filter(f{`Energy` <= 10.0 and `Radius` <= 0.5})
  ggplot(dfFlux, aes("Energy", "Radius", fill = "Flux")) +
    geom_raster() +
    scale_fill_continuous(scale = (0.0, percentile(dfFlux["Flux", float], 99))) +
    ggtitle("Flux by energy and fraction of solar radius") +
    xlab(r"Energy [$\si{keV}$]") +
    ylab(r"Relative solar radius") +
    themeLatex(fWidth = 0.5, width = 600, baseTheme = sideBySide) + # golden ratio or height = 360, ?    
    xlim(0, 10) + ylim(0, 0.5) +
    legendPosition(0.72, 0.0) + 
    ggsave("/home/basti/phd/Figs/axions/flux_by_energy_vs_radius_axion_electron.pdf", useTeX = true, standalone = true)
    #ggshow() # 640, 480)
    
main() 

Hmm, this looks a bit bizarre. But comparing with

  • sampled_radii_axion_electron.svg
  • sampled_radii_primakoff.svg

the weird structure of the axion-electron emission is actually visible a bit (in the form of sharp drops, in particular near 0.2). The question is really "why", but well. If it was some kind of sampling issue, I would assume it comes out of the solar model.

4.5.3. Black body radiation in solar core    extended

Let's compute the black body radiation for the solar core and see if it matches the energy spectrum we expect for axions.

Planck's law is defined as CITE SOMETHING:

\[ B_ν(ν, T) = \frac{2hν³}{c²} \frac{1}{e^{hν/kT} - 1} \]

where \(ν\) is the frequency of the photon and \(T\) the temperature in Kelvin. \(k\) is of course the Boltzmann constant and \(h\) the Planck constant. Let's see what this looks like for \(T = \SI{15}{\mega\kelvin}\).

blackbody_spectrum_solar_core.svg
import ggplotnim, unchained, sequtils
import ggplotnim / ggplot_sdl2

#defUnit(s⁻¹)
#defUnit(μs⁻¹)
defUnit(Watt•Steradian⁻¹•Meter⁻²•NanoMeter⁻¹)
defUnit(Joule•Meter⁻²•Steradian⁻¹)

let T_sun = 15.MegaKelvin.to(Kelvin)

proc blackBody(ν: s⁻¹, T: Kelvin): Joule•Meter⁻²•Steradian⁻¹ =
  result = (2 * hp * ν^3 / c^2 / (exp(hp * ν / (k_B * T)) - 1)).to(Joule•Meter⁻²•Steradian⁻¹)

proc xrayEnergyToFreq(E: keV): s⁻¹ = 
  ## converts the input energy in keV to a correct frequency
  result = E.to(Joule) / hp
echo 1.keV.xrayEnergyToFreq
echo "Solar core temperature ", T_Sun, " in keV : ", T_Sun.toNaturalUnit().to(keV)
echo blackBody(1.μHz.to(Hz), T_sun)
echo blackBody(1.keV.xrayEnergyToFreq, T_sun)

let energies = linspace(0.01, 15.0, 1000)
let radiance = energies.mapIt(blackBody(it.keV.xrayEnergyToFreq, T_sun).float)
let df = seqsToDf(energies, radiance)
ggplot(df, aes("energies", "radiance")) + 
  geom_line() + 
  ggtitle(r"Black body radiation @ $T = \SI{15e6}{K}$") +
  xlab(r"Energy [$\si{keV}$]") +
  ylab(r"Radiance [$\si{J.m^{-2}.sr^{-1}}$]") +
  xlim(0, 15) +
  themeLatex(fWidth = 0.5, width = 600, baseTheme = sideBySide) + # golden ratio or height = 360, ?
  #ggshow() 
  ggsave("/home/basti/phd/Figs/blackbody_spectrum_solar_core.pdf", useTeX = true, standalone = true, width = 600, height = 360)

4.6. Chameleons

Chameleons are a different type of hypothetical particle. A scalar particle arising from extensions to general relativity, which acts as a "fifth force" and can be used to model dark energy. We will not go into detail about the underlying theory here. Refer to (Waterhouse 2006) for an in-depth introduction to chameleon gravity and (Brax and Davis 2015) on how they differ from other modified gravity models.

Suffice it to say, chameleonic theory also yields a coupling to photons, \(β_γ\), which can be utilized for conversion into X-rays in transverse magnetic fields. Similarly, they can therefore be produced in the Sun. However, the chameleon field has a peculiar property in that its mass is dependent on the density of the surrounding medium. This means chameleons cannot be produced in the solar core and escape. A suitable location to produce chameleons, is the solar tachocline region. A region of differential rotation at the boundary between the inner radiative zone and outer convective zone, at around \(0.7 · R_{\odot}\) in the Sun. This leads to strong magnetic fields at low enough densities. Due to the much lower temperatures present in the solar tachocline compared to the solar core, the peak of the solar chameleon spectrum is below \(\SI{1}{keV}\). See (Brax, Lindner, and Zioutas 2012) and (Anastassopoulos et al. 2015) for details about the chameleon production in the tachocline region and resulting spectrum.

Their possible detection with CAST was proposed in 2012 (Brax, Lindner, and Zioutas 2012), with actual searches being performed initially with a Silicon Drift Detector (SDD) (Anastassopoulos et al. 2015) and later with a GridPix detector using data taken in 2014/15 (Christoph Krieger 2018; Anastassopoulos et al. 2019). In addition, the KWISP (Kinetic WISP) detector (Cuendis et al. 2019; Baier, n.d.) was deployed at CAST, attempting search for the chameleon via their chameleon-matter coupling \(β_m\), which should lead to a chameleons interacting with a force sensor (without conversion to X-rays in the CAST magnet).

As the chameleon flux is highly dependent on details of solar physics that are understood to a much lesser extent than the temperature and composition in the solar core (required for axions), uncertainty on chameleon results is much larger. See also (Zanzi and Ricci 2015) for a discussion of chameleon fields in the context of solar physics.

The conversion probability for back conversion into X-rays in a magnetic field, assuming coherent conversion, is given by (Brax, Lindner, and Zioutas 2012):

\begin{equation} \label{eq:theory:chameleon_conversion_prob} P_{c↦γ} = \left( \frac{β_γ B L}{2 m_{\text{pl}}} \right)². \end{equation}

Here, \(m_{\text{pl}}\) is the reduced Planck mass,

\[ m_{\text{pl}} = \frac{M_{\text{pl}}}{\sqrt{8 π}} = \sqrt{ \frac{ \hbar c }{ 8 π G } }, \]

i.e. using natural units with \(G = \frac{1}{8π}\) instead of \(G = 1\), used in cosmology because it removes the \(8π\) term from the Einstein field equations. Similar to the axion-photon conversion probability, it is in natural units and \(B\) and \(L\) should be converted or the missing constants reinserted. The equation is valid in the chameleon-matter coupling range \(1 \leq β_m \leq \num{1e6}\), which corresponds to non-resonant production in the Sun as well as evades significant chameleon interaction with materials on its path to the CAST magnet.

Note that in many cases \(M_γ = \frac{m_{\text{pl}}}{β_γ}\) is introduced. If however, one defines its inverse, \(g_{βγ} = \frac{β_γ}{m_{\text{pl}}}\), eq. \eqref{eq:theory:chameleon_conversion_prob} takes the form of eq. \eqref{eq:theory:conversion_prob} for the axion, including an effective coupling constant of units \(\si{GeV⁻¹}\).

4.7. Current bounds on coupling constants

The field of axion searches is expanding rapidly in recent years, especially in haloscope experiments. A thorough overview of all the different possible ways to constrain axion couplings and best limits in each of them is not in the scope of this thesis. We will give a succinct overview of the general ideas and reference the best current limits on the relevant coupling constants in the regions of interest for this thesis. A great, frequently updated overview of the current best axion limits is maintained by Ciaran O'Hare, available here (O’HARE 2020). Generally, axion couplings can be probed by three main avenues:

Pure, indirect astrophysical constraints
Different astrophysical phenomena can be used to study and constrain axion couplings. One example is the cooling rate of stars. If axions were produced inside of stars and generally manage to leave the star without interaction, they carry energy away. Similar to neutrinos they would therefore contribute to stellar cooling. From observed cooling rates and knowledge of solar models, constraints can be set on potential axion contributions. Many other astrophysical sources can be probed in similar ways. In all cases these constraints are indirect in nature. Which coupling can be constrained depends on the physical processes considered.
Direct astrophysical constraints
Certain types of laboratory experiments attempt to measure axions directly and produce constraints due to non-detection. Solar helioscopes attempt to directly measure axions produced in the Sun; more on these in chapter 5. Haloscope experiments utilize microwave cavities in an attempt to tune to the frequency resonant with the axion mass of cold axions part of the dark matter halo. While relying on astrophysically produced axions, the intent is direct detection. Recent interest in axions also means data from WIMP experiments like the XENON collaboration (Aprile et al. 2017) is being analyzed for axion signatures. Haloscopes and helioscopes depend on the axion-photon coupling \(g_{aγ}\), while WIMP experiments may consider the axion-electron \(g_{ae}\) or axion-nucleon \(g_{aN}\) couplings. For the astrophysical production mechanism an additional dependency on the coupling producing the axion source is added, which may result in only being able to give constraints on products of different couplings.
Direct production constraints
The final approach is full laboratory experiments, which both attempt to first produce axions and then detect them. This idea is commonly done in so called 'light shining through the wall' (LSW) experiments like the ALPS experiment at DESY (Bähre et al. 2013). Here, a laser cavity in a magnetic field is intended as an axion production facility. These produced axions would leave the cavities, propagate through some kind of wall (e.g. lead) and enter a second set of equivalent cavities, just without an active laser. Produced axions could convert back into photons in the second set of cavities. The disadvantage is that one deals with the \(g_{aγ}\) coupling both in production and reconversion. 11

The bounds of interest for this thesis are the axion-photon coupling \(g_{aγ}\) for masses below around \(\SI{100}{meV}\) 12 and the product of the axion-photon and axion-electron couplings \(g_{ae}·g_{aγ}\) in similar mass ranges. The latter for the case of dominant axion-electron production \(g_{ae}\) in the Sun and detection via \(g_{aγ}\) in the CAST magnet.

The current best limit on the axion-photon coupling from laboratory experiments is from the CAST Nature paper in 2017 (Collaboration and others 2017), providing a bound of

\[ g_{aγ, \text{Nature}} < \SI{6.6e-11}{GeV^{−1}} \text{ at } \SI{95}{\%} \text{ CL}. \]

Similarly, for the product on the axion-electron and axion-photon couplings, the best limit from a helioscope experiment is also from CAST in 2013 (Barth et al. 2013). This limit is

\[ g_{ae}·g_{aγ} \lesssim \SI{8.1e-23}{GeV^{-1}} \text{ for } m_a \leq \SI{10}{meV} \]

and acts as the main comparison to the the results presented in this thesis.

From astrophysical processes the brightness of the tip of the red-giant branch (TRGB) stars is the most stringent way to restrict the axion-electron coupling \(g_{ae}\) alone. This is because axion production would induce more cooling, which would lead to a larger core mass at helium burning ignition, resulting in brighter TRGB stars. (Capozzi and Raffelt 2020) calculate a limit of \(g_{ae} = \num{1.3e-13} \text{ at } \SI{95}{\%} \text{ CL}\). A similar limit is obtained in (Straniero et al. 2020). (Bertolami et al. 2014) compute a comparable limit from the White Dwarf luminosity function. However, for purely astrophysical coupling constraints the strong assumptions needing to be made about the underlying physical processes imply these limits are by themselves not sufficient. In fact there is reason to believe at least some of the current astrohpysical bounds are overestimated. In (Dennis and Sakstein 2023) the authors use a machine learning (ML) model to predict the brightness of TRGB stars to allow for much faster simulations of the parameter space relevant for such bounds. Using Markov Chain Monte Carlo models based on the ML output, they show values up to \(g_{ae} = \num{5e-13}\) are not actually excluded if the full uncertainty of stellar parameters is included. As their calculations only went up to such values, even larger couplings may still be allowed.

Finally, observations from X-ray telescope can also be used for limits on the product of \(g_{ae}·g_{aγ}\). This has been done in (Dessert, Long, and Safdi 2019) based on data from the Suzaku mission and followed up on by the same authors in (Dessert, Long, and Safdi 2022) using data from the Chandra mission. In the latter they compute an exceptionally strong limit of

\[ g_{ae} · g_{aγ} < \SI{1.3e-25}{GeV^{-1}} \text{ at } \SI{95}{\%} \text{ CL}, \]

valid for axion masses below \(m_a \lesssim \SI{5e-6}{eV}\).

Chameleons
The current best bound on the chameleon-photon coupling was obtained using a previous GridPix detector with data taken at CAST in 2014/15. The observed limit obtained by C. Krieger in (Christoph Krieger 2018; Anastassopoulos et al. 2019) is \[ β_{γ, \text{Krieger}} = \num{5.74e10} \text{ at } \SI{95}{\%} \text{ CL}. \]

4.7.1. Potential solar axion hints

For completeness a few words about previous results which do not actually provide a limit, but rather show small hints of possible axion signals.

One of the first credible hints of such a solar axion signal comes from the XMM Newton telescope. Seasonal variation in the X-ray flux at a level of \(11σ\) is observed. The explanation provided in (Fraser et al. 2014) is axion reconversion into X-rays in Earth's magnetic field. While criticism exists (Roncadelli and Tavecchio 2015), there is anyhow recent interest in this signal (Ge et al. 2022), this time due to axion quark nuggets (AQN). An AQN explanation would produce a signal up to \(\SI{100}{keV}\) outside the XMM Newton sensitive range. The authors propose to check archival data of the Nuclear Spectroscopic Telescope Array (NuSTAR) or the Gamma-Ray Burst Monitor of the FERMI telescope for such seasonal variations.

Further, while the CAST Nature result (Collaboration and others 2017) provides the current best limit on the axion-photon coupling, its axion candidate dataset actually shows a signal excess at \SI{3}{keV} at a \(3σ\) level. Statistical effects or possibly not perfectly accounted for systematic variation resulting in slightly more argon fluorescence during tracking data is more likely.

In 2020 the XENON collaboration announced having seen an excess in their electron recoil counts at energies compatible with a solar axion signal (Aprile et al. 2020). A possible explanation not requiring solar axions was given as trace amounts of tritium below their sensitivity threshold. This garnered a lot of attention, because while only of \(3.4σ\) significance it was the first hint of a potential solar axion signal published by a large collaboration as such. However, combining the resulting axion coupling with astrophysical results indicates a more likely non axion origin for the signal (Athron et al. 2021; cite: Luzio, Fedele, et al. 2020). With the release the first XENONnT results in 2022 (Aprile et al. 2022) about new physics in which no excess is visible, the old result is ruled out. Add to that, the LUX-ZEPLIN collaboration, a similar xenon filled experiment, recently published their results on new physics. Although not as sensitive as XENONnT, but more sensitive than XENON1T, no excess was observed there either (Aalbers et al. 2023).

4.7.2. A few more words on haloscopes    extended

A haloscope is a type of axion experiment consisting of a (typically microwave) cavity placed in a magnetic field. It intends to detect axions of the dark matter halo of our galaxy. Axions that are part of the dark matter component are necessarily very low energy as they decoupled long ago and underwent cooling ever since. Thus, their energies are in the microwave range. If a cavity has a resonance frequency matching the axion mass (the kinetic energy is negligible, so the majority of the energy is in the mass), the conversion probability is enhanced by the quality factor \(Q\) of the cavity (effectively the number of reflections in the cavity). The upside of such experiments are the strong enhancements possible, which allow to reach very low coupling constants. However, a cavity has a single resonance frequency, limiting the mass range to be studied to a very narrow range. Most experiments use cavities that can be tuned to expand the mass range. At each tuned frequency data is taken for a fixed amount of time to reach a certain coupling constant. As such a tunable cavity experiment can scan a narrow band of axion masses over the course of its data taking campaign. Due to the simplicity of the setup these type of experiments are very popular nowadays.

Footnotes:

1

\(C\) refers to the discrete transformation of charge conjugation and \(P\) for parity transformation. Both refer to the idea of studying a physical system with either (or both) of these transformations applied. A \cp conserving theory (or system) would behave exactly the same under the combined transformation. The Standard Model is mathematically \cpt invariant (\(T\) being time reversal). As such, if a system exhibits different behavior under time reversal it implies a violation of \cp to achieve a combined \cpt invariance.

2

At least for me.

3

\(\mathrm{U}(1)\) refers to the "circle group", i.e. the group that describes rotations on a unit circle (consider a phase shift on the complex plane). The group operation as such can be considered as multiplication by a complex phase. \(\mathrm{SU}(n)\) is the special unitary group, which means the group of unitary matrices of rank \(n\) with determinant 1, where the group operation is matrix multiplication of these matrices (for \(\mathrm{SU}(2)\) the Pauli matrices multiplied by \(\frac{i}{2}\) are a possible set of infinitesimal generators for example).

4

The parameter \(θ\) as written in the equation is actually a compound of a pure \(θ\) from the QCD \(θ\) vacuum and an electroweak contribution. A chiral transformation is required to go to the physical mass eigenstates to diagonalize the quark mass matrix. This adds a term to \(θ\), \(\overline{θ} = θ + \arg \det M\). We simply drop the bar over \(\overline{θ}\).

5

For all practical purposes the terms interaction eigenstates and propagation eigenstates (the latter often also called mass eigenstates) are a convenient tool to work with. A field \(X\) in the interaction eigenstate corresponds to the field we can actually measure in an experiment. However, if fields interact, say with another field \(Y\), then along their space and time evolution they may mix. Therefore, it is useful (and convenient) to introduce a propagation eigenstate \(X'\) in which that new field (which is different from the physical field \(X\)!) propagates without any interaction. The nature of the field interactions have been absorbed into the time and space evolution of the field itself – it is a superposition of the \(X\) and \(Y\) interaction eigenstates. In the simplest case a field \(X'\) in the propagation state may just be oscillating between \(X\) and \(Y\), for example.

6

There is a more detailed derivation written by Biljana Lakić and Krešimir Jakovčić available internally in the IAXO collaboration. If you don't have access to it, reach out!

7

This can be done easily by treating the buffer gas as a refractive medium with complex refractive index \(n_γ = 1 - m²_γ / (2ω²) - iΓ/2ω\), which then produces attenuation via the \(Γ\) attenuation term. This is useful to know as it relates to the X-ray properties as discussed later in sec. 6.1.1 and sec. 6.1.2.

8

http://phd.vindaar.de/docs/bufferGasIAXO/v1/index.html contains the axionMass.pdf, axionMass.nim and finally axionMass.org, the document from which both other files are generated.

9

The Primakoff effect – named after Henry Primakoff – is the production of neutral pions in the presence of an atomic nucleus. Due to the pseudoscalar nature and coupling to photons of both neutral pions and axions the equivalent process is allowed for axions.

10

The difference of the peak position of the blackbody spectrum and the Primakoff flux is due to the axion production being dominantly from \(\SIrange{7.5}{17.5}{\%}\) of the solar radius (compare fig. #fig:theory:solar_axion_flux:flux_radial) where temperatures on average are closer to \(\sim\SI{12e6}{K}\).

11

While astrophysical sources of course also introduce an additional \(g²\) from their production (resulting in all experiments depending on \(g⁴\) effectively), the advantage is that in absolute terms an astrophysical source produces orders of magnitude more axion flux than a LSW experiment. In that sense LSW experiments deal with a squared suppression over those depending on astrophysical sources.

12

\(\SI{100}{meV}\) is the range in which the CAST experiment is not in the fully coherent regime of the \(\sinc\) term of the conversion probability anymore, resulting in significant sensitivity loss.

Click on any heading marked 'extended' to open it