Numerical simulations of rotating axisymmetric sunspots
G. J. J. Botha, F. H. Busse, N. E. Hurlburt, A. M. Rucklidge
aa r X i v : . [ a s t r o - ph ] A p r Mon. Not. R. Astron. Soc. , 1–21 (2007) Printed 2 December 2018 (MN L A TEX style file v2.2)
Numerical simulations of rotating axisymmetric sunspots
G. J. J. Botha ⋆ , F. H. Busse † , N. E. Hurlburt ‡ and A. M. Rucklidge § Department of Applied Mathematics, University of Leeds, Leeds, LS2 9JT, UK Institute of Geophys. Planet. Physics, UCLA, Los Angeles, CA 90024, USA Lockheed Martin Solar and Astrophysics Laboratory, Organization ADBS Building 252, Palo Alto, CA 94304, USA
Accepted 12 April 2008; First submitted 22 February 2008
ABSTRACT
A numerical model of axisymmetric convection in the presence of a vertical magneticflux bundle and rotation about the axis is presented. The model contains a compress-ible plasma described by the nonlinear MHD equations, with density and temperaturegradients simulating the upper layer of the sun’s convection zone. The solutions ex-hibit a central magnetic flux tube in a cylindrical numerical domain, with convectioncells forming collar flows around the tube. When the numerical domain is rotatedwith a constant angular velocity, the plasma forms a Rankine vortex, with the plasmarotating as a rigid body where the magnetic field is strong, as in the flux tube, whileexperiencing sheared azimuthal flow in the surrounding convection cells, forming afree vortex. As a result, the azimuthal velocity component has its maximum valueclose to the outer edge of the flux tube. The azimuthal flow inside the magnetic fluxtube and the vortex flow are prograde relative to the rotating cylindrical referenceframe. A retrograde flow appears at the outer wall. The most significant convectioncell outside the flux tube is the location for the maximum value of the azimuthal mag-netic field component. The azimuthal flow and magnetic structure are not generatedspontaneously, but decay exponentially in the absence of any imposed rotation of thecylindrical domain.
Key words:
MHD — convection — Sun: magnetic fields — sunspots
Observations of sunspots rotating around their own axis,which is perpendicular to the plane of the photosphere, havea long history throughout the twentieth century. Hale (1908)and Evershed (1910) noticed the rotation as well as a vor-tex forming around the rotating sunspot. Observations inthe photosphere and corona continued through the century,culminating in the high resolution measurements of today.(See Brown et al. (2003) and references therein.) An excit-ing development during the last decade has been the abil-ity to measure the associated flow beneath the photosphere(Gizon & Birch 2005; Zhao & Kosovichev 2003).There is a distinct radial profile associated with the az-imuthal velocity of rotating sunspots. Brown et al. (2003)found that the umbra (in which the rotation axis resides)has small average azimuthal velocities, while the fastest ro-tation occurs at some point along the radial length of thepenumbra. The rotation then tails away to a negligible value ⋆ E-mail: [email protected] † E-mail: [email protected] ‡ E-mail: [email protected] § E-mail: [email protected] outside of the sunspot. The peak azimuthal velocity in thepenumbra can be more than double that inside the umbra.Brown et al. (2003) found suggestions of rotation outsidesome sunspots, but these observations are hampered by anambiguous penumbral edge. In contrast, Yan & Qu (2007)observed a rotating sunspot where the maximum azimuthalvelocity occurred inside the umbra. The rotation persisted inthe penumbra and the area near the penumbra, with the an-gular velocity reducing as one moves radially away from theumbra. The surrounding area far removed from the penum-bra experienced a slow rotation in the opposite directionfrom that of the rotating sunspot.It is not clear if the direction of rotation has a hemi-spheric preference. Knoˇska (1975), and references therein,found that the majority of rotating sunspots in both hemi-spheres turn anticlockwise. However, Ding, Hong & Wang(1987) found a preference, with clockwise (anticlockwise)rotation predominantly in the southern (northern) hemi-sphere. This would suggest that the Sun’s differential ro-tation associated with the global flow field has an influ-ence. The Coriolis force due to the Evershed flow fieldwould cause rotation in the opposite direction. How-ever, helioseismic observations (Gizon & Birch 2005), sup-ported by numerical results (Hurlburt & Rucklidge 2000; c (cid:13) Botha, Busse, Hurlburt & Rucklidge
Botha, Rucklidge & Hurlburt 2006), show a converging hor-izontal flow below the Evershed flow, which leads to cyclonicvorticity and hence a possible contribution from the Coriolisforce.The small sample of rotating sunspots studied byBrown et al. (2003) suggests that younger sunspots rotatefaster than older ones. However, it was difficult to judge theages of the sunspots in the sample. The rotation rates aretime dependent, with all rotation eventually decreasing withtime. The peak rotation in the penumbra can be anything upto 3 degrees per hour, as observed by Brown et al. (2003).The same time evolution was observed in a rotating pore(Dorotoviˇc et al. 2002). This behaviour suggest that some ofthe rotation are caused by local events that are transitory innature. This conclusion is strengthened by an observation ofdamped oscillatory motion, which had a maximum rotationof 3.5 degrees per hour (Kuˇcera 1982).A possible mechanism causing rotating sunspots is therise of twisted flux ropes (Gibson et al. 2004). In this modelthe rotation of two flux rope poles is observed after thecentral horizontal portion of the flux rope has emergedthrough the photosphere. This implies the existence of twoco-evolving sunspots of opposite magnetic polarity. This isgenerally not in evidence in the observations, mostly be-cause leading sunspots are often followed by a more diffuseopposite polarity.It was shown by Gopasyuk & Gopasyuk (2005) thatwhen the averaged velocity and magnetic field compo-nents are subtracted from observed sunspots, the resultfits a lightly damped sinusoidal wave. This implies that aswell as the motion described so far, sunspots also experi-ence torsional oscillations. Using the thin flux tube model,Musielak et al. (2007) found that in a compressible, isother-mal field-free medium, linear torsional Alfv´en waves alongthe magnetic tube do not have a cutoff frequency.The rotation of sunspots has been linked to theformation of soft X-ray sigmoids as well as theeruption of flares (Alexander 2006; R´egnier & Canfield2006; Tian & Alexander 2006). Numerical simulations byGerrard et al. (2003) show that by adding a horizontal pho-tospheric flow to the rotation, the generation of flares isenhanced as both rotation and flow increase the complexityof the magnetic field. This is supported by the observationthat flare activity is correlated to magnetic flux and kineticvorticity (Mason et al. 2006).Evidence from helioseismic measurements show that therotation of sunspots, as observed in the photosphere, extendsinto the deeper layers of the sun. Up to a depth of approx-imately 7 Mm vortical flow in the same direction as therotation of the sunspot exists, while below 7 Mm a vorticalflow opposite to the sunspot rotation direction is observed(Kosovichev 2002; Zhao & Kosovichev 2003; Gizon & Birch2005). However, it should be noted that helioseismic mea-surements are difficult and not always consistent. For ex-ample, up to a depth of 3 Mm a converging inflow is foundwhen using p modes, while f mode measurements find onlyoutflows down to 10 Mm (Gizon & Birch 2005).In this paper an axisymmetric model is used to sim-ulate rotation around a central magnetic flux bundle. Thevalues of the physical parameters of the model are chosen todescribe the solar convection zone from a depth of approxi-mately 500 km below the visible surface of the Sun to a depth of approximately 6000 km (Botha, Rucklidge & Hurlburt2006). The numerical domain is a cylinder with an aspect ra-tio of Γ = 3, i.e. one unit deep and a radial distance of threeunits. This implies that we are simulating magnetoconvec-tion on the supergranular scale. To generate azimuthal flowand magnetic field, the whole domain is rotated at a constantangular velocity. Strictly speaking this is not equivalent tosimulating a pore or sunspot where only the magnetic fluxbundle rotates. However, in spite of driving the azimuthalflow throughout the numerical domain, we find that thesolution tends to conform to observations of Brown et al.(2003), where a maximum azimuthal flow occurs close tothe magnetic flux bundle. This means that a vortical flowforms around the flux bundle while the plasma inside thebundle rotates as a solid body. This type of flow is formallydescribed as a Rankine vortex.Our numerical results may be compared to results foundby Jones & Galloway (1993), who studied a Boussinesq fluidin an axisymmetric cylinder. They imposed two types ofboundary conditions: a stress-free outer wall as well as anexternal flow, implemented by rotating the outer wall of thefixed cylinder. A flux bundle formed at the central axis, withthe maximum angular velocity occurring near the axis. Forhigh magnetic field strength, described by the dimensionlessChandrasekhar number ( Q ), the flux bundle broadened withthe stress-free boundary conditions producing a reverse inazimuthal magnetic flux near the axis, while the imposedexternal flow produced a reverse azimuthal flow near theaxis. The reversal in the azimuthal magnetic field is ascribedto the conservation of angular momentum under stress-freeconditions, while the flow reversal obtained with the im-posed external flow is ascribed to the working of the Lorentzforce.In Section 2 the mathematical model and its numericalimplementation is described. This is followed by the numer-ical results. When no rotation of the numerical domain ispresent (Section 3.1), no azimuthal velocity and magneticfield are generated spontaneously: both quantities are smalland decay exponentially with time. This case is useful tocompare the rotating solutions against. Driven by the rotat-ing numerical domain, the solution settles into a time inde-pendent solution that shows rigid rotation of the plasma inthe magnetic flux tube and vortical rotation around it (Sec-tions 3.2 and 3.3). This solution is robust and essentiallystays the same when the magnetic field strength is increased(Section 3.4), the Prandtl number is lowered (Section 3.5),as well as when the stratification is increased (Section 3.7).The latter part of the numerical investigation explores theinfluence of the numerical domain on the solution. This wedo by changing the bottom and outside boundary condi-tions (Sections 3.6 and 3.8). We conclude the paper with asummary of the results. Partial differential equations describing compressible mag-netoconvection are solved in an axisymmetric cylindrical ge-ometry, using a numerical code developed for this purpose.A detailed description of the two dimensional (2D) model isgiven by Hurlburt & Rucklidge (2000). Here we extend themodel to 2.5D by including azimuthal components in addi- c (cid:13) , 1–21 otating axisymmetric sunspots tion to the radial and axial components. A constant angularvelocity is added that introduces the Coriolis and centrifugalforces into the Navier-Stokes equation. The initial temperature and density profiles in the vertical( z ) direction are given by T = 1 + θz, (1) ρ = (1 + θz ) m . (2)The temperature and density are scaled so that they areequal to 1 at the top of the static atmosphere. The initialtemperature gradient is given by θ , while m is the poly-tropic index. The equations for fully compressible, nonlinearaxisymmetric magnetoconvection that we use are ∂ρ∂t = −∇ · ( u ρ ); (3) ∂ u ∂t = − u · ∇ u − Ω × u + Ω (ˆ z × r ) × ˆ z − ρ ∇ P + θ ( m + 1)ˆ z + σKρ ∇ · τ − σζ K Qρ j × B ; (4) ∂T∂t = − u · ∇ T − ( γ − T ∇ · u + γKρ ∇ T + σK ( γ − ρ (cid:16) τ : τ + ζ QK j (cid:17) ; (5) ∂A φ ∂t = ( u × B ) φ − ζ Kj φ ; (6) ∂B φ ∂t = h ∇ × ( u × B ) i φ + ζ K (cid:16) ∇ B φ − B φ r (cid:17) . (7)The cylindrical reference frame is rotated about its axis at aconstant angular velocity of Ω = ( dφ/dt )ˆ z , which is respon-sible for the Coriolis and centrifugal terms in the Navier-Stokes equation (4). The vector potential A φ gives the r and z components of the magnetic field while the azimuthalcomponent is included explicitly, so that the magnetic fieldis given by B = ∇ × (ˆ φ A φ ) + ˆ φ B φ . (8)The velocity consists of three components, namely u =( u r , u φ , u z ), where u φ refers to the azimuthal velocity rel-ative to the rotating reference frame. We also use the auxil-iary equations ∇ · B = 0 , P = ρT, j = ∇ × B , (9)and the following notation: γ is the ratio of specific heats; σ the Prandtl number; ζ the magnetic diffusivity ratio at z = 0; and the Chandrasekhar number given by Q = ( Bd ) µρην , (10)where d is the depth of the domain, µ the magnetic perme-ability, η the magnetic diffusivity and ν the kinetic viscosity.The rate of strain tensor is given by τ = ∂u r ∂r ∂u φ ∂r − u φ r ∂u r ∂z + ∂u z ∂r∂u φ ∂r − u φ r u r r ∂u φ ∂z∂u r ∂z + ∂u z ∂r ∂u φ ∂z ∂u z ∂z , (11)while the dimensionless thermal conductivity K is relatedto the Rayleigh number R in the following way: R = θ ( m + 1) (cid:20) − ( m + 1)( γ − γ (cid:21) (1 + θ/ m − σK . (12) R is a measure of the importance of buoyancy forces com-pared to viscous forces in the middle of the layer, and isused to drive the convection in the model. The addition ofrotation adds an additional scaling which relates the convec-tive timescale to the rotational timescale. Following Gilman(1977) we express this ratio using the convective Rossbynumber defined as Ro = σK √ R. (13)Brummell, Hurlburt & Toomre (1996) found that this pa-rameter, which can be evaluated from the control parame-ters, is typically close to that based on the traditional def-inition of the Rossby number, namely the ratio of the rmsvorticity of the flow to the vorticity 2Ω associated with therotating cylinder. The effects of rotation are significant for Ro ≈
1, and dominate for Ro ≪
1. For supergranules andsunspot moats, Ro ≈ The computational domain is an axisymmetric cylinder ofradius Γ, situated in the ( r, z ) plane so that0 r Γ , z , (14)with z = 0 at the top of the box (Hurlburt & Rucklidge2000; Botha, Rucklidge & Hurlburt 2006). We require thatall variables be sufficiently well-behaved at the axis ( r = 0)and that the differential operators in the PDEs are non-singular. This implies that ∂ρ∂r = u r = u φ = ∂u z ∂r = ∂T∂r = 0 ,A φ = B r = B φ = ∂B z ∂r = j φ = 0 . (15)Terms like u r /r , u φ /r and B φ /r are evaluated usingl’Hˆopital’s rule, while terms containing u r /r cancel alge-braically.The outside wall ( r = Γ) is a slippery, impenetrablewall with no lateral heat flux across it (i.e. an insulator): c (cid:13) , 1–21 Botha, Busse, Hurlburt & Rucklidge ∂T∂r = u r = ∂u z ∂r = B r = B φ = j φ = 0 , ∂u φ ∂r = u φ r . (16)The magnetic potential has the value A φ = Γ / B z = 1.At the bottom boundary the magnetic field is vertical.The temperature T is chosen to be constant with value θ + 1from equation (1). The bottom boundary is impenetrableand stress free, i.e. B r = B φ = ∂B z ∂z = ∂u r ∂z = ∂u φ ∂z = u z = 0 . (17)The top of the box is treated as impenetrable for theplasma, with a radiative temperature boundary conditiongiven by Stefan’s law: ∂u r ∂z = ∂u φ ∂z = u z = 0 , ∂T∂z = θT . (18)The θ in Stefan’s law is the same as in (1), so that theequilibrium profile used as initial condition is not destroyed.The Stefan-Boltzmann constant will enter (18) only whenwe dimensionalize it. The magnetic field is matched to apotential field on top of the numerical domain, ∂A φ /∂z = M pot ( A φ ), where M pot is a linear operator, so that B r and B z are continuous across the boundary. The potential field issolved by assuming an infinitely tall conducting cylinder ofradius Γ above the domain, with the magnetic field becominguniform as z → ∞ . A more detailed description is given byHurlburt & Rucklidge (2000). No currents exist inside thepotential field and consequently we choose j z = 0 along thetop of the box. From (15) it then follows that B φ = 0 alongthe top boundary.One consequence of these boundary conditions is thatno current escapes from the numerical domain: j r = 0 on r = Γ and j z = 0 on z = 0 and 1. It follows that theboundaries do not provide any net vertical torques, and sodo not contribute to changes in the vertical component ofthe total angular momentum, L z . Nonetheless, L z is notconserved in compressible convection: the Coriolis term canlead to changes in L z , for example, when mass is transportedto larger distances from the axis. (Note that this does notoccur in incompressible convection.) A consequence of thisis that as the solution evolves through time, L z tends todrift, making the meaning of the parameter Ω less precise.Therefore, we will look for steady solutions with no net ver-tical angular momentum relative to the rotating frame. Weachieve this by calculating the drift in the value of L z af-ter each iteration and then introducing a correction in theform of an equivalent rigid body rotation in the oppositedirection. This alters the evolution of the PDEs (slightly),but steady states are correct solutions of the PDEs. In ouroscillatory solutions we remove the constraint from L z andfollow its time evolution. This is discussed in Section 3.9.A uniform, vertical magnetic field is used as initial con-dition. For a nonrotating cylinder (Ω = 0) the azimuthalmagnetic field is perturbed (Section 3.1). For a finite angu-lar velocity (Ω = 0) the evolution of the plasma is triggeredby starting the quiet, nonrotating plasma with a finite Ωand no plasma perturbation. Both these initialisations en-sure that L z = 0 at the start of the numerical simulations.The density does not in principle satisfy a boundarycondition, but we impose the value of the normal derivativeof ρ obtained from the Navier-Stokes equation (4).
1. Lines2. Lines3. Arrows4. Colour5. Lines6. Arrows7. Colour 8. Colour9. Colour10. Arrows
Figure 1.
The diagnostics used to describe the numerical results:1. Potential magnetic field lines; 2. Magnetic field lines; 3. Veloc-ity field in the ( r, z ) plane; 4. Temperature fluctuation relativeto the unperturbed state; 5. Density contour lines; 6. Magneticfield direction and strength; 7. Azimuthal current density; 8. Az-imuthal velocity field; 9. Azimuthal magnetic field; 10. Currentdensity in the ( r, z ) plane. All colour scales have red as maximumand blue as minimum. The colour green represents zero.
Figure 2.
No rotation with the azimuthal velocity and magneticfield components approaching zero asymptotically. Consequentlythe right hand side of Figure 1 is not included. The absolutetemperature variation has a maximum of max | ˜ T | = 3 .
1, whilethe azimuthal current density j φ lies in the interval (-200,350). The numerical code was developed specifically for thistype of calculation (Hurlburt & Rucklidge 2000). Sixth-order compact finite differencing is used, which reduces atthe boundaries to fifth-order accuracy for first-order deriva-tives and fourth-order accuracy for second-order derivatives.The grid intervals were chosen to be equal in the r and z directions, with 240 grid points in the horizontal and 80 inthe vertical for the majority of calculations. The time evo-lution obtained fourth-order accuracy through a modified(explicit) Bulirsch-Stoer integration scheme, with the timestep limited by the Courant condition (using the maximumsound and Alfv´en speeds, as well as the thermal diffusivelimit), multiplied by a safety factor of 0.5. Unless otherwise stated, the results shown here have beenobtained with the following parameter values: R = 10 , Q = c (cid:13) , 1–21 otating axisymmetric sunspots Figure 3.
Constant angular velocity Ω = 0 .
1. The time independent solution shows two convection cells in the radial direction instead ofthe one in the result obtained with no rotation (Figure 2). The diagnostics are described in Figure 1, with the azimuthal current density j φ ∈ ( − , r, z ) plane j r ∈ ( − , j z ∈ ( − , | ˜ T | = 2 .
88, and the measured max | u r | = 2 .
36, max | u φ | = 0 .
87, max | u z | = 2 .
35, and max | B φ | = 3 . Figure 4.
Radial profile of u r with Ω = 0 .
1. The solid line is atdepth 0.25, the dotted line at 0.5, and the dashed line at 0.75.
Figure 5.
Radial profile of azimuthal velocity u φ with Ω = 0 . Figure 6.
Radial profile of the axial or vertical velocity u z withΩ = 0 .
1. The three lines have the same meaning as in Figure 4.
Figure 7.
Radial profile of azimuthal magnetic field B φ withΩ = 0 .
1. The three lines have the same meaning as in Figure 4.c (cid:13) , 1–21
Botha, Busse, Hurlburt & Rucklidge σ = 1, ζ = 0 . θ = 10, m = 1, γ = 5 / When no is rotation present (Ω = 0 or Ro → ∞ ), the az-imuthal magnetic field is perturbed initially and the plasmais allowed to evolve through time. The solution reaches thestate depicted in Figure 2, with the explanation of the diag-nostic given by the left hand box in Figure 1. The solutiondescribed here is typical of the plasma state when the ax-isymmetric cylinder is not rotated (Hurlburt & Rucklidge2000; Botha, Rucklidge & Hurlburt 2006) and it provides aconvenient base against which to compare results obtainedwith rotation.From an initial vertical magnetic field, the convectionsweeps the magnetic field towards the central axis whereit forms a flux tube. A large anticlockwise convection cellforms next to the tube, flowing towards the flux tube at thetop of the numerical domain. This flow direction keeps themagnetic field confined to the central axis. The temperatureis time dependent in that a cold plasma blob forms at thetop next to the magnetic flux tube, is convected down theside of the tube, only to dissipate as it is convected along thebottom boundary. Figure 2 shows a new cold plasma blobforming at the top while the remnants of the previous coldblob is still visible at the bottom, moving towards the outerboundary. This temperature oscillation has a period of ap-proximately 1.275 time units, which corresponds roughly tohalf the circulation time around the convection cell. The up-per layers of the plasma are heated by the upflow next to theouter boundary and the resulting hot plasma blob is time in-dependent. The azimuthal current density has its maximumvalue (in both directions) next to the flux tube where themagnetic field gradient is the highest. Azimuthal flow andmagnetic structure are not generated spontaneously. The az-imuthal components, generated by the initial perturbation,are small and decay exponentially as the solution evolvesthrough time.Periodic oscillations, an example of which is the time de-pendent temperature next to the magnetic flux bundle, arefamiliar from Rayleigh-Benard convection (Clever & Busse1995). Jones & Galloway (1993) found periodic oscillationsfor a Boussinesq fluid in an axisymmetric cylinder. As mustbe expected, as in our case they found no spontaneous gen-eration of azimuthal velocity or magnetic field components.Given the model’s temperature and density profiles in(1) and (2), the sound speed increases and the Alfv´en speeddecreases with depth. The inflowing layer at the top of thebox is deeper than the outflowing layer at the bottom. Thisis ascribed to the fact that the total radial momentum in thesystem is zero. The higher density at the bottom leads to ashallower outflowing layer with lower radial velocities trans-porting the same momentum outwards as what the deepertop layer with higher velocities transports inwards. The sim-ulation runs with a maximum Mach number of approxi-mately 1. The time step is limited by the thermal diffusivity,with the dimensionless thermal conductivity K calculatedusing (12). This constraint on the time step is true for allthe numerical simulations in this paper. By introducing a constant angular velocity of Ω = 0 . Ro = 77 .
5) we obtain the solution in Figure 3. The oneconvection cell in the case of no rotation (Figure 2) has splitinto two cells, i.e. a finite Ω reduces the characteristic wave-length of the convection in the radial direction. This effect ofrotation on convection is well known (Chandrasekhar 1961).The convection cell next to the flux tube always has an an-ticlockwise flow direction with an inward flow at the top, sothat it forms a collar which forces the magnetic field togetherat the central axis (Botha, Rucklidge & Hurlburt 2006).By introducing a finite Ω, the centrifugal term in equa-tion (4) provides a force in the radial direction. This man-ifests as a change in density contours, which go from beingapproximately horizontal without rotation (Figure 2) to be-ing slanted at the outer wall (Figure 3). This boundary effectis localized and does not affect the solution in the domaininterior. The treatment of the outer boundary is discussedin more detail in Section 3.8.There is no time dependence in the solution. The radialprofiles of the velocities are given in Figures 4 to 6 at threedepths. All three velocity components are of the same orderof magnitude. The radial velocity (Figure 4) shows the twocells circulating in opposite directions, as well as the factthat the speed is higher in the upper part of the numericaldomain.The azimuthal velocity (Figure 5) shows that theplasma inside the strong magnetic field of the flux tube ro-tates as a solid body, with maximum rotation on the outsideedge of the tube. In the convection area of the solution therotation is in the form of a vortex with the azimuthal velocitygradually falling away with radius. This rotation pattern isuniform throughout the depth of the box and compares wellwith observations that show the largest azimuthal velocitiesare located in the penumbra (Brown et al. 2003). One canfit the profile with a Rankine vortex, described by v φ ( r ) = V rR for r R,V Rr for r > R, (19)with R the magnetic flux tube radius and V = max( u φ ).An observer in the rotating reference frame of the cylinderwill measure an azimuthal velocity profile of v ′ φ ( r ) = v φ ( r ) − Ω r. (20)The values of V and R measured for Ω = 0 . v ′ φ in Figure 8. Rank-ine vortices are used regularly to model tropical cyclones onEarth. Helioseismic measurements of flow around sunspotsin the upper convection zone show a strong resemblance tothe flow of hurricanes on Earth (Zhao & Kosovichev 2003).It is a happy coincidence that Herschel thought of sunspotsas large cyclonic storms (Thomas & Weiss 1992). In ourmodel the Rankine vortex makes physical sense. Convectionis suppressed where the magnetic field is strong. The radialdependence of the azimuthal velocity in these regions is thatof a rotating rigid body. Since the region experiencing rigidrotation corresponds to the magnetized region in all cases,we deduce that magnetic effects are responsible for the rigidbody rotation. A vortex exists around the flux tube. Angular c (cid:13) , 1–21 otating axisymmetric sunspots Figure 8.
Radial profile of a Rankine vortex inside a referenceframe rotating with Ω = 0 . v ′ φ is described by (20) and shouldbe compared with Figure 5. The lines have the same meaning asin Figure 4, and the Rankine constants are listed in Table 1. momentum mixes in axisymmetric convection, which resultsin the free vortex in the field-free convection cells where con-vection is strong. The counter flow near the outer wall is aconsequence of the treatment of L z in our solution. Sinceour solution has zero vertical angular momentum relative tothe rotating reference frame, a significant counter flow hasto occur at the edge in order to balance the peak flow nextto the flux tube.The vertical velocity (Figure 6) shows the strong down-flow at the outside of the magnetic flux tube and at the outeredge of the numerical domain, as well as the upflow betweenthe two convection cells. Comparing the two downflows, weobserve that the downflow at the outer edge is stronger thanthat next to the magnetic flux tube. It is essential to use alarge enough aspect ratio (Γ) so that the outer boundary isremoved from the physics around the magnetic flux tube. AΓ = 3 appears to be a reasonable compromise between thisand the computational limitations.The radial and axial magnetic field components are con-centrated in the magnetic flux tube, with B z three timeslarger than B r . Figure 7 shows the radial profile of B φ , thesize of which is an order of magnitude smaller than B r . B φ is confined mainly to the inner convection cell next to themagnetic flux tube. The current, obtained from the mag-netic field through equation (9), reflects the distribution ofthe magnetic field. Its azimuthal component is concentratedon the outside of the magnetic flux tube where the radialgradient in the magnetic field is the largest. The radial andvertical components are distributed in and around the innerconvection cell around the azimuthal magnetic field maxima(Figure 3). At the top and bottom boundaries the radial cur-rent density has local maxima due to the fact that no currentflows out of the box. The convective Rossby number ( Ro ) associated with the ro-tation around the central axis decreases as Ω increases. The Ro associated with the circulation around the convectioncells for the case with Ω = 0 .
1, i.e. Ro = 77 .
5, compareswell with that of supergranulation in the Sun. Ro of larger Table 1.
The measured constants in the Rankine vortex (19).Ω z (measured downward) V R Ω values correspond to even larger-scale flows. In all casesthe rotation rate is low enough that it should not signifi-cantly change the value of the critical Rayleigh number forthe onset of convection (see Brummell, Hurlburt & Toomre(1996)) and thus the cases exhibit comparable amplitudes.As the magnitude of Ω increases, the width of the mag-netic flux tube increases. Figure 9 shows a time independentsolution with Ω = 0 .
3. The magnetic field strength insidethe flux tube decreases with increasing width, allowing weakconvection to form inside the flux tube itself. Figure 9 showsthat the upflow in the flux tube heats the plasma in the toplayers of the tube, while very weak outflow forms along thetop boundary. Eventually, for Ω > .
3, the convection in-side the flux tube becomes strong enough to break it intoconcentric rings.As rotation increases and with it the width of the mag-netic flux tube, the magnetic field lines forming the flux tubestraighten. This causes the azimuthal current density j φ todecrease in both positive and negative azimuthal directions.The position of j φ stays the same: it flows around the fluxtube, created by large magnetic field gradients there.Increasing Ω also increases the size of the centrifugalforce in equation (4). For Ω = 0 . > . /r dependence in the convec-tion region corresponds to homogeneous angular momentumthat is caused by the effective mixing by the convection. Asthe width of the magnetic flux tube increases with the in-crease in Ω, the radius of the forced vortex also increases(Table 1). The radial profiles of v φ for different Ω valuescan be compared in Figures 5 and 10. They show that asΩ increases, the maxima next to the magnetic flux tubesincrease as well. (See also Table 1.) To maintain the initial L z = 0, the counterflow at the outer wall increases in sym-pathy. This is in contrast to Jones & Galloway (1993), who c (cid:13) , 1–21 Botha, Busse, Hurlburt & Rucklidge
Figure 9.
Constant angular velocity Ω = 0 .
3. The convection inside the magnetic flux tube is stronger and the flux tube wider thanin th case with Ω = 0 . | u φ | = 1 .
5, max | B φ | = 8 .
0, max | ˜ T | = 2 . j φ ∈ ( − , j r ∈ ( − , j z ∈ ( − , Figure 10.
Radial profile of azimuthal velocity u φ with Ω = 0 . found a retrograde flow at the central axis for large Q val-ues in a Boussinesq fluid. To generate a retrograde flow nearthe central axis, we had to change the temperature bound-ary condition at the lower boundary in our model (Section3.6). The rigid body rotation inside the magnetic flux tubeis perturbed when the weak convection inside the tube be-comes strong enough to influence the local magnetic field.Figure 9 shows the strength of the convection inside the fluxtube, while Figure 10 shows the deviation from rigid bodyrotation. This deviation increases deeper in the numericaldomain where convection is stronger.The azimuthal magnetic field tends to be located in theconvection cells closest to the central axis. In the case ofΩ = 0 . > . B φ located in it grows in strength Figure 11.
The behaviour with increasing rotation is presentedusing the following notation: B r is green plus signs connected bya solid line; B φ is blue stars connected by a dot-dashed line; u z iscyan diamonds connected by a dotted line; u r is magenta crossesconnected by a dashed line; u φ is black triangles connected by atriple-dot-dashed line. The peak values are plotted in each case. relative to the B φ in the collar flow outside the flux bundle.The directions of B φ in the small cell inside the flux bundleand that of B φ in the collar flow are anti-parallel. This corre-sponds to the direction of flow of the convection cell. Figure9 shows that a clockwise convection cell has a B φ point-ing in the positive φ direction (i.e. into the page), while B φ in an anticlockwise cell points in the negative φ direction(i.e. out of the page). This is caused by the interaction ofthe magnetic field with the velocity in the first term on theright hand side of equation (7). The current surroundingthe local maxima of B φ has the same direction as the localconvection, since it is calculated using equation (9).The behaviour when rotation is increased is summa-rized in Figure 11. As Ω increases the magnetic flux tube c (cid:13) , 1–21 otating axisymmetric sunspots Figure 12.
Results with Q = 128 and Ω = 0 .
3. The stronger magnetic field widens the magnetic flux tube at the expense of the sizeof the convection area, as seen when compared to Figure 9. The strength of convection and the size of the azimuthal quantities arereduced. Max | B φ | is located in the anticlockwise convection cell next to the flux tube. The measured max | u φ | = 1 .
0, max | B φ | = 2 . | ˜ T | = 2 . j φ ∈ ( − , r, z ) plane is j r ∈ ( − ,
54) and j z ∈ ( − , Figure 13.
Radial profile of u φ corresponding to Figure 12, withΩ = 0 . Q = 128. The lines have the same meaning as inFigure 4. widens and the field lines straighten and become more ver-tical, which leads to a decrease in the radial component ofthe magnetic field. At the same time the azimuthal veloc-ity and magnetic field components increase, being driven byΩ. Compared to these changes, the radial and axial velocitycomponents stay relatively stable. All velocity componentsincrease in absolute value as Ω increases. From previous numerical studies it is known that an increasein the magnetic field increases the width of the magnetic fluxtube for nonrotating solutions (Hurlburt & Rucklidge 2000).This is also true when rotation is present, which can be seenwhen Figure 9 with Ω = 0 . Q = 32 is compared withFigure 12 for Ω = 0 . Q = 128. The solution is timeindependent and the growth in flux tube width takes place at the expense of the radial dimensions of the convectioncells.The magnetic flux tube retains its configuration, withan anticlockwise convection cell holding the flux tube inplace. The stronger magnetic field suppresses the weak con-vection inside the magnetic flux tube, which was presentwhen Q = 32 (Figure 9).As the area with strong magnetic field becomes wider,the field-free convective region is compressed. The decreasein the field-free area is accompanied by lower flow veloci-ties of convection. Here the maximum Mach number of thesolution is 0.8, while max (Mach) = 0 . Q = 32. Themaximum measured azimuthal velocity for Q = 128 (Figure12) is also 2/3 of what it is for Q = 32 (Figure 9).The weaker flow in the convection cell around the mag-netic flux tube means the field lines are less compressed whencompared to the case when Q = 32 (Figure 9). This leadsto lower gradients at the flux tube’s edge, which in turn im-plies a lower azimuthal current density flowing around theflux tube, since j φ is calculated using equation (9).The azimuthal flow of this solution fits that of a Rankinevortex. By suppressing the weak convection inside the fluxtube that is present for Q = 32, the plasma flow inside theflux tube becomes more like rigid body rotation. (CompareFigures 10 and 13.) The maximum u φ next to the flux bundleis lower for Q = 128 due to the lower levels of convectionin the solution. This is also true for the counter flow at theouter wall.The azimuthal magnetic field has its maximum in theconvection cells closest to the central axis. In Figure 9 with Q = 32, there exists a small clockwise cell at the base of theflux tube. This cell is strong enough to contain a significantpart of B φ , with the anticlockwise collar flow containing ananti-parallel B φ component. Here, for Q = 128 (Figure 12),the small clockwise cell inside the flux tube is suppressed, sothat max | B φ | is mainly located in the anticlockwise collarflow. c (cid:13) , 1–21 Botha, Busse, Hurlburt & Rucklidge
Figure 14.
The time independent solution with σ = 0 . Q = 32 and Ω = 0 .
1. The lower σ value enhances convection, which leadsto a larger inner convection cell, a narrower central magnetic flux tube and more variation in the density contours. The measuredmax | u φ | = 1 .
6, max | B φ | = 5 .
5, max | ˜ T | = 3 . j φ ∈ ( − , j r ∈ ( − , j z ∈ ( − , Figure 15.
Radial profile of u φ corresponding to Figure 14, withΩ = 0 . σ = 0 . Q = 32. The lines have the same meaningas in Figure 4. Decreasing the Prandtl number σ causes the convection de-scribed by the steady solution to become more vigorous, ascan be seen when Figure 3 (with σ = 1) is compared toFigure 14 (with σ = 0 . σ = 1 the maximum Machnumber in the solution is 0.9, while for σ = 0 . .
7. The lower σ value brings the simulationcloser to the physical conditions in the upper convectionzone. The dimensionless thermal conductivity K , definedby (12), changes from 4 . × − for σ = 1 to 1 . × − for σ = 0 .
1. For σ = 0 . K = 8 . × − and withΩ = 0 . Ro = 42 . σ , withthe azimuthal component approximately a third the size ofthe radial and axial components. The azimuthal flow takesthe form of a Rankine vortex (Figure 15). The rigid bodyrotation inside the flux bundle is faster than when σ = 1(Figure 5), with a higher maximum next to the bundle. Theflux bundle is also narrower, which means that to maintainthe initial L z = 0 during the simulation, the counter flowat the outer boundary needs to be only slightly larger thanwhen σ = 1.The stronger convection pushes the magnetic flux into athinner flux bundle at the central axis, so that the strengthof B z increases for lower values of σ . The size of B φ relativeto B r stays approximately the same for all values of σ . As inthe case for σ = 1 (Figure 3), the azimuthal magnetic field ismostly located in the inner convection cell, with its maximanext to the magnetic flux bundle. The stronger convectionalso causes the curvature and gradients of the magnetic fieldlines in the ( r, z ) plane to increase, which increases the sizeof the azimuthal current density obtained from equation (9).The position of max | j φ | stays the same.The effects of the enhanced convection are visible inthe density contours. For σ = 1 the density contours inthe convection cells are approximately horizontal (Figure 3)while for σ = 0 . σ = 1 and its lowerconvection strengths.There are changes in the temperature profile of the so-lution. The stronger upflow between the two convection cellsleads to a larger variation from the original heat profile inthe upper layers of the solution, as can be seen when max | ˜ T | in Figures 3 and 14 are compared. Also, for lower σ valuesthe thermal diffusion rate becomes more significant, whichreduces the radial extent of the heated plasma above theupflow. c (cid:13) , 1–21 otating axisymmetric sunspots Figure 16.
Results with a constant ∂T/∂z used as bottom boundary condition. The parameters are Ω = 0 . Q = 32 and σ = 0 .
1. Weakconvection forms inside the wide magnetic flux tube. The overall convection is lower and the inner convection cell smaller compared tosolutions using a constant T lower boundary condition. The measured max | u φ | = 0 .
8, max | B φ | = 3 .
2, max | ˜ T | = 1 . j φ ∈ ( − , r, z ) plane is j r ∈ ( − ,
23) and j z ∈ ( − , Table 2.
Changing the lower temperature boundary condition with Q = 32, m = 1, Γ = 3, θ = 10.Ω σ Lower boundary: T constant Lower boundary: ∂T/∂z constantmax(Mach) T at z = 1 max(Mach) T range at z = 10 1.0 1.2 11 0.6 (9.34,10.19)0.1 1.0 0.9 11 0.6 (9.33,10.42)0.1 0.3 1.7 11 1.1 (9.05,10.43)0.1 0.1 – – 1.3 (8.59,10.20) To determine the extent to which the bottom boundary isinfluencing the result, we changed the temperature prescrip-tion on this boundary from a constant value T to a constant ∂T /∂z . From equation (1) we set ∂T /∂z = θ , so that theheat flux Kθ stays equivalent to the heat flux for a con-stant T boundary condition. A linear stability analysis byHurlburt, Toomre & Massaguer (1984) determined that thischange in boundary condition results in halving the criti-cal Rayleigh number for the onset of convection, and henceone would expect somewhat more vigorous convection forthe same Rayleigh number, all other aspects of the solu-tion being equal. However, the numerical results discussedin this section show that for this highly nonlinear system,the change in the lower boundary condition leads to lowerconvection levels.For low Prandtl numbers, such as σ = 0 . σ > . σ = 1, while Figures 18 to 20 give the radialprofiles of u φ for various parameter values.Figure 17 shows the solution with σ = 1. A solution inthe same parameter space but with constant temperature atthe bottom boundary is given by Figure 3. When compar-ing the two solutions, it is clear that the level of convectionis lower in the case of constant ∂T /∂z boundary condition.This is shown explicitly in Table 2: the maximum Mach num-ber is lower when a constant ∂T /∂z is used. It also showsin the fact that the maximum azimuthal velocity is weakerin Figure 17 than in Figure 3. This is true for all choices ofparameter values. The maximum Mach number in the solu-tion increases as σ decreases, in line with the discussion inSection 3.5 and the fact that more heat flows through thesystem. In Table 2 the solution for σ = 0 . ∂T /∂z does not have a counterpart when a constant tem-perature boundary condition is used, because the convectionbecomes too vigorous and large shocks form that terminatethe numerical simulation. This agrees with the conclusionthat a constant ∂T /∂z at the lower boundary leads to lowerconvection levels.The weaker convection allows the magnetic flux tube to c (cid:13) , 1–21 Botha, Busse, Hurlburt & Rucklidge
Figure 17.
Results with constant ∂T/∂z bottom boundary condition and the same parameter values as Figure 16, but with σ = 1. Oneconvection cell forms with outflow along the top boundary. This allows the magnetic flux to spread radially, with weak convection formingnear the central axis inside the flux tube. The measured max | u r | = 0 .
8, max | u φ | = 0 .
5, max | u z | = 0 .
8, max | B φ | = 3 .
7, max | ˜ T | = 1 . j φ ∈ ( − , r, z ) plane is j r ∈ ( − ,
22) and j z ∈ ( − , Figure 18.
Radial profile of u φ with Ω = 0 . σ = 0 . Q = 32and a constant ∂T/∂z as bottom boundary. The solution has aconfiguration similar to Figure 17, but the weak convection insidethe flux tube has higher amplitudes. These are responsible for thelarge perturbations in the interval 0 r
1. The lines have thesame meaning as in Figure 4. be wider, so that the value of max | B z | is lower. The differ-ent flow pattern when σ > . r, z ) plane. This leads to lower levels of az-imuthal current density, which can be seen when Figures 3and 17 are compared. In the case for σ = 0 . σ > . T lower boundary, theplasma rotates as a rigid body (or forced vortex) where the Figure 19.
Radial profile of u φ corresponding to Figure 17, withΩ = 0 . σ = 1, and Q = 32. Figure 4 defines the line notation. Figure 20.
Radial profile of u φ with Ω = 0 . σ = 0 . Q = 256and constant ∂T/∂z . The solution has a configuration similar toFigure 17. The lines have the same meaning as in Figure 4.c (cid:13) , 1–21 otating axisymmetric sunspots magnetic field is strongest, while a free vortex forms in theconvection area where the magnetic field is weaker. For σ > . σ = 0 . σ = 1.) Moving higher up in the numerical domain,the width of the magnetic flux tube increases and with itthe radius of the forced vortex, with the free vortex in theconvection area occupying less space. At the top of the boxthe magnetic field influences the azimuthal flow so muchthat most of the free vortex flow is distorted. The strongconvection associated with low σ values (Section 3.5) enablesweak convection to form inside the magnetic flux tube. Thisweak flow serves as a perturbation to the rigid body rotationinside the magnetic flux tube around the axis (Figure 18 for σ = 0 . σ one weakens the convection in the so-lution, as discussed in Section 3.5, and reduces the heat fluxthrough the system. This allows the magnetic field to be-come more uniform along the central axis with less pertur-bations in the rigid body rotation that occurs there. Figure17 shows an example of this for σ = 1. An interesting phe-nomenon occurs when σ = 1. In this case the plasma insidethe magnetic flux tube slows down to almost zero (Figure19). The plasma in the field-free convection area still formsa free vortex, which means the azimuthal flow grows froma small value to its maximum next to the flux tube overa small distance, before it tails off into the usual free vor-tex profile. The maximum flow in the vortex occurs closeto the base of the flux tube, where the convection area hasthe lowest level of magnetic field. The fact that the slowlyflowing plasma inside the flux tube is still in the same di-rection as the free vortex is a coincidence. The plasma flowinside the magnetic field area is highly sensitive to the val-ues of Q and σ . Some values give a retrograde rotation atthe central axis, in the same direction as the counter flownear the outer wall. Other values give a prograde rotation inthe direction of the free vortex, but with huge perturbationsin the rigid body rotation, while others will give a flux tubethat rotates partly prograde and partly retrograde. As oneexample of what can occur, we plot the radial profile of u φ for the values Q = 256 and σ = 0 . Q value.) Notice in Figure 20that the free vortex still rotates prograde with a counter flownext to the outer wall, similar to Figure 19 when Q = 32and σ = 1.For a constant T at the lower boundary, the azimuthalmagnetic field formed in the inner convection cell closestto the central axis. Here, for σ = 0 . B φ in it. Themuch stronger convection cell forming the collar flow aroundthe magnetic flux tube contains the rest of B φ . These twocells have opposite meridional circulations, so that the B φ components in them are anti-parallel to each other. For thecases where σ > . B φ )is situated towards the top of the numerical domain, whileit is towards the bottom of the domain for constant T lowerboundaries. This corresponds to the direction of flow in theconvection cell containing B φ in each case.The influence of the outer wall is discernible with con-stant ∂T /∂z at the lower boundary, when only one clockwiseconvection cell forms in the solution (i.e. for σ > . B φ ) forms, generating its own cur-rent around it in the ( r, z ) plane (Figure 17). The heat fluxthrough the system with σ < σ = 1 andthe convection stronger, so that the magnetic field is less ableto concentrate next to the outer wall. One can also see theinfluence of the outside wall in the azimuthal velocity pro-file; the amplitude of the counter flow next to the outer walldiminishes sharply at the wall. This effect becomes largeras the size of the local max( B φ ) at the outer wall increases.(Compare Figures 18 and 19 where σ increases, as well as 18and 20 where Q increases.) In all our results this boundaryeffect is highly localized and does not influence the solutiondeeper in the numerical domain. In Section 3.8 the treat-ment of the outer boundary is discussed.Figures 21 and 22 show that the temperatures at thetop and bottom domain boundaries are lower when con-stant ∂T /∂z is used, compared to a constant T at the lowerboundary. Inside the magnetic flux tube the temperaturenear the top of the domain (Figure 21) is largely indepen-dent of the bottom boundary condition. Where the convec-tion dominates outside the magnetic flux tube, the temper-ature variation is lower than for a constant T at the bot-tom. Figure 22 shows the temperature variation at the lowerboundary. Similar to the top of the domain, the tempera-ture inside the magnetic flux tube is least affected by theboundary condition, although the temperature is lower alongthe whole radial length for ∂T /∂z constant. Figure 22 andTable 2 show that the variation of the temperature alongthe bottom boundary is substantial. This observation mayexplain the difference between the linear stability resultsby Hurlburt, Toomre & Massaguer (1984) and the nonlin-ear behaviour, as mentioned in the beginning of this section.Although the heat flux into the domain stays the same, alower temperature along the bottom boundary will drivethe convection less vigorously. To investigate how sensitivethe velocity amplitudes are to the lower bottom tempera-ture, we increased the bottom boundary temperature andthe heat flux through it by increasing θ . This, however, hasrepercussions throughout the system, as discussed in Section3.7. Increasing the value of θ from 10 to 20 is felt through-out the system. From equations (1) and (2) we see thatthe stratification doubles. Through equation (12) the valueof the dimensionless thermal conductivity K changes from4 . × − for θ = 10 to 1 . × − for θ = 20. This meansthe heat flux through the system ( Kθ ) is raised five fold,as mentioned in Section 3.6. It also implies that the con- c (cid:13) , 1–21 Botha, Busse, Hurlburt & Rucklidge
Table 3.
Changing the stratification of the domain ( θ ) with σ = 1, m = 1, Γ = 3.Ω Q Lower boundary θ = 10 θ = 20max(Mach) T at z = 1 max(Mach) T at z = 10.1 32 ∂T/∂z constant 0.6 (9.33,10.42) 0.8 (17.81,19.53)0.3 128 T constant 0.8 11 1.0 21 Figure 21.
Top boundary temperature with constant ∂T/∂z atbottom boundary. The solid line is σ = 0 . σ = 0 . σ = 1 (Figure 17). The dottedline is σ = 0 . Figure 22.
Temperature at bottom boundary with ∂T/∂z con-stant. The lines have the same meaning as in Figure 21. The extra(triple-dot-dashed) line corresponds to σ = 1 and θ = 20 (Figure24), as discussed in Section 3.7. To fit into the graph, 10 has beensubtracted from this temperature. vective Rossby number (13) increases one order of magni-tude while Ω stays unchanged. Table 3 shows that the max-imum Mach number in the solution increases for both tem-perature prescriptions. The linear stability properties arechanged both by increasing the critical Rayleigh number forthe onset of convection, due to the change in stratification(Hurlburt, Toomre & Massaguer 1984), and by increasingthe mean value of ζ in the domain (Weiss et al. 1990).Increasing θ does not change the basic configuration of the numerical solution. The magnetic flux tube formingat the central axis remains intact, as well as the convec-tion cells in the field-free region. This is true for a constanttemperature as well as a constant ∂T /∂z bottom boundarycondition. Increasing θ increases the strength of convectionin the ( r, z ) plane, as measured by the Mach number in Ta-ble 3, as well as the width of the magnetic flux tube. Thelatter can be observed in Figure 23 for a constant T at thebottom boundary and in Figure 24 for ∂T /∂z constant. Thewider flux tube leads to lower gradients in the magnetic field,which means the azimuthal current density around the tube,calculated using equation (9), becomes weaker. All the otherazimuthal quantities ( B φ and u φ ) are also weaker when com-pared to results with θ = 10.Figure 23 shows a solution with Q = 128 and constant T at the lower boundary. When this is compared to a so-lution with the same parameter values and boundary con-ditions, but with θ = 10 (Figure 12), one observes that thewider magnetic flux tube allows weak convection cells toform inside it. This convection is strong enough to perturbthe temperature inside the flux tube in the top half of thenumerical domain. A careful inspection of the top bound-ary shows that these convection cells cause flow along theboundary, so that concentric rings start to appear at thetop. This is a consequence of the axisymmetry in our model,as one would expect cellular convection to form inside theflux tube. These flows of concentric rings around the centralaxis grow in size as the magnetic flux tube becomes wider,which occurs for higher values of Ω. The azimuthal magneticfield has its maximum value in the strong collar flow next tothe flux bundle. However, the weak convection cells insidethe flux tube are defined well enough for B φ to have sig-nificant components inside them (Figure 23). The directionof each convection cell determines the direction of the local B φ inside it: clockwise convection contains a local maximumand anticlockwise a local minimum. The azimuthal velocityforms a Rankine vortex, as shown in Figure 25. The weakconvection inside the magnetic flux tube perturbs the rigidbody rotation of the plasma inside the tube. The maximumvalue of u φ is found next to the outer edge of the flux tubeand a free vortex forms in the convection area. At the outerwall a counterflow forms, as in the case when θ = 10. A com-parison between the two sets of results ( Figure 13 for θ = 10and Figure 25 for θ = 20) shows that max( u φ ) is lower andsituated farther from the central axis for θ = 20. This meansthe counter flow at the outer boundary, generated becauseour solution has zero vertical angular momentum relativeto the rotating reference frame, is approximately the samestrength for both θ values.Figure 24 shows a solution with Q = 32 and ∂T /∂z constant at the lower boundary. This should be comparedto Figure 17 that has the same parameter values and bound-ary conditions, but with θ = 10. This comparison shows that c (cid:13) , 1–21 otating axisymmetric sunspots Figure 23.
Results with Q = 128, Ω = 0 . θ = 20 and a constant T as bottom boundary. Compared to the results for θ = 10 (Figure12), the magnetic field is more radially dispersed, allowing convection cells to form throughout the radial domain. As one moves towardsthe outer boundary, the upflows between cells grow stronger, accompanied by stronger heating above them. The measured max | u φ | = 0 . | B φ | = 1 .
3, max | ˜ T | = 4 . j φ ∈ ( − , r, z ) plane we have j r ∈ ( − ,
15) and j z ∈ ( − , Figure 24.
Results with Q = 32, Ω = 0 . θ = 20 and constant ∂T/∂z as bottom boundary. Compared to the results for θ = 10(Figure 17), the magnetic field is more radially dispersed, allowing weak convection cells to form inside the flux tube. The strong downflowat the outer edge cools the plasma in the lower layers of the domain. The measured max | u r | = 1 .
0, max | u φ | = 0 .
5, max | u z | = 1 . | B φ | = 2 .
0, max | ˜ T | = 2 . j φ ∈ ( − , j r ∈ ( − ,
7) and j z ∈ ( − , the magnetic flux tube is wider for θ = 20 and that weakconvection occurs inside the tube. This convection is notstrong enough to significantly heat the upper layers of thenumerical domain. In fact, the strong downflows in the so-lution cools the lower layers of the numerical domain muchmore than when θ = 10 (Figure 22). The azimuthal mag-netic field is located inside the large convection cell next tothe magnetic flux tube. As in the case for θ = 10, the maxi-mum value of B φ is located close to the flux tube. Inside thefield-free convection areas we observe a free vortex, with itsmaximum next to the edge of the flux tube and a counterflow next to the outer wall (Figure 26). The rotation of theplasma inside the flux tube is that of a rigid body, but inthe opposite direction from the direction of the free vortexaround the tube. We also see that the plasma inside the hor- izontal magnetic field on top of the convection zone rotatesretrograde, i.e. in the same direction as the counter flow atthe outside wall. This flow pattern is not a surprise, giventhe fact that we could generate retro flows at the centralaxis and the top of the numerical domain by playing withthe parameter values in the set of results with θ = 10. (Seethe discussion in Section 3.6.)The variation in temperature is much larger for thesecases with θ = 20 than in the comparable cases with θ = 10.For a constant T at the bottom boundary (Figure 23) theheating occurring at the top of the numerical domain due tothe strong upflow between the two large convection cells isthree times larger than the heating for θ = 10 (Figure 12).In contrast, the result with ∂T /∂z constant at the bottomboundary (Figure 24) shows that the strong downflow next c (cid:13) , 1–21 Botha, Busse, Hurlburt & Rucklidge
Figure 25.
Radial profile of u φ corresponding to Figure 23, withΩ = 0 . Q = 128 and θ = 20. This should be compared with thecase when θ = 10 in Figure 13. The lines have the same meaningas in Figure 4. Figure 26.
Radial profile of u φ corresponding to Figure 24, withΩ = 0 . Q = 32 and θ = 20. This should be compared to thecase when θ = 10 in Figure 19. The lines have the same meaningas in Figure 4. to the outer wall cools the lower part of the numerical do-main. The amount of cooling is double that which occurs for θ = 10 (Figure 17). When comparing results with a rotating cylinder, one ob-serves a collar flow next to the flux bundle and a clockwiseconvection cell at the outer boundary. The size of the cellat the outer boundary seems to be robust for the variousparameter values. The only exceptions are for ∂T /∂z con-stant at the bottom boundary and with σ > . T bottom boundary condition than when ∂T /∂z is used. In contrast, the azimuthal velocity shows a sharpdecrease in its value at the outside wall when a ∂T /∂z bot-tom boundary condition is used (Figures 18 to 20) while theoutside wall hardly register in the u φ profile for a constant T bottom boundary (Figure 10). The slanted ρ contoursand the decrease in u φ amplitude show how the influence ofthe outer boundary on the solution increases as Ω increases.Due to the formulation of the problem this is unavoidable,but it does not pose a problem as long as these effects staylocalized at the outer boundary.In order to minimize the effect of the outer boundary onthe solution, it was treated throughout as a slippery bound-ary, so that the condition on u φ is given by (16), obtainedfrom the off-diagonal elements of the rate of strain tensor(11). In this paper the boundary conditions at the outerwall were chosen so that the coupling between the numericaldomain and its outside surroundings is kept to a minimum.With boundary conditions (16) only a vertical current existsand the Lorentz force is zero at the outer wall. To measurethe influence of the outer wall on the solution, we changedits magnetic boundary condition to that of a perfect con-ductor. In this case no currents exist parallel to the wall,with a radial current moving through the outer wall. Thecondition that no vertical current exists leads to ∂B φ ∂r = − B φ r , (21)while the radial magnetic field component stays zero, as in(16). For a perfect conductor the Lorentz force has compo-nents parallel to the outer wall, but there is no force acrossthe wall. This implies that there is a torque at the outerboundary, leading to a contribution to the angular momen-tum.Changing the outer boundary conditions on B φ to (21)changes the solution slightly, but only when the azimuthalamplitudes at the boundary become significant when com-pared to the solution near the central axis, as was the case for ∂T /∂z constant at the bottom boundary and σ > . B φ on the outer boundary left the azimuthalvelocity field intact.The difference between boundary conditions (16) and(21) was thoroughly tested. We started numerical runs froma uniform magnetic field for both sets of conditions, whichlead to almost identical time independent solutions. The nu-merical results presented in this paper were started with (21)and then continued with (16). In all cases no significant dif-ference between the numerical solutions could be observed. c (cid:13) , 1–21 otating axisymmetric sunspots Table 4.
Survey of numerical solutions obtained with T constant at the lower boundary and with parameters R = 10 , ζ = 0 . m = 1, γ = 5 /
3, Γ = 3. The star superscript in the Ω column indicates time dependent solutions.Ω
Q σ θ Ro max | ˜ T | max(Mach) max | u φ | max | B r | max | B φ | max | B z | j φ range0 ∗
32 1 10 ∞ (2.9,3.1) (1.2,1.3) → → ∗
128 1 10 ∞ (2.8,2.9) (1.0,1.2) → → ∞ → → ∗
256 1 20 205.8 (3.8,3.9) 0.6 (0.1,1.5) (2.0,2.0) 0.2 (7.0,7.1) (-33,39)0.2 32 1 10 38.7 2.9 0.9 1.3 10.7 4.2 24.4 (-148,271)0.2 32 1 20 102.8 4.7 1.2 1.3 7.4 2.7 17.7 (-74,153)0.2 128 1 10 38.7 2.8 0.8 0.8 6.3 2.6 12.8 (-73,149)0.2 128 1 20 102.8 4.3 1.0 0.6 3.7 1.1 7.5 (-37,64)0.2 256 1 10 38.7 2.7 0.7 0.6 4.4 1.7 8.4 (-56,101)0.2 ∗
256 1 20 102.8 (3.7,3.9) (0.6,0.7) 0.3 (1.8,2.0) (0.3,0.4) (4.5,4.6) (-32,37)0.3 32 1 10 25.8 2.9 0.9 1.5 8.5 8.0 20.3 (-143,220)0.3 32 1 20 68.5 4.6 1.2 1.5 6.3 2.8 14.3 (-66,123)0.3 128 1 10 25.8 2.7 0.8 1.0 5.3 2.9 10.6 (-65,125)0.3 128 1 20 68.5 4.3 1.0 0.7 3.3 1.3 7.8 (-39,57)0.3 256 1 10 25.8 2.6 0.6 0.7 3.8 2.1 7.2 (-55,89)0.3 ∗
256 1 20 68.5 3.7 0.5 (1.6,1.7) (1.6,1.7) 0.4 3.8 (-31,37)
Increasing the angular velocity Ω increases the width of themagnetic flux tube and allows weak convection cells to forminside the tube, similar to those formed in Figure 23. If thevalue of Ω becomes large enough, these cells undergo peri-odic motion. For a constant T bottom boundary condition, Q = 256 and θ = 20, the weak convection cells inside theflux tube oscillate in the radial direction. As Ω increasesfrom 0.1 to 0.3, the amplitude of this oscillation increases aswell. In contrast, for a bottom boundary condition of con-stant ∂T /∂z with Ω = 0 . θ = 10, a hot blob forms dueto weak upflows inside the tube. This blob then moves to-wards the central axis where it dissipates. A new blob formsinside the flux tube and the process repeats itself. This hap-pens for Q values of 128 and 256 with Ω = 0 .
3, as well asfor Q = 256 and Ω = 0 .
2. When the value of θ is doubled to20, the solutions become time independent again.When the forced conservation of L z is lifted, its valuetends to drift, as discussed in Section 2.2. In the case ofa time dependent solution, L z oscillates in sympathy withthe oscillation in the convection, in addition to its drift. Theamplitude of the oscillation can thus be expressed in terms ofan equivalent solid body rotation, which is of order O (10 − ).This compares to a drift in the value of L z of order O (10 − )per unit time. We have investigated magnetoconvection around a magneticflux bundle in a cylinder, when the cylinder is rotated at a constant angular velocity Ω. The model uses a compressibleplasma with density and temperature gradients simulatingthe upper solar convection zone. All the numerical solutionsthat we obtained are presented in Tables 4 and 5. Through-out the calculations the maximum velocities are in the ( r, z )plane, so the the maximum Mach number in these tables area good proxy for u r and u z . For time dependent solutionswe present the range in which the different diagnostics lie.With no rotation (Ω = 0) and a constant temperatureat the lower boundary, the solution is in the form of a fluxtube situated at the central axis, surrounded by a field-freeannular convection ring that forms a collar around the fluxtube (Section 3.1). This magnetic configuration lends itselfto the description of idealized pores and sunspots. The col-lar flow has been measured in the convection around bothphenomena. (See Botha, Rucklidge & Hurlburt (2006) andreferences in it.)The introduction of a constant angular velocity Ωwidens the magnetic flux tube (Section 3.2). Other waysto increase the tube width are to increase the magnetic fieldstrength (Section 3.4) and to increase the heat flux into thenumerical domain from below (Section 3.6). If the magneticfield strength (i.e. Q ) is kept constant and the tube width isincreased by means of one of the above, then the amplitudeof the vertical magnetic field component in the flux tubeis lowered. This allows weak convection cells to form insidethe tube. As Ω increases, the flux tube widens and the weakconvection becomes stronger so that eventually concentricrings appear at the top of the numerical domain (Figures9 and 23). In a fully 3D model one would expect cellularconvection to form inside the flux tube. c (cid:13) , 1–21 Botha, Busse, Hurlburt & Rucklidge
Table 5.
Survey of numerical solutions obtained with ∂T/∂z constant at the lower boundary and with parameters R = 10 , ζ = 0 . m = 1, γ = 5 /
3, Γ = 3. The star superscript in the Ω column indicates time dependent solutions.Ω
Q σ θ Ro max | ˜ T | max(Mach) max | u φ | max | B r | max | B φ | max | B z | j φ range0 32 1 10 ∞ → → ∞ → → ∞ → → ∗
256 1 10 38.7 (1.1,1.5) 0.3 0.5 2.2 (1.8,1.9) (5.6,12.8) (-35,57)0.2 256 1 20 102.8 1.6 0.2 0.3 1.4 0.7 4.3 (-15,26)0.3 32 1 10 25.8 1.2 0.5 0.9 2.6 4.8 24.9 (-54,110)0.3 32 1 20 68.5 2.1 0.6 0.9 2.2 3.4 10.2 (-37,45)0.3 ∗
128 1 10 25.8 (1.1,1.3) (0.3,0.4) 0.7 (1.8,2.3) (2.6,2.8) (6.6,14.2) (-40,62)0.3 128 1 20 68.5 1.9 0.4 0.6 1.5 1.8 7.0 (-24,27)0.3 ∗
256 1 10 25.8 1.1 (0.2,0.3) 0.5 (1.4,1.8) (1.9,2.0) (6.8,7.4) (-33,46)0.3 256 1 20 68.5 1.6 0.2 0.3 1.1 0.9 4.1 (-18,21)
Increasing Ω also brings time dependence to the solu-tion (Section 3.9). For moderate Ω values the weak con-vection cells oscillate horizontally inside the magnetic fluxtube, while for large Ω values the weak cells push period-ically through the edge of the flux tube into the field-freeconvection area. This time dependence can be reduced byincreasing the strength of the magnetic field (Section 3.4).The collar flow around the magnetic flux tube is influ-enced by the strength of the convection and the temperatureprescription at the lower boundary. By lowering the value ofthe Prandtl number ( σ ), the convection becomes strongerand the size of the collar cell increases (Section 3.5). Thestronger convection pushes the magnetic flux tighter at thecentral axis so that the flux tube width decreases and themagnetic field strength on axis increases. For σ = 0 . ∂T /∂z (Section3.6). However, for σ > . u φ and B φ increase (Figure 11). In contrast,the amplitudes of u r and u z hardly change with Ω. For all values of Ω the azimuthal flow pattern fits that of a Rankinevortex: in areas with strong magnetic field the plasma tendsto rotate as a rigid body while around it a free vortex formsin the field-free convection areas. This means that max( u φ )is located outside the flux tube edge. A finite Ω shortensthe wavelength of convection in the radial direction, so thatthe initial convection annulus breaks up into more than oneconvection cell (Section 3.2). The vortex forming around theflux tube is not dependent on the number of convection cellsin the field-free region.The plasma inside the magnetic flux tube and the vortexaround the tube flow prograde relative to the rotating cylin-drical reference frame (Figure 5). A retrograde or counterflow appears next to the outer wall of the cylinder. Thiscounter flow is due to the fact that in our solution the ver-tical component of the angular momentum is zero relativeto the rotating reference frame. We initialize the simula-tions with L z = 0 and the counter flow appears at the outerwall to maintain the status quo. To obtain a retrograde flowat the central axis like Jones & Galloway (1993), we haveto change the bottom boundary condition on the tempera-ture from constant T to constant ∂T /∂z (Section 3.6). Thischange in boundary condition also creates a strong horizon-tal magnetic component in the top layers of the numericaldomain, which may rotate retrograde with the plasma at theaxis and the outer wall (Figure 20). Alternatively, by gen-erating weak turbulence inside the magnetic flux tube, it is c (cid:13) , 1–21 otating axisymmetric sunspots possible to perturb the rigid body rotation of the plasmainside the flux tube to such an extent that one gets pro-grade and retrograde flow inside the flux tube. This is morelikely to happen with a constant ∂T /∂z than for a constant T lower boundary condition.Unlike the azimuthal velocity, the azimuthal magneticfield is influenced by the structure of the convection cells.Max( B φ ) is confined to the strongest convection cell closestto the outer edge of the magnetic flux tube. For a constant T lower boundary condition this is usually the collar flownext to the magnetic flux tube. For constant ∂T /∂z as lowerboundary condition and σ = 0 .
1, significant parts of B φ forminside the weak convection in the flux tube as well as insidecollar cell outside the tube (Figure 16). For constant ∂T /∂z and σ > .
3, the B φ forms in the large convection cell aroundthe flux tube, with a local maximum next to the flux tube(Figure 17). The direction of B φ depends on the convectiondirection; for anticlockwise flow (as in the collar flow aroundthe flux tube) it points in the negative φ direction and viceversa for clockwise flow. When the solution has one largeconvection cell with clockwise flow, which we obtain with aconstant ∂T /∂z at the bottom boundary, a local max | B φ | forms at the outer wall, but its radial width and amplitudeis small so that it does not influence the numerical solutioninside the domain (Section 3.6).The current density in the ( r, z ) plane always formsaround the local max | B φ | and flows in the same directionas the local convection. The azimuthal current density formsaround the edge of the flux tube where the magnetic fieldlines have the largest gradients and curvature. This meansany process that widens the flux tube, i.e. straightens themagnetic field lines, will decrease j φ and vice versa. Increas-ing Ω (Section 3.3), the magnetic field strength (Section 3.4),the stratification in the domain (Section 3.7), and chang-ing the temperature lower boundary condition to constant ∂T /∂z (Section 3.6) lead to a decrease in the amplitudeof j φ , while a lower Prandtl number (Figure 14) increasesthe amplitude of j φ . When the weak convection inside themagnetic flux tube becomes strong enough to bend the fieldlines, local maxima of | j φ | starts to form.Lowering the Prandtl number ( σ ) increases the strengthof convection (Section 3.5) as well as thermal diffusivity.Thus stronger upflows lead to stronger localized heating inthe upper layer, while the radial extent of the heated plasmais reduced (Figure 14). The stronger convection also causessignificant variations in the density gradient inside the field-free convection area. In contrast, a finite Ω with σ = 1 haslittle effect on the density inside the convection area (Section3.3). Only at the outer boundary does the rotation changethe density gradient in a significant way, but the radial ex-tent of this layer is small and does not influence the restof the domain. Inside the magnetic flux tube the density isrelatively unaffected for Ω .
3. Relative large Ω values arenecessary to observe a significant influence by the centrifugalforce.To ascertain the effect of the lower boundary, wechanged the temperature boundary condition (Section 3.6)and the stratification in the numerical domain (Section 3.7).Increasing the stratification effectively increases the heatflux through the lower boundary into the domain. Thiswidens the magnetic flux tube, allowing weak convectioncells to form inside it. However, the convection in the field-
Figure 27.
Velocity field on the ( r, φ ) plane at z = 0 .
25 for Q = 32, σ = 1, Ω = 0 .
1, and constant T bottom boundary, shownin Figure 3. Arrows represent u r and u φ and colour the soundspeed perturbation ˜ c s . Max(˜ c s )=0.54 is the light shade betweenthe two convection cells at r = 1 . c s )=-0.03 the darkshade on the edge of the magnetic flux tube at r = 0 . free regions and the configuration of the magnetic field stayessentially the same (Table 3 and Figures 23 and 24). Incontrast, changing the temperature prescription from a con-stant temperature to constant ∂T /∂z drastically affected thesolution. The bottom temperature reduces slightly (Figure22), but this does not account for the changes in the solu-tion. The amplitude of the convection reduces significantly(Table 2) and for σ > . u φ ) in the photosphere is of the order 10 − km s − (Brown et al. 2003). Zhao & Kosovichev (2003) measuredmax( u φ ) ≈ . − below the photosphere at depths0-3 and 9-12 Mm. This compares well with our measuredvelocities in Tables 4 and 5, where for low angular velocity(Ω = 0 .
1) we obtain max( u φ ) ≈ O (10 − ) km s − , taking asound speed of 1.29 km s − as reference speed.We present Figures 27 and 28 to facilitate comparison ofour results with local helioseismic measurements, of whichFigures 6 to 8 in Kosovichev (2002) are examples. Figure27 shows the flow in the ( r, φ ) plane for a constant T andFigure 28 for a constant ∂T /∂z bottom boundary. Theseplanes correspond to depth z = 0 .
25 in Figures 3 and 17 re- c (cid:13) , 1–21 Botha, Busse, Hurlburt & Rucklidge
Figure 28.
Velocity field on the ( r, φ ) plane at z = 0 .
25 for Q = 32, σ = 1, Ω = 0 .
1, and constant ∂T/∂z bottom boundary,shown in Figure 17. The diagnostics has the same meaning as inFigure 27. Max(˜ c s )=0.04 is the light shade around the magneticflux tube at r = 1 and min(˜ c s )=-0.12 the dark shade at the outerwall at r = 3. spectively, so that Figure 27 represents two convection cellswith a collar flow around a well-defined magnetic flux tube,while Figure 28 represents one outflowing convection cellthat drags the magnetic field lines away from the centralaxis. The arrows show that u r dominates u φ for Ω = 0 . u φ relative to u r increases with Ω (Figure 11), so that the Rankine vortex be-comes more visible for higher values of Ω. In Figures 27 and28 the outer boundary condition u r = 0 still holds, with ar-rows chosen close to the boundary showing that the flow hasfinite size next to the boundary. The colour palette showsthe perturbed sound speed ( c s )in the plane. Where there isan upwelling the plasma is heated and vice versa. Figure 27shows the warmer plasma – and hence higher sound speed –between the two convection cells and Figure 28 at the upflownext to the magnetic flux tube. Downflow with its accompa-nied cooler plasma – and hence lower sound speed – occursaround the edge of the flux tube for Figure 27 and at theouter wall for Figure 28. The flux tube itself is cooler thanthe rest of the surrounding plasma (Figure 27), but whereweak convection inside it exists, it starts to heat up (Figure28). The difference in max( c s ) between Figures 27 and 28 isthus due to the radical different flow pattern for each case.We generate azimuthal flow in this axisymmetric modelby rotating the cylinder around its axis at a constant angularvelocity. As a result we obtain time independent solutions,in contrast to the highly time dependent observations. Forlow angular velocities the flow inside the magnetic flux tubeand the vortex flow are prograde. Due to the fact that ourmodel conserves the vertical component of the angular mo-mentum, a retrograde flow appears next to the outer wall. We find that high angular velocities tend to break the um-bra into concentric rings and introduce time dependence inthe form of periodic behaviour in the radial direction. Thesephenomena have not been observed in sunspots and may bedue to the axisymmetry in our model. It is more likely thatour numerical results obtained with low angular velocitiesare realistic models of solar observations. We would like to thank Chris Jones for informative discus-sions. GJJB and FHB would like to acknowledge financialsupport from NASA grant NNG 04GG07G. GJJB wouldalso like to acknowledge support through PPARC grantPPA/G/O/2002/00014 and STFC grant PP/E001092/1.NEH would like to acknowledge support from NASA grantsNNG06GD45G and NNM07AA01C.
REFERENCES S ., 1975, Bull. Astron. Inst. Czechosl., 26, 151Kosovichev A. G., 2002, Astronomische Nachrichten, 323,186Kuˇcera A., 1982, Bull. Astron. Inst. Czechosl., 33, 345Mason D., Komm R., Hill F., Howe R., Haber D., HindmanB. W., 2006, ApJ, 645, 1543Musielak Z. E., Routh S., Hammer R., 2007, ApJ, 659, 650R´egnier S., Canfield R. C., 2006, A&A, 451, 319 c (cid:13) , 1–21 otating axisymmetric sunspots Thomas J. H., Weiss N. O., 1992, in Sunspots: Theoryand Observations, NATO ASI Series C: Mathematical andPhysical Sciences, Volume 375, page 3, (Kluwer AcademicPublishers, Dordrecht)Tian L., Alexander D., 2006, Solar Physics, 233, 29Weiss N. O., Brownjohn D. P., Hurlburt N. E., Proctor M.R. E., 1990, MNRAS, 245, 434Yan X. L., Qu Z. Q., 2007, A&A, 468, 1083Zhao J., Kosovichev A. G., 2003, ApJ, 591, 446 c (cid:13)000