An Evolution and Eruption of the Coronal Magnetic Field through a Data-Driven MHD Simulation

We solve the following zero-$\beta$ MHD equation to make the ground-truth data and also conduct the data-driven simulation,

where $B$ is the magnetic flux density, $v$ is the velocity, $J$ is the electric current density, $\rho$ is the pseudo density, and $\phi$ is the convenient potential to remove errors derived from $\mbox{\boldmath$\nabla$}\cdot\mbox{\boldmath$B$}$ (Dedner et al. 2002), respectively. $\rho_{0}$ corresponds to the initial density. The length, magnetic field, density, velocity, time, and electric current density are normalized by $L^{*}$ , $B^{*}$, $\rho^{*}$ , $V_{\rm A}^{*}\equiv B^{*}/(\mu_{0}\rho^{*})^{1/2}$, where $\mu_{0}$ is the magnetic permeability, $\tau_{\rm A}^{*}\equiv L^{*}/V_{\rm A}^{*}$, and $J^{*}=B^{*}/\mu_{0}L^{*}$. $\nu$ and $\eta$ are viscosity and resistivity, fixed by $1.0\times 10^{-3}$ and $1.0\times 10^{-5}$, respectively, and the coefficients $c_{h}^{2}$, $c_{p}^{2}$ in Eq.(6) also fix the constant value, 0.04 and 0.1, respectively. $\zeta$ is a diffusion coefficient of the density which avoids the sudden variation of the density to improve the robustness of the simulation, where $\zeta$ is set to $1.0\times 10^{-4}$ in this study. A numerical box of 1.0 $\times$ 1.0 $\times$ 1.0, which is given in its non-dimensional value, is divided by 320 $\times$ 320 $\times$ 320 grid points.

We first make ground-truth data according to an MHD simulation done by Amari et al. (2000) and Amari et al. (2003) to verify the reproducibility of the energy storage and release processes of the coronal magnetic field by our data-driven MHD simulation. We first set a simple bipole magnetic field like sunspots as shown in Fig.1a, from which the potential field is extrapolated in the 3D space following the Green function method (Sakurai 1982) as shown in Fig.1b. The initial velocity and density are set as $|\mbox{\boldmath$v$}|=0.0$ and $\rho_{0}=1.0$ in the whole space, respectively.

Next we impose a twisting motion on the photosphere according to,

where $\psi$ is a steam function which satisfies $\mbox{\boldmath$v$}_{h}^{(b)}|_{z=0}=\mbox{\boldmath$z$}\times\mbox{\boldmath$\nabla$}_{h}\psi$. We give $\psi$ the following formula,

and

where $B_{z}^{(b)}$ and $B_{z:max}^{(b)}$ correspond to the magnetic field measured at the photosphere and the maximum value, respectively. If time $t$ is beyond the critical time set at $t_{cri}$ where we set $t_{cri}=9.0$, $\gamma(t)$ quickly falls to zero. The twisting motion is shown in Fig.1c. It is convenient to give $\psi$ as a function of $B_{z}$ because the twisting motion is imposed along the contour of $B_{z}$, i.e., $(\mbox{\boldmath$v$}_{h}\cdot\mbox{\boldmath$\nabla$})B_{z}=0$ is promised regardless of $B_{z}$. Thus the distribution of $B_{z}$ will not be changed by the twisting motion. Since we assume $v_{z}$=0 at the bottom surface as well as $(\mbox{\boldmath$v$}_{h}\cdot\mbox{\boldmath$\nabla$})B_{z}=0$, these conditions satisfy $\partial_{t}B_{z}=0$. Therefore, normal flux is not transported across the bottom surface during the evolution and our simulation can keep $\int{\bf B}\cdot d{\bf S}=0$. When $\gamma(t)$ falls to zero, the twisting motion at the photosphere comes to a complete halt. Afterward, the magnetic field is relaxed for a while, .i.e., no external motion is imposed. The horizontal magnetic components $B_{x}$ and $B_{y}$ at the photosphere are following the equations during the twisting motion and relaxation process after the twisting motion is over,

