On the origin of the solar hemispherical helicity rules: Simulations of the rise of magnetic flux concentrations in a background field

We define our SRR as the region in parameter space where one relative orientation of twist and background field rises successfully but the reverse orientation does not. Seeing as the dynamics of the selection mechanism depend on the interaction between the background field and the twist of the tube, it might be expected that the SRR is most strongly dependent on $q$ and $B_{s}$, and we concentrate on this first, fixing for now the other main parameter, the scale height of the background field configuration, at $H_{B}=0.125$ (and, as always, keeping the other parameters at their canonical values: $m=1.5,~{}R=0.125,~{}\theta=0.25$). It is perhaps easiest to think of the SRR as the parameter regime where one sign of twist rises but the other does not at fixed background field strength and orientation. Note, however, that our original simulations fixed the twist and examined the effect of reversing the field. The SRR is, therefore, really a two-dimensional region of $q-B_{s}$ parameter space (for fixed $H_{B}$) and we investigate this here.

Figure 6 shows a much broader survey of results in the $q-B_{s}$ parameter space than was shown in Paper I. The original work of Paper I would correspond to a single row at $q=2.5$ (not actually shown in this new figure). Schematic indications of the crucial twist and background field orientations are shown in the border of the table, along with their values, for ease of interpretation of the signs. The figure shows a matrix colored at each $q-B_{s}$ value according to whether a simulation at those values shows that the flux concentration clearly rises (green), clearly fails to rise (red), or something less easily determined (orange). The determination in each simulation is done using visualisations and the flux fractions, as exemplified in Figures 3 and 4. The relatively complicated intermediate dynamics represented by the orange colored cases are examined in more detail in Section 4.4. Swapping both the sign of $q$ and $B_{s}$ results in the same dynamics (equivalent to just viewing the same system from the other end of the flux tube), and so the table is diagonally-symmetric across the origin. All four quadrants of the matrix have been included for ease of using the table, but this symmetry emphasizes that the selection mechanism only depends on the relative direction and strength of $B_{s}$ and the twist, $q$. We, therefore, explain one half of the survey chart ($Bs>0$) since analogous interpretations can be drawn for the other half. Note that in our model setup, when $q$ changes sign, the current helicity of the tube changes sign, since $B_{z}$ remains fixed in the positive $z$ direction. This is irrelevant to the selection mechanism, since it depends solely on twist, not helicity. For relevance to the solar case, however, the selection mechanism can be cast in terms of the helicity, since, there, the sign of the background field and the tubes axial field are correlated, according to our understanding of the dynamics of the solar dynamo. The cases within the white dotted lines are special cases. For $q\lesssim 0.25$, all cases fail to rise coherently regardless of the background field strength because some initial twist is required to maintain the coherency of the tube (Moreno-Insertis & Emonet, 1996). For all other $q\gtrsim 0.25$, but zero background field, the tube rises successfully.

For any fixed $B_{s}$, the table in Figure 6 clearly shows that an asymmetry exists between positive and negative $q$. For $B_{s}>0$, the transition from red (unsuccessful rise) to green (successful rise) for negative $q$ occurs (if it occurs at all in the cases surveyed) at far larger $|q|$ than for positive $q$. By symmetry, for $B_{s}<0$, the roles are reversed with the $q>0$ boundary occurring at much larger $|q|$. The different $|q|$ at which the red-green (unsuccessful-successful) transition occurs (say $|q_{-}|$ for the negative $q$ and $|q_{+}|$ for the positive $q$) delineate the SRR at that $B_{s}$. The SRR can therefore be defined from the simulation data at any $B_{s}$ as the set

Within this set, $q$ of one sign rises and the other does not. It can easily be verified from Figure 6 that, for $B_{s}>0$, it is structures with $q>0$ that rise if their $q$ lies in the SRR found at that $B_{s}$, whereas for $B_{s}<0$, it is the structures with $q<0$ that are the ones that rise.

Figure 6 shows that the boundaries of the SRR, given by $|q_{-}|$ and $|q_{+}|$, vary with $B_{s}$. For higher $|B_{s}|$, stronger $q$ is clearly required to overcome the background field in order to rise, and the SRR becomes defined by higher $|q|$ values and appears to widen in extent of $|q|$. These issues are investigated in detail via a numerical model in the next section.

At this point, it is worth noting that there are some physical limits on what we might expect for reasonable values of $|q|$, despite the fact that we have no direct observational evidence from the solar interior to help us out in this matter. Firstly, in our model, structures with $|q|\gtrsim 12$ are unphysical, since then the azimuthal field is sufficiently strong that the associated magnetic pressure induces a negative density in the initial conditions for the structure ²²2This limit is actually not constant but instead a very weakly decreasing function of $B_{s}$, and is shown later as the dashed line in Fig. 8. Secondly, structures with $|q|\gtrsim 8$ are likely unstable to kink instabilities (according to Linton et al. (1996), whose three-dimensional model has a different field configuration from our two-dimensional model, but is likely a reasonable guide). These issues informed our choice of range in $q$ in these simulations. Similarly, we have no direct observational data for the strength of solar interior magnetic fields, but it seems reasonable to presume that the background poloidal field strength is only a fraction of the strength of any rising flux structure.

