eurm10 msam10

# On the construction of general large-amplitude spherically polarised Alfvén
waves

J. Squire1 Email address for correspondence<EMAIL_ADDRESS>A.
Mallet2 1Physics Department, University of Otago, Dunedin 9010, New Zealand
2Space Sciences Laboratory, University of California, Berkeley, CA 94720, USA

(2022; )

###### Abstract

In a magnetised plasma on scales well above ion kinetic scales, any constant-
magnitude magnetic field, accompanied by parallel Alfvénic flows, forms a
nonlinear solution in an isobaric, constant-density background. These
structures, which are also known as spherically polarised Alfvén waves, are
observed ubiquitously in the solar wind, presumably created by the growth of
small-amplitude fluctuations as they propagate outwards in the corona. Here,
we present a computational method to construct such solutions of arbitrary
amplitude with general multi-dimensional structure, and explore some of their
properties. The difficulty lies in computing a zero-divergence, constant-
magnitude magnetic field, which leaves a single, quasi-free function to define
the solution, while requiring strong constraints on any individual component
of the field. Motivated by the physical process of wave growth in the solar
wind, our method circumvents this issue by starting from low-amplitude
Alfvénic fluctuations dominated by a strong mean field, then “growing”
magnetic perturbations into the large-amplitude regime. We present example
solutions with nontrivial structure in one, two, and three dimensions,
demonstrating a clear tendency to form very sharp gradients or
discontinuities, unless the solution is one dimensional. As well as being
useful as an input for other calculations, particularly the study of
parametric decay, our results provide a natural explanation for the extremely
sharp field discontinuities observed across magnetic-field switchbacks in the
low solar wind.

##  1 Introduction

The existence of the incompressible Alfvén wave is perhaps the most important
distinguishing feature between the dynamics of neutral fluids and plasmas,
underlying key physics of turbulent energy dissipation and enabling nontrivial
dynamics well below the scale of the mean-free path (e.g., Schekochihin et
al., 2009). Interestingly—and unlike other wave-like perturbations of fluids
or plasmas—Alfvén waves have an unambiguous nonlinear generalisation (Barnes &
Hollweg, 1974): defining the plasma density $\rho$, pressure $P$, magnetic-
field $\bm{B}$, flow velocity $\bm{u}$, and field strength $B=|\bm{B}|$, any
perturbation that satisfies

$P={\rm const.},\quad\rho={\rm const.},\quad B^{2}={\rm
const.},\quad\delta\bm{u}=\pm\delta\bm{B}/\sqrt{4\pi\rho}$ (1)

is a nonlinear solution that propagates in the direction
$\hat{\bm{b}}=\bm{B}/B$ at speed ${v}_{\rm
A}\equiv|\overline{\bm{B}}|/\sqrt{4\pi\rho}$ (where the overline signifies a
spatial average, and $\delta\bm{B}=\bm{B}-\overline{\bm{B}}$). Such solutions
are necessarily maximally “imbalanced” (the magnitude of the fluctuation’s
cross helicity equals their energy), and, in the wave frame at speed ${v}_{\rm
A}$, involve constant plasma kinetic energy (constant $|\bm{u}|$) and zero
motional electric field $\bm{u}\times\bm{B}$ (Matteini et al., 2015).
Moreover, unlike large-amplitude sound or magnetosonic waves, spherically
polarised perturbations do not steepen into shocks even if
$|\delta\bm{B}|\gg|\overline{\bm{B}}|$, although they can be unstable (Sagdeev
& Galeev, 1969; Cohen & Dewar, 1974). They also share many properties of the
linear (small-amplitude) Alfvénic fluctuations, such as how they change in
amplitude and/or refract in a nonhomogenous background medium (e.g., Barnes &
Hollweg, 1974; Hollweg, 1974). The solution (1) is even valid in a
collisionless plasma as well as in collisional magnetohydrodynamics (MHD), so
long as $\delta\bm{B}$ varies over scales much larger than the ion gyroradius
or skin depth (Barnes & Suffolk, 1971; Kulsrud, 1983). Given these properties,
it is perhaps not surprising that perturbations close to Eq. 1 are observed
ubiquitously in our best-studied example of a natural plasma—the solar wind
(Belcher & Davis, 1971). In this context, they are known as _spherically
polarised_ and their properties likely underly a number of key aspects of
solar-wind physics, including heating and acceleration (e.g., Shoda et al.,
2021; Bale et al., 2021), turbulent spectra (e.g., Matteini et al., 2018;
Bowen et al., 2021), and properties of “switchback” field reversals (Kasper et
al., 2019; Squire et al., 2020; Johnston et al., 2022).

Despite these important features, to our knowledge, there does not exist any
general method to construct and study such nonlinear Alfvénic solutions when
$\delta\bm{B}\gtrsim|\bm{B}|$, or even any general results regarding their
existence. Previous methods (Valentini et al., 2019; Roberts, 2012, see also),
though promising, may lead to unnecessary discontinuities in the solutions,
the causes of which are discussed below. The difficulty arises from the twin
constraints of $\nabla\cdot\bm{B}=0$ and $B=|\bm{B}|={\rm const.}$, which
leaves only one degree of freedom from which to construct the 3-D vector field
$\bm{B}$. Moreover, as we show below, this degree of freedom cannot in general
be freely chosen when $\delta\bm{B}$ approaches $\overline{\bm{B}}$ in
magnitude, even though simple, smooth solutions do exist (they just require
nontrivial constraints on all components of $\bm{B}$ simultaneously). It is
the purpose of this letter to present such a method and briefly explore the
properties of the solutions (1). These are of interest for two reasons. First,
the method can be used as input or initial conditions for other numerical
calculations, such as the study of parametric decay (Del Zanna, 2001;
Primavera et al., 2019) or large-amplitude reflection-driven turbulence
(Squire et al., 2020; Johnston et al., 2022). Second, the solutions are
interesting in and of themselves—our method shares strong similarities with
the plasma-expansion process that generates large-amplitude Alfvénic states in
the solar wind, so the structures and properties that arise may be of direct
physical relevance. Of particular interest is the development of sharp field
discontinuities, which we show form much more readily 2-D or 3-D solutions
than in 1-D; this is promising, since highly discontinuous fields seem to be
observed in switchbacks by Parker Solar Probe and other spacecraft (Bale et
al., 2019; Akhavan-Tafti et al., 2021). More generally, the method provides an
extremely broad class of nonlinear solutions to the MHD equations that
complements other solutions such as force-free magnetostatic equilibria
(Marsh, 1996). Seen differently, it provides a way to construct a general
divergence-free, unit-vector field.

Below, we first discuss the basic idea of the method and equations involved,
then present its application to magnetic-field configurations that vary in
one, two, and three dimensions, respectively. The one-dimensional version,
although idealised and explored previously in Mallet et al. (2021, hereafter
M+21) and Squire et al. (2022), provides a helpful illustration of the method
and the challenges involved in higher dimensions. We finish with brief
discussion of some interesting features of the solutions, focused particularly
on the generation of sharp gradient and/or discontinuities, as well as some
possible applications of our results.

##  2 Method

The method we propose is based on splitting the magnetic field into its mean
($\overline{\bm{B}}$) and fluctuating ($\delta\bm{B}$) parts, then devising an
equation to grow $\delta\bm{B}$ in amplitude while maintaining
$\overline{\delta\bm{B}}=0$, $\nabla\cdot\delta\bm{B}=0$, fixed
$\overline{\bm{B}}$, and constant
$B^{2}=|\overline{\bm{B}}+\delta\bm{B}|^{2}$. In this way, the method “grows”
a small amplitude spherically polarised wave, which can be relatively easily
constructed from linear Alfvénic perturbations via standard optimization
methods, to any desired amplitude
$\mathcal{A}\equiv(\overline{\delta\bm{B}^{2}}/\overline{\bm{B}}^{2})^{1/2}$.
The obvious candidate to evolve $\delta\bm{B}$ is simply the induction
equation supplemented by exponential growth,

$\frac{\partial}{\partial
t}\delta\bm{B}=\delta\bm{B}+\nabla\times(\tilde{\bm{u}}\times\delta\bm{B}),$
(2)

which clearly maintains $\overline{\delta\bm{B}}=0$ and
$\nabla\cdot\delta\bm{B}=0$ by construction. The flow $\tilde{\bm{u}}$ should
be non-Alfvénic: this is not the flow associated with the Alfvén wave itself
($\delta\bm{u}$ in Eq. 1), but rather the flow needed to change the shape of
$\delta\bm{B}$ as it grows with $t$, in order to maintain constant $B$. The
clear choice is the potential/compressive flow $\tilde{\bm{u}}=\nabla\phi$,
which cannot contribute to the Alfvénic flow on a periodic domain because
$\nabla\cdot\tilde{\bm{u}}=\nabla^{2}\phi$. To proceed, we form
${\partial}/{\partial
t}|\overline{\bm{B}}+\delta\bm{B}|^{2}=2(\overline{\bm{B}}+\delta\bm{B})\cdot{\partial\delta\bm{B}}/{\partial
t}$ and assert that its fluctuating part must be zero, _viz.,_ that $\phi$
should be chosen so that $\partial/\partial t\,\delta(B^{2})=0$. Using the
fact that $\nabla B^{2}=0$ and $\nabla\overline{\bm{B}}=0$, this yields

$B^{2}\nabla_{\perp}^{2}\phi=B^{2}\nabla^{2}\phi-\sum_{ij}B_{i}B_{j}\nabla_{i}\nabla_{j}\phi=-\delta\bm{B}\cdot\overline{\bm{B}}.,$
(3)

where $B_{i}=\overline{B}_{i}+\delta B_{i}$. If Eq. 3 can be solved at each
step, with the result used to evolve $\delta\bm{B}$ through Eq. 2, the total
field will maintain spatially constant $B^{2}$ (and $\nabla\cdot\bm{B}=0$) as
$\delta\bm{B}$ grows to $\mathcal{A}\gg 1$.

