single = false, list/sort = true, cite/group = true, cite/group/cmd =
cite/group/pre = ; , EVE short = EVE, long = Extreme Ultraviolet Variability
Experiment, cite = Woods:2012 AIA short = AIA, long = Atmospheric Imaging
Assembly, cite = Lemen:2012 HMI short = HMI, long = Helioseismic Magnetic
Imager, cite = Scherrer:2012 SDO short = SDO, long = Solar Dynamics
Observatory, cite = Pesnell:2012 STEREO short = STEREO, long = Solar
Terrestrial Relations Observatory, cite = Kaiser:2008 SECCHI short = SECCHI,
long = Sun Earth Connection Coronal and Heliospheric Investigation, cite =
Howard:2008 SCIP short = SCIP, long = Sun Centered Imaging Package MGN short =
MGN, long = multi-scale Gaussian normalisation, cite = Morgan:2014 EUVI short
= EUVI, long = EUVI CME short = CME, short-plural-form = CMEs, long = coronal
mass ejection, long-plural-form = coronal mass ejections PIL short = PIL,
short-plural-form = PILs, long = polarity inversion line, long-plural-form =
polarity inversion lines LOS short = LOS, short-plural-form = LOS, long = line
of sight, long-plural-form = lines of sight FOV short = FOV, short-plural-form
= FOV, long = field of view, long-plural-form = fields of view UV short = UV,
short-plural-form = UV, long = ultraviolet, long-plural-form = ultraviolet EUV
short = EUV, short-plural-form = EUV, long = extreme ultraviolet, long-plural-
form = extreme ultraviolet GONG short = GONG, long = Global Oscillations
Network Group DST short = DST, long = Dunn Solar Telescope IBIS short = IBIS,
long = Interferometric Bidimensional Spectropolarimeter, cite =
Cavallini:2006,Reardon:2008 FIRS short = FIRS, long = Facility Infrared
Spectropolarimeter, cite = Jaeggli:2011 EIT short = EIT, long = Extreme
Ultraviolet Imaging Telescope, cite = Delaboudiniere:1995 ROSA short = ROSA,
long = Rapid Oscillations in the Solar Atmosphere, cite = Jess:2010 SPINOR
short = SPINOR, long = Spectro-Polarimeter for Infrared and Optical Regions,
cite = Socas:2006 SOHO short = SOHO, long = Solar and Heliospheric
Observatory, cite = Domingo:1995 MDI short = MDI, long = Michelson Doppler
Imager, cite = Scherrer:1995 NMSU short = NMSU, long = New Mexico State
University NSO short = NSO, long = National Solar Observatory TAC short = TAC,
long = time allocation committee DEM short = DEM, short-plural-form = DEMs,
long = Differential Emission Measure, long-plural-form = Differential Emission
Measures PI short = PI, short-plural-form = PI, long = Principle Investigator,
long-plural-form = Principle Investigators IDL short = IDL, long = Interactive
Data Language bb short = bb, long = broadband nb short = nb, long = narrowband
IPM short = IPM, long = interplanetary medium AO short = AO, long = adaptive
optics, cite = Rimmele:2004 IR short = IR, long = infrared halpha short =
H-$\alpha$, long = Hydrogen-$\alpha$ FPI short = FPI, long = Fabry-Pérot DWDM
short = DWDM, long = dense wavelength division multiplexing CCD short = CCD,
long = charge-coupled device HI short = HI, long = Heliospheric Investigation
ROI short = ROI, long = region-of-interest POV short = POV, long = point-of-
view MHS short = MHS, long = magnetohydrostatic MHD short = MHD, long =
magnetohydrodynamic RMHD short = RMHD, long = radiative magnetohydrodynamic
DOT short = DOT, long = Dutch Open Telescope, cite = Rutten:1997 BBSO short =
BBSO, long = Big Bear Solar Observatory NST short = NST, long = New Solar
Telescope, cite = Goode:2012 SOT short = SOT, long = Solar Optical Telescope,
cite = Tsuneta:2008 hinode short = Hinode, long = Hinode, cite = Kosugi:2007
GST short = GST, long = Goode Solar Telescope SST short = SST, long = Swedish
1-m Solar Telescope, cite = Scharmer:2002 TRACE short = TRACE, long =
Transition Region and Coronal Explorer, cite = Handy:1999 PCTR short = PCTR,
long = prominence-corona-transition-region DKIST short = DKIST, long = Daniel
K. Inouye Solar Telescope RTI short = RTI, long = Rayleigh-Taylor instability
SSW short = SSW, long = SolarSoftWare, cite = Freeland:1998 LTE short = LTE,
long = local thermodynamic equilibrium NLTE short = NLTE, long = non-”local
thermodynamic equilibrium” FTS short = FTS, long = Fourier Transform
Spectrometer, cite = Kurucz:1984 RTE short = RTE, long = radiative transfer
equation CRD short = CRD, long = complete frequency redistribution PRD short =
PRD, long = partial frequency redistribution BCM short = BCM, long = Beckers’
cloud model, cite=Beckers:1964 HSRA short = HSRA, long = Harvard Smithsonian
Reference Atmosphere, cite=Gingerich:1971 NICOLE short = NICOLE, long = Non-
LTE Inversion COde using the Lorien Engine, cite=Socasnavarro:2015 SIR short =
SIR, long = Stokes Inversion based on Response functions, cite=Ruizcobo:1992
HAZEL short = HAZEL, long = Hanle and Zeeman Light, cite=AsensioRamos:2008
BPSS short = BPSS, long = bald patch separatrix surface, cite=Bungey:1996 EIS
short = EIS, long = EUV Imaging Spectrometer, cite=Culhane:2007 FWHM short =
FWHM, long = full width at half maximum AMRVAC short = MPI-AMRVAC, long =
Adaptive Mesh Refinement Versatile Advection Code, cite =
Keppens:2012,Porth:2014,Xia:2018,Keppens:2020 AMR short = AMR, long = adaptive
mesh refinement CCI short = CCI, long = Convective Continuum Instability BV
short = BV, long = Brunt-Väisälä TVDLF short = TVDLF, long = Total Variation
Diminishing Lax-Friedrich TI short = TI, long = Thermal Instability TNE short
= TNE, long = thermal non-equilibrium MALI short = MALI, long = Multilevel
Accelerated Lambda Iteration yt short = yt, long = yt-project, cite =
Turk:2011 GL short = GL, long = Gauss-Legendre CM short = CM, long =
classically mounted MULTI3D short = MULTI3D, long = multi-level non-LTE 3D,
cite = Leenaarts:2009 EB short = EB, long = Eddington-Barbier KDE short = KDE,
long = kernel density estimate NLFFF short = NLFFF, long = nonlinear force-
free field

11institutetext: Instituto de Astrofísica de Canarias, 38205 La Laguna,
Tenerife, Spain 22institutetext: Universidad de La Laguna, 38206 La Laguna,
Tenerife, Spain 33institutetext: Centre for mathematical Plasma-Astrophysics,
Celestijnenlaan 200B, 3001 Leuven, KU Leuven, Belgium 44institutetext:
Rosseland Centre for Solar Physics, University of Oslo, Postboks 1029
Blindern, 0315 Oslo, Norway 55institutetext: Institute of Theoretical
Astrophysics, University of Oslo, Postboks 1029 Blindern, 0315 Oslo, Norway
∗55email<EMAIL_ADDRESS><EMAIL_ADDRESS>
# From eruption to post-flare rain: a 2.5D MHD model

Samrat Sen∗ 11 2 2 3 3 Avijeet Prasad 44 5 5 Valeriia Liakh 33 Rony Keppens 33

(Received: XXXX; accepted: XXXX)

###### Abstract

Context. Erupting magnetic flux ropes (MFRs) play an important role in
producing solar flares, whereas fine-scale, condensed coronal rain is often
found in post-flare loops. However, the formation of the MFRs in the pre-flare
stage, and how this leads to coronal rain in a post-eruption magnetic loop is
not fully understood.

Aims. We explore the formation, and eruption of MFRs, followed by the
appearance of coronal rain in the post-flare loops, to understand the magnetic
and thermodynamic properties of eruptive events and their multi-thermal
aspects in the solar atmosphere.

Methods. We perform a resistive-magnetohydrodynamic (MHD) simulation with the
open-source code MPI-AMRVAC to explore the evolution of sheared magnetic
arcades that can lead to flux rope eruptions. The system is in mechanical
imbalance at the initial state, and evolves self-consistently in a non-
adiabatic atmosphere under the influence of radiative losses, thermal
conduction, and background heating. We use an additional level of adaptive
mesh refinement to achieve the smallest cell size of $\approx 32.7$ km in each
direction to reveal the fine structures in the system.

Results. The system achieves a semi-equilibrium state after a short transient
evolution from its initial mechanically imbalanced condition. A series of
erupting MFRs is formed due to spontaneous magnetic reconnection, across
current sheets created underneath the erupting flux ropes. Gradual development
of thermal imbalance is noticed at a loop top in the post-eruption phase,
which leads to catastrophic cooling and formation of condensations. We obtain
plasma blobs which fall down along the magnetic loop in the form of coronal
rain. The dynamical and thermodynamic properties of these cool-condensations
are in good agreement with observations of post-flare coronal rain.

Conclusions. The presented simulation supports the development and eruption of
multiple MFRs, and the formation of coronal rain in post-flare loops, which is
one of the key aspects to reveal the coronal heating mystery in the solar
atmosphere.

###### Key Words.:

Magnetic reconnection – Magnetohydrodynamics (MHD) – Sun: corona – Sun:
filaments, prominences – Sun: flares

##  1 Introduction

