Cast in the form of a linear programming problem, we hold
$\mathbb{S}_{z\mathcal{I},-\frac{1}{2}}$ fixed and seek the least-squares
solution for $\bm{\mathbf{L}}_{z\mathcal{I}}$

$\mathbb{S}_{z\mathcal{I}}\cdot\bm{\mathbf{L}}_{z\mathcal{I}}=\bm{\mathbf{f}}{},$
(51)

where all terms are evaluated at the $k=-\frac{1}{2}$ interface. The solution
minimizes the residual vector $\mathcal{R}$:

$\mathcal{R}=\min_{\bm{\mathbf{L}}_{z}}\Biggl{(}\mathbb{S}_{z\mathcal{I}}\cdot\bm{\mathbf{L}}_{z\mathcal{I}}-\bm{\mathbf{f}}{}\Biggr{)}.$
(52)

Applying SVD to solve this system of equations141414As implemented in the
code, we use the LAPACK dgelsd() routine to solve the SVD problem using the
divide and conquer algorithm. See documentation at https://netlib.org/.
determines the minimum-norm vector of incoming characteristic derivatives
$\bm{\mathbf{L}}^{*}_{z\mathcal{I}}$ that results in the solution

$\bm{\mathbf{f}}^{*}=\mathbb{S}_{z\mathcal{I}}\bm{\mathbf{L}}^{*}_{z\mathcal{I}},$
(53)

where the $*$ indicates the minimizing solution to the problem described by
(51); all terms are again evaluated at the $k=-\frac{1}{2}$ interface. We
emphasize that $\bm{\mathbf{f}}^{*}$ may not equal $\bm{\mathbf{f}}$ in
general, but will if an exact solution exists. Substituting the minimum-norm
solution $\bm{\mathbf{f}}^{*}$ for the term
$\mathbb{A}_{z\mathcal{I},ij,-\frac{1}{2}}(\partial_{z}\bm{\mathbf{U}})_{ij,-\frac{1}{2}}$
into the finite-volume form of the time derivative given at Equation (45), we
find

$\partial_{t}\bm{\mathbf{U}}_{ij0}=-\bm{\mathbf{f}}^{*}-\sum_{n,l\pm}\mathbb{A}_{n,l\pm}(\partial_{n}\bm{\mathbf{U}})_{l\mp}-\bm{\mathbf{D}}_{ij0},\quad\begin{subarray}{c}n\in(x,y,z)\\\
l\pm\in(i\pm\frac{1}{2},j\pm\frac{1}{2},k+\frac{1}{2})\end{subarray}.$ (54)

The interpretation of the SVD solution is that it is the physically admissible
time derivative on the boundary, $\partial_{t}\bm{\mathbf{U}}_{0}$, that most
closely matches the target time derivative $\overline{\dot{\bm{\mathbf{U}}}}$.
To see this, formally decompose the minimum norm vector as

$\bm{\mathbf{f}}^{*}=-\sum_{n,l\pm}\mathbb{A}_{n,l\pm}(\partial_{n}\bm{\mathbf{U}})_{l\mp}-\bm{\mathbf{D}}_{ij0}-\dot{\bm{\mathbf{U}}}^{*}.$
(55)

The time derivative is then

$\displaystyle\partial_{t}\bm{\mathbf{U}}_{ij0}$
$\displaystyle=-\bm{\mathbf{f}}^{*}-\sum_{n,l\pm}\mathbb{A}_{n,l\pm}(\partial_{n}\bm{\mathbf{U}})_{l\mp}-\bm{\mathbf{D}}_{ij0}$
(56)
$\displaystyle=\Bigl{(}+\sum_{n,l\pm}\mathbb{A}_{n,l\pm}(\partial_{n}\bm{\mathbf{U}})_{l\mp}+\bm{\mathbf{D}}_{ij0}+\dot{\bm{\mathbf{U}}}^{*}\Bigr{)}-\sum_{n,l\pm}\mathbb{A}_{n,l\pm}(\partial_{n}\bm{\mathbf{U}})_{l\mp}-\bm{\mathbf{D}}_{ij0}$
(57)
$\displaystyle=\dot{\bm{\mathbf{U}}}^{*}=\bm{\mathbf{f}}{}-\bm{\mathbf{f}}^{*}+\overline{\dot{\bm{\mathbf{U}}}}.$
(58)

In essence, the SVD method produces a physically-admissible correction factor
$(\bm{\mathbf{f}}{}-\bm{\mathbf{f}}^{*})$ to the requested target time
derivative $\overline{\dot{\bm{\mathbf{U}}}}$. The correction factor is
nothing but the residual vector,
$\mathcal{R}=\bm{\mathbf{f}}{}-\bm{\mathbf{f}}^{*}$. If
$\mathcal{R}=\bm{\mathbf{0}}$ then $\bm{\mathbf{f}}^{*}=\bm{\mathbf{f}}{}$ and
the target update can be exactly matched using a linear combination of
incoming characteristic derivatives; i.e., the target lies in the image of the
matrix $\mathbb{S}_{z\mathcal{I}}$. On the other hand, any non-zero part of
the residual vector lies in the null-space of the matrix
$\mathbb{S}_{z\mathcal{I}}$, which is the span of the outgoing characteristic
modes that comprise the matrix $\mathbb{S}_{z\mathcal{O}}$.

Unlike the simple inversion method described following Equation (48), the SVD
method is a true optimization method. In particular, different weights can be
applied to each of the primitive variables to normalize their respective
orders of magnitude, or additionally, to reflect the level of certainty with
which the target value of each variable is known. A typical situation in
observational solar physics is that the magnetic fields and velocities are
inferred with more certainty than the density or temperature. SDO/HMI, for
instance, is a highly optimized instrument for determining the vector magnetic
field and Doppler velocities in certain regimes, has continuum data that can
be used to infer the temperature, but gives little information on the density
or pressure. This could be handled by providing an estimate for $\rho$ and
$\epsilon$ as constant in time, but weighting those variables much less than
the target magnetic fields. The SVD solution then matches the better-known
variables and finds a self-consistent solution for the less well known fields,
$\rho$ and $\epsilon$ in the current example. In this way, the full
information content of the incoming characteristics is harnessed to find a
physically valid solution to the data driving problem. We discuss the case of
simply holding $\rho$ and $\epsilon$ uniform in space and constant in time on
the boundary (as was done in, e.g., the data-driven simulations compared in
Toriumi et al. (2020)), in §7.2.

Uncertainties in the variables are handled by applying a weighting
transformation to the SVD problem, which simply involves left multiplying
Equation (51) by a diagonal weighting matrix $\mathbb{W}$. This transforms
$\mathbb{S}_{z\mathcal{I}}\mapsto\widetilde{\mathbb{S}_{z\mathcal{I}}}=\mathbb{W}\mathbb{S}_{z\mathcal{I}}$
and
$\bm{\mathbf{f}}{}\mapsto\widetilde{\bm{\mathbf{f}}{}}=\mathbb{W}\bm{\mathbf{f}}{}$.
SVD then produces the minimum-norm characteristic derivative vector
$\bm{\mathbf{L}}_{z\mathcal{I}}^{*}$. It can still be used directly in the
numerical integration defined above because the form of (51) is the same for
both the weighted and unweighted cases; the minimum-norm characteristic
derivative vector for the weighted system is, of course, different than that
for the unweighted system. In §7 we present results for the case with no
weighting matrix, then demonstrate use of a weighting matrix in §7.2. In both
cases, $\bm{\mathbf{f}}^{*}$ is calculated by left-multiplying
$\bm{\mathbf{L}}_{z\mathcal{I}}^{*}$ by the _unweighted_ incoming subspace
matrix $\mathbb{S}_{z\mathcal{I}}$, as written at Equation (53).

Once the minimum norm solution $\bm{\mathbf{L}}_{z\mathcal{I}}^{*}$ is known,
the derivative of primitive variables between $k=0$ and $k=-1$,
$(\partial_{z}\bm{\mathbf{U}}{})_{-\frac{1}{2}}$, can be solved for by
multiplying Equation (49) from the left by
$\mathbb{A}_{z\mathcal{I},-\frac{1}{2}}^{-1}$. Looking at Table 1, the
presence of any incoming characteristic modes involves at least 7 primitive
variables, while the eighth ($B_{z}$) can always be defined using the
solenoidal condition for magnetic fields. The derivative of the state vector
is therefore always well defined for any number of incoming characteristic
derivatives, or is unnecessary if there are no incoming characteristic
derivatives. After the spatial derivative is known the state vector at $k=-1$,
$\bm{\mathbf{U}}_{-1}$, can be set accordingly. In practice, time updating the
$k=0$ level in the _Driving_ simulation only requires the product
$\mathbb{A}_{z\mathcal{I}}\partial_{z}\bm{\mathbf{U}}=\mathbb{S}_{z\mathcal{I}}\bm{\mathbf{L}}_{z\mathcal{I}}^{*}$,
so the values of primitive variables at $k=-1$ are not needed. They are still
used, however, as part of the coupling between the _Driving_ and the _Data-
Driven_ simulation.

The best fit solution for the time derivative, (55), is then time-integrated
by any appropriate method to produce the time-updated estimate to the driving
layer. Again keeping things simple, we apply a first order forward-Euler
scheme, for which the solution is

$\bm{\mathbf{U}}(t_{O1+\eta}+\delta
t,\bm{\mathbf{x}}_{b})=\bm{\mathbf{U}}(t_{O1+\eta},\bm{\mathbf{x}}_{b})-\delta
t\biggl{[}\mathbb{S}_{z\mathcal{I}}\bm{\mathbf{L}}_{z\mathcal{I}}^{*}+\sum_{n,l\pm}\mathbb{A}_{n,l\pm}(\partial_{n}\bm{\mathbf{U}})_{l\mp}+\bm{\mathbf{D}}_{ij0}\biggr{]},\quad\begin{subarray}{c}n\in(x,y,z)\\\
l\pm\in(i\pm\frac{1}{2},j\pm\frac{1}{2},k+\frac{1}{2})\end{subarray}.$ (59)

This equation defines the time-advance of the driving layer in the _Driving_
simulation CHAR. Stability is not an issue because new information is
continually injected into the _Driving_ simulation at each time step. Recall
that this solution is still only approximately correct because we have used
the directionally split subspace defined at $k=0$ in place of the subspace at
$k=-\frac{1}{2}$. As discussed following Equation (48), a nonlinear
optimization method would determine a simultaneous solution for both the
incoming subspace and the spatial derivative across the interface state, i.e.,
by finding
$\mathbb{A}_{\mathcal{I}}\partial_{z}\bm{\mathbf{U}}=\mathbb{S}_{\mathcal{I}}\bm{\mathbf{L}}_{z\mathcal{I}}$
evaluated at $k=-\frac{1}{2}$. This would generate a fully self-consistent
solution, to which the current linear solution is a first approximation. While
we have explored using a nonlinear optimization method, thus far the solutions
are not fully stable, and we therefore retain the linear optimization for now.

Above the boundary layer (in cells $k\geq 1$), the regular numerical update to
the _Data-Driven_ simulation is carried out for time step $\eta$ using
whatever numerical scheme is implemented in that simulation. At the end of the
$\eta$ time step, the time-updated solution for the driving layer ($k=0$) and
the ghost layer ($k=-1$) are copied from the _Driving_ simulation to the
_Data-Driven_ simulation and the time-updated solution from the interior of
the _Data-Driven_ simulation is copied into the _Driving_ simulation (cells
$k=1,2$). In this way the two simulations are kept coupled together and the
_Driving_ simulation is ready to calculate the next step. A new target time
derivative is defined for the driving layer using Equation (35), and the
process repeats until the next observed state
$\widehat{\bm{\mathbf{U}}}(t_{O2},\bm{\mathbf{x}}_{b})$ is reached and new
observational data needs to be ingested.

Summarizing all the steps in the algorithm, using the SVD based optimization
method:

  1. 1.

At the beginning of a timestep $\eta$ (which updates the simulation from
$t_{O1+\eta}$ to $t_{O1{}+\eta+1}$), the MHD state is known at all cell center
locations $k\geq 0$, and is consistent between the _Driving_ simulation and
the _Data-Driven_ simulation.

  2. 2.

An observationally derived target time derivative for the primitive variables
in the $k=0$ driving layer is calculated using Equation (35). Equivalently, a
target boundary state may be specified via (36).

  3. 3.

The interface states in the driving layer are calculated using (37) at the
five faces of each cell where all information is known (everywhere but
$k=-\frac{1}{2}$).

  4. 4.

The eigensystem at each interface location is calculated through Equation (4)
and using the left and right eigenmatrices given in Appendices A, B, and C.
All quantities are calculated using the interface values of the primitive
variables.

  5. 5.

The eigensystem of each interface state is split into the grid-positive and
grid-negative subspaces using the projection operators (20) (also calculated
using the interface states), via (24).

  6. 6.

An optimization method is used to find the minimum-norm solution of (47) for
incoming characteristic derivatives through the lower face,
$\bm{\mathbf{L}}_{z\mathcal{I}}^{*}$ at $k=-\frac{1}{2}$. This solution for
the characteristic derivative vector implies an exterior state at $k=-1$ that
minimizes the difference between the target time derivative and the numerical
time derivative, i.e., the state $\bm{\mathbf{U}}_{-1}(t_{O1+\eta})$, while
also producing the physically admissible time updated state at $k=0$, i.e.,
$\bm{\mathbf{U}}_{0}(t_{O1{}+\eta+1})$; these two states are stored in
temporary variables.

  7. 7.

The primitive variables in the _Driving_ simulation at level $k=-1$ are
changed to the values just calculated ($\bm{\mathbf{U}}_{-1}(t_{O1+\eta})$),
as this is the state of the exterior plasma that will produce the desired time
update.

  8. 8.

The _Data-Driven_ simulation performs its standard numerical time integration
of the MHD equations from time $t_{O1+\eta}$ to $t_{O1{}+\eta+1}$ for cells
with $k\geq 1$. When needed, the values of primitive variables at levels $k=0$
and $k=-1$ (or spatial derivatives defined in terms of them) are pulled from
the _Driving_ simulation to be used by the _Data-Driven_ simulation during
this step (these values may be spatially interpolated by the _Data-Driven_
simulation to other grid locations, as needed).

  9. 9.

The time-updated values in the $k=0$ layer calculated in Step 6,
$\bm{\mathbf{U}}_{0}(t_{O1{}+\eta+1})$, are applied in the _Driving_
simulation and copied to the _Data-Driven_ simulation.

  10. 10.

The time-updated values in layers $k=1,2$ are copied from the _Data-Driven_
simulation to the _Driving_ simulation.

  11. 11.

The full time integration is now complete for both the _Driving_ simulation
and the _Data-Driven_ simulation, and their values are consistent with each
other for $k\geq 0$ (the $k=-1$ state can only be calculated in Step 6).

  12. 12.

Repeat Steps 1-11 until the simulation reaches time $t_{O2}$.

The major simplification in the current implementation described above comes
in Step 6, where we find the linear SVD solution to the minimization problem
as given at (53). For this, we approximated the incoming subspace at the cell
interface $k=-\frac{1}{2}$ by the incoming subspace at $k=0$; that is, we set
$\mathbb{A}_{z\mathcal{I},-\frac{1}{2}}=\mathbb{A}_{z\mathcal{I},0}$ (or the
equivalent in the eigendecomposed system). The assumption is that the cell
centered state and the interface MHD states are similar enough, and the
resulting eigendecompositions are similar enough, that this will produce an
approximately correct answer, especially since the solution is constantly
“nudged” in the right direction from timestep to timestep. However, this does
mean that the external state $\bm{\mathbf{U}}_{-1}$, interface state
$\bm{\mathbf{U}}_{-\frac{1}{2}}$, and boundary state $\bm{\mathbf{U}}_{0}$,
and the associated matrices and characteristic derivatives projected onto the
incoming subspace at the interface ($\mathbb{S}_{z\mathcal{I},-\frac{1}{2}}$,
$\mathbb{A}_{z\mathcal{I},-\frac{1}{2}}$, and
$\bm{\mathbf{L}}_{z\mathcal{I},-\frac{1}{2}}$), are not exactly consistent. A
fully nonlinear optimization approach would self-consistently determine
$\mathbb{S}_{z\mathcal{I},-\frac{1}{2}}$ and
$\bm{\mathbf{L}}_{z\mathcal{I},-\frac{1}{2}}$, and hence
$\bm{\mathbf{U}}_{-1}$ via
$\mathbb{S}_{z\mathcal{I},-\frac{1}{2}}\bm{\mathbf{L}}_{z\mathcal{I},-\frac{1}{2}}=\mathbb{A}_{z\mathcal{I},-\frac{1}{2}}(\partial_{z}\bm{\mathbf{U}})_{-\frac{1}{2}}$,
but we have not yet fully implemented and tested the nonlinear approaches.

