#  Simulating the photospheric to coronal plasma using magnetohydrodyanamic
characteristics I: data-driven boundary conditions

Lucas A. Tarr National Solar Observatory, 22 Ohi ‘ a Ku St, Makawao, HI, 96768
N. Dylan Kee National Solar Observatory, 22 Ohi ‘ a Ku St, Makawao, HI, 96768
Mark G. Linton U.S. Naval Research Laboratory, 4555 Overlook Ave SW,
Washington, DC 20375  Peter W. Schuck NASA Goddard Space Flight Center, 8800
Greenbelt Rd, Greenbelt, MD 20771  James E. Leake NASA Goddard Space Flight
Center, 8800 Greenbelt Rd, Greenbelt, MD 20771

###### Abstract

We develop a general description of how information propagates through a
magnetohydrodynamic (MHD) system based on the method of characteristics and
use that to formulate numerical boundary conditions (BCs) that are
intrinsically consistent with the MHD equations. Our formulation includes two
major advances for simulations of the Sun. First, we derive data-driven BCs
that optimally match the state of the plasma inferred from a time series of
observations of a boundary (e.g., the solar photosphere). Second, our method
directly handles random noise and systematic bias in the observations, and
finds a solution for the boundary evolution that is strictly consistent with
MHD and maximally consistent with the observations. We validate the method
against a Ground Truth (GT) simulation of an expanding spheromak. The data-
driven simulation can reproduce the GT simulation above the photosphere with
high fidelity when driven at high cadence. Errors progressively increase for
lower driving cadence until a threshold cadence is reached and the driven
simulation can no longer accurately reproduce the GT simulation. However, our
characteristic formulation of the BCs still requires adherence of the boundary
evolution to the MHD equations even when the driven solution departs from the
true solution in the driving layer. That increasing departure clearly
indicates when additional information at the boundary is needed to fully
specify the correct evolution of the system. The method functions even when no
information about the evolution of some variables on the lower boundary is
available, albeit with a further decrease in fidelity.

††software:  Correctness of the MHD eigensystem was verified using Mathematica
(Wolfram Research, 2019). The primary 3D MHD code, LaRe3D v.2.10 (Arber et
al., 2001), is written in Fortran 2003 and parallelized with MPI. 3D
renderings were generated using the VisIt visualization software (Childs et
al., 2012). The majority of the numerical analysis was carried out using the
NumPy package (Harris et al., 2020) and figures were prepared using the
Matplotlib library (Hunter, 2007) in python3. The dgelsd algorithm to solve
the SVD problem is included in LAPACK (Anderson et al., 1999).

ʻ ‘

##  1 Introduction

How is the three-dimensional (3D) state of the MHD plasma in the solar
atmosphere advanced in time given (i) the initial state of the atmosphere and
(ii) observations of the lower (photospheric) boundary at the later time? This
is a problem at the heart of solar physics research. First, all energy in the
chromosphere/corona originates in the convection zone and must transit the
photosphere. In plain language, all of the energy in photosphere-to-corona
models “originates” at the boundary! Thus, correctly solving for the MHD
evolution of the boundary is of paramount importance for correctly solving for
the evolution within the the volume of interest. Second, photospheric
emissions originate in a relatively narrow layer between optically thick and
optically thin radiative enviornments. Thus, the photospheric observables
derived from these emissions can be assumed to originate from a narrow height
range in the solar atmosphere—key for boundary driving. Third, the photosphere
represents the narrow boundary layer between the pressure dominated convection
zone physics and magnetically dominated coronal physics. Thus, photospheric
“boundary” driving of photosphere-to-corona models of the solar atmosphere is
an attractive approach for understanding the evolution of the corona because
of practicality—the photospheric observations are the most localized (and
perhaps the most reliable), and simplicity—the physics of the pressure-
dominated convection zone is lumped into the boundary conditions for
theoretical studies.

For these reasons, the rigorous time-advance of MHD models given only the
time-evolution of the boundary is of paramount importance for modeling coronal
activity such as the evolution of active regions, of coronal heating, of flare
energy release, of coronal mass ejections, and for space weather prediction.
Ideally, the solution to this problem—i.e., the values of the mass density
$\rho$, internal energy density $\epsilon$, magnetic field vector
$\bm{\mathbf{B}}$, and velocity vector $\bm{\mathbf{v}}$ at every location in
the photosphere-to-corona volume at the later time—should be exactly the same
as the case in which a more extended problem is solved wherein we model the
entire universe; that is, from the solar core, through the convection zone,
corona, heliosphere, and into interstellar space. The boundary driving must
properly encode the influence of the rest of the universe on our patch of
solar atmosphere, and the influence of our patch of solar atmosphere on the
rest of the universe. The daunting complex nature of this task is concealed by
its typical plain-language summary: we need to set the boundary conditions
“correctly.” This paper rigorously addresses “correct” in the MHD context.

The good news is that the photosphere is nearly continuously observed at
decent temporal cadence and spatial resolution by a variety of spectroscopic
and polarimetric instruments. Combined with a mix of models and assumptions,
these observations can be used to infer a good fraction of the MHD state
vector,
$\bm{\mathbf{U}}:=(\rho,\epsilon,\bm{\mathbf{v}},\bm{\mathbf{B}})^{T}$, in the
photosphere on a routine basis, especially within active regions, which are
the focus of this study. The Helioseismic and Magnetic Imager (Scherrer et
al., 2012; Schou et al., 2012) instrument onboard the Solar Dynamics
Observatory (Pesnell et al., 2012) regularly infers vector magnetic field
estimates of $\bm{\mathbf{B}}$ for the earth-facing side of the Sun at a 12
minute cadence and one arcsecond resolution ($\approx 730\,\mathrm{km}$ on the
Sun). Optical flow methods such as DAVE4VM (Schuck, 2008) or PDFI_SS (Fisher
et al., 2020) can provide estimates of plasma velocity $\bm{\mathbf{v}}$.
Thermodynamic variables such as density $\rho$ and internal energy density
$\epsilon$ (or, equivalently, the temperature $T$ or pressure $P$) can be
derived by combining model atmospheres with observations of the continuum
intensity, multiple spectral lines (Maltby et al., 1986), or using empirically
derived relations (Maltby et al., 1986; Jaeggli et al., 2012; Borrero et al.,
2019). For smaller spatial patches of the Sun, higher resolution and higher
sensitivity spectropolarimetric, multi-line observations are providing ever
better constrained estimates of the photospheric state vector (Yadav et al.,
2019; Quintero Noda et al., 2021; Vissers et al., 2021; da Silva Santos et
al., 2022), especially in combination with spatially-coupled inversion
approaches based on machine learning techniques (Asensio Ramos et al., 2017;
Asensio Ramos & Díaz Baso, 2019) and the newest generation of solar
telescopes, such as the Daniel K. Inouye Solar Telescope (Rimmele et al.,
2020).

This wealth of observational data has made photosphere-to-corona simulations
an appealing approach for studying the solar atmosphere, but such simulations
come with unique challenges in terms of their boundary conditions. Indeed,
even in the best case scenario where the entire MHD state vector may be
_estimated_ on the boundary from observations, this estimated state vector
will contain bias and noise and, in general, will not be completely consistent
with the internal state of the photosphere-to-corona model of the solar
atmosphere (Borrero et al., 2019). This same issue arises in purely
theoretical models where the model boundary is driven in a particular way
without regard to the internal state of the model, in which case the boundary
driving can over-determine the MHD model leading to an improper solution in
the volume of interest. Thus, the question arises, how best to use _estimated_
boundary data given the MHD state of the volume of interest?

Photosphere-to-corona simulations occupy a special space in plasma physics
modeling in terms of what boundary conditions are appropriate to the
situation. This is in stark contrast to many other subdisciplines of plasma
physics, where researchers are most often interested in the evolution of a
plasma that is energized in its interior, usually by instabilities, and the
boundary plays a secondary role. To give several examples:

  1. (i)

Local homogeneity and local isotropy—the volume of interest can be considered
a subvolume of a much larger turbulent system (Kolmogorov, 1941) so that
periodic boundary conditions, which pass energy out of the volume from one
boundary and inject it in on the opposite boundary, may be applicable.

  2. (ii)

Local plasma confinement—the volume of interest is either isolated from the
rest of the universe or the boundary is strictly controlled by external
electric circuits (Goedbloed & Poedts, 2004, §4.6.1) and therefore conducting
boundary conditions, which reflect energy back into the volume, may be
appropriate.

  3. (iii)

Spatially and temporally localized phenomena—the dynamics of interest are
transient, localized in space, and considered independent of the boundary
conditions (e.g., some reconnection scenarios (Nitta et al., 2001), bouyant
flux ropes (Leake et al., 2022), or interacting plasma wave packets (Verniero
et al., 2018)). Therefore, any numerically stable boundary condition may be
used provided all analysis of the interesting physics is carried out before
the boundaries can learn about the evolution of the plasma and act back on the
region of interest.

The last example (iii) makes explicit the idea that some information
propagates from a localized region of space, interacts with a distant boundary
(artificial or not), and then some other, possibly different, information
propagates back from the boundary into the localized subvolume of interest.

None of the above boundary conditions are appropriate for photosphere-to-
corona simulations because the corona is driven by energy that originates in
the convection zone and is transported through the photosphere. From the
perspective of a photosphere-to-corona simulation, the energy that ultimately
structures the corona originated in the simulation boundary. Therefore, for
photosphere-to-corona simulations, the boundary and its treatment have
elevated importance.

Despite the wealth of observational data mentioned previously, it remains an
open question how to correctly evolve an MHD simulation consistent with an
observed time series of the MHD state vector on the lower boundary. This is
true regardless of whether the entire state vector is known or only a subset.
Numerous attempts to solve related problems in the field of solar physics
stretch back decades (Schmidt, 1964; Altschuler & Newkirk, 1969; Nakagawa,
1980; Alissandrakis, 1981; Low, 1985; Bogdan & Low, 1986; Yang et al., 1986;
Nakagawa et al., 1987; Wu & Wang, 1987; Aly, 1989; Sakurai, 1989; Wiegelmann,
2004; Mackay & van Ballegooijen, 2006; Wu et al., 2006; Valori et al., 2007;
Wheatland, 2007; Jiang et al., 2011; Cheung & DeRosa, 2012; Malanushenko et
al., 2012; Aschwanden, 2013; Yeates, 2014; Fisher et al., 2015; Aschwanden,
2016; Goldstraw et al., 2018; Zhu & Wiegelmann, 2018; Pomoell et al., 2019;
Zhu & Wiegelmann, 2019; Boocock & Tsiklauri, 2019; Price et al., 2020; Lumme
et al., 2022; Mathews et al., 2022; Yeates, 2022). In brief, methods to
estimate the fully 3D state of the solar corona can be grouped into static
methods and evolving methods. Static methods construct an estimate of the
coronal state based on boundary observations at a single instance in time, and
evolving methods use a time series of observations at a boundary to
iteratively modify the volume state. Many of the methods just cited operate in
the limit of negligible pressure gradient forces and therefore only attempt to
determine the coronal magnetic field. On the other hand, the most
sophisticated methods are fully dynamic: they numerically integrate the full
MHD equations such that the state at one time is causally related to the state
at the next time through the combination of Maxwell’s equations and Newton’s
laws from which the MHD equations are derived.

In the present context we use the term “data-driven simulation” to mean a
dynamically evolving simulation that directly incorporates estimates for the
primitive variables (the elements of the MHD state vector $\bm{\mathbf{U}}$)
into the simulation boundary. We are agnostic as to where those estimates come
from: they could, for instance, be derived from observations, extracted from
another simulation, or imposed according to some analytic model. The problem
is challenging because it is possible to prescribe the boundary values in a
way that is fundamentally inconsistent with the MHD equations. In fact, this
will almost certainly be the case unless deliberate precautions are taken. For
example, most published data driving approaches hold some subset of the
primitive variables constant in time while directly setting others (Bourdin et
al., 2013; Galsgaard et al., 2015; Hayashi et al., 2018; Jiang & Toriumi,
2020; Kaneko et al., 2021; Guo et al., 2021; Inoue et al., 2023); typically,
some subset of the velocity and magnetic field vectors are time-varying while
the remaining velocity and magnetic field variables, the density, and the
pressure are held constant. Many variations on this approach exist, and while
the results may appear reasonable, in general _it cannot and does not evolve
the boundary in a way that is consistent with the MHD equations._ The reason
is that the evolution of the primitive variables are not independent from each
other but instead are coupled together by the MHD equations. Stated another
way, an arbitrary perturbation to the MHD state vector is typically not a
dynamical solution to the MHD equations; only displacements that can be
expressed as combinations of eigenmodes of the MHD equations are dynamically
allowed (see, e.g., the Lagrangian formulation of MHD in Dewar, 1970; Dewar et
al., 2015). As such, at present, none of the published data-driven MHD
approaches appear to drive their boundaries in a way that are guaranteed to be
fully consistent with the MHD equations. While this may be well recognized, it
is often not stated explicitly, which is why we have emphasized it here.

The inherent challenge of data-driven photosphere-to-corona simulations and
the shortcomings of the current methods were well illustrated in Toriumi et
al. (2020). They tested the fidelity of multiple coronal reconstruction
methods by having each method attempt to reproduce the coronal magnetic field
from an ab initio flux emergence simulation that spanned the convection zone
to the corona. Each method was given synthetic “observational” data consisting
of the time series of the MHD state vector extracted from a horizontal slice
through the emergence simulation’s photosphere. Toriumi et al. (2020) found
that all of the methods they tested for reconstructing the coronal field from
photospheric observations—static extrapolation, magnetofrictional, or data-
driven MHD—struggled to accurately reproduce the ground truth field on a
variety of metrics, including the extent of field expansion, the total
magnetic energy, and the energy in excess of the potential field (although a
magnetofrictional method was able to reproduce the relative magnetic helicity
and asymptotically match the total magnetic energy). Qualitatively and
quantitatively, the structure of field lines differed substantially between
the ground truth solution and all reproductions, even when a metric like the
total magnetic energy matched the ground truth value to a factor of two or so
(see especially Figures 3 and 5 of Toriumi et al. (2020)).

Jiang & Toriumi (2020) presented a follow-up to Toriumi et al. (2020) in which
they performed additional tests of one of the MHD based data-driven
simulations111The simulations used a variation of the code described in Jiang
et al. (2016). to study the effect of using data extracted at different
heights from ground truth simulation as a boundary for the data-driven
simulation. They tested their data driving approach using layers extracted
between the photosphere (which had a large Lorentz force imbalance) up to the
low corona (where the Lorentz force was approximately zero, a so-called “force
free” solution). They found much better agreement when driving from layers
with smaller residual Lorentz forces. However, it is unclear why an MHD-based
data driving approach was unable to reproduce the self-consistently generated
MHD state of the ground truth simulation, regardless of the size of the
Lorentz force; Lorentz forces are clearly an inherent and necessary part of
the MHD equations.

This is a puzzling and troubling result that is worth emphasizing: for the
test case of Toriumi et al. (2020) and Jiang & Toriumi (2020) the
“observational” data are perfect, the layer used for driving is pulled
directly from a self-consistent MHD simulation, and yet none of the methods
were able to fully reproduce the coronal field when driven from the
photosphere (the material plasma variables were not considered in that study).
The unavoidable conclusion is that, at a fundamental level, none of the tested
methods are fully consistent with the MHD equations _in the driven
photospheric boundary layer_. At some level these methods must not correctly
encode how information crosses the driving boundary. This emphasizes that if
the boundary is not evolved correctly then the interior will not be a valid
solution to the problem at hand, and may not be a valid solution to any MHD
system.

In the present work we derive boundary conditions for all primitive variables
that are necessarily consistent with the MHD equations, up to the level that
the numerical scheme matches the analytic equations. To do so we express the
evolution of the boundary using the characteristic form of MHD (see, e.g.,
Jeffrey & Taniuti, 1964; Goedbloed et al., 2010). The characteristic form of
MHD ultimately encodes how information (i.e., disturbances) propagate locally
through the plasma. This is possible because ideal MHD is a hyperbolic
system222MHD is not strictly hyperbolic: for certain states of the system, two
or more of the eigenmodes may become degenerate. However, those cases may
typically be resolved, see Roe & Balsara (1996). and therefore has propagating
solutions, each of which describe the directional propagation of information
(Goedbloed et al., 2010). It is therefore possible to decompose the dynamics
of the boundary into a portion that is due only to information originating
within the simulation and a portion that is due only to information
originating from outside the simulation, i.e., in the external universe. The
observed dynamics of the boundary are due to the superposition of these two
sources of propagating information through the boundary plus any information
propagating within the boundary itself, i.e., in the transverse direction. The
transverse portion, arising from transverse derivatives, can be calculated
using only data within the boundary, which is already known in the simulation.
Then, supposing the time evolution of the boundary is known, and because the
information originating inside the simulation and within the boundary is
known, it follows that information originating from outside the simulation can
be solved for and used to consistently set the boundary conditions correctly.
This is the fundamental idea behind our data-driven boundary condition (DDBC).

