# Self-trapping of active particles induced by non-reciprocal interactions in
disordered media

Rodrigo Saavedra<EMAIL_ADDRESS>Departamento de Ciencias de la
Computación, Centro de Investigación Científica y de Educación Superior de
Ensenada, Baja California, México  Fernando Peruani<EMAIL_ADDRESS>Laboratoire de Physique Theéorique et Modélisation, CNRS UMR 8089, CY Cergy
Paris Université, Cergy-Pontoise Cedex, France

###### Abstract

Disordered media is known to modify the transport properties of active matter
depending on interactions between particles and with the substrate. Here we
study systems of active particles with visual-like perception that co-align
with the positions of perceived conspecifics, and anti-align with positions of
static obstacles. We report a novel self-trapping mechanism of particles
forming closed loops that progressively shrink, surrounding one or multiple
obstacles. This mechanism corresponds to a pinning behavior preventing
particle diffusion. Increased co-alignment strength is found to reduce loop
shrinking time, although this effect reaches a plateau at higher strengths.
Loops are found to initially exhibit local polar order, but eventually they
transition to nematic states as they absorb more particles. We show a phase
diagram demonstrating self-trapping occurs within a specific range of aperture
angles of the vision cone.

Most examples of natural systems, if not all, where collective motion occurs
in the wild, take place in heterogeneous media. Examples can be found at all
scales. Microtubules driven by molecular motors form complex patterns inside
the cell where the space is filled by organelles and vesicles [1, 2].
Bacterial colonies exhibit complex collective behaviors, e.g. swarming in
heterogeneous environments such as the soil [3, 4] or highly complex tissues
such as in the gastrointestinal tract [5, 6, 7]. At a larger scale, herds of
mammals migrate long distances traversing rivers, forests, etc. [8, 9]. In
active matter, the influence of static obstacles on the collective behavior
has been a focus of research in recent years. For example, when velocity
alignment is considered, it has been found that a single small obstacle can
have a dramatic impact in the large-scale dynamics of polar flocks, leading to
flow reversals and chaotic dynamics [10]. In the presence of several randomly
placed obstacles, the system becomes disorganized overall, with some of the
particles aggregating into small independent clusters [11, 12, 12]. For
coupled arrays of active particles, e.g. active polymers, it has been found
that they can either rectify or segregate when navigating in porous media [13,
14, 15]. Non-reciprocal interactions pertains to a new class of active matter
that can be found in system interacting via chemotaxis, hydrodynamic forces,
or quorum-sensing [16]. It has been recently explored in theoretical studies
[17, 18, 19, 20, 21], as well as in particle-based numerical simulations, e.g.
in models of particles that move according to positional-based interactions
while considering conspecifics located within a restricted cone of vision [22,
23, 24, 25, 26, 27]. This mechanism, fundamentally different from the
celebrated Vicsek rule, has been successfully applied to describe the
collective dynamics of sheep [28, 29], and to induce cohesive group formation
in swarms of light-activated colloids [30, 31, 32].

In this work, we consider the impact that a heterogeneous medium, i.e.
quenched disorder [33, 34, 35, 36], may have on collective dynamics of active
particles with non-reciprocal interactions. We investigate a new mechanism of
self-trapping of active particles forming a loop around static obstacles. This
behavior is found to emerge from interactions of particles with visual-like
perception combined with positional-based co-alignment between moving
particles, as well as anti-alignment between particles and static obstacles.
We find trapped loops consist of active particles that remain tangentially
oriented in the loop with small fluctuations in the orientation given by
thermal diffusion. Once formed, loops can absorb other incoming particles in
the system, and also some of the particles can escape due to thermal
diffusion. As a result of the increasing number of particles in the loop, the
local order of the structure transitions from polar to nematic. We
characterize the trapped loops in terms of dynamic parameters like the radius
of gyration and mean-squared displacement, as well as distributions for
cluster size and displacement. Our analysis shows that self-trapping
corresponds to pinning behavior preventing particle diffusion. It is found to
induce an exponential decay in the cluster size distribution for small
clusters, with a delta peak at larger sizes, and also a quadratic dependence
in the displacement probability-density function for small displacements with
a range corresponding to the loop size. A phase diagram is also obtained for
varying values of the aperture angles of the vision cone, showing that
trapping occurs only for intermediate values of the parameters. We also study
the emergence of nematic bands that percolate while avoiding the obstacle
configuration, or they can self-trap forming a closed loop surrounding
multiple obstacles. Our study sheds light to the understanding of active
matter with non-reciprocal interactions, including the effect of crowded
environments on the collective behavior and transport properties.

_Model_.–We consider a system of $N$ motile particles and $N_{o}$ static
obstacles distributed inside of a square box of length $L$. The motion of each
motile particle $i$ is governed by

$\displaystyle\dot{\bm{r}}_{i}$ $\displaystyle=v_{0}\bm{e}_{i},$ (1a)
$\displaystyle\dot{\theta}_{i}$
$\displaystyle=\tau_{i}^{\text{att}}(\\{\bm{r}_{j}\\})+\tau_{i}^{\text{rep}}(\\{\bm{r}_{o}\\})+\xi_{i}.$
(1b)

Here $v_{0}$ and $\bm{e}_{i}=(\cos\theta,\sin\theta)^{T}$ are the self-
propulsion magnitude and direction, $\xi_{i}$ is a noise term sampled from a
normal distribution of width $\sqrt{2D_{\theta}}$, where $D_{\theta}$ is the
rotational diffusion coefficient; the terms
$\tau_{i}^{\text{att}}(\\{\bm{r}_{j}\\})$,
$\tau_{i}^{\text{rep}}(\\{\bm{r}_{o}\\})$ correspond to alignment torques
induced by neighbour particles $j$ and obstacles $o$, respectively. They are
given by

$\displaystyle\tau_{i}^{\text{att}}(\\{\bm{r}_{j}\\})$
$\displaystyle=-\frac{\gamma^{\text{att}}}{n^{\text{att}}}\sum_{j\in{V}_{i}^{\text{att}}}\sin(\beta_{ij}-\theta_{i}),$
(2a) $\displaystyle\tau_{i}^{\text{rep}}(\\{\bm{r}_{j}\\})$
$\displaystyle=-\frac{\gamma^{\text{rep}}}{n^{\text{rep}}}\sum_{o\in{V}_{i}^{\text{rep}}}\sin(\beta_{io}-\theta_{i}),$
(2b)

with $\beta_{ij}$ the polar angle of the vector
$\bm{r}_{ij}=\bm{r}_{j}-\bm{r}_{i}$, and similarly for $\beta_{io}$.
Coefficients $\gamma^{\text{att}}$, $\gamma^{\text{rep}}$ are attraction and
repulsion alignment strengths, respectively; $V_{i}^{\text{att}}$,
$V_{i}^{\text{rep}}$ are interaction zones satisfying
${V_{i}}^{\text{att}}\in{V_{i}}^{\text{rep}}$. The interaction zone
$V_{i}^{\mathrm{att}}$ corresponds to a cone of vision of aperture angle
$\alpha$, and $V_{i}^{\mathrm{rep}}$ to a cone of vision of angle $\beta$, see
Fig. 1a.

Figure 1: (a) Illustration of alignment interactions given by Eqs. (2). Each
active particle perceives neighbors within a narrow cone of vision of angle
$\alpha$, and obstacles within a wide cone of vision of aperture $\beta$. (b)
Representative snapshot for a configuration of $N=900$ particles and
$N_{o}=30$ static obstacles. Simulation parameters are $\alpha=0.3$,
$\beta=0.6$, $\sqrt{2D_{\theta}}=0.3$.

Each term is normalized by $n^{\text{att}}$, $n^{\text{rep}}$ the number of
neighbours within each interaction region.

We consider a configuration of homogeneously distributed obstacles, separated
by a distance $d$. As we increase the number of obstacles $N_{o}$, $d$
decreases as well. The effective obstacle size is given by the interaction
with the active particles, which have an interaction range cutoff
$r_{\mathrm{cut}}$. Therefore, obstacles overlap occur only when
$d<2r_{\mathrm{cut}}$. We perform particle-based numerical simulations by
solving Eq. 1 for $N=900$ particles randomly distributed within a square box
of side length $L=30$, with fixed self-propulsion speed $v_{0}=1$, alignment
strengths $\gamma^{\text{att}}=-\gamma$, $\gamma^{\text{rep}}=\gamma$ with
$\gamma=5$, and cone of vision range of $r_{\mathrm{cut}}=1.5$. We choose two
different sets of parameters for $(\alpha,D_{\theta})$. One set with
$\alpha/\pi=0.3$ and $\sqrt{2D_{\theta}}=0.3$ where particles are known to
aggregate into _polar filaments_. Another set with $\alpha=0.6\pi$ and
$\sqrt{2D_{\theta}}=0.8$ where particles form percolating _nematic bands_ ,
see Ref. [23].