Finally, we note that any optimization approach to determining the optimal
data-driven boundary condition can be extended to include other constraints.
For instance, when applied to solar data and including gravity, it may make
sense to preserve, say, the mean pressure on the boundary, the total mass flux
through the boundary, or to match the energy or helicity flux through the
boundary, as estimated by some other method. This extendable feature makes the
optimization method a robust choice for constructing physically motivated
data-driven boundary conditions.

##  6 Test implementation using CHAR-DDBC to data drive LaRe3D

Figure 5: Perspective comparison of the _GT_ (left), $\Delta t{}=20\delta t{}$
_Data-Driven_ (middle) and $\Delta t=100\delta t{}$ _Data-Driven_ (right)
simulations at four times throughout the simulation (top to bottom: 0, 125,
250, and 375$\delta t$, or times $t=0,\ 0.38,\ 0.89,\ 1.42$). The _Data-
Driven_ models are driven with data extracted from the _GT_ simulation at
$z=1.09$. The vertical cut shows vertical velocities in the $y=0$ plane, and
the horizontal cut shows vertical magnetic field in the $z=1.14$ plane (the
first non-driven cell in the _Data-Driven_ simulations). Seed points of field
lines are initiated both below the driving layer (_GT_ simulation only) and
above it (all simulations; see text in Appendix F following (F12)). An
animation of this Figure is included in the online material. The animation
shows the expansion of the spheromak through the driving layer for the _GT_
simulation and the reproduction of those field lines in the two _Data-Driven_
simulations. The low-cadence, $\Delta t=100\delta t$ simulation has notably
deformed field lines compared to the other two simulations.

To demonstrate the generality of this approach we have developed a
characteristic-based data-driven boundary condition (DDBC) add-on module
called CHAR. This module is a separate, finite-volume MHD code where the time
update is calculated in terms of the MHD characteristics. Importantly, the
CHAR module is not specific to any particular meshing, gridding, or numerical
time-advance. The interfaces to the CHAR module were written so that it can be
coupled to, and drive, any MHD code. Furthermore, CHAR can be used as a kind
of numerical glue to couple different plasma codes to an MHD code, e.g., CHAR
could couple a kinetic plasma code to an MHD code. Analogously, CHAR could
even be used to drive an MHD simulation in a layer in the interior of an MHD
domain.

As proof of concept, the characteristic based CHAR module was coupled to the
LaRe3D MHD code which utilizes a Lagrangian remap method to advance the MHD
equations (Arber et al., 2001). We emphasize that the numerical algorithm and
gridding used to advance the MHD equations is different in CHAR and LaRe3D.
CHAR runs concurrent with LaRe3D and overlaps the two LaRe3D ghost cells and
first two interior cells along the LaRe3D lower boundary, that is, CHAR is the
_Driving_ simulation and LaRe3D is the _Data-Driven_ simulation introduced at
the beginning of §5. The vertical layout of the grids are shown in Figure 3,
in which CHAR is the red simulation and LaRe3D is the green simulation.

All the primitive variables in the characteristic code are stored at cell
centers and represent cell averaged quantities, while the LaRe3D variables are
represented on a staggered grid. The CHAR simulation is updated in-step with
the normal LaRe3D simulation, i.e., at the typical LaRe3D Courant step. The
full bookkeeping is described in the algorithm summary at the end of §5, but
the cartoon version is that at the end of each numerical time step, the CHAR
variables that are located inside the LaRe3D simulation boundaries are
overwritten using values interpolated from the LaRe3D simulation values, while
the LaRe3D ghost-cell values are overwritten using interpolated CHAR values.
In this way, the two simulations are kept consistent with each other while
still using the characteristic decomposition of MHD to define the boundary
condition. The reason for taking this hybrid approach is that LaRe3D is
fundamentally a Lagrangian code, not a characteristic-based code, so the
characteristic decomposition is awkward to formulate purely in terms of the
LaRe3D numerics151515The awkwardness arises for a variety of reasons that all
stem from fundamental differences in the numerical schemes. LaRe3D treats
variables on a staggered grid with a predictor-corrector update that
incorporates mass and energy conservation by associating specific derivatives
with specific control volumes. Calculating all necessary quantities at the
intermediate steps in LaRe3D in terms of the characteristics would appear to
require evaluating the characteristic derivatives projected onto MHD states at
different locations and times and then combining them in a numerically stable
way. It is unclear how to do that. The far more simple case of a incorporating
the characteristics into a finite difference scheme failed for similar
reasons, as described in Appendix D.. Therefore, our tests below and those
described in Paper II for nonreflecting boundary conditions demonstrate that
our method can be used to robustly couple together two MHD simulations with
disparate numerical schemes or varieties of approximations.

The data driving method developed in §5 was validated using a three-step
process analogous to that described in Leake et al. (2017). First, a uniformly
meshed ground truth (_GT_) simulation was developed and run in LaRe3D. The
initial condition was an initially unbalanced spheromak that expands into a
field-free homogeneous plasma. The numerical domain in the _GT_ simulation is
larger than the total expansion of the spheromak in the portion we will
eventually use for data driving. The left hand columns of Figure 5 and Figure
6 show side-on and top-down perspectives of the _GT_ simulation, respectively,
in a cut-out domain that is horizontally the same size as the _Data-Driven_
simulations, to ease comparison.

Second, we extracted the eight primitive variables from a constant layer
$z=z_{b}$ in the _GT_ simulation that lies just outside and above the initial
spheromak. This layer is denoted by the dark plane in the left column of
Figure 5. The values of variables on the staggered LaRe3D grid are
appropriately averaged to cell-centered locations and used for the extraction.
The gray scale images in left columns of Figure 5 and Figure 6 show the
vertical magnetic field $(B_{z})$ in the extraction plane. The unbalanced
spheromak expands through this plane and, as the time-sequence in each figure
and online animations show, the evolution of the magnetic field at and above
the plane resembles many features of an emerging bipolar active region,
including separation of the opposite magnetic polarities and rotation of a
sigmoidally shaped polarity inversion line at the center of emergence.

Finally, we set up a suite of driven simulations in a reduced numerical domain
that extends from the $z_{b}$ plane selected in Step 2 to the top edge of the
_GT_ simulation. Data from the extracted plane are used to test the data
driving method. That is, time series of primitive variables from the extracted
plane are the synthetic observations that will be ingested by the driving
layer at $k=0$, as described in Section §5. This setup allows for a variety of
tests. In the following investigation we perform two types of validation
tests. The first focuses on varying the input cadence of the “observations,”
similar in spirit to the tests performed by Leake et al. (2017), except in our
data driving scheme the boundary is constrained to evolve consistently with
MHD and the interior of the simulation. The middle and right columns of Figure
5 and Figure 6 compare two examples that data drive between progressively
further spaced snapshots extracted from the _GT_ simulation, i.e., we vary the
driving cadence $\Delta t$. In a second type of test we do not provide our
algorithm any information on the time evolution of either the density or
internal energy density in the driving plane. These results are discussed in
detail in §7. Appendix F provides the full analytic description of the
spheromak and the numerical details of each simulation.

Figure 6: Same as Figure 5 but from a top-down perspective, showing _GT_
(left), $\Delta t{}=20\delta t{}$ _Data-Driven_ (middle) and $\Delta
t=100\delta t{}$ _Data-Driven_ (right) simulations at four times throughout
the spheromak expansion (top to bottom at Courant steps 0, 125, 250, and 375).
The vertical magnetic field is shown in grayscale in the $z=1.14$ plane. Seed
points of the field lines are initiated at the $z=1.25$ plane (see text in
Appendix F following (F12)). An animation of this Figure is included in the
online material. The animation shows the expansion of the field lines through
the driving layer. The vertical magnetic field in the driving layer starts at
zero. Subsequently the bipolar magnetic field appears as an expanding yin-yang
pattern that rotates as the spheromak enters the driven volume. Both the field
lines and the vertical magnetic field in the driving layer match very well
between the _GT_ and higher-cadence $\Delta t=20\delta t$ _Data-Driven_
simulations, but the lower cadence simulation, with $\Delta t=100\delta t$, is
severely distorted.

In slightly more detail, the initial condition in the _GT_ simulation produces
a blast wave with an embedded, expanding spheromak. As described in Appendix
F, the spheromak includes a shielding current sheet that magnetically isolates
it from its homogeneous, field-free surroundings; when we speak of the
“spheromak” we mean that closed magnetic volume of plasma. The blast wave
extends beyond the edge of the (magnetic) spheromak and hence the outer edge
of the blast wave is a purely hydrodynamic shock that propagates to and
through the edges of the _GT_ domain (absent numerical boundaries it would
propagate to infinity). The expansion of the blast wave is initially
nonlinear, becomes increasingly linear after $t\approx 1.5$, and is overall a
bit more rapid and extensive in the spheromak’s rotationally symmetric $y$ and
$z$ directions compared to the $x$ direction. By $t=3.0$, when we stopped the
_GT_ simulation, the blast wave had reached a distance of about $3.3$ in the
$y$ and $z$ directions and $3.2$ in the $x$ direction. At long times the blast
wave would be roughly spherically symmetric, and represents mechanically
radiated energy lost from the system.

The hydromagnetic spheromak undergoes a much more contained expansion, as
defined by the expansion of the magnetic field. Similar to the blast wave, the
spheromak’s initial expansion is nonlinear and becomes roughly linear around
$t=1.0$, by which time it has expanded to a radius of $1.5$ in the $x$
direction and $1.7$ in the $y$ and $z$ directions. By $t=3.0$ it has expanded
to $1.74$ in the $x$ direction and $2.0$ in the $y$ and $z$ directions. We
stop our analysis before any disturbance reaches the periodic boundaries into
the _Data-Driven_ domain, above $z_{b}$ (which, again, has reduced size in $x$
and $y$ compared to the _GT_ domain). Our simulations therefore fall into
category (iii) of typical boundary conditions for simulated plasmas from the
Introduction §1. Paper II explores how the structure of a passively advected
spheromak are modified by interaction with several standard nonreflecting
boundary conditions.

In order to demonstrate the capabilities of our method and explore its
limitations, we have performed two different types of data driving tests. The
first batch of tests is a series of driven simulations using progressively
lower cadence subsets $\Delta t$ of the full cadence boundary data extracted
from the _GT_ simulation at every Courant-step $\delta t$. Figure 2 provides a
schematic overview of the process. At each Courant step in the _Data-Driven_
simulation, a target boundary state is defined and boundary data is temporally
integrated from the current state to the state that is optimally close to the
target state while still fulfilling the constraints imposed by MHD, as
described in §5. At one extreme, the extracted boundary data is provided at
every Courant step from the _GT_ simulation to the _Driving_ simulation so
that $\Delta t$=$\delta t$. This best case scenario is only limited by the
accuracy of the coupling between the two different numerical schemes (i.e.,
the characteristic code CHAR and Lagrangian-remap code LaRe3D). At the other
extreme, the extracted boundary state is only provided at one out of every
hundred of the ground truth simulation time steps. This simulation uses
driving data that severely undersamples the dynamic evolution of the expanding
spheromak (see §7.1). In total, driving was tested using $\Delta
t=1,5,20,50,\text{ and }100$ timesteps $\delta t$ pulled from the _GT_
simulation. Taken together, this suite of simulations demonstrates the
transition from an accurate best case reproduction of the _GT_ by the _Data-
Driven_ simulation to a poor reproduction of the _GT_ once a dynamic sampling
threshold is crossed. This result echoes the findings of Leake et al. (2017).
However in contrast to Leake et al. (2017) the boundary of each simulation in
this suite remains consistent with MHD: our characteristic method does not
allow setting the primitive variables at the boundary in an under- or over-
determined way.

In the second batch of tests we apply our SVD-based DDBC approach to the case
where information on the thermodynamic variables is missing. Here the “target”
states for both $\rho$ and $\epsilon$ are set to their respective initial
values and small relative weights are given to those variables when setting up
the SVD optimization problem. This allows the CHAR-DDBC approach to find a
better approximation to those variables. The results are compared to the case
in which the known solution for all primitive variables are provided to the
SVD minimization and to the results of more typical boundary conditions found
throughout the literature. These latter BCs use pure linear-interpolation of
the velocity and magnetic fields between observed boundary states, and hold
the density and internal energy density uniform and constant at their initial
values.

##  7 Results

Figure 5 shows a perspective overview of the ground truth simulation in the
left column and two versions of _Data-Driven_ simulations with $\Delta t$=20
and 100$\delta t$ as the cadence between boundary observations in the middle
and right columns, respectively. Figure 6 displays the same simulations and
times as Figure 5 but from a top-down perspective. Time increases downward
along each column, and animations of both figures are available in the online
material. The vertical cut in the $y=0$ plane of Figure 5 shows the vertical
velocity $v_{z}$ in a red-blue color scale (the velocity cut-plane is omitted
in Figure 6). The horizontal cut through the $z=1.14$ plane in each figure
shows vertical magnetic field $B_{z}$ in grayscale. This plane lies just
outside of the initial spheromak, as seen in the upper left frame of Figure 5,
and is the first non-driven cell in the driven domain, i.e., one cell above
the driving layer at $z_{b}=1.109$. The tubes in Figures 5 and 6 trace out
magnetic field lines, whose seed points are described in Appendix F just
before (F13). The locations of the seed points are constant in time and the
field lines are identified by random colors (which are also constant in time
for a given seed point).

The first column in both figures show that the combined spheromak and blast
wave from the _GT_ simulation expands through the driving layer and produces a
magnetic configuration similar to a sheared arcade above the driving surface.
This process is similar to how flux emerges on the Sun where the moment of
emergence is preceded by a sheath current that forms the boundary between the
emerging active region and ambient corona. The flux then expands into the
ambient corona. This is especially easy to see towards the end of the
simulation in the bottom row of Figure 6, where the field lines originating
further from the polarity inversion line (where $B_{z}=0$) have different
orientations when they cross the inversion line compared to field lines
originating closer to the inversion line. Comparison between the left and
middle columns of each figure shows that running a _Data-Driven_ simulation
with a driving cadence of $\Delta t=20\delta t$ produces fields and flows in
the driven simulation that are very close to the _GT_ simulation (though some
slight differences are apparent). However, when the driving frequency
decreases to just one out of every hundred snapshots ($\Delta t=100\delta t$),
as shown in the right hand column of both figures, the differences relative to
the ground truth simulation are more pronounced.

The top-down perspective in Figure 6 makes it very easy to see that the
difference between the _GT_ and $20\delta t{}$ _Data-Driven_ case (left and
middle panels) are quite subtle while for the $100\delta t{}$ driven case
(right panel) the difference, even qualitatively, is drastic. In the latter
case, the distribution of magnetic flux near the driven boundary developed a
remarkably different pattern compared to the _GT_ case as a result of the
data-driven algorithm; correspondingly distorted field lines were produced
higher up in the _Data Driven_ simulation volume. The large disparity between
the true and driven solutions, even near the driving layer, is caused by the
combination of a temporal sampling of the driving layer that under-resolves
the dynamic evolution of the system and our algorithm’s prioritization of the
physics of the MHD equations over temporally-interpolated boundary evolution.

Put simply, the discrepancy with the _GT_ simulation at the largest driving
cadence of $\Delta t=100\delta t{}$ demonstrates that the linear interpolation
between two driving states on the lower boundary is inconsistent with the MHD
equations for that cadence. Instead, when a boundary is so sparsely observed
in time that a linear interpolation between two states is not a valid solution
to the MHD equations, our data driving method will instead produce an optimal
solution that _is_ consistent with the MHD equations. We discuss the issue of
a maximum allowed driving cadence further in §7.1. If one wishes to data drive
a simulation under those conditions, either additional independent information
or a different set of assumptions must be provided to give a better estimate
of the boundary evolution.

Note that any mismatch between the requested and produced boundary states
serves as a check on the validity of a _Data-Driven_ simulation regardless of
the presence of ground truth data. This check between requested versus
produced boundary states in the driving layer is an inherent feature of our
data driving method and, critically, is present in the scenario where the data
are derived from actual observations.

Figure 7: Comparison of all primitive variables in the interior of the _Data-
Driven_ simulation at a horizontal slice at $k=15$ cells above the driven
boundary $(z\approx 1.59)$, between the _GT_ simulation (columns 1, 5) and
simulations driven every 20$\delta t$ (columns 2,6), 50$\delta t$ (columns
3,7), and 100$\delta t$ (columns 4,8). The horizontal axis is $x$, the
vertical axis is $y$, and data corresponds to Courant step $320$, or $t=1.18$,
. The color scale for each primitive variable is common for all simulations.
The top row compares the density in the left four columns and internal energy
density in the right four columns. The next three rows compare components of
the velocity field (left) in a Red-Blue color map and magnetic field (right)
in a Grayscale color map.

