© ESO, 2016

1. Introduction

Sunspots are large concentrations of magnetic flux that appear on the solar surface and are due to the emergence of flux from the convection zone. Rotational motions regularly develop within sunspots in both observations and numerical simulations. This phenomena can trigger explosive events in the Sun’s atmosphere such as coronal mass ejections (CMEs) and solar flares. There has been much debate over the driver of this rotation; several theories have been introduced including photospheric flows caused by differential rotation and magneto-convective dynamics (Brown et al. 2003) and the flux emergence process (Min & Chae 2009; Brown et al. 2003).

Observations of sunspot rotation date back to the beginning of the twentieth century when Evershed (1909) first found spectral observations that exhibited signs of rotating sunspots. Since their discovery, several observational studies have been conducted (Brown et al. 2003; Yan & Qu 2007; Yan et al. 2008; Min & Chae 2009). Through the tracking of prominent features in magnetograms, these studies have shown that sunspots can undergo rotations of the order of a few hundreds degrees over a period of a few days. We direct the reader to Sturrock et al. (2015; henceforth referred to as Paper I) for a more detailed review of recent observations of sunspot rotation.

Following on from the first paper in the series, Paper I, we continue to investigate the process of flux emergence as a cause for sunspot rotation. The net torque produced by magnetic flux emergence was found to be the mechanism for the rotational movement of sunspots. In this paper, we focus on how a change in the parameters determining the magnetic structure of the flux tube affect the evolution of flux and rotation of the spots. Similar flux emergence parameter studies have been conducted in the past with a slightly different focus. Murray et al. (2006) investigated the effects of varying the magnetic field strength and twist of a sub-photospheric cylindrical flux tube with the aim of understanding how these parameters affect the evolution of the flux tube on its rise to the atmosphere. The authors found a self-similar evolution in the rise and emergence of the tube when the magnetic field strength is varied. However, they did not find such a simple self-similar evolution when varying the twist due to the non-linear interaction between the twist and the tension force acting on the tube. Nevertheless, if the field strength or twist is low enough, the flux tube cannot fully emerge into the atmosphere. Another interesting parameter study was carried out by MacTaggart & Hood (2009) where they used a different initial condition, namely a toroidal flux tube. The advantages of modelling a flux tube in the shape of a half-torus, and hence a tube rooted deeper in the interior are two-fold: the axis is able to fully emerge in cases prohibited by the cylindrical model and the sunspot pair does not continue to drift apart. We use this model as our initial condition and instead focus on how varying the above parameters affects the features of sunspot rotation.

Although sunspot rotation has been a very attractive topic to both observers and theorists in recent years, the rotation’s dependence on field strength and twist of the initial sub-photospheric field has, to the best of our knowledge, been left unexplored. Given that there are no current observations of sunspot rotation for varying initial strengths and twists, the results we find based on this simple model should be checked against future observations. We can, however, make predictions on the effect of these parameters from previous studies (Murray et al. 2006; MacTaggart & Hood 2009). For instance, we predict that the twist of the sub-photospheric tube should have a substantial effect on the rate at which the sunspots rotate as we expect a larger rotation for more highly twisted fields, as there is more twist stored in the interior field to unwind. However, as we will discuss later, the density deficit’s non-linear dependence on the twist, α, will complicate this effect. In addition, we predict that the role of the magnetic field strength may not be so important to the magnetic flux tube’s rotation at the photosphere. Nonetheless, the density deficit is directly proportional to the field strength of the tube squared, and hence the stronger fields emerge more fully in our experiments. This allows the axis to align vertically which may impact on the rotation. By considering the effects of changing the field strength we actually study the effects of a differing emerging rate.

In this paper, we present results from 3D magnetohydrodynamic (MHD) simulations of buoyant magnetic flux tubes rising through the solar interior and emerging at the photosphere. We are particularly interested in the rotational motions of the photospheric footpoints of the tube, i.e. the sunspots. We have varied two of the parameters defining the magnetic structure of the sub-photospheric flux tube, namely the magnetic field strength at the axis, B₀, and the twist of the tube, α. Our aim is to identify the effect of these parameters on a number of quantities relating to the rotation of sunspots. In addition, we seek to understand the process controlling the amount of sunspot rotation. To determine the individual effect of these parameters, we vary the field strength and twist, independently of each other.

The remainder of this paper is laid out as follows. In Sect. 2, we briefly outline the model setup and introduce the parameters we choose for the parametric study. A short outline of the general dynamics and main phases of sunspot rotation from a general experiment is given in Sect. 3. Sections 4 and 5 present the results of varying the magnetic field strength and twist of the flux tube respectively. Finally, we conclude the paper with a summary of the results and main conclusions in Sect. 6.

2. Model setup

For the experiments performed in this paper we use the same model initial set-up as defined in Paper I where we solve the 3D resistive MHD equations using the Lare3d code (Arber et al. 2001). Precisely, we model a stratified background atmosphere consisting of a sub-photospheric layer, an isothermal photosphere, a transition region, and an isothermal lower coronal layer. We then insert a twisted toroidal flux tube, introduced by Hood et al. (2009), given in Cartesian coordinates as $\begin{array}{ccc}{\mathit{B}}_{\mathit{x}}& \mathrm{=}& \mathrm{-}{\mathit{B}}_{\mathit{\theta }}\mathrm{\left(}\mathit{r}\mathrm{\right)}\frac{\mathit{R}\mathrm{-}{\mathit{R}}_{\mathrm{0}}}{\mathit{r}}\mathit{,}\\ {\mathit{B}}_{\mathit{y}}& \mathrm{=}& \mathrm{-}{\mathit{B}}_{\mathit{\phi }}\mathrm{\left(}\mathit{r}\mathrm{\right)}\frac{\mathit{z}\mathrm{-}{\mathit{z}}_{\mathrm{base}}}{\mathit{R}}\mathrm{+}{\mathit{B}}_{\mathit{\theta }}\mathrm{\left(}\mathit{r}\mathrm{\right)}\frac{\mathit{x}}{\mathit{r}}\frac{\mathit{y}}{\mathit{R}}\mathit{,}\\ {\mathit{B}}_{\mathit{z}}& \mathrm{=}& \end{array}$(1)where B_φ = B₀e^{− r2/a2} and B_θ = αrB_φ. B₀ is the axial field strength, α is the degree of twist of the fieldlines, a is the minor radius of the flux tube, and R₀ is the major radius. We define the axis of the flux tube to be the fieldline threading through the centre at r = 0.