Ultimately, the $q-B_{s}$ diagram in Figure 6 exhibits compactly the results of a selection mechanism and therefore illustrates clearly what we have referred to as the SRR, the parameter regime in which the selection occurs. The overall conclusions from the simulation data presented in Figure 6 may be that an SRR of some extent occurs for almost all reasonably physical parameter values. This seems to imply the existence of an overall bias for the rise of flux structures, even if they have significant variance in their magnetic composition. In general, a prevalence of positively twisted structures might be expected to emerge from a distribution of structures in the presence of a positive background field, and a prevalence of negatively twisted structures might be expected in a negative background field. We verify these ideas in Section 4.4, but first we turn to an analytical model that defines the SRR explicitly in terms of the parameters, allowing a more fine-grained exploration and understanding of the parameter space.

Our previous paper, Paper I, provided some brief insight into the reasons for the varying dynamics associated with different orientations of the fields and tubes. Here, we extend those ideas into a detailed model that can explain and predict the $q-B_{s}$ behaviour discussed above, and the dependence on other parameters. To formulate the model, we carry out a theoretical analysis of the key magnetic forces acting on the flux concentration in the vertical $(y)$ direction as it initially begins to rise. Our model is based on the following ideas. Firstly, we have in mind that the rise of a tube-like concentration through the overlying background field is generally hindered by tension induced in the background field as it is wrapped around the structure and gets lifted and stretched upwards. Secondly, originating from the ideas in Paper I, we expect that the introduction of the background field induces differential values of the internal tension and buoyancy forces across the tube-like concentration which can lead to net forces that can help or hinder the rise. To create a model that predicts whether rise is successful or not, we calculate the differential buoyancy and tension forces from the initial conditions at $t=0$ and then estimate the impeding tension forces from the overlying field that are induced as the structure rises. A balance of these forces should delimit the bounding case between successful and unsuccessful rise. Note that, since the selection mechanism and all the results discussed so far are robust to different formulations the magnetic initial conditions (see Section 3 and Appendix A), in order to make analytic integration of the net forces easier, we revert back to the “top-hat” magnetic initial conditions ($B_{z}=1$ for $r\leq R$ in the concentration; see Eqn 6) as used in Paper I rather than the Gaussian cross-section used in Section 3.

The buoyancy force in the y-direction, $F_{buoyancy,y}$, including its magnetic contribution, is generally thought of as the driving force for the rise of the concentration. We have deliberately set up our problem so that the presence of the magnetic field of the concentration produces a density (and temperature) perturbation (under the assumption of fast equilibrium of the total pressure; see Appendix A, Case 2) that will drive the rise. The addition of a background field that varies with height can create a perturbation to the vertical buoyancy force that is asymmetric about the center of the tube and therefore has a net contribution. Ultimately, however, we will find in the following calculation that the magnetic contribution to the buoyancy force is dominated by the axial field of the tube, and that this differential perturbation is unimportant (for the buoyancy force, although it will be crucial for the tension force). The net buoyancy force in the vertical (y) direction over the circular region of the tube-like concentration in our Cartesian coordinate system is given by

The final term in the integral is related to the cross term $(\sim B_{x}B_{back})$ and yields a difficult term in the integration. We drop this term here, expecting it to be relatively insignificant, and indeed verify that this is the case later. In order to integrate the first term in equation (16), we use Equation (5) to substitute

and then, after integrating equation (16) once with respect to $x$, we further use the definition of modified Bessel function of first order,

with $t=y/R$ and $z=R/H_{B}$. The result of the overall integration of Equation (16) (excluding the cross term) is then

where $\tau=R/(2H_{B})$. Equation (19) shows that the net buoyancy force has components that depend on the strength of the background field (through the first term and $B_{s}^{2}$) and the tube field strength (through the second term, dependent on $q^{2}$ and a constant value for the axial field). When evaluating this expression for our parameters later, we will see that the second term clearly dominates.

When a background field is included, the role of the internal tension force in the overall dynamics can be equally as important as the buoyancy force in determining the net upwards driving force. When there is no background field ($B_{s}=0$), the internal tension force in the tube is symmetric and acts towards the center of the tube. Hence, there is no net force in any direction at $t=0$ and tension serves only to maintain the coherency of the tube in this case. On the introduction of a background field that varies in the vertical, this internal tension force becomes asymmetrical around the horizontal mid-line of the tube, and, as pointed out in Paper I, the net vertical force induced can be upwards, thereby aiding the rise, or downwards, thereby hindering the rise. We denote this internal tension force in the vertical (y) direction by $F_{tension,y}$, given by

where $y$ has been translated as before. With the chosen magnetic initial conditions given in equation (6), $\dfrac{\partial B_{y}}{\partial y}=0$ and so, at $t=0$,

The third part of the force budget that we need to consider is the tension force induced in the large-scale background field by the rise of the structure upwards into this field, which we denote by $F_{back,y}$. This opposes the rise of the tube. To estimate $F_{back,y}$, we examine the tension force induced in a rise that proceeds to a general height, $H$. For a successful rise in our simulations, $H$ will be approximately equal to or greater than the maximum vertical size of the simulation box, $y_{max}=4$. Tension forces in curved field lines can be written as square of the magnitude of the field divided by the radius of curvature, pointing inwards to the curve. If we assume that a rising tube lifts all of the field it encounters during the rise and concentrates it into a narrow region wrapped around the tube, then a description of the net tension force (in the $y$-translated coordinate system used earlier) opposing the vertical rise induced during the rise might be reasonably given by the integrated square amplitude of the field divided by the radius of the tube:

For small $H$, i.e. $H<<H_{B}$, we can rewrite the exponential function using series expansion and keep only terms up to first order, giving

Clearly then the tension force induced in the background field is linear in $H$ and will be negligible when the height to which the flux tube rises, $H$, is very small, as expected. When the flux tube rises more significantly, i.e. $H>>H_{B}$, the net induced tension force tends to a constant value

The tension therefore only varies rapidly over a height of about a scale height $H_{B}$ and its role can be bounded above nicely by equation (28).

If we make the assumption that the internal buoyancy and tension forces in the tube-like concentration do not change much from their initial values, then the marginal case that defines the transition from a failed rise to a successful rise is given by a balance of the three elements described above, i.e.

Assuming some non-negligible rise, $H>>H_{B}$, using Equations (19), (28) and (23), this becomes

Table 1 shows examples of the numerical values of these roots for some selected $B_{s}$ and $H_{B}$. Figure 7 accompanies the table and shows the relative contributions of the individual components of magnetic buoyancy, internal tension and background field tension forces, plus their total, from the model equations (19), (23) and (28) respectively, as a function of $q$ for the canonical background field configuration with $H_{B}=0.125$ and the $B_{s}$ indicated. We will use these cases as illustrative examples. A more complete evaluation of the roots of the canonical case is given later in Figure 8, and other $H_{B}$ values are addressed in the next section. In Table 1, we see that for very weak values of $B_{s}$ $(B_{s}\leq 0.03~{}\text{at}~{}H_{B}=0.125)$, there exists no real roots for $q$. This is consistent with the fact that all values of twists, whether positive or negative, will rise successfully when only a very weak background field is present (assuming that $q$ is above the minimum threshold necessary to maintain coherency of the structure)³³3This is likely similar to the $q\gtrsim 0.25$ threshold for the $B_{s}=0$ case, but note that the presence of background field may adjust this.. Figure 7a illustrates this result, as $F_{total,y}$, the sum of the forces on the left hand side of equation (30), does not cross the horizontal axis for any $q$ at this $B_{s}$. For intermediate and high values of $B_{s}$, we get two real roots and these are asymmetric in $q$. For example, we find that the roots for $B_{s}=0.10$ are $q=(q_{1},q_{2})\approx(-19.99,~{}0.17)$ and for $B_{s}=0.16$ are $q=(q_{1},q_{2})\approx(-32.75,~{}1.04)$. This asymmetry in the roots defines the model definition of the SRR. For example, if $q_{1}<0$ and $q_{2}>0$, then rise occurs for structures with $q<q_{1}$ and $q>q_{2}$ and there is a region of $|q|$ between $|q_{1}|$ and $|q_{2}|$ where structures possessing $+|q|$ and $-|q|$ behave differently ($+|q|$ rises, $-|q|$ does not). The canonical examples from Paper I for $B_{s}=0.10$, showed that $q=2.5$ rises but $q=-2.5$ did not, and so $|q|=2.5$ lies in the SRR of this $B_{s}$, which we can see here is given by $0.17<|q|<19.99$. More generally, in these examples with positive $B_{s}$, positively-twisted flux tube that are fairly weak can rise whereas negatively-twisted flux structures of the same strength (or even relatively much stronger) cannot.

Figure 7 clarifies the origin of the asymmetrical roots and therefore the SRR in terms of the individual force contributions. The buoyancy force (black lines derived from Equation (19) with associated circles from simulation data – see later) is quadratic in $q$ and can now be seen to be dominated by the contribution from the initial density perturbation due to the magnetic structure (the second term on the right hand side of Equation (19)), since the buoyancy force is very close to symmetric in $q$ and only very weakly dependent on $B_{s}$. The asymmetry in the $q$ dependence of the total upwards force on the tube (green line; diamonds) instead stems from the net vertical internal tension force in the tube (red line; squares). This is linear in $q$, and significant in size because it depends on $B_{s}$ and $H_{B}$ rather than these quantities squared. For the $B_{s}>0$ shown here, negative $q$ creates a negative tension force serving to reduce the effectiveness of the buoyancy force, retarding the rise of negative $q$ structures. This effect can completely inhibit the rise unless the buoyancy force (proportional to $q^{2}$) is large enough to compensate. Where $q>0$, the internal tension force is positive and aids in the rise of the structure. The tension in the background field (blue line; no symbols) is independent of $q$ and only serves to reduce the chances of rise evenly across the range of $q$ at fixed $B_{s}$ (but is dependent on $B_{s}^{2}$). The net force, given by the left hand side of equation (30) (green line; diamonds) is therefore a shifted and tilted quadratic function. Where this function passes through zero are the roots $q_{1}$ and $q_{2}$, and now it is clear that the effect of the background field on the internal tension of the tube is responsible for introducing the major asymmetry into this function and therefore is the culprit for creating the asymmetric roots leading to the existence of a SRR.

The results have been discussed here only for $B_{s}>0$. The results for $B_{s}<0$ can be inferred by the symmetry in the $q-B_{s}$ space, as noted earlier with relation to the simulation data in Fig. 6. For example, in an equivalent of Fig. 7 for $B_{s}<0$, the buoyancy force and the tension force induced in the background field would remain essentially the same, but the internal tube tension force would switch sign. The roles of $q_{1}$ and $q_{2}$ would therefore be reversed, with the asymmetry such that $|q_{1}|<<|q_{2}|$. A general definition of the model SRR such that it is valid for any orientation of $B_{back}$ in terms of the roots of the quadratic is given similarly to Equation 12 by

