A numerical model of idealized, axisymmetric, rotating sunspots is presented. The model contains a compressible plasma described by the nonlinear MHD equations, with density and temperature gradients simulating the upper layer of the sun’s convection zone. The solution forms a central flux tube in the cylindrical numerical domain, with convection cells pushing the magnetic field to the axis. When the numerical domain is rotated with a constant angular velocity, the umbra rotates as a rigid body while the surrounding convection cells show a swirling, vortical flow. As a result, the azimuthal velocity and magnetic field have their maximum values close to the flux tube, inside the innermost convection cell

Botha, Gert

Busse, F. H.

Hurlburt, Neal

Rucklidge, Alistair

English

Northumbria Research Link

Citation:  Botha,  Gert,  Rucklidge,  Alistair,  Busse,  F.  H.  and  Hurlburt,  N.  E.  (2006) Numerical  simulations  of rotating  sunspots.  In:  Proceedings of SOHO-17 ‘10 Years of SOHO and Beyond’. ESA SP, 617 . ESA Publications Division,  ESTEC, Noordwijk,  The  Netherlands. ISBN  9789290929284Published by: ESA Publications Division, ESTECURL: http://www.esa.int/esapub/conference/toc/tocSP617.pdfThis  version  was  downloaded  from  Northumbria  Research  Link:  http://nrl.northumbria.ac.uk/10995/Northumbria  University  has  developed Northumbria  Research  Link  (NRL)  to  enable users to access the University’s research output.  Copyright  © and moral  rights  for  items on NRL  are retained by the individual  author(s) and/or other  copyright  owners.  Single  copies of full  items can be reproduced, displayed or performed, and given to third  parties  in  any  format  or  medium  for  personal  research or  study,  educational,  or  not-for-profit  purposes without  prior  permission  or  charge,  provided  the  authors,  tit le  and  full  bibliographic  details  are  given,  as  well  as  a  hyperlink  and/or  URL  to  the  original  metadata  page. The content  must  not  be changed in  any way.  Full  items must  not  be sold commercially  in  any format  or medium  without  formal  permission of the copyright  holder.  The full  policy is available online: http://nrl.northumbria.ac.uk/policies.htmlThis document  may differ  from the final,  published version of the research and has been made available online in  accordance with  publisher  policies. To read and/or cite from the  published  version  of the  research,  please visit  the  publisher’s  website  (a subscription  may be required.)NUMERICAL SIMULATIONS OF ROTATING SUNSPOTSG. J. J. Botha1, A. M. Rucklidge1, F. H. Busse2, and N. E. Hurlburt31Department of Applied Mathematics, University of Leeds, Leeds, LS2 9JT, UK, Email: gert@maths.leeds.ac.uk;A.M.Rucklidge@leeds.ac.uk2Institute of Geophys. Planet. Physics, UCLA, Los Angeles, CA 90024, USA, Email: fbusse@igpp.ucla.edu3Lockheed Martin Solar and Astrophysics Laboratory, Organization L9-41 Building 252, Palo Alto, CA 94304, USA,Email: hurlburt@lmsal.comABSTRACTA numerical model of idealized, axisymmetric, rotatingsunspots is presented. The model contains a compressibleplasma described by the nonlinear MHD equations, withdensity and temperature gradients simulating the upperlayer of the sun’s convection zone. The solution formsa central flux tube in the cylindrical numerical domain,with convection cells pushing the magnetic field to theaxis. When the numerical domain is rotated with a con-stant angular velocity, the umbra rotates as a rigid bodywhile the surrounding convection cells show a swirling,vortical flow. As a result, the azimuthal velocity and mag-netic field have their maximum values close to the fluxtube, inside the innermost convection cell.Key words: MHD; convection; Sun: magnetic fields;sunspots.1. INTRODUCTIONPhotospheric white light measurements of rotatingsunspots (Brown et al., 2003) show a rotation of up to200 degrees about the umbral centres over 3 to 5 days.Sunspots have also been observed to undergo dampedoscillations (Kucˇera, 1982). Helioseismic measurementsbelow the surface (Zhao & Kosovichev, 2003) reveal thatat a depth of 0 to 12 Mm there is evidence of structuraltwist. In addition, at a depth of 0 to 5 Mm there existsubsurface horizontal vortical flows, while there are alsosuggestions that below 9 Mm a vortical flow in the oppo-site direction to the above may exist.To investigate the physical processes associated withrotating sunspots, a numerical model of an idealizedsunspot is presented. The sunspot is placed in an axisym-metric cylinder with an aspect ratio (radius versus depth)of 3, which is then rotated at a constant angular veloc-ity. The model is described in Sections 2 and 3. Two setsof results are discussed (Sections 4 and 5) that were ob-tained with different temperature boundary conditions atthe bottom of the numerical domain. Finally, the effect ofthe amplitude of the reference frame’s angular velocity isdiscussed in Section 6.2. MATHEMATICAL MODELThe initial temperature and density profiles areT = T0(1 + θz), (1)ρ = ρ0(1 + θz)m, (2)with the 0 subscript defining the quantity at the top of thebox (z = 0), θ the initial temperature gradient, and mthe polytropic index. Throughout we have used T0 = 1,ρ0 = 1, θ = 10 and m = 1. We solve the equations forfully compressible, nonlinear axisymmetric magnetocon-vection (Hurlburt & Rucklidge, 2000; Botha et al., 2006):∂ρ∂t= −∇ · (uρ) (3)∂u∂t= −u · ∇u− 2Ω× u+Ω2r+ θ(m+ 1)zˆ−1ρ∇P + σKρ∇ · τ − σζ0K2Qρj×B (4)∂T∂t= −u · ∇T − (γ − 1)T∇ · u+ γKρ∇2T+σK(γ − 1)ρ(12τ : τ + ζ20QK2j2)(5)∂Aφ∂t= (u×B)φ − ζ0Kjφ (6)∂Bφ∂t=[∇× (u×B)]φ+ ζ0K∇2Bφ (7)The vector potential Aφ gives the r and z componentsof the magnetic field while the azimuthal component isevolved explicitly, so that the magnetic field is given byB = ∇× (φˆAφ) + φˆBφ. (8)The velocity consists of three components, namely u =u(r, φ, z). We also use the auxiliary equations∇ ·B = 0, (9)P = ρT, (10)j = ∇×B, (11)and the following notation: τ is the rate of strain tensor;γ the ratio of specific heats; σ the Prandtl number; Kthe dimensionless thermal conductivity; ζ0 the magneticdiffusivity ratio at z = 0; and Q is the Chandrasekharnumber. All the other symbols have their usual meaning.The physical quantities are made dimensionless by scal-ing length ∝ depth; time ∝ depth / (sound speed at top);and temperature, magnetic field, density, and pressure allproportional to their initial values at the top of the numer-ical domain. The cylindrical reference frame is rotated inthe plane perpendicular to the axis at a constant angularvelocity Ω = (dφ/dt)zˆ, which results in two additionalterms in the Navier-Stokes equation (4).3. NUMERICAL MODELmagnetic field: potentialtemperature: Stefan’s lawtemperature: fixedmagnetic field: verticalslippery,perfectly conducting,thermally insulatingwall1z00 r ΓFigure 1. The 2D numerical domain (r, z) of an axisym-metric cylinder with radius Γ and with z = 0 at the topof the box.The computational domain is given in Figure 1. We re-quire that all variables be sufficiently well-behaved at theaxis (r = 0) and that the differential operators in thePDEs are non-singular. This implies that∂ρ∂r= Ur =∂Uz∂r= Aφ = Br =∂Bz∂r= j =∂T∂r= 0.(12)Terms like u/r are evaluated using l’Hoˆpital’s rule. Inthe numerical simulations we use a fourth-order Bulirsch-Stoer time integration, with sixth-order compact finitedifferencing (Lele, 1992). The time step was limited bythe Courant condition (taking the maximum sound andAlfve´n speeds, as well as thermal diffusive limits into ac-count), multiplied by a safety factor of 0.5. A uniform,vertical magnetic field was used as initial condition, andno perturbation was given to the plasma. The time evolu-tion of the plasma is triggered by starting the quiet initialstate with the constant angular velocity |Ω| of the refer-ence frame.4. CONVECTION WITH CONSTANT T AT BOT-TOM BOUNDARYThese results were obtained with Rayleigh number R =105, Q = 128, Ω = 0.1, σ = 1, ζ0 = 0.2, γ = 5/3Figure 2. The radial profile of the azimuthal velocity sam-pled from Figure 3. The solid line is sampled at depth0.25, the dotted line at depth 0.5, and the dashed line at0.75.and Γ = 3. The solution is shown in Figure 3 when theplasma is almost time independent. All the quantities inthe (r, z) plane show no time evolution, while the max-ima of the azimuthal quantities are increasing at a veryslow rate.The convection forces the magnetic flux to the centralaxis where a strong flux tube forms (Hurlburt & Ruck-lidge, 2000). Outside the region with strong magneticfield two convection cells form, the inner cell rotating an-ticlockwise, i.e. at the top towards the magnetic flux tube.Inside the flux tube the density stratification stays as ini-tialized due to the low level of convection there, whileoutside the flux tube the density shows the effect of thehigh levels of convection. The radial profile of the az-imuthal velocity is presented in Figure 2, taken at dif-ferent depths from the solution as shown in Figure 3. Itshows that the magnetic flux tube rotates as a solid body,so that the maximum azimuthal velocity is reached justoutside the magnetic flux tube. Outside the flux tube thevelocity decreases in such a way that there exist a swirlingvortex around the magnetic flux tube. Figure 2 shows thatthis is true at all depths in the numerical domain. The az-imuthal magnetic field is situated mostly inside the innerconvection cell (Figure 3), so that its maximum is closeto the magnetic flux tube.5. CONVECTION WITH CONSTANT ∂T/∂Z ATBOTTOM BOUNDARYThis result was obtained with the same parameters thatwere used in the case of a fixed temperature at the bot-tom boundary (Section 4). The final state of the solutionis shown in Figure 4. Again there is little time evolutionpresent, with only the amplitudes of the azimuthal quan-tities increasing at a very slow rate. The result showsthe influence of the choice of boundary conditions on theconfiguration of the convection cells. The main differ-                     Figure 3. The diagnostic is divided into five boxes. The top box on the left hand side shows the magnetic field lines ascalculated for a potential field. The middle layer on the left hand side shows the temperature in colour (blue is cold andred is hot), the velocity field as arrows, and the magnetic field lines as contour lines. The bottom box on the left hand sideshows the azimuthal current in colour, the arrows are the magnetic field indicating size and direction, and the contourlines represent the density. The top box on the right hand side shows the azimuthal velocity in colour (blue being negativeand red positive), while the bottom box on the right hand side shows the azimuthal magnetic field in colour and the currentin the (r, z) plane as arrows. The state of the plasma shown here was obtained with a constant T at the bottom boundary.Its azimuthal velocity shows the rigid rotation of the magnetic flux bundle at the centre of the cylinder, as well as theswirling flow outside the flux tube with the maximum azimuthal flow located in the anticlockwise convection cell next tothe flux bundle. The solution shows no time dependence in the (r, z) plane, with the azimuthal quantities increasing veryslowly in amplitude.                     Figure 4. The same diagnostic as was used in Figure 3. This result was obtained with a constant ∂T/∂z at bottomboundary of the numerical domain, while all the other parameters were kept the same as in Figure 3. The convection cellclosest to the magnetic flux bundle has changed direction, resulting in the magnetic flux bundle not being well confinedby the convection. The azimuthal quantities show rigid rotation in the flux bundle and a swirling vortical flow in theconvection cell. There is no time dependence in the (r, z) plane, with the azimuthal quantities increasing very slowly inamplitude.Figure 5. The magnetic flux tube rotates as a solid body.Presented here is the angular velocity of the flux tube asa function of the angular velocity (Ω) of the referenceframe, which is constant for any given solution. The mea-surements were obtained with a constant T at the bottomboundary, as shown in Figure 3.ence compared to Section 4 (and Figure 3) is a clockwiseconvection cell that dominates outside the flux tube. As aresult, the magnetic flux is not as confined to the centralaxis as with anticlockwise convection and the azimuthalmagnetic field has changed direction.In spite of the differences between Figures 3 and 4, themagnetic flux tube rotates as a solid body while outsidethe region with strong magnetic field the azimuthal con-vection follows a swirling vortical motion as much as themagnetic field configuration allows this to happen. Asa result, the maximum of the azimuthal velocity is nextto the magnetic flux tube. The azimuthal magnetic com-ponent is distributed across the convection cell closest tothe magnetic flux bundle, as was the case in Section 5 andFigure 3.6. RIGID BODY ROTATIONThe magnetic flux tube is well defined for simulationswith a constant T at the bottom boundary. This flux tuberotates as a solid body around the axis. Although the fluxtube is less defined with a constant temperature flux asbottom boundary condition (Figure 4 and Section 5), theregions with strong field still rotate as a solid body. Inorder to investigate the influence of the size of Ω on theangular velocity of the plasma, the solid body rotation ofthe magnetic flux tube was monitored for the case of aconstant T at the bottom boundary (Figure 3 and Section4.) Figure 5 shows that there is no direct relation be-tween the angular velocity Ω of the reference frame andthe rotation of the magnetic flux tube. It seems that theconfiguration and vigour of the neighbouring convectioncells have a larger influence on the rotation of the mag-netic flux tube than the amplitude of Ω. By increasingthe size of Ω the strength of convection is increased, butthe configuration of the convection cells is dependent onthe bottom temperature boundary. An increased Ω alsowidens the magnetic flux tube at the centre. When it be-comes wide enough, convection starts to form inside theflux tube. In contrast, as the magnetic field strength (i.e.Q) increases, the solid body rotation slows down. Thiscan be explained by the fact that a strong magnetic fieldinhibits convection.7. SUMMARYRotating sunspots were simulated by evolving a plasmain an axisymmetric cylinder, with the cylinder rotating ata constant angular velocity Ω. This Ω drives the veloc-ity components of the plasma through the Navier-Stokesequation (Eq. 4) and as such no initial perturbation of theplasma was needed. Convection cells form in the solu-tion that push the magnetic field towards the central axisof the cylinder to form magnetic flux tubes. These fluxtubes display rigid body rotation while the plasma out-side the tube experiences a swirling vortical flow. Anazimuthal magnetic field forms with a maximum in theconvection cell closest to the magnetic flux tube. By in-creasing the amplitude of Ω, the width of the magneticflux tube at the central axis increases, which allows con-vection to form inside the flux tube in the lower localmagnetic field strength.ACKNOWLEDGEMENTSThis work was done with financial support from PPARCand NASA.REFERENCESBotha G.J.J, Rucklidge A.M., Hurlburt N.E., 2006, MN-RAS, accepted for publicationBrown D.S., Nightingale R.W., Alexander D., SchrijverC.J., Metcalf T.R., Shine R.A., Title A.M., Wolfson C.J.,2003, Solar Physics, 216, 79Hurlburt N.E., Rucklidge A.M., 2000, MNRAS, 314, 793Kucˇera A., 1982, Bulletin of the Astronomical Institutesof Czechoslovakia, 33, 345Lele S.K., 1992, J. Computational Physics, 103, 16Zhao J., Kosovichev A.G., 2003, ApJ, 591, 446

Numerical simulations of rotating sunspots

Abstract

Similar works

Full text

Available Versions

Northumbria Research Link