while $B_{z}$ is fixed. $E_{x}$ and $E_{y}$ represent the electric field in $x$ and $y$ components which are described as $(-\mbox{\boldmath$v$}\times\mbox{\boldmath$B$})\cdot\mbox{\boldmath$x$}$ and $(-\mbox{\boldmath$v$}\times\mbox{\boldmath$B$})\cdot{\mbox{\boldmath$y$}}$, respectively. The $v_{z}$ and density $\rho$ are fixed with initial value at the bottom boundary and $\phi$ is imposed on the Neumann condition throughout the simulation. Note that, during the relaxation process, $E_{z}$ at the bottom surface is zero, i.e., $\partial_{x}E_{z}$ = $\partial_{y}E_{y}$=0 because the evolutions of $v_{x}$ and $v_{y}$ are halted there ($v_{x}$ = $v_{y}$ = 0).

Finally, after the relaxation, we impose a diverging motion on the sunspot, which corresponds to a converging motion around the PIL, as shown in Fig.1d. The diverging motion is given as the following formula (e.g., Xia et al. 2014),

where $x_{d}$=$y_{d}$ =0.1 is given. The magnetic field at the bottom boundary follows an induction equation in each component as below,

note that the diffusion is given in the direction of converging motion (in this study, it corresponds to $x$ direction) to make the simulation more robust. This diverging motion is continuously imposed on the sunspot until the end of the simulation. In both cases, twisting motion and diverging motion, we normalized the velocity $\mbox{\boldmath$v$}^{(b)}$ by the maximum value $|v_{y:max}^{(b)}|$ and multiplied by 0.01, i.e., $\mbox{\boldmath$v$}^{(b)}\Rightarrow 0.01\mbox{\boldmath$v$}^{(b)}/|v_{y:max}^{(b)}|$, before being used in Equations (8) and (10), so that the maximum absolute value of the velocity is set as 0.01 at the bottom surface.

The same MHD equations (Eq.(2)-Eq.(6)) are applied to the data-driven MHD simulation while the electric fields, $\mbox{\boldmath$E$}^{(1)-(3)}$ are given at the bottom and a half-grid above and below as shown in Fig.2. Since the location that is one grid above the bottom ($k$=1 shown in Fig.2) corresponds to the practical bottom boundary, the derivation of the magnetic field at this position ($k$=1) is important for the data-driven simulation in this study. Especially, the derivative in z direction introduced in the induction equation should be handled with care because it includes the difference between non-MHD components ($\mbox{\boldmath$E$}^{(1)}$ and $\mbox{\boldmath$E$}^{(3)}$) which are placed at $k$=0 or 1/2 and the electric fields obtained from MHD process ($-\mbox{\boldmath$v$}\times\mbox{\boldmath$B$}$) which are placed at $k$=2.

According to a second order finite differential method, the derivative of E in the z-direction at $k$=1 in the induction equation is written as,

where $E$ is either the $x$ or $y$ component and we use the following approximation,

to obtain the last equation. When we take into account $E_{k=\frac{1}{2}}=E^{(1)}+E^{(3)}$ and $E_{k=0}=E^{(1)}$, we obtain

On the other hand, we can write down the different formula as follows,

where $E_{k=\frac{1}{2}}$ is approximated in Eq.(21). Eventually, we obtain the following equation,

Hereafter, the former written in Eq.(22) is denoted as Type-A and the later (Eq.(23)) is denoted as Type-B. We will mainly show the results obtained from Type-A and discuss the difference in results between Type-A and Type-B.

The coronal magnetic fields is driven as outlines in Fig. 3. We derive the electric field $\mbox{\boldmath$E$}^{(1)-(3)}$ at each time interval from the two output data at $t_{n}$ and $t_{n+1}$ ($n$=1,2,3$\cdots$) obtained from the referenced MHD simulation where only the magnetic field at the bottom is used. The time interval $t_{n+1}-t_{n}$ is set to 0.5625. The electric fields are derived from the three Poisson equations (14), (15) and (16). In this study, we use the Gauss-Seidel method to solve them. The performance and convergence of the Poisson equations are shown in Appendix (see Fig.A1). In the data-driven simulation, each electric field drives the upper coronal magnetic field at each interval. Note that, in this study, we drive the coronal magnetic field when the $|\mbox{\boldmath$B$}|=\sqrt{B_{x}^{2}+B_{y}^{2}+B_{z}^{2}}$ measured at the bottom is more than 0.0625, i.e., the weak magnetic field region does not change with time.

