Study of Reconnection Dynamics and Plasma Relaxation in MHD simulation of a Solar Flare

Self-organization in continuous systems is associated with dissipative processes. In particular, for magnetized plasmas, it is known as magnetic relaxation, where the magnetic energy is converted into heat and kinetic energy of flow through the process of magnetic reconnection. An example of such a system is the solar corona, where reconnection manifests as solar transients like flares and jets. Consequently, toward investigation of plasma relaxation in solar transients, we utilize a novel approach of data-constrained MHD simulation for an observed solar flare. The selected active region NOAA 12253 hosts a GOES M1.3 class flare. The investigation of extrapolated coronal magnetic field in conjunction with the spatiotemporal evolution of the flare reveals a hyperbolic flux tube (HFT), overlying the observed brightenings. MHD simulation is carried out with the EULAG-MHD numerical model to explore the corresponding reconnection dynamics. The overall simulation shows signatures of relaxation. For a detailed analysis, we consider three distinct sub-volumes. We analyze the magnetic field line dynamics along with time evolution of physically relevant quantities like magnetic energy, current density, twist, and gradients in magnetic field. In the terminal state, none of the sub-volumes are seen to reach a force-free state, thus remaining in non-equilibrium, suggesting the possibility of further relaxation. We conclude that the extent of relaxation depends on the efficacy and duration of reconnection, and hence, on the energetics and time span of the flare.


