Sunspot rotation
II. Effects of varying the field strength and twist of an emerging flux tube
School of Mathematics and Statistics, University of St Andrews, St Andrews, Fife, KY16 9SS, UK
email: zoe.sturrock@standrews.ac.uk
Received: 21 February 2016
Accepted: 12 May 2016
Context. Observations of flux emergence indicate that rotational velocities may develop within sunspots. However, the dependence of this rotation on subphotospheric field strength and twist remains largely unknown.
Aims. We investigate the effects of varying the initial field strength and twist of an emerging subphotospheric magnetic flux tube on the rotation of the sunspots at the photosphere.
Methods. We consider a simple model of a stratified domain with a subphotospheric interior layer and three overlying atmospheric layers. A twisted arched flux tube is inserted in the interior and is allowed to rise into the atmosphere. To achieve this, the magnetohydrodynamic equations are solved using the Lagrangianremap code, Lare3d. We perform a parameter study by independently varying the subphotospheric magnetic field strength and twist.
Results. Altering the initial magnetic field strength and twist of the flux tube significantly affects the tube’s evolution and the rotational motions that develop at the photosphere. The rotation angle, vorticity, and current show a direct dependence on the initial field strength. We find that an increase in field strength increases the angle through which the fieldlines rotate, the length of the fieldlines extending into the atmosphere, and the magnetic energy transported to the atmosphere. This also affects the amount of residual twist in the interior. The length of the fieldlines is crucial as we predict the twist per unit length equilibrates to a lower value on longer fieldlines. No such direct dependence is found when we modify the twist of the magnetic field owing to the complex effect this has on the tension force acting on the tube. However, there is still a clear ordering in quantities such as the rotation angle, helicity, and free energy with higher initial twist cases being related to sunspots that rotate more rapidly, transporting more helicity and magnetic energy to the atmosphere.
Key words: magnetohydrodynamics / methods: numerical / Sun: magnetic fields / Sun: interior / Sun: atmosphere
© 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 magnetoconvective 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 subphotospheric 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 selfsimilar evolution in the rise and emergence of the tube when the magnetic field strength is varied. However, they did not find such a simple selfsimilar evolution when varying the twist due to the nonlinear 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 halftorus, and hence a tube rooted deeper in the interior are twofold: 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 subphotospheric 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 subphotospheric 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 nonlinear 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 subphotospheric flux tube, namely the magnetic field strength at the axis, B_{0}, 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 setup 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 subphotospheric 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 (1)where B_{φ} = B_{0}e^{− r2/a2} and B_{θ} = αrB_{φ}. B_{0} is the axial field strength, α is the degree of twist of the fieldlines, a is the minor radius of the flux tube, and R_{0} 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 (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 . 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 overplotted. The z = 0 plane is also highlighted in grey as the solar surface. 

Open with DEXTER 
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_{0} = 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_{0}, and the degree of twist, α. The twist is assumed to be positive and constant in all experiments, ensuring that all tubes are uniformly righthand twisted. Each fieldline rotates about the axis through an angle of α radians over one unit of distance along the axis. The initial setup of a representative experiment is summarised in Fig. 1, with parameters B_{0} = 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_{0} is varied and Group 2 where B_{0} is fixed and α is varied. A summary of the B_{0} and α values under consideration is given in Table 1. It should be noted we only consider a relatively small range of B_{0} values as a consequence of our model choice. As found in previous parameter studies (see MacTaggart & Hood 2009), if we pick a lower B_{0} value, the flux tube will fail to fully emerge, and if we choose a much higher B_{0} value, we may encounter unphysical negative pressures. In fact, if we increase the initial axial field strength to B_{0} = 12, we encounter negative pressure in the particular initial setup outlined in the previous section. It is, however, important to note that the range of allowable field strengths will vary from setup to setup. 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^{19} Mx to 7.4 × 10^{19} Mx. This is typical of a small active region or large ephemeral region.
Parameter space under investigation.
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_{0} = 5 for 216 time units (90 min); B_{0} = 6 for 180 time units (75 min); B_{0} = 7 for 154 time units (64 min); B_{0} = 8 for 135 time units (56 min); B_{0} = 9 for 120 time units (50 min); and B_{0} = 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) B_{0} and b) α. 