As seen in Table-1, we ran 8 different simulations which are mainly classified into the cases of TW and the cases of CV where TW and CV mean data-driven simulations in twisting (phase energy buildup phase: $t$=0$\sim$11.25 in the MHD simulation) and the converging phase (energy released phase: t$\geq$22.5 in the MHD simulation), respectively. The initial conditions of the data-driven MHD simulations in the cases of TW and CV employ the magnetic fields at $t$=0 and $t$=22.5 obtained from the ground-truth data. Note that the data-driven simulations are not conducted continuously from $t$=0 to after the eruption. Because one of objectives of this study is that how the data-driven simulation can produce the energy release phase. If we conduct the data-driven simulation from t=0 to the eruption, it becomes difficult to isolate whether the errors are due to the energy buildup phase (t=0-11.5) or the relaxation phase (t=11.5-22.5). We focus on evaluating the data-driven simulation in the energy buildup phase without the extra errors. The Type-A and -B differential methods denote -A (A0, A1, and A2 ), and -B, respectively. The Type-A is further classified into A0, A1, and A2 where the vertical velocity ($v_{z}$) placed at one grid ($k=1$) above the bottom surface has been given through a specific update in the temporal evolution and the method is different in each case. In the cases of A1 and A2, the vertical velocity ($v_{z}$) is fixed to zero or is update through a linear interpolation with values of $k=0$ and $k=2$, respectively, while the A0 has no specific updated, i.e., the vertical velocity is derived from an equation of motion directly. In case B, we run one simulation where the vertical velocity at $k=1$ is updated using the liner interpolation.

The temporal evolution of the magnetic and kinetic energies are shown in Fig.4. Both energies are buildup by $t\sim$ 11.25 due to the twisting motion of the sunspot, then gradually decrease by $t\sim$ 22.5 in the relaxation phase. Note that the twisting motion is halted at $t\sim$ 11.25 and the magnetic fields are relaxed by $t\sim 22.5$. Within a few moments of turning on the converging motion at $t\sim$ 22.5, the kinetic energy increases dramatically, whilst the magnetic energy decreases a lot of which is dissipated and converted into the Joule heating through the flux cancellation.

The evolution of the 3D magnetic field lines, while the twisting motion is imposed, is shown in Fig.5. As time passes, the magnetic field lines are sheared gradually and eventually a sigmoidal structure is formed at $t$=10.12 which is often observed prior to flares (e.g., Savcheva et al. 2012, Inoue et al. 2012, Kawabata et al. 2018).

After the relaxation, we impose the diverging flow on the sunspots as shown in Figs.6a-c with the evolution of the sunspots. From the upper panels, a part of the magnetic flux is transported toward the polarity inversion line (PIL) at which they are canceled. The lower panels (Figs.6d-e) show the evolution of the magnetic field lines. The footpoints of the magnetic field lines are carried by the flow on the both sunspots and field lines of opposite polarity encounter each other at the PIL. Magnetic reconnection then takes place, resulting in the formation of the long and highly twisted lines as seen in Figs.6e and f. The continuous converging flow around the PIL enhances the further twisted MFR and lifts it upward as seen in Fig.7b where the magnetic field lines just before the diverging flow are imposed are shown in Fig.7a. The current sheet is formed under the lifted MFR and eventually it erupts as shown in Figs.7c and d.

We have attempted to reproduce these MHD processes by using the data-driven MHD simulation in which the electric field works as the driving source on the bottom surface, instead of the velocity field. These electric fields are derived from a time series of the magnetic field on the bottom surface of the ground-truth data produced by the MHD simulation.

In above section, we found that the case of CV-A0 does not capture the ground-truth data well, for instance, in the temporal evolution of the kinetic energy and the axis of the MFR. One other notable difference when compared to the other cases from the data-driven simulation is that the eruption occurs later in the CV-A0 case. The difference between these cases is the handling of the vertical velocity ($v_{z}$) at $k=1$, i.e., $z=\Delta z$. Therefore, we should examine the behavior of the vertical velocity at $z=\Delta z$ or around the position, in both energy build up and energy release phases, respectively.