Introduction
Continuous dissipative systems exhibit self-organization (Hasegawa, 1985) by evolving preferentially toward states characterized by long-range ordering in one physical variable and short-range disorder in other variables.In the context of magnetized plasmas, self-organization is generally termed as plasma relaxation wherein the magnetofluid relaxes toward a minimum energy state while preserving appropriate physical variables.This process is quantitatively understood by forward cascade in one variable-which gets dissipated at the smallest scales in the system while, inverse cascade of other variables generates the long range order.A qualitative analysis of the phenomenon is also possible by relying on the principle of selective decay (Matthaeus and Montgomery, 1980).The principle effectively compares the decay rates of two or more variables in the presence of dissipation and the relaxed state is obtained by a constrained minimization of the fastest decaying variable while treating the slower decaying variables as invariants (Ortolani and Schnack, 1993).
The constrained minimization can be explored by briefly revisiting the Woltjer (1958) theory.There the relaxed state is obtained by minimizing the net magnetic energy of a magnetic flux tube while keeping its volume integrated magnetic helicity: a measure of magnetic topology, invariant.The relaxed state is characterized by a volume current density (J) to be entirely parallel to the magnetic field (B) such that the Lorentz force is zero (J × B = 0)-earning the moniker "Force-Free State".Incidentally, the force-free state also describes a magnetostatic equilibrium for a low β plasma, where the magnetic tension force is balanced by the magnetic pressure force.The proportionality factor between J and B, although constant for a given flux tube, can vary over different flux tubes and hence, is a function of position.The equations representing the Woltjer state are where α(r) represents the magnetic circulation per unit flux (Parker, 2012).The second equation is necessary to impose the solenoidality of B. Together, these two equations are called nonlinear force-free equations and are used to describe various physical systems including the solar corona.J. B. Taylor further conjectured that in an isolated and slightly resistive plasma, magnetic reconnection between different flux tubes will homogenize α(r) and plasma pressure (if any), resulting in a relaxed state obtained by a minimization of the global magnetic energy while keeping the global magnetic helicity invariant (Taylor, 1974(Taylor, , 1986(Taylor, , 2000)).Notably, A is the vector potential and the volume V encompasses the whole system domain.The relaxed state is described by the linear force-free field (Chandrasekhar and Kendall, 1957), satisfying with a constant α.Notably, relative to the global magnetic energy, the treatment of global magnetic helicity as an invariant is consistent with the selective decay principle-as demonstrated in Berger (1984), Browning (1988), Choudhuri (1998), and nicely summarized in Yeates (2020).Since the conservation of global magnetic helicity is central to Taylor relaxation, it merits further discussion.Notably, the helicity quantifies the interlinking and knotting between different magnetic field lines along with twisting and kinking of a given filed line (Berger and Field, 1984).Importantly, magnetic helicity is a pseudo-scalar because of its gauge dependency.For example, the definition in equation ( 4) is gauge independent only for an isolated system with no magnetic field line cutting across its boundaries.For open systems, the definition of magnetic helicity involves reference fields along with apt boundary condition to make it gauge invariant-see Yeates (2020) for details.
Both the Woltjer and Taylor relaxed states focus only on magnetic properties of the plasma and do not include other variables like the plasma flow, kinetic pressure, or dissipation rates.Inclusions of these variables are possible within the framework of two-fluid magnetohydrodynamics (MHD).Generally in the two fluid formalism, the minimizer incorporates plasma flow whereas the invariants are either the generalized ion and electron helicities (Steinhauer and Ishida, 1997) or their derivatives (Bhattacharyya and Janaki, 2004), and can even include total (magnetic +kinetic) energy (Yoshida and Mahajan, 2002).The obtained relaxed states are always "flow-coupled", i.e the plasma flow and magnetic field are interrelated.To avoid any confusion between the Woltjer or Taylor relaxed states with the two-fluid relaxed states, hereafter the former is called magnetic relaxation which, also emphasizes on its magnetic nature.Relevantly, Yeates, Hornig, and Wilmot-Smith (2010) proposed the existence of an additional constraint along with magnetic helicity, namely the topological degree of field line mapping (also see Yeates, Russell, andHornig, 2015, Yeates, Russell, andHornig, 2021) to obtain relaxed states.The relaxed states turned out to be either linear force-free or nonlinear force-free depending on the topological degree.SOLA: main.tex;23 January 2024;3:08;p. 3 One of the essential characteristics of magnetic relaxation is the dissipation of magnetic energy through reconnection.Relevantly, in solar coronal transients such as flares, jets, and coronal mass ejections, a fraction of the stored magnetic energy is dissipated in the form of heat and kinetic energy of charged particles through the process of magnetic reconnection (Shibata and Magara, 2011, Zweibel and Yamada, 2016, Li, Priest, and Guo, 2021).The reconnection assisted decay of magnetic energy along with the wealth of multiwavelength observations of these transients from space-based observatories, make them a suitable testbed to explore magnetic relaxation in nature.In this context, Nandy et al. (2003) analyzed several flare-productive active regions and found their time evolution to be tending towards a linear force-free state.A similar result was found by Murray, Bloomfield, and Gallagher (2013), who investigated the pre-flare and post-flare coronal magnetic fields in active region NOAA 10953.The authors determined the post-flare configuration as closer to linear force-free field and suggested this to be an indicator of incomplete Taylor relaxation.Recently, Liu et al. (2023) found some evidence for Taylor relaxation in increased homogenization of α for multiple X-class flares during 2010-2017.
Along with the observational studies, numerical simulations employing analytical magnetic fields have also been carried out.The simulation by Amari and Luciani (2000) employed bipolar potential fields driven by a 2D velocity field imposed at the bottom boundary.They found the terminal state of their simulation to be far from a constant-α field.Contrarily, Browning et al. (2008) and Hood, Browning, and van der Linden (2009) investigated the nanoflare heating model by following the development of kink instability in coronal loops.The relaxed state was found to be consistent with a linear force-free configuration.In the context of topological dissipation problem (Pontin and Hornig, 2020), Pontin et al. (2011) used braided magnetic fields and found the terminal state of relaxation to be nearly nonlinear force-free.For resistive MHD simulation of a solar coronal jet, Pariat et al. (2015) analyzed the evolution of helicity for several gauge choices, and found it to be approximately conserved.Recently, Robinson, Aulanier, and Carlsson (2023) explored the formation of a magnetic flux rope in MHD simulation of the Quiet Sun, where disordered low-lying coronal magnetic field lines undergo multiple small-scale reconnections.The authors recognized the process as self-organization, where an inverse cascade of helicity occurs, making the system to tend toward Taylor relaxation.
The above simulations although explore various elements of relaxation in coronal transients but are idealized scenarios because of their well-organized analytical initial conditions.In this paper a novel approach to study the magnetic relaxation is adopted where the initial field of a flaring region is obtained through extrapolation of coronal field using photospheric magnetograms.The idea is to accommodate the field line complexity of an actual flaring region which may not get captured in analytical magnetic fields.Subsequently, MHD simulation in combination with multiwavelength analyses is carried out to explore the reconnection dynamics and its consequence on the magnetic relaxation.
The paper is organized as follows.Section 2 describes the active region, the M1.3 class flare and the rationale behind selecting the AR.It also describes the temporal development of the transient activity using Extreme Ultraviolet SOLA: main.tex;23 January 2024; 3:08; p. 4 (EUV) observations from SDO/AIA satellite.In Section 3, we present the details and setup for magnetic field extrapolation, supplemented with morphological investigation and indices that characterize the accuracy of magnetic field reconstruction.In Section 4, the numerical framework for magnetohydrodynamics simulation is presented.Section 5 presents the results and analysis of this study and in Section 6, we present the summary and discussion of this work.

Active Region and Flare
The choice of active region (AR) and solar flare in this study is governed by the following rationale (a) The flare should be GOES (Geostationary Operational Environmental Satellite) M class or higher to ensure significant dissipation in magnetic energy (b) In the post-flare phase, there should not be any other major flaring activity so that the magnetic energy buildup and decay phases are sharp and clear (c) Due to the use of magnetic field extrapolation model in this study, the chosen AR should be nearly disk centered.Such a choice pertains to low measurement error in photospheric vector magnetic field and also minimizes the projection effects due to finite curvature of the photospheric surface (e.g.see Venkatakrishnan, Hagyard, and Hathaway, 1988).Both the effects combinedly reduce error during magnetic field extrapolation.Considering these constraints, we select the AR NOAA 12253, with heliographic coordinates as S05E01 on January 4, 2015.It hosts a GOES M1.3 class flare of net duration 35 minutes (min.), having start, peak, and end time as 15:18 UT, 15:36 UT, and 15:53 UT, respectively.Importantly, in the post-flare phase, there is no flaring activity for the next six hours.Along with the aforementioned criterion's, we make sure that the selected active region complies with the condition B z = const.at the bottom boundary, used in the MHD simulation.This translates into the requirement that during the course of flaring activity, the total relative change in magnetic flux (integrated over the bottom boundary) is minimal.We use the line-ofsight magnetograms from hmi.M 45 series of the Helioseismic Magnetic Imager (SDO/HMI: Schou et al., 2012, Scherrer et al., 2012) onboard the Solar Dynamics Observatory (SDO: Pesnell, Thompson, and Chamberlin, 2012), with temporal cadence of 45 seconds to evaluate this.The original magnetogram (Panel (a) in Figure 1) having dimensions of 4096×4096 in pixel units is CEA projected (Calabretta and Greisen, 2002) and cropped to match the pre-defined dimensions of HARP active region patch (Hoeksema et al., 2014) for AR NOAA 12253, which is 877×445 in pixel units (Panel (b) in Figure 1).Using this processed magnetogram, we find that over a period of 72 min., starting from 15:00 UT up to 16:12 UT, the relative changes in positive and negative flux with respect to their initial values are 0.36 % and 0.42 %, respectively.We explore the spatiotemporal evolution of flare using observations from 1600 Å and 304 Å channels of the Atmospheric Imaging Assembly (SDO/AIA: Lemen et al., 2012).In Figure 2, Panel (a) depicts the location of a brightening (labeled B 1 ) during the beginning of flare.The location is relevant because it might host a potential reconnection site and merits attention.The subsequent evolution reveals multiple brightenings during the flare peak (labeled B 2 ), as shown in Panel (b).Notably, in the observations of 304 Å channel, we find the presence of a dome shaped structure.Its spatial location is marked by the yellow colored box in Panel (c).A zoomed in view of this boxed region, with better image contrast, is given in Panel (e).The Panel highlights the approximate edges of the structure by yellow lines.These lines depict multiple connections between the central location C and the traced, nearly circular periphery.Furthermore, the line toward the west of C indicates the association of dome structure with magnetic morphology in rest of the active region.We note that the overlaid lines are in agreement with the expected two dimensional projection of a dome structure.Further developments reveal the complete spatial extent of flare dynamics where specific chromospheric flare ribbons are recognizable, as marked in Panel (d) by white arrows.Therefore, in observations, the brightenings B 1 , B 2 and the dome shaped structure are identified to be of significance and hence, merit investigation of associated magnetic field line morphologies for an understanding of their role in the flaring activity and relaxation process.

Magnetic Field Extrapolation
In order to explore the magnetic field line morphology of the active region, we employ a non-force-free extrapolation model (Bhattacharyya and Janaki, 2004, Bhattacharyya et al., 2007and Hu et al., 2010).The rationale for using the non-force-free extrapolation follows from an order of magnitude estimate for Lorentz force and rate of change of momentum on the photosphere (see Agarwal, Bhattacharyya, and Wiegelmann (2022) for details).The estimate is obtained SOLA: main.tex;23 January 2024; 3:08; p. 7 by using equation ( 6), as follows where all the quantities have usual meanings.Since, β ≈ 1 on the photosphere, equation ( 6) then gives thus making Lorentz force a plausible driver for the photospheric motions and an apt candidate to initiate MHD simulations.The NFFF extrapolation exploits a magnetic field B satisfying an inhomogeneous double curl Beltrami equation (Bhattacharyya et al., 2007) having a and b as constants.The solenoidality of B imposes ∇ 2 ψ = 0. Notably, the double-curl equation represents a self-organized state satisfying the MDR principle, (see Bhattacharyya and Janaki, 2004, and references therein for details).To solve the double-curl equation, an auxiliary field B ′ = B − (∇ψ)/b (Hu and Dasgupta, 2008)-satisfying the corresponding homogeneous equation is constructed.The equation represents a two-fluid steady state (Mahajan and Yoshida, 1998) and has a solution The B i are Chandrasekhar-Kendall eigenfunctions (Chandrasekhar and Kendall, 1957), obeying force-free equations with constant twists α i , and form a complete orthonormal set when the eigenvalues are real (Yoshida and Giga, 1990).Straightforwardly, where B 3 = (∇ψ)/b is a potential field.Combining equations ( 10) and (11), we have where the matrix V is a Vandermonde matrix having elements α i−1 j for i, j = 1, 2, 3, and α 3 = 0 (Hu and Dasgupta, 2008).The double-curl is solved by using the technique described in Hu et al. (2010).A pair of α i are selected and B 3 is SOLA: main.tex;23 January 2024; 3:08; p. 8 set to B 3 = 0. Using B z from the observed magnetogram, along with α i , the z -components of B 1 and B 2 are obtained at the bottom boundary.Afterwards, a linear force-free solver is employed to extrapolate the transverse components of B 1 and B 2 .Subsequently, an optimal pair of α i is obtained by minimizing the average normalized deviation of the magnetogram transverse field (B t ) from its extrapolated value (b t = B 1t + B 2t ), quantified as where M =N 2 is the total number of grid-points on the transverse plane.The E n is further reduced by employing the following decomposition for B 3 where 3t , the transverse difference △b t = B t − b t is obtained and further utilized to estimate the z -component of , b t , E n are estimated, and the procedure is repeated until the value of E n approximately saturates with the number of iterations, making the solution unique.Importantly, the procedure alters the bottom boundary and a correlation with the original magnetogram is necessary to check for the accuracy.
For our purpose, the vector magnetogram at 15:12 UT, from the hmi.sharp cea 720s series (Bobra et al., 2014) of SDO/HMI is employed as the bottom boundary.Though the SHARP series accounts for projection and foreshortening effects, a disk-centered active region reduces the possibility of any distortion in the magnetogram.The magnetic field components on the photosphere are obtained as B r , B p , and B t , which satisfy (a) B z =B r (r; radial), (b) B x =B p (p; poloidal), and (c) B y =−B t (t; toroidal) in a Cartesian coordinate system.The dimensions of the observed magnetogram is 877×445 pixels (≈ 317.91 Mm×161.31Mm).To save computational cost, we suitably crop and scale the magnetogram to new dimensions of 216×110 pixels (≈ 313.2 Mm×159.5 Mm) and extrapolation is carried out in a computational box defined by 216×110×110 voxels, where a voxel represents a value on a regular grid in 3-D space.The cropping and scaling procedures render the relative changes in positive and negative magnetic fluxes to be 0.02 % and 0.84 %, respectively.These changes are minimal and approximately satisfy the condition, B z = const.in the MHD simulation.

Robustness of Extrapolated Magnetic Field
The robustness of magnetic field extrapolation is quantified by computing the following parameters.
(a) The current weighted average (σ j ) of the sine of angle (σ i ) between current density and magnetic field (Wheatland, Sturrock, and Roumeliotis, 2000), as SOLA: main.tex;23 January 2024; 3:08; p. 9 defined in equation ( 15) where i runs over all the voxels in the computational box.Afterwards, the sine inverse of σ j (denoted by θ) is computed, which represents the average angle between current density and magnetic field.In our case, θ = 63.73 • , which is expected because the model is non-force-free at the bottom boundary.
(b) The fractional flux (Gilchrist et al., 2020), which quantifies the divergence free condition of the magnetic field, as defined in equation ( 16) where ∂S i represents the surface area of any voxel and S i it's volume.We find this value to be 3.366 × 10 −9 , which is numerically small enough to justify the solenoidal property of extrapolated magnetic field.
(c) The ratio of total magnetic energy with respect to the total potential state energy, denoted by E NFFF /E P .It allows to evaluate the free magnetic energy, hence sheds light on the capability of model to account for energy released during the transient phenomenon.The ratio turns out to be 1.305 in our case, hence, suggests that the extrapolated magnetic field has ≈ 30.5% more energy than the potential field.Quantitatively, the amount of available free energy is 5.6 × 10 31 ergs, which is enough to power a GOES M class flare.Further, to characterize the extrapolated magnetic field in the solar atmosphere, we check the variation of horizontally averaged magnetic field strength (|B| H av ), current density (|J| H av ), Lorentz force (|J×B| H av ), and θ H av with height, as shown in Figure 3.The averages are defined as where l denotes the voxel index.The averages are computed over layers defined by the 2D arrangement of voxels having N = 216 × 110 voxels along the x− and y−directions, respectively.Panel (a) reveals a continuous decrease of magnetic field strength with height and with respect to the bottom boundary, the percentage decrement is 84.4 %.The profiles in Panels (b) and (c) highlight the rapid decay of respective quantities, which approach saturation asymptotically.Within a distance of nearly 3 Mm, both the quantities decay by almost 50 %, but the subsequent decrease is relatively slower.In the higher layers of solar atmosphere, magnetic field lines tend to be more potential, thus characterized analytically by zero current density and low twist.Therefore, as shown in Panel (d), θ tends to increase with height.

Morphological Investigation of Magnetic Field Lines
Focusing on brightenings B 1 , B 2 , and the dome structure, magnetic field line morphologies in the active region are explored.Cospatial with B 1 , a hyperbolic flux tube (HFT; Titov, Hornig, and Démoulin, 2002) is found, as shown in Figure 4. Panel (a) of the figure shows magnetic field line linkage of the HFT configuration -constituted by the four quasi-connectivity domains (Zhao et al., 2014) in blue, yellow, pink, and red colors.These domains comprise of two intersecting quasi separatrix layers or QSLs (Titov, Hornig, and Démoulin, 2002) -one by blue and yellow magnetic field lines (MFLs), the other by pink and red MFLs.
Notably, these configurations are preferred sites for reconnection (Démoulin, 2006) and hence are of interest.QSLs are characterized by strong but finite gradients in the magnetic field line mapping.The gradients are quantified by the estimation of squashing degree Q, defined as follows.For the two footpoints of a field line, rooted in R a (x a , y a ) and R b (x b , y b ), the Jacobian matrix for the mapping ab : R a (x a , y a ) → R b (x b , y b ) is given by from which, the squashing degree Q is defined as where B n,a (x a , y a ) and B n,b (x b , y b ) are the normal components of magnetic field at the respective footpoints.We calculated the squashing degree by using the numerical code developed in Liu et al. (2016), available at http://staff.ustc.edu.cn/∼rliu/qfactor.html.In presence of strong gradients (high Q value), magnetic field lines undergo slippage while passing through the current layers, often referred as slipping reconnection (Aulanier et al., 2006).Panel (a) shows the ln Q map in a plane perpendicular to the bottom boundary, and crossing through the HFT morphology.It reveals the characteristic X-shape (ln Q ≥ 8) along HFT, which further confirms our interpretation of the morphology.Similarly, in Panel (b), regions of high gradient in the plane of the bottom boundary are seen to be nearly cospatial with B 2 , thus suggesting a plausible scenario for slipping reconnection.In particular, as evident from Panels (c) and (d), the footpoints of yellow MFLs lie on the boundary of dome, while those at one end of pink MFLs partially cover the periphery of dome.The presence of high gradients in footpoint mapping of the MFLs is indicative of slippage, which can possibly explain parts of brightening B 2 and chromospheric flare ribbons.The robustness of the extrapolated magnetic field and it's agreement with observations suggests that it can be reliably utilized as an input for the reported magnetohydrodynamics simulation.

Magnetohydrodynamics Simulation
To successfully simulate the coronal transients, the condition of flux-freezing must hold everywhere in the computational box, while allowing for magnetic reconnection at the plausible locations.In this work, the coronal plasma is idealized to be thermodynamically inactive and incompressible.The governing MHD equations in dimensionless form are SOLA: main.tex;23 January 2024; 3:08; p. 13 where R A F = (V A L)/ν is an effective fluid Reynolds number with V A as the Alfvén speed and ν as the kinematic viscosity.Hereafter, R A F is referred as fluid Reynolds number to keep the terminology uncluttered.The dimensionless equations are obtained by the normalization listed below.
In general, B 0 and L 0 are characteristic values of the system under consideration.Importantly, although restrictive, the incompressibility may not affect magnetic reconnection directly and has been used in earlier works (Dahlburg, Antiochos, andZang, 1991, Aulanier, Pariat, andDémoulin, 2005).Moreover, utilizing the discretized incompressibility constraint, the pressure p satisfies an elliptic boundary value problem on the discrete integral form of the momentum equation.
Toward simulating relaxation physics, it is desirable to preserve flux-freezing to an appropriate fidelity by minimizing numerical diffusion and dispersion errors away from the reconnection sites.Such minimization is a signature of a class of inherently nonlinear high resolution transport methods that prevent field extrema along flow trajectories while ensuring higher order accuracy away from the steep gradients in advected fields (Bhattacharyya, Low, and Smolarkiewicz, 2010).Consequently, equations ( 23)-( 26) are solved by the numerical model EULAG-MHD (Smolarkiewicz and Charbonneau, 2013), central to which is the spatio-temporally second order accurate, nonoscillatory, and forward-in-time Multidimensional Positive-Definite Advection Transport Algorithm, MPDATA (Smolarkiewicz, 1983, Smolarkiewicz and Margolin, 1998, Smolarkiewicz, 2006).For the computations carried out in the paper, important is the widely documented dissipative property of MPDATA that mimics the action of explicit subgrid-scale turbulence models (Margolin and Rider, 2002;Margolin, Smolarkiewicz, and Wyszogrodzki, 2002;Margolin, Smolarkiewicz, and Wyszogradzki, 2006), wherever the concerned advective field is under resolvedthe property referred to as Implicit Large-Eddy Simulations (ILES; Margolin, Rider, and Grinstein, 2006;Smolarkiewicz, Margolin, and Wyszogrodzki, 2007;Grinstein, Margolin, and Rider, 2007).Therefore, the effective numerical implementation of the induction equation by EULAG-MHD is SOLA: main.tex;23 January 2024; 3:08; p. 14 where, D B represents the numerical magnetic diffusion-rendering magnetic reconnections to be solely numerically assisted.Such delegation of the entire magnetic diffusivity to ILES is advantageous but also calls for a cautious approach in analyzing and extracting simulation results.Being localized and intermittent, the magnetic reconnection in the spirit of ILES minimizes the computational effort, while tending to maximize the effective Reynolds number of simulations (Waite and Smolarkiewicz, 2008).However, the absence of physical diffusivity makes it impossible to accurately identify the relation between electric field and current density--rendering a precise estimation of magnetic Reynolds number unfeasible.Being intermittent in space and time, quantification of this numerical dissipation is strictly meaningful only in the spectral space where, analogous to the eddy viscosity of explicit subgrid-scale models for turbulent flows, it only acts on the shortest modes admissible on the grid (Domaradzki, Xiao, and Smolarkiewicz, 2003), particularly near steep gradients in simulated fields.Such a calculation is beyond the scope of this paper.Notably, earlier works (Prasad et al., 2020, Yalim et al., 2022, Bora et al., 2022, Agarwal, Bhattacharyya, and Wiegelmann, 2022) have shown the reconnections in the spirit of ILES to be consistent with the source region dynamics of coronal transients and provide the credence for adopting the same methodology here.

Numerical Setup
The MHD simulation is carried out with bottom boundary satisfying the linetied condition (Aulanier et al., 2010).In our case, we ensure this by keeping B z and v z fixed at the bottom boundary.We have kept the lateral and top boundaries of the computational box open.The simulation is initiated from a static state (zero flow) using the extrapolated non-force-free magnetic field, having dimensions 216×110×110 which is mapped on a computational grid of .This reduction in fluid Reynolds number may be envisaged as a smaller Alfvén speed, which turns out to be V A | sim.≈ 0.125× V A | cor.for L sim. = 110 × 1450 km = 159.5 Mm.The total simulation time in physical units is equivalent to nt × ∆t × (L sim./V A | sim. ) ≈ 63.6 min., where nt = 15000.Importantly, although the coronal plasma with a reduced fluid Reynolds number is not realistic, the choice does not affect the changes in field line connectivity because of reconnection, but only the rate of evolution.Additionally, it saves computational cost, as demonstrated by Jiang et al. (2016).

Results and Analysis
The initial non-zero Lorentz force pushes the magnetofluid and generates dynamics.The overall simulated dynamics pertaining to magnetic relaxation can be explored from the evolution of the following grid averaged parameters SOLA: main.tex;23 January 2024; 3:08; p. 15 where, l denotes the voxel index and the volume V encloses the volume of interest.W V av and |Γ| V av measure the grid averaged magnetic energy and twist, whereas (J/B) V av serves as a proxy to quantify gradient of magnetic field.The plot of grid averaged magnetic energy, depicted in Panel (a), shows a continuous decrease (≈ 7%) and is in alignment with the possibility of magnetic relaxation through reconnection.To support this idea, Panel (b) plots the grid averaged twist |Γ| V av .Notably, Γ is also associated with magnetic helicity.The plot shows a decay of the average twist up to ≈ 40 minutes, followed by a rise.The initial decay is in conformity with the scenario of magnetic reconnection being responsible for untwisting of global field structure (Wilmot-Smith, Pontin, and Hornig, 2010) and reducing the complexity of field lines.The scenario of reconnection assisted relaxation is further reinforced by a similar variation of (J/B) V av (Panel (c)) since, reconnection is expected to smooth out steep field gradients.The rise in both the parameters is due to a current enhancement localized near the top of the computational domain-addressed later in the paper.
In the above backdrop, understanding dynamics of magnetic field lines involved in reconnection merits further attention.For the purpose, the computational volume is partitioned into three sub-volumes, as described in Figure 1.Our selection focuses on the hyperbolic flux tube (HFT) as the principal reconnection site and on the observed extent of brightenings in the active region, as depicted in Figure 6.The two-dimensional projections of sub-volumes S 1 , S 2 , and S 3 are shown in cyan, green, and yellow color boxes.Sub-volume S 1 enshrouds brightening B 1 and is centered on the X-point of HFT, thus consisting of those regions where the development of strongest current layers is possible.S 2 encloses the HFT morphology that envelops B 1 and partly B 2 such that the field line connectivities of depicted MFLs (see Figure 4) are contained within S 2 .Lastly, S 3 covers the complete spatial extent of the observed brightening (see Figure 2) and full vertical height of the computational box.Importantly, magnetic energy in a sub-volume can change because of an interplay between dissipation, Poynting flux, and the conversion of kinetic energy SOLA: main.tex;23 January 2024; 3:08; p. 17 to magnetic energy.Consequently, the following sections analyzes magnetofluid evolution for each sub-volume.For completeness, the analyses are augmented by plotting the corresponding grid averaged current densities: in usual notations.

Sub-volume S 1
The time evolution of W V av , J V av , |Γ| V av , and (J/B) V av is depicted in Panels (a), (b), (c), and (d) ofFigure 7, respectively.To understand their dynamical evolution, the time duration of numerical simulation is partitioned into five phases, denoted by P (i) 1 , where, i = 1, 2, ..., 5.The composition of S 1 has five layers along the vertical direction, denoted by z 0 = 0, 1, ..., 4, and the contribution of each layer in the shaping of parameter profile needs to be examined.
W V av increases initially up to phase P 1 (4) , followed by a continuous decay during P 1 (5) .Auxiliary analysis (not shown here) reveals that the magnetic energy at all layers have a profile similar to W V av .Contrarily, the same is not true for J V av and |Γ| V av -an explanation for which is presented below.Owing to z 0 = 0, 1, J V av decreases sharply in the beginning phase P 1 (1) .During the rising phase P 1 (2) , while all the layers exhibit similar profile, only z 0 = 2, 3, 4 contribute most significantly because the X-point of HFT exists at these heights, which is a prominent site for development of strong currents.Lastly, from phase P 1 (3) to P 1 (5) , there is an overall decrease in J V av , again due to z 0 = 2, 3, 4 along with some wiggling-predominantly due to z 0 = 0, 1, 2.
The |Γ| V av profile during phases P 1 (1) and P 1 (2) is shaped by z 0 = 0, 1.The decline during P 1 (3) is attributed to z 0 = 2, 3, 4, while the evolution in P 1 (4) and P 1 ( 5) is strongly determined by the bottom two layers, i.e., z 0 = 0, 1.The pronounced effect of bottom two layers (z 0 = 0, 1) could be due to the fact that spatial structures near the bottom boundary are not sufficiently resolved due to reduction in resolution of the extrapolation.Consequently, the dynamics in near neighborhood of the X-point of HFT leads to fluctuations in the profile of J V av and |Γ| V av , which do not smooth out due to the small size of sub-volume S 1 .Lastly, it is seen that the evolution of (J/B) V av is qualitatively similar to |Γ| V av profile.The quantitative changes in parameters during each of the phases are summarized inFigure 2 and from the rightmost column of this table, we note that the net magnetic energy and current density have increased while the overall twist and gradients have reduced in sub-volume S 1 .The evolution of magnetic field line dynamics, as shown in Panels (a) and (b) ofFigure 8, reveals that the field lines change their connectivity as soon as the simulation is initiated.The change in connectivity occurs because of reconnection at the X-point of HFT.
SOLA: main.tex;23 January 2024; 3:08; p. 19 With reconnection being known to dissipate magnetic energy, the increase of W V av demands additional analysis.With field line twist decreasing, the energy may increase if the net energy flux entering the sub-volume S 1 supersedes the energy dissipation at the X-point of the HFT.Such an analysis requires estimations of Poynting flux and dissipation to high accuracy, which is presently beyond the scope of this article as stated upfront in the paper.Nevertheless, an attempt is made toward a coarse estimation.A variable |D| is defined to approximate the D B in equation ( 28), indicating when and where non-ideal effects can be important.

|D|
Toward evaluating energy flux entering or leaving the sub-volume, an approximate estimation of Poynting flux is attempted with only the ideal contribution of electric field satisfying because of the ILES nature of the computation.Using straightforward vector analysis, the Poynting flux (Kusano et al., 2002) across the bounding surface (B) of a sub-volume can be written as where N has usual meaning, n and t mark the normal and tangential components to the area element vector denoted by ∆a, respectively.Notably, v z remains zero at the bottom boundary throughout the computation because of the employed boundary condition and the initial static state-see section 4.1.Consequently, only the second term contributes to the Poynting flux through the bottom boundary.
In Figure 9, Panel (a) plots the two-dimensional data planes of temporally averaged (averaged over the total computation time) |D| extracted from its 3D data volume, using the Slice Renderer function of VAPOR (Li et al., 2019)   (5) 1 but in complete disagreement for P (2) 1 and briefly for P (3) 1 .An absolute reasoning for this disagreement is not viable within the employed framework of the model, nevertheless, a possibility is the model not being in strict adherence to the equation ( 24) -leading to a violation of equation (34).

Sub-volume S 2
The evolution of W V av , J V av , |Γ|   The exploration of dynamics reveals that each layer along the vertical direction of computational box (denoted by z 0 = 0, 1, ..., 19) has nearly similar profile for magnetic energy, while for J V av and |Γ| V av , this is not true.
SOLA: main.tex;23 January 2024; 3:08; p. 22 The sharp decline in J V av during phase P (1) 2 is predominantly caused by z 0 = 0, 1.The rising phase P (2) 2 is caused by the layers z 0 = 2 to 9, with dominant contributions from z 0 = 2, 3, 4 and maximum from z 0 = 3.Similarly, the declining P (3) 2 phase is shaped by layers z 0 = 0 to 4, but the most significant role is played by layers z 0 = 1, 2, 3, while the maximum contribution arises from z 0 = 2.In the later phase, i.e.P (4) 2 , J V av increases again because of z 0 = 11 to 19.Notably, during P (4) 2 , the layers z 0 = 2 to 10 display declining values of current density, thus suggesting that while current density decreases in lower layers, the overall phase is governed by the dynamical evolution in higher layers.Lastly, in the concluding phase P (5) 2 , layers from z 0 = 0 to 15 exhibit decrease of current density, thus resulting in an overall decay.
The |Γ| V av profile reveals sharp decline during phases P (1) 2 and P (2) 2 , primarily due to initial five to six layers (z 0 = 0 to 5), but as in case of J V av , z 0 = 0, 1 determine the overall profile during these phases.The subsequent rising phases P (3) 2 and P (4) 2 are seen to be governed by layers z 0 = 8 to 19 and z 0 = 12 to 19, respectively.Notably, during these two phases, the lower layers, identified by z 0 = 0 to 7 and z 0 = 0 to 11 show lowering of twist over time.This behavior is reminiscent of J V av during P (4) 2 .In the end phase P (5) 2 , all except the top four layers, show lowering of twist, thus resulting in an overall decaying profile.During the early phases, i.e., up to P (3) 2 for J V av and P (2) 2 for |Γ| V av , the lower layers (z 0 = 0, 1, ..., 5) of sub-volume S 2 are seen to be playing the major role in determining the evolution of grid averaged parameters.This is due to the fact that the non-ideal region (the X-point of HFT) is within the first five layers of bottom boundary.Since S 2 contains S 1 , the reconnection at X-point plays an important role during the beginning phase of W V av .Moreover, the energy reduction has added contribution from other sources as S 2 covers the observed brightening B 2 as well.We explored an instance of this possibility using the anticipated slipping reconnection in yellow and pink MFLs constituting the observed dome structure.Panels (a) and (b) inFigure 11 depict a situation where sudden flipping of three selective magnetic field lines occurs (more profound in the animation provided in supplementary materials), which implies slipping reconnection.To facilitate easy identification, the footpoints of the three field lines are marked with black, white, and red colored circles.The increase in twist from P  In S 2 |D| ∈ {0.0, 0.05}, which is smaller compared to that in S 1 (marked by the black colored box in Panel (a)), signifying larger values of |D| to be localized at S 1 .The Poynting flux is positive for most of the P (2) 2 , which is in conformity with the energy decay.For phases P (3) 2 to P (5) 2 , the Poynting flux is negative, which can further be visualized from FigureFigure 13, where a portion of green field lines are pushed completely inside S 2 (red colored box).The corresponding energy influx along with the increment in twist seems to overwhelm dissipation, thus resulting in the observed energy increase.SOLA: main.tex;23 January 2024; 3:08; p. 24

Sub-volume S 3
The sub-volume S 3 encompasses the complete extent of the observed brightening.For convenience, the evolution in S 3 is investigated in five phases, defined by P (i) 3 , where, i = 1, 2, ..., 5.The temporal evolution of the grid averaged parameters is shown inFigure 14.
SOLA: main.tex;23 January 2024; 3:08; p. 25 Panel (a) reveals that S 3 exhibits continuous decrease in W V av up to P (4) 3 , which is in close agreement with the end time of the flare.Such uninterrupted decrement is a prime signature of relaxation in the considered volume.From Panel (b), it is seen that after an initial drop, J V av peaks at 15:27 UT, subsequently followed by a declining profile.Panel (c) indicates that |Γ| V av decays up to 15:39 UT, which nearly corresponds to the peak time of the flare.This suggests lowering of overall twist and hence a simplification of field line complexity, which further complements the interpretation of relaxation within the sub-volume.In the later phase, there is an increase in twist while the magnetic field gradient (Panel (d)) is seen to be declining continuously with very small increment toward the end of simulation.Overall, the volume averaged MHD evolution in S 3 is similar to the overall simulated dynamics.The quantitative changes associated with the grid averaged profiles are summarized inFigure 4. Notably, in this sub-volume, the terminal state is characterized by a reduced value of all the parameters, i.e. magnetic energy, current density, twist, and magnetic field gradients.
Notably, sub-volume S 3 spans the full extent of observed brightenings and the full vertical extent of the computational box.Consequently, to understand the SOLA: main.tex;23 January 2024; 3:08; p. 26 Table 4. Summary of the quantitative changes in grid averaged profiles of magnetic energy (W V av ), current density (J V av ), twist parameter(|Γ| V av ), and magnetic field gradient ((J/B) V av ) for sub-volume S 3 , during phases P  ) is carried out.These follow the definitions given in 29, 30, and 32 but grid averaged over different z = z 0 layers, each having N = 70 × 60 voxels along the x− and y−directions, respectively.Figure 15 shows the temporal profile of W H av for selected layers.During phase P (1) 3 , all the layers exhibit decreasing W H av with significant contribution from z 0 = 0, 1, 2, thus leading to declining W V av .Further, as shown in Panels (a) and (b), the subsequent phases are seen to have increasing W H av for z 0 = 0 to 3 in P (2) 3 , z 0 = 1 to 10 in P (3) 3 , z 0 = 0 to 14 in P 3 , and z 0 = 0 to 17 in P (5) 3 , respectively.Note that the evolution of W H av differs by an order of magnitude in the two Panels.The remaining layers during these phases exhibit declining W H av , as evident from Panels (b), (c), and (d).In effect then, owing to their larger number, these remaining layers dominate the profile evolution of W V av during phases P (2) 3 , P 3 , and P (4) 3 .However, in the end phase, the dynamics in z 0 = 0 to 17 takes control, thus leading to increasing W V av during P (5) 3 .
SOLA: main.tex;23 January 2024; 3:08; p. 27 Due to larger volume of S 3 , the resulting W V av profile is jointly governed by both larger (smaller) decrements in the lower (higher) layers.However, for S 2 , whose vertical extent is restricted to z 0 = 19, layers from Panel (a) and partly from Panel (b) can be visualized to jointly reproduce an initial fall, followed by continuous rise.Next, the behavior of J H av is explored, as shown inFigure 16.

