# Comparative simulations of Kelvin-Helmholtz induced magnetic reconnection at
the Earth’s magnetospheric flanks

Silvia Ferro<EMAIL_ADDRESS>Centre for mathematical Plasma
Astrophysics, Department of Mathematics, KU Leuven, Celestijnenlaan 200B,
B-3001 Leuven, Belgium  Matteo Faganello  Aix-Marseille University, CNRS, PIIM
UMR 7345, Marseille, France  Francesco Califano  Dipartimento di Fisica “E.
Fermi”, Università di Pisa, Pisa, Italy  Fabio Bacchini  Centre for
mathematical Plasma Astrophysics, Department of Mathematics, KU Leuven,
Celestijnenlaan 200B, B-3001 Leuven, Belgium  Royal Belgian Institute for
Space Aeronomy, Solar-Terrestrial Centre of Excellence, Ringlaan 3, 1180
Uccle, Belgium

###### Abstract

This study presents 3D resistive Hall-magnetohydrodynamic numerical
simulations focusing on the Kelvin-Helmholtz instability (KHI) dynamics at
Earth’s magnetospheric flanks during northward interplanetary magnetic field
(IMF) periods. Through a comparative analysis of two simulations, we explore
the impact of distinct magnetic field orientations on plasma dynamics and
magnetic reconnection events. In one configuration, a uniform magnetic field
results in (double) Mid-Latitude Reconnection (MLR), while in the second
configuration, a magnetic shear induces both Type I Vortex Induced
Reconnection and MLR. Key findings include the symmetries and asymmetries in
the latitudinal distribution of KHI vortices and current sheets and a
quantitative comparison of reconnection events. Particularly noteworthy is the
quantification of newly closed field lines that experienced double
reconnection, ultimately becoming embedded in solar wind plasma at low
latitudes while remaining connected to magnetospheric plasma at high
latitudes. The varying abundance of such lines in the two simulations holds
implications for plasma transport at the magnetopause. This work significantly
advances our understanding of magnetospheric processes, emphasizing the
essential role of a three-dimensional perspective in accurately simulating
magnetospheric plasma phenomena.

††preprint: AIP/123-QED

##  I Introduction

Magnetic reconnection is a fundamental mechanism capable of altering the
topology of magnetic field lines within plasmas. Specifically, it plays a
significant role within the Earth’s magnetosphere (MS). Reconnection can
transfer the energy stored in the global, large-scale structure of the
magnetic field to kinetic scales, and it releases this energy locally while
also modifying the line connectivity. The classic magnetospheric example is
provided by the Dungey cycle[1] during periods of southward Interplanetary
Magnetic Field (IMF): IMF lines and magnetic flux, advected by the solar wind
(SW), pile up in front of the magnetopause, where the electric current is
strongly enhanced. At scales comparable with the electron inertial length, the
frozen-in condition breaks, and reconnection occurs[2]. Thus, magnetic
reconnection modifies the topology of the magnetospheric field and that of the
SW, creating “open” field lines connected to the Earth at one pole only. This
global topological modification of the MS allows for plasma recirculation and
plasma transport across the magnetopause. Indeed, open field lines can be
advected in the nightside region, accumulating magnetic flux along the
magnetotail, where finally magnetic reconnection occurs once again [3].
Furthermore, SW particles can easily stream along open field lines and enter
the MS[4, 5].

During northward periods, dayside reconnection is ineffective. On the
contrary, Cusp reconnection[6] can occur, causing plasma transport at dayside
and generating a low-latitude boundary layer[7]. At the same time, the Kelvin-
Helmholtz instability (KHI), active along the low-latitude MS flanks, mediates
the interaction between the SW and the MS and drives reconnection there,
boosting plasma transport. Two main mechanisms are at play in that context,
depending on the shear of the magnetic field at low latitudes.

If the IMF and the magnetospheric field are not aligned, magnetic shear
arises. In this configuration, the KHI forces reconnection to occur with a
mechanism somewhat similar to dayside reconnection: the velocity field imposed
by the growing, ideal KHI accumulates magnetic flux at the pre-existing
magnetopause current sheet until reconnection occurs. The pinching of the low-
latitude current sheet of the magnetopause is thus the peculiar signature of
Type I Vortex Induced Reconnection (VIR)[8, 9].

In the second scenario, when the magnetospheric field and the IMF are
sufficiently aligned at low latitudes, current sheets will not develop in
those regions. Consequently, the KHI will create current sheets elsewhere by
folding and/or stretching magnetic lines. In the oversimplified yet widely
considered two-dimensional picture, KH vortices can fold the in-plane (with
respect to the 2D dynamical plane) component of the initial field, thereby
leading to the development of current sheets and Type II VIR[10, 11]. In a
more realistic three-dimensional configuration, field lines are tied to high-
latitude regions, less prone to KH development. As a result, lines are
advected differently at low and high latitudes, and current sheets are created
far away from the equatorial plane, at mid-latitudes, where finally, mid-
latitude reconnection (MLR) takes place[12].

The topological modifications generated in the two scenarios above are
different. When starting from a configuration with aligned magnetic fields and
without an initial current sheet at low latitudes, simulations show that MLR
proceeds almost simultaneously in both hemispheres (double-MLR). This leads to
the formation of “twice-reconnected lines” [13, 14], half of which are “newly
closed”. Newly closed field lines are connected to the Earth at both poles but
are embedded in the SW plasma at low latitudes. Type I VIR, instead, creates
open field lines that connect the two plasma regions. Both mechanisms have
been shown to contribute substantially to the transport properties of the
flank magnetopause[15, 5, 12]. Moreover, double-MLR can explain the increase
of specific entropy inside the MS during northward periods[16]. On the
contrary, as soon as a magnetic shear exists, Type I VIR can occur, but not
alone. In fact, while Type I VIR is forced by the pinching of a low-latitude
current sheet, where the vortices are well-developed, 3D effects related to
the differential advection of field lines at low and high latitudes induce
MLR. These complex dynamics have been inferred by MMS data[17, 18] and
confirmed by numerical simulations[19, 20]. In this way, a mix of open field
lines and newly closed lines can be found at the magnetopause.

Here, with the aid of numerical simulations, we perform a comparative study to
explore the two potential 3D evolution scenarios observed during northward
periods. The two configurations vary in the orientation of the magnetic field,
with one featuring a purely northward magnetic field and the other introducing
a magnetic rotation and, consequently, a shear in the initial magnetic field.
Within these nearly identical settings, we can, for the first time,
quantitatively assess the impact of magnetic rotation on magnetic reconnection
events in such systems. Our focus is particularly on quantifying the amount of
reconnected lines generated in both configurations, determining their
implications for plasma transport at the magnetopause, and evaluating the
relative proportions of open and newly closed magnetic field lines.

This article is structured as follows: In section II, we introduce the plasma
models, starting equilibria, and the numerical method. In section III we
compare simulation results, and section IV is dedicated to conclusions.

##  II Plasma Model and Simulation Set up

###  II.1 Resistive Hall-MHD model and numerical code

The evolution of the plasma is described by a resistive Hall-MHD quasi-neutral
model[12], which accounts for different evolution for the electron and ion
pressures and diamagnetic effects for the electrons. The model equations read

$\displaystyle\frac{\partial n}{\partial t}+\nabla\cdot(n\pmb{\mathit{u}})=0,$
(1) $\displaystyle\frac{\partial(n\pmb{\mathit{u}})}{\partial
t}+\nabla\cdot[(n\pmb{\mathit{u}}\pmb{\mathit{u}})+(\Pi\overleftrightarrow{\pmb{\mathit{I}}}-\pmb{\mathit{B}}\pmb{\mathit{B}})]=\mathbf{0},$
(2)

where $n$ is the overall plasma number density, $\pmb{\mathit{u}}$ the fluid
velocity (neglecting electron inertia), $\pmb{\mathit{B}}$ the magnetic field,
and $\Pi=P_{\text{th}}+P_{\text{m}}$ is the total pressure (sum of thermal
pressure $P_{\text{th}}=P_{i}+P_{e}$ and magnetic pressure
$P_{\text{m}}=B^{2}/2$). The thermal pressures of ions and electrons evolve
following an adiabatic closure (for each species $\alpha=i,e$),

$\frac{\partial(n_{\alpha}S_{\alpha})}{\partial
t}+\nabla\cdot(n_{\alpha}S_{\alpha}\pmb{\mathit{u}}_{\alpha})=0\,\,,\quad
S_{\alpha}=P_{\text{th}}n^{-5/3}.$ (3)

To complete the set of equations, we introduce Faraday’s and Ampère’s laws to
describe the evolution of the magnetic field (neglecting the displacement
current, given the non-relativistic regime of the magnetospheric plasma),
along with the evolution of the current density $\pmb{\mathit{j}}$,

$\displaystyle\frac{\partial\pmb{\mathit{B}}}{\partial
t}=-\nabla\crossproduct\pmb{\mathit{E}}\,,\quad\pmb{\mathit{j}}=\nabla\crossproduct\pmb{\mathit{B}}=n(\pmb{\mathit{u}}_{i}-\pmb{\mathit{u}}_{e}).$
(4)

Here, $\pmb{\mathit{u}}_{i}$ and $\pmb{\mathit{u}}_{e}$ are the ion and
electron fluid velocities, with $\pmb{\mathit{u}}_{i}\simeq\pmb{\mathit{u}}$
and $\pmb{\mathit{u}}_{e}\simeq\pmb{\mathit{u}}_{i}-\pmb{\mathit{j}}/n$ in the
limit of massless electrons. All quantities are normalized using the ion mass
$m_{i}$, the ion skin depth $d_{i}$, and the ion cyclotron frequency
$\Omega_{\text{c}}$. Finally, the generalized Ohm’s law accounting for the
electric field $\pmb{\mathit{E}}$ reads

$\displaystyle\pmb{\mathit{E}}=-\pmb{\mathit{u}}\times\pmb{\mathit{B}}+\frac{\pmb{\mathit{j}}\crossproduct\pmb{\mathit{B}}}{n}-\frac{1}{n}\nabla
P_{e}+\eta\pmb{\mathit{j}}.$ (5)