Figure 7 compares the values of primitive variables directly between four sets
of simulations. Each panel is a horizontal slice extracted from 15 cells above
the driving layer ($z=1.59$) at a single time midway through the spheromak
expansion ($t=1.1822$ or 320 Courant steps). The panels are labeled by
primitive variable and simulation, and similar variables are further
distinguished by different color scales: the python standard _viridis_ for
thermodynamic variables, red-blue for velocities, and grayscale for magnetic
fields.

The first row shows slices through density ($\rho$, left 4 columns) and
internal energy density ($\epsilon$, right 4), and the next three rows show
components of the velocity field ($v_{x},v_{y},v_{z}$, left columns) and
magnetic field ($B_{x},B_{y},B_{z}$, right). Columns 1 and 5 are for the _GT_
simulation, 2 and 6 for the $20\delta t{}$ driving case, 3 and 7 for $50\delta
t{}$ driving, and 4 and 8 for the $100\delta t{}$ driving simulation. The
color bar for a given primitive variable is shared across all simulations.
Here we see that all of the primitive variables are similar between the _GT_
and _Data-Driven_ simulations up to the $50\delta t{}$ driving cadence case,
with decreasing fidelity with increased cadence. In contrast, we once again
see that the $\Delta t=100\delta t{}$ simulation differs substantially from
the _GT_ case in all primitive variables: it has a qualitatively different
appearance.

Figure 8: Weighted mean-squared-difference (wMSD) as a function of time using
different driving cadences, with $\Delta t=1,5,20,50,$ and $100\delta t$ for
the blue, green, brown, violet, and yellow lines, respectively. (a) wMSD
calculated in the driving boundary layer, i.e., cell $k=100\ (k=0)$ in the
_GT_ (_Data-Driven_) simulation. (b) Zoom-in near the middle of the
simulation. (c) Same as (a) but 15 cells higher into the simulation (i.e., in
the plane shown in Figure 7). (d) Zoom-in shortly after the spheromak edge
reaches the $k=15$ layer in the _Data-Driven_ simulation.

To demonstrate the effect of varying the driving cadence in a more
quantitative form we calculated the weighted mean-squared difference between
each of the _Data-Driven_ simulations and the _GT_ simulation. For a driven
state $\bm{\mathbf{U}}{}$ and ground truth state $\bm{\mathbf{G}}{}$ we define
the weighted mean squared difference as

$\textsf{wMSD}(\bm{\mathbf{U}},\bm{\mathbf{G}};t)=\bigl{\langle}(\mathbf{U}-\mathbf{G})\mathbb{K}^{-1}(\mathbf{U}-\mathbf{G})^{T}\bigr{\rangle},$
(60)

where $\bm{\mathbf{U}}$ and $\bm{\mathbf{G}}$ are both functions of time. The
weighting matrix $\mathbb{K}^{-1}$ (not a function of time) is defined as the
inverse of the covariance matrix $\mathbb{K}$ between each pair of primitive
variables in the _GT_ data calculated at a single reference time:

$\mathbb{K}_{\sigma\zeta}=\bigl{\langle}(G_{\sigma}-\langle
G_{\sigma}\rangle)(G_{\zeta}-\langle G_{\zeta}\rangle)\bigr{\rangle}.$ (61)

$\mathbb{K}$ is a measure of the width of each variable’s inherent
distribution and the correlations between different variables. It is therefore
a suitable normalization coefficient between the variables. We calculate
$\mathbb{K}$ at the reference time $t_{\eta}=280$ Courant steps, midway
through the spheromak expansion, and use that single covariance matrix in the
calculation of wMSD at all times.

In the above equations, the expectation value of a quantity $q$ is defined

$\langle q\rangle=\frac{1}{N_{q}}\sum_{i\in\text{domain}}q_{i}$ (62)

where the domain can be any subset of the simulation, and has $N_{q}$ cells.
In the following, we calculate expectation values separately in each
constant-$z$ plane, so that $N_{q}=N_{x}\times N_{y}$.

The weighting matrix $\mathbb{K}^{-1}$ serves as a unit of measure of the
difference between variables in the _GT_ and _Data-Driven_ simulations.
Consider the most simple scenario in which all variables are uncorrelated and
normally distributed. In that case the covariance matrix $\mathbb{K}$ is
diagonal and each element is the variance of the associated variable (i.e.,
the square of its standard deviation). Now consider the argument of the wMSD
function applied to a single cell, in several situations: if one of the
variables in the _Data-Driven_ simulation differs from its _GT_ value by one
standard deviation of that variable’s inherent distribution, and all other
variables exactly match their _GT_ values, then applied to that single cell
$\textsf{wMSD}{}=1$; if every variable differs from its _GT_ value by one
standard deviation of its inherent distribution, then $\textsf{wMSD}{}=8$; and
if each variable is off by $1/\sqrt{8}$ of the standard deviation of its
distribution in the _GT_ simulation then the wMSD will be 1. This discussion
makes clear that, when considering the full 8-variable MHD state vector, a
wMSD value of 1 represents very little difference between the _GT_ and _Data-
Driven_ simulations. While correlations between the variables will make this
simple interpretation of wMSD more complex, we have found that the metric
provides a decent measure of the relative similarity between two simulations.

Figure 8 shows the result of calculating wMSD at each simulation time and in
two different simulation layers: the driving layer $k=0$ in Panel (a) and
$k=15$ cells into the _Data-Driven_ simulation in Panel (c). The second layer
corresponds to roughly half the maximum expansion of the spheromak in the
vertical direction. The expectation value in Equation (60) is calculated over
all $N_{x}\times N_{y}$ cells in each layer161616We have also experimented
with calculating the expectation value over the subset of cells that have
either $\lvert B\rvert>10^{-10}$ or
$\lvert\rho-\rho_{0}\rvert/\rho_{0}>10^{-5}$ in either the _Data-Driven_ or
_GT_ simulation. The intent was to better penalize either over- or under-
expansion of the sphereomak relative to the _GT_ simulation. Consistently
defining the “outer edge” of the sphereomak proved problematic. The end result
is that, while calculating expectation values over all cells in a given layer
reduces the relative value of wMSD for layers further from the boundary, it
produces an easier to understand metric of comparison across many
simulations.. The inset Panels (b) and (d) show zoomed-in regions of each
parent plot. Blue, green, brown, purple, and yellow lines correspond to
driving cadences of $\Delta t=$ 1, 5, 20, 50, and 100 Courant steps $\delta
t$, respectively.

The driving cadence for each of the _Data-Driven_ simulations is clearly
visible in Figure 8(a) as the periodic rise and fall of each mean difference
line. The minima of each curve coincide simply because the longer driving
cadence simulations ingest data at integer multiplies of most of the shorter
cadence simulations. An exception to this is seen by comparing wMSD curves for
the $\Delta t=20\delta t$ and $\Delta t=50\delta t$. As the inset Panel (b)
shows, the former (brown curve) has a difference maxima at $t=150$, in the
middle of its driving cadence, while the latter has a minima, right at the
moment when it ingests fresh data.

The linear interpolation of the state vector between observation times means
that, each time new data is provided to the driving algorithm, that
information is immediately incorporated into every location across the lower
boundary (albeit with a low amplitude). For longer driving cadences, and in
the context of this expanding spheromak simulation, this results in
instantaneous transverse phase velocities across the driving layer,
encompassing the full spatial extent of what will become newly emerged field.
This is, of course, unphysical behavior, because the true evolution in this
case is a gradual expansion outward from the central region. Instead, the
instantaneously provided data then propagates away from the driven boundary
everywhere it was introduced, increasing the wMSD metric in the _Data-Driven_
simulations. The differences only decrease once the true, gradual expansion in
the GT simulation catches up with the prematurely injected data in the _Data-
Driven_ simulations, at which point the metric begins to fall. In all cadence
regimes, our optimization based SVD solution to determine boundary conditions
in terms of the characteristic derivatives does find a _physically plausible_
evolution of the boundary. For shorter driving cadences, it reproduces the
true, GT solution; for longer cadences, the solution can differ substantially.

These three phases of the driving process are labeled in Panel (a) of Figure
8—as “Ingest”, “Propagate”, “Catchup”—for the extreme case of driving at
$\Delta t=100\delta t$, where some new, and rather inconsistent, data is
introduced at the “Ingest” time. Note that, through the definition of
$\overline{\dot{\bm{\mathbf{U}}}}$, this includes information at $t+\Delta t$.
The knowledge that new information was just introduced everywhere at the lower
boundary slowly diffuses through this system during “Propgate” phase, and the
errors increase during this time as the optimization method struggles to find
a physically possible way in which such information could have been
introduced. Once the midpoint of the driving cadence is reached, the
information has propagated through enough of the system that the algorithm is
able to find a physically allowed dynamical pathway towards the next
observation point; that is, the driven system begins to “Catch-up” to the
ground truth system. At the end of the cycle, a new batch of data is
introduced, and the process repeats.

Simulations with the most frequently ingested data naturally have the lowest
error. Note that driving at the ground truth’s Courant step, $\Delta t=1\delta
t$ (blue line), does not result in $\textsf{wMSD}=0$ because LaRe3D and CHAR
use two different numerical schemes. Coupling the two codes together
essentially introduces some slight numerical diffusion due to the
interpolation of primitive variables to different grid locations. As a result,
the wMSD for the $\Delta t=1$ case represents a best-case scenario for data
driving, with as faithful a reproduction as is possible to achieve in the
present case, and has $\textsf{wMSD}\lesssim 0.01$ throughout the simulation.
If we restrict the area over which wMSD is calculated to only regions into
which the spheromak has expanded, then the maximum difference for the $\Delta
t=1$ case is $\textsf{wMSD}{}\approx 0.33$. Very roughly, this means that, on
average, the MHD state vector in a random cell in the _Data-Driven_ simulation
will differ from the GT values in that same cell by approximately
$\sqrt{0.33/8}\approx 0.2$ of a standard deviation of each variable’s inherent
distribution. Of course, some variables may adhere much more closely to their
_GT_ counterparts than others.

From the inset panel Figure 8(b) it is clear that both the $\Delta t=$ 5 and
20$\delta t$ simulations are nearly able to match the quality of the $\Delta
t=1$ case in the boundary layer at each of their respective injection times
(green and brown curves). That is to say, in between the driving times at
which the solution in the boundary is known, the _Data-Driven_ simulation
becomes progressively worse up to the midpoint of the driving cadence and then
progressively better. Increasing the time between driving data leads to
correspondingly greater excursions away from the true evolution of the driving
boundary, but for the cases $\Delta t=1,5,20\delta t{}$, our data driving
algorithm is essentially able to match the observed boundary states. For the
$\Delta t=50\delta t{}$ case, however, our MHD-constrained data driving
approach is no longer able to match the driving boundary data, as seen when
the minima of the purple curve at Courant step $150$ is nearly the same as the
maxima of the brown curve: the best 50$\delta t$ time step has the same error
as the worst 20$\delta t$ time step. The $\Delta t=100$ case (yellow line) is
correspondingly worse.

Similar properties hold in the interiors of the _Data-Driven_ simulations, as
seen in Figure 8(c) and (d). These panels show the wMSD calculated in a plane
15 cells into the _Data-Driven_ simulation ($z\approx 1.6$), about halfway
between the driving layer and the greatest extent of the expanding spheromak
in the $z$ direction. The slower driving cadences for $\Delta t=50$ and
$\Delta t=100$, and the nonphysical linearly interpolated target states
between the starting and ending states, cause signals to reach this part of
the volume earlier than they should, i.e., the purple and yellow difference
curves deviate from zero prior to the blue, green, and brown curves. The
simulations with faster driving cadences have complex, quasi-periodic
deviations from the _GT_ case as new information propagates through the
system, as seen in the inset Panel (d), especially for the $\Delta t=20$ case.
The overall difference to the _GT_ simulation remains low in the spheromak
interior in that case, however. In contrast, simulations with slower driving
cadences periodically develop much worse reproductions of the _GT_ simulation
as new driving data is ingested and (improperly) linearly interpolated in
time.

###  7.1 Condition to determine the maximum allowed driving time

Figure 9: Time-distance plot showing the horizontal expansion of the spheromak
along the $x$-direction along the line $(y=0,z_{b}=1.11)$ for four primitive
variables: (a) $\rho$, (b) $B_{z}$, (c) $v_{x}$, and (d) $v_{y}$. The
approximate average expansion rate of $dx/dt\approx
2/1.4\left[L_{N}/T_{N}\right]$ is indicated by the white dashed line.

The above results indicate that there is a critical sampling cadence $\Delta
t_{\text{crit}}$ beyond which data driving will become increasingly incorrect.
This maximum time can be of thought of as an analogue of a Courant condition
for a discretely sampled physical system: an object that moves by more than
its own physical length scale between two observations cannot be unambiguously
identified.

Leake et al. (2017) estimated the critical sampling cadence in their 2D flux
emergence simulation. Following their example, we define the actual
observation cadence $\tau$, apparent horizontal feature velocity $v_{h}$, and
feature extent $L$ (the half-width-at-half-max of the feature). An unambiguous
association of a single spatial feature observed at two times requires
$\tau<\Delta t_{\text{crit}}\equiv L/v_{h}$. Note that this physical sampling
condition is independent of the simulation’s numerical Courant condition,
which is defined with reference to the numerical grid via $\delta
t=dx/max(|v|,c_{f})$.

Figure 9 shows the application of the critical sampling cadence to our
simulation. The four panels are time-distance plots extracted from the driving
plane along the $x$-axis on the line $(y,z)=(0,z_{b}=1.109375)$ for four
primitive variables—(a) $\rho(x,t)$, (b) $B_{z}(x,t)$, (c) $v_{x}(x,t)$, and
(d) $v_{y}(x,t)$—and show the expansion of the spheromak through the driving
plane. We define the physical extent $L$ as the half-width-half-max of one of
the lobes in $B_{z}$ at time $t=0.5$, which turns out to be $L\approx 0.2$. A
by-eye estimate for the horizontal velocity is indicated by the dashed white
line, which roughly corresponds to either the fastest expansion of the
vertical magnetic field (Panels b or d) or to the average nonlinear expansion
of the fast-mode shock traveling into the surrounding field-free plasma
(Panels a or c). These velocities are approximately equal and have the value
$v_{h}=dx/dt=2/1.4\approx 1.4$. This gives a critical sampling rate for this
feature of $\Delta t_{\text{crit}}=L/v_{h}\approx 0.2/1.4=0.14$. In order to
compare with the tests of various sampling rates earlier in this section, we
note that the average Courant step across the simulations is $\langle\delta
t\rangle=0.004$. So, we should see the fidelity of the reconstruction diminish
when driving at a cadence lower than $\Delta t_{\text{crit}}/\langle\delta
t\rangle\approx 36$ (after using the actual numerical values instead of the
one decimal values above). This rough estimate is indeed consistent with the
analysis of how wMSD varies with driving cadence, reported in Figure 8.

Clearly, the estimate of a critical sampling time will depend on the type of
feature used to perform the estimate. In general, it will be given by the
smallest, fastest moving feature that is present in a system that can be
observed for a given instrumental setup. As a point of comparison, Leake et
al. (2017) found that most of their 2D flux emergence simulations produced
features with horizontal lengths and velocities that would be accurately
sampled by the $12\,\mathrm{min}$ cadence of HMI vector magnetograms, but not
always. Their simulations with the fastest emergence violated this condition,
and as they discussed in some detail, their data driving algorithm, which
linearly interpolated the MHD state vector between input states, produced poor
results in those cases.

###  7.2 Comparison to typical data-driven approaches

