Magnetohydrodynamic simulation of magnetic null-point reconnections in NOAA AR12192 initiated with an extrapolated non-force-free-field

The solar corona can be treated as a magnetized plasma having large electrical conductivity with evolution being determined by the magnetohydrodynamic (MHD) equations (Priest, 2014). The magnetic Reynolds number $R_{M}(vL/\eta$, in usual notations) for the corona is of the order of $10^{10}$ (Aschwanden, 2004), which makes the Alfvén’s theorem of flux freezing valid and ensures plasma-parcels to remain tied with magnetic field lines (MFLs) during evolution (Alfvén, 1942). The eruptive events (flares, coronal mass ejections) occurring at the corona are thought to be signatures of magnetic reconnection (MR): a process involving the topological rearrangement of MFLs with conversion of magnetic energy into heat and kinetic energy of mass motion (Shibata & Magara, 2011). Notably, the requirement to onset MRs is small $R_{M}$ which corresponds to small $L$, the length over which the magnetic field varies. The smallness of $L$ can either be pre-existing in a magnetic topology—manifested as magnetic nulls and quasi-separatrix layers (QSLs)—or can develop autonomously during the evolution of the magnetofluid. Such autonomous developments (owing to discontinuities in magnetic field) are expected from the Parker’s magnetostatic theorem (Parker, 1972, 1988, 1994) which states that for a perfect electrically conducting plasma, the conditions of flux-freezing and the equilibrium cannot be satisfied simultaneously by a magnetic field which is continuous everywhere. The reduction of $L$ and the consequent spontaneous magnetic reconnections during a quasi-static evolution of the plasma under a near-precise maintenance of the flux-freezing has been identified in contemporary MHD simulations (Kumar et al., 2015a; Kumar & Bhattacharyya, 2016; Kumar et al., 2016). However, these studies were performed using idealized scenarios of initial bipolar magnetic fields that lacked the complexities often observed in solar active regions.

Presently, the coronal field needs to be extrapolated from the photospheric magnetic field because of a lack of direct measurements. For extrapolation, the usage of the nonlinear-force-free-fields(NLFFFs), a subset of force-free-fields (Wiegelmann, 2008; Wiegelmann & Sakurai, 2012) is customary. The NLFFF can be solved analytically (in spherical polar coordinates) under the assumption of axisymmetry (Low & Lou, 1990; Prasad et al., 2014) but the analytical solution fails to effectively capture the complexity of an active region magnetogram, which is often non-axisymmetric. Such complexities are well replicated in NLFFF extrapolations (Duan et al., 2017). Recent MHD simulations based on NLFFF extrapolations were successful in simulating the coronal dynamics leading to eruptions (Jiang et al., 2013; Kliem et al., 2013; Amari et al., 2014; Inoue et al., 2014, 2015; Savcheva et al., 2015, 2016; Inoue, 2016). Importantly, only the region sandwiched between the photosphere and the upper corona is relatively force-free whereas at the photosphere—where magnetograms are obtained—the Lorentz force is non-zero (Gary, 2001). Generally, to mitigate this problem within the framework of NLFFF, a technique called ‘preprocessing’ is often performed on the photospheric data which minimizes the Lorentz force in the vector magnetograms and provides a boundary condition suitable for NLFFF extrapolations (Wiegelmann et al., 2006; Jiang & Feng, 2014).

