Turbulent convection as a significant hidden provider of magnetic helicity in solar eruptions

Solar flares and coronal mass ejections, the primary space weather disturbances affecting the entire heliosphere and near-Earth environment, mainly emanate from sunspot regions harbouring high degrees of magnetic twist. However, it is not clear how magnetic helicity, the quantity for measuring the magnetic twist, is supplied to the upper solar atmosphere via the emergence of magnetic flux from the turbulent convection zone. Here, we report state-of-the-art numerical simulations of magnetic flux emergence from the deep convection zone. By controlling the twist of emerging flux, we find that with the support of convective upflow, the untwisted emerging flux can reach the solar surface without collapsing, in contrast to previous theoretical predictions, and eventually create sunspots. Because of the turbulent twisting of magnetic flux, the produced sunspots exhibit rotation and inject magnetic helicity into the upper atmosphere, amounting to a substantial fraction of injected helicity in the twisted cases that is sufficient to produce flare eruptions. This result indicates that the turbulent convection is responsible for supplying a non-negligible amount of magnetic helicity and potentially contributes to solar flares.


Introduction
Solar flares and coronal mass ejections are the primary sources of space weather disturbances driving various plasma processes in the entire interplanetary space, including the near-Earth environment [1][2][3] .The strongest events among these eruptions emanate from the solar active regions, in which the rotation and shear motions of strongly magnetised sunspots are often observed 4 .Once a solar flare occurs and a coronal mass ejection is launched, a helical magnetic flux rope is observed in the interplanetary space 5,6 , which is the clearest evidence that solar flares are the sudden release of the excessive magnetic energy accumulated in the helical magnetic structure in the solar corona through magnetic reconnection and plasma instability [7][8][9] .
Magnetic helicity is a measure to quantify the topology of a magnetic field, such as twists, kinks, and internal linkages: where A is the vector potential of a magnetic field B, i.e., B = ∇ × A. It is well conserved even in resistive magnetohydrodynamic (MHD) processes but gauge-invariant only if the magnetic field is fully contained in a closed volume.Otherwise, the relative magnetic helicity, 10,11 is widely used because it is gauge invariant in any case, and the potential magnetic field B p is often adopted as a reference.The total amount of injected helicity, H R , and the helicity injection rate (helicity flux), where the subscripts h and z denote the horizontal and vertical directions, respectively, are widely used in analyses of flareproductive active regions.For instance, observations show that active regions with a larger amount of magnetic helicity tend to produce stronger flares [12][13][14] and that the activity level is enhanced in association with the total injected helicity or the temporal variation of helicity flux [15][16][17][18][19] .This is why magnetic helicity attracts broad attention of the solar and heliophysics community in the context of flare prediction and forecasting.
Based on a copious amount of observational evidence, it is widely believed that the injection of magnetic helicity into the corona is mainly due to the emergence of twisted magnetic flux from the convection zone and the consequential motions of sunspots, such as rotation and shearing 16,[20][21][22][23][24][25][26][27] .Therefore, investigating how magnetic helicity is accumulated in the corona as the magnetic flux emerges and builds up active regions is an important factor in understanding the occurrence of solar flares.However, it has almost never been considered how and to what extent the background convection affects the helicity injection, which is probably because we cannot probe the solar interior using direct optical observations.][30] To solve this problem, we perform numerical simulations in which a twist-free magnetic flux tube emerges from the turbulent convection zone and examine whether the helicity injection into the upper atmosphere is negligibly small or comparable to the observations.Previous theoretical studies suggested that an untwisted emerging flux cannot reach the photosphere because it experiences counteracting aerodynamic drag (see 31 and the references therein).Under the ideal condition with no background convection, even if the non-twisted flux tube reaches the photosphere, the helicity injection will be zero (as far as the emerging bipole keeps its geometrical symmetry).In this study, we explore the occurrence of helicity injection by calculating the case with convection and, if so, to what extent.
We perform our computations using radiative MHD code R2D2 32 to reproduce the realistic turbulent thermal convection in the Sun.We use the convection model that is in a statistically steady state as the initial condition (t = 0).We do not consider solar rotation in our model; therefore, the net kinetic helicity of the background flow field is negligibly small.The simulations with no solar rotation and thus infinitesimal net kinetic helicity enable investigation of the magnetic helicity injection that is caused purely by turbulence.
At t = 0, a magnetic flux tube with an axial field strength of 12 kG and a typical radius of 8 Mm is placed at a depth of 22 Mm in the rectangular computational box that covers the entire solar convection zone, with the box spanning over 0 ≤ x ≤ 98.3 Mm, 0 ≤ y ≤ 98.3 Mm, and −201 Mm ≤ z ≤ 676 km (bottom panel of Figure 1).The bottom boundary (−201 Mm) is deeper than most previous convective flux emergence simulations, which were at most −30 Mm [33][34][35] .Thus, the present model allows us to investigate the effects of large-scale convection on the emergence of magnetic flux and the resultant sunspot formation [36][37][38][39] .In the initial state, a mechanical balance between the flux tube and surroundings is achieved by lowering the entropy inside the tube.Therefore, the magnetic flux starts rising in response to the background velocity field without artificial buoyancy.For the purpose of comparison, we calculate two additional cases where the flux tubes are weakly and strongly twisted (right-handed twist), but the background convection field remains the same.The twist strengths of the weakly-and strongly-twisted tubes are q cr /4 and q cr /2, respectively, where q cr (= 0.125 Mm −1 ) is the critical twist for the kink instability 40 , namely, these tubes are stable against the instability.In all three cases, the total magnetic flux in the axial direction is on the order of 2 × 10 22 Mx.

