# AniMAIRE - A New Openly Available Tool for Calculating Atmospheric Ionising
Radiation Dose Rates and Single Event Effects During Anisotropic Conditions.

###### Abstract

AniMAIRE (Anisotropic Model for Atmospheric Ionising Radiation Effects) is a
new model and Python toolkit for calculating radiation dose rates experienced
by aircraft during anisotropic solar energetic particle events. AniMAIRE
expands the physics of the MAIRE+ model such that dose rate calculations can
be performed for anisotropic solar energetic particle conditions by supplying
a proton or alpha particle rigidity spectrum, a pitch angle distribution, and
the conditions of Earth’s magnetosphere. In this paper, we describe the
algorithm and top-level structure of AniMAIRE and showcase AniMAIRE’s
capabilities by analysing the dose rate maps that AniMAIRE produces when the
time-dependent spectra and pitch angle distribution for GLE71 are inputted. We
find that the dose rates AniMAIRE produces for the event fall between the dose
rates produced by the WASAVIES and CRAC:DOMO models. Dose rate maps that
evolve throughout the event are also shown, and it is found that each peak in
the inputted pitch angle distribution generates a dose rate hotspot in each of
the polar regions. AniMAIRE has been made available openly online so that it
can be downloaded and run freely on local machines and so that the space
weather community can easily contribute to it using Github forking.

FB FL CD ClDy Space Weather

Surrey Space Centre, University of Surrey, Guildford, U.K. CSDRadConsultancy

C. S. W<EMAIL_ADDRESS>
AniMAIRE can simulate atmospheric radiation and electronics effects from any
inputted anisotropic distribution of particles hitting Earth.

AniMAIRE predicts GLE71 dose rates between those of CRAC:DOMO and WASAVIES,
and event-integrated dose rates are found to require anisotropy.

AniMAIRE has been made openly available online, such that anyone can run and
experiment with it, or contribute to the model through forking.

## Plain Language Summary

Occasionally an eruption occurs on the surface of the Sun, ejecting high
energy particles into the solar system in what are known as solar particle
events. If this eruption of particles happens to hit Earth and have a
sufficiently high kinetic energy, the particles can penetrate Earth’s
magnetosphere and atmosphere, creating particle radiation showers that
increase radiation dose rates in Earth’s atmosphere. These situations are
known as Ground-Level Enhancements (GLEs), which occur once every few years on
average. GLEs can cause issues in aircraft electronics and cause passengers
and crew on airplanes to experience heightened radiation doses. In this paper,
we describe a new, openly available model we’ve created to calculate dose
rates in the atmosphere during GLEs, called AniMAIRE. Most models that exist
to date are only able to calculate radiation dose rates when all particles are
hitting Earth equally from all directions, however this approximation is
generally only reliable during the later stages of events, and isn’t accurate
for many GLEs. AniMAIRE has therefore been designed so that it doesn’t need to
rely on this approximation, and can calculate dose rates for situations where
it is supplied a direction-dependent particle flux. In this paper we also show
and analyse the dose rate maps AniMAIRE produces when running across a GLE in
May 2012, finding several interesting results.

##  1 Introduction

Atmospheric radiation dose rates and Single Event Effects (SEEs) induced by
Ground Level Enhancements (GLEs) are one of the numerous hazards that space
weather can create on Earth [Shea Smart (2012), Meier . (2020), Dyer . (1989),
Dyer . (1990), Tobiska . (2015)]. GLEs occur on average once every few years
[Karapetyan (2008), Usoskin . (2015)] and can cause large increases in
radiation dose rates versus normal ‘quiet’ cosmic-ray-only conditions [Dyer .
(2007), Copeland . (2008)]; for example, increases by factors of 1000 or more
versus quiet conditions could be expected on a time scale of between 1 in 40
and 1 in 70 years [Dyer . (2017)]. The UK Royal Academy of Engineering
predicted that a large GLE could cause airline crew and passengers to receive
dose rates that exceed 20 times the recommended general public exposure limit
[Cannon . (2013)]. Radiation from these events could also cause electronics
issues in aircraft systems due to Single Event Effects (SEEs) [Normand Baker
(1993), Taber Normand (1993)], potentially increasing pilot workloads if
issues were to occur in multiple systems at once. It is self-evident that
lowering the risk of these issues occurring is important to both the airline
industry as a whole and airline passengers. Therefore, simulating the dose
rates and SEEs an aircraft will experience during a large space weather event
is an important part of both understanding the dose rates aircraft will
experience, and developing dose rate mitigation strategies.

The simulation of atmospheric radiation levels is currently an active area of
research, and many models have been developed over the last few decades to
model atmospheric radiation. Many models have historically focused on
simulating cosmic ray quiet time dose rates. However, recent modeling work in
the field has moved towards simulating GLEs, although typically with
assumptions of isotropy. In this paper, we discuss the development and release
of a new and openly available Python model for simulating anisotropic dose
rates due to incoming protons - called AniMAIRE. AniMAIRE is based on the
physics behind its precursor model, MAIRE+ [Hands . (2022)], but with the
additional feature that it can simulate dose rates during anisotropic
conditions.

The MAIRE+ model [Hands . (2022)] was developed by the Space Environment and
Protection Group at Surrey Space Centre in collaboration with the UK Met
Office as part of the UK’s SWIMMR programme. MAIRE+ is a successor model to
the previous QARM [Lei . (2004), Lei . (2005)], and MAIRE [Lei (), Hands .
(2016)] models, using many of the same particle responses that were present in
these. It is an isotropic model with input conditions tailored for the UK
airspace. Other major models that can currently perform calculations during
isotropic conditions include CARI-7 [Copeland (2017)], AVIDOS [Latocha .
(2009)] and NAIRAS [C. Mertens . (2009), CJ. Mertens . (2013)].

While isotropic models are more accurate during the later stages of most GLEs,
they are unable to accurately model dose rates during more anisotropic stages
of GLEs (typically the initial 30 minutes to 60 minutes of a GLE, although the
stages and features of a GLE vary considerably from event to event [Shea Smart
(2012)]). The only currently used major models that we are aware of that can
determine dose rates during anisotropic conditions are WASAVIES [Kataoka .
(2014), Sato . (2018)], SIGLE [Lantos Fuller (2003), Lantos Fuller (2004)] and
CRAC:DOMO by mishev2015computation. SIGLE takes an empirical approach to
dealing with anisotropy, using the count rate increases at the nearest ground-
level neutron monitor to a particular flight to apply dose rate weighting for
anisotropy.

In contrast, WASAVIES takes a physics-based approach and effectively traces a
full simulation of solar particle injection profiles through the solar
heliosphere to Earth’s magnetosphere. WASAVIES then traces particle directions
through Earth’s magnetosphere using the Tsyganenko 1989 model [Tsyganenko
(1989)], and then converts the particle spectra hitting Earth’s atmosphere
into dose rates. It does all this using a combination of pre-calculated
response databases and real-time simulations. Fitting against experimental
measurements is used to determine the most appropriate injection profile for a
particular event.

The CRAC:DOMO model is similar to WASAVIES and also uses trajectory tracing
through Earth’s magnetosphere using the Tsyganenko 1989 model (implemented
through PLANETOCOSMICS [Desorgher (), Desorgher . (2006)]) to convert particle
distributions at Earth’s magnetosphere into dose rates. The key difference
between CRAC:DOMO and WASAVIES here is that CRAC:DOMO assumes nothing about
the heliosphere and instead spectra and pitch angle distributions at Earth’s
magnetosphere are directly inserted into the model, which are determined
through a separate fitting model using ground-level neutron monitor data.

Few atmospheric radiation models (isotropic or anisotropic) have been released
openly both in terms of running and for modification in general, which is
likely one of the reasons why so many models have been developed in general,
as every group in the field has so far had to create their own model to
perform research. One positive recent development has been the recent addition
of models to NASA CCMC [Community Coordinated Modelling Center ()] system, for
instance, NAIRAS [C. Mertens . (2023)], for real-time running, and the
development of ESA’s Network of Models (NoM) [European Space Agency ()] where
models can be run on-demand.

To aid the atmospheric radiation community, AniMAIRE has currently been
released openly [Davis . (2024)] on Github and PyPi, where it can be installed
with a single command using the most common Python package installer, ‘pip’
and run locally by any users with a moderate level of experience in Python.
Links to AniMAIRE and all of the software it depends on (including several
additional useful software packages that were developed as part of this work)
can be found in the Software and Data Availability section at the end of this
paper. One of the reasons why Python has become so popular in the last decade
in the scientific community as a whole is because of the ease at which outputs
from Python ‘modules’ such as AniMAIRE developed by users can be fed into each
other, and as AniMAIRE has been made open source it can be used in other
projects in the community as a dependency.

The source code of AniMAIRE is openly available for users who wish to
investigate the algorithm in more detail or for members of the community who
wish to make their own contributions and improvements to the software through
Github ‘forking’ and pull requests. We are also currently investigating adding
AniMAIRE to the CCMC [Community Coordinated Modelling Center ()] and NoM
[European Space Agency ()] so that users without Python experience can also
perform calculations with AniMAIRE. More details on how to access and install
AniMAIRE and all its dependencies can be found in the data availability
section at the end of this paper.

