\documentclass[draft]{agujournal2019}
\usepackage{url}
\usepackage{lineno}
\usepackage[inline]{trackchanges}
\usepackage{soul}
\linenumbers

\draftfalse

\journalname{JGR: Space Physics}

\begin{document}

\title{Current Closure and Joule Heating in Data-Driven 3-D Auroral Arc Simulations}

\authors{J. van Irsel\affil{1}, K. A. Lynch\affil{1}, A. Mule\affil{1}, M. D. Zettergren\affil{2}, J. K. Burchill\affil{3}, L. J. Lamarche\affil{4}, D. L. Hampton\affil{5}}

\affiliation{1}{Department of Physics and Astronomy, Dartmouth College, Hanover, NH, 03755 USA}
\affiliation{2}{Physical Sciences Department, Embry‑Riddle Aeronautical University, Daytona Beach, FL, 32114 USA}
\affiliation{3}{Department of Physics and Astronomy, University of Calgary, Calgary, AB, T2N1N4 Canada}
\affiliation{4}{SRI International, Menlo Park, CA, 94025 USA}
\affiliation{5}{Geophysical Institute, University of Alaska Fairbanks, Fairbanks, AK, 99775 USA}

\correspondingauthor{Jules van Irsel}{jules.van.irsel.gr@dartmouth.edu}

\begin{keypoints}
\item Understanding current closure in discrete auroral arc systems requires data-driven, three-dimensional ionospheric simulations
\item Large-scale convection fields play a significant role in determining auroral arc current closure morphology and associated Joule heating
\item Details of precipitating electron energy distributions can significantly affect current closure and Joule heating in auroral arc systems
\end{keypoints}

\begin{abstract}
Discrete auroral arc systems, despite many symmetries, are three-dimensional in nature, encapsulating latitude and longitude variations in precipitation and field-aligned currents combined with important altitude variations in conductivities, hence closure currents. This study presents data-driven, 3-D numerical simulations of these processes based on a coordinated campaign of heterogeneous measurements collected from the Poker Flat Research Range during a sequence of Swarm spacecraft overpasses. These measurements include field-aligned current, global-scale convection flow, and auroral emissions, which are used to create top-boundary drivers for auroral arc simulations. Six conjunctions between the spacecraft, all-sky imagers, and radars are investigated and their measurements are used to simulate auroral arcs through multiple iterations per conjunction event. We look at different estimates of the background convection flow, assumptions about the energy distributions of electron precipitation, and along-arc structures in field-aligned current, and see what effect they have on current closure and Joule heating in auroral arc systems. Across the six conjunction events, 11 comparisons of auroral arc systems are presented, covering a catalog of 17 simulations in total. These comparisons allow us to look at the sensitivity of auroral arc systems to input parameters and envelop the simulations in a qualitative confidence interval. Our results suggest that discrete aurorae should be studied in three dimensions to fully understand field-aligned current closure and, by extension, Magnetosphere-Ionosphere-Thermosphere coupling. Additionally, our results demonstrate that both large-scale convection flows and specifics about the energy distributions of auroral precipitation can significantly affect current closure and Joule heating in auroral arc systems.
\end{abstract}

\section*{Plain Language Summary}
The aurora, or northern and southern lights, are embedded within a system of interacting electric and magnetic fields, and charged particles, the more energetic of which produce the lights themselves by exciting the neutral atmosphere. This brings about a three-dimensional current system and resistive heating, known as Joule heating. These currents enter and exit the atmosphere along the Earth's magnetic field, and can only close their circuit between altitudes of 80 -- 150 km, where the current carriers collide with the atmosphere. This paper outlines the importance of simulating aurorae in three-dimensions, and looks at how sensitive these simulations are to various input choices by observing the resulting differences in current connectivity and Joule heating. We look at collections of measurements from six different events and simulate them multiple times with different combinations covering 17 simulations in total. This allows us to gain insight into how much confidence can be had in our auroral arc simulations, and, by extension, what aspects are important to get right when studying auroral arcs. We conclude that large-scale plasma motion and the distribution of energies of the light-producing electrons both significantly affect the auroral system, and that current connectivity should be studied in three dimensions.

\section{Introduction}\label{s:intro}
Laws governing the physics of auroral arc systems are intrinsically three-dimensional---the conservation of mass, momentum, and energy density, in conjunction with Maxwell's equations, outline a system whose across-arc, along-arc, and field aligned directions are coupled. In the last decade or two, interest in three-dimensional (3-D) studies of the auroral ionosphere has slowly picked up \cite{amm2008, fujii2011, fujii2012, marghitu2012, zettergren2019, clayton2019, clayton2021, lynch2022, yano2021, vanirsel2024}, and we continue this trend by investigating quiet, discrete auroral forms in 3-D. Specifically, this paper looks at how electric current closure and Joule heating are affected by global-scale electric fields, the energy distributions of precipitating electrons, and along-arc structure in field-aligned currents (FAC), to provide insight into the geophysical domain of auroral arc systems.

The conductivity of the ionospheric volume surrounding auroral arcs is highly sensitive to impact ionization from electron precipitation \cite{fang2008,fang2010}. This ionization increases with increased energy flux, varies horizontally depending on arc structure, and varies in altitude depending on the energy distribution of the precipitation. Furthermore, the overarching, large-scale convection electric field guides the current continuity solution and directly affects the Joule heating of the system. For these reasons, to better understand auroral arc system currents, it is crucial that such systems are studied in 3-D.

Auroral-arc-scale science plays an important role in interpreting magnetosphere-ionosphere-thermosphere (MIT) coupling. The ionospheric end plays a non-passive role in this coupling \cite[\& references therein]{marghitu2012} and is involved in an ongoing sequence of system science studies \cite{wolf1975,seyler1990,cowley2000,lotko2004,fujii2011,fujii2012,marghitu2012,khazanov2018,clayton2019,clayton2021,yano2021,lynch2022,enengl2023,wang2024,vanirsel2024}. Such MIT studies require \textit{F}-region ionospheric maps of FAC and electric potential to be consistent with a 3-D ionospheric conductivity volume created by sunlight and charged-particle, auroral precipitation. However, what is often looked at is the two-dimensional (2-D) perspective of auroral arc systems, whether that is north-up or east-north. In the case of the horizontal ($\perp\mathbf{B}$) perspective, high-latitude electrostatic coupling assumes the height-integrated relation between quasi-static electric field, FAC, and conductances given by \citeA[Equation 8.15]{kelley2009}:
\begin{equation}
    j_\parallel(x,y)=\Sigma_P\nabla_\perp\cdot\mathbf{E}+\mathbf{E}\cdot\nabla_\perp\Sigma_P+\left(\mathbf{b}\times\mathbf{E}\right)\cdot\nabla_\perp\Sigma_H,
    \label{e:fac}
\end{equation}
where $(x,y)$ is the plane orthogonal to the local magnetic field, $j_\parallel$ is the ionospheric topside FAC, $\Sigma_{P,H}$ are the height-integrated Pedersen and Hall conductivities, i.e. conductances, $\mathbf{E}$ is the ionospheric electric field, and $\mathbf{b}=\mathbf{B}/B$ is the magnetic field direction. \citeA{yano2021} \cite<among others,>{marghitu2012, fujii2012} however, have pointed out that integrating out altitudinal effects can hide significant information regarding polar ionospheric systems, especially in terms of current closure. They use simplified 3-D Hall-magnetohydrodynamic simulations, taking into account ion-neutral collisions, to show that 2-D FAC closure assumed by the thin-layer approximation of the ionosphere is fundamentally different from the 3-D description, if alone for the fact that current streamlines can pass underneath one another.

The electric field solution from Equation \ref{e:fac} can be separated it into a constant, large-scale electric field, $\bar{\mathbf{E}}$, and a perturbation field, $\delta\mathbf{E}$, which gives two FAC contributions: $j_\parallel=\bar{j}_{\parallel}+\delta j_\parallel$ where
\begin{equation}
    \bar{j}_\parallel(x,y)=\bar{\mathbf{E}}\cdot\nabla_\perp\Sigma_P+\left(\mathbf{b}\times\bar{\mathbf{E}}\right)\cdot\nabla_\perp\Sigma_H,
    \label{e:bar_fac}
\end{equation}
and
\begin{equation}
    \delta j_\parallel(x,y)=\Sigma_P\nabla_\perp\cdot\delta\mathbf{E}+\delta\mathbf{E}\cdot\nabla_\perp\Sigma_P+\left(\mathbf{b}\times\delta\mathbf{E}\right)\cdot\nabla_\perp\Sigma_H.
    \label{e:delta_fac}
\end{equation}
After calculating and height-integrating the conductivities at a particular point in time, one can subtract $\bar{j}_\parallel$ from a specified \textit{F}-region map of FAC, $j_\parallel$, with which $\delta\mathbf{E}$ can be determined, i.e. solving current continuity and ionospheric Ohm's law with source term $\delta j_\parallel=j_\parallel-\bar{j}_\parallel$. In this sense, the electrostatic drivers are $j_\parallel$ and $\bar{\mathbf{E}}$, and the ionosphere responds by introducing polarization fields to provide the remaining FAC. In other words, $\delta\mathbf{E}$ is a result from local polarization charge densities within the ionospheric volume, while $\bar{\mathbf{E}}$ is an electric field external to our auroral-arc-scale system. With this perspective, a constant global estimate of the background flow, $\bar{\mathbf{v}}=\bar{\mathbf{E}}\times\mathbf{b}/B$, from either SuperDARN or PFISR, is an additional current driver and thus should be accounted for when interpreting FAC observations. Both Equation \ref{e:fac} and topics discussed in this paper deal with self-consistency, not causal relationships, when finding solutions to auroral current continuity.

\citeA{marghitu2012} reviews sequentially more complex descriptions of auroral arcs, the first of which takes on a band of enhanced uniform conductance with negligible altitudinal thickness and polarization electric fields that are fully in the across-arc direction. Having no along-arc gradients whatsoever results in FAC closure which relies only on Pedersen currents (see Equation \ref{e:fac}), while the electrojet current flows underneath, but plays no part in FAC closure. The second description introduces an along-arc component in the electric field which can greatly enhance the auroral electrojet current by means of the Cowling effect \cite{cowling1932}. With a partial Cowling channel (one with some FAC blockage), \citeA{amm2008} point out that this requires taking into account the ionospheric thickness when looking at current continuity. This is because, as \citeA{yano2021} have also pointed out, divergence-free currents cannot flow through one another. \citeA{amm2011,fujii2011,fujii2012} therefore take on a finite length Cowling channel model, which includes a thin Pedersen layer on top of a thin Hall layer, allowing for primary and secondary Pedersen and Hall currents to connect. The third description by \citeA{marghitu2012} only ignores the along-arc variation in the electric field, but does take on gradients of conductance along the arc. To understand FAC closure with this description, \citeA{marghitu2012} uses 2-D (east-north) modeling given the non-linear nature of this problem. \citeA{marghitu2012} concludes, however, that, even though various one- or two-dimensional descriptions of auroral arcs capture a substantial interpretation, a complete 3-D description is necessary to fully understand, even sheet-like, auroral arc systems.

This paper builds on work done by \citeA{clayton2021}, who study auroral arc systems and, to do so, developed new methods for driving simulations with 2-D maps of auroral data to study the surrounding ionosphere in 3-D. Similar to the work presented in this paper, they use multi-spectral auroral imagery from the Poker Flat Digital All-Sky Camera to both (a) infer the electron precipitation energetics and (b) replicate one-dimensional, in situ measurements of plasma flow, creating continuous 2-D driver maps. Their plasma flow measurements are provided by the Isinglass sounding rocket campaign and the replication methods are described by \citeA{clayton2019}. In this paper, we use replication methods by \citeA{vanirsel2024}, which expand upon these ideas, yet altered slightly in order to use in situ FAC data from orbital spacecraft (Swarm) instead of plasma flow data. With these tools, and given an abundance of observational datasets from the winter months of 2023, we explore the dependencies of current closure paths and Joule heating in auroral arc systems on different values of $\bar{\mathbf{E}}$, forms of electron precipitation spectra, and top-boundary FAC structures.

In this paper, we aim to determine geophysical, self-consistent solutions to ionospheric current continuity in non-ideal discrete auroral arcs that posses structure in across-arc, along-arc, and field aligned directions. In doing so, we explore how to properly drive 3-D simulations of auroral arc systems using 2-D electrostatic, continuous top-boundary conditions from distributed, multi-platform datasets: all-sky, multi-spectral imagery, in situ FAC data, and radar-based background convection flow data. Additionally, we study the sensitivity of current continuity solutions to various driver parameters, particularly background convection flow and precipitation parameters, in order to envelop auroral arc simulations in a form of qualitative confidence estimates. This provides a better understanding of the dominant physics behind auroral current closure and Joule heating for different situations. Ancillary to this, this study provides a catalog of auroral arc simulations covering six conjunction events with multiple modeling iterations per event, as well as driver and visualization tools to facilitate future studies of auroral arc systems.

In Section \ref{s:methods} we outline the instrumentation used in this work, a brief description of the ionospheric model used to produce our simulations, along with methods for imagery inversion, the replication technique, the implementation of precipitating electron impact ionization, and our use of flux tubes for 3-D visualization of current closure. Section \ref{s:events} summarizes the 6 conjunction events and Section \ref{s:results} covers the simulation results and comparisons thereof. We conclude our findings and discuss possible improvements and future uses of our work in Section \ref{s:conclusions}. \ref{a:diffnumflux} covers the derivation of the differential hemispherical number flux of accelerated Maxwellian precipitation, and figures of simulations not included in the main body of this paper are in the Supporting Information, along with other supporting figures and descriptions.