General evolution
Figure 1 shows the vertical magnetic field B z in the photospheric surface, the emergent intensity and the magnetic field strength and field lines in the 3D space for the three runs, i.e., the non-, weakly-and strongly-twisted cases.In this study, we define the layer where the optical depth is unity (τ = 1) as the photospheric surface and measure physical quantities to make direct comparisons with observations.Figure 1 shows that the flux tube is levitated by a convective upflow located in the centre of the computational box in all three cases.Even the non-twisted flux tube successfully reaches the photosphere, in contrast to the previous prediction that a flux tube needs some twisting to maintain its cohesion against aerodynamic drag 31 .
From t = 20 hr, the positive and negative magnetic elements distributed in the photosphere gradually assemble and build up sunspots with positive and negative polarities.The positive sunspot moves beyond the side periodic boundary and at around t = 30 hr, the two sunspots collide with each other, eventually creating a strongly-packed bipolar sunspot (δ -sunspot), which is known to be highly flare-active 41 .This situation agrees with the scenario where a single flux system emerges at multiple locations to build up a colliding bipolar sunspot 36,42,43 .The magnetic energy of the flux tube at the initial state is set to be the same for the three cases.However, once the sunspots are established in the photosphere, their decay is more rapid for the weaker twist cases because the local convection cells can more easily intrude into the sunspots and break them into pieces.In the twist-free case, the sunspots disappear by around t = 60 hr, i.e., the lifetime in the photosphere is approximately 40 hr.
One remarkable feature of the developed sunspots is their continued rotation.In the strong twist case, the two sunspots rotate in the same clockwise direction because the flux tube, endowed initially with right-handed torsion, releases its twist as it appears in the photosphere 36,44,45 .Observations show that many more flaring active regions show rotations of bipolar sunspots in the same directions than in the opposite directions 46 .Of particular note, however, is that the sunspots in the no-twist case also exhibit rotations, with the negative sunspot in a clockwise direction and the positive sunspot in a counterclockwise direction (the square-framed zone in Figure 1).Because this flux tube is not given any twist at the initial state, the observed sunspot rotations are presumed to be driven by the background turbulence beneath the solar surface.

