# Nonlinear Damping and Field-Aligned Flows of Propagating Shear Alfvén Waves
with Braginskii Viscosity

Alexander J. B. Russell School of Science and Engineering,
University of Dundee,
Dundee, DD1 4HN, Scotland, UK Alexander J. B. Russell<EMAIL_ADDRESS>
###### Abstract

Braginskii MHD provides a more accurate description of many plasma
environments than classical MHD since it actively treats the stress tensor
using a closure derived from physical principles. Stress tensor effects
nonetheless remain relatively unexplored for solar MHD phenomena, especially
in nonlinear regimes. This paper analytically examines nonlinear damping and
longitudinal flows of propagating shear Alfvén waves. Most previous studies of
MHD waves in Braginskii MHD considered the strict linear limit of vanishing
wave perturbations. We show that those former linear results only apply to
Alfvén wave amplitudes in the corona that are so small as to be of little
interest, typically a wave energy less than $10^{-11}$ times the energy of the
background magnetic field. For observed wave amplitudes, the Braginskii
viscous dissipation of coronal Alfvén waves is nonlinear and a factor around
$10^{9}$ stronger than predicted by the linear theory. Furthermore, the
dominant damping occurs through the parallel viscosity coefficient $\eta_{0}$,
rather than the perpendicular viscosity coefficient $\eta_{2}$ in the
linearized solution. This paper develops the nonlinear theory, showing that
the wave energy density decays with an envelope $(1+z/L_{d})^{-1}$. The
damping length $L_{d}$ exhibits an optimal damping solution, beyond which
greater viscosity leads to lower dissipation as the viscous forces self-
organise the longitudinal flow to suppress damping. Although the nonlinear
damping greatly exceeds the linear damping, it remains negligible for many
coronal applications.

Alfvén waves (23), Solar corona (1483), Solar coronal heating (1989), Solar
coronal holes (1484), Solar wind (1534), Magnetohydrodynamics (1964), Space
plasmas (1544), Plasma astrophysics (1261), Plasma physics (2089)

##  1 Introduction

Alfvénic waves are a ubiquitous feature of natural plasmas, including the
solar corona (Tomczyk et al., 2007; De Pontieu et al., 2007; Lin et al., 2007;
Okamoto et al., 2007) and solar wind (Coleman, 1967; Belcher & Davis, 1971).
In solar physics, these waves contain sufficient energy to heat the open
corona and accelerate the fast solar wind (McIntosh et al., 2011), and they
damp significantly within a solar radius above the surface (Bemporad & Abbo,
2012; Hahn et al., 2012; Hahn & Savin, 2013; Hahn et al., 2022). How these
Alfvénic waves damp in astrophysical and space plasmas is an important
question that has remained open for almost a century (see early papers by
Alfvén 1947 and Osterbrock 1961; modern reviews by De Moortel & Browning 2015,
Arregui 2015 and Van Doorsselaere et al. 2020; and historical perspectives by
Russell 2018 and De Moortel et al. 2020).

Most theoretical knowledge about solar Alfvénic waves is based on “classical”
magnetohydrodynamics (MHD), a mathematical framework that originated from
intuitive coupling of Maxwell’s equations and Euler equations of inviscid
hydrodynamics (Hartmann, 1937; Alfvén, 1942, 1943, 1950; Batchelor, 1950) and
became widely adopted in large part due to its success providing insight into
diverse natural phenomena (see e.g. Priest, 2014). However, classical MHD is
one member of a larger family of plasma descriptions, some of which offer a
more complete description of the plasma. This paper analytically examines
Alfvén wave damping in the more general framework of Braginskii MHD, which
unlike classical MHD, retains the anisotropic viscous stress tensor.

A number of authors, including § 8 of Braginskii (1965), have previously
investigated viscous damping of Alfvén waves in the linear limit of
vanishingly small wave amplitude. When priority is given to smallness of the
wave amplitude, the problem becomes framed as a matter of how anisotropic
viscosity affects velocities that are perpendicular to the magnetic field (the
direction of which is treated as unchanging). With this approximation, damping
is determined by the “perpendicular” viscosity coefficient $\eta_{2}$, which
is extremely small in the corona. It was thus originally concluded that
viscous damping is very weak for coronal Alfvén waves unless they have very
short wavelengths. This path of reasoning is shown as the vertical branch in
Fig. 1.

There is, however, another way to view the problem. Viscous damping of Alfvén
waves can alternatively be considered with priority given to the largeness of
parallel viscosity coefficient $\eta_{0}$. Given that
$\eta_{2}/\eta_{0}\gtrsim 10^{-11}$ is typical in the corona (Hollweg, 1985),
even a very small component of $\mathbf{v}$ parallel to the total magnetic
field $\mathbf{B}$ would be expected to produce major departures from linear
theory. This path of reasoning is shown as the horizontal branch in Fig. 1.

Figure 1: Schematic paths of reasoning. The vertical branch gives priority to
smallness of the wave amplitude and concludes that damping is a linear process
governed by the perpendicular viscosity coefficient $\eta_{2}$, e.g. § 8 of
Braginskii (1965). The horizontal branch gives priority to the smallness of
$\eta_{2}/\eta_{0}$, leading to nonlinear damping via $\eta_{0}$. This paper
follows the diagonal branch, which includes deriving the validity condition
for the two outcomes. Nonlinear damping via $\eta_{0}$ is appropriate for most
coronal applications.

The second viewpoint of the problem takes impetus from the observation that
(unless wave amplitudes vanish entirely) Alfvén waves do have a non-zero
velocity component parallel to the total magnetic field. Two effects
contribute to this, which are separated if one expands
$\mathbf{V}\cdot\mathbf{B}=\mathbf{V}\cdot(\mathbf{b}+\mathbf{B}_{0})$, where
$\mathbf{B}_{0}$ is the equilibrium magnetic field and $\mathbf{b}$ is the
magnetic perturbation. First, $\mathbf{V}\cdot\mathbf{b}$ is non-zero for an
Alfvén wave, since the velocity perturbation perpendicular to $\mathbf{B}_{0}$
is aligned with the magnetic field perturbation $\mathbf{b}$. In other words,
deflection of the magnetic field from its equilibrium direction implies there
is a non-zero velocity component parallel to the total magnetic field. Second,
in compressible plasma, the magnetic pressure of the magnetic field
perturbation drives a nonlinear ponderomotive flow parallel to the equilibrium
magnetic field (e.g. Hollweg, 1971). The ponderomotive flow makes
$\mathbf{V}\cdot\mathbf{B}_{0}$ nonzero as well. Both of these effects allow
for the possibility of nonlinear viscous damping via the large parallel
viscosity coefficient $\eta_{0}$.

With the benefit of modern observations (e.g. McIntosh et al., 2011; Morton et
al., 2015), it is known that normalised wave amplitudes $b/B_{0}\sim
V/v_{A}\sim 0.1$ are typical for the base of an open coronal field region, for
example. The “smallness” of the square of this ratio is very modest in
comparison to the extreme largeness of $\eta_{0}/\eta_{2}$. Thus, it is likely
from the outset that viscous damping of Alfvén waves will be a nonlinear
process governed by $\eta_{0}$ and the wave amplitude. This paper provides
mathematical evidence that this heuristic analysis holds true, along with
detailed examination of the consequences.

Various previous studies have explored effects of Braginskii viscosity on MHD
waves since Braginskii (1965). In solar physics, the effect of linearized
Braginskii viscosity was revisited from the 1980s to the mid-1990s through the
lens of phase mixing and resonant absorption, with the aim of determining how
including the viscosity tensor modifies these scale-shortening processes and
their heating properties. At the time, it was common practice in solar MHD
wave theory to work with linearized equations. Thus, due to linearization,
Steinolfson et al. (1986), Hollweg (1987), Ruderman (1991), Ofman et al.
(1994) and Erdelyi & Goossens (1995) obtained analytical and numerical results
that strictly apply to Alfvén waves of vanishing amplitude.

In adjacent fields, the effect of anisotropic viscosity on MHD waves has also
been investigated with an eye on MHD turbulence and the solar wind. Of
particular note, Montgomery (1992) advocated that Braginskii viscosity is
important in hot tenuous plasmas, that in many circumstances it should be
treated using parallel ion viscosity, and that plasma motions may self-
organise to suppress damping. He further applied these ideas to anisotropy in
MHD turbulence, on the basis that a quasi-steady turbulence is composed of the
undamped modes. Quantitative elaboration in Montgomery (1992) was based on a
linear normal mode analysis, which captures linear damping of magnetoacoustic
waves by parallel viscosity, but excludes nonlinear viscous damping of Alfvén
waves. The conclusion that a linearized stress tensor damps Alfvén waves only
negligibly, while damping magnetoacoustic waves significantly, was further
reinforced by related work by Oughton (1996, 1997).

Similar ideas to ours regarding the importance of nonlinearity were advocated
by Nocera et al. (1986), who modelled Alfvén waves subject to the $\eta_{0}$
part of the Braginskii viscous stress tensor, retaining the leading-order
nonlinear terms in the wave perturbations. Consistent with the argument above,
their calculations found that coronal Alfvén waves damp nonlinearly by
parallel viscosity. The current paper complements and extends the previous
analysis by Nocera et al. (1986), with the goal of producing a comprehensive
understanding of the nonlinear damping and field-aligned flows of propagating
shear Alfvén waves with Braginskii viscosity.

A limitation of the mathematical techniques used in this paper is that they
exclude certain other nonlinear effects that may be important in plasmas, such
as nonlinear interactions between waves. Numerical investigations will be
required in future to verify the analytical theory presented here, compare the
relative importance of viscous damping and other nonlinear effects such as
parametric decay instability, and consider interactions between nonlinear
processes in Braginskii MHD.

This paper is organised as follows. § 2 provides scientific background on
single-fluid Braginskii MHD and its relationship to other single-fluid plasma
models. § 3 quantitatively examines Alfvén wave heating by the full Braginskii
viscous stress tensor, demonstrating the importance of nonlinear $\eta_{0}$
terms and compressibility, and obtaining the wave decay properties for the
weakly viscous limit using energy principles. In § 4, we argue that in highly
viscous limit, viscous heating is suppressed by self-organisation of the
ponderomotive flow, which implies that viscosity strongly alters the field-
aligned flow associated with Alfvén waves in this regime. § 5 further
strengthens the analysis, using multiple scale analysis to obtain the decay
properties without restrictions on the Alfvénic Reynolds number, assuming the
framework of Braginskii MHD. The paper finishes with discussion in § 6 and
summary of main conclusions in § 7.

##  2 Braginskii MHD

Braginskii MHD is an important plasma description that treats anisotropic
viscosity and thermal conduction using rigorous closure from physical
principles. This section provides a short primer on single-fluid Braginksii
MHD, its connection with pressure (or temperature) anisotropy, and its
relation to classical MHD and the CGL double-adiabatic equations.

As is described in various plasma textbooks, fluid variables can be rigorously
and robustly defined as velocity moments of the underlying particle
distribution functions. Transport equations for each particle species are then
derived by taking moments of the kinetic Boltzmann equation, and combined to
obtain the single fluid equations. Recommended presentations can be found in
Schunk & Nagy (2009) Chapter 7 and the Appendix of Spitzer (1962).

Assuming quasi-neutrality, and conservation of mass, momentum and energy, this
process yields the mass continuity equation,

$\frac{\partial\rho}{\partial t}+\nabla\cdot\left(\rho\mathbf{V}\right)=0,$
(1)

momentum equation

$\rho\frac{\mathrm{D}\mathbf{V}}{\mathrm{D}t}=-\nabla\cdot\mathbf{P}+\rho\mathbf{G}+\mathbf{j}\times\mathbf{B},$
(2)

energy equation,

$\frac{\mathrm{D}}{\mathrm{D}t}\left(\frac{3}{2}p\right)+\frac{5}{2}p\nabla\cdot\mathbf{V}=-\mathbf{\pi}:\nabla\mathbf{V}-\nabla\cdot\mathbf{q}+\mathbf{j}\cdot\left(\mathbf{E}+\mathbf{V}\times\mathbf{B}\right),$
(3)

higher-order transport equations if required, and the generalized Ohm’s law.

The pressure tensor $\mathbf{P}$ that appears in Eq. (2) is the most
fundamental representation of the internal forces associated with thermal
motions of particles. It is symmetric, so it represents six degree of freedom.
The momentum equation can also be reformulated by introducing the scalar
pressure and stress tensor as

$\displaystyle p$
$\displaystyle=\frac{1}{3}\operatorname{Trace}\left(\mathbf{P}\right)=\frac{1}{3}P_{\alpha\alpha},$
$\displaystyle\pi_{\alpha\beta}$
$\displaystyle=P_{\alpha\beta}-p\delta_{\alpha\beta},$ (4)

where $\delta_{\alpha\beta}$ is the Kronecker delta. So defined, the stress
tensor $\mathbf{\pi}$ is symmetric and traceless. These definitions gives the
replacement $-\nabla\cdot\mathbf{P}=-\nabla p-\nabla\cdot\mathbf{\pi}$.

Deriving transport equations by moment taking meets with a fundamental closure
problem: the transport equation for each fluid variable depends on a higher-
order variable, producing an infinite regress unless the system can be closed
by other considerations. The method of closure is therefore a major
distinguishing feature between different fluid models for plasmas. It is also
a major source of validity caveats. Various different methods of closure
produce governing equations that conserve mass, momentum and energy, since
these properties are already built into Eqs. (1)–(3). However, the different
models discussed below disagree on the internal forces and heating, and can
therefore produce different behaviors.

Classical MHD (Hartmann, 1937; Alfvén, 1942, 1943; Batchelor, 1950)
corresponds to a closure treatment in which the stress tensor and the heat
flow vector are dropped from Eqs. (2) and (3). Dropping the stress tensor can
be justified when particle collisions or other forms of particle scattering
such as wave-particle interactions are frequent enough that the pressure
tensor remains very close to isotropic. The resulting MHD equations are valid
for many situations, for instance modelling static equilibria, or dynamic
situations in which the divergence of the stress tensor remains small compared
to the Lorentz force. It is nonetheless a truncation since higher order
variables are set to zero rather than approximated. Furthermore,
collisionality in environments such as the solar corona is low enough that the
stress tensor can become significant for various dynamic phenomena, including
MHD waves.

