Relativistic Force-Free Electrodynamic Simulations of Neutron Star Magnetospheres

0.6 MB·English

by  Jonathan C. McKinney

Mon.Not.R.Astron.Soc.000,1–5(2006) Printed5February2008 (MNLATEXstylefilev2.2)
Relativistic Force-Free Electrodynamic Simulations of
Neutron Star Magnetospheres
Jonathan C. McKinney⋆
Institute forTheory and Computation, Harvard-Smithsonian Centerfor Astrophysics, 60 Garden Street, MS51, Cambridge, MA 02138, USA
Accepted 2006January18.Received2006January17;inoriginalform2005December8
6 ABSTRACT
0
0 The luminosity and structure of neutron star magnetospheres are crucial to our
2 understanding of pulsar and plerion emission. A solution found using the force-free

n approximation would be an interesting standard with which any model with more
a physics could be compared. Prior quasi-analytic force-free solutions may not be sta-
J ble,whilepriortime-dependentmagnetohydrodynamicmodelsusedunphysicalmodel

9 parameters.Weuseatime-dependentrelativisticforce-freeelectrodynamicscodewith
1 no free parameters to find a unique stationary solution for the axisymmetric rotat-
ing pulsar magnetosphere in a Minkowski space-time in the case of no surface cur-

1 rents on the star. The solution is similar to the force-free quasi-analytic solution
v
of Contopoulos et al. (1999) and the numerical magnetohydrodynamic solution of
1
Komissarov (2005). The magnetosphere structure and the usefulness of the classi-
1
cal y-point in the general dissipative regime are discussed. The pulsar luminosity is
4
1 found to be L≈0.99±0.01µ2Ω4⋆/c3 for dipole moment µ and angular frequency Ω⋆.
0
Key words: stars:pulsars:general– stars:winds and outflows – relativity – force-free
6
0
/
h
p
- 1 INTRODUCTION using a newly developed general relativistic force-free elec-
o trodynamics (GRFFE) extension (McKinney 2005d) to
r Neutronstarmagnetospheresaresuspectedtohavelargere- a general relativistic magnetohydrodynamics (GRMHD)
t gionsofspacethatarenearlyforce-free(Goldreich & Julian
s code called HARM (Gammie et al. 2003). The origi-
a 1969; Ruderman & Sutherland 1975). Self-consistent quasi-
nal GRMHD code has been successfully used to study
: analytic solutions exist that describe the force-free (and
v GRMHD models of accretion flows, winds, and jets
i some magnetohydrodynamic (MHD)) parts of the environ- (Gammie et al.2003;Gammie, Shapiro, & McKinney2004;
X ment of such systems (see, e.g., Contopoulos et al. 1999;
McKinney & Gammie 2004; McKinney 2005a,b,c). The
r Goodwin et al. 2004;Gruzinov 2005).Quasi-analyticmeth-
a GRFFE formulation is a simple extension of HARM, so all
ods cannot test the stability of their solutions and there is
the code developed for HARM is immediately useful and
little hope to study thegeneral nonaligned rotator.
has identical convergence properties to the GRMHD ver-
Recently,Komissarov(2005)usedatime-dependentnu-
sion. We have shown that our GRFFE formulation is at
merical code to directly integrate the force-free and MHD
leastasaccurateastheGRFFEformulation byKomissarov
equations of motion. They considered their force-free solu-
(2002,2004)(McKinney2005d).Parabolicspatialinterpola-
tionunphysicalduetouncontrolledfast reconnectioninthe
tionandfourth-orderRunge-Kuttatemporalintegrationare
current sheet that led to closed field lines far beyond the
usedtoimproveaccuracy.ComparedtoKomissarov (2005),
light cylinder. Their MHD version had no such defect, but
our force-free scheme is improved by significantly lowering
theywereforcedtointroduceunphysicalevolutionequations
the reconnection rate in current sheets by forcing numeri-
and parameters (see theirsection 4.2) in order to avoid nu-
cally induced velocities into current sheets to vanish. This
merical errors and strong unphysical features. They found
eliminates all the difficulties encountered by Komissarov
a stationary MHD solution that may be unique, but they
(2005).
were uncertain whether their solution depended upon the
unphysicalmodel parameters. Section2discussesthesolutionoftheneutronstarmag-
We investigate axisymmetric neutron star magneto- netosphere. Section 3 summarizes the results of the paper.
spheres by integrating the force-free equations of motion Thenotation follows Misner et al. (1973) and thesignature
of the metric is −+++. Tensor components are given in a
coordinate basis. The components of the tensors of interest
⋆ E-mail:[email protected] are given by Fµν for the Faraday tensor, ∗Fµν for the dual
(cid:13)c 2006RAS
2 Jonathan C. McKinney
of the Faraday, and Tµν for the stress-energy tensor. The of Ω⋆ ≈ 0.0216c3/GM such that the light cylinder is at
field angular frequency is ΩF ≡ Ftr/Frφ = Ftθ/Fθφ. The RL=46.3GM/c3.
magnetic field can be written as Bi = ∗Fit. The poloidal In practice, Br and ΩF = Ω⋆ are fixed, while Bθ and
magnetospheric structure is defined by the φ-component Bφ areparabolically extrapolatedintotheneutronstarsur-
of the vector potential (A ≡ Ψ). The current system face from the computational domain. For a stationary, ax-
φ
is defined by the current density (J) and the polar en- isymmetric force-free solution, one can show that the field
closed current (B ≡ ∗F ). The electromagnetic luminos- geometry completely determines the velocity (see equation
φ φt
ity is L ≡ −2π dθTrr2sinθ. See Gammie et al. (2003); 47 in McKinney 2005d). For arbitrary field components Bi
McKinney & GamRθmie(t2004);McKinney(2005d)fordetails. this prescription for vi always leads to a time-like velocity
within thelight “cylinders” (McKinney 2005d).
Thepolaraxisboundaryconditionissuchthattheper-
pendicular fluxes vanish. The outer boundary condition is
2 NEUTRON STAR MAGNETOSPHERE
obtainedbyextrapolatingthefieldintotheboundaryzones
Here the implementation details and the solution to the and setting the velocity as described above, although the
neutron star magnetosphere are discussed. It is most use- outerboundaryischosentobefarawaytoavoidreflections
ful to the community if a similar model to that chosen by backontothesolution.Intheeventthatthereisafluxfrom
Komissarov (2005)isstudiedandcompared.Thiswill show theouterboundaryintothecomputationalgrid,thatfluxis
that one can study pulsar magnetospheres in the force-free set tozero.
limit using time-dependent numerical models. Despite the
unphysical parameters introduced by Komissarov (2005),
2.2 Coordinates
our solutions are similar.
The only additional note is that we and they choose a Thecomputationisperformedonasetofuniformrectangu-
solution with no surface currents on the stellar surface by larcoordinatesdescribedbythevectorfieldx(µ),whereeach
our and their choice of relaxed boundary conditions. This uniformcoordinateisarbitrarilymappedto,e.g.,t,r,θ,φin
corresponds to, e.g., a neutron star with a crust in shear spherical polar coordinates. Apart from the code’s ability
equilibrium(see,e.g.,Rudermanet al.1998).Otherbound- to avoid significant numerical reconnection, the ability to
ary conditions may lead to other solutions. choose an arbitrary θ grid is crucial to obtain an accurate
solution around theequatorial current sheet.
The radial coordinate is chosen tobe
2.1 Initial and Boundary Conditions
r=R +ex(1), (2)
0
AsinKomissarov(2005),themagneticfieldisapproximated
as an aligned dipolar field and the space-time is approxi- whereR0 ischosentoconcentratethegridzonestowardthe
matedasMinkowski.Thepoloidalfieldatt=0(andforall surface (as R0 is increased from 0 to R⋆). An inner radius
time for Br) is set to be theperfect dipole with of Rin = R⋆, an outer radius of Rout = 2315GM/c2 (50
B µ light cylinders as in Komissarov 2005), and R0 = 0.903R⋆
A = poleR3sin2θ= sin2θ, (1) are used.
φ 2r ⋆ r
The θ coordinate is chosen to be
wpohlearremthaegnmetaigcnfieetlidc dstirpeonlegtmhoBmpeolnetainsdµst=ellBarporleaRdi⋆3u/s2Rfo⋆r. a θ=πx(2)+ 12(1−h(r))sin(2πx(2)), (3)
As in Komissarov (2005), the initial conditions are ar-
bitrarily chosen to have a velocity of vi = 0 and Bφ = 0, where x(2) labels an arbitrary uniform coordinate and h(r)
and the proceeding violent non-stationary evolution even- is used to concentrate grid zones toward or away from the
tually relaxestoa steady-statesolution. Thus,any solution equator.Thewidedeadzoneininner-radialregionsandthe
found must be a stable solution. As in Komissarov (2005), currentsheetinouter-radialregionsarebothresolvedusing
themodel is evolved for 55 light cylindertimes. h(r)= 1 + 1atan r−rs (h −h )+h ,(4)
A steady-state solution with no discontinuities or sur- (cid:18)2 π (cid:18) r (cid:19)(cid:19) outer inner inner
0
face currents on the stellar surface is found by choosing
boundaryconditionsdeterminedbyananalysisoftheGrad- where hinner = 1.8, houter = 0.05hinner, r0 = 5GM/c2, and
Shafranovequation(see,e.g.,Bogovalov1997;Beskin1997). rs=15GM/c2.
One is required to specify 2 constraints and to fix the Different resolutions were chosen to test convergence.
magnetic field component perpendicular to the star. As in Intermsofradialvs.θ resolution,weused64×64,80×64,
Komissarov (2005), Eφ ≡ Ftφ = 0, ΩF = Ω⋆, and Br are 80×128,80×256,80×128,160×128,160×256,andfinally
fixed in time, where Ω⋆ is the stellar angular velocity. This for the fiducial model we used 480×256 to reach a similar
assumes the particle acceleration gap is negligible, which is resolutionofKomissarov(2005),whoused496×244fortheir
not generally true (see, e.g., Mestel & Shibata 1994). The final resolution. We find that all quantities have converged
velocity obeys the frozen-in conditions (see equation 46 in towithin 1−5% at our highest resolution, where theexact
McKinney 2005d) in steady-state and axisymmetry. percentage dependson thequantity.
Komissarov (2005) set the radius of the star to be
R⋆ = 0.1RL, where RL = c/ΩF = c/Ω⋆ is the light cylin- 2.3 Magnetosphere Solution
der. For R⋆ = 10km, this means they chose a spin pe-
riod of τ ≈ 2.1ms or Ω⋆ ≈ 0.0207c3/GM for a neutron First,thecodeischeckedthatthesolutionreachesasteady-
star with M = 1.44M⊙. We choose a similar frequency state.Thesolutiondoesreachawell-definedsteady-stateby
(cid:13)c 2006RAS,MNRAS000,1–5
Relativistic Force-Free Electrodynamic Simulations of Neutron Star Magnetospheres 3
Figure1.UppertwolinesshowthevalueofΩF/Ω⋆asafunction Figure2.FluxfunctionAφinaboxwithsize3RL×3RL(RLis
of θ for r = 0.2RL (solid line) and r = 0.7RL (dashed line) for light cylinder) as in figure 3 (upper left panel) of Komissarov
64θzonesand80radialzones.LowerlinesshowΩF/Ω⋆−0.2for (2005). There are 80 contours from Aφ = 0 (polar axis) to
256 θ zones and 480 radial zones. Directly comparable to lower Aφ =µ/R⋆ =9.57µΩ⋆/c (equator on star). A single thick solid
panelsinfigure6ofKomissarov(2005). contour has beenadded fortheclosedfieldlinethattouches the
theoreticallightcylinder.Theverticalsolidlineisthetheoretical
lightcylinderRL.
aboutt∼1000GM/c3,orafterabout20lightcylinderlight
crossing times, out to about r ∼ 1000GM/c3 for all of the
domainexceptwithintheequatorialcurrentsheet.Thepart computational grid out to r ∼1000GM/c3 (large radii still
of the equatorial current sheet that connects to the closed has initial transients) is <3%.
zone undergoes slow, small-amplitude oscillations. The fre- Finally, the region where the stationary solution is ac-
quencyof theseoscillations isproportional tothereconnec- tually force-free is checked, which is tracked by our code
tion rate allowed by our reconnection model, while theam- (McKinney 2005d). Force-free is only violated in the very
plitude is inversely proportional to the reconnection rate. thin current sheet and only in thesheet for r>2.8RL. En-
Baseduponthetimingoftheoscillations,thefiducialmodel ergyisnotstrictlyconservedwhenforce-freeisviolated,but
hasareconnectionperiodof∼10msnearthelightcylinder. the volume occupied by this region is negligible and so en-
In the limit of an unlimited reconnection rate, there are no ergy loss is negligible and does not affect the luminosity as
oscillations and themagnetosphere has alarger closed zone a function of radius described below.
thatextendssomewhatbeyondthelightcylinder,similarto Figure 2 shows the magnetic flux function (Aφ). The
Komissarov (2005). theoretical light cylinder at R = RL very closely follows
Second, the code is checked for proper integration. thetrueouterlight cylinderexceptin theverythin current
Since the neutron star has Ω⋆ = ΩF over the entire sur- sheet, where the light “cylinder” extends slightly radially
face, then in axisymmetry and for a stationary flow, ΩF outwardby≈5%RL.Thismagnetosphericstructureissim-
should be the same for all field lines unless force-free (or ilar to the force-free solution of Contopoulos et al. (1999)
idealMHD)isviolated.Noticethatstrongsuper-corotation (although they have some kinks at the light cylinder) and
(sub-corotation) would artificially enhance (diminish) the theMHD solution of Komissarov (2005).
spindown luminosity and decrease (increase) the size of Figure3showsthestructureofthemagnetospherenear
the closed zone. Figure 1 shows ΩF/Ω⋆ as a function of the slightly dissipative current sheet. In our solutions and
θ for r = 0.2RL and r = 0.7RL, the same locations as in the solutions of Komissarov (2005), there are always some
Komissarov (2005) for their figure 6 (lower 2 panels). The very thin, extended closed field loops at the equator. Us-
plot is for t = 2546.3GM/c3 ≈ 55RL/c for a θ resolution ing our low-resistivity model that keeps the current sheet
of only 64 zones and 80 radial zones (upper lines) and of fromreconnecting,theheightsoftheseextendedclosedfield
256 θ zones by 480 radial zones (lower 2 lines displaced by loops are smaller for increasing resolution, which also more
0.2). Komissarov (2005)’s figureshows that evenat aθ res- sharplydefineswhatisqualitativelyassociated withtheso-
olution of 244 that the peak-to-peak fractional variation is calledy-point.Forgeneraldissipativeapplications,however,
22% near theequator and 45% overall (they haveproblems thelocationofthey-pointisill-definedandcannotbequan-
at the poles). Ourmodel with only a θ resolution of 64 has titatively determined. Notice that while Komissarov (2005)
a peak-to-peak variation of only < 10%, and there are no state that their closed zone reaches all the way to the light
significant artifacts at the poles. Over the entire domain, cylinder, it is clear from their plots and other statements
the spin of the magnetosphere deviates at worst by 5%. At they make that there is no well-defined closed-open tran-
a resolution of 480×256, the variation of ΩF/Ω⋆ over the sition, which is consistent with our results. However, one
(cid:13)c 2006RAS,MNRAS000,1–5
4 Jonathan C. McKinney
with ourresults. Sincethenatureof thedissipation isquite
differentineachcode,thissuggeststhatoursimilarlychosen
Ω⋆ and so the lack of a perfect dipole, rather than dissipa-
tion, is the reason why the first kinked field line appears
insidethetheoretical light cylinder.Otherwiseperhapsit is
intrinsic to any stable solution.
Second,inKomissarov(2005)theydefinethe“y-point”
value of the vector potential to be at the theoretical light
cylinderat θ=π/2, which we findto be
A | ≡Ψ =1.27±0.01µΩ⋆, (8)
φRL light c
which is quite close to the value of Ψseparatrix ≈1.27µΩ⋆/c
given by Gruzinov (2005), similar to Ψopen = 1.23µΩ⋆/c
given in Contopoulos (2005), similar to Ψopen =1.27µΩ⋆/c
given by Timokhin (2005), and similar to Ψ = 1.265±
open
0.005(µΩ⋆/c) given by Komissarov (2005). Notice that the
definition of the y-point value by Komissarov (2005) is
their Ψ , which we defined as Ψ . One should com-
yp light
pare our Ψ to their Ψ , which are the same measure-
light yp
ment. However, while they get 1.37, we get 1.27 in units of
Figure 3. Flux function Aφ for zoomed-in region of Fig- c=Ω⋆ =µ=1.
ure 2 around light cylinder. There are 80 contours from Aφ = Third, as in Komissarov (2005), the vector potential
1.09µΩ⋆/ctoAφ=1.80µΩ⋆/c.Therightverticallinerepresents value continues to drop along the equator until r > 5RL
the theoretical light cylinder. The left vertical line points to the where their value just oscillates around Ψ ≈ 1.265 ±
open
fieldkinkpointattheequator. 0.005(µΩ⋆/c). Wefind that such a measurement gives
A | ≡Ψ =1.226±0.005µΩ⋆, (9)
can identify at least threequantitativemeasures, similar to φopen open c
Komissarov (2005). which is quite close to the value of Ψopen = 1.23µΩ⋆/c
First, thefirst closed field line to develop a kink at the given in Contopoulos (2005) and slightly different than the
equator nearthe light cylinderhas an angle of value given by Komissarov (2005). That the value given by
θ ≈76◦±3◦ (5) Komissarov(2005)islargermeanstheyhavemoreopenflux.
kink
Thus,oneexpectstheirluminositytobelarger,whichisthe
between thefieldline andequatorat thekink.This issimi- case.
lartothesolutionfoundbyGruzinov(2005).Astheyfound, Ofsignificantinterestistheluminositycoefficientofthe
we find that this result is independent of Ω⋆ since the so- dipolar approximation. Wefind that theluminosity is
lution is smooth through the surface. For the model being
µ2Ω4
presented, the intersection of the kinked field line with the L≈0.99±0.01 ⋆. (10)
stellar surface is at θ ≈ 70◦ between a horizontal line c3
kink,⋆
and the field line at the stellar surface. The kink intersects Thisresultfor LagreeswiththeresultbyGruzinov(2005),
theequator at despite the detailed differences in the location of the kink-
point, y-point, or the amount of open flux vs. closed flux.
rkink≈43.08±0.05GM/c2 =0.931±0.001RL, (6) Sinceenergyisexplicitlyconservedinournumericalmethod,
the energy flux through each radius follows this same for-
slightly within thetheoretical light cylinder.Thekinkposi-
mula for the entire region that has reached a steady-state
tion oscillates in time by ≈±0.01RL, where theprior error
estimate is smaller since a mean position in time was iso- anddoesnotviolateforce-free, unlikeinKomissarov (2005)
where energy is not explicitly conserved. Note that for a
lated. At thekink,the vectorpotential has a value
lowerresolutionmodel,e.g.160radialzonesby128θzones,
Aφ|kink≡Ψkink=1.35±0.01µΩc⋆, (7) itshaqtuiLte≈sim1.1il0ar±t0o.0t5heinreusnuilttsbwyitKhoµmi=ssaΩr⋆ov=(2c0=051)., Gwhiviecnh
which is quite close to the value of Ψopen = 1.36µΩ⋆/c that our code with 64 θ zones generated less error in ΩF
given by Contopoulos et al. (1999) and by the value of than even their model with 244 θ zones, then likely their
Ψyp=1.375±0.005µΩ⋆/cgivenbyKomissarov(2005).That solution with 496×244 radial vs. θ zones is simply not as
thesolutionkinksbeforethetheoreticallightcylindercould converged as our solution with 480×256 zones. However,
beduetothelackofaperfectdipoleduetothelightcylinder the agreement between our force-free model and their ideal
being close to the stellar surface, the remaining dissipation MHD model is impressive given the many unphysical pa-
in the model, or could be intrinsic to our stable model. A rameterstheyhadtointroducetokeepthesolutionrealistic
comparisonofthenumericaldissipativey-pointandthethe- and thecode stable (see theirsection 4.2).
oretical dissipationless y-point (see, e.g., Uzdensky 2003) is Figure 4 shows the large-scale structure of the field at
left for futurework. thefinaltimeof55RL/coutto30RL.Thereare40contours
The first kinked field line appears around r ≈ 0.9RL from Aφ = 0 to Aφ = 1.41µΩ⋆/c. The field is essentially
in the fiducial model of Komissarov (2005), and this agrees monopolar.ThisissimilartoKomissarov(2005),althougha
(cid:13)c 2006RAS,MNRAS000,1–5
Relativistic Force-Free Electrodynamic Simulations of Neutron Star Magnetospheres 5
Contopoulos et al. (1999) and is also similar to figure 4 in
Komissarov (2005).
3 CONCLUSIONS
Astationaryforce-freesolutionisfoundfortheneutronstar
magnetosphere with a dipolar surface field in Minkowski
space-time with no stellar surface current. The luminosity
follows the standard dipolar luminosity with a coefficient
determined accurately to be k=0.99±0.01.
The difficulties encountered by Komissarov (2005) in
the force-free regime were avoided, and the unphysical pa-
rameterstheyintroducedintheMHDregimewereavoided.
Their MHD solution is similar to our force-free solution,
which suggests that the solution we and they find is accu-
rate, stable, and may be unique.
ACKNOWLEDGMENTS
ThisresearchwassupportedbyNASA-AstrophysicsTheory
Figure4.FluxfunctionAφ outto30RL.Thereare40contours ProgramgrantNAG5-10780andaHarvardCfAInstitutefor
from0to1.41µΩ⋆/c.
TheoryandComputationfellowship.IthankSergueiKomis-
sarov and Dmitri Uzdenskyfor useful comments.
REFERENCES
Beskin, V. S. 1997, UspekhiFizicheskikh Nauk,40, 659
Bogovalov, S. V.1997, A&A,323, 634
Contopoulos, I.,Kazanas, D.,&Fendt,C.1999, ApJ,511,
351
Contopoulos, I.2005, A&A,442, 579
Gammie,C.F.,McKinney,J.C.,&Ga´borTo´th2003,ApJ,
589, 444
Gammie, C. F., Shapiro, S. L., & McKinney, J. C. 2004,
ApJ,602, 312
Goldreich, P., & Julian, W. H.1969, ApJ, 157, 869
Goodwin, S.P.,Mestel, J., Mestel, L.,& Wright,G. A.E.
2004, MNRAS,349, 213
Gruzinov, A.2005, Physical Review Letters, 94, 021101
Komissarov, S. S.2002, MNRAS,336, 759
Komissarov, S. S.2004, MNRAS,350, 427
Komissarov, S. S. 2005, ArXiv Astrophysics e-prints,
arXiv:astro-ph/0510310
McKinney,J. C., & Gammie, C. F. 2004, ApJ,611, 977
Figure 5.CurrentfunctionBφ vs.fluxfunctionΨ=Aφ.
McKinney,J. C. 2005a, ApJ, 630, L5
McKinney, J. C. 2005b, ArXiv Astrophysics e-prints,
reconnectioneventandplasmoid motiondisruptstheirfield arXiv:astro-ph/0506369
at r∼23RL. McKinney, J. C. 2005c, ArXiv Astrophysics e-prints,
AcontourplotofthepulsarcurrentfunctionB shows arXiv:astro-ph/0506368
φ
that the quantity follows field lines and has a significant McKinney, J. C. 2005d, ArXiv Astrophysics e-prints,
currentincrease across theseparatrix.A plotof thecurrent arXiv:astro-ph/0601410
structure agrees with the classical model of pulsar current Mestel, L., & Shibata, S.1994, MNRAS,271, 621
closure and with Komissarov (2005). To gain quantitative Misner, C. W., Thorne, K. S.,& Wheeler, J. A.1973, San
insight, Figure 5 shows the pulsar current function vs. the Francisco: W.H.Freeman and Co., 1973,
pulsar flux function, which is a similar plot to figure 4 in Ruderman,M., Zhu,T., & Chen, K.1998, ApJ, 492, 267
Komissarov (2005). Wefind that Ruderman,M.A.,&Sutherland,P.G.1975, ApJ,196,51
Timokhin, A. 2005, ArXiv Astrophysics e-prints,
Ψ =1.07µΩ⋆ at B =±1.05µΩ2⋆, (11) arXiv:astro-ph/0507054
max c φ c2 Uzdensky,D. A.2003, ApJ,598, 446
at r = 1.1RL. This is similar to the distribution shown by
(cid:13)c 2006RAS,MNRAS000,1–5