AniMAIRE therefore provides novel features to the community of being openly
accessible, able to perform anisotropic calculations, a modular software
design, and the ability to be run locally and as part of other future
software.

##  2 The AniMAIRE Physics/Algorithm

Like MAIRE+, AniMAIRE uses altitude-dependent yield functions to convert an
energetic particle spectrum to a radiation dose rate at a given location on
Earth. However, when the SEP spectrum is anisotropic, the spectrum which can
be ‘seen’ from a given location on Earth is specific to that location.
Therefore, AniMAIRE calculates the energetic particle spectrum, $J(P)$ for a
given location on Earth as

$J(P)=F(P)\times J_{\parallel}(P)\times G(\alpha(P))$ (1)

$J_{\parallel}(P)$ is the SEP proton spectrum arriving along the ‘axis of
symmetry’, i.e. the spectrum at the direction of maximum flux, as a function
of particle rigidity $P$, which is usually close to the Interplanetary
Magnetic Field (IMF) direction [Mishev Usoskin (2016), Mishev . (2021)]. The
anisotropy is accounted for through the pitch angle distribution,
$G(\alpha(P))$. The parameter $\alpha(P)$ is the pitch angle between the
asymptotic directions corresponding to the given location and the location of
the axis of symmetry. Since the asymptotic direction at a given location is a
function of the rigidity of the particle reaching that location, the pitch
angle is also a function of particle rigidity. Finally $F(P)$ is a filtering
function, with a value of 0 for forbidden values of $P$ and 1 for allowed
values (see Cooke1991 for an overview of asymptotic directions and cut-off
values).

$J_{\parallel}(P)$ can also be set to be the cosmic ray spectrum, which can be
calculated in AniMAIRE through the DLR-modified ISO model [Matthiä . (2013)]
as implemented in the CosRayModifiedISO Python module [Davis (20242)]. In this
cosmic ray-only case, $G(\alpha(P))$ is set to be isotropic, i.e. always equal
to 1.

The general AniMAIRE process for calculating dose rates and fluxes at a single
location is as follows:

  1. 1.

Asymptotic directions are calculated for rigidities between 0.1 GV and 1010.0
GV using MAGNETOCOSMICS [Desorgher (2006), Desorgher (2004)] through the
AsympDirsCalculator Python module [Davis (20241)]. Currently, only the
Tsyganenko 1989 model [Tsyganenko (1989)] can be used, but more models could
be added in future versions of AsympDirsCalculator.

  2. 2.

The pitch angle between each asymptotic direction and the axis of symmetry is
calculated using the dot product of the asymptotic direction with the incoming
direction of the solar particle event. The pitch angle can be calculated using
the equation
$\cos(\alpha)=\sin(\Lambda)\sin(\theta_{s})+\cos(\Lambda)\cos(\theta_{s})\cos(\Phi-\phi_{s})$
given by cramp1997october. The pitch angle distribution is assumed to be
symmetric about this particular direction.

  3. 3.

The full energetic particle spectrum is calculated using eq. 1.

  4. 4.

The full energetic particle spectrum is multiplied by the relevant yield
functions for each of the output quantities and integrated with respect to
rigidity. This is achieved through the atmosphericRadiationDoseAndFlux Python
package [Davis Lei (2024)].

A diagram of this procedure is shown in fig. 1.

A standard run of AniMAIRE takes as input any generic proton and/or alpha
particle spectrum, a list of latitude/longitude pairs, a list of altitudes,
pitch angle distributions, Kp index, the date and time for the run to
correspond to, and reference pitch angle coordinates (the location on Earth
where the pitch angle is 0.0, i.e. the direction in which the GLE is hitting
Earth from). Many of these input parameters have default values so that a user
doesn’t have to specify every parameter all the time - for instance, if only a
proton spectrum is specified, AniMAIRE will automatically calculate the Kp
index using the PySpaceWeather Python module [Bender ()] using the specified
date and time. Equally, if a date and time is not specified, AniMAIRE will use
the current date and time by default.

Spectra do not always have to be specified directly. AniMAIRE can calculate
dose rates from incoming cosmic rays as stated earlier, or through supplying
parameters for several standard GLE proton spectrum distributions. Examples of
using these additional methods can be found at AniMAIRE’s Github and PyPi
pages, and AniMAIRE could be relatively trivially expanded in the future by
members of the space weather community so that it can be directly run using
spectra from other models.

Figure 1: A diagram showing the general algorithm AniMAIRE uses, as well as
the dependency packages that were created as part of the AniMAIRE development
process.

##  3 Application to Cosmic Rays and GLE71

Links to all of the data and Jupyter Notebooks used to generate the plots in
this section can be found in the Software and Data Availability section at the
end of this paper.

To showcase the outputs that AniMAIRE can produce for anisotropic situations,
AniMAIRE was run across the time-dependent spectra and pitch angle
distributions for GLE71 reported in mishev2021application. GLE71 was a complex
event, and Mishev et al. found that it was best modeled using a double
Gaussian pitch angle distribution with a modified power law rigidity spectrum.
They also found that the resulting distribution agreed well with SEP
measurements from the Payload for Antimatter Matter Exploration and Light-
nuclei Astrophysics (PAMELA)[Bruno . (2018)] when integrated across the whole
event.

The spectra and pitch angle distribution for each time interval were fed into
AniMAIRE to produce dose rate maps at each point in time throughout GLE71. The
input rigidity spectra equation corresponding to the GLE used was:

$J_{\parallel}(P)=J_{0}P^{-(\gamma+\delta\gamma(P-1))}$ (2)

and the pitch angle distribution used was:

$G(\alpha(P))=\exp{\bigg{(}-\frac{\alpha^{2}}{\sigma_{1}^{2}}\bigg{)}}+B\exp{\bigg{(}-\frac{(\alpha-\alpha^{\prime})^{2}}{\sigma_{2}^{2}}\bigg{)}}$
(3)

Here, $P$ is the rigidity in units of GigaVolts. Note that the pitch angle
itself is a function of rigidity $P$ and is defined relative to the axis of
symmetry of the incoming SEP event, given with $\Phi$ and $\Lambda$. These
parameters, as well as $\gamma$, $\delta\gamma$, $\sigma_{1}$, $\sigma_{2}$
and $\alpha^{\prime}$ were given by Mishev et al. and were determined through
their fitting method - although $\gamma$ values in the table are listed as
negative and have to be multiplied by $-1$ to insert into eq. 2.

An example of the effect of multiplying the anisotropic spectra and pitch
angle distribution together is shown in fig. 2, which shows the spectra for
three different cases which impact the latitude of 65.02° and longitude of
25.5° in the first timestamp of GLE71. The isotropic spectrum is the raw
spectrum, $J_{\parallel}(P)$, from eq. 2, and the anisotropic spectra are the
spectra multiplied by the pitch angle distribution, as calculated by AniMAIRE,
i.e. $J(P)$ from eq. 1. The latitude and longitude used were chosen somewhat
arbitrarily for this comparison but are the coordinate of the commonly used
OULU neutron monitor.

Figure 2: The GLE spectra for several different cases that are incident upon a
latitude of 65.02° and a longitude of 25.5°, the location of the OULU neutron
monitor at the first timestamp for the event. The vertical anisotropic
spectrum is a direct multiplication of pitch angle values with the isotropic
spectra, which is why the anisotropic spectra are always less than the
isotropic spectrum here (the pitch angle distribution used here was always
less than 1, with a very conservative isotropic spectrum).

The anisotropic spectra are less than the isotropic spectrum at all points as
the value of the pitch angle distribution ($G(\alpha)$) is always less than
approximately 1 in this case. The anisotropic spectra are also different in
shape to the isotropic spectrum at lower rigidities, illustrating the complex
effect of the pitch angle distribution on the spectra a coordinate receives.

Dose rate maps for the first timestamp in GLE71 were generated for both cosmic
rays and anisotropic input spectra and are shown in fig. 3, and the effect of
anisotropy on the spatial distribution of dose rates can be seen in fig. 3(b)
versus the isotropic cosmic ray dose rates. Polar regions of heightened dose
rates are the same approximate shape and location as in an isotropic situation
but contain a complex internal geographic structure. The double Gaussian pitch
angle distribution that was used in this event causes two subregions of higher
dose rates in the polar region of each hemisphere. The location of the peak
dose rate in each of these subregions corresponds to the location with the
strongest connection through the asymptotic direction to one of the peaks. The
subregions corresponding to the 0.0 radian pitch angle peak contain the most
intense dose rates, which is to be expected as the 0.0 radian pitch angle peak
represents the strongest pitch angle weighting factors. While the subregions
corresponding to the second pitch angle peak at 2.38 radian contain lower dose
rates, they encompass a wider area, which is likely because the peak is
broader than the first.

(a)

(b)

Figure 3: Dose rates at 39 kft in altitude due to galactic cosmic rays are
shown in fig. 3(a). The dose rates follow the structure of Earth’s magnetic
vertical cut-off rigidities, peaking in the polar regions. Dose rates due to
the initial anisotropic GLE spectrum are shown in fig. 3(b). The crosses on
fig. 3(b) are the locations where the connection of 1.5 GV protons to the 0.0
radian pitch angle peak (coloured in white) and 2.38 radian pitch angle peak
(coloured in yellow) was found to be maximal for the north and south
hemisphere.

