# Going Beyond the MHD Approximation: Physics-Based Numerical Solution of the
CGL Equations

Deepak Bhoriya Physics Department, University of Notre Dame, Notre Dame, USA
Dinshaw Balsara Physics Department, University of Notre Dame, Notre Dame, USA
ACMS Department, University of Notre Dame, Notre Dame, USA  Vladimir Florinski
Department of Space Science, University of Alabama in Huntsville, Huntsville,
AL, USA  Center for Space Plasma and Aeronomic Research, University of Alabama
in Huntsville, Huntsville, AL, USA  Harish Kumar Department of Mathematics,
Indian Institute of Technology, Delhi, India

###### Abstract

We present a new numerical model for solving the Chew-Goldberger-Low system of
equations describing a bi-Maxwellian plasma in a magnetic field. Heliospheric
and geospace environments are often observed to be in an anisotropic state
with distinctly different parallel and perpendicular pressure components. The
CGL system represents the simplest leading order correction to the common
isotropic MHD model that still allows to incorporate the latter’s most
desirable features. However, the CGL system presents several numerical
challenges: the system is not in conservation form, the source terms are
stiff, and unlike MHD it is prone to a loss of hyperbolicity if the parallel
and perpendicular pressures become too different. The usual cure is to bring
the parallel and perpendicular pressures closer to one another; but that has
usually been done in an ad hoc manner. We present a physics-informed method of
pressure relaxation based on the idea of pitch-angle scattering that keeps the
numerical system hyperbolic and naturally leads to zero anisotropy in the
limit of very large plasma beta. Numerical codes based on the CGL equations
can, therefore, be made to function robustly for any magnetic field strength,
including the limit where the magnetic field approaches zero. The capabilities
of our new algorithm are demonstrated using several stringent test problems
that provide a comparison of the CGL equations in the weakly and strongly
collisional limits. This includes a test problem that mimics interaction of a
shock with a magnetospheric environment in 2D.

Computational methods — Magnetohydrodynamics — Planetary magnetospheres —
Stellar winds

##  1 Introduction

###  1.1 The MHD system and its limitations

The magnetohydrodynamic (MHD) equations were first developed by Alfvén (1943)
for describing the large-scale, low-frequency motion of magnetized plasmas.
Inbuilt in the MHD single-fluid approximation is the assumption that the
plasma is neutral on the scales being considered and that the distributions of
ions and electrons that make up the plasma are isotropic, Maxwellian and in
equipartition. The assumption of charge neutrality eliminates the Debye
frequency which operates on very small scales. The assumption of an isotropic
Maxwellian distribution of particle velocities also presumes that the plasma
can draw on certain collisional processes to ensure the Maxwellian
distribution for the velocities of the electrons and ions on all the scales of
interest. The approximation also results in a single density and temperature
to describe the thermodynamics as well as the entropy of the particles that
make up the plasma. The MHD equations find wide-ranging applicability in
astrophysics, space physics, fusion research and many other fields of
engineering. Indeed, working within the MHD approximation has served us
extremely well! Thanks to the MHD equations, we have models for the magnetized
solar wind (Weber & Davis Jr, 1967; Steinolfson et al., 1975), the Earth’s
magnetosphere (Zwan & Wolf, 1976; Pudovkin & Semenov, 1977), the heliosphere
(Baranov & Zaitsev, 1995; Pogorelov & Semenov, 1997), the launching of winds
and jets (Blandford & Payne, 1982), dynamo theory (Parker, 2019; Ruzmaikin et
al., 1988), MHD turbulence (Sridhar & Goldreich, 1994; Goldreich & Sridhar,
1995, 1997), the multiphase interstellar medium (Mac Low et al., 2005), and
other applications in astrophysics and space physics that are too numerous to
mention. Numerical MHD also plays an extremely important role in plasma fusion
(Robinson et al., 2008; Hoelzl et al., 2021).

However, it is also well known that the assumption of an isotropic pressure
that is inbuilt into MHD is often violated in space plasmas where the Coulomb
collision rate is negligible. A compression or expansion of a parcel of plasma
associated with a flow can lead to a preferential increase or decrease in the
pressure component parallel to the deformation. The best known example is a
magnetic field pileup on the day side of the Earth’s magnetopause that leads
to a condition where the pressure in the directions perpendicular to the field
lines exceeds the pressure in the parallel direction – this is the so-called
plasma depletion layer (Crooker et al., 1979; Gary et al., 1993). Similar
depletion layers were subsequently discovered in front of the magnetospheres
of other planets (Gershman et al., 2013; Masters et al., 2014) and were
proposed for the heliopause (Cairns & Fuselier, 2017). Proton temperatures in
the solar wind were also observed to be anisotropic especially during periods
an enhanced wave activity (Marsch et al., 1981; Gary et al., 2001; Hellinger
et al., 2006).

Space plasmas follow a natural trend to maximize their entropy, i.e., relax to
an equilibrium state with an isotropic pressure. A collisionless MHD
description tacitly assumes strong interaction between particles and magnetic
fluctuations causing scattering in pitch angle so frequent, that any departure
from isotropy becomes instantly erased. However, in real space environments
the mean scattering time is often comparable with the dynamical timescales of
the system, and a complete relaxation does not occur. Clearly, standard MHD is
not sufficient to model such anisotropic environments numerically. The
simplest solution is to assume that the pressure tensor is gyrotropic (axially
symmetric about the direction of the magnetic field) and use two pressure
variables instead of one; this is the basis of the Chew-Goldberger-Low (CGL)
model discussed in the next subsection.

###  1.2 CGL as the simplest anisotropic MHD representations

Chew et al. (1956) proposed an extension to the MHD system based on
conservation of the first two adiabatic invariants of charged particle motion
in a magnetic field describing a preservation of magnetic flux in a
deformation normal to and parallel to the mean magnetic field, respectively.
This approximation is often known as the CGL approximation, and it is also
referred to as the double adiabatic approximation. If we take ${D}/{Dt}$ to be
the advective derivative, then the CGL theory posits that the pressure
perpendicular to the magnetic field, $p_{\perp}$, and the pressure parallel to
the magnetic field, ${p}_{\parallel}$, evolve according to

$\frac{D}{Dt}\left(\frac{{p_{\perp}}}{\rho
B}\right)=0,\quad\frac{D}{Dt}\left(\frac{{p_{\parallel}}{{B}^{2}}}{{{\rho}^{3}}}\right)=0.$
(1)

Here $\rho$ is the fluid density and “$B$” is the magnitude of the magnetic
field. The CGL equations govern large (fluid) scale phenomena, where the
physical processes take place on timescales that are much larger than the
proton gyro-period  and varies over distances much greater than the Larmor
radius (e.g., Burgers, 1969).

To illustrate the effects of the above two equations, as they relate to the
double adiabatic approximation, let us consider a simple thought experiment.
Let us consider a flux tube with length “L” and area “$\Sigma$”. Then, as the
flux tube is stretched or compressed, we have $B\propto{{\Sigma}^{-1}}$and
$\rho\propto{{\left(L\Sigma\right)}^{-1}}$. With these scalings, we can also
define the temperature perpendicular and parallel to the flux tube as
${{T}_{\parallel}}\equiv{{{p}_{\parallel}}}/{\rho}\;$ and
${{T}_{\perp}}\equiv{{{p}_{\perp}}}/{\rho}\;$. We then see from the above two
equations that ${{T}_{\parallel}}\propto{{L}^{-2}}$ and
${{T}_{\bot}}\propto{{\Sigma}^{-1}}$, showing that, as the flux tube is
stretched or compressed, different types of changes are brought about in the
parallel and perpendicular temperatures. The CGL formalism neglects kinetic
effects leading to collective interaction between particles and waves, and
therefore admits unrestricted values of the anisotropy in pressures, defined
as $A\equiv{{{p}_{\perp}}}/{{p_{\parallel}}}\;$.

In reality, particles are affected by electromagnetic fluctuations, either of
ambient origin (what is commonly known as turbulence) or self-generated as a
result of an instability (Parker, 1958). Several types of plasma instabilities
can arise in a plasma with an anisotropic pressure distribution (Chandrasekhar
et al., 1958; Barnes, 1966; Kennel & Petschek, 1966). These processes place
constraints on the range of anisotropy values attainable in space plasmas.
Thus, the physics of producing anisotropy (via the double adiabatic
approximation) as well as the physics of restoring isotropy (via inclusion of
the physics of the above-mentioned instabilities and turbulence) needs to be
put into our large-scale models for astrophysical and space plasmas. The goal
of this paper is to deliver on this need.

Formally, the only difference between the MHD and CGL system of equations is
the representation of the pressure tensor. Whereas for the former
$\mathbf{p}=p\mathbf{I}$, where $\mathbf{I}$ is the identity matrix, the
latter has a diagonal form in the coordinate frame aligned with the mean
magnetic field. While that form is convenient for particle-based models, it is
not useful in most grid-based models that rely on a fixed coordinate frame. In
that frame, the expression for the pressure tensor becomes

$\mathbf{p}=p_{\perp}\mathbf{I}+(p_{\parallel}-p_{\perp})\mathbf{bb},$ (2)

where $\mathbf{b}=\mathbf{B}/B$ is the unit vector in the direction of the
magnetic field. To describe the relative strength of thermal effects compared
with magnetization effects in the double adiabatic approximation, two plasma
beta parameters and their average are introduced,

$\beta_{\parallel}=\frac{8\pi
p_{\parallel}}{B^{2}},\quad\beta_{\perp}=\frac{8\pi
p_{\perp}}{B^{2}},\quad\bar{\beta}=\frac{\beta_{\parallel}+2\beta_{\perp}}{3}.$
(3)

Note that expression (2) is vaguely defined in the limit of a very weak
magnetic field, which has in the past led to numerical difficulties in CGL
codes. We will show in this paper (subsection 4.3) that a physics-based
approach overcomes this difficulty.

###  1.3 Previous work and the guidance it provides

The CGL system consists of five equations, three of which are scalar and two
are vector equations. Most authors use equations for the density, velocity,
and magnetic field in a conservative form (similar to MHD). However, there is
no accepted standard for incorporating equations (1) into the system, and
multiple approaches can be found in the literature. The paper of Hesse & Birn
(1992) represents an early implementation of the CGL system in a three-
dimensional numerical model of magnetotail dynamics. The authors used
evolution equations for $p_{\parallel}$ and $p_{\perp}$, while the total
energy was not explicitly conserved. It was also recognized that a relaxation
term was needed to keep the anisotropies from becoming unphysically large.
Hesse & Birn (1992) in particular used a relaxation term that drove the system
toward a marginally stable state on a timescale $\tau$ given by the linear
growth rates of the temperature anisotropy plasma instabilities (the firehose
and mirror instabilities, see below). Denton & Lyon (1996) developed a 2D
model CGL model for magnetosheath simulations that evolved the adiabatic
invariants (1) directly. The paper also introduced a bounded anisotropy model,
where instead of a relaxation term the code re-adjusted $p_{\parallel}$ and
$p_{\perp}$ to the nearest state of marginal stability (MS) without affecting
the trace of the pressure tensor at the end of each time step. An improved
version was presented in Denton & Lyon (2000).

A steady state 3D model of the day side magnetosheath was published by Erkaev
et al. (1999). These authors also used the instability threshold conditions as
a “fence” preventing the system from leaving the stability region. Pudovkin et
al. (1999) investigated pressure anisotropies in the magnetosheath with a
steady state 2D code. A different technique was demonstrated for anisotropy
reduction where the relaxation time was based on pitch angle scattering time.
In this approach, the pressure ratio approaches unity with time instead of
remaining in the state of marginal stability. To avoid the difficulties in
calculating the scattering rates from physical principles the model used a
constant scattering time as a free parameter. Samsonov et al. (2001) argued
for the existence of a magnetopause plasma depletion layer based on the
results from their 3D time-dependent anisotropic MHD model. That model solved
the energy conservation law plus a separate equation for the anisotropy $A$. A
threshold method of isotropization was used with $\tau$ again treated as a
free parameter. Note that the papers discussed so far were specifically
intended for magnetospheric simulations rather than for general use.

The first general framework for the CGL-MHD system was presented by Meng et
al. (2012b). The model evolved parallel and perpendicular pressure equations
in a non-conservative form. The properties of the characteristic matrix of the
system of 9 equations and the corresponding eigenvalues (wave speeds) were
discussed in detail. The relaxation method was also based on instabilities,
and the system was driven toward the nearest stable state. The relaxation time
was taken to be a constant for simplicity. The next revision of this model
(Meng et al., 2012a) used an improved relaxation time based on the growth rate
of the most unstable modes. Santos-Lima et al. (2014) studied turbulence on
very large (cosmological) scales with a model that combined eqns. (1) into a
single conservation law for the quantity $A(\rho/B)^{3}$. Note that this
quantity can become very large in the gas-dynamic limit ($\beta\to\infty$),
which would be undesirable for general-purpose models. Pressure relaxation was
applied separately from the main scheme by transforming to primitive
variables; the relaxation time was a free parameter. The authors pointed out
that non-Hall MHD models do not resolve gyroscales, and therefore, an attempt
to incorporate a kinetic wave growth scale would be very similar to a “hard
fence” approach. The updated model (Santos-Lima et al., 2016) incorporated
anisotropy relaxation terms derived from quasi-linear theory (Seough & Yoon,
2012; Yoon & Seough, 2012). It was found that the rate of anisotropy creation
via a flow deformation is, in most cases, an extremely small fraction of the
scattering rate, a fact lending support to the simple bounded anisotropy
model.

The above-mentioned papers used a Riemann solver to compute the interface
fluxes in their numerical codes. An alternative approach, presented in Luo et
al. (2023), developed a kinetic flux-splitting method of flux calculation by
taking moments of a bi-Maxwellian distribution. The model was tested on the
magnetic reconnection problem driven by numerical dissipation. Sharma et al.
(2006) extended the double adiabatic conservation laws with an addition of the
heat flux. In expressing the heat flux components through lower order moments
they use a Landau fluid closure, but characterized by a single wavelength that
was a multiple of the grid cell size. Unlike conventional CGL, that mode
correctly predicts the growth rate of the mirror instability. The authors used
more elaborate relaxation times that incorporated the instability growth rates
on the scales comparable with particle’s Larmor radius, but the resulting
relaxation method was conceptually similar to the hard fence. Hu et al. (2010)
compared the solutions to the CGL + Landau fluid heat flux system to a kinetic
particle-mesh simulation. They reported substantial differences between the
parallel heat fluxes in the two models; but otherwise their results suggested
that the fluid model in general is a viable alternative to the kinetic
approach for the Earths’s magnetosphere.

