{'fraction_non_alphanumeric': 0.036701860231272, 'fraction_numerical': 0.023127199597787834, 'mean_word_length': 4.587359550561798, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Spontaneous Generation of delta-sunspots in Convective Magnetohydrodynamic Simulation of Magnetic Flux Emergence

Observations reveal that strong solar flares and coronal mass ejections tend to occur in complex active regions characterized by delta-sunspots, spot rotation, sheared polarity inversion lines (PILs), and magnetic flux ropes. Here we report on the first modeling of spontaneous delta-spot generation as a result of flux emergence from the turbulent convection zone. Utilizing state-of-the-art radiative magnetohydrodynamics code R2D2, we simulate the emergence of a force-free flux tube in the convection zone that stretches down to -140 Mm. Elevated by large-scale convective upflows, the tube appears on the photosphere as two emerging bipoles. The opposite polarities collide against each other due to the subsurface connectivity, and they develop into a pair of closely-packed delta-spots. The Lorentz force drives the spot rotation and a strong counter-streaming flow of 10 km/s at the PIL in delta-spots, which, in tandem with local convection, strengthens the horizontal field to 4 kG and builds up a highly-sheared PIL. In the atmosphere above the PIL, a flux rope structure is created. All these processes follow the multi-buoyant segment theory of the delta-spot formation, and they occur as a natural consequence of interaction between magnetic flux and turbulent convection, suggesting that the generation of delta-spots and the resultant flare eruptions may be a stochastically determined process.


INTRODUCTION
Solar flares and coronal mass ejections are the plasma process through which magnetic energy is rapidly converted to heat, kinetic energy, and accelerated highenergy particles (Shibata & Magara 2011). It is known that major flares emanate from complex active regions (Toriumi & Wang 2019). For instance, statistical studies revealed that the δ-sunspots, in which umbrae of opposite polarities are in close vicinity of each other that they share a common penumbra, are prone to strongest flares (Künzel 1960;Sammis et al. 2000;. Other key features include: spot rotation (Brown et al. 2003;Yan et al. 2008); shear flows along the polarity inversion line (PIL) (Harvey & Harvey 1976;Krall et al. 1982); PIL with strong magnetic field, B z gradient, and shear (Hagyard et al. 1984;Schrijver 2007); magnetic channel (Zirin & Wang 1993;Wang et al. 2008); and magnetic flux rope (Gibson et al. 2006).
A major difficulty faced while understanding the formation of flare-productive active regions through emergence from the convection zone is that we cannot investigate the subsurface magnetic field from direct optical observations. Therefore, a number of magnetohydrodynamic (MHD) flux emergence simulations have been conducted in the past decades.
One of the suggested scenarios for δ-spot formation is the multi-segment buoyant model, where a subsurface magnetic flux rises at two locations and appears on the surface as a pair of emerging bipoles. Within this quadrupolar region, the inner polarities of opposite signs collide against each other to form a δ-structure with a sheared PIL in between. This scenario was advocated by Toriumi et al. (2014), who modeled the emergence of a single horizontal tube that is initially made buoyant at two segments. They found that the confinement of the opposite polarities on the surface occurs because the two emerging sections are connected by a dipped field beneath the surface. This situation was followed by Fang & Fan (2015), while Oi (2017) 1 and Syntelis et al. (2019) modeled flare eruptions from the quadrupolar system. Other scenarios include the emergence of a kink-unstable flux tube (Fan et al. 1998;Takasao et al. 2015;Knizhnik et al. 2018) and the collision of two emerging tubes (Murray & Hood 2007;Jouve et al. 2018;Cheung et al. 2019).
Although these simulations succeeded in reproducing some key aspects of flaring regions, many were performed in highly idealized or controlled situations. For instance, the convection zone is mimicked as a planeparallel atmosphere without including convective flows; the tube's emergence is triggered by artificially reducing the density from the tube; or the emerging flux is kinematically advected into the domain through the boundary. However, in reality, the spot formation occurs as a natural consequence of the interaction between magnetic flux and background convection, which has been difficult for modelers to accommodate.
In this Letter, we report on the first spontaneous δspot formation that follows the multi-buoyant segment scenario as a result of flux emergence from the turbulent convection zone. Utilizing the newly developed radiative MHD code, which solves thermal convection of various scales from 100 Mm-sized cells to surface granules, we are now able to overcome the above issues and address the effect of turbulence on emerging flux, the δ-spot formation, and the magnetic properties.

NUMERICAL SETUP
The numerical simulation was performed with the radiative MHD code R2D2, which stands for the Radiation and RSST (reduced speed of sound technique) for Deep Dynamics (see Hotta et al. 2019, for the details). In brief, this code solves the MHD equations with taking into account the radiative energy transfer and adopts RSST (Rempel 2005;Hotta et al. 2012Hotta et al. , 2015Iijima et al. 2019) to deal with the fast sound speed and mitigate the Courant-Friedrichs-Lewy condition. We adopt the gray approximation for the radiative transfer with the Rosseland mean opacity.
1 http://www-space.eps.s.u-tokyo.ac.jp/group/yokoyama-lab/thesis/oi2017 mastereruption.pdf The Cartesian box spans 98.3 Mm (x)×98.3 Mm (y)× 139.9 Mm (z), resolved by a 1024 × 1024 × 256 grid. The horizontal (x and y) grid spacing is 96 km (uniform), while the vertical (z) spacing increases from 48 km around the top boundary to 2950 km around the bottom boundary. The top boundary is located 700 km above the average τ = 1 surface. The bottom boundary is thus located at −139.2 Mm (i.e. 0.8R ⊙ ), which is deep enough compared to the thickness of the convection zone (∼ 200 Mm). Therefore, we can self-consistently solve large-scale convection and surface granulation at the same time.
We assume periodic boundaries for both horizontal directions. The top boundary is open for upflows and closed for downflows, while the density and entropy perturbation from the initial state are free there. The bottom is open for the flows. For mass conservation, the horizontally averaged density is fixed to the initial condition, and the perturbation from the average is free. At the bottom, the entropy in upflows is fixed to the initial value and is free in downflows. We adopt the stress-free boundary condition for the horizontal velocity at both the top and bottom boundaries. The magnetic field at the top is matched to a potential field above, whereas all three components of the field are symmetric about the bottom.
We calculated the convection without magnetic field, first for 60 solar days over a domain up to −2.35 Mm (0.997R ⊙ ) with an artificial cooling layer and then for five days over the whole domain until a statistically equilibrium state was attained. This procedure is justified because the existence of the surface does not influence the deep convection structure . Then, at t = 0 hr, we introduced a horizontal flux tube at −16.7 Mm, which is given as an x-directed force-free Lundquist field (Lundquist 1951): where r is the radial distance from the axis, B tb = 10 kG the axial field strength, J 0 and J 1 the Bessel functions, α = a 0 /R tb , a 0 = 2.404825, and R tb = 7 Mm the tube's radius. The total axial flux is 6.6 × 10 21 Mx. Unlike many of the previous δ-spot simulations where the tube was artificially made buoyant or kinematically inserted, the tube here was initially in mechanical balance and thus started moving only in response to the background convection flows. It should be noted that because of the Alfvén speed limiting (Rempel et al. 2009), with which we limit the Alfvén speed to a maximum of 40 km s −1 to accelerate the computation, the physical quantities may be affected, especially within the umbral cores in the upper photosphere. In the next section, we focus on the generation and properties of δ-spots. However, readers may consult H. Hotta & H. Iijima (in preparation) for a detailed account of flux emergence and spot formation. Figure 1 shows the temporal evolution of the emerging flux tube. The left and second columns are the emergent intensity normalized by its quiet-Sun average (I/I 0 ) and the vertical field strength at the τ = 1 surface (B z : magnetogram), respectively. The contours indicate the umbra/penumbra and penumbra/quiet-Sun boundaries, which are defined as I/I 0 = 0.45 and 0.9, respectively, measured on the intensity map smoothed using a convolution with a Gaussian kernel with a FWHM of 3 Mm (Rempel 2015). The third column shows the vertical velocity (V z ) at −16.7 Mm, i.e., the initial depth of the tube axis. The column on the right presents the normalized magnetic field strength, defined as (2) where y 1 = 24.6 Mm, y 2 = 73.8 Mm, ∆y = y 2 − y 1 , and ρ 0 (z) is the initial background density. This quantity indicates the field strength averaged in the y-direction and normalized by the square root of the local initial density and thus possesses the dimension of velocity. At t = 0 hr, a strong upflow starts to elevate the flux tube around x = 0 Mm (i.e. x = 98.3 Mm because of the periodic side boundaries) and creates an Ω-shaped loop, which develops into an emerging bipole P1-N1 on the surface (see panels for 12.0 hr). Between P1 and N1, a number of small-scale magnetic elements are scattered and exhibit a net-like structure. Through merging and cancellation, they develop into mature sunspots with penumbrae. These behaviors resemble previous observations and models of flux emergence (e.g., Strous et al. 1996;Pariat et al. 2004;Fan et al. 2003;Cheung et al. 2010).

General Evolution
Meanwhile, in the convection zone, another rising section appears around x = 40 Mm (see the V z map for 12.0 hr) and produces the secondary Ω-loop, which eventually appears as the bipole P2-N2 at the domain center between N1 and P1 (see panels for 32.0 hr). The spot area -measured as the total pixels with I/I 0 < 0.9 in the smoothed intensity map -and the total unsigned flux attain their peak values of 8.0 × 10 8 km 2 , equivalent to 260 MSH (millionths of the solar hemisphere), and 2.5 × 10 22 Mx, respectively, at 40 hr.
Because the original horizontal tube was elevated at two segments, the legs of the two Ω-loops approach each other to form the spot pairs N1-P2 and N2-P1. As time progresses, each pair collides and eventually shares a common penumbra, building up a δ-spot (see panels for 43.3 hr). Within each δ-spot, elongated dark convection cells are trapped between the umbrae of opposite polarities (δ-spot light bridge). In the accompanying movie, the spots show vigorous rotational motions in the counter-clockwise direction and the spots are connected below the photosphere by U-loops that are stretched down to about −40 Mm by strong downdrafts.

Sunspot Rotation
In order to examine the driving mechanism of the observed spot rotations, we focus our attention on N1 in Figure 2. The velocity vectors averaged over 8.3 hr in panel (a) demonstrate that the rotation is not the apparent effect but the actual plasma motion. One may find that the convection cells in the penumbra and light bridge are stretched in the direction of the spot rotation. Panel (b) shows the temporal evolution of the spot structure along the circular slit around N1, which clearly shows an apparent counter-clockwise motion with the angular velocity up to 15 • hr −1 .
The azimuthal component of the temporally averaged plasma velocity along the slit is plotted in panel (c). The velocity ranges from almost 0 km s −1 in the umbral regions to more than 5 km s −1 in the penumbra and the light bridge. Overplotted are the azimuthal components of the Lorentz force ((∇×B)×B/(4π)) and gas pressure gradient (−∇p) along the slit. In most of the regions, the Lorentz force dominates the pressure gradient, indicating that the main driving force of the spot rotation is the Lorentz force.
The above result is akin to those of the previous flux emergence simulations and is explained in the following manner (e.g., Longcope & Welsch 2000;Fan 2009;Sturrock et al. 2015). When an isolated flux tube with a left-handed twist emerges into the low-density atmosphere, it drastically expands and the field lines are therefore bent in a manner that the Lorentz force acts in the direction opposite to the original twist, i.e., counterclockwise. The observed counter-clockwise spot rotations are a consequence of this Lorentz force, and the flux tube injects its helicity into the atmosphere through this process.    two rotating umbrae. The flow runs along the PIL and the horizontal field vectors are aligned with the flow vectors in most parts, constituting a highly-sheared PIL.

Magnetic Properties of the δ-spot PIL
Cross-sectional profiles in panels (c-h) reveal that the τ = 1 surface is elevated in the PIL by a few 100 km. The vertical magnetic field inverts its sign across the PIL with a steep gradient of about 1 kG Mm −1 on average. The horizontal magnetic field is strongly sheared and intensified by the horizontal flows, with the field strength being up to about 4 kG. One may find that the strong B z concentrations in the magnetic channel structure appear at the locations of steep velocity shear.
As shown above, the velocity shear of the counterstreaming flow, driven by the Lorentz force acting on each spot umbra (Section 3.2), intensifies the magnetic shear at the PIL. In the previous ideal δ-spot model by Toriumi & Takasao (2017) that lacks thermal convection, the horizontal field at the PIL was in fact intensified, but only up to an equipartition field strength of about 1 kG. In the present model, on the contrary, the vigorous convection continues around the PIL (δ-spot light bridge), which further strengthens the field to 4 kG.
The magnetic structure above the surface is presented in Figure 4. Above the PIL, the field lines are low-lying and highly sheared (yellow to red), which is covered by the less-sheared field lines that extend from the umbrae on both sides (aqua to blue). The overall field configuration takes the form of a twisted flux rope, which is often observed in the corona in the pre-eruption phase. The sheared field at the flux rope core is due to the advection by the sheared flows, whereas the overlying arcade is more potential because the umbral cores rotate less strongly (see Figures 2(a) and 3(b)).

DISCUSSION
In this Letter, we reported on the spontaneous generation of δ-spots, which are known to be flare-productive, by performing a realistic flux emergence simulation. Thanks to the deep enough computation domain realized by the state-of-the-art R2D2 code, we can assess the effects of radiative transfer and thermal convection on emerging flux. Although some physical quantities may not be directly compared with observations because of the numerical assumptions, the general properties revealed in this study lead to a better understanding of the genesis of flare-productive regions.
The initial force-free flux tube was pushed up by largescale upflows at two separate portions as a pair of emerging bipoles. In the photosphere, umbrae of opposite polarities collided with each other and formed a pair of strongly-packed δ-spots because their legs were con- nected by deeply-anchored U-loops. 2 This result points to the possibility that the multi-buoyant segment scenario occurs on the actual Sun.
One interesting hypothesis we can derive from our results is that the δ-spot formation, and therefore the resultant flare eruptions, might be a probabilistically determined process depending on where a magnetic flux is located on the turbulent background. In fact, in some of the test cases where we only changed the location of the initial flux tube in the same background convection, we found that the tube appears as a single emerging bipole and never produces δ-spots (H. Hotta & H. Iijima, in preparation).
Another point to be discussed is the effect of deep convection structures. Helioseismology points out that convections in numerical simulations are prone to deviation from the observations, especially in large scales (> 30 Mm) (Hanasoge et al. 2012;Lord et al. 2014). However, because in the present model, we selected the small horizontal domain extent of ∼ 100 Mm, which inhibits the generation of unreasonably large cells, the cell sizes of > 30 Mm may not be critically different from the reality. Also, the spot generation may be impacted by the choice of the bottom boundary depth and thus the convective patterns in deep layers. However, again, since the present δ-spots stem from the supergranulescale emerging flux, even if we limit the bottom depth to, say, −30 Mm, the result may not differ much as far as we keep the domain that is deep enough to harbor supergranulation.
Each polarity of the generated δ-spots showed a strong rotation for hours, which was driven by the Lorentz force via unwinding of the tube's twist. The observed angular speed of 15 • hr −1 is somewhat larger than the reported values (e.g., Min & Chae 2009). The rotation speed may depend on the tube's initial twist, but the potential-field top boundary may allow for the rapid helicity release as well (Cheung et al. 2010).
The collision of two rotating spots of opposite polarity generated a counter-streaming flow at the PIL, which greatly enhanced the magnetic shear and horizontal field. The observed field strength of 4 kG was due to the velocity shear as well as the light bridge convection. A magnetic channel, some of which is suggested to be a flare-triggering field (e.g., Kusano et al. 2012), was created by the strong velocity shear, implying that small-scale local convection may play a crucial role in evoking flare eruptions (Toriumi et al. 2013).
Magnetoconvective property of a δ-spot light bridge is in many ways similar to that of a regular light bridge that separates the umbrae of the same polarity. Comparing the cross-sectional profiles in Figure 3 of this Letter and Figure 4 of Toriumi et al. (2015), one may find that in the regular bridge, the horizontal flow is dominant and the elevated iso-τ surfaces show a cusp structure that is sandwiched by the canopy field fanning out from the adjacent umbrae, whereas in the δ-spot bridge, a counter-streaming flow is confined by the low-lying sheared field that connects the neighboring umbrae.
The objective of this Letter is to reveal how the turbulent solar convection allows for the spontaneous formation of δ-spots, sheared PILs, magnetic channels, and twisted flux ropes, all of which are the key ingredients of flare-productive active regions. In the forthcoming papers, we intend to perform detailed analyses on the generation of sunspots, strong PIL horizontal fields, and coronal response to the δ-spot formation.

ACKNOWLEDGMENTS
The authors are grateful to the anonymous referee for improving the manuscript. The results were obtained using the K computer at the RIKEN (proposal nos. hp190070, hp180215, hp180042, and ra000008). This work was supported by JSPS KAK-ENHI Grant Numbers JP15H05814 (PI: K. Ichimoto) and JP18H04436 (PI: H. Hotta), by the NINS program for cross-disciplinary study (Grant Numbers 01321802 and 01311904) on Turbulence, Transport, and Heating Dynamics in Laboratory and Astrophysical Plasmas: "SoLaBo-X", and by MEXT as Exploratory Challenge on Post-K Computer (Elucidation of the Birth of Exoplanets [Second Earth] and the Environmental Variations of Planets in the Solar System).

Facilities: K computer
Software: R2D2