The introduction of an artificial, small resistivity $\eta\sim 10^{-3}$ in
Equation (5) ultimately enables the occurrence of magnetic reconnection. With
our choice of $\eta$, justified by earlier 2D studies[21, 22], the resistive
diffusion time-scale is $\tau_{R}\sim L_{\text{eq}}^{2}\eta^{-1}\approx
10^{4}$, where $L_{\text{eq}}=3d_{i}$ is the gradient length of the
equilibrium configuration. This time scale is far longer than our simulation
time. Consequently, at large scales, the magnetic field is frozen-in in the
fluid motion generated by the KH dynamics, and inside the vortices magnetic
diffusion is negligible. Furthermore, the Hall term plays a critical role in
$\beta\sim 1$ plasmas for describing the decoupling between the magnetic and
the ion dynamics at nearly the ion skin depth. Most importantly, taking into
account the Hall term enables the transition to fast reconnection
characterized by a reconnection rate weakly dependent on the value of the
resistivity[23].

We use a fourth-order Runge-Kutta scheme for the numerical integration of
Equations (1)–(5) and calculate spatial derivatives using sixth-order explicit
finite differences in the $y$\- and $z$-direction. In the $x$-direction, we
implement a sixth-order implicit compact scheme with spectral-like
resolution[24], as this is the direction in which smaller scale dynamics
emerge. We impose periodic conditions at the $y$\- and $z$-boundary. The
$x$-boundaries are open, so we apply conditions based on magnetohydrodynamic
characteristics, allowing Alfvénic and magnetosonic perturbations to freely
exit the numerical domain[25, 26].

###  II.2 Geometrical configuration and equilibrium

We adopt a 3D slab configuration to schematically model the structure of
Earth’s magnetopause flank, including 3D mechanisms such as the stabilization
of the KHI at high latitudes. In our frame of reference, the $x$-axis is
perpendicular to the undisturbed magnetopause, the $y$-axis aligns with the SW
direction, and the $z$-axis corresponds to the northward direction. In order
to track the evolution of the KH-induced vortices, the numerical box moves at
the phase velocity of the KHI along the MS flank towards the tail region,
effectively remaining fixed in the rest frame of the vortex.

For the initial configuration, we consider a 2D Grad-Shafranov equilibrium
satisfying

$\nabla^{2}\psi=-4\pi\frac{\mathrm{d}\Pi}{\mathrm{d}\psi},$ (6)

where $\psi$ is the flux function. The corresponding equilibrium magnetic
field reads

$\pmb{\mathit{B}}_{\text{eq}}(x,z)=\nabla\psi\times\hat{y}+B_{y,\text{eq}}(\psi)\hat{y}.$
(7)

In Equation (7), the first term represents the (nearly) northward MS field.
The second term $B_{y,eq}$ corresponds to the flow-aligned component of the
IMF, which reproduces the overall magnetic field rotation when combined with
the northward magnetic field component. All other physical quantities are
obtained from the flux function and are by construction constant along
magnetic field lines. Setting $\Pi$ constant with respect to $\psi$, the flux
function can be written as a harmonic function

$\displaystyle\psi$
$\displaystyle=\frac{1}{2}\biggl{[}\left(1+\frac{1}{\delta}\right){\left(x-\frac{L_{x}}{2}\right)}$
(8)
$\displaystyle+\left(1-\frac{1}{\delta}\right)\frac{L_{z}}{2\pi}\operatorname{sinh}\frac{2\pi{\left(x-L_{x}/2\right)}}{L_{z}}\operatorname{cos}\frac{2\pi(z-L_{z}/2)}{L_{z}}\biggr{]}.$

The resulting equilibrium quantities are given by

$\displaystyle B_{x,\text{eq}}$
$\displaystyle=\frac{1}{2}\left(1-\frac{1}{\delta}\right)\operatorname{sinh}\frac{2\pi{\left(x-L_{x}/2\right)}}{L_{z}}\operatorname{sin}\frac{2\pi(z-L_{z}/2)}{L_{z}},$
(9) $\displaystyle B_{y,\text{eq}}$
$\displaystyle=\operatorname{tan}\frac{\phi}{2}\operatorname{tanh}\left(\frac{\psi}{L_{\text{eq}}}\right),$
$\displaystyle B_{z,\text{eq}}$
$\displaystyle=\frac{1}{2}\biggl{[}\left(1+\frac{1}{\delta}\right)$
$\displaystyle+\left(1-\frac{1}{\delta}\right)\operatorname{cosh}\frac{2\pi{\left(x-L_{x}/2\right)}}{L_{z}}\operatorname{cos}\frac{2\pi(z-L_{z}/2)}{L_{z}}\biggr{]},$
$\displaystyle u_{y,\text{eq}}$ $\displaystyle=\frac{\Delta
u}{2}\operatorname{tanh}\frac{\psi}{L_{\text{eq}}},\qquad\text{and}\qquad
u_{x,\text{eq}}=u_{y,\text{eq}}\,=0.$

The equilibrium number density is constant, $n_{\text{eq}}=1$. The initial
thermal pressure is set so that $\Pi$ is constant, with
$P_{\text{th}}|_{\psi=0}=1$ and $P_{i}/P_{e}=4$, as inspired by earlier
works[27, 20].

Table 1:  Summary of the relevant parameters characterizing the numerical boxes and the equilibrium configurations: the box dimensions ($L_{x},L_{y},L_{z}$), the number of grid points ($N_{x},N_{y},N_{z}$), the equilibrium scale length $L_{\text{eq}}$, the magnetic shear angle of both runs $\phi_{\text{A,B}}$, the velocity jump across the two plasma regions $\Delta u$, and the parameter $\delta$ determining the ratio between $B_{z}(x=0,z=L_{z}/2)$ and $B_{z}(x=0,z=0)$. parameter | value | parameter | value
---|---|---|---
$L_{x}$ | $90$ | $N_{x}$ | $600$
$L_{y}$ | $94$ | $N_{y}$ | $512$
$L_{z}$ | $377$ | $N_{z}$ | $512$
$\delta$ | $3$ | $L_{\text{eq}}$ | $3$
$\Delta u$ | $1$ | $\phi_{\text{A,B}}$ | $0^{\circ},\ 18^{\circ}$

Furthermore, the chosen value $\delta\ll L_{x}<L_{z}/4$ ensures that the
dominant component of the magnetic field in the $(x,z)$-plane is aligned with
the northward direction. Setting $L_{\text{eq}}=3d_{i}$ allows the KHI to
develop at large, nearly MHD scales, with the wavelength of the fastest
growing mode (FGM) measuring approximately $\lambda_{\text{FGM}}\sim 47d_{i}$.
Finally, we set $L_{y}=2\lambda_{\text{FGM}}$ in order to let multiple
vortices develop in the numerical box. Since $L_{z}>8\lambda_{\text{FGM}}$,
stable high-latitude regions are far enough from the unstable equatorial one
so that KH vortices can grow and even pair[28].

In Table 1, we list the main parameters of the numerical box and the
equilibrium configuration. We have performed two identical runs differing only
by the presence (or absence) of shear in the initial magnetic field
configuration. This was achieved using the magnetic shear angle $\phi$ (see
$B_{y,\text{eq}}$ in Equation (9)), which represents the total angle between
the MS and the SW magnetic fields. In Figure 1, we show in the left panel the
results of the first simulation, “Run A”, in which $\phi_{\text{A}}=0$. This
implies no rotation in the magnetic field, hence $B_{y,\text{eq}}=0$. The
right panel shows the second simulation, “Run B”, with $\phi_{\text{B}}\simeq
18^{\circ}$. This introduces a magnetic rotation at equilibrium and,
therefore, the presence of an initial magnetic shear in the system. A
consequence of this initial setting is that at $t=0$, in Run B, there is a
small initial current sheet (absent in Run A) at the frontier between the two
plasma regions. The chosen velocity profile (see Equation (9)) is the same for
both runs and has stronger gradients at the equator and weaker at high
latitudes. This is shown in Figure 2, where we plot a slice at $y=0$ through
the numerical box for Run A. In particular, we show the hourglass-like
configuration of the initial magnetic field lines (which correspond to the
isocontours of the magnetic flux) in the $(x,z)$-plane, with lines getting
closer together at the equatorial plane ($z\sim L_{z}/2$) than at high
latitudes ($z\sim 0$ and $z\sim L_{z}$). Since the growth rate of the
instability $\gamma_{KH}$ is proportional to the velocity gradient[29], this
equilibrium drives the development of the KHI far earlier in the equatorial
plane, while high latitude regions remain “quiet” throughout the simulation.
This configuration effectively captures the MS flank dynamics. The choice of a
uniform density enables us to investigate the nonlinear evolution of the
instability in an idealized configuration, ruling out the development of
secondary Rayleigh-Taylor instabilities[30, 31, 32]. For the parameters we
adopt, the FGM at the equator has $m=2$, where $m$ is the mode number along
the $y$-direction. Except for the density, the chosen parameters (e.g., the
ion skin depth $d_{i}$, Alfvén velocity $v_{\text{A}}$, and the other
normalization quantities) represent typical values of physical quantities
observed in the outer MS, with sonic and alfvénic Mach numbers of order one.
For practical reasons, we define a scalar quantity $\zeta$ that is passively
advected by the plasma velocity, following $(\partial/\partial
t+\pmb{\mathit{u}}\cdot\nabla)\zeta=0$, with
$\zeta=0.6+0.4\operatorname{tanh}(\psi/L_{\text{eq}})$. This passive tracer
helps us visualize plasma structures even in the absence of a density jump
across the magnetopause. Indeed, the MS plasma is characterized by values
$\zeta<0.6$, while the SW plasma has $\zeta>0.6$, with the magnetopause
surface corresponding to $\zeta_{\text{fr}}=0.6$.

Figure 1:  Sketch of the initial magnetic field configuration. In Run A, the
magnetic shear angle $\phi_{\text{A}}=0$; therefore, both the MS and SW
magnetic fields are aligned northward. In Run B, the presence of
$\phi_{\text{B}}\simeq 18^{\circ}$ introduces a rotation: at the magnetopause
(the frontier) the magnetic fields are aligned with the northward direction,
but both fields are inclined by $\pm\phi_{\text{B}}/2$ for $x\rightarrow
0,L_{x}$. Figure 2:  Isocontours of the magnetic flux function $\psi(x,z)$
(white solid lines) and the spatial distribution of the ion velocity
$U_{iy}(x,z)$ at equilibrium ($t=0$). The white line at $x=L_{x}/2$ represents
the frontier between the MS and the SW plasmas at equilibrium.

##  III Results and simulations comparison

###  III.1 Asymmetrical VS symmetrical development of vortices and
reconnection events