Open with DEXTER 
Fig. 3 Radial distribution of the initial density excess, ρ_{exc}, at z = −25 for varying a) B_{0} and b) α. 

Open with DEXTER 
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_{0} 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_{0}. The buoyancy profile, too, is strongly dependent on B_{0} because the density deficit is proportional to . The B_{0} = 10 tube will therefore be 4 times more buoyant than the B_{0} = 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_{0} = 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. 

Open with DEXTER 
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^{23} Wb^{2} of helicity and 8.2 × 10^{22} 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_{0} with fixed α
In this group, we fix α at 0.3 ( turns of twist or equivalently one turn in 3.56 Mm) and vary the field strength at the axis of the tube, B_{0}, from B_{0} = 5 (6500 G) to B_{0} = 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_{0}, changes the evolution of the field in two ways. Firstly, changing B_{0} 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 B_{z} on the z = 0 plane and b) the scaled mean of B_{h} on the z = 0 plane over rescaled time for the parameters outlined in the legend for Group 1. 

Open with DEXTER 
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_{0}, 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_{0} and as such we may expect that the field will still be proportional to B_{0} 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 , 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_{0} values.
For completeness, we have also included the horizontal field strength, and find that it is proportional to B_{0}, 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_{0} 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 B_{0} over rescaled time. 

Open with DEXTER 
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 selfsimilar until approximately . The B_{0} = 5 and B_{0} = 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_{0} 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_{0} 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_{0} = 5 tube appears to plateau at a fixed length, much shorter than that of the stronger B_{0} 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_{0}. 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 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 as , given the magnetic tension force is proportional to .
Fig. 7 Torque integrals due to magnetic tension for various B_{0} cases (as defined in the key) with specifically the rescaled torque integrals, , measured over rescaled time, . 

Open with DEXTER 
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_{0} experiments and that we are not missing some phases of the evolution for lower B_{0} values. Generally, it appears that all six simulations demonstrate a selfsimilar 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_{0} cases to varying degrees. In order to calculate the rotation angle, we trace the photospheric location of a series of fieldlines using a fourthorder RungeKutta 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 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 B_{0} 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. 

Open with DEXTER 
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 and so flux tubes with larger B_{0} values appear at the photosphere first. Consequently, the time taken for the flux tubes to reach the photosphere is inversely proportional to B_{0}. 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_{0} 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.
Rotation angles for Group 1 experiments.
In Table 2, we have selected the magnitude of the final angles of rotation for the various B_{0} 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_{0} 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_{0}, the only varied quantity in this model, we are able to remove any dependency on B_{0} and reveal a selfsimilar 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_{0}, 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_{0}, 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 (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 fourthorder RungeKutta 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 (4)over 100 fieldlines.
Fig. 9 Average number of turns, N_{I}, fieldlines undergo within the interior portion of one leg of the tube, measured over rescaled time, for Group 1 cases. 

Open with DEXTER 
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_{0} 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_{0} 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_{0} experiments, we present an analysis of a proxy for the local twist referred to as the forcefree parameter, defined as 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 forcefree, α_{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_{0} 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_{0} 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. 

Open with DEXTER 
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 forcefree. The lowmagnitude 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_{0} 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_{0} = 9 and B_{0} = 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 “overrotated”. 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 twofold. Firstly, the weaker fields are equilibrating on a shorter distance and secondly, the Alfvén speed is also proportional to . Given the density deficit’s dependence on B_{0}, the density deficit of the tube is nonlinearly related to B_{0} 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 α_{i}l_{i} = α_{f}l_{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} = α_{i}l_{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 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_{0} 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_{0}, we have rescaled ⟨ ω_{z} ⟩ with respect to as well as redefining the time as . The rescaled plot is shown in Fig. 11b where again we find a selfsimilar evolution, apart from during the latter stages of the B_{0} = 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_{0} = 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 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_{0} 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, , measured over rescaled time, . 

Open with DEXTER 
Fig. 12 Maximum of j_{z} over photospheric z = 0 plane for various B_{0} cases with a) the unscaled maximum of j_{z} measured over time and b) the scaled maximum of j_{z}, Maxofj_{z}/B_{0} measured over rescaled time, . Additionally, c) shows the scaled average of j_{z}, ⟨ j_{z} ⟩, measured over rescaled time. 