In summary, the problem of data driving can be stated as this: How can the MHD
state vector at time $t$ in the volume be evolved to time $t+\Delta t$, given
a new target state vector on the lower boundary at time $t+\Delta t$ and a
requirement that the evolution is fully consistent with MHD? We stress again
that, due to incomplete or uncertain knowledge of the boundary layer, it may
not be possible to simultaneously satisfy both the estimated boundary
evolution and the MHD equations, and in such cases we side with the MHD
equations. With that in mind, the goal of this paper is to describe both the
mathematical background and numerical implementation needed to evolve an MHD
simulation from $t$ to $t+\Delta t$ in a way that is fully consistent with the
MHD equations in the volume and maximally consistent with the new, prescribed
MHD state vector on the lower boundary at $t+\Delta t$. This solution to the
problem may then be repeated $N$ times to drive the MHD system from $t$ to
$t+N\,\Delta t$.

The structure of the rest of the paper is as follows. In §2 we recast the MHD
equations into characteristic form and discuss how this describes the flow of
information, or equivalently, the propagation of specific disturbances,
through the system. In §3 we give an overview of using the characteristics to
set boundary conditions and show how the time evolution of any boundary can be
decomposed into the effect of information that passes through that boundary
from each side. In §4 we show how a time series of observations of a boundary
can be used to solve for the information propagating in one of the normal
directions, thus formally solving the problem of correctly setting boundary
conditions in the continuous case. In §5 we describe the representation of
characteristic based MHD in a finite-volume numerical scheme, CHAR, and the
implementation of data-driven boundary conditions in that code, which we call
CHAR-DDBC. In §6 we describe our validation suite, which centers on a “ground-
truth” simulation of an expanding spheromak implemented the Lagrangian Remap
3D MHD code (LaRe3D: Arber et al., 2001). We emphasize here that part of our
validation involves coupling the CHAR-DDBC code to the LaRe3D code, which have
unrelated numerical schemes; CHAR-DDBC is general-purpose and can be coupled
to any MHD code. We describe the results of the validation simulations in §7,
provide some additional discussion in §8, and summarize in §9. Several
appendices follow that serve to keep the present work self-contained.

As a final note in this Introduction, we stress that our characteristic method
can be used to formulate arbitrary boundary conditions in addition to data-
driven boundary conditions. This will be explored in subsequent papers in this
series. In particular, a problem that is intricately related to the present
one is how to set up boundary conditions that allow information to freely
leave a simulation volume, a so-called (and perhaps unfortunately named)
_nonreflecting boundary condition_ (NRBC). For that topic, we refer the reader
to Kee et al. (2023), hereafter called Paper II. Put succinctly, if the
present Paper I is about correctly getting information into a simulation
volume, Paper II is about correctly getting information out of the simulation
volume. As it turns out, the former stands on much more solid ground than the
latter.

##  2 Rewriting MHD in characteristic form

The key to formulating data-driven boundary conditions consistent with the
internal state of an MHD simulation is to transform the traditional Eulerian
description of MHD into the _characteristic form_ which directly describes the
flow of information through the plasma as specific types of propagating modes.
The characteristic form is particularly useful for writing and analyzing
boundary conditions for three reasons:

  1. 1.

By identifying each _type_ of information (each mode), the characteristics
encode precisely how the plasma at a given location responds to each mode
propagating through that location.

  2. 2.

By identifying the direction and speed at which each mode propagates relative
to the plasma velocity $\bm{\mathbf{v}}$, the characteristic form
distinguishes between information that propagates out of the simulation and
information that propagates into the simulation, through the boundary, from
the external universe.

  3. 3.

By expressing the boundary conditions in terms of the characteristics, they
are intrinsically consistent with the MHD equations themselves.
This final reason has as a corollary that:

    1. 3a.

Any physical assumptions about the dynamics in the external universe inherent
in the boundary conditions must be made explicit in terms of the allowed
characteristic modes.

The first step towards formulating boundary conditions in characteristic form
is to write the MHD equations in that form without any reference to a
boundary, i.e., for an arbitrary location in space333From a data driving
perspective, this initial formulation could be used for internal data driving
and/or data assimilation in general MHD models.. While little of the
information in this section is completely novel—among many other references,
the characteristics of perturbative MHD are detailed in the text by Jeffrey &
Taniuti (1964) and the full theory for numerical simulations was applied to a
1D Roe-type upwind differencing scheme in Brio & Wu (1988)—we include it here
to keep the present text self-contained and to highlight aspects that are
important for developing data-driven boundary conditions for 3D MHD.

The ideal MHD equations that describe conservation of mass, momentum, and
energy, plus the ideal induction equation, are:

$\displaystyle\frac{\partial\rho}{\partial
t}+\nabla\cdot(\rho\bm{\mathbf{v}})=0$ (1a)
$\displaystyle\frac{\partial\bm{\mathbf{v}}}{\partial
t}+\bm{\mathbf{v}}\cdot\nabla\bm{\mathbf{v}}+\frac{1}{\rho}\nabla
P-\frac{1}{\rho}\bm{\mathbf{j}}\times\bm{\mathbf{B}}-\bm{\mathbf{g}}=\bm{\mathbf{0}}$
(1b) $\displaystyle\frac{\partial\epsilon}{\partial
t}+\bm{\mathbf{v}}\cdot\nabla\epsilon+\frac{P}{\rho}\nabla\cdot\bm{\mathbf{v}}=0$
(1c) $\displaystyle\frac{\partial\bm{\mathbf{B}}}{\partial
t}+(\bm{\mathbf{v}}\cdot\nabla)\bm{\mathbf{B}}-(\bm{\mathbf{B}}\cdot\nabla)\bm{\mathbf{v}}+\bm{\mathbf{B}}(\nabla\cdot\bm{\mathbf{v}})=\bm{\mathbf{0}}.$
(1d) The system is closed using the ideal equation of state relating pressure
$P$ to density and internal energy density, $\displaystyle
P=(\gamma-1)\rho\epsilon.$ (1e)

A general equation of state, with pressure a general function of density and
internal energy $P(\rho,\epsilon)$, could also be used and would be better
suited to a photospheric plasma. However, the ideal equation of state leads to
slightly simplified notation and so we adopt it here. All of the following
work carries over with minor changes for the general case. The primitive
variables are density $\rho$, specific internal energy density $\epsilon$, the
velocity vector $\bm{\mathbf{v}}^{T}=(v_{x},v_{y},v_{z})$, and the magnetic
vector $\bm{\mathbf{B}}^{T}=(B_{x},B_{y},B_{z})$. The electric current density
is $\bm{\mathbf{j}}=\nabla\times\bm{\mathbf{B}}/\mu_{0}$, where $\mu_{0}$ is
the permeability of free space.

The equations are nondimensionalized by choosing normalization constants for
the magnetic field, density, and length: $B=B_{N}B^{*}$,
$\rho=\rho_{N}\rho^{*}$, and $L=L_{N}L^{*}$, where symbols subscripted with
$N$ have units and starred quantities are unitless. The normalizations for all
other quantities are defined in terms of these three, the values of which are
given in Appendix §F.2 where we describe the details of the numerical
simulations and their initial conditions. Gradients are given by
$\nabla=\nabla^{*}/L_{N}$, and for the remaining variables we have
$V_{N}=B_{N}/\sqrt{\mu_{0}\rho_{N}}$, $T_{N}=L_{N}/V_{N}$,
$P_{N}=B_{N}^{2}/\mu_{0}$, $\epsilon_{N}=B_{N}^{2}/\mu_{0}\rho_{N}$,
$g_{N}=L_{N}g_{R}/V_{N}^{2}$ (with $g_{R}$ the value of gravity at the solar
surface in SI units), and $j_{N}=B_{N}/\mu_{0}L_{N}$. The final of these means
that the current density in normalized units has the form
$\bm{\mathbf{j}}^{*}=\nabla^{*}\times\bm{\mathbf{B}}^{*}$ so that the
permeability $\mu_{0}$ drops out of the equations. Thus, the MHD equations in
Equation (1) are the same in both normalized and unnormalized (or starred and
unstarred) form, with the single difference that $\bm{\mathbf{j}}$ is replaced
with $\nabla\times\bm{\mathbf{B}}$ in the momentum equation (1b), without
$\mu_{0}$. From here on we drop reference to the starred variables and all
expressions are taken to be normalized unless stated explicitly.

The total MHD state vector written in primitive variables is
$\bm{\mathbf{U}}^{T}=(\rho,\epsilon,v_{x},v_{y},v_{z},B_{x},B_{y},B_{z})$. The
MHD equations can then be rewritten in matrix form as

$\partial_{t}\bm{\mathbf{U}}+\mathbb{A}_{x}\cdot\partial_{x}\bm{\mathbf{U}}+\mathbb{A}_{y}\cdot\partial_{y}\bm{\mathbf{U}}+\mathbb{A}_{z}\cdot\partial_{z}\bm{\mathbf{U}}+\bm{\mathbf{D}}=\bm{\mathbf{0}},$
(2)

where, for instance, the coefficient matrix $\mathbb{A}_{z}$ is given by

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

The coefficient matrices $\mathbb{A}_{x}$ and $\mathbb{A}_{y}$ are similarly
defined, and their full expressions are given in Appendices B and C. The
vector $\bm{\mathbf{D}}$ may contain any inhomogeneous terms one wishes to
include, such as gravity, resistive terms from a generalized Ohm’s law,
volumetric heating terms, and so on. We have included gravity as an example,
in Equations (1); thus, in this paper, the inhomogeneous term is simply
$\bm{\mathbf{D}}{}=\bm{\mathbf{g}}=-g\widehat{\bm{\mathbf{z}}}$, where $g=1$
is the value of gravity near the solar surface in normalized units.

The matrices $\mathbb{A}_{n}$ are individually diagonalizable (Roe & Balsara,
1996), but not simultaneously so outside of special cases.444Notation: $n$
will index the spatial directions $x,y,z$, and $\sigma$ and $\zeta$ will index
the variables of our system of equations, e.g.
$[\mathbb{A}_{z}]_{\sigma\zeta}$. Thus, the first element of the MHD state
vector is $U_{\sigma}=\rho$ with $\sigma=1$. Scalars are written in standard
font face, vectors in bold face, and matrices in a double-struck face. This is
not a concern for us: we ultimately wish to formulate a boundary condition and
are free to diagonalize the system in the direction normal to that boundary.
From here on, we will mostly focus on $\widehat{\bm{\mathbf{z}}}$ as an
example, but the analysis carries over to all directions by cyclic permutation
of $\\{x,y,z\\}$. Looking ahead towards the application of the formalism to
data driving, the diagonalization will eventually be applied individually to
each face of a computational cell to construct a discrete, finite volume
method (see Appendix D for a discussion of our choice of a finite volume over
a finite difference approach to the characteristics).

The primitive variables (and the more typical conserved variables, for that
matter) do not constitute a good set of variables for describing how
information propagates through the plasma. To find a suitable set, we can
rotate the state-space coordinate system in Equation (2) describing the MHD
equations in the standard way by diagonalizing each of the $\mathbb{A}$
matrices with a similarity transformation. As will be clear shortly, this
transformation essentially separates out the different types of information
that can propagate in a given direction. Using $\widehat{\bm{\mathbf{z}}}$ as
an example, the similarity transformation for $\mathbb{A}_{z}$ is

$\mathbb{A}_{z}=\mathbb{S}_{z}\mathbb{\Lambda}_{z}\mathbb{S}_{z}^{-1}\quad\Longleftrightarrow\quad\mathbb{\Lambda}_{z}=\mathbb{S}_{z}^{-1}\mathbb{A}_{z}\mathbb{S}_{z},$
(4)

where $\mathbb{\Lambda}_{z}$ is the diagonal matrix of eigenvalues of
$\mathbb{A}_{z}$ and the $\mathbb{S}_{z}^{-1}$ and $\mathbb{S}_{z}$ matrices
contain the eigenvectors of $\mathbb{A}_{z}$. The left eigenmatrix
$\mathbb{S}_{z}^{-1}$ has the left eigenvectors
$\bm{\mathbf{l}}_{z,\sigma}^{T}$ as each row while the right eigenmatrix
$\mathbb{S}_{z}$ has the right eigenvectors $\bm{\mathbf{r}}_{z,\sigma}$ as
each column, where the index $\sigma$ runs from 1 to 8, corresponding to the 8
primitive variables (or any transformed variables derived from them). Full
expressions for the left and right eigenvectors in the
$\widehat{\bm{\mathbf{z}}}$ direction are given in Appendix A.

The eigenvalues $\lambda_{z,\sigma}$ and left and right eigenvectors satisfy
the usual eigenequations:
$\bm{\mathbf{l}}_{z,\sigma}^{T}\cdot\mathbb{A}_{z}=\bm{\mathbf{l}}_{z,\sigma}^{T}\lambda_{z,\sigma}$
and
$\mathbb{A}_{z}\cdot\bm{\mathbf{r}}_{z,\sigma}=\lambda_{z,\sigma}\bm{\mathbf{r}}_{z,\sigma}$.
The eigenvectors are biorthonormal so that
$\bm{\mathbf{l}}_{z,\sigma}^{T}\cdot\bm{\mathbf{r}}_{z,\zeta}=\delta_{\sigma\zeta}$,
the Kroneker delta, and
$\mathbb{S}_{z}^{-1}\mathbb{S}_{z}=\mathbb{S}_{z}\mathbb{S}_{z}^{-1}=\mathbb{I}$,
the identity matrix. Complete descriptions of the left ($\mathbb{S}_{z}^{-1}$)
and right ($\mathbb{S}_{z}$) eigenmatrices are given in Appendix A; the
corresponding matrices for the $\widehat{\bm{\mathbf{x}}}$ and
$\widehat{\bm{\mathbf{y}}}$ directions are given in Appendices B and C,
respectively.

The eigenvalues of $\mathbb{A}_{z}$ are found by solving the standard equation

$\displaystyle\textsf{Det}(\mathbb{A}_{z}-\lambda\mathbb{I})=0,$ (5) or, with
the determinant expanded out,
$\displaystyle\Bigl{[}(v_{z}-\lambda)^{2}\Bigr{]}\Bigl{[}b_{z}^{2}-(v_{z}-\lambda)^{2}\Bigr{]}\Bigl{[}(v_{z}-\lambda)^{4}-(a^{2}+b^{2})(v_{z}-\lambda)^{2}+b_{z}^{2}a^{2}\Bigr{]}=0.$
(6)

This equation has eight solutions555The eight solutions are unique when
$\mathbb{A}_{z}$ is full rank, but in certain limits the eigenvalues become
degenerate. For example, when the magnetic field is reduced to zero the 8 MHD
solutions collapse onto the 5 hydrodynamic solutions. See Roe & Balsara (1996)
for a full discussion of the eigenstructure of MHD equations in various
limits.:

$\displaystyle\lambda_{1}=v_{z},$ (7a) $\displaystyle\lambda_{2}=v_{z},$ (7b)
$\displaystyle\lambda_{3}=v_{z}-c_{a},$ (7c)
$\displaystyle\lambda_{4}=v_{z}+c_{a},$ (7d)
$\displaystyle\lambda_{5}=v_{z}-c_{s},$ (7e)
$\displaystyle\lambda_{6}=v_{z}+c_{s},$ (7f)
$\displaystyle\lambda_{7}=v_{z}-c_{f},$ (7g)
$\displaystyle\lambda_{8}=v_{z}+c_{f},$ (7h) where $\displaystyle c_{a}=\lvert
b_{z}\rvert,$ (7i) $\displaystyle
c_{f}^{2}=\frac{1}{2}(a^{2}+b^{2})+\frac{1}{2}\sqrt{(a^{2}+b^{2})^{2}-4a^{2}b_{z}^{2}},$
(7j) $\displaystyle
c_{s}^{2}=\frac{1}{2}(a^{2}+b^{2})-\frac{1}{2}\sqrt{(a^{2}+b^{2})^{2}-4a^{2}b_{z}^{2}},$
(7k) $\displaystyle a^{2}=\gamma(\gamma-1)\epsilon,\quad
b_{n}=B_{n}/\sqrt{\rho}\text{ for $n\in(x,y,z)$},\quad b^{2}=\sum b_{n}^{2}.$
(7l)