A different approach to anisotropic MHD that does not involve solving the CGL
system was proposed by Hirabayashi et al. (2016). Their model evolved the ten-
moment one fluid set of equations incorporating a general non-gyrotropic
pressure tensor. Gyrotropization was achieved by a relaxation term. While this
method has the advantage of being able to readily deal with situations where
$B\to 0$, it is more cumbersome to implement and does not preserve the
familiar characteristics of the MHD system of equations. The detailed
eigenstructure for this system is also not available making it difficult to
design very stable numerical schemes that rely on reconstruction of the
eigenweights.

The present paper seeks to develop a general-purpose formalism for the CGL
system of equations that can be applied to a broad range of anisotropic space
plasma environments with (in principle) arbitrary values of the plasma beta,
including the gas-dynamic limit. The model strictly enforces the conservation
of energy and preserves the divergence-free nature of the magnetic field,
which was not addressed in earlier developments. In addition, our model
introduces an “elastic fence” approach to anisotropy relaxation wherein the
system tends to equilibrate the pressures at a rate that is driven by a
physical relaxation time for most values of the anisotropy, as long as the
anisotropy is in the physical range (i.e. inside the fence). However, for
situations where the system comes close to the hard boundaries of the fence,
it is smoothly allowed to relax back at a faster rate so that the system never
fully reaches the “hard fence”. There is a logical reason for this compromise:
If the anisotropy reaches the hard fence, the hyperbolic system loses its
physical realizability – the eigenvalues of the system change character from
real to imaginary! It is, therefore, safer to keep a small margin which
prevents the anisotropy from reaching the hard fence. If the system approaches
the hard fence, our strategy is to modify the relaxation time to smoothly move
the system away from the hard fence. This ensures that the physical variables
in the hyperbolic system always stay within the domain of hyperbolicity. In
later sections we describe how special care is taken in the implementation of
this “elastic fence” to ensure that the mathematical system remains within
physical bounds. Interestingly, we show that this also provides a salutary
cure to the CGL equations when the magnetic field tends to zero; i.e., our
strategy smoothly retrieves an isotropic pressure for the CGL equations in the
limit when the magnetic field tends to zero.

##  2 The CGL Equations and the Numerical Challenges that they Pose

The CGL equations can be written in the cgs units as a non-conservative
hyperbolic system with a stiff source term as follows:

$\displaystyle\frac{\partial\rho}{\partial t}+\nabla\cdot(\rho\mathbf{v})$
$\displaystyle=0$ (4) $\displaystyle\frac{\partial(\rho\mathbf{v})}{\partial
t}+\nabla\cdot\left[\rho\mathbf{v}\mathbf{v}+p_{\perp}\textbf{I}+(p_{\parallel}-p_{\perp})\mathbf{b}\mathbf{b}-\frac{1}{4\pi}\left(\mathbf{B}\mathbf{B}-\frac{B^{2}}{2}\textbf{I}\right)\right]$
$\displaystyle=0$ (5)
$\displaystyle\frac{\partial(p_{\parallel}-p_{\perp})}{\partial
t}+\nabla\cdot\left[(p_{\parallel}-p_{\perp})\mathbf{v}\right]+(2p_{\parallel}+p_{\perp})\mathbf{b}\cdot\nabla\mathbf{v}\cdot\mathbf{b}-p_{\perp}\nabla\cdot\mathbf{v}$
$\displaystyle=-\frac{1}{\tau}\left(p_{\parallel}-p_{\perp}\right)$ (6)
$\displaystyle\frac{\partial\mathcal{E}}{\partial
t}+\nabla\cdot\left[\mathbf{v}\left(\mathcal{E}+p_{\perp}+\frac{B^{2}}{8\pi}\right)+\mathbf{v}\cdot\left((p_{\parallel}-p_{\perp})\mathbf{b}\mathbf{b}-\frac{\mathbf{B}\mathbf{B}}{4\pi}\right)\right]$
$\displaystyle=0$ (7) $\displaystyle\frac{\partial\mathbf{B}}{\partial
t}+\nabla\times\left[-(\mathbf{v}\times\mathbf{B})\right]$ $\displaystyle=0$
(8)

In the above set of equations, $\rho$ represents the total density and
$\mathbf{v}=(v_{x},v_{y},v_{z})^{\top}$ denotes the velocity vector. The
relaxation time is denoted by $\tau$. The system is closed under the following
equation of state

$\mathcal{E}=\dfrac{1}{2}\rho{v}^{2}+\frac{B^{2}}{8\pi}+\frac{3}{2}\bar{p}$
(9)

where $\bar{p}\equiv\dfrac{2p_{\perp}+p_{\parallel}}{3}$ is the average
pressure. Our approach of solving for the pressure difference is uncommon in
the literature, but we never experienced any difficulties with this
formulation as can be seen from the results presented in §7.

We see that eqns. (4) through (7) evolve zone-centered variables; however,
eqn. (8) for the magnetic field requires the evolution of facial magnetic
fields. Since we want a general formulation for the CGL equations that can
hold at all orders of accuracy, we will have to pay attention to the non-
oscillatory reconstruction of zone-centered variables as well as face-centered
variables. The zone-centered variables can be reconstructed via well-known
Weighted Essentially Non-Oscillatory (WENO) schemes. For a very efficient
implementation of these schemes, please see the supplement in Balsara et al.
(2023). The reconstruction of face-centered variables has been described in
Balsara (2001, 2004, 2009) and Balsara et al. (2018). Therefore, while the
implementation of such methods may be slightly challenging, we will regard it
as a solved problem.

Once the solution has been reconstructed in space, it has to be evolved in
time. To keep our approach as simple as possible, we will rely on multi-stage
Runge-Kutta methods; but the next paragraph will show that the CGL system has
intricacies which require a special form of Runge-Kutta timestepping method.
Realize though that every stage in any Runge-Kutta method will, in the very
least, require the application of a one-dimensional Riemann solver at zone
boundaries. We see that the mass, momentum and total energy equations, i.e.
eqns. (4), (5) and (7) respectively, are in flux conservation form. Such PDEs
that are in flux conservation form are familiar to the astrophysics community.
However, notice that eqn. (6) for the evolution of the difference between the
parallel and perpendicular pressures, is not in conservation form. Riemann
solvers that can handle non-conservative products in the hyperbolic PDE are
not so well-known in the astrophysics community. We will, therefore, need a
Riemann solver that is designed to accommodate non-conservative products
(Dumbser & Balsara, 2016). In such circumstances, the Riemann solver has to
return fluctuations and the PDE has to be updated in fluctuation form. The
presence of a non-conservative product implies that the Riemann solver becomes
more complicated, and more expensive because it has to be iterated to
convergence. Fortunately, much of the PDE system is indeed in conservation
form, so our method should be such as to preserve conservation whenever it is
mandated by the governing equations. Even so, our method should be flexible
enough to provide proper upwinded treatment for the non-conservative products.
That is one of the challenges that we will have to face in this paper when
formulating methods for the CGL equations. A known solution has been presented
in Appendix C of Dumbser & Balsara (2016) in the form of the HLLEMNC
pseudocode. As a result, while this may be challenging to implement, it is
viewed here as a challenge that can be met, and has a good resolution in the
literature.

Notice too that the induction equation, eqn. (8), requires that facial
magnetic fields be updated with edge-centered electric fields. The electric
field will have to be obtained from a multidimensional Riemann solver.
Multidimensional Riemann solvers will therefore be needed for obtaining the
electric field at the edges of the mesh. Fortunately, the multidimensional
Riemann problem for CGL is not any different from the multidimensional Riemann
problem for MHD (Balsara, 2010, 2012, 2014; Balsara & Nkonga, 2017);
therefore, this is regarded as a challenge that is resolved.

We also see that the right-hand side of eqn. (6) has a relaxation term. When
the relaxation timescale “$\tau$” (tau) becomes small, the parallel and
perpendicular pressures must approach one another. Note however, that “$\tau$”
will not be a pre-specified constant; instead it will depend on the physics of
the problem and will change as the local flow variables change in each zone.
Since the relaxation timescale can become solution-dependent, one cannot use
simpler exponential time differencing-based Runge-Kutta methods. Observe from
the right-hand side of eqn. (6) that when the difference in pressures is non-
zero and the relaxation time is small, a stiffly stable method for treating
the source terms is not guaranteed to preserve the sign of the pressure
difference! Allowing the source terms to change the sign of the pressure
difference can cause unphysical oscillations in the pressure. Furthermore, one
class of plasma instabilities operates when the parallel pressure is greater
than the perpendicular pressure and another class of plasma instabilities
operates when the parallel pressure is less than the perpendicular pressure.
This shows us that it is important for the source term treatment to preserve
the sign of $\left(p_{\parallel}-p_{\perp}\right)$. We, therefore, need a
special type of Runge-Kutta treatment of source terms which guarantees that
the sign of the pressure difference cannot be changed by the source terms
acting by themselves. (Note though that the left-hand side of eqn. (6) can
indeed change the sign of the pressure difference.) In other words, we need a
sign-preserving Runge-Kutta treatment of source terms. That is another
challenge associated with the CGL equations. While such sign-preserving Runge-
Kutta schemes have appeared in the recent literature (Huang & Shu, 2018), this
idea has not percolated into the astrophysical literature, and we will provide
some detail here.

It should also be observed that the flux terms and other spatial gradient
terms in eqns. (4), (5) and (6) require an evaluation of the unit vector,
“$\mathbf{b}\equiv{\mathbf{B}}/B$”, that is in the direction of the magnetic
field. This requires that a well-defined magnetic field should be present at
all points within a zone. As the magnetic field tends to zero, parts of those
flux terms could become singular. Therefore, an approach has to be found to
prevent that.

Most importantly, the CGL system is prone to a loss of hyperbolicity. The loss
of hyperbolicity arises when the pressure anisotropy
$A\equiv{{{p}_{\perp}}}/{{{p}_{\parallel}}}\;$ becomes much larger than unity
or much smaller than unity. As shown in Kato et al. (1966) , the extent or the
anisotropy that can be tolerated without loss of hyperbolicity is determined
by the strength of the local magnetic field. We examine this issue further in
§3.3.

##  3 Innovations in Numerical Methods for the CGL Equations

We split this section into three sub-sections. The first subsection shows how
the timestepping is modified to yield a sign-preserving update strategy that
can accommodate the hyperbolic terms while treating the source term in a
fashion that is sign-preserving and unconditionally stable. The second
subsection shows how this is integrated into a scheme that is written in
fluctuation form; which can, therefore, accommodate non-conservative products.
The third subsection catalogues the major results from an eigenvalue analysis
of the CGL system which also reveals the limits where the system loses
hyperbolicity.

Figure 1: The arrangement of the spatial nodes in the third or fourth order
accurate RK-WENO algorithm for two scaled space dimensions. The nodes within
four abutting spatial zones are shown by the black dots. The red, double-sided
arrows indicate the application of 1D Riemann solvers at the nodal points in
the x-direction. The blue, double-sided arrows indicate the application of 1D
Riemann solvers at the nodal points in the y-direction. The green dashed
square at the right-upper vertex of zone $(i,j)$ indicates the application of
2D Riemann solvers at the vertices of the mesh.

###  3.1 The PDE System and Sign-Preserving, Stiffly Stable, Runge-Kutta Time
stepping

To keep the description tractable, we show the steps in two spatial
dimensions. Once the two-dimensional formulation is understood, the three-
dimensional extension is obvious. In 2D, we assume that we have a hyperbolic
PDE with non-conservative products of the form

$\displaystyle{{\partial}_{t}}\mathbf{U}+{{\partial}_{x}}\mathbf{F}\left(\mathbf{U}\right)+{{\partial}_{y}}\mathbf{G}\left(\mathbf{U}\right)+\mathbf{P}\left(\mathbf{U}\right){{\partial}_{x}}\mathbf{U}+\mathbf{Q}\left(\mathbf{U}\right){{\partial}_{y}}\mathbf{U}=\dfrac{1}{\tau}\mathbf{S}\left(\mathbf{U}\right)$
(10)

Here “$\mathbf{U}$” is the solution vector,
“$\mathbf{F}\left(\mathbf{U}\right)$” and
“$\mathbf{G}\left(\mathbf{U}\right)$” are the x- and y-fluxes,
“$\mathbf{P}\left(\mathbf{U}\right)$” and
“$\mathbf{Q}\left(\mathbf{U}\right)$” are matrices with non-conservative
products and “$\mathbf{S}\left(\mathbf{U}\right)$” is a source term, assumed
to be stiff. It is easy to see that Eqns. (4) to (8) have this form where the
fluxes in the x- and y-directions are easily identified and the non-
conservative products in Eqn. (6) can be written very simply as matrices
“$\mathbf{P}\left(\mathbf{U}\right)$” and
“$\mathbf{Q}\left(\mathbf{U}\right)$”. We emphasize that the source is stiff
by identifying ${\mathbf{S}\left(\mathbf{U}\right)}/{\tau}\;$ in Eqn. (10)
with the term $-{\left({{p}_{\parallel}}-{{p}_{\bot}}\right)}/{\tau}\;$ in
Eqn. (6). Notice that in order to introduce a physics-based “elastic fence”,
we will have to eventually make “$\tau$” solution-dependent. To write the
Runge-Kutta scheme compactly, we have to compress the notation so that we
write the update in Eqn. (10) as

$\displaystyle{\partial_{t}}{\mathbf{U}}=\mathcal{L}\left({\mathbf{U}}\right)+\frac{1}{\tau}{\mathbf{S}}\left({\mathbf{U}}\right)\text{\qquad
with
\qquad}\mathcal{L}\left({\mathbf{U}}\right)\equiv-{\partial_{x}}{\mathbf{F}}\left({\mathbf{U}}\right)-{\partial_{y}}{\mathbf{G}}\left({\mathbf{U}}\right)-{\mathbf{P}}\left({\mathbf{U}}\right){\partial_{x}}{\mathbf{U}}-{\mathbf{Q}}\left({\mathbf{U}}\right){\partial_{y}}{\mathbf{U}}$
(11)