Braginskii MHD uses a less restrictive method of closure. As is detailed by
Braginskii (1965), when the collisional mean free path is significantly
shorter than length scales over which fluid quantities vary, the heat flow
vector takes the form of an anisotropic thermal conduction, and the stress
tensor takes the form of an anisotropic viscosity. Closure can therefore be
achieved by expressing $\mathbf{q}$ and $\pi$ in terms of lower-order fluid
variables, which are traditionally derived using methods similar to Chapman &
Cowling (1939) or Grad (1949).

The anisotropy inherent in $\mathbf{q}$ and $\pi$ can be appreciated
heuristically, by considering the helical motion of charged particles in
magnetized plasmas. The mean free path parallel to the magnetic field is the
same as for unmagnetized plasmas, implying that transport parallel to the
magnetic field is the same as for unmagnetized plasmas. Meanwhile, the mean
free path perpendicular to the magnetic field is the gyroradius, which is
typically much less than the mean free path parallel to the magnetic field,
which supresses perpendicular transport. Hence both thermal conduction and
viscous stresses are anisotropic with respect to the magnetic field direction,
often extremely so.

The full Braginskii stress tensor, used in § 3, involves five viscosity
coefficients. A useful simplification, used in § 5, is that for strong
magnetizations, $\Omega_{i}\tau_{i}\gg 1$, the parallel $\eta_{0}$ coefficient
greatly exceeds the other viscosity coefficients. Hence, one can often
simplify by neglecting the smaller coefficients (although, as shown in § 3 it
can be necessary to retain other viscosity coefficients if length scales are
highly anisotropic). In this simplification, one has the following covariant
expressions for parallel viscosity (Lifshitz & Pitaevskii, 1981; Hollweg,
1986):

$\displaystyle\pi_{\alpha\beta}$
$\displaystyle=-3\eta_{0}\left(h_{\alpha}h_{\beta}-\frac{\delta_{\alpha\beta}}{3}\right)\left(h_{\mu}h_{\nu}-\frac{\delta_{\mu\nu}}{3}\right)\partial_{\mu}V_{\nu},$
(5) $\displaystyle Q_{visc}$
$\displaystyle=3\eta_{0}\left(\left(h_{\alpha}h_{\beta}-\frac{\delta_{\alpha\beta}}{3}\right)\partial_{\alpha}V_{\beta}\right)^{2},$

where $\mathbf{h}=\mathbf{B}/|\mathbf{B}|$ is the unit vector in the direction
of the magnetic field. These expressions are different to the isotropic
viscosity that appears in the Navier-Stokes equations, owing to the anisotropy
introduced by the magnetic field.

Parallel viscosity is closely related to pressure anisotropy. As pointed out
by Chew et al. (1956), when $\Omega_{i}\tau_{i}\gg 1$ the particle Lorentz
force makes the pressure tensor gyrotropic, giving it the form

$P_{\alpha\beta}=p_{\perp}\delta_{\alpha\beta}+(p_{||}-p_{\perp})h_{\alpha}h_{\beta}.$
(6)

This is a significant simplification, since the six degrees of freedom of a
general pressure tensor have been replaced with two variables, $p_{||}$ and
$p_{\perp}$. The definitions in Eqs. (4) then yield $p=(p_{||}+2p_{\perp})/3$
and

$\displaystyle\pi_{\alpha\beta}$
$\displaystyle=(p_{||}-p_{\perp})\left(h_{\alpha}h_{\beta}-\frac{\delta_{\alpha\beta}}{3}\right).$
(7)

Equation (7) shows that pressure anisotropy has an equivalent stress tensor,
which is proportional to $p_{||}-p_{\perp}$. Furthermore, Eqs. (5) and (7)
both have the form
$\pi_{\alpha\beta}\sim(h_{\alpha}h_{\beta}-\delta_{\alpha\beta}/3)$, so
equivalence of the stress tensors reduces to equivalence of the scalar factors
in the two equations. An illuminating analysis of the conditions under which
they converge has been written by Hollweg (1985, 1986), the most important
condition being that collisions (or other processes such wave-particle
interactions) relax the pressure anisotropy driven by velocity gradients to an
extent that the pressure is only weakly anisotropic. Classical MHD, for
comparison, assumes that pressure anisotropy can be neglected altogether.

For low collisionality, the quasistatic approximation in Braginskii MHD ceases
to be valid and strong pressure anisotropy may develop. Under these
conditions, separate evolution equations can be derived for $p_{||}$ and
$p_{\perp}$ (Chew et al., 1956; Hollweg, 1986). However, the closure problem
rears its head again, because those equations depend on the heat flow vector.
A simple approach to obtaining a closed system is to ignore the heat flow
vector, thus obtaining the CGL double adiabatic equations (Chew et al., 1956),
which are commonly used for collisionless plasma. More sophisticated
approaches also exist that solve for the evolution of the pressure anisotropy
or the evolution of the stress tensor, retaining the heat flow vector and
closing by other means. The works by Balescu (1988); Schunk & Nagy (2009);
Zank (2014); Hunana et al. (2019a, b, 2022) provide further reading on this
topic.

Summarising, there exists a family of adjacent (sometimes overlapping) single-
fluid models for plasmas. The most appropriate choice for a particular problem
and/or context depends on the collisionality. When MHD timescales are greater
than the ion collision time, Braginskii MHD provides rigorous closure and
treats the internal forces and heat flow more accurately than classical MHD.

##  3 Alfvén Wave Heating by Braginskii Viscosity

###  3.1 Model

We quantitatively examine the viscous dissipation for an Alfvén wave, which is
a transverse wave polarized so that the magnetic perturbation is perpendicular
to the equilibrium magnetic field and the wavevector. Setting the equilibrium
magnetic field in the $z$-direction, the magnetic perturbation in the
$x$-direction and the wavevector in the $yz$-plane, we consider a total
magnetic field of the form

$\mathbf{B}=b(y,z,t)\mathbf{e}_{x}+B_{0}\mathbf{e}_{z}.$ (8)

This ansatz automatically satisfies $\nabla\cdot\mathbf{B}=0$. For the
velocity field we assume the form

$\mathbf{V}=V_{x}(y,z,t)\mathbf{e}_{x}+V_{z}(y,z,t)\mathbf{e}_{z}.$ (9)

The $V_{x}$ is the dominant velocity component. In linearized theory it would
be the only component of $\mathbf{V}$. Additionally, we have explicitly
included a higher-order $V_{z}$ term that represents the nonlinear
ponderomotive flow parallel to the equilibrium magnetic field, which is driven
by gradients of the magnetic pressure perturbation $b^{2}/2\mu_{0}$ associated
with a finite-amplitude Alfvén wave (e.g. Hollweg, 1971). The $V_{z}$ term can
be dropped when the plasma is incompressible (see § 3.3), however it is
required for a nonlinear treatment of compressible plasma and affects the wave
heating via the parallel viscosity coefficient $\eta_{0}$ (as remarked in §
1). The expression for $V_{z}$ in classical MHD is given later in Eq. (29).

In a full solution, derivatives of $b^{2}/2\mu_{0}$ with respect to $y$ give
rise to an additional nonlinear $y$-component of $\mathbf{V}$, which in turn
produces a nonlinear $y$-component of $\mathbf{B}$. These terms are not shown
explicitly in Eqs. (8) and (9). Such terms were included by Nocera et al.
(1986) and appear not to affect our main conclusions, provided the
perpendicular wavelength of the Alfvén wave is sufficiently large.

The viscous force is determined from the viscous stress tensor
$\pi_{\alpha\beta}$ by

$F_{visc,\alpha}=-\frac{\partial\pi_{\alpha\beta}}{\partial x_{\beta}},$ (10)

and the viscous heating rate is determined using

$Q_{visc}=-\pi_{\alpha\beta}\frac{\partial V_{\alpha}}{\partial x_{\beta}},$
(11)

where $\alpha\in\left\\{x,y,z\right\\}$, $\beta\in\left\\{x,y,z\right\\}$, the
$x_{\beta}$ are components of the position vector, $V_{\alpha}$ are components
of $\mathbf{V}$ and repeated indices imply summation in the Einstein
convention.

A vital point is that the viscous stress tensor depends on the direction of
the magnetic field given by the unit vector
$\mathbf{h}=\mathbf{B}/|\mathbf{B}|$, which for our Alfvén wave model in Eq.
(8) has

$h_{x}=\frac{b}{\sqrt{B_{0}^{2}+b^{2}}},\quad h_{y}=0,\quad
h_{z}=\frac{B_{0}}{\sqrt{B_{0}^{2}+b^{2}}},$ (12)

with $h_{x}^{2}+h_{z}^{2}=1$. Our analysis differs from many past works by
considering ${h_{x}\neq 0}$ and identifying the dominant heating contribution
at the end, as opposed to setting $h_{x}=0$ before evaluating the damping
effect on Alfvén waves.

Applying formulas from §4 of Braginskii (1965) (equivalent matrix expressions
are given by Hogan, 1984), the stress tensor is related to five viscosity
coefficients by

$\pi_{\alpha\beta}=-\sum_{i=0}^{2}\eta_{i}W_{i\alpha\beta}+\sum_{i=3}^{4}\eta_{i}W_{i\alpha\beta}.$
(13)

The gyroviscous $\eta_{3}$ and $\eta_{4}$ terms do not contribute to heating,
so evaluating the heating rate $Q_{visc}$ requires

$\displaystyle W_{0\alpha\beta}$
$\displaystyle=\frac{3}{2}\left(h_{\alpha}h_{\beta}-\frac{1}{3}\delta_{\alpha\beta}\right)\left(h_{\mu}h_{\nu}-\frac{1}{3}\delta_{\mu\nu}\right)W_{\mu\nu},$
(14) $\displaystyle W_{1\alpha\beta}$
$\displaystyle=\left(\delta_{\alpha\mu}^{\perp}\delta_{\beta\nu}^{\perp}+\frac{1}{2}\delta_{\alpha\beta}^{\perp}h_{\mu}h_{\nu}\right)W_{\mu\nu},$
$\displaystyle W_{2\alpha\beta}$
$\displaystyle=\left(\delta_{\alpha\mu}^{\perp}h_{\beta}h_{\nu}+\delta_{\beta\nu}^{\perp}h_{\alpha}h_{\mu}\right)W_{\mu\nu},$

where $\delta_{\alpha\beta}$ is the Kronecker delta,

$\delta_{\alpha\beta}^{\perp}=\delta_{\alpha\beta}-h_{\alpha}h_{\beta},$ (15)

and the rate of strain tensor is

$W_{\alpha\beta}=\frac{\partial V_{\alpha}}{\partial x_{\beta}}+\frac{\partial
V_{\beta}}{\partial
x_{\alpha}}-\frac{2}{3}\delta_{\alpha\beta}\nabla\cdot\mathbf{V}.$ (16)

For the shear Alfvén wave geometry described by Eq. (9), the $\mathbf{W}_{i}$
tensors become

$\displaystyle\mathbf{W}_{0}$ $\displaystyle=$
$\displaystyle\left(h_{x}h_{z}\partial_{z}V_{x}+\left(\tfrac{2}{3}-h_{x}^{2}\right)\partial_{z}V_{z}\right)\begin{pmatrix}3h_{x}^{2}-1&0&3h_{x}h_{z}\\\
0&-1&0\\\ 3h_{x}h_{z}&0&3h_{z}^{2}-1\end{pmatrix},$ (17)
$\displaystyle\mathbf{W}_{1}$ $\displaystyle=$ $\displaystyle
h_{x}\left(h_{z}\partial_{z}V_{x}-h_{x}\partial_{z}V_{z}\right)\begin{pmatrix}-h_{z}^{2}&0&h_{x}h_{z}\\\
0&1&0\\\
h_{x}h_{z}&0&-h_{x}^{2}\end{pmatrix}+\left(h_{x}\partial_{y}V_{x}-h_{z}\partial_{y}V_{z}\right)\begin{pmatrix}0&h_{z}&0\\\
h_{z}&0&-h_{x}\\\ 0&-h_{x}&0\end{pmatrix},$ (18) $\displaystyle\mathbf{W}_{2}$
$\displaystyle=$
$\displaystyle\left(\left(1-2h_{x}^{2}\right)\partial_{z}V_{x}-2h_{x}h_{z}\partial_{z}V_{z}\right)\begin{pmatrix}2h_{x}h_{z}&0&1-2h_{x}^{2}\\\
0&0&0\\\
1-2h_{x}^{2}&0&-2h_{x}h_{z}\end{pmatrix}+\left(h_{x}\partial_{y}V_{x}+h_{z}\partial_{y}V_{z}\right)\begin{pmatrix}0&h_{x}&0\\\
h_{x}&0&h_{z}\\\ 0&h_{z}&0\end{pmatrix},$ (19)

The viscous heating rate with $h_{x}\neq 0$ retained is thus

$Q_{visc}=\frac{\eta_{0}}{3}\left(3h_{x}h_{z}\partial_{z}V_{x}+(2-3h_{x}^{2})\partial_{z}V_{z}\right)^{2}+\eta_{1}h_{x}^{2}\left(h_{z}\partial_{z}V_{x}-h_{x}\partial_{z}V_{z}\right)^{2}+\eta_{1}\left(h_{z}\partial_{y}V_{x}-h_{x}\partial_{y}V_{z}\right)^{2}\\\
+\eta_{2}\left(\left(1-2h_{x}^{2}\right)\partial_{z}V_{x}-2h_{x}h_{z}\partial_{z}V_{z}\right)^{2}+\eta_{2}\left(h_{x}\partial_{y}V_{x}+h_{z}\partial_{y}V_{z}\right)^{2}.$
(20)

###  3.2 Two small parameters