This method was motivated by the reduced equations of M+21, which describe the
evolution of 1-D spherically polarised Alfvénic states in the expanding solar-
wind plasma. In this case, the amplitudes of the fluctuating field and
(radial) mean field evolve as $\propto a^{-1/2}$ and $\propto a^{-1}$,
respectively, where $a$ is the plasma’s expansion factor. This leads to a
similar situation where $\delta\bm{B}$ grows compared to $\overline{\bm{B}}$.
Expansion also causes the gradient operators in Eq. 2 to differ in the
perpendicular and parallel directions and change with $a$, as well as rotating
the mean field if it has non-radial components. These effects necessitate
extra terms in Eqs. 2 and 3 (see § 3 for details), but do not fundamentally
modify the method. Such terms cause intriguingly different nonlinear solutions
to develop (see, e.g., figure 4 of Squire et al. 2022) pointing to interesting
generalisations of our method, such as making the direction of
$\overline{\bm{B}}$ change in time or modifying the gradient operators to
capture the physical effects of super-radial expansion and wind acceleration
relevant to the inner heliosphere (Tenerani & Velli, 2017). The correspondence
with M+21 also shows that, while Eqs. 2 and 3 are formulated on a purely
mathematical basis, the process of $\delta\bm{B}$ growth does share important
similarities with the real physical processes that generate large-amplitude
spherically polarised states in the solar wind.

##  3 One dimension

In one-dimensional solutions, $\delta\bm{B}$ varies only in the $\hat{\bm{p}}$
direction, which can be at an arbitrary angle to $\overline{\bm{B}}$. This
case is significantly more straightforward than two or three-dimensional
solutions and serves to illustrate some useful points. Denoting the co-
ordinate in the $\hat{\bm{p}}$ direction as $\lambda$, such that
$\nabla\phi=\hat{\bm{p}}d\phi/d\lambda=\hat{\bm{p}}\phi^{\prime}$, where
$\hat{\bm{p}}$ is a unit vector, Eqs. 2 and 3 simplify to

$\frac{\partial\delta\bm{B}}{\partial
t}=\delta\bm{B}-\frac{\partial}{\partial\lambda}\left[\phi^{\prime}(\delta\bm{B}+\overline{\bm{B}}_{{\rm
T}})\right],\quad|\delta\bm{B}+\overline{\bm{B}}_{{\rm
T}}|^{2}\phi^{\prime\prime}=-\delta\bm{B}\cdot\overline{\bm{B}}.$ (4)

Here $\overline{\bm{B}}_{{\rm
T}}=\overline{\bm{B}}-\hat{\bm{p}}(\hat{\bm{p}}\cdot\overline{\bm{B}})$ is the
part of $\overline{\bm{B}}$ that is transverse to the mean field. We see a
clear correspondence with equations (59) and (61) of M+21 (with
$\overline{\bm{B}}_{{\rm T}}$ denoted $\bm{v}_{\rm AT}$), which can actually
be equivalently derived by asserting that $\partial/\partial
t\,\delta(B^{2})=0$ (i.e., in the same way as Eqs. 2 and 3), as opposed to the
asymptotic expansion in slow expansion rate used therein. The extra terms in
M+21 result from expansion induced rotation of $\hat{\bm{p}}$ and
$\overline{\bm{B}}$, as mentioned above.

Figure 1: One-dimensional spherically polarised solutions from Eq. 4, starting
from a solution of Eq. 5 with $\mathcal{A}\approx 0.2$ constructed from a
random collection of the first 6 Fourier modes. The wavevector $\hat{\bm{p}}$
is at an angle of $30^{\circ}$ to $\overline{\bm{B}}$, which lies in the
$\hat{x}$ direction. Colors show $B_{x}$ (blue), $B_{y}$ (red), $B_{z}$
(yellow), and $B$ (black). From top to bottom, the panels show
$\mathcal{A}\approx 0.2$ (the initial conditions), $\mathcal{A}\approx 0.65$,
$\mathcal{A}\approx 5$, and $\mathcal{A}\approx 400$, illustrating how the
solution approaches the zero-mean-field limit with
$\delta\bm{B}\gg\overline{\bm{B}}$.

Solving Eq. 4 first requires an initial condition with constant $B^{2}$. This
is easily constructed by arbitrarily specifying the component of
$\delta\bm{B}$ in the
$\hat{\bm{n}}\equiv\hat{\bm{p}}\times\overline{\bm{B}}/|\hat{\bm{p}}\times\overline{\bm{B}}|$
direction—i.e., the $\delta\bm{B}$ direction for a linear Alfvén wave—then
solving for the $\delta\bm{B}$ component in the
$\hat{\bm{m}}\equiv\hat{\bm{p}}\times\hat{\bm{n}}/|\hat{\bm{p}}\times\hat{\bm{n}}|$
direction, noting also that $\hat{\bm{p}}\cdot\delta\bm{B}=0$ due to
periodicity in $\lambda$. This gives

$\hat{\bm{m}}\cdot\delta\bm{B}=-\hat{\bm{m}}\cdot\overline{\bm{B}}+\sqrt{B^{2}-(\hat{\bm{p}}\cdot\overline{\bm{B}})^{2}-(\hat{\bm{n}}\cdot\delta\bm{B})^{2}},$
(5)

with $\hat{\bm{n}}\cdot\delta\bm{B}$ an arbitrary function of $\lambda$. This
must be combined with the condition
$\overline{\hat{\bm{m}}\cdot\delta\bm{B}}=0$, which determines $B^{2}$. The
solution (5) illustrates two difficulties that also manifest in higher
dimensions: first, $B^{2}$ must be computed self-consistently for a particular
form of $\delta\bm{B}$; second, there may not exist real solutions for
arbitrary choices of $\hat{\bm{n}}\cdot\delta\bm{B}$ once $|\delta\bm{B}|$
approaches $|\overline{\bm{B}}|$. If one is not careful, either of these
difficulties will almost inevitably cause discontinuities in $\delta\bm{B}$ as
its solution is constructed. But, we reiterate that this does not signify that
spherically polarised solutions do not exist, or are non-smooth; rather, they
require constraining both components of $\delta\bm{B}$ simultaneously, as
opposed to specifying one component and solving for the other. As a concrete
example, Barnes & Hollweg (1974) show that for
$\hat{\bm{n}}\cdot\delta\bm{B}=\mathcal{A}_{n}\sin(k\lambda)$ solutions exist
only for $\mathcal{A}_{n}<\mathcal{A}_{n,{\rm crit}}=(\pi/2)\sin\vartheta$,
where
$\vartheta=\cos^{-1}(\hat{\bm{p}}\cdot\overline{\bm{B}}/|\overline{\bm{B}}|)$.
However, as we show below, if one starts with such a solution and evolves it
according to Eq. 4, there is no singular or otherwise interesting behaviour as
$\mathcal{A}$ passes through $\mathcal{A}_{n,{\rm crit}}$;
$\hat{\bm{n}}\cdot\delta\bm{B}$ simply changes shape to avoid the amplitude
limit. This example illustrates why “growing” the $\delta\bm{B}$ solution from
small amplitudes is necessary—without this, one is limited by the inability to
choose an appropriate free function that will enable the constant-$B$
constraint to be satisfied.

With a small-amplitude $\delta\bm{B}$ specified, Eq. 4 is easily solved by
standard numerical methods. Here we use sixth-order finite differences for
consistency with the 2-D and 3-D calculations, but a Fourier pseudospectral
method is equally suitable in 1-D (M+21; Squire et al., 2022). An example is
shown in Fig. 1 at several times as $\delta\bm{B}$ grows in amplitude. By the
final snapshot, with $\mathcal{A}\approx 400$, the solution is approaching the
zero-mean-field limit, which would give a purely stationary Alfvénic solution,
since the propagation velocity $\bm{v}_{\rm A}$ becomes much smaller than the
flows $\delta\bm{u}$ in the nonlinear solution. Note that once
$\delta\bm{B}\gg\overline{\bm{B}}$, $\phi$ becomes small and the shape of
$\delta\bm{B}$ remains nearly constant with growing amplitude.

##  4 Two dimensions

To explore general 2-D solutions, we stipulate that $\delta\bm{B}$ is a
function only of $\ell\equiv\hat{q}_{x}x+\hat{q}_{y}y$ and $z$, taking
$\overline{\bm{B}}=B_{0}\hat{\bm{x}}$. This geometry is the obvious
generalization of the 1-D wave described above, with $\delta\bm{B}$ varying
only in the plane angled at $\theta_{\rm
2D}\equiv\tan^{-1}(\hat{q}_{y}/\hat{q}_{x})$ to the mean field. As in 1-D, the
2-D calculation proceeds in two steps by first constructing a low-amplitude
near-linear solution, then growing this using Eqs. 2 and 3. Both steps are
significantly more complex than in 1-D.

Low-amplitude solution First, we note that the clear generalization of the
$\hat{\bm{n}}$-directed 1-D linear Alfvénic field to two (or three) dimensions
is the field $\delta\bm{B}=\nabla\times(A_{x}\hat{\bm{x}})$, with $A_{x}$
chosen arbitrarily. The goal, then, is to add additional field components to
enforce constant $B^{2}$, which is best done using the vector potential
($A_{y}$ and/or $A_{z}$) so as to maintain $\nabla\cdot\delta\bm{B}=0$ (the
vector potential is not needed in 1-D because enforcing
$\hat{\bm{p}}\cdot\delta\bm{B}=0$ ensures $\nabla\cdot\delta\bm{B}=0$). The
equation for $B^{2}$ becomes

$B^{2}={\rm
const.}=|B_{0}-\partial_{z}A_{y}+\partial_{y}A_{z}|^{2}+|\partial_{z}A_{x}-\partial_{x}A_{z}|^{2}+|\partial_{x}A_{y}-\partial_{y}A_{x}|^{2},$
(6)