An effective dose rate map for the initial phase of GLE71 has been previously
been calculated using CRAC:DOMO and presented in Figure 5 of
mishev2015computation for 39 kft. This dose rate map shows some similarities
to fig. 3(b), but also some differences, likely caused by a combination of
factors. The southern hemisphere regions show the same 0.0 radian pitch angle
peak dose rate feature between 50° and 100°; however, the map by Mishev et al.
shows no obvious sharp dose rate peaks in the northern hemisphere. A key
mapping difference here is that Mishev et al. used a smoothed 5° by 15°
resolution, in contrast to the 5° by 5° resolution in fig. 3(b), so it is
likely that many of the granular features that can be seen in fig. 3(b) have
been missed in the map by Mishev et al. Also note that as is discussed later
in this paper, the AniMAIRE dose rate yield functions likely differ from those
used by CRAC:DOMO.

Fig. 4 displays three dose rate maps calculated at three selected time periods
during the event. The three maps calculated for each time period were: an
isotropic map, an anisotropic dose rate map, and a total dose rate map, which
combines the anisotropic dose rate map with a background, cosmic ray dose rate
map. An animated map showing dose rates throughout the whole event has been
attached in the supplemental material for this paper, as well as outputted
dose rate data files.

The anisotropic dose rate maps were calculated as described above. For the
isotropic dose rate maps, $G(\alpha(P))$ was set to 1 for all pitch angles,
$\alpha$. This demonstrates the type of result which can be achieved with
isotropic-only models (which is a worst-case scenario), using the method of
Mishev et al. for isotropic dose calculation. This method is also similar to
that of copeland2008cosmic. For total dose rate maps, the anisotropic dose
rate map was added to the cosmic ray dose map, which in this case was
calculated using the previously mentioned DLR-modified ISO model [Matthiä .
(2013)] as implemented in the CosRayModifiedISO Python module [Davis (20242)],
supplied with a pre-event OULU neutron monitor count rate of 106.54 cts/s.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

Figure 4: Effective dose maps for an altitude of 39 kft for an isotropic GLE71
spectrum, the anisotropic GLE71 spectrum and pitch angle distribution, and the
total (anisotropic GLE + cosmic ray) dose for the anisotropic GLE71. As the
isotropic spectrum is a worst-case assumption, global dose rates are overall
significantly larger across Earth in the isotropic case than in the
anisotropic case but are missing any of the complex anisotropic structures
that are present in the anisotropic case.

Fig. 4 shows large differences between the isotropic and anisotropic
situations at the beginning of the event. The isotropic situation shows
broadly the same dose rates across the polar regions beyond a transition
region. In contrast, for the anisotropic event, there are two ”peaks” in the
polar regions, corresponding to the dual peaks of the double Gaussian pitch
angle distribution, but with a complex geographic curved shape formed due to
the influence of Earth’s magnetosphere and asymptotic directions. As the event
proceeds over time, the anisotropic maps become steadily more similar to the
isotropic maps due to the spreading of the Gaussian widths for the pitch angle
distribution. This convergence over time towards a more isotropic distribution
is common in GLEs [Shea Smart (2012)]. However, noticeable differences in dose
rate maps can be observed for this event even 65 minutes after it began.

In addition to investigating dose rates as a function of location on Earth,
the dependence of dose rates on altitude in the isotropic case vs the
anisotropic case has been examined. The altitude dependence of dose rates is
shown in fig. 5 for the isotropic case in the polar region and at four polar
coordinates surrounding Earth in the anisotropic case.

Figure 5: Effective dose rates as a function of altitude at different polar
latitude-longitude locations. The isotropic distributions are larger in all
cases, which is expected, and the curves are similarly shaped but do vary in
curvature, particularly towards higher altitudes.

While it is theoretically possible that the anisotropy could cause the
effective dose-altitude curve to vary in shape, the curves in fig. 5 all have
similar shapes. There is nevertheless some divergence in shape, particularly
at higher altitudes. For instance, the longitude=0.0 ° and longitude=180.0 °
cases begin to diverge at about 13 km in altitude.

One question of interest that arises is whether or not anisotropy in the event
is important for event-integrated doses, as anisotropy is typically only
significant in the beginning stages of an event. We find that for this event,
the anisotropy is indeed important, and this can be seen in the event-
integrated dose maps shown in fig. 6.

(a)

(b)

Figure 6: Integrated dose maps across the full event as defined by the
parameters given by Mishev et al., which were fed into AniMAIRE for both the
isotropic case and anisotropic case at 39 kft (excluding the final 5 minute
timestamp due to a minor numerical issue with applying parameters to the
supplied spectrum equation). Although the event is mostly anisotropic during
only the first half an hour, the anisotropy can clearly be observed in the
anisotropic map. This illustrates how the anisotropy must be considered to
accurately assess the total dose a flight will receive, even if the anisotropy
is primarily dominant specifically at the beginning of an event.

Wide hotspots of high dose rates in the polar regions can be observed in fig.
6(b), in contrast to the relatively flat polar dose rates in the isotropic
case. This indicates that the total dose experienced by an aircraft cannot be
calculated accurately using only isotropic assumptions - any aircraft
positioned in one of the polar anisotropic hot regions at the beginning of the
event will accumulate a much larger dose than an aircraft positioned at a
different polar region.

##  4 Comparisons to other models

In addition to analysing the event, it was also possible to compare isotropic
AniMAIRE GLE71 effective dose rates to dose rates produced by CRAC:DOMO and
WASAVIES. Fig. 7(a) shows isotropic GCR dose rates as a function of altitude
compared with dose rates given by CRAC:DOMO [Mishev . (2021)] and WASAVIES
[Sato . (2018)]. Comparisons with the GLE effective dose rates from CRAC:DOMO
are shown in fig. 7(b), and comparisons with WASAVIES data are shown in figs.
7(c) and 7(d).

(a)

(b)

(c)

(d)

Figure 7: Comparisons with maximum dose results from CRAC:DOMO and WASAVIES.
Comparisons between cosmic ray dose rates during GLE71 can be seen in fig.
7(a). Comparisons with GLE-only effective dose rates from CRAC:DOMO are shown
in fig. 7(b), and comparisons with WASAVIES data are shown in figs. 7(c) and
7(d), where AniMAIRE dose rates in fig. 7(d) have been arbitrarily divided by
a factor of 3 to show how the relationship between effective dose rates
calculated by AniMAIRE and WASAVIES match extremely well, except for differing
by a flat factor, indicating that this difference likely arises from a
difference in input spectra normalisation.

The GCR dose rates given by CRAC:DOMO are consistently larger than those
produced by AniMAIRE, particularly in the 20 kft to 40 kft range. The GLE
effective dose rates of CRAC:DOMO over time also follow a different profile to
those produced by AniMAIRE, as the dose rates of CRAC:DOMO decrease at a much
slower rate over time. Interestingly, the opposite conclusions are true for
the comparisons with WASAVIES dose rates. WASAVIES dose rates are
significantly lower than those produced by AniMAIRE; however, the variation in
dose rates over time matches extremely well, as does the ratio between dose
rates at different altitudes. This is illustrated in figs. 7(c) and 7(d) where
the AniMAIRE dose vs time curves match the WASAVIES curves remarkably well
when divided flatly by an arbitrarily chosen value of three.

A plausible explanation for the differences in shape and structure between
AniMAIRE and CRAC:DOMO is provided by the plot of cosmic ray comparisons shown
in fig. 7(a). AniMAIRE’s cosmic ray-induced dose rates agree well with the
WASAVIES cosmic ray-induced dose rates above 15 kft. In contrast, the
CRAC:DOMO cosmic ray-induced dose rates are much larger than the AniMAIRE dose
rates between about 20 kft and 35 kft. This implies that AniMAIRE and WASAVIES
use similar yield functions to calculate dose rates from the incoming spectra
but that the yield functions used by CRAC:DOMO are much larger in the 20 kft
to 35 kft region. This explains why the differences between CRAC:DOMO and
AniMAIRE’s dose rates in fig. 7(b) are largest in the 20 kft to 35 kft region
and it also explains why the AniMAIRE dose rate curves match the shape and
structure of WASAVIES but not CRAC:DOMO’s.

The fact that dose rates produced by AniMAIRE in this instance are
significantly larger than those predicted by WASAVIES is unsurprising as a
previous study by sato2018real using WASAVIES on GLE69 found that the
effective dose rates produced were a factor of between 2 to 3 lower than prior
published dose rates, including a factor of 2 lower than Mishev et al.’s model
(now named CRAC:DOMO) for that GLE. Sato et al. explain this as WASAVIES
predicting lower spectral indices for some GLEs than many other models, which
also causes other parameters in WASAVIES to vary to compensate, leading to
lower predicted dose rates. Note that WASAVIES here calculated spectral
indices to a precision of 0.5, and that a spectral index of exactly 4.5 was
used for almost the entirety of the first hour of GLE71.

Overall, these comparisons show that the AniMAIRE results for GLE71 fall
between those of CRAC:DOMO and WASAVIES, but agree better with WASAVIES in
terms of altitude structure and time variations. In future work, it would be
interesting for AniMAIRE to be run with other spectra and pitch angle
distributions and compared to other anisotropic models to see if these
differences remain, and it would also be interesting for comparisons to be
performed using other spectrum calculation algorithms and magnetospheric
models, which could be added to future versions of AniMAIRE.

