# Energetic Proton Propagation and Acceleration Simulated for the Bastille Day
Event of July 14, 2000

Matthew A. Young University of New Hampshire
Morse Hall
8 College Road
Durham, NH, 03824, USA  Nathan A. Schwadron University of New Hampshire
Morse Hall
8 College Road
Durham, NH, 03824, USA  Matthew Gorby  University of New Hampshire
Morse Hall
8 College Road
Durham, NH, 03824, USA  Jon Linker Predictive Science Inc.
9990 Mesa Rim Road
Suite 170
San Diego, CA, 92121, USA  Ronald M. Caplan Predictive Science Inc.
9990 Mesa Rim Road
Suite 170
San Diego, CA, 92121, USA  Cooper Downs Predictive Science Inc.
9990 Mesa Rim Road
Suite 170
San Diego, CA, 92121, USA  Tibor Török Predictive Science Inc.
9990 Mesa Rim Road
Suite 170
San Diego, CA, 92121, USA  Pete Riley Predictive Science Inc.
9990 Mesa Rim Road
Suite 170
San Diego, CA, 92121, USA  Roberto Lionello Predictive Science Inc.
9990 Mesa Rim Road
Suite 170
San Diego, CA, 92121, USA  Viacheslav Titov Predictive Science Inc.
9990 Mesa Rim Road
Suite 170
San Diego, CA, 92121, USA  Richard A. Mewaldt  California Institute of
Technology
Physics Department
Pasadena, CA, 91125, USA  Christina M. S. Cohen  California Institute of
Technology
Physics Department
Pasadena, CA, 91125, USA

(Received ???; Revised ???; Accepted ???)

###### Abstract

This work presents results from simulations of the 14 July 2000 (“Bastille
Day”) solar proton event. We used the Energetic Particle Radiation Environment
Model (EPREM) and the CORona-HELiosphere (CORHEL) software suite within the
SPE Threat Assessment Tool (STAT) framework to model proton acceleration to
GeV energies due to the passage of a CME through the low solar corona, and
compared the model results to GOES-08 observations. The coupled simulation
models particle acceleration from 1 to 20 $R_{\odot}$, after which it models
only particle transport. The simulation roughly reproduces the peak event
fluxes, and timing and spatial location of the energetic particle event. While
peak fluxes and overall variation within the first few hours of the simulation
agree well with observations, the modeled CME moves beyond the inner
simulation boundary after several hours. The model therefore accurately
describes the acceleration processes in the low corona and resolves the sites
of most rapid acceleration close to the Sun. Plots of integral flux envelopes
from multiple simulated observers near Earth further improve the comparison to
observations and increase potential for predicting solar particle events.
Broken-power-law fits to fluence spectra agree with diffusive acceleration
theory over the low energy range. Over the high energy range, they demonstrate
the variability in acceleration rate and mirror the inter-event variability
observed solar-cycle 23 GLEs. We discuss ways to improve STAT predictions,
including using corrected GOES energy bins and computing fits to the seed
spectrum. This paper demonstrates a predictive tool for simulating low-coronal
SEP acceleration.

††journal: ApJ

Draft

##  1 Introduction

The acceleration of solar energetic particles (SEPs), and their transport to
Earth and elsewhere in the heliosphere, has been recognized since the 1940s
(Forbush, 1946). Periods of enhanced energetic particle flux at 1 au present
opportunities to study these fundamental processes, but also pose an
increasing threat to geospace assets. Solar particle events (SPEs) can be
hazardous to space-based crews in Earth orbit (e.g., aboard the International
Space Station) and may imperil future crewed lunar or interplanetary missions
(Cucinotta et al., 2010, 2015, and references therein). They can even threaten
aircraft communication and navigation systems, and increase long-term health
risks for airline crews and passengers on polar flights. Therefore, providing
a meaningful estimate of energetic particle flux at a particular location in
the heliosphere is an important goal of current heliophysical research
efforts.

Large SPEs are typically associated with X-class solar flares and coronal mass
ejections (CMEs) with speeds exceeding 1000 km s-1, both of which typically
arise from complex sunspot groups (also known as active regions) with strong
magnetic fields. One mechanism for accelerating SEPs to the energies observed
at Earth is via shocks — or, more generally, compressions — that form low in
the corona during the passage of a CME. After a compression forms, it
propagates outward and accelerates particles over a finite space for a finite
time, at which point magnetic connectivity and cross-field diffusion determine
whether those particles arrive at an observer’s location.

An analysis of Type-II radio bursts by Gopalswamy et al. (2005) revealed that
CME shocks responsible for many SEP events formed at a few solar radii above
the solar surface, and that particle acceleration began at less than 10 solar
radii ($R_{\odot}$). Later, Reames (2009a, b) used velocity dispersion
analysis (VDA; see Kahler & Ragot (2006)) to derive solar particle release
(SPR) times for 30 of the 45 GLEs during solar cycles 20-23. The results
suggest particle release heights of approximately 2 $R_{\odot}$ for GLEs and
above 3 $R_{\odot}$ for non-GLEs. They also suggest that shock acceleration
occurs over a spatially broad region and occurs higher in the corona on the
shock flanks. The notion of connecting an SPR time to an altitude above the
solar surface extends back to Kahler (1994), who also noted that SEP injection
in three of the four analyzed GLEs resulted from a single CME-driven shock and
likely did not occur at the center of the CME.

Zhang et al. (2001) described a three-step process in which flare-associated
CMEs are initiated in the low corona (1.3-1.5 $R_{\odot}$) and accelerated out
to a few times $R_{\odot}$. Their work relied on four events with source
regions close to the limb observed by the Large Angle and Spectrometric
Coronagraph (LASCO) and EUV Imaging Telescope (EIT) aboard the Solar and
Heliosphysics Observatory (SOHO). Detailed observations of so-called EIT
waves, named after the instrument, prompted Liu et al. (2010) to propose that
the associated sharp fronts likely reflect CME acceleration in the low corona,
in support of the work by Zhang et al. (2001). Additional remote observations
have solidified the picture of low-coronal shock formation (Veronig et al.,
2010) and their association with significant increases in proton flux at 1 au
(Kozarev et al., 2011).

Mewaldt et al. (2012) showed that charge states observed for solar-cycle 23
GLE events require that the plasma density, $n$, and the particle acceleration
time, $\tau$, obey the relation $n\tau\geq 7.5\times 10^{9}$ s/cm3. That
implies acceleration heights of around 1.5 $R_{\odot}$, which they found to be
consistent with observations. They also note that acceleration along the CME
flank, where particles have more access to higher density, would raise
$n\tau$, making it easier to meet the minimum requirement.

Despite a wealth of energetic particle data at 1 au, our understanding of the
origin of energetic particle events is limited by the historical paucity of
observations closer to the Sun. While Parker Solar Probe and Solar Orbiter
promise to provided valuable insight into the near-Sun environment, analysis
of historical events must rely more heavily on remote observations, theory,
and simulations. For example: A combination of 3-D modeling, remote
observations, and 1-au data by Kouloumvakos et al. (2016) indicated that
previous SEP events might have been due to the first of two eruptive events in
which the arrival of the shock at the STEREO-B footprint was consistent with
the calculated SPR time, while the SPR time calculated for L1 was consistent
with acceleration along the CME’s western flank.