where $\partial_{x}=\hat{q}_{x}\partial_{\ell}$ and
$\partial_{y}=\hat{q}_{y}\partial_{\ell}$. If one of $A_{y}$ or $A_{z}$ is
fixed, Eq. 6 is a first-order nonlinear partial differential equation (PDE) in
$A_{z}$ or $A_{y}$, which can in principle be solved using the method of
characteristics. However, such a method, which is effectively that used by
Valentini et al. (2019), is problematic on a periodic domain because the value
of the constant $B^{2}$ is not a-priori known, but itself depends on the
solution ($A_{z}$ or $A_{y}$). An incorrect choice of $B^{2}$ manifests as an
additional component of the mean field (as occurs for
$\hat{\bm{m}}\cdot\delta\bm{B}$ in Eq. 5), which will require $\bm{A}$ to
contain linear gradients, and thus lead to spurious discontinuities on a
periodic domain. To overcome this, we instead stipulate that $\nabla B^{2}=0$
in Eq. 6, which removes the issue of $B^{2}$ being undetermined at the cost of
increasing the order of the PDE. We also use the “mean-field Coulomb gauge”
$A_{y}=-\partial_{z}\alpha$, $A_{z}=\partial_{y}\alpha$ so that Eq. 6 involves
just one free function without causing an artificial difference between the
$\ell$ and $z$ directions. We then solve the resulting nonlinear PDE for
$\alpha(\ell,z)$ in Fourier space using MATLAB’s fsolve function with the
trust-region dogleg method. The two equations of $\nabla B^{2}=0$ are easily
combined into one by minimizing only nonzero Fourier modes with fsolve, and
the low-order approximate solution
$(\partial_{y}^{2}+\partial_{z}^{2})\alpha=-(|\partial_{z}A_{x}|^{2}+|\partial_{z}A_{x}|^{2})/(2B_{0})$
provides a good initial guess for the optimization. So long as $A_{x}$ is
chosen to be relatively smooth (e.g., the first one or two Fourier modes in
each direction), the procedure is rapid and straightforward because only a
small number of Fourier modes are needed to represent $\alpha$. We have found
no evidence for the development of discontinuities in such solutions, and a
solution for $\alpha$ can be found for most choices of smooth, low-amplitude
$A_{x}$ (at least in two dimensions; see below).

Figure 2: Two-dimensional spherically polarised solution on an $\ell,z$ plane
angled at $\theta_{\rm 2D}=30^{\circ}$ from the $\hat{\bm{x}}$ (mean-field)
direction. The top three panels show periodic traces of $B$ (black), $B_{x}$
(blue), $B_{y}$ (red), and $B_{z}$ (yellow) along the white line plotted on
the bottom-left panel (this lies at angle $\alpha\approx 11.3$ degrees from
$\hat{\bm{\ell}}$; integer $l$ values are marked to illustrate the
correspondence between the trace and the 2-D solution). The initial condition
is constructed from a random superposition of modes in $A_{x}$, scaled to give
amplitude $\mathcal{A}\approx 0.2$ (top panel). It then grows in time
according to Eqs. 2 and 3. The bottom panels show the 2-D structure of the
components of $\bm{B}$ at the time corresponding the bottom trace, when
$\mathcal{A}\approx 5$. At least to the precision of the $384^{2}$ resolution
used here, discontinuities develop in the field structure, unlike the 1-D
solutions (the most prominent is near near $l=1$ on the trace plots). Aside
from numerical ringing caused by the development of these discontinuities,
however, $B$ remains very constant throughout the domain (the colorbar of $B$
on the bottom left is scaled to $\pm 2\%$)

Large-amplitude solutions Starting from a constant-$B$ initial condition, Eqs.
2 and 3 can be solved using standard methods. For simplicity, we use an Euler
timestepper with the timestep chosen based on the maximum of $\phi$ over the
domain. The only complication arises from the non-homogenous derivative
operator in the Poisson equation (3), which must be recomputed and solved at
every timestep as $\delta\bm{B}$, and thus $\nabla_{\perp}$, change. We choose
to use a periodic sixth-order finite-difference representation111A range of
finite-difference orders were tested. Lower-order operators cause a clear
decrease in solution quality, but there was little gained for orders above
$6$. , forming $\nabla_{\perp}^{2}$ as a sparse matrix through left-
multiplication of the relevant gradient operators by $B_{i}$. A difficulty in
the interpretation and solution of Eq. 3 is that the $\nabla_{\perp}^{2}$
matrix has a rather high dimensional nullspace, which means that, depending on
the source $-\delta\bm{B}\cdot\overline{\bm{B}}$, Eq. 3 may not have a
solution222To understand why this is the case, it is helpful to consider a
simple example such as $\bm{B}=\hat{\bm{z}}$, which gives
$\nabla_{\perp}^{2}=\partial_{\ell}^{2}$. Taking arbitrary $f$ in
$\nabla_{\perp}^{2}\phi=f$, it is clear that any $\ell$-independent $f$ that
is periodic in $z$ cannot be captured by $\nabla_{\perp}^{2}\phi$ on a
periodic domain. Thus, the solution to Eq. 3 cannot capture any part of
$\delta\bm{B}\cdot\overline{\bm{B}}$ that varies along $\bm{B}$ but is
constant perpendicular to it (though, exactly what this statement means for a
complex 2-D or 3-D $\delta\bm{B}$ is not obvious). Note that this difficulty
is absent in 1-D because all qualities vary only in $\lambda$.. We circumvent
the issue by interpreting Eq. 3 in the least-squares sense, thus finding the
best approximation to the flow $\tilde{\bm{u}}=\nabla\phi$ that maintains
constant $B$. We use the least-squares conjugate-gradient method (Paige &
Saunders, 1982) with an incomplete LU preconditioner. As we show below, the
solution is generally very good, but the method certainly does not guarantee
constant $B$, at least at finite resolution, unlike in 1-D. Presumably, if
$-\delta\bm{B}\cdot\overline{\bm{B}}$ has a large component in the nullspace
of $\nabla_{\perp}^{2}$, this implies that $\delta\bm{B}$ has developed
structures that cannot continue to grow in amplitude while maintaining
constant $B$. Understanding the conditions under which this occurs requires
further study.

Figure 3: Same as Fig. 2 but starting from a different initial random
collection of Fourier modes in the low-amplitude $A_{x}$ initial conditions.
We show only the solution with $\mathcal{A}\approx 5$. In this case,
discontinuous structures do not develop and the solution is well resolved at
this resolution of $256^{2}$ (which is lower than Fig. 2)

An example solution is shown in Fig. 2. We use a periodic domain of size
$L_{\ell}=1$, $L_{z}=1$, with $384$ grid points in each direction and
$\theta_{\rm 2D}=30^{\circ}$. The initial conditions in $A_{x}$ are
constructed from a random combination of $|k_{\ell}|\leq 2\pi/L_{\ell}$ and
$|k_{z}|\leq 2\pi/L_{z}$ modes, scaled to give $\mathcal{A}\approx 0.2$. In
the top three panels, we show the field structure by taking an angled trace
through the domain, accounting for the periodicity in order to show most of
the solution (but note that some larger structures appear twice; the bottom-
left panel illustrates the correspondence to the 2-D domain). The smooth,
small-$\delta B_{x}$ initial conditions are shown in the top panel of Fig. 2
and have a very small variation in $B$, with
$({\overline{B^{2}}})^{1/2}/\overline{B}\approx 10^{-5}$. As $\delta\bm{B}$
grows, it develops quite sharp gradients by modest amplitudes
($\mathcal{A}\gtrsim 0.7$; see second panel), which are are numerically
unresolved at this resolution. The large-amplitude $\mathcal{A}\approx 5$
solution is shown both in bottom trace and in the image plots in the bottom
row, illustrating the discontinuities that have developed at several locations
in the domain. These sharp gradients represent a distinct difference compared
to 1-D solutions, which only ever steepen modestly from the initial waveform
as they grow to large amplitudes (see Fig. 1; this is discussed in more detail
below). The method clearly maintains extremely constant $B$, aside from
numerical ringing near some of the sharp gradient discontinuities that form
(excluding these regions $({\overline{B^{2}}})^{1/2}/\overline{B}\approx
0.8\%$).

We have also explored solutions on planes at other angles $\theta_{\rm 2D}$
(not shown) and with differing initial conditions. As we demonstrate in Fig.
3, which shows another example large-amplitude solution, the propensity to
form sharp gradients depends on the structure of the initial low-amplitude
$A_{x}$ (see below for further discussion). Unlike Fig. 2, the solution in
Fig. 3 is well resolved and relatively smooth, constituting a practical
demonstration that 2-D smooth, spherically polarised solutions exist (which,
as far as we are aware, was not previously known). We have not found any
simple heuristic to determine how the sharpness of the final solution relates
to the initial low-amplitude one. In general, solutions with larger
$\theta_{\rm 2D}$— _viz.,_ those which are more elongated along the magnetic
field— develop somewhat larger variation in $B$. This seems to be due to
larger inaccuracies in the solutions of Eq. 3, although whether this relates
to discontinuities remains unclear.

##  5 Three dimensions

Three-dimensional solutions are constructed through a method that is almost
identical to the 2-D case. The only additional complication is that the method
to construct the low-amplitude solution using the mean-field Coulomb gauge
($A_{y}=-\partial_{z}\alpha$, $A_{z}=\partial_{y}\alpha$) cannot be used to
construct a field that varies in $x$ but not in $y$ or $z$ (i.e., a field with
power in $k_{\perp}=\sqrt{k_{y}^{2}+k_{z}^{2}}=0$ modes). Although a different
method of constructing a low-amplitude solution may alleviate this, we opt
instead to adjust the chosen $A_{x}$ so that its associated
$B^{2}=|\partial_{y}A_{x}|^{2}+|\partial_{z}A_{x}|^{2}$ has very little power
in $k_{\perp}=0$ modes. This is easily done using MATLAB’s lsqnonlin function
after constructing $A_{x}$ from a collection of random Fourier modes. The
method seems to work well; for example, it generates smooth initial conditions
with $\mathcal{A}\approx 0.2$ and
$({\overline{B^{2}}})^{1/2}/\overline{B}\approx 3\times 10^{-5}$ (see Fig. 4).
In solving Eqs. 2 and 3, the only extra challenge compared to 2-D is the
computational expense, and we are limited to constructing solutions with
resolutions ${\lesssim}48^{3}$ due to our non-parallelised implementation
using MATLAB. However, the only significant computational difficulty—
inverting $\nabla_{\perp}^{2}$ in Eq. 3—involves standard iterative matrix-
solve methods, which have robust and efficient parallel implementations that
could allow for much higher resolutions if desired.