##  5 Conclusion, and Future Possible Extensions

AniMAIRE is a new open-source toolkit for calculating radiation dose rates due
to GLEs and cosmic rays, and the openly available nature of AniMAIRE should
provide the community with an easy way to investigate how dose rates on Earth
vary during anisotropic conditions.

The analysis of the GLE71 dose rates calculated by AniMAIRE described in this
paper found several interesting results and showcases the capabilities of the
model. Comparing the dose rates calculated with assumptions of isotropy and
anisotropy highlights the variation in predicted radiation exposure between
the two methods. While this difference is most notable at early times,
comparison of event-integrated dose rate maps shows that the proper treatment
of the anisotropy in SEP protons is key to accurately predicting dose rates
during GLEs.

The dose rates predicted by AniMAIRE for GLE71 have also been compared to
those predicted by WASAVIES and Mishev et al. (CRAC:DOMO). The time-profile of
dose rates calculated with AniMAIRE agrees well with WASAVIES, although the
dose rates predicted by AniMAIRE are 3-4 times higher than those predicted by
WASAVIES. This discrepancy is consistent with previous comparisons which found
WASAVIES predicted lower counting rates than other models. In contrast, the
dose rates predicted by CRAC:DOMO are higher than those predicted by AniMAIRE
and decay much more slowly in time. This may be explained in part by CRAC:DOMO
employing yield functions which result in a much higher dose at aviation
altitudes.

The Pythonic design of AniMAIRE should mean it is relatively easy to plug into
other projects, and other members of the space weather community could connect
AniMAIRE with their own tools. Some ideas that could be implemented in the
future include adding GLE spectral calculations to AniMAIRE, perhaps for
instance using the MAIRE+ GLE spectrum algorithm [Hands . (2022)]. Another
idea might be to add options for using magnetospheric models other than
Tsyganenko 1989 or even alternatives to MAGNETOCOSMICS such as OTSO [Larsen .
(2023)]. Finally, work is ongoing in collaboration with the Community
Coordinated Modeling Center (CCMC) [Community Coordinated Modelling Center ()]
at NASA and with ESA through their Network of Models (NoM) [European Space
Agency ()] to create web interfaces for AniMAIRE so non-Python users can
perform runs without prior programming experience.

##  6 Software and Data Availability

AniMAIRE, and all of the software that have been developed as part of this
project have been released openly online. The AniMAIRE code itself is
available at its homepage on Github at https://github.com/ssc-maire/AniMAIRE-
public . The README file in the Github repository contains information on both
how to install AniMAIRE, and how to use it, including example commands and
outputs. The Github repository also contains two Jupyter Notebooks called
AniMAIRE_examples.ipynb, which can be run by users to produce several example
dose rate maps, and
notebooks_and_data_and_figures_for_paper/GLE71_plots_for_paper.ipynb, which
was the Jupyter Notebook used to generate all of the plots and data used for
this paper. All of the data and figures used for this paper can also be found
in the notebooks_and_data_and_figures_for_paper/ directory and subdirectories
in the Github repository.

While AniMAIRE can be installed from Github directly, we would recommend users
actually install it using the command:

          pip install AniMAIRE

in a Linux terminal. This command installs AniMAIRE from its PyPi page at
https://pypi.org/project/AniMAIRE/ and almost every package it depends on
automatically, without the need for any further commands. Note that if you use
Windows rather than Linux, you can easily install Windows Subsystem for Linux
(WSL) from https://apps.microsoft.com/detail/9mttcl66cpxj, which is a Linux
terminal that runs on Windows systems (AniMAIRE was actually developed using
WSL and has been tested on the Ubuntu 20.04.6 LTS distribution of WSL).

The only additional package you need to additionally manually install at
present to run AniMAIRE is the Geant4-based MAGNETOCOSMICS package [Desorgher
(2006), Desorgher (2004)] (the AsympDirsCalculator package, which AniMAIRE
currently uses, requires MAGNETOCOSMICS to calculate asymptotic directions and
cut-off rigidities). MAGNETOCOSMICS can be found and installed from
https://cosray.unibe.ch/~laurent/magnetocosmics/, and also from several other
locations around the internet. If you are struggling to install MAGNETOCOSMICS
or have any difficulties or questions in general regarding AniMAIRE
installation and usage, please feel free to contact us through posting on the
discussion section at the AniMAIRE Github page, direct Github messaging, or by
emailing us.

Some users may also be interested in separately using some of the Python
packages AniMAIRE depends on (and which get automatically installed with
AniMAIRE) that we developed as part of this project, some of which we’ve found
very useful on their own for our research. These are:

  1. 1.

ParticleRigidityCalculationTools [Davis (20243)]: This is a Python package for
performing simple conversions between kinetic energy and rigidity, and kinetic
energy spectra and rigidity spectra. Its Github repository is at
https://github.com/ChrisSWDavis/ParticleRigidityCalculationTools and its
installation page is at
https://pypi.org/project/ParticleRigidityCalculationTools/ . Its Zenodo DOI is
https://zenodo.org/doi/10.5281/zenodo.10992232 .

  2. 2.

DoseAndFluxCalculator [Davis Lei (2024)]: This is a Python package for
calculating dose rates and single event effects due to a particular given
particle spectrum hitting Earth’s atmosphere, by applying the yield functions
from previously developed MAIRE models [Lei (), Hands . (2022)] to the
inputted spectra. Its Github repository is at https://github.com/ssc-
maire/DoseAndFluxCalculator and its installation page is at
https://pypi.org/project/atmosphericRadiationDoseAndFlux/ . Its Zenodo DOI is
https://zenodo.org/doi/10.5281/zenodo.10992476 .

  3. 3.

CosRayModifiedISO [Davis (20242)]: This is a Python package for calculating
spectra due to cosmic rays given a particular date and time, solar modulation
or OULU neutron monitor count rate. It is based on a model developed by
matthia2013ready. Its Github repository is at https://github.com/ssc-
maire/CosRayModifiedISO and its installation page is at
https://pypi.org/project/CosRayModifiedISO/ . Its Zenodo DOI is
https://zenodo.org/doi/10.5281/zenodo.10992395 .

  4. 4.

AsymptoticDirectionsCalculator [Davis (20241)]: This is a Python package for
calculating asymptotic directions and vertical cut-off rigidities. Currently
it effectively acts as a Python wrapper for running MAGNETOCOSMICS, although
there have been discussions to eventually change this to an optional
dependency, and perhaps add alternative models such as OTSO[Larsen . (2023)].
Its Github repository is at https://github.com/ssc-
maire/AsymptoticDirectionsCalculator-public and its installation page is at
https://pypi.org/project/AsympDirsCalculator/ . Its Zenodo DOI is
https://zenodo.org/doi/10.5281/zenodo.10992606 .

Several animations can also be found in the supplemental materials of this
paper and in the AniMAIRE Github repository, which show the full progression
of GLE71 maps throughout the whole event. All of AniMAIRE’s releases, its
data, and plots have been archived on Zenodo at the DOI
https://zenodo.org/doi/10.5281/zenodo.10992643, and the licences for all
bespoke software described in this paper can be found in the software
repositories.

###### Acknowledgements.

We would like to thank Alexander Mishev, Nicholas Larsen and Sergey Koldobskiy
for their help and advice while performing this research, and aiding us with
performing comparisons with their work. We also acknowledge the NMDB database
(www.nmdb.eu), founded under the European Union’s FP7 programme (contract no.
213007) for providing data. The SWARM project, that this software’s
development was part of, was funded by NERC (grant no. NE/V002899/1) under the
SWIMMR programme, under the Strategic Priorities Fund (SPF) delivered by UKRI.

