Looking at the axionic dark sector with ANITA

Abstract

The ANITA experiment has recently observed two anomalous events emerging from well below the horizon. Even though they are consistent with tau cascades, a high-energy Standard Model or Beyond the Standard Model explanation is challenging and in tension with other experiments. We study under which conditions the reflection of generic radio pulses can reproduce these signals. Furthermore, we propose that these pulses can be resonantly produced in the ionosphere via axion–photon conversion. This naturally explains the direction and polarization of the events and avoids other experimental bounds.

1 Introduction

ANITA (ANtarctic Impulsive Transient Antenna) is a flying radio antenna dedicated to measuring impulsive radio signals in the Antarctica [1,2,3]. In particular, it can trigger pulses originated by cosmic ray air showers [3]. ANITA has a very good angular resolution and is able to discern whether the events are direct or reflected in the ice by measuring the polarization and phase (so-called polarity by the ANITA collaboration) of the radio pulse. Two of the direct cosmic ray events observed in the first and third flights, which seem to be originated well below the horizon ($27^\circ$ and $35^\circ$ respectively) [4, 5], are particularly intriguing and cannot in principle be interpreted as caused by high-energy cosmic rays. The only standard model (SM) particle that could potentially traverse a large amount of Earth matter (in this case around 6000 and 7000 km) and initiate a particle cascade leading to these events is a very high-energy ($\mathcal {O}(\mathrm{EeV})$) neutrino. However, for these extremely high energies, the neutrino-nucleon cross section leads to a very small survival probability ($\lesssim 10^{-6}$) over the chord length of the events, rendering such interpretation strongly disfavored [6, 7].

Recently, two potential SM explanations have been proposed in terms of transition radiation [8] and reflection on anomalous sub-surface structures [9]. In both cases, the origin of the anomalous events would be a reflected cosmic ray shower. Unless reflection occurs on a rather tilted surface, this hypothesis is in principle in 2.5$\sigma$ tension with the observed polarization angle of the first event [4]. Interestingly, both explanations predict particular signatures. On the one hand, transition radiation predicts that events with large elevations will be anomalous, in slight tension with current data. On the other hand, sub-surface structures predict some amount of double events, which should also be generated by the calibration pulses emitted by the HiCal antenna [10, 11]. Therefore, dedicated searches and/or more exposure are required to validate these possibilities.

Several Beyond the SM (BSM) scenarios have also been proposed to explain the origin of these events in terms of high-energy particles [7, 12,13,14,15,16,17,18,19,20]. However, they are rather in tension with IceCube and Auger bounds [20].

In this Letter, we propose a novel origin for these intriguing signals. We will show that reflected radio waves tend to present the properties of the mysterious ANITA events. Furthermore, we propose that the radio signal is generated via the conversion of an axion-like pulse. For the masses suggested by the data, this transformation happens to be resonant in the Earth ionosphere. This process involves very soft ($\mathcal {O}\left( \mu \mathrm {eV}\right)$) physics, invisible to IceCube and Auger, that ANITA can potentially test with a dedicated analysis.

2 Anomalous ANITA events

Atmospheric cosmic ray showers produce radio pulses with linear polarization perpendicular to the Earth magnetic field $\mathbf {B}^\oplus$, and a well-determined phase that flips at reflection. Since the magnetic field in the Antarctica is mostly vertical, ANITA searches for high-energy cosmic ray showers looking for horizontally polarized radio signals. In particular, the anomalous events are mostly horizontally polarized, and their phase led ANITA to interpret them as generated by up-going cosmic rays.

Any reflected electromagnetic wave, though, also tends to be horizontally polarized. In addition, if its origin is not a high-energy cosmic ray shower, its phase depends on the production mechanism and can thus match the one of the anomalous events. In this section, we will thus explore generic down-going radio waves reflected in the Antarctic ice as the origin for these events.

Fig. 1

Expected angular distribution of linearly polarized events perpendicular to $\mathbf {B}^\oplus$ under different hypotheses as labelled in the legend (see the main text for further details). More information about the computation of the reflected radio wave expected distribution (solid lines) is given in the Supplemental Material. The two ANITA anomalous events are shown in gray. The black line marks the elevation of the horizon, as seen from the ANITA balloon. The elevation range plotted corresponds to ANITA’s angular acceptance [2]