Magnetic helicity injection and sunspot rotation
Figure 2 shows the temporal evolutions of injected magnetic helicity, H R , and helicity flux, dH R /dt, as compared with the total unsigned magnetic flux, all measured in the photosphere for the three cases.The initial total magnetic flux in the flux tubes in the axial direction is of the order of 2 × 10 22 Mx.Thus, if a flux tube emerges bodily as a whole to form a pair of sunspots, the total photospheric flux would be 4 × 10 22 Mx.However, in all cases it is on the order of 3 × 10 22 Mx, indicating that some fractions of the flux tubes remain below the photosphere.The total fluxes peaked around t = 40 hr and then gradually decreased as the sunspots decayed.
The middle panel of Figure 2 illustrates that positive magnetic helicity is injected in the cases of strong (black line) and weak (blue line) twists.This is reasonable considering that these flux tubes were initially given right-handed twists and that the sunspots displayed clockwise rotations.What is noteworthy about this plot is that the positive helicity injection also occurs in the non-twisted case (red line).The accumulated helicity amounts as much as to 2.9 × 10 43 Mx 2 , which is about 20% to 50% of those in the twisted cases.This result reveals that a non-negligible amount of magnetic helicity can be injected even when a twist-free flux tube emerges in the convective zone with null net kinetic helicity.It should be noted that the magnetic helicity normalised by the square of the total photospheric flux, H R / max (Φ) 2 , for the three cases is 0.029, 0.058 and 0.094 (peak values), which are comparable to the recent observations of flaring active regions 47,48 .Therefore, the present simulations provide the reasonable reproduction of solar active regions.The bottom panel of Figure 2 shows the evolution of helicity flux (see Equation ( 3)).
Then, what effect causes a finite positive helicity injection in the non-twisted case? Figure 3 shows the 3D rendering below the sunspot collision area (corresponding to the square-framed zone in Figure 1) and the plots of velocity fields at three different depths.We average the data over the period from t = 35 hr to 38 hr, when the helicity flux is peaked (bottom panel of Figure 2).Below the positive and negative sunspots, a pair of vertical magnetic fluxes (represented by yellow magnetic field lines) extend downward into the deep convection zone.The velocity plots reveal that the two magnetic pillars (indicated by contours) reside in the strong downflowing plumes, to which surrounding plasmas stream in, leaving local vortices (highlighted by cyan arrows).These structures are created because some portions of the initially horizontal flux tube are dragged into the strong downflow plumes and become vertical magnetic concentrations 35,37 .At the same time, the plumes also drive the surrounding plasmas to flow into them, accompanied by the vortices.The directions of the vortices agree with those of the sunspot rotations in the photosphere, i.e., clockwise (counterclockwise) for the negative (positive) sunspot.This indicates that the local vortices streaming into the downflow plumes spin the magnetic pillars and drive the sunspot rotations in the photosphere.If the sunspot rotation is because of the release of magnetic twists, the two magnetic pillars (sunspots in the photosphere) should rotate in the same direction 44,45 .However, this is not the case.Thus, even if the flux tube is initially endowed with no magnetic twist, it is possible that the surrounding turbulent flow, in which the local vortices reside, exerts a spinning effect on the flux tube and, accordingly, the sunspot rotation occurs in the photosphere, leading to the significant injection of magnetic helicity into the upper atmosphere.

Flare productivity
One of the most promising methodologies that have been suggested to predict the flare eruptions is to analyse the spatial distribution of high free-energy regions (HiFERs), where the non-potential magnetic field, B np = |B − B p |, exceeds the critical value of 1 kG, and examine the occurrence of the double-arc instability (DAI 50 ) for each HiFER patch near the polarity inversion line 49 .Top and middle panels of Figure 4 show the distributions of B np and the magnetic twist flux density τ twist at the representative timings for the non-twisted and strongly-twisted cases.The bottom panels are the flare phase diagram, which is the scatterplot of the minimum critical radius of the circular reconnection region required to satisfy the DAI condition, r c , versus the magnetic energy that can be released, E r , for the identified HiFERs.Observations show that large X-class flares likely occur when there is a HiFER with r c < 1 Mm and E r > 4 × 10 31 erg on the flare phase diagram 49 .
In the strongly-twisted case, the HiFER patches (B np > 1 kG) are coherent and the positive twist is distributed at the polarity inversion line between the two major sunspots, which is in agreement with the analysis results in the previous subsections.In the flare phase diagram, several locations are very close to the X-flare criticality.Therefore, we can predict that M-class flares are possible with a slight chance of X flares for the strong twist case.In the non-twisted case, however, although there are regions where B np exceeds 2 kG, the HiFERs are more fragmented.The twists around the polarity inversion lines are a mixture of positive and negative signs, suggesting that the twists of both signs can be generated, probably randomly, because of the nature of turbulent convection, when a twist-less magnetic flux emerges.No HiFERs satisfy the X-flare threshold in the flare phase diagram, but lower-class (e.g.C-class) flares may still occur.