We must also highlight the difference in pressure and density of the flux tube from its surroundings as this has important consequences later. In order to reach an equilibrium we must establish a pressure excess, p_exc. We make the flux tube buoyant by introducing a density excess by maintaining the temperature and pressure excess as given by ${\mathit{\rho }}_{\mathrm{exc}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\frac{{\mathit{p}}_{\mathrm{exc}}}{\mathit{T}\mathrm{\left(}\mathit{z}\mathrm{\right)}}\mathrm{=}\frac{{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}}{\mathrm{4}\mathit{T}\mathrm{\left(}\mathit{z}\mathrm{\right)}}{\mathrm{e}}^{\mathrm{-}\mathrm{2}{\mathit{r}}^{\mathrm{2}}\mathit{/}{\mathit{a}}^{\mathrm{2}}}\mathrm{(}{\mathit{\alpha }}^{\mathrm{2}}{\mathit{a}}^{\mathrm{2}}\mathrm{-}\mathrm{2}{\mathit{\alpha }}^{\mathrm{2}}{\mathit{r}}^{\mathrm{2}}\mathrm{-}\mathrm{2}\mathrm{)}\mathit{.}$(2)We direct the reader to Paper I for further details. We note that the pressure excess, and in turn density excess, is negative for all r if $\mathit{\alpha a}\mathit{<}\sqrt{\mathrm{2}}$. In all of our parameter choices this is satisfied and so the density excess is negative and is instead a density deficit. Equivalently, this means that the outwardly directed magnetic pressure force is larger than the inwardly directed magnetic tension force, and hence the gas pressure gradient acts inwardly to balance the forces.

Fig. 1 Summary of initial setup. The background density distribution is shown on the right wall, the temperature distribution on the back wall, a selection of fieldlines are shown in red, and an isosurface of the magnetic field (| B | = 1) is over-plotted. The z = 0 plane is also highlighted in grey as the solar surface.

2.1. Parameter choice

In this study, we fix the base of the computational domain at z_base = −25. The major radius of the torus is R₀ = 15 and the minor radius is a = 2.5 for all experiments. We choose to vary the magnetic field strength at the axis of the interior tube, B₀, and the degree of twist, α. The twist is assumed to be positive and constant in all experiments, ensuring that all tubes are uniformly right-hand twisted. Each fieldline rotates about the axis through an angle of α radians over one unit of distance along the axis. The initial set-up of a representative experiment is summarised in Fig. 1, with parameters B₀ = 7 (axial field strength of 9100 G) and α = 0.4 (three full turns of twist in interior tube).

The experiments are split into two groups: Group 1 where α is kept fixed and B₀ is varied and Group 2 where B₀ is fixed and α is varied. A summary of the B₀ and α values under consideration is given in Table 1. It should be noted we only consider a relatively small range of B₀ values as a consequence of our model choice. As found in previous parameter studies (see MacTaggart & Hood 2009), if we pick a lower B₀ value, the flux tube will fail to fully emerge, and if we choose a much higher B₀ value, we may encounter unphysical negative pressures. In fact, if we increase the initial axial field strength to B₀ = 12, we encounter negative pressure in the particular initial set-up outlined in the previous section. It is, however, important to note that the range of allowable field strengths will vary from set-up to set-up. The twist values range from α = 0.2 corresponding to a turn and a half of twist to α = 0.4 that corresponds to three full turns of twist in the initial field. The total flux threading a cross section of the tubes we study range from 3.7 × 10¹⁹ Mx to 7.4 × 10¹⁹ Mx. This is typical of a small active region or large ephemeral region.

All experiments have been performed for different numbers of normalised time units, and hence different final times. The Group 1 experiments are carried out as follows: B₀ = 5 for 216 time units (90 min); B₀ = 6 for 180 time units (75 min); B₀ = 7 for 154 time units (64 min); B₀ = 8 for 135 time units (56 min); B₀ = 9 for 120 time units (50 min); and B₀ = 10 for 108 time units (45 min). This ensures all experiments are performed for the same rescaled time, i.e. . More details of rescaling the time follow in Sect. 4. The Group 2 experiments, on the other hand, are all executed for the same number of time units, 120 time units or equivalently 50 min.

Fig. 2 Radial distribution of the initial magnetic field strength, | B |, at z = −25 for varying a) B0 and b) α.

Fig. 3 Radial distribution of the initial density excess, ρexc, at z = −25 for varying a) B0 and b) α.

Before we proceed to discuss the results of our experiments, we analyse how the parameter choices affect the structure of the magnetic field and resulting density. From the density excess in Eq. (2), it is clear that both the twist and magnetic field strength play an important role in controlling the buoyancy of the flux tube. For Group 1, the variation in B₀ significantly affects the magnetic field strength for all radii until a significant distance from the axis is reached, i.e. the edge of the tube, as shown in Fig. 2a. The initial field strength of the tube is, of course, directly proportional to B₀. The buoyancy profile, too, is strongly dependent on B₀ because the density deficit is proportional to ${\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}$. The B₀ = 10 tube will therefore be 4 times more buoyant than the B₀ = 5 tube (see Fig. 3a).

In Group 2, the variation in α has a small effect on the field strength for all radii, leaving the field strength at the axis of the tube unchanged, as displayed in Fig. 2b. However, the more strongly twisted case has a slightly larger field strength as we move away from the axis. The buoyancy profile, on the other hand, is influenced by the value of α as increasing the amount of twist increases the inward acting magnetic tension force more than the outward acting magnetic pressure force, therefore altering the Lorentz force and in turn the density deficit. This results in the higher α cases having a smaller inwardly acting gas pressure gradient, and hence being less buoyant at the centre of the flux tube, as demonstrated in Fig. 3b. At the outside of the flux tube the larger field strength of the higher α case causes the tubes to be more buoyant.

3. General dynamics

In this section, we discuss the general dynamics and features of an experiment with axial field strength B₀ = 9 and twist α = 0.4. For a full analysis of the properties of sunspot rotation in this particular experiment, we refer the reader to Paper I. The flux tube begins to rise buoyantly when the experiment begins due to the density deficit introduced. The decreasing temperature profile in the solar interior permits the flux tube to continue to rise until it reaches the photosphere. Due to the photosphere’s stable stratification, the flux tube is no longer buoyant past this point. In order for the flux tube to fully emerge, the initiation of a secondary instability is necessary. If the field is strong enough, as in the case of this experiment, the flux tube then rises into the corona by means of the magnetic buoyancy instability. Further details of this instability can be found in the recent flux emergence review paper, Hood et al. (2012).

Fig. 4 Analysis of sunspot rotation based on an experiment from Paper I (Sturrock et al. 2015). a) The trajectories of two fieldlines relative to the axis (centre of sunspot) with time coloured in the key. b) Coloured contour of vertical vorticity, ωz, at a time soon after the intersecting field becomes vertical.

Once the field emerges and expands into the low density atmosphere, both the magnetic field and plasma exhibit signs of rotation within the two photospheric polarity sources. In this particular experiment, both sunspots undergo rotations through angles of up to 353°. Particular fieldlines are traced in order to follow the trajectories through the photospheric plane as shown in Fig. 4a. Following the approach introduced in Paper I, this figure shows the position of the fieldlines relative to the central axial fieldline and clearly shows the fieldlines follow an almost circular path as they rotate around the sunspot centre. This is also evidenced by vortical motions that develop within both sunspot centres as displayed in Fig. 4b, indicating that the plasma within the sunspots are rotating clockwise. This has important ramifications for both the interior and atmospheric field. We find a clear depletion in magnetic helicity and energy in the interior portion of the volume as the rotational photospheric motions untwist the interior field. In the atmosphere, on the other hand, there is a steady increase in magnetic helicity and energy owing to both the emergence of magnetic flux and rotational motions at the photosphere. In total, 3.6 × 10²³ Wb² of helicity and 8.2 × 10²² J of energy are transported to the atmosphere. We believe that this is related to the propagation of a torsional Alfvén wave at the instance of emergence that untwists the twisted interior flux tube in an effort to equilibrate the twisted interior with the stretched corona (Longcope & Welsch 2000). Now that the basic model and general behaviour of sunspot rotation is established, we consider two parameter groups introduced in Table 1.

4. Varying B₀ with fixed α

In this group, we fix α at 0.3 ($\mathrm{2}{\mathrm{1}}^{}\mathit{/}\mathrm{4}$ turns of twist or equivalently one turn in 3.56 Mm) and vary the field strength at the axis of the tube, B₀, from B₀ = 5 (6500 G) to B₀ = 10 (13 000 G) in steps of 1 (1300 G). This allows us to understand the sole effect that the magnetic field strength has on the rotation of the sunspots. Altering the initial field strength, B₀, changes the evolution of the field in two ways. Firstly, changing B₀ alters the initial density deficit thereby controlling the speed at which the flux tube rises through the solar interior (see Fig. 3a). Secondly, the tube’s evolution is altered on its journey from the photosphere. In order for the flux tube to emerge, the magnetic buoyancy instability must be triggered which occurs when the plasma β is lowered to one. This occurs more quickly for stronger fields due to their higher magnetic pressure. For a weaker field, on the other hand, the magnetic pressure is built up more slowly as the flux tube is squashed and the field is spread at the photosphere. Rotational motions are manifested in several different ways, and hence we investigate a variety of different quantities. Before delving into the rotational properties, we first analyse the general evolution of the magnetic field.

Fig. 5 The evolution of a) the scaled mean of Bz on the z = 0 plane and b) the scaled mean of Bh on the z = 0 plane over rescaled time for the parameters outlined in the legend for Group 1.

4.1. General evolution

