Magnetohydrodynamic modeling of solar coronal dynamics with an initial non-force-free magnetic field

TL;DR

This paper presents a magnetohydrodynamic simulation of solar coronal magnetic fields initialized with a non-force-free, non-axisymmetric field, revealing complex dynamics, magnetic reconnections, and flux rope behavior consistent with observed solar flares.

Contribution

It introduces a novel MHD modeling approach starting from a realistic, non-force-free magnetic field to simulate solar active region dynamics.

Findings

Simulated flux rope undergoes asymmetric rise matching observations.

Reconnection occurs at specific X-type null lines, driving dynamics.

Spontaneous mass flow development is observed in the simulation.

Abstract

The magnetic fields in the solar corona are generally neither force-free nor axisymmetric and have complex dynamics that are difficult to characterize. Here we simulate the topological evolution of solar coronal magnetic field lines (MFLs) using a magnetohydrodynamic model. The simulation is initialized with a non-axisymmetric non-force-free magnetic field that best correlates with the observed vector magnetograms of solar active regions (ARs). To focus on these ideas, simulations are performed for the flaring AR 11283 noted for its complexity and well-documented dynamics. The simulated dynamics develops as the initial Lorentz force pushes the plasma and facilitates successive magnetic reconnections at the two X-type null lines present in the initial field. Importantly, the simulation allows for the spontaneous development of mass flow, unique among contemporary works, that preferentially reconnects field lines at one of the X-type null lines. Consequently, a flux rope consisting of low-lying twisted MFLs, which approximately traces the major polarity inversion line, undergoes an asymmetric monotonic rise. The rise is attributed to a reduction in the magnetic tension force at the region overlying the rope, resulting from the reconnection. A monotonic rise of the rope is in conformity with the standard scenario of flares. Importantly, the simulated dynamics leads to bifurcations of the flux rope, which, being akin to the observed filament bifurcation in AR 11283, establishes the appropriateness of the initial field in describing ARs.

Tables

Table 1. Table 1 : Parameters used in construction of the initial non-axisymmetric non-force-free field. Here, the parameters α 1 subscript 𝛼 1 \alpha_{1} , α 2 subscript 𝛼 2 \alpha_{2} , l 𝑙 l , m 𝑚 m , κ 𝜅 \kappa and γ 𝛾 \gamma determine a particular mode while r 0 subscript 𝑟 0 r_{0} , r ′ superscript 𝑟 ′ r^{\prime} r 1 subscript 𝑟 1 r_{1} , θ ′ superscript 𝜃 ′ \theta^{\prime} and ψ ′ superscript 𝜓 ′ \psi^{\prime} determine the extent of the computation domain.

$c$	${\alpha }_1,{\alpha }_2$	$l$	$m$	$\kappa$	$\gamma$	$r_0$	$r^{\prime }$	$r_1$	${\theta }^{\prime }$	${\psi }^{\prime }$
0.54	9.67, 13.23	3	2	0.125	0.375	0.1	0.34	1	1.83	0.13

Equations

$$\frac{\partial \mathbf{v}}{\partial t}+\left(\mathbf{v}\cdot \nabla \right)\mathbf{v}=-\nabla p+\left(\nabla \times \mathbf{B}\right)\times \mathbf{B}+\frac{{\tau }_a}{{\tau }_{\nu }}{\nabla }^2\mathbf{v},$$

$$\nabla \cdot \mathbf{v}=0,$$

$$\frac{\partial \mathbf{B}}{\partial t}=\nabla \times (\mathbf{v}\times \mathbf{B}),$$

$$\nabla \cdot \mathbf{B}=0,$$

$$\mathbf{B}⟶\frac{\mathbf{B}}{B_0},$$

$$\mathbf{v}⟶\frac{\mathbf{v}}{v_a},$$

$$L⟶\frac{L}{L_0},$$

$$t⟶\frac{t}{{\tau }_a},$$

$$p⟶\frac{p}{\rho {v_a}^2}.$$