Figure 8 shows an example of the delineation of the SRR in $q-B_{s}$ space via this method (at the canonical parameter values, in particular $H_{B}=0.125$; other $H_{B}$ are done shortly). The red diamond points show $\text{min}(|q_{1}|,|q_{2}|)$ values and the black circular points show $\text{max}(|q_{1}|,|q_{2}|)$. For these cases with $B_{s}>0$, the root of smaller modulus is $q_{2}$. The SRR is the region between the lines. Structures with values of $q$ of either sign such that $|q|$ is above the upper limit of the SRR (line of black points) all rise. For $B_{s}\gtrsim 0.09$, where $q_{1}$ and $q_{2}$ are of opposite sign, structures with both signs of $q$ would fail to rise if $|q|$ is below the lower limit (red diamonds). For $0<B_{s}\lesssim 0.09$, where $q_{1}$ and $q_{2}$ are both negative, both signs would rise for $|q|$ below the red line. Only within the SRR range do structures of one sign rise whereas structures with the opposite sign do not. Positive twist rises preferentially for the cases shown here with $B_{s}>0$, but negative twist would rise preferentially if $B_{s}<0$. The SRR definition given in Equation (31) is presented as a range of $q$ at fixed $B_{s}$ and is therefore a function of $B_{s}$. The SRR is actually a region in $q-B_{s}$ space and could equally well be presented as a range of $B_{s}$ that is a function of $q$, as is clear from both the data and the model, as shown in Figs. 6 and 8.

The model proposed therefore seems to explain well the nature of the dynamics that have been observed. However, a number of simplifications were introduced, namely: the assumption that the background temperature variation across the tube was reasonably represented by its average; that it is reasonable to omit the $q$ dependent term in equation (16); and, that it is a good approximation to estimate the internal tube forces (buoyancy and internal tension) only at $t=0$ and then use that to deduce the subsequent $t>0$ rise behavior. Since we have all the data from the simulations, it is constructive to compare what we can in the model with the numerical simulations in order to verify our assumptions. We therefore start by crudely calculating the net internal forces in the structure at $t=0$ by summing over the relevant quantities inside the tube area as follows:

where $r_{i,j}=x_{i}^{2}+(y_{j}-y_{c})^{2}$ with $i,j$ the indices of the discrete simulation grid in the $x$ and $y$ direction respectively. For each relevant simulation that we have performed, we plot the values of these net internal forces as the appropriately-colored symbols in Figure 7. Internal magnetic tension (red squares) calculated using the simulation data exactly matches with the model estimate. The black circles representing the internal tube buoyancy from the simulations do have a slight linear dependence on $q$ added to the symmetrical quadratic term, but so slight that it is not visually apparent in the figure, vindicating our omission of this term in the model.

The major source of uncertainty in our model then comes from the estimate of $F_{back,y}$. This quantity is a very rough approximation to one aspect of more complex dynamics. Omitted effects are the deformed geometry of the structure, effects of the wake and the overlaying field on the wake, any drainage from the main structure, and diffusive effects, at least. These other effects are hard to model and would be hard to check against data even if they were modeled, since they are not compactly spatially located to the structure and are time-evolving. The discrepancy associated with these errors shows up in comparisons of the roots that we obtain by solving the quadratic Equation (30) and the SRR range deduced from observing the simulations. The upper and lower bounds of the SRR deduced from the simulations are shown for a few sample $B_{s}$ as the green circles in Figure 8. Although the simulation SRR boundaries are only visually extracted from a coarse set of simulations and are therefore relatively inaccurate, it is clear from this comparison that the simulated SRR is substantially narrower than the model predicts. It appears that our estimate of the rise-quenching effects from the background field is indeed quite a severe underestimation, as we imagined. We could attempt ad hoc changes to the model (e.g. increasing the radius of curvature of the retarding external tension force in the background field to reduce its effect and to account for the topology of the rising structure later in its rise), or try to calibrate the model to the data, but since the simulations are fairly coarsely spaced in parameter ($q-B_{s}$ space) and the boundaries between successful and failed rise are sometimes not entirely obvious, then this seems unlikely to be particularly accurate. Furthermore, the omitted effects are probably dependent on the parameters (such as $B_{s}$ and $H_{B}$) too, and therefore a simple calibration is unlikely to work universally. However, even with some significant uncertainty, the model still provides a good explanation of the dynamics behind the SRR, and gives useful predictions of the trends seen in the simulations, at least. In particular, it appears to be a robust result that the extent of the SRR grows and the root of lower modulus grows (in signed value) with increasing $B_{s}$ (see also Fig. 12 later). These facts have a significant influence on the expected distribution of the signs of the emerging helicity, as will be discussed in detail in Section 4.4.