To try and understand the influence of the interior magnetic field strength on the evolution of the tube as it rises, we consider how this affects the magnetic field strength at the photosphere. We consider two proxies for the magnetic field strength at the photosphere, namely the mean vertical field strength, ⟨ B_z ⟩ _{z = 0}, and the mean horizontal field strength, ⟨ B_h ⟩ _{z = 0}. Both expressions have been plotted in Figs. 5a and b over rescaled time for all six Group 1 simulations as coloured by the key. We calculate the mean by averaging over the photospheric region where B_z> 3/4max(B_z). Hence, we assume that contributions from the main sunspot field are included and weaker undular field regions outside the spots are excluded. This has been compared with other proxies for the magnetic field, such as the maximum field strength and we find the same general behaviour in this case. In order to take into account the density deficit’s dependence on B₀, we rescale the horizontal axis by redefining time as . This is equivalent to measuring time on an Alfvén timescale rather than a sound timescale. Before we proceed, we note that the symbols on this plot and all subsequent plots do not reflect the spatial or temporal resolution of the experiments. Unless stated otherwise, symbols are only plotted every five grid points (or time units) for visualisation purposes.

All three components of the initial magnetic field are proportional to B₀ and as such we may expect that the field will still be proportional to B₀ when the tube reaches the photosphere. Interestingly, from Fig. 5a, we find that the vertical field strength at the photosphere can instead be scaled by ${\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}$, suggesting that stronger initial fields tend to concentrate and strengthen in the vertical direction at the photosphere. We find that stronger fields emerge more fully with a vertical axis, and hence possess a larger vertical field, B_z, at the photosphere. Flux tubes with weaker fields, on the other hand, tend to spread at the photosphere before the magnetic buoyancy instability is initiated. This could be responsible for a smaller than expected B_z at the photosphere for lower B₀ values.

For completeness, we have also included the horizontal field strength, ${\mathit{B}}_{\mathrm{h}}\mathrm{=}\sqrt{{\mathit{B}}_{\mathit{x}}^{\mathrm{2}}\mathrm{+}{\mathit{B}}_{\mathit{y}}^{\mathrm{2}}}$ and find that it is proportional to B₀, as we predict. As we are averaging over the region where B_z> 3/4max(B_z), we do not see the effects of the horizontal expansion of the field for weaker B₀ values. It is important that we bear these scalings in mind when analysing later results as the magnetic field is altered on its journey to the photosphere, and hence we may not find the scalings we expect. Observations often only consider the line of sight magnetic field, the vertical magnetic field in our case, and so caution must be taken when making deductions about the interior magnetic field.

Fig. 6 The evolution of a) the height of the leading edge of the flux system and b) the length of the axis fieldline in the atmosphere above z = 0 for the varying B0 over rescaled time.

In order to investigate the magnetic field’s journey to the atmosphere, we also consider the leading edge of the flux system over rescaled time for the Group 1 cases in Fig. 6a. We find the evolution to be self-similar until approximately . The B₀ = 5 and B₀ = 6 cases appear to plateau at a fixed height, not reaching the top of the box. The height that the flux tubes reach is determined by pressure balance on the boundary, i.e. where the total pressure within the tube equals the background gas pressure. Emergence slows for the weaker B₀ experiments and there is not enough magnetic pressure to push the boundary upward. Hence, the maximum height reached is lower for weaker field experiments. It is hard to determine where stronger fields reach their pressure balance boundary as they reach the top of the box during the experiment. Consequently, the fieldlines in the stronger B₀ experiments extend further into the atmosphere and so their axis is longer as shown in Fig. 6b. This plot shows the length of the axis fieldline as measured from the point the fieldline enters the z = 0 plane to the point where it leaves through z = 0 after passing through the atmosphere. The axis of the B₀ = 5 tube appears to plateau at a fixed length, much shorter than that of the stronger B₀ tubes given that they extend higher into the corona. The length of the axis fieldline has an approximately linear dependence on the initial field strength, B₀. We return to this concept later as it has important consequences when calculating the twist.

4.2. Torque

In Paper I, a full analysis was performed of the unbalanced torque caused by magnetic forces. This was concluded to be the driver of the rotational motion at the photosphere. The torque, τ_F, is the tendency of a force to rotate an object about an axis and is given by r × F where r is the displacement vector from the centre and F is the force of interest. If we consider a closed circular curve surrounding the maximum of B_z on the photospheric plane, and calculate a surface integral of the torque within this integral, we find the magnetic tension force to be the only contributor. Explicitly, we find the torque due to magnetic forces, τ_F, to be equal to the torque due to magnetic tension, τ_t, as follows $\begin{array}{ccc}{\mathit{\tau }}_{\mathrm{F}}\mathrm{=}{\mathit{\tau }}_{\mathit{t}}\mathrm{=}”{r}\mathrm{\times }\left(\mathrm{(}{B}\mathrm{·}\mathrm{\nabla }\mathrm{)}{B}\right)\mathrm{·}\mathrm{d}{S}\mathit{,}& & \end{array}$where S is the surface contained within a circular contour of radius a = 2.5 surrounding the maximum of B_z. Hence, we speculate that the unbalanced torque produced by the magnetic tension force drives the rotation. This has been plotted in Fig. 7 for all of the Group 1 cases where we have rescaled both quantities as follows. We have again rescaled time as and we also redefine the torque integral by scaling it with respect to ${\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}$ as $\overline{{\mathit{\tau }}_{\mathrm{T}}}\mathrm{=}{\mathit{\tau }}_{\mathrm{T}}\mathit{/}{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}$, given the magnetic tension force is proportional to ${\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}$.

Fig. 7 Torque integrals due to magnetic tension for various B0 cases (as defined in the key) with specifically the rescaled torque integrals, $\overline{{\mathit{\tau }}_{\mathit{T}}}\mathrm{=}{\mathit{\tau }}_{\mathit{T}}\mathit{/}{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}$, measured over rescaled time, .

As stated earlier, all experiments have been executed for different final times, in order to ensure they have the same final . This should ensure that we capture the same period of evolution for all of the varying B₀ experiments and that we are not missing some phases of the evolution for lower B₀ values. Generally, it appears that all six simulations demonstrate a self-similar torque imbalance which drives the rotation. The largest magnitude of torque is found between and , after which the torque damps suggesting a slowing of the rotation after this point.

4.3. Rotation angle

As discussed in Sect. 3, both sunspots experienced significant rotations in the general experiment. This is true for all B₀ cases to varying degrees. In order to calculate the rotation angle, we trace the photospheric location of a series of fieldlines using a fourth-order Runge-Kutta method from the base of the computational domain. In particular, we trace the axis fieldline from the centre of the negative flux source at (0,−15,−25) and, as we expect, it follows the centre of the sunspot. In order to calculate the rotation angle, we also trace a selection of fieldlines from the base within a radius of one around the axis. Given the x and y coordinates of the intersections of selected fieldlines through the photosphere, we can calculate the angle of rotation using $\begin{array}{ccc}\mathit{\phi }\mathrm{=}{\tan }^{-1}\left(\frac{\mathit{y}\mathrm{-}{\mathit{y}}_{\mathrm{axis}}}{\mathit{x}\mathrm{-}{\mathit{x}}_{\mathrm{axis}}}\right)\mathit{,}& & \end{array}$where (x_axis,y_axis) is the location of the axis and (x,y) is the location of another fieldline we have traced. To calculate the angle, we trace 100 fieldlines from a circular footpoint of radius one on the base. We can then calculate the mean rotation angle by averaging the rotation angle over the traced fieldlines within the footpoint of interest. However, as all traced fieldlines intersect at different locations on the photosphere, and hence have different initial rotation angles, we must first subtract off the initial rotation angle, and therefore the average begins at φ = 0 for all cases. The resulting mean rotation angles are shown in Fig. 8a for the six cases we are investigating.

Fig. 8 Rotation angles for various B0 cases with a) the unscaled rotation angles measured over time; b) the rescaled rotation angles, , measured over rescaled time, and; c) the unscaled rotation rate, dφ/ dt measured over rescaled time.

After the emergence of the fields at the photosphere, there is a short period with little change in rotation angle while the sunspots drift apart. From Eq. (2), the buoyancy force is proportional to ${\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}$ and so flux tubes with larger B₀ values appear at the photosphere first. Consequently, the time taken for the flux tubes to reach the photosphere is inversely proportional to B₀. To incorporate this, we again redefine time as and rescale the horizontal axis as shown in Fig. 8b. We also notice a direct relationship between φ and B₀ so we redefine and find that the scaled rotation angles are approximately similar to each other. The scaled time evolution of the scaled mean rotation angles is shown in Fig. 8b.

