Formation of super-strong horizontal magnetic field in delta-type sunspot in radiation magnetohydrodynamic simulations

Our calculation procedure is shown in Fig. 1. We start our calculations from the low-resolution deep convection zone run that lacks the photosphere with the top boundary at 0.99$R_{\odot}$, resolved by the number of grid points of $256^{3}$ (Case DEEP). The computation speed is accelerated by omitting the photosphere because we can circumvent the calculation of small-scale, short-lived, granular convection. Case DEEP continues for 90 days. The typical convection time scale in the deep convection zone is 30 days and 90 days corresponds to three turn over times. However, Hotta et al. (2019) show that calculations for 75 days and 360 days yield almost the same stratification. Thus we conclude that 90 days is enough for relaxation.

The calculation domain is extended to include the photosphere (Case HYDRO). From here, the top boundary is at $700~{}\mathrm{km}$ above the average $\tau=1$ surface. The resolution is upgraded to $384\times 1024^{2}$. The horizontal grid spacing is 96 km, which is acceptable for resolving the photosphere. We use non-uniform grid spacing in the vertical direction, which is 48 km and 1.5 Mm at the photosphere and the base of the convection zone, respectively. This grid spacing of 1.5 Mm is acceptable for resolving the deep convection zone, where the pressure scale height is 60 Mm. Case HYDRO continues for five days. Then, statistically steady convection covering the whole convection zone is prepared. This approach is justified because the near-surface layer does not have a significant influence on the deep convection (Hotta et al., 2019).

From Case FE1, we start the flux emergence simulations. The force-free twisted flux tube is inserted at a depth of 22 Mm. We use the Bessel function for the force-free flux tube (see HI20 for the details). The axial magnetic field at the centre of the flux tube and the radius of the flux tube are $10^{4}~{}\mathrm{G}$ and 10 Mm, respectively. The total magnetic flux is $1.35\times 10^{22}~{}\mathrm{Mx}$, which is twice larger than that in TH19 to have stronger magnetic field in the photosphere. We define the beginning of Case FE1 as $t=0$. In this case, we only solve two rays for the radiation transfer; thus, only upward and downward energy transfers are allowed. Case FE1 continues for 92.7 h. The initial location of the flux tube is shown in Fig. 2. The coherent large-scale downflow at the centre of the calculation domain pins down the middle of the flux tube, whereas the two ends emerge to the photosphere. As a result, the sunspots of opposite polarities collide with each other to generate a delta-type spot in a self-consistent manner (the multi-buoyant segment scenario: see Toriumi et al., 2014, and TH19 for the details).

After the initial sunspots are created at $t=33.3~{}\mathrm{h}$, we upgrade the calculation to Case FEM, in which the number of grids is $256\times 2048^{2}$. The horizontal grid spacing is 48 km. We keep the 48 km vertical grid spacing in Case FEM in the photosphere, but the grid spacing at the base of the convection zone is made coarse by using a spacing of 3 Mm because the details of the convection in the deep layer do not have a significant influence on the sunspot evolution in the photosphere. From here, we use the multi-ray radiation transfer. We solve 24 rays (see Appendix A for the details of the calculation scheme). The one-ray radiation transfer prohibits the horizontal radiation energy transfer, and the convection pattern tends to show unrealistic small-scale features in high-resolution calculations. In addition, we adopt the potential magnetic field at the top boundary ($\alpha=1$). In this study, we mainly analyse the Case FEM.

By using the data at $t=43.3~{}\mathrm{h}$ of the Case FEM, we restart Cases FEM80, FEM160, and FEM$\alpha$. In Cases FEM80 and FEM160, we change the Alfven speed limit to $C_{\mathrm{max}}=80$ and $160~{}\mathrm{km~{}s^{-1}}$, respectively. In Case FEM$\alpha$, we adopt $\alpha=2$, in which the magnetic field at the top boundary is more horizontal than the potential magnetic field. Cases FEM80 and FEM160 continue for 10 h, whereas Case FEM$\alpha$ continues for 14 h. In these cases, the other settings are the same as those of Case FEM.