$$\mathbf{B}={\mathbf{B}}^1+\kappa {\mathbf{B}}^2.$$

$$\nabla \times {\mathbf{B}}^i={\alpha }_i{\mathbf{B}}^i,$$

$${\mathbf{B}}^i=\frac{1}{{\alpha }_i}\nabla \times \nabla \times {\psi }^i\mathbf{r}+\nabla \times {\psi }^i\mathbf{r},$$

$${\psi }_{lm}^i(r,\theta ,\varphi )=[C_{lm}^1{j}_l({\alpha }_{i}r)+C_{lm}^2{y}_l({\alpha }_{i}r)]P_l^m(\cos \theta )\exp (im\varphi ),$$

$$B_r^i(r,\theta ,\varphi )$$

$$B_{\theta }^i(r,\theta ,\varphi )$$

$$B_{\varphi }^i(r,\theta ,\varphi )$$

$${\psi }_{lm}^i(r,\theta ,\varphi )=$$

$$c=\frac{⟨(\mathbf{B}\cdot {\mathbf{B}}_O)|{\mathbf{B}}_O|⟩}{⟨|\mathbf{B}{|}^3{⟩}^{1/3}⟨|{\mathbf{B}}_O{|}^3{⟩}^{2/3}},$$

Full text

magnetohydrodynamics (MHD) – Sun: activity – Sun: corona – Sun: flares – Sun: magnetic fields – Sun: photosphere

1 Introduction

The solar corona is generally treated as a magnetized plasma with high electrical conductivity whose evolution is determined by the equations of magnetohydrodynamics (MHD) (Priest, 2014). The large magnetic Reynolds number $R_{M}(vL/\eta$, in usual notations) of the corona allows the plasma to achieve its non-diffusive limit (Aschwanden, 2004). Consequently, Alfvén’s theorem of flux freezing is satisfied, ensuring magnetic field lines (MFLs) to remain tied to fluid parcels during evolution (Alfvén, 1942). In contrast, instances of various eruptive events (flares, coronal mass ejections) occurring in the corona are believed to be signatures of magnetic reconnection (MR): the topological rearrangement of MFLs along with the conversion of magnetic energy into heat and the kinetic energy of mass motion (Shibata & Magara, 2011). To trigger MRs, the value of $R_{M}$ needs to be lowered, which can be attained by a reduction in $L$, the length scale determining the gradient of the magnetic field. A vanishing $L$ can either pre-exist in a magnetic topology in the form of various two- and three-dimensional nulls and quasi-separatrix layers (QSLs), or develop autonomously. The autonomous development is assured by the Parker’s magnetostatic theorem (Parker, 1972, 1988, 1994), which claims that discontinuities of the magnetic field are ubiquitous in an infinite electrically conducting plasma in equilibrium, because the flux-freezing condition and the equilibrium generally cannot be satisfied simultaneously by a magnetic field that is continuous everywhere. In an evolving fluid, MHD simulations establish the autonomous reduction of $L$ and the occurrence of consequent spontaneous reconnections as the fluid relaxes toward a quasi-steady state under a nearly precise maintenance of flux freezing (Kumar et al., 2015a; Kumar & Bhattacharyya, 2016; Kumar et al., 2016). However, in these studies, the initial bipolar magnetic fields were constructed analytically, reducing the natural complexity often observed in active regions (ARs).

Without a direct measurement (Wiegelmann & Sakurai, 2012), the coronal fields are numerically extrapolated from the photospheric magnetic fields based on a force-free assumption; i.e. exert no Lorentz force is exerted (Régnier, 2013). The corresponding analytical expressions of the force-free field (in spherical polar coordinates) are often obtained under the assumption of axisymmetry (Low & Lou, 1990; Prasad et al., 2014). However, such a symmetry generally does not exist in the real case, as many vector magnetograms depict the photospheric magnetic fields to be morphologically complex.