Eruptive events such as flares and coronal mass ejections (CMEs) occur due to
the conversion of the stored magnetic energy of twisted magnetic field lines
into the abrupt release of kinetic, thermal, and non-thermal energies
(Antiochos et al., 1999; Fleishman et al., 2020). Astrophysical plasmas, in
general, have very high Lundquist numbers, $S_{L}=v_{a}L/\eta_{m}$ (where,
$v_{a}$, $L$, and $\eta_{m}$ are the Alfvén velocity, length scale, and
magnetic diffusivity, respectively), and particularly for the solar corona,
$S_{L}\sim 10^{10}$ (Aschwanden, 2004). Under such conditions, Alfvén’s flux
freezing theorem (Alfvén, 1942) is valid throughout most of the coronal
volume, except in localized current sheets and in regions with strong magnetic
gradients. Magnetic structures associated with the eruptive events can
sometimes be described as magnetic flux ropes (MFRs), which are a bundle of
twisted magnetic field lines winding about a common axis (Chen, 2011; Priest,
2014; He et al., 2022). These MFRs can rise under favourable conditions
through ideal magnetohydrodynamic (MHD) instabilities such as kink and torus
instabilities (Fan & Gibson, 2007; Kliem & Török, 2006), or through magnetic
reconnection, where the magnetic field lines (MFLs) change their topology when
they come close to each other in a resistive medium. The reconnection process
converts the magnetic energy into thermal and kinetic energies and can
generate acceleration of the charged particles below the flux ropes (Amari et
al., 2003; Aulanier et al., 2012). A large fraction of this released energy
can be transported downward from the corona to the denser chromosphere through
thermal conduction and energetic electrons (Antiochos & Sturrock, 1978; Ruan
et al., 2020; Druett et al., 2023) in an arcade-like magnetic loop system.
This energy deposition can lead to sudden heating of the local plasma in the
chromosphere and causes upward evaporation into the corona to increase the
density and temperature of the magnetic loop system (Druett et al., 2023). The
loops gradually cool down to their usual coronal condition mainly through
radiative cooling and field-aligned thermal conduction (Cargill et al., 1995;
Aschwanden & Alexander, 2001). These loops can be visible in extreme
ultraviolet (EUV) (Scullion et al., 2014, 2016; Mason et al., 2019, and
references therein) and also in H$\alpha$ (Bruzek, 1964) bands, and are called
post-flare loops (PFLs).

One of the most fascinating plasma processes in the solar corona is its multi-
thermal behavior, where cool ($\sim 10^{4}$ K), and dense ($\sim 10^{10}$
cm-3) materials appear in a background hot ($\sim 1$ MK) medium. This
phenomenon mainly manifests itself in two distinct forms:
prominence/filaments, and coronal rain, which are basically end results of the
same underlying physics. Radiative processes are always at the origin of
coronal rain or coronal condensation in general. When there is an abrupt
cutoff in heating, which can happen after flares, the radiative cooling in a
localized domain can dominate over heating. Temperature in that domain can
drop down to $\sim 10^{4}$ K forming cool-condensed plasma that falls down
towards the chromosphere in the form of rain (Antiochos, 1980). Coronal rain
has been observed in PFLs in the classic “loop prominence model” by Bruzek
(1964). It is believed that condensations form in a PFL as a consequence of
flare heating (Antiochos, 1980), and would lead to a single appearance and
subsequent fall of coronal rain. This is supported by the observations of
Scullion et al. (2014, 2016) and Liu et al. (2015). A related model studies
loops undergoing “thermal non-equilibrium” (TNE) as proposed by Antiochos &
Klimchuk (1991) as the origin of long-lived prominence condensation. The idea
is when the heating in the coronal loops is sufficiently localized at the loop
foot points, the materials at the loop top cannot achieve thermal equilibrium.
Then the loop materials can undergo a quasi-cyclic evolution of heating and
chromospheric evaporation followed by catastrophic radiative cooling to form
coronal rains. Active region loops can undergo this process to form many
cycles of coronal rain, which is also supported by the observations of Antolin
et al. (2015); Froment et al. (2015) and Froment et al. (2017). More recently,
many 2D and 3D models have emerged (Claes & Keppens, 2019; Claes et al.,
2020), also of purely coronal volumes, where the TNE scenario is excluded by
construction, and where the local trigger due to thermal instability is
emphasized instead. Simulations of coronal rain in multidimensional arcade
setups (Fang et al., 2015; Xia et al., 2017; Li et al., 2022) show a fair
agreement with actual observations of rain in loops, and in essence combine
the quasi-cyclic rain appearance with the thermal instability trigger process.

Towards realizing the aspect of coronal condensation after flux rope
eruptions, we use a 2.5D resistive MHD simulation in a non-adiabatic
description including optically thin radiative loss, (field-aligned) thermal
conduction and steady (but spatially varying) background heating. To mimick
coronal loop structures, we use a bipolar magnetic field configuration of
sheared arcades, which are non-force-free and embedded in a stratified solar
atmosphere with gravity. The Lorentz force (LF) of the system from the initial
state drives a rapid evolution where one forms flux ropes from initial bipolar
ﬁeld lines which ascend continuously towards eruption, mediated via the common
process of spontaneous magnetic reconnections. The non-adiabatic effects of
the system are crucial to trigger the gradual development of thermal
instability, and the onset of the catastrophic cooling phase, which may lead
to the formation of coronal condensation in general, and coronal rain in
particular.

Recent studies of reconnection-driven coronal condensations in solar current
sheets in non-adiabatic regimes have been reported by Sen & Keppens (2022) in
2D geometry, and in a more complex topology of magnetic flux ropes in 3D
geometry (Sen et al., 2023). Those studies did not include gravitational
stratification, and highlighted the interplay between tearing and thermal
instabilities in coronal current sheet setups. The formation of coronal rain
in a post-eruption coronal loop has received less attention to date. A study
by Ruan et al. (2021) shows the formation of coronal rain demonstrating
multiple post-flare rain cycles in a post-flare loop. However, the eruption of
MFRs is not captured in their study. In this work, we include the self-
consistent development and eruption of a series of MFRs, followed by the
formation of coronal rain in a loop top that gradually evolves with time,
splits into different rain blobs, and falls downwards. We investigate the
magnetic and thermodynamic properties of this event in detail to understand
their formation mechanisms and the fine-scale structures to reveal the multi-
thermal aspect of the solar atmosphere.

The rest of the paper is organized as follows. In Section 2, we describe the
numerical model with the precise initial condition, the governing equations
with the numerical schemes used, and the boundary conditions. In Section 3, we
report the main results, and present a detailed analysis of the underlying
physics. Section 4 summarizes the key results of this study, discusses the
applicability of the model to the solar atmosphere, makes comparison with
existing observations, and addresses the novelty of the work. We conclude how
this work can be useful in a broader aspect to understand the solar atmosphere
in future studies.

##  2 Numerical setup

The MHD processes in the solar atmosphere are extremely dynamic and non-
linear. To explore the small-scale dynamics, magnetic, and thermodynamic
properties of eruptions, and associated condensations in the solar corona, we
construct a non-adiabatic and resistive-MHD simulation in a 2.5D cartesian
geometry using MPI-parallelised Adaptive Mesh Refinement Versatile Advection
code (MPI-AMRVAC) 111Open source at: http://amrvac.org/ (Keppens et al., 2012;
Porth et al., 2014; Xia et al., 2018; Keppens et al., 2021, 2023).

###  2.1 Model setup and governing equations

To construct a magnetic field topology similar to sheared coronal loops, we
use an initial magnetic field configuration motivated by Priest (2014) and
Kumar et al. (2016),

$\displaystyle B_{x}$ $\displaystyle=B_{0}\ \text{sin}(k_{x}(x-x_{0}))\
\text{exp}(-yk_{x}/\sigma)\,,$ (1) $\displaystyle B_{y}$
$\displaystyle=B_{0}\sigma\ \text{cos}(k_{x}(x-x_{0}))\
\text{exp}(-yk_{x}/\sigma)\,,$ (2) $\displaystyle B_{z}$ $\displaystyle=B_{0}\
\text{tan}(\gamma)\ \text{sin}(k_{x}(x-x_{0}))\ \text{exp}(-yk_{x}/\sigma)\,,$
(3)

where, the spatial domain of the simulation extends from $x=0$ to $2\pi$ Mm,
along the horizontal direction, and $y=0$ to 25 Mm along the vertical
direction, whereas the $z$-direction has translational symmetry. The
parameter, $k_{x}=\frac{2\pi}{L_{x}}$ (where, $L_{x}=2\pi$ Mm is the
horizontal span of the simulation domain), which gives $k_{x}=1$ Mm-1. The
arcade is centered at the numerical domain by setting $x_{0}=\pi/2$ Mm. The
shearing angle of the projected field lines at the $y=0$ plane with the $+x$
axis is determined by

$\displaystyle\gamma=\text{tan}^{-1}(B_{z}/B_{x})\,.$ (4)

Note that, for $\gamma=0$, the guide field $B_{z}$ vanishes, and the field
lines are tangential to the constant $z$ plane which leads to no shearing. The
associated Lorentz force is given by,

$\displaystyle{\bf F_{L}}=$ $\displaystyle\
(\sigma^{2}-\text{sec}^{2}\gamma)B_{0}^{2}k_{x}\text{sin}\big{(}k_{x}(x-x_{0})\big{)}\
\text{exp}\bigg{(}-\frac{2k_{x}y}{\sigma}\bigg{)}$ (5)
$\displaystyle\times\bigg{(}\text{cos}\big{(}k_{x}(x-x_{0})\big{)}\
{\bf\hat{x}}-\frac{1}{\sigma}\text{sin}\big{(}k_{x}(x-x_{0})\big{)}\
{\bf\hat{y}}\bigg{)},$