During P
(1) 3 , other than the top two, all layers exhibit decreasing J H av , thereby causing the sharp decline of J V av in this phase.Notably, the dominant role is played by the bottom layers z 0 = 0, 1, as may be seen from Panel (a).Subsequently, in the next two phases, the process of current formation and dissipation within the HFT governs the evolution.As depicted in Panels (a) and (b), the increase in J V av during P (2) 3 is essentially due rising J H av in layers z 0 = 2 to 11.Similarly, the decreasing J H av in z 0 = 0 to 7 causes the decline of J V av during phase P (3) 3 despite the increasing J H av in z 0 = 8 to 20.In the remaining two phases P (4) 3 and P (5) 3 , the segregation of any dominant contribution from z = z 0 layers was found to be difficult.However, the J V av profile is understood from the finding that J H av decreases significantly in layers z 0 = 23 to 107 and z 0 = 26 to 109, respectively, as evident from Panels (c) and (d).Interestingly, the two topmost layers reveal an abrupt increase in J H av , an understanding of which requires investigation of field line dynamics.Lastly, the behavior of |Γ| H av is investigated, as depicted inFigure 17.
SOLA: main.tex;23 January 2024; 3:08; p. 30   Such localization of |D| is compatible with the general idea of ILES.On an average the Poynting flux is ≈ 30% of its value for S 2 and is predominantly negative, implying energy influx.As a consequence, the decrease in magnetic energy can be attributed to the overall decrease in twist conjointly with nonideal effects, contributed primarily from S 1 and further augmented by slipping reconnections in S 2 .