In order to have an “elastic fence” that keeps
$\left({{p_{\parallel}}-{p_{\bot}}}\right)$ within the physically allowed
bounds, we will need a method that can treat the source term in an
unconditionally stable fashion, even in the limit where the relaxation time
$\tau\to 0$. This property is called stiff stability. While several approaches
for treating the source term can be stiffly stable, many of them will not
preserve the sign of $\left({{p_{\parallel}}-{p_{\bot}}}\right)$ when the
source term, and only the source term, is applied to the PDE system. If the
source term integrator is not sign-preserving, the sign of
$\left({{p_{\parallel}}-{p_{\bot}}}\right)$ can change even when the source is
treated by itself. The unhappy consequence of that is that the source term
treatment can, by itself, cause oscillations in
$\left({{p_{\parallel}}-{p_{\bot}}}\right)$, which runs counter to the physics
of the CGL equations. We want a stiffly stable treatment of the source terms
that is sign-preserving. Technically, such schemes are known as bound-
preserving schemes. Such bound-preserving schemes were designed with second
and third-order temporal accuracy by Huang & Shu (2018). The trick is to make
the timestepping as close as possible to an exponential time differencing
scheme, which is guaranteed to be bounds preserving because of the very nature
of the exponential function. (In the next section we also describe how the
relaxation time “$\tau$” is made solution-dependent in order to construct a
“soft fence” that is physics-based.) We make a small modification of Huang &
Shu (2018) to obtain a bounds-preserving, two-stage, second-order, Strong
Stability Preserving Runge-Kutta scheme that takes the solution from a state
vector ${{\mathbf{U}}^{n}}$ at time ${t^{n}}$ to a state
${{\mathbf{U}}^{n+1}}$ at time ${t^{n+1}}={t^{n}}+\Delta t$, can be written
as:

$\displaystyle{\mathbf{U}}^{(1)}$
$\displaystyle=\dfrac{1}{1+z+\dfrac{1}{2}z^{2}}\left({\mathbf{U}}^{n}+\Delta
t\mathcal{L}(\mathbf{U}^{n})\right),$ (12a) $\displaystyle{\mathbf{U}}^{n+1}$
$\displaystyle=\dfrac{1}{2\left(1+z+\dfrac{1}{2}z^{2}\right)}\textbf{U}^{n}+\dfrac{1}{2}\left({\mathbf{U}}^{(1)}+\Delta
t\mathcal{L}(\mathbf{U}^{(1)})\right),$ (12b) where $z=\dfrac{\Delta
t}{\tau}$.

We are assuming that we have a linear source term within each zone; however
the relaxation time can vary from one zone to the next. Notice that the scheme
is explicit in the treatment of the source term, a favorable property that it
derives from the original exponential time differentiation. The corresponding
three-stage, third-order, Strong Stability Preserving Runge-Kutta scheme can
be written as:-

$\displaystyle{\mathbf{U}}^{(1)}$
$\displaystyle=\dfrac{1}{1+z+\dfrac{1}{2}z^{2}+\dfrac{1}{6}z^{3}}\left({\mathbf{U}}^{n}+\Delta
t\mathcal{L}(\mathbf{U}^{n})\right),$ (13a) $\displaystyle{\mathbf{U}}^{(2)}$
$\displaystyle=\dfrac{3}{4\left(1+\dfrac{1}{2}z+\dfrac{1}{8}z^{2}+\dfrac{1}{48}z^{3}\right)}\textbf{U}^{n}+\dfrac{1+z+\dfrac{1}{2}z^{2}+\dfrac{1}{6}z^{3}}{4\left(1+\dfrac{1}{2}z+\dfrac{1}{8}z^{2}+\dfrac{1}{48}z^{3}\right)}\left({\mathbf{U}}^{(1)}+\Delta
t\mathcal{L}(\mathbf{U}^{(1)})\right),$ (13b)
$\displaystyle{\mathbf{U}}^{n+1}$
$\displaystyle=\dfrac{1}{3\left(1+z+\dfrac{1}{2}z^{2}+\dfrac{1}{6}z^{3}\right)}\textbf{U}^{n}+\dfrac{2\left(1+\dfrac{1}{2}z+\dfrac{1}{8}z^{2}+\dfrac{1}{48}z^{3}\right)}{3\left(1+z+\dfrac{1}{2}z^{2}+\dfrac{1}{6}z^{3}\right)}\left({\mathbf{U}}^{(2)}+\Delta
t\mathcal{L}(\mathbf{U}^{(2)})\right),$ (13c) where $z\equiv\dfrac{\Delta
t}{\tau}$.

###  3.2 Detailed Description of a Single Stage in the Update

We realize that a higher-order WENO reconstruction can be constructed for the
zone-centered variables as well as the face-centered variables using the WENO
methods that are cited soon after eqn. (9). This gives us a spatially high-
order representation of the solution within each zone of the computational
domain. Now please see Figure 1, which shows four contiguous zones along with
the fourth-order accurate nodal points that are used in each zone. The nodal
points can be integrated in any one dimension using the Simpson rule and in
multiple dimensions using a tensor product of the Simpson rule. Within each
physical zone, we map $\left({x,y}\right)\in\left[{-\Delta x/2,\Delta
x/2}\right]\times\left[{-\Delta y/2,\Delta y/2}\right]$ to the reference
element
$\left({\xi,\psi}\right)\in\left[{-1/2,1/2}\right]\times\left[{-1/2,1/2}\right]$
so that all points within a zone can be accessed by ranging from $-1/2$ to
$1/2$ in each of the local coordinate directions for that zone. In eqns. (12)
and (13) we have already described the temporal update. Our task in this
section is, therefore, to describe the spatial update by describing how
$\mathcal{L}\left({\mathbf{U}}\right)$ is obtained in any zone.

Now imagine four zones like the ones in Figure 1 being placed side by side. We
realize that at each interface between any two neighboring zones we have two
sets of states, one from one side of the zone boundary and the other from the
other side of the zone boundary. These two states can be sent as input states
into a one-dimensional Riemann solver. The Riemann solver gives us information
about the wave contributions that move to the left of the zone boundary and to
the right of the zone boundary. Appendix C of Dumbser & Balsara (2016) gives
us one example of such a Riemann solver that can work with PDEs in
conservation form as well as PDEs with non-conservative products. The
fluctuation formulation then gives us an equivalent view of how to extract the
wave contributions that move to the left of the boundary and to the right of
the boundary. The Riemann solver applied in the x-direction at a zone boundary
gives us two sets of fluctuations. At each x-face of the mesh we get left-
going x-directional fluctuations at the nodal points in that face, which we
denote as ${\mathbf{D}^{-x}}$, as well as right-going fluctuations in the same
face, which we denote as ${\mathbf{D}^{+x}}$. These fluctuations in either
x-direction are shown by the red double-sided arrows in Figure 1. Likewise, at
each y-face of the mesh we get downward-going y-directional fluctuations at
the nodal points in that face, which we denote as ${\mathbf{D}^{-y}}$, as well
as upward-going y-directional fluctuations in the same face, which we denote
as ${\mathbf{D}^{+y}}$. These fluctuations in either y-direction are shown by
the blue double-sided arrows in Figure 1. We can now write the update for the
zone-centered variables as

$\displaystyle\mathcal{L}\left({\mathbf{U}}\right)=$
$\displaystyle-\frac{1}{{\Delta
x}}\left\\{{\int\limits_{-1/2}^{1/2}{\left[{{\mathbf{D}}^{-x}\left({\xi=1/2,\psi}\right)}\right]d\psi}+\int\limits_{-1/2}^{1/2}{\left[{{\mathbf{D}}^{+x}\left({\xi=-1/2,\psi}\right)}\right]d\psi}}\right\\}$
(14) $\displaystyle-\frac{1}{{\Delta
y}}\left\\{{\int\limits_{-1/2}^{1/2}{\left[{{\mathbf{D}}^{-y}\left({\xi,\psi=1/2}\right)}\right]d\xi}+\int\limits_{-1/2}^{1/2}{\left[{{\mathbf{D}}^{+y}\left({\xi,\psi=-1/2}\right)}\right]d\psi}}\right\\}$
$\displaystyle-\frac{1}{{\Delta
x}}\left\\{{\int\limits_{-1/2}^{1/2}{\left[{{\mathbf{F}}\left({\xi={{\left({1/2}\right)}^{-}},\psi}\right)}\right]d\psi}-\int\limits_{-1/2}^{1/2}{\left[{{\mathbf{F}}\left({\xi=-{{\left({1/2}\right)}^{+}},\psi}\right)}\right]d\psi}}\right\\}$
$\displaystyle-\frac{1}{{\Delta
y}}\left\\{{\int\limits_{-1/2}^{1/2}{\left[{{\mathbf{G}}\left({\xi,\psi={{\left({1/2}\right)}^{-}}}\right)}\right]d\xi}-\int\limits_{-1/2}^{1/2}{\left[{{\mathbf{G}}\left({\xi,\psi=-{{\left({1/2}\right)}^{+}}}\right)}\right]d\psi}}\right\\}$
$\displaystyle-\frac{1}{{\Delta
x}}\int\limits_{-1/2}^{1/2}{\int\limits_{-1/2}^{1/2}{{\mathbf{P}}\left({\xi,\psi}\right){\partial_{\xi}}{\mathbf{U}}\left({\xi,\psi}\right)d\xi}d\psi}-\frac{1}{{\Delta
y}}\int\limits_{-1/2}^{1/2}{\int\limits_{-1/2}^{1/2}{{\mathbf{Q}}\left({\xi,\psi}\right){\partial_{\psi}}{\mathbf{U}}\left({\xi,\psi}\right)d\xi}d\psi}$

Also notice from the right upper vertex of zone $\left({i,j}\right)$ in Figure
1 that we have four states surrounding that vertex. As a result, those four
states can be sent into a multidimensional Riemann solver Balsara (2010, 2012,
2014); Balsara & Nkonga (2017) and they produce a multidimensionally upwinded
representation of the electric field at the vertex. An update for the facial
magnetic fields can, therefore, be similarly designed and this is well-
described in the above-mentioned references as well as Balsara & Spicer
(1999). Having described all the details for a single stage in the Runge-Kutta
update in this subsection, as well as having described how the different
stages have to be coupled together in the previous subsection, we can now
conclude that the full spatially and temporally high-order scheme has been
described here. In this Sub-section we have described the method at a
conceptual level. In Sub-section 5.3 we provide a step-by-step implementation
plan.

###  3.3 Eigenvalue analysis for the CGL equations

The set of eigenvalues plays a crucial role in the design of efficient
numerical schemes. For that reason, we provide the set of eigenvalues (denoted
by $\Lambda_{d}$, $d=x,y,z$) for the $d-$directional characteristic matrix
below.

$\Lambda_{d}=\\{v_{d}-m_{f},v_{d}-m_{s},v_{d}-c_{a},v_{d},v_{d},0,v_{d}+c_{a},v_{d}+m_{s},v_{d}+m_{f}\\},$
(15)

where

$\displaystyle c_{a}$
$\displaystyle=\sqrt{\frac{{B_{d}}^{2}}{4\,\rho\,\pi}-\frac{{b_{d}}^{2}\,{\left(p_{\parallel}-p_{\perp}\right)}}{\rho}},\
(\text{where}\ b_{d}=B_{d}/|\mathbf{B}|)$ $\displaystyle m_{s}$
$\displaystyle=\sqrt{\frac{-\mathrm{b}-\sqrt{{\mathrm{b}}^{2}-4\,\mathrm{a}\,\mathrm{c}}}{2\,\rho}},$
$\displaystyle m_{f}$
$\displaystyle=\sqrt{\frac{-\mathrm{b}+\sqrt{{\mathrm{b}}^{2}-4\,\mathrm{a}\,\mathrm{c}}}{2\,\rho}}.$

The expressions for the unknowns ‘a’, ‘b’ and ‘c’ in the above equations are
given by

$\displaystyle\mathrm{a}$ $\displaystyle=2\rho$ $\displaystyle\mathrm{b}$
$\displaystyle=-2\,p_{\perp}-{b_{d}}^{2}\,{\left(2\,p_{\parallel}-p_{\perp}\right)}-\frac{B^{2}}{4\,\pi}$
$\displaystyle\mathrm{c}$
$\displaystyle=-\frac{3\,{b_{d}}^{4}\,{p_{\parallel}}^{2}-{b_{d}}^{2}\,{p_{\perp}}^{2}\,{\left({b_{d}}^{2}-1\right)}-\dfrac{3\,{B_{d}}^{2}\,p_{\parallel}}{4\,\pi}+3\,{b_{d}}^{2}\,p_{\parallel}\,p_{\perp}\,{\left({b_{d}}^{2}-2\right)}}{2\,\rho}.$

Note that the CGL wave speeds do not fully reduce to the MHD wave speeds when
the parallel and perpendicular pressures are equal because of different
adiabatic indices for the plasma in those two directions. The above set of
eigenvalues ($\Lambda_{d}$) is valid within a specific admissible domain. To
understand this physical domain, we first need to introduce two limits, which
we call the lower limit (represented as $p_{\mathrm{lo}}$ ) and the upper
limit (represented as $p_{\mathrm{hi}}$ ). These limits are defined as
follows:

$p_{\mathrm{lo}}=\dfrac{p_{\perp}^{2}}{6p_{\perp}+\frac{3B^{2}}{4\pi}}\ \ ;\ \
p_{\mathrm{hi}}=\dfrac{B^{2}}{4\pi}+p_{\perp}.$ (16)

The admissible domain, $\Omega$, to maintain hyperbolicity in the CGL system
is given by

$\Omega=\\{\mathbf{W}\in\mathbb{R}^{9}:~{}\rho>0,~{}p_{\parallel},p_{\perp}>0,~{}p_{\mathrm{lo}}\leq
p_{\parallel}\leq p_{\mathrm{hi}}\\},$ (17)

where $\mathbf{W}=(\rho,\mathbf{v},p_{\parallel},p_{\perp},\mathbf{B})^{\top}$
is the vector of primitive variables. Unlike ordinary MHD, in the CGL system
the slow wave is not necessarily slower than the Alfvén wave. Therefore, the
admissible region for $p_{\parallel}$ can be divided into three parts (see
Kato et al., 1966) based on the values of Alfvén wave ($c_{A}$), slow wave
($c_{s}$) and fast wave ($c_{f}$).

  1. 1.

$p_{\mathrm{lo}}\leq p_{\parallel}\leq\dfrac{p_{\mathrm{hi}}}{4}$, $c_{s}\leq
c_{A}\leq c_{f}$

  2. 2.

$\dfrac{p_{\mathrm{hi}}}{4}\leq
p_{\parallel}\leq\dfrac{p_{\mathrm{hi}}}{4}+\dfrac{3p_{\mathrm{lo}}}{4}$,
$c_{s}\leq c_{A}<c_{f}$

  3. 3.

$\dfrac{p_{\mathrm{hi}}}{4}+\dfrac{3p_{\mathrm{lo}}}{4}\leq p_{\parallel}\leq
p_{\mathrm{hi}}$, $c_{A}\leq c_{s}<c_{f}$.

The complexity of the wave speeds and the realization that violating the
bounds in the previous three equations could cause a code crash, leads us to
an important realization: we will have to use physics to always keep guiding
the numerics in such a way that the numerics can stay within the physically
realizable domain. Eqn. (6) gives us an important clue that the relaxation
time always drives the system towards isotropic pressure distributions. But
for a given constant physical relaxation time, the system may not always stay
out of the danger zone. In this paper, the above-mentioned goal is achieved by
allowing the relaxation time to be modified when the pressure anisotropy
becomes so large as to bring the system dangerously close to the limits of
physical realizability. This is what we mean by an ”elastic fence”. The next
section provides more detail.