In this paper, we present results from simulations of the 14 July 2000
(“Bastille Day”) SPE using a software suite, called the SPE Threat Assessment
Tool (STAT) framework (Linker et al., 2019), that combines the Earth-Moon-Mars
Radiation Environment Model (EMMREM; Schwadron et al., 2010) and the CORona-
HELiosphere (CORHEL; Riley et al., 2012) software suite. The tool enables us
to model the proton acceleration process to GeV energies when a CME erupts,
expands, and propagates through the low solar corona. This helps us understand
how the acceleration process develops very close to the Sun — a region where
the interplay between CME initiation and expansion, the complex magnetic
structures of the corona, and the expansion of the solar wind make energetic
particle acceleration extremely difficult to accurately describe. Section 2
briefly describes the STAT framework and the components used in this
investigations. Section 3 describes the SPE event within the context of a
multi-day period of solar activity; section 4 presents the results of our
study; section 5 discusses the implications of our results and describes
future avenues of research; and section 6 concludes the paper.

##  2 STAT Framework

Modeling particle acceleration in the low coronal environment is especially
challenging. The plasma and magnetic properties can vary by orders of
magnitude, which leads to significant variation in the local Alfvén ($V_{A}$)
and sound speeds ($C_{S}$). To investigate this complex but crucial region,
STAT couples magnetohydrodynamic (MHD) simulations of CMEs in the low corona
with 3D solutions of the focused transport equation for SEPs. It allows users
to run the Energetic Particle Radiation Environment Module (EPREM, a component
of EMMREM) for precomputed Magnetohydrodynamic Algorithm outside a Sphere
(MAS, a component of CORHEL) simulations of real CME events to simulate SEP
events and provide diagnostics that can be compared with observations. STAT
and its diagnostics have been recently described by Linker et al. (2019); here
we briefly describe its primary components.

###  2.1 MAS Simulations

The properties of compressional regions such as shocks that drive SEP
acceleration depend critically on the properties of the local plasma
environment. To model a specific event, the simulation must realistically
capture these properties for the time period under study. The MAS model has a
long history of continued development and applications to this problem. While
models with a simple energy equation can qualitatively reproduce coronal
properties (Mikić & Linker, 1996; Linker et al., 1999; Mikić et al., 1999) and
are sufficient for exploring some dynamical aspects of boundary evolution
(Linker et al., 2011), so-called thermodynamic MHD models (Lionello et al.,
2009; Riley et al., 2011, 2012; Downs et al., 2013; Titov et al., 2017; Linker
et al., 2017; Mikić et al., 2018) are necessary to compute the plasma density
and temperature with sufficient accuracy to simulate EUV and X-ray emission
observed from space. In this approach, the energy equation accounts for
anisotropic thermal conduction, radiative losses, and coronal heating.
Inclusion of these extra physical terms is vital for obtaining realistic
$V_{A}$ and $C_{S}$.

To model a specific time period, a full-sun map of the photospheric magnetic
field is obtained from an observatory or flux transport model and processed to
create a boundary condition for the radial magnetic field (e.g., Linker et
al., 2017). For this event, we developed a thermodynamic MHD simulation of the
global corona using the procedure, equations, and coronal heating
specification described by Lionello et al. (2009); more recent MAS simulations
use a Wave-Turbulence-Driven description of coronal heating (e.g., Mikić et
al., 2018). In the thermodynamic model, the temperature at the lower boundary
is set to 20,000 K, similar to the upper chromosphere, and the upper boundary
is at 20 solar radii ($R_{\odot}$), beyond the sonic and Alfvén critical
points.

CORHEL computes MAS solutions in the coronal and heliospheric domains
separately. Coronal solutions are used to provide the inner boundary condition
for the heliospheric solutions for both steady-state background and dynamic
CME simulations (Lionello et al., 2013). In principle, MHD solutions for both
the coronal and heliospheric domains can be included in the EPREM calculation
(described in section 2.2). At the present time, STAT employs only the MHD
coronal domain within EPREM, and the remainder of the heliosphere is modeled
with a simple spiral magnetic field created with a radially constant solar
wind speed. This restricts us to modeling the first few hours of an SEP event;
once the CME leaves the coronal portion of the domain, possible SEP
acceleration from the CME propagation in the heliosphere is not modeled. STAT
is presently being modified to incorporate both the coronal and heliospheric
solutions in the EPREM simulations; these results will be the subject of
future papers.

###  2.2 EPREM Focused Transport Simulations

EPREM models energetic particle acceleration and transport using a Lagrangian
system, which co-moves with the plasma. To accomplish this, EPREM creates a
spherical shell of simulation nodes at each time step and advances each node
according to $\Delta\vec{r}=\vec{V}\Delta t$, where $\vec{r}$ is the 3-D node
displacement, $\vec{V}$ is the 3-D flow velocity, and $\Delta t$ is the time
step duration. It then calculates the distribution function for species $s$,
$f_{s}(t,r,p,\mu)$, according to the focused transport equation,

$\displaystyle\left[1-\frac{\left(\vec{V}\cdot\hat{b}\right)v\mu}{c^{2}}\right]\frac{df_{s}}{dt}$
(convection) $\displaystyle+v\mu\hat{b}\cdot\nabla f_{s}$ (streaming)
$\displaystyle+\frac{\left(1-\mu^{2}\right)}{2}\left[v\hat{b}\cdot\nabla\ln
B-\frac{2}{v}\hat{b}\cdot\frac{d\vec{V}}{dt}+\mu\frac{d\ln\left(n^{2}/B^{3}\right)}{dt}\right]\frac{\partial
f_{s}}{\partial\mu}$ (adiabatic focusing)
$\displaystyle+\left[-\frac{\mu}{v}\hat{b}\cdot\frac{d\vec{V}}{dt}+\mu^{2}\frac{d\ln\left(n/B\right)}{dt}+\frac{\left(1-\mu^{2}\right)}{2}\frac{d\ln
B}{dt}\right]\frac{\partial f_{s}}{\partial\ln p}$ (cooling)
$\displaystyle=\frac{\partial}{\partial\mu}\left(\frac{D_{\mu\mu}}{2}\frac{\partial
f_{s}}{\partial\mu}\right)+q\left(\vec{r},p,t\right)$ (pitch-angle scattering
and injection) (1)

where $t$ is time, $r$ is distance along the stream line, $p$ is momentum,
$\mu$ is the pitch-angle cosine, $\hat{b}$ is a unit vector parallel to the
magnetic field $\vec{B}$, $v$ is the particle velocity (distinct from the flow
velocity $\vec{V}$), $c$ is the speed of light, $q$ is the elementary charge,
and $D$ is the diffusion tensor.

Each node advances outward with the solar wind flow and is linked to nodes on
the neighboring shells. Each linked sequence of nodes defines a simulation
stream representing a velocity path line — the trajectory of fluid particles.
In steady-state (i.e., in the frame rotating with the Sun) these are also
streamlines; in places where the frozen-in assumption of ideal
magnetohydrodynamics (MHD) holds, these lines also represent magnetic field
lines.