We have now established the importance of the presence of the large-scale background field on the dynamics of a rising flux tube, concentrating so far on the relationship between the background field strength and the twist strength, and their relative orientations. We now examine the dependence on the configuration of the background field, which we have characterized as an exponential with scale height $H_{B}$. Figure 9 shows the variation of the background field, $B_{Back}$, as a function of height, $y$, for some example values of the scale-height, $H_{B}$, when $B_{s}=0.1$. The dotted horizontal lines in the figure show the initial vertical location of the flux concentration. Recall that $B_{s}$ specifies the field strength at the center of the tube so all these lines intersect at that point with value $B_{s}=0.1$. This figure is intended to emphasize the fact that, when $H_{B}$ is small, the background field variation across the flux tube is substantial, whereas for high values of $H_{B}$ the background field changes very little with height overall and hence its variation across the tube is small. Since this differential is key to the new dynamics discussed above, we might expect $H_{B}$ to play a role.

Figure 10 shows intensity plots of the (normalized) axial field $B_{z}$ at three different times in the evolution of system for three values of scale height, $H_{B}=0.125,~{}0.175,~{}0.300$, whilst keeping the background field strength and twist of the flux tube constant (at $B_{s}=0.10$ and $q=2.5$). Overall, Figure 10 shows that the rise of the flux structure is quenched for higher $H_{B}$. Figure 10$a$, exhibiting the evolution for the canonical $H_{B}=0.125$, shows that a typical flux structure – a “squashed” head with trailing vortices – rises through the background field towards the top of the simulation domain. Figure 10$b$, for $H_{B}=0.175$, shows another overall successful rise, but this time there is considerable drainage of the axial flux $B_{z}$ away from the tube head down the the overlaying background fieldlines, plus some significant change to the geometry of the head and the trailing vortices. Note that the flux tube is further from the top of the simulation domain at $t=20$ in Fig. 10$b$ and is therefore rising more slowly compared to the previous case, even though it appears to be still rising as a fairly independent structure. Figure 10$c$ shows the case with an even larger scale height, $H_{B}=0.3$. In this case, there does not appear to be any coherent rising flux structure.

The quenching of the rise of the flux structure with the increase in scale height of the background field can be understood from a physical perspective, with confirmation from our mathematical model. One can imagine that increasing the scale height of the background field could lead to three effects that could potentially affect the dynamics. Firstly, a larger $H_{B}$ redistributes the background field in the tube and this affects the magnetic pressure that drives the rise of the tube. The analytical approximation in Equation (19) reveals that this dependence is quite complex, being proportional to $H_{B}^{2}I_{1}(R/H_{B})$. This is a function that increases with $H_{B}$ and therefore its contribution to the magnetic buoyancy is enhanced at larger $H_{B}$. However, as we discussed earlier, this magnetic pressure contribution is much smaller than the contribution from the initial density perturbation in the buoyancy term proportional to $g=\theta(m+1)$. This first effect is therefore likely not very significant.

The second potential effect is that an increase in $H_{B}$ leads to a decrease in the differential of the amplitude of the field across the tube vertically. It is this differential that drives the selection mechanism via its effect on the magnetic tension in the tube. Any effect (either increased likelihood of rise when $B_{s}$ and $q$ are of the same sign, or decreased likelihood when $B_{s}$ and $q$ are of opposite signs) is reduced for larger $H_{B}$. In the case shown in Figure 10, we might expect the tension-driven enhancement of the rise to be reduced by larger $H_{B}$. This expectation can again be examined with the analytical model, this time via Equation (23). There, it can be seen that the direction of the force depends on the sign of $qB_{s}$ and also that the amplitude again has a complicated dependence on $H_{B}$: $F_{tension,y}\sim H_{B}I_{1}(R/2H_{B})$. This dependence is almost constant, only very slightly decreasing over a significant range of $H_{B}$. This second effect is therefore again probably not very important.

The third effect remains as the likely major dynamical factor. Increasing $H_{B}$ also serves to increase the integrated field that contributes to the retarding force of the overlaying background field. One might expect larger $H_{B}$ to imply a larger value of the integrated field and therefore a large force countering the rise. Since this is an integral over the distance travelled upwards by the tube, this is probably a substantial effect. These ideas can again be verified by the analytical model, where this time the dependence can be seen from Equation (28) to be straightforwardly linear: $F_{back,y}\sim H_{B}$. As mentioned previously, this estimate is also potentially an underestimate, with some significant factors related to other effects omitted. This also lends some credence to this being the dominant physical term here, although we do not know the $H_{B}$ dependence of the other effects.

Having established some physical intuition, we can further employ the analytical model developed in Section 4.2 to analyze the broader dependence of the SRR on the scale height of the background field, $H_{B}$. Figure 11 plots the SRR from Equation (31) in $q-B_{s}$ space for fixed $H_{b}$ (Figure 11$a$) and in $q-H_{B}$ space for fixed $B_{s}$ (Figure 11$b$). The blue vertical dashed lines in each of the subplots show the values where numerical simulations exist for comparison with the model. For the values of $H_{B}$ surveyed, by comparing the panels of 11$a$ or examining an individual panel of Figure 11$b$ at fixed $B_{s}$, we see that the demarcation (and therefore the extent) of the SRR (in $|q|$) only varies a very moderate amount. This is commensurate with our earlier conclusion that the tension effect that drives the existence of the SRR is only very weakly dependent on $H_{B}$. On the other hand, individual panels of Figure 11$a$ or comparison of the panels in Figure 11$b$ confirms that varying the background field strength, $B_{s}$, can significantly change the extent of the SRR for fixed values of $H_{B}$. This is all again clarified by the analytical model where the tube internal tension term in Equation (23) that establishes the SRR is linearly dependent on $B_{s}$, but is roughly constant with $H_{B}$ (via the term $H_{B}I_{1}(R/2H_{B})$). An important, and perhaps unexpected, result overall here is that an SRR still exists even for high $H_{B}$, although the $q$ values at which it occurs become less physical.