Figure 3:  Panels (A1) and (B1): 3D representation of the current sheets and
KH vortices forming respectively at $t=460\Omega_{\text{c}}^{-1}$ for Run A
and at $t=500\Omega_{\text{c}}^{-1}$ for Run B. (A1): Slices at latitudes
$z=L_{z}$, $z=L_{z}/2$, and $z=0$, showing the latitudinal evolution of
vortices through the passive tracer. Green current sheets form around mid-
latitudes. Two slices ($z=3L_{z}/4$ top, $z=L_{z}/4$ bottom) in panels A1 and
A2 show the magnitude of the current density $|\pmb{\mathit{j}}|$. Panel (B1):
Slices at latitudes $z=L_{z}$, $z=115d_{i}$, and $z=0$ show the latitudinal
evolution of vortices through the passive tracer. Green current sheets form in
the northern hemisphere. In panels B1 and B2, we show the slice at
$z=240d_{i}$ with the current density magnitude $|\pmb{\mathit{j}}|$. (A3–4)
and (B3–4): latitudinal distribution of the maximum of the current density
$\mathrm{max}(|\pmb{\mathit{j}}|)$ and the ion velocity
$\mathrm{max}(|\pmb{\mathit{u}}_{i}|)$ respectively at
$t=200,380,460\Omega_{\text{c}}^{-1}$ for Run A and
$t=320,400,500\Omega_{\text{c}}^{-1}$ for Run B. The maximum values of
$|\pmb{\mathit{j}}|$ and $|\pmb{\mathit{u}}_{i}|$ are taken over the whole
$(x,y)$-plane at each $z$.

In Figure 3, we show a 3D rendering of the nonlinear evolution of Run A and B
(panels A1 and B1, respectively). The 3D boxes in this and the following
figures have been doubled along the $y$ periodic direction ($L_{y}\rightarrow
2L_{y}$) for the sake of clearness. First, we note that high latitudes remain
stable, as shown by the lack of structures at the planes $z=0,L_{z}$, while
folded vortices have formed around the equatorial region in both simulations.
Second, while in Run A vortices develop almost symmetrically with respect to
the equatorial plane, in particular, they develop at around $z=L_{z}/2$, Run B
exhibits an asymmetric evolution with more developed vortices in the southern
hemisphere, in particular at around at $z\simeq 115d_{i}<L_{z}/2$. Data from
Run A and Run B are plotted at different times, $t=460\Omega_{\text{c}}^{-1}$
and $t=500\Omega_{\text{c}}^{-1}$ respectively, but corresponding to the same
phase of the instability’s evolution. This apparent mismatch between the
instability stage and the simulation time is due to the different growth rates
for the two runs because of the stabilizing effect of the magnetic shear on
the KHI in Run B. To compute the growth rate of the instability at the
equatorial plane, we can approximate our system to a slab configuration of two
uniform inviscid magnetized fluids in relative motion. The dispersion
relation[33, 34] suggests that the growth $\gamma_{\text{B}}$ for Run B should
be slightly smaller than $\gamma_{\text{A}}$, the one for Run A. Indeed,
during the early development of the KHI, we infer from simulation data
$\gamma_{\text{A}}=0.0112$ and $\gamma_{\text{B}}=0.0088$. Given that the
folding time of the vortices is inversely proportional to the growth rate of
the KHI[35], Run B takes longer to reach the folding phase of the instability.

The symmetric and asymmetric evolution in the two runs is also visible in the
development of the current sheets. In the symmetric case, the formation of
current sheets is driven by the sole differential advection of magnetic field
lines [12]. In Run A, we show that magnetic field lines anchored at high
latitudes in the unperturbed plasma undergo counterstreaming flows and,
simultaneously, are advected at low latitudes by the growing KH vortices.
Since vortices are nearly at rest in our frame, lines belonging to different
plasma regions are bent in opposite directions. This process results in the
mostly symmetric formation of current sheets, highlighted as green sheets with
peaks (reddish regions) near mid-latitudes in panel (A1). At those latitudes
($z=L_{z}/4,3L_{z}/4$), the slices in panel (A2) display the magnitude of the
current density, showing the current sheets between vortices.

Panels (B1–2) of Figure 3 refer to Run B at $t=500\Omega_{\text{c}}^{-1}$. In
panel (B1), current sheets lie primarily in the northern hemisphere, with the
green highlighted region expanding across the upper part of the numerical box
asymmetrically with respect to the equatorial plane. The peaks in the current
sheets (reddish regions) are located around $z=240d_{i}>L_{z}/2$, where the
horizontal slice in panel (B2) displays the magnitude of the current density.
In this case, current formation is due to a combined action. Current pinching,
imposed by the growing KHI, raises the original current intensity. At the same
time, differential advection modifies the original magnetic rotation by
arching magnetic field lines, thus increasing the current in the northern
hemisphere while decreasing rotation and current in the southern one. The
hemisphere where the current is enhanced is determined by the scalar product
between the equilibrium current and the equilibrium vorticity[36]:
${\bm{\omega}}_{\text{eq}}\cdot{\pmb{\mathit{j}}}_{\text{eq}}$, where
${\bm{\omega}}_{\text{eq}}=\nabla\times{\pmb{\mathit{u}}}_{\text{eq}}$ and
${\pmb{\mathit{j}}}_{\text{eq}}=\nabla\times{\pmb{\mathit{B}}}_{\text{eq}}$
are the equilibrium vorticity and current, respectively.

The symmetrical development of current sheets and vortices with respect to low
latitudes for Run A is further displayed in Figure 3 panels (A3) and (A4). The
latitudinal distribution of the current density magnitude in panel (A3)
exhibits two peaks, which develop with the beginning of the nonlinear phase of
the KHI and grow as the simulation proceeds. Instead, the latitudinal
distribution of the ion velocity in panel (A4) has a single peak around the
equator ($z=L_{z}/2$), where the KH vortices fully develop. On the other hand,
the simulation with an initial magnetic shear displays an asymmetric evolution
with respect to the equatorial plane, as shown in Figure 3 panels (B3) and
(B4). The latitudinal distribution of the magnitude of the current density
(B3) starts with a single peak at $t=320Omega_{\text{c}}^{-1}$ due to the
pinching of the equilibrium current sheet induced by the velocity field
imposed by the growing KHI. This peak then shifts toward the northern
hemisphere ($z>L_{z}/2$) since differential advection increases the magnetic
rotation there. A second, smaller peak develops at later times
($t=500\Omega_{\text{c}}^{-1}$) in the southern hemisphere at $z\simeq
80d_{i}<L_{z}/2$. Indeed, after differential advection has lowered the
magnetic angle $\phi_{\text{B}}$, the latter changes sign and creates local
current sheets with $\pmb{\mathit{j}}$ oppositely directed with respect to the
equilibrium current. Looking at the evolution of KH vortices, the ion velocity
latitudinal distribution in panel (B4) shows a single peak that moves toward
the southern hemisphere. Consequently, there is a southward latitudinal shift
of the vortices as the KHI is inhibited in the northern hemisphere. Indeed,
the larger the magnetic rotation ($\phi_{\text{B}}$ increases in the northern
hemisphere), the larger the impact of the magnetic tension on the KHI[33, 34].

Figure 4: Slices of Run A at the equator (top row) and Run B below the equator
in the southern hemisphere (bottom row) showing the temporal evolution of the
KH vortices through the passive tracer $\zeta$ that highlights the
magnetosheath plasma (in red) and the magnetospheric plasma (in blue).

As the development of the KHI dominates the plasma dynamics in both
simulations, large-scale vortex structures form, as shown in Figure 4. The
passive tracer enables us to follow the morphological evolution of the
vortices for Run A (top row) and Run B (bottom row). For both configurations,
we show the different phases of the evolution: several vortices fully form
(left column) and start pairing (middle column) but get disrupted by secondary
instabilities before the full pairing (right column). In Figure 3, we have
shown that vortices are enhanced at different latitudes, depending on the
initial magnetic rotation. For this reason, in Figure 4, the horizontal slices
($(x,y)$-planes) are taken at $z=L_{z}/2\simeq 188d_{i}$ for Run A and at
$z\simeq 105d_{i}$ for Run B, i.e. at the equator and slightly below it,
respectively.

###  III.2 Quantitative analysis of reconnecting lines

To obtain quantitative insight into the reconnection process brought about by
KH mixing, we measure the number of reconnected lines that are generated in
both configurations by tracing $\sim 22000$ magnetic field lines in both runs
at each time of the simulations. We start by identifying a region in the
$(x,y)$-plane at $z=0$ of size $[-L_{x}/6,L_{x}/6]\times[0,L_{y}]$ and proceed
by selecting the position of the magnetic lines’ origin evenly across this
region, integrating the lines along the whole box[37]. At $t=0$, being the
equilibrium field lines tangent to $\psi$-isosurfaces, and thus to
$\zeta$-isosurfaces too, the passive tracer $\zeta$ is constant along each
line. In an ideal evolution for the magnetic field, even if field lines are
stretched and distorted, $\zeta$ would remain constant along each line, i.e.,
each line would maintain its original connectivity to the sole blue/red (in
terms of the passive tracer color scale) plasma. This is what we observe
during the early evolution of our simulations. As soon as the non-ideal terms
in Ohm’s law start to increase in importance, allowing for reconnection, a
variation in $\zeta$ along a field line would underline the fact that this
line has changed connectivity. For this reason, we compute the value of the
passive tracer along each line, which enables us to determine the latitude at
which a given line crosses the magnetopause and changes its connection and to
keep track of how many times a line undergoes such a process. In order to
ensure that a line is crossing the frontier between the two plasma regions, we
select a threshold of $\zeta_{\text{jump}}=0.05$, and we define a line as
“reconnected” if it exhibits a variation of $\zeta$ of, at least,
$2\zeta_{\text{jump}}$ across the magnetopause. In this way, if, for example,
a line originates at the bottom of the numerical box at
$\zeta_{\text{in}}<\zeta_{\text{fr}}-\zeta_{\text{jump}}$ (i.e., in the MS
region), then it is counted as reconnected only if the value of the passive
tracer computed along the line itself, as it goes up the box, reaches anywhere
a value $\zeta>\zeta_{\text{fr}}+\zeta_{\text{jump}}$. Hence, we define a
Connection Change (CC) as the crossing of the magnetopause (satisfying the
threshold condition) by a reconnected line. For a given line, the number of
CCs, from the bottom to the top of the box, defines it as a once-, twice-, or
more than twice-reconnected.