As anticipated in § 1 (e.g. Fig. 1), two parameters determine the relative
importance of individual terms in Eq. (20). The first small parameter is
$\left(b/B_{0}\right)^{2}$, the ratio of the wave’s magnetic energy density to
the energy density of the background magnetic field, which enters through
$h_{x}$ and $h_{z}$. In the modern era, extensive observations of coronal MHD
waves (Nakariakov & Verwichte, 2005; De Moortel & Nakariakov, 2012) allow
$\left(b/B_{0}\right)^{2}$ to be quantified with good certainty, directly from
resolved wave observations or indirectly from spectral line widths. For
example, Morton et al. (2015) studied waves at the base of a coronal open
field region using both approaches and reported a wave speed $v_{A}=400\mbox{
km s}^{-1}$ and wave motions at $v=35\mbox{ km s}^{-1}$. Both measurements are
consistent with earlier findings for coronal holes and the quiet Sun (e.g.
McIntosh et al., 2011). From observations like these,
$h_{x}^{2}\sim\left(b/B_{0}\right)^{2}\sim\left(v/v_{A}\right)^{2}\sim
10^{-2}$ .

The second small parameter is $\left(\Omega_{i}\tau_{i}\right)^{-2}$, which
sets the viscosity coefficients $\eta_{1}$ and $\eta_{2}$ relative to
$\eta_{0}$. The value of $\Omega_{i}\tau_{i}$ can vary significantly in the
corona, but if magnetic null points are excluded one obtains values similar to
the estimates made by Hollweg (1985), who found $3.4\times 10^{5}$ for a solar
active region and $7.2\times 10^{5}$ near the base of a coronal hole. We
therefore expect $(\Omega_{i}\tau_{i})^{-2}\lesssim 10^{-11}$ under common
conditions, and $\eta_{1}$ and $\eta_{2}$ simplify to

$\eta_{2}=\frac{6408}{5125}\left(\Omega_{i}\tau_{i}\right)^{-2}\eta_{0},\quad\eta_{1}=\frac{1}{4}\eta_{2}.$
(21)

The numerical coefficients in Eq. (21) are obtained in the limit
$(\Omega_{i}\tau_{i})^{-2}\to 0$, e.g. from Eq. (73) of Hunana et al. (2022).
They are approximate for finite $(\Omega_{i}\tau_{i})^{-2}$ but have a high
degree of accuracy because the corrections to the coefficients are of the
order of $(\Omega_{i}\tau_{i})^{-2}\lesssim 10^{-11}$. Inspecting Eq. (21) and
considering $(\Omega_{i}\tau_{i})^{-2}\lesssim 10^{-11}$, the $\eta_{2}$ and
$\eta_{1}$ coefficients are both vastly smaller than $\eta_{0}$.

The smallness of $\left(b/B_{0}\right)^{2}$ and the smallness of
$\left(\Omega_{i}\tau_{i}\right)^{-2}$ compete to make different terms
dominate the viscous heating. If one tries to simplify Eq. (20) by setting
$\left(b/B_{0}\right)^{2}$ to zero, then $h_{x}=0$ and $V_{z}=0$ gives
$Q_{visc}=\eta_{2}\left(\partial_{z}V_{x}\right)^{2}+\eta_{1}\left(\partial_{y}V_{x}\right)^{2}$,
as obtained by Braginskii (1965). On the other hand, if one tries to simplify
by first taking $\left(\Omega_{i}\tau_{i}\right)^{-2}$ to zero then only
$\eta_{0}$ terms remain, suggesting a different conclusion. Thus, the
quantitative results recover the two branches shown in Fig. 1. To correctly
determine the damping under coronal conditions, one must carefully compare
terms in the full Eq. (20), bearing in mind that there are two small
parameters, which we do now (diagonal branch in Fig. 1).

###  3.3 Heating rate for incompressible plasma

The analysis for incompressible plasma is relatively straightforward, which
makes it a natural starting point for discussion. The assumption of
incompressibility is appropriate for liquid metals or high-beta plasmas, but
not, we note, for the corona. The use of coronal wave amplitudes and
magnetizations in this section is therefore intended to be instructive only,
with the compressible finite-beta treatment that follows later in this paper
being required to treat the corona.

In the incompressible case, $\nabla\cdot\mathbf{V}=0$ applied to our Alfvén
wave geometry implies $V_{z}=0$. Thus Eq. (20) with
$\eta_{1}=\frac{1}{4}\eta_{2}$ simplifies to

$Q_{visc}=\left\\{3\eta_{0}h_{x}^{2}h_{z}^{2}+\eta_{2}\left(\left(1-2h_{x}^{2}\right)^{2}+\frac{1}{4}h_{x}^{2}h_{z}^{2}\right)\right\\}\left(\partial_{z}V_{x}\right)^{2}\\\
+\eta_{2}\left(\frac{1}{4}h_{z}^{2}+h_{x}^{2}\right)\left(\partial_{y}V_{x}\right)^{2}.$
(22)

The terms involving $\partial_{z}V_{x}$ set the viscous dissipation due to
wavelengths parallel to the equilibrium magnetic field, and we first ask
whether dissipation due to parallel wavelengths is dominated by the linear
$\eta_{2}$ contribution that has been widely recognised since Braginskii
(1965), or the nonlinear $\eta_{0}$ contribution. The ratio of the two terms
inside the curly brackets in Eq. (22) is

$3h_{x}^{2}\frac{\eta_{0}}{\eta_{2}}\left[\frac{1-h_{x}^{2}}{(1-2h_{x}^{2})^{2}+\frac{1}{4}h_{x}^{2}(1-h_{x}^{2})}\right]\approx
3h_{x}^{2}\frac{\eta_{0}}{\eta_{2}}\approx
2.4h_{x}^{2}(\Omega_{i}\tau_{i})^{2},$ (23)

where the first step simplifies using $h_{x}^{2}\ll 1$ (for the observed value
of $h_{x}^{2}\approx 10^{-2}$, retaining the terms in the square bracket
increases the ratio by 2.8%, so this approximation is both accurate and
conservative) and the substitution for $\eta_{0}/\eta_{2}$ is by Eq. (21). For
the coronal wave amplitudes and $\Omega_{i}\tau_{i}$ values noted in § 3.2,
this ratio exceeds $10^{9}$, with the nonlinear damping via $\eta_{0}$
dominating the heating rate by that factor. In other words, the viscous
dissipation of Alfvén waves via derivatives aligned with the equilibrium
magnetic field is a factor $10^{9}$ stronger than predicted by linear theory.

We now evaluate the role of derivatives perpendicular to the equilibrium
magnetic field by comparing the nonlinear $\eta_{0}$ term in Eq. (22) to the
term involving $\partial_{y}V_{x}$. The ratio of these heating rate terms is

$12h_{x}^{2}\frac{\eta_{0}}{\eta_{2}}\left(\frac{\lambda_{\perp}}{\lambda_{||}}\right)^{2}\left[\frac{1-h_{x}^{2}}{1+3h_{x}^{2}}\right]\approx
12h_{x}^{2}\frac{\eta_{0}}{\eta_{2}}\left(\frac{\lambda_{\perp}}{\lambda_{||}}\right)^{2}\approx
9.6h_{x}^{2}\left(\Omega_{i}\tau_{i}\right)^{2}\left(\frac{\lambda_{\perp}}{\lambda_{||}}\right)^{2},$
(24)

where $\lambda_{\perp}$ and $\lambda_{||}$ are the wavelengths perpendicular
and parallel to the equilibrium magnetic field (for the observed value of
$h_{x}^{2}\approx 10^{-2}$, the approximation of the terms in $h_{x}^{2}$ is
accurate to 3.9%). For the coronal parameters noted above, if
$\lambda_{\perp}\approx\lambda_{||}$ then the nonlinear $\eta_{0}$ term again
dominates by a factor that exceeds $10^{9}$. For smaller transverse
wavelengths, the nonlinear $\eta_{0}$ term dominates whenever
$\lambda_{\perp}\gtrsim 10^{-5}\lambda_{||}$. If one considers a wave speed of
$400\mbox{ km s}^{-1}$ and a frequency of 3 mHz, consistent with the
observations by Morton et al. (2015), the condition that the nonlinear
$\eta_{0}$ term dominates becomes $\lambda_{\perp}\gtrsim 800\mbox{ m}$. Given
that CoMP has imaged Alfvénic waves using 3 Mm pixels, this condition appears
to be met by a very large margin, making the nonlinear $\eta_{0}$ dissipation
dominant over the $\eta_{2}$ linear dissipation.

The purpose of deriving Eq. (22) and the ratios on the left hand sides of Eq.
(23) and (24) such that they include all appearances of $h_{x}$ and $h_{z}$ is
that they can be evaluated exactly for a given value of $h_{x}$. This makes it
explicit that our conclusions are insensitive to the precise value of $h_{x}$,
only that the value of $h_{x}$ is broadly consistent with coronal
observations. While that approach is most comprehensive, the same conclusions
can also be reached by separately simplifying each term in Eq. (22) using
$h_{x}^{2}\ll 1$ and $h_{z}^{2}=1-h_{x}^{2}\approx 1$ to obtain the less
cumbersome formula

$Q_{visc}=\left\\{3\eta_{0}h_{x}^{2}+\eta_{2}\right\\}\left(\partial_{z}V_{x}\right)^{2}+\eta_{1}\left(\partial_{y}V_{x}\right)^{2},$
(25)

and comparing terms to reach the same conclusions.

###  3.4 Compressible plasma with large Re

Under typical coronal conditions, the thermal pressure is too small to prevent
compression of the plasma by nonlinear magnetic pressure forces, thus a
nonlinear $V_{z}$ develops that is known as the ponderomotive flow (Hollweg,
1971). This flow component affects the viscous heating rate via the parallel
viscosity coefficient $\eta_{0}$, hence compressible theory is required for
nonlinear viscous damping of Alfvén waves in plasma.

We define the Alfvénic Reynolds number as

$\mathrm{Re}=\frac{\rho v_{A}}{k_{||}\eta_{0}}.$ (26)

This dimensionless parameter differs from the traditional Reynolds number
since it refers to the Alfvén speed $v_{A}=B/\sqrt{\mu_{0}\rho}$ instead of a
typical fluid velocity. This distinction mirrors that between the Lundquist
number and magnetic Reynolds number in resistive MHD. Justification for
defining $\mathrm{Re}$ according to Eq. (26) will be found in the detailed
mathematical solutions in § 5, in which it is found to be a natural parameter
of the system (also see Nocera et al., 1986).

In this section, the ponderomotive $V_{z}$ will be related to $V_{x}$ using
expansions in the amplitude of the primary wave fields. Several assumptions
are used to accomplish this. First, we make use of the result that a
travelling wave solution propagating in the positive $z$-direction has

$\frac{\partial}{\partial t}\equiv-v_{A}\frac{\partial}{\partial z},$ (27)

where $v_{A}$ is the wave speed. For simplicity it is assumed that derivatives
of background quantities are sufficiently weak to play a higher-order role on
the dynamics. We also simplify here by replacing full treatment of thermal
conduction with two thermodynamic cases: adiabatic and isothermal. Finally, it
is assumed that $\mathrm{Re}$ is large enough that viscous forces can be
neglected at leading order when evaluating $V_{z}$, which makes it possible to
obtain an algebraic relationship between $V_{z}$ and $b^{2}$. This assumption
will be removed for § 5, in which the effect of viscous forces on $V_{x}$ and
$V_{z}$ is included.

The $x$-components of the momentum and induction equations are unaffected by
the ponderomotive flow at leading order in the wave amplitude. From them one
recovers the Walén relation for propagating Alfvén waves,
$b/B_{0}=-V_{x}/v_{A}$.

At leading order, the $z$-component of the momentum equation is

$\rho_{0}\frac{\partial V_{z}}{\partial t}+\frac{\partial}{\partial
z}\left(\delta p+\frac{b^{2}}{2\mu_{0}}\right)=0,$ (28)

where the viscous force has been neglected since we currently consider the
limit of large $\mathrm{Re}$. Using Eq (27) and integrating yields an
algebraic relationship between $V_{z}$, $\delta p$ and $b^{2}$. In an
adiabatic treatment, the energy equation yields $\delta p=\gamma
p_{0}V_{z}/v_{A}$, hence we obtain

$\frac{V_{z}}{v_{A}}=\frac{1}{2(1-\beta)}\left(\frac{b}{B_{0}}\right)^{2},$
(29)

where

$\beta=\left(\frac{c_{s}}{v_{A}}\right)^{2}=\frac{\gamma}{2}\left(\frac{p}{B_{0}^{2}/2\mu_{0}}\right).$
(30)

In an isothermal treatment, the ideal gas law $p=\rho RT$ yields $\delta
p/p_{0}=\delta\rho/\rho_{0}=V_{z}/v_{A}$. This does not change the form of Eq.
(29); instead, the isothermal case is recovered simply by setting $\gamma=1$
in the definition of $\beta$.

Defining $\beta$ as the square of the ratio of the sound speed
($c_{s}=\sqrt{\gamma p_{0}/\rho_{0}}$) to the Alfvén speed differs slightly
from the convention of defining $\beta$ as the ratio of thermal pressure to
magnetic pressure, due to the factor $\gamma/2$, which is 5/6 for an adiabatic
monoatomic gas, and $1/2$ for an isothermal model. Defining $\beta$ as the
speed ratio squared leads to cleaner mathematics for many MHD wave problems,
including this one, and it has therefore become established practice in MHD
wave theory.

The $\beta=1$ singularity in Eq. (29) arises because $c_{s}=v_{A}$ implies
resonance between the Alfvén wave and an acoustic wave, which resonantly
transfers energy between the waves. In this specific case, Eq. (27) does not
apply because it does not account for evolution due to resonance. Similarly,
Eq. (28) assumes that $V_{z}\ll v_{A}$ to simplify the convective derivative,
and the solution in Eq. (29) does not satisfy this condition in the immediate
vicinity of $\beta=1$. The $\beta=1$ resonance and the singularity in Eq. (29)
are not of concern for most coronal applications, which typically have
$\beta<0.2$, but there are special cases in which it is of interest, such as
waves propagating towards coronal magnetic nulls or across the $\beta=1$ layer
in the lower solar atmosphere. Russell et al. (2016) have previously applied
such nonlinear resonant coupling to the problem of sunquake generation by
magnetic field changes during solar flares.