This hypothesis is illustrated in Fig. 1, where we show in solid the expected angular distribution of reflected events perpendicular to $\mathbf {B}^\oplus$ as a function of their elevation $\varepsilon$. We have assumed an incident isotropic flux, linearly polarized with random polarization angles and a given degree of polarization $P_i$. The reflected flux peaks at elevations where light reflects close to the Brewster angle $\theta _B \sim 53^\circ$ (corresponding to $\varepsilon \sim -37^\circ$),¹ defined as the angle at which the reflected signal is polarized exactly in the horizontal direction. The elevations of the observed events, within $1 \sigma$, are shown in gray; they are both close to the peak.

In Fig. 1 we also show the expected fluxes associated to other relevant hypotheses for the origin of the anomalous events. First, a SM tau neutrino flux (dotted line) which strongly peaks at the horizon and is therefore highly disfavored. Next, a generic BSM high-energy particle with a nucleon interaction cross section 10 times weaker than that of the SM neutrino (dot-dashed). The latter hypothesis partially alleviates the tension in the angular distribution, but tension with IceCube and Auger data remains [20]. Another possibility is that ANITA misidentified reflected events originated by ultra-high-energy cosmic ray (UHECR) air showers, classifying them instead as direct events. This hypothesis, shown by the dashed line, is disfavored since it would require phase misidentification by ANITA, which is excluded at $\sim 4.5\sigma$ [21]. However, alternative scenarios have been recently proposed to support this possibility [8, 9].

According to Fig. 1, our hypothesis of reflected radio waves, linearly polarized in random initial directions, explains the up-going direction of the two events. The other alternative scenarios shown in the figure predict more events close to the horizon, and so they could be discriminated from our proposal as ANITA accumulates more exposure.

Fig. 2

Observed (orange) and incident (red) polarization angles. We also show the projection of the Earth magnetic field [22, 23], in the reflected (light green) and incident (dark green) polarization planes. The shaded regions correspond to $1\sigma$ uncertainties assuming a $4.6^\circ$ uncertainty in the determination of the polarization direction [5], and a $2^\circ$ uncertainty in the orientation of the Earth magnetic field [21]

Another key observable in ANITA, not shown in Fig. 1, is the polarization angle. In Fig. 2 we quantify under which initial conditions our proposed hypothesis is able to reproduce not only the angular distribution, but also the observed polarization angle of the anomalous events. We show the reflected (orange) and incident (red) polarization angles for both events assuming that the signal is fully polarized. For the signal to be identified as an UHECR, the reflected (orange) electric field should be orthogonal to the Earth magnetic field (green), i.e, it should essentially be along the horizontal direction (H) as it is shown in the figure.

Since both events emerge at elevations close to the Brewster angle $\theta _B$, the vertical (V) component of the reflected electric field becomes quite suppressed. This is the reason why the uncertainty on the incident polarization angle is larger than the uncertainty on the reflected polarization angle.² This effect is particularly relevant for the second event, since it is closer to $\theta _B$ and has a reflected vertical component compatible with 0. The first event, on the other hand, has a non-zero vertical component, which significantly tightens the range of allowed incident polarization angles.

3 Axion-like origin

An incoming isotropic flux of linearly polarized radio waves with the required initial conditions to reproduce the observed ANITA signals, though, cannot in principle be explained in the SM. Below, we will construct a BSM explanation based on the following requirements:

• The source must generate a flux of impulsive radio signals, spatially isolated with a linear polarization and phase consistent with Fig. 2.

• Due to the tension with IceCube and Auger data, the production process must not involve high-energy particle cascades.

The first requirement comes from the ANITA triggering system [21, 25, 26], that requires a source of isolated impulsive signals. Notice that, even though it is not a necessary requirement, astrophysical sources are expected to be isotropically distributed, as assumed in Sect. 2.