##  4 The physics of pressure relaxation

This section is split into three subsections. The first subsection makes a
physics-based case for having an “elastic fence”. The second subsection
provides the numerical details for implementing these ideas in code. The third
subsection shows that this concept of an “elastic fence” also makes it
possible for the CGL equations to attain an isotropic pressure in the non-
magnetized limit.

###  4.1 The physics-based case for an “elastic fence”

It was already recognized by Parker (1958) that a plasma with an anisotropic
pressure distribution could become unstable and generate waves that would lead
to a scattering of the particles and consequently drive the plasma toward
isotropy. The CGL-MHD system predicts two kinds of instability; one operating
in the regime $A<1$ (firehose) and the other requiring $A>1$ (mirror). The
firehose instability appears on the Alfven wave branch. The region of
stability is given by (Kato et al., 1966)

$\beta_{\perp}>\beta_{\parallel}-2.$ (18)

The mirror instability, on the other hand, appears on the slow magnetosonic
branch, for which the CGL stability criterion is (Kato et al., 1966)

$\beta_{\parallel}>\frac{\beta_{\perp}^{2}}{6(1+\beta_{\perp})}.$ (19)

Linear collisionless kinetic theory identifies three instabilities in a plasma
with a bi-Maxwellian proton distribution. The firehose instability affects the
Alfven-ion cyclotron (AIC) branch of the linear dispersion relation. It is
strongest for parallel propagating low-frequency ($\omega\ll\Omega$) modes
(Yoon, 1990; Seough & Yoon, 2012). The ion cyclotron instability is also
associated with the AIC mode, but has higher frequencies ($\omega\sim\Omega$);
its growth rate is maximized for parallel propagation. Finally, the mirror
instability (Yoon & Seough, 2012) is a non-propagating compressive mode with
maximum growth at oblique propagation. The mirror and AIC modes operate in a
similar range of the $A$–$\beta$ space, for which reason they are sometimes
grouped together in the literature despite being different in origin.

The kinetic instability thresholds are usually simplified to the limit of cold
plasma, in which case the firehose criterion is identical to (18). The kinetic
mirror instability has no factor of 6 in the denominator. If we delete that
factor from (19) and solve for $\beta_{\perp}$, the expression becomes

$\beta_{\perp}<\frac{\beta_{\parallel}}{2}+\sqrt{\beta_{\parallel}\left(1+\frac{\beta_{\parallel}}{4}\right)}$
(20)

Pressure anisotropies measured in the solar wind typically lie well inside the
stability range and far from the marginal stability limits (Hellinger et al.,
2006).

Expressions (18) and (19) correspond to the hyperbolicity conditions for the
CGL system presented in Section 2. As such, they cannot be violated, and a
numerical model must perform an immediate correction to keep the system in the
state of marginal stability. Note that there is no point in applying a
relaxation term if the physical realizability criteria for the CGL system are
not met because the model breaks down when the eigenvalues of the governing
system of equations become complex. Physically, the imaginary part of the
complex frequency (the growth rate) is comparable to the cyclotron frequency
$\Omega$. The quasi-linear scattering rates are typically smaller, but still
well outside of the time scales accessible to MHD models, so from that
perspective the relaxation should be effectively instantaneous.

If the micro-instabilities were the only physical processes leading to a
reduction in the anisotropy, then one would expect the plasma to often be in
the state of marginal stability. However, the measured anisotropies are
typically smaller than suggested by the stability criteria (Marsch et al.,
2004; Hellinger et al., 2006). In a large astrophysical system, such as the
solar wind, the instabilities contribute to the turbulent spectrum of magnetic
fluctuations at frequencies resonant with the particle’s gyromotion. However,
wavepower at relevant frequencies is also provided by the turbulent cascade
from larger spatial scales. If the fluctuating power is not too large, the
scattering frequencies could be much smaller than the cyclotron frequency and,
therefore, comparable to the MHD timescales. The scattering rate can be
estimated from the quasi-linear theory of interactions between particles and
waves (Jokipii, 1966; Hall & Sturrock, 1967; Achatz et al., 1991). The pitch-
angle scattering coefficient is a complicated expression that depends on the
pitch angle cosine $\mu$ and includes contributions from waves with multiple
polarizations and directions of propagation. Very crudely, the expression for
the turbulent scattering rate $\nu_{\mathrm{scat}}$ can be written as

$\nu_{\mathrm{scat}}\simeq\frac{\Omega^{2}P(k_{\mathrm{res}})}{v^{*}B^{2}},$
(21)

where $v^{*}$ is the characteristic speed (the larger of the wave phase
velocity, typically the Alfvén speed, and the particle thermal velocity), and
the $P$ is the turbulent power spectral density at the resonant wavenumber
$k_{\mathrm{res}}=\Omega/v^{*}$.

A typical core plasma particle is resonant in the dissipation turbulent range,
where the spectral power is very small. Assuming a power law spectrum
$P(k)\sim\langle\delta B^{2}\rangle l_{c}(kl_{c})^{-\gamma}$, where
$\langle\delta B^{2}\rangle$ is the mean square turbulent magnetic field
magnitude, $l_{c}$ is the turbulent outer scale, and $v^{*}\sim v_{t}$, the
thermal speed, one obtains

$\nu_{\mathrm{scat}}\simeq\Omega\frac{\langle\delta
B^{2}\rangle}{B^{2}}\left(\frac{r_{L}}{l_{c}}\right)^{\gamma-1},$ (22)

where $\langle\delta B^{2}\rangle/B^{2}$ is typically of the order of 0.1–1
(depending on the environment), and $r_{L}$ is the thermal Larmor radius. The
spectral power index is typically between $-3/2$ and $-2$ in the inertial
range and $-3$ in the dissipation range (Bruno & Carbone, 2016). It is
expected that as the mean field $B\to 0$, the factor $\langle\delta
B^{2}\rangle/B^{2}$ in expression (22) becomes very large leading to a very
efficient isotropization.

Figure 2: Domains of realizability of the CGL equations (shaded) for two
values of $\epsilon$. The domain is most compactly expressed in terms of the
plasma betas. The plasma betas for the parallel and perpendicular pressures
form the horizontal and vertical axes of the plot. The magenta line shows the
situation when the plasma pressures are isotropic. The mirror instability,
shown by the thick blue line, forms the upper boundary of the domain of
realizability for the CGL equations. The firehose instability, shown by the
thick red line (as well as part of the x-axis), forms the lower boundary of
the domain of realizability for the CGL equations. The dotted lines show
slightly inner margins within which the plasma pressures should lie in a
computation. The dashed red lines are lines of constant total pressure, and a
small patch of plasma can only move along those lines if it is to conserve
energy. (Of course, the energy of a patch of plasma is allowed to evolve with
a simulation.) The arrows in Figure 2a indicate directions along surfaces of
constant total pressure.

###  4.2 Numerical design of an “elastic fence”

Equation (22) above gives us a physics-based measure of the relaxation time so
that we can define a physical relaxation time as $\tau_{phys}\equiv
1/\nu_{scat}$. Usually, for any physical system one has an estimate for the
turbulent properties of the medium so that a physically-motivated relaxation
time can be found. Even when the pressures are very close to isotropic, i.e.
when ${p_{\parallel}}\simeq{p_{\bot}}$, this relaxation time will be
operational and will keep driving the pressures towards isotropization.
However, in a computer code for modeling collisionless plasmas, it cannot be
guaranteed that this rate of isotropization will always be sufficient to keep
the system hyperbolic. Therefore, we develop a complimentary phenomenology
based on the instability threshold expressions (18) and (20) for the firehose
and mirror instabilities, respectively. Figure 2a illustrates the situation in
the ${\beta_{\parallel}}$ versus ${\beta_{\bot}}$ plane; and expressing the
plasma pressures in terms of plasma beta is indeed the most compact way to
show the results. The line of isotropy runs at a $45^{\text{o}}$ angle in
Figure 2a and is shown by the solid magenta line in Figure 2a. The limit for
the development of a firehose instability runs parallel to the magenta line
and is shown by the red line in Figure 2a. The development of the mirror
instability is parallel to the $45^{\text{o}}$ line for large values of
${\beta_{\parallel}}$, with the result that asymptotically, i.e., for large
values of ${\beta_{\parallel}}$, eqn. (20) becomes
${\beta_{\bot}}\leq{\beta_{\parallel}}+1$. Of course, for smaller values of
${\beta_{\parallel}}$, which can indeed occur when the magnetic field is
strong, we will pay more careful attention to the permissible boundary for the
mirror instability. The limit for the mirror instability is shown by the blue
curve in Figure 2a.

Now realize that the energy equation, i.e. eqns. (7) and (9), will always
permit us to use the code variables to extract the total pressure $\bar{p}$
and the mean plasma beta $\bar{\beta}$ which can be written as:-

$\displaystyle\bar{p}\equiv{{\left({{p_{\parallel}}+2{p_{\bot}}}\right)}/3}\Rightarrow\bar{\beta}\equiv{{\bar{p}}/{\left({{{\bf{B}}^{2}}/(8\pi)}\right)}}.$
(23)

To find independent values of ${p_{\parallel}}$ and ${p_{\bot}}$ we have to
resort to eqn. (6) which gives us the evolution of
${p_{\parallel}}-{p_{\bot}}$. This gives us two equations from which to
disentangle ${p_{\parallel}}$ and ${p_{\bot}}$. In practice, the two pressures
should lie slightly inside the domain bounded by the mirror and firehose
instabilities, as shown by the dotted curves in Figure 2a. If they do not lie
within the physically realizable domain at the end of a timestep in the code,
energy should still be conserved with the result that we propose that the
parallel and perpendicular pressures should move along the lines of constant
total pressure till they do lie within the realizable domain. In other words
the pressures move along the line
$\bar{p}\equiv{{\left({{p_{\parallel}}+2{p_{\bot}}}\right)}/3}$, which is
equivalent to the dashed black lines
${\beta_{\parallel}}+2{\beta_{\bot}}=3\bar{\beta}$ in Figure 2a, until they
lie within the realizable domain. Implicit in this realization is the fact
that once the mirror or firehose instabilities are activated, they operate on
plasma timescales that are shorter than the MHD timescales that govern the
numerical timestep in the code. These fast-acting plasma instabilities re-
establish a modicum of pressure isotropization, which will at least push the
pressures within the physically realizable domain. Given that it is easier to
work with plasma $\beta^{\prime}s$, for any given value of $\bar{\beta}$ from
the energy equation, we can find the limiting values of
$\beta_{\parallel}^{m}$ and $\beta_{\parallel}^{f}$ for the mirror and
firehose instabilities. They are given by

$\displaystyle\beta_{\parallel}^{m}=2\left({\bar{\beta}+1/3}\right)-\sqrt{{{(\bar{\beta})}^{2}}+8\bar{\beta}/3+4/9}$
(24)

and