## References

  * Bender () PySpaceWeatherBender, S.  .  PySpaceWeather PySpaceWeather [].  https://github.com/st-bender/pyspaceweather version 0.2.4
  * Bruno . (2018) bruno2018solarBruno, A., Bazilevskaya, G., Boezio, M., Christian, ER., De Nolfo, G., Martucci, M.others  2018\.  Solar energetic particle events observed by the PAMELA mission Solar energetic particle events observed by the pamela mission.  The Astrophysical Journal862297.
  * Cannon . (2013) cannon2013extremeCannon, P., Angling, M., Barclay, L., Curry, C., Dyer, C., Edwards, R.others  2013\.  Extreme space weather: impacts on engineered systems and infrastructure Extreme space weather: impacts on engineered systems and infrastructure.  Royal Academy of Engineering.
  * Community Coordinated Modelling Center () CCMCCommunity Coordinated Modelling Center.  .  CCMC. CCMC.  https://ccmc.gsfc.nasa.gov
  * Cooke . (1991) Cooke1991Cooke, DJ., Humble, JE., Shea, MA., Smart, DF., Lund, N., Rasmussen, IL.Petrou, N.  19915\.  On cosmic-ray cut-off terminology On cosmic-ray cut-off terminology.  Il Nuovo Cimento C14213-234.  https://link.springer.com/article/10.1007/BF02509357 10.1007/BF02509357/METRICS
  * Copeland (2017) copeland2017cariCopeland, K.  2017\.  CARI-7A: Development and validation Cari-7a: Development and validation.  Radiation Protection Dosimetry1754419–431.
  * Copeland . (2008) copeland2008cosmicCopeland, K., Sauer, HH., Duke, FE. Friedberg, W.  2008\.  Cosmic radiation exposure of aircraft occupants on simulated high-latitude flights during solar proton events from 1 January 1986 through 1 January 2008 Cosmic radiation exposure of aircraft occupants on simulated high-latitude flights during solar proton events from 1 january 1986 through 1 january 2008.  Advances in Space Research4261008–1029.
  * Cramp . (1997) cramp1997octoberCramp, J., Duldig, M., Flückiger, E., Humble, J., Shea, M. Smart, D.  1997\.  The October 22, 1989, solar cosmic ray enhancement: An analysis of the anisotropy and spectral characteristics The october 22, 1989, solar cosmic ray enhancement: An analysis of the anisotropy and spectral characteristics.  Journal of Geophysical Research: Space Physics102A1124237–24248.
  * Davis (20241) AsympDirsCalculatorDavis, CSW.  2024104\.  AsymptoticDirectionsCalculator-public AsymptoticDirectionsCalculator-public [].  https://github.com/ssc-maire/AsymptoticDirectionsCalculator-public version 1.0.9  10.5281/zenodo.10992606
  * Davis (20242) cosraymodifiedisoDavis, CSW.  2024204\.  CosRayModifiedISO CosRayModifiedISO [].  https://github.com/ssc-maire/CosRayModifiedISO version 1.2.4  10.5281/zenodo.10992395
  * Davis (20243) ParticleRigidityCalculationToolsDavis, CSW.  2024304\.  ParticleRigidityCalculationTools ParticleRigidityCalculationTools [].  https://github.com/ChrisSWDavis/ParticleRigidityCalculationTools version 1.5.7  10.5281/zenodo.10992232
  * Davis Lei (2024) atmosphericRadiationDoseAndFluxDavis, CSW. Lei, F.  202404\.  DoseAndFluxCalculator DoseAndFluxCalculator [].  https://github.com/ssc-maire/DoseAndFluxCalculator version 1.0.4  10.5281/zenodo.10992476
  * Davis . (2024) AniMAIREDavis, CSW., Lei, F., Baird, F., Ryden, K. Dyer, C.  202404\.  AniMAIRE AniMAIRE [].  https://github.com/ssc-maire/AniMAIRE-public version 1.1.5  10.5281/zenodo.10992643
  * Desorgher () PLANETOCOSMICSDesorgher, L.  .  PLANETOCOSMICS. Planetocosmics.  http://cosray.unibe.ch/∼laurent/planetocosmics/
  * Desorgher (2004) desorgher2004magnetocosmicsDesorgher, L.  2004\.  MAGNETOCOSMICS Software User Manual Magnetocosmics software user manual [].  Physikalisches Institut University of Bern, http://reat.space. qinetiq.com/septimess/magcos/magnetocosmics_sum.pdf, Bern
  * Desorgher (2006) MAGNETOCOSMICSDesorgher, L.  200609\.  MAGNETOCOSMICS MAGNETOCOSMICS [].  http://cosray.unibe.ch/~laurent/magnetocosmics/ version 2.0
  * Desorgher . (2006) desorgher2006planetocosmicsDesorgher, L., Flückiger, EO. Gurtner, M.  2006\.  The planetocosmics geant4 application The planetocosmics geant4 application.  36th COSPAR Scientific Assembly 36th cospar scientific assembly ( 36, 2361).
  * Dyer . (2017) dyer2017extremeDyer, C., Hands, A., Ryden, K. Lei, F.  2017\.  Extreme atmospheric radiation environments and single event effects Extreme atmospheric radiation environments and single event effects.  IEEE Transactions on Nuclear Science651432–438.
  * Dyer . (2007) dyer2007solarDyer, C., Lei, F., Hands, A. Truscott, P.  2007\.  Solar particle events in the QinetiQ atmospheric radiation model Solar particle events in the qinetiq atmospheric radiation model.  IEEE Transactions on Nuclear Science5441071–1075.
  * Dyer . (1989) dyer1989measurementsDyer, C., Sims, A., Farren, J. Stephen, J.  1989\.  Measurements of the SEU environment in the upper atmosphere Measurements of the seu environment in the upper atmosphere.  IEEE Transactions on Nuclear Science3662275–2280.
  * Dyer . (1990) dyer1990measurementsDyer, C., Sims, A., Farren, J. Stephen, J.  1990\.  Measurements of solar flare enhancements to the single event upset environment in the upper atmosphere (avionics) Measurements of solar flare enhancements to the single event upset environment in the upper atmosphere (avionics).  IEEE Transactions on Nuclear Science3761929–1937.
  * European Space Agency () NoMEuropean Space Agency.  .  ESA Network of Models. ESA Network of Models.  https://nom.esa.int/
  * Hands . (2022) hands2022newHands, A., Lei, F., Davis, C., Clewer, B., Dyer, C. Ryden, K.  2022\.  A New Model for Nowcasting the Aviation Radiation Environment With Comparisons to In Situ Measurements During GLEs A new model for nowcasting the aviation radiation environment with comparisons to in situ measurements during gles.  Space Weather208e2022SW003155.
  * Hands . (2016) hands2016newHands, A., Lei, F., Ryden, K., Dyer, C., Underwood, C. Mertens, C.  2016\.  New data and modelling for single event effects in the stratospheric radiation environment New data and modelling for single event effects in the stratospheric radiation environment.  IEEE Transactions on Nuclear Science641587–595.
  * Karapetyan (2008) karapetyan2008detectionKarapetyan, G.  2008\.  Detection of high energy solar protons during ground level enhancements Detection of high energy solar protons during ground level enhancements.  Astroparticle Physics305234–238.
  * Kataoka . (2014) kataoka2014radiationKataoka, R., Sato, T., Kubo, Y., Shiota, D., Kuwabara, T., Yashiro, S. Yasuda, H.  2014\.  Radiation dose forecast of WASAVIES during ground-level enhancement Radiation dose forecast of wasavies during ground-level enhancement.  Space Weather126380–386.
  * Lantos Fuller (2003) lantos2003historyLantos, P. Fuller, N.  2003\.  History of the solar particle event radiation doses on-board aeroplanes using a semi-empirical model and Concorde measurements History of the solar particle event radiation doses on-board aeroplanes using a semi-empirical model and concorde measurements.  Radiation Protection Dosimetry1043199–210.
  * Lantos Fuller (2004) lantos2004semiLantos, P. Fuller, N.  2004\.  Semi-empirical model to calculate potential radiation exposure on board airplane during solar particle events Semi-empirical model to calculate potential radiation exposure on board airplane during solar particle events.  IEEE transactions on plasma science3241468–1477.
  * Larsen . (2023) larsen2023newLarsen, N., Mishev, A. Usoskin, I.  2023\.  A New Open-Source Geomagnetosphere Propagation Tool (OTSO) and Its Applications A new open-source geomagnetosphere propagation tool (otso) and its applications.  Journal of Geophysical Research: Space Physics1283e2022JA031061.
  * Latocha . (2009) latocha2009avidosLatocha, M., Beck, P. Rollet, S.  2009\.  AVIDOS—A software package for European accredited aviation dosimetry Avidos—a software package for european accredited aviation dosimetry.  Radiation Protection Dosimetry1364286–290.
  * Lei () MAIREUKLei, F.  .  Models for Atmospheric Ionising Radiation Effects - MAIRE. Models for Atmospheric Ionising Radiation Effects - MAIRE.  http://maire.uk/maire/
  * Lei . (2004) lei2004atmosphericLei, F., Clucas, S., Dyer, C. Truscott, P.  2004\.  An atmospheric radiation model based on response matrices generated by detailed Monte Carlo simulations of cosmic ray interactions An atmospheric radiation model based on response matrices generated by detailed monte carlo simulations of cosmic ray interactions.  IEEE Transactions on Nuclear Science5163442–3451.
  * Lei . (2005) lei2005improvementsLei, F., Hands, A., Clucas, S., Dyer, C. Truscott, P.  2005\.  Improvements to and validations of the QinetiQ Atmospheric Radiation Model (QARM) Improvements to and validations of the qinetiq atmospheric radiation model (qarm).  2005 8th European Conference on Radiation and Its Effects on Components and Systems 2005 8th european conference on radiation and its effects on components and systems ( D3–1).
  * Matthiä . (2013) matthia2013readyMatthiä, D., Berger, T., Mrigakshi, AI. Reitz, G.  2013\.  A ready-to-use galactic cosmic ray model A ready-to-use galactic cosmic ray model.  Advances in Space Research513329–338.
  * Meier . (2020) meier2020radiationMeier, MM., Copeland, K., Klöble, KE., Matthiä, D., Plettenberg, MC., Schennetten, K.Hellweg, CE.  2020\.  Radiation in the atmosphere—A hazard to aviation safety? Radiation in the atmosphere—a hazard to aviation safety?  Atmosphere11121358.
  * C. Mertens . (2023) mertens2023nairasMertens, C., Gronoff, G., Zheng, Y., Petrenko, M., Buhler, J., Phoenix, D.Minow, J.  2023\.  NAIRAS Model Run-On-Request Service at CCMC Nairas model run-on-request service at ccmc.  Space Weather215e2023SW003473.
  * C. Mertens . (2009) mertens2009developmentMertens, C., Tobiska, WK., Bouwer, D., Kress, B., Wiltberger, M., Solomon, S. Murray, J.  2009\.  Development of nowcast of atmospheric ionizing radiation for aviation safety (nairas) model Development of nowcast of atmospheric ionizing radiation for aviation safety (nairas) model.  1st AIAA Atmospheric and Space Environments Conference 1st aiaa atmospheric and space environments conference ( 3633).
  * CJ. Mertens . (2013) mertens2013nairasMertens, CJ., Meier, MM., Brown, S., Norman, RB. Xu, X.  2013\.  NAIRAS aircraft radiation model development, dose climatology, and initial validation. Nairas aircraft radiation model development, dose climatology, and initial validation.  Wiley Online Library.
  * Mishev . (2015) mishev2015computationMishev, A., Adibpour, F., Usoskin, I. Felsberger, E.  2015\.  Computation of dose rate at flight altitudes during ground level enhancements no. 69, 70 and 71 Computation of dose rate at flight altitudes during ground level enhancements no. 69, 70 and 71.  Advances in Space Research551354–362.
  * Mishev . (2021) mishev2021applicationMishev, A., Koldobskiy, S., Usoskin, I., Kocharov, L. Kovaltsov, G.  2021\.  Application of the verified neutron monitor yield function for an extended analysis of the GLE# 71 on 17 May 2012 Application of the verified neutron monitor yield function for an extended analysis of the gle# 71 on 17 may 2012.  Space Weather192e2020SW002626.
  * Mishev Usoskin (2016) mishev2016analysisMishev, A. Usoskin, I.  2016\.  Analysis of the ground-level enhancements on 14 July 2000 and 13 December 2006 using neutron monitor data Analysis of the ground-level enhancements on 14 july 2000 and 13 december 2006 using neutron monitor data.  Solar Physics2911225–1239.
  * Normand Baker (1993) normand1993altitudeNormand, E. Baker, T.  1993\.  Altitude and latitude variations in avionics SEU and atmospheric neutron flux Altitude and latitude variations in avionics seu and atmospheric neutron flux.  IEEE transactions on Nuclear Science4061484–1490.
  * Sato . (2018) sato2018realSato, T., Kataoka, R., Shiota, D., Kubo, Y., Ishii, M., Yasuda, H.Miyoshi, Y.  2018\.  Real time and automatic analysis program for WASAVIES: Warning system for aviation exposure to solar energetic particles Real time and automatic analysis program for wasavies: Warning system for aviation exposure to solar energetic particles.  Space Weather167924–936.
  * Shea Smart (2012) shea2012spaceShea, M. Smart, D.  2012\.  Space weather and the ground-level solar proton events of the 23rd solar cycle Space weather and the ground-level solar proton events of the 23rd solar cycle.  Space science reviews1711161–188.
  * Taber Normand (1993) taber1993singleTaber, A. Normand, E.  1993\.  Single event upset in avionics Single event upset in avionics.  IEEE Transactions on Nuclear Science402120–126.
  * Tobiska . (2015) tobiska2015advancesTobiska, WK., Atwell, W., Beck, P., Benton, E., Copeland, K., Dyer, C.others  2015\.  Advances in atmospheric radiation measurements and modeling needed to improve air safety Advances in atmospheric radiation measurements and modeling needed to improve air safety.  Space Weather134202–210.
  * Tsyganenko (1989) tsyganenko1989magnetosphericTsyganenko, NA.  1989\.  A magnetospheric magnetic field model with a warped tail current sheet A magnetospheric magnetic field model with a warped tail current sheet.  Planetary and Space Science3715–20.
  * Usoskin . (2015) usoskin2015databaseUsoskin, I., Ibragimov, A., Shea, M. Smart, D.  2015\.  Database of ground level enhancements (GLE) of high energy solar proton events Database of ground level enhancements (GLE) of high energy solar proton events.  Proceedings of Science, Proc. of 34th ICRC Hague, Netherlands30.


