# Theory of the skyrmion, meron, anti-skyrmion and anti-meron in chiral
magnets

Sandip Bera1 and Sudhansu S. Mandal1,2 1Department of Physics, Indian
Institute of Technology, Kharagpur 721302, India
2Centre for Theoretical Studies, Indian Institute of Technology, Kharagpur
721302, India

###### Abstract

We find closed-form solution of the Euler equation for a chiral magnet in
terms of a skyrmion or a meron depending on the relative strengths of magnetic
anisotropy and magnetic field. We show that the relevant length scales for
these solutions primarily depend on the strengths of Dzyaloshinskii-Moriya
interaction through its ratios, respectively, with magnetic field and magnetic
anisotropy. We thus unambiguously determine the parameter dependencies on the
radius of the topological structures particularly of the skyrmions, showing an
excellent agreement with experiments and first-principle studies. An
anisotropic Dzyaloshinskii-Moriya interaction suitable for thin films made
with $C_{nv}$ symmetric materials is found to stabilize anti-skyrmion and
anti-meron, which are prototypical for $D_{2d}$ symmetric systems, depending
on the degree of anisotropy. Based on these solutions, we obtain phase diagram
by comparing the energies of various collinear and non-collinear competing
phases.

The chiral Dzyaloshinskii-Moriya interaction (DMI) D57 ; M60  for broken
inversion symmetric systems is one of the most important mechanisms including
frustrated exchange interactions and long-ranged dipolar interaction for
producing one-dimensional modulation in magnetization known as spin-spiral
Bogdanov89 ; Bogdanov06 ; Meckler09 ; Ezawa ; Leonov16 ; Zhang17 ; Okubo12  in
ferromagnetic systems. An application of magnetic field in such a system
stabilizes Bogdanov94 ; Leonov16  skyrmions (Sks) having topologically
protected quasiparticle-like spin structure Muhlbauer09 ; Tokura10a ; Heinze11
; Tokura12 ; Romming13 ; Tanigaki15  in a ferromagnetic background. Néel and
Bloch type Sks are generally realized respectively in $C_{nv}$ and $D_{n}$
symmetric Bogdanov89 ; Bogdanov94  bulk and thin-film materials for wide range
of magnetic fields and temperatures Muhlbauer09 ; Tokura10a ; Heinze11 ;
Tokura12 ; Romming13 ; Tanigaki15 ; Kezsmarki15 ; Fujima17 ; Tokura10 ; Chacon
; Das19 .

Recent observations Nayak17 ; Nayak19  of anti-skyrmions (ASks) in Heusler
alloys with $D_{2d}$ crystal symmetry have raised an issue about the
microscopic environment which will stabilize a Sk or an ASk. While a Sk has
either Néel or Bloch type of orientation of magnetization vector governed by
the respective transverse and longitudinal DMI, an ASk displays a combination
of both. It is thus tempting to think that an anti-skyrmion may be produced in
a crystal whose symmetry gives rise to both types of DMI. Numerical
simulations, on the contrary, Koshibae14 ; Koshibae16 ; Camosi18  indicate
that the ASks do stabilize only in the presence of dipolar interaction. A
micromagnetic study Hoffmann17  suggests that Sks and ASks can, however,
coexist and this coexistence is predicted by electronic structure calculation
at interfaces due to anisotropic DMI. These ASks even take part in current-
induced motion Huang17 .

Hoffmann et al Hoffmann17  have recently observed ASks in $C_{2v}$ symmetric
systems grown on semiconductor or heavy-metal substrates, while $C_{nv}$
symmetric systems are known to stabilize Sks only Bogdanov89 ; Leonov16 . This
motivates us to study a system of thin film chiral magnet that may be
fabricated with $C_{nv}$ symmetric crystals with an anisotropic DMI in a
continuum model in search of ASk solution. Camosi et al Camosi17  have
recently reported that the epitaxially grown thin Co films on W(110) brings
anisotropy in DMI along two orthogonal growth directions of a $C_{2v}$
symmetric bulk system. Although this reported anisotropy does not correspond
to two opposite signs along two orthogonal directions, a micromagnetic
simulation seems to suggest anisotropy in thin films not only in magnitude but
also in sign Huang17 . We introduce a model DMI with such an anisotropy.

Moreover, recent observation of another topological spin structure, viz, meron
Nych17 ; Tokura18  have further raised the theoretical issue on the parameter
regimes on which all these different kinds of topological structures emerge.
Further, definite parameter dependencies on the radius Romming15 ; Malottki ;
Thiaville ; Wang  and appropriate length scale Leonov16  of a SK are not yet
settled. Our focus is thus solving basic Euler equation for angular variables
representing magnetization with isotropic DMI for Sks and merons and then
study the consequences of its anisotropy followed by the determination of
phase diagram by comparing energies of different possible solutions for
thermodynamically stable magnetic structures.

In this letter, we solve the Euler equation in a continuum model Bogdanov89 ;
Leonov16  with ferromagnetic exchange coupling, $J$, DMI strength, $D$,
strength of magnetic anisotropy, $A$, and net Zeeman energy due to magnetic
field, $H$. For moderate to high $HJ/D^{2}$ and $\gamma=2A/H<1$, we find that
the relevant length scale of the corresponding skyrmion solution is
$r_{s}=D/H$, contrasting the belief NT  of the relevant length scale
$r_{d}=J/D$. This enables us to determine the magnetic field and anisotropy
dependencies of the radius of a skyrmion and find that it is in excellent
agreement with experiments Romming15 ; Vousden16  and first-principle studies
Malottki . The meron solution at zero magnetic field is obtained for $A>0$
(easy-plane anisotropy) by minimizing energy and the relevant length scale is
found to be $D/A$. We show the formation of meron lattice and argue how a
symmetric Sk is evolved from a meron via an asymmetric skyrmion, explaining a
recent experiment as well as simulation resultsTokura18 . Further, our model
with an anisotropic DMI is shown to stabilize ASks and antimerons in $C_{nv}$
symmetric systems, as evident in recent realization Hoffmann17  of ASks in
$C_{2v}$ symmetric systems. We finally determine phase diagram for $\gamma<1$
by comparing energies of the skyrmion solution with other collinear and non-
collinear competing phases.

We begin with considering a two-dimensional chiral magnet having energy
$E=\int d^{2}r({\cal E}_{{\rm EX}}+{\cal E}^{\pm}_{{\rm DM}}+{\cal E}_{{\rm
AH}})$ with respect to an overall ferromagnet orienting along perpendicular to
the plane of the system, described by exchange energy density ${\cal E}_{{\rm
EX}}=\frac{J}{2}\left(\nabla\hat{m}\right)^{2}$, DMI energy density ${\cal
E}^{\pm}_{{\rm DM}}=-D\left(L_{xz}^{(x)}\pm L_{yz}^{(y)}\right)$. Here
$L_{ij}^{(k)}=\hat{m}_{i}\partial_{x_{k}}\hat{m}_{j}-\hat{m}_{j}\partial_{x_{k}}\hat{m}_{i}$,
$\hat{m}$ is unit magnetization vector, $\pm$ signs, respectively, refer to
the systems with $C_{nv}$ and $D_{2d}$ symmetries when the Dzyaloshinskii-
Moriya vector is transverse fnote1  to the lattice-bond. (While the former
supports SKs the later is suitable for stabilizing ASks.) The energy density
for magnetic anisotropy and applied magnetic field along $\hat{z}$ direction
given by ${\cal E}_{{\rm AH}}=-A\,(1-\hat{m}_{z}^{2})+H\,(1-\hat{m}_{z})$,
where $A>0\,(<0)$ refers to easy-plane (easy-axis) anisotropy. In spherical
polar representation,

$\hat{m}({\bm{r}})=\left[\cos\Phi(\bm{r})\sin\Theta(\bm{r}),\,\sin\Phi(\bm{r})\sin\Theta(\bm{r}),\,\cos\Theta(\bm{r})\right]$

with $\bm{r}=(r\cos\phi,r\sin\phi)$ in polar coordinate system.

A topological structure defined by its topological quantum number $N_{{\rm
sk}}=\frac{1}{4\pi}\int
d^{2}r\,\hat{m}\cdot\left(\partial_{x}\hat{m}\times\partial_{y}\hat{m}\right)\equiv\pm
N_{0}$, where positive sign refers to a Sk or meron and negative sign refers
to an ASk or anti-meron. The solutions Leonov16  of a Sk/meron and an
ASk/anti-meron correspond to $\Theta(\bm{r})=\Theta(r)$ and, respectively,
$\Phi(\bm{r})=\pm\phi+\eta$. Here $\eta$ determines a constant extra planar
rotation of magnetic moment at all points; $\eta=0\,(\pi/2)$ for Neél(Bloch)
type topological structures. Here $N_{0}$ represents the winding number NT :
its positive (negative) sign determines inward (outward) spin orientation with
respect the origin, corresponding to negative (positive) sign of $D$, and its
magnitude is $1$ for Sks and ASks, and $1/2$ for merons and anti-merons. The
boundary condition, $\hat{m}=(0,0,+1)$ for $r\to\infty$, i.e.,
$\Theta(r\to\infty)=0$ and $\hat{m}=(0,0,-1)$ at $r=0$, i.e.,
$\Theta(r=0)=\pi$ is for both SK and aSK. Meron and anti-meron correspond to
the boundary condition $\Theta(r=\infty)=\pi/2$ and $\Theta(r=0)=\pi\,(0)$ for
inward (outward) helicity.

Figure 1: (color online) Skyrmion or anti-Skyrmion solution: Numerical
solution of the Euler equation (1), i.e., $\Theta(r)$ vs. $r/r_{s}$ for
different values of $\gamma$ in the panels (a)–(e) for same set of $H_{0}/H$,
viz, 0.36, 0.49, 0.64 and 0.81. Solid lines represent the best fit solution
with the form given by Eq. (2). Inset in the panel (c) is a representative of
the panels (a)–(e) for showing the crossing of the curves for different values
of $H_{0}/H$. Also see supplemental material Suppli  for wider range of
parameters. This crossing point is identified as the radius, $R_{s}$, of a
skyrmion. The dependence of $R_{s}$ on $\gamma$ along with its approximate
fitted form is shown in the panel (f).

No matter, be it $C_{nv}$, $D_{2d}$ or $D_{n}$ systems, the Euler equation for
$\Theta(r)$ is identical Suppli . By introducing a length scale $r_{s}=D/H$
and rescaling $r\to r_{s}\rho$, we obtain Suppli  the Euler equation

$\displaystyle\frac{d^{2}\Theta}{d\rho^{2}}+\frac{1}{\rho}\frac{d\Theta}{d\rho}-\frac{\sin\Theta\cos\Theta}{\rho^{2}}$
$\displaystyle=\frac{H_{0}}{H}\left(-\frac{2}{\rho}\sin^{2}\Theta+\sin\Theta-\gamma\sin\Theta\cos\Theta\right)$
(1)

where $H_{0}=D^{2}/J$ and $\gamma=2A/H$. The numerical solutions of the Eq.
(1) with the boundary conditions $\Theta(0)=\pi$ and $\Theta(\infty)=0$ for
different values of $H_{0}/H$ and $\gamma$ are shown in Fig. 1. The length
scale $r_{s}$ which is independent of exchange energy $J$ defines the relevant
length scale as for a fixed value of $\gamma$, the deviation of the curves of
$\Theta(r)$ for different values of $H_{0}/H$ are almost negligible; the
complete solution of Eq.(1) is thus best approximated by

$\Theta(r)=4\arctan\left(\exp\left(-\beta(\gamma)r/r_{s}\right)\right)$ (2)

with $\beta(\gamma)\approx 0.91-0.55\gamma$. We note that all the curves for a
fixed $\gamma$ cross (see inset of Fig. 1(c)) at a particular $r$ and we
identify that to be the radius, $R_{s}$, of a Sk. We find its dependency on
$\gamma$ as $R_{s}=r_{s}w(\gamma)$ with $w(\gamma)\approx
0.26+\frac{2.09}{1-0.36\gamma}$. Therefore, the magnetic field dependence of
the radius of a Sk may be parametrized as

$R_{s}=\frac{C_{1}}{H}+\frac{C_{2}}{H-C_{3}}$ (3)

where the coefficients $C_{1}$ and $C_{2}$ are proportional to $|D|$ and
$C_{3}$ is proportional to $A$. We note that for a fixed $H$, radius of a Sk
increases with positive $A$, in agreement with an experiment Vousden16 .
However, an increase of easy-axis anisotropy will reduce the size of an Sk.
Figure 2 shows that the skyrmion radius obtained in an experiment Romming15
and first-principle studies Malottki  obey the relation (3) very well and the
sign of the corresponding fitted $C_{3}$ are consistent with the sign of the
reported $A$. For the systems with positive $A$, lower bound of the magnetic
field needed for producing a Sk is $H_{{\rm lb}}=2A$ and thereafter the radius
monotonically decreases with increasing $H$.

