Extrapolation of Three-dimensional Magnetic Field Structure in Flare-Productive Active Regions with Different Initial Conditions

Solar active regions contain the strong magnetic field and sometimes produce explosive events, such as solar flares and coronal mass ejections (CMEs) by releasing the magnetic energy. The magnetic reconnection releases the non-potential magnetic energy, heats and accelerates the plasmas. What is the condition on the onset of the solar flares and CMEs is one of the important open questions for solar flares and CMEs. Inferring correct 3D magnetic field is a crucial task to answer the question.

The magnetic field in the solar atmosphere is usually measured through the spectropolarimetric observations with slit-based or filter-based instruments. The state-of-art spacecrafts, such as Hinode (Kosugi et al., 2007) and Solar Dynamics Observatory (SDO; Pesnell et al., 2012) have enabled us to observe magnetic field in the photosphere with high spatial resolution and high polarimetric accuracy. While photospheric magnetic field can be inferred with high accuracy through spectropolarimetric observations, it is difficult to obtain magnetic field in the corona. The coronal lines have low brightness, are optically thin, and suffer from larger Doppler broadening.

The force-free field modeling is one of the alternative methods to infer the 3D magnetic field in the solar corona. The main concept of the force-free field modeling is to extrapolate the coronal magnetic field from the spatial map of the magnetic field in the photosphere based on two assumptions (Wiegelmann & Sakurai, 2012). First assumption is the mechanical equilibrium of the plasmas in the solar corona. Second assumption is the domination of the Lorentz force in the solar corona. In the solar corona, the plasma $\beta(=8\pi p/B^{2})$, which is the ratio between the plasma pressure and the magnetic pressure, is thought to be sufficiently small, $\beta\ll 1$ (Gary, 2001). The above two assumptions lead to the condition that the Lorentz force vanishes in the solar corona,

where $j$ is the current density. The current density follows the Ampére’s law,

Equation (1) can be written

where $\alpha$ is called the force-free parameter. Generally, $\alpha$ is not a constant in space and such magnetic field distribution is called nonlinear force-free field (NLFFF).

There exist mainly two approaches to solve NLFFF equations. First approach is to solve the NLFFF equations with the boundary conditions, which is a mathematically well-posed problem. The example of this approach was proposed by Grad & Rubin (1958), and so is often called Grad-Rubin method. In this approach, the distribution of $\alpha$ in only one polarity and the vertical magnetic field are prescribed to the bottom boundary. One weak point of this approach is that the vector magnetic field on the bottom boundary is not consistent with the observed photospheric magnetic field. Second approach is to find the closest force-free equilibrium field matching the observed vector magnetic field in the photosphere, which is prescribed to the bottom boundary. While this approach keeps the bottom boundary consistent with the observed magnetic vector in the photosphere, this approach is an ill-posed problem and there is no proof of unique and stable solutions. The examples of this approach are the optimization methods (Wheatland et al., 2000; Wiegelmann & Neukirch, 2006) and the magnetohydrodynamics (MHD) relaxation methods (Chodura & Schlueter, 1981; Mikic & McClymont, 1994; Inoue et al., 2014).

The NLFFF applications to the solar observations have some uncertainties, e.g., observational errors and numerical methods. The question is, in practical way, how large the uncertainty of the NLFFF modeling is. In other words, is there a possibility that completely different NLFFF results are obtained based on the same bottom boundary? Schrijver et al. (2006) applied six NLFFF algorithm to analytical force-free field solutions and compared solutions with each other. The best method reproduced the total energy in the magnetic field within 2% error. Metcalf et al. (2008) also compared NLFFF methods by using solar-like reference model. They applied NLFFF methods to both forced (not force-free) and force-free bottom boundaries. They show that while the NLFFF with the force-free bottom boundary reproduced helical flux bundle, that with forced bottom boundary did not reproduce reference model well. De Rosa et al. (2009) applied various kinds of NLFFF algorithm (Grad-Rubin, optimization, and MHD relaxation) to the photospheric magnetic field observed with Hinode and investigated how different the solutions from different algorithm are. They found that field lines from NLFFF do not agree with EUV images and the average misalignment is 20-40 degrees. Thalmann et al. (2013) investigated the effect of the observational instrument on the solution of NLFFF. They used a photospheric magnetic field observed with different instrument, i.e. , the Hinode/SOT SP and SDO/Helioseismic and Magnetic Imager (HMI; Schou et al., 2012). They concluded that the relative estimates such as normalized magnetic energy and the overall structure of the magnetic fields might be reliable, although there were remarkable deviations in the absolute value between the instruments. DeRosa et al. (2015) investigated the influence of spatial resolution on NLFFF with three kinds of methods similar to De Rosa et al. (2009). They showed that the free energy tends to be higher with increasing resolutions and magnetic helicity values vary significantly among different resolutions.