By far the most common variety of data-driven boundary conditions for solar
research currently found in the literature are those that linearly interpolate
in time a set of observed magnetic field variables while holding the density
and energy fixed to their initial (often uniform) values (Bourdin et al.,
2013; Hayashi et al., 2018; Warnecke & Peter, 2019; Jiang & Toriumi, 2020;
Kaneko et al., 2021; Jiang et al., 2021; Guo et al., 2021; Inoue et al.,
2023). This is because the most consistently available synoptic data set for
solar observations over the last 13 years are provided by the SDO/HMI
instrument, which is optimized for magnetic field and Doppler velocity
measurements at moderate spatial resolution ($\approx 1^{\prime\prime}$) and
temporal cadence (typically 12 minutes and 45 seconds for vector and line-of-
sight magnetograms, respectively); other synoptic programs, such as SDO’s
predecessor SOHO/MDI or the ground-based NSO/GONG network, are similar but
have lower temporal- and spatial-resolution. The plasma density and
temperature (and thus pressure or internal energy density) can be inferred by
matching observed continuum intensity to a model solar atmosphere, however,
one typically desires multiple spectral regions and lines from different
elements to reduce uncertainties and biases (Maltby et al., 1986). Even the
most sophisticated inversions that combine multiple spectral regions, high
spectral resolution, sensitive polarimetry, and advanced model atmospheres can
contain order of magnitude errors in some quantities and regions (Borrero et
al., 2019; Borrero & Pastor Yabar, 2023). Beyond those issues, such data are
not available from the typical synoptic data sets, and we are not aware of any
data driving approaches that currently attempt to incorporate those derived
quantities.

Of the data driving approaches cited above, some use the entire magnetic field
vector while others use only the line-of-sight (or, approximately, radial)
component, and this choice varies from group to group and even paper to paper
within a single research group, depending on the scientific goals. In addition
to the magnetic field and Doppler velocity, transverse velocities can be
derived by applying feature tracking algorithms to magnetograms or other data
series (Welsch et al., 2004; Schuck, 2008; Fisher et al., 2015), but this
extra step is not used in all of the above references. Instead, different
authors again employ many different velocity boundary conditions, some
linearly interpolating the results of feature tracking algorithms, some using
symmetric boundary conditions, and some simply setting all the velocities to
zero on the boundary, even though this choice is physically inconsistent with
the induction equation. The through-line of this discussion is that the simple
approximation of holding $\rho$ and $T$ (or its equivalent) fixed while
linearly interpolating all or part of the magnetic and velocity fields remains
a popular data driving approach, even for groups that have previously used a
version of characteristic boundary conditions (e.g., compare Jiang & Toriumi
(2020) or Jiang et al. (2021) to Jiang et al. (2016)).

We therefore perform a new test to compare our CHAR-DDBC algorithm to these
more typical approaches. In the following, we refer to simulations that use
these typical (interpolation-only) boundary conditions as the “baseline”
simulations; note that the baseline simulations do not use the characteristics
at all. We develop an apples-to-apples comparison between our characterstic
approach and the baseline simulations by considering the case where no
information about either the density or internal energy density is supplied to
any driving simulation, i.e., these variables are estimated as being constant
and uniform at their initial values. All other primitive variables (velocity
and magnetic field) are assumed to be perfectly known at the driving cadence.
As such, this new test still represents a greater degree of certainty in the
MHD state vector at each location than one typically has from synoptic
observations. Specifically, the baseline data-driven simulations use the
typical approach of linearly interpolating $\bm{\mathbf{B}}$ and
$\bm{\mathbf{v}}$ in the $k=0$ layer between injection times, using a zero-
gradient assumption for the $k=-\frac{1}{2}$ layer (that is, setting
$\bm{\mathbf{U}}_{k=-1}=\bm{\mathbf{U}}_{k=0}$), and then modifying the
magnetic field normal to the driving plane ($\widehat{\bm{\mathbf{z}}}$) as
needed to enforce $\nabla\cdot\bm{\mathbf{B}}=0$. The mass density and
internal energy density are held constant and uniform at their initial values
in both the $k=0$ and $k=-1$ cells. As described, this procedure roughly
mimics typical approaches found in the literature cited above and related
works by many of the same authors.

Figure 10: Comparison of wMSD for our data driving algorithm versus linear
“baseline” interpolation algorithms. For each line type, the red (blue) line
shows driving at $\Delta t=20(50)\delta t$. The solid curves reproduce the SVD
based CHAR-DDBC results from Figure 8, for reference, the dashed curves show
the difference metric for SVD CHAR-DDBC when no information is known about
either the density $\rho$ or internal energy density $\epsilon$, and the
dotted curves show the metric when $\bm{\mathbf{B}}$ and $\bm{\mathbf{v}}$ are
linearly interpolated while $\rho$ and $\epsilon$ are held fixed to their
initial values. Panel (a) gives wMSD calculated in the driving boundary layer,
i.e., the $k=100\ (k=0)$ cell in the _GT_ (_Data-Driven_) simulation. (b) Same
as (a) but 15 cells higher into the simulation, i.e., the layer from Figure 7.

We then compare those baseline simulations to our data driving approach by
running another set of the SVD-based CHAR-DDBC simulations using the same
boundary data as provided to the baseline simulations (i.e., with constant and
uniform $\rho$ and $\epsilon$ on the boundary). This new set of CHAR-DDBC
simulations incorporates the weighting matrix for the SVD problem, as
described following Equation (58). The diagonal weighting matrix $\mathbb{W}$
has values for $\rho$ and $\epsilon$ of $1\times 10^{-6}$, while the weights
for magnetic and velocity variables have values of $1.0$. After applying the
weighting transformation the CHAR-DDBC algorithm runs as described in §5.

The results of both sets of simulations are summarized in Figure 10, which
compares the difference metric wMSD for the baseline data driving methods that
use only linear interpolation (dotted lines, labeled
$\text{{INT}}[\bm{\mathbf{U}}\setminus\\{\rho,\epsilon\\}]$) and our approach
based on the characteristics (dashed lines, labeled
$\text{{SVD}}[\bm{\mathbf{U}}\setminus\\{\rho,\epsilon\\}]$); the notation
$\setminus\\{\rho,\epsilon\\}$ emphasizes that information about the evolution
of density and internal energy density are not provided to either the baseline
or our characteristic boundary conditions. For reference, the solid lines
reproduce curves from Figure 8, i.e., they are the results of providing the
full MHD state vector to the CHAR-DDBC algorithm at every driving time step.
Shown in Panel (a) are results for the driving layer and in Panel (b) for the
interior of the driven region, i.e., $k=15$ cells above the driving layer, as
in Figure 10. Red (blue) curves correspond to driving at a cadence of $\Delta
t=20(50)\delta t$.

Comparing the dashed and dotted curves clearly shows that our CHAR-DDBC
algorithm performs better than the baseline approach to data driving, given
identical input. By the end of the simulation, the mean difference to the _GT_
simulation is reduced by more than a factor of two for our approach over the
baseline case. Comparison to the solid curves shows the effect of withholding
any information about the evolution of the thermodynamic variables. The large
jump in wMSD from the case where the correct evolution of mass density and
internal energy density are given to the case where they are not (solid curves
to broken curves) is due, in large part, to a reduced expansion of the
spheromak and a decrease in amplitude of the blast wave relative to the _GT_
simulation. These systematic errors are greater when the mass and internal
energy density are held fixed along the lower boundary (dotted lines), hence
the corresponding larger values in wMSD of the baseline simulations compared
to our characteristic-based (dashed). This demonstrates the ability of our
approach, based on the characteristic form of MHD, to fill in missing
information about the evolution of state vector at the driving boundary (upper
panel), leading to improved results in the simulation interior (lower panel).

In the current implementation of our approach, the weighted SVD algorithm
tends to find solutions for which the mass and internal energy densities
evolve in tandem, whereas the true solution has a low plasma density and high
internal energy density central region (see Figure 7, top row). This behavior
likely stems from the fact that SVD produces a minimum norm solution.
Additionally, as currently set up, the SVD problem attempts to find a solution
that will drive the density and internal energy density back to their initial
values. The result is that the weighted version of CHAR-DDBC produces low
density and low energy in the spheromak core, and so the wMSD metric is
dominated by an inaccurate (but still physically realizable!) evolution of the
internal energy density. This is a primary reason why the wMSD curves in Panel
(a) are only improved by a factor of 2 over the base baseline approach that
uses linear interpolation for $\bm{\mathbf{B}}$ and $\bm{\mathbf{v}}$.
Improvements to our method are expected to more closely match the case where
the full MHD state vector is known on the boundary at the driving times (solid
lines), and these will be explored in future work.

As mentioned above, there exist in the literature many variations on (non-
characteristic based) data driven boundaries that linearly interpolate some
subset of the primitive variables. In addition to the baseline simulation
described above, we have tested several additional versions of these driven
boundary conditions. We omit the full results to save space, but do summarize
a few important points here that help provide broader context for the
comparison of these other versions of baseline (interpolation-only) driven
boundaries to our CHAR-DDBC method. First, we explored several variations of
holding components of the velocity fixed at zero and linearly interpolating
the other variables. In all cases this generates similar results to the
baseline linear interpolation case presented in Figure 10 (dotted) but with
correspondingly larger values of the difference metric. Second, we ran a set
of simulations that linearly interpolate all of the primitive variables in
time using different driving cadences. This approach produces wMSDs that are
close to those based on the characteristics, i.e., as presented in Figure 8,
and outperform it in some instances. This is to be expected: a simple linear
interpolation is how the target update is defined in CHAR-DDBC to begin with
at (35), so at high cadence the two approaches should give nearly identical
results.

This latter test demonstrates that a linear interpolation is sufficient to
produce high fidelity results _provided the dynamic evolution of the system is
adequately sampled and accurate for every primitive variable._ However, the
utility of a simple linear interpolation decreases when information about the
system’s evolution is missing, either by low cadence observations or by
variables with large observational bias and uncertainty, in which case the
observed boundary evolution is unlikely to be consistent with the evolution of
the simulation interior. We would expect that a similar challenge would arise
when there is inconsistency between well observed variables, for example, when
the spectral lines that are used to determine temperature and magnetic field
are sensitive to different heights in the solar atmosphere. In that case, even
though the evolution of each variable may be well determined, they may not
relate to each other in a physically meaningful way, at least when used to
directly define a boundary condition. Indeed, Carlsson & Stein (1994) already
encountered this issue in the relatively simple case of a piston driver at the
lower boundary of a 1D radiative hydrodynamic simulation. The piston velocity
was to be given by the observed Doppler velocity of an Fe I line, but the
mismatch between the formation height of the line and their actual lower
boundary caused significant errors in the resulting synthesized Fe I line from
the simulation. These arose, essentially, because they had not properly
encoded the physics of wave propagation through the solar atmosphere. While
they were able to mitigate the problem in a clever way in that case, it well
illustrates the challenges posed by even slight mismatches between the
observations and simulations. Such inherent inconsistencies are expected
whenever primitive variables are inferred from observations.

The above tests show that the baseline approach of holding the density and
energy to fixed values at the boundary while linearly interpolating other
variables over time can significantly alter the expansion of the spheromak
into the driven simulation. This is likely because the spheromak expansion is
primarily a fast-mode type disturbance that requires gradients in both the
density and energy. Instead, when density and energy are held fixed, the first
phase of expansion away from the boundary has to be driven purely by the
magnetic field. If we analyze the baseline simulation’s behavior in terms of
the characteristics then the linear interpolation boundary driving is most
similar to one of the Alfvén modes (however, we again stress that there is no
guarantee that the behavior of the driven cells, when simply set using
linearly interpolated primitive variables, represents any allowed motion under
MHD). Just above the driven boundary—one cell into the driven simulation—the
Alfvénic disturbance at the boundary must couple to a combination of fast- and
slow-modes to drive further expansion away from the boundary. This results in
less expansion, weaker field strengths, and thus a smaller modeled spheromak
away from the boundary compared to the _GT_ simulation. In contrast, our
approach based on the characteristics allows the mass and internal energy
density to self-consistently evolve with the other variables and therefore
mitigates these issues, but (so far) does not fully eliminate them. We hope to
address remaining issues in future work.

Finally, we note that the magnetic topology of all the driven simulations
appear similar to the _GT_ simulation. However, we caution that the spheromak
test case has a rather simple topology to begin with. Our conclusion from this
final suite of simulations is that using simple linear interpolations of
primitive variables when holding some number of those variables fixed can
produce inaccurate results. The output may appear reasonable in the simulated
domain for some of the variables, but it is unclear how one could tell that
the recovered solution is far from the truth, especially if the linear driving
at the boundary is inconsistent with the actual equations being solved in the
rest of the simulation. In contrast, our use of the characteristics to ensure
an MHD-allowed evolution of the boundary, combined with an optimization
approach to find a solution that matches the observations as closely as
possible, produces a substantial improvement over the baseline approach.
Additionally and critically, for the characteristics, one can tell from the
difference between the requested and produced boundary evolution when the
requested evolution can no longer be accommodated by the MHD equations.
Finally, the optimization method can be extended to incorporate additional
constraints, such as estimated relative helicity, mass flux, or Poynting flux,
that could help to provide a better solution when limited, uncertain, or
inconsistent observations are available.

##  8 Discussion

We have presented a self contained description of the theory, numerical
implementation, and testing of a solution to the rather nuanced problem of
data-driven boundary conditions for magnetohydrodynamic simulations. The
previous sections already contain a good deal of discussion, so here we
summarize only at a high level and try to weave together some thoughts that
are only possible after all of the foregoing exposition.

The main problem we are trying to solve is that, in order to study the
dynamics of the solar atmosphere, we need an accurate representation of the
corona, especially (_but not solely!_) the magnetic field. Direct inference of
the complete 3D coronal state from observations, for instance using polarized
spectral lines that form under coronal conditions, is at least difficult if
not practically impossible on a routine basis with current technology. This is
despite very real and ongoing improvements in that field in terms of both
instrumentation (Tomczyk et al., 2008, 2016; Rimmele et al., 2020) and
reconstruction techniques based on observations (Kramar et al., 2016; Casini
et al., 2017; Dima & Schad, 2020; Li et al., 2021; Schad et al., 2023). On the
other hand, photospheric observations, even with low spectral resolution,
allow routine inference of much of the MHD state vector at the solar surface,
for instance using GONG (Harvey et al., 1996) or SDO/HMI (Pesnell et al.,
2012) for magnetic fields and Doppler velocities, optical flow techniques for
transverse velocities (Schuck, 2008). Inferring thermodynamic quantities
typically requires simultaneous observations of multiple spectral lines whose
ratios are sensitive to density and temperature. When such observations are
not available, as is the case for synoptic solar observations, estimates based
on empirically- and model-derived relations can be used (see Jaeggli et al.,
2012, and references in the Introduction of that work). The time evolution of
the photospheric state vector then provides a dynamic constraint to the plasma
above the photosphere. That fact raises the question of how best to
incorporate the observational data into the boundary of an MHD simulation.

The main goals of this paper were developed to answer that question: first, to
develop general boundary conditions that require all the primitive variables
on the boundary to be set in a way that is consistent with the MHD equations,
and second, to apply that formulation of boundary conditions to the data
driving problem. This is far from a trivial task because the dynamics of all
the variables are coupled together in a limited number of physically
admissible ways, a constraint that is very much at the heart of our approach
to data driving.

Setting one primitive variable to some value makes implicit assumptions about
how the system is evolving and thus implies specific values for a subset of
the other variables. As a simple example, suppose a boundary cell is updated
in time by $B_{x}\rightarrow B_{x}+\delta B_{x}$, which is meant to represent
a propagating Alfvén wave. In order to be self-consistent, this requires that
the velocity field should also updated by $v_{x}\rightarrow v_{x}-\delta
B_{x}/\sqrt{\rho}$ (for propagation in the same direction as the background
field $B_{z}$). Setting the velocity to anything else will not accurately
represent a propagating Alfvén wave, and may not represent _any_ physical way
the plasma could move. This freedom to incorrectly prescribe the boundary is
permitted when the boundary conditions are set arbitrarily in terms of the
primitive variables.

For more complicated scenarios (say, a multitude of overlapping waves or a
dynamically emerging flux rope) it is much more difficult, if not impossible,
to directly determine how the variables should be set in a self-consistent
manner. Depending on the state of the system, certain types of usual MHD
dynamics may not even be allowed. Again using a simple wave as an example,
suppose we set the boundary updates for each primitive variable,
$\bm{\mathbf{U}}\rightarrow\bm{\mathbf{U}}+\delta\bm{\mathbf{U}}$, to
represent an upward propagating slow mode wave, but did so in a location where
the plasma bulk velocity is directed downward at greater than the local Alfvén
speed (row 2 of Table 1). This is not a physically admissible situation, as
the plasma cannot propagate information in that way. Those boundary conditions
are therefore nonsense! _However, almost any MHD code will happily apply
nonsensical boundary conditions anyway_ , and this includes our own workhorse
code, LaRe3D. It is unknown how that mismatch will affect the results of the
simulation, and almost impossible to determine after the fact that nonsense
information has been injected into the simulation. In general, setting any of
the primitive variables independently will not represent a physically
admissible solution to the MHD equations in the boundary, and it is unclear
what implications this may have for the ongoing simulation. What is clear is
that, because MHD is a hyperbolic system that fundamentally describes the
propagation of information, bad information injected though the boundary in
this way will eventually corrupt the dynamics of the entire system. _Bad
boundary conditions lead to bad solar simulations!_