Extent of Magnetic Relaxation
As describe earlier, force-free states are characterized by field aligned current density.Further, the distribution of α distinguishes between the nonlinear and linear force-free states.Consequently, in order to check the extent of relaxation in our simulation, we compare the histogram of angle (θ) between J and B at the beginning and end of simulation for each of the sub-volumes.The θ plots inFigure 20 utilize the transformation 180 • − θ to map θ ≥ 90 • in the range 0 • ≤ θ ≤ 90 • .Panels (a) and (b) for sub-volumes S 1 and S 2 reveal wide distributions extending over the entire range of angles for both the time instants.On the other hand, Panel (c) for S 3 shows comparatively narrow distributions peaking around 90 • .Presumably, this is due to the fact that S 3 spans the full vertical extent of the computational box and the variation of θ along height in non-force-free extrapolation model exhibits an increasing trend up to 90 • (seeFigure 3).Due to small size of S 1 and hence limited number of voxels, we could not identify any trend in the variation of θ with time except that the distribution is wide, which does not support the presence of field aligned current in terminal state.However, careful comparison of the blue and red profiles in Panels (b) and (c) suggests that during simulation, fraction of voxels with θ ≥ 60 • in S 2 and S 3 decrease, which we estimate to be 20% and 24%, respectively.This suggests that the magnetic configuration tends to relax towards a force-free state.However, in the present simulation, neither the wide distribution in S 2 nor the narrow distribution centered around θ = 90 • in S 3 support a strictly field aligned current density.Therefore, in our case, the terminal state of the simulation remains in nonequilibrium, suggesting that further magnetic relaxation is possible.To check this, we carried out an auxiliary simulation, extending its time duration to twice that of the original one.When integrated over the whole computational domain, the grid averaged angle drops by 5.7 • (64.32 • to 58.62 • ) in the auxiliary simulation, as compared to 4.3 • (64.32 • to 60.01 • ) in the original simulation, validating the possibility of further relaxation.Furthermore, we explored the time evolution of the twist parameter (Γ) for each of the sub-volumes, as shown inFigure 21.The red and blue colors represent the negative and positive values, respectively.We note that at the initial time instant, the distribution is dominated by positive Γ for each sub-volume, as shown in Panels (a),(c), and (e).As the simulation progresses, negative Γ begins to increase and terminal state consists of both positively and negatively signed values.The noteworthy aspect is the progressively increasing intermixing of the blue and red colors, exhibiting gradual fragmentation (better visualized in the movie provided with supplementary materials) into smaller structures, as evident from Panels (d) and (f).Such fragmentation is indicative of development of turbulence (e.g.Pontin et al., 2011, also see Veltri et al., 2009) but since, a quantitative investigation regarding the extent of developed turbulence is presently beyond the scope of this work, it is difficult to comment on this aspect further.