Figure 5:  Histograms of the CC counts (for $\zeta_{\text{jump}}=0.05$) for
both Run A (top) and B (bottom) during the nonlinear phase of the instability
at $t=460$, $500$, $580\Omega_{\text{c}}^{-1}$. The grey dotted line in both
panels indicates the $z$-position of the equatorial plane.

Figure 5 presents the cumulative count of CCs across both simulations, as a
function of $z$, for three different times. Panel (A) shows how, in Run A, the
latitudinal distribution of CCs develops two peaks at mid-latitudes in both
hemispheres, exhibiting, therefore, a symmetrical evolution with respect to
the equator $z=L_{z}$. This pattern aligns with the results depicted in Figure
3 panel (A3), where current sheets manifest at mid-latitudes. Consequently,
the development of current sheets at mid-latitudes initiates magnetic
reconnection, resulting in a concentration of CCs for the analyzed field
lines, precisely at the latitudes where the earlier current sheets formed.
This configuration for the current sheets (at mid-latitudes, away from the
vortices), combined with the distribution of reconnection events, suggests
that only the double-MLR process is active in this simulation.

In Figure 5 (B), the profiles of the latitudinal distributions in Run B
develop in time asymmetrically with respect to the equatorial plane. At
$t=460\Omega_{\text{c}}^{-1}$, there is only one main peak, which then
migrates towards the northern hemisphere as time progresses. Already at
$t=500\Omega_{\text{c}}^{-1}$, a second, smaller peak forms in the southern
hemisphere and proceeds to grow and shift towards the mid-latitudes of that
hemisphere ($z\sim 80d_{i}$). In this configuration, the position of the
vortices with respect to the most prominent current sheets suggests that both
MLR and Type-I VIR are active. The main current sheet in the northern
hemisphere (see Figure 3 panel (B3)) is the original current sheet already
present at the equilibrium that has shifted northward because of the
differential advection of the field lines, becoming squeezed by the KH
vortices, thus leading to the corresponding peak in the distribution of
reconnection events (see Figure 5 (B)), compatibly with Type-I VIR events.
Instead, the second, smaller peak in the distribution of reconnection events
is consistent with MLR occurring in the smaller current sheets forming away
from the equator. Interestingly, when comparing the two runs at similar stages
of the KHI evolution, the sum of CCs is of the same order of magnitude. For
example, at $t=460\Omega_{\text{c}}^{-1}$ for Run A and at
$t=500\Omega_{\text{c}}^{-1}$ for Run B, the sum of CCs is of the order of
$8\times 10^{4}$ and $7.7\times 10^{4}$, respectively for Run A and B.

Figure 6:  The temporal evolution represents the count of CC as a percentage
of the $\sim 22000$ integrated magnetic field lines at each time step. Both in
(A) and (B), we show the percentage of magnetic field lines that reconnected
once (black solid line), twice (dashed line), and more than twice (dotted
line). The cumulative percentage is shown in red.

Figure 6 illustrates the temporal progression of CC counts for the reconnected
magnetic field lines, showcasing the topological evolution during the
nonlinear phase of both simulations. In particular, it shows how many lines
have reconnected once (solid black line), twice (dashed black line), or more
than twice (dotted black line) during the nonlinear phase. As a consequence of
the different growth rates of the instability in the two runs, the nonlinear
phase starts at slightly different times for the two simulations: around
$t=380\Omega_{\text{c}}^{-1}$ for Run A and $t=420\Omega_{\text{c}}^{-1}$ for
Run B. As previously mentioned, we can identify three primary stages in the
vortex evolution for both runs: vortex formation, pairing, and disruption.
These stages aid in interpreting the data in Figure 6. Consequently, we define
the start of the pairing phase at $t=480\Omega_{\text{c}}^{-1}$ and
$t=580\Omega_{\text{c}}^{-1}$ for Run A and B, respectively. The disruption
phase, instead, starts around $t=580\Omega_{\text{c}}^{-1}$ and
$t=660\Omega_{\text{c}}^{-1}$ for Run A and B, respectively.

In Figure 6 (A), for a short initial period, the percentage of twice-
reconnected lines in Run A is of the same order as the percentage of once-
reconnected lines. However, as the vortices fully form and right before they
start pairing ($t=460\Omega_{\text{c}}^{-1}$), we see that the percentage of
twice-reconnected field lines becomes dominant. This evolution that sees a
dominance in twice-reconnected lines also continues during the pairing phase
and is compatible with what is expected from a symmetric evolution of current
sheets, which develop at mid-latitudes as shown in Figure 3 (A1) and (A2).
Thus, this information on the twice-reconnected lines further supports the
idea that only double-MLR is active when an initial magnetic rotation is
absent. In Figure 6 (B), instead, the percentage of once-reconnected lines in
Run B is dominant throughout the evolution until the start of the disruption
stage, where at $t\sim 660\Omega_{\text{c}}^{-1}$ we observe a sudden drop and
the percentage of twice- and more than twice-reconnected lines becomes more
important. The dominance of only once-reconnected field lines is also in
accordance with what was earlier discussed on the formation of main current
sheets in the northern hemisphere, with the consequent dominance of Type-I
VIR. As the vortex disruption phase starts ($t=580\Omega_{\text{c}}^{-1}$ in
Run A, and $t=660\Omega_{\text{c}}^{-1}$ in Run B), the number of field lines
that have reconnected more than twice gradually increases, overtakes the
number of once-reconnected lines, and finally approaches the number of twice-
reconnected lines. The increase in reconnection events along each line can be
justified by the start of secondary instabilities, like the development of
secondary KHIs. In Figure 7, we show indeed that for both runs, the large
vortex structures are unstable to secondary KHI, again at different times
because of the different growth rates of the primary instability. Small-scale
KH vortices form along the arms of the primary vortices at latitudes at which
they are most enhanced, specifically the equatorial plane in Run A and the
southern hemisphere in Run B. There, local velocity shears have developed
around the vortex arms and are sufficiently intense to be unstable. Since the
current is not zero, the secondary vortices pinch it and boost the production
of reconnected lines[19]. In particular, in Figure 6 (A) and (B), the dotted
line, which represents the percentage of lines that reconnect many times,
displays a sudden increase around the times at which the secondary KHI
develops (see Figure 7). In Run A, the difference in percentage between the
once- and twice-reconnected lines remains mostly constant at $\sim
20\mbox{--}25\%$ from the start of the pairing phase
($t=480\Omega_{\text{c}}^{-1}$) until the end of the simulation. On the
contrary, the number of more than twice-reconnected lines suddenly increases
as secondary KHI starts. In Run B instead, this difference between once- and
twice-reconnected lines (with, initially, once-reconnected lines dominating
the dynamics) is not only non-constant but is even inverted. We see indeed
that as the pairing stage starts ($t=580\Omega_{\text{c}}^{-1}$), this
difference is about $\sim 20\%$, but it decreases until there is an inversion
and twice-reconnected lines become dominant (and far from negligible) at
$t=680\Omega_{\text{c}}^{-1}$, after the start of the disruption phase of the
vortices. While the symmetric nature of Run A naturally leads to (and
confirms) a higher number of twice-reconnected lines compared to once- and
more than twice-reconnected lines, it is noteworthy that their overall count
is approximately the same as in Run B at the end of the simulations.

Figure 7:  Zoom-in onto a region of size $60\times 60d_{i}^{2}$ in the
$(x,y)$-plane, highlighting secondary KHI developing in both runs. Left: Run A
at time $t=580\Omega_{\text{c}}^{-1}$ in the equatorial plane ($z\simeq
188d_{i}$). Right: Run B at time $t=660\Omega_{\text{c}}^{-1}$ in the southern
hemisphere ($z=105d_{i}$).