$\displaystyle\beta_{\parallel}^{f}=\left\\{{\begin{array}[]{*{20}{c}}{\bar{\beta}+4/3{\qquad\rm{for}}\qquad\bar{\beta}\geq
2/3}\\\
{3\bar{\beta}{\qquad\qquad\rm{for}}\qquad\bar{\beta}<2/3}\end{array}}\right.$
(27)

The physically permissible values of ${\beta_{\parallel}}$ are given by
$\beta_{\parallel}^{m}\leq{\beta_{\parallel}}\leq\beta_{\parallel}^{f}$;
please see Figure 2a. This explains how the parallel and perpendicular
pressures can always be brought within the physically realizable limits. The
limits shown in Figure 2a will also ensure that the system retains
hyperbolicity as described in §3.3.

Since we understand how ${p_{\parallel}}$ and ${p_{\bot}}$ can be kept within
the bounds that ensure that the hyperbolic system is physically realizable, we
now focus on obtaining a numerical relaxation time, ${\tau_{num}}$ which is
suitable for numerical work. When it is numerically safe, ${\tau_{num}}$
should approach ${\tau_{phys}}$ quite closely. However, ${\tau_{num}}$ should
be engineered so that it nudges the system away from the boundaries that are
determined by the mirror and firehose instabilities in Figure 2a. In other
words, we seek a function $f\left({{\beta_{\parallel}},\bar{\beta}}\right)$
such that

$\displaystyle{\tau_{num}}={\tau_{phys}}{\rm{}}f\left({{\beta_{\parallel}},\bar{\beta}}\right)$
(28)

When ${\beta_{\parallel}}\to\beta_{\parallel}^{m}$ or when
${\beta_{\parallel}}\to\beta_{\parallel}^{f}$ we want the relaxation time to
become very small; i.e., we want
$f\left({{\beta_{\parallel}},\bar{\beta}}\right)\to 0$. In all other
situations we want $f\left({{\beta_{\parallel}},\bar{\beta}}\right)\to 1$.
Numerous choices for such a function are possible but the function should be
such that when ${\beta_{\parallel}}<\bar{\beta}$, the plasma pressures may not
hit the boundary formed by the mirror instability. Likewise, when
${\beta_{\parallel}}>\bar{\beta}$, the plasma pressures may not hit the
boundary formed by the firehose instability. When
${\beta_{\parallel}}\to\bar{\beta}$ we want
$f\left({{\beta_{\parallel}},\bar{\beta}}\right)\to 1$ so that the physics-
based relaxation is achieved. To have maximal continuity, we want
$f\left({{\beta_{\parallel}},\bar{\beta}}\right)$ to be a continuous function
with zero slope at ${\beta_{\parallel}}=\bar{\beta}$. We see that it is
acceptable to define $f\left({{\beta_{\parallel}},\bar{\beta}}\right)$ in a
piecewise fashion as:-

$\displaystyle
f\left({{\beta_{\parallel}},\bar{\beta}}\right)=\left\\{{\begin{array}[]{*{20}{c}}{{{\cos}^{\kappa}}\left[{\dfrac{\pi}{2}{{\left|{\dfrac{{\bar{\beta}-{\beta_{\parallel}}}}{{\bar{\beta}-\beta_{\parallel}^{m}}}}\right|}^{\eta}}}\right]{\qquad\rm{when}\quad}{\beta_{\parallel}}<\bar{\beta}}\\\
{{{\cos}^{\kappa}}\left[{\dfrac{\pi}{2}{{\left|{\dfrac{{{\beta_{\parallel}}-\bar{\beta}}}{{\beta_{\parallel}^{f}-\bar{\beta}}}}\right|}^{\eta}}}\right]{\qquad\rm{when}\quad}{\beta_{\parallel}}>\bar{\beta}}\\\
{1\qquad\qquad\qquad\qquad{\rm{otherwise}}}\end{array}}\right.$ (32)

We can take $\kappa=2$ or $\kappa=4$, and also $\eta=2$ or $\eta=4$, in the
above formula to make it increasingly flat around
${\beta_{\parallel}}=\bar{\beta}$. (Even and positive values of “$\kappa$”
enable us to ensure that the derivatives of the above functions are zero at
the limits over which they operate. Increasing “$\kappa$” increases the number
of derivatives that are zero around ${\beta_{\parallel}}=\bar{\beta}$.
Increasing “$\eta$”, while keeping it an even integer helps us to increase the
range of values around $\bar{\beta}$ over which the function
$f\left({{\beta_{\parallel}},\bar{\beta}}\right)$ remains flat.) Figure 3
shows a colorized plot of the reciprocal of the above function in the
$\left({{\beta_{\parallel}},{\beta_{\bot}}}\right)$ plane. We see that as the
plasma approaches the mirror or firehose instabilities, the reciprocal of the
relaxation time becomes very large, i.e. the traversal away from the
unrealizable part of the domain becomes very rapid. Otherwise, the numerical
relaxation time should track the physical relaxation time, as evinced by unit
values around the line of pressure isotropization in Figure 3.

Figure 3: Panel shows the plot of $\log_{10}$ (min(1/$f$, 30)). We can see
that the function is smooth even at the $45^{\text{o}}$ line of isotropy. It
only assumes large values as we approach the boundaries of the mirror and
firehose instabilities. This serves to keep the simulation within the
realizable region. The choice of the number “30” in the min function above is
purely arbitrary and is intended to reveal the range of variation of “1/$f$”.

###  4.3 The “elastic fence” concept allows the CGL Equations to reach the
non-magnetic limit

One of the persistent concerns in using the CGL equations for numerical work
is that the equations seem to become badly-formed in the limit where the
magnetic field goes to zero. This was, indeed, one of the the motivating
reasons behind the paper by Hirabayashi et al. (2016). We now show that the
CGL equations, coupled with our new concept of an “elastic fence” can smoothly
produce an isotropic pressure in the limit where the magnetic field tends to
(and reaches) zero in a numerical code.

In the limit when the magnetic field becomes very small, the plasma-$\beta$
becomes very large. Therefore, we can obtain the asymptotic limits from eqns.
(24) and (27):

$\displaystyle
p_{\parallel}^{m}=\bar{p}-\frac{2}{3}\left({\frac{{{{\bf{B}}^{2}}}}{{8\pi}}}\right)\qquad;\qquad
p_{\parallel}^{f}=\bar{p}+\frac{4}{3}\left({\frac{{{{\bf{B}}^{2}}}}{{8\pi}}}\right)$
(33)

We see, therefore, that $p_{\parallel}^{m}\leq{p_{\parallel}}\leq
p_{\parallel}^{f}$ . Furthermore, we see that the interval
$p_{\parallel}^{f}-p_{\parallel}^{m}=2\left({\dfrac{{{{\bf{B}}^{2}}}}{{8\pi}}}\right)$
shrinks with decreasing magnetic field. Therefore, as the magnetic field goes
to zero, we are guaranteed to reach the limit
${p_{\parallel}}={p_{\bot}}=\bar{p}$. This ensures that we will asymptotically
reach the limit of isotropic pressure as the magnetic field tends to zero.
Notice that in that limit
${\bf{b}}={{\bf{B}}\mathord{\left/{\vphantom{{\bf{B}}{\left|{\bf{B}}\right|}}}\right.\kern-1.2pt}{\left|{\bf{B}}\right|}}$
can be any bounded unit vector that one wishes. The structure of the left-hand
side of eqn. (6) is such that any anisotropy that it induces during any
fractional timestep is washed out at the end of that fractional timestep by
our “soft fence” formulation. Furthermore, in the limit of zero magnetic
field, eqn. (5) will not depend on pressure anisotropy. Therefore, the
hydrodynamical limit will always be reached. In other words, with all the
advancements described here, the CGL equations will operate stably in the
hydrodynamic limit because they reduce exactly to the MHD equations.

##  5 Step by Step Implementation Plan

###  5.1 Nudging $\beta_{\parallel}$ and $\beta_{\perp}$ into the physical
domain

The structure of the left-hand side of eqn. (6) is such that it can always
nudge the parallel and perpendicular pressures away from one another.
Sometimes, if the terms on the left-hand side of eqn. (6) are too large, it is
always possible that the pressure anisotropy in a numerical code can become so
large as to take the system past the bounds of physical realizability. In this
section, we provide a method to nudge $\beta_{\parallel}$ and $\beta_{\perp}$
into the physical domain. We explicitly define the boundaries of the physical
domain by introducing the curves $C_{1}$, $C_{2}$ and $C_{3}$ (which are
pointed to by the black arrows in Figure 2b). To ensure a computationally
stable implementation, we maintain a positive distance ($\epsilon>0$) from the
boundaries of the actual physical domain. This practice helps to prevent
numerical errors that would arise in a numerical code if the solution were to
approach the boundaries of Figure 2b too closely. These slightly shifted
boundaries are shown with the dotted black lines in Figure 2b. To make the
dotted lines more conspicuous in Figure 2b we have used $\epsilon=0.1$. In
practice, we take $\epsilon=10^{-2}$ in all the numerical test problems. The
curve $C_{1}$ denotes the lower boundary provided by the firehose instability
with an added $\epsilon$-dependent buffer. The curve $C_{2}$ denotes the right
boundary provided by the firehose instability with an added
$\epsilon$-dependent buffer. The curve $C_{3}$ denotes the left boundary
provided by the mirror instability with an added $\epsilon$-dependent buffer.
With these notations, the required curves can be described as follows:

$\displaystyle C_{1}\qquad$
$\displaystyle:\qquad\epsilon\beta_{\parallel}=2\beta_{\perp},$ $\displaystyle
C_{2}\qquad$
$\displaystyle:\qquad\beta_{\parallel}=(2+\beta_{\perp})-\epsilon,$
$\displaystyle C_{3}\qquad$
$\displaystyle:\qquad\beta_{\parallel}=\begin{cases}\dfrac{\beta_{\perp}^{2}}{(1+\beta_{\perp})(1-\xi)}&\text{if
}\beta_{\parallel}\leq 2-2\sqrt{\frac{2}{3}}\\\
3\bar{\beta}-2\left(\beta_{\perp}^{m}-\epsilon\right)&\text{if
}\beta_{\parallel}>2-2\sqrt{\frac{2}{3}}\end{cases},$

where
$\xi=\dfrac{((-6+4\sqrt{6}-3\epsilon)\epsilon)}{(2+2(-3+\sqrt{6})\epsilon)}$,
$\beta_{\perp}^{m}=\dfrac{3\bar{\beta}-\beta_{\parallel}^{m}}{2}$ and
$\beta_{\parallel}^{m}$ is given by eqn. (24).

We call the region enclosed within the dotted curves in Figure 2b (i.e. the
region surrounded by the curves $C_{1}$, $C_{2}$ and $C_{3}$) as the
physically admissible domain, and we denote this domain by $\Omega_{\beta}$.
We desire plasma beta’s (i.e. $\beta_{\parallel}$ and $\beta_{\perp}$) to
remain within the physical region at each grid point. To achieve this
objective, we provide a step-by-step method to nudge $\beta_{\parallel}$ and
$\beta_{\perp}$ into the physical domain.

  1. 1.

Realize that the energy equation, i.e. eqns. (7) and (9), will always permit
us to extract the total pressure $\bar{p}$ and the mean plasma beta
$\bar{\beta}$ which can be written as:-

$\displaystyle\bar{p}\equiv{{\left({{p_{\parallel}}+2{p_{\bot}}}\right)}/3}\Rightarrow\bar{\beta}\equiv{{\bar{p}}/{\left({{{\bf{B}}^{2}}/(8\pi)}\right)}}.$
(34)

  2. 2.

Eqn. (6) gives us the evolution of $\Delta p\equiv{p_{\parallel}}-{p_{\bot}}$.
This, along with eqn. (34), gives us two equations which are solved to obtain
the pressure variables ${p_{\parallel}}$ and ${p_{\bot}}$.

  3. 3.

Utilizing ${p_{\parallel}}$ and ${p_{\bot}}$, we extract corresponding beta’s
(i.e. ${\beta_{\parallel}}$ and ${\beta_{\bot}}$) using the following
expressions:-

$\displaystyle\beta_{\parallel}\equiv{{p_{\parallel}}/{\left({{{\bf{B}}^{2}}/(8\pi)}\right)}}\qquad;\qquad\beta_{\perp}\equiv{{p_{\perp}}/{\left({{{\bf{B}}^{2}}/(8\pi)}\right)}}$
(35)

  4. 4.

Realize that the obtained plasma beta’s, i.e. $\beta_{\parallel}$ and
$\beta_{\perp}$ in eqns. (35), may have been nudged out of the physical domain
$\Omega_{\beta}$. We need to nudge $\beta_{\parallel}$ and $\beta_{\perp}$
back into the physical domain in a thermal energy-conserving fashion to
retrieve the physically admissible plasma beta’s $\hat{\beta}_{\parallel}$ and
$\hat{\beta}_{\perp}$.

  5. 5.

If $\beta_{\parallel}$ and $\beta_{\perp}$ have stepped out of the physical
domain, the following expressions for $\hat{\beta}_{\parallel}$ and
$\hat{\beta}_{\perp}$ will nudge them into the admissible domain
$\Omega_{\beta}$.

$\displaystyle\hat{\beta}_{\parallel}=$
$\displaystyle\max\left\\{\min\left\\{\alpha_{1},\beta_{\parallel}\right\\},3\bar{\beta}-2\alpha_{2}\right\\},\qquad;\qquad\hat{\beta}_{\perp}=\dfrac{3\bar{\beta}-\hat{\beta}_{\parallel}}{2}$
(36) where $\displaystyle\alpha_{1}$
$\displaystyle=\min\left\\{\dfrac{3\bar{\beta}}{\epsilon+1},\bar{\beta}+\dfrac{2}{3}(2-\epsilon)\right\\}\quad;\quad\alpha_{2}=\max\left\\{\dfrac{\left(-\text{b}+\sqrt{\text{b}^{2}-\text{4ac}}\right)}{\text{2
a}},\beta_{\perp}^{m}-\epsilon\right\\}\ ;$ a $\displaystyle=(3-2\xi)\ ;\
\text{b}=(1-\xi)(2-3\bar{\beta})\ ;\ \text{c}=-3(1-\xi)\bar{\beta}\quad;$
$\displaystyle\xi$
$\displaystyle=\dfrac{(-6+4\sqrt{6}-3\epsilon)\epsilon}{2+2(-3+\sqrt{6})\epsilon}\
;\ \beta_{\perp}^{m}=\dfrac{3\bar{\beta}-\beta_{\parallel}^{m}}{2}\ ;\
\beta_{\parallel}^{m}=2\left({\bar{\beta}+1/3}\right)-\sqrt{{{(\bar{\beta})}^{2}}+8\bar{\beta}/3+4/9}.$

Realize that if $\beta_{\parallel}$ and $\beta_{\perp}$ are already inside the
physical domain $\Omega_{\beta}$ then the expressions in eqn. (36) are
designed to identically retrieve the beta’s given in (35); this ensures that
the expressions in eqn. (36) do not change the beta’s in eqn. (35) if the
latter are already within the physical domain.

  6. 6.

We replace the expressions $\beta_{\parallel}$ and $\beta_{\perp}$ given in
(35) with the physically admissible beta’s $\beta_{\parallel}$ and
$\beta_{\perp}$ given in eqn. (36) (i.e., we reset
$\beta_{\parallel}=\hat{\beta}_{\parallel}$ and
$\beta_{\perp}=\hat{\beta}_{\perp}$). This in turn yields physically
realizable parallel and perpendicular pressures that can be reset in the code.
Notice that our formulation preserves the mean thermal energy in the plasma.
In other words, the $\hat{\beta}_{\parallel}$ and $\hat{\beta}_{\perp}$ lie on
the same red dashed curves in Figure 2b as the original $\beta_{\parallel}$
and $\beta_{\perp}$.

  7. 7.

When the magnetic field is below some threshold that the applications
scientist could view as zero, then we can use the logic from Sub-section 4.3
to justify the setting. The same consideration applies to our evaluation of
eqn. (32).

The above seven-step procedure can always be employed whenever one wants to go
from conserved to primitive variables in a numerical code. It ensures that the
numerical solution always lies within the domain of physical realizability.

###  5.2 Relaxation times: Physical and Numerical consideration

While the previous Sub-section shows us how to always keep the solution
physically realizable, we need to describe the “elastic fence” and explain why
it works. We will always have some observationally motivated, or numerically
generated value for the turbulence in the system. Call the timescale for those
turbulent eddies $\tau_{phys}$. In principle, we would like to insert this
value of $\tau_{phys}$ in eqn. (6) and hope that the difference in pressures
always stays within the domain of physical realizability. But that hope is
almost never respected by a numerical code. Having a bounds-preserving
timestepping strategy will help. However, when the pressure anisotropy gets
too close to the “hard fence” formed by physical realizability, it is valuable
to provide instead an “elastic fence” which keeps the pressure anisotropy from
touching, or going beyond, the “hard fence”. This is provided by eqn. (28),
which in turn relies on eqns. (24), (27) and (32). Collectively, these
equations ensure that the relaxation time is made much shorter than
$\tau_{phys}$ when the system gets too close to the boundaries of the domain
of realizability.

Figure 4: Graph of $f(\beta_{\parallel},\bar{\beta})$ for
$\bar{\beta}=6.0/3.0$.

Fig. 4 shows how this works. In this figure, we have set $\bar{\beta}=6/3$. To
conserve total thermal energy, the values of $\beta_{\parallel}$ and
$\beta_{\perp}$ must move along the red dashed line given by $\bar{\beta}=6/3$
in Fig. 2b. We now plot out $f\left({{\beta_{\parallel}},\bar{\beta}}\right)$
for that situation in Fig. 4. We see that for a large range of values,
$f\left({{\beta_{\parallel}},\bar{\beta}}\right)$ is close to unity, with the
result that the relaxation time is indeed given by $\tau_{phys}$. When
$\beta_{\parallel}$ is much greater than $\bar{\beta}$, or when it is much
smaller than $\bar{\beta}$ , i.e., when it approaches either asymptote in Fig.
4, we see that $f\left({{\beta_{\parallel}},\bar{\beta}}\right)$ approaches
zero. This causes $\tau_{num}$ to become a very small value. In turn, this
causes the numerical relaxation time in eqn. (6) to become very small and the
system is made to move away from the “hard fence” before the numerical code
loses saliency. Please also note from Fig. 4 that we want
$f\left({{\beta_{\parallel}},\bar{\beta}}\right)$ to remain close to unity for
a large range of values. This is ensured if we take $\kappa=2$ or $\kappa=4$,
and also $\eta=2$ or $\eta=4$ in eqn. (32). We also see that setting $\eta=1$
and $\kappa=1$ is a sub-optimal choice. As long as the prescription from the
previous Sub-section is used, $\beta_{\parallel}$ will never be equal to the
dashed asymptotic values shown in Fig. 4, so it will never become exactly
zero. For all the examples shown in this paper, we have used $\eta=4$ and
$\kappa=4$.

###  5.3 Overall Implementation

The previous two Sub-sections have shown us how to ensure two things. First,
Sub-section 5.1 has shown us that ${{p}_{\parallel}}$ and ${{p}_{\bot}}$ can
always be made to lie in the physical region, even if perchance the numerical
errors may take them out of the physical region. This should occur very, very
rarely in a code. Second, Sub-section 5.2 showed that we have a physics-based
relaxation strategy that always ensures that under normal operation the
numerics remain safely within the region of realizability. For the
timestepping we will use one of the two Runge-Kutta schemes described in
subsection 3.1. Here, we provide a step-by-step implementation plan for a
single stage in the Runge-Kutta timestepping:-

  1. 1.

Start with the zone-centered conserved variables and face-centered magnetic
field variables. The zone centered variables are reconstructed using finite
volume WENO methods; see for example, the Supplement for Balsara et al.
(2023). The divergence-free magnetic field can be reconstructed using the
methods described in Balsara et al. (2018). (Any reasonable TVD or WENO scheme
can be used for the reconstruction of zone-centered fluid variables at second
order. At third order, we reconstructed the zone-centered fluid variables by
using Sub-section A.1 from the Supplement in Balsara et al. (2023). For the
divergence-free magnetic field, we used the second order reconstruction from
Section III.2 of Balsara (2001). At third order, we reconstructed the magnetic
field by using Section 2 of Balsara (2009) or Section 3 of Balsara et al.
(2018).)

  2. 2.

The reconstruction from the previous step gives us non-linearly hybridized
polynomials for the variation of all the flow variables within each zone. Now
focus on Fig. 1. Using the reconstructed variables within each zone, we can
use those polynomials to obtain the nodal values within each zone. We can also
use those nodal values to evaluate the fluxes at each nodal point. This
enables us to evaluate the flux terms, which form the third and fourth rows of
eqn. (14). This is done for each zone. We can also obtain the gradient of the
solution at each nodal point. Along with the matrices for the non-conservative
products, we can use the gradients to evaluate the fifth row of eqn. (14).
This is done for each zone.

  3. 3.

As a result of the nodal evaluation in the previous step, at each zone
boundary we have two nodal values on either side of that zone boundary. These
can be fed into the one-dimensional Riemann solver from Dumbser & Balsara
(2016) to get fluctuations that propagate to either side of each zone
boundary. This allows us to evaluate the first two rows of eqn. (14). We see,
therefore, that the entire time rate of change for the zone-centered
variables, as given by eqn. (14), is now in hand.

  4. 4.

Also as a result of the nodal evaluation in step 2, we have four states at
each edge. These are fed into a two-dimensional Riemann solver to get electric
fields at each edge. Balsara (2010, 2012, 2014) or Balsara & Nkonga (2017)
describe such multidimensional Riemann solvers. The resulting electric fields
at the edges of the mesh can then be used to obtain the time rate of change
for the facial magnetic field components.

  5. 5.

This completes our description of how $\mathcal{L}(\bf{U})$ is obtained for
use in eqn. (12) or eqn. (13). (Eqn. (12) will provide a temporally second
order accurate update; eqn. (13) will provide a temporally third order
accurate update.) Note too that at the beginning of each complete timestep we
evaluate $\tau_{num}$ within each zone using eqns. (28) and (32). This
$\tau_{num}$ can differ from zone to zone, and it can indeed change from one
complete timestep to another, but it is kept frozen within each zone for the
duration of an entire timestep. This is so that it can be used in eqn. (12) or
eqn. (13), which only remain salient if $\tau_{num}$ is kept frozen within
each zone for the duration of that timestep.

##  6 Accuracy tests

In Balsara (2004), a two-dimensional vortex problem was introduced within the
framework of the MHD system. In this section, we design such a test problem
for the CGL system. This problem is especially valuable for accuracy testing
due to its unique characteristics, featuring a smoothly evolving and
dynamically stable vortex that moves diagonally within the computational
domain.

The problem is defined within a two-dimensional domain that spans the region
$[-5,5]\times[-5,5]$. Periodic boundaries have been implemented at all the
boundaries. The vortex is initiated at the center of the computational domain
through fluctuations in the velocity and magnetic fields. These fluctuations
are defined as follows:

$\displaystyle(\delta v_{x},\delta v_{y})=\zeta\ e^{0.5(1-r^{2})}(-y,x),$ (37)
$\displaystyle(\delta B_{x},\delta B_{y})=\mu\ e^{0.5(1-r^{2})}(-y,x),$ (38)

where $\zeta$ and $\mu$ are positive parameters that will be defined later.
The distance from the center of the domain is denoted by $``r"$, which is
defined as the square root of the sum of squares of the coordinates, i.e.,
$r=\sqrt{x^{2}+y^{2}}$. It is important to emphasize that the magnetic vector
potential plays a crucial role in ensuring that the magnetic field is
initialized in a divergence-free manner within the computational domain. The
expression for the magnetic vector potential, denoted as $A_{z}$, is given by

$A_{z}=\mu\ e^{0.5(1-r^{2})}.$ (39)

As in the original MHD vortex from Balsara (2004), the unperturbed density,
pressures, and velocities are taken to be constant. Moreover, for a given
constant $\sigma$, we assume the following pressure relation for the
perturbations $\delta p_{\parallel}$ and $\delta p_{\perp}$

$\delta p_{\parallel}-\delta p_{\perp}=\sigma^{2}\ e^{(1-r^{2})}r^{2}.$ (40)

For the CGL system, the equation for the dynamical balance is given by

$\displaystyle\frac{\partial}{\partial r}\left(\delta
p_{\perp}+\frac{{B}^{2}}{8\pi}\right)=\frac{1}{r}\left(\rho{v}^{2}+\delta
p_{\parallel}-\delta p_{\perp}-\frac{{B}^{2}}{4\pi}\right).$ (41)

Assuming that $\delta p_{\perp}$ is zero at the center of the domain, we solve
equation (41) to obtain

$\displaystyle\delta p_{\perp}(r)=\frac{\mu^{2}}{8\pi}\ (1-r^{2})\
e^{(1-r^{2})}-\frac{\mu^{2}}{8\pi}e-\frac{1}{2}\left(\rho\zeta^{2}+\sigma^{2}\right)\left(e^{1-r^{2}}-e\right).$
(42)

Consequently, using equation (40), we find

$\displaystyle\delta p_{\parallel}=\delta p_{\perp}+\sigma^{2}\
e^{(1-r^{2})}r^{2}.$ (43)

As in the MHD vortex, the unperturbed density and pressures ($p_{\perp}$ and
$p_{\parallel}$) are set to unity. The unperturbed $x$ and $y$ velocities are
set to unity, and the unperturbed magnetic field is set to zero. No
perturbations have been applied to the density variable. Collectively, we have

$\left(\rho,v_{x},v_{y},v_{z},p_{\parallel},p_{\perp},B_{x},B_{y},B_{z}\right)=\left(1,1+\delta
v_{x},1+\delta v_{y},0,1+\delta p_{\parallel},1+\delta p_{\perp},\delta
B_{x},\delta B_{y},0\right).$ (44)

The test problem was configured with the following parameters:
$\mu=\sqrt{4\pi}/(2\pi)$, $\zeta=1/(2\pi)$, and $\sigma=1/(4\pi)$. We remark
that setting $\sigma=0$ retrieves analytic expressions of the MHD vortex
problem that are available in Balsara (2004). We simulate the vortex problem
till time $t=10$ on different grid stacks to perform the accuracy study. Table
(1) shows the obtained $L_{1}-$ and $L_{\infty}-$errors for the
$B_{y}-$component of the magnetic field $\mathbf{B}$. The corresponding
numerical orders are also given in the same table. For the fourth-order
scheme, we double the computational domain and stopping time to minimize the
effect of small jumps in the velocity and magnetic field at the periodic
boundaries. We see that all the schemes (second, third, and fourth-order)
reach their design accuracy.

Order of the scheme | Grid size |  $L_{1}-$errors | Order |  $L_{\infty}-$errors | Order
---|---|---|---|---|---
Second order | $32\times 32$ | 1.45602E-02 | – | 2.26862E-01 | –
$64\times 64$ | 4.04014E-03 | 1.85 | 5.94939E-02 | 1.93
$128\times 128$ | 9.14568E-04 | 2.14 | 1.61236E-02 | 1.88
$256\times 256$ | 2.25164E-04 | 2.02 | 4.54362E-03 | 1.83
Third order | $32\times 32$ | 6.98221E-03 | – | 1.15328E-01 | –
$64\times 64$ | 1.42498E-03 | 2.29 | 2.28921E-02 | 2.33
$128\times 128$ | 2.04612E-04 | 2.80 | 4.73716E-03 | 2.27
$256\times 256$ | 2.74924E-05 | 2.90 | 7.50811E-04 | 2.66
Fourth order | $32\times 32$ | 5.40291E-03 | – | 2.97769E-01 | –
$64\times 64$ | 6.34923E-04 | 3.09 | 4.87071E-02 | 2.61
$128\times 128$ | 2.81249E-05 | 4.50 | 2.99674E-03 | 4.02
$256\times 256$ | 1.70317E-06 | 4.05 | 1.95570E-04 | 3.94
Table 1: Accuracy tests: $L_{1}$ errors, $L_{\infty}$ errors, and the
corresponding order of convergence for the second, third, and fourth-order
schemes.

##  7 Some example problems

In this Section, we upgrade several standard two-dimensional MHD test problems
to show their corresponding CGL variants. For schemes characterized by spatial
second-order accuracy, temporal evolution is governed by a second-order time-
stepping algorithm specified in eqn. (12). The schemes exhibiting spatial
third and fourth-order accuracy employ a third-order time-stepping strategy
detailed in eqn (13). For all the simulations presented in this section, we
used a CFL number of 0.4.

###  7.1 Rotor Problem for the CGL Equations

Figure 5: Rotor Problem for the CGL Equations: The panels show the plots of
$\rho$, $p_{\parallel}$, $p_{\perp}$, $\Delta p/\bar{p}$, $|\bf{v}|$ and
$|\bf{B}|$ that correspond to $\tau_{phys}=10^{-5}$. Figure 6: Rotor Problem
for the CGL Equations: The panels show the plots of $\rho$, $p_{\parallel}$,
$p_{\perp}$, $\Delta p/\bar{p}$, $|\bf{v}|$ and $|\bf{B}|$ that correspond to
$\tau_{phys}=10^{-2}$. Figure 7: Rotor Problem for the CGL Equations: The
panels show the plots of $\rho$, $p_{\parallel}$, $p_{\perp}$, $\Delta
p/\bar{p}$, $|\bf{v}|$ and $|\bf{B}|$ that correspond to $\tau_{phys}=0.1$.

In this sub-section, we consider a CGL version of the two-dimensional rotor
problem that was first presented for MHD in Balsara & Spicer (1999) and
Balsara (2004). The computational domain spans the domain
$[-0.5,0.5]\times[-0.5,0.5]$. A dense and rapidly spinning cylinder is set up
in the center of an initially stationary, light ambient fluid. The ambient
fluid is initially at rest. A uniform magnetic field initially threads the two
fluids. Its value is set to 2.5 units, and it initially points in the
x-direction. The total pressure in both fluids is set to unity, i.e.
$\bar{p}=1.0$. The initial parallel and perpendicular pressures are assumed to
be equal to the total pressure, i.e. $p_{\parallel}=p_{\perp}=\bar{p}$. The
density in the ambient fluid is uniformly set to unity, while the constant
density in the rotor is $10$ units out to a radius of $0.1$. A linear taper is
applied to the density between a radius of $0.1$ and $0.13$ so that the
density in the rotor decreases linearly to the value of the density in the
ambient fluid. Six zones are used for the taper to join the density of the two
fluids. That number should be kept fixed if the resolution is increased or
decreased. The initial angular velocity of the rotor is uniform out to a
radius of $0.1$. At this radius, the toroidal velocity has a value of one
unit. The toroidal velocity decreases linearly from one unit to zero between a
radius of 0.1 and 0.13 so that it joins the velocity of the ambient fluid at a
radius of 0.13. The simulations were run until a final time of $t=0.29$ using
a mesh consisting of $200\times 200$ zones. We use the fourth-order accurate
scheme for this test problem and show the results for various $\tau_{phys}$ in
Figs. 5-7.

Fig. 5 corresponds to $\tau_{phys}=10^{-5}$ with the result that our physical
relaxation time is substantially smaller than the timestep in this problem. As
a result, we see that the pressures relax to their isotropic values within one
timestep; so the MHD limit is practically reached. Consequently, we see that
this simulation looks very much like the classical MHD case with an isotropic
pressure. Fig. 5d also shows that the relative anisotropy in the pressures is
minimal. Fig. 6 corresponds to $\tau_{phys}=10^{-2}$ with the result that our
physical relaxation time is much larger than the timestep in this problem. We
expect to see substantial differences in $p_{\parallel}$ and $p_{\perp}$, and
Fig. 6d for the relative anisotropy in the pressures shows a much larger
variation than Fig. 5d. Fig. 7 corresponds to $\tau_{phys}=0.1$ with the
result that our physical relaxation time is comparable to the final simulation
time in this problem. Now Fig. 7d shows extreme anisotropies in the pressure;
please compare it to Fig. 6d and Fig. 5d.

In this problem, the parallel and perpendicular pressures are identical at the
start of the simulation and the magnetic field is initially horizontal. We see
that the pressure anisotropies have been generated spontaneously by the
governing CGL equations. In Sub-section 1.2, we established some scalings for
the parallel and perpendicular temperatures, which can also be extended to
refer to scalings for the parallel and perpendicular pressures. By looking at
the top and bottom of the rotor in Figs. 6c and 7c, we see that the outward
propagating wave is propagating into the horizontal field, causing a
substantially greater increase in the perpendicular pressure compared to the
parallel pressure. However, we can also focus on the horizontal ends of the
rotor in Figs. 6b and 7b, which show the parallel pressure. There, we see that
the rotation is stretching the field lines longitudinally with the result that
the parallel pressure has increased much faster than the perpendicular
pressure.

###  7.2 Blast Problem for the CGL Equations

Figure 8: Blast Problem for the CGL Equations: The panels show the plots of
$\rho$, $p_{\parallel}$, $p_{\perp}$, $\Delta p/\bar{p}$, $|\bf{v}|$ and
$|\bf{B}|$ that correspond to $\tau_{phys}=10^{-5}$. Figure 9: Blast Problem
for the CGL Equations: The panels show the plots of $\rho$, $p_{\parallel}$,
$p_{\perp}$, $\Delta p/\bar{p}$, $|\bf{v}|$ and $|\bf{B}|$ that correspond to
$\tau_{phys}=10^{-3}$. Figure 10: Blast Problem for the CGL Equations: The
panels show the plots of $\rho$, $p_{\parallel}$, $p_{\perp}$, $\Delta
p/\bar{p}$, $|\bf{v}|$ and $|\bf{B}|$ that correspond to
$\tau_{phys}=10^{-2}$. Figure 11: Blast Problem for the CGL Equations: Panels
show the cut plots (along $y=0$) of $\rho$, $p_{\parallel}$, and $p_{\perp}$
profiles for various values of the relaxation time
$(\tau_{phys}=10^{-2},10^{-3}\text{ and }10^{-5})$. Figure 12: Blast Problem
for the CGL Equations: The panels show the plots of $\rho$, $p_{\parallel}$,
$p_{\perp}$ and $|\bf{v}|$ for the case where the initial magnetic field is
set to zero (${\bf{B}}=0.0$).

We consider a CGL variation of the blast problem that was initially introduced
for MHD in Balsara & Spicer (1999). The problem is configured on a unit square
grid that consists of $200\times 200$ zones and spans the domain
$[-0.5,0.5]\times[-0.5,0.5]$. Initially, the density is uniformly set to unity
across the domain. The parallel and perpendicular pressures are assumed to be
equal and are uniformly set to 0.1 except within a central circle of radius
0.1, where they are elevated to 1000. Initially, the velocity is zero
throughout. Additionally, a magnetic field with a magnitude of 60 is
initialized along the $x-$direction. We simulate the problem until a final
time of $t=0.01$ using the second-order scheme and show results for various
$\tau_{phys}$ in Figs. 8–10. For a better understanding of the shock-profiles,
in Fig. 11, we plot horizontal cuts (along $y=0$) of the density, parallel
pressure, and perpendicular pressure profiles for various values of
$\tau_{phys}$.

Fig. 8 corresponds to $\tau_{phys}=10^{-5}$ with the result that our physical
relaxation time is substantially smaller than the timestep in this problem. We
see, therefore, that the parallel and perpendicular pressures appear to be
very isotropic. Fig. 9 corresponds to $\tau_{phys}=10^{-3}$ with the result
that our physical relaxation time is much larger than the timestep in this
problem. Fig. 10 corresponds to $\tau_{phys}=10^{-2}$ so that our physical
relaxation time is comparable to the final simulation time in this problem. We
see that with increasing $\tau_{phys}$ we have an increase in the pressure
anisotropy, as expected.

As in the rotor example, the parallel and perpendicular pressures are
identical at the start of the simulation, and the magnetic field is initially
horizontal. As in the rotor problem, the pressure anisotropies have been
generated spontaneously by the physics of the governing CGL equations. By
looking at the top and bottom of the rotor in Figs. 9c and 10c, we see that
the outward propagating wave is propagating into the horizontal field, causing
a substantially greater increase in the perpendicular pressure compared to the
parallel pressure. The blast has also compressed the inner part of the domain,
with the result that the center of the blast also shows an enhanced
perpendicular pressure. We can also focus on the horizontal ends of the blast
wave in Figs. 9b and 10b, to observe that the parallel pressure is much
enhanced there. The same behavior of the pressure profiles can be seen in the
cut-plots given in Figs. 11b and 11c.

In Sub-section 4.3, we showed that the CGL equations, coupled with the concept
of an ”elastic fence” can smoothly produce an isotropic pressure in the limit
where the magnetic field tends to (and reaches) zero in a code. Here, we
numerically demonstrate this by setting the magnetic field to zero
$(\mathbf{B}=0)$ in the above test case. We present the results in Fig. 12 and
observe that our scheme yields identical parallel and perpendicular pressures,
and all profiles exhibit a strong resemblance to the hydrodynamical limit.
This test problem, numerically, illustrates that the CGL equations can indeed
operate stably in the hydrodynamic limit.

###  7.3 Orszag-Tang Problem for the CGL Equations

Figure 13: Orszag-Tang Problem for the CGL Equations: The panels show the
plots of $\rho$, $p_{\parallel}$, $p_{\perp}$, $\Delta p/\bar{p}$, $|\bf{v}|$
and $|\bf{B}|$ that correspond to $\tau_{phys}=10^{-4}$. Figure 14: Orszag-
Tang Problem for the CGL Equations: The panels show the plots of $\rho$,
$p_{\parallel}$, $p_{\perp}$, $\Delta p/\bar{p}$, $|\bf{v}|$ and $|\bf{B}|$
that correspond to $\tau_{phys}=1$.

The Orszag-Tang problem is usually taken to be a model for the transition to
MHD turbulence. The problem was first introduced in Orszag & Tang (1979) (also
see Balsara (1998)). The problem is initialized on a square domain
$[0,2]\times[0,2]$ with the periodic boundaries. The density was initialized
to unity all over the domain. The parallel and perpendicular pressures are
assumed to be equal and are uniformly set to $5/3$. The velocity field is set
up as follows:

$\mathbf{v}=-\sin(\pi y)\hat{x}+\sin(\pi x)\hat{y}.$

The magnetic field is initialized in a divergence-free manner using the
following vector potential:

$A_{z}=-\dfrac{\sqrt{4\pi}}{2\pi}\left(\cos(2\pi x)+2\cos(\pi y)\right).$

We run simulations until a final time of $t=1$ on a grid of $200\times 200$
zones. We use third-order scheme, and show results for various $\tau_{phys}$
in Figs. 13–14.

Fig. 13 shows that the pressures are indeed isotropic when the relaxation time
is very small ($\tau_{phys}=10^{-4}$). Fig. 14 corresponds to a much larger
relaxation time ($\tau_{phys}=1$). In that case we see the development of very
anisotropic pressures. This gives us the perspective that in astrophysical and
space plasmas where pressure anisotropies can develop, they will indeed
contribute significantly to the structure of the turbulence. Since the
pressure sets the temperature, and since the temperature regulates the
emissivity of a plasma, this indicates that plasmas that are simulated with
the CGL equations should, under the right conditions, show larger emissivity
fluctuations compared to plasmas that are evolved with the MHD equations.

###  7.4 Field Loop Advection Problem for the CGL Equations

Figure 15: Field Loop Advection Problem for the CGL Equations: This figure
shows the plot of $|\mathbf{B}|^{2}/(4\pi)$ that corresponds to
$\tau_{phys}=10^{-1}$.

This problem is set up on a $128\times 64$ zone domain that spans
$[-1,1]\times[-0.5,0.5]$. The problem consists of advecting a two-dimensional
loop of the magnetic field with a very low magnetic pressure compared to the
gas pressure. The loop is advected along the diagonal of the computational
domain with periodic boundaries applied in both directions. It was initially
described in Gardiner & Stone (2005), and a variation of this problem in the
CGL framework has been done in Hirabayashi et al. (2016). We follow the set-up
described in Hirabayashi et al. (2016), and do not repeat the description
here. The conventional result examines the structure of the magnetic loop
after it has completed one complete orbit around the domain. The problem was
run to a stopping time of unity using fourth-order accurate scheme.

Fig. 15 shows the magnitudes of the magnetic field for the field loop
advection problem after it has executed one complete orbit along the diagonal
of the computational domain. We use an intermediate relaxation time
$\tau_{phys}=10^{-1}$ for this problem. We see that there is virtually no
diffusion of the loop’s boundaries and no oscillations in the magnetic
pressure within the loop’s interior. We, therefore, conclude that the method
presented here has adequate dissipation behavior for the advection of magnetic
fields.

###  7.5 Shock interaction with a magnetosphere

Figure 16: Shock interaction with a magnetosphere: The schematic figure
illustrates the setup of a two-dimensional test problem that very roughly
imitates a shock interacting with a magnetosphere. The domain limits are shown
in the figure. A cylinder of unit radius that is threaded on the inside with a
magnetic field with magnitude $B_{0}$ is also shown. The parameters that are
initially used for the cylinder and the field are shown in the figure. The
exterior field that threads the cylinder is intended to look like a dipolar
field in 2D. A uniform magnetic field with $B_{1}=0.1$ threads the region that
makes up the ambient medium around the cylinder. The post-shock variables are
also shown in the figure. The right-going fast magnetosonic shock eventually
overruns the magnetized cylinder.

In this sub-section, we consider a two-dimensional test problem that very
roughly imitates a shock interacting with a magnetosphere. The two-dimensional
computational domain spans
$\left(x,y\right)\in\left[-20,40\right]\times\left[-10,10\right]$. The set-up
data is schematically shown in Fig. 16. A cylinder with a unit radius is set
up at the origin. The density inside the cylinder is ${{\rho}_{0}}=300$ with
pressures within the cylinder initialized as ${{p}_{\parallel 0}}={{p}_{\bot
0}}=1$. A constant magnetic field is initialized within the cylinder in the
y-direction with ${{B}_{0}}=3.5$. This cylinder is embedded in a much lighter
ambient fluid with density ${{\rho}_{1}}=1$ and pressures ${{p}_{\parallel
1}}={{p}_{\bot 1}}=1$. The fluid that makes up the cylinder and its ambient
are initially not moving. External to the cylinder, we set up a magnetic
structure that looks like a dipole in 2D with magnetic field given by

$\mathbf{B}=(B_{x},B_{y},B_{z})^{\top}=\left(B_{0}\dfrac{2xy}{(x^{2}+y^{2})^{2}},B_{0}\dfrac{(y^{2}-x^{2})}{(x^{2}+y^{2})^{2}},0\right)^{\top}$
(45)

The corresponding vector potential, which is useful for initializing the
problem, is given by

$\displaystyle{\vec{A}}_{0}(x,y)=\begin{cases}-{{B}_{0}}x\hat{z}&\text{for}\
{{x}^{2}}+{{y}^{2}}\leq 1\\\
-\dfrac{{{B}_{0}}x}{{{x}^{2}}+{{y}^{2}}}\hat{z}&\text{for}\
{{x}^{2}}+{{y}^{2}}>1\end{cases}$ (46)

The magnetic field that is set up external to the cylinder is designed to
dovetail with the magnetic field within the cylinder. Because we intend to set
up a right-going magnetosonic shock that eventually overruns the cylinder, we
also add to the cylinder’s internal field, as well as to its ambient magnetic
field, an additional constant magnetic field ${{\vec{B}}_{1}}=0.1\hat{y}$.
This ensures that the shocked gas will have an even larger y-component of the
magnetic field. This right-going fast magnetosonic shock is initialized along
the line $x=-15$. The shocked gas has density ${{\rho}_{2}}=3.334$, pressures
${{p}_{\parallel 2}}={{p}_{\bot 2}}=18.537$, velocity
${{{\vec{v}}}_{2}}=3.5043\hat{x}$ and magnetic field
${{\vec{B}}_{2}}=0.33342\hat{y}$. The boundary conditions are periodic in the
y-direction; inflow at the left boundary and outflow at the right boundary.
The problem was run to a simulation time of 10 units, by which time the shock
has overrun the dense cloud and is close to the right boundary. A mesh of
$1200\times 400$ zones was used.

Figure 17: Shock interaction with a magnetosphere: Panel shows the plots of
$\rho$, $p_{\parallel}$, $p_{\perp}$, $\Delta p/\bar{p}$, $|\bf{v}|$ and
$|\bf{B}|$ that correspond to $\tau_{phys}=10^{-4}$. Figure 18: Shock
interaction with a magnetosphere: Panel shows the plots of $\rho$,
$p_{\parallel}$, $p_{\perp}$, $\Delta p/\bar{p}$, $|\bf{v}|$ and $|\bf{B}|$
that correspond to $\tau_{phys}=1$.

Figs. 17 and 18 show two versions of this problem with
${{\tau}_{phys}}={{10}^{-4}}$ (i.e. the MHD limit) and ${{\tau}_{phys}}=1$
(i.e. the CGL limit). Several common features are shared between the two
simulations, as expected, but we can also pick out prominent differences.
Compare Figs. 17d and 18d. We see that there is substantially more pressure
anisotropy in Fig. 18d as compared to Fig. 17d. From Fig. 18d we see that the
part of the bow shock that faces into the incoming flow displays a
${{p}_{\bot}}$ that is substantially larger than ${{p}_{\parallel}}$, which is
consistent with observations and also physical expectations. In the
magnetotail region we see parts of the flow with
${{p}_{\bot}}>{{p}_{\parallel}}$ as well as parts of the flow with
${{p}_{\bot}}<{{p}_{\parallel}}$ indicating that the pressure fluctuations in
the CGL limit could strongly contribute to the turbulence in the magnetotail.
This is just a 2D test problem, so it is not expected to produce fully-
developed turbulence. Even so, comparing the velocity plots in Figs. 17e and
18e, we see that the velocity structure in the magnetotail region shows
substantially enhanced turbulence in Fig. 18e compared to Fig. 17e. Despite
being a 2D simulation, we see several very significant trends that we will
follow up in 3D and with more realistic conditions.

##  8 Conclusions

MHD is built on assuming that there is a single isotropic pressure. Several
lines of observational evidence, stemming from the solar heliosphere and the
magnetospheres of planets in our solar system, have demonstrated that the
plasma pressures can often become anisotropic. Astrophysical systems, like
accretion disks and the intercluster medium, also bring out the need for
modeling pressure anisotropies in astrophysical plasmas. The CGL system is the
simplest example of a hyperbolic PDE system that allows plasma pressures to
become anisotropic.

We have made a detailed study of the CGL system and the challenges it poses.
At a numerical level:- 1) the system is in non-conservation form, 2) it has
stiff source terms, 3) it is prone to a loss of hyperbolicity. Furthermore,
since different physical instabilities come into play when
${{p}_{\parallel}}>{{p}_{\perp}}$, and when ${{p}_{\parallel}}<{{p}_{\perp}}$,
it is imperative that the operation of the stiff source terms also preserve
these inequalities. The numerical challenges cannot be solved without detailed
consideration of the physics of the instabilities that dominate a
collisionless plasma. For that reason, we realize that when the pressure
anisotropy becomes too large, and with ${{p}_{\parallel}}>{{p}_{\perp}}$, the
firehose instability will dominate. On the other hand, when the pressure
anisotropy becomes too large, and with ${{p}_{\parallel}}<{{p}_{\perp}}$, the
mirror instability will dominate. These instabilities can operate on the
microscales of a large-scale computation and drive the system towards
isotropy. We realize, therefore, that physics and numerics need to work hand-
in-glove in order to make the CGL equations numerically tractable. To that
end, we first identify the limits to the pressure anisotropy past which the
firehose and mirror instabilities act. To include that physics in the
numerics, we invent the idea of an “elastic fence”. It ensures that if the
system approaches the “hard fence” associated with these instabilities, it is
gently pushed away from hitting these limits. Indeed, in order to function
reliably, a numerical code must stay away from these limits because if it did
not, hyperbolicity would be violated and the code would crash immediately
because it would not be able to evaluate the physical eigenvalues. As
inability to evaluate eigenvalues translates into an inability to evaluate the
next timestep in any code for the CGL equations. We see, therefore, that any
numerical code for the CGL equations should stay safely clear of the “hard
fence” by resorting to the “elastic fence” that we propose.