Summary and Discussion
In this study, we explore the process of magnetic relaxation during an observed solar flare.For the purpose, an extrapolated magnetic field is utilized as input to carry out a data-based numerical simulation.We select a M1.3 class flare on 04 January, 2015, hosted by the active region NOAA 12253, spanning a time of nearly 35 minutes.Observations of the flare in 1600 Å and 304 Å channels of SDO/AIA (seeFigure 2) reveal the spatiotemporal evolution of identified brightenings (namely B 1 and B 2 ), the existence of chromospheric flare ribbons, and a dome shaped structure toward the eastern side.To explore the potential reconnection sites and their morphologies, the vector magnetogram at 15:12 UT from SDO/HMI is employed to extrapolate the magnetic field using a nonforce-free field (NFFF) model.We find a hyperbolic flux tube (HFT), overlying the dome structure and the brightenings B 1 , B 2 .In general, identification of SOLA: main.tex;23 January 2024;3:08;p. 35 all the individual reconnection sites within the computational volume is a nontrivial exercise and since, the HFT is spatially correlated with the observed brightenings, we envisage it as the primary reconnection site and focus on it for further analysis.To investigate reconnection dynamics and magnetic relaxation, the EULAG-MHD model is employed to execute an Implicit Large Eddy simulation (ILES), which uses the extrapolated magnetic field as initial condition.The simulation relies on the intermittent and localized dissipative property of the advection scheme MPDATA which smooths out under-resolved variables and numerically replicates magnetic reconnection.Although quantification of this dissipation is strictly meaningful only in the spectral space, nevertheless a rough estimation is carried out by analyzing |D|: the deviation of induction equation from its ideal form.
Toward exploring magnetic relaxation, we consider the temporal evolution of grid averaged parameters such as magnetic energy (W V av ), twist parameter (|Γ| V av ), magnetic field gradient ((J/B) V av ), and current density (J V av ), as defined in equations ( 29)-(32).For an overall picture, analysis of the full computational domain (seeFigure 5) reveals that W V av , |Γ| V av , and (J/B) V av decrease with time, indicating magnetic relaxation.For a more detailed analysis, we select three subvolumes of interest, namely S 1 , S 2 , and S 3 (see Figure 1), where S 1 is centered on the X-point of the HFT, S 2 focuses on the HFT configuration, and S 3 covers the full spatial extent of observed flaring region.To investigate the dynamics, we divide the simulation time into five phases, labeled by P 3 (where i = 1, 2, ..., 5) for sub-volumes S 1 , S 2 , and S 3 , respectively.Broadly, we find that in all sub-volumes the final values of |Γ| V av and (J/B) V av are smaller than their initial values, indicating a reduction in both twist and field gradient-consistent with the scenario of relaxation.Further, common to all sub-volumes, a sudden drop in J V av during the initial phase (e.g.P 1 (1) for S 1 ) is governed prominently by layers adjacent to bottom boundary, indicating a possible boundary condition effect.
Unlike the magnetic energy averaged over the whole computational grid, its grid averages over the sub-volumes do not decay monotonically.Toward understanding variations of magnetic energy in different sub-volumes, properties of |D|, Poynting flux, and magnetic twist in each sub-volume are explored.Importantly, large values of |D| are found to be localized at S 1 , particularly coinciding with the location of the X-point of HFT.This is harmonious with the spirit of ILES.In S 1 , apart from the phase P 1 (2) and briefly in P 1 (3) , the magnetic energy evolution is in conformity with physical expectations.The disagreement could be because of a failure of idealized Ohm's law as the induction equation in its ideal limit is not satisfied.Similar analyses have been carried out to explore energy variations in S 2 and S 3 also.In S 2 , the energy influx along with twist overwhelms dissipation, thus resulting in the observed energy increase from P 2 (3) to P 2 (5) .The decrease in magnetic energy in S 3 is found to be due to non-ideal effects primarily localized in S 1 .Toward an estimation of the extent of relaxation, the angle between the current density and magnetic field (θ) at every voxel is calculated.When integrated over the whole computational domain, the grid averaged angle drops by 4.3 • .For the sub-volumes, it is found that the changes in θ distribution over the course of simulation are not very clear for S 1 due to its SOLA: main.tex;23 January 2024; 3:08; p. 36 small size.In S 2 and S 3 , the peak of θ distribution becomes smaller, as realized from the decrease in fraction of voxels having θ ≥ 60 • .The decrease in higher values of θ indicates increase in alignment between current density and magnetic field.
In tandem, the above results indicate an ongoing magnetic relaxation, though it's extent remains to be explored.The angular distribution between current density and magnetic field suggests that though there is magnetic relaxation, but not enough to reach a force-free state.The terminal state of the simulation remains in non-equilibrium, suggesting the possibility for further relaxation.An auxiliary simulation with twice the computational time shows further alignment of electric current density with magnetic field, but at a slower rate.This is expected as the corresponding time span overlaps with the observed post-flare phase where reconnection plays a secondary role.Overall, the simulation suggests the extent of solar fare induced magnetic relaxation depends on the flare energetics and its duration.To further contemplate, magnetic reconnections are localized in the flaring regions.Consequently, a flaring region can exchange magnetic helicity with its surroundings.Under such circumstances, invariance of helicity is nontrivial and a complete field alignment of electric current density may not be achieved.An explicit calculation of magnetic helicity and understanding its evolution is necessary to focus on this idea-which we leave as a future exercise.