Our solution to this problem is to describe the plasma in the boundary purely
in terms of the ways it can actually move, which means implementing the
boundary conditions purely in terms of the characteristics of the system,
specifically the incoming characteristics, as described in §3. This is
wonderful because it forces all the primitive variables to be set in a self-
consistent manner, i.e., we are prevented from over-specifying the system. The
self-consistency arises because the characteristic modes are, effectively, the
allowed Lagrangian displacements of the MHD state variable, so updates to the
primitive variables calculated in terms of them are solutions to the MHD
equations. Put another way, setting the boundary conditions in terms of the
incoming characteristic derivatives implies that, in principle, one could
build a real machine in the laboratory that physically implements those
boundary conditions on the plasma.

However, this solution produces a fresh challenge, which is how to set the
incoming characteristic derivatives to represent the desired physics.
Fortunately, the data driving problem in solar physics likely represents the
best possible scenario for determining the desired physics because all the
hard work has already been done: the Sun already did it, and we observed it
happening! But typically, we have only _sparsely_ observed the Sun, either
temporally or in terms of the number of primitive variables that can be
reliably inferred from the observations. That fact lead us to the optimization
approach described in §5.

The next major task was to test all of the above developments, and for this we
chose to recreate a ground-truth simulation of an expanding spheromak. The
dynamics of the _GT_ simulation mimic some typical behaviors of emerging
active regions on the Sun: the magnetic field expands outward from the first
emergence location; the polarity inversion line in the extracted plane, where
the normal component of the field passes through zero, has a sigmoidal shape;
and the polarity inversion line appears to rotate during emergence. In this
regard, the expanding spheromak is a nice test case because it does contain
interesting and applicable dynamics while still allowing relatively simple
analysis. Notably, we did not include any inhomogeneous terms in this
validation, such as gravity or radiative transfer effects (although we did
carry those terms through the entire calculation as the variable
$\bm{\mathbf{D}}$). Future investigations will include these effects, which
are required for more realistic photosphere-to-corona simulations and when
applying our method to actual observational data.

In this investigation we performed an initial validation by varying the input
cadence of synthetic observations extracted from the _GT_ simulation. These
tests are similar in spirit to those performed by Leake et al. (2017), except
in our data driving scheme the boundary is constrained to evolve consistently
with both MHD equations and the interior of the driven simulation. Instead,
Leake et al. (2017) simply used linear interpolation in time of all of the
primitive variables extracted from their _GT_ simulation to set the values in
ghost cells of their driven simulation. As they showed, if the dynamic
evolution of the boundary is well approximated by a linear interpolation at
each location then this approach can produce good results. However, it was not
always obvious when the linear interpolation approximation was no longer
valid. In their case, they could introduce strong and completely spurious
currents into their simulation interior that dominate the system’s dynamics,
all while exactly matching the true solution at the driving times (this
produced substantially more “free magnetic energy” in the driven system than
existed in the true system). In contrast, a major strength of our optimization
approach is that it naturally makes known when the driving cadence is too
infrequent: our optimization based boundary condition will still produce a
_valid_ MHD solution for the estimated boundary evolution, but one that no
longer matches the observed evolution of the boundary. When present, this
discrepancy amounts to a “cry for help” from the _Data-Driven_ simulation: it
is being asked to conform to a dynamical trajectory disallowed by the physical
laws that govern our universe.

Our initial test cases revealed the transition from good to poor
reconstructions for driving cadences of $\Delta t\gtrsim 50\delta t$, which
was ultimately determined by length and velocity scales present in the system
dynamics (§7.1). At that low sampling rate, the driving layer was so sparsely
sampled in time that a linear interpolation between the two observed states
was no longer close to a valid solution of the MHD equations. Instead, our
data driving method produced optimal solutions that are consistent with the
MHD equations, in the sense that there could, in principle, be a series of
external states that drive the internal plasma through the optimal evolution.

As a second validation test, we directly compared the result of our data
driving method against a more common approach wherein $\rho$ and $\epsilon$
are held fixed to their initial values while the remaining primitive variables
are linearly interpolated between the observed states—we called these
“baseline” simulations. Many variations of this approach are found throughout
the literature (Bourdin et al., 2013; Galsgaard et al., 2015; Leake et al.,
2017; Jiang & Toriumi, 2020; Kaneko et al., 2021; Guo et al., 2021; Inoue et
al., 2023). In our tests, our characteristic based optimization method was
able to produce better results than the baseline approach. This test is the
closest approximation to the type of data driving that can currently be
achieved using the routine observations provided by SDO/HMI, which are the
best suited for studying active region evolution. Note that the linear
interpolation tests performed in Leake et al. (2017) took the extreme case
where all eight primitive variables were assumed to be known. We also
performed such tests, and like Leake et al. (2017), found that linear
interpolation can produce excellent results _provided that the evolution of
all the primitive variables in the driven boundary are already self-consistent
with each other and the interior evolution of the simulation_. If the
evolution is not self-consistent, either because the observational cadence is
too low and/or some variables are less well observed, then the fidelity of the
reconstruction will decrease.

The main issue with simple interpolation is that the resulting boundary
condition is over-determined: more primitive variables are independently and
inconsistently set than allowed by the MHD equations, i.e., more than the
number of incoming characteristics given in Table 1. This problem does not
arise when using high cadence input that is already self-consistent, as, for
instance, in the high cadence tests of Leake et al. (2017): even though those
authors always set the values of all primitive variables, their driving data
were pulled from a self-consistent MHD simulation. At high cadence, they
essentially provided _redundant_ information about the simultaneous evolution
of multiple primitive variables. The errors in their data driving tests only
grew significant when the information they provided was no longer self-
consistent due to interpolation times that exceeded the physically defined
critical cadence $\Delta t_{\text{crit}}$, as was the case for us.

This point is essentially a manifestation of Kirchoff’s integral theorem for
boundary conditions. If both the value and derivative of a function are set on
the boundary they are not generally consistent (for arbitrary values of both
values and derivatives). But if one sets both the values and derivatives
_self-consistently_ then the system is not actually overdetermined because the
specific combination of values and derivatives are just linear combinations of
the correct characteristics. The case of Leake et al. (2017) is an example of
this, where every variable and its derivative was set in the BC but they were
all pulled self-consistently from an MHD simulation.

Now consider how this issue of over-specifying the system applies to a case of
real data driving, for example, using primitive variables derived from
spectropolarimetric observations of the photosphere. In the limit that all
variables are well observed at high enough temporal and spatial resolution and
have co-spatial origin, a simple interpolation between observed states will
give a correct solution to the system’s evolution. It will not be over-
specified because the evolution of all the variables would be self-consistent
and include mutually redundant information. Real life seldom works that way,
but instead includes noise, bias, and other inconsistences. Dealing with such
issues is a major strength of our method based on the characteristics and
employing an optimization approach that prioritizes the validity of the
underlying equations. For the case of insufficient temporal data, where
interpolation between known states at distant times produces inconsistent
estimates for MHD states at intermediate times, our characteristic method
instead provides a plausible physical trajectory through state space between
the two known states, a trajectory constrained by the MHD equations
themselves.

Similarly, when some portion of the MHD state vector is either wholly unknown
or estimated at lower precision than others, the optimization method should
again provide an optimal trajectory in a least-squares sense between estimated
states, constrained to be as consistent as possible with the observed states.
In such cases it may be necessary to resort to an ad hoc assumption about how
a variable evolves. While this is not ideal, it can still be done in a way
that prioritizes a self-consistent evolution of the entire MHD state vector by
combining the characteristics with an optimization approach. For instance, if
we are missing information about the time evolution of the density, we could
assume a model where it exponentially relaxes back to a background stratified
value over some timescale $t_{R}:$ $\rho(t)\rightarrow\rho_{0}e^{-t/t_{R}}$.
The optimization procedure would then accommodate this prescription, but only
to the extent that the evolution of $\rho$ was consistent with that of the
other, better known variables. Our application of the weighted SVD solution,
described in §5 and tested in §7.2, demonstrated the validity of this
approach; refinements to it are forthcoming. Finally, by applying an
optimization technique to this problem, additional constraints could also be
imposed. The most immediately useful constraints for solar physics
applications are independent estimates for the energy, helicity, or mass
fluxes through the boundary. We will explore the capabilities of our method in
these cases in future work.

Other groups have performed validation studies for data-driven MHD simulations
that are similar to our own, most notably the test of several data driving
methods in Toriumi et al. (2020) and the follow-up study focused on the DARE
code in Jiang & Toriumi (2020). In this instance, the DARE code did not use a
characteristic method for the driven boundary condition, although these same
authors have done so in the past (Jiang et al., 2016). These authors tested
DARE’s ability to recreate the coronal portion of a self-consistent, ab initio
ground truth simulation of a buoyant flux rope emerging through a
gravitationally stratified atmosphere. As stated in §1, they produced a decent
reconstruction of the magnetic field when using data extracted from a height
with minimal Lorentz forces to drive the MHD simulation, i.e., data from the
coronal base. However, the technique struggled when driven using data from
closer to the photosphere. They pointed to the presence of strong Lorentz
forces in the photosphere, which are approximately balanced by pressure
gradients in the ground truth simulation but not in the data-driven
simulation. They state that these photospheric Lorentz forces are too large to
be effectively handled by DARE, but we find this to be a puzzling argument:
DARE is an MHD-based code and Lorentz forces are certainly within the scope of
MHD, so DARE should be able to handle such a case. However, we do note that
their ground truth simulation, which includes gravitational stratification and
an emerging flux rope, is a much more difficult challenge to reproduce than
our pressure-gradient-driven expanding spheromak. In a future work we will
apply our method to the Toriumi et al. (2020) test case, which will hopefully
provide another data point for when and where data-driven MHD approaches
succeed or fail.

A number of authors of coronal simulations do use a form of the
characteristics to set at least some of their boundary conditions. A long
series of work originating with Nakagawa (1980) uses a technique they name
Projected Normal Characteristics (PNC) to set boundary conditions (Nakagawa,
1980; Nakagawa et al., 1987; Wu & Wang, 1987; Hayashi, 2005; Wu et al., 2006;
Jiang et al., 2011; Yang et al., 2012; Hayashi, 2013; Feng et al., 2015; Jiang
et al., 2016; Feng et al., 2017; Hayashi et al., 2018, 2019, 2022), while
recently a separate characteristic based method was developed by Sarp Yalim et
al. (2017) and used in Singh et al. (2018) and Sarp Yalim et al. (2020). The
PNC method is used for multiple purposes, often as a form of nonreflecting
boundary condition where the evolution of the boundary is calculated by
setting the incoming characteristics to cancel the effects of transverse
terms. In our notation, this is implemented as
$\bm{\mathbf{L}}_{z\mathcal{I}}=-\mathbb{S}_{z\mathcal{I}}^{-1}\bm{\mathbf{C}}$,
and again, we present an extended discussion of this and other nonreflecting
boundary conditions in Paper II. For instance, in Appendix B of Hayashi et al.
(2018) they describe using the PNC approach for the side and top boundaries of
a driven simulation, but use a non-characteristic based method for the bottom,
driven boundary.

The PNC method has also been applied to the problem of data driving (e.g.,
Hayashi, 2005; Wu et al., 2006; Yang et al., 2012; Feng et al., 2015; Jiang et
al., 2016). Most cite Wu et al. (2006) as the originator of this approach, and
that approach seems to differ from what we have described in this paper. In
particular, in the discussion in Wu et al. (2006) following their Equation
(A9) they appear to rely on a one-to-one mapping between between each incoming
characteristic and the value of a single primitive variable on the boundary.
This one-to-one mapping is also present in the work of Sarp Yalim et al.
(2017), Singh et al. (2018), and Sarp Yalim et al. (2020); see, e.g., Equation
(11) in the first of these, where the number of primitive variables to be set
arbitrarily is assumed to be equal to the number of incoming characteristics.
This one-to-one mapping runs counter to our approach, where the values of
multiple primitive variables are coupled together through the eigenstructure
of each characteristic mode and all incoming characteristic derivatives are
solved for simultaneously, as we discussed in Section 2 in connection with
Table 1.

What is clear is that the application of the PNC method to data driving is
quite challenging. The introduction of Jiang et al. (2021) describes some of
the issues these authors have encountered, and highlight two difficulties in
particular: (1) dealing with situations where the number of variables they
wish to set using observations does not precisely match the number of incoming
characteristics; and (2) dealing with locations near a PIL, where the normal
magnetic field passes through zero. The approach we have described does not
encounter issues in either of these cases (at least for the tests we have
performed so far), which reinforces the notion that these authors’ use of the
characteristics to implement boundary conditions differs from our own.

A related use of the characteristics that has appeared in several previous
solar simulations is to use a single incoming characteristic derivative to
enforce a single specific boundary condition, for instance, constant mass flux
(Hayashi, 2005; Lionello et al., 2013; Feng et al., 2015). While this is valid
approach in the regime in which it was applied (solar wind simulations with an
inner boundary condition located outside the solar Alfvén surface, so that the
wind is always super Alfvénic and directed into the simulation), these types
of assumptions will likely conflict with observational data in the sub-
Alfvénic corona at some location on the boundary. This may partially explain
the difficulties discussed in the introduction of Jiang et al. (2021).

In contrast to data driving, most previous formulations of boundary conditions
in terms of the characteristics are mostly concerned with getting information
_out_ through the simulation boundary rather than in. The goal in these cases
is to implement a nonreflecting boundary condition that leaves the ongoing
simulation itself relatively unmolested (Hedstrom, 1979; Thompson, 1990;
Grappin et al., 2000; Gudiksen et al., 2011). Achieving this actually makes
rather strong assumptions about the state of the external universe, but those
assumptions are not always clearly specified. In particular, we have used the
characteristic formulation of general boundary conditions described in this
paper to test several standard flavors of “nonreflecting” boundaries and found
that, while they all work reasonably well for transient disturbances such as
waves and shocks, they can fail spectacularly when attempting to advect larger
scale structures out of the simulation (for instance, a spheromak). For these
results, we refer the reader to an extensive discussion of the topic in Paper
II of this series (Kee et al., 2023).

As the previous example highlights, the formalism we have developed has
application beyond data-driven boundary conditions, and indeed, beyond
boundary conditions themselves. One application has already been demonstrated,
which is the ability to couple together numerical codes that are based on
different numerical schemes, as was done to couple the finite-volume wave-
splitting CHAR code to the Lagrangian-remap LaRe3D code. A closely related
application is to use CHAR to implement data driving in the middle of a
simulation, which would represent an important step towards using multi-height
observations of the Sun for data assimilation. This idea could possibly be
extended to couple together codes that implement different physics, for
instance to create a bi-directionally coupled PIC and fluid reconnection
simulation. Coupling in terms of the characteristics would ensure that, at
least at the level of the MHD approximation, information is correctly
exchanged between the two systems. Another application of the formalism we
have developed is to the analysis of MHD systems in general and simulations in
particular. The characteristics seem under-utilized for this task, at least in
the context of solar simulations. For instance, one could use the
characteristics to identify individual wave modes (in the infinitesimal limit)
or types of information (in the general case) propagating through a
simulation, study their interactions, and the conversion of wave energy from
one type to another. Again, subsequent investigations will explore these
topics.

##  9 Conclusions

We have developed a general framework for formulating boundary conditions in
MHD simulations based on the method of characteristics and have applied it to
the problem of photosphere-to-corona simulations driven by observations in the
photospheric layer. Our philosophy has been to prioritize consistency with the
underlying physical equations of MHD over strict adherence to temporally
interpolated observational data because, in practice, such data will contain a
variety of random and systematic errors. The result of our approach is a data-
driven boundary layer that is a valid solution to the time-integrated MHD
equations which matches the observed boundary evolution as closely as
possible. This holds to the extent that the numerically discretized equations
match the true equations.

This investigation provides a detailed and practical guide for the
implementation of our data-driven boundary conditions using the combination of
a characteristic-based finite volume scheme, CHAR, and an optimization
approach to produce an optimal data-driven boundary condition, CHAR-DDBC. The
optimization approach has the effect of continually “nudging” the boundary in
an MHD-allowed direction in state-space towards the best estimate of the state
derived from observations. Critically, our approach can directly handle
missing, noisy, and/or systematically biased data.

