On rising magnetic flux tube and formation of sunspots in a deep domain

We solve the three-dimensional magnetohydrodynamic (MHD) equations in the Cartesian geometry $(x,y,z)$ with the radiation transfer, where the $z$-direction is vertical and the $x$ and $y$ directions are horizontal. In this study, $z=0$ is the solar surface, which is $R_{\odot}=696\ \mathrm{Mm}$ from the center of the sun. The equations are solved with the R2D2 code. The magnetohydrodynamic equations are expressed as:

where $\rho$, $v$, $B$, $p$, $T$, $s$, $g$, and $Q$ are the density, fluid velocity, magnetic field, gas pressure, temperature, entropy, gravitational acceleration, and radiative heating, respectively. The subscript 1 indicates the perturbation from the stationary 1D stratification indicated with subscript 0. Thus, the thermodynamic variables are expressed as:

In this study, we do not assume that the perturbation is smaller than the stationary background stratification, i.e., $\rho_{1}<<\rho_{0}$ is not assumed. The background stratification is calculated with the hydrostatic equation with the help of the Model S (Christensen-Dalsgaard et al., 1996). See Appendix A for the details of the calculation procedure.

The equations are solved with the fourth-order spatial derivative (see Appendix B for details) and the four-step Runge–Kutta method for the time integration (Vögler et al., 2005). To maintain the stability of the calculation, we adopt the slope-limited diffusion suggested by Rempel (2014).

We use the equation of state considering the partial ionization effect with the OPAL repository (Rogers et al., 1996). We switch the linear and table equations of state to address the significant change of the perturbation through the convection zone. We evaluate value

If $c_{\mathrm{eos}}$ exceeds $10^{-2}$, we use the table equation of state ; otherwise, the linear equation of state,

is used. $(\partial p/\partial\rho)_{s}$ and $(\partial p/\partial s)_{\rho}$ are prepared with the background stratification and only depend on the height ($z$). Regarding the table equation of state, we prepare $64\times 64$ grid on the density and entropy.

We adopt the RSST (Hotta et al., 2012; Hotta et al., 2015; Iijima et al., 2019) for the equation of continuity to deal with the low Mach number flow in the deep convection zone. By using the RSST, we try to keep the Mach number throughout the convection zone. To this end, we set the RSST factor,

Here, we adopt $\xi_{0}=160$. $\rho_{\mathrm{b}}=0.2\mathrm{\ g\ cm^{-3}}$ and $c_{\mathrm{b}}=2.2\times 10^{7}\ \mathrm{cm\ s^{-1}}$ are the density and the speed of sound around the base of the convection zone, respectively. $c_{\mathrm{s}}=\sqrt{\left(\partial p/\partial\rho\right)_{s}}$ is the local adiabatic speed of sound. When the energy flux is fixed, the convection velocity scales with $\rho_{0}^{-1/3}$ in the mixing length theory. This setting approximately maintains the Mach number. Fig. 1 shows the RSST factor $\xi$ (panel a) and the reduced speed of sound $c_{\mathrm{s}}/\xi$ (panel b). Hotta et al. (2012) show that if the Mach number estimated with the reduced speed of sound is smaller than 0.7, the RSST cause no side effect on the result.

We limit the Alfven velocity to 40 $\mathrm{km\ s^{-1}}$ to deal with low-$\beta$ region above the photosphere Rempel et al. (2009b). This does not affect what is occurring at the photosphere.

The radiative heating $Q$ is calculated with the radiative transfer equation; the details are shown in Appendix C.

We calculate the thermal convection with the top boundary at $z=-7\ \mathrm{Mm}$ for 90 days. The thermal convection around the photosphere is very fast, and we exclude this layer for accelerating the calculation. Then, we include the photosphere with the top boundary at 700 km above the photosphere and continue the calculation for 5 days. This procedure is justified because the existence of the photosphere does not change the deep structure (Hotta et al., 2019).

The calculation domain extends 98.304 Mm horizontally for $x$ and $y$ directions. The vertical calculation extent is from the base of the convection zone $z=-0.29R_{\odot}$ and to $z=700\ \mathrm{km}$. The number of grid points in each horizontal direction is 1024, and the grid spacing is 96 km, and this is acceptable to resolve the photosphere. In the vertical direction, we use 512 grid points and nonuniform grid spacing, and this is 48 km around the photosphere and 900 km around the base of the convection zone.

We adopt the Lundquist solution (Lundquist, 1950) for the linear force-free flux tube as follows:

where $a_{0}=2.404825$ is the first root of $J_{0}$, and $r_{\mathrm{tb}}=9.9\ \mathrm{Mm}$ is the radius of the flux tube. $J_{0}$ and $J_{1}$ are the Bessel functions. Here, we adopt $z_{\mathrm{tb}}=-0.05R_{\odot}$, i.e., the initial magnetic flux tube is located about 35 Mm below the photosphere and $B_{\mathrm{tb}}=10^{4}\ \mathrm{G}$ for the magnetic field strength at the center of the flux tube. This leads to a total flux of $1\times 10^{22}\ \mathrm{Mx}$. $y_{\mathrm{tb}}$ is an arbitrary parameter with which the horizontal location of the flux tube is determined. With our choice, the initial flux tube is located in the center of the computational domain. Fig. 2 shows the initial flux tube condition at $z=-0.05R_{\odot}$. The initial flux tube is caught by the two coherent downflows. Note that the left downflow is larger and more coherent. Because the initial magnetic flux tube is force-free, we do not change any thermodynamic variable, i.e., the density, pressure, or entropy from the prepared hydrodynamic calculation. Thus, the magnetic flux initially does not have any buoyancy by the magnetic field, but the inertia of the convective flow can distort the magnetic flux tube.

Fig. 3 shows the temporal evolution of unsigned magnetic flux on $\tau=1$ surface (panel a) and the area of the sunspot (panel b). In Fig. 3a, the black line shows the unsigned magnetic flux integrated over all horizontal domains. The red and blue lines show the unsigned magnetic flux with the emergent intensity less than 80% and 50% of the averaged quiet sun intensity $I_{\odot}$, respectively. Fig. 3b shows the area with the emergent intensity less than $0.8I_{\odot}$ (red) and $0.5I_{\odot}$ (blue). For making Fig. 3b, we use the Gaussian filter with a width of 1 Mm. The unsigned magnetic flux reaches its maximum of $2.7\times 10^{22}\ \mathrm{Mx}$ at $t=63.7$ hr. We follow the definition of the emergence time suggested by Leka et al. (2013), in which the unsigned magnetic flux reaches 10% of the maximum magnetic flux. The emergence time in this calculation is $t=37$ hr. After the magnetic flux reaches its maximum, the magnetic flux begins to decrease. Around $t=175\ \mathrm{hr}$, almost all the magnetic flux in the sunspot $I<0.5I_{\odot}$ disappears, leading to approximate emergence and decay rates ($d\Phi/dt$) of $9\times 10^{20}\ \mathrm{Mx\ hr^{-1}}$ and $-2\times 10^{20}\ \mathrm{Mx\ hr^{-1}}$, respectively. These values are similar to those in Rempel & Cheung (2014). We note that they determine the rising speed of the flux tube at the bottom boundary, while the combination of the convection and the magnetic field determines the rising property in this calculation.

Compared with the observations, Otsuji et al. (2011) and Norton et al. (2017) show the relation of the maximum sunspot flux $\Phi_{\mathrm{max}}$ and the emergence rate $\left(d\Phi/dt\right)_{\mathrm{e}}$ as $\left(d\Phi/dt\right)_{\mathrm{e}}=9.6\times 10^{7}\ \Phi_{\mathrm{max}}^{0.57}$ and $7.8\times 10^{11}\ \Phi_{\mathrm{max}}^{0.36}$, where both are in the unit of $\mathrm{Mx\ hr^{-1}}$. The maximum flux of the current simulation $\Phi_{\mathrm{max}}=2.7\times 10^{22}\ \mathrm{Mx}$ leads to $\left(d\Phi/dt\right)_{\mathrm{e}}=6\times 10^{20}\ \mathrm{Mx}$ and $9\times 10^{19}$ with using Otsuji et al. (2011) and Norton et al. (2017) relations, respectively. In addition, Namekata et al. (2019) show that the 95% confidence interval of the emergence rate reaches $\left(d\Phi/dt\right)_{\mathrm{e}}=7\times 10^{20}\ \mathrm{Mx\ hr^{-1}}$ with $\Phi_{\mathrm{max}}=2.7\times 10^{22}\ \mathrm{Mx}$. Compared with these observational results, this numerical simulation shows a slightly larger emergence rate.