Figure 2: (color online) Variation of radius of a Sk with applied magnetic
field: Experimental data (solid circles) from Ref. Romming15, , first-
principle calculations from Ref. Malottki,  for hcp lattice (open circle) and
fcc lattice (open inverted triangles), denoted respectively as I, II and III.
Solid lines are the fitted curves of these data with the functional form in
Eq. 3. Inset: coefficients obtained by fitting are tabulated where we have put
primes to distinguish from non-primes in Eq.(3) in view of dimensions.

Figure 3(a) shows phase diagram in $A$–$H$ space with $\gamma<1$. The phase
boundary between skyrmion and the polarized ferromagnet is determined by
comparing energy of a Sk,

$\displaystyle E_{{\rm
sk}}=2\pi\int_{0}^{\infty}dr\,r\left[\frac{J}{2}\left(\left(\frac{d\Theta}{dr}\right)^{2}+\frac{\sin^{2}\Theta}{r^{2}}\right)-A\sin^{2}\Theta\right.$
$\displaystyle\hskip
8.5359pt\left.+D\left(\frac{d\Theta}{dr}+\frac{\sin\Theta\cos\Theta}{r}\right)+H(1-\cos\Theta)\right]$
(4)

with the energy of the ferromagnet. Similarly by determining energy of a spin-
spiral following Ref. Bogdanov89,  in comparison to the ferromagnet, we obtain
the phase-boundary between the spin-spiral and ferromagnet. We draw phase
boundary between spin-spiral and skyrmion phases by considering maximum
possible phase-space for spin-spiral structure. The phase diagram for
$\gamma<1$ here is consistent with previously reported phase diagrams obtained
by variational and other simulations Banerjee14 ; Rowland16 . In agreement
with an experiment Herve18 , both spin-spiral and skyrmions are accessible at
zero anisotropy.

Figure 3: (color online) (a) Phase Diagram in $A/A_{0}$–$H/H_{0}$ space, where
$A_{0}=H_{0}=D^{2}/J$. Three distinct phases: Spin-spiral (SS), isolated
skyrmion (SK), and polarized ferromagnet (FM) whose magnetization direction is
along the applied $H$. The right boundary corresponds to $\gamma=1$. (b) Two
possible degenerate structures of merons (up or down spin at the center) for
$H=0$. From top, two and one-dimensional spin-structures of merons. Any finite
$H$ stabilizes meron with down-spin at the center only. (c) From top,
depiction of a symmetric meron ($H=0$), an asymmetric meron ($H>0$,
$\gamma>1$), and an asymmetric skyrmion. Spin-down at $r=0$ for all these
structures. Here the schematics (b,c) of merons and skyrmion are considered
for $D<0$.

For a sufficiently high easy-plane anisotropy $(A>0)$ and $H=0$, all the spins
will align in the plane (planar ferromagnet). This indicates a boundary
condition $\Theta(\infty)=\pi/2$ which together with another boundary
condition $\Theta(0)=0$ or $\pi$ will provide a solution of meron when $A$ is
moderate. Taking cue of the skyrmion solution, we assume the solutions of
meron Suppli  to be

$\Theta(r)=\pm\frac{\pi}{2}+2\arctan(\exp(-\zeta r/r_{a}))$ (5)

where $r_{a}=D/A$ is the characteristic length scale, positive (negative) sign
corresponds to spin down (up) at the center of the meron, and the parameter
$\zeta$ to be determined by minimizing its energy

$\displaystyle E_{{\rm meron}}$ $\displaystyle=$ $\displaystyle
2\pi\int_{0}^{\infty}r\,dr\left[\frac{J}{2}\left(\left(\frac{d\Theta}{dr}\right)^{2}+\frac{\sin^{2}\Theta}{r^{2}}\right)\right.$
(6)
$\displaystyle+A\cos^{2}\Theta+\left.D\left(\frac{d\Theta}{dr}+\frac{\sin\Theta\cos\Theta}{r}\right)\right]\,.$

We find Suppli  $\zeta=2\ln(2)/(1+2G)\approx 0.49$, where $G$ is Catalan’s
constant. These solutions of $\Theta(r)$ are degenerate and hence they occur
simultaneously and appear as neighbors to match the background of planar
ferromagnet and form a meron-lattice, as shown in Fig. 3(b). However, with the
increase of $H$, only one-kind of meron (spin-down at its core) survive as the
other will have higher energy, because the background of spin alignment will
have nonzero out of plane (up) component. For further increase of $H$, this
meron gradually converts into a skyrmion as it helps to orient more spin with
finite up component. This is in reminiscent of the recently observed merons by
Yu et al. Tokura18 . We estimate the upper-bound of $A$ for forming a meron as
$A_{{\rm ub}}\approx 2.3A_{0}$ by comparing the energy of a meron and the
planar ferromagnet.

In the presence of $H$ with $A>>0$ such that $\gamma>1$, a tilted ferromagnet
will be formed with finite amount of spin-projection along the direction of
$H$ making a tilting angle $\arccos(1/\gamma)$ with the plane. With such a
tilted ferromagnet in the background, locally formed merons with down-spin at
its core will be asymmetric as shown schematically in Fig. 3(c) when
$A<A_{{\rm ub}}$. If we look along a particular direction, a meron’s spin
alignment at one boundary will be along the tilting angle $\arccos(1/\gamma)$
and at the other boundary it will differ by an angle $\pi$. This makes the
meron asymmetric. We note that actual $A_{{\rm ub}}$ may be lower than
estimated here because of the predicted possibility of forming cone-like
structure in the intermediate regime. The cone structure Rowland16 ;
Banerjee14  and the tilted-ferromagnet are indistinguishable in our analysis
because both these structures correspond to same $\Theta$. With further
increase of $H$, some more spins will tend to align more than
$\arccos(1/\gamma)$ forming an asymmetric skyrmion (Fig. 3(c)), corroborated
with the recent numerical simulation result Leonov17 . Upon further increase
of $H$, right(left) side of the Sk becomes shorter(longer) and evolve into a
symmetric Sk at $\gamma=1$ as we enter into the Sk phase of the phase-diagram
(Fig.3(a)).

In search of ASk and anti-meron in thin films made of $C_{nv}$ symmetric
sytems fnote2 , we introduce an anisotropic DMI given by

${\cal E}_{{\rm DM}}^{+}=-D(1+\lambda\cos
2\phi)\left(L_{xz}^{(x)}+L_{yz}^{(y)}\right)\,.$ (7)

Here $\lambda$ denotes the degree of anisotropy with $\lambda=0$ representing
the symmetric DMI present in the bulk $C_{nv}$ symmetric crystals. The energy
of an ASk is then found to be

$\displaystyle E_{{\rm
ask}}=2\pi\int_{0}^{\infty}dr\,r\left[\frac{J}{2}\left(\left(\frac{d\Theta}{dr}\right)^{2}+\frac{\sin^{2}\Theta}{r^{2}}\right)-A\sin^{2}\Theta\right.$
$\displaystyle\hskip
8.5359pt\left.+\frac{\lambda}{2}D\left(\frac{d\Theta}{dr}-\frac{\sin\Theta\cos\Theta}{r}\right)+H(1-\cos\Theta)\right]$
(8)

with $\Theta(r)$ given in Eq. (2). Inset of Fig. 4(a) shows the variation of
$E_{{\rm ask}}$ with $\lambda$ for $A=0$ and we find that $E_{{\rm
ask}}<E_{{\rm sk}}$ above a critical value $\lambda_{c}\approx 1.4$ and hence
the anisotropy in DMI stabilizes an ASk. A phase diagram has been presented in
Fig. 4(a). Ferromagnet to ASk transition is also possible for
$\lambda>\lambda_{c}$, and the corresponding critical value increases with
$H$. However, ASks are not possible for lower $H/H_{0}$ where spin-spiral
phase remain unaltered for any $\lambda$. Figure 4(b) shows minimum values of
$\lambda$ above which the full phase-space of Sks and partial-phase of
ferromagnets shown in Fig. 3(a) can stabilize ASks. The outer boundary in the
ferromagnetic region is obtained with the criterion that the ratio of the
diameter of an ASk and the spin-spiral wavelength is not less than 0.4.

The energy of an anti-meron in presence of anisotropic DMI,

$\displaystyle E_{{\rm anti-
meron}}=2\pi\int_{0}^{\infty}r\,dr\left[\frac{J}{2}\left(\left(\frac{d\Theta}{dr}\right)^{2}+\frac{\sin^{2}\Theta}{r^{2}}\right)\right.$
$\displaystyle+A\cos^{2}\Theta+\left.\frac{\lambda}{2}D\left(\frac{d\Theta}{dr}-\frac{\sin\Theta\cos\Theta}{r}\right)\right]\,.$
(9)

becomes less than $E_{{\rm meron}}$ for $6.8\lesssim\lambda$. Producing anti-
meron by anisotropic DMI is less probable than producing ASk because the
former requires much higher degree of anisotropy, which is almost in the verge
of the limit of a $D_{2d}$ system.

We here have shown that the anisotropic DMI in thin films with $C_{nv}$
symmetric materials can host anti-skyrmions for wide range of phase-space of
$A$ and $H$, in comparison to hosting skyrmions. However, we do not find any
regime of the coexistence of Sks and ASks, in contrary to the numerical
simulation Hoffmann17 . Although dipolar interaction is also a suitable
mechanism Koshibae14 ; Koshibae16 ; Camosi18 for stabilizing ASks, the
anisotropic DMI is solely responsible, to the best of our knowledge, for
small-size ASks in $C_{nv}$ symmetric systems. The dipolar interaction here
may play a role in reducing Lobanov16  the effect of magnetic anisotropy. The
physics of Sk/ASk and meron/anti-meron discussed here will reverse for systems
with $D_{2d}$ symmetries. Although the structure of an anti-skyrmion is a
combination of the structures of Néel and Bloch type Sks which are
prototypical, respectively, of DMI with Dzyaloshinskii-Moriya vector
orthogonal to the neighboring bond and along the bond, their combinations do
not produce ASks. However, a pure $D_{n}$ symmetric system will stabilize
Bloch type merons and SKs, and the corresponding anti-merons and ASks may also
be produced through anisotropic DMI.

S.S.M. is supported by SRIC, IIT Kharagpur through the project code EFH.

Figure 4: (color online) (a) Phase diagram in $H/H_{0}$–$\lambda$ plane when
$A=0$. Four distinct phases are obtained by comparing their energies. While
the energies of a skyrmion and an anti-skyrmion have been calculated using the
expressions (4) and (8), respectively, the energy of spin-spiral phase has
been obtained following the procedure reported in Ref. Bogdanov89, . The inset
shows comparison of the energies of a Sk and an ASk with respect to the energy
of the ferromagnetic ground state in the back-ground as a function of
anisotropy parameter $\lambda$ of DMI. The crossing point of two lines
determines the boundary between skyrmionic and anti-skyrmionic phases. (b) The
color map indicates minimum value of $\lambda$ above which an ASk is
stabilized over a Sk or a ferromagnet in the parameter space shown in Fig.
3(a).