The data at $t=51.3~{}\mathrm{h}$ of Case FEM are upgraded to a higher resolution and restarted as Case FEMH. The number of grid points in Case FEMH is $384\times 4096^{2}$. The grid spacing in the photosphere is 24 km in all three directions. We still use the non-uniform vertical grid spacing, which is 3 Mm at the base of the convection zone. Case FEMH continues for 2 h. Fig. 3 shows the emergent intensity in Case FEMH.

The change of resolution implies the difference in the numerical diffusivity. With a higher resolution, we can capture small-scale features. Since the difference of the resolution mainly influences the small scales, which have much shorter time scales, the relaxation occurs in much smaller time scales than the typical convection time scale, which is several minutes in the photosphere.

In this subsection, the overall evolution of the flux tube and the generated sunspots are explained. Fig. 4 shows the evolution of the unsigned magnetic flux at the $\tau=1$ surface. The solid and the dashed lines indicate the results from Cases FEM and FE1, respectively. The black lines are the unsigned magnetic flux in the whole computational domain. The red and blue lines show the magnetic flux in the regions with the intensities less than 80% and 50% of the averaged photospheric intensity $I_{\odot}$, respectively, which roughly correspond to the penumbral and umbral regions. The maximum unsigned magnetic flux is almost $5\times 10^{22}~{}\mathrm{Mx}$, which is well within the solar range ($\lesssim 2\times 10^{23}$: Toriumi et al., 2017), indicating that our setup of the magnetic flux is applicable to the Sun. We note that the magnetic flux reaches the steady state after $t=60~{}\mathrm{h}$, i.e. the magnetic flux does not decrease. This plateau indicates that the amount of photospheric flux in the delta-spot simulations does not decrease as fast as in the other flux emergence simulations (see HI20). However, we do not discuss the spot decay in detail in the present study. The calculation duration of Case FEM is shorter than that of Case FE1 but covers the period of the peak unsigned magnetic flux ($t=45~{}\mathrm{h}$). The overall evolution of the magnetic flux is similar between Cases FE1 and FEM, but we see a distinct difference, especially in the blue lines ($I<0.5I_{\odot}$), corresponding to the umbra. This indicates that the resolution and/or the radiation transfer treatment can affect the umbral evolution.

Fig. 5 shows the evolution of the three-dimensional structure of the flux tube and the normalized entropy in Case FE1. The left panels show the magnetic field strength $|B|$, and the right panels are the normalized entropy $(s-\langle s\rangle)/s_{\mathrm{rms}}$, where $s$ is the specific entropy. $\langle s\rangle$ and $s_{\mathrm{rms}}$ are the horizontally averaged and root-mean-square (rms) entropy for Case HYDRO, respectively. The results show a large-scale coherent downflow in the centre of the calculation domain that extends to the deeper layer ($z<-150~{}\mathrm{Mm}$, red clump at the center of the calculation domain extending from the photosphere to the deep layer, see also the movie). Because the convection time scale in the deeper layer is long ($\sim 30~{}\mathrm{d}$), the overall structure of this downflow does not change during the calculation period of Case FE1 (Figs. 5b, d, and f). This long-lived coherent downflow is the essential origin of the delta-type spot in this study. Upflow beside the coherent downflow causes the rising of the flux tube and the formation of the sunspot in the photosphere. The coherent downflow in the deep region causes the convergent flow between the two sunspots, and they approach each other (Fig. 5c). Then, the two sunspots collide to form a delta-type sunspot. We note that the force-free state is broken as the flux tube starts to rise. The density decreases in the flux tube and the gas pressure gradient plays a significant role in maintaining the force balance of the flux tube (see also HI20)

Fig. 6 shows the emergent intensity and the magnetic field on the $\tau=0.01$ surface in Case FEM. When the two sunspots are separated from each other, the sunspot magnetic field is mostly vertical, and the horizontal magnetic field strength is only a few thousand Gauss at maximum (Figs. 6a and b). As the two sunspots approach, penumbra-like structures are constructed between them (Fig. 6c). The horizontal magnetic field perpendicular to the PIL becomes intensified (Fig. 6d). At $t=52~{}\mathrm{h}$, the two sunspots reach each other, and the horizontal magnetic field exceeds 6000 G.

