Simulating the formation and eruption of flux rope by magneto-friction model driven by time-dependent electric fields

We simulate the coronal magnetic field evolution by magneto-friction (MF) method (Yang et al., 1986), where the magnetic field evolves in response to the photospheric flux motions through the non-ideal induction equation given by

where $\mathbf{v}_{\rm MF}$ is the magnetic frictional plasma velocity, $\mathbf{A}$ is the vector potential relating the magnetic field as $\mathbf{B}=\nabla\times\mathbf{A}$. The second term is to accommodate dissipation in the corona due to electric currents $\mathbf{J}=\nabla\times\mathbf{B}/\mu_{0}$. The magnetic diffusivity can be chosen a typical value of $\eta=1\times 10^{8}\,m^{2}/s$ to be able to run the simulation stable. Here the induction equation is solved in terms of $\mathbf{A}$ to ensure that $\nabla\cdot\mathbf{B}=0$ without additional divergence-free schemes. From the assumption of static magnetic fields, the MF velocity is given by

where $\nu$ is the magneto-frictional coefficient that controls the speed of the relaxation process and $\mu_{0}$ is the magnetic permeability in vacuum. As suggested in Cheung & DeRosa (2012) a height-dependent form of the frictional coefficient $\nu$ is given by

where $\nu_{0}$ is set around $35\times 10^{-12}$ s m^-2 and $z$ is the height above the bottom boundary, and $L$ is the chosen as 15 Mm. This form of frictional coefficient gives MF velocities a smooth transition to zero towards bottom boundary $z=0$.

Along with a special driver module, Vemareddy et al. (2024) first implemented the MF model in Pencil Code (Pencil Code Collaboration et al., 2021) and tested the evolution of coronal magnetic field driven by observed photospheric magnetic fields. The Pencil Code is a highly modular physics-oriented simulation code that can be adapted to a wide range of applications. It is a finite-difference code using sixth order in space and third-order in time differentiation schemes.

We use potential field (PF) model (Gary, 1989)) as the initial field for the simulation. The PF is constructed from the normal component ($B_{n}\to B_{z}$) of the observed vector magnetic field of the AR 11429 at 06:00 UT on March 2012. As the Pencil Code operates on the vector potential ($\mathbf{A}$) rather than magnetic field, the field divergence is being satisfied without the need for extra schemes. The vector potential $\mathbf{A}_{p}$ of the PF is computed by imposing coulomb gauge condition ($\nabla\cdot\mathbf{A}_{p}=0$) and has vanishing normal component ($\mathbf{A}_{p}\cdot\hat{n}=0$) at the bottom boundary (e.g., DeVore & Antiochos 2000). With these constraints, the components of $\mathbf{A}_{p}$ are computed as,

where $k_{x}$ and $k_{y}$ are wave vectors along the x and y directions, respectively, and $k=\sqrt{k_{x}^{2}+k_{y}^{2}}$. Here, FT refers to the 2D Fourier transform operation on the observed $B_{z}(x,y)$ and $FT^{-1}$ refers to inverse Fourier transform. The observed $B_{z}$ is adjusted to satisfy flux balance condition and is inserted in a larger area by padding the dimensions. The $\mathbf{A}_{p}$ is constructed on a uniform Cartesian computational grid of $192\times 192\times 120$ representing the AR coronal volume of physical dimensions $280\times 280\times 175$ Mm³.

In the data-driven MF models, the initial magnetic field at the start time of the simulation is driven by magnetic field observations of the photospheric boundary that represents bottom of the computational domain (Mackay et al., 2011; Yardley et al., 2018). We obtained the ambiguity resolved HMI vector magnetic field components which are served as hmi.sharp.cea_720s data product in spherical coordinates. These are approximated to Cartesian components as $B_{x}\to B_{\phi}$, B${}_{y}\to-B_{\theta}$ and $B_{z}\to B_{r}$ (Bobra et al., 2014) and are smoothed with a Gaussian width of 3 pixels both spatially and temporally. Also, the flux balance condition is applied to the $B_{z}$ component. In order to reduce computation burden, the field components are rebinned by a factor four such that each pixel corresponds 2” on the computation grid. Driving the initial 3D field with photospheric magnetic fields leads to mild to moderate non-potential fields, which may not represent the twisted core field of the ARs (Yardley et al., 2018). For complex active regions such as 11429, the observations have insufficient influx of magnetic energy and helicity, as a result the formation of low lying twisted flux in the simulations has not been possible (Vemareddy et al., 2024). In this work, the initial field is driven by the electric fields ($\mathbf{E}(x,y,z=0,t)$), which are derived from time sequence vector magnetic field observations of the AR (Fisher et al., 2010). The electric field has two components: viz inductive field ($\mathbf{E}_{I}$) and non-inductive component defined by gradient of a scalar $\psi$