In Table 2, we have selected the magnitude of the final angles of rotation for the various B₀ cases. The second column of the table contains the unscaled time, t, the third column the rescaled time, , the fourth the unscaled rotation angle, φ, and the fifth the rescaled rotation angle, φ/B₀ to take into account the rotation angle’s dependence on the magnetic field strength. We have chosen to consider the rotation angles at the final rescaled time of as we expect the flux tubes to be in a similar stage of their evolution here. This is presented in the fifth column of the table and shows that the scaled rotation angles are approximately constant. A similar analysis has been conducted for the other sunspot with comparable results.

By rescaling φ with respect to B₀, the only varied quantity in this model, we are able to remove any dependency on B₀ and reveal a self-similar behaviour as the fieldlines threading the sunspot rotate around the centre fieldline. On first inspection, this result may seem surprising as initially all fieldlines have the same helical structure since the degree of twist, α, is constant in this group. Hence, we may have expected all experiments to have the same final rotation angle. This suggests that varying the field strength not only affects the timing at which key processes occur but also the amount by which the fieldlines rotate. In particular, we believe that the length of the fieldlines extending into the atmosphere, as ordered by B₀, is influencing the rotation angle at the photosphere.

Furthermore, the rotation rate, dφ/ dt, as displayed in Fig. 8c, drops off towards the end of all experiments, suggesting that the rotation may not significantly persist if the experiments were continued. By demonstrating the rotation angle’s dependence on the field strength, this corroborates the theory introduced by Min & Chae (2009) that the rotation is a consequence of the torque on the photospheric boundary rather than by apparent effects. If the rotation was due to apparent effects, altering the field strength of the flux tube would not vary the amount of rotation. The relationship between the rotation angle and B₀, and the ramifications of this are investigated in further detail in later sections.

4.4. Twist

In order to estimate the twist of the magnetic field, we investigate a number of twist related quantities. To begin, we calculate the twist of individual fieldlines. Precisely, we consider the fieldline twist of 100 fieldlines stemming from a footpoint of radius one surrounding the axis of the sunspot. This is the same approach we used when calculating the rotation angle. To determine the twist of a fieldline we calculate $\mathrm{\Psi }\mathrm{=}\mathrm{\int }\mathrm{d}\mathit{\psi }\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{L}}\frac{{\mathit{B}}_{\mathit{\psi }}}{\mathit{R}{\mathit{B}}_{\mathit{n}}}\mathrm{d}\mathit{l,}$(3)along the length of a fieldline, L, in a local cylindrical coordinate system (R,ψ,n). Consider a distance of one unit along the axis, henceforth referred to as an axial unit length. Then, B_ψ/ (RB_n) is the angle through which the fieldlines rotate over one axial unit length. Hence, by summing this quantity over the length of the axis fieldline, we calculate the angle in radians through which fieldlines rotate over the axial length. We can then deduce the number of turns, N, the fieldlines pass through over the axial length by noting Ψ = 2πN.

To transform the coordinate system from the original Cartesian system to this local cylindrical system, we trace the axial fieldline using a fourth-order Runge-Kutta method and define a plane at each step along the axis fieldline where the normal to the plane, , is directed along the axis. Next, we create a local cylindrical system within each of these planes by defining R and ψ as the respective radial and azimuthal directions within the plane. Next, the outer fieldlines are traced and the angle through which they rotate is calculated by summing over the angle found when they intersect each of the planes.

Using the method described above, the fieldline twist, and hence the number of turns the fieldlines pass through, is calculated for 100 fieldlines originating from the left footpoint within a radius of one. To gain a mean value, we average over the number of turns the fieldlines pass through for all fieldlines traced. This is plotted against scaled time in Fig. 9. More precisely, we have plotted the average number of turns fieldlines pass through within an interior section of a leg of the flux tube, i.e. up until z = 0, by averaging ${\mathit{N}}_{\mathrm{I}}\mathrm{=}\frac{{\mathrm{\Psi }}_{\mathrm{I}}}{\mathrm{2}\mathit{\pi }}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}{\mathrm{\int }}_{\mathit{y}\hspace{0.17em}\mathit{<}\hspace{0.17em}\mathrm{0}\mathit{,}\mathit{z}\hspace{0.17em}\mathit{<}\hspace{0.17em}\mathrm{0}}\frac{{\mathit{B}}_{\mathit{\psi }}}{\mathit{R}{\mathit{B}}_{\mathit{n}}}\mathrm{d}\mathit{l,}$(4)over 100 fieldlines.

Fig. 9 Average number of turns, NI, fieldlines undergo within the interior portion of one leg of the tube, measured over rescaled time, for Group 1 cases.

Notice, we are again measuring this quantity over rescaled time to take into account the time lag associated with weaker initial field strengths. Initially all fields have the same number of turns around the axis within the interior as they contain the same initial twist. The initial evolution is similar for all cases but the twist drops off more sharply for higher B₀ cases. This result may seem surprising due to the initial twist profile. However, this discrepancy is likely to be related to the expansion of the field in the atmosphere. The stronger B₀ cases expand higher into the atmosphere distributing the atmospheric twist along a longer length (see Fig. 6b) resulting in a smaller twist per unit length in the atmospheric portion of the field. This produces a larger gradient in the twist per unit length, driving the rotation and untwisting the interior field. Notice, there are large amounts of residual twist in the submerged legs of the tube for the weaker field cases. This is related to the distribution of the twist per unit length across the domain.

Unfortunately, we cannot accurately calculate the fieldline twist within the atmospheric section of the tube as the axis kinks from the transverse x direction and the assumptions we require to change into our local planar coordinate system are no longer valid. Hence, this has been excluded from our discussion of the fieldline twist.

In an attempt to explain the large amounts of twist left in the interior for lower B₀ experiments, we present an analysis of a proxy for the local twist referred to as the force-free parameter, defined as $\begin{array}{ccc}{\mathit{\alpha }}_{\mathrm{L}}\mathrm{=}\mathrm{(}\mathrm{\nabla }\mathrm{\times }{B}\mathrm{)}\mathrm{·}{B}\mathit{/}{\mathit{B}}^{\mathrm{2}}\mathit{.}& & \end{array}$This quantity has been shown to be closely linked to the twist per unit length and has been calculated as such in many previous studies, for example Fan (2009) and Sturrock et al. (2015). The parameter is also used as a measure of twist in observational studies (see Liu et al. 2014 or Hahn et al. 2005). It should be noted that this is not the α we vary in Group 2 and we have denoted this as α_L to differentiate between the two. If we assume α_L to be constant and assume the field is force-free, α_L can be shown to be equal to twice the twist per unit length. In order to visualise the distribution of α_L along different sections of the tube, we have traced α_L along the axis passing through the centre of the sunspot from the left footpoint to the apex of the fieldline. To try and compare the different B₀ cases, we have considered a snapshot of α_L along the axial fieldline at two scaled times, and as shown in Fig. 10. Notice, we have plotted α_L against the height of the tube. Hence, the stronger B₀ tubes have reached further into the atmosphere by the final snapshot at . In addition, it should be noted that we have plotted the symbols much less frequently due to the large number of steps we have taken along the axis fieldline.

Fig. 10 The quantity αL traced along the axis fieldline against height, z, of the axis for two scaled times, a) and b). The dashed black line denotes the height at the solar surface, z = 0.

Although not presented in this paper, at t = 0, α_L is constant along the axis fieldline at a value of 0.6 (twice the initial twist per unit length). However, this value drops as the field begins to expand into the atmosphere as shown in Fig. 10a at soon after the field has emerged. We see that the value of α_L drops off with increasing height for all experiments, and hence a gradient develops in α_L. Longcope & Welsch (2000) predicted that this gradient in α_L produces a torque (as we found earlier) that drives the torsional motion of the flux tube intersecting the photosphere. The authors hypothesised that the torsional motion will continue until this gradient in α_L is removed.