Figure 1 .
Figure 1.(a) The line-of-sight magnetogram at 15:12 UT on 04 January, 2015, from SDO/HMI.The corresponding dimensions in pixel units are 4096×4096 (b) CEA projected and cropped magnetogram (based on HARP active region patch) at 15:12 UT, having dimensions 877× 445 in pixel units for AR NOAA 12253.Both the figures are scaled to represent magnetic field strength within ± 1000 Gauss, with black patches representing the negative polarity and white patches representing the positive polarity.

Figure 2 .
Figure 2. Snapshots from observations of the solar flare in 1600 Å and 304 Å channels of SDO/AIA.Panels (a) and (b) reveal the brightening locations in 1600 Å during the beginning and peak phases of flare, marked by B 1 and B 2 .Panels (c) and (d) highlight the initial configuration of the identified dome-shaped structure and chromospheric flare ribbons during the peak phase of flare.A zoomed in image of the boxed region in Panel (c) is presented in Panel (e) with enhanced image contrast.The yellow color lines represent the manual tracing of the structure, while C labels the central location, where all the lines meet (The spatiotemporal evolution of flare in 1600 Å and 304 Å channels is available as high resolution movie in supplementary material).

Figure 3 .
Figure 3. Variation of horizontally averaged (a) magnetic field strength (b) current density (c) Lorentz force, and (d) θ with height for the extrapolated magnetic field.