Magnetic reconnection is expected to lead to the development of a mixing layer
across the magnetopause, in which the mixing between the two different plasmas
is dominated by the parallel streaming of particles along reconnected
lines[15, 4]. This plasma mixing, naturally present in kinetic simulation, is
not described by our fluid model. Nevertheless, we can estimate the plasma
mixing by looking at the magnetic field line mixing. For this purpose, we
define a quantity that gives us the perpendicular displacement of magnetic
field lines with respect to the position it would have during an ideal
evolution. We recall that $\psi$ at $t=0$ can be seen as a line label and, at
the same time, as a measure of the distance of a field line from the
magnetopause ($\psi=0$ isosurface). If we let $\psi$ evolve ideally (following
an ideal advection equation, as the passive tracer does, $\psi$ for $t>0$
would no longer describe the magnetic field; however, it would still measure
the distance from the perturbed magnetopause folded by the KH vortices. We
thus define the quantity $\Delta\psi=\mathrm{max}(\psi)-\mathrm{min}(\psi)$
along one line, which estimates the perpendicular displacement of the line due
to the non-ideal evolution, i.e., not taking into account the contribution of
the displacement related to the ideal evolution.

In Figure 8, we show the evolution in time of the average squared displacement
$\langle\Delta\psi^{2}\rangle$, finding similar profiles for the two
simulations. This suggests that the presence of an initial magnetic shear does
not significantly influence the movement of magnetic field lines during the
actual system evolution with respect to the positions they would have in an
ideal evolution. However, the profiles are not entirely overlapping, as we
show in Figure 8: Run B reaches a plateau by the end of the simulation, while
Run A keeps increasing and does not appear to be stabilizing around a maximum
value for the average displacement before the end of the simulation time. Even
if the relation between $\langle\Delta\psi^{2}\rangle$ and $t$ is not
perfectly linear, the slope of the curve is of order $0.2$ for Run A and of
order $0.3$ for run B. This suggests an effective magnetic and matter
diffusion coefficient of order $D_{\text{eff}}\simeq qv_{\text{A},z}d_{i}$,
with $q$ the slopes mentioned above, and $v_{\text{A},z}$ the $z$-component of
the Alfvén velocity. Taking $v_{\text{A},z}\simeq 400$ km/s and $d_{i}\simeq
250$ km [19] we obtain $D_{\text{eff}}(\text{Run~{}A})\simeq 2\times 10^{9}$
m2/s and $D_{\text{eff}}(\text{Run~{}B})\simeq 3\times 10^{9}$ m2/s. Both
values are in agreement with the diffusion coefficient needed to account for
the observed plasma mixing $D_{\text{eff}}\sim 10^{9}$ m2/s, derived in a
simplified diffusive model of the magnetopause[38].

Figure 8:  Plot of the average displacement of magnetic field lines
$\langle\Delta\psi^{2}\rangle$ during the nonlinear phase with respect to
their initial position for both simulations (red - Run B and black - Run A).
Since the two simulations evolve at different rates and enter the nonlinear
phase at different times, we shifted the profile of Run B to match that of Run
A and ease the comparison.

##  IV Conclusions

In this paper, we presented our work on 3D resistive Hall-MHD numerical
simulations that reproduce the KHI dynamics at the Earth’s magnetospheric
flanks. Our focus was on the interaction of SW and Earth’s MS, driven by the
development of KHI and induced magnetic reconnection events. We designed two
simulations with similar setups that differ mainly in the orientation of the
magnetic field of the two plasma regions, the SW and the MS. Our 3D framework
provides a more realistic representation of the magnetospheric flank,
capturing three-dimensional processes such as high-latitude stabilization of
KHI and double ML reconnection events. The two different setups enabled us to
study how a rotation in the magnetic field influences the system, allowing for
a quantitative comparison of plasma dynamics and magnetic reconnection events.

We summarize the main results we obtained from the direct comparison between
the two simulations and with previous works in the literature as follows:

  * •

In both simulations, the KHI dominates the plasma dynamics, with the creation
of large-scale vortices at low-latitude (folding phase) that merge (pairing
phase) and finally get disrupted by secondary processes (disruption phase).
The dispersion relation[33, 34] foretells that Run B should take longer to
reach the folding phase of the instability. Indeed, during the early
development of the KHI, we infer from data the growth rates
$\gamma_{\text{B}}=0.0088<\gamma_{\text{A}}=0.0112$. As a consequence, Run B
reaches each phase of the instability at a later time than Run A.

  * •

Focussing on the latitudinal evolution of KH vortices, they initially grow at
low latitudes, but for a configuration with a purely northward magnetic field
(Run A), the system exhibits a nearly symmetrical development of the plasma
dynamics with respect to the equatorial region, until the end of the
simulation. On the contrary, in Run B, vortices quickly migrate towards the
southern hemisphere. This is due to the differential advection of magnetic
field lines, tied at high latitudes to the unperturbed plasma but advected at
low latitudes by the KH dynamics. Differential advection decreases the
magnetic shear angle $\phi_{\text{B}}$ in the southern hemisphere, enhancing
the KHI there, while it increases $\phi_{\text{B}}$ in the northern
hemisphere, weakening the instability’s growth there, with the consequent
migration of vortices. This behavior is in agreement with previous works[36,
19, 20]; in particular, the asymmetrical development of the system observed in
Run B is determined by the sign of
${\bm{\omega}}_{\text{eq}}\cdot{\pmb{\mathit{j}}}_{\text{eq}}$. In our case,
${\bm{\omega}}_{\text{eq}}\cdot{\pmb{\mathit{j}}}_{\text{eq}}>0$, so vortices
develop in the southern hemisphere, as expected. Contrary to the previous
works, because of the symmetry in the rotation of the initial magnetic field
in Run B, we do not introduce a tilt in the vortex axis. This choice allows
for a one-to-one comparison between the two different configurations (with and
without the magnetic shear), having both the vortex axes aligned in the
$z$-direction.

  * •

Our analysis of the current sheet development shows a similar evolution, with
a symmetrical formation of current sheets in Run A and an asymmetrical one in
Run B with respect to the equatorial plane. The phenomenon of differential
advection induces a bending of magnetic field lines, generating current
sheets. In Run A, with zero initial current, current sheets form at mid-
latitudes. Conversely, in Run B, where an initial current was present, the
current is enhanced in the northern hemisphere by both current pinching and
differential advection. Only when the system is well into the nonlinear phase
does a second smaller current sheet develop in the southern hemisphere below
the latitude at which vortices form. Indeed, differential advection initially
lowers $\phi_{\text{B}}$ in the southern hemisphere, reducing the current, but
eventually, it is able to reverse $\phi_{\text{B}}$, creating new current
sheets there.

  * •

The latitudinal distribution of the connection changes (CCs) shows a
qualitatively similar symmetrical/asymmetrical evolution for the two
simulations, with a clear correspondence between the peaks of the latitudinal
distributions of CC count and of the current density magnitude. In Run A, the
configuration of current sheets at mid-latitudes, away from the KH vortices,
combined with the latitudinal distribution of reconnection events, shows that
only double-MLR is active. This was suggested in earlier works[12, 37], but
the CC count was never quantified. In Run B, the vortex positions relative to
the main current sheets suggest simultaneous MLR and Type-I VIR activity. The
main current sheet, which drifted towards the northern hemisphere, causes a
peak in reconnection events associated with Type-I VIR. Meanwhile, a second
smaller peak in reconnection events is consistent with MLR, occurring in less
intense current sheets away from the equator, in the southern hemisphere. The
similar characteristics of the two simulations enabled us to quantitatively
compare the number of reconnection events at similar times during the systems’
evolution. Our analysis shows that the count of all reconnection events along
the computational box for each run is of the same order of magnitude: $\sim
8\times 10^{4}$ at the folding phase, $\sim 2\times 10^{5}$ at the merging
phase, and $\sim 4\times 10^{5}$ when the onset of secondary KHI disrupts the
large-scale vortex structures.

  * •

Because of the design of the two simulations, identical in size and very
similar in initialization, for the first time, we can study and compare
quantitatively the topological evolution of magnetic field lines in these
systems. During the nonlinear phase, we found that, in Run A, there is a clear
dominance of lines that reconnect twice. The number of twice-reconnected lines
is constantly $\sim$ twice that of once-reconnected lines from the start of
the pairing phase until the end of the simulation time for Run A. Instead, in
Run B, once-reconnected lines dominate during the initial evolution of the
system, but their dominance ends when secondary processes disrupt the KH
vortices. It is essential to underline that even if once-reconnected lines are
dominant in Run B, the number of twice-reconnected lines is far from
insignificant. For the first time, we highlight that by the end of both
simulations, the proportion of twice-reconnected lines comprises nearly
identical percentages, approximately $43\%$ and $45\%$, of the total number of
reconnected lines in Runs A and B, respectively. Twice-reconnected lines,
indicating multiple reconnection events and mixing of material, play an
important role in determining not only the plasma transport across the
magnetosphere but also in allowing the specific entropy to increase when
passing from the magnetosheath to the magnetosphere[16, 13].

  * •

In earlier 2D[30, 31, 32] and 3D[5, 19, 20] works, it was observed that, after
the folding phase, the conditions were favorable for the development of
secondary KH and Rayleigh-Taylor instabilities inside the vortices. In all
these cases, the Atwood number of the equilibrium configuration was different
from zero. A density jump between the plasmas on the two sides of the
magnetopause is necessary for the development of secondary Rayleigh-Taylor
instability, for which the vortex rotation provides the effective gravity.
This difference is also needed, at least in 2D, to have strong vorticity
layers along the vortex arms and, thus, secondary KHI. We observe that, even
in the absence of a density jump in the initial setting of both simulations,
we see the formation of secondary KH vortices along the arms of primary KH
vortices. These dynamics suggest that the vortex layers where secondary KH
vortices grow are created there by 3D effects, in particular by the combined
action of magnetic tension and line tying. These secondary processes induce
magnetic reconnection at intermediate latitudes between the two peaks in the
CC count. Thus, they influence the number of times each magnetic field line
reconnects and increase the amount of line reconnected more than twice, which
reaches the total percentage of $\sim 46\%$ and $\sim 36\%$ in Run A and Run
B, respectively.

  * •

Finally, the effective magnetic diffusion coefficient $D_{\text{eff}}$
associated with reconnection is of order $\sim 10^{9}$ m2/s for both
simulations. This value can explain plasma mixing at the magnetopause[38], and
is in agreement with values inferred in earlier simulations[5, 39, 19].

In summary, our comparative analysis of these two simulations stresses the
crucial role of magnetic field orientation in shaping the Earth’s flank
magnetopause. Specifically, it influences vortex migration, current sheet
formation, and the distribution of reconnection events in latitude. These
findings highlight the necessity of considering 3D effects and line tying for
an accurate portrayal of the complex phenomena ruling the flank dynamics.

###### Acknowledgements.

F.B. acknowledges support from the FED-tWIN programme (profile Prf-2020-004,
project “ENERGY”) issued by BELSPO. This work was granted access to the HPC
resources of IDRIS under the allocation 2022-AD010413669 made by GENCI.

## References

  * Dungey [1961] J. W. Dungey, “Interplanetary magnetic field and the auroral zones,” Phys. Rev. Lett. 6, 47–48 (1961).
  * Coppi [1964] B. Coppi, ““inertial” instabilities in plasmas,” Physics Letters 11, 226–228 (1964).
  * Coppi, Laval, and Pellat [1966] B. Coppi, G. Laval, and R. Pellat, “Dynamics of the geomagnetic tail,” Phys. Rev. Lett. 16, 1207–1210 (1966).
  * Daughton  _et al._ [2014] W. Daughton, T. K. M. Nakamura, H. Karimabadi, V. Roytershteyn, and B. Loring, “Computing the reconnection rate in turbulent kinetic layers by using electron mixing to identify topology,” Physics of Plasmas 21, 052307 (2014).
  * Nakamura and Daughton [2014] T. K. M. Nakamura and W. Daughton, “Turbulent plasma transport across the Earth’s low-latitude boundary layer,” Geophysics Research Letters 41, 8704–8712 (2014).
  * Song and Russell [1992] P. Song and C. T. Russell, “Model of the formation of the low-latitude boundary layer for strongly northward interplanetary magnetic field,” Journal of Geophysical Research: Space Physics 97, 1411–1420 (1992), https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/91JA02377 .
  * Mitchell  _et al._ [1987] D. G. Mitchell, F. Kutchko, D. J. Williams, T. E. Eastman, L. A. Frank, and C. T. Russell, “An extended study of the low-latitude boundary layer on the dawn and dusk flanks of the magnetosphere,” Journal of Geophysics Research 92, 7394–7404 (1987).
  * Knoll and Chacón [2002] D. A. Knoll and L. Chacón, “Magnetic reconnection in the two-dimensional kelvin-helmholtz instability,” Phys. Rev. Lett. 88, 215003 (2002).
  * Nakamura, Fujimoto, and Otto [2008] T. K. M. Nakamura, M. Fujimoto, and A. Otto, “Structure of an mhd-scale kelvin-helmholtz vortex: Two-dimensional two-fluid simulations including finite electron inertial effects,” Journal of Geophysical Research: Space Physics 113 (2008), https://doi.org/10.1029/2007JA012803, https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2007JA012803 .
  * Otto and Fairfield [2000] A. Otto and D. H. Fairfield, “Kelvin-helmholtz instability at the magnetotail boundary: Mhd simulation and comparison with geotail observations,” Journal of Geophysical Research: Space Physics 105, 21175–21190 (2000), https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/1999JA000312 .
  * Henri  _et al._ [2012] P. Henri, F. Califano, M. Faganello, and F. Pegoraro, “Magnetised Kelvin-Helmholtz instability in the intermediate regime between subsonic and supersonic regimes,” Physics of Plasmas 19, 072908 (2012), https://pubs.aip.org/aip/pop/article-pdf/doi/10.1063/1.4739234/16115636/072908_1_online.pdf .
  * Faganello  _et al._ [2012a] M. Faganello, F. Califano, F. Pegoraro, and T. Andreussi, “Double mid-latitude dynamical reconnection at the magnetopause: An efficient mechanism allowing solar wind to enter the Earth’s magnetosphere,” EPL (Europhysics Letters) 100, 69001 (2012a).
  * Faganello  _et al._ [2012b] M. Faganello, F. Califano, F. Pegoraro, T. Andreussi, and S. Benkadda, “Magnetic reconnection and Kelvin-Helmholtz instabilities at the Earth’s magnetopause,” Plasma Physics and Controlled Fusion 54, 124037 (2012b).
  * Faganello  _et al._ [2014] M. Faganello, F. Califano, F. Pegoraro, and A. Retinò, “Kelvin-Helmholtz vortices and double mid-latitude reconnection at the Earth’s magnetopause: Comparison between observations and simulations,” EPL (Europhysics Letters) 107, 19001 (2014).
  * Nakamura  _et al._ [2013] T. K. M. Nakamura, W. Daughton, H. Karimabadi, and S. Eriksson, “Three-dimensional dynamics of vortex-induced reconnection and comparison with THEMIS observations,” Journal of Geophysical Research (Space Physics) 118, 5742–5757 (2013).
  * Johnson and Wing [2009] J. R. Johnson and S. Wing, “Northward interplanetary magnetic field plasma sheet entropies,” Journal of Geophysical Research (Space Physics) 114, A00D08 (2009).
  * Eriksson  _et al._ [2016] S. Eriksson, B. Lavraud, F. D. Wilder, J. E. Stawarz, B. L. Giles, J. L. Burch, W. Baumjohann, R. E. Ergun, P. A. Lindqvist, W. Magnes, C. J. Pollock, C. T. Russell, Y. Saito, R. J. Strangeway, R. B. Torbert, D. J. Gershman, Y. V. Khotyaintsev, J. C. Dorelli, S. J. Schwartz, L. Avanov, E. Grimes, Y. Vernisse, A. P. Sturner, T. D. Phan, G. T. Marklund, T. E. Moore, W. R. Paterson, and K. A. Goodrich, “Magnetospheric Multiscale observations of magnetic reconnection associated with Kelvin-Helmholtz waves,” Geophysics Research Letters 43, 5606–5615 (2016).
  * Vernisse  _et al._ [2016] Y. Vernisse, B. Lavraud, S. Eriksson, D. J. Gershman, J. Dorelli, C. Pollock, B. Giles, N. Aunai, L. Avanov, J. Burch, M. Chandler, V. Coffey, J. Dargent, R. E. Ergun, C. J. Farrugia, V. Génot, D. B. Graham, H. Hasegawa, C. Jacquey, I. Kacem, Y. Khotyaintsev, W. Li, W. Magnes, A. Marchaudon, T. Moore, W. Paterson, E. Penou, T. D. Phan, A. Retino, C. T. Russell, Y. Saito, J. A. Sauvaud, R. Torbert, F. D. Wilder, and S. Yokota, “Signatures of complex magnetic topologies from multiple reconnection sites induced by Kelvin-Helmholtz instability,” Journal of Geophysical Research (Space Physics) 121, 9926–9939 (2016).
  * Sisti  _et al._ [2019] M. Sisti, M. Faganello, F. Califano, and B. Lavraud, “Satellite Data-Based 3-D Simulation of Kelvin-Helmholtz Instability and Induced Magnetic Reconnection at the Earth’s Magnetopause,” Geophysical Research Letters 46, 11,597–11,605 (2019).
  * Faganello  _et al._ [2022] M. Faganello, M. Sisti, F. Califano, and B. Lavraud, “Kelvin-Helmholtz instability and induced magnetic reconnection at the Earth’s magnetopause: a 3D simulation based on satellite data,” Plasma Physics and Controlled Fusion 64, 044014 (2022).
  * Nykyri and Otto [2004] K. Nykyri and A. Otto, “Influence of the Hall term on KH instability and reconnection inside KH vortices,” Annales Geophysicae 22, 935–949 (2004).
  * Faganello, Califano, and Pegoraro [2008a] M. Faganello, F. Califano, and F. Pegoraro, “Time Window for Magnetic Reconnection in Plasma Configurations with Velocity Shear,” Physical Review Letters 101, 175003 (2008a).
  * Birn and Hesse [2001] J. Birn and M. Hesse, “Geospace Environment Modeling (GEM) magnetic reconnection challenge: Resistive tearing, anisotropic pressure and hall effects,” Journal of Geophysics Research 106, 3737–3750 (2001).
  * Lele [1992] S. K. Lele, “Compact Finite Difference Schemes with Spectral-like Resolution,” Journal of Computational Physics 103, 16–42 (1992).
  * Thompson [1990] K. W. Thompson, “Time-Dependent Boundary Conditions for Hyperbolic Systems, II,” Journal of Computational Physics 89, 439–461 (1990).
  * Faganello, Califano, and Pegoraro [2009] M. Faganello, F. Califano, and F. Pegoraro, “Being on time in magnetic reconnection,” New Journal of Physics 11, 063008 (2009).
  * Vernisse  _et al._ [2020] Y. Vernisse, B. Lavraud, M. Faganello, S. Fadanelli, M. Sisti, F. Califano, S. Eriksson, D. J. Gershman, J. Dorelli, C. Pollock, B. Giles, L. Avanov, J. Burch, J. Dargent, R. E. Ergun, C. J. Farrugia, V. Génot, H. Hasegawa, C. Jacquey, I. Kacem, R. Kieokaew, M. Kuznetsova, T. Moore, T. Nakamura, W. Paterson, E. Penou, T. D. Phan, C. T. Russell, Y. Saito, J.-A. Sauvaud, and S. Toledo-Redondo, “Latitudinal dependence of the kelvin-helmholtz instability and beta dependence of vortex-induced high-guide field magnetic reconnection,” Journal of Geophysical Research: Space Physics 125, e2019JA027333 (2020), e2019JA027333 2019JA027333, https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2019JA027333 .
  * Takagi  _et al._ [2006] K. Takagi, C. Hashimoto, H. Hasegawa, M. Fujimoto, and R. Tandokoro, “Kelvin-Helmholtz instability in a magnetotail flank-like geometry: Three-dimensional MHD simulations,” Journal of Geophysical Research (Space Physics) 111, A08202 (2006).
  * Miura and Pritchett [1982] A. Miura and P. L. Pritchett, “Nonlocal stability analysis of the MHD Kelvin-Helmholtz instability in a compressible plasma,” Journal of Geophysics Research 87, 7431–7444 (1982).
  * Matsumoto and Hoshino [2004] Y. Matsumoto and M. Hoshino, “Onset of turbulence induced by a Kelvin-Helmholtz vortex,” Geophysics Research Letters 31, L02807 (2004).
  * Faganello, Califano, and Pegoraro [2008b] M. Faganello, F. Califano, and F. Pegoraro, “Competing Mechanisms of Plasma Transport in Inhomogeneous Configurations with Velocity Shear: The Solar-Wind Interaction with Earth’s Magnetosphere,” Physical Review Letters 100, 015001 (2008b).
  * Tenerani  _et al._ [2011] A. Tenerani, M. Faganello, F. Califano, and F. Pegoraro, “Nonlinear vortex dynamics in an inhomogeneous magnetized plasma with a sheared velocity field,” Plasma Physics and Controlled Fusion 53, 015003 (2011).
  * Axford [1960] W. I. Axford, “The stability of plane current-vortex sheets,” The Quarterly Journal of Mechanics and Applied Mathematics 13, 314–324 (1960), https://academic.oup.com/qjmam/article-pdf/13/3/314/5440261/13-3-314.pdf .
  * Sen [1963] A. K. Sen, “Stability of Hydromagnetic Kelvin-Helmholtz Discontinuity,” Physics of Fluids 6, 1154–1163 (1963).
  * Faganello and Califano [2017] M. Faganello and F. Califano, “Magnetized kelvin–helmholtz instability: theory and simulations in the earth’s magnetosphere context,” Journal of Plasma Physics 83, 535830601 (2017).
  * Fadanelli  _et al._ [2018] S. Fadanelli, M. Faganello, F. Califano, S. S. Cerri, F. Pegoraro, and B. Lavraud, “North-South Asymmetric Kelvin-Helmholtz Instability and Induced Reconnection at the Earth’s Magnetospheric Flanks,” Journal of Geophysical Research (Space Physics) 123, 9340–9356 (2018), arXiv:1805.01466 [physics.space-ph] .
  * Borgogno  _et al._ [2015] D. Borgogno, F. Califano, M. Faganello, and F. Pegoraro, “Double-reconnected magnetic structures driven by Kelvin-Helmholtz vortices at the Earth’s magnetosphere,” Physics of Plasmas 22, 032301 (2015).
  * Sonnerup [1980] B. U. O. Sonnerup, “Theory of the low-latitude boundary layer,” Journal of Geophysics Research 85, 2017–2026 (1980).
  * Nakamura  _et al._ [2017] T. K. M. Nakamura, H. Hasegawa, W. Daughton, S. Eriksson, W. Y. Li, and R. Nakamura, “Turbulent mass transfer caused by vortex induced reconnection in collisionless magnetospheric plasmas,” Nature Communications 8, 1582 (2017).


# Tuning the spontaneous exchange bias effect with Ba to Sr partial
substitution in La1.5(Sr0.5-xBax)CoMnO6

M. Boldrin  Instituto de Física, Universidade Federal de Goiás, 74001-970,
Goiânia, GO, Brazil  A. G. Silva  Instituto de Física, Universidade Federal de
Goiás, 74001-970, Goiânia, GO, Brazil  L. T. Coutrim  Instituto de Física,
Universidade Federal de Goiás, 74001-970, Goiânia, GO, Brazil  J. R. Jesus
Centro Brasileiro de Pesquisas Físicas, 22290-180, Rio de Janeiro, RJ, Brazil
C. Macchiutti  Centro Brasileiro de Pesquisas Físicas, 22290-180, Rio de
Janeiro, RJ, Brazil  E. M. Bittar  Centro Brasileiro de Pesquisas Físicas,
22290-180, Rio de Janeiro, RJ, Brazil  L. Bufaiçal<EMAIL_ADDRESS>Instituto
de Física, Universidade Federal de Goiás, 74001-970, Goiânia, GO, Brazil

###### Abstract

The spontaneous exchange bias (SEB) effect is a remarkable phenomenon recently
observed in some reentrant spin-glass materials. Here we investigate the SEB
in La1.5(Sr0.5-xBax)CoMnO6 double-perovskites, a system with multifarious
magnetic phases for which a notable increase in the exchange bias field is
observed for intermediate Sr/Ba concentrations. The Ba to Sr substitution
leads to the enhancement of the crystal lattice, which is accompanied by the
raise of both the effective magnetic moment ($\mu_{eff}$) and the
antiferromagnetic (AFM) transition temperature that is observed below the
ferromagnetic ordering. Such increases are likely related to the increased
fraction of Co3+ in the high spin configuration, leading to the enhancement of
Co3+–O–Mn4+ AFM phase and to the reduction in the uncompensation of the AFM
coupling between Co and Mn. The combined effect of the increased $\mu_{eff}$
and AFM phase plausible explains the changes in the SEB effect.

Conventionally, the exchange bias (EB) effect is manifested as a shift in
hysteresis loop measurements [$M(H)$] carried after cooling the system in the
presence of an external magnetic field ($H$) Nogues . Due to its potential
application in ultra-dense magnetic recording heads, magnetoresistive devices
and spin valves, this phenomena has been extensively investigated in the last
decades. For many years it was thought that the cooling $H$ was a necessary
condition to set the unidirectional anisotropy (UA) by breaking the symmetry
at the magnetic phase interfaces. Recently, however, it was observed that some
materials exhibit such asymmetry even after being cooled in zero $H$ Wang ;
Nayak ; Maity . These are the so-called spontaneous EB (SEB) materials, for
which the EB effect observed at low $T$ seems to be related to the presence of
spin glass (SG)-like phases PRB2018 ; New_SEB .

There is an effort of the academic community to find new materials exhibiting
robust SEB at higher $T$. In this sense, the synthesis and investigation of
new materials exhibiting such effect is mandatory. Double-perovskite (DP)
compounds are prospective candidates since their intrinsic structural and
magnetic inhomogeneities usually leads to magnetic phase competition and
frustration, key ingredients to achieve glassy magnetic behavior Sami ; Mydosh
. Among the DP materials, La1.5Sr0.5CoMnO6 stands out by presenting a giant
SEB effect Murthy . The Ba substitution also leads to a very large SEB in
La1.5Ba0.5CoMnO6 Ba2020 , whilst for La1.5Ca0.5CoMnO6 the effect is very
subtle Ca2017 . A detailed investigation of these materials’ magnetic and
electronic properties by means of X-ray absorption spectroscopy (XAS) and
X-ray magnetic circular dichroism (XMCD) has shown that the UA in these
compounds results from uncompensated antiferromagnetic (AFM) coupling between
Co and Mn ions. The distinct magnitude of the SEB effect observed in each
compound is related to the proportion of Co3+ high spin (HS) and non-magnetic
low spin (LS) configuration PRB2019 .

The SEB effect in La1.5A0.5CoMnO6 (A = Ba, Sr, Ca) is directly related to the
Co3+ spin state, which in turn results from a delicate balance between the
crystal field and the interatomic exchange energies, being thus very sensitive
to structural changes Raveau . In this work we investigate the changes in the
SEB when Sr is partially replaced by Ba in the La1.5(Sr0.5-xBax)CoMnO6 series.
A detailed investigation of the structural, electronic, and magnetic
properties of the system was carried by means of X-ray powder diffraction
(XRD) and magnetometry. Our results show a pronounced enhancement of the SEB
effect for intermediate Sr/Ba concentrations. This can be causally linked to
the changes in the balance between the crystal field and the intra-orbital
Coulomb repulsion acting on the Co ions, affecting its spin state and
consequently its coupling with the Mn ions.

Polycrystalline samples of La1.5(Sr0.5-xBax)CoMnO6 ($x$ = 0, 0.2, 0.25, 0.3,
0.4 and 0.5) were prepared by conventional solid state reaction, as described
in the Supplementary Material (SM). Attempts to produce 0 $<x<$ 0.2 samples by
using the same synthesis route were not successful. Laboratory XRD data were
collected using a Bruker $D$8 $Discover$ diffractometer operating with Cu Kα
radiation, while high resolution synchrotron XRD (SXRD) patterns were recorded
at the XPD beamline of the Brazilian Synchrotron Light Laboratory (LNLS),
using a reflection geometry. The data were obtained at room $T$ by one-
dimensional Mythen-1K detector (Dectris), using $\lambda$ = 1.9074 Å. Rietveld
refinements were performed using the GSAS+EXPGUI suiteGSAS . Magnetic
measurements were carried out in both zero field cooled (ZFC) and field cooled
(FC) modes, using a Quantum Design PPMS-VSM magnetometer.

The XRD patterns revealed that all samples belong to the rhombohedral
$R\bar{3}c$ space group, as expected since both $x$ = 0 and 0.5 end members
present the same crystal symmetry Murthy ; Ba2020 . Fig. 1 shows a magnified
view of the XRD data near the ($10\bar{2}$) Bragg reflection, evidencing a
systematic shift of the peak toward smaller angles as the Ba concentration
increases. This is a manifestation of the lattice expansion, reflecting the
fact that the Ba2+ ionic radius is larger than Sr2+ Shannon . Such lattice
expansion with increasing Ba content is also manifested in the increase of the
unit cell volume, obtained from the Rietveld refinements of both XRD and SXRD
data (inset of Fig. 1). Main results obtained from the XRD Rietveld
refinements are displayed in Table 1. Results from SXRD patterns are in Table
S1 of SM. The crystal symmetry of DP compounds is directly related to the
hybridization between the $d$ orbitals, via oxygen $p$ orbitals, being thus
strongly correlated to their electronic and magnetic properties Sami ;
Woodward . Accordingly, the structural changes here observed will remarkably
affect the materials magnetic properties, in special those associated to the
Co spin state.

Figure 1: Magnified view of the ($10\bar{2}$) Bragg reflection in the normalized XRD patterns of La1.5(Sr0.5-xBax)CoMnO6. The inset shows the unit cell volume as a function of Ba concentration ($x$), obtained from XRD (squares) and SXRD (circles). Table 1: Main results obtained from the Rietveld refinements of the XRD data for the La1.5(Sr0.5-xBax)CoMnO6 series. Sample ($x$) |  $a$ (Å) |  $c$ (Å) |  $V$ (Å3) |  $R_{wp}$ (%) |  $R_{p}$ (%)
---|---|---|---|---|---
0 | 5.4733(1) | 13.2575(2) | 343.95(1) | 10.3 | 7.3
0.2 | 5.4852(1) | 13.2937(1) | 346.38(1) | 10.8 | 8.0
0.25 | 5.4886(1) | 13.3157(1) | 347.39(1) | 8.8 | 6.4
0.3 | 5.4885(1) | 13.3221(1) | 347.54(1) | 10.1 | 7.6
0.4 | 5.4952(1) | 13.3531(1) | 349.20(1) | 11.1 | 8.2
0.5 | 5.5025(1) | 13.3777(1) | 350.78(1) | 12.1 | 8.7

Figure 2(a) shows the ZFC and FC magnetization as a function of $T$ [$M(T)$]
curves of $x$ = 0.2 compound, chosen here as a representative sample for the
series. The FC curve shows a ferromagnetic (FM)-like character, while the ZFC
one highlights three anomalies associated to distinct magnetic couplings. The
same overall behavior was observed for all samples (see Fig. S2 in SM).
Previous reports of magnetization, muon spin rotation and relaxation, XAS and
XMCD measurements for the $x$ = 0 and 0.5 compounds have shown a mixed valence
state for Co (Co2+/Co3+), while for Mn it was found a nearly tetravalent
state, although some small amount of Mn3+ was detected Murthy ; Ba2020 ;
PRB2019 . In these studies the first transition temperature, $T_{C1}$, was
ascribed to the Co2+–O–Mn4+ FM coupling, the second one, $T_{C2}$, to the
Co3+–O–Mn3+ FM coupling, while the third anomaly observed at lower $T$ is most
likely related to the Co3+–O–Mn4+ AFM coupling. It is natural to consider the
same scenario for the other Ba/Sr concentrations here investigated. The
transition temperatures obtained from the first derivative of the ZFC curves
(see Figs. S4 and S5 in SM) are displayed in Table 2, where it is observed a
systematic increase of $T_{N}$ with the Ba content. Previous works have also
found SG-like behavior for both $x$ = 0 and 0.5 end members at low
temperatures Murthy ; PRB2019 . Since all samples here investigated present
the same overall magnetic behavior, SG-like features are also expected for the
intermediate Ba/Sr compounds.

Figure 2: (a) ZFC and FC $M(T)$ curves measured at $H$ = 100 Oe for $x$ = 0.2.
(b) $\theta_{CW}$ (square) and $\mu_{eff}$ (circle) as a function of $x$ for
La1.5(Sr0.5-xBax)CoMnO6.

From the fit of the magnetic susceptibility curves at the paramagnetic region,
with the Curie-Weiss (CW) law (see Fig. S3 in SM), we obtain the CW
temperature ($\theta_{CW}$) and the effective magnetic moment ($\mu_{eff}$).
As it can be seen in Table 2, $\theta_{CW}$ is positive for all investigated
samples, giving further evidence of the prominent FM coupling. Its
systematically decrease with increasing the Ba content ($x$) is in accordance
with the increase of the AFM transition temperature ($T_{N}$).

Table 2: Main results obtained from the $M(T)$ and $M(H)$ curves for the La1.5(Sr0.5-xBax)CoMnO6 series. Sample ($x$) |  $T_{C1}$ (K) |  $T_{C2}$ (K) |  $T_{N}$ (K) |  $\theta_{CW}$ (K) |  $\mu_{eff}$ ($\mu_{B}$/f.u.) | C ($K\cdot\mu_{B}$/f.u.) |  $H_{EB}$ (Oe) |  $H_{C}$ (Oe) |  $M_{S}$ ($\mu_{B}$/f.u.)
---|---|---|---|---|---|---|---|---|---
0 | 185 | 178 | 77 | 210 | 6.3 | 1323 | 5112 | 10968 | 2.01
0.2 | 175 | 165 | 84 | 197 | 6.8 | 1340 | 7119 | 6128 | 2.44
0.25 | 181 | 173 | 88 | 193 | 7.1 | 1370 | 9343 | 4629 | 2.26
0.3 | 182 | 171 | 90 | 186 | 7.2 | 1339 | 8528 | 4699 | 2.11
0.4 | 187 | 176 | 106 | 177 | 7.4 | 1310 | 7657 | 4022 | 1.96
0.5 | 181 | 165 | 110 | 165 | 7.9 | 1304 | 4261 | 5792 | 2.41

The $\mu_{eff}$ systematically increases with $x$, plausibly explaining the
increase in $T_{N}$. As aforementioned, previous studies of Co- and Mn-based
DPs have revealed the presence of Co2+/Co3+ mixed valence, with the amounts of
LS/HS Co3+ playing an important role on the SEB effect PRB2019 . Such
investigations have also shown that changes in the A site of the perovskite
structure act mainly in the Co electronic configuration, while the Mn formal
valence remains nearly close to 4+. Thus, the enhancement of $\mu_{eff}$ is
most likely related to the increase in the amount of Co3+ in high HS
configuration present in the system. This is endorsed by the decrease of
$\theta_{CW}$ accompanied by the increase of $T_{N}$, which are likely related
to the (HS)Co3+–O–Mn4+ AFM coupling PRB2019 . These changes will cause a
direct impact in the SEB effect of the La1.5(Sr0.5-xBax)CoMnO6 series as
detailed below.

Fig. 3(a) shows the ZFC $M(H)$ curve of the $x$ = 0.2 sample, measured at 5 K.
It is a closed loop and its shape results from a combined contribution of AFM
and FM phases. The shift towards the left along the $H$ axis characterizes the
SEB effect. Inset shows that the virgin magnetization curve falls outside the
hysteresis loop up to a certain field. This is a common feature of several SEB
DP compounds, being related to field induced spin transition Murthy ; Ba2020 ;
Ca2017 ; Xie . The same overall behavior was observed for the $M(H)$ curves of
all investigated samples, as shown in Fig. S6 of SM. The shift along the $H$
axis is a measure of the EB field, herein defined as $H_{EB}=|H_{1}+H_{2}|/2$,
where $H_{1}$ and $H_{2}$ are respectively the left and right coercive fields.
Fig. 3(b) shows a pronounced increase of the SEB for intermediate contents of
Ba and Sr, being maximum for the $x$ = 0.25 sample.

Figure 3: (a) ZFC $M(H)$ loop for the $x$ = 0.2 sample, measured at 5 K. Inset
shows a magnified view of the initial magnetization curve (virgin curve). (b)
$H_{EB}$ (square) and $C=\mu_{eff}\cdot\theta_{CW}$ (circle) as a function of
$x$ (see text).

The evolution of the SEB effect in La1.5(Sr0.5-xBax)CoMnO6 can be understood
within the framework of strong correlation between the structural and
electronic properties observed for this system. The partial replacement of Sr
by Ba acts to increase the crystal lattice, leading to the decrease of the
crystal field splitting energy and thus favoring the HS configuration. This is
manifested in the increase of $\mu_{eff}$ and in the decrease of
$\theta_{CW}$, which are related to the increased strength of the Co3+–O–Mn4+
AFM coupling. On the other hand, such enhancement in the Co magnetic moment
reduces the uncompensation of its coupling with Mn, which would tend to reduce
the EB field PRB2019 .

The lattice expansion is closely related to the increase in the portion of HS
Co3+, which in turn brings two main outcomes: (i) the increase in the strength
of the AFM coupling and (ii) the decrease in the uncompensated coupling
between the transition-metal (TM) ions. These opposite effects are causally
related to the changes in the SEB effect. In order to quantify this argument,
we define here a coupling parameter $C=\mu_{eff}\cdot\theta_{CW}$. Fig. 3(b)
shows the evolution of $C$ as a function of $x$. As can be noticed, it
presents the same behavior of the $H_{EB}$ curve. One can therefore conclude
that the increase in the portion of HS Co3+ with $x$, manifested in the
enhancement of $\mu_{eff}$, leads to the rise of the AFM phase, which is
necessary for the observation of the SEB effect. In addition, it also leads to
the a decrease in the uncompensation of the coupling between the Co and Mn
superstructures, which acts to decrease the SEB. The intermediate
concentration of Sr/Ba is the region of the series where this joint effect is
maxima, resulting in the pronounced increase of $H_{EB}$.

J. Krishna Murthy et al. have found a maxima in the SEB of La2-xSrxCoMnO6
(0$\leq$x$\leq$1) for intermediate $x$ concentrations Murthy2 . Such effect
was explained in terms of the antisite disorder (ASD) at the Co/Mn site, which
destroys the long-range magnetic order. The UA is set by a field-induced
metamagnetic phase transformation from canted AFM to FM during the initial
magnetization process of the magnetic hysteresis cycle, and the changes in SEB
of La2-xSrxCoMnO6 are causally linked to the ASD via the “saturated”
magnetization in the $M(H)$ curves (the magnetization at $H$ = 70 kOe,
$M_{S}$). In the case of La1.5(Sr0.5-xBax)CoMnO6, although it is not observed
a direct relation between $H_{EB}$ and the $M_{S}$ values [magnetization at
$H$ = 90 kOe in the $M(H)$ curves, see Table 2], the scenario of phase
separation is also plausible. The increment of Ba leads to the increase of the
AFM phase due to the increased portion of HS Co3+, but it also affects the FM
coupling by increasing the distance between the TM ions. Since the EB is an
interface effect, it will depend on the balance between the magnetic phases,
being maxima in the intermediate Sr/Ba region of the series.

In summary, we thoroughly investigated the structural, electronic, and
magnetic properties of the La1.5(Sr0.5-xBax)CoMnO6 system by means of XRD and
magnetometry measurements. The lattice expansion with increasing $x$ is
closely related to the increase of $\mu_{eff}$ and $T_{N}$. The ZFC $M(H)$
loops showed a pronounced increase of the giant SEB for intermediate
concentrations of Sr/Ba, where the combined effects of the enhanced AFM phase
and the decreased uncompensation in the magnetic coupling between Co and Mn,
both resulting from the increased portion of HS Co3+, result in a maximum SEB
effect for the $x$ = 0.25 compound.

## Supplementary Material

See supplementary material for details of the crystals growth and some
structural and magnetic properties of the investigated samples.

###### Acknowledgements.

This work was supported by Conselho Nacional de Desenvolvimento Científico e
Tecnológico (CNPq) (Grant No. 400134/2016-0), Fundação Carlos Chagas Filho de
Amparo à Pesquisa do Estado do Rio de Janeiro (FAPERJ) (Grant No.
E-26/202.798/2019), Fundação de Amparo à Pesquisa do Estado de Goiás (FAPEG)
and Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES).