## References

  * (1) I. E. Dzyaloshinskii, Zh. Eksp. Theo. Fiz. 32, 1547 (1957)[Sov. Phys. JETP 5, 1259 (1957)].
  * (2) T. Moriya, Phys. Rev. 120, 91 (1960).
  * (3) A. N. Bogdanov and D. A. Yabionskii, Sov. Phys. JETP 68, 101 (1989).
  * (4) U. K. Rößler, A. N. Bogdanov, and D. Pfleiderer, Nature 442, 797 (2006).
  * (5) S. Meckler, N. Mikuszeit, A. Preßler, E. Y. Vedmedenko, O. Pietzsch, and R. Wiesendanger, Phys. Rev. Lett. 103, 157201 (2009).
  * (6) M. Ezawa, Phys. Rev. Lett. 105, 197202 (2010).
  * (7) A. O. Leonov, T. L. Monchesky, N. Romming, A. Kubetzka, A. N. Bogdanov, and R. Wiesendanger, New J. Physics 18, 065003 (2016).
  * (8) X. Zhang, J. Xia, Y. Zhou, X. Liu, H. Zhang and M. Ezawa, Nature Comm. 8, 1717 (2017)
  * (9) T. Okubo, S. Chung, and H. Kawamura, Phys. Rev. Lett. 108, 017206 (2012).
  * (10) A. Bogdanov and A. Hubert, J. Mag. Mag. Mat. 138, 255 (1994).
  * (11) S. Mühlbauer, B. Binz, F. Jonietz, C. Pfleiderer, A. Rosch, A. Neubauer, R. Georgii, and P. Böni, Science 323, 915 (2009).
  * (12) X. Z. Yu, N. Kanazawa, Y. Onose, K. Kimoto, W. Z. Zhang, S. Ishiwata, Y Matsui and Y. Tokura, Nature Materials. 10, 106 (2011).
  * (13) S. Heinze, K. von Bergmann, M. Menzel, J. Brede, A. Kubetzka, R. Wiesendanger, G. Bihlmayer, and S. Blügel, Nature Phys. 7, 713 (2011).
  * (14) A. Tonomura, X. Yu, K. Yanagisawa, T. Matsuda, Y. Onose, N. Kanazawa, H. S. Park, and Y. Tokura, Nano Lett. 12, 1673 (2012).
  * (15) N. Romming, C. Hanneken, M. Menzel, J. E. Bickel, B. Wolter, K. von Bergmann, A. Kubetzka, and R. Weisendamger, Science 341 636 (2013).
  * (16) T. Tanigaki, K. Shibata, N. Kanazawa, X. Yu, Y. Onose, H. S. Park, D. Shindo, and Y. Tokura, Nano Lett. 15, 5438 (2015).
  * (17) I. Kézsmárki, S. Bordács, P. Milde, E. Neuber, L. M. Eng, J. S. White, H. M. Ronnow, C. D. Dewhurst, M. Mochizuki, K. Yanai, H. Nakamura, D. Ehlers, V. Tsurkan and A. Loidl , Nature Materials. 14, 1116 (2015).
  * (18) Y. Fujima, N. Abe, Y. Tokunaga, and T. Arima, Phys. Rev. B 95, 180410(R) (2017).
  * (19) X. Z. Yu, Y. Onose, N. Kanazawa, J. H. Park, J. H. Han, Y. Matsui, N. Nagaosa, and Y. Tokura, Nature 465, 901 (2010).
  * (20) A. Chacon, L. Heinen, M. Halder, A. Bauer, W. Simeth, S. Mühlbauer, H. Berger, M. Garst, A. Rosch and C. Pfleiderer, Nature Phys. 14, 936 (2018).
  * (21) S. Das et al., Nature 568, 368 (2019).
  * (22) A. K. Nayak, V. Kumar, T. Ma, P. Werner, E. Pippel, R. Sahoo, F. Damay, U. K. Rößler, C. Felser, S. S. P. Parkin, Nature 548, 561 (2017).
  * (23) S. Jamaluddin, S. K. Manna, B. Giri, P. V. P. Madduri, S. S. P. Parkin, and A. K. Nayak, Adv. Funct. Mater. 0, 1901776 (2019).
  * (24) W. Koshibae and N. Nagaosa, Nature Comm. 5, 5148 (2014).
  * (25) W. Koshibae and N. Nagaosa, Nature Comm. 7, 10542 (2016).
  * (26) L. Camosi, N. Rougemaille, O. Fruchart, J. Vogel, and S. Rohart, Phys. Rev. B 97, 134404 (2018).
  * (27) M. Hoffmann, B. Zimmermann, G. P. Müler, D. Schürhoff, N. S. Kiselev, C. Melcher, and S. Blügel, Nature Comm. 8, 308 (2017).
  * (28) S. Huang, C. Zhou, G. Chen, H. Shen, A. K. Schmid, K. Liu, and Y. Wu, Phys. Rev. B 96, 144412 (2017).
  * (29) L. Camosi, S. Rohart, O. Fruchart, S. Pizzini, M. Belmeguenai, Y. Roussinné, A. Stashkevich, S. M. Cherif, L. Ranno, M. Santis, J. Vogel, Phys. Rev. B 95 , 214422 (2017).
  * (30) A. Nych, Jun-ichi Fukuda, U. Ognysta, S. Zumer and I. Musevic, Nature Phys. 13, 1215 (2017).
  * (31) X. Z. Yu, W. Koshibae, Y. Tokunaga, K. Shibata, Y. Taguchi, N. Nagaosa and Y. Tokura, Nature 564, 95 (2018).
  * (32) N. Romming, A. Kubetzka, C. Hanneken, K. von Bergmann, and R. Wiesendanger, Phys. Rev. Lett. 114, 177203 (2015).
  * (33) S. von Malottki, B. Dupé, P. F. Bessarab, A. Delin and S. Heinze, Sci. Rep. 7, 12299 (2017)
  * (34) S. Rohart and A. Thiaville, Phys. Rev. B 88, 184422 (2013).
  * (35) X. S. Wang, H. Y. Yuan and X. R. Wang, Comm. Phys. 1, 31 (2018).
  * (36) N. Nagaosa and Y. Tokura, Nat. Nanotechnol. 8, 899 (2013).
  * (37) M. Vousden, M. Albert, M. Beg, Marc-Antonio Bisotti, R. Carey, D. Chernyshenko, D. Cortes-Ortuno, W. Wang, O. Hovorka, C. H. Marrows and H. Fangohr, Appl. Phys. Lett. 108, 132406 (2016).
  * (38) DMI suitable for $D_{n}$ symmetric crystals, ${\cal E}_{{\rm DM}}=-D\left(L_{zx}^{(y)}+L_{yz}^{(x)}\right)$, leads to a Bloch type skyrmionic solution with $\eta=\pi/2$.
  * (39) See supplementary material for this manuscript.
  * (40) S. Banerjee, J. Rowland,O. Erten and M. Randeria, Phys. Rev. X 4, 031045 (2014).
  * (41) J. Rowland, S. Banerjee, and M. Randeria, Phys. Rev. B 93, 020404(R) (2016).
  * (42) M. Hervé, B. Dupé, R. Lopes, M. Böttcher, M. D. Martins, T. Balashov, L. Gerhard, J. Sinova and W. Wulfhekel, Nature Comm. 9, 1015 (2018).
  * (43) A. O. Leonov and I. Kézsmárki, Phys. Rev. B 96, 014423 (2017).
  * (44) A similar anisotropy in DMI corresponding to thin films fabricated with $D_{2d}$ symmetric crystals given by energy density ${\cal E}_{{\rm DM}}^{-}=-D(1+\lambda\cos 2\phi)\left(L_{xz}^{(x)}-L_{yz}^{(y)}\right)$ will produce similar results with reversing roles of Sks and ASks.
  * (45) I. S. Lobanov, H. Jónsson and V. M. Uzdin, Phys. Rev. B 94, 174418 (2016). Supplemental Material for “Theory of skyrmion, meron, anti-skyrmion and anti-meron in chiral magnet” Euler equation and its skyrmion solution:
Energy for a two dimensional chiral magnet is given by  $E=\int
d^{2}\bm{r}\left[\frac{J}{2}\left(\nabla\hat{m}\right)^{2}+{\cal
E}_{DM}+H(1-\hat{m}_{z})-A(1-\hat{m}_{z}^{2})\right]$ (10) where energy
density for DMI, ${\cal E}_{DM}=-D(L_{xz}^{(x)}+L_{yz}^{(y)})$ for $C_{nv}$
symmetric, $-D(L_{yz}^{(x)}+L_{zx}^{(y)})$ for $D_{n}$ symmetric and
$-D(L_{xz}^{(x)}-L_{yz}^{(y)})$ for $D_{2d}$ symmetric systems and we
distinguish them by introducing a parameter $\alpha$ with respective values 1,
2 and $-3$. Unit magnetization can be parametrized through spherical variables
$(\Theta,\Phi)$ as
$\hat{m}({\bm{r}})=\left[\cos\Phi(\bm{r})\sin\Theta(\bm{r}),\,\sin\Phi(\bm{r})\sin\Theta(\bm{r}),\,\cos\Theta(\bm{r})\right]$
with $\bm{r}=(r\cos\phi,\,r\sin\phi)$ in polar coordinate. Therefore, the
expression of energy reduces to  $\displaystyle E$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}r\,dr\int_{0}^{2\pi}d\phi\left[\frac{J}{2}\left(\Theta_{r}^{2}+\frac{\sin^{2}\Theta}{r^{2}}\Phi_{\phi}^{2}\right)\right.$
(11) $\displaystyle+$ $\displaystyle H(1-\cos\Theta)-A\sin^{2}\Theta$
$\displaystyle+$
$\displaystyle\left.D\left(\Theta_{r}\pm\frac{\sin(2\Theta)}{2r}\Phi_{\phi}\right)\sin(\alpha\frac{\pi}{2}-\phi\pm\Phi)\right]$
where $\Theta_{r}=\frac{d\Theta(r)}{dr}$ and
$\Phi_{\phi}=\frac{d\Phi(\phi)}{d\phi}$ and assuming
$\Theta(\bm{r})=\Theta(r)$ and $\Phi(\bm{r})=\Phi(\phi)$. In the last term
(11), positive sign refers to $\alpha=1$ and $2$ and negative sign refers to
$\alpha=-3$. Considering $\Phi(\phi)=\zeta\phi+\eta$ with $\zeta=1$ for
$\alpha=1,\,2$ and $\zeta=-1$ for $\alpha=-3$, and $\eta=0$ for
$\alpha=1,\,-3$ and $\pi/2$ for $\alpha=2$, we find Euler equation of
$\Theta(r)$ being independent on $\alpha$ as  $\displaystyle
J\left(\Theta_{rr}+\frac{\Theta_{r}}{r}-\frac{\sin\Theta\cos\Theta}{r^{2}}\right)+\frac{2D}{r}\sin^{2}\Theta$
$\displaystyle=H\sin\Theta-A\sin(2\Theta)$ (12) where
$\Theta_{rr}=\frac{d^{2}\Theta}{dr^{2}}$. The corresponding energy is given by
$\displaystyle E_{{\rm sk}}$ $\displaystyle=$ $\displaystyle
2\pi\int_{0}^{\infty}r\,dr\left[\frac{J}{2}\left(\Theta_{r}^{2}+\frac{\sin^{2}\Theta}{r^{2}}\right)-A\sin^{2}\Theta\right.$
(13)
$\displaystyle+H(1-\cos\Theta)+\left.D\left(\Theta_{r}+\frac{\sin(2\Theta)}{2r}\right)\right]$
Short-distance singularity in exchange energy may be avoided with the boundary
condition $\Theta(r=0)=0$ or $\pi$. We look for skyrmion solution of Eq. (12)
by introducing another boundary condition $\Theta(r=\infty)=0$ (assuming the
background as a polarized ferromagnet) along with $\Theta(0)=\pi$. In the
absence of anisotropy $(A=0)$, an approximate and asymptotically $(r\to
0,\,\infty$) exact analytical solution of Eq. (1) may be obtained as the exact
solution of the simple sine-Gordon like equation
$J\frac{d^{2}\Theta}{dr^{2}}=H\sin\Theta$, i.e.,
$\Theta(r)=4\arctan\left(\exp\left(-r/r_{0}\right)\right)$ with characteristic
length scale $r_{0}=\sqrt{J/H}$. Therefore, an approximate (exact for $r\to 0$
and $\infty$) solution for $A\neq 0$ may be obtained by considering a reduced
form of Eq. (1) as $J\frac{d^{2}\Theta}{dr^{2}}=H\sin\Theta-A\sin(2\Theta)$
whose solution satisfies an integral equation
$\int\frac{d\Theta}{\sin(\Theta/2)\sqrt{1-\gamma\cos^{2}(\Theta/2)}}=-2(r/r_{{}_{0}})$
(14) with $\gamma=2A/H\,(\gamma<1)$, expressible into an algebraic equation
$\displaystyle\sqrt{1+\frac{2\gamma\cos(\Theta/2)}{1-\gamma\cos(\Theta/2)+\sqrt{1-\gamma}\sqrt{1-\gamma\cos^{2}(\Theta/2)}}}$
$\displaystyle\times\tan(\Theta/4)\hskip
28.45274pt=\exp\left(-\sqrt{1-\gamma}\,\,r/r_{0}\right)\,.$ (15) We find the
solution of Eq.(15) as
$\Theta(r)=4\arctan\left(\exp\left(-g(\gamma)\,r/r_{0}\right)\right)$ (16)
with $g(\gamma)\simeq 1-\frac{\gamma}{7}-\frac{\gamma^{2}}{30}$. However, as
the smooth change in the orientation of spin depends on $D$, it is natural
that we consider another length scale $r_{d}=J/D$ which is the appropriate
length scale for spin-spirals. By introducing $r_{d}$ and transforming $r\to
r_{d}\rho$, we find the reduced form of Eq. (12) as
$\displaystyle\frac{d^{2}\Theta}{d\rho^{2}}+\frac{1}{\rho}\frac{d\Theta}{d\rho}-\frac{\sin\Theta\cos\Theta}{\rho^{2}}+\frac{2}{\rho}\sin^{2}\Theta$
$\displaystyle=\frac{H}{H_{0}}\left(\sin\Theta-\gamma\sin(\Theta)\cos(\Theta)\right)$
(17) where $H_{0}=A_{0}=D^{2}/J$. Numerical solution of Eq. (17) shown in Fig.
5 for different values of the parameters $H/H_{0}$ and $\gamma=2A/H$. We note
that while the long-distance solution is independent on these parameters, the
short and intermediate distance behavior is strongly parameter dependent,
suggesting $r_{d}$ is not the natural length scale of the system.

Figure 5: Skyrmion solution: Numerical solution of the Euler equation (17),
i.e., $\Theta(r)$ vs. $r/r_{d}$ for different values of $\gamma$ in the panels
(a)–(e) for same set of $H_{0}/H$, viz, 0.36, 0.49, 0.64 and 0.81.