At the final scaled time, , we find that the axis of the flux tubes have reached higher into the atmosphere. The magnitude of α_L appears approximately constant in the coronal portion of the field, indicating that this section of the field is almost force-free. The low-magnitude in α_L likely arises because of the stretching and expansion of the field (see Fan 2009). Interestingly, the magnitude of α_L in the interior is ordered by B₀ such that stronger magnetic flux tubes process a lower magnitude of α_L and weaker magnetic flux tubes retain a higher magnitude of α_L in this region. As mentioned earlier, it is often conjectured that the cause of the rotation at the photosphere is related to α_L trying to equilibrate between the twisted interior and stretched coronal field. However, this is not yet the case in the weaker field experiments as a higher magnitude α_L persists in the interior at the final time. Furthermore, the strongest experiments (B₀ = 9 and B₀ = 10) appear to have distributed their twist to the stage where α_L is larger in the atmosphere. This suggests that the experiments have passed through the equilibrium state as the tubes have “over-rotated”. The experiments would need to be performed for a longer time to determine whether the twist is equilibrating on a longer timescale or whether the lower strength field cases are unable to unwind their interior twist to match their coronal twist.

It is important to note that although Fig. 10b is plotted at the same scaled time for all experiments, all experiments appear to be at slightly different stages in their evolution. The reasons behind this are two-fold. Firstly, the weaker fields are equilibrating on a shorter distance and secondly, the Alfvén speed is also proportional to $\mathrm{1}\mathit{/}\sqrt{\mathit{\rho }}$. Given the density deficit’s dependence on B₀, the density deficit of the tube is non-linearly related to B₀ with weaker tubes possessing a larger density. This in turn means weaker tubes have a slightly smaller v_A, and as such may be evolving on a slower timescale than predicted by .

Additionally, it is worth noting that as the weaker fields do not extend as far into the atmosphere, the axis fieldline is shorter. If the twist per unit length, α_L, equilibrates on a shorter fieldline, the average value of α_L will be significantly larger in both the interior and atmospheric regions. Explicitly, if we assume that α_L tends to a constant value and that the total twist is conserved, then α_il_i = α_fl_f where α_i is the initial twist per unit length, α_f the final, l_i the initial axial length, and l_f the final. In this case the predicted final twist per unit length is α_f = α_il_i/l_f. Since all experiments have the same α_i and l_i values, an increase in l_f decreases the final twist per unit length, α_f. In the weaker cases, the final rotation angle will be smaller as the tube does not need to unwind as much interior twist to achieve its larger final twist per unit length. This effect is not yet apparent but again may be seen if the experiments had been performed for a longer time.

4.5. Vorticity

As discussed in Sect. 4.3, all six simulations in Group 1 exhibit rotation in their sunspots, quantifiable in terms of an angle. To examine this further, we have also calculated the mean vertical vorticity within each sunspot, as given by $\begin{array}{ccc}\mathrm{⟨}{\mathit{\omega }}_{\mathit{z}}\mathrm{⟩}\mathrm{=}\mathrm{⟨}\mathrm{(}\mathrm{\nabla }\mathrm{\times }{v}{\mathrm{)}}_{\mathit{z}}\mathrm{⟩}\mathrm{=}⟨\frac{\mathrm{d}{\mathit{v}}_{\mathit{y}}}{\mathrm{d}\mathit{x}}\mathrm{-}\frac{\mathrm{d}{\mathit{v}}_{\mathit{x}}}{\mathrm{d}\mathit{y}}⟩\mathit{,}& & \end{array}$where we have averaged over the photospheric region where B_z> 3/4max(B_z). This quantifies the rotation of the plasma within the upper polarity source. This has been plotted for all six Group 1 experiments as shown in Fig. 11a.

The average vertical vorticity is consistently negative for all B₀ cases indicating that the dominant motion within the sunspots is a clockwise rotation, consistent with the theory suggested in Paper I. Precisely, this rotation acts to untwist the interior magnetic field and inject twist into the atmospheric field. A very clear trend develops in that tubes with a stronger initial field strength emerge more quickly and significant vortical motions develop within their sunspot centres. To try and explore the relationship between ⟨ ω_z ⟩ and B₀, we have rescaled ⟨ ω_z ⟩ with respect to ${\mathit{B}}_{\mathrm{0}}^{\mathrm{3}}$ as well as redefining the time as . The rescaled plot is shown in Fig. 11b where again we find a self-similar evolution, apart from during the latter stages of the B₀ = 5 case. The difference in this case likely arises because we are not capturing the correct area within which the vorticity lies. In the weak B₀ = 5 case, the field spreads over the photosphere in order to build up enough field strength to initiate the magnetic buoyancy instability. This is not captured when considering the area where B_z> 3/4 maxB_z. The ${\mathit{B}}_{\mathrm{0}}^{\mathrm{3}}$ scaling is surprising and may be related to how we calculate ⟨ ω_z ⟩. Future studies should repeat this to see if this is a viable trend. In this particular case, the scaling itself is not of great significance. More importantly, we should conclude that stronger B₀ tubes tend to have stronger vortical motions developing in their polarity sources.

Fig. 11 Average vorticity for Group 1 experiments with a) the unscaled average vertical vorticity measure over time and b) the rescaled average vorticity, $\overline{\mathrm{⟨}{\mathit{\omega }}_{\mathit{z}}\mathrm{⟩}}\mathrm{=}\mathrm{⟨}{\mathit{\omega }}_{\mathit{z}}\mathrm{⟩}\mathit{/}{\mathit{B}}_{\mathrm{0}}^{\mathrm{3}}$, measured over rescaled time, .

Fig. 12 Maximum of jz over photospheric z = 0 plane for various B0 cases with a) the unscaled maximum of jz measured over time and b) the scaled maximum of jz, Maxofjz/B0 measured over rescaled time, . Additionally, c) shows the scaled average of jz, ⟨ jz ⟩, measured over rescaled time.

4.6. Current density

Next, we consider another estimate for the twist of the magnetic field by calculating the electric current density, specifically the z-component, $\begin{array}{ccc}{\mathit{j}}_{\mathit{z}}\mathrm{=}\frac{\mathit{\partial }{\mathit{B}}_{\mathit{y}}}{\mathit{\partial x}}\mathrm{-}\frac{\mathit{\partial }{\mathit{B}}_{\mathit{x}}}{\mathit{\partial y}}\mathit{,}& & \end{array}$as this is linked to how twisted the magnetic field is in the photospheric x−y plane. In a similar manner to Paper I, we analyse the vertical current density at two different planes, the photospheric plane, z = 0, and a plane at the centre of the interior domain, z = −12.5. There are several proxies we can consider for measuring the current. In this case, we plot the temporal evolution of the maximum of j_z for the z = 0 plane in Fig. 12a. All cases show an initial peak in the maximum of j_z at the photosphere due to the emergence of the field. The timing of this maximum is clearly dependent on the value of B₀ as the emergence timescale is proportional to B₀. The magnitude of the peak is also dependent on B₀ as we find the peak in the curve to be higher for larger B₀ values. This is as we expected given that j_z is proportional to B₀. Later, all plots show a steady decline due to the expansion of the field into the higher atmosphere. To investigate the self-similarity in this plot, we have rescaled the maximum of j_z with respect to B₀ and also rescaled time to become an Alfvén time as discussed before. The result of the rescaling is shown in Fig. 12b where we find a clear self-similarity in the curves as they lie on top of one another. From Fig. 12c, however, we note that the average of j_z does not directly scale with B₀. Stronger field experiments have larger averaged currents than predicted by the B₀ scaling, perhaps due to the larger rotation angle and vortical motions seen for stronger experiments.

Fig. 13 Maximum of jz over the interior z = −12.5 plane for various B0 cases with a) the unscaled maximum of jz measured over time and b) the scaled maximum of jz, Maxofjz/B0, measured over rescaled time, .

Similarly, we have plotted the spatial maximum of j_z half way down the solar interior at z = −12.5, as displayed in Fig. 13a. In all cases, after an initial decrease in the maximum of j_z, there is a slight increase due to the straightening of the legs of the tube. However, later there is a steady decrease in the maximum of j_z. Again, we find the timing of the maximum to be dependent on the initial field strength B₀. To take this into account, we have again rescaled the current and time to produce Fig. 13b. A clear self-similarity is seen here.

4.7. Magnetic helicity