## Data Availability Statement

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## References

  * (1) J. Nogués and I. K. Schuller, J. Magn. Magn. Mater. 192 203 (1999).
  * (2) B. M. Wang, Y. Liu, P. Ren, B. Xia, K. B. Ruan2 J. B. Yi, J. Ding, X. G. Li, and L. Wang, Phys. Rev. Lett. 106, 0077203 (2011).
  * (3) T. Maity, S. Goswami, D. Bhattacharya, and S. Roy, Phys. Rev. Lett. 110, 107201 (2013).
  * (4) A. K. Nayak, M. Nicklas, S. Chadov, C. Shekhar, Y. Skourski, J. Winterlik, and C. Felser, Phys. Rev. Lett. 110, 127204 (2013).
  * (5) L. T. Coutrim, E. M. Bittar, F. Garcia, and L. Bufaiçal, Phys. Rev. B 98, 064426 (2018).
  * (6) L. Bufaiçal, L. T. Coutrim, E. M. Bittar, and F. Garcia, J. Magn. Magn. Mater. 512 167048 (2020).
  * (7) S. Vasala and M. Karppinen, Prog. Solid State Chem. 43, 1 (2015).
  * (8) J. A. Mydosh, Spin Glasses: An Experimental Introduction (Taylor & Francis, London, 1993).
  * (9) J. Krishna Murthy and A. Venimadhav, Appl. Phys. Lett. 103, 25410 (2013).
  * (10) M. Boldrin, L. T. Coutrim and L. Bufaiçal, Braz. J. Phys. 50, 711 (2020).
  * (11) L. Bufaiçal, R. Finkler, L. T. Coutrim, P. G. Pagliuso, C. Grossi, F. Stavale, E. Baggio-Saitovitch, and E. M. Bittar, J. Magn. Magn. Mater. 433, 271 (2017).
  * (12) L. T. Coutrim, D. Rigitano, C. Macchiutti, T. J. A. Mori, R. Lora-Serrano, E. Granado, E. Sadrollahi, F. J. Litterst, M. B. Fontes, E. Baggio-Saitovitch, E. M. Bittar, and L. Bufaiçal, Phys. Rev. B 100, 054428 (2019).
  * (13) B. Raveau and Md. Motin Seikh, Cobalt Oxides: From Crystal Chemistry to Physics (Wiley-VCH, Weinheim, 2012).
  * (14) A. C. Larson and R. B. Von Dreele, Los Alamos National Laboratory Report No. LAUR 86-748, 2000; B. H. Toby, J. Appl. Crystallogr. 34, 210 (2001).
  * (15) R. D. Shannon, Acta Crystallographica A32 (1976) 751.
  * (16) C. J. Howard, B. J. Kennedy and P. M. Woodward, Acta Crystalogr. B 59, 463 (2003).
  * (17) L. Xie and H. G. Zhang, Curr. Appl. Phys. 18 261 (2018).
  * (18) J. Krishna Murthy, K. D. Chandrasekhar, H. C. Wu, H. D. Yang, J. Y. Lin, and A. Venimadhav, J. Phys.: Condens. Matter 28, 086003 (2016).