In force-free coronal field extrapolations, the use of nonlinear-force-free-field (NLFFF) models where the field line twist is a function of position but constant along a field line (Wiegelmann, 2008) is customary. The model is well known to give good fits to magnetograms, yielding a realistic coronal field line topology. For instance, Guo et al. (2013) performed NLFFF extrapolation for NOAA 10720 and found coronal flux ropes wrapped around by a particular set of QSLs, whereas Inoue et al. (2012) found that in AR 10930, the X-ray intensity associated with the sigmoid is more sensitive to the strength of the electric current rather than the twist of the field lines. Furthermore, many recent works (Jiang et al., 2013; Kliem et al., 2013; Amari et al., 2014; Inoue et al., 2014, 2015; Savcheva et al., 2015, 2016) have realistically modeled the coronal eruptive dynamics by using MHD simulations initialized with NLFFFs (see also a review by Inoue (2016)).

Although simulations based on NLFFF are capable of reproducing realistic dynamics, further issues need to be addressed. Specifically, the force-free approximation is based on the assumption of negligible plasma $\beta$ in the corona, where $\beta$ is the ratio of gas pressure to magnetic pressure. However, according to Gary (2001), it is only in the mid-corona—sandwiched between the photosphere and the upper corona—that the plasma $\beta$ is negligible. To mitigate this problem, a ‘preprocessing’ technique is often performed on the data to remove the force in the vector magnetogram and obtain a suitable boundary condition for NLFFF extrapolations (Wiegelmann et al., 2006; Jiang & Feng, 2014).

An alternative to NLFFF extrapolation is an approach based on a magnetic field model of the corona derived from the variational principle of the minimum energy dissipation rate (Hu & Dasgupta, 2008; Hu et al., 2008, 2010). It is shown that the extrapolated field lines from the resulting non-force-free-field (NFFF) model reasonably agree with the observed coronal loops and MFLs extrapolated using NLFFF. Central to the NFFF extrapolations is the double-curl Beltrami field obtained by Bhattacharyya et al. (2007), which is a dissipating self-organized state. Kumar & Bhattacharyya (2011) have also shown the same double-curl Beltrami field to be the minimum total energy state of a two-fluid magnetoplasma. One particular analytical solution of the double-curl equation can be formulated by superposing two linear-force-free fields (LFFFs) characterized by constant twists everywhere cf. Steinhauer & Ishida (1998) & Mahajan & Yoshida (1998) for details. Kumar & Bhattacharyya (2011) utilized such formalism to obtain NFFFs with twisted MFLs resembling coronal loops. A NFFF relevant to observations of AR 11283 was recently obtained by Prasad & Bhattacharyya (2016) by using the field in conjunction with semi-analytical constructs where non-axisymmetric solutions are best correlated with vector magnetograms. Interestingly, the obtained NFFF yields a magnetic morphology where regions of the sharpest magnetic field gradient were found to be cospatial with the flaring site. To further establish NFFFs in describing the corona, it is imperative to assess the credibility of the obtained magnetic field in reproducing the observed evolution of ARs. Specifically, in this paper, we perform numerical simulations with the non-axisymmetric NFFF as the initial condition and compare the evolution with the recorded dynamics of AR 11283. Importantly, unlike contemporary works that require a prescribed flow to induce dynamics (Amari et al., 2003; Aulanier et al., 2010), here the inherent non-zero Lorentz force serves that purpose.

To incorporate possible spontaneous reconnections (Kumar et al., 2016), the simulations are performed in congruence with the the magnetostatic theorem. For this purpose, an incompressible thermally homogeneous magnetofluid with infinite electrical conductivity is allowed to relax to a quasi-steady state via viscous relaxation (Bhattacharyya et al., 2010; Kumar et al., 2014, 2015b) satisfying

[TABLE]

in usual notations. The equations (1)-(4) are written in dimensionless form, with the following normalization

[TABLE]