To investigate the distribution of twist across the domain, we analyse the relative magnetic helicity within different sub-volumes of the domain. This quantitatively describes the degree of twist and shear of magnetic fieldlines. The share of magnetic helicity throughout the domain is affected by a number of factors, but mainly by vertical flows that move twisted magnetic fields into the corona and horizontal flows at the photosphere that shear and twist up magnetic fields (Berger & Field 1984).

The magnetic helicity of the magnetic field B, relative to some potential field B_p, in a volume V is defined as (Berger & Field 1984; Finn & Antonsen 1985) ${\mathit{H}}_{\mathit{r}}\mathrm{=}{\mathrm{\int }}_{\mathit{V}}\mathrm{(}{A}\mathrm{+}{{A}}_{\mathit{p}}\mathrm{)}\mathrm{·}\mathrm{(}{B}\mathrm{-}{{B}}_{\mathrm{p}}\mathrm{)}\mathrm{d}\mathit{V}\mathit{,}$(5)where A is the vector potential of B and A_p is the vector potential of B_p such that B_p is the reference potential field with the same normal flux distribution as B on all bounding surfaces of V. This form of helicity is chosen as it is independent with respect to the choice of gauge A. To calculate this quantity numerically within different sub-volumes, we use a numerical procedure first introduced in Moraitis et al. (2014). For further details of the method used to calculate the magnetic potential and vector potentials see Valori et al. (2012), and for details of how we implement the code in our particular experiments see Paper I.

Fig. 14 Relative magnetic helicity calculated within the atmospheric region z> 0 for various B0 cases. a) shows the unscaled helicity measured over time and b) details the rescaled magnetic helicity, $\overline{{\mathit{H}}_{\mathit{r}}}\mathrm{=}{\mathit{H}}_{\mathit{r}}\mathit{/}{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}$, measured over rescaled time, .

From previous work, we expect a linear increase in magnetic helicity in the atmosphere accompanied by a depletion of magnetic helicity in the interior region. This is a result of the direct emergence of flux and rotations of sunspots at the photosphere that twist and stress the atmospheric field while untwisting the interior portion of the field. Figure 14a considers the temporal evolution of atmospheric helicity for the different cases, and as expected, it increases for all cases. The injection of helicity is clearly ordered by the value of the initial field strength. Both B and B_p are directly proportional to B₀ in the initial set-up and the vector potentials are therefore also proportional to B₀ as evident by the expressions quoted in Paper I. Therefore, the helicity is initially proportional to ${\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}$. The rescaled atmospheric helicity is plotted against rescaled time in Fig. 14b. There still seems to be a larger amount of B₀ scaled helicity for larger B₀ cases. By doubling the initial magnetic field strength of the sub-photospheric tube, B₀, the helicity transported to the atmosphere is increased by over 23 times. As this is a volume integrated quantity, stronger fields occupy a larger portion of the volume, which may explain the larger helicity injection. In addition, the period of strong rotation may inject more helicity than expected by the ${\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}$ scaling.

Fig. 15 Relative magnetic helicity calculated within the interior portion z< 0 for varying B0 cases. a) considers the unscaled helicity and b) the rescaled magnetic helicity, $\overline{{\mathit{H}}_{\mathit{r}}}\mathrm{=}{\mathit{H}}_{\mathit{r}}\mathit{/}{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}$, measured over rescaled time, .

Similarly, the interior helicity is plotted in Fig. 15a, depicting a reduction in helicity in a similar way. The initial interior helictites are clearly ordered by the value of B₀ and can be directly scaled by ${\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}$. We must bear this in mind when interpreting the helicity as magnetic fields with the same amount of twist may have different amounts of helicity as scaled by their initial magnetic field strength. Furthermore, the stronger experiments clearly undergo a sharper drop-off in helicity. Despite this drop-off, the stronger B₀ experiments have considerably larger amounts of interior helicity compared with the weaker B₀ experiments during the latter stages of the experiment. For instance, in the B₀ = 10 experiment, there is 6.23 × 10²³ Wb² of helicity left in the interior. This may seem surprising as helicity is a measure of the twist of fieldlines and the fieldline twist is considerably lower for the B₀ = 10 case by the end of the experiment (see Fig. 9). However, it is important to note that although a non-zero fieldline twist implies a non-zero helicity, the opposite is not necessarily true (MacTaggart 2015). A magnetic field may have a non-zero helicity without containing a large amount of twist. For example, an untwisted magnetic field may not be exactly potential, and as such may process a non-zero relative helicity.

We again rescale the interior helicity by ${\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}$ as we have plotted in Fig. 15b. In this case, the helicity seems to demonstrate self-similar behaviour, unlike the atmospheric helicity. The change in helicity by resistive dissipation is much larger in the solar interior, so the scaling with ${\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}$ shows a better fit. For more details on the dominant contributors to the change in magnetic helicity, including the resistive dissipation term, see Sturrock et al. (2015) or Pariat et al. (2015).

4.8. Magnetic energy

Finally, an analysis of the magnetic energy is presented to assess how much free energy is produced by the rotational motions at the photosphere. To understand the distribution of energy across the domain, we calculate the free magnetic energy above z = 0, as follows, $\begin{array}{ccc}{\mathit{E}}_{\mathrm{free}}\mathrm{=}\mathrm{\int }\frac{{{B}}^{\mathrm{2}}}{\mathrm{2}}\mathrm{d}\mathit{V}\mathrm{-}\mathrm{\int }\frac{{{B}}_{\mathrm{p}}^{\mathrm{2}}}{\mathrm{2}}\mathrm{d}\mathit{V}\mathit{,}& & \end{array}$where V is defined as the volume above z = 0. This is essentially the excess energy stored in the field as we have subtracted off the energy stored in the potential field with the same normal distribution on the photospheric boundary. Both the unscaled and scaled free magnetic energy are plotted in Figs. 16a and b respectively. The self-similar evolution is followed strictly for the first 275 scaled time units. However, later the different B₀ cases deviate from the original trend suggesting a different B₀ dependence. This agrees with the trend we see in the helicity with larger than expected amounts of energy transported to the atmosphere for stronger initial fields. As the diffusion time is independent of B₀, we see greater diffusion for weaker experiments as they have run for a longer unscaled time. This may explain the drop in free magnetic energy in the latter stages of the experiment. The free magnetic energy transported to the atmosphere ranges from 7.2 × 10²¹ J to 6.9 × 10²² J over the range of B₀ values.

Fig. 16 Free magnetic energy stored in the atmosphere for Group 1 with a) the unscaled energy measured over time and b) the rescaled magnetic energy, $\overline{{\mathit{E}}_{\mathrm{free}}}\mathrm{=}{\mathit{E}}_{\mathrm{free}}\mathit{/}{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}$ measured over rescaled time, .

5. Varying α with fixed B₀

In Group 2, we fix B₀ at 7 (an axial field strength of 9100 G) and vary the initial twist of the tube from α = 0.2 (one turn in 5.34 Mm) to α = 0.4 (one turn in 2.67 Mm). Using a similar approach to the last group, we pinpoint the effect that the initial twist, α, has on the rotation of sunspots. We again investigate a variety of different features related to the rotational movements at the photosphere.

When we modify the degree of twist, α, the initial magnetic tension force acting on the tube is affected and this in turn changes the magnetic buoyancy profile of the tube. This is discussed earlier in Sect. 2.1 where Fig. 3b displays a comparison of the density excess for different values of α. This reveals that as the value of α is reduced the axial region becomes more buoyant but the surrounding plasma is less buoyant due to the smaller field strength here (see Fig. 2b). Therefore, we need to consider that the buoyancy profile is non-linearly altered in this group suggesting that α may have a more complex effect on the dynamics of the experiment. The non-linear dependence of the density deficit on α makes it difficult to rescale the time to remove this effect as we did in the Group 1 case. For brevity, we have excluded the methods behind calculating the quantities in this section. Unless stated otherwise, all quantities are calculated as introduced in Group 1.

5.1. Rotation angle

From Fig. 17, it is clear that the rotation angle, calculated in the same way as in Sect. 4.3, has some dependence on the initial twist, α. Changing the initial degree of twist changes the helical structure of fieldlines and means that fieldlines traced from the same location on the base appear at different locations at the photosphere in all three cases. To try and manage this effect, we have artificially moved all the starting angles to 0 and the subsequent evolution has been shifted. The time at which the field reaches the photosphere is also affected by α with the α = 0.2 tube reaching the photosphere first due to its larger axial buoyancy. In Fig. 17, there is clearly some trend in that the sunspot in the higher α experiment undergoes a larger rotation. If we surmise that the rotation of sunspots is due to the propagation of twist, and we expect the rotation to attempt to equilibrate the twist imbalance, we predict the rotation angles to be largest for the most highly twisted cases.