The advantage of solving the transport problem in the co-moving frame is that
it precludes the necessity of computing spatial gradients in flow velocity,
which tend to introduce numerical errors that accumulate over many time steps.
Instead, it requires the relatively simple task of computing the rates of
change in plasma number density, $n$, and $\vec{B}$ at each stream node after
being moved by a timestep. This methodology is based on the approach described
in Kóta et al. (2005), which follows from the theory developed by Skilling
(1971) and Ruffolo (1995). It was used by Kozarev et al. (2013) to study time
dependent effects of solar energetic particle (SEP) acceleration in the low
corona during CME evolution and by Schwadron et al. (2014) to model radiation
doses at 1 au during a strong SPE event. In order to solve Equation 1, EPREM
needs some model of $n$, $\vec{B}$, and $\vec{V}$ at each node. Simplified
scenarios can use analytic forms of these plasma quantities but realistic
modeling requires the use of MHD data such as that provided by CORHEL.

##  3 Event Background

The Bastille Day SPE occurred on 14 July 2000, during a four-day period of
increased solar activity, including multiple CMEs, which produced high levels
of energetic electrons and ions at 1 au (Smith et al., 2001). This four-day
activity period began on 12 July 2000 with a chromospheric H$\alpha$ flare
observed by the Advanced Composition Explorer (ACE) spacecraft and ended on 16
July 2000 when energetic particle fluxes began to return to background levels.
The solar eruption on 14 July 2000 was one of the largest of solar cycle 23
and produced one of that cycle’s sixteen ground-level enhancement (GLE) events
(Bieber et al., 2002; Mewaldt et al., 2012).

The solar eruption that precipitated the Bastille Day SPE originated in active
region (AR) NOAA 9077. This active region produced an X-class flare starting
at 10:03 UTC and peaking at 10:24 UTC on 14 July 2000. A halo CME associated
with this flare subsequently became visible in the field-of-view of the Large
Angle and Spectrometric Coronograph (LASCO) on board the Solar and
Heliophysics Observatory (SOHO) at 10:54 UTC; initial speeds derived from
second-order (i.e., constant acceleration) fits to four points on the CME’s
leading edge ranged from 1300-1700 km s-1 (Andrews, 2001). Figure 1 shows
integral proton flux, in the top row, and X-ray flux, in the bottom row,
observed by GOES-08 before and during the SPE. The left column shows data from
13-16 July 2000 and the right column shows data from 08:00-13:00 UTC on 14
July 2000. Time in both columns has been shifted so that the start of the
event, at 10:03 on 14 July 2000, corresponds to 0.0 on all horizontal axes.
The spike in X-ray flux marks the event start at 10:03 UTC, shortly after
which proton integral flux rises sharply.

Figure 1: GOES-08 integral proton flux (top) and X-ray flux (bottom) during
13-16 July 2000 (left) and 08:00-13:00 UTC on 14 July 2000 (right). All times
have been shifted to the event start at 10:03 UTC on 14 July 2000.

##  4 Results

###  4.1 MHD Results

For our SPE modeling we utilize an MHD simulation of the 14 July 2000 Bastille
Day eruption that was performed previously using CORHEL/MAS. Detailed
descriptions of the simulation can be found in Linker et al. (2016) and Török
et al. (2018), here we restrict ourselves to the properties that are relevant
to our investigation.

The simulation is conducted in several steps. First, a surface magnetogram is
constructed and used to calculate a potential field that serves as the initial
condition for the global magnetic field. This field is then subjected to a
thermodynamic MHD relaxation (see Section 2.1) until a steady-state plasma and
magnetic field solution of the corona and solar wind is obtained (Figure
2(a)), which can then be validated with observations (Figure 2(b)). Next, in
order to construct a current-carrying pre-eruptive configuration that can
store the free magnetic energy required to power an eruption, a force-free
magnetic flux rope is inserted along the polarity inversion line (PIL) of NOAA
AR 9077. The highly elongated and curved rope is constructed using seven
instances of the analytical TDm model (Titov et al., 2014), and relaxed
towards a force-free state in a separate zero-$\beta$ MHD calculation (where
thermal pressure and gravity are neglected). After the rope’s insertion into
the coronal solution and a further short relaxation, localized boundary flows
converging towards the PIL are imposed, resulting in a slow rise of the flux
rope, followed by its rapid acceleration and eruption, leading to a fast CME.

Figure 2(c) shows the core of the flux rope (white field lines) about two
minutes after eruption onset. The shock front preceding the erupting rope
(shown here via electric currents; cf. Figure 4) forms early on in the
eruption, below $1.5\,R_{\odot}$. After its rapid acceleration low in the
corona, the flux rope (or CME) reaches an almost constant propagation speed of
$\approx 1500\,\mathrm{km\,s}^{-1}$, very similar to the observed speed (see
Section 3). Figure 2(d) shows a synthetic white-light image of the simulated
halo CME, in good agreement with a SOHO/LASCO 2 observation at a corresponding
time. We note that the coronal simulation was coupled to the heliospheric
version of MAS and continued until the CME passed the Earth (see Török et al.
2018 for details). However, coupling to the heliospheric domain in EPREM is
currently under development, so the analysis described in this paper used only
the coronal portion of the MAS simulation.

Figure 2:  Snapshots of the MHD simulation. (a) Coronal configuration seen
from Earth at $\approx$ 09:30 UT on 14 July 2000. NOAA AR 9077 is located
north of disk center. Streamers are visualized by electric currents; field
lines are colored by temperature. The magnetogram, $B_{r}(R_{\odot})$, is
saturated at 60 G. (b) Simulated and observed SOHO/EIT 195 Å emission of the
corona shortly before the eruption. Emission from under-resolved ARs is not
visible in the synthetic image. (c) Flux rope shortly after eruption onset.
Electric currents show the compression regions and the shock in front of the
rope. (d) Simulated (running-ratio) and observed (difference) SOHO/LASCO C2
white-light image of the halo CME, about 40 minutes after eruption onset. The
field-of-view is is 1.5-6 $R_{\odot}$ in the synthetic image; the green circle
marks the solar surface. Adapted from Török et al. (2018).

###  4.2 SEP Results

Figure 3: Integrated flux at 1 au from STAT output, with a call-out box
showing streams (labeled by stream number) near Earth.

Figure 3 shows the distribution of EPREM nodes at 1 au, at the end of the
simulation run. Color represents the amplitude of proton integral flux at
energies above 10 MeV. Images like the top panel are part of the standard
output of a STAT run. In both panels, a black diamond represents the
approximate position of Earth at the current time step. The bottom panel zooms
in on Earth’s position and labels eight EPREM streams of interest: the four
streams closest to Earth’s position at approximately the same heliographic
longitude (350-353), a stream near the center of the region of enhanced
integral flux (326), a stream with extremely high integral flux (361), and two
streams to the west of Earth (342 and 354) that are also near the edge of the
SEP event. Stream 351 sits just behind the diamond and is therefore closest to
Earth at the end of the simulation run, but stream 352 begins closest to Earth
and remains closest during the development of the energetic particle event.
From this point forward, figures will refer to this image and these stream
numbers.

Figure 4: EPREM integral flux at 1 au on the stream nearest Earth during the
early part of the simulation (352), compared to GOES-08. The sequences of
panels in the upper rows show CME eruption and integral flux evolution at $20\
R_{\odot}$ and $3\ R_{\odot}$. Each panel corresponds to a time step as
indicated by a vertical line in the plot in integral flux. .