Discussion
We investigate the physical mechanisms that provide magnetic helicity into the solar corona, which is critically important in understanding the occurrence of flare eruptions, by comparing the flux emergence simulations with and without magnetic twists under the existence of turbulent thermal convection.It is widely believed that the emergence of twisted flux tubes and associated sunspot rotations supply helicity into the active region corona.However, this is not trivial because of the lack of optical observations below the photosphere, which is the primary motivation of this study.We summarise the key results as follows.
1.A non-twisted flux tube can reach the photosphere and develop an active region.This is contrary to the previous simulation results without thermal convection during the emergence, a rising untwisted flux tube splits into two vortex rolls and quickly loses its identity because of the counteracting aerodynamic drag 31 .With the support of convection, a non-twisted flux tube may complete its journey to the solar surface.
2. The injected magnetic helicity is non-zero even in the case of no magnetic twist.In fact, the amount of injected helicity was quite large, reaching about 50% as in the twisted cases.Our analysis reveals that in the convection zone below the developed sunspots, inflows into the strong downwelling plumes produce local vortex motions, which spin the vertically standing magnetic flux tubes that extend along the plumes.As a result, the sunspots in the photosphere show rotations in the directions determined by the subsurface vortices, and a considerable amount of helicity injection into the upper atmosphere occurs.
3. The instability analysis suggests that the generated non-potential magnetic field in the non-twisted case can produce eruptions, albeit weaker than the twisted cases.This indicates that the "twist" observed in a variety of forms in the flare-productive active regions of the Sun is not only supplied by the magnetic "twist" of the emerging flux and the associated sunspot rotations, as previously believed, but also includes a non-negligible fraction of the "twist" brought by the background turbulent convection.
In the present simulations, although the net kinetic helicity of the initial background convection was infinitesimally small, the resultant helicity injection in the photosphere is comparable to that provided by the emergence of twisted flux tubes.Considering the nature of the turbulence where we do not consider the solar rotation, the directions of the rotations of the local vortices that contributed to the spinning of the flux tubes may probably be determined by chance.Because the downflow plumes accompanied by the vortices are a coherent structure owing to the large time scale in the deeper convection zone, although the directions of the vortices are randomly determined, the sign of helicity injection may remain the same (positive in the present runs) over time.This indicates that the sign of helicity injection in simulations with different background turbulence can be negative, which will be examined in the future studies.
Asymmetry in sunspot rotation is reported in bipolar active regions 46,51 , and this could be due to the difference in vortex motions in the convection zone.Also, a considerable fraction of solar active regions do not obey the hemispheric helicity sign rule, with the fraction of violators being 20-40% 52,53 .This randomness may be the result of the stochastic imposing of magnetic helicity by the background turbulence to the twist of the emerging flux that is determined by the solar dynamo 54 .

Methods
In this study, we use the radiative MHD code R2D2 32 , which realistically simulates the thermal convection over the entire solar convection zone, to model the emergence of magnetic flux tubes and the spontaneous sunspot formation [36][37][38][39] .This code computes the 3D MHD equations with realistic radiation transfer and equation of state, and by implementing RSST [55][56][57][58] , it effectively suppresses the high sound speed in the deep solar interior and relaxes the Courant-Friedrichs-Lewy condition.
We use a rectangular Cartesian box as the computational domain, spanning over 0 ≤ x ≤ 98.3 Mm, 0 ≤ y ≤ 98.3 Mm, and −201 Mm ≤ z ≤ 676 km, resolved by a 1024 × 1024 × 384 grid.The grid spacing for the horizontal directions was ∆x = ∆y = 96 km (uniform), while that for the vertical direction was uniform at ∆z = 48 km from the top boundary to z = −5.6Mm, which linearly increased to the bottom boundary up to ∆z = 1486 km.The horizontal boundaries assumed the periodic boundary condition, while the magnetic field at the top boundary was connected to the potential field.The initial background convection was the same as that in 38 , in which a magnetic flux tube emerges to successfully produce a bipolar sunspot group.The net kinetic helicity is negligibly small because the calculation box did not consider the solar rotation.The normalised net kinetic helicity at t = 0, when the flux tube was embedded, was;

