PHYSICAL REVIEW E 99, 043206 (2019)

Objective magnetic vortex detection

Erico L. Rempel,1,2,* Tiago F. P. Gomes,2 Suzana S. A. Silva,1 and Abraham C.-L. Chian1,2,3
1
Institute of Aeronautical Technology (ITA), World Institute for Space Environment Research (WISER),
São José dos Campos, São Paulo 12228-900, Brazil
2
National Institute for Space Research (INPE), P.O. Box 515, São José dos Campos, São Paulo 12227-010, Brazil
3
University of Adelaide, School of Mathematical Sciences, Adelaide, South Australia 5005, Australia

(Received 7 November 2018; revised manuscript received 25 February 2019; published 29 April 2019)

Magnetic coherent vortical structures are ubiquitous in space and astrophysical plasmas and their detection
is key to understanding the nature of the intrinsic turbulence in those conducting fluids. A recently developed
method to detect magnetic vortices is explored in problems of two- and three-dimensional magnetohydrodynamic
simulations. The integrated averaged current deviation, the normed difference of the current density at a point
and the mean current density in the domain, integrated along a magnetic field line, is proved to be objective, i.e.,
invariant under rotations and translations of the observer. The method is shown to detect accurately the boundary
of magnetic vortices in two-dimensional simulations, as well as magnetic flux ropes in three dimensions.

DOI: 10.1103/PhysRevE.99.043206

I. INTRODUCTION conditions of [4] to facilitate the generation of small-scale

structures, and (ii) an initial state having dominant scales that
Vortices are the natural building blocks of turbulence, since
are small compared with the system size, defined by a broad
the nonlinear energy cascade where large vortices break up
spectrum of modes with random phases, to allow an inverse
into smaller vortices is thought to be responsible for the
cascade to develop. They obtained a turbulent power spectrum
scaling laws found in kinetic energy spectra in fully de-
that exhibits both Kraichnan and Kolmogorov spectral indices
veloped turbulence [1]. These persistent coherent structures
in the inertial range and a modified dissipation scale driven by
that enclose regions with high vorticity have been studied
intermittency in the dissipation range. Politano et al. [6] per-
for decades in hydrodynamical problems, but there is still
formed high-resolution simulations of the Orszag-Tang model
no standard definition of what a vortex is or how it can be
of MHD turbulence at Reynolds numbers of several thousand,
detected (see, e.g., Haller et al. [2] and references therein). In
yielding kinetic and magnetic energy spectra that display an
the context of space and astrophysical plasmas, it is usually
inertial range extending over one decade of wave numbers,
imperative to consider the presence of magnetic fields, with
and a dissipation range where resistive tearing instabilities
their own vortical structures. In two spatial dimensions, these
destabilize current sheets generated by the inertial dynamics,
structures exhibit a curling in-plane magnetic configuration,
leading to the formation of small-scale magnetic islands. They
with field lines encircling or spiraling around a center. Kinney
tested the generality of their results of the Orszag-Tang vortex
et al. [3] formally defined magnetic vortices as axisymmetric
by considering flows with random initial conditions. Wu and
monopole patterns in the current field. In three dimensions,
Chang [7] reported the onset and interactions of localized co-
magnetic vortices constitute a bundle of helical field lines in a
herent structures in two-dimensional MHD simulations. They
tubular region. This work focuses on frame invariant detection
noticed that the structures can merge locally and sporadically
of two-dimensional (2D) and three-dimensional (3D) mag-
in the absence of a background magnetic field. When there
netic vortices in magnetohydrodynamic (MHD) simulations.
is a background sheared magnetic field, the structures tend to
A number of papers have studied current sheets, mag-
align with the background current sheet, a process that was
netic islands, vortices, and magnetic reconnection in two-
related to an intermittent turbulence model for the Earth’s
dimensional incompressible MHD turbulence. For example,
magnetotail. Armstrong and Craig [8,9] studied viscous and
Orszag and Tang [4] introduced a two-dimensional prototype
resistive losses in transient magnetic reconnection driven by
model for numerical simulation of small-scale structures (sin-
the Orszag-Tang vortex and concluded that viscous dissipa-
gularities) of vortices and current sheets in an incompressible
tion in solar flaring plasmas can account for a significant
MHD flow, and showed that 2D MHD turbulence is less singu-
fraction of the flare energy release.
lar (containing less small-scale structures) than 3D hydrody-
In three spatial dimensions, magnetic vortices are often
namic turbulence, but more singular than 2D hydrodynamic
known as flux ropes [10]. Inside these tubular or cylindrical
turbulence. Biskamp and Welter [5] investigated the dynam-
regions with spiraling or helical magnetic field lines, the axial
ics of 2D decaying incompressible MHD turbulence using
magnetic field is much larger than the magnetic field outside.
two types of initial states: (i) a generalization of the initial
However, this definition does not provide a precise descrip-
tion of the vortex boundary, since a user-defined magnetic
threshold is required to determine whether field lines are
*
rempel@ita.br considered inside or outside the rope. Hence, several works

2470-0045/2019/99(4)/043206(15) 043206-1 ©2019 American Physical Society

REMPEL, GOMES, SILVA, AND CHIAN PHYSICAL REVIEW E 99, 043206 (2019)