In previous studies above, they showed how different results NLFFF produced. The questions is whether the completely different solutions are produced or not based on the same bottom boundary. In other words, how unique is the solution when focusing on one method? If completely different solutions exist, how can we obtain the result most consistent to the coronal imaging observations? Although many previous studies applied NLFFF to solar active regions and showed qualitatively agreement between NLFFF and coronal images (Régnier & Amari, 2004; Valori et al., 2012; Sun et al., 2012; Inoue et al., 2012; Jiang et al., 2014; Wang et al., 2014; Kawabata et al., 2017; Muhamad et al., 2018), there still exist some regions where NLFFF produce quite different field lines from coronal images (De Rosa et al., 2009). Therefore, the problem we have to tackle is to reveal whether completely different solution with same bottom boundary can exist or not in the practical calculation and what we should do in case that the modeled field lines do not agree with the coronal images.

To investigate whether completely different solutions with the same bottom boundary exist or not, we evaluate NLFFF extrapolation with different initial conditions. As is often the case with the nonlinear inverse problems, the different initial guesses may often produce completely different converged solutions. Therefore the uniqueness of the NLFFF calculation can be studied by giving different initial conditions.

The NLFFF extrapolation is usually performed as follows,

(1)

Set the 3D initial condition by using photospheric vertical magnetic field.

(2)

Give the information of the horizontal magnetic field to the bottom boundary.

(3)

Perform some relaxation process.

In the process (1), almost all previous studies use the potential field as an initial condition. We perform NLFFF calculation not only with the potential field but also with linear force-free field as an initial condition. Several force-free alpha values are given in the linear force-free case. Comparisons among the NLFFF results with different initial conditions will provide some insights to the uniqueness of the NLFFF extrapolation. As described above previous studies show that different results are obtained depending on the calculation methods used in the NLFFF modeling (Schrijver et al., 2006; Metcalf et al., 2008; Schrijver et al., 2008; De Rosa et al., 2009). The one of the causes of the difference between the Grad-Rubin method and the other two methods (optimization and MHD relaxation methods) is the treatment of the bottom boundary. Because the former method relax the 3D magnetic field so that the force-free alpha along the magnetic field line is constant, the horizontal component of magnetic field at the bottom boundary is different from the observational magnetic field in the photosphere. On the other hand, the latter methods keep the bottom boundary as the observational magnetic field vector in the photosphere. We give priority to keeping the magnetic field measured in the photosphere at the bottom boundary. Therefore, we focus on one method, MHD relaxation method, which is developed by Inoue et al. (2014).

The paper is organized as follows: We describe the observations and data reduction in Section 2. The method of the NLFFF extrapolation is described in Section 3. We present and discuss the results in Sections 4 and 5, respectively. We summarize our results in Section 6.

As targets to analyze, we chose two active regions, NOAA 11692 and NOAA 11967 to investigate how significantly different coronal magnetic field structures are derived depending on initial conditions. The examined active regions are two extreme examples; one has rather simple bi-pole magnetic distribution at the photosphere and the other shows a complicated field distribution at the photosphere. The photospheric vertical magnetic field distributions are shown in the left upper and lower panels of Figure 1. While the former is composed of bi-pole magnetic fields and shows a weak twist in the photosphere, the latter has multi-pole magnetic fields and shows a strong twist. The global alpha is one of the index values for expressing the degree of the the twist in the active regions. Tiwari et al. (2009a, b) defined global alpha $\alpha_{g}$ such as,