It now becomes necessary to reconcile the above statements with the dramatic quenching of the rise of a structure with increasing $H_{B}$ observed in Figure 10. First, it should be noted again that we expect realistic values of the twist $q$ to be of order $q<10$. Higher values create unphysical (negative) values of the density (see horizontal dashed line in Fig.  8) and other work has shown that high $q$ are unstable. For instance, the initial Gaussian profile for the axial field used in Section 3 would have a critical value of $q_{cr}=1/R=8$ for the flux tube to be kink unstable in three dimensions (Linton et al., 1996). This restriction means that realistic flux structures probably lie much closer to the lower limit of the SRR than the upper. If $q$ is such that it is close to the lower boundary, then even a small variation in $H_{B}$ and a weak dependence of the SRR boundary on $H_{B}$ can have a dramatic effect. This is essentially what was demonstrated through the simulations shown in Figure 10. The simulations in this figure are at $q=2.5$ and $B_{s}=0.10$, corresponding to the second blue dashed line in Figure 11$a$ and the middle panel in Figure 11$b$. In the latter, for $q=2.5$, the moderate upwards tilt of the lower boundary of the SRR (red diamonds) with increasing $H_{B}$ shifts the solution from inside the SRR (for $H_{B}=0.125$) to below the SRR (for $H_{B}=0.3$), thereby switching it from a rising solution to a quenched solution. A small change in $H_{B}$ can therefore have an effect due to the proximity of natural solutions to the lower boundary of the SRR.

To emphasize the above dependencies, Figure 12 plots the width of the SRR, $||q_{1}|-|q_{2}||$, as a function of (a) $B_{s}$, for a range of $H_{B}$, and, (b) $H_{B}$, for a wide range of $B_{s}$ (as in the previous figure). These figures solidify our interpretation that the width of the SRR is only very weakly dependent on $H_{B}$ but significantly dependent on $B_{s}$. Note again that at very low values of $B_{s}$ (e.g. $B_{s}=0.03$) there is no SRR since the differential tension force required to generate the necessary asymmetry for the SRR is negligible.

With an understanding of the dynamics of a selection mechanism established for a single flux concentration, in this section we examine the collective effect that leads to the observations that are known as the SHHR. The SHHR observations accumulate the signs and strengths of the helicity of emerging active regions over the whole solar cycle (and eventually multiple cycles) to extract any net bias. The SHHR observations have a large “scatter”, in the sense that not all individual active regions agree with the rules – overall, only 60-80% of active regions appear to concur – and the degree of agreement appears to vary over the solar cycle (Seehafer (1990); Pevtsov et al. (1995); Abramenko et al. (1997); Bao & Zhang (1998); Bao et al. (2000a, b); Hagino & Sakurai (2005); Hao & Zhang (2011); Singh et al. (2018)). The scatter and its temporal variation is not well understood. In this section, we therefore attempt to use our model to create synthetic SHHR observations with a view to potentially elucidating reasons for these effects.

To reproduce SHHR observations, we simulate the evolution of multiple (toroidal) flux structures through different background (poloidal) field strengths representing different parts of the solar cycle. We fix the values of the other parameters at the canonical values: $m=1.5,~{}H_{B}=0.125,~{}R=0.125,~{}\theta=0.25$. For each field strength, we perform a series of Monte Carlo (MC) simulations where, in each simulation, a wide domain is initialised with a number of flux tubes with randomly assigned twist strength (both positive and negative) and horizontal location. We then examine whether the selection effects shown for single flux tubes leads to a persistent bias over the whole cycle (simulated by our range of background field strengths). It may seem likely that this is the case as we have identified a selection mechanism, but this needs to be confirmed, since this setup allows the possibility of new dynamics emerging from the interaction of tubes, and the randomness may potentially swamp any bias.

In order to have reasonably significant statistics for each case in this MC study, we perform fifty independent simulations (each containing multiple tubes) for each background field strength. Since we are not addressing the origin of the flux structures, our assumption is that they are formed from a sheet of toroidal field, likely by a magnetic buoyancy instability (see for e.g. Vasil & Brummell (2008); Guerrero & Käpylä (2011)). If this were to be the case, we might then expect multiple flux structures to be formed at the same time in a turbulent environment, and such structures would most likely possess random strengths, twists and separations drawn from some distribution. We therefore choose each of our simulations in the MC series for a particular background field strength to start with five flux concentrations with random initial twists $(-5\leq q\leq 5)$ and random distance, $d$, between them (with the constraint, $d\geq 2R$, so that they do not overlap), embedded initially at the same (canonical) height in the prescribed overlying large-scale background field. The random values are drawn from a uniform distribution over the range given. Our setup certainly has some significant degree of arbitrariness because of the limitations in our quantitative understanding of the amount of twist in flux tubes, the distance between them and their strengths, at their origination in the deep solar interior. For example, our choice of drawing the random values of twist from a uniform distribution might not be very realistic and a Gaussian with a larger range might be better (or perhaps even a skewed distribution to account for the effects in the formation of tubes not included in our model, e.g. rotation: see Chatterjee et al. (2011)). However, as mentioned earlier, there are some constraints on the value of $|q|$. Firstly, $|q|\gtrsim 12$ (the dashed line in Fig. 8, which is actually a very weakly decreasing function of $B_{s}$) is unphysical, since then the azimuthal field is sufficiently strong induce a negative density in the structure. Secondly, structures with $|q|\gtrsim 8$ are unstable to kink instabilities (according to Linton et al. (1996), which does not fit our model exactly, but is likely a reasonable guide). Therefore, our chosen distribution (uniform, $|q|\leq 5$) seems reasonable, in that it makes it likely that the initial simulation setup has a mixture of weak and reasonably strong concentrations of both helicities, and, if anything, overemphasizes the stronger twists (that are more likely to violate the SHHR). An MC series of this type (fifty simulations, five random tubes each) is reported for four different background field strengths, ranging from relatively weak to relatively influential ($0.01\leq B_{s}\leq 0.10$). We identify the weak background field with the beginning and end of the cycle and the stronger background with the peak of the cycle, although we have no direct observational evidence from the deep solar interior of the relative strengths of these fields from which to set their values.