The second requirement calls for a new physics mechanism able to coherently generate electromagnetic waves. This phenomenon should produce waves with frequencies $\sim \mathcal {O}(1 \,\mathrm {GHz})$, i.e., it can be associated with very low energies $\sim \mathcal {O}(10^{-7}\,\mathrm {eV})$. An archetypal BSM example are axion-like particles (ALPs), that convert into photons in the presence of an external electromagnetic field.

The ALPs, first proposed to solve the Strong CP problem [27,28,29], arise in many extensions of the SM and can constitute the dark matter in our universe [30,31,32]. Furthermore, they present self-interactions that produce a very rich and complex phenomenology. In particular, different phenomena, like condensation into a bosonic soliton or instabilities leading to scalar field bursts, may produce impulsive, spatially localized configurations of the scalar field [33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48]. If an axion pulse, with a macroscopic occupation number, reaches us, it can transform into the electromagnetic pulses observed by ANITA via the interaction with the Earth magnetic field.

For this phenomenon to explain the ANITA anomalous events, its rate should be $\sim \mathrm {month}^{-1}$. Any calculation of such rate is highly model dependent, as the phenomenology of non-linear ALP interactions is complex and still under study. Nevertheless, it is estimated that in a local neighbourhood $\sim 1 \, \mathrm {pc}^3$ there can be between $10^{10}$ and $10^4$ ALP overdensities [49]. Each overdensity would develop unstable bosonic solitons within time scales $\sim 10^{-2}$–$10^7$ years [41, 50]. Thus, the required rate to explain the anomalous events could plausibly be attained.

An ALP is a pseudo-scalar field, a, that interacts with photons via a Lagrangian density $\frac{1}{4} g_{a \gamma \gamma } a F_{\mu \nu } \tilde{F}^{\mu \nu }$, where $F_{\mu \nu }$ is the electromagnetic field tensor and $\tilde{F}_{\mu \nu }$ its dual. In presence of an external magnetic field $\mathbf {B}^\oplus$, the classical equations of motion for the axion a and electric field $\mathbf {E}$ in a plasma with free electron density $n_e$ are given by [51, 52]

$$\begin{aligned} \left[ \partial _z^2\! +\! \omega ^2\! +\! \begin{pmatrix} -\omega _p^2 &{}\quad -i \omega _p^2 \frac{\varOmega _z}{\omega } &{}\quad 0 \\ i \omega _p^2 \frac{\varOmega _z}{\omega } &{}\quad - \omega _p^2 &{}\quad - g_{a \gamma \gamma } B^\oplus _y \omega ^2 \\ 0 &{}\quad -g_{a \gamma \gamma } B^\oplus _y &{}\quad -m_a^2 \end{pmatrix}\right] \! \begin{pmatrix} E_x \\ E_y \\ a \end{pmatrix}\!\! =\! 0 \, , \end{aligned}$$

(1)

where we have assumed waves propagating in the Z direction, a magnetic field in the YZ plane, and we have taken the Fourier transform in time. Here $\omega _p = \sqrt{\frac{e^2}{m_e} n_e}$ is the plasma frequency, $\varOmega _z = \frac{e B^\oplus _z}{m_e}$ is the cyclotron frequency of the plasma, $m_e$ is the electron mass, and e is the electron charge in natural units. Notice that the axion mass $m_a$ in the above equations can include extra contributions depending on the particular scalar self interactions considered.

When the Faraday rotation effects (non-diagonal terms proportional to $\varOmega _z$) are switched off, the solution of these equations gives an electromagnetic wave linearly polarized parallel to the external magnetic field. The Faraday rotation then leads to the rotation of this polarization in the XY plane.

4 Resonance in the ionosphere