Figure 15 shows the temporal evolutions of the velocity field obtained from the ground-truth data and the data-driven MHD simulation (TW-A0 and TW-A2), respectively, at the surface at $z=3\Delta z$. Although the horizontal velocity plotted by the arrows forms a twisting motion in each case, the vertical velocity ($v_{z}$) distributions are quite different. The most striking difference is that the negative vertical velocity appears in late phase of the data-driven MHD simulations. So, we hereafter discuss the spacial and temporal evolutions of the vertical velocity.

Figure 16 shows the temporal evolutions of the vertical velocity along the center of the numerical box, i.e., the vertical dashed line as shown in Fig.14a, which are obtained from ground-truth data and the data-driven simulations, respectively, in the energy buildup phase (TW-A0 to TW-A2). For the ground-truth data, the velocity temporally increases according to the twisting motion given on the photosphere. Although the vertical velocity in TW-A1 and TW-A2 increases as times goes on, it enhances the negative value with time in the region close to the solar surface, which is different to the ground-truth data. The typical result is the case of TW-A0 in which the velocity is negatively increased close to the solar surface with time. The small inset is an enlarged view in the height range 0 to 0.01. We found that the velocity located at one grid above the bottom surface, i.e., $z=\Delta z$, is fastest way to enhance the negative value. The time variation of the vertical velocity in the case of TW-A0 is determined by the equation of motion while in TW-A1 it is set of zero and TW-A2 it is updated by a linear interpolation in the vertical direction, respectively. Either the convective term $\mbox{\boldmath$v$}\cdot\mbox{\boldmath$\nabla$}\mbox{\boldmath$v$}$ or Lorentz force $\mbox{\boldmath$J$}\times\mbox{\boldmath$B$}$, or both would have a negative influence on the vertical velocity. Despite each velocity profile begins different in each case of energy buildup phase, magnetic energy is found to follow almost the same behavior in each case (Fig.8) . Since the magnetic energy is dominant over the kinetic energy, the evolution of the magnetic field is almost unaffected by the velocity field, conversely a small fluctuation of the magnetic field would greatly influence the evolution of the velocity. Therefore, the Lorentz force would be a major factor inhibits the vertical velocity from leading the correct solution. Note that, in the real situation, since no one can say that the magnetic energy dominates and the velocity field does not affect the magnetic evolution much on the photosphere, the situation is expected to become more complex.

Figure 17 plots the temporal evolution of the vertical velocity in the energy release phase, i.e., the erupting phase of the MFR in the same format as in Fig.16. The data-driven simulations, CA-V1 and CA-V2 capture the ground-truth data well while the results obtained from CV-A0 show the different behavior. From the small inset follows the same format as in Fig.16b, the value at $z=\Delta z$ suddenly becomes the negative value. This behavior is the same as seen in the case of TW-A0 in the energy buildup phase. However, this case differs from TW-A0 in that the velocity recovers from its negative value during the late phase.

Figures 18a and b compare the temporal evolutions of the vertical velocity located at $z=\Delta z$ for the ground-truth data and the results obtained from between the data-driven simulation in energy buildup phase (TW-A0 and TW-A2) and energy release phase (CV-A0 and CV-A2). Note that the vertical velocity at $z=\Delta z$ for the cases of TW-A1 and CV-A1 is set to zero, so we exclude these plots. The TW-A2 and CV-A2 cases somewhat capture the ground-truth data, while cases TW-A0 and CV-A0 show large deviations from the ground-truth by exhibiting a steep drop in the negative value from the initial values. CV-A0 however is found to increase the velocity back towards that seen in the ground-truth data. It is likely that when the magnetic field is converted into the dynamics phase from the static phase, the strong positive enhancement of the velocity is associated with the eruption, which returns to the ground-truth data even at $z=\Delta z$.

As seen in these results, we need a careful treatment of the vertical velocity at $z=\Delta z$ in the static phase, which is found to largely deviate from the ground-truth data without the treatments. Since the kinetic energy in the energy buildup phase is much smaller than the magnetic energy, the evolution of the magnetic field is unaffected by the velocity field even if the velocity strongly deviates from the ground-truth data as seen in Fig. 8. On the other hand, it affects the initiation of the eruption, which causes delays because the down flow, which is an unlikely result in the ground-truth data, inhibits local reconnection at the photosphere and therefore the formation of the MFR. Thus, some treatments (simple linear interpolation was applicable in this study) would be required.