Fig. 13 shows the results of the MC simulations for each of the four different background field strengths. Each circle in the plot represents an accumulation over a specified time interval $[t_{0},t_{1}]$ in the simulation of the absolute value of the z-component of the net horizontal current helicity measured along a line at $y=3$. We call this quantity $H_{3}$:

This quantity is evaluated from the simulations as a sum over the $x$ grid points at height $y=3$ at a particular time, and over the data output in the time interval from $t_{0}=0$ to $t_{1}=70$. This quantity represents the magnitude of the net helicity of all the flux concentrations that successfully rise through $75\%$ of the model solar interior in the simulation setup. The open and filled circles indicate the net sign, corresponding to net negative and positive helicity respectively. Each panel in the plot shows the 50 MC simulations for a given field strength, $B_{s}$, representing a given time in the solar cycle, plotted against the simulation number, which could be loosely interpreted as short time intervals at that point in the cycle, but are in reality just different realizations whose ensemble average becomes meaningful. The overall bias in the emerging helicity at any specific time in the solar cycle is therefore given by any statistically-significant deviation of the percentage of each sign of helicity away from a $50\%-50\%$ balance. These percentage values are noted on the plots. Note that owing to the orientation of the tube’s axial field along the positive z direction and the background field along the the positive x direction, our setup corresponds to the the southern hemisphere in the first half of the cycle as shown in Fig. 5a, where positive helicity is expected to dominate by the SHHR. Furthermore, to relate to observations, we are assuming that some three dimensional process leads to arching of our simulated flux concentrations for emergence so that the axial helicity that we measure horizontally in 2D is strongly correlated with the vertical helicity component that emerges in the observations. We will test these assumptions in 3D in the future.

Fig. 13 displays a systematic trend in the relative percentages of the two signs of accumulated helicity with increasing background field strength. For a very weak background field strength of $B_{s}=0.01$ in Figure 13a, we find that the surviving and accumulating helicity at $y=3$ is roughly equally divided amongst both signs of helicity. There is a slight bias towards an anti-SHHR result (more negative helicity) as has actually been observed in the solar case, but we do not dwell on this as this seems fairly statistically insignificant here, since the number of MC simulations is only 50. At higher background field strengths (which we relate to closer to solar maximum in the the solar cycle), we find results that are definitely commensurate with the SHHR. At $B_{s}=0.03$, we find that around $66\%$ of the successfully rising flux concentrations possess positive helicity; at $B_{s}=0.05$ we find that $78\%$ are positive; and at $B_{s}=0.1$, $86\%$ are positive. These results are in agreement with the overall SHHR trend, which is that positive helicity should generally be preferred in this Southern hemisphere setup (and, of course, the results for the Northern hemisphere can be inferred by symmetry). Furthermore, these results also agree with the observed trend that there are fewer violations to the SHHR (i.e. better agreement) around solar maximum and more violations (worse agreement) around solar minimum (Hagino & Sakurai, 2005). Note that we correlate a stronger deep interior background poloidal field in our model with solar maximum and a weaker poloidal field with solar minima. Since we do not know the deep interior fields, we rely upon mean-field dynamo models to justify this (see e.g. Charbonneau, 2020, and references there-in). Many such models exhibit that the poloidal and toroidal fields are reasonably in phase, and the latter are widely used as a proxy for the butterfly diagram.