STAT provides a global picture of integral proton flux in the standard
reference bins of the NOAA GOES satellites at 1 au. These are commonly used to
assess energetic particle threats to space-based assets; STAT diagnostics
allow inferences to be drawn about the physical drivers of these quantities.
Figure 4 shows integral proton flux in three energy channels — $E\geq$ 10, 50,
and 100 MeV — on 14 July 2000. Solid lines represent GOES-08 proton data from
the Energetic Particle Sensors (EPS) and dotted lines represent STAT output
calculated from EPREM stream 352. The horizontal axis of the bottom panel
represents time in hours since the start of the event at 10:03:00 UTC and the
vertical axis represents integral proton flux in counts cm-2 s-1 sr-1 . The
two rows of color images above the line plot show meridional cuts at 3 and 20
$R_{\odot}$ from the CORHEL run, at the times indicated by vertical lines on
the bottom plot.

In each image, the gray sphere represents the solar surface, thin lines
represent magnetic field lines, a thick white line indicates the peak in
velocity compression ($-\nabla\cdot\vec{V}$), and the color scale represents
the base-10 logarithm of integral flux at $E\geq 10$ MeV from 2 (blue) to 7
(red). The first three pairs of images capture the initial CME eruption out to
approximately $3\ R_{\odot}$, before flux levels have risen above background
at 1 au. Integrated flux begins to rise on this simulation stream around 20
minutes into the event, as the CME continues to expand and drive the
compression outward, but the rise in GOES-08 data lags the simulation. We will
address this discrepancy below. Simulated integral flux peaks around the level
of GOES-08 between 1.0 and 1.5 hours after the event start, before flattening
and beginning to fall off. This is due to the fact that the compression has
considerably weakened by approximately 11:25:00 UTC (fifth panel) and has
begun to leave the coronal domain. After three hours, the compression region
has passed $20\ R_{\odot}$ and the simulated integral flux has noticeably
dropped below GOES-08 levels. It is important to note here that, though EPREM
tracks protons out to 1 au, the MHD coupling in this simulation run extends
only out to $20\ R_{\odot}$ (hereafter, the “coronal domain”). Therefore,
particle acceleration only occurs within the coronal domain. This paper will
show that modeling acceleration within $20\ R_{\odot}$ captures many of the
observed characteristics of this event, but the reader must bare in mind that
results past roughly three hours after the start of the event do not account
for additional acceleration between $20\ R_{\odot}$ and 1 au (hereafter, the
“heliospheric domain”).

Figure 5: 3-D projection of simulation streams during the CME eruption, with
corresponding plots of integral flux at 1 au. Stream colors represent flux
from low (dark) to high (light) along each stream, relative to the other
streams pictured. The time shown is 7 minutes and 14 seconds after the start
of the event. Streams 350-353 are each close to Earth at some point throughout
the simulation. Stream 326 is near the center of enhanced integral flux.
Stream 361 has extremely high integral flux. Streams 354 and 342 represent
streams with relatively low integral flux. Nearby uncolored streams provide
additional context. Line plots show integral flux at $E\geq 10$, 50, and 100
MeV observed by GOES-08 (solid lines) and on each stream (dotted lines). The
horizontal and vertical axis labels shown for stream 326 apply to all streams.

Figure 5 shows a subset of the EPREM simulation streams (i.e., linked EPREM
nodes) at just over 7 minutes after the nominal event start. This UTC time
corresponds to time step 44 in the simulation, during the early stage of the
CME eruption. See the second pair of panels in the top rows of Figure 4 for
images of the CME eruption in the CORHEL simulation at the same time step. The
larger, central panel of Figure 5 displays eight streams in color: Streams
350, 351, 352, and 353, are each near the position of Earth (see Figure 3) at
some point throughout the simulation, stream 326 sits near the center of the
region of enhanced integral flux at 1 au, stream 361 has some of the highest
integral flux values, and streams 354 and 342 have comparatively low integral
flux values. The color of each stream represents the relative amount of flux
it sees during the simulation, with dark colors representing relatively low
flux and light colors representing relatively high flux. Since EPREM nodes
initially trace out a Parker spiral, unperturbed streams are essentially
radial at this distance. Meanwhile, streams perturbed by the CME drape around
its path, thereby outlining its 3-D structure. The bends and kinks in these
deformed stream lines turn out to be crucial in accelerating protons to SEP
energies.

The smaller panels in Figure 5 show simulated integral flux at 1 au on one of
the highlighted streams against GOES-08 integral flux during the 14 July 2000
SEP event. The panel layouts are similar to the bottom plot in Figure 4.
Simulated integral flux matches GOES-08 data to varying degrees at different
points throughout the simulation. The onset and rise of integral flux along
stream 352, which sits just to the right (heliographically west) of Earth at
the end of the simulation, matches GOES-08 extremely well for about the first
hour of the simulation run but turns over and falls below GOES-08 after the
initial acceleration period. Stream 351, which ends up closest to Earth, shows
integral flux increasing and rising before GOES-08 but decreasing to GOES-08
levels before falling below. Stream 350, which remains heliographically east
of Earth for the duration of the simulation run, shows an even sharper rise
and slightly higher peak amplitude compared to GOES-08. Near-Earth streams at
higher and lower heliographic latitudes exhibit a longitudinal trend in
integral flux similar to that of streams 350-353, as well as a slight north-
to-south trend toward lower amplitudes. The overall northeast-to-southwest
trend toward lower amplitudes mirrors the fact that Earth is in the lower
right portion of the elevated integral flux in Figure 3; the additional east-
to-west trend toward lower amplitudes and slightly slower rise times reflects
both Earth’s position relative to the integral flux enhancement as well as the
streams’ westward drift.

Figure 6: Envelope of EPREM integral flux at 1 au, on the three streams
closest to Earth throughout the simulation run, compared to GOES-08. The
integral flux envelope is defined as the area between the smallest and largest
values (inclusive) among all three streams. The large panel shows GOES-08 and
EPREM integral flux at $E\geq 10$, 50, and 100 MeV. Each panel in to the right
shows one of the individual integral energy bins that make up the larger
panel. As in Figure 5, the time axis spans the first three hours after 10:03
UTC.

The varying agreement with observations among streams is, in part, due to
their relative motion over Earth’s position. Figure 6 provides a rough
accounting of the fact that streams move over Earth’s position by showing the
integral-flux envelope on streams 351, 352, and 353. The envelope is defined
as the area between the minimum and maximum values (inclusive) among all three
streams, at each time step. The three integral flux bins are color coded as in
Figure 5, and the three-panel column on the right shows the individual
contribution of each bin. Despite the variation in lower and upper bounds, as
implied by Figure 5, the integral-flux envelope captures the rise and turn-
over reasonably well, especially in the $E>10$ MeV bin.

Figure 7: 3-D view of nodes showing acceleration rate, $\tau_{p}^{-1}$,
velocity divergence, $\nabla\cdot\vec{V}$, differential energy flux, $J(E)$,
and cosine of the angle between the magnetic field and plasma flow,
$\cos\theta_{BV}$, along the five streams indicated. The streams are shown at
10:10:14 UTC (7 minutes and 14 seconds after the start of the event). The
values of $\tau_{p}^{-1}$ and $J(E)$ apply to protons with energy 10.09 MeV.
An interactive version of each panel is available at
https://prediccs.sr.unh.edu/sim/eprem/bastille-day-energetic-protons/