_Results_.–A polar filament consists of a chain of particles moving altogether
in the same direction. The direction of movement of the chain is given by an
incidental leader in the head of the structure. The leader is then followed by
a tail of particles. A particle $j$ in the tail reorients its direction of
motion $\bm{e}_{j}$ according to neighbors it perceives within narrow cone of
vision $V^{att}_{j}$ of aperture angle $\alpha$. The leader’s self-propulsion
orientation $\bm{e}_{i}$ is determined by rotational diffusion $D_{\theta}$.
Active particles can also scatter away from an obstacle that is perceived
within the wide cone of vision $V^{rep}_{i}$ of aperture angle $\beta$, see
Fig. 1b. Particles can escape a polar filament due to a strong fluctuation of
their self-propulsion orientation, e.g. for high values of $D_{\theta}$, on
encounter with an obstacle, or on encounter with more particles. Collisions of
several filaments can also lead to aggregation into a single bulky structure,
which does not necessarily consist of a single chain of head-to-tail
particles, but can contain many particles in a localized region.

For $N_{o}=30$ obstacles, we observe that filaments navigate while trying to
avoid obstacles, and furthermore they show to self-trap around an obstacle,
see Fig. 1b. Loop formation can occur when the incidental leader of a motile
filament starts following the tail. In the presence of obstacles, a loop can
be formed from a single polar filament that reorients due to collisions with
obstacles, thus leading to a closed polar structure. However, other particles
can also join the loop in any direction, thus leading to a final closed
nematic structure, where some of the particles in the structure rotate
clockwise with respect to the center of the group, and some others rotate
counter-clockwise. Once formed, the loop shrinks down arguably due to an
effective line tension, then when an obstacle is in the center of the
configuration the loop displays self-trapping: the obstacle prevents the loop
from shrinking down any further, meanwhile particles still align with their
neighbors keeping the loop cohesive, see Fig. 2b.

Figure 2: Loop shrinking. (a) Without an obstacle, the loop shrinks down to a
size $r<r_{o}$, then fragments into several polar filaments that propagate
radially out of the center. (b) With an obstacle the loop shrinks down to a
size $r=r_{o}$ and stabilizes around the obstacle. (c,d) Time evolution of
dynamic parameters corresponding to configurations in (a,b), namely radius of
gyration $R_{g}(t)$, and mean squared displacement $\sigma^{2}(t)$. (e)
Probability density of finding a particle with radial orientation $e_{r}$ with
respect to the center of the configuration at $t_{1}$.

In contrast, when no obstacles are present, the shrinking loop becomes
unstable when it reaches a minimum size, then fragments into several polar
filaments propagating radially outwards of the loop center, see Fig. 2a. To
quantify the shrinking behavior, we compute the radius of gyration $R_{g}$
with respect to the center of a loop formed by $N=100$. Results in Fig. 2c
show that $R_{g}$ initially decreases both with and without an obstacle. When
the loop self-traps around the obstacle $R_{g}$ reaches a constant value
proportional to the size of the obstacle, $R_{g}\approx{r_{o}}$. Without the
obstacle, the radius of gyration decreases to a minimum value of
$R_{g}\approx{0.2}$, then it diverges at longer times when the loop has lost
cohesion forming propagating polar filaments. Moreover, we quantify the mean-
squared displacement $\sigma^{2}$, which is found to transition from ballistic
to diffusive regime in the case without an obstacle, and displays a plateau in
the case with an obstacle. We conclude that self-trapping corresponds to a
pinning mechanism that drastically impacts the transport properties of the
system. Finally, to quantify the average orientation of the particles in the
loop, we obtain the radial orientation component
$e_{r}=\hat{\bm{e}}\cdot\hat{\bm{r}}$ for each particle, where $\hat{\bm{r}}$
is a unitary radial vector with origin at the center of the configuration at
time $t_{1}$, see Fig. 2a,b. Here, $e_{r}=1$ indicates particles are radially
outward oriented, and $e_{r}=-1$ are radially inward oriented, wheres
$e_{r}=0$ corresponds to tangentially oriented particles. We obtain the
probability density $\mathcal{P}(e_{r})$, see Fig. 2e. We find that, in the
case with an obstacle, $\mathcal{P}(e_{r})$ is centered around $e_{r}=0$ and
vanishes at $e_{r}=-1,1$. This shows that trapped particles have an
orientation vector mostly tangential to the obstacle with slight deviations
corresponding to the width of $\mathcal{P}$. Deviations result from
fluctuations due to rotational noise. Conversely, in the case without an
obstacle, the maximum occurs at $e_{r}=-1,1$ as here particles are mostly
propagating radially away from the center of the initial configuration.

We study the order transition of a trapped loop from polar to nematic.
Consider an initially polar closed loop. This loop consists of a certain
number of active particles that can be rotating either clockwise (CW) or
counter-clockwise (CCW). Note that, in contrast to the Vicsek model, the
positional-based alignment defined in Eq. (2) does not align orientations of
neighbouring particles to be the same. In other words, local order can be
either polar or nematic. In Fig. 3a we show a polar loop formed by
$N_{\mathrm{CCW}}<N$ particles, i.e. it contains only a fraction of the total
particles in the system. As the system evolves, more particles can join the
loop, and they can join either in the counter- or clockwise direction.
Therefore, an initially polar loop can become nematic over time with
$N_{\mathrm{CCW}}$ and $N_{\mathrm{CW}}$ particles. In Fig. 3b we show the
time evolution of the number of particles in the loop. We observe that the
number increases over time, with step increases that correspond to incoming
bulky filaments that are absorbed by the loop. The loop in this case starts
with $N_{\mathrm{CCW}}=200$, and $N_{\mathrm{CW}}=0$ particles at $t=2200$. At
a later time $t=4700$ the loop consists of $N_{\mathrm{CCW}}\approx{500}$ and
$N_{\mathrm{CW}}\approx{100}$ particles. We also compute the time evolution of
the local polar order $s_{1}$ as well as local nematic order $s_{2}$ with the
formula

$s_{n}=\bigg{|}\frac{1}{N}\sum_{i=1}^{N}e^{in\theta_{i}}\bigg{|}.$ (3)

We observe a decrease from $s_{1}\approx{0.85}$ to $s_{1}\approx{0.6}$ in the
final configuration. The local nematic order $s_{2}$ decreases only slightly
from $s_{2}\approx{0.5}$ to $s_{2}\approx{0.4}$.

Figure 3: (a) Transition of the loop from initial polar to final nematic. Time
evolution of (b) the number of clockwise $N_{\mathrm{CW}}$ and counter-
clockwise $N_{\mathrm{CCW}}$ particles in the loop, and(c) the local polar
order $s_{1}$, as well as local nematic order $s_{2}$.

We consider a nematic loop of $N=90$ particles placed at an initial distance
of $R_{0}=6$ from the center of a single obstacle. For the standard parameters
considered, namely $v_{0}=1$ and $\gamma=5$, we obtain that the loop shrinks
and gets trapped around the obstacle. We test this scenario for several values
of the turning strength $\gamma$ at $v_{0}=1$, and obtain the radius of
gyration $R_{g}$ for each realization. Results are shown in Fig. 4a.

Figure 4: Radius of gyration for several values of the alignment strength
$\gamma$ and self-propulsion velocity (a) $v_{0}=1$, and (b) $v_{0}=2$.
Horizontal dotted line indicates the minimum value $R_{g}=1.5$ given by the
obstacle size. (c) Shrinking time obtained from (a,b) when the loop reaches
minimum size.

For $v_{0}=1$, we observe that $R_{g}$ decreases and saturates at longer times
when the loop has stabilized around the obstacle, reaching the minimum size of
$R_{g}=1.5$. For larger values of $\gamma$, the radius of gyration shows to
decrease more rapidly. Furthermore, for $v_{0}=2$, we observe that with the
small turning speed $\gamma=4$ shown here, $R_{g}$ monotonously increases in
time, indicating that the loop does not shrink in this case, see Fig. 4b. On
the contrary, the spread of active particles around the obstacle grows, thus
indicating that the loop grows in size instead. For $\gamma=5$, $R_{g}$ slowly
decreases, and does not saturate for the simulation time here considered. To
further characterize the shrinking behavior, we compute the shrinking time
$\tau_{s}$ from $R_{g}(\tau_{b})=1.5$. Obtained values are shown in Fig. 4c,
indicating that the shrinking time remains mostly constant for larger values
of the turning speed, namely for $\gamma>10$. However, $\tau_{s}$
significantly increases for smaller values $\gamma<10$ at any given self-
propulsion speed $v_{0}$.