Open with DEXTER 
4.6. Current density
Next, we consider another estimate for the twist of the magnetic field by calculating the electric current density, specifically the zcomponent, 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_{0} as the emergence timescale is proportional to B_{0}. The magnitude of the peak is also dependent on B_{0} as we find the peak in the curve to be higher for larger B_{0} values. This is as we expected given that j_{z} is proportional to B_{0}. Later, all plots show a steady decline due to the expansion of the field into the higher atmosphere. To investigate the selfsimilarity in this plot, we have rescaled the maximum of j_{z} with respect to B_{0} 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 selfsimilarity 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_{0}. Stronger field experiments have larger averaged currents than predicted by the B_{0} scaling, perhaps due to the larger rotation angle and vortical motions seen for stronger experiments.
Fig. 13 Maximum of j_{z} over the interior z = −12.5 plane for various B_{0} cases with a) the unscaled maximum of j_{z} measured over time and b) the scaled maximum of j_{z}, Maxofj_{z}/B_{0}, measured over rescaled time, . 

Open with DEXTER 
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_{0}. To take this into account, we have again rescaled the current and time to produce Fig. 13b. A clear selfsimilarity is seen here.
4.7. Magnetic helicity
To investigate the distribution of twist across the domain, we analyse the relative magnetic helicity within different subvolumes 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) (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 subvolumes, 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 B_{0} cases. a) shows the unscaled helicity measured over time and b) details the rescaled magnetic helicity, , measured over rescaled time, . 

Open with DEXTER 
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_{0} in the initial setup and the vector potentials are therefore also proportional to B_{0} as evident by the expressions quoted in Paper I. Therefore, the helicity is initially proportional to . The rescaled atmospheric helicity is plotted against rescaled time in Fig. 14b. There still seems to be a larger amount of B_{0} scaled helicity for larger B_{0} cases. By doubling the initial magnetic field strength of the subphotospheric tube, B_{0}, 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 scaling.
Fig. 15 Relative magnetic helicity calculated within the interior portion z< 0 for varying B_{0} cases. a) considers the unscaled helicity and b) the rescaled magnetic helicity, , measured over rescaled time, . 

Open with DEXTER 
Final atmospheric magnetic helicity for Group 1 experiments.
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_{0} and can be directly scaled by . 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 dropoff in helicity. Despite this dropoff, the stronger B_{0} experiments have considerably larger amounts of interior helicity compared with the weaker B_{0} experiments during the latter stages of the experiment. For instance, in the B_{0} = 10 experiment, there is 6.23 × 10^{23} Wb^{2} 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_{0} = 10 case by the end of the experiment (see Fig. 9). However, it is important to note that although a nonzero fieldline twist implies a nonzero helicity, the opposite is not necessarily true (MacTaggart 2015). A magnetic field may have a nonzero 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 nonzero relative helicity.
We again rescale the interior helicity by as we have plotted in Fig. 15b. In this case, the helicity seems to demonstrate selfsimilar behaviour, unlike the atmospheric helicity. The change in helicity by resistive dissipation is much larger in the solar interior, so the scaling with 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, 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 selfsimilar evolution is followed strictly for the first 275 scaled time units. However, later the different B_{0} cases deviate from the original trend suggesting a different B_{0} 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_{0}, 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^{21} J to 6.9 × 10^{22} J over the range of B_{0} 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, measured over rescaled time, . 

