# The Equilibrium of Coronal Loops Near Separatrices

Emily I Mason Predictive Science Inc.
9990 Mesa Rim Rd, Suite 170
San Diego, CA 92121, USA  Spiro K Antiochos University of Michigan
500 S. State Street
Ann Arbor, MI 48109, USA  Stephen Bradshaw  Rice University
6100 Main St.
Houston, TX 77005-1827, USA

(Accepted December 15, 2022)

###### Abstract

We present numerical models from the field-aligned Hydrodynamics and Radiation
Code (HYDRAD) of a highly asymmetric closed coronal loop with near-singular
expansion factor. This loop was chosen to simulate a coronal magnetic flux
tube that passes close to a null point, as in the last set of closed loops
under the fan surface of a coronal jet or a pseudostreamer. The loop has a
very large cross-section localized near the coronal null. The coronal heating
was assumed to be uniform and steady. A siphon flow establishes itself within
4 hours of simulation time, flowing from the smaller-area footpoint to the
larger-area footpoint, with high initial speeds dropping rapidly as the plasma
approaches the null region. Observationally, this would translate to strong
upflows on the order of 10 km s-1 from the footpoint rooted in the localized
minority polarity, and weak downflows from the fan-surface footpoint on the
order of a few km s-1, along with near stationary plasma near the null region.
We present the model results for two heating rates. In addition, we analyzed
analogous Hinode EIS observations of null-point topologies, which show
associated Doppler shifts in the plasma that correlate well with the
simulation results in both direction and magnitude of the bulk velocity. We
discuss the implications of our results for determining observationally the
topology of the coronal magnetic field.

††journal: ApJ

##  1 Introduction

The structure and dynamics of the solar atmosphere are widely-believed to be
determined by two major effects: the magnetic field, which is due primarily to
the sources at the photosphere, and the coronal heating process that
determines the properties of the plasma contained in this field. Since the
plasma beta is measured to be low in the corona ($\beta<<1$), the field
magnitude is fairly uniform (Schrijver et al., 1999). The coronal heating, on
the other hand, must have considerable structure because it can only occur at
small scales for the high-Lundquist numbers of the corona. This in turn
structures the plasma, but only across field lines, because both mass
transport and thermal transport are highly efficient along the field. The net
result is that coronal plasma is distributed with roughly uniform temperature
and density along the field, but with clear variations across, so the corona
appears as a system of loops.

Observational analysis and numerical modeling of coronal loops have been
active fields of study for decades, since shortly after the first X-ray
observations by Skylab showed (apparently) monolithic tubes of hot plasma in
the low corona (e.g., Rosner et al. (1978)). These studies found that the loop
plasma achieves a quasi-steady equilibrium described fairly well by the
celebrated scaling laws (Craig et al., 1978; Rosner et al., 1978; Vesecky et
al., 1979), which predict that the loop plasma temperature and density are
approximately constant in the corona, and determined by only the heating rate
and the loop length, L. The expansion factor, $\Gamma$, measures the cross-
sectional variation of the loop from its chromospheric base to its apex; this
plays a more minor role on the plasma properties (Vesecky et al., 1979;
Cargill et al., 2021). For the most part plasma velocities are small,
$<10\,km\,s^{-1}$, due to asymmetries introduced by spatial variations of the
heating or by asymmetries in the loop geometry, i.e., spatial variations in
$\Gamma$. Strong dynamics in coronal loops are possible (i.e., if the heating
is strongly episodic as in certain nanoflare models Parker (1988); Klimchuk
(2006), or strongly localized as in the case of thermal non-equilibrium
Antiochos & Klimchuk (1991); Mason et al. (2019)). However, sophisticated
numerical models with a self-consistent treatment for wave heating show, for
the most part, agreement with the quasi-static scaling-law loops (Oran et al.,
2013; Lionello et al., 2009; Török et al., 2018).

Figure 1: (a) Signed $log(Q)$ (red and blue S-web map) over an example HMI
magnetogram. (b) diagram showing overall magnetic fields of a null-point
topology, with example flux tube outlined in blue. Credit: Peter MacNeice,
private communication. (c) diagram of selected field lines crossing a closed
polarity inversion line separating strong minority field from weak surrounding
field. The red arrow indicates the field line used to extrapolate flux tube
structure. Credit: Richard DeVore, private communication.

In recent years, it has been realized that the field in the corona must have
topological boundaries, since the magnetic flux distribution at the
photosphere is generally complex, exhibiting multiple polarity regions
separated by polarity inversion lines (PILs). An example of such a flux
distribution can be seen in Figure 1a, which shows the smoothed photospheric
radial field for Carrington Rotation 2099 in 2010 as measured by HMI on SDO,
(see hmi.stanford.edu/QMap/). The field strength is given by white (positive)
and black (negative). We note the presence of several active regions,
appearing as bright black and white bipoles, during this CR. It should be
emphasized that there are very many more bipoles on the surface, the so-called
magnetic carpet, but these are small-scale and are averaged out in this
smoothed magnetogram.

Even the large flux concentrations, however, produce considerable structure at
the photosphere, which gives rise to a coronal magnetic field with complex
topology. For the corona, the topology is given primarily by the connectivity
of the field, the connections between the various flux systems. In general,
the coronal field is not observed to have higher-order topological features
such as knots and disconnected rings. Figure 1a also shows the connectivity of
the magnetic field as given by the so-called squashing factor, Q (Titov,
2007). Q is a measure of the gradients in the mapping defined by the magnetic
field lines. In the particular case of Figure 1a, the mapping is from the
spherical surface at 1.1 RS down to the photosphere, and the contours shown
are for large values of $log(Q)$, colored by the polarity of the field, red
(negative) and blue (positive). The intricate web of strong Q contours shown
in the figure is referred to as the S-Web (Antiochos et al., 2011). This S-Web
indicates locations where the magnetic connectivity is discontinuous, as at a
flux boundary, or near-discontinuous, as in the case of so-called quasi-
separatrix layers, (e.g., Demoulin et al., 1996).

There are two key points to the S-Web shown in Fig. 1a. First, wherever there
is an arc of large Q, electric currents can form easily, leading to magnetic
reconnection and eruptive solar activity. This result is the physical
underpinning for essentially all models of coronal jets and flares (Shibata et
al., 1994; Demoulin et al., 1996; Antiochos et al., 1999; Wyper et al., 2018).
The second key point is that unlike the actual magnetic flux at the
photosphere, the Q-map shown in the figure is derived from a model – in
particular, a potential field source surface model. Given that the coronal
field is never potential, (there are always filament channels and
sheared/twisted structures present), we conclude that the locations of all the
arcs and lines in Fig. 1a are likely to be incorrect. It would be extremely
valuable, therefore, to have some method for determining observationally the
locations of the separatrix curves as in Fig. 1a. We propose a possible method
in this work.