# An improved synthetic signal injection routine for HAYSTAC

Yuqi Zhu  Department of Physics, Yale University, New Haven, Connecticut
06520, USA  Wright Laboratory, Department of Physics, Yale University, New
Haven, Connecticut 06520, USA  M. J. Jewell  Department of Physics, Yale
University, New Haven, Connecticut 06520, USA  Wright Laboratory, Department
of Physics, Yale University, New Haven, Connecticut 06520, USA  Claire Laffan
Department of Physics, Yale University, New Haven, Connecticut 06520, USA
Wright Laboratory, Department of Physics, Yale University, New Haven,
Connecticut 06520, USA  Xiran Bai  Department of Physics, Yale University, New
Haven, Connecticut 06520, USA  Wright Laboratory, Department of Physics, Yale
University, New Haven, Connecticut 06520, USA  Sumita Ghosh  Wright
Laboratory, Department of Physics, Yale University, New Haven, Connecticut
06520, USA  Department of Applied Physics, Yale University, New Haven,
Connecticut 06520, USA  Eleanor Graham  Department of Physics, Yale
University, New Haven, Connecticut 06520, USA  Wright Laboratory, Department
of Physics, Yale University, New Haven, Connecticut 06520, USA  S. B. Cahn
Department of Physics, Yale University, New Haven, Connecticut 06520, USA
Wright Laboratory, Department of Physics, Yale University, New Haven,
Connecticut 06520, USA  Reina H. Maruyama<EMAIL_ADDRESS>Department
of Physics, Yale University, New Haven, Connecticut 06520, USA  Wright
Laboratory, Department of Physics, Yale University, New Haven, Connecticut
06520, USA  S. K. Lamoreaux  Department of Physics, Yale University, New
Haven, Connecticut 06520, USA  Wright Laboratory, Department of Physics, Yale
University, New Haven, Connecticut 06520, USA

###### Abstract

Microwave cavity haloscopes are among the most sensitive direct detection
experiments searching for dark matter axions via their coupling to photons.
The power of the expected microwave signal due to axion-photon conversion is
on the order of $10^{-24}$ W, so having the ability to inject synthetic axion
signals is helpful. Here we present a method based on frequency hopping spread
spectrum for synthesizing axion signals in a microwave cavity haloscope
experiment. It allows us to generate a narrow and asymmetric shape in
frequency space that mimics an axion’s spectral distribution, which is derived
from a Maxwell-Boltzmann distribution. In addition, we show that the synthetic
axion’s power can be calibrated with reference to the system noise. Compared
to the synthetic axion injection in HAYSTAC phase I, we demonstrated synthetic
signal injection with a more realistic lineshape and calibrated power.

##  I Introduction

The axion is a well-motivated solution to the strong charge–parity problem in
quantum chromodynamics and is also a dark matter candidate [1, 2, 3, 4]. The
most sensitive direct detection experiments by far are based on the microwave
cavity haloscope technique [5]. In a haloscope experiment, the axions would
convert into photons of equal energy inside a resonant microwave cavity
permeated by a strong magnetic field. The energy of each axion-converted
photon would be equal to the axion’s total energy—the sum of its rest-mass
energy $m_{a}$ and kinetic energy. Therefore, the lineshape of the resulting
microwave signal appearing in a haloscope experiment would inherit the axion’s
kinetic energy distribution. Suppose $\nu_{a}\\!=\\!m_{a}/(2\pi)$ is the axion
mass in SI-frequency units, we denote the spectral distribution of an axion
with mass $m_{a}$ as $f_{\nu_{a}}(\nu)$. The distribution $f_{\nu_{a}}(\nu)$
encodes axion’s properties derived from the pseudo-isothermal halo model [6],
and accounts for the modulation due to the Earth’s rotation around the center
of the galaxy [7, 8]. Here we describe a method to synthesize a microwave
signal whose lineshape resembles that from an expected axion signal. Injecting
such synthesized signals into the detector allows us to characterize the
detector’s response, validate the analysis procedure [9, 8, 10], and perform a
blind analysis [11].

Our method for synthesizing axion signals is inspired by a patent for radio
frequency (rf) hopping by Markey and Antheil [12]; it is related to the
frequency hopping spread spectrum (FHSS) technique, which has been applied in
military and wireless communication to prevent interception and reduce
interference. In our experiment, this technique allows us to generate an
axion’s spectrum by hopping between a large number of random rf tones sampled
from the axion’s spectral distribution. This way we can produce the spectral
spread associated with the axion’s kinetic energy distribution. As will become
more apparent later, the axion’s lineshape function is asymmetric and narrow.
The frequency hopping method allows us to overcome the technical difficulties
associated with signal shaping in the frequency domain. In principle, this
technique can be used to synthesize any generic spectral shapes. It is also a
realistic simulation: each axion’s energy is a random variable that follows
$f_{\nu_{a}}(\nu)$ and the observed lineshape is an ensemble-averaged result.

In this note, we describe the methods involved in signal injection in Sec. II
and provide a summary and outlook in Sec. III.

##  II Methods

PC2 PSG2 rf switch power splitter USB DAC Cavity AMPs VTS PSG1 $\bm{\otimes}$
LPF ADC SQZ circulatorVISATTLother tx in PC1 RFcalibrationLOweak
portreceiver’s port