4/13
The magnetic flux tube was embedded at a depth of 22 Mm.Unlike the previous R2D2 simulations, which used the force-free magnetic flux tubes, we adopt commonly used Gaussian-type flux tubes, where the longitudinal and azimuthal components of the magnetic field are given by and where B tb , r, R tb and q are the axial field strength, radial distance from the tube axis, typical radius and twist intensity, respectively.We calculate three cases with varying q (Table 1).Here, B tb was also varied so that the magnetic energies of the flux tubes, E mag (= B 2 /(8π) dV ), are the same in all three cases.The total axial magnetic flux Φ x (= B x dS) is of the order of 2 × 10 22 Mx.To achieve the mechanical balance with no magnetic buoyancy, we adjust the internal pressure and density by reducing the entropy by . Therefore, the flux tubes started emergence only because of advection exerted by the background turbulence.In all three cases, we set q to be 0 or below the critical value for the kink instability q cr (= 1/R tb ) 40 , where q > 0 indicates that the tube has a right-handed twist.1. Summary of the simulation cases.The critical twist for the kink instability corresponds to q cr = 1/R tb = 0.125 Mm −1 .

Case
We measure several physical quantities including the relative magnetic helicity (H R ), helicity flux (dH R /dt) and total unsigned magnetic flux (Φ) at the photospheric surface, here defined as where the optical depth is unity (τ = 1).For measuring the helicity, there is a degree of freedom in choosing the vector potential A p and, in this study, we select vector potentials satisfying A p • ẑ = 0 and calculate it using the method by 59 , where r = x − x ′ .In this process, the grid number was reduced from 1024 × 1024 to 256 × 256 to accelerate the computation speed.The helicity flux was then calculated using Equation (3), and the injected magnetic helicity over the course of flux emergence (Equation ( 2)) was obtained by integrating the helicity flux over time.
To investigate the flare productivity of the generated active regions, the instability analysis is performed based on the DAI theory, which states that a flux rope becomes unstable if it is sufficiently twisted and has a sufficient amount of magnetic flux against the overlying confinement field 50 .This instability is characterised by the parameter or equivalently, where T W is the amount of twist integrated over each magnetic field line in a flux rope, S rec is the footpoint area of the magnetic field lines that reconnect to form the flux rope, Φ over is the magnetic flux of the overlying field, and τ twist = T W |B z |.In DAI, κ is usually larger for larger S rec ; therefore, there is a critical value S c at which the instability occurs: κ > κ 0 ∼ 0.1.It is shown that the ratio of the magnetic helicity of current-carrying magnetic field (|H j |) to the total relative helicity (|H R |) well discriminates whether a flare event becomes eruptive or not 60 .While the κ parameter is the critical parameter which determines the onset of DAI (i.e., whether a flare occur), the helicity ratio |H j |/|H R | may be capable of distinguishing the eruptivity when the flare occurs.
In the κ-scheme 49 , the coronal magnetic field is first calculated by the non-linear force-free field extrapolation based on the vector magnetic field of the photosphere.For the computational resources, the grid number was reduced from 1024 × 1024 to 512 × 256 in this process.Then, HiFERs are identified as the regions where the non-potential field B np (= |B − B p |) exceeds the threshold value of 1 kG, which is based on the observational result that B np = 1 kG sufficiently encompasses the the distribution of non-potential magnetic fields that drive large flares 49 .For each HiFER, the critical area S c is measured as the minimum circular area that satisfies the DAI condition (κ > κ 0 ), and r c is obtained as the radius of S c , i.e., S c = πr 2 c The releasable energy for each HiFER is estimated as where S r is the area of the footpoint of the magnetic flux that pass over the circular area S c .

Figure 3 .Figure 4 .
Figure 3. 3D view of the magnetic structures below the two rotating sunspots, whose field of view is indicated in Figure 1.The reddish volume rendering shows the total field strength (|B|), while the yellow lines are the selected magnetic field lines.The 2D slice at the top presents the photospheric vertical field (B z ) at the τ = 1 surface, with the grayscale saturating at ±2 kG.The three panels on the right show the velocity fields (V z and V h vector) at −5, −10 and −20 Mm.Contours (black or green depending on the background) show the 50% and 80% levels of the maximum total field strength (|B|) of the entire horizontal (i.e.98.3 Mm × 98.3 Mm) plane at each depth, representing the two magnetic pillars.Cyan arrows highlight the directions of the horizontal flows.All physical values are obtained by averaging the data sets over a period between t = 35.0hr and 38.0 hr.