Fig. 17 Rotation angles measured over time for various α cases, as depicted in the key.

Summarised in Table 4, at the end of the section, are the final rotational angles at t = 120. Notice, in this case all experiments have been run for the same amount of time and we do not rescale t. This clearly shows some ordering of the rotation angle with the parameter α. By doubling the initial twist, we more than double the final rotation angle. However, the relationship is not as clear as we find in Group 1 due to the more complicated effect of α on the tube’s evolution.

5.2. Twist

As discussed in the previous section, the rotational motions at the photosphere extract twist from the interior. The number of turns the field takes around the axis within an interior section of the tube, N_I, is presented in Fig. 18 as calculated in Eq. (4). Given that all experiments start with a different initial twist, α, they all contain differing amounts of twist within the interior when they first intersect the photosphere. In addition, weakly twisted tubes reach the photosphere first as they are more buoyant. Due to the larger rotation seen for more highly twisted fields, the interior twist is extracted more efficiently.

Fig. 18 Average number of turns, NI, fieldlines undergo within the interior portion of one leg of the tube, measured over time for Group 2 cases.

5.3. Vorticity

As the flux tubes reach the photosphere, vortical motions develop on the sunspot centres, as shown in Fig. 19. As B₀ is constant throughout this group of simulations, the region where B_z> 3/4max(B_z) is approximately constant due to B_z’s weak dependence on α. There is, however, a clear trend indicating that tubes with a higher initial degree of twist have larger vortical motions developing on their sunspot centres. This is expected as we predict vortical motions at the photosphere untwist interior field in an attempt to equilibrate the twisted interior with the stretched atmospheric field. If the initial flux tube is highly twisted, the fieldlines threading through the sunspot must rotate through a larger angle (see Fig. 17 and Table 4), producing a higher magnitude of vorticity.

Fig. 19 Average vorticity over time for various α cases.

Fig. 20 Relative magnetic helicity calculated above z = 0 in the atmospheric portion of the volume for varying α cases.

5.4. Magnetic helicity

To complete our analysis of the Group 2 simulations, we present an analysis of the magnetic helicity within two distinct sections of the domain, namely the solar interior and atmosphere, as separated by the z = 0 photospheric boundary. The temporal evolution of this quantity, as calculated using Eq. (5), is shown in Figs. 20 and 21 for the atmosphere and interior respectively. As expected, there is a linear increase in the atmospheric helicity in all cases due to the emergence of flux into the atmosphere and the twisting of the atmospheric field caused by photospheric horizontal flows. The degree of twist clearly alters the amount of magnetic helicity transported to the atmosphere. The amount of helicity in the atmosphere tends to saturate much more quickly for the lower twist cases (red and green). The highly twisted (blue) case, on the other hand, continues to increase over the whole experiment. The final helicity values for the three cases are summarised in Table 4 in physical units. There is some evidence that the helicity may be proportional to α². However, there is not a direct relationship given the non-linear dependence of the magnetic field on α, made clear by the initial field, B, outlined in Eq. (1).

Fig. 21 Relative magnetic helicity calculated within the interior portion below z = 0 for varying α cases.

The initial magnetic helicity stored in the interior is clearly altered by the degree of initial twist of the field as evidenced in Fig. 21. Initially, the total magnetic helicity within the volume is given solely by the magnetic helicity stored within the interior as the magnetic flux tube is yet to emerge and enter the atmosphere. The starting interior helicities range from 3.25 × 10²³ Wb² for the α = 0.2 case to 6.5 × 10²³ Wb² for the α = 0.4 case. The rate at which the helicity in the interior decreases also appears to be dependent on the initial twist with a steeper decline for the highest α experiment. This directly links to the more rapid rotation observed in the highly twisted case.

6. Conclusions

In this paper, we presented results from 3D MHD simulations of buoyant twisted toroidal flux tubes as they rise through the solar interior and emerge into the atmosphere. Our primary aim was to investigate the rotation of the photospheric footpoints. We varied the magnitudes of two parameters governing the magnetic structure of the tube, namely the axial magnetic field strength, B₀, and the twist, α. Our focus was to identify the distinct effect of each of these parameters on the rotational motion at the photosphere and the many ramifications of this photospheric velocity.

To investigate this effect, we analysed various quantities relating to the plasma and magnetic field. To directly measure the rotation, we calculated the rotation angle based on the axis of the flux tube as the centre of the sunspot. This allowed us to make a direct comparison of how magnetic field strength and twist affect rotation rates. Similarly, we looked at how the field strength and twist affect the plasma vorticity within the sunspots. In addition, we analysed the twist of individual fieldlines, magnetic energy, and helicity to study the twist and writhe contained within the different sub-volumes of the domain. This allowed us to understand the distribution of twist across the system, and the transport of twist from the interior to atmosphere of the model.

Many interesting relationships were found for Group 1 in which we kept the twist, α, constant and varied the axial magnetic field strength, B₀. This parameter investigation provides us with an insight into how the initial magnetic field strength affects the amount by which the flux tube rotates. Surprisingly, we found the vertical photospheric magnetic field strength to scale with ${\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}$ when we varied the initial axial field strength, B₀, of the sub-photospheric field. All components of the magnetic field were initially proportional to B₀, but by the time the tube reached the photosphere the magnetic field’s magnitude and direction were adapted as governed by the initial B₀. Stronger fields tended to emerge more fully with a vertically directed axis, whereas weaker fields tended to spread horizontally at the photosphere to allow the magnetic buoyancy instability to initiate. In addition, stronger fields extended higher into the atmosphere and possessed a greater axial length. In brief, the magnetic field was altered on its journey to the photosphere in all experiments, and hence the results we found may be surprising and unforeseen.

Another particularly interesting result we found in Group 1 was that the rotation angle is dependent on B₀. The timescale over which the rotation occurred is dependent on B₀ owing to the density deficit’s dependence on B₀. Hence, in a fixed time, a larger rotation angle was passed through by sunspots in higher B₀ experiments. To remove this time dependence, we scaled the time as , but this did not reveal a self-similarity. Instead, we discovered that the rescaled rotation angle, , is self-similar. This result is surprising as we may have expected the final rotation angle to be the same for varying magnetic field strength as the initial fields shared the same twist and helical structure. The basis for this relationship was difficult to ascertain from this result alone, but became clearer when we investigated the twist. It is conceivable that if we performed the experiments for a longer time, the rotation would cease for stronger experiments and continue for weaker experiments until a plateau was reached for all cases. Owing to the diffusion timescale and computational expense, we are currently unable to check this. However, this seems unlikely as the rotation rate dropped off significantly to almost zero by the end of the weaker B₀ experiments. If the rotation rate does cease for the weaker field experiments and there is not any latter difference in rotation angle for later times, it is possible that the magnetic fields in the weaker experiments are unable to extract as much interior twist as the stronger experiments.

The investigations of twist were also in agreement with the rotation angle results. By considering the fieldline twist and the force-free parameter α_L, we found a considerable amount of twist left in the interior for the weaker experiments. Based on an idealised analytic model, Longcope & Welsch (2000) suggested that the photospheric footpoints of an emerging tube will rotate until the twist per unit length equilibrates along the length of the field, with the interior α_L matching the coronal α_L. In our experiments we do see some evidence of α_L heading to a constant, but found a higher interior α_L for weaker experiments and a lower interior α_L for stronger experiments compared to the constant coronal value. In order to test Longcope’s theory we performed a separate idealised experiment (not presented in this paper), composed of a twisted sub-photospheric flux tube connected to a straight coronal field in a simple stratified domain, and allowed it to evolve. We found that a torsional Alfvén wave propagated twist from the interior to the corona, and importantly found that this process continued until the twist per unit length equilibrated. Interestingly, we found that the photospheric footpoints underwent an over-rotation as the rate of twist in the interior dropped below that of the coronal value. This simple experiment helps us to understand how the flux tubes in our parameter study would have behaved if we had been able to run the experiments for a longer time. Based on this, we expect the sunspots in stronger experiments to rotate in the opposite sense and to return to a constant twist rate, and expect the sunspots in weaker experiment to continue to rotate. However, as the axis fieldline is shorter for weaker fields, the final twist per unit length is larger so it is conceivable that even once α_L has equilibrated, there may still be a considerable amount of twist left in the interior. This could help us explain why the rotation angles are smaller for weaker experiments. In addition, we found the helicity and magnetic energy to be ordered by B₀. In stronger field experiments, we noticed a larger transport of magnetic energy and helicity to the atmosphere.