where, ${\bf\hat{x}}$, and ${\bf\hat{y}}$ are the unit vectors along $x$, and
$y$ directions respectively. The advantage of choosing this magnetic field
configuration is that the system becomes force-free for
$\displaystyle{|\sigma|=\text{sec}(\gamma)}$, and non-force-free otherwise.
The initial magnetic field configuration (Equations (1)-(3)) obeys the
solenoidal condition, $\displaystyle{\nabla\cdot{\bf B}=0}$. In order to keep
this condition throughout the simulation, we use the parabolic diffusion
method (Keppens et al., 2003, 2023). To achieve a sufficient shearing in the
field lines at the initial stage, we set the value of $\gamma=72^{\circ}.5$.
Furthermore, we set $\sigma=10$, which ensures the presence of non-zero
Lorentz force at the initial state of the system. The parameter, $B_{0}$,
decides the strength of the magnetic field, where the maximum field strength
(at $y=0$), $B_{max}=B_{0}\sigma$. To set the magnetic field strength similar
to typical coronal loops, we use $B_{0}=9$ G, that gives $B_{max}=90$ G, which
is comparable with observation by Brooks et al. (2021), where they have
reported a field strength in active region loops in the range of $60-150$ G.
To obtain a typical solar coronal condition, we set the unit density,
$\bar{\rho}=2.34\times 10^{-15}$ g cm-3, and unit temperature,
$\bar{T}=10^{6}$ K.

We use as base resolution of the simulation 96 $\times$ 384 along the $x$ and
$y$ directions respectively, with only one additional adaptive mesh refinement
(AMR) level. This achieves the smallest cell size of $\approx$32.7 km in each
direction. The (de-)refinement is based on the errors estimated by the
gradients of the instantaneous density and of magnetic field components at
each time step (Lohner, 1987). We follow the time evolution of the system up
to 83.37 minutes.

To explore the evolution of the system, the following MHD equations are solved
numerically,

$\displaystyle\frac{\partial\rho}{\partial t}+\nabla\cdot(\rho\,{\bf v})=0,$
(6) $\displaystyle\frac{\partial(\rho\,{\bf v})}{\partial
t}+\nabla\cdot(\rho\,{\bf v}\otimes{\bf v}+p_{tot}\,{\bf I}-{\bf B}\otimes{\bf
B})=\rho\,{\bf g},$ (7) $\displaystyle\frac{\partial\mathcal{E}}{\partial
t}+\nabla\cdot\left[\mathcal{E}\,{\bf v}+p_{tot}\,{\bf v}-({\bf B}\cdot{\bf
v})\,{\bf B}\right]=\rho\,{\bf g}\cdot{\bf v}+\eta\,\textbf{J}^{2}-{\bf
B}\cdot\left[\nabla\times(\eta{\bf J})\right]$ (8) $\displaystyle\hskip
56.9055pt-\rho^{2}\Lambda(T)+H_{bgr}+\nabla\cdot\left[\kappa_{||}\,{\bf
b}\,({\bf b}\cdot\nabla)\,T\right],$ $\displaystyle\frac{\partial{\bf
B}}{\partial t}+\nabla\cdot({\bf v}\otimes{\bf B}-{\bf B}\otimes{\bf
v})+\nabla\times(\eta\,{\bf J})=\mathbf{0}\,,$ (9)
$\displaystyle\nabla\cdot\textbf{B}=0\,,$ (10)
$\displaystyle\textbf{J}=\nabla\times\textbf{B}\,.$ (11)

Here, we use magnetic units with unit magnetic permeability, I is the unit
tensor, and $\rho$, $T$, v, $\eta$ and ${\bf b}$ represent plasma density,
temperature, velocity vector, resistivity, and unit magnetic field vector,
respectively, and $\otimes$ is the symbol for the dyadic product of two
vectors. A uniform resistivity of $\eta=5\times 10^{-5}$ (in code unit), which
corresponds to $5.82\times 10^{10}$ cm2 s-1 is used throughout the simulation
domain. The associated Lundquist number, $S_{L}=Lv_{a}/\eta$ (where, $L$ is
the length scale and $v_{a}$ is the Alfvénic speed) $\sim 10^{5}$. We adopt
the Spitzer-type thermal conductivity, purely aligned along the magnetic
fields, ${\bf\kappa}_{||}=10^{-6}\,T^{5/2}$ erg cm-1 s-1 K-1. The total
pressure, $p_{tot}$ is the sum of the plasma and magnetic pressure,

$\displaystyle p_{tot}=p+\frac{B^{2}}{2},$ (12)

where $p$ is the plasma pressure determined by the thermodynamic quantities
through the ideal gas law of the system. The total energy density per unit
volume, ${\mathcal{E}}$, is obtained as the sum of the internal, kinetic and
magnetic energy densities,

$\displaystyle\mathcal{E}=\frac{p}{\gamma-1}+\frac{\rho
v^{2}}{2}+\frac{B^{2}}{2},$ (13)

where $\gamma=5/3$ is the ratio of the specific heats.

The density along the vertical domain is stratified according to hydrostatic
equilibrium with gravity where ${\bf g}=-g(y)\mathbf{e}_{y}$ and,

$\displaystyle g(y)=g_{0}\frac{R_{\odot}^{2}}{(R_{\odot}+y)^{2}},$ (14)

where, $g_{0}=274$ m s-2 is the gravitational acceleration at the solar
surface, and $R_{\odot}=695.7$ Mm is the radius of the Sun. We assume an
isothermal atmosphere, $T_{0}=1$ MK at the initial state of the system. The
isothermal condition is relevant as the vertical extension of the domain is
$\approx 25$ Mm. The initial density variation follows from

$\displaystyle\rho_{i}(y)=\rho_{0}\ \text{exp}\left[-y/H(y)\right],$ (15)

where, $\displaystyle{H(y)=\frac{\mathcal{R}T_{0}}{\mu
g_{0}}\frac{R_{\odot}+y}{R_{\odot}}}$ is the density scale height, which is
based on the gas constant $\mathcal{R}$, and mean molecular mass, $\mu$. The
scale height at the base ($y=0$) of our simulation domain is $\approx 50$ Mm,
which increases only up to $\approx 51.8$ Mm at the top boundary ($y=25$ Mm).
To use a typical coronal condition, we set $\rho_{0}=3.2\times 10^{-15}$ g
cm-3.

To achieve a coronal condensation phase through thermal runaway process, we
incorporate the non-adiabatic effects of optically thin radiation (fourth term
at the RHS of equation (8)), background heating (fifth term at the RHS of
equation (8)), and field-aligned thermal conduction (last term at the RHS of
equation (8)). To compensate for the initial radiative cooling exactly with
the background heating ($H_{bgr}$), we use

$\displaystyle H_{bgr}=\rho_{i}^{2}(y)\Lambda(T_{0})\,,$ (16)

where, $\Lambda(T)$ is estimated from the cooling function model developed by
Colgan et al. (2008) and Dalgarno & McCray (1972). This initial choice ensures
thermal balance between radiative cooling and background heat at the start,
but the deviations will occur as this height-dependent heating is kept
constant throughout the simulation. Thermal conduction is zero at the initial
state of the system due to our isothermal atmosphere assumption, but plays a
role in the posterior evolution away from isothermal conditions.

###  2.2 Numerical schemes and Boundary conditions

The equations (6-11) are solved numerically using a three-step, third-order
Runge-Kutta time integration method with the Van Leer flux limiter van Leer
(1974) and Total Variation Diminishing Lax-Friedrichs (TVDLF) flux scheme. We
use the ‘splitting’ of the magnetic field variable, ${\bf B}$, where it is
decomposed into a steady background term and a time-dependent deviation part
(Xia et al., 2018).

We use periodic boundary conditions at the side boundaries. Magnetic fields
are extrapolated at the top and bottom boundaries by third-order and second-
order zero-gradient methods respectively. The pressure and density values at
the bottom boundary are fixed according to their local initial values, as
prescribed in Jenkins & Keppens (2021). At the top boundary, they are set
according to the hydrostatic assumption, as prescribed in Zhao et al. (2017).
The velocity components are set anti-symmetric in the form of mirror
reflections at the bottom boundary to fill the ghost cells. At the top
boundary, we use third order extrapolation with zero-gradient condition for
all the velocity components, additionally with no inflows.

(a)

(b)

Figure 1: (a) Temporal evolution of spatial averaged absolute velocity
components $\langle|v_{x}(t)|\rangle$, $\langle|v_{y}(t)|\rangle$,
$\langle|v_{z}(t)|\rangle$ and ratio of free ($E_{f}$) to potential
($E_{pot}$) energy (dashed curve in magenta). The shaded region between
$t=5.44$ and 7 min (shown by two vertical dashed lines) represents the
dynamical semi-equilibrium phase of the system, where
$\langle|v_{y}(t)|\rangle\ \approx 19$ km s-1, $\langle|v_{x}(t)|\rangle\
\approx 2.2$ km s-1, $\langle|v_{z}(t)|\rangle\ \approx 5.4$ km s-1, and
$E_{f}/E_{pot}\approx 2.1$. (b) Evolution of the spatial average,
$\alpha=\langle\frac{|{\bf J}\times{\bf B}|}{|{\bf J}|\ |{\bf B}|}\rangle$.

(a)

(b)

(c)

Figure 2: Discrete probability densities of the characteristic time scales of
(a) background heating ($\tau_{bh}$), (b) radiative cooling ($\tau_{rc}$), and
(c) thermal conduction ($\tau_{con}$) at $t=5.44$ min.

Figure 3: Top row: Spatial distribution of plasma density for (a) semi-
equilibrium stage. (b)-(d) shows the different MFR eruption stages, where the
over-plotted solid curves represent the magnetic field lines projected in the
$x-y$ plane. Middle row: Same as the top panel for temperature. Bottom row:
Distribution of Lorentz force (LF) over the entire simulation domain, where
the over-plotted arrows represent the LF direction in the $x-y$ plane. The
saturation level of the colour bar is taken up to $10^{-5}$ G2 cm-1 for a
better visualization of the enhanced LF regions. Animations of the figures are
available online (height and width aspect ratio of the figures are not in
scale).

(a) $t=7.66$ min

(b) $t=12.67$ min

(c) $x=3.14$ Mm

(d) $x=3.14$ Mm