An alternative is the extrapolation using non-force-free-fields (NFFFs) described by the double-curl Beltrami equation for the magnetic field ${\bf{B}}$ (Hu & Dasgupta, 2008). The equation has been analytically solved for an idealized corona (Bhattacharyya et al., 2007; Kumar & Bhattacharyya, 2011) to obtain MFLs resembling coronal loops. Recently, a semi-analytical construction based on maximizing correlations of non-axisymmetric NFFFs with photospheric vector magnetograms of NOAA AR11283 successfully mimicked an event of filament bifurcation by tracking MHD evolution of a pre-existing flux -rope (Prasad & Bhattacharyya, 2016; Prasad et al., 2017) However, missing from the simulation were the small scale magnetic features and their influence on the MFL dynamics—which cannot be captured by analytical/semi-analytical models. To include these magnetic features and determine their role in overall magnetofluid evolution, here we numerically simulate evolution of AR 12192 initiated with the NFFF extrapolation model developed by Hu & Dasgupta (2008); Hu et al. (2008); Gary (2009). The focus is to assess the viability of the NFFF extrapolation in generating MFLs, the evolution of which can reliably imitate the observed dynamics. For the purpose, we follow the evolution of AR 12192, starting at 20:46 UT on 24th October 2014 and study the X3.1 confined flare occurring at 21:10 UT. Importantly, a non-zero Lorentz force, instead of prescribed flows (Amari et al., 2003; Aulanier et al., 2010), is envisaged here to initiate dynamics.

In the rest of the paper, Section 2 discusses the flare-event and the observations required for the NFFF extrapolation. In Section 3, we present the details of the initial extrapolated field. The MHD model is discussed in Section 4. The results of the simulation are presented in Section 5 and the Section 6 summarizes important results.

The AR12192 was the largest of all active regions appearing in the solar cycle 24 which produced a series of X-class flares (Chen et al., 2015). The X3.1 flare on October 24, 2014 around 21:15 UT was the strongest in a series which did not lead to any coronal mass ejection (CME) (Sun et al., 2015; Sarkar & Srivastava, 2018). Since there is a very strong correlation between flare intensity and occurrence of CMEs (Yashiro et al., 2005), this event has been extensively studied. An absence of flux rope was suggested in Jiang et al. (2016) for explaining the confined nature whereas the onset of the flare was attributed to tether-cutting (TC) MRs (Moore et al., 2001) between sheared arcades. Further studies of successive strong X-class flares triggered by TC reconnections, in the same AR, were also reported in Chen et al. (2015). Contrarily, using NLFFF extrapolation, Inoue et al. (2016) found a multiple-flux-tube system located near a Polarity Inversion Line (PIL) to be favorable for the TC reconnections. They attributed the stability of the flux-tube-system to the overlying strong tethering MFLs. Similar results were also documented in Chen et al. (2015), where the mean decay index of the horizontal background field was found to be less than the typical threshold required for the torus instability (Kliem & Török, 2006) to set in. An alternative explanation was provided by Zhang et al. (2017) who attributed the confined nature to the complexity of the involved magnetic field structures.

The confined X3.1 flare was of long duration, lasting for 6 to 7 hours as shown in Figure 1(a). The figure shows the GOES 15 X-ray flux observed during this event in the 1-8 Å  and 0.5-4 Å  channels. It should be noted that no appreciable change in the vertical magnetic field flux was recorded during this period at the photospheric boundary. This is shown in Figure 1(b) which depicts the evolution of negative (dashed line in red) and positive magnetic fluxes (continuous line in blue), calculated by using the photospheric vector magnetograms from the Heliospheric Magnetic Imager (HMI)(Schou et al., 2012) on board the Solar Dynamics Observatory (SDO)(Pesnell et al., 2012). The magnetograms are taken from the ‘hmi.sharp_cea_720s data series’ that provides full-disk vector magnetograms of the Sun with a temporal cadence of 12 minutes and a spatial resolution of $0^{\prime\prime}.5$. In order to obtain the magnetic field on a Cartesian grid, the magnetogram is initially remapped on to a Lambert cylindrical equal-area (CEA) projection and then transformed into the heliographic coordinates (Gary & Hagyard, 1990). The dotted vertical lines mark the beginning and peak phase of the flare. Hence, to a good approximation, the vertical magnetic field $B_{z}$ at the bottom boundary remains constant during the interval. Accordingly, the photosphere can be approximated to be line tied — a boundary condition used in the simulations discussed later in the paper.