Figure 4: Three-dimensional spherically polarised solution in a cubic box with
a resolution of $48^{3}$. As in Fig. 2, the top three panels show periodic
traces of $B$ (black), $B_{x}$ (blue), $B_{y}$ (red), and $B_{z}$ (yellow)
along a line in the direction $(\cos\theta_{\rm 3D},\sin\theta_{\rm
3D}\,\cos\varphi_{\rm 3D},\sin\theta_{\rm 3D}\,\sin\varphi_{\rm 3D})$, with
$\theta_{\rm 3D}\approx 30^{\circ}$ and $\varphi_{\rm 3D}\approx
11.3^{\circ}$, with $l=0$ at the center of the box (the units are scaled to
the size of the box). The initial condition is constructed from random modes
in $A_{x}$ with an amplitude such that $\mathcal{A}\approx 0.2$ (top panel),
then growing in time according to Eqs. 2 and 3. The bottom panels show the 3-D
structure of the components of $\bm{B}$ at the time corresponding the bottom
trace, when $\mathcal{A}\approx 5$.

An example 3-D solution in a cubic box is shown in Fig. 4. As in Fig. 2, we
show the solution structure using a trace along an angled, periodic line, at
$\mathcal{A}\approx 0.2$ (the initial conditions), $\mathcal{A}\approx 0.7$,
and $\mathcal{A}\approx 5$. Like in 2-D, the solutions appear to become
discontinuous, although it is more difficult to diagnose in detail because of
the limited resolution. The solutions have somewhat larger variation in
$B^{2}$ compared to Fig. 2, with
$({\overline{B^{2}}})^{1/2}/\overline{B}\approx 3\%$ by $\mathcal{A}\approx
5$, although this is at least partially a result of the lower resolution. We
have also explored solutions in boxes that are elongated along
$\overline{\bm{B}}$ with $L_{x}=4$, finding similar structures and properties,
at least to the accuracy achievable here. We also note that we have confirmed
that $\nabla\cdot\bm{B}$ remains zero within the tolerances of the finite
difference representation (with sixth-order finite differences at $48^{3}$ the
root-mean-square of $\nabla\cdot\delta\bm{B}$ is ${\approx}1.5\times
10^{-6}$).

##  6 Discussion and conclusions

Above, we have presented a method to construct large-amplitude spherically
polarised Alfvénic structures, which form a broad class of nonlinear solutions
to the compressible MHD (or kinetic MHD) equations. The method works by
“growing” a low-amplitude, nearly linear ($\delta\bm{B}\ll\overline{\bm{B}}$)
wave, thus allowing the exploration of nonlinear solutions in the large-
amplitude $\delta\bm{B}\sim\overline{\bm{B}}$ regime applicable to the solar
wind, or even zero-mean-field limiting solutions with
$\delta\bm{B}\gg\overline{\bm{B}}$. We have presented some examples of such
solutions in one, two, and three dimensions, albeit with somewhat limited
resolution due to computational challenges. Others—including, in principle,
solutions with turbulent-like spectra333Reliably doing this would require a
more refined optimization method for generating the low-amplitude solution, as
described in § 4; our current implementation can work reliably only with a
relatively smooth choice for $A_{x}$. —could be generated using our method by
choosing a different form for the initial, small-amplitude wave. These
solutions exhibit some important features observed in the solar wind,
particularly the asymmetry (one-sidedness) of $\delta
B_{\|}=\delta\bm{B}\cdot\overline{\bm{B}}/|\overline{\bm{B}}|=\delta B_{x}$
fluctuations at modest $\mathcal{A}$, which is a consequence of maintaining
constant $B$ through large changes in $\delta\bm{B}$ (Gosling et al., 2009).
At larger $\mathcal{A}$, all solutions cause magnetic-field reversals
($|\delta B_{\|}|>|\overline{\bm{B}}|$, a commonly used definition of
switchbacks), although these effects would be stronger if we chose more
oblique solutions (e.g., larger $\theta_{\rm 2D}$ or a more elongated box in
3-D; Mallet et al., 2021). However, it is also worth noting that despite its
apparent success in nearly all cases we have explored, the method is not
guaranteed to produce perfectly constant-$B$ structures because the Poisson-
like equation (3), which determines how $\delta\bm{B}$ changes shape as it
grows, lacks exact solutions in general, at least at finite resolution. These
mathematical issues should be explored in more detail in future work. Other
possible ideas for future studies include using a time-dependent mean field or
gradient operators, which would change the large-amplitude $\delta\bm{B}$ that
results from chosen small-amplitude initial conditions (M+21; Squire et al.,
2022), or the extension to the relativistic regime (Mallet & Chandran, 2021),
which may have interesting applications to a number of high-energy processes
such as disk coronae (Chandran et al., 2018) or pulsar magnetospheres (Kumar &
Bošnjak, 2020; Zhang, 2020).

One interesting use case for these solutions is as input for nonlinear MHD or
kinetic simulations to study processes such as parametric decay (instability
of Alfvén waves; Sagdeev & Galeev, 1969) or large-amplitude reflection-driven
turbulence (Johnston et al., 2022). For example, a 3-D MHD simulation could be
initialised with a 1-D, 2-D, or 3-D solution $\delta\bm{B}$ from our method by
first extending it via symmetry to three dimensions (e.g., $\delta\bm{B}_{\rm
3D}(x,y,z)=\delta\bm{B}(\hat{\bm{p}}\cdot\bm{x})$ for a 1-D solution), then
setting $\delta\bm{u}=\pm\delta\bm{B}/\sqrt{4\pi\rho}$. Study of parametric
decay especially has been somewhat limited by a lack of general, large-
amplitude Alfvénic solutions other than circularly polarized waves that can be
evolved numerically, so our method may enable important progress in this area
(the modest-amplitude case with $\delta\bm{B}\lesssim\overline{\bm{B}}$ has
been studied in 1-D by Del Zanna 2001; Del Zanna et al. 2015; Tenerani et al.
2020 considered a larger-amplitude 2-D solution). A particularly interesting
regime, which could in principle be studied for solutions of any
dimensionality, would be the zero-mean-field limit, where the nonlinear
Alfvénic solutions become a stationary, non-propagating tangle of magnetic
field lines and Alfvénic flows. The system is reminiscent of the tangled
fields possible in magnetostatic force-free equilibria (those with $\bm{u}=0$
and $\bm{B}\times(\nabla\times\bm{B})=0$; e.g., Chandrasekhar & Woltjer 1958;
Marsh 1996; Hosking et al. 2020), but with a different class of equilibria
that includes flows and has seen comparatively little study.

Figure 5: Measurement of discontinuity formation in 1-D solutions (panel a)
and 2-D solutions starting from two different initial conditions (panels b and
c). We plot the normalised infinity norm of the gradient of the solutions,
$||\nabla\bm{B}||^{\infty}/||\bm{B}||^{\infty}=(1/3)\sum_{i}\sum_{j}\max(\nabla_{j}B_{i})/\max(B_{i})$
(with the $\nabla_{j}$ taken along either $\lambda$ or $\ell$ and $z$), as a
function of $\mathcal{A}$ for a scan in resolution in each case (an $N\times
N$ grid is used in 2-D; we list only every second $N$ in the legend for
clarity). In 1-D (left panel), we initialise with a linear combination of
modes with $k\leq 4\pi/L$; in 2-D, we initialise with a linear combination of
modes with $k_{\ell}\leq 2\pi/L_{\ell}$ and $k_{z}\leq 2\pi/L_{z}$ (the
$N=384$ case of panel b is that shown in Fig. 2). Dotted lines in each case
show a scaling $N^{\chi}$, with $\chi$ chosen to match the scaling of
unconverged solutions ($\chi\approx 0.7$ in 1-D and $\chi\approx 0.8$ in 2-D).
Clearly, the 1-D solution converges at very low resolution ($N\approx 64$),
showing that Eq. 2 does not lead to particularly small-scale features. In
contrast, in 2-D, sharp field structures form much more readily: in the first
example in panel (b), which is that from Fig. 2, there is no convergence even
at the highest resolution that is feasible using our current computational
implementation ($N=384$); but, the second example in panel (c), which is that
from Fig. 3 and simply starts from a different random initial condition,
achieves convergences around $N=128$.