Figure 1: Schematic. New components added for hardware injection are shown as
gray nodes; the rest already exists as part of the experiment. Blue arrows and
black dashed arrows indicate analog and digital signals respectively. PC2
controls PSG2 to inject a number of rf tones whose frequencies are sampled
from $f_{\nu_{a}}(\nu)$ to simulate axions with mass $\nu_{a}$. The synthetic
axion signals are injected through the weak port of the cavity and then
measured and amplified by the receiver chain consisting of two amplifiers
(AMPs): a JPA followed by a High-electron-mobility transistor. Subsequently,
the amplified signals are down-converted, mixed with the local oscillator (LO)
signal provided by PSG1, and then filtered by the low-pass filter (LPF). The
time-series signal collected by the Analog/Digital Converter (ADC) is then
Fourier transformed into an intermediate-frequency (IF) spectrum. The VTS
helps us characterize the system’s total noise by comparing the cavity’s
spectrum to blackbody spectra at controlled temperatures. The signal from the
VTS is only coupled in during a calibration measurement, indicated by the
dotted blue arrows.

The synthetic axion signal is injected as a transmission (tx) input through
the weak port of the cavity and subsequently coupled out from the receiver’s
port. Fig. 1 shows the schematic of synthetic signal injection. The main
device that synthesizes axion signals is an analog signal generator (model:
Agilent E8257D), labeled as PSG2. In between PSG2 and the cavity’s weak port,
an rf switch (model: ZFSWA2-63DR+) and a power splitter/combiner (model:
ZX10-2-71-S+) are added in this order. The rf switch allows us to interrupt
signal injection from PC1 via a USB DAC (model: USB-6009-NI). For instance,
the injection signals are blocked during calibration measurements between
axion datasets, as they would bias the noise measurement. The power splitter
allows us to simultaneously send in a weak tone near the JPA resonance as a
way to monitor and stabilize the JPA gain [9]. The variable temperature stage
(VTS) allows us to carry out in situ noise calibration [8] as part of the
power calibration for the axion signal injected. HAYSTAC has been using
squeezed-state receivers since 2021 [13]: the squeezing JPA (SQZ) is also
included in this diagram, as it affects the total noise, even though the
injected signal does not go through it; the full scheme with additional
details for the squeezed-state receiver chain can be found in Refs. [13, 14].

To synthesize the microwave signal for injection, we use the “Step (digital)
sweep” mode of PSG2 by providing it with a list of frequencies, amplitudes,
and dwell times. In aggregate, the list of rf tones would reproduce the
axion’s lineshape. In this experiment, by default, we fix the amplitude and
dwell time and only change the frequency of each point. The frequency
precision is set to 0.1 kHz or 8 significant digits for a GHz frequency,
chosen based on the precision of data acquisition. Due to the device’s output-
byte limit during transmission, the maximum length for such a list is about
$n_{8}\\!=\\!35$. For sub-kHz frequency precision, $N\\!\gg$2,000 samples are
needed to reveal the shape of the sampling distribution. Therefore, to iterate
through a list with $N\\!\gg\\!n_{8}$ rf tones, it is necessary to transmit
$\lceil N/n_{8}\rceil$ sublists. To begin, we establish the connection between
PSG2 and PC2 via virtual instrument software architecture (VISA) and enable
the rf output on PSG2. The steps for sending $N$ rf tones are as follows:

  1. 1.

Send $n_{8}$ frequencies, dwell time $\tau_{d}$, and amplitude
$P_{\text{syn}}$ to PSG2 from PC2 in one command.

  2. 2.

Wait for $\tau_{w}$ before sending another command to PSG2. Here
$\tau_{w}=n_{8}\tau_{d}+\tau_{r}$ includes the time required for the sweep
$n_{8}\tau_{d}$ and an additional uniformly random delay between 2 and 5
seconds $\tau_{r}\sim U(2,5)$ required for VISA communication.

  3. 3.

Check PSG2’s status to confirm that the sweep completed. If not, generate an
error message reporting PSG2’s status (settling or sweeping) and wait 2 more
seconds before checking the status again.

  4. 4.

Repeat steps 1-3 until all $N$ items are iterated over.

Based on the timing parameters, we can derive the average duty cycle of this
signal injection routine as $\bar{\eta}=n_{8}\tau_{d}/\bar{\tau}_{w}$, after
averaging over $\lceil N/n_{8}\rceil\gg$1 repetitions. By default, the dwell
time is $\tau_{d}=\\!5$ ms, resulting in $\eta=\\!4.8\%\pm\frac{1\%}{\lceil
N/n_{8}\rceil}$. In step 2, $\tau_{r}$ is randomized to avoid introducing any
additional frequency patterns into the signal.

###  II.1 Rejection sampling

To generate a list of random frequencies for signal injection, the probability
density function of the sampling distribution is the axion’s spectral function

$\displaystyle f_{\nu_{a}}(\nu)$
$\displaystyle=\frac{2}{\sqrt{\pi}}\left(\sqrt{\frac{3}{2}}\frac{1}{r}\frac{1}{\nu_{a}\braket{\beta^{2}}}\right)\sinh\left(3r\sqrt{\frac{2(\nu-\nu_{a})}{\nu_{a}\braket{\beta^{2}}}}\right)e^{-\frac{3(\nu-\nu_{a})}{\nu_{a}\braket{\beta^{2}}}-3r^{2}/2},$
(II.1)

where $\sqrt{\braket{v^{2}}}\approx$270 km/s is the virial velocity,
$\braket{\beta^{2}}=\braket{v^{2}}/c^{2}\approx\\!8\times 10^{-7}$ ($c$ is the
speed of light), $v_{s}\approx\\!220$ km/s the orbital velocity of the solar
system about the center of the galaxy, and
$r=v_{s}/\sqrt{\braket{v^{2}}}\approx\\!\sqrt{2/3}$ [8]. Given the function
form of $f_{\nu_{a}}$, it can be seen that applying inverse transform sampling
to sample from it is not obviously trivial. Here we use rejection sampling
[15] as an alternative. To apply it, we find another distribution whose
probability density function $y(\nu)$ satisfies $\xi\cdot
y(\nu)>f_{\nu_{a}}(\nu)\>\forall\nu$ for some $\xi>$0\. In our case, the
Cauchy distribution

$\displaystyle y(\nu)=\frac{1}{\sigma\pi}\frac{1}{1+((\nu-\mu)/\sigma)^{2}}$
(II.2)

with $\mu=\\!\nu_{a}+1$ kHz and $\sigma$=3 kHz, is sufficient because as
$\nu\rightarrow\infty$, $f_{\nu_{a}}/y\rightarrow 0$. Furthermore, the
cumulative function of a Cauchy distribution is invertible; therefore it is
straightforward to apply inverse transform sampling to sample from $y$. Let us
consider rejection sampling as an iterative process. In each iteration, we
sample a random number from $y$, noting it as $\nu^{\prime}$, and then decide
whether to reject or accept it based on another random number, $u$, chosen
from the uniform random distribution between 0 and 1. If $u\leq
f_{\nu_{a}}(\nu^{\prime})/(\xi\cdot y(\nu^{\prime}))$, we accept
$\nu^{\prime}$; otherwise reject. This iteration continues until we have $N$
samples. Fig. 2 shows the outcome of this method when the sampling
distribution is $f_{\nu_{a}}$ with $\nu_{a}$=4.7216300 GHz as an example. It
can be seen that the overall acceptance rate depends on the ratio of $\int
d\nu f_{\nu_{a}}$ to $\int d\nu\>\xi\cdot y$, i.e. the ratio of areas under
$f_{\nu_{a}}$ and $\xi\cdot y$; with $\xi=$3, $\int d\nu f_{\nu_{a}}/\int
d\nu\>\xi\cdot y\approx 0.3$.

Figure 2: Applying rejection sampling to create a list of $N\\!=\\!40,000$
random frequencies from an axion’s spectral distribution with
$\nu_{a}$=4.7216300 GHz. The envelope distribution (dash-dotted gray line) for
the sampling distribution $f_{\nu_{a}}$ (solid gray line) is obtained by
scaling the Cauchy distribution $y$ by $\xi=$3\. The histogram (blue) shows
the samples obtained using this method conform to the sampling distribution
$f_{\nu_{a}}$.

###  II.2 Power calibration

(a)

(b)

Figure 3: (a) Power excess in a rescaled spectrum overlaid with $f_{\nu_{a}}$
(dashed line). (b) Expected signal power $P_{\text{KSVZ}}$ (gray) and system
noise $h\nu_{c}N_{\text{sys}}$ integrated over an analysis binwidth
$\Delta\nu_{b}$ (orange) at fixed cavity frequency $\nu_{c}$ as functions of
detuning $|\nu-\nu_{c}|$. The green dashed line indicates the frequency of the
injected signal. The cavity presents a noise source to the receiver that is
frequency dependent and varies with $|\nu-\nu_{c}|$. The response to an axion
signal has a different frequency dependence that is more sharply peaked at the
cavity resonance. The plot shows the detuning-dependent variation of the axion
signal and system noise.