Figure 4: Variation of $B_{y}$ field in the $x-y$ plane is shown at (a)
$t=7.66$ and (b) 12.67 min, where the over-plotted black arrows represent the
velocity field $\\{v_{x},v_{y}\\}$ in the $x-y$ plane. The region inside the
green ellipses (in Figures (a) and (b)) contain the reconnection point
locations, where we see the plasma inflow and outflow due to reconnection
(height and width aspect ratio in the figures are not in scale). (c) $v_{y}$
and $B_{y}$ distributions along $y$-direction at $x=3.14$ Mm (which covers the
reconnection point location). (d) Absolute velocity of $v_{x}$ and $B_{x}$
distributions along $y$-direction at $x=3.14$ Mm for the eruption phases. The
velocities are scaled in units of instantaneous Alfvén speed ($v_{a}$). The
solid and dashed curves in Figures (c) and (d) represent the velocities and
magnetic field strengths respectively. Figure 5: Temporal evolution of
(normalized) mean kinetic and magnetic energy densities during the eruption
phases, where the black arrows represent the time when first (Erp 1) and
second (Erp 2) eruptions from the central region have the peak KE densities.
The other peaks between the two eruptions (Erp 1 and 2) show the signature due
to eruptions from the side boundaries.

(a)

(b)

Figure 6: Time-distance map of LF (in arbitrary unit) along the vertical cut
at $x=3.14$ Mm. (a) and (b) represent the time domain around the first and
second flux rope eruption phases, respectively. The slope(s) at the selected
part(s) in the maps are highlighted by the black lines to estimate the
eruption speed of the flux ropes.

(a)

(b)

(c)

Figure 7: Velocity map along the projected field lines on the vertical $x-y$
plane, $v_{fl}$, for three different times.

Figure 8: Top panel: The temperature evolution at four different horizontal
locations at $x=1,2,3.14$, and 4.7 Mm all at the same height $y=4$ Mm. The
shaded region between $t=7$ to 38 min represent the eruption phase, where
several MFRs leave the domain. The vertical dashed line at $t=42$ min marks
the time when the upper coronal part again settles down to around 1 MK. The
dashed line at $t=64$ min represents the onset of a thermal runaway process,
which happens first in the central $x=3.14$ Mm location. Bottom panel:
Temperature evolution in the central $(x,y)=(3.14,4)$ Mm location. The
temperature drops between $t=42$ and 50 min with a decay rate of 156 K s-1 as
represented by the red dashed line. Figure 9: The development of thermal
imbalance between heating and cooling at the loop top location
$(x,y)=(3.14,4)$ Mm in the late phase of our simulation. The red curve denotes
the heating rate due to thermal conduction ($H_{con}$) and background heat
($H_{bgr}$), and the blue curve denotes the cooling rate due to radiation
($H_{rc}$). The units for all the quantities are in erg cm-3 s-1. The vertical
dashed line marks the time after which the radiative cooling starts to
dominate over the heating due to conduction and background heat.

##  3 Results and analysis

###  3.1 Achieving semi-equilibrium phase of the system

Due to the initial non-force-free nature of the magnetic field configuration,
and the corresponding mechanical imbalance, the system is forced to evolve
away from this state. We follow the evolution of the spatial average of the
absolute velocity components,

$\displaystyle\langle v_{x,y,z}\rangle=\frac{1}{V}\int_{V}|v_{x,y,z}|\ {\rm
d}V$ (17)

over the entire simulation domain, $V$ as shown in Figure 1(a). The average
velocity components of the system take a jump due to the large Lorentz force
from the initial state ($t=0$). The maximum velocity jump of the $v_{y}$
component, which is $\approx 100$ km s-1, whereas the other two components are
smaller. Then the system relaxes achieving a minimum velocity
$\langle|v_{y}|\rangle\approx$ 19 km s-1 at $t=5.44$ min which remains almost
constant till $t=7$ min. The other components, $\langle|v_{x}|\rangle$, and
$\langle|v_{y}|\rangle$ also remain constant to $\approx 2.2$ and 5.4 km s-1
respectively within this time domain. This is marked by the shaded region in
Figure 1(a). To estimate the energetics of this rapid initial relaxation phase
of the evolution, we calculate the ratio of the free and potential field
energies at each time steps between $t=0$ to 15 min. The potential field
energy ($E_{pot}$) is calculated as

$\displaystyle E_{pot}=\int_{V}\frac{B_{pot}^{2}}{2}{\rm d}V,$ (18)

where $\@vec{B}_{pot}$ is the potential field extrapolated using the
instantaneous $B_{y}$ component at the base of the simulation domain ($y=0$).
The free energy $E_{f}$ is the excess magnetic energy of the system from its
corresponding potential field energy, given by

$\displaystyle E_{f}=\int_{V}\frac{B^{2}}{2}{\rm d}V-E_{pot}\,.$ (19)

We notice from Figure 1(a) that the value of $E_{f}/E_{pot}$ decreases sharply
from 4.38 at the initial state to about 2.1 within $\approx 1$ min, which
remains constant till 7 min. The reason why we see an evolution of the ratio
$E_{f}/E_{pot}$ towards the value two can be argued from changes in the
magnetic energy observed in active regions during major flares (Gupta et al.,
2021; Liu et al., 2023). To estimate the evolution of the Lorentz force, we
calculate the spatial average of the parameter,

$\displaystyle\alpha=\frac{1}{V}\int_{V}\frac{|{\bf J}\times{\bf B}|}{|{\bf
J}|\ |{\bf B}|}\ {\rm d}V,$ (20)

within the full simulation domain at each time step as shown in Figure 1(b).
This parameter represents the average relative angle formed by the current and
the magnetic field which is useful for determining the (non-)force-free nature
of the magnetic field. The initial value of the average angle between ${\bf
J}$ and $\bf{B}$ is $\arcsin(\alpha)\approx 18^{\circ}$ which reduces to
$\approx 5^{\circ}$ at $t\approx 5.44$ min, and has a minimum value of
$\approx 4^{\circ}$ at $t=7$ min. This implies that the system achieves a
nearly force-free state within minutes, and that the period between 5 and 7
minutes corresponds to a nearly semi-equilibrium phase. We will explain
further on how the changes after $t=7$ min are accompanied by large flux rope
expulsions and flare-loop configurations, which again cause large velocity and
Lorentz force variations. We first characterize the thermodynamical timescales
for the nearly-force-free state between $t=5.44$ and 7 min.

We estimate the characteristic time scales of background heating
($\tau_{bh}$), radiative cooling ($\tau_{rc}$), and field-aligned thermal
conduction ($\tau_{con}$) at the starting time of the dynamical semi-
equilibrium state at $t=5.44$ min, which are calculated by

$\displaystyle\tau_{bh}$ $\displaystyle=\frac{\epsilon}{H_{bgr}}\,,$ (21)
$\displaystyle\tau_{rc}$
$\displaystyle=\frac{\epsilon}{\rho^{2}\Lambda(T)}\,,$ (22)
$\displaystyle\tau_{con}$
$\displaystyle=\frac{\epsilon}{|\nabla\cdot(\kappa_{||}{\bf b}({\bf
b}\cdot\nabla)T)|}\,,$ (23)

respectively, where $\epsilon=p/(\gamma-1)$ is the internal energy density per
unit volume. We calculate these characteristic time scales at each pixel of
the simulation domain, and obtain probability densities for them, as shown in
Figure 2. We notice that $\tau_{bh}$ and $\tau_{rc}$ are more than 1000 s for
the entire simulation domain (see Figures 2(a) and 2(b)). For estimating the
(magnetic field aligned) thermal conduction time scales $\tau_{con}$, we use a
convolution filter using a uniform kernel of size $5\times 5$ to smooth out
numerical discretization bias that arises in specific regions. We notice that
$\approx 76\%$ of the domain corresponds to $\tau_{con}>100$ s, and only 1.05%
of the domain has less than 1 s. However, the dynamical semi-equilibrium phase
of the system lasts for $\approx 93$ s (the shaded region in Figure 1(a)).
Taking all these aspects into consideration, where cooling and heating times
are long and of similar magnitude, while thermal conduction quickly
equilibrates temperatures along loops, we argue that the semi-equilibrium
phase between 5.44 and 7 min is a relevant starting point to study the later
evolution of the system.

###  3.2 Development and eruption of flux ropes

The plasma density, temperature, Lorentz force, and magnetic field line
distribution of the semi-equilibrium state (at $t=5.44$ min) are shown in the
first column of Figure 3. The magnetic field lines resemble a coronal loop
structure embedded in a bipolar magnetic field configuration. The evolution
during the semi-equilibrium phase is slow in comparison to the violent initial
phases and the thermal processes are in no position to substantially modify
the plasma during that phase; yet it is not a pure equilibrium; rather, there
is a small imbalance in the forces which leads to the events described in the
following.

The LF tries to squeeze the magnetic arcades along the horizontal direction
near the central region at $x=3.14$ Mm (and also towards the periodic
horizontal boundaries). This is clearly noticeable in Figure 3(i) for the
upper $y\gtrapprox 7$ Mm region. The squeezing of opposite polarity field
lines creates a localized region with small (of order one Gauss) $B_{x}$ and
$B_{y}$ components in the presence of the guide field, $B_{z}$. This location
leads to the magnetic reconnection without null point due to the presence of
(uniform) resistivity (Priest & Démoulin, 1995), which can be noticed in
Figures 3(b) and 3(f). At around $\approx 7.66$ min, we see the formation of a
vertical current sheet (CS) in the central region between $y\approx 5-10$ Mm,
where the LF is enhanced, as shown in Figure 3(j). This CS region becomes hot
up to $\approx 4$ MK as shown in Figure 3(f). Due to reconnection in the CS
plane, an erupting magnetic flux rope (MFR) is formed in the central region
(see Figure 3(b) or 3(f)). Due to the periodic condition at the horizontal
boundaries, erupting MFRs are also formed later at $\approx 10$ min at the
side boundaries. This is shown in Figs. 3(c) and 3(g). Further spontaneous
reconnections lead to the formation of a second erupting flux rope in the
central region at a lower height between $y\approx 0-5$ Mm at around $\approx
12.67$ min (see Figs 3(d) and 3(h)), again accompanied by the formation of a
CS underneath it. The temperature of this CS rises to $\approx 7$ MK due to
resistive heating. Cusp like structures are formed beneath the vertical CSs as
shown in Figures 3(b), (d), (f), and (h). This is similar to the classic CSHKP
solar flare model (Carmichael, 1964; Sturrock, 1966; Hirayama, 1974; Kopp &
Pneuman, 1976), where hot CSs are formed below erupting MFRs, and cusp like
magnetic structures form underneath CSs.

