Probing millisecond pulsar emission geometry using light curves from the Fermi/Large Area Telescope

TL;DR: In this article, the authors present results from 3D emission modeling, including the Special Relativistic effects of aberration and time-of-flight delays and also rotational sweepback of 13-field lines, in the geometric context of polar cap (PC), slot gap (SG), outer gap (OG), and two-pole caustic (TPC) models.

Abstract: An interesting new high-energy pulsar sub-population is emerging following early discoveries of gamma-ray millisecond pulsars (MSPs) by the Fermi Large Area Telescope (LAT). We present results from 3D emission modeling, including the Special Relativistic effects of aberration and time-of-flight delays and also rotational sweepback of 13-field lines, in the geometric context of polar cap (PC), slot gap (SG), outer gap (OG), and two-pole caustic (TPC) pulsar models. In contrast to the general belief that these very old, rapidly-rotating neutron stars (NSs) should have largely pair-starved magnetospheres due to the absence of significant pair production, we find that most of the light curves are best fit by SG and OG models, which indicates the presence of narrow accelerating gaps limited by robust pair production -- even in these pulsars with very low spin-down luminosities. The gamma-ray pulse shapes and relative phase lags with respect to the radio pulses point to high-altitude emission being dominant for all geometries. We also find exclusive differentiation of the current gamma-ray MSP population into two MSP sub-classes: light curve shapes and lags across wavebands impose either pair-starved PC (PSPC) or SG / OG-type geometries. In the first case, the radio pulse has a small lag with respect to the single gamma-ray pulse, while the (first) gamma-ray peak usually trails the radio by a large phase offset in the latter case. Finally, we find that the flux correction factor as a function of magnetic inclination and observer angles is typically of order unity for all models. Our calculation of light curves and flux correction factor f(_, _, P) for the case of MSPs is therefore complementary to the "ATLAS paper" of Watters et al. for younger pulsars.

Summary

1. INTRODUCTION

• The field of gamma-ray pulsars has already benefited profoundly from discoveries made during the first year of operation of the Fermi/Large Area Telescope (LAT).
• The pairs screen the accelerating E-field, and limit both the parallel and transverse gap size (Takata et al. 2004).
• The authors calculation of light curves and flux correction factors fΩ(α, ζ, P ) for the case of MSPs is therefore complementary to the work of Watters et al. (2009) which focuses on younger pulsars, although their TPC and OG models include non-zero emission width.

2.1. B-field and SR Effects

• Deutsch (1955) found the solution of the B- and E-fields exterior to a perfectly conducting sphere which rotates in vacuum as an inclined rotator.
• Following Chiang & Romani (1992), Cheng et al. (2000) and Dyks & Rudak (2003), the authors assume constant emissivity along the Blines in the gap regions of the geometric PC, OG, and TPC models (but not for the PSPC model), so that they do not need to include any particular E-field (or calculate dS explicitly) for these.
• The authors calculate the position and direction of photon propagation (assumed to be initially tangent to the local B-line) in the co-rotating frame, and then aberrate this direction using a Lorentz transformation, transforming from the instantaneously co-moving frame to the inertial observer frame (IOF).
• Watters et al. (2009) assume that the emission occurs only along the inner edge of both the TPC and OG gaps (rovc,1 in their notation), and so their radiation width is confined to a single field line and not equal to their gap width (w in their notation).

2.2. Particle Transport and PSPC E-field

• The authors only consider CR losses suffered by electron primaries moving along the B-field lines when modeling the HE emission.
• The authors tested their full solution of E||, which incorporates E (1) || through E(3)|| , for conservation of energy when solving the transport equation (Equation (1)) for relativistic electron primaries.

2.3. Generation of Phaseplots

• (Even though dṄ varies across the PC, the authors assume that it stays constant along B-lines, i.e., that it has no η-dependence.).
• The authors next followed the motion of electron primaries along these lines (Section 2.2), collecting HE radiation and binning as described below.
• (The “weighting factor” W therefore scales the total power to the power radiated in the particular energy bin.).

2.4. Radio Beam Model

• The authors model the radio emission beam using an empirical cone model that has been developed over the years through detailed study of pulse morphology and polarization characteristics of the average-pulse profile.
• (We do not assume a longitudinal extension of the radio emission region, but only use a single emission altitude.)the authors.the authors.
• For the current study, the authors are only interested in pulse shapes and phase shifts between the radio and gamma-ray pulses.
• The authors therefore use relative units for the cone beam luminosity.

2.5. Flux Correction Factor

• Such a flux correction factor (fΩ) is necessarily model dependent, as any observer only sees a small part of the total radiation: that coming from a slice through the emission beam, determined by the line-of-sight ζ .
• Venter (2008) defined the total gamma-ray luminosity using Lγ = Λd2Gobs, (31) with Λ = εΔΩbeam/βobs, ε = βobsGbeam/Gobs, βobs being the duty cycle, ΔΩ beam the average beaming angle, and Gbeam the all-sky total energy flux.

3. LIGHT CURVE DATA