Differentiating Eq. (29) and employing the relation $b/B_{0}=-V_{x}/v_{A}$
yields

$\partial_{\alpha}V_{z}=-\frac{1}{(1-\beta)}\frac{h_{x}}{h_{z}}\partial_{\alpha}V_{x},$
(31)

which can be used to eliminate $V_{z}$ from Eq. (20). To leading order in
$h_{x}^{2}$ in each viscosity coefficient, we find

$Q_{visc}=\left\\{C\eta_{0}h_{x}^{2}+\eta_{2}\right\\}\left(\partial_{z}V_{x}\right)^{2}+\eta_{1}(\partial_{y}V_{x})^{2},$
(32)

where

$C=\frac{1}{3}\left(\frac{1-3\beta}{1-\beta}\right)^{2}.$ (33)

Some special cases are noteworthy. The incompressible results of § 3.3 are
recovered for $\beta\to\infty$, which gives $V_{z}\to 0$ and $C\to 3$.
Similarly, the cold plasma solution is recovered by setting $\beta=0$, which
gives $V_{z}=\left(b/B_{0}\right)^{2}/2$ and $C=1/3$.

The nonlinear heating rate for cold plasma ($\beta=0$) is a factor nine
smaller than for incompressible plasma ($\beta\to\infty$), which demonstrates
the importance of compressibility for this problem. Furthermore, $C(\beta)$ is
monotonically decreasing between $\beta=0$ and $\beta=1/3$. Since
$0<\beta<1/3$ for most coronal applications, the dissipation rate due to
nonlinear Braginskii viscosity in these environments is reduced compared to
the cold plasma solution. For example, given $\beta=0.1$, the heating rate is
approximately 60% of the value for cold plasma. It is therefore evident that
compressibility and finite-beta effects must be treated when assessing viscous
dissipation of Alfvén waves.

Another important feature is that $C$ has a zero for $\beta=1/3$. This is one
circumstance in which

$\frac{V_{z}}{v_{A}}=\frac{3}{4}\left(\frac{V_{x}}{v_{A}}\right)^{2}=\frac{3}{4}\left(\frac{b}{B_{0}}\right)^{2},$
(34)

which causes cancellation within the $\eta_{0}$ contribution to $Q_{visc}$.
That a particular organisation of $V_{z}/v_{A}$ can suppress nonlinear viscous
dissipation is an important novel finding that § 4 explores further in the
context of low $\mathrm{Re}$.

The final feature of $C$ is the singularity at $\beta=1$. As noted earlier in
this section, $c_{s}=v_{A}$ implies that the Alfvén wave is in resonance with
a sound wave, which transfers energy between the Alfvén wave and the sound
wave. Caution is needed around the resonance, since resonant energy transfer
cannot be described using Eq. (27), which was used to derive Eq. (33).

Evaluation of ratios of heating terms from Eq. (32) proceeds as for the
comparison in Sec. 3.3, but with $\eta_{0}$ multiplied by $C/3$. The top-level
conclusions remain intact: heating by the Braginskii viscous stress tensor is
dominated by an $\eta_{0}$ term that is nonlinear in the wave amplitude, and
for coronal values, heating due to the nonlinear $\eta_{0}$ term is many
orders of magnitude larger than the heating due to the linear $\eta_{1}$ and
$\eta_{2}$ terms.

###  3.5 What wave amplitude is linear?

An important implication of the preceding analysis is that nonlinear effects
become significant for anisotropic viscosity at far lower wave amplitudes than
they do for other terms in the MHD equations. Linearizing the Braginskii
viscous stress tensor is only appropriate when
$h_{x}^{2}\sim(b/B_{0})^{2}\ll\left(\Omega_{i}\tau_{i}\right)^{-2}$, which in
the corona corresponds to a requirement that the wave energy density is less
than $10^{-11}$ times the energy density of the background magnetic field, far
too small to be relevant to coronal energetics. Waves that have small enough
amplitudes to be governed by linear viscous damping theory would be
unobservable and have no effect on the coronal energy balance. Thus, for
coronal Alfvén waves, viscosity must be treated nonlinearly in the wave
amplitude, as well as anisotropically due to the magnetic field.
Interestingly, this linearization condition is far more stringent than the
linearization condition for other terms in the MHD equations, whereby
$h_{x}^{2}$ is normally compared to unity. The extreme difference in these
linearization conditions is due to the large $\eta_{0}/\eta_{2}$ ratio
produced by the strong magnetization.

###  3.6 Damping scales for large Re (energy derivation)

It is of major interest to know the time and length scales over which waves
damp. This section provides a relatively simple derivation of the decay scales
for nonlinear viscous damping of propagating shear Alfvén waves for large
$\mathrm{Re}$, using energy principles.

Dropping the $\eta_{1}$ and $\eta_{2}$ terms from Eq. (32), the heating rate
due to the $\eta_{0}$ parallel viscosity coefficient for large $\mathrm{Re}$
is

$Q_{visc}=\frac{\eta_{0}}{3}\left(\frac{1-3\beta}{1-\beta}\right)^{2}\left(\frac{b}{B_{0}}\right)^{2}(\partial_{z}V_{x})^{2}.$
(35)

A wave energy decay time can be defined according to
$\tau_{d}=\left<E_{w}\right>/\left<Q_{visc}\right>$, where $\left<.\right>$
denotes the time average over a wave period, and $E_{w}$ is the wave energy
density. The corresponding decay length is $L_{d}=v_{A}\tau_{d}$.

For forward propagating Alfvén waves, $E_{w}\approx\rho_{0}V_{x}^{2}$, and
$V_{x}=a\cos(\phi)$ where $\phi=k_{||}(z-v_{A}t)$. Hence,

$E_{w}=\rho_{0}a^{2}\frac{(1+\cos(2\phi))}{2}.$ (36)

Similarly, using Eq. (35) with $(b/B_{0})^{2}\approx(V_{x}/V_{A})^{2}$,

$Q_{visc}=\frac{\eta_{0}k_{||}^{2}}{3v_{A}^{2}}\left(\frac{1-3\beta}{1-\beta}\right)^{2}a^{4}\frac{(1-\cos(4\phi))}{8},$
(37)

which give the fast time averages

$\displaystyle\left<E_{w}\right>$ $\displaystyle=$
$\displaystyle\frac{\rho_{0}a^{2}}{2},$ (38) $\displaystyle\left<Q\right>$
$\displaystyle=$
$\displaystyle\frac{\eta_{0}k_{||}^{2}a^{4}}{24v_{A}^{2}}\left(\frac{1-3\beta}{1-\beta}\right)^{2}.$
(39)

The wave energy decay scales are therefore

$\displaystyle\tau_{d}=\frac{12\rho_{0}}{\eta_{0}k_{||}^{2}(a/v_{A})^{2}}\left(\frac{1-\beta}{1-3\beta}\right)^{2},$
(40) $\displaystyle
L_{d}=\frac{12\rho_{0}v_{A}}{\eta_{0}k_{||}^{2}(a/v_{A})^{2}}\left(\frac{1-\beta}{1-3\beta}\right)^{2}.$
(41)

Equations (40) and (41) show that waves with larger $k_{||}$ (equivalently,
higher frequencies) are damped on shorter scales. We also remark that since
$L_{d}$ depends on the amplitude of $V_{x}$ (the constant $a$), the decay
envelope is non-exponential. The damping properties are elaborated on more
fully in § 5, in which the assumption of large $\mathrm{Re}$ is removed and
the functional form of the wave envelope is determined.

##  4 Self-organised viscous flow

In the limit $\mathrm{Re}\to 0$, the viscous force in the $z$-component of the
momentum equation risk becoming extremely large, unless the flow self-
organises to prevent this. Correspondingly, in the limit $\mathrm{Re}\to 0$,
strong dissipation will prevent waves from propagating, unless $V_{z}$ is
determined by viscosity. One can therefore expect self-organisation of the
flow pattern for Alfvén waves in highly viscous plasma (small values of
$\mathrm{Re}$), which is a concept previously advanced by Montgomery (1992).

To investigate quantitatively, we analyze the highly magnetised regime
$\Omega_{i}\tau_{i}\gg 1$, simplifying the stress tensor and heating rate by
retaining only the $\eta_{0}$ parallel viscosity coefficient. Inspecting Eq.
(5), components of $\pi_{\alpha\beta}$ are proportional to
$(h_{\mu}h_{\nu}-\delta_{\mu\nu}/3)\partial_{\mu}V_{\nu}$, and $Q_{visc}$ is
proportional to the square of this expression. Applying the shear Alfvén wave
geometry of Eqs. (8) and (9) and simplifying by $h_{x}^{2}\ll 1$,

$\left(h_{\mu}h_{\nu}-\frac{\delta_{\mu\nu}}{3}\right)\frac{\partial
V_{\nu}}{\partial x_{\mu}}\approx h_{x}\frac{\partial V_{x}}{\partial
z}+\frac{2}{3}\frac{\partial V_{z}}{\partial z}.$ (42)

Viscous forces and heating can be suppressed, allowing Alfvén wave
propagation, if the flow self-organises to keep this expression close to zero.
Using the Alfvén wave relation to substitute $h_{x}\approx b/B_{0}\approx-
V_{x}/v_{A}$ and integrating, we find that for small $\mathrm{Re}$

$\frac{V_{z}}{v_{A}}=\frac{3}{4}\left(\frac{V_{x}}{v_{A}}\right)^{2}=\frac{3}{4}\left(\frac{b}{B_{0}}\right)^{2}.$
(43)

The relation specified by Eq. (43) appeared previously in the different
context of § 3.4, where it was seen that viscous dissipation of Alfvén waves
in high $\mathrm{Re}$ plasma is suppressed for the special case of
$\beta=1/3$. The flow pattern required to produce cancellation within the
$\eta_{0}$ part of $Q_{visc}$ is independent of $\mathrm{Re}$ and $\beta$, but
it occurs for different reasons in the two cases: in § 3.4 it arose as a
special case of ponderomotive flow with finite $\beta$; when $\mathrm{Re}$ is
small, it occurs because of self-organisation through viscous forces.

This novel result demonstrates that decay scales and other properties derived
in § 3 should not be extrapolated to small $\mathrm{Re}$. Instead, we expect
that as $\mathrm{Re}\to 0$, the viscous force organises the flow such that
$V_{z}$ obeys Eq. (43), for which dissipation is suppressed by cancellation
within the $\eta_{0}$ part of $Q_{visc}$.

##  5 Multiple scale analysis

Section 3 used methods of analysis based on heating rates and energy
principles. Section 5 now takes a complementary approach of solving the full
set of governing equations using multiple scale analysis, to reinforce the
results of § 3, extend to general $\mathrm{Re}$ by including the effect of the
viscous force on $\mathbf{V}$, and obtain additional results including the
functional form of the nonlinear decay.

###  5.1 Comparison to Nocera et al. (1986)

We preface the multiple scale analysis part of this paper with some remarks
about related calculations by Nocera et al. (1986). Their work and ours both
concentrate on $\eta_{0}$ viscosity as the main source of wave damping,
treating this nonlinearly in the wave amplitude (the horizontal branch of
Figure 1). Also in common, both treat ponderomotive and finite $\beta$
effects.

The previous work of Nocera et al. (1986) derived a version of the viscous
stress tensor that includes the leading order effect of $h_{x}\neq 0$ in the
$\eta_{0}$ term. Terms in the viscosity tensor were then compared, concluding
like our § 3 (but by different arguments) that the nonlinear $\eta_{0}$ term
exceeds contributions from other viscosity coefficients when
$(b/B_{0})^{2}\gg(\Omega_{i}\tau_{i})^{-2}$ (their Eq. (3.13)). The two
studies thus agree on the dominance of nonlinear $\eta_{0}$ viscosity.

Nocera et al. (1986) then found a decay length using the following strategy. A
self-consistent perturbation ordering was introduced, then the $x$-components
of the momentum and induction equations were combined to obtain a single
equation for $V_{x}$, which at linear order is a wave equation. Next, all
variables apart from $V_{x}$ were eliminated from the leading-order nonlinear
term. Finally, they concluded from a stability analysis that waves with
$k_{\perp}=0$ are damped nonlinearly, with a decay time that has the same form
as our Eq. (40) (their Eq. (5.7), given in terms of normalised variables).

The detailed derivation that follows in § 5.2 draws inspiration from the
framework developed by Nocera et al. (1986). We have also taken the
opportunity to make several changes that we regard as improvements, most
importantly:

  1. 1.

Nocera et al. (1986) assumed that the fast time average of $V_{z}$ is zero,
which necessitated adding a non-zero constant of integration to $V_{z}$. By
contrast, we will set the constant of integration to zero, which is the only
choice for which an Alfvén wave driver switching on at one boundary does not
unphysically send an instantaneous signal to infinity. Additional support for
our choice comes from simulations of nonlinear longitudinal flows produced by
Alfvén waves (e.g. McLaughlin et al., 2011), which are consistent with the
constraint used in our work.

  2. 2.

The stability analysis in §5 Nocera et al. (1986) is replaced with a multiple
scale analysis of the type covered in Chapter 11 of Bender & Orszag (1978).

  3. 3.

Nocera et al. (1986) made their wave envelope a function of $z+v_{A}t$. We
treat the envelope as time-independent and thus explicitly investigate damping
of a propagating wave with respect to distance.

  4. 4.

Our derivation provides the envelope of $V_{x}$ as well as the decay length.

  5. 5.

Our solution is valid for general $\mathrm{Re}$, whereas Nocera et al. (1986)
solved for the decay scales in the low-viscosity limit of high $\mathrm{Re}$
only.

Equally, Nocera et al. (1986) treated cases that we do not, including the
possibility of $k_{\perp}$ large enough for coupling between the Alfvén and
fast modes to alter the wave properties (referred to in their paper as the
case of phase mixed waves).

###  5.2 Detailed solution

####  5.2.1 Geometry and perturbations