Why does the vertical velocity ($v_{z}$) located at $z=\Delta z$ show different behaviors in each case? From the above results, when the difference to the ground-truth data is striking, the magnetic field evolves in the quasi statical energy buildup phase toward the pre-erupting stage. Since the magnitude of the kinetic energy is much lower than the magnetic energy in the energy buildup phase, $\mbox{\boldmath$J$}\times\mbox{\boldmath$B$}$ is the major factor to lowering velocity. Each component of the Lorentz force ($F$) is described as

where $B_{z}$ at $z=\Delta z$ is derived from the rotation of the electric fields located at $z=\Delta z$ while a derivative of $E_{y}$ and $E_{x}$ in z component is included to derive the $B_{x}$ and $B_{y}$, respectively, according to equation (22). $B_{z}$ is determined by $E_{k=1}$ which is derived from MHD electric field ($-\mbox{\boldmath$v$}\times\mbox{\boldmath$B$}$) while $B_{x}$ and $B_{y}$ are derived from $E_{k=2}$, $E_{k=1}$ and $E^{(3)}$ too (see Equation (22)). $E_{k=2}$ and $E_{k=1}$ are different from $E^{(3)}$ because the $E^{(3)}$ is derived from a non-MHD process. Therefore, the z-derivative is generally not allowed at $z=\Delta z$ and the physical consistency regarding of the obtained $B_{x}$ and $B_{y}$ is not guaranteed. $F_{z}$ includes more $B_{x}$ and $B_{y}$ than the other components $F_{x}$ and $F_{y}$, so the value of $v_{z}$ that is derived from $F_{z}$ will have a lesser accuracy.

We now discuss results obtained from the data-driven simulation for the case of TW-B that implements a differential method type-B described in equation (23). We first compare with the results obtained from the ground-truth data. Figures.19a and b plot the temporal evolutions of the magnetic and kinetic energies in the energy buildup phase, i.e., the magnetic field is twisted due to the photospheric motion. Both evolutions obtained from the data-driven simulation shows large deviations from the ground-truth data. Figures 19c, d, and e show the distribution of $|\mbox{\boldmath$J$}|/|\mbox{\boldmath$B$}|$ plotted on the vertical cross section for the ground-truth data and the data-driven MHD simulations (TW-A2 and TW-B), respectively, from which the TW-A2 case reproduces the ground-truth data well, while the case of TW-B shows a totally different solution.

The difference between the TW-A2 and TW-B is the application of differing numerical differential methods, especially, at $z=\Delta z$ at $k=1$. The derivative forms in the z direction for each are as follows, with type-A written as

and type-B is written as

Type-A is an average value of the central differential method of electric fields obtained from the MHD process (the first term in the right hand side of Eq.(25)) and non-MHD process (the second term). Although each derivative has physical meaning, there is a discrepancy between the MHD term and non-MHD term. The use of the average process in this case smoothly connects them, making a more robust calculation. Type-B is a central differential method of the electric fields obtained from the MHD and non-MHD processes. In principle, this derivative is not allowed as discussed in above section because the origin of those electric fields is different, i.e., there are not connected smoothly. Therefore, this discrepancy might negatively impact the simulation result.

On the other hand, Hayashi et al. (2018) and Hayashi et al. (2019) implement type-B and the simulation works well. One reason for the difference in results from this study is the difference in the numerical schemes used. Hayashi et al. (2018) and Hayashi et al. (2019) are designed with the finite volume method, and numerical viscosity and resistivity might efficiently work on the negative impact and make the simulation stable. However, such viscosity and resistivity do not work efficiently in the central differential method as used in this study despite these being included explicitly in the induction equation and the equation of motion. Therefore, numerical instability is not suppressed in this study and the solution is not reproduced correctly as seen in Fig.19e.