\section{Observational Data, Instrumentation, \& Methodologies}\label{s:methods}
The data products we use are of six conjunction events that are part of the \textit{Swarm-over-Poker-2023} campaign. This campaign facilitated simultaneous observations in February -- March, 2023, of a variety of auroral arcs during times when the European Space Agency's (ESA) Swarm spacecraft orbited overhead of the Poker Flat Research Range in Alaska. These observations are of key ionospheric electromagnetic parameters including, but not limited to, (1) the ESA Swarm mission's ion flow data from the Thermal Ion Imagers \cite<TII,>{knudsen2017} and (2) FAC data derived from its magnetometers \cite{swarmfacdata}, (3) convection flow data from AMISR's Poker Flat Incoherent Scatter Radar \cite<PFISR,>{pfisrdata1,pfisrdata2,pfisrdata3}, (4) global convection flow maps from the Super Dual Auroral Radar Network \cite<SuperDARN,>{superdarn1}, and (5) multi-spectral, all-sky imagery from the Poker Flat Digital All-Sky Camera \cite<DASC,>{dascdata1}. The six events are chosen based on them (a) occurring during February or March of 2023, (b) including at least one Swarm crossing with imagery, (c) having an unobstructed night sky, and (d) having one or more arcs that fit the quiet discrete arc description \cite{karlsson2020}. Figure \ref{f:context} shows the geographical context of the February 10, 2023 conjunction event. In this section we cover the details surrounding these data products and any methodologies applied to them, as well as the model used to create our auroral arc simulations.

\begin{figure}
    \centering
    \includegraphics[width=0.9\textwidth]{01_grid_context.pdf}
    \caption{Geographical context of our simulations, using the February 10, 2023 conjunction as an example, showing the model space (black), the Swarm A and C crossings (yellow), the PFISR track (green), the top-boundary for the driver maps (red), the approximate location of the imagery from below (blue), and a symbolic depiction of some flow vectors from the SuperDARN global data map (orange) on top of Alaska.}
    \label{f:context}
\end{figure}

\subsection{Poker Flat Digital All-Sky Cameras \& Imagery Inversion}\label{s:inversion}
The all-sky, multi-spectral auroral imagery we use comes from the University of Alaska Fairbanks Geophysical Institute's Poker Flat Digital All-Sky Cameras \cite<DASC,>{dascdata1} located at  212.57° east and 65.12° north (geographic). From this imagery we use a Python-based routine and the GLobal airglOW model \cite<GLOW,>{solomon2017} to produce estimated maps of both total precipitating energy flux, $Q_p$, and expected energy, $\langle E\rangle$. In this work, the expected energy is either the characteristic energy, $E_0$, or acceleration potential, $U_a$ (see Section \ref{s:precip}).

As shown by \citeA{rees1974}, and later expanded on by several others \cite{strickland1989,janhunen2001,hecht2006,grubbs2018a,grubbs2018b}, the ratio of green line (558 nm) to red line (630 nm) intensity for emissions driven by electron precipitation mostly depends on $\langle E\rangle$, while the blue line (428 nm) intensity mostly depends on $Q_p$. Roughly following \citeA{grubbs2018b}, we use GLOW, driven with ionospheric background conditions, to generate a lookup table of emission line intensities for a variety of driving precipitation energy spectra. Each energy spectrum in the table is parameterized by its values of $Q_p$ and $\langle E\rangle$, and GLOW simulates emission line intensities separately for each.

After denoising and calibrating the imagery, mapping each color to its rough emission altitude, and removing background brightness, we apply a simple Python routine \cite{asilink} that uses the lookup tables to invert each usable pixel of the image to a value of $Q_p$ and $\langle E\rangle$, along with rough error bars associated with the inversion. After inversion, all precipitation maps are Gaussian smoothed in the magnetic northward direction with a window size of 32 km ($\sigma\approx5.3$ km).

\subsection{Swarm Spacecraft}\label{s:swarm}
The European Space Agency's Swarm mission consists of three satellites which were launched into nearly polar, low Earth orbits on 22 November, 2013, with the goal of providing highly detailed measurements of variations in the Earth's magnetic field. We use their version 0401, level 2 FAC data derived down to 1 Hz from the Vector Field Magnetometer \cite<VFM,>{swarmfacdata} data, along with their version 0302, level 1B Electric Field Instruments data, specifically the 2 Hz TII ion drift measurements \cite{knudsen2017,burchill2022}. The TII data, like the precipitation maps, are Gaussian smoothed to 32 km, while the FAC data are smoothed to 16 km ($\sigma\approx2.7$ km) to account for their inherently finer structure. The ion drift measurements have a 100 -- 200 m/s one-sigma accuracy, and are used only in our discussions (Section \ref{s:discussions}) as a point of comparison with our simulation results.

We note that the choice of smoothing window, an important and carefully deliberated choice, strongly affects the science scales we can investigate. The specific smoothing window is chosen to match and align the available input data scales; we know that driving the model with inconsistent drivers (i.e., fine-scale fields data and large-scale imagery) leads to spurious signatures. For this study, therefore, we have not fully characterized the dependence on this scale choice. Instead we focus our studies on permutations of input parameters at these scales (i.e. on/off or from instrument A versus instrument B and so on).

\subsection{Poker Flat Incoherent Scatter Radar}
The Poker Flat Incoherent Scatter Radar \cite<PFISR,>{pfisrdata1,pfisrdata2,pfisrdata3} is an Advanced Modular Incoherent Scatter Radar facility and has been operational since 2007. PFISR is located at the Poker Flat Research Range (212.53° E, 65.13° N), which is owned by the University of Alaska Fairbanks Geophysical Institute, and the radar is maintained for the US National Science Foundation by SRI International. The antenna boresight points at an azimuth of 15° east-of-north and elevation of 74°. In this paper, we take single-value, uniform averages of plasma drift velocity within the latitude ranges of our simulation regions, and use these averages as large-scale background flow estimates. We use their resolved vector velocity (``vvels'') data based on long pulse experiments with a five minute integration time. These data products are produced by Python scripts found at \cite{vvelslink}. We use these data to provide one plasma drift velocity average per conjunction event.

\subsection{Super Dual Auroral Radar Network}
The Super Dual Auroral Radar Network (SuperDARN) is comprised of 35+ HF and VHF radars located across the northern and southern hemispheres and is operated by 20 institutions across 10 nations. This paper uses plasma convection flow estimates over Poker Flat, AK---one global estimate per conjunction event---that are interpolated by the pyDARN open-source python library. \citeA{superdarn1} describe the SuperDARN global-scale network and the pyDARN repository can be found at \cite{pydarnlink}. SuperDARN convection map data shown in this paper were processed using the FITACF3 algorithm with a spectral width-based Heppner-Maynard Boundary. Both the order and degree of the fit was 6.

\subsection{FAC Replication}\label{s:replication}
Our simulations require spatially continuous, topside ionospheric FAC maps. \citeA{vanirsel2024} outline how this can be done for electrostatic plasma convection maps. Here we have adjusted their methods for FAC maps instead. The replications can be done using distributed optical data, provided by all-sky, multi-spectral imagery, combined with FAC data tracks, provided by spacecraft or sounding rockets. We first invert the imagery using methods outlined in Section \ref{s:inversion}, from which preliminary estimates of the height-integrated conductivities (conductances) are gathered using Equations 3 and 4 in \citeA{robinson1987}. The conductance maps are then queried for two isocontours at user-defined conductance values which are the primary and secondary arc boundaries. With these boundaries, the replication process is as follows:
\begin{enumerate}
    \item The original FAC data track is translated in the east-north plane by some amount following the primary arc boundary such that the original and replicated data are equal at the primary boundary-track intersections.
    \item The replicated data track is scaled in the along-track direction such that the original and replicated data are equal at the secondary boundary-track intersections.
    \item This replication is repeated for multiple translations along the arc until the top-boundary is filled with FAC values at a sufficient replication density.
    \item The replicated FAC data map is then interpolated onto the simulation grid, providing the top-boundary simulation driver.
\end{enumerate}
For replications whose data lie just outside of the simulation region, the arc boundaries are extrapolated, ensuring sensible matching between FAC and precipitation. For example uses of the (weighted) replication method, see Sections 2.2 and 2.3 in \citeA{vanirsel2024}. For use of this method with FAC data, the electrostatic enforcement (Section 2.2.3) is omitted and the FAC data is treated as one of the components of the flow data.

\subsection{GEMINI Simulations}
Simulations for this study use the Geospace Environment Model of Ion-Neutral Interactions \cite<GEMINI,>{zettergren2012,zettergren2019}. GEMINI solves for 3-D electrostatic current continuity and ionospheric Ohm's law, accounting for changes in state parameters which affect conductivities as it steps forward in time; it calculates the electric field that is consistent with how the top-boundary FAC requirements connect through the ionospheric volume---one whose conductivity is highly sensitive to impact ionization from electron precipitation, which is implemented into GEMINI using methods by \citeA{fang2008,fang2010}.

GEMINI is a multi-fluid (electrons and six ion species), quasi-electrostatic model with its calculations of particle continuity consisting of chemical production/loss and photo/impact ionization. Calculations of local densities, plasma flows, and temperatures are treated self-consistently and the model includes thermal conduction heat flux, collisional heating, thermoelectric electron heat flux, and inelastic cooling/heating from photoelectrons. This is supplemented with Maxwell's equations and, at the time of writing, includes no displacement current or magnetic induction. With this, the system is solved through enforcing divergence-free currents, curl-free electric fields, and invoking Ohm's law. GEMINI can be driven with (aside from maps of precipitation energetics handling impact ionization) a map of FAC \textit{or} electric potential at the top-boundary. When driving GEMINI with a top-boundary map of FAC, a user-specified background electric field, $\bar{\mathbf{E}}$, is input separately. GEMINI assumes equipotential magnetic field lines, providing horizontal electric fields that are constant in altitude \cite{farley1959}. For a full description of the governing equations solved by GEMINI, see \citeA[Appendix A]{zettergren2015a}.

\subsection{Electron Precipitation Methods}\label{s:precip}
\subsubsection{Electron Precipitation Energy Spectra}
For auroral arc systems, electron precipitation energy spectra, $\phi(E)$, are often assumed to be of a standard unaccelerated Maxwellian form \cite{fang2008} whose differential number flux, $\phi_u(E)$, is
\begin{equation}
    \phi_u(E)=\frac{Q_p}{2E_0^2}\frac{E}{E_0}\exp\left(-\frac{E}{E_0}\right),
    \label{e:unacc-max}
\end{equation}
where $Q_p$ is the total precipitating energy flux, $E_0$ is the characteristic energy, and $E$ is the precipitation energy. This has its flux peak at an energy of $E_0$, representing the arc energy, however, it also incurs an energy spread of
\begin{equation}
    \sqrt{\langle(E-E_0)^2}\rangle=\sqrt{\frac{\int_0^\infty\left(E-E_0\right)^2\phi_u(E)\mathrm{d}E}{\int_0^\infty\phi_u(E)\mathrm{d}E}}=\sqrt{3}E_0.
\end{equation}
In contrast to this formulation, in auroral situations, there is often an accelerated signature \cite{evans1968,issi2003}, where the energy spread is related to the source region thermal motions, while the peak energy is related to the auroral acceleration region \cite{evans1974}. Therefore, we look at an alternative $\phi(E)$; that of an accelerated Maxwellian whose differential number flux, $\phi_a(E)$, is (see \ref{a:diffnumflux})
\begin{equation}
    \phi_a(E)=\frac{Q_p}{T_s^2+\left(T_s+U_a\right)^2}\frac{E}{T_s}\exp\left(-\frac{E-U_a}{T_s}\right),\ E\ge U_a,
    \label{e:acc_max}
\end{equation}
where $T_s$ is now the source region characteristic energy, and $U_a$ is the auroral acceleration region potential drop. The energy spread, with $U_a>T_s$, is $\sqrt{\langle(E-U_a)^2\rangle}=T_s\sqrt{(6+2U_a/T_s)/(1+U_a/T_s)}$. With $U_a/T_s\sim3$, which is not untypical, this has an energy spread of $\sqrt{3}T_s$. This choice for $\phi(E)$ has decoupled the energy spread and peak energy, which in this case is $U_a$ when $U_a>T_s$, which is the case for all our conjunction events.

Relationships between the acceleration potential and the source region/ionospheric characteristic energy exists via the FAC this system holds \cite{knight1973, ronnmark2002}, but these are not the focus of this paper. Equation \ref{e:acc_max} is implemented into GEMINI using methods described by \citeA{fang2010}. Both the GLOW model and the methods described by \citeA{fang2008,fang2010} take into account secondary and back-scattering electrons \cite{evans1974}.

Figure \ref{f:spectra} shows examples of $\phi_u(E)$ and $\phi_a(E)$ (Equation \ref{e:unacc-max} and \ref{e:acc_max}) with $U_a=E_0=3$ keV and $T_s=490$ eV. Both these spectra have the same integrated energy flux, $Q_p$, and both peak at 3 keV, yet the accelerated Maxwellian has a significantly lower energy spread: 0.8 keV compared to 5.2 keV in the unaccelerated case. Along with this, their respective electron density altitude profiles are shown, determined using the GLOW transport model \cite{solomon2017}. It is evident that the assumption of $\phi_u(E)$ can overestimate the electron density at lower altitudes given the high energy tail of these spectra. It is noted that with $U_a=0$, i.e. no auroral acceleration, we have $\phi_a/\phi_u=1$, which covers the relatively low energy background precipitation surrounding auroral arcs. This fact is what we use to determine $T_s$.

\begin{figure}
    \centering
    \includegraphics[width=\textwidth]{02_spectra_comparison.pdf}
    \caption{Comparison between unaccelerated and accelerated Maxwellian electron precipitation spectra. (a) Normalized energy spectra of $\phi_u(E)/Q_p$ (red) and $\phi_a(E)/Q_p$ (blue). Note that both spectra peak at 3 keV. (b) Electron density altitude profiles modeled by GLOW \cite{solomon2017} with the same color scheme.}
    \label{f:spectra}
\end{figure}

\subsubsection{Determining Source Region Characteristic Energies, $T_s$}\label{s:ts}
The differential number flux for an accelerated Maxwellian population approaches that of the unaccelerated population as $U_a$ approaches zero. In this limit $T_s$ becomes analogous to $E_0$, thus, in order to find an estimate for $T_s$, we first invert the imagery (see Section \ref{s:inversion}) assuming an unaccelerated population, which provides a map of $E_0$. Figure \ref{f:srce}, panels a -- b, show this map of $E_0$ and the total energy flux, $Q_p$, for our February 10, 2023 event.