Varying the initial degree of twist also had an effect on the amount of rotation we calculated within the sunspots. However, it was difficult to find direct relationships with α, given the non-linear dependence of the initial field on α and that the magnetic buoyancy profile is also altered by the degree of twist. We found the rotation angle, twist, helicity, and vorticity to be ordered by the degree of twist, α. Larger vortical motions developed in the highly twisted experiment, transporting more helicity into the atmosphere. Although not included in this paper, we also found a larger increase in free magnetic energy in the atmosphere for the highly twisted field as depicted in Table 4. Work must be done to understand the non-linear effect of the twist on the evolution of the tube. However as the twist and magnetic tension force are inherently linked non-linearly it is difficult to scale quantities in a simple linear manner.

As mentioned in Paper I, it should be noted that we analysed particularly small active regions in our experiments and, hence, the rotation angles we observed occur over much shorter timescales than those found in observations. If we scaled up our experiments, we expect the timescales of rotation, and hence the rotation rate, to be comparable with observations. However, as the size of active region was kept constant across all experiments, this did not impact on this parameter study.

Varying the magnetic field strength and twist of the interior flux tube has a profound effect on the evolution of the flux tube in our experiments as well as the rotational properties at the photosphere. Although this paper has provided insight into this effect, further questions have been raised. Are the trends we find in our numerical parameter investigation seen in observations? For instance, are weaker fields unable to release the twist stored within the interior? If so, what are the reasons behind this? We predict that sunspots stop rotating once the interior α_L balances the coronal α_L (Longcope & Welsch 2000) and believe it is the height of the axis and in turn length of the fieldlines that hinder the transport of twist in weaker magnetic fields. We would like to run experiments for longer in order to study the final rotation angle and distribution of twist across the system. There is much work to be done here from both an observational and modelling point of view. Furthermore, we assumed that the atmosphere is unmagnetised in our experiments but this is certainly not the case within the Sun’s corona. It would be interesting to investigate whether the addition of a magnetised corona affects the rate and amount of rotation within these experiments. The insights found in this paper should be tested and compared with future observations as a model to predict how both the initial magnetic field strength and twist of a sub-photospheric flux tube affect the level of rotation within sunspots.

Acknowledgments

Z.S. acknowledges the financial support of the Carnegie Trust for Scotland. This work used the DIRAC 1, UKMHD Consortium machine at the University of St Andrews and the DiRAC Data Centric system at Durham University, operated by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). This equipment was funded by BIS National E-infrastructure capital grant ST/K00042X/1, STFC capital grant ST/H008519/1, and STFC DiRAC Operations grant ST/K003267/1 and Durham University. DiRAC is part of the National E-Infrastructure. The authors thank the unknown referee for their helpful comments.

References

All Tables

All Figures

	Fig. 1 Summary of initial setup. The background density distribution is shown on the right wall, the temperature distribution on the back wall, a selection of fieldlines are shown in red, and an isosurface of the magnetic field (| B | = 1) is over-plotted. The z = 0 plane is also highlighted in grey as the solar surface.
In the text

	Fig. 2 Radial distribution of the initial magnetic field strength, | B |, at z = −25 for varying a) B0 and b) α.
In the text

	Fig. 3 Radial distribution of the initial density excess, ρexc, at z = −25 for varying a) B0 and b) α.
In the text

	Fig. 4 Analysis of sunspot rotation based on an experiment from Paper I (Sturrock et al. 2015). a) The trajectories of two fieldlines relative to the axis (centre of sunspot) with time coloured in the key. b) Coloured contour of vertical vorticity, ωz, at a time soon after the intersecting field becomes vertical.
In the text

	Fig. 5 The evolution of a) the scaled mean of Bz on the z = 0 plane and b) the scaled mean of Bh on the z = 0 plane over rescaled time for the parameters outlined in the legend for Group 1.
In the text

	Fig. 6 The evolution of a) the height of the leading edge of the flux system and b) the length of the axis fieldline in the atmosphere above z = 0 for the varying B0 over rescaled time.
In the text

	Fig. 7 Torque integrals due to magnetic tension for various B0 cases (as defined in the key) with specifically the rescaled torque integrals, $\overline{{\mathit{\tau }}_{\mathit{T}}}\mathrm{=}{\mathit{\tau }}_{\mathit{T}}\mathit{/}{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}$, measured over rescaled time, .
In the text

	Fig. 8 Rotation angles for various B0 cases with a) the unscaled rotation angles measured over time; b) the rescaled rotation angles, , measured over rescaled time, and; c) the unscaled rotation rate, dφ/ dt measured over rescaled time.
In the text

	Fig. 9 Average number of turns, NI, fieldlines undergo within the interior portion of one leg of the tube, measured over rescaled time, for Group 1 cases.
In the text

	Fig. 10 The quantity αL traced along the axis fieldline against height, z, of the axis for two scaled times, a) and b). The dashed black line denotes the height at the solar surface, z = 0.
In the text

	Fig. 11 Average vorticity for Group 1 experiments with a) the unscaled average vertical vorticity measure over time and b) the rescaled average vorticity, $\overline{\mathrm{⟨}{\mathit{\omega }}_{\mathit{z}}\mathrm{⟩}}\mathrm{=}\mathrm{⟨}{\mathit{\omega }}_{\mathit{z}}\mathrm{⟩}\mathit{/}{\mathit{B}}_{\mathrm{0}}^{\mathrm{3}}$, measured over rescaled time, .
In the text

	Fig. 12 Maximum of jz over photospheric z = 0 plane for various B0 cases with a) the unscaled maximum of jz measured over time and b) the scaled maximum of jz, Maxofjz/B0 measured over rescaled time, . Additionally, c) shows the scaled average of jz, ⟨ jz ⟩, measured over rescaled time.
In the text

	Fig. 13 Maximum of jz over the interior z = −12.5 plane for various B0 cases with a) the unscaled maximum of jz measured over time and b) the scaled maximum of jz, Maxofjz/B0, measured over rescaled time, .
In the text

	Fig. 14 Relative magnetic helicity calculated within the atmospheric region z> 0 for various B0 cases. a) shows the unscaled helicity measured over time and b) details the rescaled magnetic helicity, $\overline{{\mathit{H}}_{\mathit{r}}}\mathrm{=}{\mathit{H}}_{\mathit{r}}\mathit{/}{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}$, measured over rescaled time, .
In the text

	Fig. 15 Relative magnetic helicity calculated within the interior portion z< 0 for varying B0 cases. a) considers the unscaled helicity and b) the rescaled magnetic helicity, $\overline{{\mathit{H}}_{\mathit{r}}}\mathrm{=}{\mathit{H}}_{\mathit{r}}\mathit{/}{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}$, measured over rescaled time, .
In the text

	Fig. 16 Free magnetic energy stored in the atmosphere for Group 1 with a) the unscaled energy measured over time and b) the rescaled magnetic energy, $\overline{{\mathit{E}}_{\mathrm{free}}}\mathrm{=}{\mathit{E}}_{\mathrm{free}}\mathit{/}{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}$ measured over rescaled time, .
In the text

	Fig. 17 Rotation angles measured over time for various α cases, as depicted in the key.
In the text

	Fig. 18 Average number of turns, NI, fieldlines undergo within the interior portion of one leg of the tube, measured over time for Group 2 cases.
In the text

	Fig. 19 Average vorticity over time for various α cases.
In the text

	Fig. 20 Relative magnetic helicity calculated above z = 0 in the atmospheric portion of the volume for varying α cases.
In the text

	Fig. 21 Relative magnetic helicity calculated within the interior portion below z = 0 for varying α cases.
In the text