Importantly, a circular brightening was observed in the chromospheric flare ribbons at the ultra-violet (UV) 1600 Å  channel preceded by a brightening of the flaring loops in the extreme-ultra-violet (EUV) channel 131 Å  of the Atmospheric Imaging Assembly (AIA) on board SDO(Lemen et al., 2012)). The brightenings occur in the interval 21:20 to 21:35 UT in the 1600 Å  channel (Figure 2 (a)) and is co-located with the brightening in the 131 Å  channel as seen around 20:58 UT, which is just before the X-class flare. The circular flare ribbons are known to map MFLs constituting the fan plane of a 3D null on the photosphere (Masson et al., 2009). To our knowledge, the generation of the circular ribbon was not reported in the earlier works and is the main focus of the paper.

To simulate evolution of such MFLs, we select the vector magnetogram at 20:46 UT, roughly 30 minutes prior to the flare. The Figure 3(a) shows the magnetogram of the active region where the positive and the negative polarities of the longitudinal component of the magnetic field are depicted in white and black, and the gray represents the background. The transverse components of the positive and negative fields are shown by blue and red arrows respectively. The PIL is represented in the figure by green lines. The AR is visibly complex, with two main polarities and multiple small-scale features. The MFL topology can be inferred using the extreme ultra-violet (EUV) channel data as observed in the 171 Å , shown in Figure 3(b), which plotted on the same CEA spatial grid as in Figure 3(a). The EUV coronal loops near the PIL are markedly sheared and twisted, indicating a high degree of complexity in the initial magnetic field topology.

The evolution is governed by the incompressible Navier-Stokes MHD equations under the assumption of thermal homogeneity and perfect electrical conductivity (Bhattacharyya et al., 2010; Kumar et al., 2014, 2015b):

written in the usual notations in dimensionless form. The normalizations for various terms in Equation (3) are as follows

The constants $B_{0}$ and $L_{0}$ are fixed using the average magnetic field strength and length-scale of the vector magnetogram respectively. Here, $v_{a}\equiv B_{0}/\sqrt{4\pi\rho_{0}}$ is the Alfvén speed and $\rho_{0}$ is the constant mass density. The constants $\tau_{a}$ and $\tau_{\nu}$, having dimensions of time, represent the Alfvén transit time ($\tau_{a}=L_{0}/v_{a}$) and viscous diffusion time scale ($\tau_{\nu}=L_{0}^{2}/\nu$), respectively. The kinematic viscosity is denoted by $\nu$. The ratio $\tau_{a}/\tau_{\nu}$ represents an effective viscosity of the system which, along with the other forces, influences the dynamics.

To solve the MHD Equations (3a)-(3d), we utilize the well established magnetohydrodynamic numerical model EULAG-MHD (Smolarkiewicz & Charbonneau, 2013), which is an extension of the hydrodynamic model EULAG predominantly used in atmospheric and climate research (Prusa et al., 2008). The pressure perturbation, denoted by $p$, about a thermodynamically uniform ambient state satisfies an elliptic boundary value problem, which is generated by imposing the discretized incompressibility constraint (Equation 3b) on the discrete integral form of the momentum equation (Equation 3a); cf.(Bhattacharyya et al., 2010) and the references therein. An identical procedure involving the gradient of an auxiliary potential in the induction equation (Equation 3c) is employed to keep ${\mathbf{B}}$ solenoidal, see Ghizaru et al. (2010) and Smolarkiewicz & Charbonneau (2013) for details. For the completeness, here we mention only important features of the EULAG-MHD and refer the readers to Smolarkiewicz & Charbonneau (2013) and references therein for detailed discussions. The model is based on the spatio-temporally second-order accurate non-oscillatory forward-in-time multidimensional positive definite advection transport algorithm, MPDATA (Smolarkiewicz, 2006). Important is the proven dissipative property of the MPDATA which, intermittently and adaptively, regularizes the under-resolved scales by simulating MRs and mimicking the action of explicit subgrid-scale turbulence models (Margolin et al., 2006) in the spirit of Implicit Large Eddy Simulations (ILES) (Grinstein et al., 2007). Such ILESs performed with the model have already been successfully utilized to simulate magnetic reconnections (MRs) to understand their role in the development of various magnetic structures in the solar corona (Kumar et al., 2015a, 2016; Prasad et al., 2017). The simulations presented continue to rely on the effectiveness of ILES in regularizing the onset of MRs.