The constants $B_{0}$ and $L_{0}$ are generally arbitrary, but can be fixed by the magnetic field strength and size of the system. Also, $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}$ have dimensions of time, and represent the Alfvén transit time ($\tau_{a}=L_{0}/v_{a}$) and viscous diffusion time scale ($\tau_{\nu}=L_{0}^{2}/\nu$), respectively, with $\nu$ being the kinematic viscosity. The ratio $\tau_{a}/\tau_{\nu}$ can be interpreted as an effective viscosity, which, along with other forces, determines the dynamics. The pressure perturbation, $p$, about a thermodynamically uniform ambient state satisfies the elliptic boundary value problem, which is generated by imposing the discretized incompressibility constraint (equation (2)) on the discrete integral form of the momentum equation (equation (1)); cf.(Bhattacharyya et al., 2010) and the references therein. An identical procedure involving the gradient of an auxiliary potential in the induction equation (equation (3)) is employed to keep ${\mathbf{B}}$ solenoidal, see Ghizaru et al. (2010) and Smolarkiewicz & Charbonneau (2013) for details.

If released from non-equilibrium, the initial Lorentz force will push the plasma and generate mass motion which is dissipated by the viscosity. With the flux freezing condition satisfied, the magnetic field cannot decay out completely but the terminal kinetic energy is expected to be zero because of viscous dissipation. In the corresponding numerical simulation scales inevitably become under-resolved and numerical artifacts such as spurious oscillations are generated through the employed numerical techniques. These underresolved scales can be removed by utilizing an appropriate numerical diffusivity of the non-oscillatory finite difference algorithm. In the literature, such calculations relying on the non-oscillatory numerical diffusivity are referred to as the Implicit Large Eddy Simulations (ILESs), (Grinstein et al., 2007). The computations performed in this work are in the spirit of ILESs and the simulated MRs, being intermittent in space and time, mimic physical reconnections in high-$R_{M}$ fluids.

The rest of the paper is organized as follows. Section 2 highlights the important steps in constructing the non-axisymmetric non-force-free magnetic field. In Section 3, the numerical model is briefly discussed. The results of the simulation are presented in Section 4. Finally, in Section 5 we summarize the results and highlight key findings.

2 Initial magnetic field

The initial NFFF $\mathbf{B}$ is generated by superposing two linear force-free fields, ${\mathbf{B}}^{1}$ and ${\mathbf{B}}^{2}$ (Bhattacharyya et al., 2007; Kumar & Bhattacharyya, 2011):

[TABLE]

The constant $\kappa$ is the relative amplitude of the superposing fields. Being linear-force-free, the $\mathbf{B}^{1}$ and $\mathbf{B}^{2}$ satisfy (Chandrasekhar, 1956; Chandrasekhar & Kendall, 1957)

[TABLE]

where $i=1,2$. The constant $\alpha_{i}$ represents the torsion coefficient (Parker, 2012) for the field ${\mathbf{B}}^{i}$ and $\alpha_{2}=(1+\gamma)\alpha_{1}$ with $\gamma$ measuring the deviation of $\mathbf{B}$ from being linear-force-free. A convenient expression of $\mathbf{B}^{i}$ in spherical polar coordinates is

[TABLE]

where ${\bf{r}}$ is the radial vector. The scalar eigenfunction $\psi^{i}$ satisfies the Helmholtz equation $(\nabla^{2}+\alpha_{i}^{2})\psi^{i}=0$ and can be expressed as

[TABLE]

where $j_{l}(r,\theta,\phi)$ and $y_{l}(r,\theta,\phi)$ are the spherical Bessel functions of the first and second kinds respectively and $P_{l}^{m}$ is an associated Legendre polynomial of degree $l$ and order $m$; $C_{lm}^{1}$ and $C_{lm}^{2}$ are constants. The individual components of $\mathbf{B}^{i}$ are

[TABLE]

Further, to make $B^{i}$ physically meaningful without having infinite energies at $r=\infty$ (Seehafer, 1978), we bound it above by confining the MFLs in a spherical shell of radius $r_{0}\leq r\leq r_{1}$ (Prasad & Bhattacharyya, 2016). Consequently,

[TABLE]