• The authors compare the light curves generated with the different theoretical models (by making constant-ζ cuts through the respective phaseplots of gamma-ray and radio emission) to the light curves of the eight MSPs recently discovered by FermiLAT (Atwood et al.
• In order to reduce the contamination of the gamma-ray signal by the Galactic and extragalactic diffuse emission or nearby sources, and thereby maximize the signal-to-noise ratio, photons were selected in narrow regions of interest, with radii of 0.◦5–1◦ around the pulsar locations.
• Therefore, the light curves shown in Figures 16–19 are biased toward energies above 1 GeV and may not reflect the actual profile shape obtained using larger regions of interest.

4. RESULTS

• Table 1 summarizes some of the properties of the MSPs discovered by Fermi-LAT (Abdo et al. 2009d).
• Differences of the geometric TPC, OG, PC, and also the PSPC models are graphically presented in Figures 7 and 8 (the geometric PC models do not provide particularly good fits to the observed light curves, and the authors will therefore not concentrate on their detailed properties in what follows).
• In the OG model, no emission originates below the null charge surface (note that the null charge surface is at ζ = 90◦ in these plots), so that an observer can only see emission from one magnetic pole, in contrast to the SG/TPC models where an observer sees emission from both poles.
• Double-peaked radio profiles occur for both large α and ζ .
• The authors bestfit model light curves allow us to infer values for α and ζ for each MSP.

5. DISCUSSION AND CONCLUSIONS

• The authors presented results from three-dimensional emission modeling of gamma-ray and radio radiation in the framework of geometric PC, OG, and TPC pulsar models, and also for the full-radiation PSPC model.
• The authors present results complementary to those obtained by Watters et al. (2009) for young pulsars.
• Another possible origin for higher surface fields is the movement of magnetic poles toward the spin axis during the spin-up phase of an MSP (Lamb et al. 2009).
• For the remaining two MSPs, the radio light curves slightly lag the gamma-ray light curves, and these are fit by the PSPC model (and not by the TPC/OG models).
• With more statistics, this effect may possibly become a discriminator between the TPC and OG models.

Figures

Table 2 MSP Model Descriptions

Figure 4. log10 of gain (acceleration) rate γ̇gain (solid line), loss rate γ̇loss (dashdotted line), curvature radius ρc (dash-dot-dotted line), and the Lorentz factor γ (short-dashed line) as a function of normalized radial distance η. We used φpc = 360◦, ξ = 0.7, α = 40◦, and P = 5 ms in this plot.

Figure 10. Sample light curves for a PC2 model with P = 5 ms. (A color version of this figure is available in the online journal.)

Figure 14. Sample light curves for an OG1 model with P = 2 ms. (A color version of this figure is available in the online journal.)

Figure 21. Flux correction factor fΩ(α, ζ, P ) vs. α and ζ for an OG1 model, with panels (a) and (b) for P = 2 ms and P = 5 ms, respectively.

Figure 20. Flux correction factor fΩ(α, ζ, P ) vs. α and ζ for a TPC2 model, with panels (a) and (b) for P = 2 ms and P = 5 ms, respectively. The same color scale is used throughout, and values of fΩ > 4 were set to zero.

Figure 1. P –Ṗ diagram, indicating contours of constant Ėrot (dashed lines) and rotational age (solid lines), as well as pulsars from the ATNF catalog (Manchester et al. 2005). We used values of Ṗ > 0 corrected for the Shklovskii effect (Shklovskii 1970), and removed pulsars in globular clusters. The squares are the eight newly discovered Fermi MSPs (Abdo et al. 2009d). All except PSR J0218+4232 lie below the ICS deathline, and all eight lie below the CR deathline (modeled by Harding & Muslimov 2002).

Figure 23. Flux correction factor fΩ(α, ζ, P ) vs. α and ζ for a PC2 model, with panels (a) and (b) for P = 2 ms and P = 5 ms, respectively.

Figure 22. Flux correction factor fΩ(α, ζ, P ) vs. α and ζ for a PC1 model, with panels (a) and (b) for P = 2 ms and P = 5 ms, respectively.

Figure 9. Sample light curves (black: gamma-ray; gray/magenta: radio at 1.4 GHz) for the PC2 model, with P = 2 ms. The observer angle ζ changes along the columns, and the inclination angle α along the rows. All pulse shape maxima are normalized to unity, and the phase range goes from φ = 0–1 in each case (and similar for subsequent figures).

Figure 13. Sample light curves for a TPC2 model with P = 5 ms. (A color version of this figure is available in the online journal.)

Figure 2. Examples of the final E-field we obtain after matching E(1)|| through E (3) || for different parameters: −E|| vs. log10 of the height above the PC, normalized by the PC radius RPC = (ΩR3/c)1/2. These plots were obtained for P = 5.75 × 10−3 s, Ṗ = 10−20, R = 106 cm, and M = 1.4M . In panel (a), we chose α = 20◦, ξ = 0.3, φpc = 45◦, in panel (b), α = 35◦, ξ = 0.7, φpc = 150◦, and in panel (c), α = 80◦, ξ = 0.8, φpc = 200◦. In the last panel, the final −E|| is negative, so no solution of ηc is obtained. In each case, we label E (1) || through E (3) || (thick solid lines), indicate potential solutions (which vary with ηc) by thin gray (cyan) lines, and the final solution by thick (red) dashed lines. Also, we indicate ηb where we match E (1) || and E (2) || , and ηc where we match E(2)|| and E (3) || , by thin vertical dashed lines. (Although E (3) || does slightly vary with ηc, we only indicate the E (3) || -solution corresponding to the ηc found for the final solution. For panel (c), we show a typical E(3)|| -solution.) (A color version of this figure is available in the online journal.)

Figure 17. Same as Figure 16, but for PSR J0437−4715 (panels (a)–(d), P = 5 ms) and PSR J0613−0200 (panels (e)–(h), P = 3 ms). Panels (a) and (e) are for a TPC2 model with (α, ζ ) = (30◦, 60◦), while (b) and (f) are for an OG1 model with (α, ζ ) = (30◦, 60◦). For the cases where we use the same (α, ζ )-combination for both the TPC and OG fits, we only have a single radio light curve fit.

Figure 12. Sample light curves for a TPC2 model with P = 3 ms. (A color version of this figure is available in the online journal.)

Table 3 Model Fits for α, ζ , and fΩ(α, ζ, P )

Table 1 Parameters of MSPs Discovered by Fermi-LAT (Abdo et al. 2009d)

Figure 6. Example phaseplots of the radio conal beam, for α = 70◦ and at a frequency of 1.4 GHz. Panel (a) is for P= 2 ms and panel (b) for P= 5 ms. Beam and annulus widths become increasingly larger for shorter periods. Increased intensity for the leading part is due to bunching of the B-field lines around the notch. (Note that in this and following phaseplots, the color scales are not the same for the different panels, but are chosen to show the most details for each case.)

Figure 18. Same as Figure 16, but for PSR J0751+1807 (panels (a)–(d), P = 3 ms) and PSR J16142−230 (panels (e)–(h), P = 3 ms). Panel (a) is for a TPC2 model with (α, ζ ) = (50◦, 50◦), (b) for an OG1 model with (α, ζ ) = (50◦, 50◦), (e) for a TPC2 model with (α, ζ ) = (40◦, 80◦), and (f) for an OG1 model with (α, ζ ) = (40◦, 80◦).

Figure 16. Gamma-ray phaseplots (left panels) and observed and fitted gamma-ray and radio light curves (right panels) for PSR J0030+0451 (panels (a)–(d),P = 5 ms) and PSR J0218+4232 (panels (e)–(h), P = 2 ms). Panel (a) is for a TPC1 model with (α, ζ ) = (70◦, 80◦), (b) for an OG1 model with (α, ζ ) = (80◦, 70◦), (e) for a TPC1 model with (α, ζ ) = (60◦, 60◦), and (f) for an OG1 model with (α, ζ ) = (50◦, 70◦). For the gamma-ray light curves (panels (c) and (g)), the histograms represent the Fermi-LAT data (Guillemot 2009), the horizontal dashed line the estimated background level, the dashed (online: magenta) lines are TPC fits, and dash-dotted (online: green) lines are OG fits (see Table 3). For the radio light curves (panels (d) and (h)) the solid (blue) line represents the radio data, while the dashed (magenta) and dash-dotted (green) lines correspond to the same (α, ζ )-combinations as those of the respective TPC and OG fits.

Figure 11. Sample light curves for a TPC1 model with P = 2 ms. (A color version of this figure is available in the online journal.)

Figure 8. Example phaseplots for the PC1 and PC2 models (panels (a) and (b), respectively), for α = 70◦ and P = 5 ms.

Full text

The Astrophysical Journal, 707:800–822, 2009 December 10 doi:10.1088/0004-637X/707/1/800
C

2009. The American Astronomical Society. All rights reserved. Printed in the U.S.A.
PROBING MILLISECOND PULSAR EMISSION GEOMETRY USING LIGHT CURVES FROM THE
FERMI/LARGE AREA TELESCOPE
C. Venter
1,2,5
, A. K. Harding
1
, and L. Guillemot
3,4,6
1
Astrophysics Science Division, NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA
2
Unit for Space Physics, North-West University, Potchefstroom Campus, Private Bag X6001, Potchefstroom 2520, South Africa
3
Universit
´
e de Bordeaux, Centr
´
ed’
´
Etudes Nucl
´
eaires Bordeaux Gradignan, UMR 5797, Gradignan, 33175, France
4
CNRS/IN2P3, Centre d’
´
Etudes Nucl
´
eaires Bordeaux Gradignan, UMR 5797, Gradignan, 33175, France
Received 2009 August 28; accepted 2009 October 27; published 2009 November 25
ABSTRACT
An interesting new high-energy pulsar sub-population is emerging following early discoveries of gamma-ray
millisecond pulsars (MSPs) by the Fermi Large Area Telescope (LAT). We present results from three-dimensional
emission modeling, including the special relativistic effects of aberration and time-of-flight delays and also
rotational sweepback of B-field lines, in the geometric context of polar cap (PC), outer gap (OG), and two-
pole caustic (TPC) pulsar models. In contrast to the general belief that these very old, rapidly rotating neutron
stars (NSs) should have largely pair-starved magnetospheres due to the absence of significant pair production,
we find that most of the light curves are best fit by TPC and OG models, which indicates the presence of
narrow accelerating gaps limited by robust pair production—even in these pulsars with very low spin-down
luminosities. The gamma-ray pulse shapes and relative phase lags with respect to the radio pulses point to
high-altitude emission being dominant for all geometries. We also find exclusive differentiation of the current
gamma-ray MSP population into two MSP sub-classes: light curve shapes and lags across wavebands impose
either pair-starved PC (PSPC) or TPC/OG-type geometries. In the first case, the radio pulse has a small lag with
respect to the single gamma-ray pulse, while the (first) gamma-ray peak usually trails the radio by a large phase
offset in the latter case. Finally, we find that the flux correction factor as a function of magnetic inclination and
observer angles is typically of order unity for all models. Our calculation of light curves and flux correction factor
for the case of MSPs is therefore complementary to the “ATLAS paper” of Watters et al. for younger pulsars.
Key words: acceleration of particles – gamma rays: theory – pulsars: general – radiation mechanisms: non-thermal
– stars: neutron
Online-only material: color figures
1. INTRODUCTION
The field of gamma-ray pulsars has already benefited pro-
foundly from discoveries made during the first year of operation
of the Fermi/Large Area Telescope (LAT). These include detec-
tions of the radio-quiet gamma-ray pulsar inside the supernova
remnant CTA 1 (Abdo et al. 2008), the second gamma-ray mil-
lisecond pulsar (MSP; Abdo et al. 2009a) following the EGRET
4.9σ detection of PSR J0218+4232 (Kuiper et al. 2004), the
six high-confidence EGRET pulsars (Thompson et al. 1999;
Thompson 2004), and discovery of 16 radio-quiet pulsars using
blind searches (Abdo et al. 2009b). In addition, eight MSPs have
now been unveiled (Abdo et al. 2009d; see Table 1), confirm-
ing expectations prior to Fermi’s launch in 2008 June (Harding
et al. 2005; Venter & De Jager 2005b). A Fermi 6 month pulsar
catalog is expected to be released shortly (Abdo et al. 2009e).
AGILE has also reported the discovery of four new gamma-ray
pulsars, and marginal detection of four more (Halpern et al.
2008; Pellizzoni et al. 2009b), in addition to the detection of
four of the EGRET pulsars (Pellizzoni et al. 2009a). Except for
the detection of the Crab at energies above 25 GeV (Aliu et al.
2008), no other pulsed emission from pulsars has as yet been
detected by ground-based Cherenkov telescopes (Schmidt et al.
2005;Albertetal.2007, 2008; Aharonian et al. 2007;F
¨
ussling
et al.
2008; Kildea 2008; Konopelko 2008; Celik et al. 2008;
De los Reyes 2009).
5
NASA Postdoctoral Program Fellow.
6
Now at Max-Planck-Institut f
¨
ur Radioastronomie, Auf dem H
¨
ugel 69,
53121 Bonn, Germany.
Millisecond pulsars (MSPs) are characterized by relatively
short periods P  30 ms and low surface magnetic fields
B
0
∼ 10
8
–10
9
G, and appear in the lower left corner of
the P –
˙
P diagram (with
˙
P being the time-derivative of P;
see Figure 1, where the newly discovered Fermi MSPs are
indicated by squares). MSPs are thought to have been spun-up to
millisecond periods by transfer of mass and angular momentum
from a binary companion during an accretion phase (Alpar et al.
1982). This follows an evolutionary phase of cessation of radio
emission from their mature pulsar progenitors, after these have
spun down to long periods and crossed the “death line” for
radio emission. These “radio-silent” progenitors (Glendenning
& Weber 2000) are thought to reside in the “death valley” of the
P –
˙
P diagram, which lies below the inverse Compton scattering
(ICS) pair death line (Harding & Muslimov 2002).
The standard “recycling scenario” (Bhattacharya & Van den
Heuvel 1991) hypothesizing that MSP birth is connected to low-
mass X-ray binaries (LMXRBs) might have been confirmed re-
cently by the detection of radio pulsations from a nearby MSP in
an LMXRB system, with an optical companion star (Archibald
et al. 2009). Optical observations indicate the presence of an
accretion disk within the past decade, but none today, raising
the possibility that the radio MSP has “turned on” after termina-
tion of recent accretion activity, thus providing a link between
LMXRBs and the birth of radio MSPs.
High-energy (HE) radiation from pulsars has mainly been
explained as originating from two emission regions. Polar cap
(PC) models (Harding et al. 1978; Daugherty & Harding 1982,
1996; Sturner et al. 1995) assume extraction of primaries from
800

No. 1, 2009 PROBING MSP LIGHT CURVES USING FERMI-LAT 801
Figure 1. P –
˙
P diagram, indicating contours of constant
˙
E
rot
(dashed lines)
and rotational age (solid lines), as well as pulsars from the ATNF catalog
(Manchester et al. 2005). We used values of
˙
P>0 corrected for the Shklovskii
effect (Shklovskii 1970), and removed pulsars in globular clusters. The squares
are the eight newly discovered Fermi MSPs (Abdo et al. 2009d). All except
PSR J0218+4232 lie below the ICS deathline, and all eight lie below the CR
deathline (modeled by Harding & Muslimov 2002).
(A color version of this figure is available in the online journal.)
the stellar surface and magnetic pair production of ensuing
HE curvature radiation (CR) or ICS gamma rays, leading to
low-altitude pair formation fronts (PFFs) which screen the
accelerating electric field (Harding & Muslimov 1998, 2001,
2002). These space-charge-limited-flow (SCLF) models have
since been extended to allow for the variation of the CR PFF
altitude across the PC and therefore acceleration of primaries
along the last open magnetic field lines in a slot gap (SG)
scenario (Arons & Scharlemann 1979; Arons 1983; Muslimov &
Harding 2003, 2004a; Harding et al. 2008). The SG results from
the absence of pair creation along these field lines, forming a
narrow acceleration gap that extends from the neutron star (NS)
surface to near the light cylinder. The SG model is thus a possible
physical realization of the two-pole caustic (TPC) geometry
(Dyks & Rudak 2003), developed to explain pulsar HE light
curves. On the other hand, outer gap (OG) models (Cheng et al.
1986a, 1986b, 2000; Chiang & Romani 1992, 1994; Zhang et al.
2004) assume that HE radiation is produced by photon–photon
pair production-induced cascades along the last open field lines
above the null-charge surfaces (
Ω · B = 0, with Ω = 2π/P ),
where the Goldreich–Julian charge density (Goldreich & Julian
1969) changes sign. The pairs screen the accelerating E-field,
and limit both the parallel and transverse gap size (Takata et al.
2004). Classical OG models may be categorized as “one-pole
caustic models,” as the assumed geometry prevents observation
of radiation from gaps (caustics) associated with both magnetic
poles (Harding 2005). More recently, however, Hirotani (2006,
2007) found and applied a two-dimensional, and subsequently a
three-dimensional (Hirotani 2008b) OG solution which extends
toward the NS surface, where a small acceleration field extracts
ions from the stellar surface in an SCLF-regime (see also Takata
et al. 2004, 2006, and in particular Takata et al. 2008 for
application to Vela). Lastly, Takata & Chang (2009) modeled
Geminga using an OG residing between a “critical” B-field line
(perpendicular to the rotational axis at the light cylinder) and
the last open field line.
Current models using the dipole field structure to model MSPs
predict largely unscreened magnetospheres due to the relatively
low B-fields inhibiting copious magnetic pair production. Such
pulsars may be described by a variation of the PC model
(applicable for younger pulsars), which we will refer to as a
“pair-starved polar cap” (PSPC) model (Muslimov & Harding
2004b, 2009; Harding et al. 2005). In a PSPC model, the pair
multiplicity is not high enough to screen the accelerating electric
field, and charges are continually accelerated up to high altitudes
over the full open-field-line region. The formation of a PSPC
“gap” is furthermore naturally understood in the context of an
SG accelerator progressively increasing in size with pulsar age,
which, in the limit of no electric field screening, relaxes to a
PSPC structure.
Several authors have modeled MSP gamma-ray fluxes, spec-
tra, and light curves in both the PSPC (Frackowiak & Rudak
2005a, 2005b; Harding et al. 2005; Venter & De Jager 2005a;
Venter 2008; Zajczyk 2008) and OG (Zhang & Cheng 2003;
Zhang et al. 2007) cases. Collective emission from a popula-
tion of MSPs in globular clusters (Harding et al. 2005; Zhang
et al. 2007; Bednarek & Sitarek 2007; Venter & De Jager 2008;
Venter et al. 2009) and in the Galactic Center (Wang 2006)have
also been considered. Watters et al. (2009) recently calculated
beaming patterns and light curves from a population of canonical
pulsars with spin-down luminosities
˙
E
rot
> 10
34
erg s
−1
using
geometric PC, TPC, and OG models. They obtained predictions
of peak multiplicity, peak separation, and flux correction factor
f
Ω
as functions of magnetic inclination and observer angles α
and ζ , and gap width w. The latter factor f
Ω
is used for con-
verting observed phase-averaged energy flux G
obs
to the total
radiated (gamma-ray) luminosity L
γ
, which is important for
calculating the efficiency of converting
˙
E
rot
into L
γ
. A good ex-
ample is the inference of the conversion efficiencies of globular-
cluster MSPs which may be collectively responsible for the HE
radiation observed from 47 Tucanae by Fermi-LAT (Abdo et al.
2009c).
In this paper, we present results from three-dimensional emis-
sion modeling, including special relativistic (SR) effects of aber-
ration and time-of-flight delays, and rotational sweepback of
B-field lines, in the geometric context of OG, TPC, and PSPC
pulsar models. We study the newly discovered gamma-ray MSP
population (Abdo et al. 2009d), and obtain fits for gamma-ray
and radio light curves. Our calculation of light curves and flux
correction factors f
Ω
(α, ζ, P ) for the case of MSPs is therefore
complementary to the work of Watters et al. (2009) which fo-
cuses on younger pulsars, although our TPC and OG models
include non-zero emission width. Section 2 deals with details
of the various models we have applied. We discuss light curves
from both observational and theoretical perspectives in Sec-
tion 3, and present our results and conclusions in Sections 4
and 5.
2. MODEL DESCRIPTION
2.1. B-field and SR Effects
Deutsch (1955) found the solution of the B- and E-fields
exterior to a perfectly conducting sphere which rotates in
vacuum as an inclined rotator. We assume that this retarded
vacuum dipolar B-field is representative of the magnetospheric
structure, and we use the implementation by Dyks et al. (2004a)
and Dyks & Harding (2004), following earlier work by Romani
& Yadigaroglu (1995), Higgins & Henriksen (1997), Arendt
& Eilek (1998), and Cheng et al. (2000). For this B-field, the

802 VENTER, HARDING, & GUILLEMOT Vol. 707
(a)
(b)
(c)
Figure 2. Examples of the final E-field we obtain after matching E
(1)
||
through
E
(3)
||
for different parameters: −E
||
vs. log
10
of the height above the PC,
normalized by the PC radius R
PC
= (ΩR
3
/c)
1/2
. These plots were obtained
for P = 5.75 × 10
−3
s,
˙
P = 10
−20
, R = 10
6
cm, and M = 1.4 M

. In panel
(a), we chose α = 20
◦
, ξ = 0.3, φ
pc
= 45
◦
, in panel (b), α = 35
◦
, ξ = 0.7,
φ
pc
= 150
◦
, and in panel (c), α = 80
◦
, ξ = 0.8, φ
pc
= 200
◦
. In the last panel,
the final −E
||
is negative, so no solution of η
c
is obtained. In each case, we label
E
(1)
||
through E
(3)
||
(thick solid lines), indicate potential solutions (which vary
with η
c
) by thin gray (cyan) lines, and the final solution by thick (red) dashed
lines. Also, we indicate η
b
where we match E
(1)
||
and E
(2)
||
,andη
c
where we
match E
(2)
||
and E
(3)
||
, by thin vertical dashed lines. (Although E
(3)
||
does slightly
vary with η
c
, we only indicate the E
(3)
||
-solution corresponding to the η
c
found
for the final solution. For panel (c), we show a typical E
(3)
||
-solution.)
(A color version of this figure is available in the online journal.)
PC shape is distorted asymmetrically by rotational sweepback
of field lines. Each field line’s footpoint is labeled by the
open volume coordinates (r
ovc
,l
ovc
) as defined by Dyks et al.
(2004a), with r
ovc
labeling self-similar contours or “rings” (r
ovc
is normalized to the PC radius R
PC
), and l
ovc
giving the arclength
along a ring (analogous to azimuthal angle; also refer to Harding
et al. 2008 for more details).
We calculate the rim of the PC by tracing field lines which
close at the light cylinder back to the stellar surface, and then
divide this PC into rings (see, e.g., Figure 2 of Dyks & Harding
2004) and azimuthal bins, with each surface patch dS associated
with a particular B-field line. We follow primary electrons
moving along each field line, and collect radiation (corrected
for SR-effects) in a phaseplot map (Section 2.3). Following
Chiang & Romani (1992), Cheng et al. (2000) and Dyks &
Rudak (2003), we assume constant emissivity along the B-
lines in the gap regions of the geometric PC, OG, and TPC
models (but not for the PSPC model), so that we do not need
to include any particular E-field (or calculate dS explicitly) for
these. In the case of the PSPC model, we use the approximation
ξ ≈ r
ovc
(with ξ ≡ θ/θ
pc
the normalized polar angle, and
θ
pc
≈ (ΩR/c)
1/2
the PC angle), and include the full E-field up
to high altitudes (Section 2.2).
In addition to the rotational sweepback (retardation) of the B-
lines, we include the effects of aberration and time-of-flight
delays. We calculate the position and direction of photon
propagation (assumed to be initially tangent to the local B-line)
in the co-rotating frame, and then aberrate this direction using a
Lorentz transformation, transforming from the instantaneously
co-moving frame to the inertial observer frame (IOF). Lastly,
we correct the phase at which the photon reaches the observer
for time delays due to the finite speed of light. More details
about calculation of these SR effects may be found in Dyks
et al. (2004b) and Dyks & Harding (2004), following previous
work by, e.g., Morini (1983) and Romani & Yadigaroglu (1995).
We furthermore explicitly use the curvature radius of the B-field
lines as calculated in the IOF, and not in the co-rotating frame,
when performing particle transport calculations (Section 2.2).
Such a model has also recently been applied to the Crab by
Harding et al. (2008).
We have calculated TPC and OG models assuming gaps that
are confined between two B-field lines with footpoints at r
ovc,1
and r
ovc,2
. We therefore activated only a small number of rings
near the rim (r
ovc
∼ 1) with r
ovc
∈ [r
ovc,1
,r
ovc,2
], and binning
radiation from these, assuming constant emissivity over the
emitting volume. For TPC models, we used r
ovc
∈ [0.80, 1.00],
[0.60, 1.00], [0.90, 1.00], [0.95, 1.00], and [1.00, 1.00] (see
Table 2) corresponding to gap widths of w ≡ r
ovc,2
− r
ovc,1
=
0.20, 0.40, 0.10, 0.05, and 0.00. Similarly, we investigated
OG models with r
ovc
∈ [0.90, 0.90], [1.00, 1.00], [0.95, 1.00]
(widths of w = 0.00, 0.00, and 0.05). These widths are smaller
than, e.g., the value of ∼ 0.14 used by Hirotani (2008a).
We did not find good light curve fits for TPC models with
large w. In the case of OG models, one should consider non-
uniform emission when choosing large w, which is beyond
the scope of this paper. The assumption of constant emissivity
in the emitting volume is a simplification, as OG models are
expected to produce the bulk of the gamma-radiation along the
inner edge (r
ovc,inner
<r
ovc
<r
ovc,PFF
)ofthegap(r
ovc,PFF
<
r
ovc
< 1), with r
ovc,PFF
indicating the position of the PFF, and
r
ovc,inner
some smaller radius depending on the radiation surface
thickness (Watters et al. 2009). We lastly modeled the PC and
PSPC cases with r
ovc
∈ [0.00, 1.00] (i.e., the full open-field-
line volume, for both constant emissivity and full radiation
codes). We used 180 colatitude (ζ ) and phase (φ) bins and
individual ring separations of δr
ovc
= 0.005, while collecting
all photons with energies above 100 MeV (in the case of the
PSPC model) when producing phaseplots and subsequent light
curves.
It is important to note a critical difference between the
radiation distribution in our TPC and OG models and that of
Watters et al. (2009). We assume that emission is distributed
uniformly throughout the gaps between r
ovc,1
and r
ovc,2
,sothe
radiation originates from a volume with non-zero width across
field lines and the radiation and gap widths are the same. For the
TPC model, this geometry is similar to that adopted by Dyks
et al. (2004a), although Dyks et al. (2004a) assumed a Gaussian
distribution of emission centered at the gap midpoint while we
simply assume a constant emissivity across the gap, both of
which crudely approximate the radiation pattern expected in
the SG. Watters et al. (2009) assume that the emission occurs
only along the inner edge of both the TPC and OG gaps (r
ovc,1
in our notation), and so their radiation width is confined to
a single field line and not equal to their gap width (w in
their notation). In the case of the OG, the physically realistic
emission pattern would have a non-zero width lying somewhere

No. 1, 2009 PROBING MSP LIGHT CURVES USING FERMI-LAT 803
between infinitely thin and uniform assumptions (see Hirotani
2008b).
2.2. Particle Transport and PSPC E-field
We only consider CR losses suffered by electron primaries
moving along the B-field lines when modeling the HE emission.
In this case, the (single electron) transport equation is given by
(e.g., Sturner 1995; Daugherty & Harding 1996)
˙
E
e
=
˙
E
e,gain
+ ˙γ
CR
m
e
c
2
= eβ
r
cE
||
−
2e
2
c
3ρ
2
c
β
4
r
γ
4
, (1)
where c is the speed of light in vacuum, β
r
= v
r
/c ∼ 1is
the particle velocity, e is the electron charge, γ is the electron
Lorentz factor, ˙γ
CR
m
e
c
2
is the frequency-integrated (total) CR
loss rate per particle (Bulik et al. 2000), ρ
c
is the curvature
radius (as calculated in the IOF; see Section 2.1), and E
||
is the
accelerating E-field parallel to the B-field. The acceleration and
loss terms balance at a particular γ
RR
in the radiation reaction
regime (Luo et al. 2000):
γ
RR
=

3E
||
ρ
2
c
2eβ
3
r

1/4
. (2)
Previous studies (e.g., Venter & De Jager 2005a, 2008;
Frackowiak & Rudak 2005a, 2005b; Harding et al. 2005;
Zajczyk 2008) have used the solutions of Muslimov & Harding
(1997) and Harding & Muslimov (1998) for the PSPC E-field:
E
(1)
||
=−
Φ
0
R

θ
GR
0

2

12κ

s
1
cos α
+6s
2
θ
GR
0
H (1)δ

(1) sin α cos φ
pc

, (3)
E
(2)
||
=−
Φ
0
R

θ
GR
0

2

3
2
κ

η
4
cos α
+
3
8
θ
GR
(η)H (η)δ

(η)ξ sin α cos φ
pc

(1 − ξ
2
), (4)
with
Φ
0
≡
B
0
ΩR
2
c
, (5)
 ≡
2GM
c
2
R
, (6)
θ
GR
(η) ≈

ΩR
c
η
f (η)

1/2
≈ θ
pc
, (7)
s
1
=
∞

i=1
J
0
(k
i
ξ)
k
3
i
J
1
(k
i
)
F
1
(γ
i
(1),η), (8)
s
2
=
∞

i=1
J
1
(
˜
k
i
ξ)
˜
k
3
i
J
2
(
˜
k
i
)
F
1
( ˜γ
i
(1),η), (9)
γ
i
(η) =
k
i
ηθ
GR
(η)(1 − /η)
1/2
, (10)
˜γ
i
(η) =
˜
k
i
ηθ
GR
(η)(1 − /η)
1/2
, (11)
F
1
(γ,η) = 1 − e
−γ (1)(η−1)
, (12)
and k
i
and
˜
k
i
are the positive roots of the Bessel functions J
0
and
J
1
(with k
i+1
>k
i
and
˜
k
i+1
>
˜
k
i
); θ
GR
0
≡ θ
GR
(1); γ (1) may be
γ
i
(1) or ˜γ
i
(1) in the expression for F
1
. The functions H (η), f (η),
and δ

(η) are all of order unity, and are defined in Muslimov &
Tsygan (1992). The first solution E
(1)
||
is valid for η − 1  1,
and E
(2)
||
for θ
GR
0
 η − 1  c/(ΩR); R is the stellar radius,
η = r/R, α is the angle between the rotation and magnetic axes,
φ
pc
is the magnetic azimuthal angle, κ

= 2GI /(c
2
R
3
)isthe
general relativistic (GR) inertial frame-dragging factor (distinct
from the κ(x) function to be defined later), and I is the moment
of inertia.
Muslimov & Harding (2004b) found the solution of E
||
for altitudes close to the light cylinder in the small-angle
approximation (small α, ξ , and high altitude):
E
(3)
||
≈−
3
16

ΩR
c

3
B
0
f (1)

κ


1 −
1
η
3
c

(1 + ξ
2
) cos α
+
1
2

√
η
c
− 1


ΩR
c

1/2
λ(1 + 2ξ
2
)
× ξ sin α cos φ
pc

(1 − ξ
2
), (13)
and λ is defined after Equation (35) of Muslimov & Harding
(2004b). They proposed that one should employ the following
formula to match the last two solutions:
E
||
≈ E
(2)
||
exp[−(η − 1)/(η
c
− 1)] + E
(3)
||
, (14)
with η
c
being a radial parameter to be determined using a
matching procedure. Muslimov & Harding (2004b) estimated
that η
c
∼ 3–4 for MSPs when ξ = θ/θ
GR
0
∼ 0.5.
It is important to include the high-altitude solution E
(3)
||
,
as Fermi results seem to indicate that the HE radiation is
originating in the outer magnetosphere (e.g., Abdo et al. 2009d).
Beaming properties and spectral characteristics of the emission
may therefore be quite different in comparison to calculations
which only employ E
(1)
||
and E
(2)
||
. In addition, while we use
E-field expressions derived in the small-angle approximation,
it is preferable to use the full solution of the Poisson equation,
particularly in the case of MSPs which have relatively small
magnetospheres and therefore much larger PC angles compared
to canonical pulsars.
In this paper, we calculate η
c
(P,
˙
P,α,ξ,φ
pc
) explicitly for
each B-field line according to the following criteria (we use
˙
P = 10
−20
, M = 1.4 M

, R = 10
6
cm, and I = 0.4 MR
2
throughout). We require that the resulting E-field should
1. be negative for all 1  η  c/(ΩR);
2. match the part of the E
(2)
||
-solution which exceeds E
(3)
||
in
absolute magnitude (i.e., where −E
(2)
||
> −E
(3)
||
) as closely
as possible; and
3. tend toward E
(3)
||
for large η.
The first criterion is required to mitigate the problem of
particle oscillations which occurs when the E-field reverses sign
beyond some altitude. Instead of this happening ∼ 40% of the

804 VENTER, HARDING, & GUILLEMOT Vol. 707
Figure 3. Contour plots of our solutions of η
c
for P = 5ms,andforα = 10
◦
, 20
◦
,...,90
◦
; ξ is the radial and φ
pc
the azimuthal coordinate in each case. The
magnetic dipole axis μ is situated at the origin, pointing outward normal to the plane of the page, in each case. The rotation axis
Ω is in the direction of φ
pc
= 0, while
the leading (trailing) edge of the pulse profile originates on B-field lines with footpoints around φ
pc
∼ 90
◦
(φ
pc
∼ 270
◦
). The η
c
-solutions get progressively smaller
for φ
pc
∼ 180
◦
, and for large α, until no solution is found which satisfies our solution matching criteria (denoted by zero values or no values at all on the plots above).
We ignore the emission from those particular field lines. We expect the η
c
-distribution to reflect the symmetry of the cos φ
pc
function which is found in E
||
; the small
irregularities stem from the fact that we used interpolation on a non-uniform (ξ,φ
pc
)-grid when preparing the contour plots.
(A color version of this figure is available in the online journal.)
time (Venter & De Jager 2005a; Venter 2008), we now only have
to ignore solutions where −E
||
< 0forη>1.1for∼ 8% of the
time. The two lower-altitude solutions E
(1)
||
and E
(2)
||
have been
matched at η = η
b
, using (Venter 2008)
η
b
≈ 1+0.0123P
−0.333
. (15)
Example fits for E
||
are shown in Figure 2 for different
parameters, as noted in the caption. The top two panels show
fits for two different η
c
, while the bottom panel is an example
where no solution for η
c
is found (according to the first criterion
above).
For illustration, Figure 3 shows contour plots of η
c
∼ 1–6
for different α, ξ, and φ
pc
, and for P = 5ms;ξ is the “radial”
and φ
pc
the azimuthal coordinate for these polar plots. From
these plots, one may infer that the “oscillatory solutions” are
encountered when φ
pc
∼ 180
◦
, and for large α (which is where
the second term of E
(2)
||
becomes negative and dominates the
first positive term inside the square brackets of Equation (4)).
The η
c
-solutions become progressively smaller for these cases,
until no solution is found which satisfies the above criteria; we
ignore emission from those particular field lines.
We tested our full solution of E
||
, which incorporates E
(1)
||
through E
(3)
||
, for conservation of energy when solving the
transport equation (Equation (1)) for relativistic electron pri-
maries. Figure 4 indicates the log
10
of acceleration rate ˙γ
gain
=
˙
E
e,gain
/m
e
c
2
, loss rate ˙γ
loss
=˙γ
CR
, curvature radius ρ
c
, and
the Lorentz factor γ as functions of distance. Although we did
not find perfect radiation reaction where the acceleration and
loss terms are equal in magnitude (similar to the findings of
Venter 2008), integration of these terms along different B-field
lines yielded energy balance (i.e., conversion of electric poten-
tial energy into gamma-radiation and particle kinetic energy)
for each integration step of the particle trajectory. An example
of this is shown in Figure 5, where the graph of the cumula-
tive energy gain (

η
η=1
dγ
gain
) coincides with that of the sum
of the cumulative energy losses and the acquired particle en-
ergy (

η
η=1
dγ
loss
+ γ (η) −γ
0
) for all η (to within ∼0.3%), with
γ
0
= γ (η = 1) the initial Lorentz factor at the stellar surface.
We used γ
0
= 100, but the calculation is quite insensitive to this
assumption, as γ quickly reaches values of ∼ 10
7
(Figure 4).
The quantities in Figure 5 are plotted in units of m
e
c
2
.
2.3. Generation of Phaseplots
In the case of the PSPC model, we normalize the particle
outflow along each B-line according to
d
˙
N(ξ,φ
pc
) =−
ρ
e
(η = 1,ξ,φ
pc
)
e
dSβ
0
c, (16)

Frequently asked questions

Q1. What are the contributions mentioned in the paper "C: " ?

The authors present results from three-dimensional emission modeling, including the special relativistic effects of aberration and time-of-flight delays and also rotational sweepback of B-field lines, in the geometric context of polar cap ( PC ), outer gap ( OG ), and twopole caustic ( TPC ) pulsar models. 

Q2. What are the future works in "C: " ?

Future work therefore includes using full acceleration and radiation models to study gamma-ray spectra, luminosities, and light curves, in order to constrain fundamental electrodynamical quantities, and possibly providing the opportunity of probing the emission geometry and B-field structure more deeply. Improved understanding of pulsar models will also feed back into more accurate population synthesis models ( e. g., Story et al. 2007 ).