We next introduce a length scale $r_{s}=D/H$ and rescaling $r\to r_{s}\rho$,
we obtain

$\displaystyle\frac{d^{2}\Theta}{d\rho^{2}}+\frac{1}{\rho}\frac{d\Theta}{d\rho}-\frac{\sin\Theta\cos\Theta}{\rho^{2}}$
$\displaystyle=\frac{H_{0}}{H}\left(-\frac{2}{\rho}\sin^{2}\Theta+\sin\Theta-\gamma\sin\Theta\cos\Theta\right)$
(18)

whose numerical solution (Fig. 1 of the article and Fig. 6 below) is almost
$H_{0}/H$ independent for a reasonable range. We thus find natural length
scale of the system as $r_{s}$. Together with the solution of $\Theta(r)$, (i)
$\Phi=\phi$ for $\alpha=1$, (ii) $\Phi=\phi+\pi/2$ for $\alpha=2$ and (iii)
$\Phi=-\phi$ for $\alpha=-3$ respectively construct magnetic structures of
Neél type skyrmion, Bloch type skyrmion and Neél type anti-skyrmion.

Figure 6: Skyrmion solution: Numerical solution of the Euler equation (18),
i.e., $\Theta(r)$ vs. $r/r_{s}$ for different values of $\gamma$ and $H_{0}/H$
in the panels (a)–(e), covering most of the skyrmion phase shown in Fig. 3(a).

However, the nature of the solution of Eq.(18) changes with $A/A_{0}$ for low
field regime. For a moderate to large $H/H_{0}$, any amount of negative (easy-
axis) anisotropy provide normal skyrmion solution (Fig. 7(a)) as shown in Fig.
1 and Fig. 6. When magnetic filed is low and $|A|/A_{0}\lesssim\pi/2\sqrt{2}$
(Dzyaloshinskii criterion S (1) for non-collinear state at zero $H$), the
nature of the solution is ‘chiral-bubble’ like S (2)–where very slow change of
$\Theta$ occurs near $r=0$, as shown in Fig. 7(b). The nature of the solution
changes for $|A|/A_{0}\gg\pi/2\sqrt{2}$ at low $H/H_{0}$ from chiral bubble to
‘meta-stable’ skyrmion, as shown in Fig. 7(c), for which $\Theta$ sharply
falls near $r=0$. The behavior of metastable skyrmions are, however,
fundamentally different from the normal skyrmions presented in the main text.
While the normal skyrmions are appropriately scaled with $r_{s}$, the
metastable skyrmions are better suited with the length scale $r_{d}$ (see Fig.
8). Further in contrary to the metastable skyrmions, the radius of a normal
skyrmion increases rapidly with the decrease of $H$. The metastable skyrmions
are energetically unfavorable to polarized ferromagnet.

Figure 7: Numerical solution of the Euler equation (18), i.e., $\Theta(r)$ vs.
$r/r_{s}$ for different values of $H/H_{0}$ and $A/A_{0}$. The nature of the
solutions changes with the parameters: These are (a) skyrmions (b) chiral
bubbles (c) meta-stable skyrmions (also chiral bubbles).  Figure 8: Numerical
solution of the Euler equations (18) and (17), i.e., $\Theta(r)$ vs. $r/r_{s}$
and $r/r_{d}$ in 9(a) and (b) respectively for different values of $H/H_{0}$
and $A/A_{0}$. Clearly, these metastable skyrmion solutions are better suited
with the length scale $r_{d}$.

Euler equation and its meron solution:

The Euler equation (12) will reduce to

$J\left(\Theta_{rr}+\frac{\Theta_{r}}{r}-\frac{\sin\Theta\cos\Theta}{r^{2}}\right)+\frac{2D}{r}\sin^{2}\Theta+A\sin(2\Theta)=0$
(19)

in the absence of magnetic field. For a sufficiently high easy-plane
anisotropy $(A>0)$, all the spins will align in the plane (planar
ferromagnet). This indicates a boundary condition $\Theta(\infty)=\pi/2$ which
together with another boundary condition $\Theta(0)=0$ or $\pi$ will provide a
solution for meron. By introducing a length scale $r_{a}=D/A$ and rescaling
$r\to r_{a}\rho$, Eq. (19) will reduce to

$\frac{d^{2}\Theta}{d\rho^{2}}+\frac{1}{\rho}\frac{d\Theta}{d\rho}-\frac{\sin\Theta\cos\Theta}{\rho^{2}}=-\frac{A_{0}}{A}\left(\frac{2}{\rho}\sin^{2}\Theta+\sin(2\Theta)\right)$
(20)

Taking cue of the skyrmion solution, we assume the solutions of meron is in
the form

$\Theta(r)=\pm\frac{\pi}{2}+2\arctan(\exp(-\zeta r/r_{a}))$ (21)

where positive (negative) sign corresponds to spin down (up) at the center of
the meron, and the parameter $\zeta$ to be determined by minimizing the
corresponding energy. The energy of a meron is given by

$\displaystyle E_{{\rm meron}}$ $\displaystyle=$ $\displaystyle
2\pi\int_{0}^{\infty}r\,dr\left[\frac{J}{2}\left(\Theta_{r}^{2}+\frac{\sin^{2}\Theta}{r^{2}}\right)+A\cos^{2}\Theta\right.$
(22)
$\displaystyle+\left.D\left(\Theta_{r}+\frac{\sin(2\Theta)}{2r}\right)\right]$

which may be simplified to

$E_{{\rm
meron}}=\frac{J}{2}\left[\ln(2)+I_{1}\right]+\frac{D^{2}}{A}\left[\frac{\ln(2)}{\zeta^{2}}-\frac{I_{2}}{\zeta}\right]$
(23)

where $I_{1}=\int_{0}^{\infty}(1/r)\tanh^{2}(r)dr$,
$I_{2}=\int_{0}^{\infty}\left(\frac{\tanh r}{\cosh r}+\frac{r}{\cosh
r}\right)dr=1+2G$, and Catalan’s constant
$G=\sum_{k=0}^{\infty}\frac{(-1)^{k}}{(2k+1)^{2}}\approx 0.91$.

## References

  * S (1) I. E. Dzyaloshinskii, Sov. Phys. JETP 20, 665 (1965).
  * S (2) See, for example, A. Fert, N. Reyren, and V. Cros, Nature Review Materials 2, 17031 (2017).


# Numerical MHD simulations of solar flares and their associated small-scale
structures

Mauricio González-Servín 1 and J. J. González-Avilés 2
1 Universidad Nacional Autónoma de México. Instituto de Radioastronomía y
Astrofísica. 58089. Morelia, Michoacán, México.
2 Universidad Nacional Autónoma de México. Escuela Nacional de Estudios
Superiores, Unidad Morelia. 58190. Morelia, Michoacán, México E-mail:
<EMAIL_ADDRESS>
###### Abstract

Using numerical simulations, we study the formation and dynamics of solar
flares in a local region of the solar atmosphere. The magnetohydrodynamics
(MHD) equations describe the dynamic evolution of flares, including space-
dependent and anomalous magnetic resistivity and highly anisotropic thermal
conduction on a 2.5 D slice. We adopt an initial solar atmospheric model in
magnetohydrostatic equilibrium, with a magnetic configuration consisting of a
vertical current sheet, which helps trigger the magnetic reconnection process.
Specifically, we study three scenarios, two with only resistivity and the
third with resistivity plus thermal conduction. The main results of the
numerical simulations show differences in the global morphology of the flares,
including the post-flare loops and the current sheet in three cases. In
particular, localized resistivity produces more substructure around the post-
flare loops that could be related to the Ritchmyer-Meshkov Instability (RMI).
Furthermore, in the scenario of anomalous resistivity, we identify the
formation of a plasmoid and a jet at coronal heights. On the other hand, in
the scenario with resistivity plus thermal conduction, the post-flare loops
are smooth, and no apparent substructures develop. Besides, in the $z-$
component of the current density for the Res+TC case, we observe the
development of multiple magnetic islands generated due to the Tearing
instability in the non-linear regime.

###### keywords:

(magnetohydrodynamics) MHD – Sun: flares – Sun: corona – methods: numerical

††pubyear: 2023††pagerange: Numerical MHD simulations of solar flares and
their associated small-scale structures–Numerical MHD simulations of solar
flares and their associated small-scale structures

##  1 Introduction

Solar flares are the most energetic explosions in the heliosphere, and they
power space weather events that can affect the magnetic environment of our
Earth. Solar flares efficiently convert magnetic energy into internal and
kinetic energy, producing emissions in different wavelengths. One of the
standard models of the solar flares has been discussed in Tajima & Shibata,
2002; Forbes, 2010; Priest, 2014a, where the authors extended the cartoon view
of magnetic reconnection using observational facts and, therefore, identify
processes such as fast evaporation upflows, X-ray loops, and hard X-ray (HRX)
source regions (Masuda et al., 1994; Su et al., 2013). One of the most
relevant 2D MHD simulations of a solar flare based on a magnetic reconnection
model with anisotropic thermal conduction was performed by Yokoyama & Shibata
(2001), who studied the thermal evolution in the post-flare loops to determine
the flare temperature. Similar MHD simulations of magnetic reconnection have
been performed to study coronal mass ejections and associated giant arcades
(Shiota et al., 2005), and chromospheric evaporation jets (Miyagoshi &
Yokoyama, 2004). More recently, Ruan et al. (2020) presented a self-consistent
solar MHD flare model that demonstrated the observationally suggested
relationship between flux swept out by the hard X-ray footpoint regions and
the actual reconnection rate at the X-point, which is a significant unknown in
flaring scenarios. In solar flares, various plasma flows above the post-
reconnection flare arcade have been reported by Liu et al. (2013); Lin et al.
(2005); Warren et al. (2018). In particular, flare current sheets have rapidly
descended, dark finger-like features commonly referred to as supra-arcade
downflows (SADs), often observed in extreme ultraviolet (EUV) soft X-ray
images. The slow speeds of SADs have been discussed in numerical models
containing the bursty jet (Cécere et al., 2015) or Rayleigh-Taylor-type (RT)
instabilities in reconnection downflow regions (Guo et al., 2014). Other
features observed in post-flare loops are the formation of magnetic islands at
the top of the loops and along the current sheet (CS). These magnetic islands
are believed to be closely related to fine structures CS in the wake of the
eruption (Shen et al., 2011; Lin et al., 2005). Moreover, the magnetic
islands’ coalescence can further produce secondary CSs and islands (Bárta et
al., 2011). More recently, Wang et al. (2021) performed 2.5D MHD simulations
to study the dynamics of magnetic reconnection in the solar CS, and they
suggest that annihilation of magnetic islands at the flare loop top, which is
not included in the standard flare model, plays a non-negligible role in
releasing magnetic energy to heat flare plasma and accelerate particles.

In this paper, we study the formation and evolution of a solar flare,
including the post-flare loops, the CS, and the development of associated
substructures such as shocks, MHD instabilities, and small-scale magnetic
islands, including the effects of localized and anomalous resistivity, and
highly anisotropic thermal conduction on the formation and dynamics of a solar
flare following the main ideas of the classical model developed by Yokoyama &
Shibata (2001). In particular, we base our simulations on three cases: i) Res
case, where we only consider the effects of a localized resistivity; ii)
anomalous resistivity case; and iii) Res+TC, where we include the effect of
the highly anisotropic thermal conduction additional to the localized
resistivity. We briefly describe the results of case ii) corresponding to the
anomalous resistivity profile.

We organize the paper as follows. First, in section 2, we describe the system
of MHD equations, the solar atmospheric model, the magnetic field
configuration, the resistivity and thermal conduction models, and the
numerical methods used to solve the MHD equations. Next, section 3 contains
the results of the two study cases, including the analysis of MHD
instabilities, shock development, and magnetic island formation. Finally, we
draw our concluding remarks in section 4.

We organize the paper as follows. First, in section 2, we describe the system
of MHD equations, the solar atmospheric model, the magnetic field
configuration, the resistivity and thermal conduction models, and the
numerical methods used to solve the MHD equations. Next, section 3 contains
the results of the two study cases, including the analysis of MHD
instabilities, shock development, and magnetic island formation. Finally, we
draw our concluding remarks in section 4.

##  2 Model and numerical methods

###  2.1 The system of MHD equations

We consider a gravitationally stratified solar atmosphere described by a
plasma obeying the resistive MHD equations with anisotropic thermal conduction
to model the plasma that constitutes the solar flare. In particular, we write
the equations in a conservative dimensionless form that helps apply the
numerical methods easily:

$\displaystyle\frac{\partial\varrho}{\partial t}+\nabla\cdot(\varrho{\bf v})$
$\displaystyle=0,$ (1) $\displaystyle\frac{\partial(\varrho{\bf v})}{\partial
t}+\nabla\cdot(\varrho{\bf v}{\bf v}-{\bf B}{\bf B}+p_{t}{\bf I})$
$\displaystyle=\varrho{\bf g},$ (2) $\displaystyle\frac{\partial E}{\partial
t}+\nabla\cdot((E+p_{t}){\bf v}-{\bf B}({\bf v}\cdot{\bf B}))$
$\displaystyle={\varrho}{\bf v}\cdot{\bf g}$
$\displaystyle-\nabla\cdot(\eta{\bf J}\times{\bf B})+\nabla\cdot{\bf F_{c}},$
(3) $\displaystyle\frac{\partial{\bf B}}{\partial t}+\nabla\cdot({\bf v}{\bf
B}-{\bf B}{\bf v})$ $\displaystyle=-\nabla\times({\bf\eta J}),$ (4)
$\displaystyle\nabla\cdot{\bf B}$ $\displaystyle=0,$ (5)

where $\varrho$ is the mass density, ${\bf v}$ represents the velocity field,
${\bf B}$ is the magnetic field, $T$ is the temperature, $\eta$ is the
magnetic resistivity, and ${\bf F_{c}}$ is the thermal conduction flux.

Besides, $p_{t}=p+{\bf B}^{2}/2$ is the total (thermal + magnetic) pressure,
${\bf I}$ is the unit matrix, and $E$ is the total energy density, that is,
the sum of the internal, kinetic, and magnetic energy densities,

$E=\frac{p}{\gamma-1}+\frac{\varrho{\bf v}^{2}}{2}+\frac{{\bf B}^{2}}{2}.$ (6)

For the fluid, we consider the adiabatic index $\gamma=5/3$. The system of
equations (1)-(5) is closed with the ideal gas law

$p=\frac{k_{B}}{\bar{m}}\varrho T,$ (7)

where $T$ is the temperature of the plasma, $\bar{m}=\mu m_{H}$ is the
particle mass specified by a mean molecular weight value $\mu=$0.6, for a
fully ionized gas, and $m_{H}$ is hydrogen’s mass, and $k_{B}$ is Boltzmann’s
constant. The gravitational source term on the right-hand side of equations
(2) and (3) is given by ${\bf g}=[0,-g]$ with magnitude $g=2.74\times 10^{4}$
cm s-2, which represents an average over the solar surface. In this paper, we
use a 2.5D model, meaning all the state variables depend on $x$ and $y$, where
$x$ is a horizontal coordinate, and $y$ represents height. Besides, a
nontrivial magnetic field $z-$component $B_{z}$ varies with $x$ is allowed. We
solve the 2.5D version of the resistive MHD equations with thermal conduction
on the $xy$ plane with resolution $\Delta x$ and $\Delta y$.

####  2.1.1 Model of the solar atmosphere

At the initial time of the simulations, we assume the solar atmosphere is in
hydrostatic equilibrium. In particular, we choose the C7 semiempirical model
to describe the temperature field of the chromosphere-transition region
(Avrett & Loeser, 2008), smoothly extended to the solar corona as in González-
Avilés et al. (2017). In Fig. 1, we display the equilibrium temperature and
mass density profiles as functions of height $y$, where the steep gradient at
the transition region $y\sim 2.1$ Mm is discernible.

Figure 1: Equilibrium temperature in Kelvin (left) and mass density in gr cm-3
(right) as functions of height $y$.

####  2.1.2 The magnetic field

Initially, we consider a force-free magnetic field ($(\nabla\times{\bf
B})\times{\bf B}={\bf 0}$) similar to that in Yokoyama & Shibata (2001); Ruan
et al. (2020), which essentially represents a vertical CS defined as follows:

$\displaystyle B_{x}$ $\displaystyle=$ $\displaystyle 0,$ (8) $\displaystyle
B_{y}$ $\displaystyle=$ $\displaystyle B_{0}\tanh(x/w),$ (9) $\displaystyle
B_{z}$ $\displaystyle=$ $\displaystyle B_{0}/\cosh(x/w),$ (10)

here $B_{0}=22.42$ G is the magnetic field strength and $w=5\times 10^{7}$ cm
is the width of the CS. We show the 2D maps of the $B_{y}$ and $B_{z}$
components in Fig. 2, where it is notable that both components only vary in
$x$ and are constant along the vertical direction $y$. However, initially,
there is a non-negligible $z-$ component of the current density, which
produces a CS that helps trigger the magnetic reconnection at late times in
the simulation. Additionally, we plot the plasma $\beta$, which is a parameter
that estimates the ratio between gas pressure and magnetic pressure, i.e.,

$\beta(x,y)=p(y)/B^{2}/2.$ (11)

Here, the pressure $p(y)$ is given by the hydrostatic equilibrium, and
$B^{2}=(B_{y}^{2}+B_{z}^{2})$. We display the plasma $\beta$ at the initial
time in the right of Fig. 2, where it is evident that $\beta$ is about $10$ in
the bottom of the photosphere $y=0$ Mm, while in the solar corona ($y>2.1$ Mm)
is approximately $10^{-3}$. This parameter will help us explain the dominance
of pressures in regions where the magnetic islands form, as described in
Section 3. Moreover, the solar atmospheric model, which is in hydrostatic
equilibrium, and the magnetic field configuration, which is force-free,
satisfying the magnetohydrostatic equation, ${\bf J}\times{\bf B}-\nabla
p-\rho g={\bf 0}$. Therefore, the whole system achieves the magnetohydrostatic
equilibrium at the initial time of the simulations.

Figure 2: Magnetic field components $B_{y}$ (left) and $B_{z}$ (middle) in
Gauss, and plasma $\beta$ (right) at the initial time of the simulations.

####  2.1.3 Spatially dependent and anomalous resistivity profiles

To perturb the atmosphere in magnetohydrostatic equilibrium and trigger the
magnetic reconnection, we set a spatial dependent function for the
resistivity, $\eta(x,y)$. In particular, we use a similar localized
resistivity function as in Takasao et al. (2015), which is defined as follows:

$\eta(x,y)=\eta_{0}\exp\left[{-(r/w_{\eta})^{2}}\right],$ (12)

where $\eta_{0}=10^{16}$ cm2 s-1 is the strength of the resistivity,
$r=\sqrt{x^{2}-(y-h_{p})^{2}}$, with $h_{p}=6\times 10^{8}$ cm (6 Mm) is the
vertical location of the Gaussian profile, and $w_{\eta}=1\times 10^{8}$ cm (1
Mm) is its width. This profile is time-independent, but it can induce a fast
and quasi-steady magnetic reconnection with a single X-point (see, e.g., Ugai,
1992).

We also include an anomalous resistivity function. Microscopic instabilities
physically cause this resistivity. Specifically, we adopt the following
profile:

$\eta(x,y)=\begin{cases}\text{min}(\alpha_{\eta}(\nu_{d}/\nu_{c})-1,\eta_{min}),\
\text{if}\ t<t_{\eta}\\\ 0,\ \text{if}\ t>t_{\eta}\end{cases}$ (13)

where, $\alpha_{\eta}=1\times 10^{-4}$, $\eta_{min}=$0.1, and $t_{\eta}=10$ s.
In the expression above, $\nu_{d}(x,y)={\bf J}/(eN_{e})$ represents an
instantaneous distribution that quantifies the relative ion-electron drift
velocity, where ${\bf J}$ is the current density magnitude, $e$ is the
electron charge and $N_{e}$ is the electron number density. Profile of
equation 13 is similar to the used in Ruan et al. (2020); however, we do not
couple with the localized profile of equation 12 and consider a small
threshold drift velocity $\nu_{c}$.

####  2.1.4 Anisotropic thermal conduction

In the equation (3), the thermal conduction flux ${\bf F_{c}}$ defines a flux-
limited expression that smoothly varies between the classical and saturated
thermal conduction regimes $F_{class}$ and $F_{sat}$, respectively, that is:

${\bf F_{c}}=\frac{F_{sat}}{F_{sat}+|{\bf F}_{class}|}{\bf F}_{class}.$ (14)

The above equation represents highly anisotropic thermal conduction suppressed
mainly in the direction transverse to the magnetic field lines. For example,
in this paper, we denote $\hat{{\bf b}}={\bf B}/|{\bf B}|$ as the unit vector
in the direction of the magnetic field; therefore, the classical thermal
conduction flux is as follows:

$\displaystyle{\bf F}_{class}$ $\displaystyle=\kappa_{\parallel}\hat{{\bf
b}}(\hat{{\bf b}}\cdot\nabla T)+\kappa_{\perp}[\nabla T-\hat{{\bf
b}}(\hat{{\bf b}}\cdot\nabla T)],$ (15) $\displaystyle|{\bf F}_{class}|$
$\displaystyle=\sqrt{({\bf\hat{b}}\cdot\nabla
T)^{2}(\kappa_{\parallel}^{2}-\kappa_{\perp}^{2})+\kappa_{\perp}^{2}\nabla
T^{2}}$ (16)

where the subscripts $\parallel$ and $\perp$ denote, respectively, the
parallel and normal components to the magnetic field, $T$ is the temperature,
$\kappa_{\parallel}$ and $\kappa_{\perp}$ are the thermal conduction
coefficients along and across the field. In most astrophysical scenarios, such
as in the Sun, $\kappa_{\perp}/\kappa_{\parallel}=2\times
10^{-31}n_{H}^{2}/(T^{3}|{\bf B}|^{2})$ is negligible (Priest, 2014b), and the
thermal conduction is mainly along the field lines. In this paper, we chose
$\kappa_{\parallel}=8\times 10^{-7}$ erg s-1 cm-1 K-1, which is a typical
value in the solar corona (Ruan et al., 2020). Saturated effects are accounted
for, making the flux independent of $\nabla T$ for substantial temperature
gradients, such as the transition region to solar corona in the solar
atmosphere, which is a novel ingredient of this paper compared to the similar
works where include the thermal conduction. In this limit, the flux magnitude
approaches $F_{sat}=5\phi\rho c_{iso}^{3}$, where $c_{iso}$ is the isothermal
speed of sound and $\phi<1$ is a free parameter that helps to account for the
uncertain numerical coefficient in $F_{sat}$ (Giuliani, 1984), which we take
equal to 0.3 in this paper. The equations 1-5 were normalized by the scale
factors in CGS units as shown in table 1, which are typical scales of the
solar atmosphere. The CGS units are used throughout the manuscript to be
consistent with the normalization factors.

Variable | Quantity | Unit | Value
---|---|---|---
$x$,$y$ | Length | $l_{0}$ |  $10^{8}$ cm
$\varrho$ | Density | $\varrho_{0}$ |  $10^{-15}$ gr cm-3
${\bf v}$ | Velocity | $v_{0}$ |  $10^{8}$ cm s-1
$t$ | Time | $t_{0}=l_{0}/v_{0}$ | 1 s
$p$ | Gas pressure | $p_{0}=\varrho_{0}v_{0}^{2}$ | 10 gr cm-1 s-2 (dyn cm-2)
${\bf B}$ | Magnetic field | $B_{0}=\sqrt{4\pi\varrho_{0}v_{0}^{2}}$ | 11.21 G
${\bf J}$ | Current density | $J_{0}=cB_{0}/l_{0}$ | 3360.5 gr1/2 cm-1/2 s-2 (statA cm-2)
$\eta$ | Resistivity | $\eta_{0}=v_{0}l_{0}$ |  $10^{16}$ cm2 s-1
$\kappa$ | Thermal conductivity | $\kappa_{0}=\mu m_{u}/\varrho_{0}v_{0}l_{0}k_{B}$ |  $7.265\times 10^{-10}$ erg s-1 cm-1 K-1
Table 1: Normalization scale factor in CGS units

###  2.2 Numerical methods

We numerically solve the equations (1)-(5) using the PLUTO code (Mignone et
al., 2007). In all simulations, we set the Courant-Friedrichs-Levy (CFL)
number equal to 0.4 and choose the second-order total variation diminishing
(TVD) Runge-Kutta time integrator. Additionally, we use the Harten-Lax-van
Leer discontinuities (HLLD) approximate Riemann solver (Miyoshi & Kusano,
2005) in combination with a linear reconstructor and the minmod limiter.

The numerical evolution of the MHD equations can lead to the violation of the
divergence-free constraint equation given by (5), developing unphysical
results like the presence of a net magnetic charge. We use the extended
generalized Lagrange multiplier method (Dedner et al., 2002) to control the
constraint violation’s growth. This method is robust in highly magnetized
regions or low plasma beta scenarios, as implied in the solar corona on the
left panel of Fig. 2. PLUTO evolves the magnetic resistivity and thermal
conduction separately from advection through operator splitting, using the
super-time-stepping technique (Alexiades et al., 1996). This method is robust
enough to handle the highly parabolic (diffusion) behavior of space-dependent
resistivity and anisotropic thermal conduction in the MHD equations adopted in
this paper. Our simulations are carried out in a domain $x\in[-10,10]$,
$y\in[0,20]$, in units of Mm ($10^{8}$ cm), covered by 800$\times$1000 grid
cells. Here, $y=0$ Mm represents the bottom of the photosphere. We used the
following boundary conditions: outflow boundary conditions at $x_{min}=-10$
Mm, $x_{max}=10$ Mm, and $y_{max}=20$ Mm. In contrast, at $y_{min}=0$ Mm, we
employ asymmetric boundary conditions in the $y-$component of the magnetic
field $B_{y}$, outflow in the $z-$component of the magnetic field $B_{z}$,
mass density and pressure, and zero boundary conditions in the component of
the velocity field. These boundary conditions are similar to those used in
Yokoyama & Shibata (2001).