In the absence of external disturbances, the stability of a particle loop
surrounding an obstacle is not guaranteed. For example, large values of
$D_{\theta}$ can also trigger de-trapping, as it directly influences the
particle orientation $\bm{e}$, allowing them to point radially out of the loop
configuration. When the aperture angle $\alpha$ of the cone of vision is
narrow, particles pointing out of the loop configuration do not perceive other
neighbors and do not co-align to join the loop anymore, instead they escape
the trap. This can trigger de-trapping of the whole structure, as the escaping
particle can serve as an incidental leader which will be followed by the rest
of the particles in the loop. To test particle trapping for different
parameters of the cone of vision as well as noise strength, we perform
simulations of a single loop around one obstacle. We consider several values
of $\alpha$ and $\beta$, both at low and also high values of the noise
strength $\sqrt{2D_{\theta}}$. See results in Fig. 5. We observe that at low
noise, $\sqrt{2D_{\theta}}=0.3$, trapping only occurs for large values of the
obstacle perception angle $\beta\geq{0.5}$, and for intermediate values of the
neighbor perception angle $\alpha$. For high noise, $\sqrt{2D_{\theta}}=0.8$,
we observe trapping only occurs for $\alpha=0.6$ and large $\beta=0.5$.

Figure 5: Single loop trapping diagram for several values of $\alpha$ and
$\beta$ at (a) low noise $\sqrt{2D_{\theta}}=0.2$, and (b) high noise
$\sqrt{2D_{\theta}}=0.8$.

We obtain the cluster size distribution $P(m)$ by performing a clustering
analysis with a cutoff radius of $r_{c}=1.5$, such that neighbors are
considered to be part of a single cluster only when they are close together a
distance $r<r_{c}$. The clustering analysis is averaged over a time interval
of $\Delta{\tau}=100$, during which particles are aggregated into a large
filament in the case without obstacles $N_{o}=0$, or into a dense nematic loop
in the case with obstacles $N_{o}=30$. Results are shown in 6a. Without
obstacles, $P(m)$ shows as a delta distribution around $N=900$, as well an
exponential distribution for clusters of small size $m<10$. With obstacles, a
delta distribution is shown around $N_{o}=600$, as well as a power law
distribution for clusters of small size $m<10^{2}$.

Figure 6: (a) Polar filament found in the case without obstacles, and closed
loops for self-trapping around 1, 2, and 3 obstacles. (b) Cluster size
distribution $P(m)$, and (c) marginal displacement probability density
$\mathcal{P}(x)$ for $N_{o}=0,30$. Cases considered correspond to those in
panel (a). For $P(m)$ a clustering analysis is performed with a cutoff radius
of $r_{c}=r_{o}$. For $\mathcal{P}(x)$, displacements are calculated during a
time interval of $\Delta{t}=L/v_{0}$.

We also obtain the marginal displacement probability density $\mathcal{P}(x)$,
where each displacement is measured during a time interval of $\Delta{t}=15$.
Results are shown in 6b. Without obstacles $N_{o}=0$, $\mathcal{P}$ shows to
be a uniform distribution. With obstacles $N_{o}=30$, the marginal
displacement PDF shows a quadratic dependence $\mathcal{P}\sim{-x^{2}}$ for
small displacements which correspond to the trapped particles in the nematic
loop. Moreover, for larger displacements a Gaussian-like dependence is shown.

A percolating _nematic band_ is formed due to a reorientation mechanism
similar to the filaments case, however here the particles can reorient and
reverse their direction due to the broader vision cone aperture angle
$\alpha=0.6$. Over time the band tends to rectify and percolate, as described
in Ref. [23]. However, in the presence of obstacles, a band can find its way
around an obstacle. In some cases, the bands won’t percolate, but the
structure might close within itself, forming a nematic loop enclosing multiple
obstacles, see Fig. 7.

Figure 7: Representative snapshots of nematic bands formed for two different
obstacle configurations where the band is (a) percolated, and (b) forming a
closed loop around several obstacles.

We obtain local polar and local nematic order parameters, $s_{1}$ and $s_{2}$,
respectively. For parameters ($\alpha,D$)=($0.3\pi$,0.3) corresponding to
filaments, we observe that $s_{1}$ monotonously decreases with the number of
obstacles $N_{o}$, see Fig. 8a. Without obstacles particles aggregate into a
single dense filament, and with increasing number of obstacles the filament
has a larger probability of fragmenting, which explains the decrease in polar
order. For ($\alpha,\sqrt{2D_{\theta}}$)=($0.3\pi$,0.3) corresponding to
bands, $s_{1}$ remains zero at any $N_{o}$, as in this case the bands are
nematic due to the local interaction of between particles. Furthermore, the
local nematic order parameter $s_{2}$ shows monotonously decrease both for
filaments and bands, see Fig. 8b. Note that local polar order contributes to
the overall value of $s_{2}$. For this reason, filaments show the same
behavior for both $s_{1}$ and $s_{2}$, with slightly larger values of $s_{2}$
at larger $N_{o}$, as nematic loops appear for higher $N_{o}$. For bands,
$s_{2}$ is non-vanishing and shows to decrease with increasing $N_{o}$,
indicating the effect of obstacles is to diminish local nematic order due to
the bending around obstacles.

Figure 8: Steady state local order parameters for bands and filaments with
varying number of obstacles $N$. (a) Local polar order. (b) Local nematic
order.

_Discussion_.–Self-trapping of polar filaments depends on the parameters of
the cone of vision $\alpha$, and $\beta$. Such perception parameters are
intrinsic to the particle and independent of the obstacle configuration. In a
configuration of homogeneously distributed obstacles like the one employed
here, particles have the same probability to self-trap around any of the
obstacles. However, a non-homogeneous configuration can be considered in order
to fabricate a filter [11]. When trapping occurs, the loops become absorbent
and a pinning behavior is induced. Interactions between conspecifics lead to
localized nematic alignment within the loops. Note that these interactions do
not favor any type of local order (polar or nematic), instead they tend to
influence the orientations of the particles towards regions of higher particle
density. Trapped loops can destabilize due to different mechanisms, like
external filaments passing by or joining the loop, as well as fluctuations
given by rotational diffusion $D_{\theta}$. Stability can depend on the size
of the loops, as well as the number of surrounded obstacles. Noise is also an
important parameter, which is already known to change the collective behaviour
of the active particles to transition from polar filaments to nematic bands,
and even to form aggregates at higher values [23]. Trapping around multiple
obstacles is more likely to be observed for intermediate values of the noise,
which corresponds to nematic bands. However, bands are also likely to
percolate, and it is still to be verified whether this is a finite-size
effect, or truly long-range order.

Our results corresponds to a first examination of the effects of crowded
environments in active particles with non-reciprocal interactions. Self-
trapping is a novel effect emerging from perception parameters, which could be
present in animals that interact through visual perception, or bacteria with
non-reciprocal interactions.

This work was financially supported by the CONACYT-DAAD scholarship program.