We validated our approach by reproducing a ground-truth LaRe3D simulation of
an expanding spheromak. This test problem was chosen because it is as simple
as possible while still retaining many dynamical properties in common with
emerging active regions on the Sun. The results were satisfyingly unexciting.
When high cadence synthetic observations were provided to the CHAR-DDBC, the
output matched the ground truth simulation at the expected level of accuracy,
given the different numerical schemes of the two codes. On the other hand,
when the input cadence was lowered, so that a great deal of (physically
incorrect) interpolation was used between observed boundary states, our
approach lost the ability to accurately replicate the ground truth simulation.
Further, the loss of fidelity occurred at a previously predicted critical
cadence determined by physical properties of the system dynamics.

Our validation also included the case where we removed all information about
the evolution of some variables, namely, the mass and internal energy
densities. This mirrors the practical application to SDO/HMI data, which will
remain the most promising observational source for data driving for the
foreseeable future. In this case, setting up an SVD problem to find an
optimized solution for incoming characteristic derivatives was able to produce
a solution with improved fidelity compared to holding those variables constant
in the lower boundary. This test demonstrates that the characteristics can be
used in practical applications to constrain variables about which one has
little reliable observational information.

Our approach to boundary conditions, and their implementation in the CHAR
code, has uses beyond data-driven boundary conditions. For one, the use of
characteristics for analyzing simulation interiors seems underappreciated.
Because the characteristics are associated with the distinct types of
information that propagate through a plasma, they could be used, for example,
to track energy associated with different modes as they propagate through a
system, and thus identify locations of mode conversion or dissipation.

As a second example, we have already demonstrated that MHD codes with
different numerical schemes can be coupled together with CHAR. This ability
could be taken further, for instance to join together PIC codes with MHD
codes, or to incorporate multi-height observations of the solar atmosphere
into a single driven simulation.

Third, as we have stated, _any_ valid boundary conditions can be represented
in characteristic form. In particular, arbitrary open boundary conditions can
be represented in characteristic form. Obviously, this has long been
recognized (Hedstrom, 1979), but mostly applied to cases where one wishes
waves to leave a simulation with minimal back-reaction on the internal
dynamics. But for photosphere-to-corona simulations, we often want to study
the eruptive phenomena of flares and coronal mass ejections that produce great
changes both locally and in the far-field regions. It is unclear what boundary
conditions for the side and top boundaries are best suited to that task. We
refer the reader to Paper II for an initial exploration of that topic, in
which we implement and compare two common versions of nonreflecting boundary
conditions using CHAR. In this light, data-driven open boundary conditions are
far easier to defend than other open boundary conditions, which essentially
impose ad hoc constraints at each location in space and time. In contrast, for
data-driven boundary conditions, the data themselves guide what the spatially
and temporally varying boundary condition must be in order to accommodate the
ongoing simulation, the observations at the boundary, and the propagation of
information governed by the MHD equations.

Finally, for data driving itself, the framework we have developed in this
paper can be greatly extended. The current tests are the initial proof-of-
concept, but have already demonstrated the ability of combining the
characteristics with an optimization approach to back out the evolution of
primitive variables about which one nominally knows nothing. Future work will
apply the method to actual observations of solar active regions, extend to a
nonlinear optimization method, and include other physical constraints.

This work is supported by the National Solar Observatory, the Office of Naval
Research, and the NASA HSR and LWS Programs. MGL acknowledges support from the
Office of Naval Research. JEL, MGL, PWS, and LAT acknowledge support from NASA
grant NNH16ZDA001N-LWS “Implementing and Evaluating a Vector-Magnetogram-
Driven Magnetohydrodynamic Model of the Magnetic Field in the Low Solar
Atmosphere.” MGL acknowledges support from NASA grants NNH17ZDA001N-LWS
“Investigating Magnetic Flux Emergence with Modeling and Observations to
Understand the Onset of Major Solar Eruptions,” NNH18ZDA001N-HSR
“Investigating Magnetic Flux Rope Emergence as the Source of Flaring Activity
in Delta-Spot Active Regions,” and NNH20ZDA001N-HSR “Investigating the
Influence of Coronal Magnetic Geometry on the Acceleration of the Solar Wind.”
PWS and JEL acknowledge support from the NASA Internal Science Funding Model
(H-ISFM) program “Magnetic Energy Buildup and Explosive Release in the Solar
Atmosphere.” This research has made use of NASA’s Astrophysics Data System.

## References

  * Abramowitz & Stegun (1964) Abramowitz, M., & Stegun, I. A. 1964, Handbook of mathematical functions : with formulas, graphs, and mathematical tables (U.S. Dept. of Commerce, National Bureau of Standards)
  * Alissandrakis (1981) Alissandrakis, C. E. 1981, A&A, 100, 197
  * Altschuler & Newkirk (1969) Altschuler, M. D., & Newkirk, G. 1969, Sol. Phys., 9, 131, doi: 10.1007/BF00145734
  * Aly (1989) Aly, J. J. 1989, Sol. Phys., 120, 19, doi: 10.1007/BF00148533
  * Anderson et al. (1999) Anderson, E., Bai, Z., Bischof, C., et al. 1999, LAPACK Users’ Guide, 3rd edn. (Philadelphia, PA: Society for Industrial and Applied Mathematics)
  * Arber et al. (2001) Arber, T. D., Longbottom, A. W., Gerrard, C. L., & Milne, A. M. 2001, Journal of Computational Physics, 171, 151, doi: 10.1006/jcph.2001.6780
  * Aschwanden (2013) Aschwanden, M. J. 2013, Sol. Phys., 287, 323, doi: 10.1007/s11207-012-0069-7
  * Aschwanden (2016) —. 2016, ApJS, 224, 25, doi: 10.3847/0067-0049/224/2/25
  * Asensio Ramos & Díaz Baso (2019) Asensio Ramos, A., & Díaz Baso, C. J. 2019, A&A, 626, A102, doi: 10.1051/0004-6361/201935628
  * Asensio Ramos et al. (2017) Asensio Ramos, A., Requerey, I. S., & Vitas, N. 2017, A&A, 604, A11, doi: 10.1051/0004-6361/201730783
  * Bogdan & Low (1986) Bogdan, T. J., & Low, B. C. 1986, ApJ, 306, 271, doi: 10.1086/164341
  * Boocock & Tsiklauri (2019) Boocock, C. M., & Tsiklauri, D. 2019, A&A, 625, A47, doi: 10.1051/0004-6361/201834684
  * Borrero & Pastor Yabar (2023) Borrero, J. M., & Pastor Yabar, A. 2023, A&A, 669, A122, doi: 10.1051/0004-6361/202244716
  * Borrero et al. (2019) Borrero, J. M., Pastor Yabar, A., Rempel, M., & Ruiz Cobo, B. 2019, A&A, 632, A111, doi: 10.1051/0004-6361/201936367
  * Bourdin et al. (2013) Bourdin, P. A., Bingert, S., & Peter, H. 2013, A&A, 555, A123, doi: 10.1051/0004-6361/201321185
  * Brio & Wu (1988) Brio, M., & Wu, C. C. 1988, Journal of Computational Physics, 75, 400, doi: 10.1016/0021-9991(88)90120-9
  * Carlsson & Stein (1994) Carlsson, M., & Stein, R. F. 1994, in Chromospheric Dynamics, ed. M. Carlsson, 47
  * Casini et al. (2017) Casini, R., White, S. M., & Judge, P. G. 2017, Space Sci. Rev., 210, 145, doi: 10.1007/s11214-017-0400-6
  * Cheung & DeRosa (2012) Cheung, M. C. M., & DeRosa, M. L. 2012, ApJ, 757, 147, doi: 10.1088/0004-637X/757/2/147
  * Childs et al. (2012) Childs, H., Brugger, E., Whitlock, B., et al. 2012, in High Performance Visualization–Enabling Extreme-Scale Scientific Insight, 357–372, doi: 10.1201/b12985
  * da Silva Santos et al. (2022) da Silva Santos, J. M., Danilovic, S., Leenaarts, J., et al. 2022, A&A, 661, A59, doi: 10.1051/0004-6361/202243191
  * De Rosa et al. (2009) De Rosa, M. L., Schrijver, C. J., Barnes, G., et al. 2009, ApJ, 696, 1780, doi: 10.1088/0004-637X/696/2/1780
  * Dewar (1970) Dewar, R. L. 1970, Physics of Fluids, 13, 2710, doi: 10.1063/1.1692854
  * Dewar et al. (2015) Dewar, R. L., Yoshida, Z., Bhattacharjee, A., & Hudson, S. R. 2015, Journal of Plasma Physics, 81, doi: 10.1017/S0022377815001336
  * Dima & Schad (2020) Dima, G. I., & Schad, T. A. 2020, ApJ, 889, 109, doi: 10.3847/1538-4357/ab616f
  * Feng et al. (2017) Feng, X., Li, C., Xiang, C., et al. 2017, ApJS, 233, 10, doi: 10.3847/1538-4365/aa957a
  * Feng et al. (2015) Feng, X., Ma, X., & Xiang, C. 2015, Journal of Geophysical Research (Space Physics), 120, 10, doi: 10.1002/2015JA021911
  * Fisher et al. (2020) Fisher, G. H., Kazachenko, M. D., Welsch, B. T., et al. 2020, ApJS, 248, 2, doi: 10.3847/1538-4365/ab8303
  * Fisher et al. (2015) Fisher, G. H., Abbett, W. P., Bercik, D. J., et al. 2015, Space Weather, 13, 369, doi: 10.1002/2015SW001191
  * Galsgaard et al. (2015) Galsgaard, K., Madjarska, M. S., Vanninathan, K., Huang, Z., & Presmann, M. 2015, A&A, 584, A39, doi: 10.1051/0004-6361/201526339
  * Goedbloed et al. (2010) Goedbloed, J. P., Keppens, R., & Poedts, S. 2010, Advanced Magnetohydrodynamics (Cambridge University Press), doi: 10.1017/CBO9781139195560
  * Goedbloed & Poedts (2004) Goedbloed, J. P. H., & Poedts, S. 2004, Principles of Magnetohydrodynamics: With Applications to Laboratory and Astrophysical Plasmas (Cambridge University Press), doi: 10.1017/CBO9780511616945
  * Goldstraw et al. (2018) Goldstraw, E. E., Hood, A. W., Browning, P. K., & Cargill, P. J. 2018, A&A, 610, A48, doi: 10.1051/0004-6361/201731069
  * Grappin et al. (2008) Grappin, R., Aulanier, G., & Pinto, R. 2008, A&A, 490, 353, doi: 10.1051/0004-6361:200810386
  * Grappin et al. (2000) Grappin, R., Léorat, J., & Buttighoffer, A. 2000, A&A, 362, 342
  * Gudiksen et al. (2011) Gudiksen, B. V., Carlsson, M., Hansteen, V. H., et al. 2011, A&A, 531, A154, doi: 10.1051/0004-6361/201116520
  * Guo et al. (2021) Guo, Y., Zhong, Z., Ding, M. D., et al. 2021, ApJ, 919, 39, doi: 10.3847/1538-4357/ac10c8
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * Harvey et al. (1996) Harvey, J. W., Hill, F., Hubbard, R. P., et al. 1996, Science, 272, 1284, doi: 10.1126/science.272.5266.1284
  * Hayashi (2005) Hayashi, K. 2005, ApJS, 161, 480, doi: 10.1086/491791
  * Hayashi (2013) —. 2013, Journal of Geophysical Research (Space Physics), 118, 6889, doi: 10.1002/2013JA018991
  * Hayashi et al. (2018) Hayashi, K., Feng, X., Xiong, M., & Jiang, C. 2018, ApJ, 855, 11, doi: 10.3847/1538-4357/aaacd8
  * Hayashi et al. (2019) —. 2019, ApJ, 871, L28, doi: 10.3847/2041-8213/aaffcf
  * Hayashi et al. (2022) Hayashi, K., Wu, C.-C., & Liou, K. 2022, ApJ, 930, 60, doi: 10.3847/1538-4357/ac6173
  * Hedstrom (1979) Hedstrom, G. W. 1979, Journal of Computational Physics, 30, 222, doi: 10.1016/0021-9991(79)90100-1
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Inoue et al. (2023) Inoue, S., Hayashi, K., & Miyoshi, T. 2023, ApJ, 946, 46, doi: 10.3847/1538-4357/ac9eaa
  * Jaeggli et al. (2012) Jaeggli, S. A., Lin, H., & Uitenbroek, H. 2012, ApJ, 745, 133, doi: 10.1088/0004-637X/745/2/133
  * Jeffrey & Taniuti (1964) Jeffrey, A., & Taniuti, T. 1964, Non-linear wave propagation (Academic Press Inc.)
  * Jiang et al. (2021) Jiang, C., Bian, X., Sun, T., & Feng, X. 2021, Frontiers in Physics, 9, 224, doi: 10.3389/fphy.2021.646750
  * Jiang et al. (2011) Jiang, C., Feng, X., Fan, Y., & Xiang, C. 2011, ApJ, 727, 101, doi: 10.1088/0004-637X/727/2/101
  * Jiang & Toriumi (2020) Jiang, C., & Toriumi, S. 2020, ApJ, 903, 11, doi: 10.3847/1538-4357/abb5ac
  * Jiang et al. (2016) Jiang, C., Wu, S. T., Feng, X., & Hu, Q. 2016, Nature Communications, 7, 11522, doi: 10.1038/ncomms11522
  * Kaneko et al. (2021) Kaneko, T., Park, S.-H., & Kusano, K. 2021, ApJ, 909, 155, doi: 10.3847/1538-4357/abe414
  * Kee et al. (2023) Kee, N. D., Tarr, L. A., Linton, M., Schuck, P., & Leake, J. 2023, ApJ
  * Kolmogorov (1941) Kolmogorov, A. 1941, Akademiia Nauk SSSR Doklady, 30, 301
  * Kramar et al. (2016) Kramar, M., Lin, H., & Tomczyk, S. 2016, ApJ, 819, L36, doi: 10.3847/2041-8205/819/2/L36
  * Landi et al. (2005) Landi, S., Velli, M., & Einaudi, G. 2005, ApJ, 624, 392, doi: 10.1086/428822
  * Leake et al. (2022) Leake, J. E., Linton, M. G., & Antiochos, S. K. 2022, ApJ, 934, 10, doi: 10.3847/1538-4357/ac74b7
  * Leake et al. (2017) Leake, J. E., Linton, M. G., & Schuck, P. W. 2017, ApJ, 838, 113, doi: 10.3847/1538-4357/aa6578
  * Li et al. (2021) Li, W., Li, M., Wang, K., et al. 2021, ApJ, 913, 135, doi: 10.3847/1538-4357/abfa97
  * Lionello et al. (2013) Lionello, R., Downs, C., Linker, J. A., et al. 2013, ApJ, 777, 76, doi: 10.1088/0004-637X/777/1/76
  * Low (1985) Low, B. C. 1985, ApJ, 293, 31, doi: 10.1086/163211
  * Lumme et al. (2022) Lumme, E., Pomoell, J., Price, D. J., et al. 2022, A&A, 658, A200, doi: 10.1051/0004-6361/202038744
  * Mackay & van Ballegooijen (2006) Mackay, D. H., & van Ballegooijen, A. A. 2006, ApJ, 641, 577, doi: 10.1086/500425
  * Malanushenko et al. (2012) Malanushenko, A., Schrijver, C. J., DeRosa, M. L., Wheatland, M. S., & Gilchrist, S. A. 2012, ApJ, 756, 153, doi: 10.1088/0004-637X/756/2/153
  * Maltby et al. (1986) Maltby, P., Avrett, E. H., Carlsson, M., et al. 1986, ApJ, 306, 284, doi: 10.1086/164342
  * Mathews et al. (2022) Mathews, N. H., Flyer, N., & Gibson, S. E. 2022, Journal of Computational Physics, 462, 111214, doi: 10.1016/j.jcp.2022.111214
  * Morse & Feshbach (1953) Morse, P. M., & Feshbach, H. 1953, Methods of theoretical physics, Part I (New York: McGraw-Hill)
  * Morse & Feshbach (1953) Morse, P. M., & Feshbach, H. 1953, Methods of Theoretical Physics. Part II (New York: McGraw-Hill)
  * Nakagawa (1980) Nakagawa, Y. 1980, ApJ, 240, 275, doi: 10.1086/158232
  * Nakagawa et al. (1987) Nakagawa, Y., HU, Y. Q., & Wu, S. T. 1987, å, 179, 354
  * Nitta et al. (2001) Nitta, S., Tanuma, S., Shibata, K., & Maezawa, K. 2001, ApJ, 550, 1119, doi: 10.1086/319774
  * Pesnell et al. (2012) Pesnell, W. D., Thompson, B. J., & Chamberlin, P. C. 2012, Sol. Phys., 275, 3, doi: 10.1007/s11207-011-9841-3
  * Peter et al. (2015) Peter, H., Warnecke, J., Chitta, L. P., & Cameron, R. H. 2015, A&A, 584, A68, doi: 10.1051/0004-6361/201527057
  * Poinsot & Lele (1992) Poinsot, T. J., & Lele, S. K. 1992, Journal of Computational Physics, 101, 104, doi: 10.1016/0021-9991(92)90046-2
  * Pomoell et al. (2019) Pomoell, J., Lumme, E., & Kilpua, E. 2019, Sol. Phys., 294, 41, doi: 10.1007/s11207-019-1430-x
  * Price et al. (2020) Price, D. J., Pomoell, J., & Kilpua, E. K. J. 2020, A&A, 644, A28, doi: 10.1051/0004-6361/202038925
  * Quintero Noda et al. (2021) Quintero Noda, C., Barklem, P. S., Gafeira, R., et al. 2021, A&A, 652, A161, doi: 10.1051/0004-6361/202037735
  * Rimmele et al. (2020) Rimmele, T. R., Warner, M., Keil, S. L., et al. 2020, Sol. Phys., 295, 172, doi: 10.1007/s11207-020-01736-7
  * Roe & Balsara (1996) Roe, P. L., & Balsara, D. S. 1996, SIAM J. of Applied Mathematics, 56, 57, doi: 10.1137/S003613999427084X
  * Rosenbluth & Bussac (1979) Rosenbluth, M. N., & Bussac, M. N. 1979, Nuclear Fusion, 19, 489, doi: 10.1088/0029-5515/19/4/007
  * Sakurai (1989) Sakurai, T. 1989, Space Sci. Rev., 51, 11, doi: 10.1007/BF00226267
  * Sarp Yalim et al. (2017) Sarp Yalim, M., Pogorelov, N., & Liu, Y. 2017, in Journal of Physics Conference Series, Vol. 837, Journal of Physics Conference Series, 012015, doi: 10.1088/1742-6596/837/1/012015
  * Sarp Yalim et al. (2020) Sarp Yalim, M., Prasad, A., Pogorelov, N., Zank, G., & Hu, Q. 2020, in Journal of Physics Conference Series, Vol. 1620, Journal of Physics Conference Series, 012026, doi: 10.1088/1742-6596/1620/1/012026
  * Schad et al. (2023) Schad, T. A., Kuhn, J. R., Fehlmann, A., et al. 2023, ApJ, 943, 59, doi: 10.3847/1538-4357/acabbd
  * Scherrer et al. (2012) Scherrer, P. H., Schou, J., Bush, R. I., et al. 2012, Sol. Phys., 275, 207, doi: 10.1007/s11207-011-9834-2
  * Schmidt (1964) Schmidt, H. U. 1964, in NASA Special Publication, Vol. 50 (National Aeronautics and Space Administration, Science and Technical Information Division), 107
  * Schou et al. (2012) Schou, J., Scherrer, P. H., Bush, R. I., et al. 2012, Sol. Phys., 275, 229, doi: 10.1007/s11207-011-9842-2
  * Schuck (2008) Schuck, P. W. 2008, ApJ, 683, 1134, doi: 10.1086/589434
  * Singh et al. (2018) Singh, T., Yalim, M. S., & Pogorelov, N. V. 2018, ApJ, 864, 18, doi: 10.3847/1538-4357/aad3b4
  * Tarr et al. (2023) Tarr, L. A., Kee, N. D., Linton, M., Schuck, P., & Leake, J. 2023, ApJ
  * Thompson (1987) Thompson, K. W. 1987, Journal of Computational Physics, 68, 1, doi: 10.1016/0021-9991(87)90041-6
  * Thompson (1990) —. 1990, Journal of Computational Physics, 89, 439, doi: 10.1016/0021-9991(90)90152-Q
  * Tomczyk et al. (2008) Tomczyk, S., Card, G. L., Darnell, T., et al. 2008, Sol. Phys., 247, 411, doi: 10.1007/s11207-007-9103-6
  * Tomczyk et al. (2016) Tomczyk, S., Landi, E., Burkepile, J. T., et al. 2016, Journal of Geophysical Research (Space Physics), 121, 7470, doi: 10.1002/2016JA022871
  * Toriumi et al. (2020) Toriumi, S., Takasao, S., Cheung, M. C. M., et al. 2020, ApJ, 890, 103, doi: 10.3847/1538-4357/ab6b1f
  * Valori et al. (2007) Valori, G., Kliem, B., & Fuhrmann, M. 2007, Sol. Phys., 245, 263, doi: 10.1007/s11207-007-9046-y
  * Verniero et al. (2018) Verniero, J. L., Howes, G. G., & Klein, K. G. 2018, Journal of Plasma Physics, 84, 905840103, doi: 10.1017/S0022377817001003
  * Vissers et al. (2021) Vissers, G. J. M., Danilovic, S., de la Cruz Rodríguez, J., et al. 2021, A&A, 645, A1, doi: 10.1051/0004-6361/202038900
  * Warnecke & Peter (2019) Warnecke, J., & Peter, H. 2019, A&A, 624, L12, doi: 10.1051/0004-6361/201935385
  * Welsch et al. (2004) Welsch, B. T., Fisher, G. H., Abbett, W. P., & Regnier, S. 2004, ApJ, 610, 1148, doi: 10.1086/421767
  * Wheatland (2007) Wheatland, M. S. 2007, Sol. Phys., 245, 251, doi: 10.1007/s11207-007-9054-y
  * Wiegelmann (2004) Wiegelmann, T. 2004, Sol. Phys., 219, 87, doi: 10.1023/B:SOLA.0000021799.39465.36
  * Wiegelmann & Inhester (2010) Wiegelmann, T., & Inhester, B. 2010, A&A, 516, A107, doi: 10.1051/0004-6361/201014391
  * Wiegelmann et al. (2012) Wiegelmann, T., Thalmann, J. K., Inhester, B., et al. 2012, Sol. Phys., 281, 37, doi: 10.1007/s11207-012-9966-z
  * Wolfram Research (2019) Wolfram Research, I. 2019, Mathematica, Version 12.0, Champaign, IL: Wolfram Research, Inc.
  * Wu et al. (2006) Wu, S. T., Wang, A. H., Liu, Y., & Hoeksema, J. T. 2006, ApJ, 652, 800, doi: 10.1086/507864
  * Wu & Wang (1987) Wu, S. T., & Wang, J. F. 1987, Computer Methods in Applied Mechanics and Engineering, 64, 267, doi: 10.1016/0045-7825(87)90043-0
  * Yadav et al. (2019) Yadav, R., de la Cruz Rodríguez, J., Díaz Baso, C. J., et al. 2019, A&A, 632, A112, doi: 10.1051/0004-6361/201936790
  * Yang et al. (2012) Yang, L. P., Feng, X. S., Xiang, C. Q., et al. 2012, Journal of Geophysical Research (Space Physics), 117, A08110, doi: 10.1029/2011JA017494
  * Yang et al. (1986) Yang, W. H., Sturrock, P. A., & Antiochos, S. K. 1986, ApJ, 309, 383, doi: 10.1086/164610
  * Yeates (2014) Yeates, A. R. 2014, Sol. Phys., 289, 631, doi: 10.1007/s11207-013-0301-0
  * Yeates (2022) —. 2022, Geophysical and Astrophysical Fluid Dynamics, 116, 305, doi: 10.1080/03091929.2021.2021197
  * Zhu & Wiegelmann (2018) Zhu, X., & Wiegelmann, T. 2018, ApJ, 866, 130, doi: 10.3847/1538-4357/aadf7f
  * Zhu & Wiegelmann (2019) —. 2019, A&A, 631, A162, doi: 10.1051/0004-6361/201936433