Equation (1) is a system of coupled wave equations with two characteristic frequencies, $\omega ^2 - \omega _p^2$ and $\omega ^2 - m_a^2$. When both frequencies are equal, axions convert resonantly into photons [51,52,53]. The smallest observed frequency of the anomalous events, $\omega _\mathrm {min} \sim 0.25 \, \mathrm {GHz}$, requires $m_a \lesssim \omega _\mathrm {min} \lesssim 10^{-7} \, \mathrm {eV}$. These masses happen to be in the range of typical $\omega _p$ values of the Earth ionosphere [54], and thus an incoming axion pulse with $m_a \lesssim 10^{-7} \,\mathrm {eV}$ would cross a region where $m_a \simeq \omega _p$, resonantly transforming into electromagnetic waves. This resonant conversion takes place in a relatively narrow region $\mathcal {O}(10\, \mathrm {km})$. The produced radio pulse will initially be polarized parallel to the Earth magnetic field, with a phase that depends on the axion wave phase and the sign of the axion–photon coupling $g_{a\gamma \gamma }$. This signal will later traverse the rest of the ionosphere, rotating its polarization vector due to the Faraday effect.

Generically, the axion burst will cross the resonant region twice, producing two radio pulses. Since the propagation through the ionosphere partially unpolarizes and decoheres the radio pulse generated during the resonance, there is a trade-off between having enough ionosphere to generate the Faraday rotation required to match the signals (see Fig. 2) and keeping the pulse coherent.³ These effects are mainly controlled by the density of free electrons $n_e$ in the ionosphere, which fluctuates in time in more than one order of magnitude [54]. Therefore, we can distinguish different phenomenological scenarios depending on the values of $n_e$. For small values, the pulse produced in the second resonant region will not undergo enough Faraday rotation, producing a mostly vertically polarized (parallel to $\mathbf {B}^\oplus$) signal. This pulse would be triggered out by the ANITA analysis or strongly suppressed due to the reflection close to $\theta _B$. However, the pulse generated during the first resonance would potentially be detectable. For high values of $n_e$, the first resonantly produced pulse would become incoherent, but the one generated in the second resonance would still be coherent and also experience enough Faraday rotation to pass ANITA’s triggers.

To illustrate these effects, in Fig. 3 we show the numerical solution of Eq. (1) for the second resonant conversion, which takes place at the lowest altitudes (details on how Eq. 1 is numerically solved can be found in Appendix C).

The top panel shows the plasma frequency $\omega _p$ and the cyclotron frequency $\varOmega _z$⁴ as a function of the propagated distance in the ionosphere, together with the axion mass $m_a$ value considered; the resonance occurs when both lines coincide. The central panel shows the total squared amplitude of the electric field due to resonant axion–photon conversion for different frequencies, normalized to the vacuum squared amplitude $|\mathcal {A}_\mathrm {vac}|^2=\left( \frac{2 a_0 g_{a \gamma \gamma } B^\oplus _y \omega ^2}{m_a^2}\right) ^2$, where $a_0$ is the amplitude of the incoming axion field. Finally, the bottom panel shows the squared projection of the electric field in the direction perpendicular to the Earth magnetic field, generated due to Faraday rotation.

Fig. 3

Top panel: $\omega _p$ and $\varOmega _z$ profile in the ionosphere, assuming a Chapman layer profile [54, 55] with a plausible maximum free electron density $n_e^\mathrm {max} = 2 \times 10^6 \, \mathrm {cm}^{-3}$, along with the axion mass chosen in the simulation. Central panel: electric field squared amplitude, normalized to the vacuum axion–photon conversion squared amplitude. In dashed, we show the result using the analytical approximation given by Eq. (2). The dashed gray line corresponds to the vacuum conversion squared amplitude. In the shaded gray region, we have switched off the axion–photon coupling in order to show only the propagation of the second resonant burst. Bottom panel: squared component of the electric field perpendicular to the Earth magnetic field, generated via Faraday rotation. We have normalized it to the total squared amplitude as given by Eq. (2)

The enhancement factor with respect to the vacuum axion–photon transition observed in Fig. 3 can be qualitatively understood using the WKB and stationary phase approximations [56] to solve Eq. (1)

$$\begin{aligned} |E|^2 = |a_0|^2\left( \frac{2 g_{a \gamma \gamma } B^\oplus _y \omega ^2}{m_a^2}\right) ^2 \left[ \frac{\pi }{4} \frac{m_a^2/k^2}{\left. \frac{\mathrm {d} \omega _p^2/m_a^2}{k \, \mathrm {d} z} \right| _\mathrm {res}} \right] \, , \end{aligned}$$