The eigenvalues $\lambda_{\sigma}$ in Equations (7) constitute the components
of the diagonal eigenvalue matrix in Equation (4),
$\mathbb{\Lambda}_{z}=\textsf{diag}(\bm{\mathbf{\lambda}}_{z})$. The
eigenvalues have units of velocity and are expressed in terms of the bulk
velocity in the diagonalization direction $v_{z}$, the sound speed $a$, the
total and projected Alfvén speeds $b$ and $c_{a}$ (with projected components
$b_{n}$ in each direction $n\in(x,y,z)$), the fast magnetosonic speed $c_{f}$,
and the slow magnetosonic speed $c_{s}$. Note that the Alfvén, fast, and slow
speeds are modified from the standard 1D expressions by incorporating a
projection onto the diagonalization direction: only $b_{z}$ appears in the
definitions of the projected Alfvén speed $c_{a}$, or in the $4a^{2}b_{z}^{2}$
term under the radical in Equations (7j) and (7k).

The eigenvalues themselves are propagation speeds, i.e., the speed at which
each possible type of information described by the MHD equations flows through
the plasma in the diagonalization direction: $\lambda_{1,2}$ are advection in
$\widehat{\bm{\mathbf{z}}}{}$, $\lambda_{3,4}$ the Alfvén speeds,
$\lambda_{5,6}$ the slow magnetosonic speeds, and $\lambda_{7,8}$ the fast
magnetosonic speeds, all calculated relative to the background velocity in the
diagonalization direction, in this example $v_{z}$. These characteristic
speeds are equivalent to those reported by, e.g., Roe & Balsara (1996),
although ours has one extra eigenvalue $v_{z}$; this is because Roe & Balsara
(1996) solved the 1-D system with $\nabla\cdot\bm{\mathbf{B}}=0$ enforced by
construction, while we choose to enforce the solenoidal constraint with an
additional eigenmode; we discuss this further at the end of this section,
after Equations (18).

The eigendecomposition of the $\mathbb{A}_{z}$ matrix described in Equation
(4) can also be carried out for $\mathbb{A}_{x}$ and $\mathbb{A}_{y}$, as
given in in Appendices B and C, respectively. Doing so leads to essentially
the same set of expressions for the eigenvalues and eigenvectors in each case,
up to cyclic permutations of the variables in $x$, $y$, and $z$. The fact that
the eigenvalues and vectors for each direction only match up to a cyclic
permutation of the projection direction is the reason that the MHD equations
(2) are not diagonalizable in all directions simultaneously: the eigensystem
of one matrix does not match those of the others. Nonetheless, the
eigendecomposition can be carried out in each direction separately.

Substituting the diagonalized forms for all the coefficient matrices
$\mathbb{A}_{x},\ \mathbb{A}_{y},\text{ and }\mathbb{A}_{z}$ into Equation
(2), we rewrite the MHD equations as

$\displaystyle\partial_{t}\bm{\mathbf{U}}+\mathbb{S}_{x}\mathbb{\Lambda}_{x}\mathbb{S}_{x}^{-1}\cdot\partial_{x}\bm{\mathbf{U}}+\mathbb{S}_{y}\mathbb{\Lambda}_{y}\mathbb{S}_{y}^{-1}\cdot\partial_{y}\bm{\mathbf{U}}+\mathbb{S}_{z}\mathbb{\Lambda}_{z}\mathbb{S}_{z}^{-1}\cdot\partial_{z}\bm{\mathbf{U}}+\bm{\mathbf{D}}{}=\bm{\mathbf{0}}.$
(8)

Grouping the derivative term in each direction
$(\partial_{n}\bm{\mathbf{U}}{})$ with its respective eigenvalue matrix
$(\mathbb{\Lambda}_{n})$ and left eigenvector matrix $(\mathbb{S}^{-1}_{n})$
defines the characteristic derivative vector $\bm{\mathbf{L}}_{n}$ in each
direction:

$\displaystyle\bm{\mathbf{L}}_{x}=\mathbb{\Lambda}_{x}{}\mathbb{S}^{-1}_{x}\cdot\partial_{x}\bm{\mathbf{U}},$
(9a)
$\displaystyle\bm{\mathbf{L}}_{y}=\mathbb{\Lambda}_{y}{}\mathbb{S}^{-1}_{y}\cdot\partial_{y}\bm{\mathbf{U}},$
(9b)
$\displaystyle\bm{\mathbf{L}}_{z}=\mathbb{\Lambda}_{z}{}\mathbb{S}^{-1}_{z}\cdot\partial_{z}\bm{\mathbf{U}}.$
(9c)

Substituting these $\bm{\mathbf{L}}_{n}$ into Equation (8), the MHD equations
become

$\displaystyle\partial_{t}\bm{\mathbf{U}}+\mathbb{S}_{x}\cdot\bm{\mathbf{L}}_{x}+\mathbb{S}_{y}\cdot\bm{\mathbf{L}}_{y}+\mathbb{S}_{z}\cdot\bm{\mathbf{L}}_{z}+\bm{\mathbf{D}}=\bm{\mathbf{0}}.$
(10)

Equation (10) is equivalent to the full system of MHD equations (1) or (2),
but rewritten in characteristic form. Note that we have not yet chosen a
boundary direction; this is still a description of a plasma whose properties
are known at all spatial locations.

The physical meaning of each element in Equation (10) can be understood even
before the full expressions are worked out. Each component ($\sigma$) of one
of the characteristic derivatives, $L_{n,\sigma}$, corresponds to a specific
combination of information (derivatives of primitive variables) that
propagates with a given speed and direction relative to the diagonalization
direction, $\widehat{\bm{\mathbf{n}}}$; for example, a fast mode disturbance
propagating in $\widehat{\bm{\mathbf{z}}}{}$. To see this, observe that:

  1. (i)

the derivatives in Equation (9) are changes in the state of the plasma from
point to point in space,

  2. (ii)

a single row $\sigma$ of a left eigenvalue matrix $\mathbb{S}^{-1}_{n}$, which
is the left eigenvector $\bm{\mathbf{l}}^{T}_{\sigma}$, couples together a set
of those changes through the dot product,
and

  3. (iii)

the coupled set of changes is multiplied by a single eigenvector speed
$\lambda_{n,\sigma}$ (recall the matrix $\mathbb{\Lambda}{}_{n}$ is diagonal).

The result of the calculation is an object which represents a specific set of
coupled changes in the primitive variables that propagate at a particular
speed in a particular direction, i.e., it represents a specific type of
information propagating through the plasma.

We emphasize that one component $\sigma$ of a characteristic derivative vector
in one direction (e.g., $L_{n,\sigma}$) is a _scalar_ quantity that encodes
how _all_ components of the MHD state vector are simultaneously modified due
to the $\sigma$ mode propagating in the $\widehat{\bm{\mathbf{n}}}$ direction.
That coupled change between all the primitive variables is the embodiment of
the dynamical constraints imposed by the MHD equations: out of all arbitrary
variations to the MHD state vector $\bm{\mathbf{U}}{}$ one might consider,
only the eight described by the eight components of the characteristic
derivative $\bm{\mathbf{L}}_{n}$ in the $\widehat{\bm{\mathbf{n}}}$ direction
are allowed by, and consistent with, the MHD equations.

A superposition of multiple modes is, of course, allowed, and the right
eigenvectors $\bm{\mathbf{r}}_{\sigma}$ describe the relative contribution of
each mode to the time update of each primitive variable. To see this, look
again to Equation (10). The right eigenmatrix $\mathbb{S}_{n}$ multiplies its
associated characteristic derivative vector. Each column of the right
eigenmatrix $\mathbb{S}_{n}$ is a right eigenvector
$\bm{\mathbf{r}}_{\sigma}$, and, carrying though the multiplication between
$\mathbb{S}_{n}$ and $\bm{\mathbf{L}}_{n}$, the elements of a right
eigenvector give the relative contribution of each mode to the time-update of
the primitive variables.

Summarizing the above few paragraphs, the specific combination of weights of
primitive variable derivatives prescribed by a left eigenvector
$\bm{\mathbf{l}}_{n,\sigma}$ defines a mode, the scalar value of the
characteristic derivative $L_{n,\sigma}$ defines the strength of the mode, and
each right eigenvector $\bm{\mathbf{r}}_{n,\zeta}$ describes the response of a
primitive variable to each mode666Throughout, we will often refer to “the
$L_{n,\sigma}$ mode” or “$L_{n,\sigma}$ characteristic” when, strictly
speaking, the characteristic derivative $L_{n,\sigma}$ is the value of the
weighted sum of primitive variable derivatives associated with the propagation
of the $\sigma$ mode evaluated at a certain spatial location. In practice, the
modes and the derivatives associated with them are matched 1-to-1, so we will
typically use the less verbose terminology.. The interpretation of Equation
(10) that follows from this discussion is that the temporal update to the MHD
state vector $\bm{\mathbf{U}}{}$ at some location $\bm{\mathbf{x}}$ is due to
the combined effect of all information propagating through that location from
each type of eigenmode in each direction. The eigendecomposition performed at
Equation (4) separated the various types of information that can propagate
through the plasma at $\bm{\mathbf{x}}$ from the effect that each type of
propagating information has on the plasma at location $\bm{\mathbf{x}}$. This
decomposition will eventually allow us to identify the information propagating
into a subvolume of plasma from the external universe, and thereby enable us
to enforce only allowed changes to the plasma at a simulation boundary. Note
that, if the supplied information is not exactly self-consistent then the
boundary condition will not satisfy the Cauchy conditions and hence will be
over-determined (Morse & Feshbach, 1953, Chapter 6). In that case, we must
either depart from the known physics of our universe or search for the closest
boundary state that _does_ obey MHD; our solution for the latter case is in
§5.

We turn now to the details of analyzing how information propagates in a
specific direction with the eventual goal of determining the way in which
information can cross a boundary and what effect it has on the plasma at the
boundary. To make contact with photosphere-to-corona simulations, we take the
$\widehat{\bm{\mathbf{z}}}$ direction as an example and group the transverse
terms (involving $\partial_{x},\partial_{y}$) and inhomogeneous terms
($\bm{\mathbf{D}}$) in Equation (10) into a single vector $\bm{\mathbf{C}}$:

$\displaystyle\bm{\mathbf{C}}=\mathbb{A}_{x}\cdot\partial_{x}\bm{\mathbf{U}}+\mathbb{A}_{y}\cdot\partial_{y}\bm{\mathbf{U}}+\bm{\mathbf{D}}=\mathbb{S}_{x}\cdot\bm{\mathbf{L}}_{x}+\mathbb{S}_{y}\cdot\bm{\mathbf{L}}_{y}+\bm{\mathbf{D}}.$
(11)

Full expressions for the vector $\bm{\mathbf{C}}$ are given in Appendix D.
Note that $\bm{\mathbf{C}}$ lives in the space of the MHD primitive variables,
i.e., each element $C_{\sigma}$ contributes to the update of one of the
primitive variables. Substituting $\bm{\mathbf{C}}$ into Equation (10), we
finally arrive at

$\partial_{t}\bm{\mathbf{U}}{}+\mathbb{S}_{z}\cdot\bm{\mathbf{L}}_{z}+\bm{\mathbf{C}}=\bm{\mathbf{0}}.$
(12)

Equation (12) is the primary equation we will focus on for the remainder of
this paper. It is still just an equivalent form of the MHD equations but now
written in a way that is useful for describing boundary conditions in the
$\widehat{\bm{\mathbf{z}}}$ direction. In particular, it consolidates all our
knowledge about the _flow of information in the $\widehat{\bm{\mathbf{z}}}$
direction_ into the vector $\bm{\mathbf{L}}_{z}$.

Consider an arbitrary plane at constant $z=z_{b}$, which we suggestively label
$\bm{\mathbf{x}}_{b}=(x,y,z_{b})$. If the state of the system is known
everywhere in the plane $\bm{\mathbf{x}}_{b}$ at some time $t$ then every
component of the combined transverse and inhomogeneous term
$\bm{\mathbf{C}}(t;\bm{\mathbf{x}}_{b})$ can be calculated in the plane at
that time. On the other hand, the characteristic derivative vector
$\bm{\mathbf{L}}_{z}$ defined in Equation (9c) contains derivatives in
direction normal to the plane and therefore inherently requires information
from outside of the plane in order to fully specify it. At the moment, the
information propagating in each direction through the plane is still bundled
together in the vector $\bm{\mathbf{L}}_{z}$, but these two pieces will soon
be separated out. If the arbitrary plane is actually a boundary, so that the
MHD state vector is known everywhere on one side (the interior) but not the
other (the exterior), then the components of $\bm{\mathbf{L}}_{z}$ that
describe inward propagating modes will represent independent information
coming into the system.

Having discussed the meaning of each element of the MHD equations when written
in characteristic form, either as Equation (10) or (12), we turn now to the
mathematical details of those elements. The following discussion still applies
to the continuous case and makes no assumptions about boundaries or numerical
discretization of the system. We will return to the application of the
characteristic derivatives to boundary conditions (and data-driven boundary
conditions in particular) in Section 3.

Each component of the characteristic derivative777To reduce notational
clutter, for the rest of this section elements of an eigensystem that do not
explicitly refer to a direction will refer to the $\widehat{\bm{\mathbf{z}}}$
direction. When either the $\widehat{\bm{\mathbf{x}}}$ or
$\widehat{\bm{\mathbf{y}}}$ directions are needed, they will be explicitly
labeled with a subscript. For example, $L_{1}$ is taken to mean $L_{z,1}$, and
the fast mode speed labeled in the most complete fashion is $c_{z,f}$.
$\bm{\mathbf{L}}$, $L_{\sigma}$, is a weighted sum of the normal derivatives
of the primitive variables multiplied by a characteristic speed
$\lambda_{\sigma}$. The components of the left eigenvector are the weights for
each primitive variable derivative for the $\sigma$ mode:

$\displaystyle
L_{\sigma}=\lambda_{\sigma}(\bm{\mathbf{l}}_{\sigma}^{T}\cdot\partial_{z}\bm{\mathbf{U}})=\lambda_{\sigma}\sum_{\zeta}S_{\sigma\zeta}^{-1}\partial_{z}U_{\zeta}.$
(13)

The eight components of $\bm{\mathbf{L}}$ (from which the left eigenvectors
can be read off, cf. Appendix A) are

divB $\displaystyle L_{1}$
$\displaystyle=v_{z}\biggl{[}B_{z}^{\prime}\biggr{]}$ (14a) entropy
$\displaystyle L_{2}$
$\displaystyle=v_{z}\biggl{[}\frac{1-\gamma}{\gamma\rho}\rho^{\prime}+\frac{1}{\gamma\epsilon}\epsilon^{\prime}\biggr{]}$
(14b) reverse Alfvén $\displaystyle L_{3}$
$\displaystyle=\frac{v_{z}-c_{a}}{2}\biggl{[}-\beta_{y}v_{x}^{\prime}+\beta_{x}v_{y}^{\prime}-\frac{\beta_{y}s_{z}}{\sqrt{\rho}}B_{x}^{\prime}+\frac{\beta_{x}s_{z}}{\sqrt{\rho}}B_{y}^{\prime}\biggr{]}$
(14c) forward Alfvén $\displaystyle L_{4}$
$\displaystyle=\frac{v_{z}+c_{a}}{2}\biggl{[}\beta_{y}v_{x}^{\prime}-\beta_{x}v_{y}^{\prime}-\frac{\beta_{y}s_{z}}{\sqrt{\rho}}B_{x}^{\prime}+\frac{\beta_{x}s_{z}}{\sqrt{\rho}}B_{y}^{\prime}\biggr{]}$
(14d) reverse slow $\displaystyle L_{5}$
$\displaystyle=\frac{v_{z}-c_{s}}{2}\biggl{[}\frac{\alpha_{s}}{\gamma\rho}\rho^{\prime}+\frac{\alpha_{s}}{\gamma\epsilon}\epsilon^{\prime}-\frac{\beta_{x}\alpha_{f}c_{f}}{a^{2}}s_{z}v_{x}^{\prime}-\frac{\beta_{y}\alpha_{f}c_{f}}{a^{2}}s_{z}v_{y}^{\prime}-\frac{\alpha_{s}c_{s}}{a^{2}}v_{z}^{\prime}-\frac{\beta_{x}\alpha_{f}}{\sqrt{\rho}a}B_{x}^{\prime}-\frac{\beta_{y}\alpha_{f}}{\sqrt{\rho}a}B_{y}^{\prime}\biggr{]}$
(14e) forward slow $\displaystyle L_{6}$
$\displaystyle=\frac{v_{z}+c_{s}}{2}\biggl{[}\frac{\alpha_{s}}{\gamma\rho}\rho^{\prime}+\frac{\alpha_{s}}{\gamma\epsilon}\epsilon^{\prime}+\frac{\beta_{x}\alpha_{f}c_{f}}{a^{2}}s_{z}v_{x}^{\prime}+\frac{\beta_{y}\alpha_{f}c_{f}}{a^{2}}s_{z}v_{y}^{\prime}+\frac{\alpha_{s}c_{s}}{a^{2}}v_{z}^{\prime}-\frac{\beta_{x}\alpha_{f}}{\sqrt{\rho}a}B_{x}^{\prime}-\frac{\beta_{y}\alpha_{f}}{\sqrt{\rho}a}B_{y}^{\prime}\biggr{]}$
(14f) reverse fast $\displaystyle L_{7}$
$\displaystyle=\frac{v_{z}-c_{f}}{2}\biggl{[}\frac{\alpha_{f}}{\gamma\rho}\rho^{\prime}+\frac{\alpha_{f}}{\gamma\epsilon}\epsilon^{\prime}+\frac{\beta_{x}\alpha_{s}c_{s}}{a^{2}}s_{z}v_{x}^{\prime}+\frac{\beta_{y}\alpha_{s}c_{s}}{a^{2}}s_{z}v_{y}^{\prime}-\frac{\alpha_{f}c_{f}}{a^{2}}v_{z}^{\prime}+\frac{\beta_{x}\alpha_{s}}{\sqrt{\rho}a}B_{x}^{\prime}+\frac{\beta_{y}\alpha_{s}}{\sqrt{\rho}a}B_{y}^{\prime}\biggr{]}$
(14g) forward fast $\displaystyle L_{8}$
$\displaystyle=\frac{v_{z}+c_{f}}{2}\biggl{[}\frac{\alpha_{f}}{\gamma\rho}\rho^{\prime}+\frac{\alpha_{f}}{\gamma\epsilon}\epsilon^{\prime}-\frac{\beta_{x}\alpha_{s}c_{s}}{a^{2}}s_{z}v_{x}^{\prime}-\frac{\beta_{y}\alpha_{s}c_{s}}{a^{2}}s_{z}v_{y}^{\prime}+\frac{\alpha_{f}c_{f}}{a^{2}}v_{z}^{\prime}+\frac{\beta_{x}\alpha_{s}}{\sqrt{\rho}a}B_{x}^{\prime}+\frac{\beta_{y}\alpha_{s}}{\sqrt{\rho}a}B_{y}^{\prime}\biggr{]},$
(14h)