##  Appendix A Left and Right Eigenvectors of $\mathbb{A}_{z}$

The right eigenvectors, which satisify
$\mathbb{A}_{z}\cdot\bm{\mathbf{r}}_{\sigma}=\lambda_{\sigma}\bm{\mathbf{r}}_{\sigma}$,
are as follows:

$\bm{\mathbf{r}}_{1}=\left.\begin{aligned} \begin{pmatrix}0\\\ 0\\\ 0\\\ 0\\\
0\\\ 0\\\ 0\\\
1\end{pmatrix}\end{aligned}\right.,\bm{\mathbf{r}}_{2}=\left.\begin{aligned}
\begin{pmatrix}-\rho\\\ \epsilon\\\ 0\\\ 0\\\ 0\\\ 0\\\ 0\\\
0\end{pmatrix}\end{aligned}\right.,\bm{\mathbf{r}}_{3,4}^{\pm}=\left.\begin{aligned}
\begin{pmatrix}0\\\ 0\\\ \pm\beta_{y}\\\ \mp\beta_{x}\\\ 0\\\
-\beta_{y}\sqrt{\rho}s_{z}\\\ \beta_{x}\sqrt{\rho}s_{z}\\\
0\end{pmatrix}\end{aligned}\right.,\bm{\mathbf{r}}_{5,6}^{\pm}=\left.\begin{aligned}
\begin{pmatrix}\rho\alpha_{s}\\\ \frac{a^{2}}{\gamma}\alpha_{s}\\\
\pm\beta_{x}\alpha_{f}c_{f}s_{z}\\\ \pm\beta_{y}\alpha_{f}c_{f}s_{z}\\\ \pm
c_{s}\alpha_{s}\\\ -\beta_{x}\sqrt{\rho}a\alpha_{f}\\\
-\beta_{y}\sqrt{\rho}a\alpha_{f}\\\
0\end{pmatrix}\end{aligned}\right.,\bm{\mathbf{r}}_{7,8}^{\pm}=\left.\begin{aligned}
\begin{pmatrix}\rho\alpha_{f}\\\ \frac{a^{2}}{\gamma}\alpha_{f}\\\
\mp\beta_{x}\alpha_{s}c_{s}s_{z}\\\ \mp\beta_{y}\alpha_{s}c_{s}s_{z}\\\ \pm
c_{f}\alpha_{f}\\\ \beta_{x}\sqrt{\rho}a\alpha_{s}\\\
\beta_{y}\sqrt{\rho}a\alpha_{s}\\\ 0\end{pmatrix}\end{aligned}\right..$ (A1)

Likewise, the left eigenvectors, which satisfy
$\bm{\mathbf{l}}_{\sigma}^{T}\cdot\mathbb{A}_{z}=\lambda_{\sigma}\bm{\mathbf{l}}_{\sigma}^{T}$,
are

$\bm{\mathbf{l}}_{1}=\left.\begin{aligned} \begin{pmatrix}0\\\ 0\\\ 0\\\ 0\\\
0\\\ 0\\\ 0\\\
1\end{pmatrix}\end{aligned}\right.,\bm{\mathbf{l}}_{2}=\left.\begin{aligned}
\begin{pmatrix}\frac{1-\gamma}{\gamma\rho}\\\ \frac{1}{\gamma\epsilon}\\\ 0\\\
0\\\ 0\\\ 0\\\ 0\\\
0\end{pmatrix}\end{aligned}\right.,(\bm{\mathbf{l}}_{3,4}^{\pm})=\frac{1}{2}\left.\begin{aligned}
\begin{pmatrix}0\\\ 0\\\ \pm\beta_{y}\\\ \mp\beta_{x}\\\ 0\\\
-\frac{\beta_{y}s_{z}}{\sqrt{\rho}}\\\ \frac{\beta_{x}s_{z}}{\sqrt{\rho}}\\\
0\\\
\end{pmatrix}\end{aligned}\right.,(\bm{\mathbf{l}}_{5,6}^{\pm})=\frac{1}{2a^{2}}\left.\begin{aligned}
\begin{pmatrix}\frac{\alpha_{s}a^{2}}{\gamma\rho}\\\
\frac{\alpha_{s}a^{2}}{\gamma\epsilon}\\\ \pm\beta_{x}\alpha_{f}c_{f}s_{z}\\\
\pm\beta_{y}\alpha_{f}c_{f}s_{z}\\\ \pm c_{s}\alpha_{s}\\\
-\beta_{x}\frac{a\alpha_{f}}{\sqrt{\rho}}\\\
-\beta_{y}\frac{a\alpha_{f}}{\sqrt{\rho}}\\\
0\end{pmatrix}\end{aligned}\right.,(\bm{\mathbf{l}}_{7,8}^{\pm})=\frac{1}{2a^{2}}\left.\begin{aligned}
\begin{pmatrix}\frac{\alpha_{f}a^{2}}{\gamma\rho}\\\
\frac{\alpha_{f}a^{2}}{\gamma\epsilon}\\\ \mp\beta_{x}\alpha_{s}c_{s}s_{z}\\\
\mp\beta_{y}\alpha_{s}c_{s}s_{z}\\\ \pm c_{f}\alpha_{f}\\\
\beta_{x}\frac{a\alpha_{s}}{\sqrt{\rho}}\\\
\beta_{y}\frac{a\alpha_{s}}{\sqrt{\rho}}\\\
0\end{pmatrix}\end{aligned}\right..$ (A2)

In all these expressions, the first subscript refers to the lower sign, and
the second to the upper sign. The eigenvectors are biorthonormal in the sense
that
$\bm{\mathbf{l}}_{\sigma}^{T}\cdot\bm{\mathbf{r}}_{\zeta}=\delta_{\sigma,\zeta}$.
Note that, because $\mathbb{A}_{z}$ is not orthogonal matrix, the left and
right eigenvectors are elements of two distinct vector spaces.

The right eigenmatrix has columns given by the right eigenvectors, while the
left eigenmatrix has rows given by the left eigenvectors:

$\mathbb{S}_{z}=\left.\begin{bmatrix}\begin{array}[]{cccccccc}|&|&|&|&|&|&|&|\\\
\bm{\mathbf{r}}_{1}&\bm{\mathbf{r}}_{2}&\bm{\mathbf{r}}_{3}&\bm{\mathbf{r}}_{4}&\bm{\mathbf{r}}_{5}&\bm{\mathbf{r}}_{6}&\bm{\mathbf{r}}_{7}&\bm{\mathbf{r}}_{8}\\\
|&|&|&|&|&|&|&|\end{array}\end{bmatrix}\right.,\quad\mathbb{S}^{-1}_{z}=\left.\begin{bmatrix}\begin{array}[]{ccc}\rule[2.15277pt]{10.76385pt}{0.5pt}&\bm{\mathbf{l}}_{1}^{T}&\rule[2.15277pt]{10.76385pt}{0.5pt}\\\
\rule[2.15277pt]{10.76385pt}{0.5pt}&\bm{\mathbf{l}}_{2}^{T}&\rule[2.15277pt]{10.76385pt}{0.5pt}\\\
\rule[2.15277pt]{10.76385pt}{0.5pt}&\bm{\mathbf{l}}_{3}^{T}&\rule[2.15277pt]{10.76385pt}{0.5pt}\\\
\rule[2.15277pt]{10.76385pt}{0.5pt}&\bm{\mathbf{l}}_{4}^{T}&\rule[2.15277pt]{10.76385pt}{0.5pt}\\\
\rule[2.15277pt]{10.76385pt}{0.5pt}&\bm{\mathbf{l}}_{5}^{T}&\rule[2.15277pt]{10.76385pt}{0.5pt}\\\
\rule[2.15277pt]{10.76385pt}{0.5pt}&\bm{\mathbf{l}}_{6}^{T}&\rule[2.15277pt]{10.76385pt}{0.5pt}\\\
\rule[2.15277pt]{10.76385pt}{0.5pt}&\bm{\mathbf{l}}_{7}^{T}&\rule[2.15277pt]{10.76385pt}{0.5pt}\\\
\rule[2.15277pt]{10.76385pt}{0.5pt}&\bm{\mathbf{l}}_{8}^{T}&\rule[2.15277pt]{10.76385pt}{0.5pt}\end{array}\end{bmatrix}\right..$
(A3)

Fully written out, the matrices are