The S-Web of Fig. 1a consists of a variety of separatrix and quasi-separatrix
curves, but the simplest such curve (and generally the most common) is that
due to a single polarity region surrounded by an opposite-polarity background,
giving rise to the so-called embedded bipole topology as seen in Figure 1b
(Antiochos, 1998; Pariat et al., 2009). This topology consists of a dome-like
separatrix surface in the corona, the fan surface that separates the flux of
the embedded bipole from surrounding flux. Somewhere on the fan the field
vanishes to form a 3D null point, and two singular field lines – the “spines”
– emanate from the null (e.g., Shibata et al., 1994; Young & Muglach, 2014;
Karpen et al., 2017; Wyper et al., 2018). Such a magnetic topology occurs for
a single polarity region that has fairly circular geometry, but the actual
photosphere contains many polarity regions, sometimes with nested embedded
polarities and/or with extremely structured shapes. As a result, the magnetic
topology in the corona consists of a complex web of interconnected
separatrices and numerous null points/spines, as implied by the figure.

The S-Web has been postulated to have major effects on the structure and
dynamics of the solar wind (Antiochos et al., 2011; Titov et al., 2011; Linker
et al., 2011; Higginson et al., 2017), especially the slow wind, but has
generally not been considered in modeling the closed field plasma. The primary
reason for this is that the topology of the coronal field is very difficult to
observe directly in XUV or X-ray images. Only certain topological features,
such as domes and spines have been observed, and only on the limb (e.g., Mason
et al., 2021). Note also that the S-Web and topological features are generally
global quantities and cannot be determined from local measurements of the
magnetic field, even if such were available. Consequently, the existence of
the S-web itself is currently only inferred from theoretical models.

In this paper we propose using loop models to determine observational proxies
of coronal topology, specifically, the connectivity of the magnetic field. The
key point for our studies is that flux tubes near the features of the S-Web
must exhibit extreme variations in their geometry (such as that indicated by
the highlighted field line in Figure 1c), because they must pass close to a
null point, and consequently, will have near-singular expansion factor
$\Gamma$ there. Such variations are likely to produce observable signatures in
the coronal plasma, such as rapid flows from one end to the other. By
determining the characteristics expected for these signatures, it should be
possible to map out the S-Web using high-resolution spectroscopic observations
expected from missions such as the upcoming EUVST (Shimizu et al., 2019) and
MUSE (De Pontieu et al., 2022), for example.

Furthermore, investigating flux tubes near topological singularities may well
result into new insights into coronal loops, themselves. It is not clear that
equilibrium states are even possible for extreme cross-section variation, even
if the heating is spatially and temporally constant. It is well known that
extreme localization of the heating process can lead to a lack of equilibrium,
a so-called thermal non-equilibrium (Antiochos & Klimchuk, 1991); the same may
occur if the expansion is extremely localized.

For this first study of the properties of S-Web plasma, we adopt the classic
assumption of a uniform and constant coronal heating and consider the most
common near-singular magnetic structure, a flux tube that lies close to the
separatrix surface of an embedded bipole and, therefore, passes very near to a
coronal null. Such flux tubes must undergo a dramatic, localized increase in
cross-sectional area, in order to satisfy the conservation of magnetic flux.
These flux tubes will be asymmetric in two ways: the aforementioned varying
cross-sectional area, and a left-to-right asymmetry due to the fact that one
leg of the flux tube lies near the inner spine (associated with stronger
field) whereas the other leg lies near the (weaker-field) fan surface (Mason
et al., 2019). This results in a strong global area factor variation in
addition to the localized maximum near the null. This type of coronal loop,
modeled here, can be seen in blue in Figure 1b.

In the following section we describe the field-aligned hydrodynamic code,
HYDRAD, that we use and provide the details of the loop model: geometry,
heating rate, etc. Section 3 presents the results of the simulation, including
an analysis of expected observational signatures. For comparison with these
results, the next section shows Hinode EUV Imaging Spectrometer (EIS, Kosugi
et al. (2007); Culhane et al. (2007)) observations of similar loops in several
medium-sized embedded bipoles like those discussed above, called null-point
topologies (NPTs) near disk center. The final section discusses the
implications of this study for observational determination of the S-web in
closed field regions, and for the fundamental physics of coronal loops.
Furthermore, we discuss the likely effects of important extensions to the
model, such as a spatially or temporally varying heating and effects of
magnetic field dynamics.

##  2 HYDRAD

###  2.1 Equations

For the simulations presented here, we use HYDRAD (Bradshaw & Mason, 2003;
Bradshaw & Cargill, 2013), a field-aligned, adaptive mesh, hydrodynamic code
that can accommodate arbitrary loop geometries via a polynomial fit to the
cross-sectional area factor and parallel component of gravity along the field.
HYDRAD is well-documented and has been widely used for calculating coronal
loop plasma dynamics in a wide variety of situations, including nanoflare
heating and thermal nonequilibrium (Bradshaw et al., 2004; Bradshaw & Cargill,
2005; Bradshaw et al., 2011; Cargill et al., 2012; Brooks et al., 2016; Reep &
Russell, 2016; Johnston et al., 2019). We do not utilize the partial
ionization state capabilities of the code, since these generally are more
relevant in impulsive heating scenarios. The loops we model here have only
steady background heating, rendering the special ionization rates less likely
to be necessary. It should be noted, however, that a fast steady flow through
the transition region can affect the radiative losses in that region, but the
qualitative structure remains the same (e.g., Spadaro et al. (1990)).

HYDRAD solves the conservative form of the standard hydrodynamic loop
equations, excerpted partially here from Bradshaw & Mason (2003); Bradshaw &
Cargill (2013):

$\frac{\partial\rho}{\partial t}+\frac{1}{A(s)}\frac{\partial}{\partial
s}(A(s)\rho v)=0$ (1)

$\frac{\partial}{\partial t}(\rho v)+\frac{1}{A(s)}\frac{\partial}{\partial
s}(A(s)\rho v^{2})=\rho g_{||}-\frac{\partial P}{\partial
s}+\frac{\partial}{\partial s}(\frac{4}{3}mu_{i}\frac{\partial v}{\partial
s})$ (2)