which can be calculated from the photospheric vector magnetic field, $B_{x}$, $B_{y}$, and $B_{z}$, where z-axis is defined as the vertical direction to the solar surface. A global alpha of NOAA 11692 is $\alpha_{g}=-1.0\times 10^{-8}{\rm m^{-1}}$ and NOAA 11967 has $\alpha_{g}=-5.0\times 10^{-8}{\rm m^{-1}}$. Absolute values of $\alpha_{g}$ in various active regions usually vary from $0$ to $5.0\times 10^{-8}{\rm m^{-1}}$ (Pevtsov et al., 1995, 1997). Therefore, the two active regions we analyzed are contrasting examples in terms of the twist.

The nonlinear force-free field extrapolation is performed by the MHD relaxation method (Inoue et al., 2014; Inoue, 2016), which uses the following equations,

where $\rho$ is the pseudo density, which is assumed to be equal to $|\mbox{\boldmath$B$}|$ to ease the relaxation by equalizing the Alfvén speed in space, $\phi$ is the convenient potential for $\nabla\cdot\mbox{\boldmath$B$}$ cleaning and $\nu$ is the viscosity, which is set to a constant ($1.0\times 10^{-3}$). The length, magnetic field, velocity, and time were normalized by $L_{0}=157$ Mm for NOAA 11967 and $L_{0}=314$ Mm for NOAA 11692 and $B_{0}=4000$ G, $V_{\rm A}\equiv B_{0}/(4\pi\rho_{0})^{1/2}$, and $\tau_{A}\equiv L_{0}/V_{A}$, where $V_{A}$ is the Alfv$\rm{\acute{e}}$n velocity. Equations (5), (6), (7) and (8) are the equation of motion, the induction equation, the ${\rm Amp\grave{e}re^{\prime}s}$ law, and $\nabla\cdot\mbox{\boldmath$B$}$ cleaning introduced by Dedner et al. (2002), respectively. The parameters $c_{p}^{2}$ and $c_{h}^{2}$ are the advection and diffusion coefficients, respectively and fixed at $0.1$ and $0.04$. The non-dimensional resistivity $\eta$ is given by

where $\eta_{0}$ and $\eta_{1}$ are fixed at $5.0\times 10^{-5}$ and $1.0\times 10^{-3}$ in non-dimensional units. The second term is introduced to accelerate the relaxation to the force-free state.

The velocity field at each grid was adjusted at each time step below in order to avoid becoming large value. When the value of $v^{*}$ becomes larger than the value of $v_{\rm{max}}$,

where $v^{*}=|\mbox{\boldmath$v$}|/|\mbox{\boldmath$v_{A}$}|$, and $v_{{\rm max}}=0.1$.

In previous NLFFF calculations, the potential field has been used as initial guess for the 3D magnetic field structure. We chose linear force free field as an initial condition including a potential field, which satisfies $\nabla\times\mbox{\boldmath$B$}=\alpha_{0}\mbox{\boldmath$B$}$, where $\alpha_{0}$ is a constant. We examined different 5 initial conditions $\alpha_{0}=[0,\pm 1.2,\pm 2.3]\ \times 10^{-8}{\rm m^{-1}}$ for NOAA 11692 and 12 initial conditions $\alpha_{0}=[0,\pm 0.70,\pm 1.2,\pm 2.3,\pm 4.6,\pm 7.0,-12]\ \times 10^{-8}{\rm m^{-1}}$ for NOAA11967. The values of $\alpha_{0}$ are chosen by the pixel size. The $\alpha_{0}$ values of $[0,\pm 0.70,\pm 1.2,\pm 2.3,\pm 4.6,\pm 7.0,-12]\times 10^{-8}$ $m^{-1}$ correspond to $[0,\pm 0.003,\pm 0.005,\pm 0.01,\pm 0.02,\pm 0.03,-0.05]$ pix^-1 for NOAA11967. Because the magnetic field of NOAA 11967 is less affected by the initial condition than that of NOAA 11692, we have more initial conditions for NOAA 11967.