We assume the Alfvén wave geometry of Eqs. (8) and (9), set $\partial/\partial
y\equiv 0$ to concentrate on waves without short perpendicular scales, and
introduce density and pressure perturbations $\delta\rho$ and $\delta p$
together with a self-consistent perturbation ordering that has
$V_{x}/v_{A}\sim b/B_{0}\sim\epsilon^{1/2}$ and
$V_{z}/v_{A}\sim\delta\rho/\rho_{0}\sim\delta p/p_{0}\sim\epsilon$. The
viscosity $\eta_{0}$ and background quantities $B_{0}$, $\rho_{0}$ and $p_{0}$
are treated as locally homogeneous for simplicity.

####  5.2.2 Nonlinear wave equation

Starting from the ideal induction equation,

$\frac{\partial\mathbf{B}}{\partial
t}=\nabla\times\left(\mathbf{V}\times\mathbf{B}\right),$ (44)

we have

$\frac{\partial b}{\partial t}-B_{0}\frac{\partial V_{x}}{\partial
z}=-\frac{\partial}{\partial z}\left(bV_{z}\right)\quad\text{(exact)},$ (45)

where the linear terms have been grouped on the left hand side and the
nonlinear term on the right hand side.

The momentum equation is

$\rho\left(\frac{\partial V_{\alpha}}{\partial
t}+(\mathbf{V}\cdot\nabla)V_{\alpha}\right)=-\frac{\partial}{\partial
x_{\alpha}}\left(p+\frac{B^{2}}{2\mu_{0}}\right)+\frac{1}{\mu_{0}}(\mathbf{B}\cdot\nabla)B_{\alpha}-\frac{\partial\pi_{\alpha\beta}}{\partial
x_{\beta}}.$ (46)

When $\eta_{0}$ contributions dominate the viscous force, Eqs. (13) and (17)
give

$-\pi_{xz}=3\eta_{0}\frac{b}{B_{0}}\left(\frac{b}{B_{0}}\frac{\partial
V_{x}}{\partial z}+\frac{2}{3}\frac{\partial V_{z}}{\partial
z}\right)+O(\epsilon^{5/2})$ (47)

so the $x$-component of Eq. (46) becomes

$\frac{\partial V_{x}}{\partial t}-\frac{B_{0}}{\mu_{0}\rho_{0}}\frac{\partial
b}{\partial z}=-\frac{\delta\rho}{\rho_{0}}\frac{\partial V_{x}}{\partial
t}-V_{z}\frac{\partial V_{x}}{\partial
z}+\frac{3\eta_{0}}{\rho_{0}}\frac{\partial}{\partial
z}\left(\frac{b}{B_{0}}\left[\frac{b}{B_{0}}\frac{\partial V_{x}}{\partial
z}+\frac{2}{3}\frac{\partial V_{z}}{\partial
z}\right]\right)+O(\epsilon^{5/2}),$ (48)

where linear terms and nonlinear terms have again been placed on opposite
sides of the equation.

Taking the time derivative of Eq. (48) and using Eq. (45) to eliminate $b$
from the linear terms,

$\left(\frac{\partial^{2}}{\partial
t^{2}}-v_{A}^{2}\frac{\partial^{2}}{\partial
z^{2}}\right)V_{x}=-v_{A}^{2}\frac{\partial^{2}}{\partial
z^{2}}\left(\frac{b}{B_{0}}V_{z}\right)-\frac{\partial}{\partial
t}\left(\frac{\delta\rho}{\rho_{0}}\frac{\partial V_{x}}{\partial
t}+V_{z}\frac{\partial V_{x}}{\partial
z}\right)+\frac{3\eta_{0}}{\rho_{0}}\frac{\partial^{2}}{\partial t\partial
z}\left(\frac{b}{B_{0}}\left[\frac{b}{B_{0}}\frac{\partial V_{x}}{\partial
z}+\frac{2}{3}\frac{\partial V_{z}}{\partial
z}\right]\right)+O(\epsilon^{5/2}).$ (49)

Interpreting Eq. (49), the linear terms (on the left hand side) correspond to
a wave equation with wave speed $v_{A}$. The leading nonlinear terms (those
shown explicitly on the right hand side) include the leading-order effect of
the anisotropic viscosity, which enters at the same order as the leading
nonlinear terms that appear in perturbative nonlinear theory of ideal Alfvén
waves.

Next, we eliminate $b$ and $\delta\rho$ from the $O(\epsilon^{3/2})$ nonlinear
terms in Eq. (49). Equations (45) and (48) are solved at linear order by the
Alfvén wave relation

$\frac{b}{B_{0}}=\pm\frac{V_{x}}{v_{A}}+O(\epsilon^{3/2}).$ (50)

We choose the negative sign so waves travel in the positive $z$ direction,
giving

$\frac{b}{B_{0}}=-\frac{V_{x}}{v_{A}}+O(\epsilon^{3/2}).$ (51)

The travelling wave behaviour of the linear solution together with assumption
that the wave envelope changes over a distance controlled by the leading order
nonlinear terms in Eq. (49) allows replacement

$\frac{\partial}{\partial t}=-v_{A}\frac{\partial}{\partial z}+O(\epsilon).$
(52)

The density perturbation is governed by the mass continuity equation

$\frac{\partial\rho}{\partial t}+\nabla\cdot\left(\rho\mathbf{V}\right)=0,$
(53)

which gives for our shear Alfvén wave

$\frac{\partial\delta\rho}{\partial t}=-\rho_{0}\frac{\partial V_{z}}{\partial
z}+O(\epsilon^{2}).$ (54)

Then, using Eq. (52) and integrating,

$\frac{\delta\rho}{\rho_{0}}=\frac{V_{z}}{v_{A}}+O(\epsilon^{2}).$ (55)

The constant of integration has been set to zero, for reasons discussed in §
5.1.

Using these results, Eq. (49) becomes

$\left(\frac{\partial^{2}}{\partial
t^{2}}-v_{A}^{2}\frac{\partial^{2}}{\partial
z^{2}}\right)V_{x}=v_{A}\frac{\partial^{2}}{\partial
z^{2}}\left(V_{x}V_{z}\right)+\frac{3\eta_{0}}{\rho_{0}}\frac{\partial^{2}}{\partial
t\partial z}\left(\frac{V_{x}}{v_{A}}\left[\frac{V_{x}}{v_{A}}\frac{\partial
V_{x}}{\partial z}-\frac{2}{3}\frac{\partial V_{z}}{\partial
z}\right]\right)+O(\epsilon^{5/2}).$ (56)

Now that the problem has been reduced to the two variables $V_{x}$ and
$V_{z}$, it is convenient to make the $\epsilon$ dependence explicit by
introducing dimensionless variables $v$ and $w$ defined by

$V_{x}(z,t)=\epsilon^{1/2}v_{A}v(z,t),\quad V_{z}(z,t)=\epsilon v_{A}w(z,t).$
(57)

Expressing Eq. (56) in the dimensionless variables $v$ and $w$, and dropping
the non-explicit higher order terms from the right hand side, we seek
solutions to

$\left(\frac{\partial^{2}}{\partial
t^{2}}-v_{A}^{2}\frac{\partial^{2}}{\partial
z^{2}}\right)v=\epsilon\left(v_{A}^{2}\frac{\partial^{2}}{\partial
z^{2}}(vw)+\frac{\eta_{0}}{\rho_{0}}\frac{\partial^{2}}{\partial t\partial
z}\left(\frac{\partial v^{3}}{\partial z}-2v\frac{\partial w}{\partial
z}\right)\right).$ (58)

####  5.2.3 Multiple scale analysis

Equation (58) is now solved using multiple scale analysis (e.g. Bender &
Orszag, 1978). Applying this technique, one introduces a new variable
$Z=\epsilon z$ that defines a long length scale, and the perturbation
expansions

$\displaystyle v(z,t)$ $\displaystyle=$ $\displaystyle v_{0}(z,Z,t)+\epsilon
v_{1}(z,Z,t)+\ldots$ (59) $\displaystyle w(z,t)$ $\displaystyle=$
$\displaystyle w_{0}(z,Z,t)+\epsilon w_{1}(z,Z,t)+\ldots.$ (60)

Derivatives are treated using the chain rule as though $z$ and $Z$ are
independent variables and setting $\mathrm{d}Z/\mathrm{d}z=\epsilon$. Thus,

$\displaystyle\frac{\partial v}{\partial z}$ $\displaystyle=$
$\displaystyle\frac{\partial v_{0}}{\partial z}+\epsilon\left(\frac{\partial
v_{0}}{\partial Z}+\frac{\partial v_{1}}{\partial z}\right)+O(\epsilon^{2}),$
(61) $\displaystyle\frac{\partial^{2}v}{\partial z^{2}}$ $\displaystyle=$
$\displaystyle\frac{\partial^{2}v_{0}}{\partial
z^{2}}+\epsilon\left(2\frac{\partial^{2}v_{0}}{\partial Z\partial
z}+\frac{\partial^{2}v_{1}}{\partial z^{2}}\right)+O(\epsilon^{2}),$ (62)

with equivalent expressions for derivatives of $w$.

Substituting into Eq. (58), collecting $\epsilon^{0}$ terms, and thus solving
the homogeneous wave equation

$\left(\frac{\partial^{2}}{\partial
t^{2}}-v_{A}^{2}\frac{\partial^{2}}{\partial z^{2}}\right)v_{0}=0,$ (63)

obtains d’Alembert’s solution

$v_{0}(z,Z,t)=f(z-v_{A}t,Z)+g(z+v_{A}t,Z).$ (64)

For forward propagating waves, the function $g$ is zero.

The corresponding $w_{0}$ is obtained by integrating the $z$-component of the
momentum equation, Eq. (46), which gives

$\frac{V_{z}}{v_{A}}=\frac{1}{\rho_{0}v_{A}^{2}}\left(\delta
p+\frac{b^{2}}{2\mu_{0}}+\pi_{zz}\right)+O(\epsilon^{2}).$ (65)

From Eqs. (13) and (17), we have

$-\pi_{zz}=2\eta_{0}\left(\frac{b}{B_{0}}\frac{\partial V_{x}}{\partial
z}+\frac{2}{3}\frac{\partial V_{z}}{\partial z}\right)+O(\epsilon^{2}).$ (66)

A substitution for the pressure perturbation $\delta p$ is obtained by
integrating the energy equation,

$\frac{\partial p}{\partial t}+\mathbf{V}\cdot\nabla p+\gamma
p\nabla\cdot\mathbf{V}=(\gamma-1)Q_{visc}.$ (67)

The viscous heating $Q_{visc}$ is of order $O(\epsilon^{2})$, so integration
gives the adiabatic relation

$\frac{\delta p}{p_{0}}=\gamma\frac{V_{z}}{v_{A}}+O(\epsilon^{2}).$ (68)

Alternatively, one can consider isothermal conditions using $\delta
p/p_{0}=V_{z}/v_{A}$ from the ideal gas law, which is recovered from Eq. (68)
by setting $\gamma=1$.

Using Eqs. (66) and (68), and eliminating $b$ terms using Eq. (51), Eq. (65)
can be expressed as

$\left(1-\beta+\frac{4}{3}\frac{\eta_{0}}{\rho_{0}v_{A}}\frac{\partial}{\partial
z}\right)w=\left(\frac{1}{2}+\frac{\eta_{0}}{\rho_{0}v_{A}}\frac{\partial}{\partial
z}\right)v^{2},$ (69)

where $\beta=(c_{s}/v_{A})^{2}$ and dropped terms are $O(\epsilon)$. When $v$
and $w$ are expanded according to Eqs. (59) and (60), an equation identical to
(69) connects $w_{0}$ and $v_{0}$.

Inspecting Eq. (69), it is evident that obtaining $w$ for a known $v$ in
general requires solving a first order linear partial differential equation.
In the limit where the viscous terms can be neglected, the problem simplifies
to the algebraic $w=v^{2}/2(1-\beta)$ relation used in Sec. 3.4. Similarly,
when the viscous terms dominate, one obtains the $w=(3/4)v^{2}$ relation for
viscously self-organised parallel flow discussed in Sec. 3.4. For the detailed
solution in this section, we retain the complete set of forces that determine
$V_{z}$, solving the full Eq. (69).

Solution is facilitated by considering the special case where $v_{0}$
oscillates sinusoidally in time. For the rest of this derivation we therefore
set

$\displaystyle v_{0}$
$\displaystyle=A(Z)\mathrm{e}^{i\phi}+A^{*}(Z)\mathrm{e}^{-i\phi},$
$\displaystyle\phi$ $\displaystyle=k_{||}(z-v_{A}t),$ (70)

where $A(Z)\in\mathbb{C}$ and $*$ denotes the complex conjugate. Representing
$A$ in polar form,

$A(Z)=R(Z)\mathrm{e}^{i\theta(Z)},$ (71)

Eq. (70) is equivalent to

$v_{0}(z,Z,t)=2R(Z)\cos(k_{||}(z-v_{A}t)+\theta(Z)).$ (72)

From inspection, $2R(Z)$ is the local amplitude and $\theta(Z)$ is a phase
shift. We have found the complex form in Eq. (70) more convenient to work with
in the following.

We now solve for the corresponding $w_{0}$. Noting that

$v_{0}^{2}=A^{2}\mathrm{e}^{2i\phi}+2|A|^{2}+{A^{*}}^{2}\mathrm{e}^{-2i\phi},$
(73)

where $|A|^{2}=AA^{*}$, we seek a solution of the form

$w_{0}=D\mathrm{e}^{2i\phi}+D^{*}\mathrm{e}^{-2i\phi}+\hat{w}_{0}.$ (74)

Substituting into Eq. (69), terms in $\mathrm{e}^{0}$ give

$\hat{w}_{0}=\frac{|A|^{2}}{1-\beta},$ (75)

while terms in $\mathrm{e}^{2i\phi}$ and $\mathrm{e}^{-2i\phi}$ independently
give

$D=\alpha
A^{2},\quad\alpha=\frac{1}{2}\frac{(1+4ik_{||}\eta_{0}/\rho_{0}v_{A})}{(1-\beta+(8/3)ik_{||}\eta_{0}/\rho_{0}v_{A})}.$
(76)