We show the flow patterns and the $B_{y}$ distributions for the instants with
clear central MFR eruptions in Fig. 4(a)-4(b). The reconnection location in
the bipolar magnetic structure is highlighted by the green ellipses in Figures
4(a), and 4(b), where the reconnection-driven inflows and outflows are
evident. To appreciate the reconnection location, we quantify the distribution
of $v_{y}$ and $B_{y}$ shown in Figure 4(c), and $v_{x}$ and $B_{x}$ in Fig.
4(d) by taking a vertical cut at $x=3.14$ Mm at $t=7.66$ and $12.67$ min
respectively. The regions where the outflow velocities flip directions are
spatially co-located with the small (of order one Gauss) magnetic field
strength locations of $B_{x}$ and $B_{y}$. The reconnection points are located
in the central region ($x=3.14$ Mm) at $y\approx 7.5$ Mm for the first, and
$y\approx 3$ Mm for the second (central) eruptions. The maximum outflow
velocities remain sub-Alfvénic, where the Alfvén velocity ($v_{a}$) is
calculated based on the instantaneous local plasma density and magnetic field
strength at a fixed location in the vicinity of the inflow side near the
reconnection points. To estimate the reconnection rate at those locations we
consider the ratio of the magnitude of inflow velocity with respect to the
instantaneous Alfvén velocity, $\mathcal{M}=|v_{x}/v_{a}|$. The variation of
$|v_{x}/v_{a}|$ in the central region is shown in Figure 4(d), where we see
the reconnection rate is $\approx 0.05$ for the first eruption, and $\approx
0.15$ for the second eruption near the reconnection points. It implies that
the second eruption is more violent than the first one.

This MFR eruption phase, where successive flux ropes are shed from the domain
and escape through the top boundary, which is open to outflows, is
energetically quantified in Fig. 5. Temporal distributions of the mean kinetic
(KE) and magnetic (ME) energy densities, which are calculated by the spatial
average of the entire simulation domain,

KE $\displaystyle=\frac{1}{V}\int_{V}\frac{\rho v^{2}}{2}{\rm d}V,$ (24) ME
$\displaystyle=\frac{1}{V}\int_{V}\frac{B^{2}}{2}{\rm d}V,$ (25)

are shown in Figure 5. The KE exhibits strong peaks two of which occur at
around $\approx 8$ and 13 min and are associated with the first and second
flux rope eruptions (marked as ‘Erp 1’ and ‘Erp 2’ respectively) from the
central region. The other peaks (one strong and one relatively weaker) in
between the two central eruptions are due to the secondary eruptions occurring
from the side boundaries (the first one is shown in Figure 3(c); the second
one is not shown in the figure, but can be found in the associated animation
available online). During the eruptions, the magnetic energy density decreases
as free magnetic energy is converted into kinetic (and thermal) energy, and
thus we see the corresponding drop of the magnetic energy density. The drops
are sharper during the peak eruption times when the KE density has peaks as
shown in Figure 5. We see clear signatures in Figures 3(j), and 3(l) that the
enhanced LF regions are co-located with the MFR locations of Figures 3(b), and
3(d) respectively.

To explore the dynamics of the eruptions themselves, we produce time-distance
(TD) maps of LF along a vertical cut at $x=3.14$ Mm, as shown in Figure 6.
Note that the choice of LF instead of density (or temperature) is preferred in
this case as the LF maps show a clear signature of the MFRs location. For
analysing the first eruption, we select the temporal domain of the TD map
between $\approx 6.3$ to 9.33 min, as shown in Figure 6(a), where we track the
enhanced LF region in the appropriate spatial and temporal domain
corresponding to the flux rope eruption, and estimate the eruption speed by
calculating the slope (shown by the black solid line) of that region in the
$y-t$ plane. Similarly, for the second eruption (from the central region), we
select the temporal domain between $\approx 11.4$ to 14.33 min, as shown in
Figure 6(b), and track the enhanced LF regions corresponding to the eruption
to estimate the eruption speeds. The speed of the first eruption is estimated
to be $\approx 167$ km s-1. For the second eruption, the flux rope speed is
$\approx 441$ km s-1 at the early phase and then decelerates to $\approx 217$
km s-1. This is in fair agreement with observational evidence in solar flares,
where the average rising velocity of the flux rope is reported $\approx 224$
km s-1 (Yan et al., 2018). The density is decreased in the ambient medium at
around $y\approx 10-20$ Mm after the first eruption (see Figure 3(d)). Also,
the formation of the second flux rope occurs at a lower height with a higher
LF. These two effects make the eruption of the second MFR faster than the
first one. This is also reflected in Figure 5, where we see that the mean
kinetic energy peak for the second eruption (Erp 2) is more than the peak
value of the first eruption (Erp 1).

(a) $t=64.83$ min

(b) $t=80.10$ min

(c) $t=83.01$ min

(d) $t=64.83$ min

(e) $t=80.10$ min

(f) $t=83.01$ min

(g) $t=64.83$ min

(h) $t=80.10$ min

(i) $t=83.01$ min

Figure 10: Density (top panel), temperature (middle panel), and gas pressure
(bottom panel) in different stages when the cool ($\sim 10^{4}$ K), and
condensed ($\sim 10^{-13}$ g cm-3) coronal rain forms. The location of three
separated density blobs is marked by the white boxes in (b). Overplotted black
solid curves represent the magnetic field lines, and the arrows in the bottom
panels represent the velocity field, $\\{v_{x},v_{y}\\}$ in the $x-y$ plane.

(a) $t=64.83$ min

(b) $t=80.1$ min

(c) $t=83.01$ min

Figure 11: Optically thin radiative cooling (RC) for the same time snapshots
as in Figure 10.

(a)

(b)

Figure 12: (a) Time series of instantaneous maximal and minimal temperature,
and peak density. (b) represents the temporal variation of instantaneous
minimum temperature highlighting the catastrophic cooling phase, where the
maximum cooling rate is estimated by the slope $\approx-6200$ K s-1 marked by
the red line. Figure 13: Variation of coronal rain mass and its area with
time.

(a)

(b)

Figure 14: (a) shows the trajectories of the individual rain blobs on $x-y$
plane. (b) shows the instantaneous vertical velocities of the individual
falling rain clumps.

###  3.3 Development of thermal instability

As a consequence of post-eruption of the MFRs, the loop footpoints are heated
due to the combined effect of (field-aligned) thermal conduction, steady
background heat (which is based on the initial density and temperature), and
radiative cooling. This leads to the evaporation of plasma materials along the
field lines. To investigate the flows of plasma materials, we estimate the
velocity along the projection of the magnetic field lines on the vertical $XY$
plane, namely, $\displaystyle{v_{fl}=\text{sign}(v_{y})\
\frac{|v_{x}B_{x}+v_{y}B_{y}|}{\sqrt{B_{x}^{2}+B_{y}^{2}}}}$, as shown in
Figure 7. The eruptions with several MFRs that escape upwards are followed by
the flow of plasma materials along the field lines from the loop foot points
to the central loop top. These materials accumulate at the loop top resulting
in a gradual density enhancement. The radiative cooling rate in the optically
thin corona, $H_{rc}\sim\rho^{2}\Lambda(T)$, has derivatives
$\displaystyle{\frac{\partial\,\text{log}\,H_{rc}}{\partial\,\text{log}\,\rho}=2}$,
and
$\displaystyle{\frac{\partial\,\text{log}\,H_{rc}}{\partial\,\text{log}\,T}\approx-0.4}$
at $T\approx 1$ MK for the $\Lambda(T)$ cooling curve we use in this model.
The temperature at the loop top location $(x,y)=(3.14,4)$ Mm remains $\sim 1$
MK between $t\approx 20$ to $40$ min (see Figure 8, top panel); this implies
that the radiative cooling rate at $T\approx 1$ MK is more sensitive to the
density perturbations than to those of temperature. The gradual mass
accumulation at the loop top at $x=3.14$ and around $y=4$ Mm leads to an
extended phase where a thermal instability slowly sets in and causes a local
condensation in the post-flare loop top. We follow the temporal evolution of
temperature at four different locations at $x=1,2,3.14$, and $4.7$ Mm all
taken at the same height namely $y=4$ Mm, as shown in Figure 8 (top panel).
The temperature at those locations can reach values up to several million
degrees during the eruption phase. The maximum temperature at $(x,y)=(3.14,4)$
Mm rises up to $\approx 7.6$ MK at the CS location beneath the MFR during the
second central eruption. After the second MFR eruption from the central
region, there is a secondary eruption from the side boundaries (as can be seen
in the animations associated with Figure 3), which leaves the simulation
domain from the top boundary at $\approx 38$ min. Indeed, about half an hour
into our simulation, no more MFRs erupt and the magnetic topology is seen to
have changed into an arcade configuration with an overlying mostly horizontal
field. The temperature evolution at the different selected $x$-values for
$y=4$ Mm show minor thermal fluctuations around a typical coronal $1$ MK
degree for $t=25$ min onward (see top panel of Figure 8). These fluctuations
are post-eruption effects. These perturbations settle down at $\approx 42$
min. The temperature at the loop top location $(x,y)=(3.14,4)$ Mm gradually
starts falling down right after that (as shown by the green curve in the top
panel of Figure 8). This signals the gradual development of thermal
instability at the post-eruption loop top. The temperature drop occurs as the
cooling rate, which is a combination of optically-thin radiation and field-
aligned thermal conduction dominates over the (steady but spatially varying)
background heating given by equation (16). The temperature evolution at this
$(x,y)=(3.14,4)$ Mm location in between $t=42$ and $50$ min represents the
(nearly) linear decay phase of temperature as shown by the bottom panel of
Figure 8. The temperature decay rate within this period is estimated to be
$\approx 156$ K s-1 by the slope marked by the red dashed line. The steady
background heating rate at the loop top location is $\approx 4.33\times
10^{-3}$ erg cm-3 s-1. Figure 9 shows that the net heating rate due to
conduction and background heat dominates over the radiative cooling rate at
the loop top location at the earlier stage of post eruption (between $t\approx
20$ and $40$ min). Later, radiative cooling gradually starts to dominate over
the heating from around $t=40$ min onwards, and at $t=64$ min the radiative
cooling rate increases sharply, when a catastrophic thermal runaway process
starts. This leads to a pressure drop at the loop top location compared to its
surroundings (see Figure 10(g)), which drives pressure gradient flows into
that location, and leads to a rapid density enhancement and drop in
temperature (shown by the right most vertical dashed line in the top panel of
Fig. 8).