The inductive component is obtained by solving the Faradays’ law of induction equation, $\partial_{t}\mathbf{B}=-\nabla\times\mathbf{E}_{I}$, which involves decomposing the poloidal-toroidal decomposition of the time derivative of the vector magnetic fields (Fisher et al., 2010). The non-inductive component requires a functional form for $\psi$ which is unknown and therefore invokes the following three different assumptions (Cheung & DeRosa, 2012; Cheung et al., 2015):

where $\Omega$ and $U$ are the free scalar parameters referring to the extent to which the axisymmetric vertical flux tube rotates and emerges, respectively. The first assumption is equivalent to zero contribution of non-inductive components, whereas the latter two cases are suggested by Cheung & DeRosa (2012); Cheung et al. (2015) in order to ensure sufficient injection of magnetic helicity and free magnetic energy. The above chosen functional forms are linked to photospheric magnetic field observations and their horizontal spatial derivatives, therefore the non-inductive part of $\mathbf{E}$ has horizontal components (x and y) that are being added to inductive part of $\mathbf{E}$. This form for non-inductive $\mathbf{E}$ have vanishing horizontal divergence, as a result satisfies the induction equation such that $\partial B_{z}/\partial t=-\hat{z}\cdot(\nabla\times\mathbf{E}_{h})$ without modifying the $B_{z}$.

The derived photospheric electric fields at an interval of 12 minutes, are being used as driver field $\mathbf{E}$ at $z=0$, instead of magnetic fields or vector potentials. At each instant of time, we add $\mathbf{E}$ at the bottom layer of the vector potential $\mathbf{A}$ in the computational domain and the time evolution of the magnetic field is determined by the time-integration of the Faraday’s law

In the Pencil Code , we developed a driver module to read the time-series of $\mathbf{E}$ at 12 minute intervals and interpolate them at the simulation time-step.

The lateral sides of simulation domain are specified with periodic boundary conditions and the top one is set to open boundary condition.

We drive the initial PF with the time-dependent electric field derived from the ad hoc assumptions given in equations 7-9. With the electric fields from the three assumptions, we perform essentially three simulation runs. The run R1 is based on the inductive electric field ($E_{I}$), whereas runs R2 and R3 are included with non-inductive contributions in order ensure sufficient injection of energy and helicity fluxes. For the later two runs, one has to specify appropriate values for free parameters $U$ and $\Omega$. As suggested by Pomoell et al. (2019), these values are constrained by comparing the Poynting flux deduced from the DAVE4VM vector velocity and magnetic field (Figure 1e, and equation 9) with that of the derived electric field ($dE/dt=\int\mathbf{E}\times\mathbf{B}\cdot d\mathbf{A}$) from these assumptions (Pomoell et al., 2019).

We use the values of $U=120$ m s^-1 and $\Omega=0.096$ turns/day for R2 and R3, respectively. As shown in Figure 2, these values give a time-integrated poynting flux of $18\times 10^{32}$ ergs till the eruption time (vertical dotted line), which is roughly a factor of two higher when compared with the observed one (Figure 1e). This difference is admitted for the reasons that i) the AR could be in a non-potential state even after the first eruption, ii) the Poynting flux derived from DAVE4VM may be underestimated because the velocities represent averaged values over an apodising window typically $19\times 19$ pixel², iii) observational sensitivity of magnetic field measurements. With these existing difficulties, the chosen values of U and $\Omega$ are broadly constrained with an uncertainty of upto 40%. It should be noted that the coronal free energy is small (about 18 times) and the formation of the twisted flux is quite implausible with $\mathbf{E}_{I}$ alone, therefore higher values of these free parameters are deemed in order to pump sufficient free magnetic energy to generate magnetic structure that is comparable to observations

The above said three types of MF simulations are performed by driving the initial PF with a time dependent bottom boundary electric field for 72 hours duration starting from 6:00 UT on March 7, 2012, and further 24 hours without changing the lower boundary.

As explained earlier, the $\mathbf{E}$-field satisfies $\partial B_{z}/\partial t=-\hat{z}\cdot(\nabla\times\mathbf{E}_{h})$ reproducing the observed evolution of $B_{z}$ self-consistently during the simulations. To check this , in Figure 3, we compare observed and simulated $B_{z}$ at $z=0$ plane at 10h and 47h time instants of simulation R2. As can be noticed that simulated $B_{z}$ has a high degree of spatial correlation including small magnetic elements in the vicinity of sunspot regions. The scatter plots refer to a correlation of $>99.95\%$ which implies the observed evolution of photospheric magnetic fields during the simulation. However, such a correlation is not true for horizontal components of $\mathbf{B}$ which are scaled by the free-parameters in non-inductive contribution of $\mathbf{E}$ in order to inject additional free-energy and helicity.