Regarding the decay rate of the sunspot, Hathaway & Choudhary (2008) show the observational relation between the sunspot area $A_{\mathrm{spot}}$ and the dissipation rate $dA_{\mathrm{spot}}/dt$ as $dA_{\mathrm{spot}}/dt=-3\times 10^{16}-5\times 10^{-3}A_{\mathrm{spot}}\ [\mathrm{cm^{2}\ hr^{-1}}]$. We assume that the mean magnetic field strength is 2000 G, which is used in Namekata et al. (2019). This relation leads to $d\Phi/dt=-1.9\times 10^{20}\ \mathrm{Mx\ hr^{-1}}$, which is consistent with the simulation in this study. Fig. 3b shows that the area of the sunspot reaches its maximum $A_{\mathrm{spot}}=5.9\times 10^{18}\ \mathrm{cm^{2}}$ at $t=65$ hr and that the area decreases to zero at $t=154$ hr. This leads to the decay rate of $dA_{\mathrm{spot}}/dt=-6.6\times 10^{16}\ \mathrm{cm^{2}\ s^{-1}}$. The Hathaway & Choudhary (2008) relation leads to $dA_{\mathrm{spot}}/dt=-6.0\times 10^{16}\ \mathrm{cm^{2}\ s^{-1}}$ for our case. Also, in terms of area, the calculation result is consistent with the observation.

Figs. 4 and 5 show the emergent intensity and line-of-sight magnetic field $B_{z}$ at $\tau=1$ surface, respectively. Around the emergence time $t=37$ hr, we begin to see small pores in the intensity map (Fig. 4b) and a diffused pattern in the magnetic field map (Fig. 5b). As time progresses, the small pores merge and construct the large-scale structure. At $t=48\ \mathrm{hr}$, we observe elongated granules between the positive and negative spots (Fig. 4c and 5c). When the photospheric magnetic flux reaches its maximum, a coherent sunspot appears (Fig. 4d and 5d). While we see some evidence of the penumbra around the sunspot, the reproduced penumbra is much less prominent than the observation. Rempel (2012) shows that the existence of the penumbra is significantly affected by the top boundary condition. Here, we use the potential magnetic field condition, while Rempel (2012) argues that the horizontally inclined magnetic field at the top boundary is required for a prominent penumbra. We also note that Rempel et al. (2009a) show that the penumbra can be observed between two sunspots because the horizontal magnetic field is expected between them. In this study, the matured sunspots are far apart from each other, and even the penumbra between the sunspot pair cannot be observed.

After the magnetic flux reaches its maximum, the sunspots begin to lose their magnetic flux. At $t=80\ \mathrm{hr}$, the right sunspot loses significant flux, while the left sunspot keeps a coherent shape with some bright features in the umbra (Fig. 5e). At $t=180$ hr, almost all the magnetic flux disappears from the umbra (Fig. 4f and 5f).

Fig. 6 shows the three-dimensional structure of the magnetic field deep in the convection zone. At $t=20\ \mathrm{hr}$, the magnetic field shows an $\Omega$-shape structure with two anchoring downflows and a broad upflow in the center region (Fig. 6b). At $t=40\ \mathrm{hr}$ (Fig. 6c), which is around the emergence time, a significant fraction of the magnetic flux reaches the near-surface layer ($z>-10\ \mathrm{Mm}$). At $t=60\ \mathrm{hr}$ (Fig. 6d), which is around the time of the maximum photospheric magnetic flux, the root of the left sunspot reaches a depth of around 80 Mm, while the root of the right sunspot remains at a depth of around 30 Mm. Below the photosphere, the magnetic field of the sunspot is mostly vertical. At $t=80$ hr (Fig.6e), the sunspots begin to decay; at the same time, the subsurface structure is destroyed. At $t=180$ hr (Fig. 6f), when the sunspot disappears almost completely, most of the flux is transported downward, and coherent features can still be seen in the deep region below a depth of 60 Mm.

The Coriolis force is an important factor for the tilt angle (Wang & Sheeley, 1991) and the asymmetry of the sunspot pair (D’Silva & Choudhuri, 1993). Cheung et al. (2010) show an almost symmetric sunspot pair because of their almost symmetric initial condition and setting. Rempel & Cheung (2014) show an asymmetric sunspot pair resulting from the horizontal flow mimicking the Coriolis force-induced flow. Chen et al. (2017) also show asymmetry resulting from the boundary condition motivated by the dynamo calculation influenced by the Coriolis force. While the Coriolis force must be importanat to generate a statistical trend of the sunspot pair feature (Weber et al., 2011), our calculation without the Coriolis force can also cause some asymmetry in the sunspot pair. The asymmetry in this study reflects the convection flow structure in the deep convection zone. The initial condition (Fig. 2) shows that the left downflow has a more circular shape than does the right downflow, which shows an elongated feature. Interestingly, the emerging sunspot follows this morphology. Fig. 7 shows the flow and magnetic structure in deep layers at $t=60$ hr. Color and line contours show the vertical velocity ($v_{z}$) and magnetic field ($B_{z}$), respectively. At a depth of at least 30 Mm (Fig. 7b), we can still observe the asymmetry of the magnetic feature, i.e., the positive magnetic feature (solid line) shows the circular shape, while the negative feature (dashed line) is elongated. In the deeper layer, the positive feature maintains a circular shape. The negative feature crosses the boundary and reaches the other side indicated with a black arrow in Fig. 7c because of the periodic boundary condition. At a depth of 50 Mm, the negative feature has a circular shape (black arrow in Fig. 7d). We note that the time scale of the convection at 50 Mm depth is about 5 days, and the structure of the magnetic flux tube has not been distorted well at this depth. This result indicates that the shape of the sunspot is influenced by the deep convection structure at least 30 Mm below the photosphere.