A second reason for interest in these solutions relates directly to their
properties and structure. As far as we are aware, it was previously not known
whether smooth, spherically polarised solutions with two- or three-dimensional
structure even existed; our method has provided near-perfect examples in 2-D,
one of which is shown in Fig. 3. In 3-D, we are resolution limited due to our
non-parallelised numerical implementation (the example in Fig. 4 is not smooth
at this resolution), but it seems unlikely that it would differ fundamentally
from 2-D. Interestingly, in some other solutions, an example of which is shown
in Fig. 2, $\bm{B}$ develops extremely sharp gradients that remain
discontinuous at our highest resolution of $384^{2}$. We demonstrate this in
Fig. 5, which quantifies the development of sharp structures by plotting how
the maximum gradient of $\bm{B}$ evolves with solution amplitude, over a scan
in resolution. For reference, we plot with dotted lines the empirically
determined power-law behaviour of an unconverged solution
($||\nabla\bm{B}||^{\infty}/||\bm{B}||^{\infty}\propto N^{\chi}$; the exponent
$\chi$ differs between 1-D and 2-D). Convergence occurs at low resolution in
1-D solutions (Fig. 5a), which only ever develop modestly sharper structures
compared to the initial wave, and are converged by $N\approx 48$ for this
example ($b_{n}$ initialised with $k\leq 4\pi/L$ modes). In 2-D, both
solutions become significantly sharper, but panel (b) (that from Fig. 2) shows
no sign of convergence at all by $384^{2}$, while panel (c) (the case in Fig.
3) converges at around $128^{2}$. Given there does not seem to be any simple
distinguishing feature(s) that determine whether discontinuities develop444The
solution in Fig. 3 could be considered “more 1-D” than that in Fig. 2, in that
in varies predominantly along a diagonal line across the domain. However, some
other cases explored do not share this property and still do not develop sharp
gradients, while other “nearly 1-D” solutions do seem to become extremely
sharp. , the most likely scenario seems to be that continuous solutions do
generally exist, but they require extremely sharp gradients in some cases. If
true, this suggests 1-D waves are simply a particular special case of more
general 2-D or 3-D fields that happens to allow large-amplitude solutions with
an especially smooth structure.

The conclusions of the previous paragraph and Fig. 5 may have interesting
consequences for switchback formation in the solar wind. Specifically, our
method based on Eqs. 2 and 3 shares clear similarities with the _physical_
processes that occur in the solar wind, where Alfvénic structures grow in
normalised amplitude due to plasma expansion. Although the direct effect of
expansion is not included in Eqs. 2 and 3 (this would cause time dependence of
the gradient operators), the general idea—whereby, in order to maintain
constant $B^{2}$ as it grows, $\delta\bm{B}$ changes shape by means of a
compressive flow—is very similar to physical expansion, and would presumably
produce similar $\delta\bm{B}$ evolution555This ignores the effect of wave
reflection in creating backwards propagating waves that seed turbulence, which
creates smaller-scale structures in a very different way. In ignoring such
effects, we assume the system remains strongly dominated by outwards-
propagating fluctuations and thus can remain nearly Alfvénic.. In support of
this idea, in 1-D with expansion effects added, Eqs. 2 and 3 become exactly
the system of M+21, which was derived asymptotically from the expanding MHD
equations and produces similar large-amplitude structures to Eqs. 2 and 3. If
this correspondence holds in 3-D, our results imply that small-amplitude
Alfvénic perturbations released from the low corona would often develop very
sharp gradients or discontinuities in the process of growing, unless they are
in some sense one dimensional (unlikely if they are created by turbulence in
the chromosphere; van Ballegooijen et al., 2011). This is quite promising for
the _in-situ_ Alfvénic scenario of switchback formation (Squire et al., 2020;
Shoda et al., 2021), in which switchbacks are simply Alfvénic fluctuations
that have grown to large amplitudes by expansion. While M+21 suggested, based
on 1-D solutions, that the model may struggle to explain the extremely sharp
switchback structures seen in some observations (Kasper et al., 2019; Bale et
al., 2019; Akhavan-Tafti et al., 2021), the answer may simply be that near-
discontinuities naturally develop when starting from low-amplitude
fluctuations with more general (non 1-D) structure. Whether this idea has
direct observable consequences—for example, in the relative populations of
rotational versus tangential discontinuities (e.g., Neugebauer et al.,
1984)—requires better understanding of why and how sharp gradients develop, so
is left to future work.

###### Acknowledgements.

The authors acknowledge useful discussion with J. Burby, R. Meyrand, B.
Chandran, Z. Johnston, and J. Nättilä in relation to this work. Support for
J.S. was provided by Rutherford Discovery Fellowship RDF-U001804, which is
managed through the Royal Society Te Apārangi. A.M. acknowledges the support
of NASA through grant 80NSSC21K0462.

Declaration of Interests: The authors report no conflict of interest.

## References

  * Akhavan-Tafti et al. (2021) Akhavan-Tafti, M., Kasper, J., Huang, J. & Bale, S. 2021 Discontinuity analysis of the leading switchback transition regions. Astron. Astro. 650, A4.
  * Bale et al. (2019) Bale, S. D., Badman, S. T., Bonnell, J. W., Bowen, T. A., Burgess, D., Case, A. W., Cattell, C. A., Chandran, B. D. G., Chaston, C. C., Chen, C. H. K., Drake, J. F., de Wit, T. D. et al. 2019 Highly structured slow solar wind emerging from an equatorial coronal hole. Nature 576 (7786), 237–242.
  * Bale et al. (2021) Bale, S. D., Horbury, T. S., Velli, M., Desai, M. I., Halekas, J. S., McManus, M. D., Panasenco, O., Badman, S. T., Bowen, T. A., Chandran, B. D. G., Drake, J. F., Kasper, J. C., Laker, R., Mallet, A., Matteini, L., Phan, T. D., Raouafi, N. E., Squire, J., Woodham, L. D. & Woolley, T. 2021 A Solar Source of Alfvénic Magnetic Field Switchbacks: In Situ Remnants of Magnetic Funnels on Supergranulation Scales. Astrophys. J. 923 (2), 174.
  * Barnes & Hollweg (1974) Barnes, A. & Hollweg, J. V. 1974 Large-amplitude hydromagnetic waves. J. Geophys. Res. 79 (16), 2302.
  * Barnes & Suffolk (1971) Barnes, A. & Suffolk, G. C. J. 1971 Relativistic kinetic theory of the large-amplitude transverse Alfvén wave. J. Plasma Phys. 5 (3), 315–329.
  * Belcher & Davis (1971) Belcher, J. W. & Davis, Leverett, J. 1971 Large-amplitude Alfvén waves in the interplanetary medium, 2. J. Geophys. Res. 76 (16), 3534.
  * Bowen et al. (2021) Bowen, T. A., Badman, S. T., Bale, S. D., Dudok de Wit, T., Horbury, T. S., Klein, K. G., Larson, D., Mallet, A., Matteini, L., McManus, M. D. & Squire, J. 2021 Nonlinear Interactions in Spherically Polarized Alfvénic Turbulence. arXiv e-prints p. 2110.11454.
  * Chandran et al. (2018) Chandran, B. D. G., Foucart, F. & Tchekhovskoy, A. 2018 Heating of accretion-disk coronae and jets by general relativistic magnetohydrodynamic turbulence. J. Plasma Phys. 84 (3), 905840310\.
  * Chandrasekhar & Woltjer (1958) Chandrasekhar, S. & Woltjer, L. 1958 On force-free magnetic fields. Proc. Nat. Acad. Sci. 44 (4), 285–289.
  * Cohen & Dewar (1974) Cohen, R. H. & Dewar, R. L. 1974 On the backscatter instability of solar wind Alfvén waves. J. Geophys. Res. 79 (28), 4174.
  * Del Zanna (2001) Del Zanna, L. 2001 Parametric decay of oblique arc-polarized Alfvén waves. Geophys. Res. Lett. 28 (13), 2585–2588.
  * Del Zanna et al. (2015) Del Zanna, L., Matteini, L., Landi, S., Verdini, A. & Velli, M. 2015 Parametric decay of parallel and oblique Alfvén waves in the expanding solar wind. J. Plasma Phys. 81 (1), 325810102.
  * Gosling et al. (2009) Gosling, J. T., McComas, D. J., Roberts, D. A. & Skoug, R. M. 2009 A One-Sided Aspect of Alfvenic Fluctuations in the Solar Wind. Astrophys. J. Lett. 695 (2), L213–L216.
  * Hollweg (1974) Hollweg, J. V. 1974 Transverse Alfvén waves in the solar wind: Arbitrary $k$, $v_{0}$, $B_{0}$, and $|\delta B|$. J. Geophys. Res. 79 (10), 1539.
  * Hosking et al. (2020) Hosking, D. N., Schekochihin, A. A. & Balbus, S. A. 2020 Elasticity of tangled magnetic fields. J. Plasma Phys. 86 (5), 905860511.
  * Johnston et al. (2022) Johnston, Z., Squire, J., Mallet, A. & Meyrand, R. 2022 On the properties of Alfvénic switchbacks in the expanding solar wind: Three-dimensional numerical simulations. Phys. Plasmas 29 (7), 072902\.
  * Kasper et al. (2019) Kasper, J. C., Bale, S. D., Belcher, J. W., Berthomier, M., Case, A. W., Chandran, B. D. G., Curtis, D. W., Gallagher, D., Gary, S. P., Golub, L., Halekas, J. S. et al. 2019 Alfvénic velocity spikes and rotational flows in the near-Sun solar wind. Nature 576, 228–231.
  * Kulsrud (1983) Kulsrud, R. M. 1983 MHD description of plasma. In Handbook of Plasma Physics (ed. R N Sagdeev & M N Rosenbluth). Princeton University.
  * Kumar & Bošnjak (2020) Kumar, P. & Bošnjak, Ž. 2020 FRB coherent emission from decay of Alfvén waves. Mon. Not. R. Astron. Soc. 494 (2), 2385–2395.
  * Mallet & Chandran (2021) Mallet, A. & Chandran, B. D. G. 2021 Exact nonlinear solutions for three-dimensional Alfvén-wave packets in relativistic magnetohydrodynamics. J. Plasma Phys. 87 (6), 175870601.
  * Mallet et al. (2021) Mallet, A., Squire, J., Chandran, B. D. G., Bowen, T. & Bale, S. D. 2021 Evolution of Large-amplitude Alfvén Waves and Generation of Switchbacks in the Expanding Solar Wind. Astrophys. J. 918 (2), 62.
  * Marsh (1996) Marsh, G. E. 1996 Force-free magnetic fields: solutions, topology and applications. World Scientific.
  * Matteini et al. (2015) Matteini, L., Horbury, T. S., Pantellini, F., Velli, M. & Schwartz, S. J. 2015 Ion Kinetic Energy Conservation and Magnetic Field Strength Constancy in Multi-fluid Solar Wind Alfvénic Turbulence. Astrophys. J. 802 (1), 11.
  * Matteini et al. (2018) Matteini, L., Stansby, D., Horbury, T. S. & Chen, C. H. K. 2018 On the 1/f Spectrum in the Solar Wind and Its Connection with Magnetic Compressibility. Astrophys. J. Lett. 869 (2), L32.
  * Neugebauer et al. (1984) Neugebauer, M., Clay, D. R., Goldstein, B. E., Tsurutani, B. T. & Zwickl, R. D. 1984 A reexamination of rotational and tangential discontinuities in the solar wind. J. Geophys. Res. 89 (A7), 5395–5408.
  * Paige & Saunders (1982) Paige, C. C. & Saunders, M. A. 1982 Lsqr: An algorithm for sparse linear equations and sparse least squares. ACM Trans. Math. Software pp. 43–71.
  * Primavera et al. (2019) Primavera, L., Malara, F., Servidio, S., Nigro, G. & Veltri, P. 2019 Parametric Instability in Two-dimensional Alfvénic Turbulence. Astrophys. J. 880 (2), 156.
  * Roberts (2012) Roberts, D. A. 2012 Construction of Solar-Wind-Like Magnetic Fields. Phys. Rev. Lett. 109 (23), 231102.
  * Sagdeev & Galeev (1969) Sagdeev, R. Z. & Galeev, A. A. 1969 Nonlinear Plasma Theory. New York: Benjamin.
  * Schekochihin et al. (2009) Schekochihin, A. A., Cowley, S. C., Dorland, W., Hammett, G. W., Howes, G. G., Quataert, E. & Tatsuno, T. 2009 Astrophysical gyrokinetics: Kinetic and fluid turbulent cascades in magnetized weakly collisional plasmas. Astrophys. J. Supp. 182 (1), 310.
  * Shoda et al. (2021) Shoda, M., Chandran, B. D. G. & Cranmer, S. R. 2021 Turbulent Generation of Magnetic Switchbacks in the Alfvénic Solar Wind. Astrophys. J. 915 (1), 52.
  * Squire et al. (2020) Squire, J., Chandran, B. D. G. & Meyrand, R. 2020 In-situ Switchback Formation in the Expanding Solar Wind. Astrophys. J. Lett. 891 (1), L2.
  * Squire et al. (2022) Squire, J., Meyrand, R., Kunz, M. W., Arzamasskiy, L., Schekochihin, A. A. & Quataert, E. 2022 High-frequency heating of the solar wind triggered by low-frequency turbulence. Nature Astronomy 6, 715–723.
  * Tenerani & Velli (2017) Tenerani, A. & Velli, M. 2017 Evolving Waves and Turbulence in the Outer Corona and Inner Heliosphere: The Accelerating Expanding Box. Astrophys. J. 843 (1), 26.
  * Tenerani et al. (2020) Tenerani, A., Velli, M., Matteini, L., Réville, V., Shi, C., Bale, S. D., Kasper, J. C., Bonnell, J. W., Case, A. W., de Wit, T. D., Goetz, K., Harvey, P. R., Klein, K. G., Korreck, K., Larson, D., Livi, R., MacDowall, R. J., Malaspina, D. M., Pulupa, M., Stevens, M. & Whittlesey, P. 2020 Magnetic field kinks and folds in the solar wind. Astrophys. J. Supp. 246 (2), 32.
  * Valentini et al. (2019) Valentini, F., Malara, F., Sorriso-Valvo, L., Bruno, R. & Primavera, L. 2019 Building Up Solar-wind-like 3D Uniform-intensity Magnetic Fields. Astrophys. J. Lett. 881 (1), L5.
  * van Ballegooijen et al. (2011) van Ballegooijen, A. A., Asgari-Targhi, M., Cranmer, S. R. & DeLuca, E. E. 2011 Heating of the Solar Chromosphere and Corona by Alfvén Wave Turbulence. Astrophys. J. 736 (1), 3.
  * Zhang (2020) Zhang, B. 2020 The physical mechanisms of fast radio bursts. Nature 587 (7832), 45–53.