(2)

where E is the electric field after the resonance, $k=\sqrt{\omega ^2-m_a^2}$ is the wave number, and the derivative is evaluated at the resonance. The term inside the square brackets gives the enhancement due to the resonance in the plasma. This approximate solution is shown by the dashed lines in Fig. 3. In the ionosphere, $\omega _p$ varies over distances $\mathcal {O}(10\, \mathrm {km})$, much longer than the wavelengths that can be observed by ANITA, $\mathcal {O}(\mathrm {m})$. Therefore, the denominator in Eq. (2) is rather small, leading to the $\mathcal {O}(10^2{-}10^3)$ global enhancement observed in Fig. 3. As can be observed in the central panel of Fig. 3, before reaching the resonance the squared amplitude of the electric field is basically equal to its value when the axion–photon conversion takes place in vacuum. That is, outside the resonant region, the axion–photon conversion rate is essentially the one in vacuum, different from zero but negligible compared with the resonant value.

The projection of the electric field shown in the bottom panels of Fig. 3 is directly related to the polarization angle of the radio pulse. When the pulse leaves the ionosphere, this angle $\psi _i$ must be consistent with the red regions in Fig. 2 for the ANITA anomalous events to be reproduced. In particular, for $\omega =1 \,\mathrm {GHz}$ we obtain a polarization angle of $\psi _i \approx 41^\circ$, in agreement with Fig. 2 and thus consistent with ANITA-I and ANITA-III. For $\omega =1.5\, \mathrm {GHz}$ we have $\psi _i \approx 74^\circ$, consistent with ANITA-III but in disagreement with ANITA-I. This is because in this case the Faraday rotation effect does not generate enough horizontal (essentially orthogonal to $\mathbf {B}^\oplus$) component of the electric field to reproduce the first event. A different ionospheric profile and/or $m_a$, though, would in general give a different result. Thus, pulses experiencing different levels of Faraday rotation, with diverse polarization directions after the ionosphere, can be generated.

In summary, the ANITA anomalous events could be due to an axion burst that resonantly converts into photons in the Earth ionosphere matching the required conditions of the isotropic, linearly polarized flux previously discussed above and shown in Figs. 1 and 2. Our proposal is schematically summarized in Fig. 4.

Fig. 4

Sketch of an axion burst arriving to the ionosphere, undergoing resonant conversion, Faraday rotation, and reflecting on the ice surface before reaching ANITA

5 Spectral properties and ALP scenario

The characteristics of the observed ANITA events can also be used to extract properties of the ALP burst.

On the one hand, we can infer information about the frequency of the burst. Both anomalous events show a large correlation with a cosmic-ray template, and [5] shows the Amplitude Spectral Density of the ANITA-III anomalous event. We have checked that both spectral requirements can be satisfied within experimental uncertainties with a Gaussian pulse with central frequency $\omega \lesssim 2.5\, \mathrm {GHz}$ and width $\sigma _\omega = 1.5-4\, \mathrm {GHz}$.

Fig. 5

In gray, estimated axion mass $m_a$ and axion–photon coupling $g_{a \gamma \gamma }$ consistent with the ANITA events for different densities $\rho$ of the incoming axion burst. The minimum observed frequency $\omega \gtrsim 0.25 \, \mathrm {GHz}$ forces $m_a \lesssim 1.6 \cdot 10^{-7} \, \mathrm {eV}$. We also show current experimental constraints [57,58,59,60,61,62]. The yellow region is compatible with the QCD axion models [63]