Open with DEXTER 
5. Varying α with fixed B_{0}
In Group 2, we fix B_{0} 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 nonlinearly altered in this group suggesting that α may have a more complex effect on the dynamics of the experiment. The nonlinear 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. 

Open with DEXTER 
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, N_{I}, fieldlines undergo within the interior portion of one leg of the tube, measured over time for Group 2 cases. 

Open with DEXTER 
Summary of key results at t = 120 for Group 2 experiments.
5.3. Vorticity
As the flux tubes reach the photosphere, vortical motions develop on the sunspot centres, as shown in Fig. 19. As B_{0} 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. 

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

Open with DEXTER 
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 α^{2}. However, there is not a direct relationship given the nonlinear 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. 

Open with DEXTER 
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^{23} Wb^{2} for the α = 0.2 case to 6.5 × 10^{23} Wb^{2} 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_{0}, 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 subvolumes 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_{0}. 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 when we varied the initial axial field strength, B_{0}, of the subphotospheric field. All components of the magnetic field were initially proportional to B_{0}, but by the time the tube reached the photosphere the magnetic field’s magnitude and direction were adapted as governed by the initial B_{0}. 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_{0}. The timescale over which the rotation occurred is dependent on B_{0} owing to the density deficit’s dependence on B_{0}. Hence, in a fixed time, a larger rotation angle was passed through by sunspots in higher B_{0} experiments. To remove this time dependence, we scaled the time as , but this did not reveal a selfsimilarity. Instead, we discovered that the rescaled rotation angle, , is selfsimilar. 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_{0} 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 forcefree 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 subphotospheric 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 overrotation 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_{0}. 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 nonlinear 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 nonlinear effect of the twist on the evolution of the tube. However as the twist and magnetic tension force are inherently linked nonlinearly 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 subphotospheric 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 Einfrastructure 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 EInfrastructure. The authors thank the unknown referee for their helpful comments.
References
 Arber, T. D., Longbottom, A. W., Gerrard, C. L., & Milne, A. M. 2001, J. Comput. Phys., 171, 151 [NASA ADS] [CrossRef] [Google Scholar]
 Berger, M. A., & Field, G. B. 1984, J. Fluid Mech., 147, 133 [NASA ADS] [CrossRef] [Google Scholar]
 Brown, D. S., Nightingale, D., Alexander, D., et al. 2003, Sol. Phys., 216, 79 [NASA ADS] [CrossRef] [Google Scholar]
 Evershed, J. 1909, MNRAS, 69, 454 [NASA ADS] [CrossRef] [Google Scholar]
 Fan, Y. 2009, ApJ, 697, 1529 [NASA ADS] [CrossRef] [Google Scholar]
 Finn, J. M., & Antonsen, T. M. 1985, Comments Plasma Phys. Control. Fusion, 111, 1529 [Google Scholar]
 Hahn, M., Gaard, S., Jibben, P., Canfield, R. C., & Nandy, D. 2005, ApJ, 629, 1135 [NASA ADS] [CrossRef] [Google Scholar]
 Hood, A. W., Archontis, V., Galsgaard, K., & MorenoInsertis, F. 2009, A&A, 503, 999 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hood, A. W., Archontis, V., & MacTaggart, D. 2012, Sol. Phys., 278, 3 [NASA ADS] [CrossRef] [Google Scholar]
 Liu, Y., Hoeksema, J. T., & Sun, X. 2014, ApJ, 783, L1 [NASA ADS] [CrossRef] [Google Scholar]
 Longcope, D. W., & Welsch, B. T. 2000, ApJ, 545, 1089 [NASA ADS] [CrossRef] [Google Scholar]
 MacTaggart, D., & Hood, A. W. 2009, A&A, 507, 995 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Min, S., & Chae, J. 2009, Sol. Phys., 258, 203 [NASA ADS] [CrossRef] [Google Scholar]
 Moraitis, K., Tziotziou, K., Georgoulis, M. K., & Archontis, V. 2014, Sol. Phys., 289, 4453 [NASA ADS] [CrossRef] [Google Scholar]
 Murray, M. J., Hood, A. W., MorenoInsertis, F., Galsgaard, K., & Archontis, V. 2006, A&A, 460, 909 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pariat, E., Valori, G., Démoulin, P., & Dalmasse, K. 2015, A&A, 580, A128 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sturrock, Z., Hood, A. W., Archontis, V., & McNeill, C. M. 2015, A&A, 582, A76 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Valori, G., Démoulin, P., & Pariat, E. 2012, Sol. Phys, 278, 347 [NASA ADS] [CrossRef] [Google Scholar]
 Yan, X. L., & Qu, Z. Q. 2007, A&A, 468, 1083 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Yan, X. L., Qu, Z. Q., & Xu, C. L. 2008, ApJ, 682, L65 [NASA ADS] [CrossRef] [Google Scholar]
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 overplotted. The z = 0 plane is also highlighted in grey as the solar surface. 