$\frac{\partial E_{e}}{\partial t}+\frac{1}{A(s)}\frac{\partial}{\partial
s}[(E_{e}+P_{e})A(s)v]=\frac{1}{A(s)}\frac{\partial}{\partial
s}(A(s)\kappa_{0e}T^{\frac{5}{2}}_{e}\frac{\partial T_{e}}{\partial
s})+v\frac{\partial P_{e}}{\partial
s}+\frac{k_{B}n}{\gamma-1}\nu_{ie}(T_{i}-T_{e})-R+H$ (3)

$\frac{\partial E_{i}}{\partial t}+\frac{1}{A(s)}\frac{\partial}{\partial
s}[(E_{i}+P_{i})A(s)v]=\frac{1}{A(s)}\frac{\partial}{\partial
s}(A(s)\kappa_{0i}T^{\frac{5}{2}}_{i}\frac{\partial T_{i}}{\partial
s})+v\frac{\partial P_{i}}{\partial
s}+\frac{k_{B}n}{\gamma-1}\nu_{ie}(T_{e}-T_{i})+\frac{\partial}{\partial
s}(\frac{4}{3}mu_{i}v\frac{\partial v}{\partial s})+\rho vg_{||}$ (4)

$E_{e}=\frac{P_{e}}{\gamma-1},E_{i}=\frac{P_{i}}{\gamma-1}+\frac{1}{2}\rho
v^{2}$ (5)

$P_{e,i}=2k_{B}nT_{e,i}$ (6)

where $\rho$ is the plasma density, $A(s)$ is the cross-sectional area, $s$ is
the dimension along the loop, $v$ is velocity in the field-aligned direction,
$g_{||}$ is the component of gravity along the field, $P$ is pressure, $E$ is
energy, $\kappa$ is the conductivity, $T$ is temperature, $H$ is the heating
function term, $R$ is the radiative loss function, $\nu_{ie}$ is the Coulomb
collision frequency, $u_{i}$ is the ion viscosity, and $n$ is the particle
density. HYDRAD solves the electron and ion temperatures separately, hence the
subscripts on energy, pressure, and temperature above, and the presence of two
energy equations; however, these deviate appreciably only in highly impulsive,
low-density situations. For the simulations presented below, the deviations
between the electron and ion temperatures remained negligible.

###  2.2 Geometry & Boundary Conditions

To obtain the magnetic structure of our near-null loop we use the well-known
embedded bipole field consisting of a uniform background vertical field and a
vertically-oriented point dipole located some distance d below the
photospheric surface. This produces an axisymmetric system and has been used
by many authors for studies of coronal jets (e.g., Pariat et al. (2009); Wyper
et al. (2018)). Note that the embedded bipole field is the simplest possible
topology that exhibits a separatrix and is very frequently present in the
actual corona. In cylindrical coordinates the field takes the simple form:

$\vec{B}_{total}=\vec{B}_{0}+\vec{B}_{dipole}$ (7)

where

$\vec{B}_{0}=B_{0}\hat{z}$ (8)

and

$\vec{B}_{dipole}=\frac{\mu_{0}d}{4\pi}(\rho^{2}+(z+d)^{2})^{-\frac{3}{2}}[3\rho(z+d)\hat{\rho}+(2(z+d)^{2}-\rho^{2})\hat{z}]$
(9)

The resulting field is shown in Figure 1, exhibiting a null point along with a
dome separatrix and spine lines, as expected. In the corona structures of this
form have long been observed and identified as jets (Bohlin et al., 1975),
anemones (Shibata et al., 1994), or fountains (Tousey et al., 1973), depending
upon the scale size at which they form. Several of the closed field lines
under the fan surface are shown in Figure 1b. To define the loop geometry for
the simulations presented here, we extracted the coordinates and arrays of
$g_{||}$ and $B$ values of a single magnetic field line from a coronal jet
simulation utilizing the field above (Karpen et al., 2017).

From the embedded bipole field above, and shown in Figure 1, we selected a
field line that passes very near to the coronal null seen in Figure 1b. Figure
2a shows the cross-sectional area profile for the loop, and Figure 2b shows
the parallel gravity ($g_{||}$) for both the original loop and the polynomial
fit which is the input for HYDRAD. The gravity profile is not significantly
different than the coronal loops of the many previous studies, and while the
polynomial fit is not as accurate as the magnetic field fit for the same loop,
the relatively short length of the loop precludes significant impact to the
simulation. The area factor, on the other hand, shows some important new
features. First, we note that the left footpoint of the loop is located in a
strong field region, while the right footpoint is located in field about one-
third the strength and, consequently, the area on the right is approximately 3
times that on the left. More important, we note that there is a region of very
strong expansion localized near $s=28$ Mm. The maximum area there is roughly a
factor of 40 greater than that at the left footpoint. We applied an
$8^{th}$-order polynomial fit to the magnetic field to model the loop geometry
in HYDRAD. The maximum standard fit in HYDRAD is $6^{th}$-order; however, this
did not sufficiently capture the sharp apex of the loop, requiring us to
increase the order. This again emphasizes the unique properties of near-null
loops. The polynomial equations utilized for the loop in the simulations are

$B(s)=7802s^{8}-34,218s^{7}+62,460s^{6}-61,617s^{5}+35,866s^{4}-12,786s^{3}+2881s^{2}-417s+32.37$
(10)

$g_{||}=-(6.958\cdot 10^{6})s^{6}+(2.081\cdot 10^{7})s^{5}-(2.288\cdot
10^{7})s^{4}+(1.100\cdot 10^{7})s^{3}-(2.014\cdot 10^{6})s^{2}+(1.127\cdot
10^{5})s-2.771\cdot 10^{4}$ (11)

The maximum refinement level in an adaptive grid simulation determines the
size of the smallest grid cell. For the results presented here, a maximum
refinement level of 10 was selected, resulting in a minimum grid cell size of
$1.22\cdot 10^{5}$ cm. We additionally constrained the simulation to have a
maximum grid cell width of $1.25\cdot 10^{8}$ cm. Note that both footpoints
are several scale heights deep in the chromosphere, so that there is
essentially no evolution there.

Figure 2: (a) graph showing the variation of the cross-sectional area along
the length of the loop used in the HYDRAD code. The specific values at the
footpoints and near the null are shown on the graph. (b) graph showing the
variation of the parallel component of gravity along the length of the
original loop (in black) and along the length of the loop used in the HYDRAD
code (blue).