Figure 4 .
Figure 4. Regions of high squashing degree and morphology of the hyperbolic flux tube (HFT).Panels (a), (b) and Panels (c), (d) are overlaid with observations in 1600 Å and 304 Å channel of SDO/AIA at the bottom boundary, corresponding to the same time instants as in Figure 2. In all the Panels, the map of squashing degree ln Q is given with color table.In Panel (a), the map is perpendicular to the bottom boundary, crossing through the HFT, while in Panels (b), (c) and (d), the ln Q map is in the plane of the bottom boundary.Panel (a) shows HFT from a side view, while Panels (c) and (d) show HFT from a top-down view.Panels (a) and (c) use a zoomed in viewpoint while Panels (b) and (d) use a zoomed out viewpoint.In Panel (c), the dome shaped structure is marked with a white color box.

Figure 6 .
Figure 6.Visual representation of sub-volumes S 1 , S 2 , and S 3 .The hyperbolic flux tube (HFT) configuration, overlaid with the vertical component of magnetic field and observation of the flaring event in 304 Å channel of SDO/AIA at 15:35:52 UT is shown.A zoomed in viewpoint is used, corresponding to a cutout of 150 × 90 pixels.The extent and spatial position of sub-volumes S 1 , S 2 , and S 3 are marked by the cyan, green, and yellow colored boxes.The arrows indicate the extent of sub-volumes along the z-direction and are drawn in proportion to the actual vertical sizes of sub-volumes given inFigure 1.

Figure 7 .
Figure 7. Time evolution of grid averaged (a) magnetic energy (W V av ) (b) current density (J V av ) (c) twist parameter(|Γ| V av ), and (d) (J/B) V av in sub-volume S 1 , during phases P (1) 1 (marked by the black arrow), P (2) 1 , P (3) 1 , P (4) 1 , and P . The dashed lines in blue, green, orange, and red colors separate the different phases in each of the profiles.The origin of the time scale maps to 15:12 UT.

Figure 8 .
Figure 8. Panels (a) and (b): Illustration of changes in the field line connectivity of yellow, blue, and red MFLs due to reconnection at the X-point of HFT configuration.The red colored box marks the edges of S 1 while the bottom boundary is overlaid with an image in 304 Å channel of SDO/AIA (The spatiotemporal evolution of the magnetic field line dynamics in sub-volume S 1 is available as movie in supplementary materials).
and |D| ∈ {0.01, 0.11}.Notably, |D| is largest in the neighborhood of the X-point of HFT depicted in Panel (a) of Figure 8 and decreases away from it.The Panel (b) plots S B av .A positive value of S B av indicates outflow of magnetic energy whereas negative value means energy influx.