\begin{figure}
    \centering
    \includegraphics[width=\linewidth]{03_determining_srce.pdf}
    \caption{Steps in determining the source region characteristic energy. (a) The total precipitating electron energy flux, $Q_p$, inverted assuming unaccelerated Maxwellian energy spectra. (b) The characteristic energy, $E_0$, inverted assuming unaccelerated Maxwellian energy spectra. (c) $E_0$ filtered by removing the top 40th percentile of $Q_p$ and the lower 30th percentile of the red line emissions. (d) Histograms of data in panels b (orange) and c (light blue) along with Gaussian magnitude fits (black and red respectively) and their peaks (dashed). Data source: \citeA{dasclink}.}
    \label{f:srce}
\end{figure}

Next, assuming that $U_a$ vanishes outside of discrete auroral arcs, we filter the arc region out of this map of $E_0$ by removing pixels corresponding to the top 40th percentile of $Q_p$. We also remove the lower 30th percentile of the red emissions, as the inversion to $E_0$ performs sub-optimally for lower red intensities. This is shown in Figure \ref{f:srce}c. We then look at the histogram of the remaining $E_0$ values and fit a Gaussian magnitude distribution to it, the peak of which is selected as the source region characteristic energy. In this case, we have $T_s=490$ eV, as is shown in panel d. This panel also shows the unfiltered distribution which shows two distinct populations, suggesting different physics behind them---presumably that of the accelerated electrons and that of the unaccelerated precipitation.

The choice of top 40th and lower 30th percentiles for the energy flux and red emissions are determined by balancing low 95\% confidence ranges and high adjusted R-squared values when fitting the filtered data in Figure \ref{f:srce}d. The percentiles investigated for $Q_p$ range from the top 10th -- 90th and for the red emissions ranged from the lower 10th -- 40th. In this process, we find that the $T_s$ estimates we report have around a $\pm$10 -- 20\% precision.

We assume this value of $T_s$ to be constant over the relevant source region and use it in Equation \ref{e:acc_max}, with which we perform the inversion described in Section \ref{s:inversion}. This inversion now happens over a $(Q_p,U_a)$ parameter space, for a given $T_s$, instead of $(Q_p,E_0)$, when creating lookup tables, providing 2-D maps of $Q_p$ and $U_a$. This is all done for each of the six conjunction events. Reassuringly, we find that inversions of these six events done with either the $\phi_u(E)$ or $\phi_a(E)$ assumptions provide nearly identical maps of $Q_p$; however, as we will show, they imply quite different conductivity and current density distributions through the ionosphere.

The methodology described in this section, at the time of writing, has not been validated experimentally because the measurement data available in our conjunctions exclude electron spectrometry. We do analyze the different spectral assumptions, and the effect they have on auroral system science, to great lengths in Section \ref{s:precipresults}. That said, this methodology invites future studies applying these techniques with multi-spectral imagery data and data used in auroral acceleration region spectrometry \cite{marghitu2006,imajo2024}.

\subsection{Current Flux Tube Visualization}\label{s:fluxtubes}
In order to visualize current closure in GEMINI output data, we show flux tubes of electric current. GEMINI enforces $\nabla\cdot\mathbf{j}=0$, where $\mathbf{j}$ is the current density, which makes the usage of flux tubes as a visualization tool sensible. We have developed tools to generate current flux tubes starting at user-defined ellipses contained inside the GEMINI simulation volume. From these ellipses, a number of current vector streamlines are sourced, which, by definition, are tangent to $\mathbf{j}$ throughout the simulation volume. This ensures the current flux through such ellipses is equal to the flux through the orientable surface enclosed by the curve connected by the streamline endpoints. Current fluxes are calculated for tubes that meet flat exit surfaces and are compared against entry fluxes as a check for numerical error. Streamline endpoints that are too far apart, or that meet at a corner of the simulation volume, are locations where the flux tube splits into multiple tubes. In this case, the fluxes of each tube are provided separately. This method of visualization is part of the toolset available online \cite{aurogemlink}.

Figure \ref{f:fluxtubes} shows three example current flux tubes. This $425\times288\times384$ cell (up, east, north) magnetically aligned volume contains a GEMINI calculated 3-D current density from which the flux tubes are derived. In this paper, simulations are all located in the northern hemisphere and magnetic east, north, and up refer to a locally orthonormal basis with up being anti-parallel to the local magnetic field, east in the direction of increasing modified apex longitudes \cite{laundal2017}, and north completing the set. The simulation in Figure \ref{f:fluxtubes} is driven by a top-boundary map of FAC which is plotted at the bottom for visualization purposes. The colormap of FAC has red associated with the downward, parallel-to-$\mathbf{B}$ (in the northern hemisphere) current vector, also referred to hereinafter as return current (i.e. “red is return”). The blue represents the upward current (downward-moving electrons in the Northern hemisphere) where, often, the accelerated auroral electron precipitation is found. On the eastern wall, a central cut of electron density is plotted. The density perturbations, which are in most part the result of the top-boundary precipitation driver maps, govern the 3-D conductivity volume and thus, in part, the current closure. The black arrows plotted on the FAC map are a sparse sample of the GEMINI calculated electric field---the second aspect governing the current closure---and the yellow arrow is the imposed constant, background convection electric field, $\bar{\mathbf{E}}$. The pink lines indicate the FAC current data from, in this example, Swarm A and C, that are footpointed down to the top-boundary and plotted at the bottom as well (these form the basis of the replicated FAC map in red and blue).

\begin{figure}
    \centering
    \includegraphics[width=\textwidth]{04_fluxtube_example.pdf}
    \caption{Example of a current flux tube plot using a February 10 simulation. The top-boundary FAC driver is plotted at the bottom for visualization purposes. Similarly, a central cut of electron density is plotted at the eastern wall along with flux tube projections. The current flux tubes are color-coded for distinction purposes and start/end at the bold/thin black solid curves. The black dashed lines are their counterparts projected on top of the FAC map. The pink lines indicate FAC data from Swarm A (right track) and C (left track) with parallel being right. The black arrows are a sparse sample of the electric field calculated by GEMINI and the yellow arrow indicates the constant background electric field. Data sources: \citeA{swarmlink}, \citeA{superdarnlink}, and \citeA{simulationslink}.}
    \label{f:fluxtubes}
\end{figure}

The current flux tubes are color-coded for easy distinction. In this example, the red flux tube originates from an ellipse at the top-boundary inside the southernmost downward, return current sheet. It carries 1.4 kA down through the volume, splitting in three, finding its way out through the top-, south-, and east-boundary. The influx and outflux regions are outlined by bold and thin closed black curves, and shadows of these curves are projected to the bottom to visualize which portion of the FAC map they embody. For select flux tubes, a projection of the entire tube is also plotted on the eastern wall. The green flux tube has its user-defined ellipse in the upward current and is calculated in reverse. It carries around 0.2 kA from two sources on the western wall, combines into a single tube, and connects with the top-boundary. Lastly, the orange flux tube (also calculated in reverse) is sourced at the northern boundary and also connects to the upward FAC. Figure \ref{f:fluxtubes}, and similar figures in the remainder of this paper, display in- and out-fluxes to two decimal places and illustrate the degree of precision of the flux tube calculations. Most current flux tubes in this paper are precise up to one decimal place, with a few exceptions of more complex current flux tubes or ones with higher amperage ($>10$ kA).

\section{Conjunction Events}\label{s:events}
This study uses a total of six conjunction events ranging from February 10 to March 19, 2023, from the Swarm-over-Poker-2023 campaign (Poker Flat Research Range, AK). As a summary of the conjunctions used in this work, Figure \ref{f:event_summary} shows selected top-boundary simulation data-drivers for each of the six events: the total energy flux of the precipitating electrons, $Q_p$, the acceleration potential, $U_a$, and the replicated FAC maps, $j_\parallel$. Driver maps of $E_0$ or those of $j_\parallel$ using fewer than all available spacecraft are not shown. Also plotted are the primary and secondary boundaries used in the replication process (see section \ref{s:replication}) and the FAC data tracks themselves. In addition, Table \ref{t:events} displays information regarding which Swarm spacecraft are part of the conjunction, the activity levels, the PFISR and SuperDARN background flow estimates, and the rough peak values of the simulation top-boundary drivers for each event. The distance from Poker Flat to the nearest SuperDARN plasma flow estimate, $\bar{\mathbf{v}}_{SD}$, is denoted $d_{SD}$.

\begin{figure}
    \centering
    \includegraphics[width=\textwidth]{05_event_summary.pdf}
    \caption{Top-boundary drivers of conjunction events. (a) The total precipitating electron energy flux, $Q_p$, for the February 10, 9:51 UT event. (b) The acceleration potential, $U_a$, for the same event. (c) The replicated FAC map, $j_\parallel$, for the same event. (d-r) Same format for remaining events. Note that the respective colorbars change per event. The solid black feather plot indicates the Swarm FAC data tracks with right being parallel. Not all Swarm data tracks are within the simulation volume and are thus not shown. Data sources: \citeA{swarmlink} and \citeA{dasclink}.}
    \label{f:event_summary}
\end{figure}

\begin{table}
\caption{Summary of conjunction events with input map values$^a$.}
\centering
\begin{tabular}{llllllll}
\hline
\textbf{Event ID}                          & 1         & 2         & 3         & 4         & 5         & 6         \\
\textbf{Date}                              & Feb 10    & Feb 12    & Mar 4     & Mar 4     & Mar 14    & Mar 19    \\
\textbf{Time} (UT)$^b$                     & 9:51:27   & 10:22:11  & 7:30:12   & 10:13:49  & 6:49:07   & 8:23:30   \\ \hline
\textbf{MLT} (Hours)                       & 23.1      & 23.3      & 20.7      & 22.9      & 20.1      & 21.4      \\
\textbf{Region} (km)$^c$                   & $290\times182$ & $290\times189$ & $290\times126$ & $290\times225$ & $220\times126$ & $432\times291$ \\
\textbf{Swarm S/C}                         & A + C     & C         & C         & B         & A + C     & B         \\ \hline
\textbf{Ap}                                & 15        & 7         & 16        & 16        & 18        & 9         \\
\textbf{F10.7 (a)}$^d$ (s.f.u.)            & 208 (175) & 200 (175) & 182 (161) & 182 (161) & 138 (162) & 143 (162) \\ \hline
\textbf{$\bar{\mathbf{v}}_{SD}$} (m/s)$^c$ & -14, 29   & -170, -31 & -323, 269 & -45, 0    & -200, -9  & -494, 96  \\
\textbf{$d_{SD}$} (km)$^e$                 & 51        & 51        & 184       & 373       & 51        & 375       \\
\textbf{$\bar{\mathbf{v}}_{PF}$} (m/s)$^c$ & -343, 2   & -237, -17 & -         & -         & -418, -44 & 178, -68  \\ \hline
\textbf{$Q_p$} (mW/m$^2$)$^f$              & 10.0      & 2.3       & 32.3      & 4.1       & 5.8       & 31.3      \\
\textbf{$U_a$} (keV)$^f$                   & 5.8       & 1.9       & 5.4       & 2.9       & 3.0       & 8.5       \\
\textbf{$T_s$} (eV)                        & 490       & 580       & 800       & 860       & 240       & 680       \\
\textbf{$E_0$} (keV)$^f$                   & 4.2       & 1.4       & 4.0       & 2.3       &  -        & -         \\ \hline
\textbf{$j_\parallel$} ($\mu$A/m$^2$)$^f$  & -2.3, 2.0 & -0.7, 1.9 & -4.5, 3.8 & -1.1, 1.0 & -1.2, 2.8 & -1.9, 1.4 \\ \hline
\multicolumn{7}{l}{$^a$Variables $\bar{\mathbf{v}}_{SD}$, $\bar{\mathbf{v}}_{PF}$, $Q_p$, $U_a$, $T_s$, $E_0$, and $j_\parallel$ are defined in-text.}                   \\
\multicolumn{7}{l}{$^b$Times indicate the spacecraft crossing approximately through the simulation center.}\\
\multicolumn{7}{l}{$^c$Regions and flows are in GEMINI magnetic coordinates/components.}\\
\multicolumn{7}{l}{$^d$F10.7a is the 81-day F10.7 average centered around the event dates.}\\
\multicolumn{7}{l}{$^e$Distances from Poker Flat to nearest SuperDARN data point.}\\
\multicolumn{7}{l}{$^f$Values for $Q_p$, $U_a$, $E_0$, and $j_\parallel$ are 99\% quantiles of maps within a 10 cell border.}\\
\end{tabular}
\label{t:events}
\end{table}

\begin{figure}
    \centering
    \includegraphics[width=\textwidth]{06_superdarn_summary.pdf}
    \caption{SuperDARN convection maps of conjunction events. Panels a -- f represent event IDs 1 -- 6 (see Table \ref{t:events}). Purple boxes are approximately centered on Poker Flat, AK and are on the order of the simulation sizes. The bold black line is the Heppner-Maynard Boundary. The colormap shows the electric potential and the ``$+$'' and ``$-$'' symbols indicate the maximum and minimum potential points. Local magnetic midnight is at the bottom and the dusk side is left. Data source: \citeA{superdarnlink}.}
    \label{f:superdarn_summary}
\end{figure}

Not all events have PFISR data tracks available because either they are too far from their respective, chosen simulation regions, or the data are considered inadequate for our purposes. Also, not all of the events have a simulation using the unaccelerated assumption for $\phi(E)$. Determining plausible arc boundaries requires meticulous care and controls where the simulation boundaries are, which is why, for several conjunction events, the FAC data track(s) lie(s) just east or west of simulation region. In such cases, the arc boundaries are extrapolated to the data tracks. Following are brief synopses of each of the six conjunction events after which, in Section \ref{s:results}, we cover their simulation results.

\subsection{February 10, 9:51:27 UT}
Figure \ref{f:event_summary}a -- c: This event includes both Swarm A and C cutting through the center of the simulation around 47 km apart. It has a curved double arc precipitation pattern with each peaking around a total energy flux of $Q_p=10.0$ mW/m$^2$ and acceleration potential of $U_a=5.8$ keV. The precipitation is collocated with the FAC replication where the precipitating and return current sheets are between $j_\parallel=-2.3$ to $2.0$ $\mu$A/m$^2$. The PFISR convection flow data are positioned at the western edge of the simulation space and estimate a strong magnetic westward flow of $\bar{\mathbf{v}}_{PF}=(-343, 2)$ m/s. In contrast, SuperDARN estimates a nearly stagnant flow of $\bar{\mathbf{v}}_{SD}=(-14, 29)$ m/s. The Magnetic Local Time (MLT) is 23.1, however, as is shown in Figure \ref{f:superdarn_summary}a, the event occurs 3 -- 4 hours duskside of the two-cell split.