###  3.4 Evolution of coronal rain

The location of the first condensation seed that appears in the central loop
top is shown in Figure 10(a). What happens next is shown visually in Figure
10, which contains three rows, for density, temperature and pressure, and
three columns, for $t=64.83$ min, $t=80.10$ min and $t=83.01$ min. Figures
10(a) and 10(d), shows the maps of the density and temperature at $t=64.83$
min. The density is enhanced up to $\sim 10^{-13}$ g cm-3 which is around
$100$ times more than the initial density ($\sim 10^{-15}$ g cm-3), and the
temperature drops to $\sim 10$ kK, which is around $100$ times lower than the
initial temperature (1 MK). We see that this thermal runaway process initiates
only at the magnetic loop top. The dips in the temperature curves in Figure 8
(top panel) for both the $(x,y)=(2,4)$ and (3.14, 4) Mm locations (shown by
the blue and green curves respectively) show what happens in the temperature
evolution at the fixed points, and must be understood from the evolution as
shown in Fig. 10(a)-10(c). The condensation accumulates surrounding plasma,
which is driven by the pressure gradient flows due to the pressure drop in the
condensation site. The pressure is as shown in the bottom panels of Figure 10,
where we notice the inflow of matter from higher to lower gas pressure regions
towards the condensation sites. In Figure 10(b), we see that the loop top
condensation has split into three isolated density clumps (or ‘blobs’) as
marked by the white boxes. These clumps fall downward along the magnetic field
lines due to the combination of gravity and the Lorentz force, and they
finally merge together again as shown in Figure 10(c), at $t=83.01$ min. These
cool-condensed structures share similarities with thermodynamic properties of
coronal rain in the solar atmosphere as observed by (Müller et al., 2005;
Antolin et al., 2015; Scullion et al., 2016; Mason et al., 2019; Şahin et al.,
2023, and references therein), where the density and temperature of the rain
clumps are reported between $\sim 10^{-14}-10^{-13}$ g cm-3, and $\sim
10^{4}-10^{5}$ K respectively.

From Figures 10(d)-10(f), we notice that the cool ($\sim 10^{4}$ K) material
is spatially co-located with the condensed structures in Figures 10(a)-10(c).
The corresponding radiative cooling losses are shown in Figure 11. The maximum
heat loss rate due to (optically thin) radiation in Figure 11(a) is $\approx
9.06\times 10^{-1}$ erg cm-3 s-1 which dominates over the background heating
of $\approx 4.33\times 10^{-3}$ erg cm-3 s-1. The maximum cooling rate
increases at later times to $3.65$ and $5.47$ erg cm-3 s-1, as shown in
Figures 11(b) and 11(c), respectively. The regions with high cooling rates
correspond to the cool-condensed regions, as shown in the top and middle rows
of Figure 10. Our 2.5D assumption influences the radiative cooling appearance,
showing the condensations as dark structures with bright edges, akin to the
prominence-corona-transition-region (PCTR) transition obtained in other recent
2.5D coronal rain simulations (Li et al., 2022).

The instantaneous maximum and minimum temperature evolution within the entire
simulation domain are shown in Figure 12(a), covering the entire period from
eruptions to the late condensation stages. The minimum temperature
distribution has two strong dips between $\approx 10-15$ min. The first dip
appears when the plasma blobs which are denser than the surrounding medium,
are trapped inside secondary erupting flux ropes from the side boundaries,
which also correspond to lower temperature than the background medium (see
Figures 3(c), and 3(g)). The second dip at about 14 minutes appears due to a
trapped plasma blob inside the rope from the second central eruption (see
Figures 3(d), and 3(h)). The corresponding signatures are also evident in the
red curve in Figure 12(a), which quantifies the instantaneous maximum density.
The same figure shows later on a sharp rise in this peak density (from $\sim
10^{-15}-10^{-13}$ g cm-3), and drop in minimum temperature from $\sim$ MK to
$\sim 10$ kK at $\approx 64$ min. This signals the onset of the thermal
runaway process to form the cool condensation at the loop top location
$(x,y)=(3.14,4)$ Mm. This catastrophic drop in (minimal) temperature happens
between $\approx 64$ to 65 min, and the sudden temperature change during that
short time interval is highlighted in Figure 12(b). The maximum rate of change
of the temperature during this minute is $\approx 6.2$ kK s-1. This
catastrophic cooling rate is in good agreement with a 2.5D MHD simulation
(Ruan et al., 2021) of a post-flare coronal rain scenario, where they found
the average catastrophic cooling rate to be 9 kK s-1, and in fair agreement
with observations for post-flare arcades (Scullion et al., 2016), where it is
reported to be $22.7$ kK s-1. The estimation of the cooling rates sometimes
can be challenging in the observations, since it requires measurements of
different spectral lines with different instruments and their responses,
sampling different parts of the flare. Note that the details of the cooling
curve prescription, $\Lambda(T)$, will also influence the condensation rates,
as was shown in a pure coronal volume study which compared many popular
cooling curves in Hermans & Keppens (2021). The temporal evolution of the rain
mass and its area are shown in Figure 13. Here, we define the localized
condensations as coronal rain clumps when the density of those regions is more
than or equal to a density threshold value, which is 10 times the maximum
density of the initial state. To calculate the total rain mass at an
instantaneous time, we collect all the cells with densities above the
threshold, then multiply the local densities with the area of the individual
cells to obtain the mass per unit length of those individual cells, and
finally aggregate the quantity of all those cells. The evolution of the rain
mass and area increase sharply from the onset time of thermal runaway at
$\approx 64$ min. This is also consistent with Figure 12(a) since the
catastrophic temperature drop occurs at the same time. The rain mass and area
reach maximum values at $\approx 78.5$ min. At about this time, the condensed
region is separated into three individual density blobs (shown in Figure
10(b)). Thereafter, we see a decrease in rain mass and area. This is because
the blobs start to coalesce while falling down, finally forming a single blob
again (see Figure 10(c)). The decrease of the rain mass and area after 79
minutes is also influenced by the fact that the smaller blobs are more easily
evaporated again, as they radiate energy in their PCTR and that they attain
lower heights where the background heating is higher. To track the motion of
individual rain blobs, we use an automated method using binary masking from a
python based module, OpenCV222more information can be found at
https://opencv.org/ to locate the instantaneous centroid locations of the
density blobs. We start the tracking when it first separates into three blobs
(around $t\approx 78$ min), and continue until the end of the simulation when
they merge together upon reaching the bottom ($t\approx 83.37$ min). Our ‘blob
1’ first falls down at the bottom, followed by ‘blob 2’ coalescing with ‘blob
1’, and finally the ‘blob 3’ falls down and merges to form a single blob. The
trajectories of the individual blobs are shown in Figure 14(a), where the
falling blobs follow trajectories (nearly) along magnetic field lines (see the
animation associated with Figure 3). Instantaneous vertical velocities of the
individual blobs are calculated based on the instantaneous distance traversed
along the vertical direction ($y$-direction) by the blobs between two
consecutive simulation frames, taken at a temporal cadence of 4.29 s. The
evolution of the falling rain velocity is shown in Figure 14(b), where we
notice that the velocity of the rain clumps is $\approx 1-40$ km s-1. It
increases as they accelerate while falling downward due to the combined effect
of gravity, gas pressure gradients, and LF. This is clearly noticeable,
especially for ‘blob 2’ and ‘blob 3’ (in Figure 14(b)), where they obtain the
maximum velocities when they reach the bottom just before coalescing with
other. They are decelerated suddenly, which is due to our lower boundary
condition and pressure gradient. The velocity of the downward falling rain
clumps in our study is in good agreement with other simulation models.
Examples include Fang et al. (2015), where the maximum velocity of the blobs
in a rain shower are $\approx 60$ km s-1, with most of the rain blobs having
velocities $\approx 30$ km s-1, and Li et al. (2022), where the velocity
varies from a few to 120 km s-1. It is also in reasonable agreement with the
the observations of Antolin et al. (2015) and Mason et al. (2019), who
reported velocities of the rain clumps in the range $\approx 50-100$ km s-1.

##  4 Summary and Discussion

In this work, we use a 2.5D resistive MHD simulation incorporating non-
adiabatic effects of optically thin radiative cooling, field-aligned thermal
conduction, and steady (but spatially varying) background heating to explore a
series of flux rope eruptions, followed by gradual development of thermal
instability, ultimately leading to a catastrophic cooling phase in the solar
corona. The initial magnetic field configuration used in this simulation is
motivated from Kumar et al. (2016), where they study the magnetic and
dynamical properties of an erupting flux rope in an incompressible medium.
Here, we use the similar initial magnetic field prescription, but extend the
study using a gravitationally stratified solar atmosphere in a compressible
medium, and the multi-thermal behaviour in the solar corona.