Though the power of each rf tone injected can be programmed on PSG2, neither
the weak port’s coupling efficiency nor the transmission line loss is known,
so it is necessary to carry out a power calibration measurement to determine
the time-averaged power of a synthetic axion signal. At fixed cavity
frequency, suppose we inject axions, each with power $P_{\text{syn}}$, near
the cavity resonance for duration $\tau_{\text{int}}$. The observed signal-to-
noise ratio (SNR) depends on the system noise power at the injected axion’s
detuning with respect to the cavity frequency $\nu_{c}$,
$h\nu_{c}N_{\text{sys}}(|\nu_{a}-\nu_{c}|)$ ($h$ is the Planck constant),
$P_{\text{syn}}$, $\tau_{\text{int}}$, $\eta$, and the measurement bandwidth
which is chosen to match the axion linewidth
$\Delta\nu_{a}\sim\nu_{a}/10^{6}$. This then allows us to estimate
$P_{\text{syn}}$ from a standard axion measurement

$P_{\text{syn}}=\text{SNR}\cdot(h\nu_{c}N_{\text{sys}}(|\nu_{a}-\nu_{c}|))\cdot\sqrt{\frac{\Delta\nu_{a}}{\eta\tau_{\text{int}}}}.$
(II.3)

As an example, Fig. 3 shows some results from a 12-hour measurement with the
cavity and JPA at fixed frequencies. From this dataset, SNR and
$N_{\text{sys}}$ are derived following the analysis routine outlined in Ref.
[8]. Each raw spectrum is first normalized by dividing out the average
baseline. This step allows us to identify bins contaminated by IF noise. After
removing contaminated bins from the first normalized spectrum, the spectrum is
normalized again by dividing out the Savitzky-Golay fit to itself and
subtracting 1. Then the resulting spectrum is rescaled to account for the
detuning-dependent sensitivity, producing the rescaled spectrum in Fig. 3(b).
To derive the SNR of the injected signal, we find the total power of the
injected signal by convolving the spectrum with the axion lineshape function
$f_{\nu_{a}}$. In this case, the resulting SNR is about 17.6. The scale factor
used in producing the rescaled spectrum is the ratio between the noise power
per analysis bin $\Delta\nu_{b}$=100 Hz and expected signal power, both of
which are plotted in Fig. 3(b). At detuning $|\nu_{a}-\nu_{c}|\sim$390 kHz,
the power of an axion-converted microwave signal is expected to be
$P_{\text{KSVZ}}=\\!5.6(1)\times 10^{-25}$ W, assuming the axion-photon
coupling strength $g_{a\gamma\gamma}$ (appearing in the axion-photon
Lagrangian) in the Kim-Shifman-Vainshtein-Zakharov (KSVZ) model [16, 17],
$g_{a\gamma\gamma}^{\text{KSVZ}}=\\!(-3.70\times 10^{-7}\text{/MeV}^{2})\cdot
m_{a}$ nominally [9]. This is calculated from various experimental parameters,
including the magnetic field strength and cavity volume [18, 13]. The system
noise in units of photon numbers $N_{\text{sys}}(|\nu_{a}-\nu_{c}|)$=0.26(2)
is derived from in situ calibration measurements [13, 14]. Evaluating Eq. II.3
with these values, we find $P_{\text{syn}}\\!=\\!2.2(2)\times 10^{-23}$ W, and
that the injected axions have a coupling strength at
$6.2(5)\>g_{a\gamma\gamma}^{\text{KSVZ}}$ level.

##  III Summary

We demonstrated synthetic axion injection using the rf hopping method. The
injected axion has the spectral distribution of a cold dark matter axion and
its power is calibrated using the total noise power as a reference. This is a
more realistic synthetic signal, in terms of spectral shape and power, as
compared to the prior demonstration in HAYSTAC Phase I (cf. Appendix F in Ref.
[8]). The code for implementing synthetic axion injection is made available as
an open-source project at
https://github.com/yuqizhuyqz/syntheticaxioninjection.git.

###### Acknowledgements.

We acknowledge the full HAYSTAC collaboration for the apparatus we used in
this experiment, including the microwave cavity and JPAs. HAYSTAC is supported
by the National Science 834 Foundation under grant numbers PHY-1701396,
PHY-835 1607223, PHY-1734006 and PHY-1914199 and the 836 Heising-Simons
Foundation under grants 2014-0904 and 837 2016-044. Y. Z., S. G., E. G., S.
C., M. J., and R. M. are supported in part by the Department of Energy under
Grant No. DE-AC02-07CH11359. S. G. is supported in part by the National
Science Foundation under Grant No. DMR-1747426.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

Yuqi Zhu: Investigation (equal); Data Curation (equal); Formal analysis
(equal); Writing–original draft (lead). M. J. Jewell: Investigation (equal);
Data Curation (equal); Formal analysis (equal); Writing—Review & Editing
(equal). Claire Laffan: Data Curation (equal); Formal analysis (equal). Xiran
Bai: Writing—Review & Editing (equal). Sumita Ghosh: Writing—Review & Editing
(equal). Eleanor Graham: Writing—Review & Editing (equal). S. B. Cahn:
Writing—Review & Editing (equal). Reina H. Maruyama: Project Administration
(equal); Funding acquisition (equal); Writing—Review & Editing (equal). S. K.
Lamoreaux: Project Administration (equal); Funding acquisition (equal);
Conceptualization(lead); Writing—Review & Editing (equal).

## Data Availability Statement

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## References

  * Peccei and Quinn [1977a] R. D. Peccei and H. R. Quinn, “$\mathrm{CP}$ conservation in the presence of pseudoparticles,” Phys. Rev. Lett. 38, 1440–1443 (1977a).
  * Peccei and Quinn [1977b] R. D. Peccei and H. R. Quinn, “Constraints imposed by $\mathrm{CP}$ conservation in the presence of pseudoparticles,” Phys. Rev. D 16, 1791–1797 (1977b).
  * Weinberg [1978] S. Weinberg, “A new light boson?” Phys. Rev. Lett. 40, 223–226 (1978).
  * Wilczek [1978] F. Wilczek, “Problem of strong $P$ and $T$ invariance in the presence of instantons,” Phys. Rev. Lett. 40, 279–282 (1978).
  * Sikivie [1983] P. Sikivie, “Experimental tests of the “invisible” axion,” Phys. Rev. Lett. 51, 1415–1417 (1983).
  * Jimenez, Verde, and Oh [2003] R. Jimenez, L. Verde, and S. P. Oh, “Dark halo properties from rotation curves,” Monthly Notices of the Royal Astronomical Society 339, 243–259 (2003).
  * Turner [1990] M. S. Turner, “Periodic signatures for the detection of cosmic axions,” Phys. Rev. D 42, 3572–3575 (1990).
  * Brubaker  _et al._ [2017a] B. M. Brubaker, L. Zhong, S. K. Lamoreaux, K. W. Lehnert, and K. A. van Bibber, “HAYSTAC axion search analysis procedure,” Phys. Rev. D 96, 123008 (2017a).
  * Brubaker  _et al._ [2017b] B. Brubaker, L. Zhong, Y. Gurevich, S. Cahn, S. Lamoreaux, M. Simanovskaia, J. Root, S. Lewis, S. Al Kenany, K. Backes, _et al._ , “First results from a microwave cavity axion search at 24 $\mu$eV,” Phys. Rev. Lett. 118, 061302 (2017b).
  * Al Kenany  _et al._ [2017] S. Al Kenany, M. Anil, K. Backes, B. Brubaker, S. Cahn, G. Carosi, Y. Gurevich, W. Kindel, S. Lamoreaux, K. Lehnert, _et al._ , “Design and operational experience of a microwave cavity axion detector for the 20–100 $\mathrm{\mu}$ev range,” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 854, 11–24 (2017).
  * Bartram  _et al._ [2021] C. Bartram, T. Braine, E. Burns, R. Cervantes, N. Crisosto, N. Du, H. Korandla, G. Leum, P. Mohapatra, T. Nitta, _et al._ (ADMX Collaboration), “Search for ”invisible” axion dark matter in the $3.3\text{-}4.2~{}{\mu}$ev mass range,” Phys. Rev. Lett. 127, 261803 (2021).
  * Markey and Antheil [1941] H. K. Markey and G. Antheil, “Secret communication system,”  (1941).
  * Backes  _et al._ [2021] K. M. Backes, D. A. Palken, S. A. Kenany, B. M. Brubaker, S. B. Cahn, A. Droster, G. C. Hilton, S. Ghosh, H. Jackson, S. K. Lamoreaux, _et al._ , “A quantum enhanced search for dark matter axions,” Nature 590, 238–242 (2021).
  * Backes [2021] K. M. Backes,  _A Quantum-Enhanced Search for Dark 1222 Matter Axions_ , Ph.D. thesis (2021).
  * Flury [1990] B. Flury, “Acceptance–rejection sampling made easy,” SIAM Review 32, 474–476 (1990).
  * Dine, Fischler, and Srednicki [1981] M. Dine, W. Fischler, and M. Srednicki, “A simple solution to the strong CP problem with a harmless axion,” Physics Letters B 104, 199–202 (1981).
  * Kim [1979] J. E. Kim, “Weak-interaction singlet and strong $\mathrm{CP}$ invariance,” Phys. Rev. Lett. 43, 103–107 (1979).
  * Sikivie [1985] P. Sikivie, “Detection rates for “invisible”-axion searches,” Phys. Rev. D 32, 2988–2991 (1985).