where $r_{0}$ is the inner radius of the shell selected to avoid singularity at the origin.

To simulate the solar AR field using $\mathbf{B}$, we follow the mathematical steps described in Prasad et al. (2014), which provide fast and reasonably good fits to vector magnetograms. The procedure involves generating of two-dimensional vector magnetogram templates from three-dimensional analytical solutions computed over the volume spanned by the spherical shell. To achieve this, a plane tangential to the shell at a location $r^{\prime}$ is constructed, and the components of $\mathbf{B}$ are computed on this plane. The plane is identified by a specific AR cutout whose orientation can further be changed by varying the two Euler angles $\theta^{\prime}$ and $\psi^{\prime}$; see Figure 1. Finally, the best correlated $\mathbf{B}$ at the cutout is obtained by maximizing

[TABLE]

which correlates $\mathbf{B}$ with the field $\mathbf{B}_{O}$ obtained from the magnetogram. Explicitly, the parameters $\alpha_{1}$, $l$, $m$, $\kappa$ and $\gamma$, and the variables $r_{0}$, $r_{1}$, $r^{\prime}$, $\theta^{\prime}$, $\psi^{\prime}$ are obtained by maximizing $c$. Using this technique, Prasad et al. (2014) analyzed vector magnetograms for AR 10930 spanning a period of three days with two X-class flares. Their estimates for the changes in free energy and relative helicity of the AR before and after the flares were found to be in good agreement with those obtained from other NLFFF extrapolations. In this work, we use vector magnetograms from the Heliospheric Magnetic Imager (HMI)(Schou et al., 2012) on board the Solar Dynamics Observatory (SDO)(Pesnell et al., 2012). The HMI (hmi.sharp_cea_720s data series) provides full-disk vector magnetograms of the Sun with a spatial resolution of $0^{\prime\prime}.5$ and a temporal cadence of 12 minutes. The magnetogram $\mathbf{B}_{O}$ is initially remapped to a Lambert cylindrical equal-area projection and then transformed into heliographic coordinates (Gary & Hagyard, 1990). The selected vector magnetogram is for the AR 11283 on 2011 September 6 at 21:00 (Figure 2 a), the same AR as discussed in Prasad & Bhattacharyya (2016) but observed at an earlier time. The AR produced an X- class flare on 2011 September 6, peaked around 22:20 hr, and was followed by a coronal mass ejection (Feng et al., 2013; Janvier et al., 2016). The initial magnetic field is constructed by using the magnetogram data at 21:00 hr, which documents that the magnetic field is morphologically complex with a strong bipolar region and other weaker diffuse field regions (Figure 2).

The details of the parameters for the best correlated non-axisymmetric NFFF are given in Table 1 and the corresponding MFLs (top and side views) are shown in Figure 2 b and Figure 3.

The figures are overlaid with the contours of $\mathbf{B}_{z}$ on the $z=0$ plane, which confirm the overall morphological similarity of the constructed field with the magnetogram in terms of having strong bipolar and weaker diffuse field regions. The MFLs connecting the strong polarity regions are twisted loops, having the appearance of a forward S, and hence resembles sigmoid emission in AR 11283 seen by the Atmospheric Imaging Assembly (AIA) 94 Å  channel (Jiang et al., 2013). The field ${\mathbf{B}}$ has two-dimensional X-type nulls (Figure 4) constituting the two null lines, marked by X1 and X2. The location of the nulls is determined using the technique reported in Kumar et al. (2015b). To further analyze the initial field line geometry, in Figure 5(a), we identify four sets of twisted MFLs (colored red, yellow, green and gray), structurally resembling a magnetic flux rope, which is a stack of magnetic flux surfaces (MFSs) spanned by twisted field lines (Schrijver & Siscoe, 2009). To closely examine the resemblance, Figure 5(b) plots a group of MFLs (in color red) with circularly distributed initial footpoints. The circle is marked by the curve P with its center located at $(x,y,z)=(0.23,0.66,0)$. Notably, the figure confirms the MFLs to be to be twisted by approximately a half turn (angle $\pi$) about the central axis. Such visual confirmation of twist in field lines has also been done in Jiang et al. (2014). After tracing the loop trajectories in the computational domain, the locus of the endpoints of the MFLs generate a second closed curve P1 on the $z=0$ plane. Since field lines cannot intersect each other, the group of MFLs generating the two closed curves must have a single common axis, which is identified as the central axis of the rope. Further, MFLs join the closed curves and are tangential to a surface, making it a MFS (Bhattacharyya et al., 2010). By changing the radius of the circle, a co-axial stack of MFSs is generated (Figure 5(c)) which, within the accuracy of the field line integration scheme, is identified as a magnetic flux rope. Further analysis (not reported) finds a threshold for the radius of the circularly distributed initial points, above which the loci of the endpoints are no longer closed curves and do not give rise to a flux rope.