##  3 Results of numerical simulations

In this section, we show the most representative results of the numerical
simulations performed for the following three scenarios: (i) resistivity case
(Res), ii) anomalous resistivity case, and (iii) resistivity plus thermal
conduction case (Res+TC). We only include the most representative results for
the case of the anomalous resistivity profile.

###  3.1 Resistivity case

In Fig. 3, we show the temperature in Kelvin and mass density in gr cm-3 with
the magnetic field lines for the Res case. For example, in panel (a), we show
the temperature at $t=25$ s. There, we note the following plasma structures
that are consistent with the standard flare model (Yokoyama & Shibata, 2001):
i) formation of a hot flare loop ($\sim 10^{8}$ K) due to the magnetic
reconnection triggered by the localized resistivity, ii) the reconnection
outflow collides with the flare loop, and iii) production of evaporation flows
of about $10^{6}$ K in the flare loop. In panel (c), the temperature at $t=50$
s shows the formation of a small hot region of about $10^{6}$ K at the top of
the loops, and it is also evident that inside the loop region, the
chromospheric evaporation develops. In panel (e), we show the mass density at
$t=25$ s, consistently proving the behavior demonstrated on temperature maps.
In particular, chromospheric evaporation appears in dense regions of about
$10^{-14}$ gr cm-3. Finally, in panel (g), the mass density at $t=50$ s shows
how dense material ($\sim 10^{-12}$) moves upwards while it forms dense region
($\sim 10^{-17}$ gr cm-3) around it, and it is also visible the small region
(but hot) at the top of it, as shown by the temperature maps.

blackRes blackRes+TC

black(a) black(b)

black(c) black(d)

black(e) black(f)

black(g) black(h)

Figure 3: Temperature maps in kelvin overlap with the magnetic field lines at
$t=25$ s (panels (a) and (b)) and $t=50$ s (panels (c) and (d)). In panels
(e), (f), we show the mass density in gr cm-3 for the time $t=25$ s, while in
panels (g) and (h), we display mass density at $t=50$ s.

In Fig. 4, we show snapshots of the temperature in kelvin and mass density in
gr cm-3 at times $t=25,50$ s for the anomalous resistivity case. In this
scenario, we only turn on the resistivity given by the equation (13) and let
the system evolve. In panels (a) and (b), we display the temperature, where it
is discernible that any post-flare structure is formed, especially since we do
not identify the post-flare loop structures and the development of a thin CS
along the vertical direction. Instead, we can see the formation of a plasma
structure that becomes wider and moves upwards while the magnetic file lines
reconnect. This kind of plasma feature could be related to a plasmoid. In
panels (c) and (d), we note that in the mass density, there are no clues of
the formation of post-flare loops, but interestingly, we identify the
development of a small jet (3 Mm in length) moving upwards from the transition
region. This feature is interesting since it could be related to the magnetic
reconnection process at microscopic scales, which is triggered by anomalous
resistivity at coronal heights, which indicates a relevant result despite the
limitations of the MHD to capture the microscopic mechanism associated with
the magnetic reconnection process. Finally, it is discernible that these plots
are different from the plots of panels (a), (c), (e), and (g) shown in Fig. 3
for the Res case.

blackAnomalous resistivity case

black(a) black(b)

black(c) black(d)

Figure 4: Temperature maps in kelvin and mass density in gr cm-3 overlay with
the magnetic field lines at $t=25$ s (panels (a) and (c)) and at $t=50$ s
(panels (b) and (d)) for the anomalous resistivity case.

In Fig. 5, we show the $z-$component of current density $J_{z}$ in statA cm-2
and plasma pressure in dyn cm-2 with the magnetic field lines overlap for the
two simulation cases. Specifically, in panels (a), (c), (e), and (g), we
display the results for the Res case at times $t=25$ s and $t=50$ s. For
instance, in panel (a) at $t=25$ s, we identify substructures in $J_{z}$. In
particular, there are regions with high current density inside the post-flare
loops, and a thin region of high $J_{z}$ is visible that is part of the CS.
While, in panel (c), at $t=50$ s, we observe substructures at about $y\approx
12$ Mm and $y\approx 17$ Mm. In panel (e), at $t=25$ s, we identify regions of
high plasma pressure ($\sim 18$ dyn cm-2) inside the post-flare loops, and
$t=50$ s (g), in the plasma pressure is also evident the structure shown in
$J_{z}$ of panel (c), in this particular region, the plasma pressure is around
1.5 dyn cm-2.

blackRes blackRes+TC

black(a) black(b)

black(c) black(d)

black(e) black(f)

black(g) black(h)

Figure 5: The $z$ component of current density in statA cm-2 overlap with the
magnetic field lines for the time $t=25$ s (panels (a) and (b)) and $t=50$ s
(panels (c) and (d)). In panels (e) (f), we show plasma pressure in dyn cm-2
for the time $t=25$ s, while in panels (g) and (h), we display mass density in
gr cm-3 for $t=50$ s, for the Res case (right panels) and the Res+TC (left
panels).

###  3.2 Resistivity with thermal conduction case

In the right panels of Fig. 3, we display the results for the Res+TC case. For
example, in panel (b), we identify the formation of a typical flare structure
and its associated post-flare loops; however, in this case, the thermal
conduction redistributes the plasma inside the flare loops, which makes the
evaporation flow less visible than in the Res case. Furthermore, the maximum
temperature reaches $\sim 10^{7}$ K, smaller than the maximum temperature in
the Res case, as shown in panel (a). In panel (d), we do not observe any
structured region on the top of loops, nor the increase of evaporation of the
plasma inside it, as shown in panel (c). This feature could be related to the
effect of thermal conduction along the magnetic field lines, which in this
case would be high along the loops, redistributing the heat inside the loops
and dominating the resistivity. Panel (f) shows the mass density at $t=25$ s,
which is consistent with the behavior shown in the temperature maps of panels
(b) and (d). In particular, inside the loops, the density is of the order of
$10^{-15}$ gr cm-3. Finally, in panel (h), the mass density at $t=50$ s shows
a smooth distribution predominantly inside the loops.

In panels (b), (d), (f), and (h) of Fig. 5, we show the z-component of current
density $J_{z}$ in statA cm-2 and plasma pressure in dyn cm-2 for the Res+TC
case. For instance, in panel (b), $J_{z}$ at $t=25$ s shows similarities to
the obtained in the Res case. However, at $t=50$ s, as displayed in panel (d),
we do not identify the substructure at the top of the post-flares as in the
Res case. Instead, the CS becomes thin. In panels (f) and (h), the plasma
pressure shows the same smooth behavior as $J_{z}$ at both times, $t=25$ and
$t=50$ s. Notably, the regions with high plasma pressure ($\sim 18$ dyn cm-2)
are mainly around the post-flare loops.

To show how the thermal conduction behaves on the flare and post-flare loop
structures, we estimate the magnitude of the heat flux q. In Fig. 6, we
display heat flux maps for two times $t=25$ and $t=50$ s. We see that the heat
flux is high within the post-flare loops and the CS, which indicates that heat
conduction could determine the thermal behavior of the post-flare loops.
Therefore, this makes the structure smooth, as shown in panels (b), (d), (f),
and (g) in the temperature and density maps of Fig. 3.

blackRes+TC

black(a) black(b)

Figure 6: Maps of the local heat flux $|{\bf q}|$ in erg s-1 cm-2.

In Fig. 7, we show 1D cuts along the vertical distance $y$ and at $x=0$ Mm of
$J_{z}$ in statA cm-2 (top panels) and the temperature in Kelvin (bottom
panels), at times $t=25,50$ s. These plots help identify the current density
and temperature variations with the vertical distance and highlight the
behavior between these variables in the Res and Res+TC cases. For example, in
the top-left panel, we show that $J_{z}$ at $t=25$ s behaves similarly in both
scenarios; however, in the Res case, there is an abrupt growth at around the
distance $y=3$ Mm, which is related to a high current density region, as shown
in panel (a) of Fig. 5. On the other hand, at $t=50$ s, the current density
variations are remarkably different for the Res+TC case in the interval of
10-20 Mm than the behavior in the Res case. These variations could be related
to the formation of multiple small-scale magnetic islands along the CS, as we
will analyze in detail in the following subsections. Finally, in the bottom
panels, we show the variations of temperature as a function of the distance,
which behave consistently with the results shown in Fig. 3. In particular, we
observe a temperature growth at about $y\approx 3$ Mm, in both cases at the
two times. Nevertheless, in the Res case, the plasma temperature reaches
higher values ($\sim 10^{9}$ K). In contrast, in the Res+TC case, the
temperature of the plasma maintains an order of magnitude of about $10^{7}$ K.
We can interpret these results in terms of the effects that thermal conduction
has over the plasma temperature, i.e., in the Res case, the heat continues to
affect the fluid due to high dynamics. In contrast, in the Res+TC case, the
heat is dissipated and redistributed, so the temperature does not rise
considerably.

Figure 7: 1D cuts along the vertical direction $y$ of the $z-$component of the
current density $J_{z}$ in statA cm-2 (top) and temperature in Kelvin (bottom)
at times $t=25$ s (left panels) and $t=50$ s (right panels) the Res case (blue
curves) and the Res+TC case (orange curves).

####  3.2.1 Instabilities

In general, in a fluid the hydrodynamic instabilities such as Rayleigh-Taylor
(RT) and Ritchmyer-Meshkov (RM) frequently appear when turbulence takes place.
In particular, the RT instability occurs at an interface between different
fluids when the lighter fluid is accelerated into the heavy. While the RM
instability may be considered a particular case of RT instability, it develops
when the acceleration provided is impulsive, resulting from a shock wave (Zhou
et al., 2021).