on flux tubes or ropes display those structures by plotting the random vortices as initial conditions, and three-dimensional
isosurfaces of the magnetic field or current density, where compressible simulations of interlocked flux rings as initial
the user chooses some energy level to define the surface to conditions. Section V gives the conclusions.
be plotted. For example, in a 3D data set assembled using
the relaxation scaling experiment device, Sears et al. [11] II. OBJECTIVE KINEMATIC AND MAGNETIC VORTICES
detected the signature of flux ropes transversal to a horizontal
plane as regions where the contours of flux rope current A. Objectivity
density J z (x, y) fall below a certain threshold. Isosurfaces The concept of objectivity refers to invariance under
and contour plots of different quantities related to the parallel time-dependent transformations of a reference frame or
current and helical magnetic field components were used to observer [17]. The invariance of a quantity is important in the
locate flux ropes in MHD simulations of the Pegasus Toroidal context of vortex identification because a detection technique
Experiment in [12]. In [13] the thresholding is applied to the must be able to correctly classify a flow region as a vortex
magnetic helicity mapping, a quantity defined as the magnetic independent of the observer. The topic has been extensively
helicity normalized by the magnitude of the magnetic field, discussed in the literature (see, e.g., [18]).
integrated along the field lines. The difficulty in defining The invariance of a physical quantity is also important
vortex boundaries also occurs in two dimensions, where the when there is no time evolution involved, i.e., in steady
vortex centers can be identified as local extrema of the current flows. For example, the description of coherent structures in
[3], but their boundaries must be inferred in some way from an instantaneous magnetic field might depend on the frame
contour plots. used, since an observer O translated and rotated in relation to
Rempel et al. [14] proposed a new and simple way to another observer Ō may see the magnetic field lines with a
define magnetic vortices based on the integrated averaged different topology. As an example, suppose that observer O
current deviation (IACD), which is the integral of the normed measures the following magnetic field at time t0 :
difference between the current density along a magnetic field ⎛ ⎞
line and the current density averaged in the spatial domain. 0 1 0
The vortex boundary is defined as the outermost convex B(x) = ⎝1 0 0⎠x, x ∈ R3 . (1)
contour surrounding a local maximum of the IACD field. The 0 0 1
method is based on an adaptation of the Lagrangian averaged The magnetic field lines observed by O at time t0 are
vorticity deviation (LAVD), introduced by Haller et al. [2] for shown in Fig. 1(a) and reveal diverging lines surrounding a
kinematic vortices and recently employed in the detection of saddle fixed point in the middle of the box. The field lines are
coherent swirls in an ocean general circulation model simula- solutions of the initial value problem
tion of the Gulf of Mexico [15]. A few preliminary results on
IACD were shown in [14] using numerical simulations of a dx
= B(x(s)), x(s0 ) = x0 , (2)
nonlinear dynamo. In the present paper, we apply the IACD ds
to magnetic vortex detection in a series of two- and three- where x(s) is the position vector and s is a scalar parameter
dimensional MHD simulations to test its robustness in differ- related to the distance l along the field line by dl = |B|ds.
ent circumstances, both in compressible and incompressible Next, suppose that observer Ō is rotated and translated in
flows. The IACD method is shown to be effective in detecting relation to O, such that position vectors are related by the
vortex boundaries in which the magnetic field lines spiral following Euclidean transformation:
together without reconnecting and spreading away. In three
dimensions, the field lines from boundaries found on a plane x̄ = c(s) + Q(s)x, s̄ = s + a, (3)
section keep their coherence in a tubular region, showing that
the method can be employed to detect flux ropes. Possible where Q(s) is an arbitrary orthogonal rotation matrix, c(s)
applications in the context of solar plasmas are discussed. We is a translation vector, and a is an arbitrary real constant.
also prove that IACD vortices are objective, which means that Following Eq. (2), the magnetic field observed by Ō can be
they remain unchanged under translations and rotations of the obtained by deriving Eq. (3):
coordinate frame. This point is fundamental in a definition of d x̄ dc dQ dx
a magnetic vortex, and is a hallmark of the IACD method. B̄ = = + x+Q
ds ds ds ds
If a detection method is not objective it means that a perfectly
= c + Q x + QB
straight magnetic field line could be classified as a twisted line
when measured with respect to a rotating frame (see Fig. 5 in = c + Q Q (x̄ − c) + QB, (4)
[16]).
In Sec. II, we discuss the important concept of objectivity, where  denotes the transpose operator.
which guarantees that an operator is invariant under contin- Now, suppose that Q and c are given by
⎛ ⎞ ⎛ ⎞
uous changes of an observer. The concept is applied to the cos(s) 0 sin(s) 0
LAVD and IACD operators. A precise definition of a magnetic Q=⎝ 0 1 0 ⎠, c = ⎝0⎠. (5)
vortex is provided. Section III contains a detailed proof of −sin(s) 0 cos(s) 0
the objectivity of IACD. Section IV presents the results of
MHD simulations in three situations: two-dimensional in- Then, for s = 2.1, the field lines computed from B̄ as observed
compressible simulations with Orszag-Tang vortices as initial by Ō behave as in Fig. 1(b). The field lines reveal swirling
conditions, two-dimensional incompressible simulations with patterns that could be classified as a magnetic vortex. This

043206-2
OBJECTIVE MAGNETIC VORTEX DETECTION PHYSICAL REVIEW E 99, 043206 (2019)

discuss objectivity in terms of position and s. Suppose that an

event is characterized by observers O and Ō, respectively, by
(x, s) and (x̄, s̄), which are related by a Euclidean transforma-
tion as in Eq. (3). Then, a point difference vector v is defined
as objective or frame indifferent if it conserves its magnitude
in different frames; thus, v is objective if it transforms as
v̄ = Q(s)v, (6)
since the rotation matrix preserves the length of a vector. Also,
in the context of steady fields, a tensor is said to be objective
if
T̄ = Q(s)TQ(s) . (7)
Equations (3), (6), and (7) can be applied when describing
the properties of a steady flow or a magnetic vector field at a
given time.

B. Objective vortex detection

1. Lagrangian averaged vorticity deviation
Most of the existing methods for defining or detecting
vortices are not objective. Haller et al. [2] developed a method
to define coherent kinematic vortices objectively. Consider a
velocity field u in a spatial domain U . A fluid particle at initial
position x(t0 ) = x0 ∈ U follows a trajectory in space given by
solutions of the initial value problem
dx
= u(x, t ), x(t0 ) = x0 . (8)
dt
The Lagrangian averaged vorticity deviation is defined as
 t0 +τ
t0 +τ
LAVDt0 (x0 ) = |ω(x(t ), t ) − ω(t )|dt, (9)
t0

where τ is a given time interval related to the typical life-

time of the vortices that are to be found, ω = ∇ × u is the
vorticity, and ω(t ) is the spatial mean of the vorticity at
time t. An objective kinematic Lagrangian vortex is, then,
defined as an evolving convex material domain VL (t ) such
FIG. 1. Magnetic field lines computed from a vector field mea- that VL (t0 ) is filled with a set of level curves (for two-
sured by two different observers: (a) an observer O detects a saddle dimensional flows) or surfaces (for three-dimensional flows)
point structure and (b) another observer Ō rotated and translated in of LAVDtt00 +τ (x) with outward nonincreasing LAVD values.
relation to O detects a magnetic vortex. In practice, one can search for kinematic Lagrangian vortices
in a two-dimensional flow by computing the LAVD field
shows that magnetic field lines are not an objective quantity, on the plane, extracting the initial positions of vortices as
since the same magnetic structure could be classified as either the local maxima of LAVD, and then computing a vortex
a saddle point or a vortex core, depending on the observer. boundary as the outermost convex closed contour of LAVD
A practical situation is the study of magnetic structures in surrounding a vortex center, i.e., a local LAVD maximum.
space and solar plasmas as observed by different satellites In three-dimensional flows, the above procedure can be used
at a given time. One would like them both to agree on the to find the intersection of a vortex with a plane in the flow
kind of structure observed, and correcting the orientation of domain.
one satellite to match the other does not solve the problem, In order to take into account small deviations from convex-
since they could both accuse the detection of a vortex when in ity that are always present in data with finite resolution, the
fact they are observing a rotated saddle structure. An example convexity deficiency was introduced by [2]
similar to the one in Fig. 1 is provided by [19] in the context Ac − Ach
of continuous rotations of a frame in two-dimensional time- ε := , (10)
dependent velocity fields. Ac
The problem described above is essentially the same as where Ac is the area enclosed by the contour extracted and Ach
computing the streamlines traced by a passive scalar in a is the area enclosed by its convex hull. Thus, an outermost
steady velocity field as t evolves. In the case of magnetic field closed contour is classified as a vortex boundary if ε is
lines, the line parameter s is used in place of t, and one can less than a given threshold. The threshold value is problem