where the prime denotes differentiation in the $z$-direction, i.e.,
$\partial_{z}q=q^{\prime}$, for $q$ an arbitrary variable. Each subequation is
labeled by the type of mode (or disturbance) it represents. The
$\textsl{div}B$ and entropy modes each stand alone, while the Alfvén-, slow-,
and fast-mode families each include a pair of related modes. The related modes
in a given family have the same coefficients on each primitive variable
derivative, but unique signs.

The above equations make use of the auxiliary variables

$\displaystyle s_{z}$ $\displaystyle=\operatorname{sign}(b_{z})$ (15a)
$\displaystyle b_{\perp}^{2}$ $\displaystyle=b_{x}^{2}+b_{y}^{2}$ (15b)
$\displaystyle\beta_{x}$ $\displaystyle=\frac{b_{x}}{b_{\perp}}$ (15c)
$\displaystyle\beta_{y}$ $\displaystyle=\frac{b_{y}}{b_{\perp}}$ (15d)
$\displaystyle\alpha_{f}^{2}$
$\displaystyle=\frac{a^{2}-c_{s}^{2}}{c_{f}^{2}-c_{s}^{2}}$ (15e)
$\displaystyle\alpha_{s}^{2}$
$\displaystyle=\frac{c_{f}^{2}-a^{2}}{c_{f}^{2}-c_{s}^{2}},$ (15f)

where $\alpha_{f}^{2}+\alpha_{s}^{2}=1$. The prefactor in front of the
bracketed term on the right hand side of each of Equations (14) is
proportional to the eigenvalue $\lambda_{\sigma}$ of the eigenmode and has
units of velocity (see Eq. (7)).

Figure 1: Diagram summarizing the relationship between the three types of
propagation, defined relative to: (i) the bulk velocity, (ii) the coordinate
axes, and (iii) the boundary of a subvolume of the plasma. Two examples are
shown for the velocity $v_{z}$ in the $\widehat{\bm{\mathbf{z}}}$ direction.
The modes are labeled by their characteristic derivatives $L_{\sigma}$. (a)
$v_{z}=0$, for which the forward-propagating modes (relative to the bulk
velocity $v_{z}$) coincide with the positive propagating modes (relative to
unit direction $\widehat{\bm{\mathbf{z}}}$), and the reverse-propagating modes
coincide with the negative propagating modes. In example (b), the bulk
velocity is directed downward with a magnitude exceeding all but the fast mode
speed $c_{f}$. The forward and reverse propagating modes are the same as in
case (a), but now the identification of the (grid-) positive and negative
modes has changed, due to the different bulk velocity. Finally, the
identification of _inward_ versus _outward_ propagating modes, defined
relative the boundary-normal directions, are also indicated.

There are three ways of defining the direction of propagation for each mode
that need to be considered, as depicted diagrammatically in Figure 1. First,
the sign of the eigenvalue of the $\sigma$ mode,
$\textsf{sign}(\lambda_{\sigma})$, determines the direction of propagation for
that mode relative to the diagonalization direction (i.e., the coordinate
axes), with _positive_ (+) eigenvalues propagating in
$+\widehat{\bm{\mathbf{z}}}$ and _negative_ $(-)$ eigenvalues in
$-\widehat{\bm{\mathbf{z}}}$. Second, the $+$ or $-$ sign in the expression
for each eigenvalue in Equations (7a)-(7h) determines if that mode propagates
in the _forward_ or _reverse_ direction relative to the bulk velocity $v_{z}$.
Eigenmodes 1 and 2 are purely advective (with $\lambda=v_{z}$), and hence do
not propagate relative to the bulk velocity. The final way of defining the
propagation direction is relative to the boundary of a subvolume of interest.
For instance, at the bottom of a simulation, the positive modes propagate
_inward_ ($\mathcal{I}$) towards the simulation interior and the negative
modes propagate _outward_ ($\mathcal{O}$) from the simulation towards the
exterior universe. Thus, at the top boundary of a simulation, the roles of the
positive (negative) propagating modes are instead identified as outward
(inward) propagating.

Figure 1 presents two examples that help to clarify the three types of
propagation that we will repeatedly refer to in this paper. Panel (a) presents
the case of a static plasma with $v_{z}=0$, and Panel (b) the case where the
bulk velocity is downward with a magnitude in between the fast magnetosonic
speed and the Alfvén speed, that is, $-c_{f}<v_{z}<-c_{a}$. The _forward_ and
_reverse_ propagating modes are labeled at the top of each panel in the
Figure. Because there is a 1-1 correspondence between characteristic
derivatives and the modes that define them, we simply use the names of the
derivatives to refer to the modes. Note that the characteristic modes
identified as forward- or reverse-propagating are always the same regardless
of the normal bulk velocity, so they do not change between Panels (a) and (b).
This is because the characteristics are ultimately defined in a co-moving
frame with the plasma: all the eigenvalues in (7) are Doppler shifted relative
to $v_{z}$.

Next, the terms _positive_ and _negative_ refer to propagation relative to the
coordinate axes, or equivalently, the grid axes in a simulation, and the
number and identification of these modes changes depending on the bulk plasma
velocity. For the static case $v_{z}=0$ in Panel (a), the positive and
negative modes precisely coincide with the forward and reverse modes. In this
case, the purely advective modes $L_{1}$ and $L_{2}$ are non-propagating. In
contrast, for the example of the large magnitude downward velocity in Panel
(b), only the $\sigma=8$ mode can propagate in the positive direction while
all the other modes are negatively propagating relative to the normal
coordinate axis. In the extreme limit where the bulk velocity exceeds the
fastest characteristic speed, $\lvert v_{z}\rvert>c_{f}$, all modes will
propagate in a single coordinate direction, e.g.,
$\pm\widehat{\bm{\mathbf{z}}}$.

Finally, the identification of the _incoming_ and _outgoing_ propagating modes
in each panel is simply a renaming of the positive and negative propagating
modes with the direction of a boundary taken into account. However, this is
still necessary nomenclature for properly describing boundary conditions,
because the outward propagating modes represent calculable information
communicated from the known plasma properties to the external universe while
the inward propagating modes represent the injection of new, previously
unknown information into the region of space where the state of the plasma is
known.

The type of information propagated through the plasma by each mode is labeled
in each of Equations (14). In order, the characteristic derivatives describe:
the conservation of magnetic flux with the fluid flow ($L_{1}$), advection of
entropy with the fluid flow ($L_{2}$), reverse and forward propagation of
Alfvénic disturbances (relative to the fluid flow; $L_{3}$ and $L_{4}$),
reverse and forward propagation of slow disturbances (relative to the fluid
flow; $L_{4}$ and $L_{6}$), and reverse and forward propagation of fast
disturbances (relative to the fluid flow; $L_{7}$ and $L_{8}$). Note that,
while each of these modes is related to the waves in linearized MHD, they are
much more general, but also local (waves are less general but non-local).
Because the characteristic description is exactly equivalent to other forms
for full MHD, any static configuration or dynamic process allowed by the full,
nonlinear MHD equations is expressible in terms of the characteristic modes,
although carrying out such a description may be difficult.

As expressed in Equations (9) and (10), the characteristic derivatives are
defined in terms of left eigenvectors and then multiplied by the right
eigenvectors in order to determine the time update to the primitive variables
in the vector $\bm{\mathbf{U}}$ . Full expressions for the left and right
eigenvectors are given in Appendix A. These are equivalent to those given by
Brio & Wu (1988) and Roe & Balsara (1996), both of whom used the same set of
auxiliary variables (15) but different primitive variables. These auxiliary
variables basically define the normalization of the eigenvectors. The
eigenvectors of any system can be arbitrarily normalized, and the
normalization given above in terms of the auxiliary variables is preferable
for numerical computation because, as Roe & Balsara (1996) show in their §2,
the eigenvectors of an MHD coefficient matrix $\mathbb{A}$ remain well behaved
in various limits when written in terms of the auxiliary variables (15). In
particular, the apparent divergence in (15c)-(15d) as $b_{\perp}\rightarrow 0$
remains bounded, and setting $\beta_{x}=\beta_{y}\rightarrow 1/\sqrt{2}$ in
that case preserves the orthonormality of the eigenvectors.

The characteristic derivatives in a direction (e.g., $\bm{\mathbf{L}}_{z}$)
have a complicated relationship to the spatial derivatives of the primitive
variables in that same direction ($\partial_{z}\bm{\mathbf{U}}$). The
coefficients of the normal derivatives in Equations (14) do define a
transformation between the two through the transformation matrix
$\mathbb{E}$888Appendix E provides the inverse matrix $\mathbb{E}^{-1}$ as
well as the derivatives of the primitive variables in terms of the
characteristic derivatives.

$\displaystyle\bm{\mathbf{L}}_{z}=\mathbb{\Lambda}_{z}\mathbb{S}^{-1}_{z}\partial_{z}\bm{\mathbf{U}}\equiv\mathbb{E}\partial_{z}\bm{\mathbf{U}},$
(16) which has an inverse transformation
$\displaystyle\partial_{z}\bm{\mathbf{U}}=\mathbb{E}^{-1}\bm{\mathbf{L}}_{z}=\mathbb{S}_{z}\mathbb{\Lambda}_{z}^{-1}\bm{\mathbf{L}}_{z},$
(17)

where $\mathbb{\Lambda}^{-1}_{z,\sigma\sigma}=\frac{1}{\lambda_{\sigma}}$ for
$\lambda_{\sigma}\neq 0$ and 0 otherwise (i.e., it is the pseudoinverse). For
any valid 3D MHD field the transformation matrix $\mathbb{E}$ should be full
rank999Alternatively, the transformation matrix may gracefully collapse onto a
subspace of the equations when any eigenvalues and eigenvectors overlap. This
happens, for instance, in the hydrodynamic limit of
$\bm{\mathbf{B}}\rightarrow\bm{\mathbf{0}}$. Again, see Roe & Balsara (1996)
for a full description of these limits.; hence, if either the characteristic
derivatives or primitive variable derivatives are given then the other set can
be calculated from the first. The fact that the transformation matrix
$\mathbb{E}$ is not a diagonal matrix makes clear that varying a single
component of a characteristic derivative vector produces corresponding
variations in multiple primitive variables, and vice versa. This
interdependence is why the characteristics are useful for setting boundary
conditions: each mode encodes a coupled, physically admissible variation of
the entire MHD state vector; or in the other direction, out of all possible
variations of the MHD state vector, those expressible in terms of
characteristic derivatives are the physically allowed variations. However,
because the characteristics represent the directional flow of information, it
is not immediately obvious how one should set up a discrete version of the
transformation. The characteristic derivatives should presumably be cast in
terms of one-sided derivatives (i.e., upwind and downwind), while this is not
necessary for derivatives of the primitive variables. We will return to that
discussion in §5.

Carrying out the matrix multiplication in Equation (12) between the right
eigenvector matrix $\mathbb{S}_{z}$ and the characteristic derivatives
$\bm{\mathbf{L}}_{z}$ allows the MHD equations to be expressed directly in
terms of the characteristic derivatives $L_{\sigma}$:

$\displaystyle\partial_{t}\rho-\rho[L_{2}]+\alpha_{s}\rho[L_{5}+L_{6}]+\alpha_{f}\rho[L_{7}+L_{8}]+C_{\rho}=0$
(18a)
$\displaystyle\partial_{t}\epsilon+\epsilon[L_{2}]+\alpha_{s}\frac{a^{2}}{\gamma}[L_{5}+L_{6}]+\alpha_{f}\frac{a^{2}}{\gamma}[L_{7}+L_{8}]+C_{\epsilon}=0$
(18b)
$\displaystyle\partial_{t}v_{x}-\beta_{y}[L_{3}-L_{4}]-\alpha_{f}\beta_{x}c_{f}s_{z}[L_{5}-L_{6}]+\alpha_{s}\beta_{x}c_{s}s_{z}[L_{7}-L_{8}]+C_{v_{x}}=0$
(18c)
$\displaystyle\partial_{t}v_{y}+\beta_{x}[L_{3}-L_{4}]-\alpha_{f}\beta_{y}c_{f}s_{z}[L_{5}-L_{6}]+\alpha_{s}\beta_{y}c_{s}s_{z}[L_{7}-L_{8}]+C_{v_{y}}=0$
(18d)
$\displaystyle\partial_{t}v_{z}-\alpha_{s}c_{s}[L_{5}-L_{6}]-\alpha_{f}c_{f}[L_{7}-L_{8}]+C_{v_{z}}=0$
(18e)
$\displaystyle\partial_{t}B_{x}-\beta_{y}\sqrt{\rho}s_{z}[L_{3}+L_{4}]-\alpha_{f}\beta_{x}\sqrt{\rho}a[L_{5}+L_{6}]+\alpha_{s}\beta_{x}\sqrt{\rho}a[L_{7}+L_{8}]+C_{B_{x}}=0$
(18f)
$\displaystyle\partial_{t}B_{y}+\beta_{x}\sqrt{\rho}s_{z}[L_{3}+L_{4}]-\alpha_{f}\beta_{y}\sqrt{\rho}a[L_{5}+L_{6}]+\alpha_{s}\beta_{y}\sqrt{\rho}a[L_{7}+L_{8}]+C_{B_{y}}=0$
(18g) $\displaystyle\partial_{t}B_{z}+[L_{1}]+C_{B_{z}}=0.$ (18h)

This expanded form of the MHD equations is exactly equivalent to any set of
Equations (1), (2), (8), (10), or (12). The subscripted $C$ in Equations (18)
are the components of the vector of transverse and inhomogeneous terms
$\bm{\mathbf{C}}$, defined in Equation (11). Full expressions for each
component are given in Appendix D (the subscript notation for $C$ emphasizes
that each component contributes to the time update of a single primitive
variable).

The Equation (18) form of the MHD equations makes it abundantly clear how the
time variation of each primitive variable depends on the characteristic
derivatives and the flow of information they imply. For instance, the
continuity equation (18a) shows that the density at some location will change
because fluid with a different entropy is advected to that location (via
$L_{2}$), or because slow or fast disturbances propagate to that location (via
any of $L_{5,6,7,8}$). Similar interpretations hold for the other primitive
variables. From the Eulerian perspective of a primitive variable at a given
location, it does not matter in which direction any of the modes propagate:
the variable simply responds to all the modes that contribute to its update,
with the strength of the response given by the value of each mode’s
characteristic derivative multiplied by the corresponding component of the
right eigenvector. On the other hand, calculating the value of a
characteristic derivative does depend on the direction of propagation for the
mode it represents: the derivative should be calculated in the upwind
direction so that the changes in the plasma represented by a characteristic
derivative are propagated at the characteristic speed in the correct direction
for the associated mode.