Figure 7 presents four versions of the same low-coronal view during the CME
eruption as shown in Figure 5. Here, nodes along each stream are colored
according to one of four physical quantities: The top left panel shows the
local theoretical acceleration rate, $\tau_{p}^{-1}$ (Schwadron et al., 2015),
of protons with $E\approx 10$ MeV, derived from MHD quantities along each
stream. The top right panel shows velocity divergence $\nabla\cdot\vec{V}$ of
the MHD flow. The bottom left panel shows the pitch-angle-averaged flux of
protons with $E\approx 10$ MeV, computed from the simulated distribution
function. The bottom right panel shows the cosine of the angle between the MHD
magnetic and velocity fields.

Most streams with strong proton flux also have distinct bends where they drape
around the flank of the CME (e.g., stream 326). Conversely, streams with
little to no bend show weaker flux (e.g., stream 342). The regions of
relatively large $-\nabla\cdot\vec{V}$, corresponding to the CME-driven shock,
occur ahead of these stream-line bends whereas the peak in $\tau_{p}^{-1}$
occurs behind them. The shape of the CME suggested by the five stream lines
suggests that there is considerable acceleration along the CME flank and that
this acceleration is producing the flux observed at Earth. The change in
$\cos\theta_{BV}$ from 1 to near or slightly below 0 indicates a change from
field-aligned flow to cross-field flow and is approximately co-spatial with
the region of peak $\tau_{p}^{-1}$. One exception to the correlation between
stream-line bending and strong flux is stream 361, which traces more of a
helical structure at this time step. Section 5 will resolve this discrepancy
and further discuss the implications of these results.

Figure 8: Low-energy power-law index, $\gamma_{L}$, versus high-energy power-
law index, $\gamma_{H}$, for simulated and observed fluence spectra at 1 au.
The bottom panel expands the lower left corner of the top panel. In both
panels, filled circles represent spectral indices from all 16 GLEs observed by
GOES during solar-cycle 23. The gray circle highlights the Bastille Day event.
The gray square also corresponds to the Bastille Day GLE event, but was
derived using corrected data (see Section 5.3.1). The colored markers
represent Bastille Day EPREM streams, with the color scale indicating distance
from the center of elevated flux. The empty red circle corresponds to stream
351, which was near Earth. The top axis of the bottom panel lists the
equivalent compression ratio, $r_{c}$, from diffusive shock acceleration
theory for a given value of $\gamma_{L}$.

Broken power laws are characterized by a power law with spectral index
$\gamma_{L}$ over a low-energy domain, a power law with spectral index
$\gamma_{H}>\gamma_{L}$ over a high-energy domain, and a break energy,
$E_{0}$, separating the two domains. Proton fluence spectra observed at 1 au
tend to exhibit broken power-law shapes due to the interplay between particle
acceleration in the shock or compression ahead of a CME and the finite scales
of the accelerator region. The former process controls the lower spectral
index and the latter process controls the upper spectral index. Schwadron et
al. (2015) showed theoretical fluence spectra with and without loss of
particles due to a given flux bundle moving out of the acceleration region.
They demonstrated that the spectrum above the break energy steepens (i.e.,
softens) when incorporating particle loss from the accelerator region since
the effects of diffusion and convection can remove particles before they reach
high energies. In their theoretical framework, $\gamma_{L}$ follows the
prediction of diffusive shock acceleration (DSA), while $\gamma_{H}$ depends
on magnetic rigidity, $R_{g}$. The reason for the former is that DSA
efficiently accelerates particles to energies in the MeV range at quasi-
parallel shocks until the particles gain enough energy to escape along the
magnetic field. The reason for the latter is that parallel diffusion depends
on rigidity as $\kappa_{\parallel}=\kappa_{\parallel
0}\left(R_{g}/R_{g0}\right)^{\chi}$, where $R_{g0}$ is a reference rigidity.

Figure 8 shows spectral indices in the high-energy subdomain versus spectral
indices in the low-energy subdomain from power-law fits to simulated and
observed fluence at 1 au. The color scale represents distance from the center
of the region of enhanced integral flux, taken here to be stream 326, which is
at $\left(\theta,\phi\right)$ = $\left(70.4^{\circ},291.9^{\circ}\right)$ (see
Figure 3). Stream 351, which is approximately co-located with Earth at the
final time step, is $18^{\circ}$ from stream 326. The upper horizontal axis on
the bottom panel lists the compression ratio predicted by DSA for the
corresponding lower power-law index, $r_{c}=2(\gamma_{L}+1)/(2\gamma_{L}-1)$.
A dashed line indicates $\gamma_{L}=\gamma_{H}$.

The top panel of Figure 8 shows values up to 30 for each spectral index. The
dark-to-light trend in color as $\gamma_{H}$ increases implies that spectra
quickly steepen in the high-energy subdomain as distance increases away from
the region of enhanced integral flux. Note that this vertical trend continues
above $\gamma_{H}=30$. The spectral-index pairs shown in the top panel fall
into two groups, labeled “soft, gradual events” and “hard, prompt events”.
Section 5 will further discuss these groups. The bottom panel corresponds to
the bottom-left region of the left panel, outlined in orange, extending to 6
in both axes. A substantial fraction of these points fall within $45^{\circ}$
of the reference point at stream 326.

Both panels highlight additional relevant points in $\gamma_{H}-\gamma_{L}$
space (see legend on bottom panel): Filled circles mark the spectral indices
observed by GOES-08 during the 16 ground-level events (GLEs) of solar cycle 23
as reported by Mewaldt et al. (2012), with the Bastille Day (GLE) shown in
gray and all others shown in white. A gray filled square marks spectral
indices computed for this work using corrected GOES-08 fluence data – see
Section 5.3.1 for an explanation of this corrected data. Finally, an empty red
circle indicates the EPREM stream nearest Earth.

##  5 Discussion

Here we discuss the implications of the results presented in Section 4.
Specifically, the physical quantities presented in Figure 7 and what they
indicate about low-coronal acceleration, the distribution of spectral-index
pairs in Figure 8, and potential sources of uncertainty in predictions based
on our results.

###  5.1 Acceleration rates

Kóta et al. (2005) noted that the expanding CME causes magnetic field lines to
bend and compress, leading to two types of possible acceleration sites:
density jumps at the parallel shock and magnetic field amplification just
behind the shock. The images of stream 351 shown in Figure 7 suggest the
latter form of acceleration, in which magnetic-field amplification due to
localized deflections along the CME flank, rather than the strong plasma
compression in front of the CME, drives the high flux of energetic protons at
1 au. The highest values of $\tau_{p}^{-1}$ along this stream are $\sim 0.1$
s-1.

The analysis presented in Schwadron et al. (2015), which used the results of
previous MAS simulations of a CME, anticipated the results from the coupled
simulations shown here. They predicted that compressions, including (but not
limited to) shocks, can readily accelerate particles to tens of MeV or more in
a matter of minutes or less, especially as the geometry approaches
perpendicular. Figure 7 clearly shows large acceleration rates (small
acceleration times) near the westward flank of the simulated CME, where the
flow is far from parallel to the magnetic field. These extremely short
acceleration times – down to a few seconds on the stream with the highest
proton flux – rival those typically found in flare acceleration. The
relationship between velocity and magnetic field in the MHD simulation, as
captured by $\cos\theta_{BV}$, indicates that significant particle
acceleration occurs where the field lines drape across the CME flank, rather
than solely at the nose of the CME. The fact that $\tau_{p}^{-1}$ and $J(E)$
mirror the behavior of $\cos\theta_{BV}$ more closely than that of
$\nabla\cdot\vec{V}$ further suggest that the relatively weak compressions on
the CME flank are capable of accelerating protons to energies of tens to
hundreds of MeV.

