Issue 
A&A
Volume 582, October 2015



Article Number  A76  
Number of page(s)  14  
Section  The Sun  
DOI  https://doi.org/10.1051/00046361/201526521  
Published online  12 October 2015 
Sunspot rotation
I. A consequence of flux emergence ^{⋆}
School of Mathematics and Statistics, University of St. Andrews, St. Andrews, Fife, KY16 9SS, UK
email: zoe.sturrock@standrews.ac.uk
Received: 11 May 2015
Accepted: 2 August 2015
Context. Solar eruptions and high flare activity often accompany the rapid rotation of sunspots. The study of sunspot rotation and the mechanisms driving this motion are therefore key to our understanding of how the solar atmosphere attains the conditions necessary for large energy release.
Aims. We aim to demonstrate and investigate the rotation of sunspots in a 3D numerical experiment of the emergence of a magnetic flux tube as it rises through the solar interior and emerges into the atmosphere. Furthermore, we seek to show that the subphotospheric twist stored in the interior is injected into the solar atmosphere by means of a definitive rotation of the sunspots.
Methods. A numerical experiment is performed to solve the 3D resistive magnetohydrodynamic equations using a LagrangianRemap code. We track the emergence of a toroidal flux tube as it rises through the solar interior and emerges into the atmosphere investigating various quantities related to both the magnetic field and plasma.
Results. Through detailed analysis of the numerical experiment, we find clear evidence that the photospheric footprints or sunspots of the flux tube undergo a rotation. Significant vertical vortical motions are found to develop within the two polarity sources after the field emerges. These rotational motions are found to leave the interior portion of the field untwisted and twist up the atmospheric portion of the field. This is shown by our analysis of the relative magnetic helicity as a significant portion of the interior helicity is transported to the atmosphere. In addition, there is a substantial transport of magnetic energy to the atmosphere. Rotation angles are also calculated by tracing selected fieldlines; the fieldlines threading through the sunspot are found to rotate through angles of up to 353° over the course of the experiment. We explain the rotation by an unbalanced torque produced by the magnetic tension force, rather than an apparent effect.
Key words: magnetohydrodynamics (MHD) / Sun: magnetic fields / sunspots / methods: numerical
The movies associated to Figs. 3, 5, and 11 are available in electronic form at http://www.aanda.org
© ESO, 2015
1. Introduction
Sunspot rotation has attracted the interest of many researchers over the years from both an observational and modelling perspective. The wide scope in observations of sunspot rotation merits a study of the mechanisms driving this motion. Coronal mass ejections (CMEs) and solar flares are often related to the rapid rotation of sunspots. Hence, the study of sunspot rotation is crucial to our understanding of such explosive events on the Sun as a possible mechanism for allowing the corona to achieve the necessary conditions for high energy release.
Just over a century ago, whilst working at the Kodaikanal Solar Observatory, Evershed (1909) first discovered evidence of the rotation of sunspots based on spectral observations. Since this initial evidence, sunspot rotations have been analysed in numerous observational studies. In recent years, detailed case studies such as Yan & Qu (2007), Yan et al. (2009), and Min & Chae (2009) and statistical analyses from Brown et al. (2003) and Yan et al. (2008) have investigated this phenomena. These studies, along with several others, have shown that sunspots can exhibit significant rotation of the order of several hundreds of degrees over a few days. Yan et al. (2008) conducted a statistical study of rotating sunspots using Transition Region And Coronal Explorer (TRACE), Hinode, and Michelson Doppler Imager (MDI) magnetograms from 1996 to 2007, in which they individually analysed 2959 active regions. They found 182 significantly rotating sunspots within 153 active regions. This is equivalent to approximately 5% of active regions harbouring rotating sunspots. On the other hand, Brown et al. (2003) conducted a more detailed study of seven sunspots using white light TRACE data to find sunspot centres and track notable features over time to calculate rotation rates. Rotation angles of between 40° and 200° were observed over periods of three to five days, resulting in an average rotation rate of a few degrees per hour. Six of these rotating spots resulted in subsequent flaring activity and the energisation of the corona.
Interestingly, measurements of sunspot rotation have given variable results depending on the methodology employed. For example, Min & Chae (2009) and Yan et al. (2009) analysed the same active region, NOAA 10930, and found notably different results. Min & Chae (2009) noted a counterclockwise rotation of 540 degrees over five days, whereas Yan et al. (2009) focused on the X3.2 flare that followed the rapid rotation of 259 degrees over four days. Furthermore, although Brown et al. (2003) and Yan & Qu (2007) both concluded that different parts of sunspots often rotate at differing speeds, Brown et al. (2003) noted that the highest rotation rate was found in the penumbra while Yan & Qu (2007) concluded that the highest rotation rate was found in the umbra. This suggests that sunspots do not necessarily rotate as a rigid body. Yan & Qu (2007) concluded that twist can be created by a variation in rotation rate with distance from the centre of a sunspot. This twist can then be injected into the corona kinking the magnetic loops and driving flare activity.
Eruptions and flares are, in fact, often correlated with rotational motions of sunspots. Observations have shown direct evidence of the energisation of the corona by these rotations. The initiation of CMEs by sunspot rotation was studied in detail by Török et al. (2014) from both an observational and modelling viewpoint. They found that the rotation of sunspots can significantly weaken the magnetic tension of the overlying field in active regions and can then trigger an eruption in this region. This is an alternative explanation for eruptions caused by rotation of sunspots to the common theory that eruptions are triggered by a direct injection of twist (Yan & Qu 2007).
This paper aims to discuss and simulate a possible mechanism for the rotation of sunspots. Two mechanisms were discussed in Brown et al. (2003), namely, photospheric flows and magnetic flux emergence. Photospheric flows are primarily due to the large scale effect of differential rotation and localised motions resulting from magnetoconvective dynamics. The effects of differential rotation are kept to a minimum in Brown et al. (2003) as the images were derotated prior to measuring the velocities. The second mechanism discussed was the emergence of magnetic flux. Brown et al. (2003) suggested that the photospheric footpoints of a flux tube are observed to rotate as the tube emerges. This is the proposed mechanism for the rotation of sunspots which we will investigate. This mechanism seems the most appropriate according to the findings of Brown et al. (2003) as well as a case study from Min & Chae (2009). The latter study supports flux emergence as a viable mechanism as they discovered that the rotation speed increases in close relation to the growth of the sunspot of interest, which we attribute to the emergence of flux. The next logical question is, how can flux emergence drive sunspot rotation? Two possible explanations were suggested by Min & Chae (2009). They conjectured that the rotation may be an apparent motion caused by a twisted flux tube rising vertically and the fieldlines successively crossing the photospheric boundary at different locations as a result of the twisted structure of the field. In this case there is no real horizontal motion of the plasma; instead this rotation is a virtual effect caused by continued displacement of the footpoints. Alternatively, Min & Chae (2009) proposed that the observed rotation may represent the real horizontal motion of the plasma due to a net torque originating in the interior driving the plasma to rotate on the photospheric boundary. The torque will be examined later in an effort to determine its origin.
The process of magnetic flux emergence has been considered numerically in various experiments in recent years (for instance Fan 2001, Archontis et al. 2004, Hood et al. 2009 amongst many others). The widely accepted picture of sunspot formation is that an Ωshaped flux tube rises from the base of the convection zone until its apex intersects the photosphere to form a pair of sunspots. If a magnetic flux tube is in pressure balance and thermal equilibrium with its surroundings, the tube will be less dense than its surroundings, and will therefore be buoyant. This mechanism allows a flux tube to rise to the stably stratified photosphere where it remains until the tube is able to enter the atmosphere by initiation of a second instability, namely the magnetic buoyancy instability. A pair of concentrations of opposite polarities, known as bipolar sunspots, mark the intersection of the field at the surface. In these experiments, the computational domain typically models a region from the top of the solar interior to the lower corona. A magnetic flux tube is then placed in the solar interior and is made buoyant by either introducing a density deficit or imposing an initial upward velocity. A recent review of simulations relating to this emergence process was undertaken by Hood et al. (2012).
Investigations of magnetic flux emergence as a possible mechanism for sunspot rotation have been conducted in the past. Longcope & Welsch (2000) suggested that the rotation of sunspots is a consequence of the transport of helicity from the convection zone to the corona as a twisted flux tube emerges. They also demonstrated that a torsional Alfvén wave will propagate downwards at the instance of emergence due to a mismatch in current between the highly twisted interior and stretched coronal portion of the field. In addition, Gibson et al. (2004) explained the rotation as an observational consequence of the emergence of a flux tube through the photosphere. An investigation of the transport of magnetic energy and helicity in an emerging flux model was conducted by Magara & Longcope (2003). They found that emergence generates horizontal flows as well as vertical flows, both of which contribute to the injection of energy and helicity to the atmosphere. The contributions by vertical flows are dominant initially but horizontal flows are the primary source at a later stage. A more comprehensive study of the horizontal flows at the photosphere during emergence was later performed by Magara (2006) in which they found that rotational flows formed in each of the polarity concentrations soon after the intersecting flux tubes became vertical. Furthermore, the rotational movement of sunspots in a 3D magnetohydrodynamic (MHD) simulation has been examined by Fan (2009) where significant vortical motions developed as a torsional Alfvén wave propagated along the tube. Fan (2009) noted that the rotation of the two polarities twisted up the inner fieldlines of the emerged field, thereby transporting twist from the interior portion of the flux tube to the stretched coronal loop. In this simulation, a cylindrical flux tube is inserted into the solar interior using the cylindrical model developed by Fan (2001).
Hood et al. (2009) perform a complementary simulation to that of Fan (2001), replacing the initial cylindrical flux tube with a toroidal tube. A common shortcoming of the cylindrical model in simulations without convective flows is that the axis of the tube never fully emerges. Altering the geometry of the flux tube to a curved shape allows for the axis of the tube to rise into the corona. Current theories suggest that the Sun’s magnetic field is created in the solar interior by dynamo action. Hence, a halftorus shaped flux tube imitates a field anchored deeper within the solar interior. The rotation of sunspots has not yet been investigated using this toroidal model; this is what we aim to study.
In the present paper, we perform a resistive 3D MHD simulation of a toroidal flux tube placed in the solar interior and track its emergence through the photosphere and lower atmosphere. We study the role of rotational flows at the photosphere while also investigating the transport of magnetic helicity and energy. In addition, we explicitly calculate the angles of rotation in our experiment to directly compare with observations. Our main aim is to demonstrate that the interior magnetic field is untwisting as the tube emerges causing an injection of twist into the atmosphere as well as a rotation of the sunspots. This will be accomplished by performing a detailed study investigating quantities relating to both the magnetic field and plasma.
The remainder of the paper is structured as follows. In Sect. 2, we describe the MHD equations used in our model as well as outlining the numerical approach employed. We also describe the initial setup of our experiment, detailing both the initial atmosphere and the model of the subphotospheric magnetic field inserted in the solar interior. In Sect. 3, the simulation results are presented, with emphasis on the rotational motions that develop within the two polarities on the photospheric plane. Our analysis also includes an investigation of the sunspot rotation angle, the flow vorticity at the photosphere, the current, relative magnetic helicity, magnetic energy and twist. Finally, in Sect. 4 we summarise the conclusions and discuss the implications of this simulation and future projects that stem from this work.
2. Model setup
In this section, we outline the numerical setup of our experiment, detailing the initial atmospheric and magnetic field configuration.
2.1. Model equations and numerical approach
For the experiment performed in this paper, we numerically solve the 3D timedependent resistive MHD equations, as described below in nondimensionalised form: with specific energy density, (5)and the electric current density, (6)The basic dimensionless quantities used in these equations are the density ρ, time t, velocity v, magnetic field B, pressure p, gravity g, electric field E, specific energy density ϵ, resistivity η and temperature T. The ratio of specific heats, γ, is taken to be 5/3. The viscosity tensor is denoted by T and the contribution of viscous heating to the energy equation is represented by Q_{visc}. The difference between simulations with and without a small viscosity are found to be negligible. The variables are made dimensionless against photospheric values. These values are: pressure, p_{ph} = 1.4 × 10^{4} Pa; density, ρ_{ph} = 3 × 10^{4} kg m^{3}; temperature T_{ph} = 5.6 × 10^{3} K; pressure scale height, H_{ph} = 170 km; velocity, v_{ph} = 6.8 km s^{1}; time, t_{ph} = 25 s and magnetic field, B_{ph} = 1300 Gauss. The dimensionless resistivity η has been taken as uniform set at a value of 0.005 in our experiment. All quantities from now on will be dimensionless, unless units are provided. In order to reach the true quantities with physical units, the dimensionless quantities should be multiplied by their photospheric values as given above.
The code used to simulate the emergence process is a 3D Lagrangian remap code, Lare3d (Arber et al. 2001). The code uses a staggered grid and is second order accurate in both space and time. The LARE code can be divided into two main steps; the Lagrangian step and the remap step. In the Lagrangian step the equations are solved in a frame that moves with the fluid. This causes the grid to be distorted, so in order to put the variables back on to the original (Eulerian) grid, a remap step is used. The code accurately resolves shocks by using a combination of artificial viscosity and Van Leer flux limiters. This code also includes a small shock viscosity to resolve shocks and the associated shock heating term in the energy equation.
The equations are solved on a uniform Cartesian grid (x, y, z) comprised of 512^{3} grid points. The box spans from −50 to 50 in the x and y directions and from −25 to 75 in the z direction in dimensionless units. This corresponds to a physical size of 17 Mm^{3}. The boundary conditions are periodic on the side walls of the computational domain and the top and bottom boundaries are closed with v = 0. For all other variables we set the normal derivatives to zero. As a consequence, the magnetic field is fixed on the bottom boundary.
2.2. Initial atmosphere
The initial stratification of the atmosphere is the same as many previous flux emergence studies, for example Hood et al. (2009). The computational domain is split into four regions: the solar interior (z< 0); the photosphere/chromosphere 0 ≤ z< 10; the transition region 10 ≤ z< 20 and the lower corona z ≥ 20. The solar surface is taken to be the plane at z = 0. The stratification is uniform across the horizontal plane and varies only with height z. The solar interior is taken to be marginally stable to convection by assuming constant entropy in this region. This assumption seems appropriate as the focus of this experiment is the evolution of the magnetic field. The photosphere/chromosphere is taken as an isothermal region with a dimensionless temperature of unity by design. The temperature distribution in the transition region is a powerlaw profile to model the steep temperature gradient observed here. Lastly, the isothermal corona is set at a temperature 150 times that of the photosphere (T_{cor} = 150). To summarise, the nondimensionalised temperature profile is specified as
The pressure and density are then calculated by numerically solving the dimensionless hydrostatic balance equation dp/ dz = −ρg and the dimensionless ideal gas law. The resulting logarithms of the initial temperature, density and pressure of the stratified background are shown in Fig. 1.
Fig. 1
Initial stratification of model atmosphere. The initial profiles are plotted on a log scale against height where red denotes the temperature distribution, black denotes the density and blue represents the magnetic pressure in the solar interior and the gas pressure throughout the domain. 
2.3. Initial magnetic field
We choose to leave the atmosphere unmagnetised by neglecting an ambient field and concentrating on the subphotospheric field. We insert a magnetic flux tube into the hydrostatic solar interior. As we want the setup to remain in equilibrium, and the tube to be in force balance, we require Given that the background environment is in hydrostatic pressure balance, this reduces to (7)where p_{exc} is the pressure excess such that the gas pressure in the tube, p_{t}, is defined to be p_{b} + p_{exc} where p_{b} is the background gas pressure. Explicitly, we have split the gas pressure into a background pressure component that balances gravity and a second component that balances the Lorentz force. The next obstacle is to choose the form of magnetic field to prescribe in the solar interior. The interior magnetic field cannot be observed therefore we choose simple models for the subphotospheric field to initiate emergence. We choose to place a toroidal flux tube with twisted fieldlines in the solar interior. The toroidal model we utilise was derived fully in Hood et al. (2009). Compared with the standard choice of a cylindrical tube, this models the emergence of the top part of an Ωshaped loop that is rooted much deeper in the solar interior. Here, we simply outline the main steps of the derivation.
First, we express our original Cartesian coordinates (x,y,z) in terms of cylindrical coordinates (R,φ,x), given that y is the coordinate describing the direction of the axis of the tube, and z denotes the height from the solar surface as previously. We note that where z_{base} is the value of z at the base of the computational domain. The magnetic field is then expressed, in terms of the flux function A = A(R,x), as where A is constant along magnetic fieldlines. In this derivation, we assume that the magnetic field is rotationally invariant, i.e. independent of φ. We note that this form of magnetic field automatically satisfies the solenoidal constraint (∇·B = 0). Inserting the magnetic field, B, in terms of the flux function, A, into Eq. (7), noting p_{exc} = F(A(R,x)) and RB_{φ} = G(A(R,x)) along with some algebraic manipulation yields, where b_{φ} = RB_{φ}. We note that all the terms are in the direction of ∇A, therefore the GradShafranov equation is of the form (8)Following the strategy used by Hood et al. (2009), we now convert our system to a local toroidal coordinate system (r,θ,φ), such that where R_{0} is the major axis of the toroidal loop and a is the minor radius of the toroidal loop. Eq. (8)can be rewritten in these local coordinates by assuming a ≪ R_{0} and in turn r ≪ R_{0}, i.e. assuming that the minor radius of the torus is much smaller than the major radius. We can expand the solution in powers of a/R_{0} such that to leading order, Eq. (8)becomes
This has exactly the same form as the standard cylindrical equation found in Archontis et al. (2004). We can therefore choose the solutions to be the same as the cylindrical flux tube. Specifically, where B_{0} is the axial field strength and α is the twist. The local toroidal field can also be approximated to O(a/R_{0}). Again using R = R_{0} + rcosθ ~ R_{0} + O(a/R_{0}) this approximates the local toroidal field to With these approximations to the local toroidal field, the pressure can be calculated by exact comparison with the cylindrical model, This balances Eq. (7)and forces the flux tube to sit in equilibrium. However, to initiate the emergence we introduce a density deficit as given by using the temperature profile specified in Sect. 2.2. This makes the flux tube lighter than its surroundings and allows it to begin rising.
Lastly, we need to reexpress the magnetic field in terms of Cartesian coordinates as this is the form that we require for the input of our MHD simulations. Before we can do this directly, we must first express the field in cylindrical coordinates (R,φ,x) as Similarly, this cylindrical field can then be converted into Cartesian coordinates to set the magnetic field as where We note that this corrects an error in Hood et al. (2009) where the sign of B_{R} and B_{x} were interchanged. Fortunately, this error was not significant as it is equivalent to using an initial twist, α, of the opposite sign.
2.4. Parameter choice
Fig. 2
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 purple and an isosurface of the magnetic field ( B  = 1) is overplotted. 
In this experiment, we set the magnetic field strength at the apex of the tube as B_{0} = 9 (11700 G) and the twist as α = 0.4 (righthand twisted). This means that each fieldline undergoes an angle of 0.4 rad per unit distance, and is equivalent to approximately three full turns of twist in the initial subphotospheric field. The base of the computational domain is set at z = −25. The major radius of the torus is R_{0} = 15 and the minor radius is a = 2.5. The total flux threading a cross section of the tube is 6.6 × 10^{19} Mx, typical of a small active region or a large ephemeral region. The initial setup of the experiment is summarised in Fig. 2.
3. Analysis
In order to examine and quantify the rotation of the sunspots in our simulation, we undertake a detailed investigation of a number of quantities relating to both the plasma and magnetic field in an attempt to demonstrate the untwisting of the interior field and rotational flows at the photosphere. First of all, we briefly describe the overall evolution of the field detailing the main phases of emergence.
3.1. Evolution of magnetic field
The flux tube begins to rise buoyantly to the photosphere due to the initial density deficit introduced and continues to rise because of the decreasing temperature stratification in the interior. However, at the photosphere, the plasma is stably stratified with a constant temperature and the flux tube is no longer buoyant. Therefore, the magnetic field finds another way to rise and expand into the corona, namely the magnetic buoyancy instability. In order to initiate this instability, a criterion must be satisfied as derived by Acheson (1979). Typically, the onset of this instability occurs when the plasma β drops to one. At this stage, the magnetic pressure exceeds the gas pressure and the field expands into the atmosphere. See Archontis et al. (2004) and Hood et al. (2012) for a full description of this instability.
Fig. 3
Visualisation of the field in the interior at times t = 40 and t = 100 respectively as traced from the lower negative footpoint (left). A movie of this figure is available online. 
In Fig. 3, we have included two figures illustrating the evolution of the magnetic fieldlines as the flux tube emerges concentrating on the interior portion of the field. We have traced three fieldlines from (0,−14,−25) (blue), (0,−15,−25) (black) and (0,−16,−25) (red) respectively. Initially the flux tube has 3 full turns of twist in total. At t = 40, the flux tube reaches the photosphere and the legs start to straighten. At this time, there is an emerged section of flux tube with about half a turn of twist that will subsequently expand into the corona. However, there is still a considerable amount of twist submerged, approximately a full twist in each leg. In the subsequent evolution there is a definitive unwinding of the submerged twist leading to a final state in which there is virtually no twist in the interior portion of the tube as evidenced by the visualisation of the field at t = 100. We suggest that this may be governed by torsional Alfvén waves. The release of subphotospheric twist as a mechanism for the propagation of a torsional Alfvén wave is investigated later.
3.2. Driver of rotational motion
Before we examine various quantities at the photosphere in an effort to demonstrate the untwisting of the interior field, we discuss the underlying reason for a rotational movement at the photosphere when a twisted magnetic structure emerges. The basic cause for this rotational motion is the behaviour of the Lorentz force. Following the explanation given by Cheung & Isobe (2014), we consider a circular closed curve lying on the photospheric plane enclosing some point P denoting the location of the maximum of B_{z}. We note that the torque is the tendency of a force to rotate an object about an axis, and is given by T = r × F where r is the displacement vector and F is the given force. If we consider the torque due to the various forces acting on the plasma and magnetic field through the surface confined by this closed contour, we find that the magnetic tension is the only force that provides a net torque. That is, the torque contributions from the magnetic pressure and gas pressure forces through this surface vanish and any nonzero torque is a result of the magnetic tension. Explicitly, where r is the displacement vector of a point on the curve from P. This is in fact a general result that can be proved for any force of the form F = ∇f. To demonstrate this we have calculated the torque due to the magnetic tension and magnetic pressure within a circular contour of radius a (2.5) surrounding the location of the maximum of B_{z}, as displayed in Fig. 4. In this case, it is clear that there is no contribution from the magnetic pressure force. Hence, we speculate that the driving motion of the rotation at the photosphere may be governed by the unbalanced torque produced by the magnetic tension force. This is characteristic of an Alfvén wave.
Fig. 4
Surface integral of torque due to the magnetic tension force (red) and the magnetic pressure force (blue) within a circular contour of radius 2.5 around a point P corresponding to the maximum of the vertical magnetic field. 
3.3. Rotation angle
As our previous analyses suggest that the interior portion of the flux tube is untwisting and significant rotational motions develop within the sunspots, we have again traced three fieldlines from the base to the photosphere in an attempt to visualise this motion and quantify the amount by which the fieldlines have rotated. The axis of the flux tube has been traced from the lower negative footpoint as well as two fieldlines either side of the axis in the y −direction, i.e. we have traced fieldlines from (0,−14.5,−25) (blue), (0,−15,−25) (black) and (0,−15.5,−25) (red). A schematic of the traced fieldlines is shown in Figs. 5a and b for times t = 40 and times t = 80 respectively.
Fig. 5
Visualisation of the axis of the flux tube (black fieldline) as well as two fieldlines (red and blue) spaced either side of the axis for comparison at selected times. A movie of this figure is available online. 
Fig. 6
Trajectories of fieldlines as they pass through the photospheric plane coloured with increasing intensity as time progresses. The colour scale on the right shows the times during the evolution. 
These figures clearly demonstrate the rotation of the fieldlines threading the sunspots in a visual manner. Examining Fig. 5, it appears that both the red and blue fieldlines have rotated through an angle of at least π radians over 40 time steps.
In order to follow the fieldlines undergoing this rotation, we have traced the x and y coordinates of the locations of the red, black and blue fieldlines as they pass through the photospheric plane. The trajectories of these fieldlines are shown in Fig. 6a. Initially, the three fieldlines drift outwards in a line as the sunspots separate. Later, the motions slow down and rotation develops. This is not particularly clear given the translational aspect of the motion as the sunspots separate. To remove this feature of the motion, we have subtracted off the location of the black axis from the position of the blue and red fieldlines in Fig. 6b. This gives us the relative position of the blue and red fieldlines and indicates that the blue and red fieldlines rotate clockwise around the black axis by an angle of almost 360°.
With the x and y coordinates of the intersections of selected fieldlines through the photosphere, we can calculate the angle of rotation using (9)where x_{axis} and y_{axis} are the x and y coordinates of the axis of the tube (black fieldline) and x_{0} and y_{0} are the coordinates of the fieldline we are investigating, i.e. the red or blue fieldline. In Fig. 7 a schematic has been included to help us visualise the meaning of the angle φ. Through the use of Eq. (9), the angle φ has been calculated for both the red and blue fieldline as displayed in Fig. 8.
Fig. 7
Representation of angle φ. 
As the two chosen fieldlines were initially equally spaced on either side of the axis, the rotation angles are π out of phase at the beginning. Both fieldlines undergo a rotation of between 7π/ 4 and 2π over 90 time steps. More precisely, the red fieldline undergoes a rotation of 340° and the blue fieldline undergoes a rotation of 353°.
Fig. 8
Evolution of angle of rotation φ for both the red and blue fieldlines as depicted in Figs. 5a and b. 
Next, we consider the temporal rate of change of the angle of rotation, i.e. dφ/ dt, as shown in Fig. 9. This illustrates that different regions of the sunspot are rotating at slightly different rates. There is an initial peak in the size of the rotation rate of both fieldlines of between −0.15 and −0.2. This peak in rotation rate occurs at about t = 44 for the red fieldline and at about t = 62 for the blue fieldline. The rate of rotation diminishes as the experiment proceeds until it reaches zero indicating that the fieldlines have essentially stopped rotating.
Fig. 9
Evolution of rate of change of the angle of rotation φ for both the red and blue fieldlines as shown in Figs. 5a and b. 
Fig. 10
Comparison of terms (solid line) and (dashed line) for both the blue and red fieldlines respectively. 
This rotation is investigated further by checking if the assumption of solid body rotation is reasonable, i.e. that the rotation angle does not depend on the radius of the fieldline. We have concluded from above that different areas of the sunspot are rotating at slightly different rates. If we assume the rotation is a solid body rotation then it follows that the velocity in the φ direction, at radius R, is given by and the zcomponent of the vorticity is given by where we have assumed that φ does not depend on R. We can therefore relate the vertical vorticity and the rate of change of the angle φ by (10)To check if our assumption is valid, we can investigate Eq. (10)by plotting dφ/ dt and ω_{z}/ 2 for both the red and blue fieldlines. In both panels of Fig. 10, the two terms balance each other well suggesting that Eq. (10)is valid and the rotation angle may not have a large dependence on the radius from the axis of the tube.
3.4. Vorticity
In another attempt to demonstrate the rotational movement at the photosphere, we analyse the plasma motion on the photospheric plane. Hence, the vorticity, calculated as the curl of the plasma velocity, is examined as this quantifies the rotation of the plasma. As we are concerned with the rotation of sunspots, the vertical component of the vorticity is the quantity of interest, as it measures the rotation in the x−y plane. This is expressed as, A positive ω_{z} refers to a counterclockwise motion while a negative ω_{z} refers to a clockwise motion. In our experiment, the initial field is righthand twisted; that is, if we were to create this field from a straight field, both footpoints would have been rotated in a clockwise motion. In other words, the fieldlines are wound counterclockwise in the positive B_{z} footpoint and clockwise in the negative B_{z} footpoint.
In order to visualise the field, let us consider the schematic of the fieldlines soon after the field has intersected the photosphere as shown in Fig. 3a. At this point, the legs of the tube have started to straighten out and almost resemble a vertical cylindrical tube originating at z = −25 and intersecting the photosphere at z = 0. As the magnetic field is fixed at the base of the computational domain, a clockwise rotation at the photosphere is required in order to unwind the interior portion of the field. The atmospheric field, on the other hand, is twisted up by a clockwise rotational motion on the two footpoints at the photosphere.
Fig. 11
Coloured contours of ω_{z} accompanied by line contours of B_{z} at z = 0, the base of the photosphere, for the specified times. A movie of this figure is available online. 
A series of coloured contour plots of the vertical vorticity are displayed in Fig. 11 at times t = 40 and t = 80. For visualisation purposes, line contours of B_{z} have been overplotted to show the location of the sunspots. At the centre of each of the polarities, there is a concentration of negative ω_{z} corresponding to a clockwise rotation. This concentration of ω_{z} appears when the legs of the tube straighten and builds with time until it peaks at about t = 50 before decaying as the simulation progresses. This hints that there is some bulk rotation of the sunspots, similar to the sunspot rotations observed by Brown et al. (2003) and Yan et al. (2009). This result complements our previous investigation of the rotation angle where we found fieldlines undergoing significant rotations. Interestingly the same sign of vertical vorticity is present in both concentrations of opposite polarity. As discussed earlier, this has important consequences for the evolution of the field in both the atmosphere and interior. We suggest that these sunspot rotations are due to the untwisting of the field in the interior injecting twist into the atmosphere. Another interesting feature of the contour plots is the red streak of negative vorticity between the sunspots. Elongated streaks of vorticity most likely correspond to shearing motions rather than rotation suggesting that these are due to shear flows between the sunspots. In addition, there are blue tails of positive vorticity located on the inner side of each of the sunspots again indicative of shearing motions.
Now that we have established a clear rotational velocity at the photosphere, the evolution of the vertical vorticity at the photosphere is expressed in a more quantifiable manner. Following the method used by Fan (2009), we achieve this by plotting the time variation of the mean vertical vorticity, ⟨ ω_{z} ⟩, averaged over the area of the upper sunspot where B_{z} is greater than 3/4 of its peak value in Fig. 12. Explicitly, (11)where x_{k} and y_{k} are the x and y coordinates of the region satisfying and N is the number of points that satisfy this criteria. This has been compared to the lower sunspot with no notable difference in result.
Fig. 12
Evolution of the mean vertical vorticity ⟨ ω_{z} ⟩ averaged over the area of the upper polarity concentration where B_{z} is above 75% of the max of B_{z} on z = 0 as described in Eq. (11). 
Considering the average vertical vorticity at the photosphere in Fig. 12, the vorticity is consistently negative suggesting that the dominant motion is a clockwise rotation. We find that a clockwise vortical motion appears in each polarity source when the field reaches the photosphere. The vortical motion quickly rises to a peak at roughly t = 50 soon after the emerged field becomes vertical and the photospheric footprints have reached their maximum separation. At this time rapid rotation commences. The horizontal velocity at this time is approximately 0.5 in magnitude which corresponds to a physical velocity of 3.4 km s^{1}. Soon after, the vorticity steadily begins to decline. This significant clockwise rotation twists up the emerged fieldlines in the atmosphere transporting twist from the tube’s interior portion to its stretched coronal portion. As noted earlier, we speculate that this transport of twist is due to some form of torsional Alfvén wave.
As discussed earlier, Min & Chae (2009) conjectured that the observed rotation of sunspots due to flux emergence may be an apparent effect due to a twisted field rising and each fieldline appearing in a different position at the photosphere. However, in our experiment we have established a clear rotation in the plasma velocity suggesting that this may not be the case. To estimate the contribution to the rotation by apparent effects, we quantify the vertical advection of the flux tube by averaging the vertical velocity in a similar fashion to the way in which we averaged the vorticity. To achieve an upper bound for our estimate, we take the vertical speed of the tube to be the maximum of ⟨ v_{z} ⟩ and assume that the vertical leg has a full turn of twist at t = 40 when the field intersects. This is equivalent to the field being advected vertically by 2.4 units by the end of the experiment, corresponding to a 34.6° apparent rotation. This is an overestimate for the apparent rotation angle as we have assumed the maximum velocity for all time, and yet, this is still significantly smaller than the calculated rotation angle. To conclude, this helps us disregard this theory and explain the rotation solely by a dynamical consequence of the emergence of flux.
Though the vorticity plots do give us a useful insight into the rotational properties of the plasma at the photosphere, further examination is necessary in order to quantify the untwisting of the interior field, and the transport of twist into the atmosphere.
3.5. Current density
A useful quantity when estimating the twist of the magnetic field is the current density, specifically the zcomponent, as denoted by as this is related to how twisted the field is in the x−y plane. We consider coloured contours of j_{z} at a height half way down the solar interior and at the base of the photosphere, as shown in Fig. 13.
Fig. 13
Coloured contours of j_{z} at the plane in the middle of the interior (z = −12.5) in the top panel and at the solar surface (z = 0) in the bottom panel for the specified times, as well as line contours of B_{z} for comparison of the size of sunspots. 
From the coloured contours of j_{z} in Fig. 13, it is clear that although the majority of each sunspot is positive or negative, the periphery of each sunspot is dominated by the opposite sign of j_{z}. As we insist that our initial subphotospheric flux tube is isolated and therefore surrounded by unmagnetised plasma, Faraday’s law demands that the flux tube must carry no net current. As the flux tube carries current inside due to its twist, a reverse current surrounds the sunspot to ensure a zero net current in this region. Focussing on the top panel of Fig. 13, there is some evidence that the two concentrations of strong j_{z} centred on the sunspots are depleting with time. Similarly, considering the bottom panel of Fig. 13, the concentrations of j_{z} at the photosphere intensify when the field first emerges then diminish as the experiment proceeds. As j_{z} signifies the twist of the field in the x−y direction and is related to the azimuthal magnetic field, a decrease in j_{z} could indicate a decline in the amount of twist. The mechanism responsible for this decrease in j_{z} requires further investigation.
Fig. 14
Maximum value of j_{z} against time for varying heights depicted in the key. 
The temporal variation of the maximum value of j_{z} for specific heights below the photosphere is displayed in Fig. 14. There is an initial increase in the maximum of j_{z} for all heights due to the emergence of the field before a steady decline as the experiment proceeds. There are two possible explanations for this steady decrease in the vertical current. This could be caused by the expansion and stretching of the field as a result of emergence or by a decline in the amount of twist stored in the field. The stretching of the field results in a decline in the gradients of B_{x} and B_{y} and lowers j_{z}. To evaluate the extent of the expansion of the field, we estimate the diameter of a contour of B_{z} as shown in Fig. 15. This plot shows the maximum yseparation of the contour of B_{z} = 1. Initially, the separation increases for heights near the photosphere as the flux tube buoyantly rises and the legs of the tube straighten. Later, there is very little change in the separation of the legs for heights deep in the solar interior. There is however some expansion of the magnetic field for the photospheric height z = 0 and 5 units below the boundary as expected as the magnetic field expands into the low density atmosphere. This helps us disregard this cause and explain the decrease in j_{z} in the solar interior solely by an untwisting of the field. Specifically, a decrease in j_{z} implies a decline in the azimuthal field which corresponds to the interior portion of the field untwisting. However, a decrease in j_{z} at the photosphere may be explained by the expansion of the field in this region.
Fig. 15
Line plot of the diameter of the B_{z} = 1.0 contour as a function of time for varying heights depicted in the key. 
Fig. 16
Time evolution of the vertical current ⟨ j_{z} ⟩ averaged over the area of the positive polarity flux source where B_{z} is greater than 75% of its maximum. 
Finally, to analyse the current further, an examination of the total j_{z} in each sunspot is necessary. As we noted from the coloured contours, although the centre of each sunspot is dominated by one sign of current, the outer boundary consists of reverse current. Hence, we cannot estimate the current in each sunspot by measuring the total positive or negative current. Instead, we estimate the current in the centre of the upper sunspot by averaging the vertical current over the area where the vertical magnetic field is greater than 3/4 of its maximum in a similar fashion to our calculation of the average vorticity. The temporal evolution of this quantity is depicted in Fig. 16 for the interior plane (z = −12.5) and the photospheric plane respectively.
At the photosphere, the mean vertical current increases as the magnetic field emerges then steadily declines as the field expands into the atmosphere. Lower down in the solar interior, the mean vertical current generally declines in a linear manner. After an initial drop in ⟨ j_{z} ⟩, there is a small increase as the field straightens out and the negative outerboundary plays a less dominant role. Overall, there is a general decrease in ⟨ j_{z} ⟩ suggesting a drop in the twist stored in the interior portion of the field as the field untwists.
3.6. Magnetic helicity
In order to quantify the transport of twist from the solar interior to the solar atmosphere due to the emergence and untwisting of the field in the interior, we have calculated the evolution of the relative magnetic helicity both above and below the photosphere. The magnetic helicity is a topological quantity describing how much a magnetic structure is twisted, sheared or braided.
The relative magnetic helicity (to the reference field B_{p}) of the field B in a given volume V is given by (Berger & Field 1984, Finn & Antonsen 1985) (12)where A is the vector potential of B (B = ∇ × A), B_{p} is the reference potential field with the same normal flux distribution as B on all bounding surfaces and A_{p} is the vector potential of B_{p} (B_{p} = ∇ × A_{p}). The relative magnetic helicity is favoured over the magnetic helicity as it has been shown to be gaugeindependent with respect to the gauges for A and A_{p}. This is a necessary condition for a physically meaningful definition of helicity.
In order to calculate the relative magnetic helicity numerically, we have tested and compared two approaches; the calculation of A and A_{p} using DeVore’s method (DeVore 2000) and the numerical procedure from Moraitis et al. (2014). As the results were comparatively similar, we will discuss the calculation of the relative magnetic helicity from Moraitis et al. (2014).
In calculating the potential field within the volume V = [x_{1},x_{2}] × [y_{1},y_{2}] × [z_{1},z_{2}], the numerical procedure utilised in Moraitis et al. (2014) takes into account all boundaries within the finite volume. This is advantageous over DeVore’s method which is only valid for a semiinfinite space above a lower boundary. The potential field satisfies j_{p} = ∇ × B_{p} = 0 within V, thus implying B_{p} = −∇φ where φ is a scalar potential. The solenoidal constraint ∇·B_{p} = 0 then implies that the scalar potential is a solution of Laplace’s equation ∇^{2}φ = 0 in V. The condition that B and B_{p} have the same normal components along the boundaries of the volume translates to Neumann boundary conditions for φ, i.e. . The Laplace equation is then solved numerically using a standard FORTRAN routine included in the FISHPACK library (Swarztrauber & Sweet 1979).
The original and potential fields are now stored for the given time step and desired volume. The next step is to calculate the corresponding vector potentials given the method proposed by Valori et al. (2012). As the relative magnetic helicity is gaugeindependent, we are free to choose the gauge A·ẑ = 0 throughout V so that the x and y components of B = ∇ × A are integrated over the interval (z_{1},z) to where A_{0} = A(x,y,z = z_{1}) = (A_{0x},A_{0y},0) is a solution to the zcomponent of B = ∇ × A, i.e. Following Valori’s method we choose the simplest solution to the above equation, given by Similarly, the vector potential of the potential field is calculated using where we have noted that A_{p0} = A_{0} as B and B_{p} share the same normal component on the boundary at z = z_{1}.
An alternative solution for the vector potentials can be obtained if we use the top boundary, i.e. integrating over the interval (z,z_{2}) as This has been checked for comparison and there is no notable difference between the two solutions.
Fig. 17
Evolution of the relative helicity H_{r} when calculated above z = 0 in the solar atmosphere. 
With the corresponding potential and vector potentials for the magnetic field, we can calculate the relative magnetic helicity within different subvolumes of the total simulation domain. First, we calculate the atmospheric magnetic helicity above the photosphere, as shown in Fig. 17. As expected, the atmospheric helicity is zero until t = 20 when the magnetic field first emerges. Subsequently, there is a linear increase in helicity as twist is steadily injected into the atmosphere by the emergence of flux and rotational motions at the photosphere. By the end of the experiment, the magnetic helicity injected into the atmosphere has reached 3.6 × 10^{23} Wb^{2} (3.6 × 10^{39} Mx^{2}), typical of a small event. This can be compared with observations where Min & Chae (2009) quote a helicity transport of 4 × 10^{42} Mx^{2} when considering a much larger active region. To investigate the normalised magnetic helicity that is transported to the atmosphere, we divide by and find that the evolution follows the same linear monotonic shape and reaches a value of 0.83 by t = 120 corresponding to almost one full twist of the flux in the original tube.
To further investigate the magnetic helicity in the atmosphere we have calculated the time derivative of the magnetic helicity above the photosphere, both numerically using finite differencing and analytically. This helps us to understand the main sources contributing to the production and depletion of helicity. Magnetic helicity is mainly contributed to by vertical flows that advect twisted magnetic flux into the corona and by surface flows that shear and twist magnetic fields (Berger & Field 1984). The rate of change of relative magnetic helicity, H_{r}, can be evaluated analytically by differentiating the expression in Eq. (12) (Berger & Field 1984), (13)The first term on the righthand side of Eq. (13)relates to the depletion of helicity by internal dissipation (dissipation term), the second corresponds to a surface correction to the resistive dissipation (surface correction term), the third relates to the generation of helicity by horizontal motions of the boundary (shear term) and the last corresponds to the injection of helicity by direct emergence (emergence term). Let us consider the rate of change of atmospheric helicity in Fig. 18. From the point the field emerges until t = 45, the helicity flux due to emergence (blue solid line) dominates the change in helicity. Later, the horizontal shearing and rotational motions at the photospheric footpoints (red solid line) are the primary sources of helicity change, in agreement with (Fan 2009). The change in helicity due to internal helicity dissipation (purple) and the surface correction term (yellow) are much less significant and do not contribute to the overall change in helicity in the atmosphere. The derivative of H_{r} from Fig. 17 is over plotted for comparison. The two approaches agree very well indicating that numerical effects are kept to a minimum. Care must be taken when considering the rate of change of helicity given in Eq. (13). See Pariat et al. (2015) for the full derivative including additional terms. Clearly, from Fig. 13, the additional terms are not important in this particular case as the flux through the surfaces closely follows the time derivative of the helicity. Pariat et al. (2015) also notes that precaution must also be taken when dividing the helicity flux into individual terms as although their sum is gaugeindependent the individual terms are not, hence, limiting their physical meaning.
Now that we have established a clear increase in the helicity in the atmosphere, we analyse the helicity in the solar interior. As the toroidal flux tube has shown clear signs of untwisting, we expect that the magnetic helicity in the interior will decrease due to both the emergence of magnetic flux and the rotational movements on the photospheric boundary. The calculation of the relative magnetic helicity and the change in helicity below the photosphere are shown in Fig. 19.
Fig. 18
Evolution of the time derivative of the relative helicity H_{r} when calculated above z = 0 in the solar atmosphere using Eq. (13). The total time derivative (black solid line) is split into the dissipation term (purple solid line), the surface correction term (yellow solid line), the shear term (red solid line) and emergence term (blue solid line). The dashed black line is the derivative of the curve from Fig. 17 calculated using finite differencing. 
Fig. 19
Evolution of the relative helicity and rate of change of helicity when calculated below z = 0 in the solar interior. 
As expected, the interior relative magnetic helicity decreases throughout the experiment. Initially, the decrease in magnetic helicity is dominated by internal helicity dissipation. Later, at t = 20, the field emerges from the interior through to the atmosphere resulting in a sharper decrease due to both the emergence of magnetic flux and the rotational movements on the photosphere extracting helicity from the interior as the flux tube untwists.
3.7. Magnetic energy
Now that we have demonstrated the transport of magnetic helicity from the interior portion of the domain to the coronal portion, we consider the magnetic energy and its distribution across the domain. In order to understand the amount of magnetic energy available for solar eruptive events, we consider the free magnetic energy relative to the potential field. That is, we calculate the excess magnetic energy after subtracting off the energy stored in the potential field, i.e. . The evolution of the free magnetic energy above z = 0 is shown in Fig. 20a.
Fig. 20
Evolution of a) the free energy when calculated above z = 0 and b) the vertical Poynting flux through the surface z = 0. The total Poynting flux is split into the emergence term (blue), the shear term (red) and the resistive term (purple) as defined in Eq. (14). 
The excess energy above z = 0 builds from the time the field first emerges. To investigate the contributions to magnetic energy by flux through the photospheric boundary we consider the Poynting flux through z = 0 as given by, (14)The first term contributing to the Poynting flux in Eq. (14)corresponds to the contribution by vertical flows owing to emergence, the second denotes the generation of magnetic energy by shearing/rotational flows and the third term is a result of resistive effects. The rate of increase of energy is largest during the initial stages of emergence due primarily to the emergence term. However, later the shear term (attributed to by rotational motions) is the main contributor to magnetic energy increase in the atmosphere. This pattern corroborates the trend that appeared in the helicity flux whereby vertical flows dominate the flux initially and horizontal flows become important later. Keeping in mind the precaution in Sect. 3.6 about the individual terms of helicity flux, the behaviour of the Poynting flux helps us to believe that the helicity flux trend may have physical meaning. At the end of the experiment, the free magnetic energy transported to the atmosphere has reached 8.2 × 10^{22} J (8.2 × 10^{29} erg).
3.8. Forcefree parameter
Fig. 21
Visualisation of magnetic fieldlines traced from both footpoints coloured by the parameter α such that red represents a strong twist (0.2 <α< 0.4) and blue denotes a weaker twist (0 <α< 0.2). 
Now that we have considered the behaviour of the plasma flows, current and magnetic helicity and energy, a proxy for the local twist is presented. Consider the quantity, α, normally referred to as the forcefree parameter or sometimes the fieldline torsion parameter,
In order to investigate this expression, we trace α along fieldlines as a tool to help us visualise the distribution of twist along the field as shown in Fig. 21. Initially the field is coloured red to signify the field is highly twisted with a value of α greater than 0.2. Rapid emergence results in a coronal magnetic flux that is initially quite untwisted (Longcope & Welsch 2000) as shown by the low coronal α. A clear gradient in α develops along the fieldlines from a high magnitude α in the interior to a lower magnitude of α in the atmosphere as displayed in Fig. 21. Fan (2009) suggests that this gradient produces a torque that drives the rotational motions observed. This result has been proven by Longcope & Klapper (1997) and it is suggested that this motion will continue until the gradient in α is removed (Longcope & Welsch 2000). This agrees with our results.
Later, after the rotational motions have set in at the photosphere, the interior field is left coloured blue with weak twist (α< 0.2) and the atmospheric field is threaded with a mixture of blue, red and white fieldlines. Thus, we can conclude that the field in the legs of the tube undergo a definitive untwisting motion as the field is transformed from highly twisted to weakly twisted by the end of the simulation. Although a large portion of the atmospheric field is coloured blue due to the expansion and stretching of the field, some fieldlines are coloured red indicating that some highly twisted field has been transported into the atmosphere. This is demonstrated in Fig. 21c where the twisted structure of the atmospheric field is shown. As α ~ 1 /L, an increase in the length scale results in a decrease in α. The length scales in the corona are much larger so a decrease in α in the atmosphere can be explained by the expansion and stretching of magnetic fieldlines. The length scales in the interior, on the other hand, are much smaller and so we explain a decrease in α by an untwisting of the field.
3.9. Propagation of torsional Alfvén wave
As touched on earlier, our results indicate that the rotational motions we observe may be governed by some form of torsional Alfvén wave propagating downwards as a consequence of the transport of twist from the interior to the corona. The travel time for an Alfvén wave to propagate from the photosphere to the base of the domain is approximately 20 time steps. This suggests that an Alfvén wave would take approximately 40 time steps to travel down to the base, reflect and return to the photospheric plane. We propose that this downward propagating Alfvén wave, first proposed by Longcope & Welsch (2000), untwists the magnetic field in the interior. The rotation will only slow down once the reflected wave has returned to the photosphere. This appears to be in fairly good agreement with Fig. 12 where the rapid rotation and large  ω_{z}  occurs from about t = 50 to t = 90.
4. Discussions and conclusions
We have presented a 3D MHD numerical experiment of the emergence of a toroidal flux tube from the solar interior through the photosphere and into the solar corona. Based on our detailed analysis, there is evidence that the interior magnetic field untwists and the photospheric footprints undergo a rotation. This rotational motion acts to untwist the interior field fixed at the base and injects twist into the emerged atmospheric portion. Our analysis of the plasma vorticity at the photospheric plane reveals that significant vortical motions develop in the centre of the sunspots. A definitive rotation of the sunspots is also demonstrated by tracking the fieldlines and calculating the rotation rate of the fieldlines threading the sunspot. Rotations of the order of one full rotation (360°) are observed in our experiment. This is similar in magnitude to the angles of rotation reported in studies that concluded a direct relationship between swift sunspot rotation and enhanced eruptive activity (Brown et al. 2003; Yan & Qu 2007 and Yan et al. 2009, etc.). However, the sunspots in our experiment rotate by one full rotation over the course of about 40 minutes, whereas Brown et al. (2003) found sunspots to rotate through angles of the order of 200° over a period of days. The timescales in this experiment are clearly not in line with observations. However, this is related to the size of our emerging active region (~10 Mm). If we scale up our experiment to a more typical active region, we expect the timescales to be much larger and in line with observations. This requires further investigation.
The direct emergence of flux paired with the continual rotational motions at the photosphere transport magnetic energy and helicity from the solar interior to the corona. The magnetic energy in the atmosphere reaches a value of 8.2 × 10^{22}J by the latter stage of the experiment. The Poynting flux of energy is split into contributions from horizontal shearing and vertical emergence terms. The initial flux of energy across the photospheric boundary is dominated by emergence but latterly the primary source is the horizontal shear. The rate of change of relative magnetic helicity in the solar atmosphere also has two predominant sources; namely helicity flux due to emergence and helicity flux by rotational motions. In a similar trend to the energy, initially the predominant source of helicity is the emergence of the magnetic flux tube but later this is replaced by the flux of helicity due to rotational motions at the photospheric level. The magnetic helicity transported to the atmosphere reaches a value of 3.6 × 10^{19} Mx^{2}. As well as the production of helicity in the atmosphere, we find a clear decrease in the magnetic helicity in the interior, supporting our understanding that this portion of the field undergoes an untwisting motion as also evidenced by a clear decrease in j_{z} in the interior.
The appearance of rotational motions centred on both sunspots has been found before in other MHD simulations including Magara (2006) and Fan (2009). Fan (2009) also explained these rotations as a consequence of torsional Alfvén wave propagation and established an increase in helicity in the atmosphere. Our work, however, explicitly discusses the effect that this rotational motion has on the interior portion of the field by establishing a depletion in the magnetic helicity stored in the interior segment of the domain and a drop in the vertical current in this region. We also show that the magnetic tension force may govern this rotational motion as it appears to produce an unbalanced torque that drives the rotation. Furthermore, we explicitly rule out that the rotation observed may be an apparent effect, helping us to explain the rotation primarily as a dynamical result of the emergence of magnetic flux. In addition, we trace fieldlines from the base of the domain as they pass through the photosphere and explicitly calculate their angles of rotation which appear to be approximately in line with the angles calculated in observations. By considering the trajectories of these selected fieldlines, we find a remarkable visual representation of two fieldlines rotating in a circular motion around the axis (the centre of the sunspot), as shown in Fig. 6b.
The simple hypothesis from Longcope & Welsch (2000) has been confirmed in this 3D resistive MHD model, namely that only a fraction of the current carried by a twisted flux tube will pass into the corona and that the propagation of a torsional Alfvén wave at the time of emergence will transport twist from the highly twisted interior to the stretched coronal loops. The rotational motions we observe are a manifestation of the propagation of these waves.
In a future paper we plan to perform a parameter study where we investigate the effects of the field strength B_{0} and twist α on our analysis of the rotation of sunspots. This will allow us to determine if there is a relationship between the strength or twist of the subphotospheric flux tube and the amount of vorticity in the sunspots or the rate at which the sunspots rotate. Another possible avenue for future research is to add in an ambient magnetic field to the coronal portion of the simulation domain. This would add in a further realism and allow us to understand the effect of an ambient field to the rate of rotation, as well as understanding whether rotation can lead to an eruption in our experiments.
Furthermore, we would like to try and understand why the rotation rate we calculate in our simulation is much larger than those calculated in observations. There are many possibilities for this discrepancy, including the size of active region which was discussed earlier. In addition, varying the strength or twist of the tube may change the time it takes for the flux tube to rise to the photosphere and hence govern the rotation rate of the tube. The model presented by Longcope & Welsch (2000) predicts that the level of rotation will depend on the rapidity of flux emergence so we plan to investigate how this affects the rotation. The length of time for the rotation may also be related to the depth at which the flux tube is anchored; this is another approach that requires investigation.
Acknowledgments
Z.S. acknowledges the financial support of the Carnegie Trust for Scotland and CMM the support of the Royal Society of Edinburgh. 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 are grateful to Dr. M. Berger for useful comments concerning magnetic helicity and very thankful to Dr. A. Wright for his helpful discussions and input concerning the propagation of a torsional Alfvén wave. We also thank the unknown referee for the helpful comments.
References
 Acheson, D. J. 1979, Sol. Phys., 62, 23 [NASA ADS] [CrossRef] [Google Scholar]
 Arber, T. D., Longbottom, A. W., Gerrard, C. L., & Milne, A. M. 2001, J. Comp. Phys., 171, 151 [Google Scholar]
 Archontis, V., MoernoInsertis, F., Galsgaard, K., Hood, A., & O’Shea, E. 2004, A&A, 426, 1047 [NASA ADS] [CrossRef] [EDP Sciences] [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]
 Cheung, M. C. M., & Isobe, H. 2014, Liv. Rev. Sol. Phys., 11 [Google Scholar]
 DeVore, C. R. 2000, ApJ, 539, 944 [NASA ADS] [CrossRef] [Google Scholar]
 Evershed, J. 1909, Mon. Not. R. Astron. Soc., 69, 454 [Google Scholar]
 Fan, Y. 2001, ApJ, 554, 111 [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]
 Gibson, S. E., Fan, Y., Mandrini, C., Fisher, G., & Demoulin, P. 2004, ApJ, 617, 600 [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]
 Longcope, D. W., & Klapper, I. 1997, ApJ, 488, 443 [NASA ADS] [CrossRef] [Google Scholar]
 Longcope, D. W., & Welsch, B. T. 2000, ApJ, 545, 1089 [NASA ADS] [CrossRef] [Google Scholar]
 Magara, T. 2006, ApJ, 653, 1499 [NASA ADS] [CrossRef] [Google Scholar]
 Magara, T., & Longcope, D. W. 2003, ApJ, 586, 630 [NASA ADS] [CrossRef] [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 [Google Scholar]
 Pariat, E., Valori, G., Démoulin, P., & Dalmasse, K. 2015, A&A, 580, A128 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Swarztrauber, P. N., & Sweet, R. A. 1979, ACM Trans. Math. Software, 5, 352 [CrossRef] [Google Scholar]
 Török, T., Temmer, M., Valori, G., et al. 2014, Sol. Phys., 286, 453 [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 [Google Scholar]
 Yan, X. L., Qu, Z. Q., Xu, C. L., Xue, Z. K., & Kong, D. F. 2009, Res. Astron. Astrophys., 9, 596 [NASA ADS] [CrossRef] [Google Scholar]
Online material
Movie of Fig. 3 (Access here)
Movie of Fig. 5 (Access here)
Movie of Fig. 11 (Access here)
All Figures
Fig. 1
Initial stratification of model atmosphere. The initial profiles are plotted on a log scale against height where red denotes the temperature distribution, black denotes the density and blue represents the magnetic pressure in the solar interior and the gas pressure throughout the domain. 

In the text 
Fig. 2
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 purple and an isosurface of the magnetic field ( B  = 1) is overplotted. 

In the text 
Fig. 3
Visualisation of the field in the interior at times t = 40 and t = 100 respectively as traced from the lower negative footpoint (left). A movie of this figure is available online. 

In the text 
Fig. 4
Surface integral of torque due to the magnetic tension force (red) and the magnetic pressure force (blue) within a circular contour of radius 2.5 around a point P corresponding to the maximum of the vertical magnetic field. 

In the text 
Fig. 5
Visualisation of the axis of the flux tube (black fieldline) as well as two fieldlines (red and blue) spaced either side of the axis for comparison at selected times. A movie of this figure is available online. 

In the text 
Fig. 6
Trajectories of fieldlines as they pass through the photospheric plane coloured with increasing intensity as time progresses. The colour scale on the right shows the times during the evolution. 

In the text 
Fig. 7
Representation of angle φ. 

In the text 
Fig. 8
Evolution of angle of rotation φ for both the red and blue fieldlines as depicted in Figs. 5a and b. 

In the text 
Fig. 9
Evolution of rate of change of the angle of rotation φ for both the red and blue fieldlines as shown in Figs. 5a and b. 

In the text 
Fig. 10
Comparison of terms (solid line) and (dashed line) for both the blue and red fieldlines respectively. 

In the text 
Fig. 11
Coloured contours of ω_{z} accompanied by line contours of B_{z} at z = 0, the base of the photosphere, for the specified times. A movie of this figure is available online. 

In the text 
Fig. 12
Evolution of the mean vertical vorticity ⟨ ω_{z} ⟩ averaged over the area of the upper polarity concentration where B_{z} is above 75% of the max of B_{z} on z = 0 as described in Eq. (11). 

In the text 
Fig. 13
Coloured contours of j_{z} at the plane in the middle of the interior (z = −12.5) in the top panel and at the solar surface (z = 0) in the bottom panel for the specified times, as well as line contours of B_{z} for comparison of the size of sunspots. 

In the text 
Fig. 14
Maximum value of j_{z} against time for varying heights depicted in the key. 

In the text 
Fig. 15
Line plot of the diameter of the B_{z} = 1.0 contour as a function of time for varying heights depicted in the key. 

In the text 
Fig. 16
Time evolution of the vertical current ⟨ j_{z} ⟩ averaged over the area of the positive polarity flux source where B_{z} is greater than 75% of its maximum. 

In the text 
Fig. 17
Evolution of the relative helicity H_{r} when calculated above z = 0 in the solar atmosphere. 

In the text 
Fig. 18
Evolution of the time derivative of the relative helicity H_{r} when calculated above z = 0 in the solar atmosphere using Eq. (13). The total time derivative (black solid line) is split into the dissipation term (purple solid line), the surface correction term (yellow solid line), the shear term (red solid line) and emergence term (blue solid line). The dashed black line is the derivative of the curve from Fig. 17 calculated using finite differencing. 

In the text 
Fig. 19
Evolution of the relative helicity and rate of change of helicity when calculated below z = 0 in the solar interior. 

In the text 
Fig. 20
Evolution of a) the free energy when calculated above z = 0 and b) the vertical Poynting flux through the surface z = 0. The total Poynting flux is split into the emergence term (blue), the shear term (red) and the resistive term (purple) as defined in Eq. (14). 

In the text 
Fig. 21
Visualisation of magnetic fieldlines traced from both footpoints coloured by the parameter α such that red represents a strong twist (0.2 <α< 0.4) and blue denotes a weaker twist (0 <α< 0.2). 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.