Recall that we defined three flavors of propagation direction when discussing
Figure 1: forward/reverse, positive/negative ($+/-$), and inward/outward
($\mathcal{I}/\mathcal{O}$), which refer to propagation relative to the bulk
velocity, the coordinates (or numerical grid), and a boundary, respectively.
In each of Equations (18) the forward– and reverse–propagating pair of modes
in each family (i.e., the slow, Alfvén, and fast modes relative to the bulk
velocity) are grouped within square brackets101010The
$\nabla\cdot\bm{\mathbf{B}}$ preserving mode, $L_{1}$, and the entropy mode,
$L_{2}$, are purely advective modes and do not come in pairs.. Referring to
the definitions of eigenvalue speeds in (7), depending on the value of the
bulk velocity $v_{z}$, either or both of the modes grouped in a single set of
brackets can propagate in the positive or negative direction (i.e., relative
to the coordinate axes). Then, once a boundary $\bm{\mathbf{x}}_{b}$ is
specified, for instance taking the $z=0$ plane as the $z_{\text{min}}$
boundary of our volume of interest, the positive and negative propagating
modes in turn define the inward and outward propagating modes relative to that
boundary.

Table 1: Characteristic derivatives (second and fourth columns) in the $\widehat{\bm{\mathbf{z}}}$ direction ordered by boundary–normal velocity $v_{z}$ (first column), for the bottom boundary ($z_{\text{min}}$) of some volume. The horizontal line marks $v_{z}=0$, with $v_{z}\leq 0$ above the line and $0\leq v_{z}$ below the line. For reference, the fourth column shows the outgoing characteristic derivatives, which are simply the reverse ordering of the inward characteristics. $P$ refers to the number of inward propagating modes and $Q$ to the number subequations in (18) that the incoming characteristic derivatives participate in. The slow, Alfvén, and fast speeds are all defined to be non-negative and are ordered $0\leq c_{s}\leq c_{a}\leq c_{f}$; see (7). Note that whenever the bulk velocity matches a mode’s propagation velocity the corresponding mode is nonpropagating and, technically, neither incoming nor outgoing. However, in that case the mode has zero amplitude (see Eqns. (14)) and thus no effect on the system. Practically, this means we can include the equals sign on both sides of each division (‘$\lesssim$’ instead of ‘$<$’) without loss of generality.  Velocity | Incoming ($\mathcal{I}$) $L_{\sigma}$ | Dependent Derivatives | Outgoing ($\mathcal{O}$) $L_{\sigma}$ | P | Q
---|---|---|---|---|---
$v_{z}\leq-c_{f}$ | - | - | $L_{7},L_{3},L_{5},{L_{1}},{L_{2}},L_{6},L_{4},L_{8}$ | 0 | 0
$-c_{f}\leq v_{z}\leq-c_{a}$ | $L_{8}$ | $\rho^{\prime},\epsilon^{\prime},v_{x}^{\prime},v_{y}^{\prime},v_{z}^{\prime},B_{x}^{\prime},B_{y}^{\prime}$ | $L_{7},L_{3},L_{5},{L_{1}},{L_{2}},L_{6},L_{4}$ | 1 | 7
$-c_{a}\leq v_{z}\leq-c_{s}$ | $L_{4},L_{8}$ | $\rho^{\prime},\epsilon^{\prime},v_{x}^{\prime},v_{y}^{\prime},v_{z}^{\prime},B_{x}^{\prime},B_{y}^{\prime}$ | $L_{7},L_{3},L_{5},{L_{1}},{L_{2}},L_{6}$ | 2 | 7
$-c_{s}\leq v_{z}\leq 0$ | $L_{6},L_{4},L_{8}$ | $\rho^{\prime},\epsilon^{\prime},v_{x}^{\prime},v_{y}^{\prime},v_{z}^{\prime},B_{x}^{\prime},B_{y}^{\prime}$ | $L_{7},L_{3},L_{5},{L_{1}},{L_{2}}$ | 3 | 7
$0\leq v_{z}\leq c_{s}$ | ${L_{1}},{L_{2}},L_{6},L_{4},L_{8}$ | ${\rho^{\prime},\epsilon^{\prime}},v_{x}^{\prime},v_{y}^{\prime},v_{z}^{\prime},B_{x}^{\prime},B_{y}^{\prime},{B_{z}^{\prime}}$ | $L_{7},L_{3},L_{5}$ | 5 | 8
$c_{s}\leq v_{z}\leq c_{a}$ | $L_{5},{L_{1}},{L_{2}},L_{6},L_{4},L_{8}$ | ${\rho^{\prime},\epsilon^{\prime}},v_{x}^{\prime},v_{y}^{\prime},v_{z}^{\prime},B_{x}^{\prime},B_{y}^{\prime},{B_{z}^{\prime}}$ | $L_{7},L_{3}$ | 6 | 8
$c_{a}\leq v_{z}\leq c_{f}$ | $L_{3},L_{5},{L_{1}},{L_{2}},L_{6},L_{4},L_{8}$ | ${\rho^{\prime},\epsilon^{\prime}},v_{x}^{\prime},v_{y}^{\prime},v_{z}^{\prime},B_{x}^{\prime},B_{y}^{\prime},{B_{z}^{\prime}}$ | $L_{7}$ | 7 | 8
$c_{f}\leq v_{z}$ | $L_{7},L_{3},L_{5},{L_{1}},{L_{2}},L_{6},L_{4},L_{8}$ | ${\rho^{\prime},\epsilon^{\prime}},v_{x}^{\prime},v_{y}^{\prime},v_{z}^{\prime},B_{x}^{\prime},B_{y}^{\prime},{B_{z}^{\prime}}$ | - | 8 | 8

Table 1 lists all the possibilities for inward and outward propagating modes
as $v_{z}$ is varied at a single location on a $z_{\text{min}}$ boundary. When
the bulk velocity $v_{z}$ is directed outward (in the
-$\widehat{\bm{\mathbf{z}}}{}$ direction) and has magnitude greater than the
fast mode speed ($v_{z}<-c_{f}$) then every eigenvalue is negative, and all
the modes propagate outward from the volume of interest, as indicated in the
first row of the Table. In this case, whatever dynamics are happening inside
the volume completely determine the behavior of the plasma in the boundary,
while anything happening outside the boundary cannot influence the boundary or
interior in any way. Any independent estimate for the evolution of the
boundary data cannot be matched unless it is already consistent with the
outgoing information.

If the magnitude of the outflowing bulk velocity is decreased to less than the
fast speed but still greater than the Alfvén speed ($-c_{f}<v_{z}<-c_{a}$)
then the forward propagating fast mode $L_{8}$ has a positive eigenvalue
$\lambda_{8}$, and is therefore inward propagating relative to the
$z_{\text{min}}$ boundary, as indicated in the second row of Table 1. From
Equation 14h and A2, we see that the $L_{8}$ characteristic derivative,
constructed from the $\bm{\mathbf{l}}_{8}$ eigenvector, couples together the
derivatives of seven of the eight primitive variables: everything but the
derivative of the normal component of the magnetic field,
$\partial_{z}B_{z}\equiv B_{z}^{\prime}$. As discussed previously and
reemphasized here, the _single_ inward propagating fast mode with
characteristic derivative $L_{8}$ will simultaneously modify seven of the
eight primitive variables in the boundary, _but not arbitrarily so._ The seven
variables move in lock-step with fixed relative amplitudes given by the
components of the eighth right eigenvector $\bm{\mathbf{r}}_{8}$, all scaled
by the (scalar!) value of the eighth element of characteristic derivative,
$L_{8}$. The external universe dictates the value of that scalar.

At the next level of Table 1, the magnitude of $v_{z}$ is greater than the
slow speed and less than the Alfvén speed, but all three are still negative
valued $(-c_{a}<v_{z}<-c_{s})$. Then, the forward-propagating Alfvén mode
$L_{4}$ joins the forward propagating fast mode $L_{8}$ as a positive mode,
which we identify as an incoming mode through the $z_{\text{min}}$ boundary.
As seen in Equations (14d) and (A2), the Alfvén mode couples together the
derivatives of four primitive variables, but still not the derivative of the
normal magnetic field $B_{z}$. As indicated in the third column of Table 1,
there is no new dependence on the primitive variable derivatives compared to
the previous row. However, because there are now two inward propagating modes,
the external universe can modify the boundary in new ways through the joint
action of both incoming modes. Thus, the space of allowed changes to the
boundary caused by the external universe has increased, while the control the
internal plasma exerts on the boundary has concomitantly decreased. Stepping
the bulk velocity across each of the eigenvalue speeds in turn, this process
continues until the bulk velocity is inward directed at a magnitude exceeding
the local fast mode velocity $c_{f}$, in which case all characteristic modes
are inward propagating. In that case, the external universe completely
dictates the evolution of the plasma in the boundary, and (consistent with
that) the internal dynamics of the plasma have no effect on the boundary
whatsoever.

As a final note about propagation directions, the decomposition into
positive/negative and outward/inward propagating information happens on a
point-by-point and time-by-time basis throughout the plasma. In particular, in
our example of a constant-$z$ plane $\bm{\mathbf{x}}_{b}$, the number of
inward propagating modes is some function $f(x,y,t)$ in that plane.

Returning to the expanded characteristic form of the MHD equations given at
(18), it is interesting to note that each component of the momentum equation
depends only on the difference between forward– and reverse–propagating pairs
of characteristics in each family (e.g., the difference between fast modes
$L_{7}-L_{8}$). In contrast, the continuity, energy, and induction equations
only involve summations between the pairs in each family (e.g.,
$L_{7}+L_{8}$). This same property holds when the transverse terms in
$\bm{\mathbf{C}}$ are fully expanded in characteristic form, that is, when
expanding out the $x$ and $y$ terms in Equation (10) as we did the $z$ terms
in (18). This pairing of sums and differences within each family of
characteristics allows their combinations to describe both the propagation of
disturbances within a plasma as well as the spatial structure of the plasma
itself.

A simple example illustrates this point best. Hydrostatic equilibrium (where
$\partial_{t}v_{z}=0$) can be interpreted as a carefully balanced set of
propagating disturbances that sum to support large scale gradients in the
density and energy (i.e., in the energy and continuity equations) while
simultaneously canceling to zero out any forces, i.e., the terms in the
momentum equation. From the perspective of the characteristics, the
stratification is due to two counter-propagating modes whose net transfer of
momentum is precisely balanced by gravity to produce a static structure. This
behavior is analogous to a zero-frequency standing wave.

Finally, as mentioned above, the characteristic derivative $L_{1}$ is a
somewhat special case. Its role in the present formulation is to preserve the
value of $\nabla\cdot\bm{\mathbf{B}}$ at each spatial location. We included
the empirical lack of magnetic monopoles ($\nabla\cdot\bm{\mathbf{B}}=0$)
implicitly in the induction equation (1d). Thus, if a problem’s initial
condition is solenoidal then $L_{1}$ ensures that the solution remains
solenoidal over time.

##  3 Characteristic boundary conditions: overview and notation

Setting boundary conditions using the characteristics amounts to setting the
value of each characteristic derivative $L_{\sigma}$ associated with an
incoming mode $\sigma$. Thompson (1990) created a catalog of standard
hydrodynamic BCs written in terms of characteristics (no–slip walls, hard
walls, uniform inflow, non–reflecting, constant pressure, etc.) while Poinsot
& Lele (1992) extended the method to include nonhyperbolic terms, e.g.,
viscous terms in the Navier–Stokes equations. Typically, these boundary
conditions are derived by first assuming a physical constraint (such as
constant pressure), calculating the values for the boundary–normal spatial
derivatives of primitive variables implied by that constraint, and then
solving for the incoming characteristic derivatives in terms of those normal
derivatives. The examples in Thompson (1990) essentially amount to applying
the transformation at (17) in situations where an easily calculable constraint
can be provided, but the situation becomes much more complicated for our use
case of data-driven MHD.

Numerous other authors have used the characteristics to formulate boundary
conditions that are approximately transparent to waves, so-called
“nonreflecting boundary conditions.” Hedstrom (1979) appears to be the first
to have done this for a general system of nonlinear hyperbolic equations
(which includes MHD) by setting the amplitudes of incoming characteristic
derivatives to zero (e.g, setting $L_{8}=0$ when the system is at level 1 of
Table 1); Grappin et al. (2000), Landi et al. (2005), Grappin et al. (2008),
or Gudiksen et al. (2011) provide more recent examples in the context of
heliophysics. These (perhaps misnamed) nonreflecting boundary conditions are
only truly nonreflecting for simple waves, but as Hedstrom (1979) already
pointed out, slightly more complex situations such as a weak shock or the
interaction of multiple simple waves begin to present a problem. Those cases
change the external system in ways that should eventually produce reflections.
However, a fundamental assumption of the nonreflecting boundary conditions is
that the external system is unchanged by the information flowing out of the
boundary: it acts as a passive information reservoir that does not directly
influence the behavior inside the simulation.

Of course, this is exactly the _opposite_ behavior we want for data driving,
where the external system explicitly forces the dynamics of the internal
system, often on a rapid timescale. In fact, when applied to photosphere-to-
corona simulations whose purpose is to accurately model processes like the
emergence of active regions, a wave picture is fundamentally inadequate. We
need to model the structural change of the magnetic field, including changes
in the topology, bulk material flows, and so on. Since the region just inside
the boundary is also evolving, the boundary layer experiences a strong forcing
from both sides. For data driving, we therefore need to be able to isolate and
decompose the evolution of the observed boundary into the portion that is due
to internal disturbances propagating outward and the portion due to external
disturbances propagating inward. This requires very clearly identifying the
types of information flowing in all directions through a given point in space.

Some additional notation is necessary to explicitly describe the flow of
information throughout the system. We start again with the MHD equations
written in characteristic form and focusing on the $\widehat{\bm{\mathbf{z}}}$
direction,

$\partial_{t}\bm{\mathbf{U}}+\underbrace{\mathbb{S}\cdot\bm{\mathbf{L}}}_{\mathbb{A}_{z}\cdot\partial_{z}\bm{\mathbf{U}}}+\bm{\mathbf{C}}=\bm{\mathbf{0}},$
(12 revisited)

where we have also included the equivalent expression in terms of the
primitive variables, for reference. The decomposition defined at Equation 4
can be used to split the eigensystem into two subspaces that represent the
positive $(+)$ and negative $(-)$ propagation of information relative to the
diagonalization direction $\widehat{\bm{\mathbf{z}}}$. The split can be
represented in multiple equivalent ways, all based on separating the positive
and negative eigenvalues in the diagonal $\mathbb{\Lambda}$ matrix:

$\mathbb{\Lambda}=\mathbb{\Lambda}_{+}+\mathbb{\Lambda}_{-}\quad\text{with}\quad\mathbb{\Lambda}_{\pm}=\frac{1}{2}(\mathbb{\Lambda}\pm\lvert\mathbb{\Lambda}\rvert).$
(19)

The $\mathbb{\Lambda}_{\pm}$ are two $8\times 8$ diagonal matrices with only
the positive or negative eigenvalues at appropriate locations, and zeros
elsewhere. The “$+$” matrix selects modes for which $\lambda_{i}>0$ and “$-$”
matrix selects modes for which $\lambda_{i}<0$. Characteristics for which
$\lambda_{i}=0$ have zero amplitude, propagate no information, and therefore
require no special treatment in this formalism; see, e.g., Equations (14) and
(18).

From Equation (19) we can also define the positive and negative projection
matrices, or equivalently, the projection operators

$\mathbb{P}_{\pm}=\mathbb{\Lambda}_{\pm}\mathbb{\Lambda}^{-1},$ (20)

where $\mathbb{\Lambda}^{-1}$ is again defined using the pseudoinverse. The
$\mathbb{P}_{+}$ and $\mathbb{P}_{-}$ operators project any object in MHD
statespace onto the positive and negative orthogonal subspaces, respectively.
Like all orthogonal projection operators, these matrices have the property
that repeated projection has no effect,

$\displaystyle\mathbb{P}_{+}\mathbb{P}_{+}=\mathbb{P}_{+},\quad\mathbb{P}_{-}\mathbb{P}_{-}=\mathbb{P}_{-},$
(21) while orthogonality ensures that
$\displaystyle\mathbb{P}_{+}\mathbb{P}_{-}=\mathbb{P}_{-}\mathbb{P}_{+}=\mathbb{0}.$
(22)