To accommodate this blend of physics-based and numerically-driven
requirements, we invent a bounds-preserving time-integration strategy that is
stable even in the presence of ultra-stiff source terms. The bounds-preserving
aspect of the scheme ensures that the sign of
$\left({{p}_{\parallel}}-{{p}_{\perp}}\right)$ is preserved by an application
of the stiff source terms. (Note, however, that the non-conservative products
in the equation for
${\partial\left({{p}_{\parallel}}-{{p}_{\perp}}\right)}/{\partial t}\;$ can
indeed drive the system away from isotropy and we design accurate methods that
allow that to happen if the physical problem has to go in that direction.) We
also recast the system in fluctuation form so that the system retains
conservation when that is warranted while accommodating the presence of non-
conservative products. The “elastic fence” is, therefore, a way of emulating
in a numerical code the same thing that nature chooses to do via the operation
of the firehose and mirror instabilities.

Several practical ideas are presented for conserving thermal energy while
keeping the pressure anisotropy within physically acceptable bounds. In
practice, it means that the CGL equations by themselves do not have the
requisite “knowledge” to stay within the domain of physical realizability, and
our innovations provide the extra “guidance” to ensure that this is achieved.
Details on the bounds-preserving timestepping, the fluctuation form, the
Riemann solvers that support this fluctuation form, and the design of the
“elastic fence” are presented in this work. Section 5 brings all these
innovations together by explaining how the different steps are integrated into
numerical code.