The simulations are initialized from a motionless state with the initial magnetic field given by the NFFF extrapolation and the magnetofluid idealized to be thermally homogeneous and having perfect electrical conductivity. The flow is generated as the initial Lorentz force pushes the plasma. To ensure the net magnetic flux to be zero in the computational domain, all components of volume ${\bf{B}}$ except for $B_{z}$, are continued to the boundaries for a given time step (Kumar et al., 2015a). At the bottom boundary, $B_{z}$ is kept constant (line-tied boundary) since the change of magnetic flux at the boundary is minimal (see Figure 1(b)). For the simulation, we set the dimensionless constant $\tau_{a}/\tau_{\nu}\approx 7\times 10^{-3}$, which is roughly two orders of magnitude larger than its coronal value. The higher value of $\tau_{a}/\tau_{\nu}$ speeds up the relaxation because of a more efficient viscous dissipation without affecting magnetic topologies. The density is set to $\rho_{0}=1$ and kinematic viscosity to $\nu=0.002$, in scaled units. The spatial unit step $\Delta x=0.0078$, while the time step is taken as $\Delta t=5\times 10^{-3}$ to satisfy the Courant-Friedrichs-Lewy (CFL) stability condition (Courant et al., 1967). The results presented here pertain to a run for 1000 $\Delta t$ which roughly corresponds to an observation time of one and half hour. Due to the constant mass density, the flow generated in the computation is incompressible—an assumption also used in earlier works (Dahlburg et al., 1991; Aulanier et al., 2005). Although the compressibility of the fluid is important for the thermodynamics of coronal loops (Ruderman & Roberts, 2002), our focus for the present is on their magnetic topology only. Notably, the $R_{M}$ throughout the simulation is infinity expect during MRs facilitated by the MPDATA driven dissipation.

Figures 9 and 10 depict MFL evolution in the neighborhood of the 3D null at two different viewing angles. The $B_{z}$ contours are plotted on the bottom boundary. Four sets of MFLs are highlighted. The fan and the spine of the null are made by the yellow MFLs whereas the red MFLs are overlying the null. The blue MFLs are located inside the dome whereas the arcade below the null is formed by the green MFLs. With evolution, the null and the constituent yellow MFLs do not sustain an appreciable ascent whereas the red MFLs expand significantly to a threshold height ($\approx 78$Mm), after which they contract. To explore the underlying physics, we note the arcade MFLs (in green) and the dome (yellow) constitute an X-type geometry cf. panel (b) of Figure 7. As reconnection occurs at the X-type null, blue MFLs come out of the dome and overlays it. The consequent increase in local magnetic pressure pushes the red MFLs upward, resulting in their overall rise. Furthermore, the red MFLs get stretched as they rise and at a threshold generate enough magnetic tension to stop additional upward motion. The threshold corresponds to a critical value of $\approx 1.5$ for the decay index, where the decay index is defined as $\displaystyle{n=-\frac{\partial(\log|\mathbf{B}|)}{\partial(\log z)}}$ (Kliem & Török, 2006). Figure 10 confirms, in their maximal rise, the MFLs can only attain $n\simeq 1.3-1.4$ which is in conformity with the confined nature of the X3.1 flare. Moreover, like (Jiang et al., 2016), we also fail to identify a flux rope, which further agrees with the confined nature of the flare.