Thus defined, the projection operators can be used as needed.

Equation (20) can be inverted to read

$\mathbb{\Lambda}_{\pm}=\mathbb{\Lambda}\mathbb{P}_{\pm}=\mathbb{P}_{\pm}\mathbb{\Lambda}.$
(23)

Using these projection operators, the original coefficient matrices
$\mathbb{A}$ of the MHD Equations (2) can be directionally split into the
positive and negative coefficient matrices that govern the flow of information
in the grid-positive and grid-negative directions:

$\displaystyle\mathbb{A}=\mathbb{A}_{+}+\mathbb{A}_{-}=\mathbb{S}(\mathbb{\Lambda}_{+}+\mathbb{\Lambda}_{-})\mathbb{S}^{-1}.$
(24)

Equivalently, the left and right eigenvector matrices can be split into
positive and negative portions:

$\displaystyle\mathbb{S}=\mathbb{S}_{+}+\mathbb{S}_{-}=\mathbb{S}\mathbb{P}_{+}+\mathbb{S}\mathbb{P}_{-}$
(25)
$\displaystyle\mathbb{S}^{-1}=\mathbb{S}_{+}^{-1}+\mathbb{S}_{-}^{-1}=\mathbb{P}_{+}\mathbb{S}^{-1}+\mathbb{P}_{-}\mathbb{S}^{-1},$
(26)

as well as the characteristic vectors:

$\displaystyle\bm{\mathbf{L}}_{z,\pm}=\mathbb{P}_{z,\pm}\bm{\mathbf{L}}_{z}=\mathbb{P}_{z,\pm}\mathbb{S}^{-1}_{z}\mathbb{A}_{z}\partial_{z}\bm{\mathbf{U}}{},$
(27)

where we have explicitly labeled a normalization direction ($z$) in this last
equation. Similar expressions for all of the above equations apply in the $x$
and $y$ directions, as well.

The eigensystem describes a basis for information propagation, and therefore,
once a subspace has been selected for any element of the eigensystem it
automatically applies to all elements of the system by the associative
property (this follows from repeated projections via (21)). That is to say,
$\mathbb{S}\bm{\mathbf{L}}_{\pm}=\mathbb{S}_{\pm}\bm{\mathbf{L}}=\mathbb{S}_{\pm}\bm{\mathbf{L}}_{\pm}=\mathbb{P}_{\pm}\mathbb{S}\bm{\mathbf{L}}$.

The orthogonality of the positive and negative subspaces means that the
elements of either system do not interact, for example,
$\mathbb{S}_{\mp}\bm{\mathbf{L}}_{\pm}=\bm{\mathbf{0}}$. What this means in
terms of the MHD equations is that the spatial derivatives calculated inside
the characteristic vectors (the $\partial_{z}U_{\zeta}$ terms in (13)) are all
taken in the upwind direction: the positive propagating modes have
characteristic derivatives defined from below while negative propagating modes
have characteristic derivatives defined from above111111More generally, for a
right handed coordinate system the positive (negative) propagating modes have
derivatives defined from the left (right) in the standard sense along each of
the coordinate axes, respectively.. For smooth systems these derivatives will
be equal, but the characteristics intrinsically describe discontinuous
solutions to the MHD equations, as well.

As described just after Equation (15), we use subscripts to represent the
subpaces associated with inward propagating ($\mathcal{I}$) and outward
propagating ($\mathcal{O}$) information through a boundary. For instance, as
indicated in Figure 1, at a $z_{\text{min}}$-boundary we identify inward with
the positive subspace and outward with the negative subspace,
$\mathbb{S}_{\mathcal{I}}=\mathbb{S}_{+}$ and
$\mathbb{S}_{\mathcal{O}}=\mathbb{S}_{-}$, respectively. In general, the $+$
and $-$ notation applies anywhere in space, while we reserve the $\mathcal{I}$
and $\mathcal{O}$ notation for a specified boundary.

As an example of how the projection onto subspaces works in practice, suppose
a location on the $z_{\text{min}}$ boundary has an outward $(v_{z}<0)$ bulk
velocity with amplitude less than the slow-mode speed ($-c_{s}\leq v_{z}\leq
0$). Then, Table 1 shows that $L_{6},\ L_{4},$ and $L_{8}$ correspond to
inward propagating modes, meaning that rows 6, 4, and 8 of the matrix
$\mathbb{S}^{-1}_{+}=\mathbb{S}^{-1}_{\mathcal{I}}$ are given as in Equation
(A5); the remaining rows contain zeros, because the projection
$\mathbb{P}_{+}\mathbb{S}^{-1}$ gives zero for all those elements. Similarly,
$L_{7},\ L_{3},\ L_{5},\ L_{1},$ and $L_{2}$ correspond to outward propagating
modes, the corresponding rows in $\mathbb{S}^{-1}_{\mathcal{O}}$ are given as
in Equation (A5), and the remaining three rows contain only zeros. With the
above notation in hand we are well equipped to describe the flow of
information anywhere within the system, including any boundary.

##  4 Data-driven boundary conditions: continuous case

The data-driven boundary condition (DDBC) problem can be stated concisely as
follows: How is an MHD system evolved _consistently_ from $t$ to $t+\Delta t$
given the MHD state vector everywhere in the system at the _initial_ time $t$,
i.e., $\bm{\mathbf{U}}\left(t,\bm{\mathbf{x}}\right)$, and the MHD state
vector at the boundary(ies) $\bm{\mathbf{x}}_{b}$ at the _new_ time $t+\Delta
t$, i.e., $\bm{\mathbf{U}}\left(t+\Delta t,\bm{\mathbf{x}}_{b}\right)$?
Traditional boundary conditions set values and/or derivatives of the MHD
primitive variables on the boundary, i.e., Dirichlet or Neumann boundary
conditions. These standard boundary conditions effectively determine the
characteristic derivatives $\bm{\mathbf{L}}$ from the values in the domain and
boundary at time $t$. The boundary state evolves purely in terms of values
determined121212Here we have ignored the details and stability of the
numerical time advance. For example, predictor/corrector or implicit schemes
_appear_ to formally use information from future times, e.g.,
$\bm{\mathbf{U}}(\bm{\mathbf{x}}_{b},t+\Delta t/2)$, to advance the MHD state,
but in fact these estimates of intermediate MHD states are calculated
deterministically from the values of $\bm{\mathbf{U}}(\bm{\mathbf{x}}_{b},t)$.
at time $t$

$\partial_{t}\bm{\mathbf{U}}(\bm{\mathbf{x}}_{b},t)=-\mathbb{S}\cdot\bm{\mathbf{L}}(\bm{\mathbf{x}}_{b},t)-\bm{\mathbf{C}}(\bm{\mathbf{x}}_{b},t)\Longrightarrow\bm{\mathbf{U}}(\bm{\mathbf{x}}_{b},t+\Delta
t).$ (28)

The philosophy here is that $\bm{\mathbf{U}}(\bm{\mathbf{x}},t)$,
$\mathbb{S}$, $\bm{\mathbf{L}}(\bm{\mathbf{x}}_{b},t)$, and
$\bm{\mathbf{C}}(\bm{\mathbf{x}}_{b},t)$ are known and used to calculate the
final unknown boundary state $\bm{\mathbf{U}}(\bm{\mathbf{x}}_{b},t+\Delta
t)$.

From the perspective of the method of characteristics, these standard boundary
conditions determine the incoming characteristics purely from the evolution of
the interior of the domain. For example, homogeneous Neumann boundary
conditions could be used to determine the values of the primitive variables
outside the simulation in terms of primitive variables inside the simulation
by requiring that $\partial\bm{\mathbf{v}}/\partial
z|_{\bm{\mathbf{x}}_{b}}=0$. Traditional boundary conditions therefore
correspond to an assumption (often ad hoc) about the evolution of the external
universe in response to the evolution of the MHD system.

In contrast to traditional boundary conditions described above, for DDBCs we
do not want to determine the incoming characteristics
$\bm{\mathbf{L}}_{\mathcal{I}}(\bm{\mathbf{x}}_{b},t)$ from the evolution of
the interior of the domain because we want the external universe to inject
_independent information_ into the MHD system through the boundary. Instead we
want to choose the incoming characteristics at the boundary such that the
superposition of the incoming and outgoing characteristics leads to evolution
of the boundary from $\bm{\mathbf{U}}\left(t,\bm{\mathbf{x}}_{b}\right)$
towards the prescribed target state $\bm{\mathbf{U}}\left(t+\Delta
t,\bm{\mathbf{x}}_{b}\right)$:

$\mathbb{S}\cdot\bm{\mathbf{L}}(t,\bm{\mathbf{x}}_{b})=-\partial_{t}\bm{\mathbf{U}}(t,\bm{\mathbf{x}}_{b})-\bm{\mathbf{C}}(t\bm{\mathbf{x}}_{b})\Longrightarrow\bm{\mathbf{U}}(t+\Delta
t,\bm{\mathbf{x}}_{b}).$ (29)

The philosophy here is that
$\partial_{t}\bm{\mathbf{U}}(t,\bm{\mathbf{x}}_{b})$, and through it the final
state of the boundary $\bm{\mathbf{U}}(t+\Delta t,\bm{\mathbf{x}}_{b})$, and
$\bm{\mathbf{C}}(t,\bm{\mathbf{x}}_{b})$ are known and used to calculate the
best estimate of the incoming characteristic derivatives
$\bm{\mathbf{L}}(t,\bm{\mathbf{x}}_{b})$ to evolve the boundary state
_towards_ $\bm{\mathbf{U}}(t+\Delta t,\bm{\mathbf{x}}_{b})$. Put another way,
the time evolution of the boundary
$\partial_{t}\bm{\mathbf{U}}(t,\bm{\mathbf{x}}_{b})$ is known independently of
(28) and (conveniently) the MHD equations as represented in Equation (12) in
terms of the characteristic derivatives $\bm{\mathbf{L}}$ in Equation (13)
give no preference to either the temporal or the spatial derivatives, which
means the unknown spatial derivatives can be solved for. The result is that,
by carefully and correctly choosing the spatial derivatives one can evolve the
system to match the estimated boundary evolution. We again stress that, by
using the characteristic derivatives rather than the primitive variables
directly, we are restricting ourselves to the subset of boundary evolutions
allowed by MHD. This will result in exact matching of the estimated boundary
evolution provided the estimated evolution is already self-consistent with the
MHD equations, but will only _approximately_ match when the estimated
evolution contains errors.

The projection operator formalism introduced in the preceding §3 can be
applied to make these arguments mathematically rigorous. Equation (28) can be
rearranged to determine the spatial derivatives (via characteristic
derivatives $\bm{\mathbf{L}}$) necessary to achieve the time derivative
$\partial_{t}\bm{\mathbf{U}}(t,\bm{\mathbf{x}}_{b})$ which will produce the
(possibly approximate) transformation
$\bm{\mathbf{U}}\left(t,\bm{\mathbf{x}}_{b}\right)\rightarrow\bm{\mathbf{U}}\left(t+\Delta
t,\bm{\mathbf{x}}_{b}\right)$. To achieve this we substitute the decomposition
(25) into (28), multiply through by $\mathbb{S}^{-1}$, and use the fact that
$\mathbb{S}^{-1}\mathbb{S}_{\mathcal{I},\mathcal{O}}=\mathbb{P}_{\mathcal{I},\mathcal{O}}$
to get

$\displaystyle\mathbb{S}^{-1}\partial_{t}\bm{\mathbf{U}}+\mathbb{P}_{\mathcal{I}}\bm{\mathbf{L}}+\mathbb{P}_{\mathcal{O}}\bm{\mathbf{L}}+\mathbb{S}^{-1}\bm{\mathbf{C}}=\bm{\mathbf{0}},$
(30)

where each term in the expression is evaluated at the boundary
$\bm{\mathbf{x}}_{b}$. Multiplying through by $\mathbb{P}_{\mathcal{I}}$ the
above equation decomposes completely for the incoming system

$\displaystyle\bm{\mathbf{L}}_{\mathcal{I}}=-\mathbb{S}^{-1}_{\mathcal{I}}\bigl{(}\partial_{t}\bm{\mathbf{U}}+\bm{\mathbf{C}}\bigr{)};$
(31) for completeness we could multiply by $\mathbb{P}_{\mathcal{O}}$ to do
the same for the outgoing system
$\displaystyle\bm{\mathbf{L}}_{\mathcal{O}}=-\mathbb{S}^{-1}_{\mathcal{O}}\bigl{(}\partial_{t}\bm{\mathbf{U}}+\bm{\mathbf{C}}\bigr{)},$
(32)

but the latter has no explicit dependence on anything outside of the domain
and therefore it cannot be used to inject independent information into the
domain. Of course, Equation (31) is only rigorously valid when our estimate
for $\partial_{t}\bm{\mathbf{U}}(t,\bm{\mathbf{x}}_{b})$ is completely
consistent with the evolution of the MHD system. This is unlikely when
$\bm{\mathbf{U}}\left(t,\bm{\mathbf{x}}_{b}\right)$ corresponds to
observations. We discuss this point in the next section, which focuses on the
numerical integration in the discrete case.

We can also transform (31) into the space of spatial derivatives of primitive
variables by substituting in for
$\bm{\mathbf{L}}_{z\mathcal{I}}=\mathbb{\Lambda}\mathbb{P}_{\mathcal{I}}\mathbb{S}^{-1}\partial_{z}\bm{\mathbf{U}}$
and then solving for those derivatives, after which we find (written in
several equivalent forms)

$\displaystyle\partial_{z}\bm{\mathbf{U}}$
$\displaystyle=-(\mathbb{S}\mathbb{\Lambda}\mathbb{P}_{\mathcal{I}}\mathbb{S}^{-1})^{-1}\bigl{(}\partial_{t}\bm{\mathbf{U}}+\bm{\mathbf{C}}\bigr{)}$
(33a)
$\displaystyle=-\mathbb{S}\mathbb{\Lambda}_{\mathcal{I}}^{-1}\mathbb{S}^{-1}\bigl{(}\partial_{t}\bm{\mathbf{U}}+\bm{\mathbf{C}}\bigr{)}$
(33b)
$\displaystyle=-\mathbb{A}_{\mathcal{I}}^{-1}\bigl{(}\partial_{t}\bm{\mathbf{U}}+\bm{\mathbf{C}}\bigr{)}.$
(33c)

Formally, the continuous problem is now solved: if we know the time derivative
$\partial_{t}\bm{\mathbf{U}}(t,\bm{\mathbf{x}}_{b})$ and the transverse and
inhomogeneous terms $\bm{\mathbf{C}}(t,\bm{\mathbf{x}}_{b})$ then we can
directly solve for the (external, one-sided) spatial derivatives
$\partial_{z}\bm{\mathbf{U}}(t,\bm{\mathbf{x}}_{b})$. What’s more, we can do
so either in terms of the characteristic derivatives of inward propagating
modes through Equation (31) or directly in terms of spatial derivatives of the
primitive variables via Equation (33c). In the latter case, one cannot
arbitrarily set derivatives of primitive variables but is instead restricted
to the subspace of spatial derivatives that are linear combinations of the
characteristic derivatives of inward propagating modes. This property is
enforced by the inverse of the positive subspace matrix,
$\mathbb{A}_{\mathcal{I}}^{-1}=\mathbb{S}\mathbb{P}_{\mathcal{I}}\mathbb{\Lambda}^{-1}\mathbb{S}^{-1}$.
On the other hand, if a given primitive variable is not a function of any
incoming characteristic then its exterior derivative is unconstrained: it has
no influence on the dynamics at or within the boundary and can therefore be
set arbitrarily. The most extreme example of this is when material leaves a
boundary faster than the fast mode speed, so that all the characteristics are
outward propagating. In that case, it doesn’t matter how the variables are set
outside the simulation because there is no incoming information. This is a
common situation, for instance, at the outer boundary of solar wind
simulations (Grappin et al., 2000).

To complete the description of a characteristics-based boundary condition, the
incoming characteristic derivatives are calculated using Equation (31) (or
equivalently for many simulations, the values for simulation ghost cells are
calculated using (33c)) while the outgoing characteristic derivatives
$L_{\mathcal{O}{},\sigma}$ (given in Table 1) are still calculated using
equations (14). Then, all the $L_{\sigma}$ are substituted into the MHD
equations (28) and the boundary values are advanced in time by standard
time–integration methods. Because the MHD state vector on the boundary
$\bm{\mathbf{U}}(t,\bm{\mathbf{x}}_{b})$ varies as a function of space,
different locations on the boundary will simultaneously sit at different rows
of Table 1 and therefore have different numbers of incoming and outgoing
characteristics. In the same way, the number of incoming characteristics at a
given location will vary as a function of time as the system evolves.