The case of stream 361 in Figure 7 may appear, at first, to contradict the
argument that stream-line (or field-line) draping along the flank of a CME
produces strong particle acceleration in the low corona. In fact, it simply
serves as a reminder that there are multiple ways to accelerate solar-wind
particles to high energies. The lower apparent acceleration rate along stream
361 is partially due to the effect of projecting dynamic 3-D data onto a
static 2-D image but it is true that its peak values of $\tau_{p}^{-1}$ at the
time step shown are smaller than even those of stream 353. However, it
nevertheless produces significant low-coronal flux and extremely high 1-au
integral flux. The resolution of this discrepancy lies in the fact that the
region of large $-\nabla\cdot\vec{V}$ coincides with a sharp transition from
$\cos\theta_{BV}<0$ to $\cos\theta_{BV}>0$, indicating the location of a
quasi-parallel shock. Analysis of the density compression ratio along this
stream (not shown) suggests that a strong shock is directly responsible for
accelerating protons along this stream via classical DSA

###  5.2 Broken power laws

A second result from Schwadron et al. (2015) was that low-coronal particle
acceleration by compressions leads to broken power-law fluence spectra at 1
au. Our present work shows the consistency between simulated fluence at 1 au
and broken power laws formed by diffusive acceleration in shocks or
compressions, followed by trapping and further acceleration for a fraction of
the originally accelerated distribution.

A key feature of the bottom panel of Figure 8 is that values of $\gamma_{L}$
are constrained to within $1\leq\gamma_{L}\leq 2$, or $4\geq r_{c}\geq 2$. The
upper bound in $r_{c}$ represents the theoretical upper limit for strong
shocks; the lower bound represents a weaker yet still relatively strong shock
with respect to fast-forward shocks observed in the heliosphere (Kilpua et
al., 2015). Values of $\gamma_{H}$, on the other hand, exhibit a much wider
spread, implying that there is a high degree of variability in the spatial and
temporal scales of the coronal accelerator regions along the streams that
produce the SEPs observed at 1 au. Furthermore, the spread in fluence spectral
indices measured by GOES-08 during the 16 GLEs of solar cycle 23 is consistent
with the spread in EPREM fluences for this single simulated event. This latter
point suggests that inter-event variability among strong SEP events may be
due, at least in part, to intra-event variability as a function of
heliographic latitude and longitude.

Desai et al. (2016) found that low-energy spectral slopes in H-Fe spectra from
46 SEP events had values in the range 0.1-3 while high-energy spectral slopes
had values in the range 0.5-9, and that the high-energy slope was typically
larger than the low-energy slope. The values of $\gamma_{L}$ and $\gamma_{H}$
in Figure 8 exhibit a similar relationship to their findings. They also showed
that the spectral slopes were independent of species during a given SEP event.
This suggests that the results shown in Figure 8 should apply to other
species.

As noted in Section 4, the spectral-index pairs in the top panel of Figure 8
fall into two groups: soft, gradual spectra, with $\gamma_{H}\gtrsim 6.5$, and
hard, prompt spectra, with $\gamma_{H}\lesssim 6.5$. The spectral-index pairs
in the bottom panel correspond to hard, prompt spectra; they naturally include
all the GLEs since only sufficiently energetic particles produce the
atmospheric secondaries necessary for ground-level detection (Lopate, 2006).
This clustering into two groups is a result of the processes that create
different spectral forms in fluence at 1 au. Many of the soft, gradual spectra
occur far from the region of enhanced flux and essentially represent
background fluences. They develop extremely soft high-energy spectra simply
because there are few or no particles with energies above $E_{0}\approx 30$
MeV. The asymptotic trend in soft, gradual spectra toward $\gamma_{L}\approx
7$ suggests particle acceleration by compressions with $r_{c}$ as low as
$16/13\approx 1.23$. The hard, prompt spectra likewise cluster because of the
physical processes that produce measurable energetic particle events: The
field-line draping shown in Figure 7 accelerates a significant fraction of
available particles to high energies in general, but the specific variation in
$\tau_{p}$ among streams along the CME flank leads to the observed variability
in $\gamma_{H}$ at 1 au. At the same time, diffusive acceleration due to
shocks and compressions produces relatively consistent values of $\gamma_{L}$.

###  5.3 Uncertainty in predictions

####  5.3.1 GOES energy bins

The reliability of STAT predictions requires good agreement between simulation
results and data. In most cases, one assumes that observations represent true
values against which to judge simulation output. However, properly converting
raw data into higher-level products requires care and an appropriate
interpretation of the raw data. Sandberg et al. (2014) used a novel
intercalibration scheme to derive effective energy values for GOES/EPS bins,
using the IMP-8/Goddard Medium Energy Experiment (GME) data set as a
reference. See Cohen & Mewaldt (2018) for additional motivation behind
deriving effective energy values and ranges. So far, this work has compared
simulation results to the publicly available (uncorrected) GOES-08 data.
Figure 9 reproduces the integral flux envelope from Figure 6 (“Original”)
alongside a similar plot that uses re-calibrated GOES-08 energy bins
(“Corrected”). Due to the increased importance of corrections in higher energy
bins, the effect is more pronounced in the $\geq 50$ and $\geq 100$ MeV
channels. This effectively spreads out both the relative rise times and
amplitudes of the three traces. When compared to corrected integral flux, the
simulated integral flux envelope captures more of each GOES-08 trace,
especially in the highest energy bin.

Figure 9: GOES-08 integral proton flux using original (left) and corrected
(right) bins. The axis spans are identical to those in Figure 6.

####  5.3.2 Seed-particle spectrum

The observed fluence depends on the form of the seed population through which
the CME propagates. EPREM assumes a seed-particle spectrum that depends on
energy, $E$, and radial distance from the Sun, $r$, as

$f_{seed}(E,r)=\frac{m_{p}^{2}}{2E}\left(\frac{J_{0}}{\xi}\right)\left(\frac{r_{1}}{r}\right)^{\beta}\left(\frac{E}{E_{r}}\right)^{-\gamma}\exp\left(-\frac{E}{E_{0}}\right),$
(2)

where $m$ is the particle mass, $J_{0}$ is the amplitude of the 4He fluence
spectrum, $\xi$ is the ratio of 4He to H, $r_{1}$ is a fixed reference
distance, $\beta$ is the radial scaling dependence, $E_{r}$ is a reference
energy, $\gamma$ is the energy power-law dependence, and $E_{0}$ is the roll-
over energy. Table 1 lists the parameter values used for this simulation run.
The values of $\xi$, $\gamma$, and $E_{0}$ were chosen based on the work
described in Dayeh et al. (2009) (See also Kozarev et al., 2013). The value of
$J_{0}$ was a rough estimate based on averaged quiet-time GOES-08 proton flux
before the event; the value of $E_{r}$ simply sets the reference point for
this estimate. Future simulation runs will use physical estimates of seed-
population parameters derived from fits to 1-au satellite data (e.g., from
GOES-08, ACE, Wind, and STEREO), then compare the results to simulated fluence
and integral flux to improve predictions. Observations near the Sun by Parker
Solar Probe and Solar Orbiter will also provide direct estimates of the SEP
seed spectrum.