The simulated 3D null appears to rotate with evolution (Figure 11). For aiding visualization, the MFLs have been color-coded based on their distance along the $y$ axis. We have also shown the volume wherein the seed points of MFLs are located. When viewed from the top, an anti-clockwise rotation of the MFLs is quite prominent which matched very well with the similar dynamics seen in the AIA 131 Å  channel. This correspondence with observations gives more credibility to the simulation. The Figure 12 is also overlaid with streamlines (green) and $|\mathbf{J}|/|\mathbf{B}|$. Noticeable is the initial high value of $|\mathbf{J}|/|\mathbf{B}|$ near the null. The value increases with time, becoming maximum at $t=400$, decaying subsequently. The peaking of $|\mathbf{J}|/|\mathbf{B}|$ is indicative of magnetic reconnections occurring near the null. The resultant outflow is shown by the red streamlines. For further investigation, Figure 13 plots the Q-map where the squashing factor Q is calculated by following Demoulin et al. (1996); Liu et al. (2016) and ascertains the dome to have high gradient of magnetic connectivity which results in slipping reconnections (Aulanier et al., 2007). The subsequent change in magnetic connectivity manifests as the seeming MFL rotation. For validation, we note the co-located flow (in green) is not along the rotation and hence, cannot cause it.

In panels (a) and (b) of Figure 14, we overlay intensity structures in wavelengths 1600 Å  at 21:25 UT and 131 Å  at 20:50 UT with corresponding MFLs. Importantly, the almost exact match of the footpoints for both wavelengths with brightenings not only establishes the importance of the 3D null in the circular flare ribbon but also being in agreement with the contemporary understanding, validates the effectiveness of the NFFF extrapolation in constructing a valid coronal field model.

The paper presents simulated dynamics of AR12192 from 20:48 UT. The plasma is idealized to have perfect electrical conductivity while being viscid, thermally homogeneous and incompressible. The simulations are initialized with magnetic field lines extrapolated from SDO/HMI vector magnetograms using a new technique which employs a model where the corona is not strictly force-free and has some Lorentz force. Nevertheless, the Lorentz force decreases rapidly with height making the corona to be force-free in an asymptotic limit—agreeing with the standard scenario of the coronal field. Advantageously, this non-force-free-field extrapolation model self-consistently initialize the coronal dynamics without requiring prescribed plasma flows which are somewhat custom-made.

The extrapolated magnetic field is found to have a 3D null located approximately at a height of 3 Mm from the photosphere and has clearly distinguishable spine and a dome shaped fan. Importantly, a magnetic arcade is found to be located within the dome and making an X-type null with it. A Q-map of the initial field identifies the dome with a region where the gradient of the field line connectivity is large.

The simulation focuses on a circular brightening recorded in 1600 Å  channel at around 21:21 UT. The absence of any flux emergence in the window of 19:00 UT-24:00 UT allows the vertical field at the bottom boundary to be assumed as line-tied. To optimize the computation cost, the simulation is performed on $256\times 128\times 128$ grids along the $x$, $y$ and, $z$ respectively, resolving a physical domain of $360\times 180\times 180$ Mm³. It is initiated not by a prescribed flow, but by the initial Lorentz force which onset the evolution autonomously. Subsequently, the favorable forces bring non-parallel field lines in close proximity of the 3D null, which ultimately leads to unresolved scales. In the spirit of ILES, the MPDATA then generates locally adaptive residual dissipation to regularize the underresolved scales with simulated MRs. Further, the MRs are found to be consistent with the idea of slipping reconnection, standard at a 3D null, and imparts a sense of rotation to the footpoints of the dome. Such a rotation is also observed in the channel Å, corroborating the observed circular brightening to be caused by MRs at the 3D null.

Magnetic reconnections also occur at the X-type null formed by the MFLs belonging to the arcade and the spine. Interestingly, reconnections enable MFLs contained within the dome to come out of it and overlay the spine. The consequent increase in magnetic pressure raises the overlying MFLs further up and in principle, can cause the X3.1 flare observed at 21:15 UT. The flare was confined in nature and resulted in no CMEs. In the simulation, the MFLs are found to never reach a height where the decay index becomes more the critical value required for the torus instability to set in, confirming further the efficacy of the simulation in replicating the observation.