# The Importance of Data Preprocessing in Malware Classification
Explainability

Elyes Manai  Department of Computer Science and Software Engineering
Laval University
Quebec, Canada
<EMAIL_ADDRESS>Mohamed Mejri  Department of Computer Science and
Software Engineering
Laval University
Quebec, Canada
<EMAIL_ADDRESS>Jaouhar Fattahi  Department of Computer Science
and Software Engineering
Laval University
Quebec, Canada
<EMAIL_ADDRESS>
###### Abstract

This paper investigates the impact of preprocessing techniques on the
explainability of XAI (Explainable Artificial Intelligence) algorithms. Using
a malware classification dataset, we trained an XGBoost model and compared the
performance of two feature encoding methods: Label Encoding (LE) and One Hot
Encoding (OHE). Our findings reveal a marginal performance loss when using OHE
instead of LE. However, the more detailed explanations provided by OHE
compensated for this loss. We observed that OHE enables deeper exploration of
details in both global and local contexts, facilitating more comprehensive
answers. Additionally, we observed that using OHE resulted in smaller
explanation files and reduced analysis time for human analysts. These findings
emphasize the significance of considering feature encoding techniques in XAI
research and suggest potential for further exploration by incorporating
additional encoding methods and innovative visualization approaches.

###### Index Terms:

Explainability, XAI, preprocessing, feature encoding

##  I Introduction

Machine learning has witnessed remarkable advancements in recent years,
enabling the development of sophisticated models that achieve impressive
performance on various tasks. As these tasks and the data they are trained on
become more complex, so does the model complexity. This often causes the
decision-making process to lack transparency, making it difficult to
understand the reasons behind their predictions. In a society that uses AI for
an ever-growing number of use cases, however, that lack of understanding can
pose serious risks to the users. Averting these risks and allowing more
control over what our AI is doing, thus allowing more responsible AIs, is the
goal behind the Explainable Artificial Intelligence (XAI) subfield. This
subdomain of AI focuses on making black-box models transparent by providing
understandable explanations for their decisions. XAI also allows us to combine
the powerful pattern-recognition learning capabilities of AI with human-
readable explanations that humans can instinctively understand and explain.
the algorithms used in XAI usually work by finding out what parts of the input
and of the model weights most affect the model’s predictions. The end result
will be a summary of each feature’s contribution to the model. How helpful are
these summaries, however, which we can call the quality of the generated
explanations, depends on several parameters such as the chosen algorithm, the
model architecture, and the data preprocessing technique. This last parameter,
however, is not as popular as the others. While most XAI research focuses on
algorithms, use cases, and the quality of explanations generated, there is a
lack of research on the impact of preprocessing on generated explanations. We
think that the preprocessing technique has a sizable impact on the quality of
generated explanations and should be more explored. More specifically, we are
interested in the feature encoding step of the preprocessing pipeline. Since
XAI methods summarize feature contribution, the way we encode our models will
directly affect the understandability of the generated explanations. Since
preprocessing directly affects model performance, considerations must be taken
to not trade off too much performance for better explanations, as better
explanations on an unprecise model are not useful. Nonetheless, we think that
a minor performance loss for a major boost in explainability is worth it, as
it also opens up the door for better model and data understanding, bias
discovery, robustness tests, and overall higher quality assurance. This is
especially important in critical industries such as Medicine, Finance, and
Cyber Security. To showcase the added value of our idea in a real use case, we
will apply Machine Learning and Explainability on a common problem in Cyber
security: Malware Classification. It is one of the most common tasks that
Machine Learning is applied to in modern antiviruses and Intrusion Detection
Systems. We will train a model on a publicly available malware dataset, apply
the XAI algorithm, switch the preprocessing technique and compare the
generated explanations. We will show that new rules and pain points can be
detected and further explored by just changing the preprocessing technique. To
the best of our knowledge, no prior studies specifically addressed the subject
of the direct impact of preprocessing on explanation quality in the field of
XAI have been identified in the existing literature. Our comprehensive review
of the literature revealed that research in XAI is more geared towards XAI
algorithms [4, 3, 5, 6], the generated explanations[7, 8], alternative ways to
bake explainability into the input features[9, 10] and other related
problems[11, 12]. Our focus in this paper can be summarized as follows: Given
that XAI algorithms use the input features as the key components for the
generated explanations, it is safe to assume that the type of feature encoding
used will directly affect the clarity of the explanations. The more explicit
the feature, the more detailed should be the explanation we get. With that in
mind, we will study two main questions in this paper:

  1. 1.

Does feature encoding affect explainability?

  2. 2.

If yes, what encoding yields better explainability and why?

##  II Concepts

###  II-A Machine Learning Modeling

Machine learning modeling refers to the process of creating a predictive
mathematical model using machine learning algorithms. It involves training a
model on a labeled dataset, evaluating its performance, and using it to make
predictions or gain insights from new, unseen data. Here is a brief overview
of the steps involved in machine learning modeling:

  1. 1.

Data Preparation: collecting and preparing the dataset for modeling. It
includes tasks such as data cleaning, handling missing values, dealing with
outliers, and feature engineering (creating new features from existing ones).

  2. 2.

Splitting the Dataset: splitting it into two or more subsets: a training set,
a validation set, and a test set. The training set is used to train the model,
the validation set is used for hyperparameter tuning and model selection, and
the test set is used to evaluate the final model’s performance.

  3. 3.

Model Selection: Choosing the most adequate model(s) among various machine
learning algorithms such as linear regression, decision trees, random forests,
support vector machines, neural networks, etc. The choice of model depends on
the type of problem, the nature of the data, and the desired outcome.

  4. 4.

Model Training: training the model on the training dataset by fitting it to
the input features and the corresponding target variables. During the training
process, the model learns the underlying patterns and relationships in the
data.

  5. 5.