We used the equations of Alissandrakis (1981) for the calculations of the linear force-free field and the resulting initial conditions are shown in Figures 2 and 3. Although the case $\alpha_{0}=+12\ \times 10^{-8}{\rm m^{-1}}$ for NOAA 11967 was also calculated, the calculation did not converge. Therefore, we do not include it in the results in this paper. The numerical domain is set to $(0,0,0)<(x,y,z)<(1.0,0.7,0.7)$ resolved by $360\times 252\times 252$ nodes for NOAA11692 and $(0,0,0)<(x,y,z)<(1.5,1.0,0.5)$ resolved by $540\times 360\times 180$ nodes for NOAA 11967. In order to set the same top boundary for all calculations, the initial conditions above $z=0.417$ are set to the potential field.

The magnetic field at the top was fixed with the initial state (potential field) and the normal component of the magnetic field on the bottom boundary was also fixed. The side boundary was periodic. We varied the transverse component on the bottom boundary $\mbox{\boldmath$B$}_{\rm BC}$ as follows,

where $\mbox{\boldmath$B$}_{\rm obs}$ and $\mbox{\boldmath$B$}_{\rm initial}$ are the transverse component of the observational and initial bottom boundary, respectively. We increased $\gamma=\gamma+d\gamma$ when $\int|\mbox{\boldmath$j$}\times\mbox{\boldmath$B$}|^{2}dV$ dropped below a critical value. In this study we set $d\gamma=0.1$. When $\gamma$ becomes equal to $1$, $\mbox{\boldmath$B$}_{BC}$ is consistent with the observed field. The number of calculation steps were set to 24000 steps for NOAA 11692 and 25000 steps for NOAA 11967

Spatial derivatives are calculated by the second-order central differences and temporal derivatives are integrated by the Runge-Kutta-Gill method to fourth order accuracy.

We calculated NLFFF with 5 and 12 different initial conditions for NOAA 11692 and 11967, respectively. Summary of our results is as follows.

• •

  According to the comparison with the soft X-ray image, the NLFFF shows better correspondence in the simple active region NOAA 11692, when the initial constant force-free alpha is the same sign with the global alpha calculated from the photospheric magnetic field, as shown in Figure 5. On the other hand, in the complex multi-pole active region NOAA 11967, the results of the NLFFF extrapolations are less dependent on the initial condition, as shown in Figure 7.

• •

  Total magnetic energy of the NLFFF extrapolation does not strongly depend on the initial condition as shown in Tables 2 and 3. The dependence of the free energy is larger compared to the total magnetic energy.

• •

  The solution of NLFFF at the region where the strong magnetic field exists, e.g., magnetic field in the lower height, tends to be less affected by the initial condition, as shown in Figure 9. On the other hand, the region in the weak magnetic field around the polarity inversion line, tends to be affected by the initial condition.

• •

  Except for the case $\alpha_{0}=-12\times 10^{-8}{\rm m^{-1}}$ in NOAA 11967, the NLFFF extrapolation is considered to be converged, as shown in Figure 18. The increase of the calculation steps may not affect the results in this study.

From Figure 9, the magnetic field in the lower height region tends to be less affected by the initial condition, while the magnetic field at the higher region is strongly affected by the initial condition. There is a possibility that this result is caused by the convergence problem of the NLFFF modeling, however this issue is ruled out by an analysis of the Lorentz forces as a function of height, as shown in the Appendix.

Our results suggest that NLFFF models calculated using relaxation methods, including the method used in this study, permit different results (depending on the initialization field used) in the upper regions of the computational domain. In other words, there exist several local minima for the force-free equilibrium in the higher regions. Because the initial condition is force-free everywhere except at the bottom boundary, a larger degree of freedom is allowed at the higher region compared to the lower region. For example, in our MHD relaxation method, the information of the bottom boundary is transported by the propagation of the psuedo-Alfvén wave as mentioned above. The region which is far from the boundary tends to be less affected by the psuedo-Alfvén wave and remain almost unchanged from the initial equilibrium state (potential-like or LFFF-like). Therefore, the NLFFF extrapolation can not reproduce the magnetic field of active regions well such as have twisted structures in the corona, but have potential-like photospheric magnetic field. Note that this possibility does not prove that the degree of freedom of the NLFFF solution is mathematically larger at the higher region because certain amount of the Lorentz force remains in our calculation. To solve this problem, we suggest that additional observational limitation should be given to the current NLFFF modeling such that the magnetic field at the higher region can converge to the NLFFF result, which is consistent with X-ray and/or EUV imaging observations.