Several points are worth emphasizing at this juncture:

  * •

The MHD equations give no preference to either spatial or temporal
derivatives. If some subset of those derivatives are known, the others can be
solved for.

  * •

There is no one-to-one mapping between characteristics and primitive
variables, so arbitrary constraints on the primitive variables could lead to a
consistent, under-, or over-determined system of equations when applied to
MHD.

  * •

Boundary conditions on the primitive variables that are also consistent with
the MHD equations are necessarily expressible in terms of incoming
characteristics only.

  * •

Possibly the only guaranteed safe way to apply arbitrary boundary conditions
is via the method of characteristics.

The first point means we can use an estimate for the temporal evolution of the
primitive variables on a boundary—for instance from a time series of
observations—to find a set of incoming characteristic derivatives that are
maximally consistent with the observed evolution at the boundary; this amounts
to exchanging known temporal derivatives for unknown spatial derivatives in
the MHD equations. However, the second point cautions that the spatial
derivatives of multiple primitive variables are simultaneously constrained by
a single incoming characteristic, and the derivative of a single primitive
variable simultaneously contributes to multiple characteristics: the
derivatives are coupled together. While the primitive variables themselves
make no reference to the amount of independent information that can be self-
consistently supplied to the system, this is exactly what the characteristics
provide. The number of incoming characteristics for a given boundary state
$\bm{\mathbf{U}}(t,\bm{\mathbf{x}}_{b})$ is equivalent to the amount of
independent information that can be supplied to the system. The incoming
eigvenvectors then give the precise, coupled forms that incoming information
can take for an MHD system. This leads to our third point, that any self-
consistent boundary condition must be expressible in terms of the incoming
characteristics. Conversely, if a boundary condition is not expressible in
terms of only incoming characteristics then the system being described is
either under- or over- determined. Finally, based on the above chain of logic,
we posit that the characteristics likely provide the only generic method to
guarantee that an arbitrary boundary condition will satisfy the MHD equations.

##  5 Data-driven boundary conditions: finite volume discretization, time-
integration, and data-optimization

Figure 2: Schematic diagram of the data driving algorithm. The upper row shows
several times from a _Ground Truth (GT)_ simulation and highlights the
vertical magnetic field in the layer that will be used for data driving (gray
scale). Time flows left to right. The time between each state O in the _GT_
simulation is $\delta t$ while the time between the extracted data series used
for driving ($t_{O1},t_{O2},\ldots$) is $\Delta t$. Blue arrows represent the
intermittent injection of data extracted (or “observed”) from the _GT_
simulation and used in the _Driving_ simulation (Dr) The bottom row shows
several times from the resulting _Data-Driven_ simulation from a side-on 3D
perspective which shows vertical velocities in the vertical cut-plane and the
magnetic field in the driving plane. Field lines highlight the 3D magnetic
structure for each simulation. See text for the full relation between the _GT_
, _Driving_ , and _Data-Driven_ simulations.

Thus far we have described the problem of data-driven simulations with only
minimal reference to the specifics of either the observations or the
simulations. Those details must now be brought to the forefront, but we keep
the discussion as general as possible so that the technique may later be
applied to diverse observational data sets or numerical schemes. Our goal in
the present investigation is to validate our data driving approach, and
therefore an independently run “ground truth” simulation with a larger domain
will serve as a stand-in for observations used to drive the boundary. Figure 2
provides a cartoon-overview of our data driving algorithm and its validation
using the independent Ground Truth simulation, _GT_.

In our numerical formulation, data driving results from performing a special
time integration of the MHD equations in a single layer of driving cells. In
this section we describe the details of that time integration, where the
driving layer is taken to be the $k=0$ cell index in the $z$-direction. The
integration extends temporally between one observation of the boundary at,
say, time $t_{O1}$ and the next observation at time $t_{O2}$, as indicated in
Figure 2. It is thus the numerical implementation of our data-driven boundary
condition, DDBC, which solves for the set of incoming characteristic
derivatives at each space-time point on the driven boundary,
$\bm{\mathbf{L}}_{z\mathcal{I}}(t,\bm{\mathbf{x}}_{b})$. At the heart of this
task is an optimization approach that finds the values of incoming
characteristic derivatives that achieve the best possible consistency between
the observations and the MHD equations. Those incoming characteristic
derivatives then determine the optimal evolution of primitive variables in the
$k=0$ driving layer given the limited set of temporal observations of the MHD
variables at that layer, i.e., the intermittently available observations (or
the intermittently extracted driving layer from the _GT_ simulation). The use
of intermittent data is depicted as blue arrows in Figure 2. As will become
apparent, the boundary condition is effectively applied at the
$k=-\frac{1}{2}$ interface, so the $k=-1$ cell is a ghost cell that does not
directly come into play in our data driving algorithm; it is, however, used to
couple separate simulations together.

In more detail, our implementation of data driving uses two simulations that
partially overlap in space: the _Driving_ simulation and the _Data-Driven_
simulation. The _Driving_ simulation spans a narrow band in height and is used
exclusively to update the driving layer. The _Data-Driven_ simulation spans
the entire coronal volume above the driving layer and is numerically coupled
to the _Driving_ simulation. In addition to these, the present work also
includes a third simulation, the _GT_ simulation, that is a stand-in for the
observations. This latter simulation does not exist when applying our method
to the real Sun.

Figure 2 sketches the relation between the Ground Truth simulation (top row),
the time-integration of the MHD equations in the $k=0$ driving layer in the
_Driving_ simulation (middle region), and the resulting _Data-Driven_
simulation (bottom row). Full details of the simulations are given in §6. For
now, the important points are that: (1) the driving layer sits at cell center
index $k=0$; (2) the logically separate _Driving_ simulation (represented by
the Dr: O……O row of the Figure), spanning $-1\leq k\leq 2$ in the
$z$-direction, is used to numerically update in the $k=0$ driving layer
($k=-1$ is a ghost-cell for this simulation); (3) the _Data-Driven_ simulation
(bottom portion of the Figure), spanning cells $1\leq k\leq N_{z}$ in the $z$
direction with numerical lower boundary at $k=\frac{1}{2}$, represents the
simulated photosphere-to-corona MHD system above the driving layer; and (4)
the _Ground Truth (GT)_ simulation (top portion of the Figure) spans the much
larger range $k_{\text{min, GT}}<k<N_{z}$ with $k_{\text{min, GT}}\ll-1$, from
which the driving layer at $k=0$ is occasionally extracted and supplied to the
_Driving_ simulation. A more realistic _GT_ simulation would span from the
convection zone to the corona, similar to what was used in Toriumi et al.
(2020).

Figure 3 describes the vertical overlap of numerical grids in the vicinity of
the driving layer. All three simulations use a common set of indices, i.e., a
layer $k$ refers to the same height in all three simulations. We use integer
valued indices $(i,j,k)$ for cell centered locations and half-integer values
for cell faces (e.g., the left and right $x$-faces $i\pm\frac{1}{2}$, front
and back $y$-faces $j\pm\frac{1}{2}$, and top and bottom $z$-faces
$k\pm\frac{1}{2}$). In the figure, blue regions represent the extent of the
_GT_ simulation, red the _Driving_ simulation, and green the _Data-Driven_
simulation. Horizontal black dashed lines show the cell interface locations
and are labeled by their half-integer indices in the $z$ direction.

Our numerical scheme to implement DDBCs is ultimately cast in finite-volume
form, but this specific choice is not required for characteristic based
boundary conditions in general or data-driven boundary conditions in
particular131313Indeed, we originally implemented a finite-difference scheme,
which worked well in some instances, but not others.. We call the
characteristic-based finite-volume code CHAR, and our implementation of data-
driven boundary conditions in this code CHAR-DDBC.

In brief, a small-domain _Driving_ MHD simulation (represented by the dots in
the Dr line in Figure 2 and by the red region in Figure 3) runs currently with
and coupled to the main _Data-Driven_ MHD simulation (green region in Figure
/2). The _Driving_ simulation is where the DDBC is actually implemented.
Appropriate averaging of spatially overlapping cells couples the two
simulations together. In a real-world scenario these are the only two
simulations that would run, and observational data would be ingested into cell
$k=0$ of the _Driving_ simulation. For testing purposes, a third MHD
simulation, (_GT_ ; blue), is used to extract synthetic observations that
serve as input to the _Driving_ simulation. Direct comparison between the _GT_
simulation and the _Data-Driven_ simulation in the cells $k\geq 0$ is used to
assess how accurately the DDBC performs. In practical application to data-
driven simulations of the Sun, the extracted layer from the _GT_ will be
replaced with a time series of plasma parameters derived from observations of
the solar photosphere. The main benefit of the logical separation into three
simulations is that it keeps the key elements of our data driving method
distinct from a specific numerical code and makes the method adaptable to a
wide variety of numerical techniques.

Figure 3: Relationship between the vertical grids ($z$ direction, indexed by
$k$) in the three overlapping simulations. The Ground Truth simulation (_GT_ ;
blue) extends over a large domain. The Driven Simulation (green) extends over
the upper portion of the _GT_ domain. The _Driving_ simulation (red) extends
for four full computational cells that span the boundary of the Driven
simulation. Cell faces are shown as horizontal dashed lines. Cell center
locations are given integer indices and cell interfaces are given half-integer
indices. The $k=0$ cell-center layer is the driving layer and $k=\frac{1}{2}$
is the interface boundary of the Driven domain. The desired incoming
characteristics are solved for at the $k=-\frac{1}{2}$ interface,
$\bm{\mathbf{L}}_{z,\mathcal{I}}(t,i,j,k=-\frac{1}{2})$.

In addition to the spatial dimensions, there are two time sequences to
consider. First, there is the time step between one observation and the next
observation, $\Delta t$. These are the only times at which boundary data
(i.e., the driving layer) is known,
$\bm{\mathbf{U}}(t,\bm{\mathbf{x}}_{b})\ldots\bm{\mathbf{U}}(t+\Delta
t,\bm{\mathbf{x}}_{b})$. In Figure 2 these are represented, for example, by
$t_{O1},\ t_{O2},\ldots$ with $\Delta t=t_{O2}-t_{O1}$. Second, there is the
simulation time step $\delta t$, given by the simulation’s Courant condition,
which is typically much shorter than the cadence of boundary observations. For
example, the standard HMI data series hmi.ME_720s_fd10 provides vector
magnetic field data at a cadence of 720 seconds, and the less frequently
available high-cadence HMI data series provides data at 130 seconds, while a
typical MHD code may have a Courant limited timestep of 1 second or less,
depending on the simulation parameters. The upper portion of Figure 2 shows
several times from the _GT_ simulation, the middle portion represents the
occasional extraction (at cadence $\Delta t$) of the $k=0$ layer from the _GT_
simulation and its use in the driving simulation, and the lower portion shows
the result of using the driving layer to produce the driven simulation in the
reduced spatial domain.

In general there will be $\Gamma$ simulation steps in the time $\Delta t$
between one observation and the next. The driving simulation performs the
optimization procedure described below $\Gamma$ times to determine the
incoming characteristic derivatives at each time and integrates the MHD
equations using the calculated derivatives. This process will drive the _Data-
Driven_ simulation, in a self–consistent way, between the two observation
times. Those $\Gamma$ steps comprise one observation cycle, and we repeat this
process for as many observational cycles as exist.

In the following, we only need to describe a single observational cycle, which
has an observation at the start of the cycle at time $t_{O1}$ and an
observation at the end of the cycle at time $t_{O2}$. The total duration of
the cycle between two observations is $\Delta t=t_{O2}-t_{O1}$. The
$\eta^{\text{th}}$ simulation step from one observed time to the next we
designate by $t_{O1+\eta}$, where $0\leq\eta\leq\Gamma-1$, so that
$t_{O1{}+\Gamma}$=$t_{O2}$. A single application of the optimization method is
used to update the driving layer at $k=0$ between $t_{O1+\eta}$ and
$t_{O1{}+\eta+1}$, i.e., for the Courant-limited simulation time step $\delta
t=t_{O1{}+\eta+1}-t_{O1+\eta}$. The optimization method is repeated until the
_Data-Driven_ simulation has been numerically integrated from $t_{O1}$ to
$t_{O2}$. Note that neither $\delta t$ nor $\Gamma$ is known in advance, but
this is unimportant.

We designate the observed state in the driving layer with a hat,
$\widehat{\bm{\mathbf{U}}}$, so a single observational cycle is bounded by
$\widehat{\bm{\mathbf{U}}}(t_{O1},\bm{\mathbf{x}}_{b})$ and
$\widehat{\bm{\mathbf{U}}}(t_{O2},\bm{\mathbf{x}}_{b})$. The _Data-Driven_
simulation’s boundary state may not exactly match any of the observed states,
so that in general
$\bm{\mathbf{U}}(t_{O1},\bm{\mathbf{x}}_{b})\neq\widehat{\bm{\mathbf{U}}}(t_{O1},\bm{\mathbf{x}}_{b})$.
We wish to drive the current state
$\bm{\mathbf{U}}(t_{O1},\bm{\mathbf{x}}_{b})$ toward
$\widehat{\bm{\mathbf{U}}}(t_{O2},\bm{\mathbf{x}}_{b})$ as best we can in a
“hound chasing the hare” fashion.

The observations in the driving layer allow us to estimate the time derivative
we are attempting to reproduce by the boundary driving, i.e., this is the
_independent estimate_ discussed in §4 in connection with Equation (29). We
designate the target time derivative by an over-bar,
$\overline{\dot{\bm{\mathbf{U}}}}$, where the over-dot denotes a generic time
derivative. The simplest estimate is a linear interpolation over the two
observed states,

$\dot{\bm{\mathbf{U}}}\approx\dot{\bm{\mathbf{U}}}{}_{\text{simple}}=\frac{\widehat{\bm{\mathbf{U}}}(t_{O2})-\widehat{\bm{\mathbf{U}}}(t_{O1})}{\Delta
t},$ (34)

so that the estimated time derivative is constant over the cycle. A somewhat
better approximation is also a linear estimate but uses the current value of
the state variable (at step $\eta$) instead of the value at the beginning of
the observation cycle:

$\dot{\bm{\mathbf{U}}}\approx\dot{\bm{\mathbf{U}}}{}_{\text{current
time}}=\frac{\widehat{\bm{\mathbf{U}}}(t_{O2})-\bm{\mathbf{U}}_{O1{}+\eta}}{t_{O2}-t_{O1+\eta}}.$
(35)

This has the effect of constantly nudging the current solution in the best
direction to achieve the target final state
$\widehat{\bm{\mathbf{U}}}(t_{O2})$. Higher order estimates could be explored,
but we will take
$\overline{\dot{\bm{\mathbf{U}}}}{}=\dot{\bm{\mathbf{U}}}{}_{\text{current
time}}$ from here on.

Regardless of how it is estimated, $\overline{\dot{\bm{\mathbf{U}}}}$ provides
a numerical target value for the boundary state vector at the next simulation
time, $t=t_{O1+\eta}+\delta t$, for instance using a forward Euler update

$\overline{\bm{\mathbf{U}}}(t_{O1+\eta}+\delta
t,\bm{\mathbf{x}}_{b})=\bm{\mathbf{U}}(t_{O1+\eta},\bm{\mathbf{x}}_{b})+\delta
t\overline{\dot{\bm{\mathbf{U}}}}(t_{O1+\eta},\bm{\mathbf{x}}_{b}).$ (36)

Just as for the time derivative, we could also use a higher order integration
method to define a higher order estimate for the target boundary state. But
again, in this investigation we simply use Equation (36) to define the target
boundary state.

Figure 4: Key elements of finite volume scheme to update cell-centered
primitive variables $\bm{\mathbf{U}}_{ij0}$ in the driving layer $k=0$. The
directionally split coefficient matrices, $\mathbb{A}_{\pm}$, are depicted at
each interface (the y-faces are omitted for clarity). Integer and half-integer
indices refer to cell center and cell face locations, respectively, and
several sets of vertex indices are included as examples. Matrices associated
with green arrows represent information coming into to the $ij0$ cell and
those associated with blue arrows represent information leaving that cell. The
matrix at the lower face ($k=-\frac{1}{2}$) is given a red arrow to indicate
that it needs to be calculated using independent information; i.e., the
incoming information represented by
$\mathbb{A}_{\mathcal{I},(i,j,-\frac{1}{2})}$ is the boundary condition. The
$k=-1$ cell is not shown as it is not used explicitly in the current numerical
scheme (see discussion at the end of §5).