On the other hand, we can also extract information on the mass $m_a$ and coupling $g_{a \gamma \gamma }$ of the ALP. Using Eq. (2), we can estimate the relation between the observed electric field amplitude $\sim 1 \,\mathrm {mV/m}$, $m_a$, $g_{a \gamma \gamma }$, and the amplitude of the axion field burst $a_0$. The latter can, in turn, be determined by the energy density of the burst $\rho \sim |a_0|^2 \omega ^2$. Considering $B^\oplus \sim 0.45 \, \mathrm {G}$ [22, 23], a wave frequency $\omega \sim 1.5\,\mathrm {GHz}$, attenuation due to reflection on the Antarctic ice, and typical resonance-enhancement factors $\frac{\mathrm {d}\omega _p^2/m_a^2}{\mathrm {d}z} \sim 10^{-2} \, \mathrm {km}^{-1}$, we show in Fig. 5 an estimation for the values of $g_{a \gamma \gamma }$ and $m_a$ consistent with ANITA for different energy densities $\rho$ of the incoming axion burst. The gray region shows the mass-coupling range compatible with the minimum observed frequency $\omega \sim 0.25 \, \mathrm {GHz}$  (corresponding to $m_a \lesssim 10^{-7} \, \mathrm {eV}$), and the dashed lines label different densities of the incoming axionic burst. For densities $\gtrsim 10^{-12} \, \mathrm {g}/\mathrm {cm}^2$ all present experimental bounds are evaded [57,58,59,60,61,62].

6 Conclusions

In this work, we have explored the directional and polarization properties of the anomalous ANITA events 3985267 (ANITA-I) and 15717147 (ANITA-III). We have found that the reflection of an isotropic flux, linearly polarized in arbitrary directions, can naturally accommodate both observables. This is mostly due to the triggering of ANITA, that favors horizontally polarized events, together with reflection close to the Brewster angle.

Requiring a polarized flux not produced via high-energy cascades in order to avoid the IceCube and Auger bounds, we have proposed a generation mechanism based on the axion–photon conversion in the Earth magnetic field of a classical, high occupation number, axion burst.

Interestingly, we have also found that this conversion is dominated by a resonance naturally occurring in the ionosphere for the radio frequencies observed by ANITA. After traversing the remaining part of the ionosphere, the signal will undergo Faraday rotation, providing pulses polarized in different directions that can explain the mysterious events.

Our proposal can already be tested reanalyzing the data collected by ANITA, including the fourth flight data currently under analysis. If the hypothesis presented in the second section is correct, relaxing the triggering that requires geomagnetically correlated events should reveal events with a non-suppressed vertical component emerging from angles different from $\theta _B$. On the other hand, an axion-like origin generically predicts two consecutive events with different polarizations and/or coherence, since the axion burst experiences two resonant transitions into photons along its propagation through the ionosphere. Extra signals could thus be observed by searching for doubled events which may require decreasing the coherence threshold.

Appendices

Appendix A: Expected elevation distribution of reflected events

The ANITA anomalous events are compatible with being 100% linearly polarized perpendicular to the Earth magnetic field [4, 5]. Therefore, in Fig. 1 we have computed the expected distribution of events that satisfy these properties within experimental uncertainties.

To do so, we have considered a distribution of incoming radio pulses, characterized by their incident angle $\theta _i$, polarization angle $\psi _i$ and degree of polarization $P_i$. The pulses are assumed to be uniformly distributed in $\theta _i$ and $\psi _i$, with constant $P_i$. The Stokes formalism [24] relates the reflected and incident polarization states, and so it allows to obtain the expected distribution of events as a function of the corresponding reflected polarization angle, reflected degree of polarization, and angle of reflection.

The reflected radio pulses should match the polarization properties observed by the ANITA collaboration [4, 5]. They have to be, within experimental uncertainties, compatible with being observed as 100% linearly polarized perpendicular to the Earth magnetic field. Imposing this requirement, we obtain the expected distribution of events as a function of the reflected angle (which is directly related to the elevation $\varepsilon$). The resulting distribution is what we show in solid lines in Fig. 1.

In the computation, we have assumed an index of refraction for the Antarctic surface of $n=1.35$ [2, 10]. In addition, in order to impose the expected signal to be compatible with being observed as 100% linearly polarized perpendicular to the Earth magnetic field, information about the experimental uncertainties on the polarization angle $\psi$ and degree of polarization P is required. The former is reported by the ANITA collaboration to be 4.6$^\circ$ [5], whereas to our knowledge the latter is not available. In order to estimate this uncertainty, we have considered that the degree of polarization is given by

$$\begin{aligned} P = \frac{\sqrt{Q^2 + U^2 + V^2}}{I} \, , \end{aligned}$$