In this section, we investigate the rising process of the flux tube mainly around the center region ($x=41$ Mm). Fig. 8 shows the temporal evolution of the rising flux tube. The left and right panels show the vertical velocity $v_{z}$ and Alfven velocity of the axial field $c_{\mathrm{A}}=B_{x}/\sqrt{4\pi\rho}$, respectively. To investigate the motion of the flux tube, we adopt the following definition of the center of the flux tube. Basically, we use a clumping method. We detect clumps of the region having an Alfven velocity beyond a given threshold. With iteration, we survey the threshold $c_{\mathrm{A0}}$ with which the largest clump has a magnetic flux of $6\times 10^{21}\ \mathrm{Mx}$, where the magnetic flux is estimated with $B_{x}$. The largest clump ($S_{\mathrm{c}}$) in each time step is defined as the rising flux tube. The center of the flux tube $(y_{\mathrm{c}},z_{\mathrm{c}})$ is defined as follows:

The boundary and center of the flux tube are shown with black lines and black crosses in Fig. 8, respectively. The time evolution of the center of the flux tube is shown in Fig. 9. Fig. 9a shows the vertical position of the center of flux tube $z_{\mathrm{c}}$. In the beginning, the rising speed increases with time, when the flux tube reaches the solar surface, the speed is reduced. The time profile of the vertical position is not smooth because of sudden merging and splitting. To evaluate the rising speed, we use a Savitzky–Golay filter for the vertical position, which is shown as the red line in Fig. 9a. By using the filtered profile of the vertical position, we estimate the rising speed of the center of the flux tube. The black line in Fig. 9b shows the rising speed of the flux tube. The red line shows the root-mean-square (RMS) upflow convection velocity in a hydrodynamic run, i.e., without the magnetic field. We note that the downflow convection velocity is typically larger than the upflow velocity. The blue line shows the mean Alfven velocity in the flux tube. In the beginning, the rising speed is almost the same as the upward convection velocity. This is natural because the initial flux tube is force-free and the flux tube needs to obey the convective flow. As time progresses, the rising speed exceeds the convection velocity, indicating some contribution of the magnetic field to the rising process. Even at the end of the rising process, the rising speed does not reach the local Alfven velocity. This result indicates that the rising process is neither a pure magnetic nor a pure convective process.

To understand the rising mechanism, we investigate the vertical equation of motion. Fig. 10a shows the vertical forces averaged in the flux tube. Before evaluating the equation of motion, we filtered out the sound waves in the phase space (Georgobiani et al., 2007). The equation of motion in the vertical direction is the following:

The red and blue lines in Fig. 10a show the buoyancy ($-dp/dz-\rho g$) and Lorentz force ($[(\nabla\times{\mbox{\boldmath$B$}}\times{\mbox{\boldmath$B$}})]_{z}/4\pi$), respectively. The moving flux tube can be roughly regarded as a Lagrangean parcel in the $y-z$ plane. Only the contribution from the inertia term ($-\rho\left({\mbox{\boldmath$v$}}\cdot\nabla\right)v_{z}$) to the motion of the rising flux tube is ($-\rho v_{x}dv_{z}/dx$). The black line is the sum of these terms. The buoyancy and total force direct upward during the rising process, while the Lorentz force directs downward. The main contribution of the Lorentz force is magnetic tension. The roots of the magnetic flux are anchored in the deep layer; thus, the magnetic tension in the rising part directs downward. Fig. 10b shows the mean normalized density, pressure, and entropy in the flux tube. The normalization procedure, i.e., the tilde, for quantities are defined as follows:

where $\langle q\rangle_{\mathrm{hd}}$ and $q_{\mathrm{rms(hd)}}$ are the horizontal average and RMS value of a quantity $q$ in the calculation without the magnetic field, respectively. After normalization, the density, pressure, and entropy are related as $\tilde{\rho}=\tilde{s}+\tilde{p}$. We also note that because $(\partial\rho/\partial s)_{p}$ is negative, the negative value of $\tilde{s}$ corresponds to a positive perturbation of the entropy $s_{1}-\langle s_{1}\rangle_{\mathrm{hd}}$ and vice versa. The result shows that the mean normalized density in the flux tube is significantly low (the perturbation is 1.7 times larger than the RMS value); this is the primary driver of the rising flux tube. There are two contributions to the low density in the flux tube. The first is the low pressure caused by the magnetic pressure. After the rising starts, the twist becomes relaxed and the force-free balance is destroyed. Then, the gas pressure gradient force plays a role in maintaining the flux tube and decreasing the gas pressure in the flux tube. Because the initial flux tube is force-free, the magnetic tension plays a role in balancing the magnetic pressure, even during the rising process; thus, the gas pressure is not significantly low in this case. The other contribution to the low density is the high entropy maintained by the suppression of the mixing. The initial magnetic flux tube is located in the upflow region and has higher entropy than does the surrounding plasma. As seen in Fig. 8c and d, the magnetic field suppresses the mixing between upflows and downflows. In an ordinary medium in the convection zone, the upflow warm medium is mixed in a mixing length and loses its high entropy, but this magnetic flux tube can avoid this process and maintain high entropy, leading to continuously low density in the flux tube and further acceleration. Fig. 10b shows that the high entropy is the main contribution to the low density in the flux tube.

In this section, we investigate the divergent flow prior to emergence time to compare the results with Birch et al. (2016), in which clear divergent flow is not seen at 3 hr before the flux emergence.

To this end, we follow a similar procedure to that of Birch et al. (2016). We take the divergence of the horizontal flow at $\tau=1$ surface and the Gaussian filter with a 6-Mm window. Fig. 11 shows the horizontal divergence $\nabla_{\mathrm{h}}\cdot\boldsymbol{v}$ at $\tau=1$ surface at $t=34$ and 37, and 42 hr. The standard error of the distribution is $3\times 10^{-5}\ \mathrm{s^{-1}}$, and the black line in Fig. 11 shows the 3$\sigma$ value. We again note that $t=37$ hr is the emergence time. Panels a and c are 5 hr after and 3 hr before the emergence time. While we see a clear divergent flow at and after the emergence time (panels b and c), we do not see it 3 hr before the emergence time (panel a). This result is consistent with the observation of Birch et al. (2016). By contrast, the rising speed of the flux tube exceeds $250\ \mathrm{m\ s^{-1}}$ at a depth of 18 Mm, while Birch et al. (2016) suggest that the rising speed is no larger than $150\ \mathrm{m\ s^{-1}}$ to avoid divergent flow prior to the emergence time.

Fig. 12 shows the temporal evolution of the horizontal divergence at $\tau=1$ surface. The black and red lines show the maximum and top 1% value of the divergent flow, respectively. This also indicates that the clear evidence of the divergent flow is seen at only 1 hr before the emergence time.

In this section, we investigate the flow in the flux tube during the flux emergence process. First, we discuss the vertical flow at $z=-30$ Mm in detail. Then, we discuss the flow at different heights. Fig. 13 shows the temporal evolution of the vertical flow and the density in the flux tube. From this section, we focus on the left coherent positive sunspot. Because the root of the sunspot is created with the downflow from the horizontal magnetic flux tube, the root is initially filled with the downflow (Fig. 13a: $t=22.67$ hr). This is somewhat consistent with Chen et al. (2017), in which the strong magnetic concentration occurs at a coherent downflow in the deep layer. A difference is seen in the following evolution. The flow in the flux tube changes signs around $t=30$ hr in Fig. 13c. After the sunspot is formed (Fig. 13e: $t=40$ hr), the flux tube is filled with the upflow. During this process, the flux tube is filled with the low-density medium (Figs. 13b, d, and f).