Open with DEXTER  
In the text 
Fig. 2 Radial distribution of the initial magnetic field strength,  B , at z = −25 for varying a) B_{0} and b) α. 

Open with DEXTER  
In the text 
Fig. 3 Radial distribution of the initial density excess, ρ_{exc}, at z = −25 for varying a) B_{0} and b) α. 

Open with DEXTER  
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. 

Open with DEXTER  
In the text 
Fig. 5 The evolution of a) the scaled mean of B_{z} on the z = 0 plane and b) the scaled mean of B_{h} on the z = 0 plane over rescaled time for the parameters outlined in the legend for Group 1. 

Open with DEXTER  
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 B_{0} over rescaled time. 

Open with DEXTER  
In the text 
Fig. 7 Torque integrals due to magnetic tension for various B_{0} cases (as defined in the key) with specifically the rescaled torque integrals, , measured over rescaled time, . 

Open with DEXTER  
In the text 
Fig. 8 Rotation angles for various B_{0} 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. 

Open with DEXTER  
In the text 
Fig. 9 Average number of turns, N_{I}, fieldlines undergo within the interior portion of one leg of the tube, measured over rescaled time, for Group 1 cases. 

Open with DEXTER  
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. 

Open with DEXTER  
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, , measured over rescaled time, . 

Open with DEXTER  
In the text 
Fig. 12 Maximum of j_{z} over photospheric z = 0 plane for various B_{0} cases with a) the unscaled maximum of j_{z} measured over time and b) the scaled maximum of j_{z}, Maxofj_{z}/B_{0} measured over rescaled time, . Additionally, c) shows the scaled average of j_{z}, ⟨ j_{z} ⟩, measured over rescaled time. 

Open with DEXTER  
In the text 
Fig. 13 Maximum of j_{z} over the interior z = −12.5 plane for various B_{0} cases with a) the unscaled maximum of j_{z} measured over time and b) the scaled maximum of j_{z}, Maxofj_{z}/B_{0}, measured over rescaled time, . 

Open with DEXTER  
In the text 
Fig. 14 Relative magnetic helicity calculated within the atmospheric region z> 0 for various B_{0} cases. a) shows the unscaled helicity measured over time and b) details the rescaled magnetic helicity, , measured over rescaled time, . 

Open with DEXTER  
In the text 
Fig. 15 Relative magnetic helicity calculated within the interior portion z< 0 for varying B_{0} cases. a) considers the unscaled helicity and b) the rescaled magnetic helicity, , measured over rescaled time, . 

Open with DEXTER  
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, measured over rescaled time, . 

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

Open with DEXTER  
In the text 
Fig. 18 Average number of turns, N_{I}, fieldlines undergo within the interior portion of one leg of the tube, measured over time for Group 2 cases. 

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

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

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

Open with DEXTER  
In the text 