(A.1)

which is determined by the Stokes parameters of the incoming radio pulse,

$$\begin{aligned} I&= \frac{1}{N} \sum _{i=1}^N |\varepsilon _H^i|^2 + |\varepsilon _V^i|^2 \, , \end{aligned}$$

(A.2)

$$\begin{aligned} Q&= \frac{1}{N} \sum _{i=1}^N |\varepsilon _H^i|^2 - |\varepsilon _V^i|^2 \, , \end{aligned}$$

(A.3)

$$\begin{aligned} U&= \frac{1}{N} 2 \mathrm {Re} \left[ \sum _{i=1}^N \varepsilon _H^i \left( \varepsilon _V^i\right) ^* \right] \, , \end{aligned}$$

(A.4)

$$\begin{aligned} V&= \frac{1}{N} 2 \mathrm {Im} \left[ \sum _{i=1}^N \varepsilon _H^i \left( \varepsilon _V^i\right) ^* \right] \, , \end{aligned}$$

(A.5)

where $\varepsilon \equiv E + i \hat{E}$, $\hat{E}$ is the Hilbert transform of the electric field E and the subindices H and V denote the horizontal and vertical components of the electric field at $N = 500$ discrete time instants [21].

According to Ref. [21], the uncertainty on the Stokes parameters (and so the uncertainty on P) is dominated by random Gaussian white noise. We have thus assumed $\varepsilon ^i$ to be uncorrelated Gaussian random variables with a mean value given by the true electric field and a standard deviation $\delta$. The experimental uncertainties on the Stokes parameters (and therefore on P) have been estimated by randomly generating $\varepsilon ^i$. The only free parameter, $\delta$, has been adjusted to reproduce ANITA’s uncertainty on the polarization angle reported in [5].

Our estimation leads to an uncertainty on the degree of polarization of $\sim 0.13$. In any case, we have checked that considering a different input for this uncertainty does not change significantly our conclusions.

Finally, for the SM $\nu _\tau$ hypothesis and the flux coming from a generic particle with an interaction cross section 10 times smaller than the SM $\nu _\tau$, we have simulated propagation through the Earth with the numerical library $\nu$-SQuIDS [64,65,66].

Appendix B: Initial conditions for a partially polarized pulse

The incident pulse does not need to be fully polarized. In the following, we generalize Fig. 2 by relaxing this hypothesis. To visualize its effect, in the top and bottom panels of Fig. 6 we show the relation among the reflected degree of polarization and polarization angle, $P_r$ and $\psi _r$; and the incident degree of polarization and polarization angle, $P_i$ and $\psi _i$. The corresponding $1\sigma$ allowed region for $\psi _r$ extracted from [5] is shown in gray in the bottom panels. The gray region in the top panels corresponds to our estimation for the degree of polarization’s allowed region at $1\sigma$.

Fig. 6

Reflected polarization angle $\psi _r$ and degree of polarization $P_r$ as a function of the incident polarization angle $\psi _i$ and degree of polarization $P_i$ (labeled by the color). The $1\sigma$ allowed region for the polarization angle [5] (bottom panels) and our corresponding estimation for the degree of polarization (top panels) are shown in gray. All angles are measured with respect to the horizontal

Comparing Fig. 2 to the bottom panels in Fig. 6, we conclude that the range of allowed incident polarization angles increases once we slightly relax the assumption of a fully polarized signal. This is because the unpolarized part of the signal leads to a horizontally polarized component after reflection, which tends to tilt the polarization angle closer to the horizontal, i.e., closer to being perpendicular to $\mathbf {B}^\oplus$. Comparing the ANITA-I and ANITA-III panels, we notice that the level of initial polarization and range of values of the incident polarization angle is less stringent for ANITA-III, as expected since its incident angle is closer to $\theta _B$.

Appendix C: Numerical integration