The initial magnetic field is non-force-free, and its strength depends on
three parameters, $\sigma$, $\gamma$ (which is related to the shearing angle),
and the field strength $B_{0}$. These initial forces instantly induce flows
that drive the plasma evolution. We showed how this imbalance transiently
evolves to a more force-balanced, energetically settled semi-equilibrium stage
after 5.44 minutes, which we use as a reference point for the rest of the
simulation. The Lorentz force-driven flows cause convergence of plasma along
the central axis and trigger magnetic reconnection leading to flux rope
formation and eruptions. The shearing angle of the magnetic arcades is also
one of the important parameters for the formation of the flux ropes, and our
simulation shows clear successive MFR formations with a shearing angle of
$\gamma=72^{\circ}.5$. To investigate the effect of a lower $\gamma$, we also
ran the simulation up to $\approx 90$ minutes with $\gamma=25^{\circ}$ keeping
all the other parameters the same (not shown in our figures), and notice that
there is no formation (and eruption) of flux ropes. Successive MFR eruptions
are observed in a multiple flux rope system in a solar flare (Hou et al.,
2018), where two consecutive “double-decker flux rope” eruptions happen within
an interval of five minutes due to strong shearing motions. This is very
similar with our findings, where the interval between two (central) eruptions
in our study is $\approx 4.9$ minutes. Multiple MFR eruptions are also found
in other simulations, such as in three-arcade configurations by Ding et al.
(2006) in a resistive medium, and Soenen et al. (2009) in an ideal MHD regime.
These studies show homologous eruptions of multiple flux ropes from the same
sites, which agrees with our findings where the two MFRs are ejected from the
central arcade region. However, the main focus of our study is to investigate
the multiple eruption of flux ropes and the multi-thermal behaviour in a
localized domain in corona. This present setup is not motivated to portray a
global MHD model, where it fits in the larger picture of CME evolutions, e.g.
due to breakout scenario, and how the large-scale magnetic topology (at
further out radii, i.e. in the top part of our simulation domain) evolves in
solar corona. We believe that the use of periodic boundary condition at the
side be adopted meaningfully and in fact a standard procedure for a local-box
simulation. Examples include, Kumar et al. (2016), where they have studied the
eruption, and Robinson et al. (2022), where the formation of a pre-flare flux
rope is studied in periodic side boundaries. In our simulation, the magnetic
fields above the central arcade switch to a primary horizontal orientation
when the eruptive flux ropes expand laterally, and it is of course mediated by
the choice of periodic boundary condition at sides. We witness a flipping of
the horizontal magnetic field ($B_{x}$ component) at the top part of the
domain after each flux rope eruption. This is because the magnetic field (and
mainly its horizontal component as introduced by the flux rope formations) due
to each eruption dominates over the remnant magnetic field of the previous
eruption. However, we perform an additional simulation with a modified
boundary condition, where we use a three-arcade system, and by extending the
lateral and vertical domain of our simulation box (see Appendix A). Changing
into this setup also gives the similar evolution: multiple eruption of flux
ropes from an initial non-force-free setup, and the late stages of the post-
flare arcade system that shows the formation of coronal condensation in a
localized domain in corona. It does change (obviously) the way the eruptions
proceed at the top part of the domain, but qualitatively the evolution are
common in both cases simulated. We also notice that the central flux ropes
expand laterally while eruption, and create horizontal field overlying the
arcades. We admit that the details of the magnetic field topology in the top
part of our domain, which drastically changes as successive flux ropes form
and pass through, is not fully realistic in neither of the two simulated
scenarios. We understand the importance of a more realistic setup which needs
to be handled better in future, where e.g. the actual expansion radially is
allowed, and the neigbouring magnetic flux evolutions are incorporated, such
as when doing a full spherical model, and portray a large-scale (global)
magnetic field of solar corona. Nevertheless, our main findings (multiple
eruptions and a post-flare multi-thermal evolution) is not affected by this.
The temperature evolution at the CSs underneath the MFRs during the second
central eruptions at $(x,y)=(3.14,4)$ Mm goes up to $\approx 7$ MK (from an
initial 1 MK temperature). This high temperature would show clear brightenings
in the soft X-ray band, as relevant for flares.

The development of the thermal instability is noticed to start after eruption
of the flux ropes (see Figure 8). It first leads to a linear decay of
temperature followed by a later catastrophic temperature drop and formation of
a condensation that evolves, splits and falls down due to gravity. The cool
material basically falls along the magnetic field lines of the loop. The
thermodynamic behaviour of these cool-condensed structures is in good
agreement with coronal rain observed in the solar atmosphere. We notice the
formation of the condensation in the post-flare loop after $\sim 30$ min of
the last flux rope eruption. This timescale is in agreement with the
observation for a C8.2-class post-flare coronal rain, where the time
difference between the flare ribbon formation, and the first appearance of
coronal rain in the post-flare loops is $\approx 26$ min (Scullion et al.,
2016). The observational evidence for coronal rain formation at the loop top
regions that fall downwards along the magnetic field lines is similar to what
we find in our simulation. The condensation time scale is also in a reasonable
agreement with an $X2.1$ class post-flare coronal rain observation, where the
time gap between the peak phase of the flare, and the first appearance of the
condensation at the loop top is $\approx 50$ min (Brooks et al., 2024). This
indicates that our model is in close correspondence with actual post-flare
coronal rain scenarios. The clumpy structure of the plasma blobs in our
simulation is very much a feature of flare-driven rain, more so than in the
quiescent (active region loop) rain, because the densities are much larger in
the flare case than in quiescent rain. The density and temperature structuring
of the rain, forming a coronal sheath around the cooler and denser cores of
the rain clumps we found in our study is one of the unique aspects of the
observations reported in Scullion et al. (2016). They detected the source of
the formation of the rain first in the hotter, and progressively later in the
cooler Extreme-ultraviolet (EUV) channels.

The role of resistivity $\eta$ in our model is an important aspect in the
formation of MFRs and the CSs they involve, but also regulates current-driven
Ohmic heating in all later phases where coronal condensations can form. The
uniform resistivity we use in this model is $\eta=5\times 10^{-5}$ (in
dimensionless units). To explore the effect of enhanced resistivity in the
condensation process, we have also run the experiment with a higher value,
$\eta=5\times 10^{-4}$ (in dimensionless units), keeping all the other
parameters unchanged. We do not see any occurrence of a catastrophic cooling
phase in the system till $90$ min. A more detailed parametric survey of
$\gamma$, $\eta$, $\sigma$, and $B_{0}$ for this model may be interesting to
study in the future to understand the dependence of the eruption and
condensation mechanisms on those parameters.

The maximum temperature in the upper half of the simulation domain reaches
high values of up to $\approx 25$ MK (see Figure 12(a)). This is influenced by
the fact that the successive MFRs also tend to evacuate this upper region,
causing it to be heated more easily. The heating occurs due to our steady
background heat prescription, which is based on the initial density and
temperature values. The mass outflow from the top boundary due to the
eruptions, and the prescription of no inflow through it leads to a density
drop in that upper region, and therefore the background heat dominates over
the cooling, leading to a high temperature. The use of the background heat
prescription in this model is motivated by the fact that it is exactly
compensated by the radiative cooling term in the initial state (there is no
cooling due to thermal conduction at the initial stage because of our initial
isothermal condition). Future extensions of this work should consider the
detailed coupling to the lower chromosphere and photospheric regions. A
localized heating prescription in the loop foot points can also trigger
thermal non-equilibrium (TNE) to promote and regulate the appearence of
coronal condensations, as proposed by Antiochos & Klimchuk (1991). A recent
study by Li et al. (2022) demonstrates the formation of coronal rain showers
due to turbulent heating at the chromospheric foot points. Advancement of our
model by extending the spatial domain from the chromosphere to the corona and
incorporating localized heating at the loop foot points can be interesting to
carry out, which we plan to explore in future simulations. This may produce a
quasi-cyclic formation of a larger number of rain clumps, or rain showers.

To conclude, our current model demonstrates the ability to reproduce the
formation of post-flare coronal rain, where multiple flux rope eruptions
happen from the same sites in the form of homologous ejections. It also sets
direction for a better understanding of the the multi-thermal ($\sim$ 10 kK -
10 MK) behaviour in the solar corona, which is important to reveal the broader
aspect of coronal heating.

###### Acknowledgements.

Data visualization and analysis were performed using Visit and yt-project. The
authors thank the anonymous referee for the constructive suggestions, and also
express their thanks to F. Moreno-Insertis, Eamon Scullion, Eric Priest, and
Manuel Luna for the useful discussions and suggestions, that improved the
manuscript considerably. SS and AP acknowledge support by the European
Research Council through the Synergy Grant #810218 (“The Whole Sun”,
ERC-2018-SyG). SS and RK acknowledge support by the C1 project TRACESpace
funded by KU Leuven, the ERC grant agreement 833251 Prominent ERC-ADG 2018,
and FWO grant G0B4521N. RK acknowledges ISSI in Bern, Team project #545. AP
would like to acknowledge the support from the Research Council of Norway
through its Centres of Excellence scheme, project number 262622. The
computational resources and services used in this work were provided by the
VSC (Flemish Supercomputer Center), funded by the Research Foundation Flanders
(FWO) and the Flemish Government - department EWI; La Palma supercomputer,
provided by Instituto de Astrofísica de Canarias (IAC).