## References

  * Bruno  _et al._ [2011] L. Bruno, M. Salierno, D. E. Wetzler, M. A. Despósito, and V. Levi, PLOS ONE 6, e18332 (2011).
  * Barlan and Gelfand [2017] K. Barlan and V. I. Gelfand, Cold Spring Harb. Perspect. Biol. 9, a025817 (2017).
  * Partridge and Harshey [2013] J. D. Partridge and R. M. Harshey, J. Bacteriol. 195, 909 (2013).
  * Yang  _et al._ [2017] A. Yang, W. S. Tang, T. Si, and J. X. Tang, Biophys. J. 112, 1462 (2017).
  * Altshuler  _et al._ [2013] E. Altshuler, G. Miño, C. Pérez-Penichet, L. del Río, A. Lindner, A. Rousselet, and E. Clément, Soft Matter 9, 1864 (2013).
  * Figueroa-Morales  _et al._ [2015] N. Figueroa-Morales, G. L. Miño, A. Rivera, R. Caballero, E. Clément, E. Altshuler, and A. Lindner, Soft Matter 11, 6284 (2015).
  * Nguyen  _et al._ [2020] J. Nguyen, J. Lara-Gutiérrez, and R. Stocker, FEMS Microbiol. Rev. 45, fuaa068 (2020).
  * Seidler  _et al._ [2015] R. G. Seidler, R. A. Long, J. Berger, S. Bergen, and J. P. Beckmann, Biol. Conserv. 29, 99 (2015).
  * Leclerc  _et al._ [2021] M. Leclerc, M. Leblond, M. Le Corre, C. Dussault, and S. D. Côté, J. Mammal. 102, 1342 (2021).
  * Codina  _et al._ [2022] J. Codina, B. Mahault, H. Chaté, J. Dobnikar, I. Pagonabarraga, and X.-q. Shi, Phys. Rev. Lett. 128, 218001 (2022).
  * Chepizhko and Peruani [2013] O. Chepizhko and F. Peruani, Phys. Rev. Lett. 111, 160604 (2013).
  * Chepizhko  _et al._ [2013] O. Chepizhko, E. G. Altmann, and F. Peruani, Phys. Rev. Lett. 110, 238101 (2013).
  * Mokhtari and Zippelius [2019] Z. Mokhtari and A. Zippelius, Phys. Rev. Lett. 123, 028001 (2019).
  * Zhang  _et al._ [2021] B. Zhang, F. Tan, and N. Zhao, Soft Matter 17, 10334 (2021).
  * Theeyancheri  _et al._ [2022] L. Theeyancheri, S. Chaki, T. Bhattacharjee, and R. Chakrabarti, Phys. Rev. E 106, 014504 (2022).
  * Dinelli  _et al._ [2023] A. Dinelli, J. O’Byrne, A. Curatolo, Y. Zhao, P. Sollich, and J. Tailleur, Nat. Commun. 14, 7035 (2023).
  * Durve  _et al._ [2018] M. Durve, A. Saha, and A. Sayeed, Eur. Phys. J. E 41, 49 (2018).
  * You  _et al._ [2020] Z. You, A. Baskaran, and M. C. Marchetti, PNAS 117, 19767 (2020).
  * Fruchart  _et al._ [2021] M. Fruchart, R. Hanai, P. B. Littlewood, and V. Vitelli, Nature 592, 363 (2021).
  * Kreienkamp and Klapp [2022] K. L. Kreienkamp and S. H. L. Klapp, New J. Phys. 24, 123009 (2022).
  * Knežević  _et al._ [2022] M. Knežević, T. Welker, and H. Stark, Sci. Rep. 12, 19437 (2022).
  * Newman and Sayama [2008] J. P. Newman and H. Sayama, Phys. Rev. E 78, 011913 (2008).
  * Barberis and Peruani [2016] L. Barberis and F. Peruani, Phys. Rev. Lett. 117, 248001 (2016).
  * Stengele  _et al._ [2022] P. Stengele, A. Lüders, and P. Nielaba, Phys. Rev. E 106, 014603 (2022).
  * Qi  _et al._ [2022] J. Qi, L. Bai, Y. Xiao, Y. Wei, and W. Wu, Inf. Sci. 582, 850 (2022).
  * Negi  _et al._ [2022] R. S. Negi, R. G. Winkler, and G. Gompper, Soft Matter 18, 6167 (2022).
  * Negi  _et al._ [2024] R. S. Negi, R. G. Winkler, and G. Gompper, Phys. Rev. Research 6, 013118 (2024).
  * Ginelli  _et al._ [2015] F. Ginelli, F. Peruani, M.-H. Pillot, H. Chaté, G. Theraulaz, and R. Bon, Proc. Nat. Acad. Sci. 112, 12729 (2015).
  * Martinez  _et al._ [2018] R. Martinez, F. Alarcon, D. R. Rodriguez, J. L. Aragones, and C. Valeriani, Eur. Phys. J. E 41, 91 (2018).
  * Lavergne  _et al._ [2019] F. A. Lavergne, H. Wendehenne, T. Bäuerle, and C. Bechinger, Science 364, 70 (2019).
  * Qi  _et al._ [2020] J. Qi, L. Bai, Y. Xiao, W. Wu, and L. Liu, IEEE Access 8, 160490 (2020).
  * Chen and Bechinger [2022] C.-J. Chen and C. Bechinger, New J. Phys. 24, 033001 (2022).
  * Sándor  _et al._ [2017] Cs. Sándor, A. Libál, C. Reichhardt, and C. J. Olson Reichhardt, Phys. Rev. E 95, 032606 (2017).
  * Toner  _et al._ [2018] J. Toner, N. Guttenberg, and Y. Tu, Phys. Rev. Lett. 121, 248002 (2018).
  * Kumar and Mishra [2022] S. Kumar and S. Mishra, Phys. Rev. E 106, 044603 (2022).
  * Duan  _et al._ [2021] Y. Duan, B. Mahault, Y.-q. Ma, X.-q. Shi, and H. Chaté, Phys. Rev. Lett. 126, 178001 (2021).


# Data-driven MHD simulation of successive solar plasma eruptions

Takafumi Kaneko Institute for Space-Earth Environmental Research, Nagoya
University,
Furo-cho, Chikusa-ku, Nagoya, Aichi, 464-8601, Japan  Sung-Hong Park Institute
for Space-Earth Environmental Research, Nagoya University,
Furo-cho, Chikusa-ku, Nagoya, Aichi, 464-8601, Japan  Kanya Kusano Institute
for Space-Earth Environmental Research, Nagoya University,
Furo-cho, Chikusa-ku, Nagoya, Aichi, 464-8601, Japan Takafumi Kaneko
<EMAIL_ADDRESS>
(Received November, 12, 2020; Revised January 10, 2021; Accepted )

###### Abstract

Solar flares and plasma eruptions are sudden releases of magnetic energy
stored in the plasma atmosphere. To understand the physical mechanisms
governing their occurrences, three-dimensional magnetic fields from the
photosphere up to the corona must be studied. The solar photospheric magnetic
fields are observable, whereas the coronal magnetic fields cannot be measured.
One method for inferring coronal magnetic fields is performing data-driven
simulations, which involves time-series observational data of the photospheric
magnetic fields with the bottom boundary of magnetohydrodynamic simulations.
We developed a data-driven method in which temporal evolutions of the
observational vector magnetic field can be reproduced at the bottom boundary
in the simulation by introducing an inverted velocity field. This velocity
field is obtained by inversely solving the induction equation and applying an
appropriate gauge transformation. Using this method, we performed a data-
driven simulation of successive small eruptions observed by the Solar Dynamics
Observatory and the Solar Magnetic Activity Telescope in November 2017. The
simulation well reproduced the converging motion between opposite-polarity
magnetic patches, demonstrating successive formation and eruptions of helical
flux ropes.

Sun: flares—Sun: filaments, prominences—Sun: corona—Sun: photosphere

††journal: ApJ

##  1 Introduction

Solar flares are the sudden releases of energy from the sun. Flares are often
accompanied by plasma eruptions such as prominence eruptions and coronal mass
ejections. Flares and plasma eruptions are caused by the release of magnetic
energy stored in the plasma atmosphere. Evidences of flares and plasma
eruptions have also been found in other sun-like stars (Osten et al., 2005;
Pandey & Singh, 2008; Maehara et al., 2012; Notsu et al., 2019; Namekata et
al., 2020). The sun is the only star of which the photospheric magnetic fields
can be observed with a high spatio-temporal resolution. From the solar
observations and theoretical studies based on magnetohydrodynamic (MHD)
theories, we can infer the detailed magnetic activities leading to the sudden
energy release, which can be common in other sun-like stars. In the current
understanding of solar physics, magnetic reconnection and MHD instabilities
(Hood & Priest, 1979; Kliem & Török, 2006; Ishiguro & Kusano, 2017) are the
essential mechanisms leading to the magnetic energy release.

In the solar observations, the photospheric magnetic fields are temporally
changed via advection by convective flows and magnetic fluxes emerging from
the deeper convection zone. To reveal the mechanisms of the explosive events
and develop methodologies to predict them, previous studies attempted to
evaluate the possibility of magnetic reconnection and the critical conditions
of MHD instabilities (Amari et al., 2014; Kusano et al., 2020). For this, the
information of three-dimensional magnetic fields from the photosphere up to
the corona is required. The photospheric magnetic fields can be observed,
whereas the coronal magnetic fields cannot be measured directly. Previous
studies developed numerical methods to extrapolate three-dimensional coronal
magnetic fields from the two-dimensional observational vector magnetic fields
in the photosphere, e.g., nonlinear force-free field (NLFFF) approximation
(reviewed by Inoue, 2016). There have been attempts of data-constrained
simulations, where the NLFFF approximation was used as the initial condition
of MHD simulation (Amari et al., 2014; Muhamad et al., 2017). In these
simulations, the photospheric magnetic fields after temporal integration did
not always reproduce the observed ones. Another attempt was the data-driven
simulation in which time-series photospheric magnetic data were involved in
the bottom boundary of MHD simulations. The expected advantage of the data-
driven methods, compared with the NLFFF or data-constrained model, is that the
results are free from the assumption of force-free field. We can follow more
realistic temporal evolution of coronal magnetic fields as a response of
temporal change of the observational photospheric magnetic fields. Several
data-driven MHD simulations have been performed, and their results agree with
some aspects in the observations, e.g., morphology of the coronal magnetic
loops (Cheung & DeRosa, 2012; Cheung et al., 2015; Jiang et al., 2016; Hayashi
et al., 2018, 2019; Pomoell et al., 2019; Guo et al., 2019; He et al., 2020).
In contrast, a recent comparative study by Toriumi et al. (2020) reported that
the numerical solutions obtained from the different data-driven simulations
using the same time-series magnetic data were different from each other. The
data-driven methods must be improved further to resolve these discrepancies.