Fig. 7 presents a comparison between Cases FEM, FEM$\alpha$, and FEMH. Figs. 7a, c, and e indicate the emergent intensity; and b, d, and f depict the vertical velocity on the $\tau=1$ surface. All panels show the quantities at $t=53.3~{}\mathrm{h}$. While the photospheric surface is filled with fine-scale structures in Case FEMH (panels e and f), the overall structure, such as the width of the light bridge between the two sunspots, is similar to that of the lower-resolution run (Case FEM: panels a and b). The typical width of the filamentary structure across the PIL is 1 Mm for both Cases FEM and FEMH. One notable difference between Cases FEM and FEMH is the intensity in the umbra. Case FEMH (Fig. 7e) shows lower intensity in the umbra than that in Case FEM (Fig. 7a). Rempel (2012) suggests that the numerical diffusivity causes a mass diffusion and the resulting increase of the umbral dot. Case FEMH has less numerical diffusivity than Case FEM, which probably is the cause of the reduced umbral intensity in Case FEMH. This result implies that the actual Sun would show smaller umbral intensity than the simulations shown in this paper since the actual Sun has much smaller diffusivity. The vertical velocity structures in Cases FEM and FEMH are almost identical. Along the light bridge, one may find the region of weak upflows at the edges of the sunspot on the right. The corresponding region of the left sunspot is filled with downflows. Some fluctuations may be found in this flow structure, but the upflow–downflow relation does not change along the light bridge.

Case FEM$\alpha$ (Figs. 7c and d) exhibits significantly different features from Cases FEM and FEMH. Specifically, the width of the light bridge in Cases FEM and FEMH is 2–3 Mm, which is larger than 5 Mm in Case FEM$\alpha$. In response to this change, the area of the umbrae is decreased. Rempel (2012) reveals that the increase of $\alpha$ extends the length of the penumbra, but the umbral area is not influenced significantly in his comparison (see Figs. 2 and 3 of Rempel, 2012). We expect that the difference in topology, i.e. whether the sunspots are isolated or delta type, causes this decrease in the umbral area. Another difference in Case FEM$\alpha$ from Cases FEM and FEMH is found with regard to the vertical velocity, especially in the light bridge (Fig. 7d). While the flow structure is coherent in Cases FEM and FEMH, the flow in Case FEM$\alpha$ is somewhat messy. The left (right) edge of the light bridge is mostly occupied by the downflows (upflows). However, in some segments, this rule is violated and the upflow (downflow) is observed.

Apart from the light bridge, coherent downflows are also observed at the immediate edges of the umbrae (rather than at the penumbra/quiet-Sun boundaries, which is regularly seen). Rempel (2011) shows that the downflows at the umbral edges are created when the spot lacks a penumbra. Without a penumbra, the deeper layer around the umbra is exposed and the radiative energy loss becomes significant, leading to the production of such downflows. In our simulations, however, the downflows at the umbral edges are observed even in the case with prominent penumbrae (Case FEM$\alpha$). Therefore, we conjecture that the spot rotation causes these downflows. Because of the Lorentz force, the horizontal rotational flow is bent to the vertical direction and the downflow is driven there.