The solution for $w_{0}$ can also be expressed without complex numbers. Making
explicit the real and imaginary parts of $\alpha=\alpha_{r}+i\alpha_{i}$, we
have the real constants

$\displaystyle\alpha_{r}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\frac{(1-\beta+(32/3)(\mathrm{Re})^{-2})}{((1-\beta)^{2}+(64/9)(\mathrm{Re})^{-2})},$
(77) $\displaystyle\alpha_{i}$ $\displaystyle=$
$\displaystyle\frac{2}{3}\frac{(1-3\beta)(\mathrm{Re})^{-1}}{((1-\beta)^{2}+(64/9)(\mathrm{Re})^{-2})},$
(78)

where $\mathrm{Re}=\rho_{0}v_{A}/k_{||}\eta_{0}$ consistent with Eq. (26). It
is then easily shown that

$\frac{w_{0}}{2R(Z)^{2}}=\alpha_{r}\cos(2(k_{||}(z-v_{A}t)+\theta(Z)))\\\
-\alpha_{i}\sin(2(k_{||}(z-v_{A}t)+\theta(Z)))+\frac{1}{2(1-\beta)}.$ (79)

To deduce $R(Z)$ and $\theta(Z)$, we return to analysing Eq. (58). The
$\epsilon^{1}$ terms in Eq. (58) give the inhomogeneous partial differential
equation

$\left(\frac{\partial^{2}}{\partial
t^{2}}-v_{A}^{2}\frac{\partial^{2}}{\partial
z^{2}}\right)v_{1}=2v_{A}^{2}\frac{\partial^{2}v_{0}}{\partial Z\partial
z}+v_{A}^{2}\frac{\partial^{2}}{\partial
z^{2}}\left(v_{0}w_{0}\right)+\frac{\eta_{0}}{\rho_{0}}\frac{\partial^{2}}{\partial
t\partial z}\left(\frac{\partial v_{0}^{3}}{\partial z}-2v_{0}\frac{\partial
w_{0}}{\partial z}\right).$ (80)

The $v_{0}$ and $w_{0}$ terms drive $v_{1}$, and the solution for $v_{1}$ will
have a secular contribution (i.e. one or more terms that grow relative to
corresponding solutions of the homogeneous equation) if terms on the right
hand side resonate with the solution to the undriven wave equation. In the
specific case where $v_{0}$ is given by Eq. (70), secular terms in the
solution for $v_{1}$ will restrict the domain for which $v_{0}$ is a valid
approximation if the right hand side of Eq. (80) contains $\mathrm{e}^{i\phi}$
or $\mathrm{e}^{-i\phi}$ terms. The central idea in multiple scale analysis is
to solve for the $A(Z)$ that makes the resonance disappear, making $v_{0}$ a
durable approximation for $v$.

Using Eqs. (70) and (74)–(76), $\mathrm{e}^{i\phi}$ terms vanish from the
right hand side of Eq. (80) if and only if

$\frac{1}{A}\frac{\mathrm{d}A}{\mathrm{d}Z}=-\frac{k_{||}|A|^{2}}{2}\left[i\left(\alpha+\frac{1}{1-\beta}\right)+\frac{k_{||}\eta_{0}}{\rho_{0}v_{A}}\left(3-4\alpha\right)\right].$
(81)

The same condition also removes the $\mathrm{e}^{-i\phi}$ terms.

Changing to polar form, Eq. (71) implies

$\frac{1}{A}\frac{\mathrm{d}A}{\mathrm{d}Z}=\frac{1}{R}\frac{\mathrm{d}R}{\mathrm{d}Z}+i\frac{\mathrm{d}\theta}{\mathrm{d}Z}.$
(82)

Hence, the real and imaginary parts of Eq. (81) yield the real ordinary
differential equations

$\displaystyle\frac{\mathrm{d}R}{\mathrm{d}Z}$
$\displaystyle=-\frac{\kappa_{1}}{2}R^{3},$ (83)
$\displaystyle\frac{\mathrm{d}\theta}{\mathrm{d}Z}$
$\displaystyle=\kappa_{2}R^{2},$ (84)

where

$\displaystyle\kappa_{1}$
$\displaystyle=k_{||}\left(\frac{k_{||}\eta_{0}}{\rho_{0}v_{A}}\left(3-4\alpha_{r}\right)-\alpha_{i}\right),$
(85) $\displaystyle\kappa_{2}$
$\displaystyle=-\frac{k_{||}}{2}\left(\alpha_{r}+\frac{1}{1-\beta}-\frac{k_{||}\eta_{0}}{\rho_{0}v_{A}}4\alpha_{i}\right).$
(86)

Eqs. (83) and (84) govern the local amplitude and phase drift of the Alfvén
wave respectively (c.f. Eq. (72)).

Our main interest is in $R(Z)$, which determines how the waves decay. Equation
(83) is a separable first order differential equation. The solution is

$R(Z)=\frac{R(0)}{\sqrt{1+\kappa_{1}R(0)^{2}Z}}.$ (87)

For $\kappa_{1}>0$ the wave envelope decays non-exponentially, over a damping
length that is inversely proportional to the square of the initial wave
amplitude. Having obtained $R(Z)$, the solution for $\theta(Z)$ is obtained by
directly integrating Eq. (84). Using Eq. (87),

$\theta(Z)=\theta(0)+\frac{\kappa_{2}}{\kappa_{1}}\ln\left|1+\kappa_{1}R(0)^{2}Z\right|.$
(88)

####  5.2.4 Solution in original variables

Having ensured corrections to $v\approx v_{0}$ remain of order
$\epsilon\sim(b/B_{0})^{2}\ll 1$, the multiple scale analysis is concluded by
using $v_{0}$ as the approximation for $v$. Returning to the original
variables,

$V_{x}(z,t)=\frac{a_{0}}{\sqrt{1+z/L_{d}}}\cos\left(k_{||}(z-v_{A}t)+(\kappa_{2}/\kappa_{1})\ln\left|1+z/L_{d}\right|+\theta_{0}\right),$
(89)

where $a_{0}$ is the amplitude of $V_{x}(0,t)$, $\theta_{0}$ sets the initial
phase of the wave (at $z=0$, $t=0$) and

$L_{d}=\frac{4}{\kappa_{1}(a_{0}/v_{A})^{2}}$ (90)

is the decay length. Using Eqs. (77), (78) and (85),

$\kappa_{1}=\frac{k_{||}}{3\mathrm{Re}}\frac{(1-3\beta)^{2}}{(1-\beta)^{2}+(64/9)(\mathrm{Re})^{-2}},$
(91)

where $\mathrm{Re}$ is the Reynolds number for the wave, as defined by Eq.
(26). Thus,

$L_{d}=\frac{12\mathrm{Re}}{k_{||}(a_{0}/v_{A})^{2}}\frac{(1-\beta)^{2}+(64/9)(\mathrm{Re})^{-2}}{(1-3\beta)^{2}}.$
(92)

If one neglects the $\mathrm{Re}^{-2}$ term in the numerator of Eq. (92), then
$L_{d}$ agrees exactly with the formula in Eq. (41) that we derived from
energy principles. The formula for $L_{d}$ in Eq. (92) is more general since
it was derived without direct assumptions about the value of
$\mathrm{Re}=k_{||}\eta_{0}/\rho_{0}v_{A}$, although the multiple scale
analysis requires that the combination of parameters $k_{||}$, $a_{0}^{2}$ and
$\mathrm{Re}$ are such that waves damp over a significantly longer scale than
the wavelength.

###  5.3 Non-exponential decay and interpretation of damping length

As a general principle, the Alfvén wave energy density $E_{w}=\rho V_{x}^{2}$
decays more rapidly than the perturbation $V_{x}$, due to the quadratic power.
For exponential decay this is reflected in a factor two difference in the
respective $\mathrm{e}^{-1}$ decay lengths. For the non-exponential decay
produced by nonlinear viscous damping, the situation is handled differently.
The same $L_{d}$ describes $V_{x}$ and $E_{w}$, however they have different
functional forms. The velocity amplitude decays as $(1+z/L_{d})^{-1/2}$ (see
Eq. (89)), while the wave energy density decays as $(1+z/L_{d})^{-1}$.
Therefore, over a distance $L_{d}$, the velocity amplitude reduces by a factor
$\sqrt{2}$ and the energy density halves.

###  5.4 Inclusion of thermal conduction

The multiple scale analysis can also be modified to include explicit thermal
conduction. Since thermal conduction is highly anisotropic, we include the
parallel thermal conduction, setting the heat flow vector to

$\mathbf{q}=-\mathcal{K}_{||}(\mathbf{h}\cdot\nabla T)\mathbf{h},$ (93)

where $\mathcal{K}_{||}$ is the coefficient of parallel thermal conduction,
and temperature $T=p/\rho\mathcal{R}$ where $\mathcal{R}$ is the gas constant.
The energy equation with heat flow is

$\frac{\partial p}{\partial t}+\mathbf{V}\cdot\nabla p+\gamma
p\nabla\cdot\mathbf{V}=(\gamma-1)(Q_{visc}-\nabla\cdot\mathbf{q}),$ (94)

which replaces Eq. (67).

It follows that $\delta p/p_{0}$ is related to $v_{z}/v_{A}$ by the partial
differential equation

$\left(1+\Lambda\frac{\partial}{\partial z}\right)\frac{\delta
p}{p_{0}}=\left(\gamma+\Lambda\frac{\partial}{\partial
z}\right)\frac{V_{z}}{v_{A}}+O(\epsilon^{2}).$ (95)

where

$\Lambda=\frac{(\gamma-1)\mathcal{K}_{||}}{\rho_{0}\mathcal{R}v_{A}}$ (96)

is a conductive length scale. In the limit of weak thermal conduction,
$\Lambda\to 0$ gives $\delta p/p_{0}=\gamma V_{z}/v_{A}$, recovering the
adiabatic case treated above. Similarly, for strong thermal conduction,
$\Lambda\to\infty$ gives $\delta p/p_{0}=V_{z}/v_{A}$, recovering the
isothermal case.

Introducing a dimensionless pressure variable $c$ defined by $\delta
p=\epsilon p_{0}c(z,t)$, expanding $c(z,t)=c_{0}(z,Z,t)+\epsilon
c_{1}(z,Z,t)+\ldots$ and setting
$c_{0}=C\mathrm{e}^{2i\phi}+C^{*}\mathrm{e}^{2i\phi}+\hat{c}_{0}$, terms in
$\mathrm{e}^{0}$ in Eq. (95) yield $\hat{c}_{0}=\gamma\hat{w}_{0}$, and terms
in $\mathrm{e}^{2i\phi}$ yield $C=\Gamma D$, where

$\Gamma=\frac{\gamma+2ik_{||}\Lambda}{1+2ik_{||}\Lambda}.$ (97)

Solving further, an equation $D=\alpha A^{2}$ analogous to Eq. (76) is
obtained but with $\beta$ replaced by the complex-valued $\Gamma
p_{0}/\rho_{0}v_{A}^{2}$ in the formula for $\alpha$. Meanwhile, (75) and (81)
are unchanged, retaining the real-valued $\beta=\gamma
p_{0}/\rho_{0}v_{A}^{2}$. The wave amplitude is therefore governed by results
identical to Eqs. (83) and (85), with the aforementioned change in the
definition of $\alpha$.

##  6 Discussion

###  6.1 Optimum damping

Inspecting Eq. (92), the formula for $L_{d}k_{||}$ has a minimum with respect
to the Alfvénic Reynolds number at $\mathrm{Re}=8/(3|1-\beta|)$. Thus, shear
Alfvén waves with $k_{||}=3\rho_{0}v_{A}|1-\beta|/3\eta_{0}$ are damped in the
fewest number of wavelengths, which we refer to as optimum damping. The
optimally damped waves have

$\frac{L_{d}}{\lambda_{||}}=\frac{32}{\pi}\frac{|1-\beta|}{(a_{0}/v_{A})^{2}(1-3\beta)^{2}}.$
(98)

When $\beta\ll 1$, the right hand side of Eq. (98) is approximately ten
divided by the square of the normalised wave amplitude. Hence, while nonlinear
viscous damping can in principle damp Alfvén waves in a small number of
wavelengths, this requires large amplitudes $a/v_{A}\sim 1$, or $\beta\sim 1$.
For a more typically encountered amplitudes $a/v_{A}\sim 10^{-1}$ and low
$\beta$, one finds $L_{d}/\lambda_{||}\gtrsim 1000$, making nonlinear viscous
damping negligible for many coronal situations.

###  6.2 Viscous self-organisation

The suppression of nonlinear viscous damping for small $\mathrm{Re}$ (highly
viscous plasma) does not mean that viscous effects are unimportant in this
regime. On the contrary, nonlinear damping is suppressed for small
$\mathrm{R_{e}}$ because viscous forces organise the parallel flow associated
with the Alfvén wave to approach the relationship
$V_{z}/v_{A}=(3/4)(V_{x}/v_{A})^{2}$. This modification of the parallel flow
plays a crucial role in avoiding significant nonlinear damping in highly
viscous plasma, when modelled using Braginskii MHD.

###  6.3 Validity constraints

Throughout this paper, we have assumed that $\beta\neq 1$ to avoid resonant
wave coupling. This condition holds throughout most of the corona, so it is
appropriate for our primary applications. Additionally, the multiple scale
analysis in § 5 uses $\epsilon\sim(V_{x}/v_{A})^{2}\sim(b/B_{0})^{2}$ as a
small parameter, one consequence of which is that the nonlinearly damping
occurs over a distance considerably greater than the parallel wavelength. As
noted in § 3.2, transverse coronal waves are observed in open-field regions
with $\epsilon\sim 10^{-2}$, making weakly nonlinear theory appropriate for
such situations. Obtaining nonlinear viscous solutions in the resonant and
strongly nonlinear regimes nonetheless remain interesting future challenges
for plasma theory.