In this study, we focus on the velocity fields in the bottom boundary of MHD
simulation. In several data-driven methods, the velocity fields at the bottom
boundary were set to be zero, leading to physical inconsistency between
velocity fields and electric or magnetic fields in terms of the induction
equation. A recent study by Hayashi et al. (2019) combined the velocity fields
derived from differential affine velocity estimator for vector magnetograms
(DAVE4VM; Schuck, 2006) with their own data-driven method (denoted as the
$v$-driven method in their paper). They confirmed that the frozen-in condition
between plasmas and magnetic fields was well established. This is because the
DAVE4VM-inferred velocity works as the bottom boundary condition to the
equation of motion in MHD simulation, providing the motion of plasmas coherent
with the time evolution of the magnetic fields in the observation. Another
recent study by Guo et al. (2019) reported that the numerical results of data-
driven simulations with and without velocity fields by DAVE4VM were similar in
terms of morphology and propagation path of the erupted flux ropes. They
argued that eruption inevitably happens if the initial condition of MHD
simulation is already close to the dynamic eruptive phase. Their conclusion
was that the change of the bottom boundary condition had subtle effect to the
onset mechanism, while it would affect magnetic energy build-up before
eruptive phase. Since the observational targets and many aspects of numerical
techniques were different in Hayashi et al. (2019) and Guo et al. (2019), it
is fairly difficult to compare their results. One issue we concern about their
numerical techniques is that the DAVE4VM-inferred velocity is not always
consistent with the inverted electric fields or the time evolution of the
observational magnetic fields used as the bottom boundary condition of MHD
simulations. In the present study, we attempted to implement an inversion
technique of the induction equation directly in our simulation code, and
proposed a method to derive the velocity fields reproducing the observed time
evolution of magnetic field as a numerical solution of MHD equations. To
confirm the feasibility of the method, we applied it to the successive small
eruptive events that occurred in November 2017.

The observation of the eruptive events are described in Section 2. The
numerical method including velocity inversion is described in Section 3. The
numerical results are shown in Section 4. We summarize and discuss the results
in Section 5

##  2 Observation

The Solar Dynamics Doppler Imager (SDDI; Ichimoto et al., 2017) installed on
the Solar Magnetic Activity Research Telescope (SMART; UeNo et al., 2004) at
Hida Observatory of Kyoto University provides full-disk solar images at
multiple wavelengths around the H$\mathrm{\alpha}$ 6563 Å line with a 0.25 Å
bandpass. The top and bottom panels of Figure 1 show H$\mathrm{\alpha}$ blue
wing images at $-$0.5 Å from the line center and H$\mathrm{\alpha}$ line
center images, respectively, in six snapshots taken from SMART/SDDI
observations on November 4–5, 2017, which demonstrate two successive
eruptions. As indicated by the arrow in panel (a2) of Figure 1, the first
eruption event started at 23:40 UT on November 4 and appeared as a compact
dark feature with a size of $\sim$10″ in H$\mathrm{\alpha}-$0.5 Å. This dark
feature (i.e., a so-called H$\mathrm{\alpha}$ upflow event), which is only
visible in the H$\mathrm{\alpha}$ blue wing, displays an upward motion and is
known to be often associated with magnetic reconnection (Wang et al., 1998;
Chae et al., 1998). In a sequence of H$\mathrm{\alpha}-$0.5 Å images, it was
found that the H$\mathrm{\alpha}-$0.5 Å upflow features increase in size as
they erupt in the south-west direction (see panel (a3)). During the eruption,
an enhanced brightening in H$\mathrm{\alpha}$ was observed near the magnetic
polarity inversion line (PIL), where two opposite-polarity magnetic patches
approached each other. Approximately 3 h after the first eruption, another
eruption event began at 02:50 UT on November 5, exhibiting characteristics
similar to those of the first eruption event in the context of the south-west
eruption direction and H$\mathrm{\alpha}$ brightening. In the case of the
second eruption, contrarily, we note that an inverse S-shaped structure is
clearly seen in the H$\mathrm{\alpha}$ line center (refer to panel (b6)).

Figure 1: Two successive eruption events observed in H$\mathrm{\alpha}-$0.5 Å
images (top panels) and H$\mathrm{\alpha}$ line center images (bottom panels).
The arrows in panels (a2) and (a3) indicate the first event, while those in
panels (a5) and (a6) indicate the second event. In each panel, the red and
blue contours represent $\pm$50 G of the vertical magnetic field $B_{z}$. The
yellow box in panel (b1) represents the region used as the bottom boundary in
our MHD simulation.

In this study, we used a sequence of photospheric vector magnetograms obtained
at 12-min cadence by the Helioseismic and Magnetic Imager (HMI; Schou et al.,
2012) onboard the Solar Dynamics Observatory (SDO; Pesnell et al., 2012). The
pixel size of the HMI vector magnetograms was $\sim$360 km. Figure 2 shows two
co-aligned images of the vertical ($B_{z}$) and horizontal ($B_{x}$ and
$B_{y}$) components of the photospheric magnetic field at 22:58 UT on November
4, 2017 (left column) and at 00:58 UT on November 5, 2017 (right column). The
field-of-view of the co-aligned magnetic field images is marked by the yellow
box in panel (b1) of Figure 1, which contains the magnetic source region that
produced the two eruptions. The source region consists of two main opposite-
polarity magnetic patches that, in general, showed a converging motion as well
as a decrease in the magnetic flux of both polarities over a 5 h interval
around the times that the two eruptions occurred. Moreover, as shown in panels
(d) and (f) of Figure 2, the strengths of the horizontal components $B_{x}$
and $B_{y}$ are found to increase after the first eruption.

Figure 2: Snapshots of time-series data of magnetic field observed by SDO/HMI.
Panels (a), (c), and (e) show snapshots of 22:58 UT on November 4, 2017
(corresponding to $t=0$ in the simulation). Panels (b), (d), and (f) show
snapshots of 00:58 UT on November 5, 2017 (corresponding to
$t=120~{}\mathrm{min}$ in the simulation). The field-of-view of these figures
is represented by the yellow box in Fig. 1 (b1)

##  3 Numerical Method

We numerically solved the zero-beta MHD equations as follows.

$\frac{\partial\rho}{\partial
t}+\nabla\cdot\left(\rho\mbox{\boldmath$v$}\right)=0,$ (1)
$\frac{\partial\left(\rho\mbox{\boldmath$v$}\right)}{\partial
t}+\nabla\cdot\left(\rho\mbox{\boldmath$v$}\mbox{\boldmath$v$}+\frac{B^{2}}{\
8\pi}\mbox{\boldmath$I$}-\frac{\mbox{\boldmath$B$}\mbox{\boldmath$B$}}{4\pi}\right)=0,$
(2) $\frac{\partial\mbox{\boldmath$B$}}{\partial
t}=\nabla\times(\mbox{\boldmath$v$}\times\mbox{\boldmath$B$}-\eta\mbox{\boldmath$J$}),$
(3) $\mbox{\boldmath$J$}=\frac{1}{4\pi}\nabla\times\mbox{\boldmath$B$},$ (4)

where $t$, $\rho$, $v$, $B$, $J$, $\eta$, and $I$ denote time, mass density,
velocity fields, magnetic fields, current density, resistivity, and unit
vector, respectively. We used the anomalous resistivity in the following form:

$\eta=0,~{}~{}(J<J_{c}),$ (5) $\eta=\eta_{0}(J/J_{c}-1)^{2},~{}~{}(J\geq
J_{c}),$ (6)

where $J_{c}=10^{-9}~{}\mathrm{G/cm}$, $\eta_{0}=10^{11}~{}\mathrm{cm^{2}/s}$,
and we restrict $\eta\leq\eta_{\mathrm{max}}=10^{11}~{}\mathrm{cm^{2}/s}$.