Fig. 8 shows the horizontal magnetic field between two sunspots at $t=53.3~{}\mathrm{h}$ for Case FEM. Figs. 8a and b indicate the emergent intensity and the horizontal magnetic field strength $B_{\mathrm{h}}=\sqrt{B_{x}^{2}+B_{y}^{2}}$ at $\tau=0.01$, respectively. The maximum horizontal strength in this period at $\tau=0.01$ and 1 exceeds 6600 G and 7200 G, respectively. Figs. 8c and d present the profiles of the horizontal magnetic field strength along the red dashed line in Figs. 8a and b, showing the results at $\tau=0.01$ and $1$, respectively. Because the density is larger on the $\tau=1$ surface than that at $\tau=0.01$, the magnetic field strength is slightly stronger on the $\tau=1$ surface. Even at $\tau=0.01$, the horizontal magnetic field strength exceeds 6000 G in all cases. Cases FEM (black), FEM80 (yellow), FEM160 (blue), and FEMH (purple) show almost the same profile of the horizontal magnetic field strength. Case FEM$\alpha$ (green) shows a slightly stronger field around $45~{}\mathrm{Mm}<y<60~{}\mathrm{Mm}$, but the peak magnetic field around $y=40~{}\mathrm{Mm}$ is almost the same as for the other cases except for Case FE1. This finding makes a marked contrast with the result in Rempel (2012). Rempel (2012) shows that a larger $\alpha$ causes a large horizontal magnetic field. While this is true in our calculations, the horizontal field of the light bridge does not depend on the top boundary condition. Case FE1 displays a slightly different magnetic field strength in some positions. In Case FE1, the calculation cannot properly resolve the fine-scale filamentary structures, potentially leading to this slight difference.

The consistency of the results with different numerical conditions indicates the robustness of our result, especially of the exceptionally strong magnetic fields. The fact that the different Alfven speed limits do not influence the peak strength of the horizontal magnetic field indicates that these magnetic fields are in a force-free state. Fig. 9 shows a quantity:

on $z=0$. Because the Lorentz force $\bm{F}_{\mathrm{L}}$ is expressed as $\bm{F}_{\mathrm{L}}=\left(\nabla\times\bm{B}\right)\times\bm{B}/4\pi$, the force is not acting when the electric current and the magnetic field are parallel. If $\cos\theta=1$, the current and the magnetic field are parallel. Fig. 9 shows that the force-free state is achieved within the light bridge, i.e. the strong horizontal field region.

In this subsection, we investigate the amplification and maintenance mechanism of the strong horizontal magnetic field in the light bridge.

Figs. 10a, b, and c show the flow structures $v_{x}$, $v_{y}$, and $v_{z}$ at $z=0$, respectively. $z=0$ is at approximately the same level as the $\tau=0.01$ surface at the light bridge. Because the initial flux tube is twisted in a right-handed manner to achieve the force-free state, the sunspots in this phase are rotating clockwise by relaxing the twist. The flow structure in the light bridge reflects this rotation. At the edge of the left (right) sunspot, the horizontal flow is directed to the lower left (upper right). This means that the light bridge is filled with a divergent and sheared flow. The vertical flow in the centre of the light bridge is a weak upflow and is surrounded by the downflow. These vertical flows are also caused by the rotation.

To understand the evolution of the magnetic field, we investigate the induction equation. The induction equation is written as

where $i=x,y,z$. The terms $P_{\mathrm{adv}(i)}$, $P_{\mathrm{cmp}(i)}$, and $P_{\mathrm{str}(i)}$ denote the advection, compression, and stretching, respectively. These are expressed as

In this study, we omit $B_{i}\partial v_{i}/\partial i$ from $P_{\mathrm{cmp}(i)}$ and $P_{\mathrm{str}(i)}$ because $\partial v_{i}/\partial i$ has no influence on the evolution of the $i$-component magnetic field ($B_{i}$). We note that these terms are cancelled out when we discuss the sum of the terms.
Fig. 11 shows the terms in the $x$- and $y$-component induction equation along the green dashed line in Fig. 10. The black, red, and blue lines represent $P_{\mathrm{adv}(i)}$, $P_{\mathrm{cmp}(i)}$, and $P_{\mathrm{str}(i)}$, respectively. The dashed line indicates the sum of these terms. The values are averaged over the period of $t=48.3$ to $53.3~{}\mathrm{h}$. While the magnetic field is evolving during this period, the sum of the terms is much smaller than the individual terms. The amplification of the magnetic field takes about 10 h (from Fig. 6d to Fig. 6h). Since the amplified magnetic field strength is 6000 G, the amplification rate of the magnetic field is about 0.2 $\mathrm{G\ s^{-1}}$, which is two orders of magnitude smaller than each term in the induction equations. Therefore, the best we can do is to determine the contributing term(s) for the evolution (see also Siu-Tapia et al., 2019). Regarding the $x$-component induction equation, whose individual terms are plotted in Fig. 11a, the contributing factor is not clear. We cannot determine which term is important for the amplification and maintenance of the $x$-component field. This is discussed further in the next paragraph. Meanwhile, Fig. 11b shows that the $y$-component magnetic field $B_{y}$ is amplified by the stretching. We note that the amplified $y$-component magnetic field $B_{y}$ at the light bridge is negative and thus the negative values for the terms correspond to the field amplification. Interestingly, the compression term weakens the horizontal field $B_{y}$, as this term is positive, i.e. the flow is diverging.