Furthermore, in Fig. 8, we show a zoom region of the temperature and vertical
velocity overlap with the velocity vector field at time $t=50$ s, for the Res
case (panels (a) and (c)) and for the Res+TC case (panels (b) and (d). In
panels (a) and (c) of the figure, we identify plasma flowing downward and
upward, shown in the velocity vector field; therefore, these flows could
generate the ’spike’ structures inside the post-flare loop regions. In
contrast, in panels (b) and (d), we see that plasma moves predominantly
upwards, and the post-flare loop looks homogeneous; this is consistent with
the behavior of the heat flow that is high inside the flare loops, as
schematized by the maps of $|{\bf q}|$ in Fig. 6. In the four figures is an
orange vertical dashed line, which represents the distance where the pressure,
density, and velocity gradients are estimated for the instabilities analysis
shown below.

blackRes blackRes+TC

black(a) black(b)

black(c) black(d)

Figure 8: Zoom of the temperature in kelvin (top panels), and vertical
velocity $v_{y}$ in cm s-1 (bottom panels) overlap with vector velocity field
at $t=50$ s, corresponding to the Res case (panels (a) and (c)), and the
Res+TC case (panels (b) and (d)).

To investigate the possible development of RTI/RMI, we estimate the variations
of the pressure gradient, density gradient, vertical velocity gradient, and
plasma $\beta$ along the orange vertical dashed line drawn in the panels of
Fig. 8, that is similar to the analysis performed by Shen et al. (2022).
Notably, we estimate how the above-mentioned variable changes when crossing
the interface around the transition region ($y=\sim$2.1 Mm). In Fig. 9, we
show 1D cuts of pressure gradient in dyn cm-3, density gradient in gr cm-4,
vertical velocity ($v_{y}$) gradient in s-1 and plasma $\beta$, corresponding
to the Res case in panels (a), (b), (c), and (d), and to the Res+TC case in
panels (e), (f), (g) and (h) for times $t=25$ s and $t=50$ s.

The RMI generally occurs when the sign of the density and pressure gradients
are opposite, i.e., $\nabla\rho\cdot\nabla p<0$, if a shock propagates from
the light to the dense gas. Nevertheless, if a shock propagates from the heavy
to the light gas, which is the scenario of our simulations, the RMI should
occur when $\nabla\rho\cdot\nabla p>0$. However, in the Res case, we see that
pressure and density gradients have the same signs just below the interface
($y=\sim$2.1 Mm) at $t=25$ s; therefore, the condition $\nabla\rho\cdot\nabla
p>0$ it seems to be satisfied globally. If we analyze locally, let us say at
$y=2$ Mm, there $\nabla p\approx-1.335\times 10^{-6}$, and
$\nabla\rho\approx-4.736\times 10^{-20}$, so the condition
$\nabla\rho\cdot\nabla p>0$ it is satisfied. If we estimate the condition over
the interface, at about $y=2.5$ Mm, we get $\nabla p\approx-2.993\times
10^{-7}$, and $\nabla\rho\approx-4.326\times 10^{-22}$, and the condition is
also satisfied. If we go higher, i.e., at about $y=4.5$ Mm, $\nabla
p\approx-8.741\times 10^{-8}$, and $\nabla\rho\approx-3.753\times 10^{-24}$,
the condition remains. Now, if we estimate the condition for the same point,
but at $t=50$ s, we obtain that at $y=2$ Mm, $\nabla p\approx-2.014\times
10^{-8}$, and $\nabla\rho\approx-3.121\times 10^{-20}$, so the condition
$\nabla\rho\cdot\nabla p>0$ is present there. If we go over the interface, at
about $y=2.5$ Mm, we get, $\nabla p\approx 6.095\times 10^{-8}$, and
$\nabla\rho\approx-4.730\times 10^{-20}$, which does not satisfy the
condition. At $y=4.5$ Mm, we get, $\nabla p\approx-1.659\times 10^{-8}$, and
$\nabla\rho\approx-1.716\times 10^{-22}$, and again, $\nabla\rho\cdot\nabla
p>0$. The same happen at about $y=5$ Mm, where $\nabla p\approx-7.299\times
10^{-8}$, and $\nabla\rho\approx-7.308\times 10^{-23}$. According to the above
estimations, we identify that the condition $\nabla\rho\cdot\nabla p>0$ takes
place for most of the points, so it is likely the development of an RMI for
the Res case.

The vertical velocity gradient (panel (c)) shows a positive value near the
interface but higher up, and this gradient varies from negative to positive
values, which means the presence of upward and downward flows. In panel (c),
plasma $\beta$ at $t=50$ s shows an evident variation from a magnetically
dominant regime ($\beta<1$) to a fluid-dominated regime ($\beta>1$). In
particular, the RMI features appear in regions where $\beta\approx 1$, i.e.,
where the fluid and magnetic pressure are approximately equal. The latter
results are consistent with the found in Shen et al. (2022). On the other
hand, in the Res+TC case, we obtain that for $t=25$ s, $\nabla
p\approx-8.534\times 10^{-8}$, and $\nabla\rho\approx-8.263\times 10^{-22}$ at
$y=2$ Mm, while at $y=2.5$ Mm, $t=25$ s, $\nabla p\approx-3.094\times
10^{-8}$, and $\nabla\rho\approx-8.179\times 10^{-24}$. At approximately
$y=4.5$ Mm, $t=25$ s, $\nabla p\approx-8.978\times 10^{-8}$, and
$\nabla\rho\approx-1.189\times 10^{-22}$. Finally, at $y=5$ Mm, $\nabla
p\approx-5.298\times 10^{-8}$, and $\nabla\rho\approx-1.189\times 10^{-22}$,
and the condition $\nabla\rho\cdot\nabla p>0$ satisfies for all the points.
But, if we estimate for $t=50$ s, we obtain $\nabla p\approx-6.398\times
10^{-7}$, and $\nabla\rho\approx 3.380\times 10^{-20}$ at $y=2$ Mm, while at
$y=2.5$ Mm, $\nabla p\approx 4.872\times 10^{-7}$, and
$\nabla\rho\approx-1.114\times 10^{-20}$. At approximately $y=4.5$ Mm, $\nabla
p\approx-7.993\times 10^{-8}$, and $\nabla\rho\approx 2.145\times 10^{-22}$.
Finally, at $y=5$ Mm, $\nabla p\approx-9.791\times 10^{-8}$, and
$\nabla\rho\approx-8.918\times 10^{-23}$. In most of the points, we get
$\nabla\rho\cdot\nabla p>0$. However, we do not see instabilities in all the
maps for the Res+TC case; consequently, this analysis does not apply to this
scenario. In panel (g), the vertical velocity gradient oscillates between
positive and negative values below and over the interface up to about $y=2.5$
Mm. However, at a higher distance, the gradient is positive, which indicates a
predominantly upward flow, consistent with the behavior shown in the 2D map of
the vertical velocity in panel (f) of Fig. 8. Finally, in panel (h), we note
that plasma $\beta$ is always smaller than 1, which means that magnetic
pressure dominates over fluid pressure, and therefore, the plasma is flowing
in this for the two times $t=25,50$ s. This behavior is consistent with any
observed RMI since there is no transition between magnetically-dominated and
fluid-dominated regimes.

According to the previous analysis, we can conclude that the internal
structure of the loop region is related to an RMI in the Res case. In the
Res+TC case, the thermal conductivity determines the plasma dynamics in the
loop region, so the plasma moves through the loops and dissipates the heat,
and therefore, any RMI develops.

blackRes case

black(a) black(b)

black(c) black(d)

blackRes+TC case

black(e) black(f)

black(g) black(h)

Figure 9: 1D cuts of pressure gradient in dyn cm-3, density gradient in gr
cm-4, vertical velocity ($v_{y}$) gradient in s-1 and plasma $\beta$ along the
vertical distance $y$ at $x=0$ Mm, corresponding to the Res case in panels
(a), (b), (c), and (d), and to the Res+TC case in panels (e), (f), (g) and (h)
for times $t=25$ s and $t=50$ s. The black vertical dashed line represents the
interface between the dense and tenuous plasma around the transition region at
$y\sim$2.1 Mm.

####  3.2.2 Shock development

To investigate if there are the physical conditions so that a shock could
develop in the post-flare loops structures of any of the two simulation cases,
we estimate the sonic Mach number, $M_{s}=V_{y}(C_{s}$ number, and the
Alfvénic Mach number, $M_{A}=V_{y}/C_{A}$, where, $V_{y}$ is the local
vertical velocity of the plasma, $C_{s}=\sqrt{\gamma p/\rho}$, is the local
speed of sound, and $C_{A}=|{\bf B}|/\sqrt{4\pi\rho}$, is the local Alfvén
speed. For example, in panel (a) of Fig. 10, we display the sonic Mach number
and the Alfvénic Mach number for the Res case at $t=50$ s. There, we see that
the sonic Mach number is greater than one, i.e., the plasma is supersonic in
the bottom regions of the loops, where the plasma propagates upwards. We also
note in panel (b) that the Alfvénic Mach number is greater than one in regions
where the plasma propagates upwards, which means that plasma flow is
superAlfvénic, and as a consequence of the two later conditions, a shock
develops in the Res case. On the other hand, in panels (c) and (d), we show
the sonic and Alfvénic Mach numbers for the Res+TC case. In this scenario, it
is discernible that plasma is neither supersonic nor superAlfvénic; therefore,
the environment is unsuitable for the shock to develop. The results shown for
the Res case in Fig. 10 are consistent with those obtained in Chen et al.
(2015) in the scenario of particle acceleration by a solar flare termination
shock.

blackRes case

black(a) black(b)

blackRes+TC case

black(c) black(d)

Figure 10: Sonic and Alfvénic Mach numbers, $M_{s}$ (panels (a) and (c)) and
$M_{A}$ (panels (b) and (d)) correspond to the Res case (top panels) and to
the Res+TC case (bottom panels) at $t=50$ s.

####  3.2.3 Magnetic islands formation

In the Res+TC case, we can zoom in, for example, the $J_{z}$ and observe the
CS closely to identify multiple magnetic islands, as shown in Fig. 11. In
panel (a), we display a zoom region that covers $x\in[-2,2]$ Mm and
$y\in[9,13]$ Mm, while in panel (b), we show a zoom region that covers
$x\in[-1,1]$ Mm and $y\in[15,20]$ Mm. For instance, in panel (a), we note the
development of three magnetic islands along the CS, i.e., in regions where
$J_{z}$ is high. These structures are also observed along the CS at higher
altitudes as schematized panels (b). In both panels, we can identify
separatrices regions that indicate the presence of the magnetic islands. These
regions represent the most recent reconnected field line that separates
regions with magnetic fields of different topology; these magnetic fields may
be anti-parallel, while the reconnected fields are found where the oppositely
directed fields have been interconnected. Generally, the separatrices meet at
the X-line, surrounded by a diffusion region, i.e., a high current density
region, as shown in the small plots indicated by the blue arrows in panel (a).
These magnetic islands can be formed due to a Tearing-instability in the non-
linear regime. The tearing-instability is a linear stability of a CS of
infinite (Furth et al., 1963), which is characterized by making a system in
equilibrium unstable due to the tearing mode, leading to the formation of
X-points and plasmoids during the reconnection process. This instability is
part of a basic resistive instability that is linked with a long-wave tearing
mode, corresponding to the breakup of the layer along current-flow lines (see,
e.g., Furth et al., 1963). It can occur on timescales $\tau$, such that
$\tau_{A}<\tau<\tau_{d}$. The latter times are the time last the instability
to traverse the CS at the diffusion speed $v_{d}$ and the Alfvén speed
$v_{A}$, respectively, and it has the effect of creating many small-scale
magnetic islands in the CS (Shen et al., 2011). Besides, these magnetic
islands are likely to develop in solar flares, (see, e.g., Kliem et al., 2000;
Wang et al., 2021).

blackRes+TC case

black(a) black(b)

Figure 11: Zoom of the $z$ component of the current density in statA cm-2 at
$t=50$ s that schematize the magnetic islands corresponding to the Res+TC
case.

To complement Fig. 11 and go into more detail in the formation of the magnetic
islands due to the Tearing-instability, in Fig. 12, we show a zoom region that
covers $x\in[-2,2]$ Mm and $y\in[9,20]$ Mm of the magnetic field lines (a) and
the plasma $\beta$ (b) at time $t=51$ s. For example, in panel (a), we note
the formation of various small-scale ($\sim 0.5$ Mm of width) magnetic
islands, which can be related to regions of high ($\sim 10^{4}$) plasma
$\beta$ values. We can also see that magnetic islands are souring by low
($\sim 10^{-1}$) plasma $\beta$ values, which means that overall magnetic
islands are localized in a mixed region. The latter means that magnetic
islands are more likely to form in regions where gas pressure dominates over
magnetic pressure, i.e., in regions of high plasma $\beta$. This result makes
sense since the magnetic islands form in X-points or O-points, where, by
definition, the magnetic field is zero. Besides, it is interesting that
magnetic islands are surrounded by regions where the magnetic pressure
dominates over gas pressure ($\beta\sim 10^{-1}$), which implies that the
magnetic field confines the material that forms the islands. Specifically, in
the Res+TC case, the thermal conduction helps distribute the plasma along the
magnetic field lines, leading to the generation of a thinner CS compared to
the Res case.

blackRes+TC case

black(a) black(b)

Figure 12: Zoom of the magnetic field lines (a) and the plasma $\beta$ (b) at
time $t=51$ s for the Res+TC case.

Moreover, it is possible to indicate the evolution of the CS’s width, which
could help identify the formation of multiple X- and O-points in the CS, (see,
e.g., Shen et al., 2011). For instance, in Fig. 13, we plot the half-width
$(w)$ of the CS in Mm and the $z-$ component of the current density $J_{z}$ at
the point $(x,y)=$(0,10.4) Mm as functions of time (a) and the vertical
component of the plasma velocity along the $y-$axis at two representative
times $t=43$ s and $t=51$ s (b) for the Res+TC case. Interestingly, the
formation of magnetic islands is related to the width of the CS, and these
structures are expected to form for small values of $w$ (see, e.g., Shen et
al., 2011; MacTaggart & Hillier, 2020). For instance, in Fig. 13(a), we show
the time evolution of the $w$ at the point $(x,y)=(0,10.4)$ Mm represented by
the magenta-colored curve. This figure shows that the CS’s width reduces over
time. Also, in Fig. 13(a), we display the evolution of $J_{z}$ at the same
point as $w$, which is shown by the Mm black curve. One can notice how $J_{z}$
increases with time until around $t\approx 43$ s, and then it experiences a
sudden decrease up to the final time of the simulation, which is $t=51$ s.
This increase is related to the location of an O-type (neutral) point at $y=$
10.4 Mm, which in turn triggers the magnetic island (see, e.g., Shen et al.,
2011). In this scenario, plasma pressure dominates over magnetic pressures
inside the magnetic islands, as shown in the plasma $\beta$ map of Fig. 12.
Since magnetic effects are small, a sudden drop in current density is observed
in the last seconds of the simulation. The appearance of magnetic islands in
the last seconds of the simulation is also consistent with $w$ minimum values
($\sim$ 0.1 Mm). Besides, to identify the formation of six O-points along the
CS, in Fig. 13(b), we display the plasma velocity along the y-direction for
two times $t=43$ s and $t=51$ s. According to our simulation results, these
times cover a period before and after the formation of magnetic islands. Each
arrow in Fig. 13(b) indicates the location of each magnetic island and its
respective vertical velocity, which is normalized with the velocity scaled
factor $v_{0}=10^{8}$ cm s-1. Before forming the magnetic islands, the plasma
vertical velocity increases with time and does not show significant
variations. For example, at $t=51$ s, which is the time when we observe the
formation of magnetic islands, and from $10$ Mm onwards, the plasma vertical
velocity fluctuates. These fluctuations indicate the presence of the magnetic
islands, known as plasma blobs (see, e.g., Shen et al., 2011). In Fig. 13(b),
we also specify the location and vertical velocities of the magnetic islands
themselves. We see an island that moves downstream and is close to $y=10.4$ Mm
with a vertical velocity of about $v_{y}=-0.17$. The other five magnetic
islands move upstream at $y=11.5$ Mm, $y=14.8$ Mm, $y=16.3$ Mm, $y=17.8$ Mm,
and $y=19.6$ Mm. They correspond to vertical velocities of approximately
$v_{y}=0.08$, $v_{y}=0.79$, $v_{y}=0.7$, $v_{y}=0.71$, and $v_{y}=0.76$,
respectively. The plasma vertical velocity around these islands is high,
although for the magnetic island of approximately $y=$16.3 Mm, the vertical
velocity of the plasma around and the island are similar. For other magnetic
islands, the vertical velocity of the surrounding plasma is considerably
higher, which is consistent with the found in, for example, Shen et al.
(2011).

black(a) black(b)

Figure 13: Variations of $w$ (magenta curve), and the $z-$ component of the
current density $J_{z}$ (black curve) at $y=$10.4 Mm as functions of time (a)
and the vertical velocity along the vertical distance $y$ at two times $t=43$
s and $t=51$ s (b) for the Res+TC case.

Additionally, to indicate how the variations of the vertical velocity shown in
Fig. 13(b) are related to the magnetic island formation, in Fig. 14, we show
the temporal evolution of the $z-$ component of the current density $J_{z}$ at
the locations of the first and second magnetic island as schematized by the
arrows of Fig. 13(b). Essentially, these plots show that the behavior of
$J_{z}$ in the two magnetic islands is similar, i.e., all reach their maximum
value of around $t\approx 43$ s and then suddenly decrease to the final of the
simulation. The increase in $J_{z}$ is due to the appearance of magnetic
islands. In the bottom of Fig. 14, we show $J_{z}$ (colored in blue)
calculated at $y=$13.2 Mm, where it is a point without magnetic island. Here,
we note that the behavior of $J_{z}$ is different from the plots where
magnetic islands form since, in the final times of the simulation, $J_{z}$
grows instead of suddenly decreasing.

Figure 14: Temporal evolution of the $z-$component of the current density
$J_{z}$ estimated at the locations of the first and second magnetic islands of
Fig. 14 (black curves) and at $y=$13.2 Mm (blue curve).

##  4 Concluding remarks

In this paper, we study the formation and dynamics of solar flares, including
the post-flare loops in a local region of the solar atmosphere, for two
scenarios, according to the classical 2D model. In the Res case, localized
resistivity produces more substructure around the post-flare loops, related to
the plasma moving upwards and downwards. In addition, we identify the
development of a complex structure, which could be associated with a magnetic
island at the top of post-flare loops. On the other hand, in the Res+TC case,
the whole solar flare structure looks smooth, including the post-flare loops
and no apparent substructures are developed, which means that thermal
conductivity redistributes the plasma homogeneously. However, we observe the
formation of multiple small-scale magnetic islands from the top of the post-
flare loops and along the CS in the vertical direction. Additionally, we
include the case of an anomalous resistivity without thermal conduction. In
that scenario, any post-flare loops and thin CS are formed. Instead, we see a
hot plasma structure forming that becomes wider and moves upwards while
magnetic field lines reconnect at the coronal level. This kind of plasma
feature could be related to a plasmoid. Furthermore, in the mass density maps,
we identify the development of a small jet ($\sim 3$ Mm of height) moving
upwards from the transition region. This feature is interesting since it could
be related to the magnetic reconnection process at microscopic scales, which
is triggered by anomalous resistivity, given a relevant result despite the
limitations of the MHD to capture the microscopic mechanism associated with
the magnetic reconnection process. However, it is clear that in the case with
anomalous resistivity, we do not observe the formation of a typical solar
flare and associated post-flare loops as indicated by the classical 2D model
introduced by Yokoyama & Shibata (2001).

The development of the multiple magnetic islands in the Res+TC case is related
to the Tearing instability in the nonlinear regime. This instability is
characterized by making a system in equilibrium unstable due to the tearing
mode, leading to the formation of X-points and plasmoids during the
reconnection process. This instability is part of a fundamental resistive
instability linked with a long-wave tearing mode, corresponding to the breakup
of the layer along current-flow lines. Besides, in the Res+TC case, the
magnetic islands form along the CS in the vertical direction, consistent with
the results obtained by Wang et al. (2021). Moreover, we observe that magnetic
islands tend to form in regions with high plasma $\beta$ ($\sim 10^{4}$) but
surrounded by low plasma $\beta$ regions ($\sim 10^{-2}$), which suggest that
on those regions plasma pressure and magnetic pressure must contribute. That
is, the magnetic pressure confines the plasma, and the plasma pressure makes
the structure of the islands. In the magnetic islands, the current density
tends to increase, and the CS’s half-width becomes thinner when forming.
Finally, when magnetic islands form, we identify the generation of plasma bobs
in the vertical velocity, which appears due to the flow generated by the
magnetic reconnection process.

In the case of the analysis for the instabilities performed to identify the
development of RTI and RMI, we conclude that in the Res case, the internal
spike-shaped structures are likely related to an RMI since the general
physical conditions that state this instability are satisfied. In the Res+TC
case, we show that any condition for RMI or RTI is satisfied since the post-
flare loops look smooth. Although our simulation setup for the Res case is
more straightforward than the one used in Shen et al. (2022), we identify
desirable conditions for developing the RMI/RTI instabilities since we use
only a localized resistivity profile. This result is crucial since it could
indicate that a simple model with magnetic resistivity can solely produce such
instabilities observed in underdense plasma downflows associated with magnetic
reconnection in solar flares.

A difference of the present study regarding, for example, Yokoyama & Shibata
(2001); Takasao et al. (2015); Ruan et al. (2020); Shen et al. (2022), is the
inclusion of a highly anisotropic thermal conductivity that smoothly varies
between the classical and saturated thermal conduction regimes, which
according to our knowledge, it had not been used in numerical simulations of
solar flares, so far. Additionally, we adopt a model for the solar atmosphere
that includes the minimum temperature of the chromosphere and the sharp
gradient of the transition region, which is more realistic than the one
adopted by the abovementioned works since they typically consider a simple
tangent function to describe a smooth variation of temperature from the
photosphere to the solar corona. Furthermore, in comparison to Mondal et al.
(2024), where the authors studied the formation and evolution of plasmoids in
coronal CS, subjected to an external velocity perturbation under the condition
of uniform resistivity, in our paper, we employed a localized resistivity that
triggers the magnetic reconnection and leads the formation of the CS. Besides,
Mondal et al. (2024) indicates that CS initially goes in thinning and
Petschek-type magnetic reconnection followed by tearing instability and
plasmoid formation. In contrast, one of the main objectives of our paper
focuses on studying the associated small-scale structures of solar flares
without entering into detail about the reconnection mechanism. However, both
analyses could complement each other and engage in the context of the highly
dynamic behavior observed along the CS in flare-like structures.

Finally, the results of this paper could help to continue persuading the
dynamics of the formation of solar flares based on the standard model. In
particular, it opens the gate to investigate if the anomalous resistivity
itself could generate plasmoids and jets, which complete the picture of the
magnetic reconnection process using the MHD approximations. Additionally, the
current results could help better determine that solar flares can develop
substructures even in simple scenarios where, for example, turbulence and
particle acceleration are not considered.

## Acknowledgements

We gratefully acknowledge the constructive remarks of the referee that highly
improved our manuscript. J.J.G.-A. acknowledges to "Consejo Nacional de
Humanidades Ciencias y Tecnologías (CONAHCYT)" 319216 project "Modalidad:
Paradigmas y Controversias de la Ciencia 2022." M.G.-S. is grateful to
“Programa de Becas para Estudios de Posgrado en México” granted by CONAHCYT.

## Data Availability

The data underlying this article will be shared on reasonable request to the
corresponding author.

## References

  * Alexiades et al. (1996) Alexiades V., Amiez G., Gremaud P.-A., 1996, Communications in Numerical Methods in Engineering, 12, 31
  * Avrett & Loeser (2008) Avrett E. H., Loeser R., 2008, The Astrophysical Journal Supplement Series, 175, 229
  * Bárta et al. (2011) Bárta M., Büchner J., Karlický M., Skála J., 2011, The Astrophysical Journal, 737, 24
  * Cécere et al. (2015) Cécere M., Zurbriggen E., Costa A., Schneiter M., 2015, The Astrophysical Journal, 807, 6
  * Chen et al. (2015) Chen B., Bastian T. S., Shen C., Gary D. E., Krucker S., Glesener L., 2015, Science, 350, 1238
  * Dedner et al. (2002) Dedner A., Kemm F., Kröner D., Munz C. D., Schnitzer T., Wesenberg M., 2002, Journal of Computational Physics, 175, 645
  * Forbes (2010) Forbes T., 2010, in Schrijver C. J., Siscoe G. L., eds, , Heliophysics: Space Storms and Radiation: Causes and Effects. p. 159
  * Furth et al. (1963) Furth H. P., Killeen J., Rosenbluth M. N., 1963, The Physics of Fluids, 6, 459
  * Giuliani (1984) Giuliani J. L. J., 1984, The Astrophysical Journal, 277, 605
  * González-Avilés et al. (2017) González-Avilés J. J., Guzmán F. S., Fedun V., 2017, The Astrophysical Journal, 836, 24
  * Guo et al. (2014) Guo L.-J., Huang Y.-M., Bhattacharjee A., Innes D. E., 2014, The Astrophysical Journal, 796, L29
  * Kliem et al. (2000) Kliem B., Karlický M., Benz A. O., 2000, A%A, 360, 715
  * Lin et al. (2005) Lin J., Ko Y.-K., Sui L., Raymond J. C., Stenborg G. A., Jiang Y., Zhao S., Mancuso S., 2005, The Astrophysical Journal, 622, 1251
  * Liu et al. (2013) Liu W., Chen Q., Petrosian V., 2013, The Astrophysical Journal, 767, 168
  * MacTaggart & Hillier (2020) MacTaggart D., Hillier A., 2020, Topics in Magnetohydrodynamic Topology, Reconnection and Stability Theory. Springer
  * Masuda et al. (1994) Masuda S., Kosugi T., Hara H., Tsuneta S., Ogawara Y., 1994, Nature, 371, 495
  * Mignone et al. (2007) Mignone A., Bodo G., Massaglia S., Matsakos T., Tesileanu O., Zanni C., Ferrari A., 2007, The Astrophysical Journal Supplement Series, 170, 228
  * Miyagoshi & Yokoyama (2004) Miyagoshi T., Yokoyama T., 2004, The Astrophysical Journal, 614, 1042
  * Miyoshi & Kusano (2005) Miyoshi T., Kusano K., 2005, Journal of Computational Physics, 208, 315
  * Mondal et al. (2024) Mondal S., Srivastava A. K., Pontin D. I., Yuan D., Priest E. R., 2024, arXiv e-prints, p. arXiv:2401.07048
  * Priest (2014a) Priest E., 2014a, Magnetohydrodynamics of the Sun, doi:10.1017/CBO9781139020732.
  * Priest (2014b) Priest E., 2014b, References. Cambridge University Press, p. 493–542, doi:10.1017/CBO9781139020732.017
  * Ruan et al. (2020) Ruan W., Xia C., Keppens R., 2020, The Astrophysical Journal, 896, 97
  * Shen et al. (2011) Shen C., Lin J., Murphy N. A., 2011, The Astrophysical Journal, 737, 14
  * Shen et al. (2022) Shen C., Chen B., Reeves K. K., Yu S., Polito V., Xie X., 2022, Nature Astronomy, 6, 317
  * Shiota et al. (2005) Shiota D., Isobe H., Chen P. F., Yamamoto T. T., Sakajiri T., Shibata K., 2005, The Astrophysical Journal, 634, 663
  * Su et al. (2013) Su Y., Veronig A. M., Holman G. D., Dennis B. R., Wang T., Temmer M., Gan W., 2013, Nature, 9, 489
  * Tajima & Shibata (2002) Tajima T., Shibata K., 2002, Plasma astrophysics
  * Takasao et al. (2015) Takasao S., Matsumoto T., Nakamura N., Shibata K., 2015, The Astrophysical Journal, 805, 135
  * Ugai (1992) Ugai M., 1992, Physics of Fluids B: Plasma Physics, 4, 2953
  * Wang et al. (2021) Wang Y., Cheng X., Ding M., Lu Q., 2021, The Astrophysical Journal, 923, 227
  * Warren et al. (2018) Warren H. P., Brooks D. H., Ugarte-Urra I., Reep J. W., Crump N. A., Doschek G. A., 2018, The Astrophysical Journal, 854, 122
  * Yokoyama & Shibata (2001) Yokoyama T., Shibata K., 2001, The Astrophysical Journal, 549, 1160
  * Zhou et al. (2021) Zhou Y., et al., 2021, Physica D: Nonlinear Phenomena, 423, 132838