We inverted velocity fields which reproduce the observational photospheric
magnetic fields by solving Eq. (3), and implemented them in the bottom
boundary layer of the MHD simulation. The inverted velocity fields were
computed by the following three steps:

  1. 1.

Inversion of the induction equation
We solved an inverse problem of the induction equation as, in principle,

$\frac{\mbox{\boldmath$B$}_{\mathrm{obs}}^{n+1}-\mbox{\boldmath$B$}_{\mathrm{obs}}^{n}}{\tau}=-\nabla\times\mbox{\boldmath$E$}^{I},$
(7)

where $\mbox{\boldmath$B$}_{\mathrm{obs}}^{n}$, $\tau$, and
$\mbox{\boldmath$E$}^{I}$ represent the $n$-th snapshot in the time-series
data of the observational magnetic fields, the temporal cadence of the HMI
observation, and an inverted electric field, respectively. As pointed out in a
previous study (Kusano et al., 2002), we cannot solve this inverse problem
completely because Eq. (7) includes the derivative in the $z$-direction (the
direction normal to the photosphere), whereas the observational magnetic data
have only two-dimensional information in the $xy$-plane (corresponding to the
solar surface). Several methods have been proposed to resolve this problem. In
this study, we adopt the poloidal-toloidal decomposition method (Fisher et
al., 2010) and obtain $\mbox{\boldmath$E$}^{I}$. The advantage of this method
is that we can estimate the vertical derivative of electric fields to some
extent. However, the complete solution cannot be obtained even by this method.
We carried out the inversion of the electric fields between the simulated
magnetic fields and the observational magnetic fields during the observational
time cadence:

$\frac{\mbox{\boldmath$B$}_{\mathrm{obs}}^{n+1}-\mbox{\boldmath$B$}_{\mathrm{sim}}^{n+m/M}}{(1-m/M)\tau}=-\nabla\times\mbox{\boldmath$E$}^{I},$
(8)

where $\mbox{\boldmath$B$}_{\mathrm{sim}}$ denotes the simulated magnetic
fields, $m=1,2,...,M-1$ represents the $m$-th sub-snapshot between the $n$-th
and the $(n+1)$-th observational snapshots. We adopted $M=6$ in this study,
hence, the inversion was performed every 2 min during 12-min observational
cadence of HMI. This piecewise inversion technique increases feasibility of
the observational magnetic fields compared with the case that electric fields
are inverted only once between $B_{\mathrm{obs}}^{n}$ and
$B_{\mathrm{obs}}^{n+1}$.

  2. 2.

Gauge transformation
Electric field is mathematically gauge-invariant to the induction equation; we
can add an arbitrary scalar potential $\phi$ in the following form.

$\mbox{\boldmath$E$}=\mbox{\boldmath$E$}^{I}-\nabla\phi,$ (9)

where $E$ represents the electric field after gauge transformation. In
contrast, as demonstrated by Pomoell et al. (2019), the results of data-driven
simulations are influenced by gauge transformation. We adopted a gauge
transformation that satisfied $\mbox{\boldmath$E$}\cdot\mbox{\boldmath$B$}=0$
using the iterative approach in Fisher et al. (2010). The motivation to use
this gauge transformation is as follows: the electric fields defined by
$\mbox{\boldmath$E$}=-\mbox{\boldmath$v$}\times\mbox{\boldmath$B$}$ are always
perpendicular to magnetic fields. In contrast, the inverted electric fields
$\mbox{\boldmath$E$}^{I}$ before the gauge transformation usually contain
nonzero $\mbox{\boldmath$E$}_{\parallel}$ (parallel component to $B$). The
nonzero $\mbox{\boldmath$E$}_{\parallel}$ can cause mismatch of the electric
fields between the bottom boundary layer and the main simulation domain
because the electric fields in the main simulation domain are computed as
$\mbox{\boldmath$E$}=-\mbox{\boldmath$v$}\times\mbox{\boldmath$B$}$. Thus, we
assumed that $\mbox{\boldmath$E$}\cdot\mbox{\boldmath$B$}=0$ is a necessary
condition for the boundary electric fields in data-driven MHD simulations.
Note that this assumption is valid even if resistive term
$\eta\mbox{\boldmath$J$}$ was introduced because the magnitude of
$\eta\mbox{\boldmath$J$}$ is constrained much smaller than that of
$-\mbox{\boldmath$v$}\times\mbox{\boldmath$B$}$ in MHD simulations.

  3. 3.

Derivation of velocity fields
We compute velocity fields as follows:

$\mbox{\boldmath$v$}^{I}=\frac{\mbox{\boldmath$E$}\times\mbox{\boldmath$B$}}{B^{2}},$
(10)

where $\mbox{\boldmath$v$}^{I}$ represents the inverted velocity field. We
substituted $\mbox{\boldmath$B$}_{\mathrm{sim}}$ to $B$ in Eq. (10).
$\mbox{\boldmath$v$}^{I}$ was updated every numerical time step.

In Step 3, in the case that $E$ contains $\mbox{\boldmath$E$}_{\parallel}$,
$\mbox{\boldmath$v$}^{I}$ loses the information of
$\mbox{\boldmath$E$}_{\parallel}$ (because
$\mbox{\boldmath$E$}_{\parallel}\times\mbox{\boldmath$B$}=0$). In our
manipulations, in Step 2, $\mbox{\boldmath$E$}_{\parallel}$ has already been
eliminated by the gauge transformation. The inductive electric fields were
calculated in a part of the right-hand side of Eq. (3) as follows:

$\mbox{\boldmath$v$}^{I}\times\mbox{\boldmath$B$}=\frac{(\mbox{\boldmath$E$}\times\mbox{\boldmath$B$})\times\mbox{\boldmath$B$}}{B^{2}}=-\mbox{\boldmath$E$}+\frac{(\mbox{\boldmath$E$}\cdot\mbox{\boldmath$B$})\mbox{\boldmath$B$}}{B^{2}}=-\mbox{\boldmath$E$}.$
(11)

Thus, we expect that the observed photospheric magnetic fields are reproduced
as a self-consistent numerical solution of Eq. (3) only by introducing
$\mbox{\boldmath$v$}^{I}$ in the bottom boundary. To reduce the observational
noise in the area of weak magnetic fields, which can damage the inversion of
electric fields and the gauge transformation, we applied a low-pass filter
using FFT to the original magnetic data. The practical spatial resolution was
8 times lower than the original spatial resolution of the HMI.

The simulation domain is a rectangular box. Its Cartesian coordinates
$(x,y,z)$ are extended to $0<x<89.6~{}\mathrm{Mm}$, $0<y<89.6~{}\mathrm{Mm}$,
and $-1.44~{}\mathrm{Mm}<z<73.8~{}\mathrm{Mm}$, respectively, where the
$xy$-plane is the horizontal plane parallel to the solar surface, and the
$z$-direction represents the height. We adopted uniform grid spacing in every
direction, and the grid size was $360~{}\mathrm{km}$ corresponding to the
spatial resolution of the HMI. Below the $z=0$ plane, we set 5 grids in
$z$-direction where Eq. (3) was numerically solved by introducing
$\mbox{\boldmath$v$}^{I}$. Note that only the horizontal derivatives were
calculated in the lowest 2 grids. The $z=-360~{}\mathrm{km}$ plane (one grid
below $z=0$) is at the height where the observational magnetic fields were
expected to be reproduced. The method of Fisher et al. (2010) derives
$\partial_{z}E_{x}$ and $\partial_{z}E_{y}$. Assuming $\partial_{z}E_{z}=0$,
we linearly extrapolated the inverted electric fields in the $z$-direction
below $z=0$ and computed $\mbox{\boldmath$v$}^{I}$ with the local magnetic
fields using Eq. (10). The density was fixed to the initial values below
$z=0$. We adopted free boundary condition to the top boundary and fixed to the
side boundaries. Our simulation only included the corona with typical density
of $10^{9}~{}\mathrm{cm^{-3}}$, which is much smaller than the typical
photospheric density $10^{17}~{}\mathrm{cm^{-3}}$. To suppress the
unrealistically fast Alfvèn speed, we reduce the magnetic field strength to be
10 times smaller than the original observed values. The same modification was
also adopted in Jiang et al. (2016).

The initial condition was a potential field computed by the Fourier expansion
method (Priest, 2014) from the vertical magnetic field at 22:58 UT on November
4, 2017, observed by HMI. The initial density was given by
$\rho=\rho_{0}\exp[-z/H]$, where $\rho_{0}=3.2\times
10^{-15}~{}\mathrm{g/cm^{-3}}$ and $H=3.0\times 10^{4}~{}\mathrm{km}$. The
numerical scheme used was a four-step Runge-Kutta method (Jameson, 2017) and a
fourth order central finite difference method with an artificial viscosity
(Rempel, 2014).