Table 1:  Seed spectrum parameters used in the simulation  Symbol | Value | Units | Description
---|---|---|---
$J_{0}$ | 20.0 | cm-2 s-1 sr-1 MeV-1 | Fluence amplitude
$\xi$ | 0.10 | – |  4He/H
$r_{1}$ | 1 | au | Reference distance from Sun
$\beta$ | 1.7 | – | Radial scaling dependence
$E_{r}$ | 1.0 | MeV | Power-law reference energy
$\gamma$ | 2.0 | – | Power-law index
$E_{0}$ | 1.0 | MeV | Roll-over energy

##  6 Conclusion

This paper has presented results from the Bastille Day energetic particle
event on 14 July 2000. We used two core elements of the STAT simulation
framework — EPREM and CORHEL — to model particle acceleration due to a CME-
driven shock, and subsequent transport to 1 au. The Bastille Day event was one
of the GLE events of solar cycle 23 and is an excellent example of low-coronal
proton acceleration. Comparison between simulated and observed integral flux
shows that, though acceleration in the coupled simulation occurs only within
20-30 $R_{\odot}$, simulated and observed integral flux at a single near-Earth
observer show good agreement early in the event. Including multiple simulated
near-Earth observers in a min-max envelope further improves the predictive
capabilities of this tool. Simulated particle arrival times match observations
very well despite the fact that the CME launch time was chosen based on EUV
observations, rather than a fit to particle arrival times.

Broken power-law fits to simulated fluence spectra at 1 au yield a trend in
low-energy spectral index versus high-energy spectral index that is consistent
with observations: The low-energy spectral index is confined to a narrow range
while the high-energy spectral index varies over a wide range. Values of the
low-energy spectral index match predictions of diffusive acceleration theory
for density compression ratios from 2-4, which are also observed in the low-
coronal portion of simulated field lines. Values of the high-energy spectral
index are determined by the spatial and temporal extent of low-coronal
particle-acceleration regions; they exhibit much greater variation due to the
low-coronal variation in density, velocity, and magnetic field during CME
passage. The variability in high-energy spectral index of simulated fluence is
consistent with variability in observed fluence during all 16 GLEs of solar
cycle 23, suggesting that at least some historical _inter_ -event variability
may be due to point observations of _intra_ -event variability.

The STAT framework is in active development, and we discussed two ways in
which we may improve future predictions: comparing simulated integral flux to
GOES data using corrected energy bins, and fitting the energetic-proton seed
spectrum to fluence data. Motivation for the first improvement comes from
published works that have demonstrated how novel intercalibration schemes for
GOES energy bins can increase agreement with the IMP-8 reference data set. We
found that 1-au integral flux from our simulation agreed better with the
corrected GOES-08 data during the early part of the event, and no worse after
the period of initial acceleration, than with the original GOES-08 data. The
second proposed improvement seeks to refine our method for estimating EPREM
input parameters related to the energetic-proton seed spectrum, which is
currently based on historical estimates of the 4He spectrum and does not vary
from event to event. We have begun by fitting the seed-spectrum functional
form to averaged GOES-08 proton fluence data from a few hours before each
event and will continue by incorporating additional data sets (e.g., ACE,
Wind, and STEREO). Results will be the subject of future publications. An
additional improvement to STAT’s predictive capabilities that is in progress
is a framework for coupling heliospheric MHD data from CORHEL into EPREM, to
allow acceleration beyond the current 20-30 $R_{\odot}$ domain. This is
critical for modeling a full event, which is a requirement in calculating
radiation dose levels.

We have thus provided results of the SPE Threat Assessment Tool (STAT)
framework for the onset of the Bastille Day Event of 2000. The simulation
itself roughly reproduces the timing of the event, the peak event fluxes, and
the spatial location of peak flux. The coupled EPREM-CORHEL model accurately
describes the acceleration processes active in the low corona and resolves the
sites of the most rapid acceleration close to the Sun. Disagreement between
simulated and observed particle fluxes after about three hours post-onset are
to be expected, since the coupling only extends to the outer edge of the
coronal domain. This tool nonetheless demonstrates predictive capabilities
that may be a key component for understanding the conditions associated with
space radiation events and the hazards they pose to future space-based
missions.