043206-3
REMPEL, GOMES, SILVA, AND CHIAN PHYSICAL REVIEW E 99, 043206 (2019)

dependent and must be chosen after some trial and error.

It must be small enough to avoid classifying the numerical P z

approximation of nonconvex structures as convex and large

Z
enough to account for the ragged edges resulting from the P x
r
line
discretization of the LAVD field. Previous experiments with field
)
observational and simulated data in two and three dimensions
R C(s
using a LAVD field with 5122 [2] or 5123 [14] grid points ^
y

revealed good results using ε between 10−3 and 10−5 . Z

^ ^ Y
X Y
2. Integrated averaged current deviation
In order to detect magnetic vortices, Rempel et al. [14]
introduced a variation of the LAVD technique based on the
X
magnetic field and the current density. Consider the magnetic
field line equation [Eq. (2)] at t = t0 . Then, the IACD is FIG. 2. “Lagrangian” and spatial coordinates of a deformable
defined as body.
 s0 +ξ
s0 +ξ
IACDs0 (x0 ) = |J(x(s), t0 ) − J(t0 )|ds, (11)
s0 tensor components can be written as [20]
⎛ ⎞
where J = ∇ × B/μ0 is the current density and J(t0 ) is ∂x Bx ∂y Bx ∂z Bx
the mean current density. Note that, differently from LAVD, ∇B = ⎝∂x By ∂y By ∂z By ⎠.
the calculations of IACD are performed for a fixed time
∂x Bz ∂y Bz ∂z Bz
t0 , due to the fact that the field line equation, Eq. (2), is
defined for a fixed time. The IACD is an objective quantity, The magnetic gradient tensor can be written in terms of its
which guarantees its invariance under arbitrary, continually symmetric (S) and antisymmetric (A) parts,
varying observer changes along magnetic lines. Thus, given
arbitrary rotation matrix Q(s) and translation vector c(s), if ∇B = S + A, (13)
J̄ = Q(s)J + c(s), then where
|J̄ − J̄| = |J − J|, (12)
∇B + ∇B
as demonstrated in Sec. III. S= , (14)
2
For a given parameter interval [s0 , s0 + ξ ], a magnetic
∇B − ∇B
vortex is defined as a convex domain VL (t0 ) such that VL (t0 ) A= . (15)
is filled with a set of tubular level curves or surfaces of 2
IACDss00 +ξ (x0 ) with outward nonincreasing IACD values; the Substituting the components of ∇B in A, we obtain the
boundary of VL (t0 ) is the outermost convex level curve or sur- following matrix:
face of IACDss00 +ξ (x0 ) in VL (t0 ). The numerical computation of ⎛ ⎞
magnetic vortices in two- and three-dimensional fields follows 0 ∂y Bx − ∂x By ∂z Bx − ∂x Bz
1⎜ ⎟
the same procedure described above for objective kinematic A = ⎝∂x By − ∂y Bx 0 ∂z By − ∂y Bz ⎠.
vortices, just changing LAVD by IACD. 2
∂x Bz − ∂z Bx ∂y Bz − ∂z By 0
The parameter ξ determines the length along the field lines
inside the structure to be considered as a coherent magnetic (16)
vortex and is usually chosen after some trial and error, de-
From Ampère’s law,
pending on each problem. In three-dimensional systems, ξ is
related to the length of magnetic flux ropes, so larger values of ∇ × B = (∂y Bz − ∂z By )x̂ + (∂z Bx − ∂x Bz )ŷ
ξ should be used to detect longer ropes, which will typically
be thinner than the ropes found with smaller ξ , as illustrated + (∂x By − ∂y Bx )ẑ = μ0 J. (17)
in Sec. IV B. Thus, matrix A can be written as
Solving Eq. (2) requires integrating field lines in a con- ⎛ ⎞
tinuous space domain. Since we are working with discrete 0 −Jz Jy
μ0 ⎝ μ0
data for B obtained from numerical simulations of the MHD A= Jz 0 −Jx ⎠ ⇒ A = J, (18)
equations, we use cubic splines interpolation for the B field 2 2
−Jy Jx 0
and solve Eq. (2) with a fourth-order Runge-Kutta integrator.
where J is defined as the current density tensor [20]. This
tensor has the property that J = −J .
III. PROOF OF THE OBJECTIVITY OF IACD
Now, let r = r(R, s) be a spatial vector position and R =
The objectivity of LAVD was demonstrated in [2] using the R(r, s) the associated material vector position, where s is a pa-
velocity gradient tensor and the definition of the spin tensor. rameter. The components of r and R provide the “Lagrangian”
We follow a similar approach for IACD based on the magnetic and spatial coordinates, respectively. Figure 2 illustrates the
gradient and current density tensors. The magnetic gradient Eulerian and Lagrangian position vectors of a point P in a

043206-4
OBJECTIVE MAGNETIC VORTEX DETECTION PHYSICAL REVIEW E 99, 043206 (2019)

(a) (a)

(b) (b)

(c)
FIG. 4. Orszag-Tang vortices at t = 1: (a) the magnetic field
vectors and (b) line integral convolution plot showing the magnetic
field lines.

description of the deformation gradient tensor is given by

⎛ ∂x ∂x ∂x
⎞
∂X ∂Y ∂Z
∂r ⎜ ∂y ∂y ∂y ⎟
F= = ⎝ ∂X ∂Y ∂Z ⎠
. (19)
∂R ∂z ∂z ∂z
∂X ∂Y ∂Z

The transformed deformation gradient tensor is given by

∂r̄ ∂ (Q(s)r + c(s)) ∂r
F̄ = = =Q ⇒ F̄ = QF . (20)
∂R ∂R ∂R
The derivative of the deformation gradient tensor with
FIG. 3. Orszag-Tang vortices at t0 : (a) the magnetic field vectors,
(b) an enlargement of the central region of (a), and (c) line integral
respect to the parameter s is defined by
convolution plot showing the magnetic field lines. The large arrows ∂F ∂ ∂r ∂ ∂r ∂B
indicate the location of an X point. F = = = ⇒ F = , (21)
∂s ∂s ∂R ∂R ∂s ∂R
∂r
where we used the field line equation ( ∂s = B) and the fact
body represented by the red line. In a frame following the field that the material coordinate R does not depend on the parame-
line parametrized by s (the “Lagrangian” frame), the material ter s. From Eq. (21), we can write the magnetic gradient tensor
coordinates are fixed for a given s. in terms of the deformation gradient tensor as
The material coordinates are used to define the undeformed
reference configuration and the spatial coordinates define the ∂B ∂B ∂r
F = = ⇒ F = ∇BF ⇒ ∇B = F F−1 . (22)
deformed current configuration. In this sense, the Lagrangian ∂R ∂r ∂R

043206-5
REMPEL, GOMES, SILVA, AND CHIAN PHYSICAL REVIEW E 99, 043206 (2019)

(a) (a)

(b) (b)

FIG. 5. Orszag-Tang current density at t0 : (a) the normalized FIG. 6. Orszag-Tang IACD at t0 with ξ = 100: (a) the normal-
magnitude of the current density and the vortices detected as convex ized IACD field and the vortices detected as convex closed curves
closed curves (black) from the contour plot of |J| and (b) the line from the contour plot of IACD and (b) the line integral convolution
integral convolution plot of the magnetic field lines (dark gray) and plot of the magnetic field lines (dark gray) nearby one of the vortices.
the field lines starting from initial conditions (light gray [yellow The circles are initial conditions on the vortex detected by IACD and
online] circles) on the J vortex boundary shown in (a) (red). the yellow lines are the corresponding magnetic field lines obtained
by solving Eq. (2).
Using Eqs. (20) and (22) we can show that the magnetic
gradient tensor is not objective,
where the relation Q Q = −QQ  was used [21]. The cur-
∇¯ B̄ = F̄ F̄−1 = (QF) (QF)−1 rent density tensor J, being proportional to tensor A, Eq. (18),
= (Q F + QF )(F−1 Q−1 ) is also not objective.
Although the antisymmetrical tensor A is not objective, the
= Q FF−1 Q + QF F−1 Q
averaged deviation is an objective quantity:
= Q Q + Q∇BQ , (23)
Ā − Ā = QAQ + Q Q − QAQ + Q Q 
where the fact that Q is an orthogonal matrix was used
(Q = Q−1 ). = QAQ + Q Q − QAQ − Q Q
We can use the results from Eq. (23) to prove that the = Q(A − A)Q . (25)
antisymmetrical part of ∇B, Eq. (15), is not objective either:
 Although we have shown that the averaged deviation of the
∇¯ B̄ − ∇¯ B̄
Ā = antisymmetrical tensor is objective, we are interested in the
2
objectivity of the current density tensor. For this purpose, let
1   
= [Q Q + Q∇BQ − (Q Q + Q∇BQ ) ] us write the tensor A in Eq. (18) in index notation,
2
1   μ0
=

[Q Q + Q∇BQ − QQ − Q∇B Q ] Ai j = − εi jk Jk ,
2 2
μ0 μ0
1 Ai j e j = − εi jk Jk e j = εik j Jk e j , (26)
= [Q(∇B − ∇B )Q + 2Q Q ] 2 2
2 μ0
Ā = QAQ + Q Q , (24) Ae = (J × e),
2

043206-6
OBJECTIVE MAGNETIC VORTEX DETECTION PHYSICAL REVIEW E 99, 043206 (2019)

(a)

FIG. 7. Orszag-Tang current density at t = 1: the normalized

(b)
magnitude of the current density and the small vortices detected as
convex closed curves (rings) from the contour plot of |J|.

where e is an arbitrary objective vector in R3 , εi jk is the

Levi-Civita tensor [22], and J is the current density vector
associated to tensor A.
From Eq. (24),
Ae = Q ĀQe − Q Q e. (27)
Now, using the objective transformation ē = Qe we can
rewrite Eq. (27) as
Ae = Q Āē − Q Q Q ē. (28)
FIG. 8. Orszag-Tang IACD at t = 1, with ξ = 25: (a) the nor-
Analogously, we can define a vector α associated with the malized IACD field and the vortices detected as convex closed curves
tensor Q Q . Following in Eq. (26), since Q Q is an an- from the contour plot of IACD and (b) the line integral convolution
tisymmetrical tensor, its dual vector α is a non-null vector plot of the magnetic field lines (dark gray) and the field lines starting
uniquely defined by [23] from initial conditions (light gray [yellow online] circles) on an
μ0 IACD vortex boundary shown in (a).
Q Q ē = (α × ē). (29)
2
By substituting Eqs. (26) and (29) into Eq. (28) we obtain Note that the mean value J is also not an objective
quantity, since by taking averages in Eq. (30) we have
μ0 μ0 μ0
(J × e) = Q (J̄ × ē) − Q (α × ē) J = Q J̄ − Q α. (33)
2 2 2
In the particular case when J = 0, J̄ = α, which is not
J × e = Q (J̄ × ē) − Q (α × ē) null, since α is uniquely defined by Eq. (29), where Q is an
J × e = [(Q J̄) × (Q ē)] − [(Q α) × (Q ē)] arbitrary rotation matrix.
J × e = (Q J̄) × e − (Q α) × e
IV. RESULTS
J × e = [Q J̄ − Q α] × e
A. Two-dimensional incompressible MHD
J = Q J̄ − Q α. (30)
In dimensionless conservative form, the ideal MHD equa-
The transformation described in Eq. (30) demonstrates that the tions for an incompressible fluid are
current density vector is not an objective quantity. However,
∂u
using Eq. (30), we obtain that the averaged current deviation = −∇ p − ∇ · (uu − BB), (34)
is an objective vector: ∂t
∂B
J − J = Q J̄ − Q α − Q J̄ − Q α = −∇ · (uB − Bu), (35)
   
∂t
= Q J̄ − Q α − Q J̄ + Q α
∇ · u = 0, (36)
= Q (J̄ − J̄). (31)
∇ · B = 0, (37)
Therefore, the IACD method, Eq. (11), is objective, for
where u is the velocity, B is the magnetic field, and p is
|J̄ − J̄| = |Q(J − J)| = |J − J|. (32) the total pressure. Equations (34)–(37) are solved in a square

043206-7
REMPEL, GOMES, SILVA, AND CHIAN PHYSICAL REVIEW E 99, 043206 (2019)

(a) (b) (c)

(d) (e) (f)

FIG. 9. Random vortices Eulerian fields at t = 300 (upper panel) and t = 386 (lower panel): (a, d) magnetic vectors, (b, e) magnetic field
streamlines, and (c, f) current density.

domain with sides L = 2π and periodic boundary conditions the magnetic shear and consequently enhancing the current
in both directions. Numerical integration is performed with a density in that region.
finite difference method proposed by Wu [24], where a fifth- As indicated by Eq. (11), the IACD is a measure of the
order weighted essentially non-oscillatory (WENO) scheme deviation of the current density from its mean value, inte-
is used in the spatial discretization, a variable-step fourth- grated along magnetic field lines. In order to distinguish the
order Runge-Kutta method is employed for time integration, qualities of IACD from the current density itself, we compare
and the incompressibility condition is enforced through the the vortex detection using J and IACD. Figure 5 illustrates
fractional-step method of Kim and Moin [25]. The numerical the detection using J at t0 . Figure 5(a) shows a color plot
resolution is 512×512 grid points. In the following sections, of |J| (normalized by max(|J|)) with two vortices indicated
Eqs. (34)–(37) are solved with two different sets of initial by black lines. Those vortices were found by employing the
conditions. algorithm described in Sec. II B, where first a local maximum
of the magnitude of the current density field is found, then
the outermost closed convex contour line surrounding that
1. Deterministic initial conditions: Orszag-Tang vortex
maximum is extracted using MATLAB’s contourc function.
The Orszag-Tang (OT) vortex is defined by the following Figure 5(b) shows the magnetic field streamlines superposed
stream function ψ and magnetic potential a: by one of the two vortices (red line). The yellow spots are
initial conditions from where the magnetic field lines are
ψ (x, y) = 2(cos x + cos y), (38) computed using Eq. (2). It can be seen that most lines spread
away from the vortex region; therefore, the coherent structure
a(x, y) = cos 2x + 2 cos y. (39) identified from J is not a good representation of a magnetic
vortex.
Figure 3 shows the state of the magnetic field at t0 = 0.01, Next, we compute the IACD field at t0 using Eq. (11) with
where Fig. 3(a) depicts the magnetic field vectors and ξ = 100 and a convexity deficiency of ε = 10−3 ; the result
Fig. 3(b) is an enlargement showing the X point formed at is shown in Fig. 6. The top panel is the IACD field (normal-
the vortex interface. The magnetic field lines are exhibited ized by max(IACD)) and the vortices centered on two local
in Fig. 3(c) as a line integral convolution (LIC) plot, which maxima, denoted by the yellow lines. When initial conditions
displays the integral curves of (Bx , By ) in different shades of for Eq. (2) are chosen on the outermost convex contour sur-
gray. After one time unit, the magnetic field structures are rounding an IACD maximum (red spots in the lower panel),
distorted as in Fig. 4. The vectors near the interface between the resulting magnetic field lines (yellow) remain confined
the two central vortices have become larger, thus increasing to the vortical area, just inside the region surrounded by

043206-8
OBJECTIVE MAGNETIC VORTEX DETECTION PHYSICAL REVIEW E 99, 043206 (2019)

(a) (a)

(b) (b)

FIG. 10. Random vortices current density field at t = 300: FIG. 11. Random vortices IACD at t = 300 from s0 = 0 to s0 +ξ ,
(a) surface plot of the |J| field with rings denoting the outermost con- with ξ = 100: (a) surface plot of the IACD field with rings denoting
vex closed level curves surrounding selected local maxima; (b) the the outermost convex closed level curves surrounding selected local
same as (a) but in a 2D color plot. maxima; (b) the same as (a) but in a 2D color plot.

2. Random initial conditions: Wu-Chang model of magnetotail

the separatrix lines starting at the two neighboring X points,
which can be seen in the LIC plot. In this section, the MHD model for vortex interactions
When the current density and IACD are applied to vortex presented by Wu and Chang [7] is adopted to test the IACD
detection at t = 1, the results are as shown in Figs. 7 and 8, method. The initial conditions are given by random velocity
respectively. The current density failed to identify large vor- and magnetic fields, where the magnetic field is specified by
tices in the domain, as the main coherent structures are the a flux function. In the Fourier space, the two components of
current sheets formed at the vortex interfaces due to the the velocity and the magnetic flux function are proportional
strong magnetic shear. The shear also impacts the vortex to k 2 exp (−αk 2 ), where α = 0.01 is adopted, so the initial
detection using IACD, since its computation is based on kinetic energy spectrum is peaked at k = 10. The initial
|J|. Nonetheless, at t = 1 it still manages to extract the condition is such that the maximum amplitude is 0.2 for the
main vortical regions, as seen in Fig. 8. For larger t, the velocity field and 0.3 for the magnetic field.
shear increases even further, which impacts the numerical Figure 9 depicts the system state at t = 300 (upper panels)
computation of the closed contours; thus thresholds in |J| and t = 386 (lower panels). Figures 9(a) and 9(d) plot the
need to be imposed if one wants to use these techniques magnetic field vectors, Figs. 9(b) and 9(e) are line integral
to find the vortical structures. Since IACD is based on the convolution (LIC) plots, showing the instantaneous stream-
curl of B, it will be large on regions with rotation of mag- lines of the magnetic field, and Figs. 9(c) and 9(f) show the
netic field lines as well as shear regions, a problem shared current density plots. The arrows in the upper panel point to a
with the LAVD operator in velocity fields. For a discus- couple of magnetic vortices that move with the flow until they
sion on problems related to kinematic vortex detection in coalesce and merge into a single vortex, indicated by arrows
strong shear regions of observational velocity fields using in the lower panel.
LAVD, as well as possible solutions to these problems, In order to show this vortex-vortex interaction in more de-
see Silva et al. [26]. tail, we proceed to extract the vortex boundaries by applying
We now proceed with a different example where the initial the algorithm described in Sec. II B. First, we perform the
conditions produce a result that may be more relevant to analysis using the current density alone. Thus, a local max-
turbulence studies. imum of the magnitude of the current density field is found

043206-9
REMPEL, GOMES, SILVA, AND CHIAN PHYSICAL REVIEW E 99, 043206 (2019)

detect the magnetic vortex boundaries, outside which the mag-

netic field lines will tend to disperse away. This can be seen
in Fig. 12, which shows an enlargement of the rectangular
region in Fig. 11. The white circles denote a set of initial
conditions on the vortex boundary computed from the IACD
field; the white lines represent the trajectories of these initial
conditions under the magnetic field line equation [Eq. (2)].
Note that all magnetic field lines stay together in a coherent
bundle surrounding the vortex center for ξ = 100 s units.
However, when we start from initial conditions outside the
IACD vortex (green circles), the resulting magnetic field lines
spread away in different directions (green lines). Therefore,
the IACD method can detect magnetic vortex boundaries more
FIG. 12. Random vortices magnetic field lines at t = 300 from precisely than the |J| method in this example.
s0 = 0 to s0 + ξ , with ξ = 100. The white lines denote the field lines Back to the previously mentioned vortex merging, Fig. 13
starting from initial conditions (white circles) on an IACD vortex; shows an enlargement of the upper left corner of the domain
the light gray (green online) lines denote field lines starting from box, where the merging pointed out in Fig. 9 occurs. The up-
surrounding initial conditions. per panels represent the magnetic vectors and the lower panels
the LIC plots of the magnetic streamlines at different times,
until the two vortices coalesce into a larger vortex. These
and the outermost closed convex contour line surrounding coalescence events are important because they form local
that maximum is extracted. Figure 10 shows two views of sites of enstrophy production, conversion, and dissipation [3].
the magnitude of the current density at t = 300 with a set of The boundaries of these vortices were computed with IACD
vortices found with the given algorithm, where a convexity using ξ = 100 and are shown in Fig. 14, in both color scales
deficiency ε = 10−3 was adopted. A comparison between [Figs. 14(a)–14(d)] and surface plots [Figs. 14(e)–14(h)]. The
Figs. 9(b) and 10(b) seems to suggest that the numerically moment when the merging is effected can be defined by
computed vortex boundaries from the |J| field are underes- the presence of only one peak in the IACD surface plot in
timated and wider vortical regions should be expected. The Fig. 14(h).
choice of larger values of ε (within reasonable limits) does
not significantly change the results, which are related to the
intrinsic ruggedness of the |J| field. This problem is not B. Three-dimensional compressible MHD
present in the IACD vortices, shown in Fig. 11, where the The previous section showed that the IACD can be used to
IACD field is computed at time t = 300 from Eq. (11) with accurately detect magnetic vortex boundaries in planar fields,
s0 = 0 and ξ = 100. The IACD method managed to correctly where the structures detected were shown to match the LIC

(a) (b) (c) (d)

(e) (f) (g) (h)

FIG. 13. Random magnetic vortices merging: magnetic field vectors at (a) t = 362.34, (b) t = 372.32, (c) t = 384.54, and (d) t = 386.82,
respectively; (e–h) the line integral convolution plots corresponding to the magnetic field lines of (a–d). The arrow points to two vortices that
coalesce at t = 386.82.

043206-10
OBJECTIVE MAGNETIC VORTEX DETECTION PHYSICAL REVIEW E 99, 043206 (2019)

(a) (b)

(c) (d)

(e) (f)

FIG. 14. Random magnetic vortices merging in the IACD field: magnetic vortices detected from the IACD field at (a) t = 362.34,
(b) t = 372.32, (c) t = 384.54, and (d) t = 386.82, respectively; (e–h) surface plots of IACD corresponding to (a–d).

plots of the magnetic streamlines. However, the qualities of

IACD can be better appreciated in 3D fields. In this section, ∂ ln ρ
+ u · ∇ ln ρ = −∇ · u, (42)
the IACD operator is employed in detecting coherent mag- ∂t
netic structures in three-dimensional numerical simulations of
the resistive MHD equations for a compressible isothermal where u is the velocity, B is the magnetic field, ρ is the
gas. The model equations are the same as used in [27] for density, cs is the isothermal sound speed, p is the pressure,
simulations of hydromagnetic turbulence, J = ∇ × B/μ0 is the current density, with μ0 the constant
vacuum permeability, A is the magnetic vector potential, and
∂u ∇p J × B η is the constant magnetic diffusivity. The pressure gradient
+ u · ∇u = − + + F visc , (40)
∂t ρ ρ is written as ∇ p/ρ = cs2 ∇ ln ρ, where cs2 = p/ρ is constant.
∂A The viscous term is given by F visc = ν(∇ 2 u + ∇∇ · u/3 +
= u × B + η∇ 2 A, (41) 2S · ∇ ln ρ), where ν is the constant kinematic viscosity and S
∂t

043206-11
REMPEL, GOMES, SILVA, AND CHIAN PHYSICAL REVIEW E 99, 043206 (2019)

(a)

(b)

FIG. 15. Interlocked magnetic flux rings at t = 0.

is the traceless rate of strain tensor,

1 ∂ui ∂u j 2
Si j = + − δi j ∇ · u , (43)
2 ∂x j ∂xi 3

(c)

FIG. 17. Horizontal magnetic field components on z = 0 at

t = 100: (a) horizontal components (Bx , By ) of the magnetic field
vectors, (b) two-dimensional streamlines of the magnetic field com-
puted using only the horizontal components, and (c) the IACD field
with two magnetic vortices (black lines).

where δi j is the Kronecker delta. The domain is a square

box with sides L = 2π and periodic boundary conditions,
which implies that mass is conserved, ρ = ρ0 , where ρ0 is
FIG. 16. Magnetic field lines from random initial conditions at the initial uniform density and the brackets denote volume
t = 100: (a, b) two views of the same data from different angles. average. Nondimensional quantities are adopted where length

043206-12
OBJECTIVE MAGNETIC VORTEX DETECTION PHYSICAL REVIEW E 99, 043206 (2019)

is in units of L/2π , u is in units of cs , ρ is in units of ρ0 , and on the plane which sections the middle of the box at z = 0.
B is in units of (μ0 ρ0 cs2 )1/2 . Figures 17(a) and 17(b) show the horizontal components of
Equations (40)–(42) are solved with sixth-order centered the magnetic field vectors (Bx , By ) and their two-dimensional
finite differences in space and third-order variable-step streamlines, respectively. The presence of two large swirling
Runge-Kutta in time, using the PENCIL CODE [28] in a square regions interacting via the periodic boundary conditions can
box with sides L = 2π and periodic boundary conditions. The be seen. However, this view of the field on the plane does not
initial conditions are u = 0, ρ = ρ0 = 1. The magnetic initial take into account the torsion of the lines in the z direction, thus
condition consists of a pair of interlocked magnetic tori, or concealing the true complexity of the structures in the box.
flux rings, both with the same flux and with inner radius Figure 17(c) shows the IACDξs0=2000
=0 (x0 ) scalar field, where x0
Ri = 0.3 and outer radius Ro = 1, as seen in Fig. 15; inside belongs to a grid of points on the plane and the closed black
the tori ∇ × B = 0 and the initial B profile across each ring is curves delineate two magnetic vortices located in the central
a Gaussian. One torus is aligned with the plane z = 0 and the region of the box. Note that, although the large swirling
other with the plane x = 0. The initial tori setup is provided structure seen in Fig. 17(b) certainly has an impact on the
with the PENCIL CODE and described in the code’s manual; coherence of the field lines in the short term, it is not the main
variations of this problem with two or three interlocked rings convex coherent structure for a long integration of field lines
were previously studied by Kerr and Brandenburg [29] and with ξ = 2000.
Del Sordo et al. [30], respectively. We set η = ν = 0.003 and The coherence of the structures identified by IACD in
adopt a numerical resolution of 1283 grid points. Fig. 17(c) can be seen in Fig. 18, where we plot the magnetic
The tension due to the curvature of the flux tubes initially field lines (red) emanating from the two vortex contours
shrinks the rings and induces a velocity [29]. Eventually, (black). The left panels refer to the smaller vortex of Fig. 17(c)
magnetic reconnection takes place and the topology of the and the right panels to the larger vortex. The lower panels
field lines is changed. The configuration of the magnetic field are the same as the upper panels, but with the inclusion of
lines after t = 100 time units is shown in Fig. 16. The complex a few magnetic field lines from random initial conditions. It
entanglement plotted from random initial conditions hides the is remarkable that IACD is able to detect the magnetic flux
inherent order present in the magnetic coherent structures ropes among such a complex entanglement of field lines, espe-
that are still embedded in the field. We start searching for cially considering that the boundary was calculated only on a
those structures by considering their horizontal components plane.

(a) (b)

(c) (d)

FIG. 18. Magnetic flux-rope reconstruction at t = 100 with ξ = 2000: (a) magnetic field lines (dark gray [red online] lines) emanating
from the small vortex boundary (black line), (b) magnetic field lines from the large vortex boundary; (c, d) the same as in (a) and (b), but with
the inclusion of magnetic field lines emanating from random initial conditions (light gray [green online]).

043206-13
REMPEL, GOMES, SILVA, AND CHIAN PHYSICAL REVIEW E 99, 043206 (2019)

(a) has employed LCS techniques to detect structures in plasmas

from numerical simulations of the MHD equations [32–34] as
well as observational data of the solar photosphere [26,35,36].
In [14], the first definition of a magnetic vortex based on the
theory of LCS was introduced in the form of the IACD opera-
tor, which was explored in the present paper by means of three
MHD examples, namely, a two-dimensional incompressible
Orszag-Tang vortex, two-dimensional incompressible random
vortices, and three-dimensional compressible interlocked flux
rings. The IACD was shown to precisely detect the outer
contour of the coherent magnetic structures that keep the mag-
netic field lines bundled together, even among surrounding
turbulence. A comparison with detection using the current
density alone proved the superiority of the IACD, which is
also invariant to spatially dependent rotations and translations
of the observer, which the current density is not (see Sec. III).
The method can be used in conjunction with the squashing
(b)
factor Q [35–37] to investigate magnetic reconnection in
simulated and observational plasmas, since the IACD can
identify the main structures involved in the reconnection while
the Q operator detects the quasi-separatrix layers, which are
the main magnetic reconnection sites between those struc-
tures [12]. A possible application is the reconnection of the
solar emerging flux ropes with the coronal field [38]. Detec-
tion of the boundary layers in interacting magnetic flux ropes
is also crucial for understanding the origin of intermittent
turbulence in the solar wind, with resulting non-Gaussian
statistics and cross-scale coupling [39,40]. Additionally, flux
ropes are found in coronal mass ejections (CMEs), giant
releases of plasma from the Sun. These flux ropes are formed
in the low corona and carried away through the interplanetary
medium, where they can reconnect with magnetospheres [41].
The vortex merging observed in our two-dimensional ran-
FIG. 19. Magnetic flux-rope reconstruction at t = 100 with dom vortex simulations (Sec. IV A 2) was previously applied
ξ = 500: (a) the IACD field on the midplane z = 0, with one by Wu and Chang [7] to the observation of bursty bulk flows
magnetic vortex (black line), and (b) magnetic field lines (red lines) in the Earth’s magnetotail. In fact, sequential releases of tilted
emanating from the vortex boundary (black line). flux ropes due to reconnection were observed in the Earth’s
magnetotail by Hietala et al. [42], where the magnetic flux
As previously mentioned in Sec. II B 2, the parameter ξ in ropes were reconstructed by a combination of two-spacecraft
Eq. (11) basically controls the length and shape of the vortices timing for their global orientation and Grad-Shafranov re-
detected by IACD. In Fig. 19, we illustrate what happens construction (GSR) for their local orientation. The GSR
when a smaller value of ξ is adopted for the same magnetic method models the magnetic field structure traversed by a
field used in Figs. 17 and 18. The IACDξs0=500 spacecraft [43] and is the same technique employed in [39]
=0 (x0 ) field in
Fig. 19(a) is much simpler than in Fig. 17(a), since the field to reconstruct interplanetary magnetic flux ropes. Another
lines are short enough to stay away from the boundaries. A example of the use of the GSR method is the reconstruction
single wide vortex is found and a few magnetic field lines of flux ropes from ACE and Wind spacecraft data, which was
from its boundary are plotted in Fig. 19(b). The structure will recently used to study flux ropes merging in the solar wind by
lose coherence and the field lines will spread out if longer Zheng et al. [44]. This renders support for the use of IACD not
integration is conducted to compute the field lines. only in numerical simulations, but also in observational data.
Perhaps an even better natural laboratory for this technique
is the Sun, where the vector magnetic field can be measured
V. DISCUSSION AND CONCLUSIONS
by space instruments such as the Helioseismic and Mag-
In the past decade, objective (observer-independent) netic Imager (HMI) onboard the Solar Dynamics Observatory
Lagrangian techniques have prevailed in the detection of (SDO) [45]; alternatively the horizontal magnetic field in the
coherent structures in hydrodynamic flows by means of the solar surface can be estimated from the photospheric flow,
theory of Lagrangian coherent structures (LCS), with a set of as proposed by Yeates et al. [35]. In conclusion, we believe
mathematical tools being developed to precisely define hyper- that IACD is a novel and simple objective tool that can be
bolic (manifold lines), parabolic (jet cores), and elliptic (vor- used to detect magnetic vortices in numerical simulations and
tices) LCS [2,31]. Following that tendency, a series of works observations of space and astrophysical plasmas.

043206-14
OBJECTIVE MAGNETIC VORTEX DETECTION PHYSICAL REVIEW E 99, 043206 (2019)

ACKNOWLEDGMENTS 001. E.L.R. acknowledges Brazilian agencies CAPES, CNPq,

and FAPESP (Grants No. 2013/26258-4 and No. 16/24970-7)
We are indebted to Dr. C. C. Wu of UCLA for kindly pro-
for their financial support. S.S.A.S. and T.F.P.G. acknowledge
viding us his numerical codes for our 2D MHD simulations,
Brazilian agency CAPES for the financial support. Research
without which this paper would have not been possible. We
was carried out using the computational resources of the Cen-
thank Dr. F. J. Beron-Vera for helpful discussions. This study
ter for Mathematical Sciences Applied to Industry (CeMEAI),
was financed in part by the Coordenação de Aperfeiçoamento
funded by FAPESP (Grant No. 2013/07375-0).
de Pessoal de Nível Superior - Brasil (CAPES), Finance Code

[1] U. Frisch, Turbulence: The Legacy of A. N. Kolmogorov [23] C. Fang, An Introduction to Fluid Mechanics (Springer Interna-
(Cambridge University Press, Cambridge, UK, 1995). tional Publishing, Cham, 2019).
[2] G. Haller, A. Hadjighasem, M. Farazmand, and F. Huhn, [24] C.-C. Wu, Geophys. Astrophys. Fluid Dyn. 101, 37 (2007).
J. Fluid Mech. 795, 136 (2016). [25] J. Kim and P. Moin, J. Comput. Phys. 59, 308 (1985).
[3] R. Kinney, J. C. McWilliams, and T. Tajima, Phys. Plasmas 2, [26] S. S. A. Silva et al., Astrophys. J. 863, L2 (2018).
3623 (1995). [27] N. E. L. Haugen, A. Brandenburg, and W. Dobler, Phys. Rev. E
[4] S. A. Orszag and C.-M. Tang, J. Fluid Mech. 90, 129 (1979). 70, 016308 (2004).
[5] D. Biskamp and H. Welter, Phys. Fluids B 1, 1964 (1989). [28] http://pencil-code.nordita.org/
[6] H. Politano, A. Pouquet, and P. L. Sulem, Phys. Fluids B 1, [29] R. M. Kerr and A. Brandenburg, Phys. Rev. Lett. 83, 1155
2330 (1989). (1999).
[7] C.-C. Wu and T. Chang, J. Atmos. Terr. Phys. 63, 1447 (2001). [30] F. DelSordo, S. Candelaresi, and A. Brandenburg, Phys. Rev. E
[8] C. K. Armstrong and I. J. D. Craig, Solar Phys. 283, 463 81, 036401 (2010).
(2013). [31] G. Haller, Annu. Rev. Fluid Mech. 47, 137 (2015).
[9] C. K. Armstrong and I. J. D. Craig, Solar Phys. 289, 869 [32] E. L. Rempel, A. C.-L. Chian, and A. Brandenburg, Astrophys.
(2014). J. 735, L9 (2011).
[10] V. S. Lukin, Plasma Phys. Controlled Fusion 56, 060301 (2014). [33] E. L. Rempel, A. C.-L. Chian, and A. Brandenburg, Phys. Scr.
[11] J. Sears et al., Plasma Phys. Controlled Fusion 56, 095022 86, 018405 (2012).
(2014). [34] E. L. Rempel et al., J. Fluid Mech. 729, 309 (2013).
[12] J. B. O’Bryan and C. R. Sovinec, Plasma Phys. Controlled [35] A. R. Yeates, G. Hornig, and B. T. Welsch, Astron. Astrophys.
Fusion 56, 064005 (2014). 539, A1 (2012).
[13] C. Lawder and A. Yeates, Astrophys. J. 846, 106 (2017). [36] A. C.-L. Chian et al., Astrophys. J. 786, 51 (2014).
[14] E. L. Rempel et al., Mon. Not. R. Astron. Soc. 466, L108 [37] V. S. Titov, G. Hornig, and P. Démoulin, J. Geophys. Res. 107,
(2017). 1164 (2002).
[15] F. J. Beron-Vera, A. Hadjighasem, Q. Xia, M. J. Olascoaga, [38] D. MacTaggart and A. L. Haynes, Mon. Not. R. Astron. Soc.
and G. Haller, Proc. Natl. Acad. Sci. USA (2018), 43856, 1500 (2014).
doi:10.1073/pnas.1701392115. [39] A. C.-L. Chian et al., Astrophys. J. 832, 179 (2016).
[16] C. Prior and A. R. Yeates, Astrophys. J. 787, 100 (2014). [40] R. A. Miranda, A. B. Schelin, A. C.-L. Chian, and J. L. Ferreira,
[17] C. Truesdell and W. Noll, The Nonlinear Field Theories of Ann. Geophys. 36, 497 (2018).
Mechanics (Springer, Berlin, 2004). [41] A. Vourlidas, Plasma Phys. Controlled Fusion 56, 064001
[18] T. Günther and H. Theisel, Comput. Graph. Forum 37, 149 (2014).
(2018). [42] H. Hietala, J. P. Eastwood, and A. Isavnin, Plasma Phys.
[19] G. Haller, J. Fluid Mech. 525, 1 (2005). Controlled Fusion 56, 064011 (2014).
[20] S. Kotsiaros and N. Olsen, Int. J. Geomath. 3, 297 (2012). [43] Q. Hu, C. W. Smith, N. F. Ness, and R. M. Skoug, J. Geophys.

[21] QQ = QQ−1 = 1 ⇒ (QQ ) = Q Q + QQ  = 0 ⇒ Q Q = Res. 109, A03102 (2004).

−QQ . [44] J. Zheng, Q. Hu, Y. Chen, and J. le Roux, J. Phys.: Conf. Ser.
[22] ε123 = ε312 = ε231 = 1, ε213 = ε312 = ε321 = −1, and εi jk = 0 900, 012024 (2017).
for repeated indices. [45] J. T. Hoeksema et al., Solar Phys. 289, 3483 (2014).

043206-15