Figure 4 presents the snapshots of the simulated magnetic structure of all three runs. Field lines are traced from the foot point locations where the horizontal magnetic field and total electric current are strong. From R1, one can see that the initial potential field does not change significantly over several hours of evolution as the boundary electric field has inductive components derived from $B_{z}$ components alone. In the runs R2 and R3, the initial magnetic structure evolved progressively to a highly twisted structure at the core of the AR surrounded by a less sheared field which is evidently noticed in the snapshots at 25th hour. Further evolution till 45th hour, leads to transforming neighbouring less sheared fields into twisted fields akin to a coherent flux rope along the PIL. Following this, the flux rope and the surrounding field grow and appear to rise in height, which is a typical signature of initiation of the eruption. In the presence of magnetic diffusion and supply of magnetic helicity through a time-dependent bottom boundary, the MF relaxation method thus captures the formation and eruption of MFR in the AR. Using the other codes, earlier studies (eg., Cheung & DeRosa 2012; Pomoell et al. 2019) reported the simulated coronal field evolution capturing the flux rope formation in the emerging ARs.

It is worth pointing out that the electric fields in R2 are derived with $J_{z}$ distribution at the photosphere, as a result, one can expect a spatially varying twist about the PIL scaled by the free parameter $U$. And in R3, the non-inductive contributions to the derived electric fields are based on $B_{z}$ distribution scaled by the $\Omega$, thus the spatial variation of magnetic twist along the PIL and its vicinity do not reflect the same observations as those of electric fields in R2. Given this fact, the simulated 3D magnetic structure in R2 better reflects the coronal EUV observations than that in R3, as will be convinced from the following comparative analysis with the observations.

The time evolution of computed relative magnetic helicity (H) and magnetic energy in the simulated domain are plotted in Figure 5. In the MF relaxation of R1, the volume helicity (H) increases in magnitude from $-4\times 10^{-9}$ Mx² to $-0.37\times 10^{42}$ Mx². This is solely due to shear motions of fluxes along PIL. Although there is no appreciable increase of total magnetic energy (E), the free magnetic energy ($E-E_{p}$) increases continuously, and it is $1.31\times 10^{32}$ ergs at the time of the eruption. During the simulation time of R1, the fractional free energy increases upto 10%.

On the other hand, the H accumulates in the corona continuously (blue, red curves) and amounts to $-39.2\times 10^{42}$ Mx², $-81.6\times 10^{42}$ Mx² for R2 and R3 respectively. Similarly, the total and free magnetic energies exhibit increasing behaviour as the coronal magnetic field becomes twisted. In these runs, the fractional free energy corresponds to 47% (R2), 57% (R3) at the time of the eruption, which is probably significant enough to initiate the rise motion of the flux rope. Numerical simulations by Amari et al. (2003) reported a fractional free energy of 5% to initiate the eruption, while TMF simulations by Pomoell et al. (2019) refers to 14-50% of fractional free energy that resulted in erupting flux rope, therefore, our results are more aligned with later one as the simulation setup is same. It is important to note that the coronal estimates of helicity ($H$) for R3 is higher by a factor of two than that for R2, even though the same amounts of energy fluxes are injected for these runs. As discussed earlier, the higher value of $H$ is due to uniform distribution of twist nature arising from equation 8, correspondingly, the coronal field is highly twisted in R3.

We performed another simulation of R2 driven by electric fields with free parameter $U=180\,m/s$. This run is to understand the flux rope formation time scale and rise motion during a given time evolution. From the comparison, the twisted flux formation in this run happened early and later its rise motion to higher heights. Figure 6a displays the rendered magnetic structure of R2 ($U=180\,m/s$) with $B_{z}$ as background image. The modelled magnetic structure consisted of a low-lying, twisted core field along the PIL overlaid by potential arcades and J-shaped field lines (lobes or elbow field lines) sheared past each other in SW and NE regions, which together manifest an inverse-S sigmoidal structure. Although the global AR magnetic configuration is inverse S-shaped, the magnetic field in the SW region alone mimics another inverse-S sigmoid. As shown in Figure 6(b-c), the magnetic structure qualitatively resembles the morphology of coronal features captured in AIA images. Owing to high temperatures, sigmoids are often seen in hot EUV channels. AIA 335 Å image ($T\approx 2.5$MK) displays coronal sigmoid resembling the simulated magnetic structure and plasma loops as the tracers of magnetic field in AIA 171 Åsnapshot. The AIA 304 Å waveband is a cool channel ($T\approx 5\times 10^{4}$K) image exhibiting a filament channel in the SW region, which is plasma embedded in the dips of twisted field of the sigmoid configuration.