Here, we define the flux tube as the largest clump in which the vertical magnetic field exceeds the threshold value of 7000 G. Fig. 14a shows the temporal evolution of the mean vertical flow in the flux tube. In the beginning, the flux tube is filled with the downflow. Around $t=23$ hr, the downflow stops increasing in amplitude. Around $t=40$ hr, the vertical flow changes signs. Afterward, the amplitude of the upflow increases continuously. Fig. 14b shows the force balance in the flux tube. The red and blue lines show the buoyancy and the Lorentz force, respectively. The flux tube roughly acts as the Lagrangean parcel in the $x-y$ plane. Thus, for the total dynamics, the contribution from the inertia term is the vertical inertia term $-\rho v_{z}\partial v_{z}/\partial z$. The black line in Fig. 14b shows the sum of the buoyancy, Lorentz force, and vertical inertia terms. The total force is consistent with the temporal evolution of the vertical velocity in the flux tube. Around $t=23$ hr, the total force changes its sign and keeps the positive value after that. The main driver of the upflow in the flux tube is the buoyancy. Fig. 14c shows the normalized density, pressure, and entropy averaged in the flux tube (see eqs. (21)–(23)). In this case, the main contribution of the low density is the low gas pressure caused by the adiabatic expansion from the Lorentz force.

Fig. 15 shows the azimuthally averaged values in the polar coordinate. Because the generated sunspot and its flux tube are tangling along the vertical direction, different heights cannot share the same origin as the polar coordinate. Thus, we define the origin of the polar coordinate at every depth separately. We adopt the following procedure.

1. 1.

   We estimate the magnetic flux of the positive vertical magnetic field in the left half of the computational domain ($\Phi_{+}$).

2. 2.

   The threshold value of the magnetic flux ($\Phi_{\mathrm{s}}$) is $25\%$ of $\Phi_{+}$.

3. 3.

   We search the critical value $B_{\mathrm{c}}$ for the vertical magnetic field $B_{z}$ with which the largest clump has a vertical magnetic flux of $\Phi_{\mathrm{s}}$, where the definition of the clump is the continuous region in which the vertical magnetic field exceeds the critical value. The flux tube is identified with $S_{\mathrm{o}}$.

4. 4.

   The origin of the polar coordinate $(x_{\mathrm{o}},y_{\mathrm{o}})$ at depth is defined with the following:

   	$\displaystyle x_{\mathrm{o}}$	$\displaystyle=$	$\displaystyle\frac{\int_{S_{\mathrm{o}}}xB_{z}dS}{\int_{S_{\mathrm{o}}}B_{z}dS},$		(24)
   	$\displaystyle y_{\mathrm{o}}$	$\displaystyle=$	$\displaystyle\frac{\int_{S_{\mathrm{o}}}yB_{z}dS}{\int_{S_{\mathrm{o}}}B_{z}dS}.$		(25)

5. 5.

   If $\Phi_{+}$ is smaller than $2.5\times 10^{21}\ \mathrm{Mx}$, we use the origin of the polar coordinate defined at a place one grid below.

6. 6.

   The vertical profile of the origin of the polar coordinate is not smooth. We obtain the smooth profile of the origin of the polar coordinate with the Savitzky– Golay filter. The filtered profile is used for the origin of the polar coordinate in the following analyses.

Fig. 15 shows the azimuthally averaged profile of the vertical velocity (Fig. 15a) and the normalized density $\tilde{\rho}$ (Fig. 15b) at $t=66$ hr, where the parenthesis $\langle\rangle$ shows the azimuthal average at constant radius $r$. The finding at a depth of 30 Mm is applicable to different depths. Except for the near-surface layer ($z>-10$ Mm), the center of the flux tube is filled with the upflow, surrounded by the coherent downflow in the outer side of the flux tube. The coherent upflow corresponds to low density (Fig. 15b), indicating that the driving mechanism of the upflow is the buoyancy.

In this section, we discuss the formation mechanism of sunspots. The most important process for the formation of sunspots is the collection of the magnetic flux toward the center region. As discussed in Cheung et al. (2010), the temporal evolution of the magnetic flux is expressed as follows:

where area $S$ is a circular area with fixed radius $r$ at height $z$, and $\langle\rangle$ shows the azimuthal average. Cheung et al. (2010) suggests that the essential term for collecting the magnetic flux at the photosphere is $-\langle v^{\prime}_{r}B^{\prime}_{z}\rangle$. This is caused by the mass loss at the highly inclined field (see Fig. 8 in Cheung et al. (2010)).