3 Numerical model

To solve the MHD equations (1)-(4), 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). 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). Importantly, the MPDATA has a proven dissipative property, which is intermittent and adaptive, allowing the formation of underresolved scales in field variables for a selected grid resolution. In the magnetofluid undergoing viscous relaxation, sharpening of the magnetic field gradient is expected to be unbounded, ultimately leading to MRs at locations where the separation of non-parallel field lines approaches the grid resolution. The MR process per se is underresolved, but effectively regularized by the local second-order residual dissipation of MPDATA, which is sufficient to sustain the monotonic nature of the solution. Being intermittent and adaptive, the residual dissipation facilitates the model performing ILESs that mimic the action of explicit subgrid-scale turbulence models (Margolin et al., 2006) whenever the concerned advective field is underresolved. 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., 2016, 2015a). The simulations presented continue to rely on the effectiveness of ILES in regularizing the onset of MRs.

4 Simulation results

The simulations are performed on a $128\times 128\times 128$ grid, resolving the normalized computational domain $\Gamma$ spanning $[0,1]\times[0,1]\times[0,0.25]$ (respectively in $x$, $y$, and $z$); see Figure 1, where a unit length is equivalent to 50 Mm. The vertical extent of the domain is constrained by the condition that the magnetic field is to be contained within the flux surface at $r_{1}$ in Figure 1. The simulations start from a motionless state with the NFFF obtained in Section 2. To ensure net zero magnetic flux through $\Gamma$, the volume magnetic field $\mathbf{B}$ is continued to the domain boundaries (Kumar et al., 2015a). The dimensionless constant $\tau_{a}/\tau_{\nu}\approx 7\times 10^{-3}$, which is two orders of magnitude larger than its coronal value. To calculate ${\tau_{a}}$, we set $\rho_{0}=1$. Notably, the larger value of $\tau_{a}/\tau_{\nu}$ only speeds up the computation because of a more efficient viscous dissipation by the increased effective drag. With time being normalized by $\tau_{a}$, a unit of the Alfvén transit time can be interpreted as $\approx 15$ minutes of HMI data. With a constant mass density, flow in the computation is incompressible—an assumption also used in other works (Dahlburg et al., 1991; Aulanier et al., 2005). Although compressibility is important for the thermodynamics of coronal loops (Ruderman & Roberts, 2002), our focus is on elucidating the changes in magnetic topology idealized with a thermally homogeneous magnetofluid. Throughout the simulation, the maximum plasma $\beta$ varies in the range $\{0.003,0.007\}$ and the fluid Reynolds number $R_{F}\approx 10^{2}$, which reasonably approximate the coronal values. The non-force-free nature of the initial field can be assessed from the histogram and the direct volume rendering of $|\mathbf{J}\times\mathbf{B}|/|\nabla(B^{2}/2)|$ for the initial field (Jiang & Feng, 2012) as shown in Figure 6, where we find most of the values populated in the range $\{0.005,0.01\}$. Additionally, with zero resistivity $R_{M}=\infty$ and the flux freezing is satisfied throughout evolution, except during reconnections. Consequently, the MFLs evolve in response to unbalanced forces (including Lorentz force) pushing the plasma. The plasma motion then naturally incorporates line-tied boundaries, which is conventional in coronal simulations (Schrijver & Siscoe, 2009). In addition, coronal simulations may require more sophisticated boundary treatment to replicate the observed dynamics (Jiang et al., 2013). An evolution of field lines by forcing the plasma inevitably makes the footpoints to move appreciably. A recent work concurs by suggesting the rapid evolution of photospheric magnetic field during major eruptive events (Sun et al., 2017). Nevertheless, we have also performed simulations where the bottom boundary of the domain is kept fixed (not shown). A different flux rope, co-located with the gray colored rope (Figure 5a) and having similar connectivity, replicates the reported evolution.