##  4 Results

Figure 3 shows snapshots of magnetic fields in the bottom boundary at the
height $z=-360~{}\mathrm{km}$, where the observed magnetic fields were
expected to be reproduced. Note that the magnetic fields shown in Fig. 3 are
the numerically obtained solutions, not merely smoothed observational data.
Compared with Fig. 2, we confirmed that the converging motion of the opposite-
polarity magnetic patches and the intrusion of the negative patch to the
positive patch were well reproduced in our simulation. The small structures
were smoothed out by the low-pass filter used for the inversion and the
anomalous resistivity during the temporal integration of the MHD simulation.
The structural similarity (SSIM) values (Wang et al., 2004) between the raw
observational data and the low-pass filtered observational data at
$t=120~{}\mathrm{min}$ were 0.22, 0.12, and 0.58 for the components $B_{x}$,
$B_{y}$, and $B_{z}$, respectively, and the SSIM values between the low-pass
filtered observational data and the simulated data were 0.82, 0.61, and 0.93
for the components $B_{x}$, $B_{y}$, and $B_{z}$, respectively. We used the
low-pass filtered data for calculation of $\mbox{\boldmath$v$}^{I}$. We
confirmed that the inverted velocities work well because the given magnetic
fields were reproduced with high accuracy.

Figure 3: Snapshots of time evolution of magnetic fields at bottom boundary.
Panels (a), (c), and (e) show snapshots of $t=0$ (corresponding to Nov. 4.
2017, 22:58 UT in the observation). Panels (b), (d), and (f) show snapshots of
$t=120~{}\mathrm{min}$ (corresponding to Nov. 5. 2017, 0:58 UT in the
observation).

As the opposite-polarity patches converged with each other and the negative
magnetic patch further trespassed into the positive magnetic patch, the
formation and eruption of flux ropes via reconnection successively occurred in
our simulation. Figure 4 shows the temporal evolution of the three-dimensional
magnetic field in the corona. Figures 4 (a) and (b) show snapshots of when the
horizontal magnetic fields were shifted from the potential fields to the
observed ones at 22:58 UT on November 4, 2017. Figures 4 (c) and (d) show
snapshots of the first eruption. Figures 4 (e) and (f) show snapshots of the
second eruption. We succeeded in reproducing the successive eruptions of flux
ropes. In both cases, the flux ropes erupted in the south-west direction. We
can interpret that the erupting filamentary structures in the observation were
manifestations of the erupting flux ropes. Compared with H$\mathrm{\alpha}$
blue wing images in the observation (see Fig. 1 (a3) and (a6)), the direction
of eruptions in the simulation were in agreement with the observational
results.

Figure 4: Temporal evolution of the coronal magnetic fields. Lines and colors
on the lines represent the magnetic field lines and vertical velocity,
respectively. Blue and red represent upward and downward velocities,
respectively. The grayscale on the bottom surface represents the vertical
magnetic fields.

Figure 5 (a) shows the temporal evolution of the kinetic energy integrated in
the simulation domain over the $z=0$ plane. The rapid increase in kinetic
energy at $t\sim 100~{}\mathrm{min}$ and $t\sim 200~{}\mathrm{min}$ represents
the eruptions. The onset time of the first eruption in the simulation was
delayed $40~{}\mathrm{min}$ compared with the observational results, and that
of the second eruption was $40~{}\mathrm{min}$ earlier. Figure 5 (b) shows the
temporal evolution of the nonpotential magnetic energy computed as difference
of the total magnetic energy and the potential magnetic energy. The rapid
increase of kinetic energy temporally coincided with the reduction of
nonpotential magnetic energy. The kinetic energy was approximately $1\times
10^{25}~{}\mathrm{erg}$ and the released magnetic energy was $2-3\times
10^{25}~{}\mathrm{erg}$. Note that the magnetic energy can change also due to
the energy flux at the bottom and the top boundaries. Because we reduced the
magnetic field strength in the simulation to 10 times smaller than the
original observational values, we can speculate that the actual energy release
was of the order of $10^{27}~{}\mathrm{erg}$ for these eruptive events (the
right axis of Fig. 5 (a)). This is because magnetic energy is proportional to
$B^{2}$, and we solved scale-free MHD equations.

Figure 5: The solid, dashed, and dash-dotted lines in panel (a) represent the
kinetic energy of all velocity components
$\int_{z>0}\frac{1}{2}\rho(v_{x}^{2}+v_{y}^{2}+v_{z}^{2})dV$, of the vertical
component $\int_{z>0}\frac{1}{2}\rho v_{z}^{2}dV$, and of the horizontal
components $\int_{z>0}\frac{1}{2}\rho(v_{x}^{2}+v_{y}^{2})dV$, respectively.
The vertical dotted lines indicate the actual onset times of the eruptions in
the observation. The right axis in panel (a) represents the estimated released
energy in reality. The solid and dashed lines in panel (b) represent the
nonpotential magnetic energy and the kinetic energy of all velocity
components, respectively.

##  5 Summary and Discussion

We developed a numerical methodology that reproduces the temporal evolution of
observational vector magnetic fields in the bottom boundary of MHD simulation
by introducing the velocity fields inverted from the time-series observational
magnetic data ($\mbox{\boldmath$E$}\times\mbox{\boldmath$B$}$-driven method).
The inverted velocity fields were computed by the formula of
$\mbox{\boldmath$E$}\times\mbox{\boldmath$B$}$ drift. The gauge transformation
of electric field satisfying $\mbox{\boldmath$E$}\cdot\mbox{\boldmath$B$}=0$
enables this simple formulation. In the previous data-driven simulations, the
velocity fields inferred by DAVE4VM were introduced in addition with the
inverted electric fields or the observational magnetic fields (Hayashi et al.,
2019; Guo et al., 2019). The DAVE4VM inversely solves the induction equation
as well, whereas the vertical derivative of the horizontal components of
electric fields is assumed negligible. In practice, the observed time
evolution of the magnetic fields is not always reproduced completely even if
the induction equation was integrated along with the DAVE4VM-inferred velocity
fields. Therefore, the previous data-driven simulations used the DAVE4VM-
inferred velocity as the bottom boundary condition of the equation of motion,
and the observational magnetic fields or the inverted electric fields for the
induction equation. In the present simulation, the observational magnetic
fields were reproduced only by introducing the inverted velocity fields. The
bottom boundary condition for the equation of motion and the induction
equation is both physically and numerically consistent in our method. We
applied our method to the successive eruptive events. Our simulation succeeded
in reproducing the successive formation and eruption of flux ropes as a
response of temporal evolution of the observational magnetic fields. The
inversion technique of electric fields by Fisher et al. (2010) used in Step 1,
described in Section 3, was widely used in the previous studies (Pomoell et
al., 2019). The computation of the inverted velocity fields in Step 3 is
straightforward. Our method is simple yet feasible to reproduce magnetic
activities in the solar atmosphere.

The energy release of flares in the solar active regions, the typical field
strength of which is several thousand gauss, is in the range of
$10^{28}-10^{32}~{}\mathrm{erg}$. The field strength of the magnetic patches
in our study was approximately one hundred gauss. The estimation of the
released energy of $10^{27}~{}\mathrm{erg}$ from our simulation result is
plausible for small eruptive events.

The successive flares and eruptions from active regions are often reported.
The largest flare in Solar Cycle 24, which marked X9.3 in GOES X-ray
classification, also had the preceding X2.2 flare. The flares were triggered
by the continuous intrusion of the opposite-polarity magnetic fluxes,
according to the analyses by Bamba et al. (2020). In our case, although the
spatial size and magnetic field strength were much smaller, it is common that
the continuous convergence of the opposite magnetic flux and the subsequent
partial deformation of the PIL were the triggers of the successive eruptions.
It is worthy to note that the triggering mechanism of successive eruptions was
similar over a wide range of different spatial and temporal scales. The
previous theoretical studies (Kusano et al., 2012, 2020) also support that the
partial deformation of PILs can trigger eruption. Kusano et al. (2020)
discussed how a small area of reconnection can lead to MHD instability. We
speculate that the local converging motion can create small reconnetion-favor
(opposite-polarity or reversed-shear type) regions along the PILs in the
active regions. In contrast, the origin of continuous converging motion is
still unclear. It is also unclear whether the converging motion is
concentrated nearby PILs or ubiquitous in the photosphere. The ultra high
resolution observations by the Daniel K. Inouye Solar Telescope may reveal
this issue. A comprehensive understanding of the photospheric motion coupling
with magnetic activity in the convection zone is also required (Cheung et al.,
2019; Hotta et al., 2019; Toriumi & Hotta, 2019).