The simulated magnetic structure is rendered in perspective, as depicted in Figure 7 at different times. The initial potential field becomes sheared progressively forming inverse J-shaped field lines surrounding the low lying twisted core field along the PIL. In a time scale of two days a well developed twisted flux builds up along the PIL. Thereafter the twisted magnetic structure ascends in height with time, especially the field in the SW region, from 40 Mm to 120 Mm (See the attached animation). This rise motion corroborates the observations of the eruption from the SW region of the AR (Dhakal et al., 2020). We emphasise that the rise motion does not turns into an eventual eruption like the observed onset of eruptions, rather it is a progressive in time even after driving stopped from 72 hours. Since the velocity is controlled by $\nu$ parameter, the TMF simulations implicitly lacks the short term dynamic evolution, although they capture the build up of free magnetic energy during slow evolution over long time.

Since the MF model does not include the thermodynamic evolution of the plasma, we cannot synthesize images of coronal emission in EUV or X-ray wavelengths. To produce synthetic maps of coronal loops similar to EUV images, Cheung & DeRosa (2012) introduced a method to generate proxy emissivities based on the modelled magnetic field alone. Proxy emissivities are derived with the values of square of the total current density ($J^{2}$) averaged along a magnetic field line. This technique is useful for a qualitative visualization of coronal loops in a magnetic model such as MF, however it is not a replacement for more sophisticated techniques that use the thermodynamic variables from MHD models (e.g., Lionello et al. 2009; Hansteen et al. 2010). Following their method, we computed the emissivity ($\varepsilon(x,y,z)$) at each grid point in the computational domain. This volume distribution of emissivity is integrated in the vertical direction to generate a 2D distribution, as seen in coronal EUV observations. Such proxy emission maps at two epochs are compared with the respective AIA 304, 335 Å images in Figure 8. We can notice a striking morphological similarity of the modelled emission with the diffuse plasma emission of the sigmoid captured in 335 Åimage and trace of dark filament present in 304 Å images. This filament is regarded as flux rope in the models and it is expanding and rising in time as observed in the vertical cross section planes.

For a topological study, we analyse the quasi separatrix layers (QSLs, Titov et al. 2002) where the magnetic field line linkage changes drastically in the volume. The strength of QSLs is measured by squashing factor Q. From the simulated magnetic structure, we computed Q with the procedures of Liu et al. (2016) and its map is displayed in Figure 9a. In this map, QSLs with large Q values ($>10^{6}$) are identified by intense white traces in strong field regions. Continuous trace of Q can be noticed along the main PIL.

To compare QSL locations with the flare ribbon emission, the AIA 1600 Å,  304Å, 131 Å observations during the impulsive phase of the flare are displayed in the panels of Figure 9(b-d) respectively. On these images, contours of Q at $10^{5}$, $10^{6}$ levels are over-plotted by removing the irrelevant QSLs with a mask on the Q-map. The ribbon emission is related to the erupting MFR in the SW region, which is co-spatial with the QSL section of the simulated magnetic configuration. Unlike the QSLs determined in the NLFFF models (Vemareddy, 2021), the QSLs in this simulation better represent the twisted core flux along the PIL and the observed ribbon emission.

From the theoretical studies, it was predicted that the flare ribbons are the photospheric/chromospheric foot prints of QSLs that encloses a twisted FR (Demoulin et al., 1996). The extremities of the ribbons are found to be hook shaped for weakly twisted FRs and are spiral shaped for highly twisted FRs (Zhao et al., 2016). The observed ribbon morphology in our case delineates inverse-S shape with co-spatial hooks in the extreme ends, therefore the enclosed flux rope is indicated to be moderately twisted. This finding is consistent with previously reported pre-eruptive NLFFF configurations (Bobra et al., 2008; Su & van Ballegooijen, 2013; Liu et al., 2016; Vemareddy, 2021).

In Figure 10, we display Q and twist number ($tw$) maps computed from the vertical cross section of the erupting flux rope at three different times. In these maps, the QSLs of large Q values well distinguish two closed domains belonging to twisted core sections above the PIL and the surrounding less sheared arcade. With the twisted flux (flux rope) at the core, the QSLs in the cross section resembles an inverse tear-drop shape as evidenced by the theoretical/numerical models (Jiang et al., 2018; Vemareddy, 2019). The twist of the field lines at the core is negative; therefore, field lines are left helical and the flux rope has an inverse-S morphology. Within the flux rope cross-section, the total twist of the field lines varies in the range of 2 turns; however, average twist number is well below 1 turn, so that the twisted flux is in stable equilibrium. After the formation in the first few hours, the twisted flux rope expands and rises up in time, as delineated in Figure 10.