We perform further analysis on the evolution of $x$-component magnetic field $B_{x}$ with Fig. 12. Figs. 12a, b, and c show $B_{x}$, $P_{\mathrm{tra}(xz)}=P_{\mathrm{adv}(xz)}+P_{\mathrm{cmp}(xz)}$, and $P_{\mathrm{str}(xz)}$, respectively, where

These terms correspond to the vertical transport of $B_{x}$ and the stretching of the vertical field $B_{z}$ in the direction of $B_{x}$, respectively. We find that these terms play the most important roles in enhancing the magnitude of $B_{x}$. In the centre of the light bridge, the vertical transport term contributes to increase $B_{x}$, while at the edge of the light bridge, the vertical stretching creates $B_{x}$. Both effects indicate that the vertical magnetic fields from the two sunspots reconnect with each other in the upper atmosphere above the light bridge (i.e. above the PIL) and create arcade-shaped field lines. Then, the amplified $B_{x}$ is transported downward and observed in the lower atmosphere at $\tau=0.01$ or $z=0$. It should be noted that the reconnection takes place probably above the top boundary because of the potential boundary condition and thus we do not directly observe the site of reconnection.

Fig. 13 shows the contribution to the $y$-component induction equation. Figs. 13a, b, and c describe $B_{y}$, $P_{\mathrm{str}(yx)}$, and $P_{\mathrm{cmp}(yx)}$, respectively, where

These terms correspond to the stretching of the $B_{x}$ field in the direction of $B_{y}$ and the compression of the $x$-component flow, respectively. Again, these terms are found to be most critical to the evolution of $B_{y}$. Fig. 13b shows that $B_{y}$ in the centre of the light bridge is mainly amplified by the shear of $v_{y}$ (see also Fig. 10b). Because there is a divergent flow in the centre of the light bridge, the compression term decreases $B_{y}$. In contrast, at both edges, the compression plays a primary role in amplifying the magnetic field.

Our finding is confirmed by the vertical slice at $y=45~{}\mathrm{Mm}$ (Fig. 14). The red arrows show the magnetic field on the $x$–$z$ plane ($B_{x}$ and $B_{z}$). The greyscale represents $B_{y}$. The solid and dashed lines are the $\tau=1$ and $\tau=0.01$ surfaces, respectively. As discussed, the arcade-shaped magnetic field across the PIL is seen in the upper layer above $\tau=0.01$. The strong $B_{y}$ is created around or below the photosphere down to $-5$ Mm. As a result, a sheared arcade field is created above the PIL. The established force-free state in this region (Fig. 9) indicates that the magnetic pressure gradient of the axial field $B_{y}$ and the magnetic tension of the other components $B_{x}$ and $B_{z}$ are balanced.

Fig. 15 shows the horizontal flow ($v_{x}$ and $v_{y}$: the red arrows) and the vertical magnetic field ($B_{z}$: colour contour) at $z=-10~{}\mathrm{Mm}$ (panel a) and $z=-20~{}\mathrm{Mm}$ (panel b). At $z=-10~{}\mathrm{Mm}$, both the left and right sunspots show clockwise rotations, which is consistent with the photosphere. In the deeper layer ($z=-20~{}\mathrm{Mm}$), only the right sunspot appears to rotate.