Fig. 16 shows the azimuthally averaged radial velocity ($\langle v_{r}\rangle$: panel a) and one of the turbulent induction terms ($\langle v^{\prime}_{z}B^{\prime}_{r}\rangle$: panel b). The definition of the origin of the polar coordinate is the same as that in Section 3.5. The values are averaged between $t=33$ and 37 hr, and this range is around the emergence time, which is before the formation of the sunspot. At this time, the sunspot location in the near-surface is filled with the outflow from the center to the outside, which is the negative effect for the generation of the sunspot. Among the induction terms, the most important contribution is from $\langle v^{\prime}_{z}B^{\prime}_{r}\rangle$. In Fig. 16b, this term is positive in the magnetic flux. This can be explained with Fig. 17. The horizontal magnetic field with positive $B_{x}$ can be regarded as a negative radial magnetic field ($B^{\prime}_{r}<0$) in the downflow side ($v^{\prime}_{z}<0$). In the upflow region ($v^{\prime}_{z}>0$), the radial magnetic field is positive ($B^{\prime}_{r}>0$). Both cause a positive correlation, meaning that before the formation of the sunspot, the shear of the downflow and the upflow are the most important mechanisms for the generation of the vertical magnetic flux.

Fig. 18 shows the azimuthally averaged flows and induction terms averaged between $t=47$ and 53 hr, i.e., during the formation of sunspots. We find a strong coherent downflow and inflow at the near-surface layer ($z>-15$ Mm) in the center (Figs. 18a and b, indicated with a dash-dot line). This makes an important contribution to the collection of the vertical magnetic flux (Fig. 18c). This flow is transient and continues less than 10 hr. While the existence of this downflow and inflow is not reported in Cheung et al. (2010), Rempel & Cheung (2014) show its importance (see their Fig. 11). The treatment of the bottom boundary would be important for the flow for calculations without the deep layer.

In the outer side of the sunspot, we see the Evershed-like outflow in the near-surface layer ($-10$ Mm, indicated with a dashed line in Fig. 18b), which does not contribute to the loss of the vertical magnetic flux in this phase, because the Evershed flow and the magnetic field line are roughly aligned (Fig. 18c). In the deeper layer, the sunspot core is filled with weak outflow, reducing the vertical magnetic flux in this phase (indicated with a doted line in Fig. 18b). In the outer area, we also see weak inflow that does not play any role in the collection of the vertical flux (indicated with a dash-dot-dot line in Fig. 18b).

The most important contribution of the formation of the sunspot in the outer region is a turbulent induction term ($-\langle v^{\prime}_{r}B^{\prime}_{z}\rangle$), which is consistent with Cheung et al. (2010) and Rempel & Cheung (2014). The turbulent induction term ($-\langle v^{\prime}_{r}B^{\prime}_{z}\rangle$) has a negative value in the deeper layer ($z<-15$ Mm). Even the sum of all the induction terms has a negative value in the deep layer during the formation of the sunspot at the solar surface. Consequently, the decay of the sunspot starts during its formation.

In this section, we report the overall flow structure of the matured sunspot in comparison with the sunspot during formation.

Fig. 19 shows the azimuthally averaged flows and induction terms in the same manner as Fig. 18, but the values are averaged between $t=63$ to 69 hr, in which the unsigned magnetic flux has its maximum value. Compared with the sunspot during formation, the strong coherent downflow and inflow in the near-surface layer disappear (Figs. 19a and b), whereas, weak downflow and inflow are seen in the near-surface layer ($z>-5\ \mathrm{Mm}$). Nonetheless, the weak inflow still contribute to collecting the vertical magnetic flux in the near-surface layer (Fig. 19c) . In the deep layer ($z=-5$ to $-20$ Mm), the root of the sunspot is filled with the outflow, which is connected to the Evershed-like flow in the outer area in the near-surface. This outflow in the middle layer makes an important contribution to the decay of the sunspot. In the deeper layer ($z<-20$ Mm), the center of the flux tube is filled with the upflow, as discussed in Section 3.5. This upflow is also connected to the outflow in the middle layer and the Evershed-like flow in the near-surface. We see coherent inflow in the outside the flux tube (Fig. 19b), which is connected to the downflow in the outer layer of the flux tube in the deeper layer. This coherent inflow makes a significant contribution to the collection of the magnetic flux, i.e., the maintenance of the sunspot. The turbulent induction term $-\langle v^{\prime}_{r}B^{\prime}_{z}\rangle$ still makes a significant contribution to the maintenance of the sunspot in the near-surface layer, even if the sunspot is matured.

Fig. 20 shows the azimuthally averaged normalized temperature $\tilde{T}$ averaged between $t=47$ and 53 hr (panel a, during the formation) and $t=63$ and 69 hr (panel b, matured sunspot), where $\tilde{T}=(T-\langle T\rangle_{\mathrm{hd}})/T_{\mathrm{rms(hd)}}$. As the formation process of the sunspot proceeds, the low-temperature region expands because of the suppression of the convective energy transport and the radiative cooling at the photosphere. In both cases, we see a high-temperature region in the center of the flux tube below a depth of $-15$ Mm. While the phase of the sunspot formation is different in panels a and b, the depth where the temperature changes from high to low does not change. In the matured sunspot, the coherent low-temperature region continues to a depth of 40 Mm.