With the target state and target time derivative defined, we turn to the
details of integrating the cells in the driving layer forward in time in a way
that is consistent with both the estimated target states and the directional
flow of information given by the characteristics. To achieve this, we apply
the characteristic eigendecomposition and projection formalism introduced in
§2 and §3 to each face of a computational cell and thereby develop a finite
volume scheme based on the characteristics, CHAR. This essentially treats each
cell face as a boundary and we cyclically apply our formalism at each face in
order to update the cell center values of the primitive variables. By itself,
this produces a fully functional 3D MHD code that works anywhere in a
simulation interior. Data driving results from applying the formalism
developed in §4 to the one cell face through which new information enters the
system: the lower $z$ face in the driving layer. That new information is, of
course, the target time derivative $\overline{\dot{\bm{\mathbf{U}}}}$ just
described (or, equivalently, the target state $\overline{\bm{\mathbf{U}}}$).

Figure 4 depicts the CHAR numerical scheme applied to a computational cell in
driving layer $k=0$. The MHD variables in CHAR are cell averaged quantities
that live at cell centers,
$\bm{\mathbf{U}}_{ijk}=\bm{\mathbf{U}}(\bm{\mathbf{x}}_{i,j,k})=\bm{\mathbf{U}}(x_{i},y_{j},z_{k})$.
The highlighted cell in the figure is an element of the driving layer at
$k=0$, so
$\bm{\mathbf{U}}_{ij0}=\bm{\mathbf{U}}(\bm{\mathbf{x}}_{b})=\bm{\mathbf{U}}(x_{i},y_{j},z_{b})$.
The time update to a given cell uses information from the six surrounding cell
centers, as well as the average state at each cell interface (the $y$
direction is omitted, for clarity). As mentioned above, the cell faces are
identified by half-integer index labels, e.g., the right-side face of cell
$(i,j,k)$ is labeled $(i+\frac{1}{2},j,k)$. The objects in the Figure provide
several more examples; for example, the right-back-top vertex of the
highlighted $ij0$ cell is labeled $(i+\frac{1}{2},j+\frac{1}{2},\frac{1}{2})$.

The flow of information across each cell interface is defined in terms of the
eigendecomposition described in §3 that splits the system into positive and
negative subspaces at each interface in the face-normal
($\widehat{\bm{\mathbf{n}}}$) direction. Each eigendecomposition is calculated
using the _interface_ MHD state, defined as the average of the two surrounding
cell centered states:

$\bm{\mathbf{U}}_{n\pm\frac{1}{2}}=\frac{1}{2}(\bm{\mathbf{U}}_{n}+\bm{\mathbf{U}}_{n\pm
1}),\qquad n\in\\{i,j,k\\}.$ (37)

For example, at the $i+\frac{1}{2}$ face, whose normal direction is
$\widehat{\bm{\mathbf{x}}}$, the eigendecomposition is

$\displaystyle\mathbb{A}_{x,i+\frac{1}{2}}=\mathbb{S}_{x,i+\frac{1}{2}}\mathbb{\Lambda}_{x,i+\frac{1}{2}}\mathbb{S}^{-1}_{x,i+\frac{1}{2}}=\mathbb{A}_{x+,i+\frac{1}{2}}+\mathbb{A}_{x-,i+\frac{1}{2}},$
(38)

where all versions of the $\mathbb{A}_{x}$, $\mathbb{S}_{x}$, and
$\mathbb{\Lambda}_{x}$ matrices are constructed using the $x$-interface-
averaged MHD state $\bm{\mathbf{U}}{}_{i+\frac{1}{2},j,k}$. The second
equality shows the coefficient matrix
$\mathbb{A}_{x}=\mathbb{A}_{x+}+\mathbb{A}_{x-}$ split by projecting onto the
positive and negative subspaces as in Equation (24) (recall from the
discussion following Equation (27) that the projection of any of $\mathbb{A}$,
$\mathbb{\Lambda}$, $\mathbb{S}$, or $\bm{\mathbf{L}}$ onto a subspace
necessarily projects the other terms onto the same subspace in the same
expression).

In Figure 4, the split at each cell face is represented by different colored
arrows (for clarity we have only shown the $x$ and $z$ faces, but the front
and back faces of the cell in the $y$ direction are to be included, as well).
The colors distinguish objects in _incoming_ subspaces (green) from _outgoing_
subspaces (blue), from the perspective of the $ij0$ cell. Thus, at the right
cell face $(i+1/2,j,0)$ the incoming subspace in the $x$ direction to cell
$ij0$ is the negative subspace, $\mathbb{A}_{x-,(i+1/2,j,0)}$, while at the
left cell face $(i-1/2,j,0)$ the incoming subspace is the positive subspace,
$\mathbb{A}_{x+,(i-1/2,j,0)}$. The lower face of a $k=0$ cell is a special
case, because new information from outside the simulation must be supplied at
that location, hence we label it red in the figure.

The characteristics are incorporated into the calculation by projecting the
normal derivative of the MHD state vector onto the split eigensystem at the
interface state, e.g., in the $\widehat{\bm{\mathbf{n}}}$ direction,

$\displaystyle\mathbb{A}_{n}\partial_{n}\bm{\mathbf{U}}{}$
$\displaystyle=(\mathbb{A}_{n+}+\mathbb{A}_{n-})\partial_{n}\bm{\mathbf{U}}{}$
(39)
$\displaystyle=\mathbb{S}_{n}(\mathbb{P}_{n+}+\mathbb{P}_{n-})\mathbb{\Lambda}_{n}\mathbb{S}_{n}^{-1}\partial_{n}\bm{\mathbf{U}}$
(40) $\displaystyle=\mathbb{S}_{n}(\bm{\mathbf{L}}_{n+}+\bm{\mathbf{L}}_{n-})$
(41)
$\displaystyle=\mathbb{S}_{n+}\bm{\mathbf{L}}_{n+}+\mathbb{S}_{n-}\bm{\mathbf{L}}_{n-}.$
(42)

The normal derivatives of the primitive variables are centered on the cell
interface,

$\bigl{(}\partial_{n}\bm{\mathbf{U}}\bigr{)}_{l+\frac{1}{2}}=\frac{\bm{\mathbf{U}}_{l+1}-\bm{\mathbf{U}}_{l}}{\Delta
x_{l}},$ (43)

and these derivatives are projected onto the positive and negative
eigensystems of the interface. Here, $\Delta x_{l}=x_{l+1}-x_{l}$ is the
distance between cell center locations on either side of the cell face, and
the index $l$ is the index associated with the normal direction
$\widehat{\bm{\mathbf{n}}}$ (as it will be throughout the following
discussion). The positive subspace $\mathbb{A}_{+}$ (and associated
derivatives $\bm{\mathbf{L}}_{+}$) contributes to the time-update of the cell
on the right side of an interface, while the negative subspace
$\mathbb{A}_{-}$ (and associated $\bm{\mathbf{L}}_{-}$) contribute to the
update of the cell on the left side of the interface.

Taking the $\widehat{\bm{\mathbf{x}}}$ direction as an example, cell $i$ is
updated according to its incoming characteristics, which means the positive
subspace from the interface at $i-\frac{1}{2}$ and the negative subspace from
the interface at $i+\frac{1}{2}$. Each of these terms are labeled at the
appropriate faces in Figure 4. Analogous expressions hold in
$\widehat{\bm{\mathbf{y}}}$ and $\widehat{\bm{\mathbf{z}}}$ so that the total
time update to the cell-centered state variable at cell $(i,j,k)$ is the sum
of incoming characteristics over all six cell faces, plus the inhomogeneous
terms:

$\displaystyle\partial_{t}\bm{\mathbf{U}}_{ijk}+\sum_{(n,l)\in(x,i),(y,j),(z,k)}\biggl{(}\mathbb{A}_{n+,l-\frac{1}{2}}\bigl{(}\partial_{n}\bm{\mathbf{U}}\bigr{)}_{l-\frac{1}{2}}+\mathbb{A}_{n-,l+\frac{1}{2}}\bigl{(}\partial_{n}\bm{\mathbf{U}}\bigr{)}_{l+\frac{1}{2}}\biggr{)}+\bm{\mathbf{D}}_{ijk}=\bm{\mathbf{0}},$
(44)

where the derivative direction $n$ is taken in the direction normal to each
face, paired with each combination of $l\pm\frac{1}{2}$, so that $x$ is paired
with $i$, $y$ is paired with $j$, and $z$ is paired with $k$. Equation (44) is
a finite volume version of the characteristic form of the discretized MHD
equations. For $k>0$ (away from the driving layer) every term can be
calculated and so the time derivative can be calculated and numerically
integrated—this is just a standard expression for the numerical update
anywhere in a simulation interior. Applied to the $k=0$ driving layer, the
only unknown terms in (44) are associated with the incoming subspace at the
$k=-\frac{1}{2}$ face, that is,
$\mathbb{A}_{z\mathcal{I},-\frac{1}{2}}\bigl{(}\partial_{z}\bm{\mathbf{U}}\bigr{)}_{-\frac{1}{2}}$
where we now make the notational switch
$\mathbb{A}_{z+}\rightarrow\mathbb{A}_{z\mathcal{I}}$ to reinforce that the
positive subspace at face $k=-\frac{1}{2}$ is incoming. This incoming matrix
is marked by the red arrow in Figure 4.

To make progress we apply the finite volume update to the $k=0$ cells,
separate out the unknown terms from the summation, and solve for the time
derivative to write

$\displaystyle\partial_{t}\bm{\mathbf{U}}_{ij0}=-\mathbb{A}_{z\mathcal{I},ij,-\frac{1}{2}}(\partial_{z}\bm{\mathbf{U}})_{ij,-\frac{1}{2}}-\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}$ (45)

where the $l\pm$ notation indicates that summation is taken over, and
derivatives are evaluated at, all cell faces except the $k=-\frac{1}{2}$ face.
Next we define a residual vector $\mathcal{R}$ as the difference between the
finite volume calculation of the time derivative and our observationally
estimated time derivative from (35):

$\displaystyle\mathcal{R}$
$\displaystyle=\overbrace{\partial_{t}\bm{\mathbf{U}}}^{\text{calculation}}-\overbrace{\overline{\dot{\bm{\mathbf{U}}}}}^{\text{estimate}}$
(46) $\displaystyle=-\hbox
to0.0pt{$\displaystyle\overbrace{\phantom{\mathbb{A}_{z\mathcal{I},ij,-\frac{1}{2}}(\partial_{z}\bm{\mathbf{U}})_{ij,-\frac{1}{2}}-\sum_{n,l\pm}\mathbb{A}_{n,l\pm}(\partial_{n}\bm{\mathbf{U}})_{l\mp}-\bm{\mathbf{D}}_{ij0}}}^{\text{calculation}}$\hss}\underbrace{\mathbb{A}_{z\mathcal{I},ij,-\frac{1}{2}}(\partial_{z}\bm{\mathbf{U}})_{ij,-\frac{1}{2}}}_{\text{unknown}}-\underbrace{\sum_{n,l\pm}\mathbb{A}_{n,l\pm}(\partial_{n}\bm{\mathbf{U}})_{l\mp}-\bm{\mathbf{D}}_{ij0}-\overbrace{\overline{\dot{\bm{\mathbf{U}}}}}^{\text{estimate}}}_{\text{known}},\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}.$ (47)

Each term is labeled according to whether it is a known or unknown quantity,
and by whether it is part of the finite volume calculation or the
observational estimate. Our goal is to find the unknown quantities, that is,
the $k=-1$ cell centered state $\bm{\mathbf{U}}_{-1}$ and the $k=-\frac{1}{2}$
incoming subspace $\mathbb{A}_{\mathcal{I},-\frac{1}{2}}$, that together
minimize $\mathcal{R}$. Note that the first spatial derivative term on the
right hand side of the equation implicitly includes the unknown external state
via $(\partial_{z}\bm{\mathbf{U}})_{ij,-\frac{1}{2}}\approx\frac{1}{\Delta
z}(\bm{\mathbf{U}}_{ij0}-\bm{\mathbf{U}}_{ij,-1})$. Therefore, the problem is
generally nonlinear because the explicit $\bm{\mathbf{U}}_{-1}$ term from the
derivative multiplies the coefficient matrix
$\mathbb{A}_{\mathcal{I},-\frac{1}{2}}$, which itself depends on
$\bm{\mathbf{U}}_{-1}$ though the $k=-\frac{1}{2}$ interface state as defined
in (37).

An optimal solution can be found by minimizing $\mathcal{R}^{2}$ via a
nonlinear optimization algorithm such as Gauss-Newton, Levenberg-Marquardt,
conjugate gradient, etc. In that case, the optimal result is the state
$\bm{\mathbf{U}}_{-1}$ that combines with $\bm{\mathbf{U}}_{0}$ to
simultaneously generate $\mathbb{A}_{\mathcal{I},-\frac{1}{2}}$ and produce
$\partial_{t}\bm{\mathbf{U}}_{0}\rightarrow\overline{\dot{\bm{\mathbf{U}}}}$.
Instead, let us seek a solution for the $\mathcal{R}=0$ case by naively taking
the interface state as a given, inserting the finite difference version of the
vertical derivative at the $k=-\frac{1}{2}$ face, and formally solving for the
$\bm{\mathbf{U}}_{-1}$ by matrix inversion:

$\displaystyle\bm{\mathbf{U}}_{ij,-1}=\bm{\mathbf{U}}_{ij0}+\Delta
z\biggl{[}\mathbb{A}_{z\mathcal{I},ij,-\frac{1}{2}}\biggr{]}^{-1}\biggl{(}\overline{\dot{\bm{\mathbf{U}}}}_{ij0}+\sum_{n,l\pm}\mathbb{A}_{n,l\pm}(\partial_{n}\bm{\mathbf{U}})_{l\pm}+\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}.$ (48)

The incoming coefficient matrix is invertible (see (33c)), so provided it is
known (48) can be used to numerically calculate the state
$\bm{\mathbf{U}}_{ij,-1}$.

Equation (48) can in principle be used to find an iterative solution though
the following process. Starting with an initial guess for
$\mathbb{A}_{\mathcal{I},ij,-\frac{1}{2}}$, use (48) to solve for
$\bm{\mathbf{U}}{}_{ij,-1}$, use that state and $\bm{\mathbf{U}}{}_{ij0}$ to
calculate a new interface state, recalculate the eigendecomposition
$\mathbb{A}_{\mathcal{I},ij,-\frac{1}{2}}$, and then iterate. It can be shown
that this iterative state is stable, so that the correct solution reproduces
itself on iteration. However, convergence to the correct state is not
guaranteed.

We tested the iterative approach above by using the positive subspace of the
eigensystem at the $k=0$ cell center location in place of
$\mathbb{A}_{\mathcal{I},ij,-\frac{1}{2}}$ in (48) and then solving for the
$k=-1$ state. This simplification did work well for a variety of tests of our
data-driven boundary condition described in §7. A major drawback, however, was
that it could not be adapted to give different weights to estimates of
different primitive variables, and could not be extended to deal with missing
data; i.e., it was not a true optimization method.

A true optimization method, and the one we focus on for the remainder of this
paper, was implemented using a Singular Value Decomposition (SVD) approach. It
still only provides a solution to the linear problem by approximating the
incoming subspace at $k=-\frac{1}{2}$ using the values of primitive variables
at $k=0$; that is, we will set
$\mathbb{S}_{z\mathcal{I},-\frac{1}{2}}=\mathbb{S}_{z\mathcal{I},0}$ (or,
equivalently,
$\mathbb{A}_{z\mathcal{I},-\frac{1}{2}}=\mathbb{A}_{z\mathcal{I},0}$; see (42)
or (49)) and not consider additional updates to the $k=-\frac{1}{2}$ interface
state. This approximation can be relaxed by solving the nonlinear optimization
problem in future work.

To derive the SVD solution, the “unknown” term in (47) is transformed from the
space of primitive variable derivatives to the space of characteristic
derivatives via

$\displaystyle\mathbb{A}_{z\mathcal{I}}\cdot\partial_{z}\bm{\mathbf{U}}=\mathbb{S}_{z}\mathbb{\Lambda}_{z\mathcal{I}}\mathbb{S}^{-1}_{z}\cdot\partial_{z}\bm{\mathbf{U}}=\mathbb{S}_{z\mathcal{I}}\cdot\bm{\mathbf{L}}_{z\mathcal{I}},$
(49)

while the “known” term is collected into a single vector $\bm{\mathbf{f}}$:

$\bm{\mathbf{f}}=-\sum_{n,l\pm}\mathbb{A}_{n,l\pm}(\partial_{n}\bm{\mathbf{U}})_{l\mp}-\bm{\mathbf{D}}_{ij0}-\overline{\dot{\bm{\mathbf{U}}}},\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}.$ (50)