These major results can be interpreted fairly straightforwardly in the light of the work in Sections 4.1 and 4.2, and, in particular, the $q-B_{s}$ diagrams from simulations (e.g. Fig. 6) and the model (e.g. Figs. 8, 11 and see also Fig. 12). In general, these plots show that there is only a narrow SRR at weak $B_{s}$ but that the SRR grows wider in extent in $|q|$ as $B_{s}$ increases (a robust trend, even though the data and the model do not agree precisely on location of the boundaries of the SRR due to the inaccuracies of the model, as discussed earlier). Violations of the SRR would mainly be detected when the twist of a structure is strong enough for it to lie above the SRR (although this is not the only possibility, as discussed shortly). Since there are good physical reasons for the value of $|q|$ to be capped above, clearly, where the SRR is narrower at lower $B_{s}$, more violations might be expected to occur, since our distribution of $q$ in the MC simulations is more likely to place them above the SRR region (and below the cap). In the limit of the SRR vanishing, all structures lie outside the SRR and either don’t rise, or are equally likely to rise, therefore yielding percentage values for the helicity sign distribution close to 50%-50%. As the SRR becomes wider at higher $B_{s}$, the number of violations drops as $q$ values drawn from the distribution become less likely to lie above the SRR. This leads to more complete compliance with the SRR and a ratio of the emerging helicity signs that would tend towards 100%-0% in favor of the sign of helicity expected from the SHHR (positive in our example cases). Indeed, for higher $B_{s}$ ($B_{s}\gtrsim 0.06$ in Fig. 8), the cap of our distribution ($|q|=5$) lies below the upper bound of the SRR from both the model and data, and almost 100% agreement would be expected, potentially violated only by very low $|q|$ lying in the small region below the lower bound of the SRR (red diamonds in the model). Structures with such sub-SRR twist are only expected to rise when $B_{s}\lesssim 0.09$ (where both roots $q_{1}$ and $q_{2}$ are negative); violations are not expected for sub-SRR values of twist for $B_{s}\gtrsim 0.09$ since both signs should not rise, but interestingly in this region another method of violation (discussed shortly) operates, reducing the agreement with the SHHR from 100%. Regardless of these details, the trend to fewer violations for stronger $B_{s}$ seems very robust and explainable by the selection mechanism found here.

The above result, related to the expected degree of scatter in the observations and its dependence on the solar cycle, is a reflection of how violations to the SHHR are achieved in the MC simulations. We therefore investigate the causes of these violations in more details now. To exhibit the mechanisms we find, we focus on the cases with $B_{s}=0.05$ and $0.10$. We find three distinct, relevant reasons for violations which are discussed below:

1. 1.

   Strongly twisted flux tube: The first, and main, reason, mentioned previously, is simply the existence of strongly twisted flux concentrations of the “wrong” sign (where, here, “wrong” means that the sign of helicity of the structure violates the SHHR; in the cases shown here, the “wrong” sign is negative since the SHHR would predict predominantly positive helicity). If the structures are sufficiently strongly twisted to lie above the SRR, then these structures will rise. As proof of this from the simulations, Figures 14$a,b$ show examples of the time evolution of normalized axial magnetic field strength, $B_{z}$, from the MC simulations performed with $B_{s}=0.05$. The twist sign and strength of individual flux concentrations are shown, numbered by their relative position along the horizontal axis at $t=t_{1}$. Fig. 14$a$ shows a case where the randomly-assigned twists of all the flux concentrations in the simulation are negative (clockwise) and therefore not what would be expected from the SHHR. However, tube 1 has almost the maximum twist allowable from our distribution and we see that at $t=t_{6}$, this tube has risen significantly through the model solar interior and is approaching the $y=3$ level where we measure the net helicity. Note that tubes 3 and 5 are behaving similarly, but to a lesser extent. From Figures 6 and 8 it can be seen that these high values of $|q|$ lie outside the numerically-identified SRR (but inside the potentially inaccurate upper boundary of the model SRR) and therefore might be expected to rise. This is the main source of scatter, and was predicted from our single structure simulations.

2. 2.

   Tube-tube interaction: Fig. 14$b$ exhibits the second class of violation, which is a case where interactions between multiple structures do play an important role. Of the five flux concentrations in this figure, four of them have $q<0$ and the remaining one ($q_{1}$) has $q>0$. The first three, $q_{1},q_{2},q_{3}$, behave as one might expect: $q_{1}$ is moderately strong and positive and therefore from being in the SRR, it is expected to rise; $q_{2}$ and $q_{3}$ are strong but negative and therefore may still rise since they are above the SRR (by the previous method of violation). The other two, $q_{4}$ and $q_{5}$, are very close together, thanks to the use of random distances between concentrations in these MC simulations. The dynamics of the tubes then are no longer those of individual single tubes as interactions between structures become important. Since both these concentrations have the same twist, they coalesce very quickly and form a single, more strongly buoyant region. This structure is still negative but has effectively more twist and more buoyancy, and therefore rises more quickly, becoming a violation of the previous kind due to its strong effective twist. This shows that separation of the structures can be an important element of the dynamics, as the original twist distribution alone would not have predicted the resulting outcome.

3. 3.

   Flux tube twist rearrangement: A third scenario that can lead to violations of the SHHR is particularly peculiar and results from the conversion of weak “correct” signed structures into the “wrong” sign during the rise. An example is shown from one of the MC simulations for the case of $B_{s}=0.1$ in Fig. 15. These types of violations are quite common at higher $B_{s}$ and account for the majority of the violations shown in Fig. 13$d$. Fig. 15 shows the evolution of $B_{x}$, a clear indicator of twist of the tube, in this example. The first three of the structures initially have negative twist, and the remaining two are positive. The first three tubes do not rise successfully since they have only moderate twist strength lying in the SRR and are therefore filtered out by the selection mechanism. Tube 4 has positive twist but is too weak to rise (lies below the SRR). Tube 5 is the peculiar case. This structure has an intermediate positive twist and is expected to rise from consideration of the SRR. However, this structure undergoes a topological arrangement early in its transit, flipping the sign of its twist (visible in the plot as a flip from blue over red in the $B_{x}$ field to red over blue). The exact mechanism for the flip is a little difficult to discern, but appears to be facilitated by the interaction of the strong background field wrapping around into the trailing vortices of the original structure. The structure continues to rise despite this rearrangement, and therefore emerges with the “wrong” sign.