As discussed in Hotta et al. (2012), we need to keep the Mach number smaller than 0.7 in order to avoid the influence of the RSST. In this study, we keep the Mach number smaller than 0.1 when the RSST is used. Fig. 22 shows the RMS velocity in the hydrodynamic calculation without the magnetic field (panel a) and the Mach number estimated with the reduced speed of sound (panel b). The Mach number is always smaller than 0.1 at locations where the RSST is used (indicated with the dashed line in Fig 22b). Phenomena related with the magnetic field, i.e., the flux emergence and sunspot formation, occurs in a similar time scale as the thermal convection. For example, Fig. 13 shows that the up- and downflow are typically 200–300 $\mathrm{m\ s^{-1}}$, where the reduced speed of sound at z=-30 Mm is about 4 $\mathrm{km\ s^{-1}}$. Thus the Mach number is smaller than 0.1 and we do not expect any significant influence of the RSST on the result shown in this paper.

In this study, we prepare an unprecedentedly deep computational domain to minimize the bottom boundary influence. This also has a possible drawback on the convection structure. It is known that numerical simulations tend to overestimate the large-scale (> 30 Mm) convection energy in the deep convection zone (Hanasoge et al., 2012; Lord et al., 2014). To reduce the power on a larger scale than the supergranulation, we intentionally prepare the small box in the horizontal direction of 100 Mm (see also Toriumi & Hotta, 2019). This setting causes another side effect. Because of the insufficient horizontal box size and the periodic boundary condition, a generated sunspot is close to the sunspot on the other side. We expect that the calculation results in this study are not an evolution of the single isolated sunspot pair, but rather, the sequence of the sunspot pairs. This fact mainly influences the decay of the sunspot, and we do not discuss this effect in detail in this study.

In this study, the emergence rate of the magnetic flux is $9\times 10^{20}\ \mathrm{Mx\ hr^{-1}}$, which is slightly larger than the 95% confidence interval ($7\times 10^{20}\ \mathrm{Mx\ hr^{-1}}$, see Section 3.2) suggested by Namekata et al. (2019). There are several possible reasons for this discrepancy. The emergence rate would be determined by the rising speed and the structure of the magnetic flux tube. A fast-rising and intense structure of the flux tube tends to cause a high emergence rate. In this study, we find that the rising speed of the flux tube is not purely determined by the thermal convection; the magnetic effect is also important. Thus, there would be an appropriate initial setting of the magnetic flux tube for reproducing the emergence rate consistent with the observation. Parameter surveys for the initial setting of the flux tube should therefore be carried out in a future study.

By contrast, the number of observations is insufficient to conclude that our emergence rate is larger than reality. There are only several estimations of the emergence rate for the magnetic flux larger than $10^{22}\ \mathrm{Mx}$ (Toriumi et al., 2014). For example, Sun & Norton (2017) show $1.12\times 10^{21}\ \mathrm{Mx\ hr^{-1}}$ for an instantaneous emergence rate. More observations are needed to identify the most probable flux tube structure in the deep convection zone.

Previously, flux emergence simulations in the deep convection zone with a thin-flux tube (D’Silva & Choudhuri, 1993; Moreno-Insertis et al., 1995), anelastic approximation (Fan, 2008; Weber et al., 2011), RSST (Hotta & Yokoyama, 2012), and simulations around the solar surface (Cheung et al., 2010; Rempel & Cheung, 2014; Chen et al., 2017) were almost completely separated. Our approach using the R2D2 code sheds light on connecting these two regions. This is an integral part of both the solar dynamo and the formation of the sunspot. Even in the simulation in this study with the computational domain of the whole convection zone, we still assume the initial magnetic flux tube condition because we have a very low vertical resolution in the deep convection zone ($\sim 900\ \mathrm{km}$). To maintain the magnetic flux of the magnetic flux tube during the rising process from the base of the convection zone to the surface, we need to prepare a high resolution even in the deep convection zone (Cheung et al., 2006). By using the current numerical resources, we are unable to achieve this purpose; however, we will soon be able to perform these types of calculations using next-generation supercomputers, such as Fugaku in Japan. In future studies, we plan to connect dynamo calculation and sunspot formation simulation in a calculation. This approach is expected to reveal the formation process of sunspots in a more self-consistent manner.