Figure 9 .
Figure 9. (a) Two-dimensional data planes of temporally averaged |D| at three different heights in sub-volume S 1 .The mapping of data values is shown in the color bar.(b) Temporal evolution of S B av for sub-volume S 1 .The dashed lines in blue, green, orange, and red colors separate the different phases.The origin of the time scale maps to 15:12 UT.
2 are presented in Panels (a), (b), (c), and (d) ofFigure 10.Again, five phases are considered, denoted by P continuously.The quantitative changes corresponding to different phases are summarized inFigure 3, from which, a comparison of the terminal and initial states of the simulation reveals that the net magnetic energy in S 2 increases, while the other parameters decrease.

Figure 10 .
Figure 10.Time evolution of grid averaged (a) magnetic energy (W V av ) (b) current density (J V av ) (c) twist parameter(|Γ| V av ), and (d) (J/B) V av in sub-volume S 2 , during phases P (1) 2 (marked by the black arrow), P (2) 2 , P (3) 2 , P (4) 2 , and P (5) 2 (also marked by black arrow ), respectively.The dashed lines in blue, green, orange, and red colors separate the different phases in each of the profiles.The origin of the time scale maps to 15:12 UT.

Figure 11 .
Figure 11.Panels (a) and (b): Illustration of sudden shift in the footpoints of yellow and pink magnetic field lines within the sub-volume S 2 due to slipping reconnection.Panel (a) highlights the initial footpoints of three selective MFLs in black, white, and red colored circles.Panel (b) depicts the sudden movement of these footpoints, which is not along the direction of plasma flow (shown in white arrows).The bottom boundary is overlaid with squashing degree map and an image in 304 Å channel of SDO/AIA (The spatiotemporal evolution of the magnetic field line dynamics undergoing slipping reconnection in sub-volume S 2 is available as movie in supplementary materials).
accordance with the magnetic energy increase.To gain further insight, Figure 12 plots the time averaged deviation |D| and Poynting flux in Panels (a) and (b), respectively.

Figure 12 .
Figure 12.(a) Two-dimensional data planes of temporally averaged |D| at two different heights in sub-volume S 2 .The mapping of data values is shown in the color bar.The black box marks the sub-volume S 1 (b) Temporal evolution of S B av for sub-volume S 2 .The dashed lines in blue, green, orange, and red colors separate the different phases.The origin of the time scale maps to 15:12 UT.

Figure 13 .
Figure 13.Panels (a) and (b): Illustration of magnetic flux transfer within the sub-volume S 2 .The field lines comprising the HFT are shown, along with an additional set of green colored MFLs, which are pushed completely inside the sub-volume S 2 during simulation.The red colored box marks the edges of S 2 while the bottom boundary is overlaid with an image in 304 Å channel of SDO/AIA (The spatiotemporal evolution of the magnetic field line dynamics in sub-volume S 2 is available as movie in supplementary materials).

Figure 14 .
Figure 14.Time evolution of grid averaged (a) magnetic energy (W V av ) (b) current density (J V av ) (c) twist parameter(|Γ| V av ), and (d) (J/B) V av in sub-volume S 3 during phases P (1) 3 (marked by the black arrow), P (2) 3 , P (3) 3 , P (4) 3 , and P (5) 3 , respectively.The dashed lines in blue, green, orange, and red colors separate the different phases in each of the profiles.The origin of the time scale maps to 15:12 UT.
negative values indicate the rising and declining phases, while the net value in the rightmost column tells about the difference between terminal and initial states.P

Figure 15 .
Figure 15.Time evolution of W H av for sub-volume S 3 at different z = z 0 layers, shown by the black, magenta, cyan, and indigo color solid lines.The labels indicate the chosen z 0 value in each Panel.The different phases P (i) 3 , where i = 0, 1, ..., 5, are marked only in Panel (a) to avoid clutter, while the dashed lines separating the phases are marked in each of the Panels.The y-scale in Panels (b), (c), and (d) differs by an order of magnitude (10 −1 ) than in Panel (a).The origin of time on x-axis maps to 15:12 UT.

Figure 16 .
Figure 16.Time evolution of J H av for sub-volume S 3 at different z = z 0 layers, shown by the black, magenta, cyan, and indigo color solid lines.The labels indicate the chosen z 0 value in each Panel.The different phases P (i) 3 , where i = 0, 1, ..., 5, are marked only in Panel (a) to avoid clutter, while the dashed lines separating the phases are marked in each of the Panels.The origin of time on x-axis maps to 15:12 UT.

Figure 17 .
Figure 17.Time evolution of |Γ| H av for sub-volume S 3 at different z = z 0 layers, shown by the black, magenta, cyan, and indigo color solid lines.The labels indicate the chosen z 0 value in each Panel.The different phases P (i) 3 , where i = 0, 1, ..., 5, are marked only in Panel (a) to avoid clutter, while the dashed lines separating the phases are marked in each of the Panels.The origin of time on x-axis maps to 15:12 UT.

Figure 18 .
Figure 18.Illustration of magnetic field line dynamics responsible for abrupt rise of J H av and |Γ| H av in the top two layers of the computational box.The blue MFLs and red arrows in Panel (a) depict the bipolar potential field lines and direction of Lorentz force in the beginning of simulation.Panel (b) depicts the deformation in MFLs due to action of Lorentz force over the course of simulation.The bottom boundary is overlaid with image in 304 Å channel of SDO/AIA (The spatiotemporal evolution of these blue MFLs is available as movie in supplementary material).

Figure 19
Figure 19 plots slice rendering of time averaged |D| along with the Poynting flux.The |D| ∈ {0, 0.003}, which is one and two orders less than its values in S 2 and S 1 (marked in Panel (a) with arrows), respectively.Comparison of |D| in all the three sub-volumes indicates localization of maximal |D| at S 1 and specifically, at the neighborhood of the X-point-the primary reconnection site.Such localization of |D| is compatible with the general idea of ILES.On an average the Poynting flux is ≈ 30% of its value for S 2 and is predominantly negative, implying energy influx.As a consequence, the decrease in magnetic energy can be attributed to the overall decrease in twist conjointly with nonideal effects, contributed primarily from S 1 and further augmented by slipping reconnections in S 2 .

Figure 19 .
Figure 19.(a) Two-dimensional data planes of temporally averaged |D| at different heights in sub-volume S 3 .The mapping of data values is shown in the color bar.The black boxes mark the sub-volumes S 1 and S 2 (b) Temporal evolution of S B av for sub-volume S 3 .The dashed lines in blue, green, orange, and red colors separate the different phases.The origin of the time scale maps to 15:12 UT.

Figure 20 .
Figure 20.Distribution of angles between current density (J) and magnetic field vectors (B) in sub-volumes S 1 , S 2 , and S 3 at the beginning (blue) and end (red) of simulation.

Figure 21 .
Figure 21.Direct Volume Rendering (DVR) of the twist parameter (Γ) for sub-volumes S 1 , S 2 , and S 3 at the initial and terminal state of numerical simulation.The blue and red colors represent positive and negative values of Γ (The spatiotemporal evolution of distribution of twist in each sub-volume is available as a movie in supplementary materials).

Table 2 .
Summary of the quantitative changes in grid averaged profiles of magnetic energy (W V av ), current density (J V av ), twist parameter(|Γ| V av ), and magnetic field gradient ((J/B) V av ) for sub-volume S 1 , during phases P

Table 3 .
Summary of the quantitative changes in grid averaged profiles of magnetic energy (W V av ), current density (J V av ), twist parameter(|Γ| V av ), and magnetic field gradient ((J/B) V av ) for sub-volume S 2 , during phases P