Computing nonlinear force free coronal magnetic fields

0.41 MB·English

by  T. Wiegelmann

8
0
0
2
Computing nonlinear force free coronal magnetic fields.

n
a
J

1
2 T. Wiegelmann1,2 and T.Neukirch1

] 1School ofMathematicsand Statistics,UniversityofSt. Andrews,St. Andrews,KY169SS,
h
p UnitedKingdom
- 2Max-Planck-Institutfu¨rAeronomie,Max-Planck-Strasse2, 37191Katlenburg-Lindau,
o
r Germany
t
s
a
[

1
v
5
1
2 Camera-ready Copy for
3
. NonlinearProcesses inGeophysics
1
0 Manuscript-No. NonlinearProcesses inGeophysics,Volume10,Issue4/5,2003,pp.313-322
8
0
:
v
Offsetrequests to:
i
X
T. Wiegelmann
r
a
NonlinearProcessesinGeophysics(2002)25:1–10
Nonlinear Processes
in Geophysics
(cid:13)cEuropeanGeophysicalSociety2002
Computing nonlinear force free coronal magnetic fields.
T.Wiegelmann1,2andT.Neukirch1
1SchoolofMathematicsandStatistics,UniversityofSt. Andrews,St. Andrews,KY169SS,UnitedKingdom
2Max-Planck-Institutfu¨rAeronomie,Max-Planck-Strasse2,37191Katlenburg-Lindau,Germany
Received:12. August2002–Revised:22.November2002–Accepted:27. November2002
Abstract. Knowledge of the structure of the coronal mag- measurement of the coronal magnetic field is not feasible
netic field is importantforourunderstandingof manysolar with the observationalmethods presently available. There-
activity phenomena, e.g. flares and CMEs. However, the fore the extrapolation of the magnetic field into the corona
direct measurementof coronalmagnetic fields is not possi- frommeasurementsofthemagneticfieldatthephotospheric
blewithpresentmethods,andthereforethecoronalfieldhas levelisextremelyimportant.
tobeextrapolatedfromphotosphericmeasurements. Dueto It is generally assumed that the magnetic pressure in the
the low plasma beta the coronal magnetic field can usually coronaismuchhigherthantheplasmapressure(smallplasma
beassumedtobeapproximatelyforcefree,withelectriccur- β) and that therefore the magnetic field is nearly force-free
rentsflowingalongthemagneticfieldlines. Thereareboth (foracriticalviewofthisassumptionseeGary,2001). The
observational and theoretical reasons which suggest that at extrapolationmethodsbasedonthisassumptionincludepo-
least prior to an eruption the coronal magnetic field is in a tentialfieldextrapolation(e.g.,Schmidt,1964;Semel,1967),
nonlinearforcefreestate. Unfortunatelythecomputationof linear force-free field extrapolation (e.g., Chiu and Hilton
nonlinearforcefree fields is way more difficultthan poten- , 1977; Seehafer , 1978, 1982; Semel , 1988) and nonlin-
tial or linear force free fields and analytic solutions are not ear force-freefieldextrapolation(e.g.,Amarietal. , 1997).
generallyavailable. Wediscussseveralmethodswhichhave Methodsfornon-force-freefieldextrapolationhavealsobeen
beenproposedtocomputenonlinearforcefreefieldsandfo- developed(PetrieandNeukirch,2000)butarenotusedrou-
cus particularly on an optimizationmethod which has been tinely.
suggestedrecently. Wecomparethenumericalperformance Whereaspotentialandlinearforce-freefieldscanbeused
ofanewlydevelopednumericalcodebasedontheoptimiza- as a first step to model the general structure of magnetic
tionmethodwiththeperformanceofanothercodebasedon fields in the solar corona, the use of non-linear force free
an MHD relaxationmethodif bothcodesare appliedto the fieldsisessentialtounderstanderuptivephenomenaasthere
reconstruction of a semi-analytic nonlinear force-free solu- arebothobservationalandtheoreticalreasonswhichsuggest
tion. Theoptimizationmethodhasalsobeentestedforcases thatthepre-eruptivemagneticfieldsarenon-linearforce-free
where we add random noise to the perfect boundary con- fields (see below). The calculation of nonlinear force free
ditions of the analytic solution, in this way mimicking the fieldsis complicatedby theintrinsic nonlinearityof theun-
morerealisticcasewheretheboundaryconditionsaregiven derlyingmathematicalproblem.Anotherproblemwhichbe-
by vector magnetogramdata. We find that the convergence comes especially important when nonlinear force free field
propertiesoftheoptimizationmethodareaffectedbyadding areusedformagneticfieldextrapolationisthecorrectformu-
noise to the boundary data and we discuss possibilities to lation of the problemwith respectto boundaryvalues. The
overcomethisdifficulty. available magnetographobservationprovideeither the line-
of-sightmagneticfield(B ), whichissufficientforpoten-
los
tial and linear force free fields, or all three components of
thephotosphericmagneticfield. Onlythelatterinformation
1 Introduction
is sufficient to determine non-linear force free fields com-
pletely. There are, however, tremendous difficulties to be
Themagneticfieldisthemostimportantquantityforunder-
overcomebothinextractingthenecessaryinformationabout
standingtheplasmastructureofthesolarcoronaandtheac-
themagneticfieldcomponentsontheboundaryfromthedata
tivityphenomenaitshows. Unfortunatelyasystematicdirect
inanunambiguouswayandintheuseofthisinformationin
Correspondenceto: T.Wiegelmann the computationof the correspondingmagnetic field. Even
1
2 Coronalmagneticfields.
thoughalotofresearchhasbeendevotedtotheseproblems oflinear(orconstant-α)forcefreefields. Standardmethods
in the past so far none of the proposed methods has been forcalculatingpotentialandlinearforcefreefieldareavail-
found to be outstanding in combining simplicity, usability able(seee.g.Gary,1989)andtheycanatleastgivearough
androbustness. impressionofthecoronalmagneticfieldstructure.
Thepurposeofthepresentcontributionistoassess some Thereare,however,bothobservationalandtheoreticalar-
ofthemethodspresentlyavailableforcalculatingnon-linear gumentsthatatleastthemagneticfieldpriortoeruptivepro-
force-free fields, with a view to implementing them into a cesses in the corona is not a linear force free (or potential)
generalextrapolationcodewhichalsotakesstereoscopicin- field. If the normal componentor the line-of-sightcompo-
formationinto accountin the extrapolationprocess. A first nentofthe magneticfieldonthe boundaryis giventhecor-
versionofsuchamethodbasedonlinearforcefreeextrapola- respondingpotentialfield is uniquelydeterminedand is ac-
tionhasrecentlybeenproposed(WiegelmannandNeukirch tuallythemagneticfieldwiththelowestmagneticenergyfor
,2002). AsstatedinWiegelmannandNeukirch (2002)and this boundary condition. Since magnetic energy is needed
bearing in mind the remarksabove, it is would be very de- to drive an eruptive process pre-eruptive coronal magnetic
sirable to generalize the extrapolation method to nonlinear fields cannot be expected to be in a potential state. This is
forcefreefields. alsocorroboratedbyobservations(e.g.Hagyard,1990;Hag-
As a first step in this directionwe investigatethe robust- yardetal. , 1990;Falconeretal. , 2002). Similarobserva-
nessoftheoptimizationmethodrecentlyproposedbyWheat- tionsalsoindicatethatitwillusuallybedifficulttomatchthe
landetal. (2000),andcompareitsperformancewithanother boundarymagneticfieldwithalinearforce-freefield. Also,
methodbasedonMHDrelaxation.Thepaperisorganizedas from a more generalpointof view it seems highly unlikely
follows. InSect. 2wediscussthebasicequationsandback- that the complexpre-eruptioncoronalmagneticfield which
ground theory. Section 3 describes various methods which isslowlybutconstantlystressedbychangingboundarycon-
havebeenproposedfornonlinearforcefreefield extrapola- ditions will have the current density distribution of a linear
tion. Outofthesewefocusontheoptimisationmethodand, force-freefield,i.e. aconstantα. A(strongly)localizedcur-
forthecaseofanalyticallygivenboundarydata,ontheMHD rent density in the erupting configuration seems to be a lot
relaxationmethod.Thesemethodsaretestedusingananalyt- morenatural.Anotherargumentwhichissometimesputfor-
icalsolutioninSect. 4.Theoptimisationmethodisalsosub- wardtoruleoutlinearforce-freemagneticfieldsformodels
jectedtotestswithmoregeneralboundaryconditionswhere ofpre-eruptivestatesisthatlinearforcefreefieldsminimize
weaddnoisetothegivenanalyticalboundarydata. Thepa- themagneticenergyundertheassumptionofglobalhelicity
percloseswiththeconclusionsinSect. 5. conservation (e.g., Taylor , 1974), and therefore if helicity
conservation could be assumed to hold for a coronal erup-
tion,alinearforce-freefieldwouldnotbeabletoprovidethe
2 Basicequations: Force-FreeEquilibria
energyfortheeruption.Eventhoughithasbeenproposedto
applyvariantsofTaylor’shypothesistothestateofthecoro-
Forforcefreefieldstheequationstosolveare
nalmagneticfieldonlargerscales(e.g.,HeyvaertsandPriest
j×B = 0 (1) , 1984), it generally does not seem to apply for models of
eruptiveprocesses(e.g.,AmariandLuciani,2000).Further-
∇×B = µ j, (2)
0
more,thereisthepossibilityofhelicitytransportintoandout
∇·B = 0. (3)
oftheregionofinterestandforthesereasons,wedonotwant
Equation(1)impliesthatforforcefreefieldsthecurrentden- toputtoomuchemphasisonthisargument.
sityandthemagneticfieldareparallel,i.e. However, we believe that the other argumentsare strong
enoughtoassumethatgenerallyonehastotakeintoaccount
µ0j=αB. (4) that for the force free magnetic fields of pre-eruptivestates
thevalueofαchangesfromfieldlinetofieldline. Thisau-
Hereαisafunctionofspace.Thisfunctionhastosatisfythe
tomaticallyleadsustoconsidernonlinearforcefreefields.
equation
B·∇α=0 (5)
3 ExtrapolationMethods
which is obtained by taking the divergence of Eq. (4) and
makinguseofEqs. (2)and(3).SubstitutingEq. (4)intoEq.
The problem we discuss in this section is to compute non-
(2)wegetasecondequationwhichdeterminesB:
linear force free magnetic fields in Cartesian coordinatesif
givenallthreecomponentsofthemagneticfieldonthebound-
∇×B=αB (6)
ary, i.e. the photosphere. Various methods have been pro-
Equation (5) implies that α is constant along magnetic posed(seee.g.Amarietal. ,1997;De´moulinetal. ,1997;
fieldlines,butthatitcanvaryacrossthemagneticfield. Ob- McClymontetal. ,1997;Semel,1998;Amarietal. ,2000;
viously,Eq. (5)issolvedbyα = 0implyingthatj = 0by YanandSakurai,2000)ofwhichwewilldiscussthree. For
Eq. (4). Inthiscaseweobtainpotentialfields. Anotherob- the first method we assume that the magnetic field is given
viouspossibilityisαconstant,butnonzero. Thisisthecase on the lower boundary (photosphere), whereas for the sec-
Coronalmagneticfields. 3
ondandthirdmethodweassumethatBisprescribedonthe 3.2 MHD-Relaxation
surfacesofavolumeV.
Anotherpossibilitytocalculatenonlinearforcefreefieldsis
3.1 DirectExtrapolation bythemethodofMHDrelaxation(e.g.,ChoduraandSchlu¨ter
, 1981). The idea is to start with a suitable magnetic field
Directextrapolation(Wuetal. ,1990)isaconceptuallysim- which is not in equilibrium and to relax it into a force free
ple method, in which the Eqs. (6) and (5) are reformulated state. This is done by using the MHD equationsin the fol-
in such way that they can be used to extrapolatethe photo- lowingform:
sphericboundaryconditionsgivenbyavectormagnetogram
νv = (∇×B)×B (14)
intothesolarcorona.
E+v×B = 0 (15)
If we suppose that the boundary condition on the lower
∂B
boundaryis given by B0(x,y,0) then we can calculate the = −∇×E (16)
z-componentofthecurrentdensityatz =0byusingEq.(2) ∂t
∇·B = 0 (17)
∂B ∂B
µ j = y0 − x0. (7) Theequationofmotion(14)hasbeenmodifiedinsuchaway
0 z0
∂x ∂y
that it ensures that the (artificial) velocity field is reduced.
Knowingthez-componentofthephotosphericmagneticfield Equation(15)ensuresthatthemagneticconnectivityremains
(B ),wecanuseEqs. (7)and(4)todetermineα(x,y,0)= unchangedduringthe relaxation. Therelaxationcoefficient
z0
α : canbechoseninsuchwaythatitacceleratestheapproachto
0
theequilibriumstate. Weuse
j
α = z0 (8) 1
0 B ν = |B|2 (18)
z0 µ
We remarkthatspecialcarehasto betakenatphotospheric withµ=constant.Choosingtherelaxationcoefficientνpro-
polarityinversionlines,i.e. linesalongwhichBz0 = 0(see portionalto B2 speeds up the relaxation process in regions
e.g.,Cupermanetal. ,1991). Equation(8)allowsustocal- of weak magnetic field (see Roumeliotis , 1996). Combin-
culate the x and y-component of the current density using ingEqs. (14),(15),(16)and(18)wegetanequationforthe
Eq. (4) evolutionofthemagneticfieldduringtherelaxationprocess:
∂B [(∇×B)×B]×B
jx0 = α0Bx0, (9) =µ∇× . (19)
∂t (cid:18) B2 (cid:19)
j = α B . (10)
y0 0 y0
Thisequationisthensolvednumericallystartingwithagiven
WenowuseEq.(3)andthexandy-componentofEq.(2) initialconditionforB. Equation(19)ensuresthatEq. (17)
to obtain expressionsfor the z-derivativesof all three mag- is satisfied duringthe relaxationif the initialmagneticfield
neticfieldcomponentsintheform satisfiesit.
The difficulty with the relaxation method in the present
∂B ∂B
x0 = j + z0, (11) form is that it cannot be guaranteed that for the boundary
y0
∂z ∂x conditionswe impose and for a given initial magnetic field
∂B ∂B
y0 = z0 −j , (12) (i.e. givenconnectivity),asmoothforce-freeequilibriumex-
x0
∂z ∂y ists to which the system can relax. If such a smooth equi-
∂Bz0 = −∂Bx0 − ∂By0. (13) libriumdoesnotexisttheformationofcurrentsheetsistobe
∂z ∂x ∂y expectedwhichwillleadtonumericaldifficulties.Therefore,
carehastobetakenwhenchoosinganinitialmagneticfield.
Theideaistointegratethissetofequations(numerically)in
Thisclearlylimitstheapplicabilityoftherelaxationmethod,
thedirectionofincreasingz,basicallyrepeatingtheprevious
and further work will be needed to overcome this obstacle.
stepsateachheight.
We therefore only apply the relaxation method for compar-
Unfortunatelyitcanbeshownthattheformulationofthe
ison with the optimizationmethodto cases wherewe know
force free equationsin this way is unstable because it is an
theboundaryconditionstobecompatible.
ill-posedproblem(e.g.,Cupermanetal. ,1990;Amarietal.
One can show, however, that the relaxation method con-
,1997).Inparticularonefindsthatexponentialgrowthofthe
vergesunderlessrestrictiveboundaryconditions. Bymulti-
magneticfieldwithincreasingheightisatypicalbehaviour.
plyingEq. (16)byB,integratingoverthecompletecompu-
Thereasonforthisisthatthemethodtransportsinformation
tationalvolumeV andusingEq. (15),wefindthat
onlyfromthephotosphereupwards. Otherboundarycondi-
d
tions. e.g. at an upperboundary,either at a finite heightor W = [(v·B)B−B2v]·ndS− (j×B)·vdV,(20)
B
atinfinitycannotbetakenintoaccount. Attemptshavebeen dt I∂V ZV
madetoregularizethemethod(e.g.,Cupermanetal. ,1991; where
De´moulin et al. , 1992), but cannot be considered as fully 1
W = B2dV. (21)
successful. B
2Z
V
4 Coronalmagneticfields.
Underline-tyingboundaryconditionswehavev =0onthe Here µ is a constant which can be chosen to speed up the
boundary,andonlythevolumeintegralsurvives. UsingEq. convergenceof the iterationprocess. Equation(19)and the
(14)andthefactthatν >0onecanshowthatinthiscase leadingtermsofEq.(31)areidentical,butEq. (31)contains
d µ additionalterms.
W =− 0 (j×B)2dV <0, (22)
B ForthismethodthevectorfieldBisnotnecessarilysole-
dt Z ν
V
noidalduringthecomputation,butwillbedivergence-freeif
and an equilibrium is achieved if and only if the magnetic
theoptimalstate withL = 0isreached. A disadvantageof
fieldisforcefree.
the method is that it cannot be guaranteedthat this optimal
We emphasizethattheseboundaryconditionsareusually
stateisindeedreachedforagiveninitialfieldandboundary
notcompatiblewithprescribingallthreemagneticfieldcom-
conditions. If this is not the case then the resulting B will
ponentsonthe boundaries,butonlythe normalcomponent.
eitherbenotforcefreeornotsolenoidalorboth.
Therefore,inthefollowingwediscusstherelaxationmethod
onlyfortestcasesinwhichweknowthattheboundarycon-
ditionsarecompatible.Thisisdoneforcomparisonwiththe 4 Tests
othermethodwediscussnow,theoptimizationmethod.
We encodedthe MHD-relaxationmethodand the optimiza-
3.3 Optimization
tion method in one program. The code is written in C and
uses4thorderfinitedifferencesonanequalspacedgrid.The
Recently Wheatland et al. (2000) have proposed an opti-
timeiterationiscomputedwiththemethodofsteepestgradi-
mization method which refines a proposal by Roumeliotis
ent(seee.g.textbookslikeGeigerandKanzow,1999).The
(1996).Inthismethodthefunctional
programisparallelizedwithOpenMP.
L= B−2|(∇×B)×B|2+|∇·B|2 d3V (23)
Z 4.1 Thesemi-analyticaltestfield
V (cid:2) (cid:3)
is minimized. Obviously, L is bounded from below by 0.
To test the reconstruction methods, we try to reconstruct a
Thisboundisattainedifthemagneticfieldsatisfiestheforce
semi-analyticnonlinearforcefreesolutionfoundbyLowand
freeequations
Lou (1990).Wheatlandetal. (2000)haveusedsimilartests.
(∇×B)×B = 0 (24) The main difference between their paper and ours is in the
∇·B = 0. (25) diagnosticquantitiesusedandinthecomparisonwithi)the
relaxationmethodandii)thenoisyboundarydata(whichwe
The variation of L with respect to an iteration parameter t
thinkismorerepresentativeofarealisticsituation).Lowand
leadsto(seeWheatlandetal. (2000)fordetails)
Lou (1990) solved the Grad-Shafranov equation for axis-
1 dL ∂B ∂B
=− ·Fd3V − ·Gd2S (26) symmetric force free fields in sphericalcoordinatesr, θ, φ.
2 dt Z ∂t Z ∂t
V S For axis-symmetry the magnetic field can be written in the
where
form
F = ∇×(Ω×B)−Ω×(∇×B) 1 1 ∂A ∂A
−∇(Ω·B)+Ω(∇·B)+Ω2 B (27) B= rsinθ (cid:18)r ∂θer− ∂reθ+Qeφ(cid:19) (32)
G = nˆ×(Ω×B)−nˆ(Ω·B) (28) whereAisthefluxfunction,andQrepresentstheφ-component
Ω = B−2 [(∇×B)×B−(∇·B)B] (29) ofB,dependingonlyonA. ThefluxfunctionAsatisfiesthe
Grad-Shafranovequation
The surface term in (26) vanishes for ∂B/∂t = 0 1 on the
boundaryandLdecreasesfor ∂2A 1−µ2 ∂2A dQ
+ +Q =0 (33)
∂B ∂r2 r2 ∂µ2 dA
=µF (30)
∂t whereµ=cosθ.Amongothers,LowandLou (1990)derive
Thisleadstoaniterationprocedureforthemagneticfield solutionsfor
whichisbasedontheequation
dQ
∂B [(∇×B)×B]×B =α=constant (34)
= µ∇× dA
∂t (cid:18) B2 (cid:19)
bylookingforseparablesolutionsoftheform
+µ −Ω×(∇×B)−∇(Ω·B) P(µ)
n A(r,θ)= . (35)
rn
+Ω(∇·B)+Ω2 B (31)
o The solutions are axis-symmetric in spherical coordinates
1Thisconditionmakesitnecessarythatallthreecomponentsofthemag- withapointsourceattheorigin,butifusedfortestingforce
neticfieldhavetobeprescribedonthesixboundariesofthecomputational freecodesinCartesiangeometrythesymmetryisnolonger
box.Insection4.2.3wediscusshowthesurfaceintegralin(26)canbeused obviousafteratranslationwhichplacesthepointsourceout-
toupdatethelateral andtopboundaryconditions duringtheoptimization
sidethecomputationaldomainandarotationofthesymme-
process. Underthisconditiononlythebottomboundaryisprescribedtime
tryaxiswithrespecttotheCartesiancoordinateaxis.
independentlywiththehelpofphotosphericvectormagnetograms.
Coronalmagneticfields. 5
4.2 Testruns
Table 1. Details ofruns toreconstruct aLow andLou (1990)solution.
ThevaluesoftheparametersusedbyLowandLouarel = 0.3andΦ =
Forthe testwe havefirst calculatedBz(x,y,0)onthe pho- π/4.Thefirstcolumncontainstheusedgridsizeandcommentswhichcode
tosphericboundaryz = 0 fortheappropriateLowandLou (relaxation oroptimization) and boundary conditions have been used. If
solution.UsingthisB asaboundaryconditionwecalculate notspecifiedtheboundaryconditionshavebeenextractedfromtheanalytic
z
solution. We specify if noise is added to the boundary conditions. The
thepotentialfieldcorrespondingtothisboundaryconditions
remarkPot-boundarymeansthatthelateralandtopboundaryconditionsare
inside our computational box. The potential field is com-
notextractedfromtheanalyticsolution,butprescribedasaglobalpotential
putedwithhelpofamethoddevelopedbySeehafer (1978). field. The second column contains the iteration step. The third column
This method gives the componentsof the magnetic field in containsthevalueofthefunctionalL,thefourthcolumntheLorentzforce
terms of a Fourier series. The observed magnetogram (or (averaged over the numerical grid) and the last column the relative error
comparedwiththeanalyticsolution.
here the extractedmagnetogramfromthe Low and Lou so-
lution) covers a rectangular region extending from 0 to Lx nx×ny×nz Step [T2Lm] [n|Nj×mB−|3] RelativeError
in x and 0 to L in y is artificially extendedonto a rectan- 40×40×20 discret.error 3.2 2.0 Reference
y Start 0 46573 264 0.42
gularregioncovering−L to L and−L to L bytaking
x x y y Relaxation 50 3092 38.2 0.26
anantisymmetricmirrorimageoftheoriginalmagnetogram 500 769 7.5 0.095
in the extended region, i.e. Bz(−x,y) = −Bz(x,y) and 5,000 25.0 1.4 0.0016
B (x,−y) = −B (x,y). We use a Fast Fourier Transfor- Optimization 50 832 44.1 0.19
z z
mation (FFT) scheme (see Alissandrakis (1981)) to deter- 500 23.1 7.0 0.03
5,000 3.15 2.05 0.0001
minethecoefficientsoftheFourierseries. Formoredetails
1%noise 5,000 9.5 3.7 0.00017
regardingthismethodseeSeehafer (1978,1982). 10%noise 5,000 645.0 29.2 0.0055
Thepotentialfieldisusedasstartingfieldinsidethecom- 20%noise 5,000 2584 58.6 0.028
putational box but the Low and Lou B field is imposed on Pot-boundary 500 468 9.1 0.39
the boundaries. We then use either the MHD relaxation or Pot-boundary 5,000 458 6.3 0.37
the optimization method to calculate the correct force free 80×80×40 discret.error 0.6 0.65 Reference
Start 0 109203 297 0.48
equilibriumfield. Duringthecomputationswecalculatethe
Optimization 50 1590 52 0.26
quantities L/[T2m] (for both relaxation and the optimiza- 500 81 12 0.079
tion method), the absolute value of the Lorentz force |J× 5,000 0.79 1.02 0.0013
B|/[nN m−3](averagedoverthenumericalgrid),thevalue 25,000 0.6 0.65 0.000003
of|∇·B|/[mGaussm−1](averagedoverthenumericalgrid), 10%noise 25,000 1876 37.6 0.0078
Pot-boundary 1,000 563 12.8 0.42
andthedifferencebetweenthenumericalmagneticfieldand
Pot-boundary 10,000 529 5.3 0.36
the known analytical solution |B(t)−Bana|2 (averaged over
|Bana|2
the numerical grid) at each time step. In section 4.2.1 and
4.2.2allcomponentsofthemagneticfieldarefixedonallsix
canbeseen.
boundaries. Insection4.2.3onlythebottomboundarycon-
ToquantifythisstatementweshowinFig.2acomparison
ditionisprescribedtime-independentandthelateralandtop
oftheevolutionofthefourdiagnosticquantitiesforthetwo
boundariesareupdatedduringthe iteration. Thedetails are
methodsduringthecomputation. Forbothmethodsallfour
describedin4.2.3.
diagnosticquantitiesdecreaseduringthecomputationbutaf-
4.2.1 Standardtests ter 5,000 steps the optimization method shows significant
smaller valuesforthe quantitiesL, ∇·B and thecompari-
InFig. 1weshowthree-dimensionalplotsofselectedmag- sonwiththeanalyticsolution.Thereforewecanstatethatin
netic field linesfor the Low and Lou solution, the potential thiscasetheoptimizationmethodseemstobeslightlymore
fieldcalculatedbytakingtheLowandLouB onthelower efficientthantheMHD relaxationmethod. Thisiscorrobo-
z
boundary,the field ofthe MHD relaxationmethodafter50, rated by a look a Table 1 in which we summarize the main
500and5,000relaxationtimesteps, andthesameplotsfor resultsofthevarioustestrunswehavemade. Thefirstrow
the optimization method. In these runs we used a grid size shows the discretisation errors for L and the Lorentz force
of40×40×20. Thecolourcodingofthebottomboundary if the knownsolution is discretised ona 40×40×20-grid
indicatestheB distributiononthatboundary. and usedto calculate the valuesof thisquantities. Thesec-
z
It can bee seen that the potential field which is used as ondrowcontainsthevaluesofL, theLorentzforceandthe
startingfieldofthecomputationsforbothmethodsisclearly relativeerroraftertheinteriorgridpointshavebeenreplaced
differentfrom the Low and Lou test field. The state of the bythepotentialfieldcalculatedfromthephotosphericBz of
system after 50 steps still shows some resemblance to the theLowandLousolution. Therelativelylargevaluesofthe
initialpotentialfieldforbothmethodsbutthefieldlineshave threediagnosticquantitiesshowthedeviationfromtheequi-
started to evolve away from the potential field. One can librium.
also see small differences between the two methods. Af- Thenexttworowsshowhowthethreediagnosticquanti-
ter5,000stepsnoobviousdifferencesbetweenthemagnetic tiesevolveduringtheMHDrelaxationmethodforthe40×
fieldreachedwitheithermethodandtheLowandLoufield 40×20-grid. One can see that after 5,000 steps the value
6 Coronalmagneticfields.
Fig.1.Toprow,leftpanel:AsetofselectedfieldlinesoftheLowandLousolution.Toprow,rightpanel:Thesamefieldlinesforthecorrespondingpotential
fieldcalculatedasdescribedinthemaintext. Thedifferencebetweenthetwofieldsisobvious. Middlerow: Thefieldlinesafter50(left)and5,000(right)
stepsoftheMHDrelaxationmethod. Bottomrow: Thefieldlinesafter50(left)and5,000(right)stepsoftheoptimizationmethod. Theboxdrawnshows
thespatialextensionofthenumericalbox. Thecolourcodingatthebottomoftheboxesshowsthenormalcomponentofthephotosphericmagneticfieldin
Gauss.
Coronalmagneticfields. 7
ofLhasdroppedbyalmostthreeordersofmagnitude,butis
stilloneorderofmagnitudeabovethevaluecalculatedwith
thediscretisedexactsolution.TheLorentzforcehasdropped
tothelevelofthedicretisationerror,andtherelativeerrorhas
droppedmorethantwoordersofmagnitude.
We are giving the same quantities for the optimization
methodinthefollowingthreerows. Itcanbeseenthatespe-
ciallythevaluesofLarealwayswaybelowthecorrespond-
ing values of the MHD relaxation method. This is not sur-
prisingastheoptimizationmethodreliesontheminimization
ofthefunctionalL. ThefinalvalueofLisactuallybelowthe
discretisationerror. Itcanalsobeseenthattherelativeerror
after 5,000 iterations is more than one order of magnitude
smallerthanthecorrespondingvalueoftheMHDrelaxation
method.
Wehavemadeasteptowardscheckingnumericalconver-
gencefortheoptimizationmethodbyrepeatingthecalcula-
tiononagridwith doubledresolution(80×80×40). The
resultsareshowninthelowerpartofTable1. Thefirstthing
tonoticeisthatthediscretisationerrorisalmostanorderof
magnitude smaller than for the previousgrid. For this grid
25,000 iteration steps have been carried out, and after that
thevaluesofLandtheLorentzforcehavereachedthelevel
ofthediscretisationerror.Therelativeerrorismorethantwo
ordersofmagnitudesmallerthanforthecoarsergrid.
4.2.2 Effectofaddingnoise
The previous calculations have been carried out under the
assumption that the magnetic field on the boundary of the
computationalboxisknownexactly. Suchanidealizedsitu-
ationwillnotbefoundwhenrealdataareused.Vectormag-
netogramswillhavefinite resolutionandsufferfromobser-
vational uncertainties making the reconstruction potentially
moredifficult. Especiallytheoptimizationmethodhasbeen
proposed in view of coping with these difficulties (Wheat-
landetal. ,2000).
To investigate how the optimization method works if the
boundary conditions are not given by the exact analytical
field, but to keep control over the amount of uncertainty,
we have carried out test runs with the optimization method
adding random noise with to the boundary conditions. We
add the noise by multiplyingthe exactboundaryconditions
withanumber1+δwhereδisarandomnumberintherange
−n ≤δ ≤n andn isthenoiselevel.
l l l
Fig. 2. Evolution of the diagnostic quantities. Top panel: functional L
To study the effect of the noise we have done runs with
(definedinEq. (23)),secondpanel: Lorentzforce|j×B|(averagedover
differentamplitudesofnoiseonthe40×40×20grid. The
thenumericalgrid),thirdpanel:|∇·B|(averagedoverthenumericalgrid)
bottom panel: relative error |B(t)−Bana|2 (averagedoverthenumerical evolutionofLandtheLorentzforcewiththenumbersofit-
|Bana|2 eration for variousnoise levels are shown in Fig. 3. It can
grid). The quantities are drawn for the MHD relaxation method (dotted
lines)andfortheoptimizationmethod(solidlines). Thegridsizeusedwas clearlybeseenthatthemethodconvergeslessandlesswell
40×40×20. withincreasingnoise.Fornoiselevelsof10%and20%,the
correspondingvaluesafter5,000stepsforthe40×40×20
gridaregiveninTable1. Wehavealsocarriedoutarunon
the80×80×40gridwithanoiselevelof10%. Forthisrun
thevaluesofthediagnosticquantitiesafter25,000stepsare
stillhigherthanthevaluesofthecorrespondingquantitieson
thecoarsergridafter5,000steps. Wehavetoconcludethat
8 Coronalmagneticfields.
wellposedproblem. Apopularwellposedboundarycondi-
tion for non linear force free fields is to impose the normal
component B of the magnetic field and the normal com-
n
ponentofthecurrentdensityj inregionsforapositive(or
n
negative)B . (seeAly,1989, regardingthecompatibilityof
n
photospheric vector magnetograph data). These conditions
have been derived under the very strict constraint of a flux
balanced magnetogram where all magnetic flux is closed.
Thisimpliesthateachmagneticfieldlinestartingatonepoint
onthephotospherealsohastoendonthephotosphere(ina
regionofoppositemagneticflux,butwiththesamevalueof
α). Suchstrictlyisolatedactiveregionsseemtobeatoore-
strictiveconstraintforrealvectormagnetograms. Realvec-
tormagnetogramsprovideall three componentsof the mag-
neticfieldforeithersignofB andthusalsoj onthecom-
n n
pletebottomboundary. Itseemsreasonabletoimposethese
observed data on the bottom boundary. The price we have
to pay is that field lines starting on the photosphere might
pass the lateraland topboundary. 2 Consequentlywe need
toupdatethelateralandtopboundaryduringtheiteration.
Fig.3. EvolutionofLand|j×B|fortheoptimization codeforvarious
Asafirststeptowardsconsistentboundaryconditionsfor
levels of boundary value noise. The lowest lines correspond to an ideal
vectormagnetogramandthehigherlinestovectormagnetogramswith1%, the optimizationmethod, we choosea potentialfield on the
10%and20%noise,respectively.Thegridsizeusedhereis40×40×20. lateral and top boundaries. Firstly, a potential field can be
Theconvergenceofthemethodgetssmallerwithincreasingnoiselevel. easily computedjust from the measurednormalcomponent
of the photospheric magnetic field and secondly, it can be
assumedthatthesolarmagneticfieldisreasonablywellap-
evena10%uncertaintyintheboundaryconditionscouldaf- proximatedbya potentialfield outsideactiveregions. Here
fecttheconvergenceoftheoptimizationmethodquitebadly. we compute the potential field on the boundaryas a global
Themainproblemisthatwithnoiseandfixedboundarycon- potential field computed from the B distribution at z = 0
z
ditionsonallsixboundariestheboundaryconditionareover- alone.Table1showsthatthevalueofLforthe40×40×20
imposedandwillnolongerbecompatible. Wediscusshow griddropstoaslightlylowervaluethanforthecorrespond-
thisproblemcanbesolvedinthenextsection. ingrunwith10%noise,whiletheremaining|j×B|forces
areafactorofthreelower. Therelativeerrorcomparedwith
4.2.3 The problem of the lateral and top boundary condi- theexactsolutionisquitehighhere,whichisnowonderbe-
tions causethemagneticfieldisforcedtostaypotentialclosetothe
side and top boundaries. One mightnotice that the error in
Untilnowtherunshavebeencarriedoutundertheassump- theforcesisafactorof3largerthanthediscretisationerror,
tionthatthemagneticfieldboundaryconditionsarewellknown whileLismorethantwoordersofmagnitudelargerthanthe
onall six boundariesofthe box. Unfortunatelyrealvector- discretisation error. The reason is that an inconsistencyoc-
magnetogramsonly provideinformationregardingthe pho- cursattheedgesoftheboxbetweenthephotosphericbound-
tospheric magnetic field. The top and side boundary con- aryandthesideboundaries.TheLowandLousolutionisnot
ditions are unknown and have to be fixed somehow. Here potentialonthephotosphereclosetothesideboundaries,but
we want to investigate the influence of the choice of these thechosensideboundaryconditionsarepotential.
boundaryconditions. An additionalproblemoccursby fix- To overcome this difficulty we cannot fix the lateral and
ingthevectormagneticfieldonallsixboundaries(asdonein top boundaries during the iteration, but have to relax also
the previoussections) because one overimposesthe bound- these boundaries. As the iteration equation (30) has been
ary conditions. For the test cases in the previous sections derived under the condition ∂B = 0 on all boundaries we
∂t
the boundary conditions have been extracted from an ana- havetoextendtheoptimizationapproach.Ifweallow ∂B 6=
lyticsolutionandconsequentlytheboundariesareautomati- ∂t
0onsomeboundariesthesurfaceterminequation(26)does
callycompatible(thisimpliesthattherearenoproblemswith
over-imposed boundary conditions). Section 4.2.2 showed 2Ifthefluxisnotbalancedintheregion,itimpliesthatpartoftheflux
that this compatibility mightget lost if the boundariescon- distributionismissingfromthelimited-sizeofthemagnetogramorthatthe
tain noise. To overcome this difficulty and to take into ac- calibrationofthemagnetographisbad.Inbothcases,anymagneticextrap-
count that measured data are only available for the bottom olationmethodwillonlyderiveanapproximatefield.Unfortunatelythesize
boundary we describe how the strict condition ∂B = 0 on of observed vectormagnetograms is limited and a real magnetogram will
∂t usuallynotbeexactlyfluxbalanced.Fieldlinespassingtheboundaryofthe
allboundariescanbeavoidedwithintheoptimizationproce-
computationalboxcanbeeitheropenfieldlinesorcloseonthephotosphere
dure. First we havetoimposethe boundaryconditionsasa outsidetheobservedmagnetogram.
Coronalmagneticfields. 9
not necessarily vanish. It is straightforward to extend the aknownsemi-analyticsolution(LowandLou,1990)froma
iterationby givennon-equilibriuminitial condition. Both methodscon-
vergetotheexactsolution,buttheoptimizationmethodhas
∂B
=µG (36) higheraccuracy.TheMHDrelaxationmethodisonlyapplied
∂t tothiscasesinceinitspresentnumericalimplementation,the
on the open boundaries. Equation (36) changesthe bound- boundaryconditionscouldgiverisetoinconsistencies.
aryvaluesinsuchwaythatLdecreases.Inthefollowingwe To simulate a more realistic situation which is closer to
choose the lateral and top boundariesas open (theyare ini- workingwithobservationaldatafromvectormagnetograms
tializedwithapotentialfield)andapplyEq. (36)here. The we added noise to the boundary conditions. We have only
bottom boundaryremainstime independentduring the iter- appliedtheoptimizationmethodtothistypeofproblemand
ation. Let us remark that one might as well use the above found that already a relatively small noise level of 10 %
described photospheric boundary condition for isolated ac- canaffecttheconvergenceofthemethodquiteconsiderably.
tive regions(j fixed on the photosphereonly in regionsof Moreworkisneededtosee howthisdifficultycanbeover-
n
positiveB ). TodosoonehastoapplyEq. (36)alsoforthe come. Wealsointendtogeneralizetherelaxationmethodin
n
transversalmagneticfieldonthephotosphereinregionswith suchwaythatisableofcopingwithmorerealisticboundary
negativeB . data.
n
Ifwedonotkeepthelateralandtopboundaryconditions Theultimateaimforthefutureistoapplythesemethodsto
fixed,butallowthemtorelaxwithhelpofEq. (36),Ldrops datafromvectormagnetograms. Unfortunatelyvectormag-
to36.2whichisonlyoneorderofmagnitudeabovethedis- netogramsarelessaccuratethanline-of-sightmagnetograms
cretisationerror. Thecorrespondingvalueofthe totalforce (e.g. fromMDIonSOHO).Therefore,inthelightofourre-
decreasesslightly(|j×B|=5.3)comparedtofixedbound- sultsinthecaseswherenoisewasaddedtotheboundarycon-
ary conditions. For observed vectormagnetogramsit might ditionsitwillbeinterestingtoseewhetherandhowquickly
beusefultochooseasufficientlylargeareaaroundanactive the methodsconverge. We also planto use stereoscopicin-
region,i.e. withthesideboundariesrelativelyfarawayfrom formation as a further constraint in the reconstruction pro-
thenon-potentialactiveregion,sothatthemagneticfieldcan cess. This will become especially important for the analy-
be assumed to be approximately potential close to the side sis of data from the STEREO mission. A method which is
boundaries. We remark that it might be hard to find cor- based onlinearforcefree fieldshasrecentlybeenproposed
responding vector magnetogramsas these instruments (e.g. by Wiegelmannand Neukirch (2002), but linear force free
IVMinHawaii)donotobservethecompletesolardisk,but fieldsseemtobetoorestrictivetodescribecoronalphenom-
onlyrestricted areas. We havealso carriedouta runonthe enaappropriately.Therefore,thenaturalnextstepwillbethe
80×80×40gridwithpotentialfieldboundaryconditions. useofnonlinearforcefreefieldsinsuchacode.
Forthisruna stationarystate isreachedafter10,000steps,
and the forces are less than one order of magnitude higher Acknowledgements. TheauthorsthankBerndInhesterandHardiPeterfor
thanthediscretisationerror,whileLisnearlythreeordersof useful discussions. TWacknowledges supportbytheEuropean Commu-
nitysHumanPotentialProgrammethroughaMarie-Curie-Fellowship. TN
magnitudeabovethediscretisationerror.Pleasenotethatthe
thanksPPARCforsupportbyanAdvancedFellowship. Thecomputations
absoluteerrorinthe|j×B|hasthesameorderofmagnitude
havebeencarriedoutonthePPARCfundedMHDClusterinSt. Andrews.
hereasforthecorresponding40×40×20run. Ifweallow WethankNorbertSeehaferandanunknownrefereeforeusefulcomments.
relaxationonthesideandtopboundaries,Ldropsto96and
|j×B|=8.7isslightlyhigherthanforfixedsideboundary
conditions. Weconcludethattheuncertaintyofthesideand References
topboundaryconditionsaffectstheconvergenceoftheopti-
mizationmethod.Itsinfluencecanbereducedsignificantlyif Alissandrakis, C.E.: Onthecomputation ofaconstant-α force-free mag-
neticfield,Astron.Astrophys.,100,197-200,1981.
weallowforarelaxationoftheseboundariesandonlyfixthe
Aly,J.J:Onthereconstructionofthenonlinearforce-freecoronalmagnetic
photosphericboundaryconditions,butfurtherimprovements
fieldfromboundarydata,SolarPhys.,120,19-48,1989.
wouldstillbewelcome. Amari, T. and Luciani, J. F.: Helicity redistribution during relaxation of
astrophysicalplasmas,Phys.Rev.Lett.,84,1196-1199,2000.
Amari, T., Aly, J. J., Luciani, J. F., Boulmezaoud, T. Z., and Mikic, Z.:
5 Conclusions Reconstructingthesolarcoronalmagneticfieldasaforce-freemagnetic
field,SolarPhys.,174,129-149,1997.
Amari,T.,Boulmezaoud,T.Z.,andMikic,Z.:Aniterativemethodforthe
In the present paper we have presented an assessment of
reconstructionofthesolarmagneticfield.I.Methodforregularsolutions,
thepropertiesofanoptimizationmethodwhichhasrecently
Astron.Astrophys.,350,1051-1059,2000.
beenproposedforthecalculationofnonlinearforce-freefields Chiu,Y.T.andHilton,H.H.:ExactGreen’sfunctionmethodofsolarforce-
fromgivenboundarydata(Wheatlandetal.,2000).Onepart freemagnetic-fieldcomputationswithconstantalpha.I.Theoryandba-
oftheassessmentwasa comparisontheperformanceofthe sictestcases,Astrophys.J.,212,873-885,1977.
Chodura,R.,andSchlu¨ter,A.:A3-DcodeforMHDequilibriumandstabil-
optimizationmethodwiththeperformanceofanMHDrelax-
ity,J.Comp.Phys.,41,68-88,1981.
ationmethod.Forbothmethodsnewparallelizedcodeshave
Cuperman,S.,Ofman,L.,andSemel,M.:Determinationofforce-freemag-
beendeveloped. Bothmethodshavebeenappliedtofinding neticfieldsabovethephotosphereusingthree-componentboundarycon-