\subsection{February 12, 10:22:11 UT}
Figure \ref{f:event_summary}d -- f: This is a low flux, low energy, and generally inactive event with a Swarm A conjunction roughly 153 km west of the simulation space and with a PFISR data cut through the center. It has a single, blurry but straight arc of around $Q_p=2.3$ mW/m$^2$ and $U_a=1.9$ keV, with the FAC sheets ranging from $j_\parallel=-0.7$ to $1.9$ $\mu$A/m$^2$. Both PFISR and SuperDARN suggest a large westward flow of $\bar{\mathbf{v}}_{PF}=(-237, -17)$ and $\bar{\mathbf{v}}_{SD}=(-170, -31)$ m/s respectively. The MLT is 23.3---roughly 1 hour prior to the two-cell split.

\subsection{March 4, 7:30:12 UT}
Figure \ref{f:event_summary}g -- i: In contrast to the previous event, this one has an intense arc of $Q_p=32.3$ mW/m$^2$ and $U_a=5.4$ keV with a Swarm C crossing around 141 km eastward of the simulation space and FAC data of $j_\parallel=-4.5$ to $3.8$ $\mu$A/m$^2$. This arc has reasonable along-arc structure; the total energy flux ranges from its peak to around 20 mW/m$^2$ going from west to east. Unfortunately, this event does not have usable PFISR data, but SuperDARN shows a very strong northwestern flow of $\bar{\mathbf{v}}_{SD}=(-323, 269)$ m/s. This strong, skewed flow is the result of a skewed two-cell convection pattern determined by pyDARN v4.1 \cite{superdarn1} as shown in Figure \ref{f:superdarn_summary}c. The event's MLT is 20.7, but this convection pattern places it around 5 -- 7 hours before the two-cell split.

\subsection{March 4, 10:13:49 UT}
Figure \ref{f:event_summary}j -- l: This event, just under three hours later than the previous at an MLT of 22.9, has a straight double arc pattern at $Q_p=4.1$ mW/m$^2$ and $U_a=2.9$ keV with Swarm B an average of 94 km westward of the simulation. This event has $T_s=860$ eV, which is $60$ eV higher than 2.75 hours earlier (this change is within the $\pm$10 -- 20\% precision of estimated $T_s$ values), and the currents have now subsided down to $j_\parallel=-1.1$ to $1.0$ $\mu$A/m$^2$. Again, this event includes no PFISR data, while SuperDARN now estimates a stagnant flow of $\bar{\mathbf{v}}_{SD}=(-45,0)$ m/s. Compared to the previous event, Figure \ref{f:superdarn_summary}d shows a much subdued convection pattern with the two-cell split sitting right around local magnetic midnight.

\subsection{March 14, 6:49:07 UT}
Figure \ref{f:event_summary}m -- o: This event is distinct in that it has its precipitation collocated with downward, rather than upward, FAC. There is a down-up-down FAC sheet set ranging from $j_\parallel=2.8$ to $-1.2$ to $2.0$ $\mu$A/m$^2$ centered around a $Q_p=5.8$ mW/m$^2$, $U_a=3.0$ keV precipitation pattern. It is also the second event with both Swarm A and C conjunctions. Swarm A sits around 44 km east of the model space, while the Swarm C crossing is just inside at the northeastern corner, and the southernmost PFISR data point is located around 100 km west of the simulation. The direction of both the PFISR and SuperDARN convection flow estimates are very similar, however, the PFISR flow estimate of $\bar{\mathbf{v}}_{PF}=(-418, -44)$ m/s is around twice as strong as the SuperDARN estimate of $\bar{\mathbf{v}}_{SD}=(-200,-9)$ m/s. This 20.1 MLT event sits at around 2 hours duskside to the two-cell split.

\subsection{March 19, 8:23:30 UT}
Figure \ref{f:event_summary}p -- r: The last event, and the second Swarm B conjunction, is unaligned to magnetic latitudes and has strong precipitation with along-arc structure; the energy flux peaks at around $Q_p=31.3$ mW/m$^2$ and subsides to around 20 mW/m$^2$ at the eastern and western boundaries. The acceleration potential is the highest among our events, peaking at around $U_a=8.5$ keV and the FAC data range from around $j_\parallel=-1.9$ to $1.4$ $\mu$A/m$^2$. PFISR cuts through the center and estimates a flow of $\bar{\mathbf{v}}_{PF}=(178, -68)$, where SuperDARN estimates $\bar{\mathbf{v}}_{SD}=(-494, 96)$ m/s. The MLT is 21.4, however, Figure \ref{f:superdarn_summary}f shows a multi-cell convection pattern which gives a relatively nonstandard context.

\section{Simulation Results}\label{s:results}
The six conjunction events are each simulated multiple times, iterating through different parameters, allowing the simulations to be systematically compared. This highlights and isolates the relevant physics involved and allows us to study sensitivities to these parameters. Table \ref{t:comparisons} provides the list of simulation comparisons covered in this paper (and its Supporting Information), labeled IDs I-XI, where individual simulations are denoted Ia, Ib, IIa, and so on. The comparisons are divided into three categories of feature permutations: (1) background convection flow and its source, (2) the assumption of unaccelerated versus accelerated Maxwellian precipitation spectra, and (3) single versus double spacecraft replications, highlighting along-arc FAC structure.

\begin{table}
\caption{Summary of event comparisons$^a$}
\centering
\begin{tabular}{ccccccccccc}
\textbf{Category} & \textbf{ID} & \textbf{Datetime} (UT)$^b$ & \multicolumn{2}{c}{\textbf{BG flow} (m/s)$^c$} & \multicolumn{2}{c}{\textbf{BG source}} & \multicolumn{2}{c}{\textbf{Acc.}} & \multicolumn{2}{c}{\textbf{Swarm}} \\ \hline
                         &               &               & \textbf{a}  & \textbf{b}  & \textbf{a} & \textbf{b} &            &            &            &            \\ \hline
                         & \textbf{I}    & Feb 10, 9:51  & (-14, 29)   & (-343, 2)   & SD         & PFISR      & Y          & -          & AC         & -          \\
\textbf{Background}      & \textbf{II}   & Feb 12, 10:22 & (-170, -31) & (-237, -17) & SD         & PFISR      & Y          & -          & C          & -          \\
\textbf{flow}            & \textbf{III}  & Mar 4, 7:30   & (-323, 269) & (0, 0)      & SD         & None       & Y          & -          & C          & -          \\
                         & \textbf{IV}   & Mar 14, 6:49  & (-200, -9)  & (-418, -44) & SD         & PFISR      & Y          & -          & AC         & -          \\
                         & \textbf{V}    & Mar 19, 8:23  & (-494, 96)  & (178, -68)  & SD         & PFISR      & Y          & -          & B          & -          \\ \hline
                         &               &               &             &             &            &            & \textbf{a} & \textbf{b} &            &            \\ \hline
                         & \textbf{VI}   & Feb 10, 9:51  & (-14, 29)   & -           & SD         & -          & Y          & N          & AC         & -          \\
\textbf{Accelerated vs.} & \textbf{VII}  & Feb 12, 10:22 & (-170, -31) & -           & SD         & -          & Y          & N          & C          & -          \\
\textbf{unaccelerated}   & \textbf{VIII} & Mar 4, 7:30   & (-323, 269) & -           & SD         & -          & Y          & N          & C          & -          \\
                         & \textbf{IX}   & Mar 4, 10:14  & (-45, 0)    & -           & SD         & -          & Y          & N          & B          & -          \\ \hline
                         &               &               &             &             &            &            &            &            & \textbf{a} & \textbf{b} \\ \hline
\textbf{Along-arc}       & \textbf{X}    & Feb 10, 9:51  & (-14, 29)   & -           & SD         & -          & Y          & -          & AC         & A          \\
\textbf{structure}       & \textbf{XI}   & Mar 14, 6:49  & (-200, -9)  & -           & SD         & -          & Y          & -          & AC         & A          \\ \hline
\multicolumn{11}{l}{$^a$Comparisons are labeled I-XI with individual simulations labeled Ia, Ib, IIa, etc.}                                                        \\
\multicolumn{11}{l}{$^b$Times indicate the spacecraft crossing approximately through the simulation center.}                                                       \\
\multicolumn{11}{l}{$^c$PFISR and SuperDARN background flows are in GEMINI magnetic east/north components.}                                                        \\
\end{tabular}
\label{t:comparisons}
\end{table}

Each simulation has $425\times288\times384$ cells in the magnetic up, east, and north directions respectively and a duration of 60 seconds with static drivers. The altitudinal extent is 80 -- 507 km, with cell heights of 0.3 -- 10 km respectively, and the magnetic east/north extents are given in Table \ref{t:events} and Figure \ref{f:event_summary}. Horizontal cell dimensions settle at 700 -- 1400 m in the magnetic east direction, and 238 -- 700 m in the magnetic north direction. Unless otherwise stated, all simulations default to SuperDARN background flow estimates, accelerated Maxwellian precipitation, and FAC replication using maximal data tracks. The simulations can be found online \cite{simulationslink}.

\subsection{Background Flow \& Electric Field}\label{s:bef}
There are two factors which dictate the existence of closure currents: (1) the Pedersen and Hall conductivities, and (2) the strength of the electric field. The conductivities require enhanced ionization at closure altitudes which is largely dictated by electron precipitation---enhanced energy fluxes, $Q_p$, increase the conductivity overall, while stronger acceleration potentials, $U_a$, give preference to Hall over Pedersen closure. Adding to this, spatial structure in the precipitation means that these conductivities have 3-D structure, affecting current closure in all directions. The magnitude of the electric field, however, dictates the magnitude of closure currents overall. We argue that strong electric fields can render the need for Hall closure to be negligible. We therefore begin by looking at comparisons of simulations that have different background electric field assumptions.

Figure \ref{f:feb10} shows three view angles of the results for Simulation Ia, referenced in Table \ref{t:comparisons}, where Section \ref{s:fluxtubes} explains the format of this figure. It uses FAC data from Swarm A and C, the accelerated Maxwellian precipitation assumption, and a background plasma flow estimate from SuperDARN. In this first example, the background flow of $\bar{\mathbf{v}}_{SD}=(-14, 29)$ m/s amounts to a constant background electric field of 1.6 mV/m directed roughly 26 degrees north-of-east (geomagnetic).

\begin{figure}
    \centering
    \includegraphics[width=\textwidth]{07_gemini_output_20230210_35487_09_SD_AM_AC_1.pdf}
    \caption{Isometric (a), side (b), and top (c) view of the GEMINI results for Simulation Ia. For plot details, see Section \ref{s:fluxtubes}. Data sources: \citeA{swarmlink}, \citeA{superdarnlink}, and \citeA{simulationslink}.}
    \label{f:feb10}
\end{figure}

What follows are descriptions of three of our five comparisons (see Table \ref{t:comparisons}) that outline the sensitivity of auroral current closure to the constant background electric field, $\bar{\mathbf{E}}$, around which GEMINI solves current continuity and Ohm's law for $\mathbf{E}=\bar{\mathbf{E}}+\delta\mathbf{E}$. The remaining comparisons, along with their associated figures and descriptions, can be found in the Supporting Information. For a more complete 3-D visualization of current closure, the Supporting Information also contains isometric views, along with the side and top views, for all simulations.

\subsubsection{Comparison I: Background Flow}
In Comparison I, we compare and contrast the use of SuperDARN derived background flow against using the PFISR observed background flow. Figure \ref{f:feb10} illustrates three current closure paths of Simulation Ia, which assumes the SuperDARN background flow, and shows the complexity of current closure in a reasonably typical discrete auroral arc system. The red current flux tube, carrying 1.4 kA, starts at the center of southernmost return current sheet and rotates to closure currents at an altitude range of 110 -- 150 km. The bulk of the current continues northward, however, 0.2 kA exits through the southern boundary and $<0.1$ kA exits through the eastern boundary. Focusing on the remaining 1.2 kA, panel c shows that this segment opens up to the northeast, aligning the tube with the electric field at first, i.e. Pedersen closure. The weak strength of the electric field, however, renders the Pedersen closure ineffective and requires the tube to traverse through lower altitudes to find sufficient paths for closure. At these lower altitudes, the Hall currents dominate and thus this portion of the tube rotates perpendicular to the electric field. This increases the length it has to travel while crossing into the upward FAC region and stretches the overall current closure morphology in the along-arc direction. The portion which exits through the eastern wall, presumably, would follow this same pattern somewhere outside the simulation volume, but this is speculative. More notably, however, the remaining unclosed portion on the southern part of the tube traverses southward, but this is for the same reason: the tube rotates in the Hall layer looking for upward FAC somewhere outside the simulation.

The green tube is sourced from the western boundary with two ends, both carrying around 0.1 kA, which combine into a single, 0.2 kA upward segment of the tube closing in between the two precipitation current sheets. Panel a shows how they cling to the higher density, i.e. higher conductivity, regions caused by the double-arc precipitation; they wrap around these density enhancements in the northward direction following the local electric field.

The orange tube is sourced from the northern boundary with 0.5 kA and travels southward, somewhat aligned to the electric field, before it hits an electric field convergence. Thus, to avoid going against the electric field, the flux tube lowers in altitude, in search of Hall conductivity, and abruptly turns to the east. This outlines the self-consistency aspect of the nature of auroral current closure---the flux tube (a) lowers in altitude where (b) the density is higher, (c) the electric field converges, and (d) the Hall conductivity increases allowing for an eastward turn, all spatially coincident. Finally, the tube further rotates to gain just enough Pedersen current, and hence altitude, to allow for a connection with the upward FAC sheet. This current flux tube, along with the previous two, highlights a set of 3-D considerations needed when trying to understand current closure morphology, and thus MIT coupling. This is especially true when Hall currents are required in this closure, which is the case for Simulation Ia, given its weaker electric field.