A completely non-trivial test problem is constructed in Section 6 to show that
the proposed numerical methods for the CGL equations meet their design
accuracy. We demonstrate that our numerical code can reach its designed
second, third, and fourth-order accuracy. We also present several stringent
numerical examples. The test problems illustrate that when the physical
relaxation time is small, pressure isotropy is indeed very closely achieved by
the CGL equations – in other words, the CGL equations can access the MHD limit
if our proposed methods are used. However, when the physical relaxation time
is large, the CGL equations show us the different types of physics that arise
when pressure isotropy is violated. It has also been claimed in the literature
that the CGL equations cannot reclaim the hydrodynamic limit because
$\mathbf{b}={\mathbf{B}}/{\left|\mathbf{B}\right|}\;$ becomes undefined. But
we show that this claim is specious in light of the fact that the physics of
the “elastic fence” restores pressure isotropy when the plasma beta becomes
extremely large. As a result, we demonstrate, perhaps for the first time, that
a CGL code will stably operate in the hydrodynamical limit. We also include a
test problem that mimics the interaction of a shock with a magnetospheric
environment in 2D.

##  9 Acknowledgments

This work was supported by NSF grant 2009871 and by NASA grant 80NSSC20K0786.

## References

  * Achatz et al. (1991) Achatz, U., Steinacker, J., & Schlickeiser, R. 1991, A&A, 250, 266
  * Alfvén (1943) Alfvén, H. 1943, Arkiv. Mat. Astron. Fys
  * Balsara (1998) Balsara, D. S. 1998, The Astrophysical Journal Supplement Series, 116, 133
  * Balsara (2001) —. 2001, Journal of Computational Physics, 174, 614
  * Balsara (2004) —. 2004, The Astrophysical Journal Supplement Series, 151, 149
  * Balsara (2009) —. 2009, Journal of Computational Physics, 228, 5040
  * Balsara (2010) —. 2010, Journal of Computational Physics, 229, 1970
  * Balsara (2012) —. 2012, Journal of Computational Physics, 231, 7476
  * Balsara (2014) —. 2014, Journal of Computational Physics, 277, 163
  * Balsara et al. (2018) Balsara, D. S., Garain, S., Taflove, A., & Montecinos, G. 2018, Journal of Computational Physics, 354, 613
  * Balsara & Nkonga (2017) Balsara, D. S., & Nkonga, B. 2017, Journal of Computational Physics, 346, 25
  * Balsara et al. (2023) Balsara, D. S., Samantaray, S., & Subramanian, S. 2023, Communications on Applied Mathematics and Computation, 5, 428
  * Balsara & Spicer (1999) Balsara, D. S., & Spicer, D. S. 1999, Journal of Computational Physics, 149, 270
  * Baranov & Zaitsev (1995) Baranov, V. B., & Zaitsev, N. A. 1995, A&A, 304, 631
  * Barnes (1966) Barnes, A. 1966, The Physics of Fluids, 9, 1483
  * Blandford & Payne (1982) Blandford, R. D., & Payne, D. 1982, Monthly Notices of the Royal Astronomical Society, 199, 883
  * Bruno & Carbone (2016) Bruno, R., & Carbone, V. 2016, Turbulence in the Solar Wind, Vol. 928, doi: 10.1007/978-3-319-43440-7
  * Burgers (1969) Burgers, J. M. 1969, Flow Equations for Composite Gases
  * Cairns & Fuselier (2017) Cairns, I. H., & Fuselier, S. A. 2017, ApJ, 834, 197, doi: 10.3847/1538-4357/834/2/197
  * Chandrasekhar et al. (1958) Chandrasekhar, S., Kaufman, A., & Watson, K. 1958, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences, 245, 435
  * Chew et al. (1956) Chew, G. F., Goldberger, M. L., & Low, F. E. 1956, Proceedings of the Royal Society of London Series A, 236, 112, doi: 10.1098/rspa.1956.0116
  * Crooker et al. (1979) Crooker, N. U., Eastman, T. E., & Stiles, G. S. 1979, J. Geophys. Res., 84, 869, doi: 10.1029/JA084iA03p00869
  * Denton & Lyon (1996) Denton, R. E., & Lyon, J. G. 1996, Geophys. Res. Lett., 23, 2891, doi: 10.1029/96GL01590
  * Denton & Lyon (2000) —. 2000, J. Geophys. Res., 105, 7545, doi: 10.1029/1999JA000360
  * Dumbser & Balsara (2016) Dumbser, M., & Balsara, D. S. 2016, Journal of Computational Physics, 304, 275
  * Erkaev et al. (1999) Erkaev, N. V., Farrugia, C. J., & Biernat, H. K. 1999, J. Geophys. Res., 104, 6877, doi: 10.1029/1998JA900134
  * Gardiner & Stone (2005) Gardiner, T. A., & Stone, J. M. 2005, Journal of Computational Physics, 205, 509
  * Gary et al. (1993) Gary, S. P., Fuselier, S. A., & Anderson, B. J. 1993, J. Geophys. Res., 98, 1481, doi: 10.1029/92JA01844
  * Gary et al. (2001) Gary, S. P., Skoug, R. M., Steinberg, J. T., & Smith, C. W. 2001, Geophys. Res. Lett., 28, 2759, doi: 10.1029/2001GL013165
  * Gershman et al. (2013) Gershman, D. J., Slavin, J. A., Raines, J. M., et al. 2013, Journal of Geophysical Research (Space Physics), 118, 7181, doi: 10.1002/2013JA019244
  * Goldreich & Sridhar (1995) Goldreich, P., & Sridhar, S. 1995, Astrophysical Journal, Part 1 (ISSN 0004-637X), vol. 438, no. 2, p. 763-775, 438, 763
  * Goldreich & Sridhar (1997) —. 1997, The Astrophysical Journal, 485, 680
  * Hall & Sturrock (1967) Hall, D. E., & Sturrock, P. A. 1967, Physics of Fluids, 10, 2620, doi: 10.1063/1.1762084
  * Hellinger et al. (2006) Hellinger, P., Trávníček, P., Kasper, J. C., & Lazarus, A. J. 2006, Geophys. Res. Lett., 33, L09101, doi: 10.1029/2006GL025925
  * Hesse & Birn (1992) Hesse, M., & Birn, J. 1992, J. Geophys. Res., 97, 10643, doi: 10.1029/92JA00793
  * Hirabayashi et al. (2016) Hirabayashi, K., Hoshino, M., & Amano, T. 2016, Journal of Computational Physics, 327, 851, doi: 10.1016/j.jcp.2016.09.064
  * Hirabayashi et al. (2016) Hirabayashi, K., Hoshino, M., & Amano, T. 2016, Journal of Computational Physics, 327, 851
  * Hoelzl et al. (2021) Hoelzl, M., Huijsmans, G. T., Pamela, S., et al. 2021, Nuclear Fusion, 61, 065001
  * Hu et al. (2010) Hu, Y., Denton, R. E., & Lin, Y. 2010, Journal of Atmospheric and Solar-Terrestrial Physics, 72, 1155–1162, doi: 10.1016/j.jastp.2010.07.007
  * Huang & Shu (2018) Huang, J., & Shu, C.-W. 2018, Journal of Computational Physics, 361, 111
  * Jokipii (1966) Jokipii, J. R. 1966, ApJ, 146, 480, doi: 10.1086/148912
  * Kato et al. (1966) Kato, Y., Tajiri, M., & Taniuti, T. 1966, Journal of the Physical Society of Japan, 21, 765
  * Kennel & Petschek (1966) Kennel, C. F., & Petschek, H. 1966, Journal of Geophysical Research, 71, 1
  * Luo et al. (2023) Luo, H., Lyon, J. G., & Zhang, B. 2023, Journal of Computational Physics, 490, 112311, doi: 10.1016/j.jcp.2023.112311
  * Mac Low et al. (2005) Mac Low, M.-M., Balsara, D. S., Kim, J., & de Avillez, M. A. 2005, The Astrophysical Journal, 626, 864
  * Marsch et al. (2004) Marsch, E., Ao, X. Z., & Tu, C. Y. 2004, Journal of Geophysical Research (Space Physics), 109, A04102, doi: 10.1029/2003JA010330
  * Marsch et al. (1981) Marsch, E., Rosenbauer, H., Schwenn, R., Muehlhaeuser, K. H., & Denskat, K. U. 1981, J. Geophys. Res., 86, 9199, doi: 10.1029/JA086iA11p09199
  * Masters et al. (2014) Masters, A., Phan, T. D., Badman, S. V., et al. 2014, Journal of Geophysical Research (Space Physics), 119, 121, doi: 10.1002/2013JA019516
  * Meng et al. (2012a) Meng, X., Tóth, G., Liemohn, M. W., Gombosi, T. I., & Runov, A. 2012a, Journal of Geophysical Research (Space Physics), 117, A08216, doi: 10.1029/2012JA017791
  * Meng et al. (2012b) Meng, X., Tóth, G., Sokolov, I. V., & Gombosi, T. I. 2012b, Journal of Computational Physics, 231, 3610, doi: 10.1016/j.jcp.2011.12.042
  * Orszag & Tang (1979) Orszag, S. A., & Tang, C.-M. 1979, Journal of Fluid Mechanics, 90, 129
  * Parker (1958) Parker, E. N. 1958, Physical Review, 109, 1874, doi: 10.1103/PhysRev.109.1874
  * Parker (2019) Parker, E. N. 2019, Cosmical magnetic fields: Their origin and their activity (Oxford university press)
  * Pogorelov & Semenov (1997) Pogorelov, N. V., & Semenov, A. Y. 1997, A&A, 321, 330
  * Pudovkin & Semenov (1977) Pudovkin, M., & Semenov, V. 1977, in Annales de Geophysique, Vol. 33, 423–427
  * Pudovkin et al. (1999) Pudovkin, M. I., Besser, B. P., Lebedeva, V. V., Zaitseva, S. A., & Meister, C.-V. 1999, Physics of Plasmas, 6, 2887, doi: 10.1063/1.873246
  * Robinson et al. (2008) Robinson, A., Brunner, T., Carroll, S., et al. 2008, in 46th aiaa aerospace sciences meeting and exhibit, 1235
  * Ruzmaikin et al. (1988) Ruzmaikin, A., Shukurov, A., & Sokoloff, D. 1988, Magnetic fields of galaxies Kluwer Academic Press, Dordrecht
  * Samsonov et al. (2001) Samsonov, A. A., Pudovkin, M. I., Gary, S. P., & Hubert, D. 2001, J. Geophys. Res., 106, 21689, doi: 10.1029/2000JA900150
  * Santos-Lima et al. (2014) Santos-Lima, R., de Gouveia Dal Pino, E. M., Kowal, G., et al. 2014, ApJ, 781, 84, doi: 10.1088/0004-637X/781/2/84
  * Santos-Lima et al. (2016) Santos-Lima, R., Yan, H., de Gouveia Dal Pino, E. M., & Lazarian, A. 2016, MNRAS, 460, 2492, doi: 10.1093/mnras/stw1079
  * Seough & Yoon (2012) Seough, J., & Yoon, P. H. 2012, Journal of Geophysical Research (Space Physics), 117, A08101, doi: 10.1029/2012JA017645
  * Sharma et al. (2006) Sharma, P., Hammett, G. W., Quataert, E., & Stone, J. M. 2006, The Astrophysical Journal, 637, 952–967, doi: 10.1086/498405
  * Sridhar & Goldreich (1994) Sridhar, S., & Goldreich, P. 1994, Astrophysical Journal, 432, 612
  * Steinolfson et al. (1975) Steinolfson, R. S., Dryer, M., & Nakagawa, Y. 1975, J. Geophys. Res., 80, 1223, doi: 10.1029/JA080i010p01223
  * Weber & Davis Jr (1967) Weber, E. J., & Davis Jr, L. 1967, Astrophysical Journal, vol. 148, p. 217-227, 148, 217
  * Yoon (1990) Yoon, P. H. 1990, Physics of Fluids B, 2, 842, doi: 10.1063/1.859568
  * Yoon & Seough (2012) Yoon, P. H., & Seough, J. 2012, Journal of Geophysical Research (Space Physics), 117, A08102, doi: 10.1029/2012JA017697
  * Zwan & Wolf (1976) Zwan, B., & Wolf, R. 1976, Journal of Geophysical Research, 81, 1636