## References

  * Alfvén (1942) Alfvén, H. 1942, Nature, 150, 405
  * Amari et al. (2003) Amari, T., Luciani, J. F., Aly, J. J., Mikic, Z., & Linker, J. 2003, ApJ, 585, 1073
  * Antiochos (1980) Antiochos, S. K. 1980, ApJ, 241, 385
  * Antiochos et al. (1999) Antiochos, S. K., DeVore, C. R., & Klimchuk, J. A. 1999, ApJ, 510, 485
  * Antiochos & Klimchuk (1991) Antiochos, S. K. & Klimchuk, J. A. 1991, ApJ, 378, 372
  * Antiochos & Sturrock (1978) Antiochos, S. K. & Sturrock, P. A. 1978, ApJ, 220, 1137
  * Antolin et al. (2015) Antolin, P., Vissers, G., Pereira, T. M. D., Rouppe van der Voort, L., & Scullion, E. 2015, ApJ, 806, 81
  * Aschwanden (2004) Aschwanden, M. J. 2004, Physics of the Solar Corona. An Introduction
  * Aschwanden & Alexander (2001) Aschwanden, M. J. & Alexander, D. 2001, Sol. Phys., 204, 91
  * Aulanier et al. (2012) Aulanier, G., Janvier, M., & Schmieder, B. 2012, A&A, 543, A110
  * Brooks et al. (2024) Brooks, D. H., Reep, J. W., Ugarte-Urra, I., Unverferth, J. E., & Warren, H. P. 2024, arXiv e-prints, arXiv:2401.04537
  * Brooks et al. (2021) Brooks, D. H., Warren, H. P., & Landi, E. 2021, ApJ, 915, L24
  * Bruzek (1964) Bruzek, A. 1964, ApJ, 140, 746
  * Cargill et al. (1995) Cargill, P. J., Mariska, J. T., & Antiochos, S. K. 1995, ApJ, 439, 1034
  * Carmichael (1964) Carmichael, H. 1964, in NASA Special Publication, Vol. 50, 451
  * Chen (2011) Chen, P. F. 2011, Living Reviews in Solar Physics, 8, 1
  * Claes & Keppens (2019) Claes, N. & Keppens, R. 2019, A&A, 624, A96
  * Claes et al. (2020) Claes, N., Keppens, R., & Xia, C. 2020, A&A, 636, A112
  * Colgan et al. (2008) Colgan, J., Abdallah, J., J., Sherrill, M. E., et al. 2008, ApJ, 689, 585
  * Şahin et al. (2023) Şahin, S., Antolin, P., Froment, C., & Schad, T. A. 2023, ApJ, 950, 171
  * Dalgarno & McCray (1972) Dalgarno, A. & McCray, R. A. 1972, ARA&A, 10, 375
  * Ding et al. (2006) Ding, J. Y., Hu, Y. Q., & Wang, J. X. 2006, Sol. Phys., 235, 223
  * Druett et al. (2023) Druett, M. K., Ruan, W., & Keppens, R. 2023, Sol. Phys., 298, 134
  * Fan & Gibson (2007) Fan, Y. & Gibson, S. E. 2007, ApJ, 668, 1232
  * Fang et al. (2015) Fang, X., Xia, C., Keppens, R., & Van Doorsselaere, T. 2015, ApJ, 807, 142
  * Fleishman et al. (2020) Fleishman, G. D., Gary, D. E., Chen, B., et al. 2020, Science, 367, 278
  * Froment et al. (2017) Froment, C., Auchère, F., Aulanier, G., et al. 2017, ApJ, 835, 272
  * Froment et al. (2015) Froment, C., Auchère, F., Bocchialini, K., et al. 2015, ApJ, 807, 158
  * Gupta et al. (2021) Gupta, M., Thalmann, J. K., & Veronig, A. M. 2021, A&A, 653, A69
  * He et al. (2022) He, W., Hu, Q., Jiang, C., Qiu, J., & Prasad, A. 2022, ApJ, 934, 103
  * Hermans & Keppens (2021) Hermans, J. & Keppens, R. 2021, A&A, 655, A36
  * Hirayama (1974) Hirayama, T. 1974, Sol. Phys., 34, 323
  * Hou et al. (2018) Hou, Y. J., Zhang, J., Li, T., Yang, S. H., & Li, X. H. 2018, A&A, 619, A100
  * Jenkins & Keppens (2021) Jenkins, J. M. & Keppens, R. 2021, A&A, 646, A134
  * Keppens et al. (2012) Keppens, R., Meliani, Z., van Marle, A. J., et al. 2012, Journal of Computational Physics, 231, 718
  * Keppens et al. (2003) Keppens, R., Nool, M., Tóth, G., & Goedbloed, J. P. 2003, Computer Physics Communications, 153, 317
  * Keppens et al. (2023) Keppens, R., Popescu Braileanu, B., Zhou, Y., et al. 2023, A&A, 673, A66
  * Keppens et al. (2021) Keppens, R., Teunissen, J., Xia, C., & Porth, O. 2021, Computers & Mathematics with Applications, 81, 316, development and Application of Open-source Software for Problems with Numerical PDEs
  * Kliem & Török (2006) Kliem, B. & Török, T. 2006, Phys. Rev. Lett., 96, 255002
  * Kopp & Pneuman (1976) Kopp, R. A. & Pneuman, G. W. 1976, Sol. Phys., 50, 85
  * Kumar et al. (2016) Kumar, S., Bhattacharyya, R., Joshi, B., & Smolarkiewicz, P. K. 2016, ApJ, 830, 80
  * Li et al. (2022) Li, X., Keppens, R., & Zhou, Y. 2022, ApJ, 926, 216
  * Liu et al. (2015) Liu, W., De Pontieu, B., Vial, J.-C., et al. 2015, ApJ, 803, 85
  * Liu et al. (2023) Liu, Y., Welsch, B. T., Valori, G., et al. 2023, ApJ, 942, 27
  * Lohner (1987) Lohner, R. 1987, Computer Methods in Applied Mechanics and Engineering, 61, 323
  * Mason et al. (2019) Mason, E. I., Antiochos, S. K., & Viall, N. M. 2019, ApJ, 874, L33
  * Müller et al. (2005) Müller, D. A. N., de Groof, A., de Pontieu, B., & Hansteen, V. H. 2005, in ESA Special Publication, Vol. 596, Chromospheric and Coronal Magnetic Fields, ed. D. E. Innes, A. Lagg, & S. A. Solanki, 37.1
  * Porth et al. (2014) Porth, O., Xia, C., Hendrix, T., Moschou, S. P., & Keppens, R. 2014, ApJS, 214, 4
  * Priest (2014) Priest, E. 2014, Magnetohydrodynamics of the Sun
  * Priest & Démoulin (1995) Priest, E. R. & Démoulin, P. 1995, J. Geophys. Res., 100, 23443
  * Robinson et al. (2022) Robinson, R. A., Carlsson, M., & Aulanier, G. 2022, A&A, 668, A177
  * Ruan et al. (2020) Ruan, W., Xia, C., & Keppens, R. 2020, ApJ, 896, 97
  * Ruan et al. (2021) Ruan, W., Zhou, Y., & Keppens, R. 2021, ApJ, 920, L15
  * Scullion et al. (2016) Scullion, E., Rouppe van der Voort, L., Antolin, P., et al. 2016, ApJ, 833, 184
  * Scullion et al. (2014) Scullion, E., Rouppe van der Voort, L., Wedemeyer, S., & Antolin, P. 2014, ApJ, 797, 36
  * Sen et al. (2023) Sen, S., Jenkins, J., & Keppens, R. 2023, A&A, 678, A132
  * Sen & Keppens (2022) Sen, S. & Keppens, R. 2022, A&A, 666, A28
  * Soenen et al. (2009) Soenen, A., Zuccarello, F. P., Jacobs, C., et al. 2009, A&A, 501, 1123
  * Sturrock (1966) Sturrock, P. A. 1966, Nature, 211, 695
  * van Leer (1974) van Leer, B. 1974, Journal of Computational Physics, 14, 361
  * Xia et al. (2017) Xia, C., Keppens, R., & Fang, X. 2017, A&A, 603, A42
  * Xia et al. (2018) Xia, C., Teunissen, J., El Mellah, I., Chané, E., & Keppens, R. 2018, ApJS, 234, 30
  * Yan et al. (2018) Yan, X. L., Yang, L. H., Xue, Z. K., et al. 2018, ApJ, 853, L18
  * Zhao et al. (2017) Zhao, X., Xia, C., Keppens, R., & Gan, W. 2017, ApJ, 841, 106

##  Appendix A Evolution of the system with modified setup

Figure 15: Spatial distribution of plasma density (top panel), and temperature
(bottom panel) for $t=13.23$ min, $21.11$ min, and $33.5$ min. The over-
plotted solid curves represent the magnetic field lines projected in the $x-y$
plane (height and width aspect ratio of the figures are not in scale). Figure
16: Spatial distribution of density (left) and temperature (right)
distribution at $t=35.78$ min, showing the splitting of the coronal rain blobs
that forms at the loop top. The overplotted black curves represent the
magnetic field lines projected in the $x-y$ plane.

We perform an additional simulation with a modified setup, where we use a
three-arcade system by extending the lateral domain of our boundary from $x=0$
to $3\pi$ Mm (9.42 Mm), and the top boundary from y=0 to 40 Mm. We use the
highest resolution of $65$ km and $54$ km along $x$ and $y$ directions
respectively. Here, we modify the side boundary condition from periodic to the
prescription according to Jenkins & Keppens (2021), where the density and
pressure are symmetric; $v_{x}$, $v_{y}$, $v_{z}$ are anti-symmetric,
symmetric, and anti-symmetric respectively; $B_{x}$, $B_{y}$, and $B_{z}$ are
anti-symmetric, symmetric, and anti-symmetric respectively at left and right
boundaries. The top and bottom boundary conditions follow the same
prescription as described in section 2. The density, temperature, and the
magnetic field line evolution for the eruptions are shown in Figure 15. We
notice the formation of multiple flux ropes from the central arcade which
expand laterally while eruption, and create horizontal field overlying the
arcades. The magnetic field topology above the arcades changes drastically as
successive flux ropes form and pass through. Finally, we notice the formation
of coronal condensation at the central loop top at $t\approx 33$ min, which
splits into different rain blobs and fall down, as shown in Figure 16. This
time scale of condensation is earlier than the other setup (which is around 64
min). This difference in time scale appears due to the initial magnetic field
structure with three arcades, modified boundary condition at the side,
different spatial domain, and spatial resolution.