This work was supported by NASA under the LWS (grant # 80NSSC19K0067), O2R
(grant # 80NSSC20K0285), and SBIR/STTR programs, by the NSF PREEVENTS program
(ICER1854790), and AFOSR (contract # FA9550-15-C-0001). Computational
resources were provided by NASA’s NAS (Pleiades) and NSF’s XSEDE (TACC &
SDSC).

## References

  * Andrews (2001) Andrews, M. D. 2001, Solar Phys., 204, 179, doi: 10.1023/A:1014215923912
  * Bieber et al. (2002) Bieber, J. W., Dröge, W., Evenson, P. A., et al. 2002, Astrophys. J., 567, 622, doi: 10.1086/338246
  * Cohen & Mewaldt (2018) Cohen, C. M. S., & Mewaldt, R. A. 2018, Space Weather, 16, 1616, doi: 10.1029/2018SW002006
  * Cucinotta et al. (2015) Cucinotta, F. A., Alp, M., Rowedder, B., & Kim, M.-H. Y. 2015, Life Sciences and Space Research, 5, 31, doi: 10.1016/j.lssr.2015.04.002
  * Cucinotta et al. (2010) Cucinotta, F. A., Hu, S., Schwadron, N. A., et al. 2010, Space Weather, 8, S00E09, doi: 10.1029/2010SW000572
  * Dayeh et al. (2009) Dayeh, M. A., Desai, M. I., Dwyer, J. R., et al. 2009, Astrophys. J., 693, 1588, doi: 10.1088/0004-637X/693/2/1588
  * Desai et al. (2016) Desai, M. I., Mason, G. M., Dayeh, M. A., et al. 2016, in Journal of Physics Conference Series, Vol. 767, Journal of Physics Conference Series, 012004, doi: 10.1088/1742-6596/767/1/012004
  * Downs et al. (2013) Downs, C., Linker, J. A., Mikić, Z., et al. 2013, Science, 340, 1196, doi: 10.1126/science.1236550
  * Forbush (1946) Forbush, S. E. 1946, Physical Review, 70, 771, doi: 10.1103/PhysRev.70.771
  * Gopalswamy et al. (2005) Gopalswamy, N., Aguilar-Rodriguez, E., Yashiro, S., et al. 2005, J. Geophys. Res. Space Physics, 110, A12S07, doi: 10.1029/2005JA011158
  * Gopalswamy et al. (2014) Gopalswamy, N., Xie, H., Akiyama, S., Mäkelä, P. A., & Yashiro, S. 2014, Earth, Planets, and Space, 66, 104, doi: 10.1186/1880-5981-66-104
  * Kahler (1994) Kahler, S. 1994, Astrophys. J., 428, 837, doi: 10.1086/174292
  * Kahler & Ragot (2006) Kahler, S., & Ragot, B. R. 2006, Astrophys. J., 646, 634, doi: 10.1086/504674
  * Kilpua et al. (2015) Kilpua, E. K. J., Lumme, E., Andreeova, K., Isavnin, A., & Koskinen, H. E. J. 2015, J. Geophys. Res. Space Physics, 120, 4112, doi: 10.1002/2015JA021138
  * Kóta et al. (2005) Kóta, J., Manchester, W. B., Jokipii, J. R., de Zeeuw, D. L., & Gombosi, T. I. 2005, in American Institute of Physics Conference Series, Vol. 781, The Physics of Collisionless Shocks: 4th Annual IGPP International Astrophysics Conference, ed. G. Li, G. P. Zank, & C. T. Russell, 201–206, doi: 10.1063/1.2032697
  * Kouloumvakos et al. (2016) Kouloumvakos, A., Patsourakos, S., Nindos, A., et al. 2016, Astrophys. J., 821, 31, doi: 10.3847/0004-637X/821/1/31
  * Kozarev et al. (2013) Kozarev, K. A., Evans, R. M., Schwadron, N. A., et al. 2013, Astrophys. J., 778, 43, doi: 10.1088/0004-637x/778/1/43
  * Kozarev et al. (2011) Kozarev, K. A., Korreck, K. E., Lobzin, V. V., Weber, M. A., & Schwadron, N. A. 2011, Astrophys. J. Lett., 733, L25, doi: 10.1088/2041-8205/733/2/L25
  * Linker et al. (2016) Linker, J., Török, T., Downs, C., et al. 2016, in American Institute of Physics Conference Series, Vol. 1720, American Institute of Physics Conference Series, 020002, doi: 10.1063/1.4943803
  * Linker et al. (2019) Linker, J. A., Caplan, R. M., Schwadron, N., et al. 2019, in Journal of Physics Conference Series, Vol. 1225, Journal of Physics Conference Series, 012007, doi: 10.1088/1742-6596/1225/1/012007
  * Linker et al. (2011) Linker, J. A., Lionello, R., Mikić, Z., Titov, V. S., & Antiochos, S. K. 2011, Astrophys. J., 731, 110, doi: 10.1088/0004-637X/731/2/110
  * Linker et al. (1999) Linker, J. A., Mikić, Z., Biesecker, D. A., et al. 1999, J. Geophys. Res., 104, 9809, doi: 10.1029/1998JA900159
  * Linker et al. (2017) Linker, J. A., Caplan, R. M., Downs, C., et al. 2017, ApJ, 848, 70, doi: 10.3847/1538-4357/aa8a70
  * Lionello et al. (2013) Lionello, R., Downs, C., Linker, J. A., et al. 2013, Astrophys. J., 777, 76, doi: 10.1088/0004-637X/777/1/76
  * Lionello et al. (2009) Lionello, R., Linker, J. A., & Mikić, Z. 2009, Astrophys. J., 690, 902, doi: 10.1088/0004-637X/690/1/902
  * Liu et al. (2010) Liu, W., Nitta, N. V., Schrijver, C. J., Title, A. M., & Tarbell, T. D. 2010, Astrophys. J. Lett., 723, L53, doi: 10.1088/2041-8205/723/1/L53
  * Lopate (2006) Lopate, C. 2006, Washington DC American Geophysical Union Geophysical Monograph Series, 165, 283, doi: 10.1029/165GM27
  * Mewaldt et al. (2012) Mewaldt, R. A., Looper, M. D., Cohen, C. M. S., et al. 2012, Space Sci. Rev., 171, 97, doi: 10.1007/s11214-012-9884-2
  * Mikić & Linker (1996) Mikić, Z., & Linker, J. A. 1996, in American Institute of Physics Conference Series, Vol. 382, American Institute of Physics Conference Series, ed. D. Winterhalter, J. T. Gosling, S. R. Habbal, W. S. Kurth, & M. Neugebauer, 104–107, doi: 10.1063/1.51370
  * Mikić et al. (1999) Mikić, Z., Linker, J. A., Schnack, D. D., Lionello, R., & Tarditi, A. 1999, Physics of Plasmas, 6, 2217, doi: 10.1063/1.873474
  * Mikić et al. (2018) Mikić, Z., Downs, C., Linker, J. A., et al. 2018, Nature Astronomy, 2, 913, doi: 10.1038/s41550-018-0562-5
  * Reames (2009a) Reames, D. V. 2009a, Astrophys. J., 706, 844, doi: 10.1088/0004-637X/706/1/844
  * Reames (2009b) —. 2009b, Astrophys. J., 693, 812, doi: 10.1088/0004-637X/693/1/812
  * Riley et al. (2012) Riley, P., Linker, J. A., Lionello, R., & Mikic, Z. 2012, Journal of Atmospheric and Solar-Terrestrial Physics, 83, 1, doi: 10.1016/j.jastp.2011.12.013
  * Riley et al. (2011) Riley, P., Lionello, R., Linker, J. A., et al. 2011, Solar Phys., 274, 361, doi: 10.1007/s11207-010-9698-x
  * Ruffolo (1995) Ruffolo, D. 1995, Astrophys. J., 442, 861, doi: 10.1086/175489
  * Sandberg et al. (2014) Sandberg, I., Jiggens, P., Heynderickx, D., & Daglis, I. A. 2014, Geophysical Research Letters, 41, 4435, doi: 10.1002/2014GL060469
  * Schwadron et al. (2010) Schwadron, N. A., Townsend, L., Kozarev, K., et al. 2010, Space Weather, 8, S00E02, doi: 10.1029/2009SW000523
  * Schwadron et al. (2014) Schwadron, N. A., Gorby, M., Török, T., et al. 2014, Space Weather, 12, 323, doi: 10.1002/2014SW001086
  * Schwadron et al. (2015) Schwadron, N. A., Lee, M. A., Gorby, M., et al. 2015, Astrophys. J., 810, 97, doi: 10.1088/0004-637X/810/2/97
  * Skilling (1971) Skilling, J. 1971, Astrophys. J., 170, 265, doi: 10.1086/151210
  * Smith et al. (2001) Smith, C. W., Ness, N. F., Burlaga, L. F., et al. 2001, Solar Phys., 204, 227, doi: 10.1023/A:1014265108171
  * Titov et al. (2017) Titov, V. S., Mikić, Z., Török, T., Linker, J. A., & Panasenco, O. 2017, ApJ, 845, 141, doi: 10.3847/1538-4357/aa81ce
  * Titov et al. (2014) Titov, V. S., Török, T., Mikic, Z., & Linker, J. A. 2014, ApJ, 790, 163, doi: 10.1088/0004-637X/790/2/163
  * Török et al. (2018) Török, T., Downs, C., Linker, J. A., et al. 2018, Astrophys. J., 856, 75, doi: 10.3847/1538-4357/aab36d
  * Veronig et al. (2010) Veronig, A. M., Muhr, N., Kienreich, I. W., Temmer, M., & Vršnak, B. 2010, Astrophys. J. Lett., 716, L57, doi: 10.1088/2041-8205/716/1/L57
  * Zhang et al. (2001) Zhang, J., Dere, K. P., Howard, R. A., Kundu, M. R., & White, S. M. 2001, Astrophys. J., 559, 452, doi: 10.1086/322405