##  3 Simulation Results

We executed two steady-state simulations, differentiated only by disparate
background heating rates. The background heating rate in HYDRAD is determined
by the initial conditions calculation, for which the user inputs the loop
length, the distance between the tops of the constant-temperature
chromospheric regions at each end (i.e., the coronal loop length), and the
plasma density at these points. The top of the chromosphere, where the
temperature gradient vanishes, is generally defined to be the loop “base” in
coronal loop models. The initial-conditions code then computes the steady
background heating rate which would be necessary to sustain a symmetric,
uniform cross-section loop in hydrostatic equilibrium with the given footpoint
density, (see (Bradshaw et al., 2004, 2011) for full details). We then apply
the geometry seen in Figure 2 and allow the loop to find a new steady-state
equilibrium, if one exists, with the previously-determined background heating.

For the simulation that we hereafter refer to as low-heat (LH), the code’s
standard base (top of the chromosphere) plasma density of $10^{11}\,cm^{-3}$
resulted in a background heating rate of $4.882\cdot
10^{-4}\,erg\,s^{-1}\,cm^{-3}$. The high-heat (HH) run was set to a base
density of $4\cdot 10^{11}\,cm^{-3}$, resulting in a background heating rate
of $2.417\cdot 10^{-3}\,erg\,s^{-1}\,cm^{-3}$. Both of these values are
typical for coronal loops.

Figure 3 shows the final velocities that we find from the two simulations, and
is the primary result of our calculations. First, and most important, we find
that these velocities are essentially steady. Both runs reached their steady
state quickly after the beginning of the simulation; the LH steady state was
established around 2.5 hours of elapsed solar time, while the HH converged
even faster, after only 2.2 hours. We do find some “jitter” in the solution
due to numerical noise and waves bounding back and forth, but these time-
average out to zero. Therefore, even for this extreme variation in the loop
cross-section, a steady solution is possible – at least for a near-uniform
heating. The velocity profile, however, is distinctly different than that of
previous loop models due to the unique nature of the area variation.

The effect of a cross-sectional area variation $A(s)$ on loop plasma structure
has been studied for many years (e.g., Vesecky et al., 1979; Mikić et al.,
2013a). In particular, recent studies by Klimchuk & Luna (2019) and Cargill et
al. (2022) present a thorough analysis on the effect of an area variation on
loop plasma structure and dynamics, including detailed numerical simulations
for spatially- and/or temporally-varying coronal heating. Here we only
summarize the underlying physical understanding, but encourage the reader to
refer to those papers for the rigorous modeling. Consider the classic loop
model that has a uniform cross-section, $A(s)\equiv 1$, and is symmetric about
the apex. Since the heat flux at the loop base vanishes, the total heating
into the loop, $E$, must be balanced by the radiation losses from the coronal
and transition-region sections of the loop $E=R_{c}+R_{tr}$. It has been
established from many simulations that the losses from these two regions are
of the same order of magnitude, but with the transition region emission
somewhat larger, $R_{tr}\approx 2-3\times R_{c}$ (e.g., Klimchuk et al.,
2008). Note that for this loop the average area in the corona is equal to that
in the transition region, since the area is constant everywhere. Now consider
the effect of an area variation $A(s)\neq const$, as for example, in a typical
dipole magnetic field, and assume that the area at the loop base $A_{b}=0.1$
whereas that at the apex remains at $A_{c}=1$. For this loop the area factor
$\Gamma=10$, which is a typical value. The ratio of the radiation losses from
the transition region versus that from the coronal section of the loop now
changes due to the relative average area of these sections. Since the
transition region lies at the very base of the loop and is extremely narrow,
its average area $\langle A_{tr}\rangle\approx A_{b}$. As a result, the losses
from this region are reduced by this factor, $R_{tr}\times A_{b}$, whereas the
coronal losses and total heat input are reduced by the value of the area
factor averaged along the whole loop $\langle A(s)\rangle$. But $\langle
A(s)\rangle>A_{b}$ so as a result, the loop is no longer in energy balance;
the heating is larger than the radiative cooling,
$E>R_{c}+R_{tr}\times\frac{A_{b}}{\langle A(s)\rangle}$. To establish a new
balance the coronal losses must increase from their previous value, which
requires that the density and hence the pressure increases. We conclude that
an area constriction causes the coronal pressure to increase (see (Klimchuk &
Luna, 2019; Cargill et al., 2022) for quantitative details on how this
increase varies with $A(s)$). If the area factor (effective constriction) is
larger in the left leg than in the right leg, then the pressure will be larger
on the left so that no static equilibrium is possible and a flow from left to
right must occur. For example, if the left leg occurs in a high-field strength
$B\approx 100$ G active region, whereas the right leg is in the surrounding
photosphere where typically, $B\approx 10$ G, then an area expansion of 10
will be present from left to right. This will result in a relatively-slow,
$V<10$ km/s, steady siphon flow (e.g. Antiochos, 1984) that accelerates from
the left leg, reaches a maximum in the corona and decelerates on the right.

Figure 3: Time-averaged velocity graphs for LH (blue) and HH (red) steady
state runs. The LH steady state velocity peaks on the left leg at 9.78 km s-1,
while the HH run produced a maximum velocity near the same point at 11.22 km
s-1. Both runs experienced near-null velocities of less than 1 km s-1, and the
right leg maxima were 3.43 and 3.88 km s-1 for the LH and HH simulations,
respectively.

These straightforward arguments describe the physical origin of the flows seen
in Fig. 3. Defining each loop leg from the apex, which coincides very closely
with the point of maximum area, then the left leg has an effective $\Gamma$
factor much larger than that on the right, which accounts for the left-to-
right flow. Such a flow would inevitably be present in any loop that lies very
near a separatrix boundary, and would be relatively large compared to loops
far from the separatrix, because the area variation is much larger at the
separatrix, becoming near-singular there.

A unique feature of the flow is that it becomes very small near the apex
(vicinity of null). The reason for this is straightforward. Since the loop is
in a quasi-steady state, the flow must be steady so that the mass flux $nvA$
is constant. As in all coronal loop models, thermal conduction keeps the
temperature in the coronal section of the loop near-constant; consequently,
the density must also be nearly constant in the corona, otherwise there would
be large pressure imbalances. The area, however, becomes very large near the
null, so the velocity must become small there in order to maintain a near-
constant mass flux. It should be emphasized that such a velocity profile would
be a feature of only a separatrix loop, because it requires a very large and
very localized expansion factor; in other words, a null point region.