Equation (1) has fast oscillating solutions with very different time-scales: the frequency of the waves, $\omega$, is much larger than all the other scales in the system such as the plasma frequency, $\omega _p$; the axion mass, $m_a$; and the axion–photon coupling term, $g_{a \gamma \gamma } B^\oplus _y \omega ^2$. For this reason, the numerical integration is not straightforward: we have first written all fields C(z) as $C(z) = \tilde{C}(z) e^{i k z}$, with $k = \sqrt{\omega ^2 - m_a^2}$. In this way, the fast oscillations driven by $\omega$ are effectively separated from the axion–photon conversion and Faraday rotation effects, which take place with much longer characteristic times.

In our computation, the initial condition for Eq. (1) is a burst of axion field and no electromagnetic field. Given the smallness of $g_{a\gamma \gamma }$, the feedback of the generated electromagnetic wave on the axion field can be neglected to a good approximation. Thus, the axion field amplitude can be considered constant. This allows us to simplify the problem and solve Eq. (1) for $E_i$ ($i=\{x,y\}$) considering the axion field as a constant source.

Furthermore, any set of coupled second order differential equations can be written as a set of first order differential equations (at the cost of doubling the number of equations) and thus, separating the real and imaginary parts of the electric field $\tilde{E}_i=\mathfrak {R}(\tilde{E}_i) + i\mathfrak {I}(\tilde{E}_i)$ and denoting $\tilde{D}_i = \partial _z \tilde{E}_i$, we can rewrite Eq. (1) as

$$\begin{aligned}&\partial _z \begin{pmatrix} \mathfrak {R}(\tilde{D}_x) \\ \mathfrak {I}(\tilde{D}_x)\\ \mathfrak {R}(\tilde{D}_y) \\ \mathfrak {I}(\tilde{D}_y) \\ \mathfrak {R}(\tilde{E}_x)\\ \mathfrak {I}(\tilde{E}_x)\\ \mathfrak {R}(\tilde{E}_y)\\ \mathfrak {I}(\tilde{E}_y) \end{pmatrix} = \begin{pmatrix} 0&{}2k&{}0&{}0&{} \omega _p^2 - m_a^2 &{} 0 &{} 0 &{} -\omega _p^2 \frac{\varOmega _z}{\omega } \\ -2k&{}0&{}0&{}0&{}0&{} \omega _p^2 - m_a^2 &{} \omega _p^2 \frac{\varOmega _z}{\omega } &{} 0 \\ 0&{}0&{}0&{}2k&{} 0 &{} \omega _p^2\frac{\varOmega _z}{\omega } &{} \omega _p^2 - m_a^2 &{} 0 \\ 0&{}0&{}-2k&{}0&{} - \omega _p^2\frac{\varOmega _z}{\omega }&{} 0 &{} 0 &{} \omega _p^2 - m_a^2 \\ 1&{}0&{}0&{}0&{}0&{}0&{}0&{}0\\ 0&{}1&{}0&{}0&{}0&{}0&{}0&{}0\\ 0&{}0&{}1&{}0&{}0&{}0&{}0&{}0\\ 0&{}0&{}0&{}1&{}0&{}0&{}0&{}0 \end{pmatrix}\nonumber \\&\qquad \begin{pmatrix} \mathfrak {R}(\tilde{D}_x) \\ \mathfrak {I}(\tilde{D}_x\\ \mathfrak {R}(\tilde{D}_y) \\ \mathfrak {I}(\tilde{D}_y) \\ \mathfrak {R}(\tilde{E}_x\\ \mathfrak {I}(\tilde{E}_x\\ \mathfrak {R}(\tilde{E}_y\\ \mathfrak {I}(\tilde{E}_y) \end{pmatrix} + \begin{pmatrix} 0 \\ 0\\ g_{a \gamma \gamma } B^\oplus _y \omega ^2 \tilde{a} \\ 0 \\ 0\\ 0 \\ 0\\ 0 \end{pmatrix}. \end{aligned}$$

(C.6)

where a is the amplitude of the axion pulse arriving to the Earth. This differential equation can be easily solved with a standard numerical integrator. In particular, we have used the one from the scipy python package, which is essentially a wrapper of the FORTRAN library odepack [67]. We used the default integrator: an explicit Runge–Kutta method of order 5(4) featuring an adaptive step size, setting the precision goal to be unnoticeable in Fig. 3.