Figures 7 and 8 show the top and side views of a sequence of snapshots from the evolution of the overlying MFLs to the flux rope colored red in Figure 5(a). Until reconnection at $t=2.2$, the rope is located mostly above the main polarity inversion line (PIL) which separates the two strong opposite polarity regions R1 and R2. The figure is further overlaid with the X-type null lines X1 and X2 (in gray). Evidently, the initial MFLs of the rope are twisted while being anchored at the $z=0$ plane. The sequence indicates an overall anti-clockwise rotation of the overlying field lines, the rope, and the strong polarity regions. The flux rope also expands along with preferential reconnections at the null line X1. To assess the effect of rotation on the preferential reconnection, we plot streamlines at $t=0.4$ near X1 and X2, as shown in Figure 9. As evident from the figure, the streamlines favor a clumping of oppositely directed field lines in the neighborhood of X2, leading to reconnections. However, near X1, the streamlines bring field lines with the same direction toward the null line and hence no reconnection occurs. The net result is an asymmetric evolution that favors reconnection at X2 rather than X1. Being situated closer to the X2, the MFLs overlying the rope reconnects, first resulting in a local decrease of the magnetic tension force.

In response, the flux rope rises vertically while simultaneously getting rotated toward X2. Figure 10 shows the evolution of the maximum height of the innermost flux surface of the rope (marked in black) which confirms that the rise is monotonic. Furthermore, the rise is not symmetric above the two legs of the rope due to the asymmetry of the initial field and the preferential reconnection. Notably, such a result is essentially similar to a more realistic MHD simulation initialized with an NLFFF performed by Jiang et al. (2013), who show that the expansion of a low-lying flux rope is assisted by the reconnection of its overlying flux. Furthermore, MFLs on the outermost surface of the rope (depicted in color red) reconnect, which partially opens and bifurcates the rope, making it connect region R1 with regions R2 and R3 (Figure 7). Afterwards, the inner surfaces of the rope (marked in yellow and black) also bifurcate in a similar way, leading to the breakup of the rope. The concept of flux-rope bifurcation is well-established and responsible for the possible drainage of mass in an ascending filament (Gibson & Fan, 2006, 2008). Importantly, the evolution of the AR 11283 also documents two-pronged pitchfork-type bifurcations of the associated filament when observed in the 304Å  channel of the AIA (Lemen et al., 2012) at approximately UT 21:50 (Figure 11). The reasonably faithful simulation of the observed flux rope bifurcation justifies the credibility of the non-axisymmetric non-force-free magnetic field in exploring the coronal dynamics. However, without reproducing the bursty rise of the rope, the present computations fail to capture the subsequent X-class flare. The failure can be attributed to either the exclusion of the flux emergence occurring from 22:00 hr onward (Petrie, 2012) in the simulated dynamics, or the requirement of a more accurate method to reconstruct the initial NFFF. Notably, the vertical component of the emergent field, when favorably aligned with respect to the existing field, may reconnect and provide the trigger necessary for the impulsive rise.

An auxiliary analysis (not shown) indicates an increase in the magnetic flux of MFLs overlying the ropes, in colors white and green in Figure 5, because of reconnections at X2. The increased density results in an enhanced magnetic tension force which tethers the white and green ropes further and inhibit their ascent, whereas the yellow rope situated away from the main PIL shows an appreciable but non-monotonic rise because of similar tethering by reconnected field lines.