$\mathbb{S}_{z}=\begin{bmatrix}0&-\rho&0&0&\rho\alpha_{s}&\rho\alpha_{s}&\rho\alpha_{f}&\rho\alpha_{f}\\\
0&\epsilon&0&0&\frac{a^{2}}{\gamma}\alpha_{s}&\frac{a^{2}}{\gamma}\alpha_{s}&\frac{a^{2}}{\gamma}\alpha_{f}&\frac{a^{2}}{\gamma}\alpha_{f}\\\
0&0&-\beta_{y}&+\beta_{y}&-\beta_{x}\alpha_{f}c_{f}s_{z}&+\beta_{x}\alpha_{f}c_{f}s_{z}&+\beta_{x}\alpha_{s}c_{s}s_{z}&-\beta_{x}\alpha_{s}c_{s}s_{z}\\\
0&0&+\beta_{x}&-\beta_{x}&-\beta_{y}\alpha_{f}c_{f}s_{z}&+\beta_{y}\alpha_{f}c_{f}s_{z}&+\beta_{y}\alpha_{s}c_{s}s_{z}&-\beta_{y}\alpha_{s}c_{s}s_{z}\\\
0&0&0&0&-c_{s}\alpha_{s}&+c_{s}\alpha_{s}&-c_{f}\alpha_{f}&+c_{f}\alpha_{f}\\\
0&0&-\beta_{y}\sqrt{\rho}s_{z}&-\beta_{y}\sqrt{\rho}s_{z}&-\beta_{x}\sqrt{\rho}a\alpha_{f}&-\beta_{x}\sqrt{\rho}a\alpha_{f}&\beta_{x}\sqrt{\rho}a\alpha_{s}&\beta_{x}\sqrt{\rho}a\alpha_{s}\\\
0&0&\beta_{x}\sqrt{\rho}s_{z}&\beta_{x}\sqrt{\rho}s_{z}&-\beta_{y}\sqrt{\rho}a\alpha_{f}&-\beta_{y}\sqrt{\rho}a\alpha_{f}&\beta_{y}\sqrt{\rho}a\alpha_{s}&\beta_{y}\sqrt{\rho}a\alpha_{s}\\\
1&0&0&0&0&0&0&0\end{bmatrix}$ (A4)

and

$\mathbb{S}_{z}^{-1}=\left.\begin{aligned} \begin{bmatrix}0&0&0&0&0&0&0&1\\\
-\frac{\gamma-1}{\gamma\rho}&\frac{1}{\gamma\epsilon}&0&0&0&0&0&0\\\
0&0&-\frac{\beta_{y}}{2}&\frac{\beta_{x}}{2}&0&-\frac{\beta_{y}s_{z}}{2\sqrt{\rho}}&\frac{\beta_{x}s_{z}}{2\sqrt{\rho}}&0\\\
0&0&\frac{\beta_{y}}{2}&-\frac{\beta_{x}}{2}&0&-\frac{\beta_{y}s_{z}}{2\sqrt{\rho}}&\frac{\beta_{x}s_{z}}{2\sqrt{\rho}}&0\\\
\frac{\alpha_{s}}{2\gamma\rho}&\frac{\alpha_{s}}{2\gamma\epsilon}&-\frac{c_{f}\alpha_{f}\beta_{x}s_{z}}{2a^{2}}&-\frac{c_{f}\alpha_{f}\beta_{y}s_{z}}{2a^{2}}&-\frac{c_{s}\alpha_{s}}{2a^{2}}&-\frac{\alpha_{f}\beta_{x}}{2a\sqrt{\rho}}&-\frac{\alpha_{f}\beta_{y}}{2a\sqrt{\rho}}&0\\\
\frac{\alpha_{s}}{2\gamma\rho}&\frac{\alpha_{s}}{2\gamma\epsilon}&\frac{c_{f}\alpha_{f}\beta_{x}s_{z}}{2a^{2}}&\frac{c_{f}\alpha_{f}\beta_{y}s_{z}}{2a^{2}}&\frac{c_{s}\alpha_{s}}{2a^{2}}&-\frac{\alpha_{f}\beta_{x}}{2a\sqrt{\rho}}&-\frac{\alpha_{f}\beta_{y}}{2a\sqrt{\rho}}&0\\\
\frac{\alpha_{f}}{2\gamma\rho}&\frac{\alpha_{f}}{2\gamma\epsilon}&\frac{c_{s}\alpha_{s}\beta_{x}s_{z}}{2a^{2}}&\frac{c_{s}\alpha_{s}\beta_{y}s_{z}}{2a^{2}}&-\frac{c_{f}\alpha_{f}}{2a^{2}}&\frac{\alpha_{s}\beta_{x}}{2a\sqrt{\rho}}&\frac{\alpha_{s}\beta_{y}}{2a\sqrt{\rho}}&0\\\
\frac{\alpha_{f}}{2\gamma\rho}&\frac{\alpha_{f}}{2\gamma\epsilon}&-\frac{c_{s}\alpha_{s}\beta_{x}s_{z}}{2a^{2}}&-\frac{c_{s}\alpha_{s}\beta_{y}s_{z}}{2a^{2}}&\frac{c_{f}\alpha_{f}}{2a^{2}}&\frac{\alpha_{s}\beta_{x}}{2a\sqrt{\rho}}&\frac{\alpha_{s}\beta_{y}}{2a\sqrt{\rho}}&0\\\
\end{bmatrix}\end{aligned}\right..$ (A5)

##  Appendix B Left and Right Eigenmatrices of $\mathbb{A}_{x}$

The coefficient matrix in the x-direction, $\mathbb{A}_{x}$, from Equation (2)
is given by

$\mathbb{A}_{x}=\left.\begin{aligned} \begin{bmatrix}v_{x}&0&\rho&0&0&0&0&0\\\
0&v_{x}&(\gamma-1)\epsilon&0&0&0&0&0\\\
\frac{(\gamma-1)\epsilon}{\rho}&(\gamma-1)&v_{x}&0&0&0&B_{y}/\rho&B_{z}/\rho\\\
0&0&0&v_{x}&0&0&-B_{x}/\rho&0\\\ 0&0&0&0&v_{x}&0&0&-B_{x}/\rho\\\
0&0&0&0&0&v_{x}&0&0\\\ 0&0&B_{y}&-B_{x}&0&0&v_{x}&0\\\
0&0&B_{z}&0&-B_{x}&0&0&v_{x}\\\ \end{bmatrix}\end{aligned}\right..$ (B1)

The right and left eigenmatrices of $\mathbb{A}_{x}$ are then calculated in a
similar fashion as for $\mathbb{A}_{z}$, except the correct projections of the
Alfvén velocities onto the normal and perpendicular directions must be used in
the definitions of the eigenvalues and vectors, i.e., $v_{z}\rightarrow
v_{x}$, $c_{a}\rightarrow\lvert b_{x}\rvert$ (including inside the definitions
of $c_{s,f}$ and $\alpha_{s,f}$), $b_{\perp}^{2}\rightarrow
b_{y}^{2}+b_{z}^{2}$, and the tangential ratios are now
$\beta_{y}=b_{y}/b_{\perp}$ and $\beta_{z}=b_{z}/b_{\perp}$. With those
definitions, the right and left eigenmatrices in the x-direction are

$\mathbb{S}_{x}=\left.\begin{aligned}
\begin{bmatrix}0&-\rho&0&0&\rho\alpha_{s}&\rho\alpha_{s}&\rho\alpha_{f}&\rho\alpha_{f}\\\
0&\epsilon&0&0&\frac{a^{2}\alpha_{s}}{\gamma}&\frac{a^{2}\alpha_{s}}{\gamma}&\frac{a^{2}\alpha_{f}}{\gamma}&\frac{a^{2}\alpha_{f}}{\gamma}\\\
0&0&0&0&-c_{s}\alpha_{s}&c_{s}\alpha_{s}&-c_{f}\alpha_{f}&c_{f}\alpha_{f}\\\
0&0&-\beta_{z}&\beta_{z}&-s_{x}c_{f}\alpha_{f}\beta_{y}&s_{x}c_{f}\alpha_{f}\beta_{y}&s_{x}c_{s}\alpha_{s}\beta_{y}&-s_{x}c_{s}\alpha_{s}\beta_{y}\\\
0&0&\beta_{y}&-\beta_{y}&-s_{x}c_{f}\alpha_{f}\beta_{z}&s_{x}c_{f}\alpha_{f}\beta_{z}&s_{x}c_{s}\alpha_{s}\beta_{z}&-s_{x}c_{s}\alpha_{s}\beta_{z}\\\
1&0&0&0&0&0&0&0\\\
0&0&-\sqrt{\rho}s_{x}\beta_{z}&-\sqrt{\rho}s_{x}\beta_{z}&-a\sqrt{\rho}\alpha_{f}\beta_{y}&-a\sqrt{\rho}\alpha_{f}\beta_{y}&a\sqrt{\rho}\alpha_{s}\beta_{y}&a\sqrt{\rho}\alpha_{s}\beta_{y}\\\
0&0&\sqrt{\rho}s_{x}\beta_{y}&\sqrt{\rho}s_{x}\beta_{y}&-a\sqrt{\rho}\alpha_{f}\beta_{z}&-a\sqrt{\rho}\alpha_{f}\beta_{z}&a\sqrt{\rho}\alpha_{s}\beta_{z}&a\sqrt{\rho}\alpha_{s}\beta_{z}\\\
\end{bmatrix}\end{aligned}\right.$ (B2)

and

$\mathbb{S}_{x}^{-1}=\left.\begin{aligned} \begin{bmatrix}0&0&0&0&0&1&0&0\\\
-\frac{\gamma-1}{\gamma\rho}&\frac{1}{\gamma\epsilon}&0&0&0&0&0&0\\\
0&0&0&-\frac{\beta_{z}}{2}&\frac{\beta_{y}}{2}&0&-\frac{s_{x}\beta_{z}}{2\sqrt{\rho}}&\frac{s_{x}\beta_{y}}{2\sqrt{\rho}}\\\
0&0&0&\frac{\beta_{z}}{2}&-\frac{\beta_{y}}{2}&0&-\frac{s_{x}\beta_{z}}{2\sqrt{\rho}}&\frac{s_{x}\beta_{y}}{2\sqrt{\rho}}\\\
\frac{\alpha_{s}}{2\gamma\rho}&\frac{\alpha_{s}}{2\gamma\epsilon}&-\frac{c_{s}\alpha_{s}}{2a^{2}}&-\frac{s_{x}c_{f}\alpha_{f}\beta_{y}}{2a^{2}}&-\frac{s_{x}c_{f}\alpha_{f}\beta_{z}}{2a^{2}}&0&-\frac{\alpha_{f}\beta_{y}}{2a\sqrt{\rho}}&-\frac{\alpha_{f}\beta_{z}}{2a\sqrt{\rho}}\\\
\frac{\alpha_{s}}{2\gamma\rho}&\frac{\alpha_{s}}{2\gamma\epsilon}&\frac{c_{s}\alpha_{s}}{2a^{2}}&\frac{s_{x}c_{f}\alpha_{f}\beta_{y}}{2a^{2}}&\frac{s_{x}c_{f}\alpha_{f}\beta_{z}}{2a^{2}}&0&-\frac{\alpha_{f}\beta_{y}}{2a\sqrt{\rho}}&-\frac{\alpha_{f}\beta_{z}}{2a\sqrt{\rho}}\\\
\frac{\alpha_{f}}{2\gamma\rho}&\frac{\alpha_{f}}{2\gamma\epsilon}&-\frac{c_{f}\alpha_{f}}{2a^{2}}&\frac{s_{x}c_{s}\alpha_{s}\beta_{y}}{2a^{2}}&\frac{s_{x}c_{s}\alpha_{s}\beta_{z}}{2a^{2}}&0&\frac{\alpha_{s}\beta_{y}}{2a\sqrt{\rho}}&\frac{\alpha_{s}\beta_{z}}{2a\sqrt{\rho}}\\\
\frac{\alpha_{f}}{2\gamma\rho}&\frac{\alpha_{f}}{2\gamma\epsilon}&\frac{c_{f}\alpha_{f}}{2a^{2}}&-\frac{s_{x}c_{s}\alpha_{s}\beta_{y}}{2a^{2}}&-\frac{s_{x}c_{s}\alpha_{s}\beta_{z}}{2a^{2}}&0&\frac{\alpha_{s}\beta_{y}}{2a\sqrt{\rho}}&\frac{\alpha_{s}\beta_{z}}{2a\sqrt{\rho}}\\\
\end{bmatrix}\end{aligned}\right..$ (B3)

The left and right eigenmatrices have been verified (using Mathematica) to be
biorthonormal in the sense that
$\mathbb{S}_{x}^{-1}\mathbb{S}_{x}=\mathbb{I}$.

##  Appendix C Left and Right Eigenmatrices of $\mathbb{A}_{y}$

The coefficient matrix in the y-direction, $\mathbb{A}_{y}$, from Equation (2)
is given by

$\mathbb{A}_{y}=\left.\begin{aligned} \begin{bmatrix}v_{y}&0&0&\rho&0&0&0&0\\\
0&v_{y}&0&(\gamma-1)\epsilon&0&0&0&0\\\ 0&0&v_{y}&0&0&-B_{y}/\rho&0&0\\\
\frac{(\gamma-1)\epsilon}{\rho}&(\gamma-1)&0&v_{y}&0&B_{x}/\rho&0&B_{z}/\rho\\\
0&0&0&0&v_{y}&0&0&-B_{y}/\rho\\\ 0&0&-B_{y}&B_{x}&0&v_{y}&0&0\\\
0&0&0&0&0&0&v_{y}&0\\\ 0&0&0&B_{z}&-B_{y}&0&0&v_{y}\\\
\end{bmatrix}\end{aligned}\right..$ (C1)

The right and left eigenmatrices of $\mathbb{A}_{y}$ are then calculated in a
similar fashion as for $\mathbb{A}_{z}$, except the correct projections of the
Alfvén velocities onto the normal and perpendicular directions must be used in
the definitions of the eigenvalues and vectors, i.e., $v_{z}\rightarrow
v_{y}$, $c_{a}\rightarrow\lvert b_{y}\rvert$ (including inside the definitions
of $c_{f,s}$ and $\alpha_{s,f}$), $b_{\perp}^{2}\rightarrow
b_{x}^{2}+b_{z}^{2}$, and the tangential ratios are
$\beta_{x}=b_{x}/b_{\perp}$ and $\beta_{z}=b_{z}/b_{\perp}$. With those
definitions, the right and left eigenmatrices are

$\mathbb{S}_{y}=\left.\begin{aligned}
\begin{bmatrix}0&-\rho&0&0&\rho\alpha_{s}&\rho\alpha_{s}&\rho\alpha_{f}&\rho\alpha_{f}\\\
0&\epsilon&0&0&\frac{a^{2}\alpha_{s}}{\gamma}&\frac{a^{2}\alpha_{s}}{\gamma}&\frac{a^{2}\alpha_{f}}{\gamma}&\frac{a^{2}\alpha_{f}}{\gamma}\\\
0&0&\beta_{z}&-\beta_{z}&-s_{y}c_{f}\alpha_{f}\beta_{x}&s_{y}c_{f}\alpha_{f}\beta_{x}&s_{y}c_{s}\alpha_{s}\beta_{x}&-s_{y}c_{s}\alpha_{s}\beta_{x}\\\
0&0&0&0&-c_{s}\alpha_{s}&c_{s}\alpha_{s}&-c_{f}\alpha_{f}&c_{f}\alpha_{f}\\\
0&0&-\beta_{x}&\beta_{x}&-s_{y}c_{f}\alpha_{f}\beta_{z}&s_{y}c_{f}\alpha_{f}\beta_{z}&s_{y}c_{s}\alpha_{s}\beta_{z}&-s_{y}c_{s}\alpha_{s}\beta_{z}\\\
0&0&\sqrt{\rho}s_{y}\beta_{z}&\sqrt{\rho}s_{y}\beta_{z}&-a\sqrt{\rho}\alpha_{f}\beta_{x}&-a\sqrt{\rho}\alpha_{f}\beta_{x}&a\sqrt{\rho}\alpha_{s}\beta_{x}&a\sqrt{\rho}\alpha_{s}\beta_{x}\\\
1&0&0&0&0&0&0&0\\\
0&0&-\sqrt{\rho}s_{y}\beta_{x}&-\sqrt{\rho}s_{y}\beta_{x}&-a\sqrt{\rho}\alpha_{f}\beta_{z}&-a\sqrt{\rho}\alpha_{f}\beta_{z}&a\sqrt{\rho}\alpha_{s}\beta_{z}&a\sqrt{\rho}\alpha_{s}\beta_{z}\\\
\end{bmatrix}\end{aligned}\right.,$ (C2)

and

$\mathbb{S}_{y}^{-1}=\left.\begin{aligned} \begin{bmatrix}0&0&0&0&0&0&1&0\\\
-\frac{\gamma-1}{\gamma\rho}&\frac{1}{\gamma\epsilon}&0&0&0&0&0&0\\\
0&0&\frac{\beta_{z}}{2}&0&-\frac{\beta_{x}}{2}&\frac{s_{y}\beta_{z}}{2\sqrt{\rho}}&0&-\frac{s_{y}\beta_{x}}{2\sqrt{\rho}}\\\
0&0&-\frac{\beta_{z}}{2}&0&\frac{\beta_{x}}{2}&\frac{s_{y}\beta_{z}}{2\sqrt{\rho}}&0&-\frac{`s_{y}\beta_{x}}{2\sqrt{\rho}}\\\