Hyperparameter Tuning: selecting the optimal values for parameters to improve
the model’s performance. Techniques like grid search, random search, or
Bayesian optimization can be used for hyperparameter tuning.

  6. 6.

Model Evaluation: Assessing the performance of the trained model through
evaluation metrics such as accuracy, precision, recall, F1 score etc. The
model is evaluated on the validation set to understand its generalization
capabilities and fine-tune any parameters if needed.

Machine learning modeling is an iterative process, and it often involves
experimenting with different algorithms, feature engineering techniques, and
hyperparameter settings to find the best-performing model for a given problem.

###  II-B Feature Encoding

Feature encoding, also known as feature transformation or feature
representation, is a crucial step in data preprocessing where categorical or
textual features are converted into numerical representations that can be
effectively used by machine learning algorithms. This is a mandatory step as
ML algorithms only deal with numerical features. The choice of encoding
technique directly impacts the ML performance. Here are some common feature
encoding techniques:

  * •

One-Hot Encoding: Each category within a categorical feature is represented by
a binary feature. If a feature has n categories, it is encoded into n binary
features, where only one feature is active (1) for a particular category, and
the rest are inactive (0). One-hot encoding is useful when there is no
inherent order or relationship among the categories.

  * •

Label Encoding: Label encoding assigns a unique numerical label to each
category within a categorical feature. Each category is represented by a
distinct integer value. Label encoding is suitable when the categories have an
ordinal relationship or when using algorithms that can directly work with
integer inputs.

  * •

Ordinal Encoding: Similar to label encoding, ordinal encoding assigns
numerical labels to categories. However, ordinal encoding takes into account
the order or rank of the categories and assigns values accordingly. For
example, ”low,” ”medium,” and ”high” could be encoded as 1, 2, and 3,
respectively.

  * •

Binary Encoding: Binary encoding represents categories as binary bit patterns.
Each category is assigned a unique binary code, and each bit in the code
represents the presence or absence of a category. Binary encoding can be
efficient for high-cardinality categorical features and reduces the
dimensionality compared to one-hot encoding.

  * •

Embedding: Embedding techniques are commonly used for encoding textual or
high-dimensional categorical features. Embeddings are dense, low-dimensional
representations that capture semantic relationships between categories.
Embeddings are learned using techniques like Word2Vec
[DBLP:conf/somet/FattahiZM21, DBLP:conf/somet/FattahiMZG22] or categorical
embedding layers in deep learning models [DBLP:journals/access/DahoudaJ21].

###  II-C Explainability

Explainability in the context of machine learning[10097489, 9391727,
ZiniA239667] refers to the ability to understand and interpret the decisions
or predictions made by a machine learning model. It involves gaining insights
into how and why a model arrives at a particular output, providing
transparency and comprehensibility to the decision-making process. There are
various approaches to achieving explainability:

  * •

Model-Agnostic Approaches: These methods aim to explain any black-box machine
learning model without relying on its internal structure. They involve
techniques like feature importance analysis, partial dependence plots[2], and
surrogate models, which provide insights into the relationship between input
features and model predictions.

  * •

Rule-Based Approaches: These approaches aim to generate human-readable rules
that describe the decision-making process of the model. Rule-based models,
such as decision trees or rule lists, can provide explicit if-then statements
that explain how specific features influence predictions.

  * •

Interpretable Model Architectures: Some machine learning models, such as
linear regression, logistic regression, or decision trees, inherently provide
interpretable explanations. Their simplicity and transparency allow users to
understand the impact of each feature on the final prediction.

  * •

Local Explanations: Local explanation methods focus on explaining individual
predictions rather than the model as a whole. Techniques like LIME[3] (Local
Interpretable Model-Agnostic Explanations) or SHAP[4] (SHapley Additive
exPlanations) provide insights into which features contributed the most to a
particular prediction.

  * •

Visualizations: Visualizations play a significant role in explaining complex
models and high-dimensional data. Techniques like heatmaps, bar plots, scatter
plots, or saliency maps help in visualizing feature importance, decision
boundaries, or highlighting influential regions in the data.

###  II-D Malware Detection

To demonstrate our work, we will take the common task of detecting malware.
Malware are malicious pieces of software that are designed to infiltrate and
damage information systems without the users’ consent [AslanS20, info11030168,
Ameyed73, ZiadiaFMP20, Lou2020]. The term malware covers a lot of categories
such as viruses, ransomware, worms, trojans, backdoors, spyware, keyloggers,
adware, bots, and rootkits. Malware analysts have to discover exactly what
happened to a system and make sure that the machines damaged by malicious
software are isolated from the organization’s network. The analysis done to
single out the suspicious parts of the software can sometimes take a group of
analysts and several hours or even days. Since undetected malware can have
devastating consequences on any organization, malware detection has been
deemed one of the most important tasks in cybersecurity. Several types of
systems have been built to detect and capture malware such as Intrusion
detection systems, antiviruses and firewalls, and these systems keep getting
smarter thanks to the combined shared knowledge of the cyber security
community and the rapid advancement of technology. Current Malware detection
systems use Machine Learning and Deep Learning to detect anomalies in files
and network packets to protect the systems they’re installed on. Since Machine
learning has been known for its fantastic classification capabilities, more
and more complex architectures and models are being tested and deployed to the
current market.

##  III Related Work

To the best of our knowledge, no prior studies specifically addressed the
subject of the direct impact of preprocessing on explanation quality in the
field of XAI have been identified in the existing literature. Our
comprehensive review of the literature revealed that research in XAI is more
geared towards XAI algorithms [4, 3, 5, 6], the generated explanations[7, 8],
alternative ways to bake explainability into the input features[9, 10] and
other related problems[11, 12]. We are likely to be the first ones working on
this and hope to pave the path for more research based on this subject.

##  IV Implementation and experimental results

###  IV-A Dataset

For this project, we found a Malware classification dataset from the 2015
Microsoft Malware Classification Challenge[13]. The public variant we managed
to download contains 19611 rows and 78 features. Each row represents a single
file. The dataset is imbalanced as there are 14599 malware files and 5012 non-
malware files, so 3 times as much malware. The dataset has no missing data and
all features are numerical aside from the ”Name” one.

###  IV-B Preprocessing

The ”Name” feature has been modified by the competition organizers to include
”virus” if the file is malware and thus be removed since it does not represent
real-life data. We do not apply any other preprocessing on the data aside from
feature encoding. In this work, we apply two encoding techniques to all the
features:

  * •

Label Encoding: Each feature value is represented by a unique integer.

  * •

One Hot Encoding: Each feature value becomes a separate binary column where 1
means the file’s value of that feature is the column name, and 0 if not. This
allows for more precise knowledge of what went wrong.

###  IV-C Machine Learning Modeling

For training, we choose XGBoost[DBLP:conf/somet/FattahiMZ21,
DBLP:journals/cbm/LiuGG23] as our base model and train it using its default
parameters, namely 100 estimators, a max depth of 5 and a learning rate of
0.1. We use the free Google Colab coding environment which offers a single
sever with 12.7GB of RAM and a single NVIDIA T4 GPU with 15GB of GPU RAM. To
evaluate our model, We use four popular metrics: Accuracy, Precision, Recall
and F1. In a nutshell, accuracy measures the overall correctness of the
model’s predictions by calculating the proportion of correctly classified
instances out of the total number of instances. Precision quantifies the
proportion of true positive predictions out of all positive predictions made
by the model, indicating the model’s ability to correctly identify positive
instances and minimize false positives. Recall measures the proportion of true
positive predictions out of all actual positive instances in the dataset,
representing the model’s ability to capture positive instances and minimize
false negatives. Finally, the F1 score combines precision and recall into a
single metric by taking their harmonic mean, providing a balanced assessment
of the model’s accuracy and considering both false positives and false
negatives. We showcase the performance results of XGBoost on the label encoded
dataset in Table I.

TABLE I: XGB results on both encoding techniques

Encoding Technique | F1 | Accuracy | Precision | Recall
---|---|---|---|---
Label Encoding | 0.991 | 0.993 | 0.992 | 0.998
One Hot Encoding | 0.988 | 0.989 | 0.987 | 0.998

Although we did not preprocess our data, aside from encoding them differently,
we managed to get pretty good results. We can therefore directly go to the
explainability part.

###  IV-D Explainability

For starters, we are going to take away non useful features because one hot
encoding all 77 features created 85102 features, which kept crashing our
environment due to insufficient RAM. To do that, we will use XGBoost’s built
in feature importance function to list each feature’s impact on the model’s
decision making. In Table II, we extract the top 10 influential features and
sort them from most to least important.

TABLE II: Feature Importance according to XGBoost

Rank | Feature | Importance score
---|---|---
1 | MajorSubsystemVersion | 0.6215
2 | Subsystem | 0.1362
3 | MinorOperatingSystemVersion | 0.0575
4 | MajorLinkerVersion | 0.0454
5 | SizeOfStackReserve | 0.0198
6 | Characteristics | 0.0100
7 | SectionMaxChar | 0.0084
8 | ImageBase | 0.0081
9 | SizeOfHeapReserve | 0.0081
10 | TimeDateStamp | 0.0059

According to Table II, the combined score of the 10 most important features
are 0.9381 which means that they represent 93.81% of the model’s decision
making power. We, therefore, can just keep these 10 features and not use the
rest. Doing so, we get the results shown in Table III.

TABLE III: XGBoost results on the top 10 features

Encoding Technique | F1 | Accuracy | Precision | Recall
---|---|---|---|---
Label Encoding | 0.992 | 0.992 | 0.991 | 0.998
One Hot Encoding | 0.985 | 0.985 | 0.983 | 0.998

Comparing the results shown in Table III to those in Table I show that
although we did lose a bit of performance, the drop is marginal (less than
1%). This means that if the One Hot encoding does provide us with more
explainability power, it would be recommended to use. For the next par, we
will use a dedicated Explainability Algorithm called Shapley Additive
Explanations (SHAP) to dig deeper into the model’s inner reasoning.