In terms of the quantitative values of the velocities for our particular
parameters, the cases have different velocities, as one would expect, given
the heating differential. There is a 1.5 km s-1 difference in maximum velocity
on the left leg (9.78 and 11.22 km s-1 for LH and HH equilibria runs,
respectively); the velocities near the null are similar, and both are
approximately 0.5 km s${}^{-}1$ (0.43 and 0.56 km s-1), in essence, negligible
velocity. On the right leg, the velocity picks up again, with the peaks at
3.43 and 3.88 km s-1, the higher velocity belonging to the HH simulation.
These values would be easily observable by spectroscopic instruments, and
would be present in all coronal lines.

It should be emphasized that the flows in Figure 3 are all due to plasma
pressure gradients; there is no other force in this 1D system other than
gravity, which is always downward. Figure 4a shows the electron temperatures
and pressures for both runs. The peak loop temperature occurs to the left of
the near-null region on both runs, and there is approximately a 1.5 MK
difference between the temperatures (2.43 and 3.89 MK) in the two cases. Note
that both cases show a small asymmetry in the density, with slightly higher
density and a more gradual gradient on the right leg compared to the left (see
Figure 4b). These gradients are required to balance the flow profile, Fig. 3,
which shows a fast acceleration on the left, but a gradual deceleration on the
right. The pressure exhibits a rapid decrease on the left leg, which
accelerates the flow upward, then bottoms out to a near constant level, except
for the small variations required to balance gravity, and then increases on
the right leg in order to decelerate the flow. Again, these gradients are set
up by the need to maintain local thermal balance (heat in equals radiation
out), given that the transition region radiation from the two legs is very
different.

Figure 4: (a) Graph of the electron temperatures (solid lines) and pressures
(dashed lines) for both the LH (blue) and HH (red) equilibria runs. The LH run
reaches a maximum temperature of 2.43 MK, while the HH peaks at 3.89 MK. There
is an over four-fold increase in electron pressure between the two runs, at
0.50 and 2.39 for the LH and HH runs. (b) Electron density for the two runs.

The HH run’s minimum, coronal pressure (2.39 dyne cm-2) is over 4 times higher
than the LH run’s (0.50 dyne cm-2). This is quite close to what would be
predicted by the standard scaling laws: $P\propto H^{6/7}$ (e.g., Rosner et
al., 1978; Bradshaw et al., 2019). This result is to be expected given the
arguments above. The scaling laws basically reflect the requirement that the
coronal radiative losses must approximately balance the heating there, given
that the transition region does not dominate the losses (Rosner et al., 1978;
Bradshaw et al., 2019). This requirement is even more valid in our case,
because the transition region radiation is suppressed by the large expansion
factors. Although there are strong flows in our model, they act only to
redistribute the energy and are not an energy sink, because just like the
conductive fluxes, they vanish at both footpoints. Note also that the flows
become very small near the null region, so their enthalpy flux is negligible
there. This has important implications for the thermal stability of the loop
plasma, as will be discussed in Section 5 below.

##  4 Hinode Observations

The main conclusion from the results above is that coronal loops near
separatrix surfaces should have distinctive velocity signatures. As a
preliminary observational test of this hypothesis, we examined a well-observed
and clearly identified null-point topology region for the type of flows
predicted by our models. One caveat to our analysis – and essentially to all
observational inferences of coronal loop structure – is that there appears to
be a long-standing contradiction between the theoretical picture of loops as a
representation of the variable magnetic field of a coronal flux tube, and
observations that apparently show structures of nearly-constant cross-
sectional area (Klimchuk et al., 1992; Peter & Bingert, 2012; Klimchuk &
DeForest, 2020). Most loop simulations, like this one, and many others (e.g.,
Mikić et al., 2013b; Froment et al., 2018), consider loops with varying cross-
sectional area, as required by the conservation of magnetic flux. However,
observations by some of the same authors, including Klimchuk & DeForest (2020)
and Auchère et al. (2018), consistently show loops with cross-sectional areas
that vary very little, even in the vicinity of higher-beta regions. An
alternative approach to apparent loop geometries is presented by Malanushenko
et al. (2022), but that paper focuses mostly on more traditional active region
loops where the low-beta approximation is valid. As has been discussed in
Mason et al. (2019), which first reported the NPTs presented here, there are
frequently distinct dark regions around the null of such structures. It is
likely that these are caused by highly-expanded loop structures like the one
that we model here. It is difficult, however, to understand how a loop that
appears to have constant cross-section could sustain the pressure differential
that would attend the required variance in plasma pressure along its length,
particularly with a flow present. Without reliable estimates of the magnetic
field strength in these low-coronal regions, the mystery persists. In order to
resolve the loop-area contradiction, it is clear that much more work is
required, both modeling with a 3D MHD code with full thermodynamics as in
Mikić et al. (2018), and direct spectropolarimetric measurements of the
magnetic fields in and around null-point structures as, perhaps, with the
Daniel K. Inouye Solar Telescope.

Ignoring the contradiction above, we now present observations of several NPTs
imaged by the Hinode EIS instrument, to investigate the evidence for long-
lived end-to-end flows (Klimchuk & Luna, 2019) and thereby test our
conclusions above. The two structures shown here have already been
investigated in detail in Mason et al. (2019) and Mason et al. (2021), and the
reader is referred to those works for details on their evolution and dynamics;
here we simply state that they are NPTs formed from decaying active regions
which form on the border of (or entirely within) coronal holes.