In contrast, Figure \ref{f:compare-feb10-ebg} shows three current flux tubes for Simulation Ib (panels c -- d) that capture the same FAC regions, whether at the start or end of each tube. The only change here is that the simulation now assumes the PFISR derived constant background flow of $\bar{\mathbf{v}}_{PF}=(-343, 2)$ m/s, which amounts to 17.2 mV/m directed nearly northward compared to the northeasterly 1.6 mV/m from Simulation Ia (panels a -- b). This larger background flow drastically changes the current closure morphology of all three flux tubes. Given the tenfold increase in the electric field magnitude, on top of a more direct Pedersen pathway across the arcs, the Pedersen closure has become significantly more effective. Panels a and c show an increase in closure altitudes of 110 -- 150 to 130 -- 180 km, which means the Hall layer is virtually untouched by these Simulation Ib closure patterns. Panel d solidifies this idea, as all three tubes follow the electric fields almost directly. This outlines the ability of the background electric field, $\bar{\mathbf{E}}$, to actively drive auroral arc systems in conjunction with the top-boundary map of $j_\parallel$.

\begin{figure}
    \centering
    \includegraphics[width=\textwidth]{08_joule_heating_20230210_35487_09_SD_AM_AC_v_PF_AM_AC.pdf}
    \caption{Comparison I (February 10, 9:51 UT): Top and side views of Simulation Ia with SuperDARN derived background flow (a, b) versus Simulation Ib with PFISR derived background flow (c, d) along with height-integrated Joule heating for Simulation Ia (e) and Ib (f). For plot details, see Section \ref{s:fluxtubes}. Data sources: \citeA{swarmlink}, \citeA{superdarnlink}, \citeA{pfisrlink}, and \citeA{simulationslink}.}
    \label{f:compare-feb10-ebg}
\end{figure}

To emphasize the sensitivity to the background electric field from the perspective of energy dissipation, panels e -- f of Figure \ref{f:compare-feb10-ebg} show the height-integrated Joule heating for Simulations Ia -- b respectively. They show the extent to which this auroral arc system can be an electrostatic load, and how $\bar{\mathbf{E}}$ can change this greatly; aside from having an order-of-magnitude higher electric field strength, Simulation Ib also closes mostly in Pedersen currents---parallel to the electric field---both facts favoring higher $\mathbf{j}\cdot\mathbf{E}=\sigma_P|\mathbf{E}|^2$ values throughout. Not only does this increase the Joule heating for Simulation Ib, it also relocates a bulk portion of it equatorward of the precipitation.

The simulations in the next comparison, Comparison III, have a similar disparity in electric field strengths, yet both have higher FAC requirements, dictating a larger need for closure currents. However, they both also have more precipitation; a factor which partially fulfills this need for additional closure.

\subsubsection{Comparison III: Background Flow}
The conjunction event for Comparison III, unfortunately, occurs too far from the PFISR field-of-view and therefore has no PFISR-deduced background flow estimate. Nevertheless, Figure \ref{f:compare-mar4a-ebg} demonstrates the sensitivity to the choice of background flow by looking at Simulation IIIa, where the SuperDARN derived background convection amounts to 21.0 mV/m directed 40 degrees east of north (first row), and comparing it to Simulation IIIb, which has the background flow set to zero, as there is no estimate for it (second row). As explained in Section \ref{s:intro}, having zero background electric field amounts to assuming most of the top-boundary FAC, $j_\parallel$, comes from electric fields caused by local polarization, $\delta\mathbf{E}$, alone. This comparison shows how much such an assumption affects current closure. Note that, with $|\bar{\mathbf{E}}|=0$, for illustration purposes, the electric field label (black here) indicates the magnitude of the GEMINI calculated electric field vector shown nearest the label.

\begin{figure}
    \centering
    \includegraphics[width=\textwidth]{09_electrojet_current_20230304_27012_09_SD_AM_xC_v_NB_AM_xC.pdf}
    \caption{Comparison III (March 4, 7:30 UT): Top and side views of Simulation IIIa with SuperDARN derived background flow (a, b) versus Simulation IIIb with no derived background flow (c, d) along with height-integrated Joule heating for Simulation IIIa (e) and IIIb (f). (g, h) North-up slices of the magnetic eastward current component for Simulations IIIa -- b respectively taken at 50 km west from center with the start curves of their respective orange flux tubes (solid black). For plot details, see Section \ref{s:fluxtubes}. Data sources: \citeA{swarmlink}, \citeA{superdarnlink}, and \citeA{simulationslink}.}
    \label{f:compare-mar4a-ebg}
\end{figure}

Comparisons I and III both look at simulations with an order-of-magnitude difference in their electric field strengths and both cover conjunction events whose acceleration potentials peak at around $U_a=5$ keV. Comparison III, however, has the precipitation energy flux more than triple, and FAC requirements roughly double, with respect to Comparison I (see Tables \ref{t:events} and \ref{t:comparisons}). This creates a higher need for current closure---a need partially fulfilled by increased conductivity at all altitudes and the strong electric field strength. Hindering these needs, however, is the less direct path for Pedersen closure given the roughly 40 degree angle at which the electric field crosses the arc in Simulation IIIa. The combination of these features allows us to look at how the sensitivity to electric field strength is affected by a different arrangement of auroral arc parameters.

Simulation IIIb, with $\bar{\mathbf{E}}=0$, depicts a typical perspective of discrete aurora \cite{marghitu2012}---an arc-aligned line of diverging electric field at the downward current sheet, and a converging one at the upward current sheet, as suggested by Equations \ref{e:fac} -- \ref{e:delta_fac}. In this simulation, this is the result of the absence of a background electric field causing current continuity and Ohm's law to be solved with electric fields from local polarization alone. The red flux tube in Simulation IIIb digs deep into the Hall layer while closing and is forced to split when bottoming out. This causes 0.4 kA to exit through the southern wall, 0.7 kA through the top-boundary, and $<0.1$ kA through the eastern wall. (Note that this tube loses around 0.2 kA throughout its path which is a result of edge effects at the eastern wall). In contrast, the order-of-magnitude higher electric field in Simulation IIIa means that its red flux tube carries that 1.4 kA from the return current sheet across to the precipitation sheet all throughout Pedersen altitudes and, thus, its closure is directed almost completely in the electric field direction. Contrarily, the green flux tubes for both simulations close largely with Pedersen currents given that their ends are relatively near one-another. Even though the green Simulation IIIb flux tube finds its closing currents at lower altitudes, it is still mostly dominated by Pedersen conductivity throughout; only the bottom apex of this tube veers to the across-$\mathbf{E}$ direction.

Morphologically speaking, the most striking difference between Simulations IIIa and IIIb lies in their connection to the electrojet current. Figure \ref{f:compare-mar4a-ebg}, panels g -- h, for Simulations IIIa -- b respectively, show a slice of the magnetic eastward component of $\mathbf{j}$ taken at 50 km west-from-center, along with the intake ends of their respective orange flux tubes in panels a -- d. With its stronger electric field, Simulation IIIa has a much higher electrojet current. The $\sim13$ mV/m eastward component of the Simulation IIIa electric field gives rise to a substantial display of the Cowling mechanism  enhancing the auroral electrojet current \cite{cowling1932,amm2011} as indicated by Figure \ref{f:compare-mar4a-ebg}, panels b and g. This makes this auroral arc system closely resemble a 3-D version of the description from Section 4 by \citeA{marghitu2012}: A ``thick uniform 2-D arc'' whose current closure is separated into a thin Pedersen and Hall layer as shown by \citeA{fujii2012}. Expanding on this description, here we show how current flux tubes can navigate around each other in a coherent and self-consistent way by venturing into the 3-D perspective.

Given the more complex shape of the orange flux tube in Simulation IIIa, Figure \ref{f:mar4a} shows the isometric view of the simulation results, in addition to the side and top views in panels a -- b from Figure \ref{f:compare-mar4a-ebg}. Here we see the almost helical shape of the orange Simulation IIIa current flux tube, resembling that of Example 3 by \citeA{mallinckrodt1985} but in 3-D. This tube captures 13.4 kA of the electrojet current, while its Simulation IIIb counterpart carries around 0.5 kA. Both intake ellipses have the same dimensions and are centered on their respective peaks of magnetic eastward currents slices. The simulations both have the same relatively strong precipitation arc ($Q_p=32.3$ mW/m$^2$, $T_s=800$ keV) around 10 -- 20 km north, resulting in a high amount of impact ionization at relatively lower altitudes. This provides plenty of Hall conductivity and, thus, has both simulations susceptible to strong electrojet currents. These currents, however, are still proportional to the electric field strength which is why the order-of-magnitude increase in electric field results in a similarly increased electrojet current.

\begin{figure}
    \centering
    \includegraphics[width=\textwidth]{10_swop_20230304_27012_09_SD_AM_xC_ISO_P_000_1.pdf}
    \caption{Isometric view of the GEMINI results for Simulation IIIa. For plot details, see Section \ref{s:fluxtubes}. Data sources: \citeA{swarmlink}, \citeA{superdarnlink}, and \citeA{simulationslink}.}
    \label{f:mar4a}
\end{figure}

The enhanced electrojet current in Simulation IIIa does not partake in parallel current closure, whereas the Simulation IIIb electrojet current \textit{is} required in the coupling of magnetospheric currents. As mentioned before, the lower electric field strength overall renders all closure currents less effective, hence the FAC has to rely on enhanced conductivity---Pedersen and then Hall---to connect. Naturally, Simulation IIIa is a more energetic configuration in terms of Joule heating; the integrated Joule heating peaks are at around 26.6 mW/m$^2$ and 0.17 mW/m$^2$ for Simulations IIIa -- b respectively, as shown in panels e -- f of Figure \ref{f:compare-mar4a-ebg}. This is consistent with the order-of-magnitude difference in electric field strengths, given the $|\mathbf{E}|^2$ relationship. Given that Hall currents are dissipationless \cite{kaeppler2012}, Simulation IIIb is thus able to rely on the electrojet currents for closure instead. In Simulation IIIa, the electrojet largely is assumed to follow the global-scale convection pattern D-shaped Hall currents instead, and has negligible involvement in auroral FAC closure. This electrojet current does, however, connect with Pedersen currents as shown in Figures \ref{f:compare-mar4a-ebg}a -- b and \ref{f:mar4a}.

As with Comparison I (as well as II and IV in the Supporting Information), here, yet again, we see that a sufficiently large background electric field has FAC close with Pedersen currents, and thus in the direction of the electric field. Even with the less-direct Pedersen pathway for closure and the higher FAC requirements, the strong electric field and relatively large precipitation energy flux provides sufficient conductivity at higher altitudes and renders Pedersen closure to be the dominant method in MIT coupling for Simulation IIIa. Furthermore, the Simulation IIIb solution features a distinct $\nabla\cdot\mathbf{E}$ signature. In Simulation IIIa, however, this signature is masked by the its large background electric field (compared to no background field in Simulation IIIb). This emphasizes the dominance of the $\nabla\Sigma_{P,H}$ terms in balancing the FAC map for auroral systems with large electric fields. Next, we move onto Comparison V whose simulations both have strong electric fields, yet in severely different orientations.

\subsubsection{Comparison V: Background Flow}
Simulation Va assumes a background electric field of 25.2 mV/m directed 11 degrees east of magnetic north as estimated by SuperDARN. In the almost complete opposite direction to this, Simulation Vb has PFISR estimate 9.5 mV/m directed 21 degrees west of south. This results in drastic differences in both current closure morphology and Joule heating, as depicted by Figure \ref{f:compare-mar19-ebg}. Both the red and green flux tubes completely flip directions in their current closure. The green flux tube, in its attempt to connect to the broad, primary precipitation current sheet, changes from sourcing its roughly 5.3 kA from the southwest corner in Simulation Va, to doing so from the northern end in Simulation Vb. The red flux tube, closing the southern primary downward current sheet, simply flips direction by following the electric field, and, interestingly, in both simulations it ends up skipping over an adjacent, lesser downward current sheet when closing its 0.8 kA.

\begin{figure}
    \centering
    \includegraphics[width=\textwidth]{11_electrojet_current_20230319_30210_09_SD_AM_xB_v_PF_AM_xB.pdf}
    \caption{Comparison V (March 19, 8:23 UT): Top and side views of Simulation Va with SuperDARN derived background flow (a, b) versus Simulation Vb with PFISR derived background flow (c, d) along with height-integrated Joule heating for Simulation Va (e) and Vb (f). (g, h) Central north-up slices of the magnetic eastward current component for Simulations Va -- b respectively with the start curves of their respective orange flux tubes (solid black). For plot details, see Section \ref{s:fluxtubes}. Data sources: \citeA{swarmlink}, \citeA{superdarnlink}, \citeA{pfisrlink}, and \citeA{simulationslink}.}
    \label{f:compare-mar19-ebg}
\end{figure}

As shown in panels g -- h of Figure \ref{f:compare-mar19-ebg}, the orange flux tube, like in Comparison III, captures the electrojet current for both Simulations Va -- b. (Here, the user-defined ellipse sits at 0 km east and the tube is calculated in both directions.) As expected, this flux tube also flips its orientation, with the current flowing from west-to-east in Simulation Va, and east-to-west in Simulation Vb. However, like in Comparison III (Figure \ref{f:compare-mar4a-ebg}) but to a lesser extent, the weaker electric field strength in Simulation Vb requires the need of this electrojet current to help close some of the FAC, 0.3 kA in this case. The 62\% weaker field also has reduced this Hall current flux tube by about 52\%.

As in Comparisons I and III, the height-integrated Joule heating shown in panels e -- f of Figure \ref{f:compare-mar19-ebg} varies roughly in proportion to the electric field strength squared. One notable difference, however, lies in the tapering off of this Joule heating in simulation Vb. This indicates that the western boundary of this simulation relies more on Hall closure; an idea supported by the electrojet usage in FAC closure depicted by the orange flux tube in panels c -- d. Regardless of the reasoning behind this, Comparison V has shown that a mere directional change in the background electric field can create different dissipation characteristics of an auroral system. Moreover, Comparison V highlights how the direction of the background electric field completely changes the connectivity of a given map of FAC. It is tempting to assume that a precipitation current sheet connects with its closest adjoining return current sheet, but as we have shown here, knowledge of the global-scale convection has considerable influence when it comes to FAC connectivity.