####  IV-D1 The SHAP algorithm

SHAP[4, SHAPAleneziL21] was introduced in 2017 and provides a unified way of
explaining the contribution of each input feature to the final prediction of
the model, based on calculated values called Shapley values. A Shapley value
is a measure of the marginal contribution of a feature to the prediction,
averaged over all possible combinations of features in the dataset. To
calculate the Shapley values for a particular prediction, SHAP applies a game-
theoretic approach based on the concept of cooperative games. It considers
each feature value as a ”player” in the game and computes the contribution of
each player to the final prediction. It then calculates the average
contribution of each player across all possible coalitions of players,
weighting each coalition by its probability of occurrence. This approach
results in a set of Shapley values, which represent the relative importance of
each feature to the prediction for a specific instance. These Shapley values
can be used to generate an explanation for the prediction, showing which
features had the greatest impact and how they affected the final outcome. The
mathematical formula used by SHAP to generate the Shapley Values is presented
in Figure 1.

$\phi_{i}(x)=\sum_{S\subseteq
N\setminus\\{i\\}}\frac{{|S|!(|N|-|S|-1)!}}{{|N|!}}[f(x_{S}\cup\\{x_{i}\\})-f(x_{S})]$
(1)

Once generated, SHAP uses these values to display plots for both global
explanations and local explanations.

####  IV-D2 Global feature importance

We use the SHAP algorithm to generate global summary plots that highlight the
importance of each feature in the model’s decision-making similarly to what we
have done in Table II. Figures 1 and 2 display the importance plots for the
Label Encoded dataset and the One Hot Encoded dataset, respectively.

Figure 1: Label Encoding global importance plot

Figure 2: One Hot Encoding global importance plot

##  V Discussion

The main difference between these plots is that while we know what feature is
more important with Label Encoding, we know what exact value of that feature
is more important with One Hot Encoding. This means that we get more
specificity as a feature’s importance is the sum of the importance of its
unique values. A top ranking feature in the Label Encoding model could have
therefore reached its rank because of the importance of some of its values,
but not the others. Using One Hot Encoding, we can single out what values
exactly are the most relevant to further analyze. For example, the
”MinorOperatingSystemVersion” feature in 1 has a mean SHAP value of almost
0.6, ranking fifth. However, in 2, we can see that is actually Version 3 of
this feature that is really impactful, ranking first with a mean SHAP value of
more than 1.2. Yet, version 1 of this feature only has a score of almost 0.2,
and the rest of the version are not in the top 10 features. So using One Hot
Encoding, we can single out files with the Version 3 of
”MinorOperatingSystemVersion” and further analyze them separately in hopes of
creating an easy rule for them or see what more we can learn. One drawback of
this plot is that it is not easy to read when we have hundreds or thousands of
features. In this example, we have 16087 features. It will be unproductive to
use this plot to study feature importance. Instead, we can extract the raw
SHAP values of all one hot encoded features, group them by original feature,
and plot them side by side in another plot. We propose the plots in Figures 3
and 4 where we plot the importance of the different values of the
”MajorSubsystemVersion” feature side by side, horizontally and vertically
respectively. We chose this feature instead of the number 1 ranking
”MinorOperatingSystemVersion” feature because it has considerably fewer
distinct values making it easier to plot, wasting less space and delivering
the same message. These figures allow us to better visually grasp the
relativity in importance between the different values of a feature. This way,
we can add or remove values to and from a watchlist and also construct rules
for particular values. We can now combine this with the confidence score of
the model at inference to start a routine, a check or apply a rule when the
score doesn’t hit the certainty thresshold. At that point, we would start
investigating individual instances, thus needing different explanations called
local explanations.

Figure 3: Horizontally stacked bar plot of the ”MajorSubsystemVersion”’s
distinct values importance

Figure 4: Vertically stacked bar plot of the ”MajorSubsystemVersion”’s
distinct values importance

####  V-1 local feature importance

Local explanations focus on individual instances, displaying to the user the
step-by-step contribution of each feature on the model’s decision. Using
SHAP’s local explanation plots, we get Figures 5 and 6 which display the local
explanation of instances 2 and 3 respectively, first using Label Encoding
first and then One Hot Ecoding.

Figure 5: Local explanation for test observation number 2

Figure 6: Local explanation for test observation number 3

Again, the added refinement of the exact feature value gives us a lot more
insight into what pushed the model towards a certain classification. Although
the one hot encoding in this case may seem useless since we already know what
value of each feature the instance holds, it instead can be used as an
assertion method to make sure there are no anomalies in the decision shifting.
Finally, we can see that being trained on the individual values changes the
base value and decision shift intensity of each feature, as it has been
trained on more finegrained data and the model had the chance to learn
combinations that go together. These combinations in a tree based model such
as XGBoost can then be used extracted and used as normal conditional IF rules
or analyzed to detect vulnerabilities that went under the radar. Even then,
the feature encoding will have an impact on the generated rules.

####  V-2 IF-Rules

IF-Rules are logical statements that express conditional relationships between
input variables and output decisions and follow a simple structure: IF a
specific condition or set of conditions is satisfied, THEN a particular action
or decision should be taken. The conditions and actions are typically
expressed using logical operators, such as ”AND,” ”OR,” and ”NOT.” IF rules
provide a transparent and interpretable way to encode domain knowledge and
decision-making criteria into a system. Due to their nature, tree-based models
can be seen as a collection of IF rules combined together to form a decision-
making process. Each node in a decision tree represents an IF statement on a
specific feature or attribute, and the tree structure guides the flow of
decision-making based on these conditions. The splitting criteria at each node
determine the conditions for branching into different paths, leading to
subsequent nodes or leaves with specific outcomes or predictions. Since
XGBoost is a tree based model, we can extract the IF-Rules it learned during
the training phase and use them to build logical pipelines or to study them.
An example of the IF-Rules learned by our XGBoost model can be seen in Figures
7 and 8 for Label Encoding and One Hot Encoding respectively.

Figure 7: Example IF-Rules for Label Encoding

Figure 8: Example IF-Rules for One Hot Encoding

While there is no apparent difference between the IF-Rules of the two encoding
techniques, the difference lies in the metadata. In Table IV, we can see the
difference in the rules’ total text length in number of characters as well as
the explanation file size in KB. We can see that One Hot Encoding resulted in
less characters which means less file size. The indirect consequence of this
is less analysis time, less system complexity and less ambiguity, all of which
directly benefit analysts and systems.

TABLE IV: Difference in IF-Rules between the encoding techniques

Encoding Technique | F1 | Rules file size | Rules text length
---|---|---|---
Label Encoding | 0.991 | 221 KB | 226065 characters
One Hot Encoding | 0.982 | 175 KB | 180237 characters

##  VI Conclusion

In this work, we studied the impact of preprocessing on the explainability
power of XAI algorithms. We took a malware classification dataset as an
example on which we trained an XGBoost model. We tried two different types of
feature encoding: Label Encoding and One Hot Encoding and found there is a
marginal performance loss by the model by using OHE instead of LE. That loss
was made up with thanks to the more detailed explanations we managed to make
thanks to OHE. We found that OHE alllows us to go deeper in the details when
searching for answers, both globally and locally. We also found that using OHE
yields smaller explanation files and results in less time spent analyzing by
human analysts. We think this is an interesting aspect to be taken into
consideration when working with XAI and could be expanded by including more
feature encoding techniques and more creative plots.

## References

  * [1] T. Mikolov, K. Chen, G. Corrado, and J. Dean, “Efficient estimation of word representations in vector space,” 2013.
  * [2] J. H. Friedman, “Greedy function approximation: a gradient boosting machine,” _Annals of statistics_ , pp. 1189–1232, 2001.
  * [3] M. T. Ribeiro, S. Singh, and C. Guestrin, “”why should i trust you?”: Explaining the predictions of any classifier,” 2016.
  * [4] S. Lundberg and S.-I. Lee, “A unified approach to interpreting model predictions,” 2017.
  * [5] M. Sundararajan, A. Taly, and Q. Yan, “Axiomatic attribution for deep networks,” 2017.
  * [6] M. T. Ribeiro, S. Singh, and C. Guestrin, “Anchors: High-precision model-agnostic explanations,” _Proceedings of the AAAI Conference on Artificial Intelligence_ , vol. 32, no. 1, Apr. 2018. [Online]. Available: https://ojs.aaai.org/index.php/AAAI/article/view/11491
  * [7] K. Sokol and P. A. Flach, “One explanation does not fit all,” _KI - Künstliche Intelligenz_ , vol. 34, pp. 235–250, 2020.
  * [8] F. Cabitza, A. Campagner, G. Malgieri, C. Natali, D. Schneeberger, K. Stoeger, and A. Holzinger, “Quod erat demonstrandum? - towards a typology of the concept of explanation for the design of explainable ai,” _Expert Systems with Applications_ , vol. 213, p. 118888, 09 2022.
  * [9] S. R. Islam, W. Eberle, S. K. Ghafoor, A. Siraj, and M. Rogers, “Domain knowledge aided explainable artificial intelligence for intrusion detection and response,” 2020.
  * [10] M. Hind, D. Wei, M. Campbell, N. C. F. Codella, A. Dhurandhar, A. Mojsilović, K. N. Ramamurthy, and K. R. Varshney, “Ted: Teaching ai to explain its decisions,” 2019.
  * [11] A. Arias-Duart, F. Parés, D. Garcia-Gasulla, and V. Gimenez-Abalos, “Focus! rating xai methods and finding biases,” 2022.
  * [12] M. Setzu, R. Guidotti, A. Monreale, F. Turini, D. Pedreschi, and F. Giannotti, “GLocalX - from local to global explanations of black box AI models,” _Artificial Intelligence_ , vol. 294, p. 103457, may 2021.
  * [13] R. Ronen, M. Radu, C. Feuerstein, E. Yom-Tov, and M. Ahmadi, “Microsoft malware classification challenge,” 2018.