The onset times of the eruptions in our simulations were differed by
$40~{}\mathrm{min}$ compared with the observational ones. A possible reason is
that the small-scale structures were smoothed out in our simulation. As
mentioned in Section 4, the SSIMs of the magnetic fields after applying the
low-pass filter were already low. The smaller structures may have to be
included as much as possible to reproduce the accurate onset times. The
anomalous resistivity might also affect the results. A parameter survey on the
anomalous resistivity must be conducted in future work.

We are grateful to the anonymous referee for the constructive and thoughtful
comments. This work was supported by MEXT/JSPS KAKENHI grant number
JP15H05814, Project for Solar-Terrestrial Environmental Prediction (PSTEP),
and JSPS KAKENHI grant number JP20K14519. This work was partially supported by
MEXT as ”Program for Promoting Researches on the Supercomputer Fugaku” (Toward
a unified view of the universe: from large scale structure to planets,
Elucidation of solar and planetary dynamics and evolution). Numerical
computations were conducted on a Cray XC50 supercomputer at the Center for
Computational Astrophysics (CfCA) of the National Astronomical Observatory of
Japan. A part of this study was carried out using the computational resource
of the Center for Integrated Data Science, Institute for Space-Earth
Environmental Research, Nagoya University. HMI is an instrument on the SDO, a
mission for NASA’s Living with a Star program. We are grateful to the staff of
Hida Observatory for supporting the instrument development and daily
observations.

## References

  * Amari et al. (2014) Amari, T., Canou, A., & Aly, J.-J. 2014, Nature, 514, 465, doi: 10.1038/nature13815
  * Bamba et al. (2020) Bamba, Y., Inoue, S., & Imada, S. 2020, ApJ, 894, 29, doi: 10.3847/1538-4357/ab85ca
  * Chae et al. (1998) Chae, J., Wang, H., Lee, C.-Y., Goode, P. R., & Schühle, U. 1998, ApJ, 504, L123, doi: 10.1086/311583
  * Cheung & DeRosa (2012) Cheung, M. C. M., & DeRosa, M. L. 2012, ApJ, 757, 147, doi: 10.1088/0004-637X/757/2/147
  * Cheung et al. (2015) Cheung, M. C. M., De Pontieu, B., Tarbell, T. D., et al. 2015, ApJ, 801, 83, doi: 10.1088/0004-637X/801/2/83
  * Cheung et al. (2019) Cheung, M. C. M., Rempel, M., Chintzoglou, G., et al. 2019, Nature Astronomy, 3, 160, doi: 10.1038/s41550-018-0629-3
  * Fisher et al. (2010) Fisher, G. H., Welsch, B. T., Abbett, W. P., & Bercik, D. J. 2010, ApJ, 715, 242, doi: 10.1088/0004-637X/715/1/242
  * Guo et al. (2019) Guo, Y., Xu, Y., Ding, M. D., et al. 2019, ApJ, 884, L1, doi: 10.3847/2041-8213/ab4514
  * Hayashi et al. (2018) Hayashi, K., Feng, X., Xiong, M., & Jiang, C. 2018, ApJ, 856, 181, doi: 10.3847/1538-4357/aab787
  * Hayashi et al. (2019) —. 2019, ApJ, 871, L28, doi: 10.3847/2041-8213/aaffcf
  * He et al. (2020) He, W., Jiang, C., Zou, P., et al. 2020, ApJ, 892, 9, doi: 10.3847/1538-4357/ab75ab
  * Hood & Priest (1979) Hood, A. W., & Priest, E. R. 1979, Sol. Phys., 64, 303, doi: 10.1007/BF00151441
  * Hotta et al. (2019) Hotta, H., Iijima, H., & Kusano, K. 2019, Science Advances, 5, 2307, doi: 10.1126/sciadv.aau2307
  * Ichimoto et al. (2017) Ichimoto, K., Ishii, T. T., Otsuji, K., et al. 2017, Sol. Phys., 292, 63, doi: 10.1007/s11207-017-1082-7
  * Inoue (2016) Inoue, S. 2016, Progress in Earth and Planetary Science, 3, 19, doi: 10.1186/s40645-016-0084-7
  * Ishiguro & Kusano (2017) Ishiguro, N., & Kusano, K. 2017, ApJ, 843, 101, doi: 10.3847/1538-4357/aa799b
  * Jameson (2017) Jameson, A. 2017, AIAA Journal, 1487
  * Jiang et al. (2016) Jiang, C., Wu, S. T., Yurchyshyn, V., et al. 2016, ApJ, 828, 62, doi: 10.3847/0004-637X/828/1/62
  * Kliem & Török (2006) Kliem, B., & Török, T. 2006, Phys. Rev. Lett., 96, 255002, doi: 10.1103/PhysRevLett.96.255002
  * Kusano et al. (2012) Kusano, K., Bamba, Y., Yamamoto, T. T., et al. 2012, ApJ, 760, 31, doi: 10.1088/0004-637X/760/1/31
  * Kusano et al. (2020) Kusano, K., Iju, T., Bamba, Y., & Inoue, S. 2020, Science, 369, 587, doi: 10.1126/science.aaz2511
  * Kusano et al. (2002) Kusano, K., Maeshiro, T., Yokoyama, T., & Sakurai, T. 2002, ApJ, 577, 501, doi: 10.1086/342171
  * Maehara et al. (2012) Maehara, H., Shibayama, T., Notsu, S., et al. 2012, Nature, 485, 478, doi: 10.1038/nature11063
  * Muhamad et al. (2017) Muhamad, J., Kusano, K., Inoue, S., & Shiota, D. 2017, ApJ, 842, 86, doi: 10.3847/1538-4357/aa750e
  * Namekata et al. (2020) Namekata, K., Maehara, H., Sasaki, R., et al. 2020, PASJ, 72, 68, doi: 10.1093/pasj/psaa051
  * Notsu et al. (2019) Notsu, Y., Maehara, H., Honda, S., et al. 2019, ApJ, 876, 58, doi: 10.3847/1538-4357/ab14e6
  * Osten et al. (2005) Osten, R. A., Hawley, S. L., Allred, J. C., Johns-Krull, C. M., & Roark, C. 2005, ApJ, 621, 398, doi: 10.1086/427275
  * Pandey & Singh (2008) Pandey, J. C., & Singh, K. P. 2008, MNRAS, 387, 1627, doi: 10.1111/j.1365-2966.2008.13342.x
  * Pesnell et al. (2012) Pesnell, W. D., Thompson, B. J., & Chamberlin, P. C. 2012, Sol. Phys., 275, 3, doi: 10.1007/s11207-011-9841-3
  * Pomoell et al. (2019) Pomoell, J., Lumme, E., & Kilpua, E. 2019, Sol. Phys., 294, 41, doi: 10.1007/s11207-019-1430-x
  * Priest (2014) Priest, E. 2014, Magnetohydrodynamics of the Sun (Cambridge University Press), doi: 10.1017/CBO9781139020732
  * Rempel (2014) Rempel, M. 2014, ApJ, 789, 132, doi: 10.1088/0004-637X/789/2/132
  * Schou et al. (2012) Schou, J., Scherrer, P. H., Bush, R. I., et al. 2012, Sol. Phys., 275, 229, doi: 10.1007/s11207-011-9842-2
  * Schuck (2006) Schuck, P. W. 2006, ApJ, 646, 1358, doi: 10.1086/505015
  * Toriumi et al. (2020) Toriumi, S., Takasao, S., Cheung, M. C. M., et al. 2020, ApJ, 890, 103, doi: 10.3847/1538-4357/ab6b1f
  * Toriumi & Hotta (2019) Toriumi, S., & Hotta, H. 2019, ApJ, 886, L21, doi: 10.3847/2041-8213/ab55e7
  * UeNo et al. (2004) UeNo, S., Nagata, S.-i., Kitai, R., Kurokawa, H., & Ichimoto, K. 2004, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 5492, Ground-based Instrumentation for Astronomy, ed. A. F. M. Moorwood & M. Iye, 958–969, doi: 10.1117/12.550304
  * Wang et al. (1998) Wang, H., Johannesson, A., Stage, M., Lee, C., & Zirin, H. 1998, Sol. Phys., 178, 55, doi: 10.1023/A:1004974927114
  * Wang et al. (2004) Wang, Z., Bovik, A. C., Sheikh, H. R., & Simoncelli, E. P. 2004, IEEE transactions on image processing, 13, 600