Figure 5 shows two NPTs as seen in the Fe XII 195.12 Å wavelength; the peak
formation temperature for this line is 1.58 MK. In both 5a and 5b, there are
distinct regions around the inner spine of the NPTs that have upflows on the
order of 3-6 km s-1. Note that for a complex polarity distribution such as
found on the true photosphere – especially those that give rise to so-called
pseudostreamers – the “spine” often corresponds to a sheet-like structure
(e.g. Titov et al., 2011; Scott et al., 2021). On the opposite end of these
closed loops, around the fan surface, there are strong downflow signatures,
which peak between 10 (Figure 5b) and 16 km s-1 (Figure 5a). This general
pattern correlates with the simulation; the positive flows in the direction of
the loop on both legs correspond to Doppler shifts in the observations that
show upflows on the spine leg and downflows on the fan leg. Interestingly, the
magnitude of the upflows on the spine legs in the observations is consistently
lower than the downflows on the fan leg, in contradiction to the simulation
results. This could be due to observational affects, such as the inner spine
being at an angle from the perpendicular with respect to the instrument,
thereby lowering the vertical component of the upflow. A key point is that
overlying coronal loops that have a strong overall net flow could weaken the
upflow and strengthen the downflow signatures. Both NPTs discussed here are
adjacent to a helmet streamer, and it has been argued recently that such flows
may be common there due to the effects of thermal nonequilibrium (Schlenker et
al., 2021). Alternatively, the dynamics of the near-null expansion in creating
the strong upflow along the spine leg could be somewhat less exaggerated on
the Sun due to effects such as field-line tangling or impulsive heating that
are not included in the calculations above. Further study is required to
determine other factors that might be affecting the distinction in velocities
between the two legs, and whether this is a consistent discrepancy across NPT
closed-field loops. The basic flow pattern, however, agrees very well with our
loop modeling and demonstrates the distinctive nature of flows near
separatrices.

Figure 5: (a) This column presents the null-point topology observed by Hinode
EIS on 2015 April 29. The top image is a 171-193-211 Å composite AIA image for
context, overlaid with a potential field source surface, from The Sun Today.
The white rectangle shows the FOV for the EIS observation at the bottom, and
the yellow square shows the FOV of the SDO HMI and AIA cutouts immediately
below. The AIA and HMI cutouts are overlaid with the approximate location of
the polarity inversion line within the NPT. At the bottom are the Hinode EIS
observations using the 195.12 Å wavelength. Closed loops can be seen in the
intensity map, running from several inner spines along the long axis of the
ovate dome to the outer fan surface. The velocity map shows several locations,
marked with black arrows, where the flows are upward from the several inner
spines; the entire outline of the fan surface, however, is marked with strong
downflows. (b) Analogous observation of an NPT on 2016 February 20; locations
of upflows are again marked with black arrows on the velocity map, and
velocities are similar in range.

##  5 Discussion

The results above have a number of important implications for understanding
the structure and dynamics of coronal plasma. One key conclusion is that flux
tubes (coronal loops) near a separatrix in the closed corona should exhibit
strong persistent flows, with a direction determined by the topology of the
separatrix, spine vs fan. One caveat to this conclusion that merits some
discussion is the validity of the 1D model for near-null loops. The underlying
assumption in all coronal loop models is that the plasma beta is so low that
the field can be assumed to be rigid throughout the plasma evolution. But if a
loop has a very large expansion at some location, it implies that the field
strength becomes small there, so the low beta assumption becomes suspect. In
our loop, for example, the field strength is down by a factor of roughly 40
from its value at the chromospheric base, so the field energy density
decreases by almost three orders of magnitude over the same distance. We
conclude that simulations of observed loops near separatrices should be
checked a posteriori to verify that the low-beta assumption does remain valid.
If the beta does become significantly larger than unity, the plasma pressure
is expected to cause the field near the null to expand even further, thereby
enhancing the flows found above. It is also possible that such flows induce
interchange reconnection, which would induce additional dynamics as found by
Scott et al. (2021), for example. Furthermore, magnetic field changes could
well lead to non-steady flows and oscillations in the field. However, the
thermal conduction and mass flows will still be along the field-aligned
direction, which may maintain similar dynamics to truly low-beta regions. We
conclude, therefore, that MHD modeling of near-null loops would be highly
informative, and plan to undertake such a study in a subsequent paper.

For field-aligned models, our work proves the viability of an end-to-end flow
as a stable steady state solution in coronal loops with extreme geometry. We
find that such a flow establishes itself rapidly as a stable steady state for
a range of heating values in a loop that is highly asymmetric both in its apex
location and cross-sectional area variation. The flow is from left to right,
as expected from simple considerations of loops with asymmetric expansion
factor. It is well known that if everything else is kept constant, increasing
the expansion factor increases the coronal pressure, because the relative
amount of transition region emission decreases (Klimchuk et al., 2008; Cargill
et al., 2021). Since the left leg of our loop has a larger corona-to-footpoint
expansion factor, it has a larger pressure, which drives a flow from left to
right. This finding supports the many observations of flows in coronal loops
in structures like active regions, observed even when there is no apparent
time-dependent heating occurring. Previous work has suggested that such
apparent motion may be due to waves propagating along the loop instead of
actual plasma motion, but this work shows that stable end-to-end flows can be
maintained in a broad range of loop geometries and situations.

Although we do find stable steady state loops, it should be emphasized that
they are likely to be more susceptible to thermal instability than standard
loops. The key point is that in the vicinity of the null, the energy balance
is primarily between heating and radiative losses. As discussed above, the
transition region losses are highly suppressed by the large expansion factors.
Also, the large expansion causes the flows and, hence, enthalpy flux to be
small there. If the thermal and enthalpy fluxes were completely negligible, so
that the energy balance were only between a constant heating and coronal
radiation, then the plasma there would be more likely to be thermally unstable
(e.g. Parker, 1953; Field, 1965). We conclude that the small residual fluxes
are providing the stability in our solutions. The implication, however, is
that as one considers loops that are closer and closer to the separatrix and,
therefore, have an increasingly larger localized expansion, then the plasma
there would become thermally unstable. In fact, we have observed coronal rain
that we attributed to thermal non-equilibrium preferentially near the
separatrix surfaces of null-point topologies (Mason et al., 2019). We intend
to investigate this issue in detail with further modeling work.

The most important result from our work is that the topology of the coronal
magnetic field, in particular the imprint of the S-Web at the top of the
corona, should be observable with high-resolution spectroscopic measurements.
Although we have used the simplest possible coronal heating, constant in space
and time, our results must hold in general: fast flows (exceeding 10 km s-1)
will be induced by the extreme flux tube geometry at separatrix boundaries. In
fact, the data presented in Section 4 appear to support this conclusion. Note
that as shown above, the magnitude of the steady flows is sensitive to the
magnitude of the heating and the loop geometry. Furthermore, the flows will be
sensitive to any spatial variation in the heating, for example, if the strong
field region on the left has higher heating than the right. Precise
measurements of the plasma velocities would provide stringent constraints on
such variations. Note also that although our single-null-point loop is likely
to be the most common separatrix structure, it is not unlikely to have flux
tubes that pass by two null points as in models of large pseudostreamers (e.g.
Titov et al., 2011; Scott et al., 2021). Such loops would have multiple
localized high-expansion regions with even more complex dynamics than our case
above, and may well be even more sensitive to the details of the heating. The
wave heating models, in particular, are sensitive to the magnetic field
variations along the loop, because the heating depends strongly on wave
reflection from variations in the Alfvén speed along the loop (e.g., van der
Holst et al., 2014; Downs et al., 2016). Near-null loops would have extreme
variations in the Alfvén speed. Our results and possible extensions,
therefore, constitute a potentially powerful observational probe of not only
coronal magnetic topology, but perhaps even of the heating mechanism itself.
It is clear that a great deal more theoretical and observational work remains
to be done on this important topic of coronal loops near separatrices.