\subsubsection{Summary: Background Flow}
Auroral arc systems are very sensitive to the electric field in matters of current closure. Given that there are many self-consistent solutions for $\mathbf{E}$ in Equation \ref{e:fac} that can be considered geophysical, we have shown here that it is crucial to get a good estimate of the global background flow in order to properly interpret behavior at auroral arc scales. In terms of simulation confidence, we can have more trust in simulations whose sources for background electric field estimates agree, such as Comparisons II and IV. However, when attempting to best understand the auroral arc system pertaining to a particular conjunction event, more certainty is needed for systems like those shown in Comparisons I, III, or V. Future conjunction campaigns will therefore benefit greatly from dedicated, multi-platform observations of large-scale convection flow---observations of comparable importance to in situ FAC measurements.

Previous studies have looked at how the strength and direction of the ionospheric electric field alters auroral current closure by systematically dividing the field-aligned, Pedersen, and Hall currents into primary and secondary components \cite{amm2011,fujii2011,fujii2012}. In this work, we look at an undivided, 3-D current connectivity in which primary and secondary polarization currents blend into a set of individual current flux tubes. The reader is encouraged to compare and contrast these two perspectives and recognize the sub-current mechanisms within the flux tubes themselves.

\subsection{Electron Precipitation Spectra}\label{s:precipresults}
As discussed in Section \ref{s:precip}, the choice of precipitating electron energy spectra can affect the impact ionization rate at different altitudes; an unaccelerated Maxwellian profile, Equation \ref{e:unacc-max}, often carries an erroneous high-energy tail which overestimates the \textit{E}-region density enhancement from electron precipitation. Moreover, choosing to use unaccelerated Maxwellian spectra in inverting multi-spectral imagery results in a de facto source region characteristic energy equal to the accelerating potential drop, i.e. $T_s=U_a=E_0$. This is not unlike how a choice of $\bar{\mathbf{E}}=0$ carries hidden assumptions about $j_\parallel$ (see Equations \ref{e:bar_fac} -- \ref{e:delta_fac}). With an accelerated Maxwellian profile, Equation \ref{e:acc_max}, we estimate $T_s$ prior to multi-spectral image inversion which allows for much ``colder'' source populations and, we argue, more geophysical precipitating electron modeling.

Below are two comparisons which look at how decoupling the source region characteristic energy from the auroral acceleration potential changes auroral current closure. See Table \ref{t:comparisons} for details on these comparisons and Table \ref{t:events} for the (peak) values for $Q_p$, $U_a$, $T_s$, and $E_0$.

\subsubsection{Comparison VI: Precipitation Spectra}
Returning back to the February 10 conjunction event from Comparison I, Figure \ref{f:compare-feb10-acc} depicts Comparison VI which looks at the differences between Simulation VIa (also named Ia) with the accelerated precipitation spectra assumption given by $\phi_a(E)$ (first row), and Simulation VIb which assumes $\phi_u(E)$ instead (second row). The first feature to point out is the central, north-up electron density slices shown in panels a and c: simulation VIa has both precipitation arc induced density enhancements tucked above 100 km in altitude, while the use of $\phi_u(E)$ in Simulation VIb has these same two arcs increasing their electron density enhancements to the bottom of the simulation volume. This limits the closure paths of current flux tubes in Simulation VIa, compared to Simulation VIb.

\begin{figure}
    \centering
    \includegraphics[width=\textwidth]{12_joule_heating_20230210_35487_09_SD_AM_AC_v_SD_UM_AC.pdf}
    \caption{Comparison VI (February 10, 9:51 UT): Top and side views of Simulation VIa with accelerated Maxwellian electron precipitation (a, b) versus Simulation VIb with unaccelerated Maxwellian electron precipitation (c, d) along with height-integrated Joule heating for Simulation VIa (e) and VIb (f). For plot details, see Section \ref{s:fluxtubes}. Data sources: \citeA{swarmlink}, \citeA{superdarnlink}, and \citeA{simulationslink}.}
    \label{f:compare-feb10-acc}
\end{figure}

Panels a and c show that all three current flux tubes are squished to higher altitudes in Simulation VIa, compared to Simulation VIb, forcing them to have a preference of Pedersen, over Hall, current closure. Panels b and d show how this affects the curvature of the flux tubes from a topside view. The red flux tube in Simulation VIa extends significantly further north given its preferred direction of that of the electric field, whereas the Simulation VIb red flux tube stays more parallel to the arc, traveling perpendicular to the electric field. Note that, despite the morphology being more along-arc, this flux tube does not extend much further east compared to the one in Simulation VIa, as it is able to capture higher upward FAC densities in this direction. Subsequently, its end region has a smaller overall area needed to capture 1.5 kA of upward FAC. Similar to the red flux tubes, the Simulation VIb orange current flux tube travels more often in the direction perpendicular to $\mathbf{E}$ compared to its Simulation VIa tube, again because the flux tube is able to traverse lower altitudes.

Panels e -- f show that, in this instance, the height-integrated Joule heating increases by around 30\% with the unaccelerated, over the accelerated, Maxwellian precipitation assumption. This can be counterintuitive when considering the Pedersen closure preference of Simulation VIa. Looking at panels b and d of Figure \ref{f:compare-feb10-acc}, however, tells us that the electric fields (black arrows) surrounding the arcs are higher in strength with the unaccelerated assumption which, evidently, is consistent with an increase in Joule heating. Ultimately, along with having the same FAC and background electric field drivers, both simulations have near identical maps of total precipitating electron energy, even thought their imagery inversions assume two different spectral shapes. This implies that the altitudinal distribution of impact ionization alters the energy accounting, and thus the electric load characteristic of this auroral arc system.

Both Simulations VIa -- b assume the relatively weak, SuperDARN derived $|\bar{\mathbf{E}}|=1.6$ mV/m, which makes them more susceptible to changes in the Hall closure layer as we have shown in Section \ref{s:bef}. Adding to this susceptibility, the precipitation arcs have relatively high values of $U_a=5.8$ and $E_0=4.2$ keV respectively. This deposits the impact ionization to lower altitudes, rendering the Hall layer more important still. This, along with the altered Joule heating, puts emphasis on the energy distribution shape of precipitating electrons in such auroral systems. Up next, we look at Comparison VIII whose simulations have both much stronger background electric fields and significantly higher total precipitation energy flux, which, along this line of reasoning, implies both Simulations VIIIa -- b are less reliant on Hall closure in MIT coupling.

\subsubsection{Comparison VIII: Precipitation Spectra}\label{s:viii}
The precipitation arcs in Simulation VIIIa (also named IIIa) assume an acceleration potential peaking at around $U_a=5.4$, and the characteristic energy for the arcs in Simulation VIIIb reach around $E_0=4.0$ keV---similar to the values from Comparison VI. However, with respect to Comparison VI, the simulations in Comparison VIII have more than three times the energy flux, $Q_p=32.3$ mW/m$^2$, background electric fields that are more than 13 times stronger, $|\bar{\mathbf{E}}|=21.0$ mV/m, and FAC sheets whose magnitudes around double, $|j_\parallel|=3.8$ -- $4.5$ $\mu$A/m$^2$. Additionally, at $T_s=800$ eV, the source region characteristic energy for simulation VIIIa also nearly doubles that of Simulations VIa. Figure \ref{f:compare-mar4a-acc} shows how unaccelerated Maxwellian precipitation at these more energetic parameters compares to accelerated Maxwellian precipitation.

\begin{figure}
    \centering
    \includegraphics[width=\textwidth]{13_electrojet_current_20230304_27012_09_SD_AM_xC_v_SD_UM_xC.pdf}
    \caption{Comparison VIII (March 4, 7:30 UT): Top and side views of Simulation VIIIa with accelerated Maxwellian electron precipitation (a, b) versus Simulation VIIIb with unaccelerated Maxwellian electron precipitation (c, d) along with height-integrated Joule heating for Simulation VIIIa (e) and VIIIb (f). (g, h) North-up slices of the magnetic eastward current component for Simulations VIIIa -- b respectively taken at 50 km west from center with the start curves of their respective orange flux tubes (solid black). For plot details, see Section \ref{s:fluxtubes}. Data sources: \citeA{swarmlink}, \citeA{superdarnlink}, and \citeA{simulationslink}.}
    \label{f:compare-mar4a-acc}
\end{figure}

By proxy of the electrojet currents shown in panels g -- h, the unaccelerated Maxwellian precipitation deposits ionization to both lower altitudes---around 6 km lower compared to Simulation VIIIa---and to a larger altitudinal range given the nearly four times higher energy spread of the unaccelerated energy spectra. We focus on these electrojet currents by looking at both orange flux tubes, which capture similar values of 13.1 and 14.8 kA for Simulations VIIIa -- b respectively. As before, the density volume resulting from the accelerated Maxwellian assumption is restricted to above around 100 km, forcing the respective orange tube to take on more Pedersen current. This means the electrojet in Simulation VIIIa veers to the northeast, directed toward the electric field (see panel b). The orange flux tube in Simulation VIIIb, being overall at lower altitudes, travels more easterly, staying relatively orthogonal to the electric field.

As shown in panels a -- b, the red flux tube in Simulation VIIIa takes advantage of the energy deposition at higher altitudes and the large electric field strength, and finds closure through Pedersen alone. In Simulation VIIIb, however, only around 0.7 of the 1.4 kA is able to connect with the FAC, while the remainder exits through the eastern boundary. Interestingly, the existence of the electrojet current in Simulation VIIIb appears to push the red flux tube away from the highest densities, subsequently squeezing it to lower altitudes. The green flux tube, having to travel a shorter horizontal distance compared to the other tubes, remains at altitudes where the $\phi_u(E)$ versus $\phi_a(E)$ assumption matters much less, and so it barely changes its morphology and amperage across the two simulations.

Panels e -- f of Figure \ref{f:compare-mar4a-acc} show a band of enhanced Joule heating just equatorward of the precipitating arc in both simulations, yet Simulation VIIIa has this band peak at around 26.6 mW/m$^2$, while Simulation VIIIb peaks closer to 40.0 mW/m$^2$---around a 50\% increase. Between the two simulations, the Pedersen current density remains fairly similar; it is the significantly varying Hall current density that creates the different morphologies (see Figure \ref{f:joule_viii}, panels a -- d). This points to the electric field strength; in Simulation VIIIa there is a band of enhanced across-arc electric field collocated with the Joule heating and peaks at around 20 mV/m, while the same is true for Simulation VIIIb except that it peaks around 40 mV/m (see panels e -- f).

