NUMERICAL STUDY ON THE EMERGENCE OF KINKED FLUX TUBE FOR UNDERSTANDING OF POSSIBLE ORIGIN OF δ-SPOT REGIONS

We solved the MHD equations in the following form:

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho }{\partial t}=-{\rm{\nabla }}\cdot (\rho {\boldsymbol{v}}),\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho {\boldsymbol{v}}}{\partial t}=-{\rm{\nabla }}\cdot \left[\rho {\boldsymbol{vv}}+\left(p+\displaystyle \frac{{| {\boldsymbol{B}}| }^{2}}{8\pi }\right)\underline{{\bf{1}}}-\displaystyle \frac{{\boldsymbol{BB}}}{4\pi }\right]+\rho {\boldsymbol{g}},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{B}}}{\partial t}=-{\rm{\nabla }}\cdot ({\boldsymbol{vB}}-{\boldsymbol{Bv}})-{\rm{\nabla }}\times (\eta {\rm{\nabla }}\times {\boldsymbol{B}}),\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}\displaystyle \frac{\partial e}{\partial t} & = & -{\rm{\nabla }}\cdot \left[{\boldsymbol{v}}\left(e+p+\displaystyle \frac{{| B| }^{{\bf{2}}}}{8\pi }\right)-\displaystyle \frac{1}{4\pi }{\boldsymbol{B}}({\boldsymbol{v}}\cdot {\boldsymbol{B}})\right]\\ & & +\displaystyle \frac{1}{4\pi }{\rm{\nabla }}\cdot ({\boldsymbol{B}}\times \eta {\rm{\nabla }}\times {\boldsymbol{B}})+\rho ({\boldsymbol{g}}\cdot {\boldsymbol{v}}),\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&e=\displaystyle \frac{p}{\gamma -1}+\displaystyle \frac{\rho {| {\boldsymbol{v}}| }^{2}}{2}+\displaystyle \frac{{| {\boldsymbol{B}}| }^{2}}{8\pi },\end{eqnarray} \tag{ 5 }$$

where ρ is the density, ${\boldsymbol{v}}$ is the velocity vector, ${\boldsymbol{B}}$ is the magnetic field vector, e is the total energy density, and p is the gas pressure. γ is the specific heat ratio (here $5/3$). The equation of state for the ideal gas is used. The gravity ${\boldsymbol{g}}={-g}_{0}\hat{{\boldsymbol{z}}}$ is a constant vector and the non-dimensional form is given by ${\boldsymbol{g}}=(0,0,-1/\gamma )$. The normalization units of our simulations are summarized in Table 1. In the expanding magnetic field in the corona, the plasma-β can become very small, where the calculated gas pressure can be negative due to numerical errors. To avoid the numerical instability, the energy equation in the low-beta ($\beta \lt 0.01$) regions is replaced by the following equation:

$$\begin{eqnarray}&&\displaystyle \frac{\partial {e}_{\mathrm{int}}}{\partial t}=-{\rm{\nabla }}\cdot ({\boldsymbol{v}}{e}_{\mathrm{int}})-(\gamma -1){e}_{\mathrm{int}}{\rm{\nabla }}\cdot {\boldsymbol{v}},\end{eqnarray} \tag{ 6 }$$

where ${e}_{\mathrm{int}}=p/(\gamma -1)$. In Equation (3), η is the magnetic diffusivity, and an anomalous resistivity model is assumed:

$$\begin{eqnarray}\eta =\left\{\begin{array}{ll}0 & (J\lt {J}_{c}\;\mathrm{or}\;\rho \gt {\rho }_{c})\\ {\eta }_{0}(J/{J}_{c}-1) & (J\geqslant {J}_{c}\;\mathrm{and}\;\rho \lt {\rho }_{c}),\end{array}\right.\end{eqnarray} \tag{ 7 }$$

where ${\eta }_{0}=0.1$, J_c = 0.1 and ${\rho }_{c}=0.1$. Density is normalized such that $\rho =1$ corresponds to the density at the photospheric base (see also Table 1). Thus, the resistivity works only in the region where the current density is strong above the photosphere.

Table 1.  Normalization Units

Quantity	Unit	Value
Length	Hp0	170 km
Velocity	${C}_{s0}={\left[\gamma ({k}_{B}/m){T}_{0}\right]}^{1/2}$	6.8 km s−1
Time	$\tau ={H}_{p0}/{C}_{s0}$	25 s
Temperature	${T}_{0}={T}_{\mathrm{ph}}$	5,600 K
Density	${\rho }_{0}={\rho }_{\mathrm{ph}}$	$1.4\times {10}^{-7}$ g cm${}^{-3}$
Pressure	$\gamma ({k}_{B}/m){\rho }_{0}{T}_{0}$	$6.3\times {10}^{4}$ dyn cm−2
Magnetic field strength	B0	250 G

Download table as:  ASCIITypeset image

The numerical scheme is based on Vögler et al. (2005): Spatial derivatives are calculated by the 4th-order central differences and temporal derivatives are integrated by a four-step Runge–Kutta scheme. We also adopt the artificial diffusion term described in Rempel et al. (2009) to stabilize the numerical calculation. Errors caused in ${\rm{\nabla }}\cdot {\boldsymbol{B}}$ are controlled by using an iterative hyperbolic divergence cleaning method, where the cleaning technique is based on the method described in Dedner et al. (2002).

The initial background atmosphere consists of three regions: an adiabatically stratified static atmosphere (representing the convection zone), a cool isothermal layer (photosphere/chromosphere), and a hot isothermal layer (corona) (see Figure 1(a)). The initial distribution of temperature is given as

$$\begin{eqnarray}&&T(z)={T}_{\mathrm{ph}}-{\left|\displaystyle \frac{{dT}}{{dz}}\right|}_{\mathrm{ad}}z\end{eqnarray} \tag{ 8 }$$

for the convection zone ($0\leqslant z\leqslant {z}_{\mathrm{ph}}$), and

$$\begin{eqnarray}&&T(z)={T}_{\mathrm{ph}}+({T}_{\mathrm{co}}-{T}_{\mathrm{ph}})\left\{\displaystyle \frac{1}{2}\left[\mathrm{tanh}\left(\displaystyle \frac{z-{z}_{\mathrm{tr}}}{{w}_{\mathrm{tr}}}\right)+1\right]\right\},\end{eqnarray} \tag{ 9 }$$

for the upper atmosphere ($z\leqslant {z}_{\mathrm{ph}}$), where T_ph and T_cr are the temperatures in the photosphere and the corona, and are set to T₀ and 150T₀, respectively. $| {dT}/{{dz}| }_{\mathrm{ad}}\equiv (\gamma -1)/\gamma$ is the adiabatic temperature gradient. The photospheric height is ${z}_{\mathrm{ph}}=80{H}_{p0}$, and the transition region between the chromosphere and the corona is located at ${z}_{\mathrm{tr}}={z}_{\mathrm{ph}}+18{H}_{p0}=98{H}_{p0}$. The width of the transition w_tr is set to 2H_p0. A magnetic flux tube is initially located below the photosphere (the depth is $\sim 10$ Mm). The longitudinal and azimuthal components of the flux tube are respectively given by

$$\begin{eqnarray}&&{B}_{x}(r)={B}_{\mathrm{tube},0}\mathrm{exp}\left(-\displaystyle \frac{{r}^{2}}{{R}_{\mathrm{tube}}^{2}}\right),\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&{B}_{\theta }(r)=\alpha \displaystyle \frac{r}{{R}_{\mathrm{tube}}}{B}_{x}(r),\end{eqnarray} \tag{ 11 }$$

where $r={\left[{y}^{2}+{(z-{z}_{\mathrm{tube}})}^{2}\right]}^{1/2}$, ${z}_{\mathrm{tube}}=18{H}_{p0}$, and ${B}_{\mathrm{tube},0}=53{B}_{0}$. ${R}_{\mathrm{tube}}=5{H}_{p0}$ is the radius of the tube. These profiles are truncated at $r=3{R}_{\mathrm{tube}}$. The total magnetic flux of $3\times {10}^{20}$ Mx is assumed. The plasma-β at the center of the initial tube is approximately 20. α is a measure of the twist of the flux tube, and is set to $-1.5$ (negative sign denotes the left-handed twist). Note that α is defined as a non-dimensional constant here. $\alpha /{R}_{\mathrm{tube}}$ gives the rate of field line rotation per unit length along the tube. The absolute value of α is larger than the critical value for the kink instability, 1, and therefore the flux tube is initially kink-unstable. The vertical distributions of the density, pressure, temperature, and magnetic field strength in the initial condition are shown in Figure 1(b).

Zoom In Zoom Out Reset image size

The simulation domain is $(-124{H}_{p0},-124{H}_{p0},0)$ $\leqslant \;(x,y,z)\;\leqslant \;(124{H}_{p0},124{H}_{p0},294{H}_{p0})$ resolved by a $600\times 600\times 740$ grid. In the x and y-directions, the mesh size is ${\rm{\Delta }}x={\rm{\Delta }}y=0.25{H}_{p0}$ within $| x| ,| y| \lt 30{H}_{p0}$ and gradually increases up to 0.8H_p0 for $| x| ,| y| \gt 30{H}_{p0}$. In the z-direction, the mesh size is ${\rm{\Delta }}z=0.24{H}_{p0}$ within $0\leqslant z\leqslant {z}_{\mathrm{ph}}+40{H}_{p0}$, and for $z\gt {z}_{\mathrm{ph}}+40{H}_{p0}$, it gradually increases up to 1.6H_p0. To facilitate investigation of the coupling between the kink and Parker instabilities, we take a domain size large enough to allow the latter to develop. The critical wavelength for it is $\sim 9$ times the local pressure scale height, ∼230H_p0, which is comparable to or smaller than the calculation domain size in the x-direction, 248H_p0.

We make the flux tube buoyant by introducing a small perturbation to the density in the flux tube. The perturbation is described as

$$\begin{eqnarray}&&\rho ={\rho }_{0}\left[1-a(x)\right],\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&a(x)=\displaystyle \frac{1}{\beta }\left[(1+\epsilon )\mathrm{exp}({-x}^{2}/{\lambda }^{2})-\epsilon \right],\end{eqnarray} \tag{ 13 }$$

where ρ₀ is the unperturbed density, $\epsilon =0.2$ and $\lambda =15{H}_{p0}$. We have checked that the evolution of the rising flux tube is not sensitive to the wavelength of the initial perturbation.

The boundaries in the x-direction (the direction of the initial flux tube axis) are assumed to be periodic. The boundaries in the y and z-directions are assumed to be a perfectly conducting wall and are non-penetrating. Many previous studies adopted a fee boundary condition for the top boundary, but we applied the non-penetrating boundary for it to avoid numerical instabilities arising from the use of a free boundary condition. The top boundary is located at a much higher position ($z=294{H}_{p0}$) than the top of the emerging coronal loops ($z\sim 200{H}_{p0}$), and therefore the top boundary will not significantly affect the evolution of the flux emergence process.

The rise of the flux tube with a single buoyant segment is shown in Figure 2. Since the tube is initially kink-unstable, the knotted structure develops during the rise, as in the previous study by Fan et al. (1999). In this study, the box size in the x-direction is $\sim 9$ times larger than the pressure scale height at the initial tube axis, whereas they are similar in Fan et al. (1999). With the larger domain size, the rising tube develops a Ω-shape with a central kinked part.

Zoom In Zoom Out Reset image size

To examine the deformation and expansion of the apex of the rising tube, we looked at the relation between the plasma density and the magnetic field strength of the most buoyant part. We confirmed that the magnetic field strength of the most buoyant part can be well described by $B\propto {\rho }^{1/2}$ (see Figure 3), which means that the flux tube experiences a strong horizontal expansion during its rise to the surface (Cheung et al. 2010). The strong horizontal expansion leads to the formation of a pancake-like (flat, horizontally extended) magnetic distribution below the surface (see Figure 2). Because of the deformation by the kinking and the strong horizontal expansion, it is difficult to deduce the resulting photospheric and coronal magnetic structures only from the magnetic structure of the rising kinked flux tube in the interior of the Sun.

Zoom In Zoom Out Reset image size

The evolution of the vertical components of the photospheric ($z=80{H}_{p0}$) and chromospheric ($z=96{H}_{p0}$) magnetic field, B_z, is displayed in Figure 4. At the beginning of the emergence ($t=305\tau$), a pair of opposite polarity spots (we call them the "main pair") appear with a large inclination with respect to the initial axis direction (x-direction). As time progresses, another pair of strong magnetic regions appear at the middle of the main pair of spots (called the "middle pair"). Later, the structure of the middle pair becomes disordered, although the main spots show a coherent structure and strong twist. The main pair and middle pair have the maximum field strength of ∼10B₀ and ∼3–4B₀, respectively.

Zoom In Zoom Out Reset image size

In order to describe the motion of this quadrupole region, we measured the flux-weighted centroids of positive and negative polarities $({x}_{\pm },{y}_{\pm })$, where

$$\begin{eqnarray}&&({x}_{\pm },{y}_{\pm })=\left(\displaystyle \frac{{\rm{\Sigma }}{{xB}}_{z,\pm }}{{\rm{\Sigma }}{B}_{z,\pm }},\displaystyle \frac{{\rm{\Sigma }}{{yB}}_{z,\pm }}{{\rm{\Sigma }}{B}_{z,\pm }}\right).\end{eqnarray} \tag{ 14 }$$

Relative motion of the positive and negative polarities at the photosphere is shown in Figure 5. The centroid position for each polarity is computed including all the parts of that polarity (i.e., include both the main pair and the middle pair). The tilt angle of this region, measured from the positive x-direction, is large at the beginning of the emergence (note that the direction of the initial tube axis is in the x-axis). However, it becomes smaller later. This change may be regarded as a clockwise motion of this region, as predicted by the emergence of a knotted tube with a left-handed twist (Tanaka 1991; Linton et al. 1999). We note that the distance between the flux-weighted centers of the opposite polarities increases as time progresses.

Zoom In Zoom Out Reset image size

It should be noted that a complex quadrupole structure is formed by the emergence of a kinked tube with a single buoyant segment, not with two buoyant segments assumed in the previous studies by Toriumi et al. (2014) and Fang & Fan (2015). The narrow middle pair becomes prominent well after the emergence of the main pair, and its structure gets disordered, but the chromospheric B_z maps in Figure 4 show a more smooth quadrupole distribution than the photospheric distribution.

We investigated the formation of the middle pair. Figure 6 shows the temporal evolution of the vertical velocity v_z at the center $(x,y)=(0,0)$ during the formation. Magnetic field is rising in the early phase ($t=310\tau$), but later v_z becomes negative, which means the submergence of emerged magnetic fields. Figure 7(a) displays the 3D magnetic field evolution, where it is shown that the middle pair is formed as emerged magnetic fields submerge. For this reason, the middle pair appears after the formation of the main pair. To clarify the cause of the submergence, we measured the vertical forces at the center. Figure 7(b) indicates that the sum of the upward Lorentz force and pressure gradient becomes weaker than the downward gravitational force, which means that the submergence is caused by the downward motion of the heavy material.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figures 8(a)–(c) illustrate the evolution of the magnetic field and density distribution in the y–z plane at x = 0. Two magnetic flux concentrations are located just below the photosphere at the beginning of the emergence ($t=300\tau$) with a single arcade field emerging into the atmosphere. In this model, as in many previous numerical models of emerging flux, the expansion of the magnetic arcade into the atmosphere is enabled by a magnetic Rayleigh–Taylor instability (e.g., Shibata et al. 1989; Matsumoto et al. 1993; Archontis et al. 2004). As the expanding arcade plows through the atmosphere, plasma is compressed above the arcade, leading to a top heavy layer (see Figure 8 panel (d)). Due to the continued emergence of the two strong flux concentrations, and the comparatively weaker twist at the middle between the two concentrations, the top heavy layer is driven to accumulate plasma at a central dip (see Figure 8 panels (d) and (e)). The continued plasma accumulation results in the submergence of the heavy plasma at the dip and the formation of the two adjacent magnetic arcades.

Zoom In Zoom Out Reset image size

A snapshot of the 3D coronal magnetic structure is shown in Figure 9. The current sheet indicated in Figure 8(c) is located between the blue and yellow magnetic arcades, where magnetic reconnection takes place. Two sets of the new loops colored purple and white are interpreted as reconnected field lines (Figures 9(a)–(d)). Looking at the reconnection site, we see that the reconnection angle (the angle between the merging field lines) is not 180° (i.e., not perfectly antiparallel), and merging field lines have a large guide field (panel (e)). As shown in Figures 9(d) and (f), the lower new arcade (purple lines) is almost parallel to the PIL, and is connected to the middle pair. The upper new arcade (white lines), which was ejected upward from the reconnection site, connects the far ends of the two magnetic arcades (yellow and blue).

Zoom In Zoom Out Reset image size

We also note the magnetic connectivity above the photosphere shown in Figure 9. A fraction of the magnetic field of the main spots is connected to the middle pair. The other magnetic field connects the two main spots. A portion of the magnetic field connecting the main spots is the reconnected field (white lines in Figure 9).

We have seen that in Figure 9, both the yellow and blue arcades on two sides of the PIL as well as the (purple) reconnected loops connecting the middle pair are highly sheared along the PIL of the middle pair. When we look at the photospheric motion, we find the vortical or rotational motion within each of the two main polarities (Figure 10(a)). Figure 10(b) shows the temporal evolution of the vertical vorticity ω_z averaged over the area where $| {B}_{z}|$ is above 75% of the peak $| {B}_{z}|$ value. From the figure, we find that the counterclockwise vortical motion becomes prominent at time $t\sim 320\tau$, and persists throughout the subsequent evolution. To illustrate the development of magnetic shear, we show the temporal evolution of the vector magnetogram of the middle pair at the photosphere and middle chromosphere in Figure 11. We can clearly observe that the horizontal magnetic field becomes more parallel to the PIL at those heights as time progresses. Note that the horizontal magnetic field strength is small just on the PIL. This is because the middle pair is formed as a result of the submergence, not the emergence. As investigated by Fan (2009), the vortical motion (and the horizontal motion) of the main polarities builds-up magnetic shear in this simulation. We also note that the magnetic shear near the PIL is not due to the direct emergence of a sheared field.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The middle pair is pushed and confined in a narrow region after the formation (Figure 4). To understand the confinement mechanism, we looked at the flow structure and the Lorentz force at the photosphere. The top panel of Figure 12 shows horizontal velocity vector on the magnetogram, and we can see the converging flow to the PIL. This converging flow is a natural consequence of the development (expansion) of the two magnetic arcades. The bottom panel of Figure 12 displays the Lorentz force vector, and we can find that the horizontal Lorentz force is pushing the two polarities together. We checked the pressure balance across the PIL of the middle pair. Figure 13 displays the profiles of the total pressure, gas pressure, magnetic pressure, and dynamic pressure across the PIL. The total pressure is almost constant across the PIL, which will explain the persistent existence of the middle pair. We can see that the high gas pressure region is supported by the magnetic pressure. We also found that the magnetic pressure has its peaks at the edges of the high pressure region because of the dynamic pressure of the converging flows.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The horizontal magnetic field is important for producing the large magnetic pressure near the PIL. As shown in Figure 10, the main polarities show the vortical motion to shear the field. Figure 14(a) displays the development of the horizontal field along the PIL (see also Figure 11). Figure 14(b) shows the ratio of the horizontal field ${B}_{h}=\sqrt{{B}_{x}^{2}+{B}_{y}^{2}}$ to the vertical magnetic field B_z. It is found that the horizontal field is much stronger than the vertical field outside the middle pair and the PIL is sandwiched in the strong horizontal field regions. Therefore, the development of the magnetic shear along the PIL by the vortical motion of the main polarities is a key to confining the narrow middle pair.

Zoom In Zoom Out Reset image size

We also found persistent fast shear flows along the PIL (see the top panel of Figure 12). The maximum velocity is about 2C_s0, which is supersonic at the photospheric level. To understand the driving mechanism, we investigated the acceleration by the Lorentz force and the pressure gradient force. Considering that the PIL is roughly straight and almost parallel to the unit vector $\hat{{\boldsymbol{e}}}=(-1/\sqrt{2},1/\sqrt{2})$, we made the dot products of the two force vectors and this unit vector to see the acceleration in the direction of the PIL: ${f}_{{\rm{L}},\mathrm{PIL}}\equiv ({\boldsymbol{J}}\times {\boldsymbol{B}})\cdot \hat{{\boldsymbol{e}}}$ and ${f}_{{\rm{P}},\mathrm{PIL}}\equiv (-{\rm{\nabla }}p)\cdot \hat{{\boldsymbol{e}}}$, respectively. Figure 15 displays the relation between the two forces and horizontal flows. The color indicates ${f}_{{\rm{L}},\mathrm{PIL}}$ (Top) and ${f}_{{\rm{P}},\mathrm{PIL}}$ (Bottom). Note that positive (negative) ${f}_{{\rm{L}},\mathrm{PIL}}$ accelerates plasma in the upper-left (lower-right) direction. The same is true for ${f}_{{\rm{P}},\mathrm{PIL}}$. The horizontal velocity vector shows that the converging flows drastically change their direction to the direction of the PIL at the edges of the middle pair. From the figure, it is found that the shear flows along the PIL are driven by the Lorentz force. The pressure gradient force has no significant contribution (this implies that the hydrodynamic baroclinic vorticity generation ($-{\rm{\nabla }}(1/\rho )\times {\rm{\nabla }}p$) can be neglected). The acceleration regions coincide with the strong magnetic pressure regions at the edges of the middle pair, where the plasma-β is less than unity (see Figure 13). The direction of the shear flow is determined by the expansion of each magnetic arcade: the magnetic arcade with the positive (negative) main polarity drives the shear flows in the positive (negative) x-direction.

Zoom In Zoom Out Reset image size

We briefly compare our simulation with an observational example of δ-spot regions. We focus on Active Region NOAA 11429, which was one of the most violent δ-spot regions in solar cycle 24. The δ-spot region appeared in 2012 March, and produced three X-class flares. Because of the high activity, it have been drawing many authors' attention (e.g., Petrie 2012; Shimizu et al. 2014). Figure 16 displays snapshots of the active regions and snapshots of the magnetogram obtained from our simulation (note that the sign of B_z of the magnetogram from the simulation is reversed just for better comparison). The continuum and magnetogram images were taken by the Helioseismic and Magnetic Imager (HMI: Schou et al. 2012) on board the Solar Dynamics Observatory. The active region appeared in the northern hemisphere, and did not follow Joy's Law. The mean force-free parameter $\langle \alpha \rangle$ of this region calculated in SHARP (Spaceweather HMI Active Region Patch) data was negative, which means that the region obeyed the helicity hemispheric rule. This is consistent with the apparent magnetic structure: the largest positive and negative spots show a left-handed twist. Our simulation reproduced some observed features of this δ-spot region. The observed δ-spot region has a narrow complex polarity pair between the main polarities (indicated by the arrow in the figure), and therefore it had a complex quadrupole structure. In addition, the middle pair appeared after the emergence of the main polarities. These observed features are also found in our simulation.

Zoom In Zoom Out Reset image size