##  6 Acknowledgements

EIM’s research during the development of this paper was supported by an
appointment to the NASA Postdoctoral Program at the NASA Goddard Space Flight
Center, administered by Universities Space Research Association under contract
with NASA. EIM would like to acknowledge support from the NASA LWS program.
SKA was supported by LWS Strategic Capability Grant 80NSSC22K0892 to the
University of Michigan.

## References

  * Antiochos (1984) Antiochos, S. K. 1984, ApJ, 280, 416, doi: 10.1086/162007
  * Antiochos (1998) Antiochos, S. K. 1998, The Astrophysical Journal, 502, L181, doi: 10.1086/311507
  * Antiochos et al. (1999) Antiochos, S. K., DeVore, C. R., & Klimchuk, J. A. 1999, ApJ, 510, 485, doi: 10.1086/306563
  * Antiochos & Klimchuk (1991) Antiochos, S. K., & Klimchuk, J. A. 1991, The Astrophysical Journal, 378, 372, doi: 10.1086/170437
  * Antiochos et al. (2011) Antiochos, S. K., Mikić, Z., Titov, V. S., Lionello, R., & Linker, J. A. 2011, Astrophysical Journal, 731, doi: 10.1088/0004-637X/731/2/112
  * Auchère et al. (2018) Auchère, F., Froment, C., Soubrié, E., et al. 2018, ApJ, 853, 176, doi: 10.3847/1538-4357/aaa5a3
  * Bohlin et al. (1975) Bohlin, J. D., Vogel, S. N., Purcell, J. D., et al. 1975, The Astrophysical Journal, 197, L133, doi: 10.1086/181794
  * Bradshaw et al. (2011) Bradshaw, S. J., Aulanier, G., & Del Zanna, G. 2011, Astrophysical Journal, 743, doi: 10.1088/0004-637X/743/1/66
  * Bradshaw & Cargill (2005) Bradshaw, S. J., & Cargill, P. J. 2005, Astronomy & Astrophysics, 437, 311, doi: 10.1051/0004-6361:20042405
  * Bradshaw & Cargill (2013) —. 2013, Astrophysical Journal, 770, doi: 10.1088/0004-637X/770/1/12
  * Bradshaw et al. (2004) Bradshaw, S. J., Del Zanna, G., & Mason, H. E. 2004, Astronomy and Astrophysics, 425, 287, doi: 10.1051/0004-6361:20040521
  * Bradshaw et al. (2019) Bradshaw, S. J., Emslie, A. G., Bian, N. H., & Kontar, E. P. 2019, ApJ, 880, 80, doi: 10.3847/1538-4357/ab287f
  * Bradshaw & Mason (2003) Bradshaw, S. J., & Mason, H. E. 2003, Astronomy & Astrophysics, 401, 699, doi: 10.1051/0004-6361:20030089
  * Brooks et al. (2016) Brooks, D. H., Reep, J. W., & Warren, H. P. 2016, The Astrophysical Journal, 826, L18, doi: 10.3847/2041-8205/826/2/l18
  * Cargill et al. (2012) Cargill, P., Bradshaw, S., & Klimchuk, J. 2012, Astrophysical Journal, 752, doi: 10.1088/0004-637X/752/2/161
  * Cargill et al. (2021) Cargill, P. J., Bradshaw, S. J., Klimchuk, J. A., & Barnes, W. T. 2021, Monthly Notices of the Royal Astronomical Society, 509, 4420, doi: 10.1093/mnras/stab3163
  * Cargill et al. (2022) Cargill, P. J., Bradshaw, S. J., Klimchuk, J. A., & Barnes, W. T. 2022, MNRAS, 509, 4420, doi: 10.1093/mnras/stab3163
  * Craig et al. (1978) Craig, I. J. D., McClymont, A. N., & Underwood, J. H. 1978, A&A, 70, 1
  * Culhane et al. (2007) Culhane, J. L., Harra, L. K., James, A. M., et al. 2007, Solar Physics, 243, 19, doi: 10.1007/s01007-007-0293-1
  * De Pontieu et al. (2022) De Pontieu, B., Testa, P., Martínez-Sykora, J., et al. 2022, The Astrophysical Journal, 926, 52, doi: 10.3847/1538-4357/ac4222
  * Demoulin et al. (1996) Demoulin, P., Henoux, J. C., Priest, E. R., & Mandrini, C. H. 1996, A&A, 308, 643
  * Downs et al. (2016) Downs, C., Lionello, R., Mikić, Z., Linker, J. A., & Velli, M. 2016, ApJ, 832, 180, doi: 10.3847/0004-637X/832/2/180
  * Field (1965) Field, G. B. 1965, ApJ, 142, 531, doi: 10.1086/148317
  * Froment et al. (2018) Froment, C., Auchère, F., Mikić, Z., et al. 2018, ApJ, 855, 52, doi: 10.3847/1538-4357/aaaf1d
  * Higginson et al. (2017) Higginson, A. K., Antiochos, S. K., DeVore, C. R., Wyper, P. F., & Zurbuchen, T. H. 2017, The Astrophysical Journal, 840, L10, doi: 10.3847/2041-8213/aa6d72
  * Johnston et al. (2019) Johnston, C. D., Cargill, P. J., Antolin, P., et al. 2019, Astronomy & Astrophysics, 625, A149, doi: 10.1051/0004-6361/201834742
  * Karpen et al. (2017) Karpen, J. T., DeVore, C. R., Antiochos, S. K., & Pariat, E. 2017, The Astrophysical Journal, 834, 62, doi: 10.3847/1538-4357/834/1/62
  * Klimchuk (2006) Klimchuk, J. 2006, Solar Physics, 234, 41, doi: 10.1007/s11207-006-0055-z
  * Klimchuk & DeForest (2020) Klimchuk, J. A., & DeForest, C. E. 2020, ApJ, 900, 167, doi: 10.3847/1538-4357/abab09
  * Klimchuk et al. (1992) Klimchuk, J. A., Lemen, J. R., Feldman, U., Tsuneta, S., & Uchida, Y. 1992, PASJ, 44, L181
  * Klimchuk & Luna (2019) Klimchuk, J. A., & Luna, M. 2019, The Astrophysical Journal, 884, 68, doi: 10.3847/1538-4357/ab41f4
  * Klimchuk et al. (2008) Klimchuk, J. A., Patsourakos, S., & Cargill, P. J. 2008, The Astrophysical Journal, 682, 1351, doi: 10.1086/589426
  * Kosugi et al. (2007) Kosugi, T., Matsuzaki, K., Sakao, T., et al. 2007, Solar Physics, 243, 3, doi: 10.1007/s11207-007-9014-6
  * Linker et al. (2011) Linker, J. A., Lionello, R., Mikić, Z., Titov, V. S., & Antiochos, S. K. 2011, The Astrophysical Journal, 731, 110, doi: 10.1088/0004-637X/731/2/110
  * Lionello et al. (2009) Lionello, R., Linker, J. A., & Mikić, Z. 2009, The Astrophysical Journal, 690, 902, doi: 10.1088/0004-637X/690/1/902
  * Malanushenko et al. (2022) Malanushenko, A., Cheung, M. C. M., DeForest, C. E., Klimchuk, J. A., & Rempel, M. 2022, The Astrophysical Journal, 927, 1, doi: 10.3847/1538-4357/ac3df9
  * Mason et al. (2019) Mason, E. I., Antiochos, S. K., & Viall, N. M. 2019, The Astrophysical Journal, 874, L33, doi: 10.3847/2041-8213/ab0c5d
  * Mason et al. (2021) Mason, E. I., Antiochos, S. K., & Vourlidas, A. 2021, The Astrophysical Journal Letters, 914, L8, doi: 10.3847/2041-8213/ac0259
  * Mikić et al. (2013a) Mikić, Z., Lionello, R., Mok, Y., Linker, J. A., & Winebarger, A. R. 2013a, ApJ, 773, 94, doi: 10.1088/0004-637X/773/2/94
  * Mikić et al. (2013b) —. 2013b, ApJ, 773, 94, doi: 10.1088/0004-637X/773/2/94
  * Mikić et al. (2018) Mikić, Z., Downs, C., Linker, J. A., et al. 2018, Nature Astronomy, 2, 913, doi: 10.1038/s41550-018-0562-5
  * Oran et al. (2013) Oran, R., Van Der Holst, B., Landi, E., et al. 2013, Astrophysical Journal, 778, doi: 10.1088/0004-637X/778/2/176
  * Pariat et al. (2009) Pariat, E., Antiochos, S. K., & Devore, C. R. 2009, Astrophysical Journal, 691, 61, doi: 10.1088/0004-637X/691/1/61
  * Parker (1988) Parker, E. 1988, The Astrophysical Journal, 330, 474
  * Parker (1953) Parker, E. N. 1953, ApJ, 117, 431, doi: 10.1086/145707
  * Peter & Bingert (2012) Peter, H., & Bingert, S. 2012, A&A, 548, A1, doi: 10.1051/0004-6361/201219473
  * Reep & Russell (2016) Reep, J. W., & Russell, A. J. B. 2016, The Astrophysical Journal, 818, L20, doi: 10.3847/2041-8205/818/1/l20
  * Rosner et al. (1978) Rosner, R., Tucker, W., & Vaiana, G. 1978, The Astrophysical Journal, 220, 643, doi: 10.1017/CBO9781107415324.004
  * Schlenker et al. (2021) Schlenker, M. J., Antiochos, S. K., MacNeice, P. J., & Mason, E. I. 2021, ApJ, 916, 115, doi: 10.3847/1538-4357/ac069d
  * Schrijver et al. (1999) Schrijver, C., Title, A., Berger, T., et al. 1999, Solar Physics, 187, 261.  https://link.springer.com/article/10.1023/A:1005194519642#citeas
  * Scott et al. (2021) Scott, R. B., Pontin, D. I., Antiochos, S. K., DeVore, C. R., & Wyper, P. F. 2021, ApJ, 913, 64, doi: 10.3847/1538-4357/abec4f
  * Shibata et al. (1994) Shibata, K., Nitta, N., Strong, K., et al. 1994, The Astrophysical Journal, 431, L51
  * Shimizu et al. (2019) Shimizu, T., Imada, S., Kawate, T., et al. 2019, in UV, X-Ray, and Gamma-Ray Space Instrumentation for Astronomy XXI, ed. O. H. Siegmund (SPIE), 6, doi: 10.1117/12.2528240
  * Spadaro et al. (1990) Spadaro, D., Zappala, R. A., Antiochos, S. K., Lanzafame, G., & Noci, G. 1990, The Astrophysical Journal, 362, 370, doi: 10.1086/169273
  * Titov (2007) Titov, V. S. 2007, ApJ, 660, 863, doi: 10.1086/512671
  * Titov et al. (2011) Titov, V. S., Mikić, Z., Linker, J. A., Lionello, R., & Antiochos, S. K. 2011, Astrophysical Journal, 731, doi: 10.1088/0004-637X/731/2/111
  * Török et al. (2018) Török, T., Downs, C., Linker, J. A., et al. 2018, The Astrophysical Journal, 856, 75, doi: 10.3847/1538-4357/aab36d
  * Tousey et al. (1973) Tousey, R., Bartoe, J.-D. F., Bohlin, J., et al. 1973, Solar Physics, 33, 265
  * van der Holst et al. (2014) van der Holst, B., Sokolov, I. V., Meng, X., et al. 2014, ApJ, 782, 81, doi: 10.1088/0004-637X/782/2/81
  * Vesecky et al. (1979) Vesecky, J. F., Antiochos, S. K., & Underwood, J. H. 1979, ApJ, 233, 987, doi: 10.1086/157462
  * Wyper et al. (2018) Wyper, P. F., DeVore, C. R., & Antiochos, S. K. 2018, The Astrophysical Journal, 852, 98, doi: 10.3847/1538-4357/aa9ffc
  * Young & Muglach (2014) Young, P. R., & Muglach, K. 2014, Solar Physics, 289, 3313, doi: 10.1007/s11207-014-0484-z

HYDRADeq