\begin{figure}
    \centering
    \includegraphics[width=\textwidth]{14_joule_heating_viii.pdf}
    \caption{Factors that play a role in enhancing Joule heating for Simulation VIIIb over VIIIa. (a, b) Central up-north cuts of Pedersen current for Simulation VIIIa -- b. (c, d) Same for Hall current. (e, f) East-north plots of electric field's magnetic north component from Simulation VIIIa -- b. (g, h) Same for Hall conductance. (i, j) East-north plots of the third term in Equation \ref{e:fac} for Simulation VIIIa -- b. (k, l) Same for the first term in Equation \ref{e:fac}. Data sources: \citeA{simulationslink}.}
    \label{f:joule_viii}
\end{figure}

The bands of precipitation enhanced Hall conductance for Simulations VIIIa -- b peak at around 60 and 80 S respectively, as shown in panels g -- h. Now, since their spatial morphology comes from the same imagery, it implies that this increase in peak value also increases $\nabla_\perp\Sigma_H$, enhancing its associated FAC contributions as per Equation \ref{e:fac}. Panels i -- j show that these contributions, in this case, are in the opposing direction with respect to the total FAC driver---the third term in Equation \ref{e:fac} creates an upward current sheet where the driver map expects a downward sheet, and vice-versa. Given that the second term, $\mathbf{E}\cdot\nabla_\perp\Sigma_P$, can only help balance this by increasing $|\mathbf{E}|$, it would do so equally to that third term, $(\mathbf{E}\times\mathbf{b})\cdot\nabla_\perp\Sigma_H$. This leaves the local polarization to help balance the FAC, as is evident in panels k -- l of Figure \ref{f:joule_viii}. As before, all the input maps have nearly the same spatial morphology for both simulations. This is why, to increase $\nabla_\perp\cdot\mathbf{E}$, the simulation assuming an unaccelerated spectrum has a higher peak electric field, resulting in enhanced Joule heating despite the dissipationless Hall current enhancement.

Simulations VIIIa -- b both have intense total energy flux and strong electric field strengths. Both of these qualities aid in high-altitude FAC closure. Both simulations, however, also have large FAC requirements and higher electron energy distribution peaks---qualities which, all else held equal, suggest the need for lower-altitude-reaching FAC flux tubes. Given the self-consistent nature of our auroral arc simulations, we can only hypothesize that a more moderate FAC top-boundary driver would result in less need for Hall closure, but it could also be true that this would result in a Cowling mechanism with reduced electric field strength. In contrast, Comparisons VII and IX (see Supporting Information) both pertain to auroral arc systems whose FAC requirements, precipitation energy fluxes, and energy distribution peaks are relatively low. These combinations of parameters, even in the case of a weaker electric field in Comparison IX, results in simulations whose assumption of electron energy distributions matter less in both current closure and Joule heating as a result of FAC source term balancing.

\subsubsection{Summary: Precipitation Spectra}
We have shown that, if a particular auroral arc system requires Hall currents for FAC closure, choosing unaccelerated Maxwellian energy spectra for precipitating electrons is too restrictive when attempting to best represent the resulting impact ionization. Decoupling the energy spread from the most probable energy allows the modeling of auroral arcs whose electrons are accelerated from much colder source regions compared to their acceleration potential. Even for relatively ``hot'' accelerated precipitation, such as that from Comparisons VIII and IX ($T_s=800$ -- $860$ eV), the alternate, unaccelerated choice of $U_a=T_s=E_0$ still grossly overestimates the depth reached by the electron density enhancements. Holding FAC demands constant, this matters most when the average electric field strength is sufficiently weak, and/or the precipitation is low-reaching, i.e. any factor that puts emphasis on the Hall conductivity layer. Furthermore, unaccelerated Maxwellian electron distributions can overestimate the Hall currents as a whole, as well as the height-integrated Joule heating.

We have shown that specific assumptions of electron precipitation spectra can change the interpretation of auroral arc systems. Aptly, recent increases in the availability of multi-spectral, rather than white-light, all-sky imagery allows the community to move away from the assumption of unaccelerated Maxwellian precipitation spectra, and toward energy distributions which decouple the energy spread from the peak energy, allowing for more flexibility in modeling electron precipitation.

\subsection{Along-Arc FAC Structure}
Of our six conjunction events, two have a double-spacecraft arc crossing. This gives us an opportunity to look at two sensitivities: (1) how does along-arc structure in FAC affect current closure, and (2) how much confidence can be had in the replication technique we use. Our double replications have a weighting scale length of 50 km (roughly the distance between the orbits of Swarm A and C) when transitioning from replications of either track. This is described in more detail by \citeA[Section 2.3]{vanirsel2024}. When performing a weighted replication with plasma flow data, this can result in arbitrary along-arc gradients which affect the first term in Equation \ref{e:fac}. In our case, though the along-arc gradient in $j_\parallel$ resulting from this weighting scale length is arbitrary, these gradients do not directly enter into Equation \ref{e:fac}. I.e. the source of error these gradients introduce are the same as those introduced by the replication technique. Following are comparisons between double versus single replications of our two double-spacecraft conjunction events.

\subsubsection{Comparison X: Along-Arc FAC Structure}
Comparison X looks at our February 10 conjunction event with Simulation Xa using a weighted replication of both Swarm A (east) and C (west), and Simulation Xb which uses a replication of Swarm A data only. Their orbits cut roughly through the center of the simulation volume at about 47 km apart, which allows us to look at current closure on either side of the tracks. Looking at Figure \ref{f:compare-feb10-sat}, panels b and d, reveals an up-down FAC sheet pair that exists in Simulation Xa, but not in Xb. From the western boundary, centered around 40 km south-from-center, a roughly -2 to 1 $\mu$A/m$^2$ FAC pair follows the arc boundary up until just past the Swarm C FAC data track, from which this signature is replicated. Furthermore, the southernmost downward FAC sheet narrows and intensifies, when transitioning from the Swarm A to C tracks, from about 25 km wide and 1 -- 1.5 $\mu$A/m$^2$ in magnitude, to around 10 km and 2 $\mu$A/m$^2$. Simulation Xb has this FAC sheet remain unchanged along the arc.

\begin{figure}
    \centering
    \includegraphics[width=\textwidth]{15_joule_heating_20230210_35487_09_SD_AM_AC_v_SD_AM_xA.pdf}
    \caption{Comparison X (February 10, 9:51 UT): Top and side views of Simulation Xa with a FAC replication using both Swarm A (east) and C (west) (a, b) versus Simulation Xb with a FAC replication using only Swarm A (c, d) along with height-integrated Joule heating for Simulation Xa (e) and Xb (f). For plot details, see Section \ref{s:fluxtubes}. Data sources: \citeA{swarmlink}, \citeA{superdarnlink}, and \citeA{simulationslink}.}
    \label{f:compare-feb10-sat}
\end{figure}

With these differences in replicated FAC maps in mind, Simulation Xa (Figure \ref{f:compare-feb10-sat}a -- b), though being the same as Simulation Ia, here shows a different set of current flux tubes. They are calculated (in reverse) from ellipses placed at the southernmost upward, precipitating current sheet located east of, west of, and in between the two FAC data tracks. This helps illustrate the affect on current closure resulting from the difference between the two data tracks. Figure \ref{f:compare-feb10-sat}c -- d (Simulation Xb) shows flux tubes that are calculated from the same three ellipses.

The orange flux tube (0.5 kA) lies almost entirely east of the Swarm A track, hence it remains mostly unchanged, both in morphology and quantity. The green flux tube, however, is encroaching on the aforementioned Swarm C replicated FAC pair and thus captures around 0.1 kA more in Simulation Xb. The electric field across the arc has the flux tubes directed southwest to northeast, such that the green flux tube has its influx end entirely on the western side of Swarm C. Here, the Simulation Xa downward current sheet is stronger, but less than half the width compared to its Simulation Xb counterpart. The steeper FAC across-arc gradient in Simulation Xa pinches the downward green flux tube end into a teardrop shape (at roughly 25 km east and 50 km south), while its higher FAC density aids in capturing that additional 0.1 kA.

The red flux tube lies completely on the western side of Swarm C and captures the upward part of the FAC sheet pair introduced by Swarm C. At 0.7 kA, this gives it an additional 0.3 kA over the red flux tube in Simulation Xb. The adjacent downward current sheet helps close 0.1 kA of this added current, while the remainder is closed with a similar teardrop shaped flux tube end.

Comparison X outlines how a double versus a single FAC data track replication can introduce, albeit relatively minor, FAC signatures in the along-arc direction. We have to assume such signatures can appear and disappear over distances on the order of 50 km in every FAC replication. The major FAC structure, however, is conserved, suggesting the replication methodology holds.

\subsubsection{Comparison XI: Along-Arc FAC Structure}
Due to limitations of the all-sky imagery of the March 14 conjunction event, the simulation region for Comparison XI is almost completely west of both Swarm tracks. This prevents us from sourcing current flux tubes on either side of the data tracks, however we can still use Comparison XI to provide insight into what confidence can be had in the replication technique, and deliberate about the extent to which auroral arc FAC varies in the along-arc direction.

Panels a -- b in Figure \ref{f:compare-mar14-sat} show results from Simulation XIa, which is driven by a FAC map replicated from both Swarm A and C data. However, given the locations of the data tracks, most of this replication uses data from Swarm C, as it is the closest to the simulation region. With Simulation XIb (panels c -- d) using only Swarm A in its FAC replication, this is essentially a Swarm A versus Swarm C comparison.

\begin{figure}
    \centering
    \includegraphics[width=\textwidth]{16_joule_heating_20230314_24547_09_SD_AM_AC_v_SD_AM_xA.pdf}
    \caption{Comparison XI (March 14, 6:49 UT): Top and side views of Simulation XIa with a FAC replication using both Swarm A (east) and C (west) (a, b) versus Simulation XIb with a FAC replication using only Swarm A (c, d) along with height-integrated Joule heating for Simulation XIa (e) and XIb (f). For plot details, see Section \ref{s:fluxtubes}. Data sources: \citeA{swarmlink}, \citeA{superdarnlink}, and \citeA{simulationslink}.}
    \label{f:compare-mar14-sat}
\end{figure}

In contrast to Comparison X, here we see two FAC replications that, though varying somewhat, are structurally very similar. The southernmost return current sheets for both Simulations XIa -- b are similar in strength, width, and location, as is shown by the red flux tubes who capture around 1.2 kA in the same place for both simulations. The return current sheet just above, captured by the green flux tubes, is around half as strong in Simulation XIb and positioned $\sim$7 km southward, and the orange flux tubes carry a similar 1.5 -- 1.6 kA of Hall current at nearly the same location in both simulations. Overall, Comparison XI provides support for the extrapolation of FAC data over a distance of around 50 km.

\subsubsection{Summary: Along-Arc FAC Structure}
Two of our six conjunction events benefit from being able to use a second data track in their replications and subsequent simulations. Comparisons X and XI show to what extent the FAC map can change in just under 50 km, providing important insight into the confidence of all of our FAC replications, and consequently the resulting 3-D simulations of these auroral arc systems. Overall, contingent on the morphology indicated by the imagery and aside from minor FAC signatures, replicating the FAC data using arc boundaries defined by auroral imagery is a justifiable method for creating 2-D, continuous driver maps for 3-D simulations of auroral arc systems.

\section{Discussions \& Conclusions}\label{s:conclusions}\label{s:discussions}
Current closure morphology and Joule heating from resulting closure currents are two important aspects of ionospheric physics, particularly surrounding discrete auroral arc systems. By carefully incorporating observational data from multi-instrument conjunctions into input drivers of auroral arc simulations, we point out three aspects that the results are susceptible to: (1) the constant background flow, (2) the specifics of electron precipitation, and (3) the along-arc structure in FAC and the arc-boundary replication technique. Here, we conclude our findings and discuss possible future studies that can advance from this work.

Auroral arc systems should be studied in three dimensions to fully understand field-aligned current closure and, by extension, Magnetosphere-Ionosphere-Thermosphere coupling. We show, using several permutations of 3-D, electrostatic, data-driven, auroral arc simulations across six conjunction events, that flux tubes of electric current navigate around one another in their closure paths; something they cannot do in height-integrated (east, north), or cross-arc (north-up) two-dimensional descriptions. These current flux tubes tell the story of how FAC, ionospheric electric fields, and Pedersen and Hall conductivities interplay in a cohesive, self-consistent manner, and they do so with more detail than 2-D descriptions allow.

To produce top-boundary driver maps for our simulations, we demonstrate the use of auroral-imagery-guided FAC replication, similar to methods outlined by \citeA{clayton2019,vanirsel2024}. We show that this method can produce FAC maps that are geophysically consistent with maps of precipitation energetics, and that hold reasonably well for major arc-scale FAC structure. However, more minor FAC structure may appear or disappear when moving in the along-arc direction over distances on the order of 50 km. Along-arc variations below this scale size may well exist, but these are beyond the scope of this work given that we are limited by the 50 km Swarm spacecraft separation. Even so, this methodology uses maximal information from imagery derived precipitation maps to provide geophysically meaningful extrapolations of FAC surrounding auroral arcs. 

The 3-D auroral arc simulations covered in this paper have been shown to be very sensitive to both the magnitude and the direction of the constant, large-scale, background electric field, $\bar{\mathbf{E}}$. Equation \ref{e:bar_fac} shows what the choice of $\bar{\mathbf{E}}$ implies about the 2-D top-boundary FAC driver map, and thus, how the simulations interpret these maps. We draw the following conclusions about how $\bar{\mathbf{E}}$, in the absence of neutral winds, affects discrete auroral arc systems:
\begin{itemize}
    \item Strong background convection fields can render the use of Hall currents in FAC closure negligible, while weak background convection fields put emphasis on both local polarization fields and FAC closure through the electrojet.
    \item Across-arc electric fields provide shorter closure paths making FAC close through Pedersen current more often.
    \item FAC sheets close with adjacent ones only in the direction of the electric field.
    \item When part of the electric field is directed along the arc, it lengthens the closure paths and, as current flux tubes cannot intersect, it pushes additional tubes to Hall current altitudes.
    \item The manner in which the background electric field affects current connectivity, along with the electric field itself, significantly alters Joule heating, and thus the electrostatic load characteristics, of auroral arc systems.
\end{itemize}
These conclusions support the idea that large-scale convection flow conditions are a dominant driver of the specific morphology of auroral arc systems, with which the arc-scale ionosphere perturbs $\mathbf{E}$ in accordance with a 3-D conductivity volume.

In addition to background convection fields, auroral arc systems are also sensitive to the specifics of precipitating electron energy distributions. We show that the use of unaccelerated Maxwellian spectra can erroneously enhance impact ionization at lower altitudes, resulting in an overestimation of \textit{E}-region densities. We compare the use of such spectra against accelerated Maxwellian spectra, which decouple the peak energy from the energy spread. Even for simulations whose source region characteristic energy is estimated to be relatively large, $T_s\sim800$ eV, the unaccelerated assumption still greatly overestimates density enhancements at lower altitudes. We show the following:
\begin{itemize}
    \item The availability of multi-spectral, rather than white-light, all-sky imagery allows the community to move away from the assumption of unaccelerated Maxwellian precipitation spectra, and toward energy distributions which decouple the energy spread from the peak energy, allowing for more flexibility in modeling electron precipitation.
    \item The choice of accelerated versus unaccelerated Maxwellian electron precipitation most affects FAC closure in auroral arc systems when the Hall currents play a considerable role in this closure.
    \item Keeping FAC and total precipitating electron energy drivers constant, the choice of unaccelerated over accelerated precipitation alone can increase the calculated height-integrated Joule heating by 30 -- 50 \% in some auroral arc systems.
    \item Unaccelerated Maxwellian auroral precipitation assumptions can greatly enhance electrojet currents compared to accelerated precipitation assumptions.
\end{itemize}

This work looks at how to determine geophysical, self-consistent solutions to current continuity in auroral arc systems, and what these systems are sensitive to, thus uncovering how important various parameters can be. How then do we know which solution is correct? The existence of TII ion drift data (or other, independent data such as those from the Electrojet Zeeman Imaging Explorer mission \cite{madelaire2023}) from the Swarm spacecraft invites comparisons to the calculated GEMINI output flow maps covered in this paper. Figure \ref{f:tii_example} shows two such comparisons of the magnetic eastward TII flow (assuming no along-track component) across the model space for two of the simulations (Ib and IVb). While we have generated 17 simulations of the six events in Table \ref{t:events}, only the February 10 and March 14 conjunctions include Swarm A TII data; only the former has the crossing directly within the model space. The simulations using PFISR for the background flow for these two cases match better than the corresponding SuperDARN runs, which have smaller background flows.

\begin{figure}
    \centering
    \includegraphics[width=\textwidth]{17_tii_example.pdf}
    \caption{GEMINI versus TII flow comparisons for Simulations Ib (a) and IVb (b). The GEMINI magnetic eastward, ``v2'', and northward, ``v3'', plasma flows are interpolated through the simulation volume at the Swarm A tracks. TII magnetic eastward ion drift data, ``ViMagE'', are converted to geomagnetic coordinates assuming no along-track component. Both simulations use PFISR derived background flow, accelerated Maxwellian precipitation, and double-spacecraft replications. Data sources: \citeA{swarmlink}, \citeA{pfisrlink}, and \citeA{simulationslink}.}
    \label{f:tii_example}
\end{figure}

It is notable among the examples chosen for this study (the six events in Table \ref{t:events}) that there is not a particularly strong correlation between the magnetic and electric field signatures in the raw Swarm data---for most of these events, the $\nabla\cdot\mathbf{E}$ term in Equation \ref{e:fac} is apparently not the major player for the events in Table \ref{t:events} (except for perhaps Simulation IIIb). Thus this comparison with TII becomes mostly a question of matching the background flow to the TII value, perhaps why the nearer source (PFISR) provides the closest match. For the first example shown, there are some $\nabla\cdot\mathbf{E}$ signatures in both TII and the GEMINI results, but the GEMINI result is somewhat smoother and slightly offset. Both of these differences may well be artifacts of the image inversion process.

Finally we can consider whether the competition between the $\nabla\cdot\mathbf{E}$ and $\nabla\Sigma_{P,H}$ terms in Equation \ref{e:fac} provides a truly unique solution to the problem posed. There is a strong dependence on the chosen $\bar{\mathbf{E}}$: choosing the background electric field differently finds different situations. There may be choices, beyond what PFISR and SuperDARN provide, which more closely track the TII cross-track flow values. We do see that choosing different background flows, e.g. the no-background flow run versus the large-background flow simulations in Comparison III, generates in the GEMINI result a visible $\nabla\cdot\mathbf{E}$ signature which, though still there, is masked when the imposed background electric field is strengthened. Future work exploring these comparisons with TII should include (a) events like the one covered by \citeA{clayton2021}, with its strong $\nabla\cdot\mathbf{E}$ signature; and (b) further study of error sources stemming from matching the spacecraft data to inverted imagery, particularly for oblique camera angles which tend to blur and misplace discrete arc structures. We also note the scale of smoothing applied for these runs, as described in Subsections \ref{s:inversion} and \ref{s:swarm}: this level of smoothing may yet be hiding relevant physics, particularly at sharp arc edge boundaries.

The tools developed herein provide a means for data-driven event case study simulations to be routinely done, assuming sufficient data coverage. Upcoming iterations may consider different, incomplete combinations of input and/or adaptation of our methods into a formal physics-based assimilation scheme. A subject for further studies is the investigation of the sensitivity to other simulation choices made, such as data smoothing, conductance isocontours, etc. Also, the relevant physical gradient limit caused by recombination and collisions in the current closure altitude region: how sharp of gradients can be sustained and be relevant?

In the collective effort to try and understand the nature of aurorae, the instruments that provide our observational data are an ever-existing limitation. It would be optimal to deploy 1000s of spacecraft, radars, and imagers across the northern and southern auroral ovals \cite{nykyri2025}, but this is impractical. Hence, measurements must be targeted and focused on parameters that are most influential to the physics at hand. This work provides three such aspects to contribute to this focus and aids in making decisions as to what is important and when.

\appendix
\section{Derivation of Accelerated Bi-Maxwellian Differential Number Flux}\label{a:diffnumflux}
In order to implement the impact ionization calculations by \citeA{fang2010}, we need the differential (as a function of energy) hemispherical number flux, i.e. electrons/eV/s/cm$^2$, of precipitating energetic auroral electrons at the topside of the ionosphere for every latitude-longitude pair. To derive this flux for an accelerated population we start with a bi-Maxwellian source at the plasmasheet as is done by \citeA{fridman1980}:
\begin{equation}
g_s(v_{\parallel,s},v_{\perp,s},\varphi)\mathrm{d}^3v=n_{e,s}\left(\frac{m_e}{2\pi}\right)^{3/2}\frac{1}{E_{\parallel,s}^{1/2}E_{\perp,s}}\exp\left[-\frac{m_ev_{\parallel,s}^2}{2E_{\parallel,s}}-\frac{m_ev_{\perp,s}^2}{2E_{\perp,s}}\right]v_{\perp,s} \mathrm{d} v_\parallel \mathrm{d} v_\perp \mathrm{d}\varphi,
\end{equation}
where $n_{e,s}$ is the source region electron density, $m_e$ is the mass of an electron, $E_{\parallel,s}$ and $E_{\perp,s}$ are the parallel and perpendicular characteristic energies, $v_{\parallel,s}$ and $v_{\perp,s}$ are the source region parallel and perpendicular speeds, and $\varphi$ is the azimuthal coordinate. As electrons precipitate down towards the ionosphere they undergo no collisions---their velocities change in two ways only \cite{knight1973,fridman1980,kaeppler2013}:
\begin{enumerate}
    \item The conservation of the first adiabatic invariant, i.e. the mirror force, increases their perpendicular velocity:
    \begin{equation}
        v_{\perp,s}=\frac{1}{\sqrt{\beta}}v_{\perp,i},
    \end{equation}
    where $\beta=B_i/B_s>1$, and $B_i$ and $B_s$ are the ionospheric and source region magnetic field strengths.
    \item The conservation of energy increases the square magnitude speed as they fall through the parallel potential difference, $U_a$:
    \begin{equation}
    v_{\parallel,i}^2+v_{\perp,i}^2=v_{\parallel,s}^2+v_{\perp,s}^2+\frac{2U_a}{m_e}.
    \label{e:coe}
    \end{equation}
\end{enumerate}
This provides the parallel source region speed as a function of the ionospheric coordinates:
\begin{equation}
    v_{\parallel,s}=\pm\sqrt{v_{\parallel,i}^2+v_{\perp,i}^2\frac{\beta-1}{\beta}-\frac{2U_a}{m_e}}.
\end{equation}
From here, we use Liouville's theorem which tells us that, along a well-defined path through phase space, e.g. $(\mathbf{x},\mathbf{v})_s\to(\mathbf{x},\mathbf{v})_i$, the phase space density is held constant such that
\begin{equation}
    g_i(\mathbf{x}_i,\mathbf{v}_i)=g_s(\mathbf{x}_s,\mathbf{v}_s).
\end{equation}
A good assumption is to say that we may separate spatial and velocity coordinates, $g(\mathbf{x},\mathbf{v})=n(\mathbf{x})f(\mathbf{v})$, and that locally the densities are constants, i.e. $n_i(\mathbf{x})=n_{e,i}$, $n_s(\mathbf{x})=n_{e,s}$. This tells us
\begin{equation}
    g_i(\mathbf{v}_i)=g_s(\mathbf{v}_s)=g_s\left(\mathbf{v}_s(\mathbf{v}_i)\right),
\end{equation}
such that
\begin{equation}
    g_i(v_{\parallel,i},v_{\perp,i})\mathrm{d}^2v=n_{e,s}\frac{m_e^{3/2}/\sqrt{2\pi}}{E_{\parallel,s}^{1/2}E_{\perp,s}}\exp\left[-\frac{m_e\left(v_{\parallel,i}^2+v_{\perp,i}^2\frac{\beta-1}{\beta}-\frac{2U_a}{m_e}\right)}{2E_{\parallel,s}}-\frac{m_ev_{\perp,i}^2/\beta}{2E_{\perp,s}}\right]\frac{v_{\perp,i}}{\sqrt{\beta}} \mathrm{d} v_\parallel \mathrm{d} v_\perp
\end{equation}
where we've integrated over $\varphi$. The ionospheric density is thus
\begin{equation}
    n_{e,i}=n_{e,s}\frac{E_{\parallel,s}\sqrt{\beta}}{E_{\parallel,s}+E_{\perp,s}(\beta-1)}\exp\left[\frac{U_a}{E_{\parallel,s}}\right].
\end{equation}
Note that $U_a\to0$ and $E_{\parallel,s}\to E_{\perp,s}$ gives a familiar density relation: $n_{e,i}=n_{e,s}/\sqrt{\beta}$. Now that we have the velocity distribution function at the ionosphere, we find the differential number flux using $J_{\parallel,i}(\mathbf{v}_i)\mathrm{d}^3v=v_{\parallel,i}g_i(\mathbf{v}_i)\mathrm{d}^3v$ and then we perform the following change of coordinates:
\begin{equation}
    v_{\parallel,i}=v\cos\theta=\sqrt{2E/m_e}\cos\theta\ \mathrm{and}\ v_{\perp,i}=v\sin\theta=\sqrt{2E/m_e}\sin\theta,
\end{equation}
with $\theta$ being the pitch angle, and with Jacobian determinant $1/m_e$. The energy, $E$, has the condition
\begin{equation}
    E=\frac{m_e}{2}\left(v_{\parallel,i}^2+v_{\perp,i}^2\right)\ge U_a,
\end{equation}
as per Equation \ref{e:coe}. This gives
\begin{equation}
    J_{\parallel,i}(E,\theta)\mathrm{d} E\mathrm{d}\theta=\frac{n_{e,s}}{\sqrt{m_e}}\frac{1}{E_{\parallel,s}^{1/2}E_{\perp,s}}\frac{\sin2\theta}{\sqrt{2\pi\beta}} E\exp\left[-\frac{E-U_a}{E_{\parallel,s}}-\left(\frac{E}{E_{\perp,s}}-\frac{E}{E_{\parallel,s}}\right)\frac{\sin^2\theta}{\beta}\right]\mathrm{d}E\mathrm{d}\theta.
\end{equation}
With unit-less parameters $\varepsilon\equiv E/E_{\parallel,s}$, $U_a\equiv u_a/E_{\parallel,s}$, and $\delta\equiv E_{\perp,s}/E_{\parallel,s}$, we get
\begin{equation}
    \frac{1}{n_{e,s}}\sqrt{\frac{m_e}{E_{\parallel,s}}}J_{\parallel,i}(E,\theta)\mathrm{d}E\mathrm{d}\theta=\frac{\sin2\theta}{\sqrt{2\pi\beta}}\frac{\varepsilon}{\delta}\exp\left[-\left(\varepsilon-u_a\right)-\left(\frac{\varepsilon}{\delta}-\varepsilon\right)\frac{\sin^2\theta}{\beta}\right]\mathrm{d}\varepsilon\mathrm{d}\theta.
\end{equation}
We now integrate over $v_{\parallel,i}>0$, i.e. $0\le\theta\le\pi/2$, and find the hemispherical differential number flux.
\begin{equation}
J_{\parallel,i}(\varepsilon)\mathrm{d} \varepsilon=n_{e,s}\sqrt{\frac{E_{\parallel,s}}{m_e}}\frac{1}{\delta\sqrt{2\pi\beta}}G\left(\frac{\delta-1}{\delta\beta}\varepsilon\right)\varepsilon e^{-\varepsilon+u_a}\mathrm{d}\varepsilon,\ \mathrm{where}\ G(x)\equiv\frac{e^{x}-1}{x}.
\end{equation}
For similar parallel and perpendicular source temperatures, we have $\delta\sim1$, and we have $\beta\sim 10^3$ for a plasmasheet source region \cite{fridman1980}, where $G(x\ll1)\to1+x/2+\mathcal{O}(x^2)$ such that
\begin{equation}
    J_{\parallel,i}(\varepsilon)\mathrm{d} \varepsilon\approx n_{e,s}\sqrt{\frac{E_{\parallel,s}}{m_e}}\frac{1}{\delta\sqrt{2\pi\beta}}\left(1+\frac{\delta-1}{2\delta\beta}\varepsilon\right)\varepsilon e^{-\varepsilon+u_a}\mathrm{d}\varepsilon
\end{equation}
If we re-cast this in terms of normalized total precipitating energy flux, $q_p\equiv Q_p/E_{\parallel,s}$, where
\begin{equation}
    q_p=\int_{u_a}^\infty \varepsilon J_{\parallel,i}(\varepsilon)\mathrm{d}\varepsilon,
\end{equation}
we get
\begin{equation}
    J_{\parallel,i}(\varepsilon)\mathrm{d}\varepsilon=q_p\frac{1+\chi\varepsilon}{2+6\chi+u_a(2+u_a+(6+u_a(3+u_a))\chi)}\varepsilon e^{-\varepsilon+u_a}\mathrm{d}\varepsilon,\ \mathrm{where}\ \chi=\frac{\delta-1}{2\delta\beta}.
    \label{e:flux}
\end{equation}
We note that in our regime of $\beta\sim10^3$ we may ignore the temperature difference at the source, so if we take the limit of $\delta \to 1$ we get a familiar result
\begin{equation}
J_{||,i}(E)\mathrm{d}E=\frac{Q_p}{T_s^2+\left(T_s+U_a\right)^2}\frac{E}{T_s}\exp\left[-\frac{E-U_a}{T_s}\right]\mathrm{d}E,\ E\ge U_a
\end{equation}
where, for clarity, we have defined $T_s\equiv E_{\parallel,s}$. These results have been congregated from knowledge and derivations obtained in publications by \citeA{medicus1961,evans1974,fridman1980,strickland1989,kaeppler2013}.

\section*{Open Research Section}
All 3-D simulation data, imagery inversions, and supporting metadata are available at \url{https://rcweb.dartmouth.edu/LynchK}. The data for the Poker Flat DASC are available at \url{http://optics.gi.alaska.edu/optics/archive}, for AMISR at \url{https://data.amisr.com/database}, for SuperDARN at \url{https://superdarn.ca/data-download}, and for the Swarm at \url{https://swarm-diss.eo.esa.int}. The GEMINI source code and documentation is available at \url{https://github.com/gemini3d} and the replication/visualization tools at \url{https://github.com/317Lab/aurora_gemini}.

\section*{Conflict of Interest declaration}
The authors declare no conflicts of interest relevant to this study.

\acknowledgments
We thank NASA for funding Jules van Irsel from Grant 80NSSC23K1622 through the FINESST program and Dartmouth College for providing internal funding for Phase A of the Auroral Reconstruction CubeSwarm proposal from which a lot of this research stems, and we acknowledge NASA Grant 80GSFC21C0009 for the ARCS MIDEX CSR funding. AMISR facilities, including PFISR, are operated and maintained for the US National Science Foundation (NSF) by SRI International under a cooperative agreement \#AGS-1840962. AMISR data are produced and made publicly available for scientific and academic research purposes under the same award. This work uses version 0302, level 1B data from the Thermal Ion Imagers and version 0401, level 2 FAC data made possible by the European Space Agency's Swarm Data, Innovation, and Science Cluster and Canadian Space Agency grant 15 SUSWARM. We thank NASA for providing funding for the GEMINI model development from Grant 80NSSC24K0126. 3-D simulations supporting this work were facilitated by the Discovery cluster at Dartmouth Research Computing and the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center. The authors acknowledge the use of SuperDARN data and thank Daniel Billett for collecting SuperDARN convection flow data for use in this paper. SuperDARN is a collection of radars funded by national scientific funding agencies of Australia, Canada, China, France, Italy, Japan, Norway, South Africa, United Kingdom and the United States of America. We thank Michael Hirsch, Terry Kovacs, and John Griffin for GEMINI and computational technical support.

\bibliography{vanirsel_etal_2025}

\end{document}