5 Summary and conclusions

This paper numerically simulates the viscous relaxation of a thermally homogeneous incompressible magnetofluid with infinite electrical conductivity and having an initial non-axisymmetric non-force-free magnetic field. The initial field is constructed by forward fitting a photospheric magnetogram using an analytically non-force-free magnetic field. In our calculations, the focus is on AR 11283, with magnetogram from the HMI. Notably, as the initial field is non-force-free and non-axisymmetric, it is able to characterize the basic structural features of such a complex AR. The corresponding magnetic field consists of twisted field lines and two X-type null lines. Most importantly, there exists a sigmoidal flux rope that roughly matches the observed sigmoid emission. A direct visualization confirms a rope to be a stack of local MFSs.

As the initial magnetic field is not balanced, the Lorentz force naturally trigger evolution. Subsequently, the favorable forces bring non-parallel field lines in close proximity of the null lines and lead to an unbounded sharpening of the magnetic field gradient, which results in the development of underresolved scales. In the spirit of ILES, the MPDATA then generates locally adaptive dissipation to regularize the underresolved scales with simulated MRs.

The computation reveals, at one of the X-type null lines, that the initial flow favorably clumps non-parallel MFLs and triggers MRs. The generated outflow makes the two major polarities rotate, which brings further non-parallel field lines toward the X-line and drives secondary reconnections. In response, the magnetic tension force at regions overlying the rope decreases and with the reconnections being preferential, results in an asymmetric ascent of the rope. Although occurring in congruence with the magnetostatic theorem, the evolution leads to reconnections occurring only at the pre-existing X-type nulls. Importantly, the simulation depicts a sequential bifurcation of the rope over a range of time. The images from the AIA 304 Å  channel suggest a similar bifurcation of a filament, which adds credibility to the use of non-axisymmetric NFFF to represent a complex AR. However, the computations fail to capture the subsequent X-class flare as the required impulsive ascent of the rope is missing. We attribute the failure to not accounting for the flux emergence and other solar surface modifications between the bifurcation and the flares. For example, the recent data-driven MHD simulation of the same AR by Jiang et al. (2016), which includes not only the flux emergence but all other surface motions, accurately replicates the M-class flare at 01:50 hr.

In retrospect, in its basic utilization, the non-force-free non-axisyemmteric magnetic field is useful in modeling ARs and the evolution of the solar corona. The future goal is to develop a complete numerical scheme that will simulate the dynamics of an AR from such an initial state to the occurrence of flares, if any. In the backdrop of the successful simulation of coronal dynamics by Jiang et al. (2016), we believe that a synergetic approach exploiting the advantages of various numerical models is the key to achieve this goal. For instance, the fully numerical non-force-free extrapolation model developed by Hu & Dasgupta (2008) and its subsequent developments (Hu et al., 2008; Gary, 2009), rely on self-organized higher-curl Beltrami fields (Bhattacharyya et al., 2007; Kumar & Bhattacharyya, 2011; Janaki & Dasgupta, 2012) that are morphologically identical to the initial magnetic field $\mathbf{B}$. Therefore, it can be expected that a data-driven MHD simulation of reconnections with the initial field provided by non-force-free extrapolations will provide a more realistic dynamics of the coronal magnetic field and thus warrants further attention.

Acknowledgements: The simulations are performed using the 100 TF cluster Vikram-100 at the Physical Research Laboratory, India. We acknowledge the use of the visualization software VAPOR (www.vapor.ucar.edu) for generating relevant graphics. Data and images are courtesy of NASA/SDO and the HMI and AIA science teams. SDO/HMI is a joint effort of many teams and individuals to whom we are greatly indebted for providing the data. The authors are thankful to Dr. P. K. Smolarkiewicz for his invaluable contribution toward the better presentation of this paper. We are also thankful to the anonymous referee for providing insightful suggestions toward the overall betterment of this paper.