Applicability of this paper’s results to physical problems is also constrained
to conditions under which Braginskii MHD can be rigorously applied. As
discussed in § 2, the traditional derivation of Braginskii MHD assumes that
the collisional mean free path is less than the macroscopic scales. Comparing
the mean free path parallel to the magnetic field to the parallel wavelength,
this condition can be given as $k_{||}\lambda_{mfp}<1$, where
$\lambda_{mfp}=v_{Ti}\tau_{i}$, $v_{Ti}=\sqrt{k_{B}T/m_{i}}$ and $\tau_{i}$ is
the ion collision time. Using the formula (Braginskii, 1958, 1965; Hollweg,
1985)

$\eta_{0}=0.96nk_{B}T\tau_{i},$ (99)

and the definition of $\mathrm{Re}$ in Eq. (26), one can show that

$k_{||}\lambda_{mfp}<1\ \Leftrightarrow\ \beta^{1/2}\mathrm{Re}=\frac{\rho
c_{s}}{k_{||}\eta_{0}}>1.$ (100)

In other words, Braginskii MHD requires that the Reynolds number based on the
sound speed is greater than unity. One should therefore be cautious about
applying small Alfvénic Reynolds number results such as viscous self-
organisation to real low-$\beta$ plasmas.

###  6.4 Formulas for applications

For applications to real plasmas, the following formulas are convenient. In
cases where the parallel viscosity coefficient is determined by Coulomb
collisions,

$\eta_{0}=0.96nk_{B}T\tau_{i}=\frac{22}{\lambda_{C}}\times 10^{-17}T^{5/2},$
(101)

where this formula is stated in S.I. units with $T$ in kelvin, and
$\lambda_{C}$ is the Coulomb logarithm (e.g. Hollweg, 1985). The Reynolds
number defined in Eq. (26) can then be expressed as

$\mathrm{Re}=5.8\times 10^{20}\lambda_{C}B^{2}f^{-1}T^{-5/2},$ (102)

also in S.I. units, where $f=v_{A}k_{||}/2\pi$ is the wave frequency. This
formula makes explicit the dependences on frequency, magnetic field strength
and temperature. The Alfvénic Reynolds number is smallest when the plasma has
high temperature and low magnetic field strength, and for higher frequency
waves. Finally, we express the damping length in Eq. (92) as a function of
frequency and the mean square velocity $\left<V_{x}^{2}\right>=a_{0}^{2}/2$,
which gives

$L_{d}=\frac{3}{\pi}\frac{v_{A}^{3}\mathrm{Re}}{f\left<V_{x}^{2}\right>}\frac{(1-\beta)^{2}+(64/9)(\mathrm{Re})^{-2}}{(1-3\beta)^{2}}.$
(103)

###  6.5 Waves in a coronal open-field region

Outgoing transverse waves in the magnetically open solar corona contain
sufficient energy to heat the open corona and accelerate the fast solar wind
(McIntosh et al., 2011; Morton et al., 2015), and they are observed to damp
significantly within a solar radius above the Sun’s surface (Bemporad & Abbo,
2012; Hahn et al., 2012; Hahn & Savin, 2013; Hahn et al., 2022). Heating at
these altitudes is also thought to be important for producing the observed
rapid acceleration of the fast solar wind (Habbal et al., 1995; McKenzie et
al., 1995). The problem of how the outgoing waves damp has not been
conclusively solved, although one leading hypothesis is turbulent cascade
driven by interactions with downgoing Alfvén waves (Hollweg, 1986; Heyvaerts &
Priest, 1992; Matthaeus et al., 1999; Cranmer et al., 2007; Verdini et al.,
2010; Mikić et al., 2018) produced either by reflection from density
inhomogeneities (van Ballegooijen & Asgari-Targhi, 2016; Pascoe et al., 2022)
or by parametric decay instability (Galeev & Oraevskii, 1963; Derby, 1978;
Goldstein, 1978; Shoda et al., 2019; Hahn et al., 2022).

Here, we demonstrate that Braginskii viscosity does not cause significant
damping of Alfvén waves at the altitudes at which the traditional derivation
of Braginskii MHD holds. For concreteness, we consider the Sun’s northern
polar open-field region on 27 March 2012, using observational values reported
by Morton et al. (2015). Enhanced wave power was present around $f=5\mbox{
mHz}$, which suggests Alfvénic waves produced by p-modes (Morton et al.,
2019). We will calculate damping lengths for this frequency, noting that
$\mathrm{Re}$ and $L_{d}$ depend on $f$, with $L_{d}\sim f^{-2}$ in the limit
of high $\mathrm{Re}$. Morton et al. (2015) inferred that the Alfvén speed was
nearly constant with $v_{A}=400\mbox{ km s}^{-1}$ on their domain of $r=1.05$
to $1.20R_{\sun}$. For temperature, we set $T=1.6\times 10^{6}\mbox{ K}$, the
formation temperature of the Fe XIII lines used by the CoMP instrument, which
implies the proton thermal speed $V_{Ti}=\sqrt{k_{B}T/m_{i}}$ is $115\mbox{ km
s}^{-1}$. Hence, in for an isothermal equation of state $\beta=0.083$ and
$\beta^{1/2}=0.29$. For the wave velocity amplitude, Morton et al. (2015)
recommended that the reported non-thermal line width should be used, which
varies with altitude.

Starting with lowest altitude observed by Morton et al. (2015),
$r=1.05R_{\sun}$, we set $n=10^{14}\mbox{ m}^{-3}$, $B=2\times 10^{-4}\mbox{
T}$ and take the rms value of $V_{x}$ as $35\mbox{ km s}^{-1}$. We therefore
find $\lambda_{C}=19$ and $\mathrm{Re}=28$. Since
$\beta^{1/2}\mathrm{Re}=8>1$, Braginksii MHD applies and we evaluate
$L_{d}=4.2\times 10^{8}\mbox{ km}\equiv 600R_{\sun}$.

At $r=1.20$, we set $n=10^{13}\mbox{ m}^{-3}$, $B=6\times 10^{-5}\mbox{ T}$
and take the rms value of $V_{x}$ as $50\mbox{ km s}^{-1}$. The observed
parameters therefore give $\lambda_{C}=21$ and $\mathrm{Re}=2.7$. Since
$\beta^{1/2}\mathrm{Re}=0.8\approx 1$, this altitude is close to the maximum
at which the assumptions by which Braginskii MHD is traditionally derived
remains valid (for this particular open field region, and assuming Eq. (101)).
Evaluating the damping length here returns $L_{d}=4.2\times 10^{7}\mbox{
km}\equiv 61R_{\sun}$.

We conclude that Braginskii viscosity does not cause significant wave damping
below $r=1.2R_{\sun}$, which is consistent with observational results that
Alfvénic wave amplitudes in coronal holes follow ideal WKB scaling out to
around this altitude (Cranmer & van Ballegooijen, 2005; Hahn & Savin, 2013).

Between the altitudes we have examined, $L_{d}$ reduces by an order of
magnitude. If one were to extrapolate using high $\mathrm{Re}$ or
incompressible results, it would appear that viscous damping becomes important
near the altitudes at which the waves are observed to damp. We are cautious
about making such an assertion for two reasons. First, as discussed in § 6.1,
our results show that for $\mathrm{Re}<8/(3|1-\beta))$ the damping length in a
Braginskii MHD model increases again as the field-aligned flow self organises
to supress viscous damping. Secondly, as the plasma becomes increasingly
collisionless ($\beta^{1/2}\mathrm{Re}<1$) the traditional derivation of
Bragniskii MHD falters.

Intriguingly, it may be significant that the onset of wave damping broadly
coincides with the altitude at which Braginskii MHD can no longer be
confidently applied if one invokes the $\eta_{0}$ expression for Coulomb
collisions given in Eq. (101). This correspondence is suggestive that the wave
damping observed in coronal holes may involve collisionless and heat flow
effects not found in the most common fluid models.

###  6.6 Future work

The present types of analyses should be extended in future to other types of
propagating transverse MHD waves. The nonlinear longitudinal flow that
accompanies propagating torsional Alfvén waves differs from its counterpart
for propagating shear Alfvén waves (Vasheghani Farahani et al., 2011) and it
will be of interest to investigate how this difference affects the nonlinear
viscous damping. It is similarly desirable to determine how nonlinear
viscosity affects propagating kink waves (Edwin & Roberts, 1983).

For propagating shear Alfvén waves, viscous damping appears most promising
near the $c_{s}=v_{A}$ singularity, which must be treated using different
methods to those used in this paper. The solar wind frequently has $\beta\sim
1$, while $\beta=1$ occurs in the lower solar atmosphere and in the vicinity
of coronal nulls points. Hence, this case is of considerable physical
interest. One challenge for application to magnetic nulls is that the magnetic
field unit vector $\mathbf{h}=\mathbf{B}/|\mathbf{B}|$ is not defined at the
null itself, so one must be careful to evaluate the Braginskii stress tensor
using appropriate calculations, e.g. see recent discussion by MacTaggart et
al. (2017).

A further challenge is to develop a theory of nonlinear viscous damping
applicable to strongly nonlinear waves with amplitudes $b\sim B_{0}$ and
greater. The results of the multiple-scale analysis in § 5 are rigorous only
for the weakly nonlinear case, in which $\epsilon\sim(b/B_{0})^{2}$ can be
treated as a small parameter and it is assumed that the damping length is
significantly longer than the wavelength. Strongly nonlinear Alfvén waves with
$b\sim B_{0}$ are a feature of the solar wind, and while the low
collisionality of the solar wind means that Braginskii MHD may not be an
appropriate framework for that application, extending the current work to
strongly nonlinear waves remains an interesting problem.

There is a diverse collection of MHD wave problems beyond wave damping for
which viscous effects are likely to be significant. Prime among these are
nonlinear phenomena involving Alfvén waves, for which the nonlinear viscosity
tensor enters the equations at the same order as the effect of interest. For
example, standing Alfvén waves drive significantly stronger field-aligned
flows than occur for propagating waves because standing Alfvén waves create
inhomogeneous time-averaged magnetic pressure. There could also be significant
value in investigating how viscosity modifies wave interactions, including
Alfvén wave collisions and parametric decay instability (Galeev & Oraevskii,
1963; Derby, 1978; Goldstein, 1978), which are central to leading hypotheses
of wave heating in the magnetically open solar corona.

Finally, we point to the continuing need for basic plasma physics research to
provide increasingly rigorous derivation and validation of the appropriate
fluid equations for weakly collisional and collisionless plasma, in the face
of the closure problem summarised in § 2. As discussed in § 2 and 6.3,
Braginskii MHD breaks down at higher altitudes in the corona as the plasma
becomes increasingly collisionless (see Eqs. (100) and (102)). The CGL double-
adiabatic equations and other models that evolve the stress tensor may provide
a more suitable framework in these conditions. Hunana et al. (2019a, b, 2022)
provide recent discussions of such models and their limitations.
Alternatively, it may be necessary for the solar waves community to more
widely adopt non-fluid plasma models. However, tractability of kinetic models
remains a limiting factor, especially in light of the large separations
between kinetic and macroscopic scales that are characteristic of the Sun’s
corona. Eloquent comments on these matters can be found in Montgomery (1996).

##  7 Conclusions

This paper has investigated the properties of propagating shear Alfvén waves
subject to the nonlinear effects of the Braginskii viscous stress tensor. The
main points are as follows:

  1. 1.

For many plasma environments, including the low-altitude solar corona,
Braginskii MHD provides a more accurate description of plasma than classical
MHD does, by rigorously treating the stress tensor and thermal conduction.
Stress tensor effects nonetheless remain relatively unexplored for many solar
MHD phenomena.

  2. 2.

The dominant viscous effects for propagating shear Alfvén waves are nonlinear
in the wave amplitude and occur through the “parallel” viscosity coefficient,
$\eta_{0}$. Theoretical results based on linearizing the stress tensor with
respect to the wave amplitude are only valid for amplitudes satisfying
$(b/B_{0})^{2}\ll(\Omega_{i}\tau_{i})^{-2}$. Such waves would be energetically
insignificant under normal coronal conditions, hence nonlinear treatment is
required.

  3. 3.

Compressibility and pressure affect the nonlinear field-aligned flow
associated with shear Alfvén waves, hence they impact nonlinear wave damping.
Both must be included to produce accurate coronal results.

  4. 4.

Braginskii viscosity damps propagating shear Alfvén waves nonlinearly, such
that the primary wave fields $b$ and $V_{x}$ decay as $(1+z/L_{d})^{-1/2}$,
where the decay length

$L_{d}=\frac{12\mathrm{Re}}{k_{||}(a_{0}/v_{A})^{2}}\frac{(1-\beta)^{2}+(64/9)(\mathrm{Re})^{-2}}{(1-3\beta)^{2}}.$

Here, $a_{0}$ is the initial velocity amplitude of the wave,
$\beta=(c_{s}/v_{A})^{2}$ and $\mathrm{Re}=\rho v_{A}/k_{||}\eta_{0}$ is the
Alfvénic Reynolds number of the wave. The energy density decays as
$(1+z/L_{d})^{-1}$.

  5. 5.

Optimal damping (the minimum normalised damping length $k_{||}L_{d}$) is
obtained when $\mathrm{Re}=8/(3|1-\beta|)$. For low $\beta$ plasma and
$(a_{0}/v_{A})\lesssim 10^{-1}$, one finds $L_{d}/\lambda_{||}\gtrsim 1000$,
indicating that nonlinear viscous damping is negligible for many coronal
situations.

  6. 6.

The asymptotic behaviour that $L_{d}\to\infty$ in the highly viscous regime
$\mathrm{Re}\to 0$ is attributed to self-organisation of the parallel flow by
viscous forces such that $V_{z}/v_{A}\approx(3/4)(V_{x}/v_{A})^{2}$, which
suppresses dissipation.

  7. 7.