We compared NLFFF results with coronal loops observed with Hinode/XRT. For NOAA 11967, 3D magnetic field configuration matches well, and appears more consistent with X-ray observation when the absolute value of initial force-free alpha is smaller. On the other hand, NOAA 11692 shows strong initial condition dependence in Figure 5. The clear difference between NOAA 11967 and 11692 is the complexity of the photospheric magnetic field. Our results suggest that the NLFFF result of NOAA 11967 is less affected by the initial condition than that of NOAA 11692. Because we analyze only two active region, we can not determine the cause of this result. Therefore, we discuss the candidates of the cause of difference in terms of the property of two active regions. There are five possibilities to explain this cause.

1. 1.

   NOAA 11967 has more magnetic flux, $8.4\times 10^{22}$ Mx, than that of NOAA 11692, $3.7\times 10^{22}$ Mx. This suggests that there are more strong magnetic field region above the photosphere.

2. 2.

   Unbalanced magnetic flux at the photospheric height may also be the cause. In the NLFFF modeling, magnetic flux unbalance at the bottom boundary may produce inconsistent results and may be related to the initial condition dependence. The net vertical magnetic flux normalized by the total magnetic flux is -0.03 for NOAA 11692 and 0.17 for NOAA 11967. This shows that the magnitude of unbalance of vertical magnetic flux is larger in NOAA 11967 than in NOAA 11692. Because the NLFFF results of NOAA 11967 are less affected by the initial condition, the flux unbalance in this study may not affect the initial condition dependence.

3. 3.

   The force-freeness at the photospheric height is different between the active regions. Although force-freeness of the both active regions is relatively small, that of NOAA 11967 ($|F_{z}|/F_{p}=0.089$) is smaller than that of NOAA 11692 ($|F_{z}|/F_{p}=0.17$). The large value of the force-freeness may affect the dependence of the initial condition on the NLFFF result.

4. 4.

   The length scale of the field lines may affect the results. While the height of the loop top in NOAA 11692 is around 100 Mm as shown in Figure 6, that in NOAA 11967 is around 30 Mm as shown in Figure 8. This result means that the height of the coronal loops identified in the soft X-ray images are different between NOAA 11692 and 11967. Our result shows that the magnetic field in the higher region tends to be affected more strongly by the initial condition than that in the lower region. Therefore, in comparison with the soft X-ray image in NOAA 11692, we compared the field lines, which tend to be affected by the initial condition. However, only the difference of the length scale of the coronal loops can not explain the difference of the initial condition dependency. As shown in Figures 12 and Figures 15, NOAA 11692 is more dependent on the initial condition than NOAA 11967 even in the same height. This result indicates that only the difference of the length scale is not the cause of the difference of the dependence between the two active regions.

5. 5.

   There is a difference of the horizontal field distribution between the two active regions. As shown in Figure 1, strong horizontal magnetic field occupies a large part of the FOV of NOAA 11967 for all four polarities. On the other hand, in NOAA 11692, while the horizontal magnetic field can be seen in the negative sunspot, there is no strong horizontal magnetic field in the positive polarity. As shown in Figure 4, while 6.3 $\%$ of the horizontal magnetic field in the FOV is larger than 1000 G for NOAA 11967, 0.49 $\%$ is larger than 1000 G for NOAA 11692. This result indicates that in relaxation process, it is difficult to give sufficient electric currents (non-potential magnetic field) to the upper atmosphere in the case of NOAA 11692 because there are less currents and forces in the photosphere. As a result, there are strong dependency on initial condition in the NLFFF modeling in NOAA 11692. In NOAA 11967 case, the active region has strong horizontal field, which give strong perturbations to the upper atmosphere in the relaxation process. This may be why the NLFFF modeling of NOAA 11967 converged to a consistent solution with different initial conditions.