We test the effect on the accuracy of our results to the time cadence used in the data-driven MHD simulations used in this study. Figure 20a shows the comparison of the temporal evolution of the kinetic energies in the energy buildup phase that is obtained form the ground-truth data, the data-driven simulation (TW-A2), and a new data-driven simulation with 2.5 times the temporal resolution ($t_{n+1}-t_{n}$=0.225) of TW-A2 ( TW-A2D). These are plotted in black, red, and blue, respectively. Since the blue line almost overlaps with the red line, no significant differences due to the temporal evolution are apparent in the energy buildup phase. Thus we can say that the solution is convergent with respect to temporal resolution.

Figure 20b shows the each kinetic energy during the energy release phase, plotted in the same format as Fig. 20a. The red line, which corresponds to the data-driven simulation (CV-A2), diverges slightly from the blue line (CV-A2D) that has 2.5 times the temporal resolution of CV-A2. This is difference is from the result in the energy buildup phase. However, the difference between CV-A2 and CV-A2D is only a few percent. From these results, our solution is almost convergent in temporal resolution. According to Leake et al. (2017), using a higher temporal cadence makes the results of a data-driven simulation closer to the ground-truth data. However, we do not see this in our results. One of possible reason could be due to differences in the methods used to create the simulation. Leake et al. (2017) directly uses data that are extracted from the ground-truth data where-as our method drives the coronal magnetic field by electric fields given near the solar surface that are derived from the bottom magnetic fields of the ground-truth data. It is therefore likely that these differing methods are the reason for the differing effect to adjusting the temporal resolution of the simulation.

Furthermore, it is important to consider a ratio of the sampling time ($dt_{s}$) to the dynamical time ($dt_{d}$) proposed by Leake et al. (2017) and the ratio should be less than 1. The dynamical time is defined by $dt_{t}$=$\Delta/V_{p}$ where $\Delta$ is a length of the grid interval and $V_{p}$ is velocity of the photosphere. In this study, since $\Delta$=$1/320$ and maximum $V_{p(max)}$ = 0.01, $dt_{t(min)}$ is given as 0.3125. The sampling time $dt_{s}$ is given as 0.5625 (low time resolution cases) and 0.225 (high time resolution cases,), therefore the ratio corresponds to 1.8 and 0.72, respectively. The ratio in the low time resolution cases (other than TW-A2D and CV-A2D) is over 1. However, in this work, the bottom-boundary parameters driving the system are from a model and thoroughly given in form of functions of time and position. From this reason, there is little uncertainty in the boundary data (sampling data) and there is not much difference in energies between the low time resolution cases and high time resolution cases.

We run the data-driven simulation throughout the entire process from the energy accumulation process to the release process. The result is shown in Fig. 21(a) in which the temporal evolution of the ground-truth data and the data-driven simulation is plotted in black and red, respectively. We found that the data-driven simulation cannot cause the eruption that is different from the results discussed above. One of possible reasons is the pre-erupting magnetic field. Figures 21(b) and (c) show the pre-erupting magnetic fields of the ground-truth data and produced by the data-driven simulation, respectively. The size of the sigmoidal structure is different between them. Fig. 10(b) shows the magnetic flux of the pre-erupting magnetic fields that depends on the twist value before the relaxation process. From this result, the maximum twist value of the ground-truth data is in a rage from 0.6 to 0.7. On the other hand, the highly twisted field lines buildup in the data-driven simulation are concentrated approximately up to $T_{w}$=0.6. This may be very small to cause the eruption. For instance, Inoue et al. (2018) suggested that to cause the solar flares associated with the eruptions requires the twisted field lines with twist values from 0.5 to 1.0. Therefore, the ground-truth data does not have much of twist and it would be difficult for the magnetic field produced by the data-driven simulation to cause the eruption. This result indicates that a little difference in twist may affect the eruption. Furthermore, this is a reminder of that the 3D solution is not uniquely determined even using same boundary condition.

Although the problem may be caused by the ground-truth data that has weakly twisted field lines before the eruption in this study, we should avoid the problem of weakening of the twist of the magnetic field lines produced in the data-driven simulation when using the photospheric magnetic field. The NLFFF would be useful as the initial condition of the data-driven simulation before the eruption. Because, in some cases, the NLFFF before the eruption reproduces more highly twisted field lines than those shown in Fig. 21(c) (Inoue 2016). However more discussion is needed on the issue of whether to use the magnetic field produced in the data-driven simulation or the NLFFF as the initial condition. This would be addressed in future work.