Applicability of the Braginskii MHD solutions to real plasmas is constrained
by the traditional derivation of Braginskii MHD assuming that
$k_{||}\lambda_{mfp}<1$ which is equivalent to $\beta^{1/2}\mathrm{Re}=\rho
c_{s}/k_{||}\eta_{0}>1$. In other words, Braginskii MHD requires that the
Reynolds number based on the sound speed is greater than unity. We therefore
recommend that only the damping results for large Alfvénic Reynolds number
should be applied to real coronal plasma, using the simplified formula
$L_{d}=12\mathrm{Re}(1-\beta)^{2}/(k_{||}(a_{0}/v_{A})^{2}(1-3\beta)^{2}))$
that has been derived in this paper by two different techniques: energy
principles and multiple scale analysis.

  8. 8.

Application to transverse waves observed in a polar open-field region
concludes that nonlinear Braginskii viscosity does not cause significant
damping of the waves at the altitudes at which the assumptions by which
Braginskii MHD is traditionally derived remain valid ($r\lesssim 1.2R_{\sun}$
for the considered region and wave properties). Intriguingly, the observed
onset of wave damping broadly coincides with the altitude at which Braginskii
MHD can no longer be confidently applied if one invokes the $\eta_{0}$
expression for Coulomb collisions given in Eq. (101).

This work was prompted by and benefited from conversations with Paola Testa,
Bart De Pontieu, Vanessa Polito, Graham Kerr, Mark Cheung, Wei Liu, David
Graham, Joel Allred, Mats Carlsson, Iain Hannah and Fabio Reale during a
research visit to LMSAL funded by ESA’s support for the IRIS mission (August
2018) and meetings of International Space Science Institute (Bern)
International Team 355 (November 2018). I am grateful to Peter Cargill, Andrew
Wright, Bart De Pontieu and Paola Testa for comments on an early draft (June
2019), and Declan Diver for encouragment to explore connections with pressure
anisotropy. I thank the reviewer for considered and constructive suggestions,
and several unnamed individuals for comments that also improved the
manuscript. The research made use of the NRL Plasma Formulary, NASA’s ADS
Bibliographic Services and the British Library On Demand service with
assistance from librarians at the University of Dundee.

## References

  * Alfvén (1942) Alfvén, H. 1942, Nature, 150, 405, doi: 10.1038/150405d0
  * Alfvén (1943) —. 1943, Arkiv for Matematik, Astronomi och Fysik, 29B, 1
  * Alfvén (1947) —. 1947, MNRAS, 107, 211, doi: 10.1093/mnras/107.2.211
  * Alfvén (1950) —. 1950, Cosmical electrodynamics (Clarendon Press)
  * Arregui (2015) Arregui, I. 2015, Philosophical Transactions of the Royal Society of London Series A, 373, 20140261, doi: 10.1098/rsta.2014.0261
  * Balescu (1988) Balescu, R. 1988, Transport processes in plasmas. (Amsterdam: North-Holland)
  * Batchelor (1950) Batchelor, G. K. 1950, Proceedings of the Royal Society of London Series A, 201, 405, doi: 10.1098/rspa.1950.0069
  * Belcher & Davis (1971) Belcher, J. W., & Davis, Leverett, J. 1971, J. Geophys. Res., 76, 3534, doi: 10.1029/JA076i016p03534
  * Bemporad & Abbo (2012) Bemporad, A., & Abbo, L. 2012, ApJ, 751, 110, doi: 10.1088/0004-637X/751/2/110
  * Bender & Orszag (1978) Bender, C. M., & Orszag, S. A. 1978, Advanced Mathematical Methods for Scientists and Engineers (Springer)
  * Braginskii (1958) Braginskii, S. I. 1958, Soviet Journal of Experimental and Theoretical Physics, 6, 358
  * Braginskii (1965) —. 1965, Reviews of Plasma Physics, 1, 205
  * Chapman & Cowling (1939) Chapman, S., & Cowling, T. G. 1939, The Mathematical Theory of Non-Uniform Gases (Cambridge University Press)
  * 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
  * Coleman (1967) Coleman, Paul J., J. 1967, Planet. Space Sci., 15, 953, doi: 10.1016/0032-0633(67)90166-3
  * Cranmer & van Ballegooijen (2005) Cranmer, S. R., & van Ballegooijen, A. A. 2005, ApJS, 156, 265, doi: 10.1086/426507
  * Cranmer et al. (2007) Cranmer, S. R., van Ballegooijen, A. A., & Edgar, R. J. 2007, ApJS, 171, 520, doi: 10.1086/518001
  * De Moortel & Browning (2015) De Moortel, I., & Browning, P. 2015, Philosophical Transactions of the Royal Society of London Series A, 373, 20140269, doi: 10.1098/rsta.2014.0269
  * De Moortel et al. (2020) De Moortel, I., Falconer, I., & Stack, R. 2020, Astronomy and Geophysics, 61, 2.34, doi: 10.1093/astrogeo/ataa031
  * De Moortel & Nakariakov (2012) De Moortel, I., & Nakariakov, V. M. 2012, Philosophical Transactions of the Royal Society of London Series A, 370, 3193, doi: 10.1098/rsta.2011.0640
  * De Pontieu et al. (2007) De Pontieu, B., McIntosh, S. W., Carlsson, M., et al. 2007, Science, 318, 1574, doi: 10.1126/science.1151747
  * Derby (1978) Derby, N. F., J. 1978, ApJ, 224, 1013, doi: 10.1086/156451
  * Edwin & Roberts (1983) Edwin, P. M., & Roberts, B. 1983, Sol. Phys., 88, 179, doi: 10.1007/BF00196186
  * Erdelyi & Goossens (1995) Erdelyi, R., & Goossens, M. 1995, A&A, 294, 575
  * Galeev & Oraevskii (1963) Galeev, A. A., & Oraevskii, V. N. 1963, Soviet Physics Doklady, 7, 988
  * Goldstein (1978) Goldstein, M. L. 1978, ApJ, 219, 700, doi: 10.1086/155829
  * Grad (1949) Grad, H. 1949, Communications on Pure and Applied Mathematics, 2, 331, doi: https://doi.org/10.1002/cpa.3160020403
  * Habbal et al. (1995) Habbal, S. R., Esser, R., Guhathakurta, M., & Fisher, R. R. 1995, Geophys. Res. Lett., 22, 1465, doi: 10.1029/95GL01064
  * Hahn et al. (2022) Hahn, M., Fu, X., & Savin, D. W. 2022, ApJ, 933, 52, doi: 10.3847/1538-4357/ac7147
  * Hahn et al. (2012) Hahn, M., Landi, E., & Savin, D. W. 2012, ApJ, 753, 36, doi: 10.1088/0004-637X/753/1/36
  * Hahn & Savin (2013) Hahn, M., & Savin, D. W. 2013, ApJ, 776, 78, doi: 10.1088/0004-637X/776/2/78
  * Hartmann (1937) Hartmann, J. 1937, Mathematisk-fysiske Meddelelser. XV, 15, 6
  * Heyvaerts & Priest (1992) Heyvaerts, J., & Priest, E. R. 1992, ApJ, 390, 297, doi: 10.1086/171280
  * Hogan (1984) Hogan, J. T. 1984, Physics of Fluids, 27, 2308, doi: 10.1063/1.864886
  * Hollweg (1971) Hollweg, J. V. 1971, J. Geophys. Res., 76, 5155, doi: 10.1029/JA076i022p05155
  * Hollweg (1985) —. 1985, J. Geophys. Res., 90, 7620, doi: 10.1029/JA090iA08p07620
  * Hollweg (1986) —. 1986, J. Geophys. Res., 91, 4111, doi: 10.1029/JA091iA04p04111
  * Hollweg (1987) —. 1987, ApJ, 320, 875, doi: 10.1086/165604
  * Hunana et al. (2019a) Hunana, P., Tenerani, A., Zank, G. P., et al. 2019a, Journal of Plasma Physics, 85, 205850602, doi: 10.1017/S0022377819000801
  * Hunana et al. (2019b) —. 2019b, Journal of Plasma Physics, 85, 205850603, doi: 10.1017/S0022377819000850
  * Hunana et al. (2022) Hunana, P., Passot, T., Khomenko, E., et al. 2022, ApJS, 260, 26, doi: 10.3847/1538-4365/ac5044
  * Lifshitz & Pitaevskii (1981) Lifshitz, E. M., & Pitaevskii, L. P. 1981, Physical kinetics (Butterworth-Heinemann)
  * Lin et al. (2007) Lin, Y., Engvold, O., Rouppe van der Voort, L. H. M., & van Noort, M. 2007, Sol. Phys., 246, 65, doi: 10.1007/s11207-007-0402-8
  * MacTaggart et al. (2017) MacTaggart, D., Vergori, L., & Quinn, J. 2017, Journal of Fluid Mechanics, 826, 615, doi: 10.1017/jfm.2017.463
  * Matthaeus et al. (1999) Matthaeus, W. H., Zank, G. P., Oughton, S., Mullan, D. J., & Dmitruk, P. 1999, ApJ, 523, L93, doi: 10.1086/312259
  * McIntosh et al. (2011) McIntosh, S. W., de Pontieu, B., Carlsson, M., et al. 2011, Nature, 475, 477, doi: 10.1038/nature10235
  * McKenzie et al. (1995) McKenzie, J. F., Banaszkiewicz, M., & Axford, W. I. 1995, A&A, 303, L45
  * McLaughlin et al. (2011) McLaughlin, J. A., de Moortel, I., & Hood, A. W. 2011, A&A, 527, A149, doi: 10.1051/0004-6361/201015552
  * Mikić et al. (2018) Mikić, , Z., Downs, C., et al. 2018, Nature Astronomy, 2, 913, doi: 10.1038/s41550-018-0562-5
  * Montgomery (1992) Montgomery, D. 1992, J. Geophys. Res., 97, 4309, doi: 10.1029/92JA00144
  * Montgomery (1996) —. 1996, Journal of Plasma Physics, 56, 387, doi: 10.1017/S0022377800019358
  * Morton et al. (2015) Morton, R. J., Tomczyk, S., & Pinto, R. 2015, Nature Communications, 6, 7813, doi: 10.1038/ncomms8813
  * Morton et al. (2019) Morton, R. J., Weberg, M. J., & McLaughlin, J. A. 2019, Nature Astronomy, 3, 223, doi: 10.1038/s41550-018-0668-9
  * Nakariakov & Verwichte (2005) Nakariakov, V. M., & Verwichte, E. 2005, Living Reviews in Solar Physics, 2, 3, doi: 10.12942/lrsp-2005-3
  * Nocera et al. (1986) Nocera, L., Priest, E. R., & Hollweg, J. V. 1986, Geophysical and Astrophysical Fluid Dynamics, 35, 111, doi: 10.1080/03091928608245889
  * Ofman et al. (1994) Ofman, L., Davila, J. M., & Steinolfson, R. S. 1994, ApJ, 421, 360, doi: 10.1086/173654
  * Okamoto et al. (2007) Okamoto, T. J., Tsuneta, S., Berger, T. E., et al. 2007, Science, 318, 1577, doi: 10.1126/science.1145447
  * Osterbrock (1961) Osterbrock, D. E. 1961, ApJ, 134, 347, doi: 10.1086/147165
  * Oughton (1996) Oughton, S. 1996, Journal of Plasma Physics, 56, 641, doi: 10.1017/S0022377800019504
  * Oughton (1997) —. 1997, Journal of Plasma Physics, 58, 571, doi: 10.1017/S0022377897005928
  * Pascoe et al. (2022) Pascoe, D. J., De Moortel, I., Pagano, P., & Howson, T. A. 2022, MNRAS, 516, 2181, doi: 10.1093/mnras/stac2294
  * Priest (2014) Priest, E. 2014, Magnetohydrodynamics of the Sun (Cambridge University Press), doi: 10.1017/CBO9781139020732
  * Ruderman (1991) Ruderman, M. S. 1991, Sol. Phys., 131, 11, doi: 10.1007/BF00151740
  * Russell (2018) Russell, A. J. B. 2018, Sol. Phys., 293, 83, doi: 10.1007/s11207-018-1296-3
  * Russell et al. (2016) Russell, A. J. B., Mooney, M. K., Leake, J. E., & Hudson, H. S. 2016, ApJ, 831, 42, doi: 10.3847/0004-637X/831/1/42
  * Schunk & Nagy (2009) Schunk, R., & Nagy, A. 2009, Ionospheres: Physics, Plasma Physics, and Chemistry (Cambridge University Press), doi: 10.1017/CBO9780511635342
  * Shoda et al. (2019) Shoda, M., Suzuki, T. K., Asgari-Targhi, M., & Yokoyama, T. 2019, ApJ, 880, L2, doi: 10.3847/2041-8213/ab2b45
  * Spitzer (1962) Spitzer, L. 1962, Physics of Fully Ionized Gases, Second edn. (Interscience Publishers, reprinted by Dover Publications)
  * Steinolfson et al. (1986) Steinolfson, R. S., Priest, E. R., Poedts, S., Nocera, L., & Goossens, M. 1986, ApJ, 304, 526, doi: 10.1086/164187
  * Tomczyk et al. (2007) Tomczyk, S., McIntosh, S. W., Keil, S. L., et al. 2007, Science, 317, 1192, doi: 10.1126/science.1143304
  * van Ballegooijen & Asgari-Targhi (2016) van Ballegooijen, A. A., & Asgari-Targhi, M. 2016, ApJ, 821, 106, doi: 10.3847/0004-637X/821/2/106
  * Van Doorsselaere et al. (2020) Van Doorsselaere, T., Srivastava, A. K., Antolin, P., et al. 2020, Space Sci. Rev., 216, 140, doi: 10.1007/s11214-020-00770-y
  * Vasheghani Farahani et al. (2011) Vasheghani Farahani, S., Nakariakov, V. M., van Doorsselaere, T., & Verwichte, E. 2011, A&A, 526, A80, doi: 10.1051/0004-6361/201016063
  * Verdini et al. (2010) Verdini, A., Velli, M., Matthaeus, W. H., Oughton, S., & Dmitruk, P. 2010, ApJ, 708, L116, doi: 10.1088/2041-8205/708/2/L116
  * Zank (2014) Zank, G. P. 2014, Transport Processes in Space Physics and Astrophysics, Vol. 877 (Springer), doi: 10.1007/978-1-4614-8480-6