As shown in Tables 2 and 3, the total magnetic energy does not strongly depend on the initial condition. On the other hand, the free energy shows larger difference among each solution than the total magnetic energy. Since the free energy is defined as the deviation between the magnetic energy and potential magnetic energy, the small difference of the total magnetic energy becomes large difference in free energy. This ratio is smaller than previous studies focusing on other dependences, such as method dependence (De Rosa et al., 2009, total energy: 1.9 between ${\rm Reg^{+}}$ and ${\rm Am1^{-1}}$), instrument dependence between Hinode and SDO (Thalmann et al., 2013, total energy: $\sim 2.2$, free energy: $\sim 2.5$), and spatial resolution dependence (DeRosa et al., 2015, total energy:$\sim 1.4$, free energy: $\sim 2$ in magnetofrictional method). This means that the initial condition dependence of the total magnetic energy and free energy in our MHD relaxation method is comparatively small. The uniqueness of total energy and free energy can be explained by using our results. Since the magnetic energy and free energy concentrate in the lower height, they become less dependent to the initial condition.

Our study shows that the NLFFF results determined from relaxation methods may be strongly affected by the initial condition in the higher levels of the computational domain. We therefore recommend checking whether field lines from NLFFF models are consistent with coronal images. Our results also indicate that the sigmoidal structure can be well reproduced by changing initial condition even when the photospheric magnetic field has less twist. Regarding NOAA 11692, our NLFFF modeling produce more consistent results with X-ray observations, when we use the initial values of $\alpha_{0}=-2.3\times 10^{-8}{\rm m^{-1}}$. This value has the same sign of the global alpha in the photosphere ( $\alpha_{g}=-1.0\times 10^{-8}{\rm m^{-1}}$). Although the global alpha does not strictly correspond to the global magnetic twist when the magnetic field is not force-free, the sign of the global alpha gives the sign of global shear signed angle (Tiwari et al., 2009b). Therefore, our results offer an important suggestion that the photospheric global alpha can be used for the rough estimation of the initial condition for better NLFFF modeling.

Summarizing our findings,

(1)

The solution of NLFFF at the region where the strong magnetic field exists, e.g., magnetic field in the lower height ($<$10 Mm), tends to be less affected by the initial condition, although the Lorentz force is concentrated at the lower height.

(2)

Total magnetic energy of the NLFFF extrapolation does not strongly depend on the initial condition.

(3)

The NLFFF extrapolation of the complex active region NOAA 11967 is less dependent on the initial condition compared to that of NOAA 11692.

We proposed the problem whether completely different solutions with the same bottom boundary exist or not. We conclude that we obtain completely different 3D NLFFF structure from the different initial conditions with the same bottom boundary. However, the initial condition dependence is small (the correlation coefficient $C>0.9$) where the the magnetic field is strong, e.g., in the lower height ($<$10 Mm). We also reveal that the 10-100 times larger Lorentz force, which is normalized by the square of the magnetic field strength, remains at the lower height ($<$ 10 Mm) than that at higher region ($>$ 10 Mm). The magnitude of the dependence is also different between the two active regions.

Acknowledgements We appreciate the feedback offered by the anonymous referee. We would like to thank Tatsu-hiko Miura, Yota Watanabe, Yuichi Ike, Yohei Ema, Tomoya Kenmochi, Saiei Matsubara, Tsuyoshi Yoneda, Yusaku Tiba, and Takahito Kashiwabara for informative discussions. Hinode is a Japanese mission developed and launched by ISAS/JAXA, collaboratoing with NAOJ as a domestic partner, NASA and STFC (UK) as international partners. It is operated by these agencies in cooperation with ESA and NSC (Norway). We gratefully acknowledge the SDO/HMI team for providing data. Our calculations of NLFFF modeling were performed on JAXA Supercomputer System generation 2 (JSS2). This work was supported by MEXT/JSPS KAKENHI Grant Numbers JP15H05814.