High-energy pulsar light curves in an offset polar cap -field geometry
The light curves and spectral properties of more than 200 -ray pulsars have been measured in unsurpassed detail in the eight years since the launch of the hugely successful Fermi Large Area Telescope (LAT) -ray mission. We performed geometric pulsar light curve modelling using static, retarded vacuum, and offset polar cap (PC) dipole -fields (the latter is characterized by a parameter ), in conjunction with standard two-pole caustic (TPC) and outer gap (OG) emission geometries. In addition to constant-emissivity geometric models, we also considered a slot gap (SG) -field associated with the offset-PC dipole -field and found that its inclusion leads to qualitatively different light curves. We therefore find that the assumed -field and especially the -field structure, as well as the emission geometry (magnetic inclination and observer angles), have a great impact on the pulsar’s visibility and its high-energy pulse shape. We compared our model light curves to the superior-quality -ray light curve of the Vela pulsar (for energies MeV). Our overall optimal light curve fit (with the lowest value) is for the retarded vacuum dipole field and OG model. We found that smaller values of are favoured for the offset-PC dipole field when assuming constant emissivity, and larger values are favoured for variable emissivity, but not significantly so. When we increased the relatively low SG -fields we found improved light curve fits, with the inferred pulsar geometry being closer to best fits from independent studies in this case. In particular, we found that such a larger SG -field (leading to variable emissivity) gives a second overall best fit. This and other indications point to the fact that the actual -field may be larger than predicted by the SG model.
High-energy pulsar light curves in an offset-PC -field geometry
\FullConference4th Annual Conference on High Energy Astrophysics in Southern Africa
25-27 August, 2016
Cape Town, South Africa
The field of -ray pulsars has been revolutionised by the launch of the Fermi Large Area Telescope (LAT; ). Over the past eight years, Fermi has detected over 200 -ray pulsars and has furthermore measured their light curves and spectral characteristics in unprecedented detail. Fermi’s Second Pulsar Catalog (2PC; ) describes the properties of some 117 of these pulsars in the energy range 100 MeV100 GeV. In this paper, we will focus on the GeV band light curves of the Vela pulsar , the brightest persistent source in the -ray sky.
Physical emission models such as the slot gap (SG; ) and outer gap (OG; [10, 38]) fall short of fully explaining (global) magnetospheric characteristics, e.g., the particle acceleration and pair production, current closure, and radiation of a complex multi-wavelength spectrum. More recent developments include global magnetospheric models such as the force-free (FF) inside and dissipative outside (FIDO) model [24, 25], the wind models of, e.g., , and particle-in-cell simulations (PIC; [8, 9]). Although much progress has been made using these physical (or emission) models, geometric light curve modeling [16, 41, 42, 22, 37] still presents a crucial avenue for probing the pulsar magnetosphere in the context of traditional pulsar models. The most commonly used emission geometries include the two-pole caustic (TPC; the SG model may be its physical representation; ) and OG models and may be used to constrain the pulsar geometry (i.e., magnetic inclination angle and the observer viewing angle with respect to the spin axis ), as well as the -ray emission region’s location and extent. This may provide vital insight into the boundary conditions and help constrain the accelerator geometry of next-generation full radiation models.
The assumed -field structure is essential for predicting the light curves seen by the observer using geometric models, since photons are expected to be emitted tangentially to the local -field lines in the corotating pulsar frame . Even a small difference in the magnetospheric structure will therefore have an impact on the light curve predictions. Additionally, we have also incorporated an SG -field associated with the offset-PC dipole -field (making this latter case an emission model), which allows us to calculate the emissivity in the acceleration region in the corotating frame from first principles.
In this paper, we investigate the impact of different magnetospheric structures (i.e., the static dipole , retarded vacuum dipole (RVD; ), and an offset-PC dipole -field solution [19, 20]), as well as the SG -field on the pulsar visibility and -ray pulse shape. In combination with the different -field solutions mentioned above, we assume standard TPC and OG emission geometries. In Section 2 we briefly describe the offset-PC dipole -field and its corresponding SG -field implemented in our code [16, 4]. We also investigate the effect of increasing the -field by a factor of a 100. In Section 3, we present our phase plots and model light curves for the Vela pulsar, and we compare our results to previous multi-wavelength studies. Our conclusions follow in Section 4.
2 The Offset-PC Magnetosphere
2.1 -field structure
Several -field structures have been studied in pulsar models, including the static dipole, the RVD (a rotating vacuum magnetosphere which can in principle accelerate particles but do not contain any charges or currents), the FF (filled with charges and currents, but unable to accelerate particles since the accelerating -field is screened everywhere; ), and the offset-PC dipole. The offset-PC dipole solution analytically mimics deviations from the static dipole near the stellar surface and is azimuthally asymmetric, with field lines having a smaller curvature radius over half of the PC (in the direction of the PC offset) compared to those of the other half [19, 20]. Such small distortions in the -field structure can be due to retardation and asymmetric currents, thereby shifting the PCs by small amounts in different directions. A more realistic pulsar magnetosphere, i.e., a dissipative solution [29, 26, 28, 39, 27], would be one that is intermediate between the RVD and the FF fields.
The symmetric case involves an offset of both PCs, with respect to the magnetic () axis in the same direction and applies to neutron stars with some interior current distortions that produce multipolar components near the stellar surface [19, 20]. We study the effect of this simpler symmetric case on predicted light curves. The general expression for a symmetric offset-PC dipole -field in spherical coordinates in the magnetic frame (indicated by the primed coordinates, where ) is 
where the symbols have the same meaning as before [19, 20]. We choose the offset direction to be in the plane. The -field lines are distorted in all directions, with the distortion depending on parameters (related to the magnitude of the shift of the PC from the magnetic axis) and (we choose in what follows, with the offset being in the direction). If we set the symmetric case reduces to a symmetric static dipole.
The difference between our offset-PC field and a dipole field that is offset with respect to the stellar centre can be most clearly seen by performing a multipolar expansion of these respective fields. An offset dipolar field may be expressed (to lowest order) as the sum of a centred dipole and quadropolar terms ). Conversely, our offset-PC field may be written as
Therefore, we can see that our offset-PC model (Eq. ) consists of a centred dipole plus terms of order or . Since and , the latter terms present perturbations (e.g., poloidal and toroidal effects) to the centred dipole. These perturbed components of the distorted magnetic field were derived under the solenoidality condition [19, 20].
2.2 Incorporating a corresponding SG -field
It is important to take the accelerating -field (-field parallel to the local -field,) into account when such expressions are available, since this will modulate the emissivity in the gap as opposed to geometric models where we assume constant per unit length in the corotating frame. For the SG case we implement the full -field in the rotational frame corrected for general relativistic (GR) effects (e.g., [32, 33]).
The low-altitude solution is given by (A.K. Harding 2015, private communication)
We approximate the high-altitude SG -field by 
The critical scaled radius is where the high-altitude and low-altitude -field solutions are matched, with the critical radius, the stellar radius, , and the light cylinder radius (where the corotation speed equals the speed of light).
2.3 Increasing the relatively low -field
In the curvature radiation reaction (CRR, where the energy gain rate equals the CR loss rate) limit, we can determine the CR cutoff of the CR photon spectrum as follows 
with cm the curvature radius of the -field line and statvolt cm. Since the SG -field (see Section 2.2) is low (implying a CR cutoff around a few MeV), the phase plots for emission MeV display small caustics (Section 3.1) which result in “missing structure”. Therefore, we investigate the effect on the light curves of the offset-PC dipole -field and SG model combination when we increase the -field. As a test we multiply Eq. (5) by a factor 100. Using the above expression the estimated cutoff energy for our increased SG -field is now GeV, which is in the energy range of Fermi ( MeV).
3.1 Phase plots and light curves
As an example we show phase plots and their corresponding light curves for the offset-PC dipole, for both the TPC (assuming uniform ) and SG (assuming variable ) models. Figure 1 is for the TPC model for . For larger values of the caustics extend over a larger range in , with the emission forming a “closed loop,” which is also a feature of the static dipole -field at . The TPC model is visible at nearly all angle combinations, since some emission occurs below the null charge surface (the geometric surface across which the charge density changes sign; ) for this model, in contrast to the OG model. However, for and below no light curves are visible, i.e., no emission is observed due to the “closed loop” structure of the caustics. The TPC light curves exhibit relatively more off-pulse emission than the OG ones. In the TPC model, emission is visible from both magnetic poles, forming double peaks in some cases, whereas in the OG model emission is visible from a single pole. One does obtain double peaks in the OG case, however, when the line of sight crosses the caustic at two different phases.
If we compare Figure 1 with the static dipole case (for ; not shown), we notice that a larger PC offset results in qualitatively different phase plots and light curves, e.g., modulation at small . Also, the caustics occupy a slightly larger region of phase space and seem more pronounced for larger and values. The light curve shapes are also slightly different.
Figure 2 is for the offset-PC dipole -field and , for a variable due to using an SG -field solution (with CR the dominating process for emitting -rays; see Sections 2.2). The caustic structure and resulting light curves are qualitatively different for various compared to the constant case. The caustics appear smaller and less pronounced for larger values (since becomes lower as increases), and extend over a smaller range in . If we compare Figure 2 with the case for (for variable ; not shown) we note a new emission structure close to the PCs for small values of and . This reflects the boosted -field on the “favourably curved” -field lines. In Figure 2 a smaller region in phase space filled. The light curves generally display only one broad peak with less off-peak emission compared to Figure 1. As and increase, more peaks become visible, with emission still visible from both poles as seen for larger and values, e.g., and .
If we compare Figures 1 with 2, we notice that when we take into account, the phase plots and light curves change considerably. For example, for in the constant case, a “closed loop” emission pattern is visible in the phase plot, which is different compared to the small “wing-like” emission pattern in the variable case. Therefore, we see that both the -field and -field have an impact on the predicted light curves. This small “wing-like” caustic pattern is due to the fact that we only included photons in the phase plot with energies MeV. Given the relatively low -field, there are only a few photons with energies exceeding MeV.
In Figure 3 we present the phase plots and light curves for the SG -field (increased by a factor 100) for the offset-PC dipole and SG model solution, with . If we compare Figure 3 with Figure 2 we notice that more phase space is filled by caustics, especially at larger . At the visibility is again enhanced. The caustic structure becomes wider and more pronounced, with extra emission features arising as seen at larger and values. This leads to small changes in the light curve shapes. At smaller values, the emission around the PC forms a circular pattern that becomes smaller as increases. These rings around the PCs become visible since the low -field is boosted, leading to an increase in bridge (region between the first and second peak of a light curve) emission as well as higher signal-to-noise ratio. At low the background becomes feature-rich, but not at significant intensities, however.
3.2 Comparison of best-fit parameters for different models
We next follow the same approach as a previous study  to compare the various optimal solutions of the different models. We determine the difference between the scaled of the optimal model, , and the other models () using
with degrees of freedom . We considered two approaches: we found the best fit (i) per -field and model combination (), and (ii) overall (for all -field and model combinations, )
In Figure 4 we label the different -field structures assumed in the various models as well as the overall comparison along the -axis, and plot and on the -axis. We represent the TPC geometry with a circle, the OG with a square, and for the offset-PC dipole field we represent the various values for constant by different coloured stars, for variable by different coloured left pointing triangles, and for the case of by different coloured upright triangles, as indicated in the legend. The dashed horizontal lines indicate the confidence levels we obtained for degrees of freedom. These confidence levels are used as indicators of when to reject or accept an alternative fit compared to the optimum fit.
For the static dipole field the TPC model gives the optimum fit and the OG model lies within of this fit, implying that the OG geometry may provide an acceptable alternative fit to the data in this case. For the RVD field the TPC model is significantly rejected beyond the level (not shown on plot), and the OG model is preferred. We show three cases for the offset-PC dipole field, including the TPC model assuming constant , the SG model assuming variable , and the latter is for an -field multiplied by a factor of . The optimal fits for the offset-PC dipole field and TPC model reveal that a smaller offset is generally preferred for constant , while a larger offset is preferred for variable (but not significantly), with all alternative fits falling within of these. However, when we increase -field, a smaller offset is preferred for the SG and variable case. When we compare all model and -field combinations with the overall best fit (i.e., rescaling the values of all combinations using the optimal fit involving the RVD -field and OG model), we notice that the static dipole and TPC model falls within , whereas the static OG model lies within . We also note that the usual offset-PC dipole -field and TPC model combination (for all values) is above (with some fits ), but the offset-PC dipole -field and SG model combination (for all values) is significantly rejected (). However, the case of the offset-PC dipole field and a higher SG -field for all values leads to a recovery, since all the fits fall within or and delivers an overall optimal fit for , second only to the RVD and OG model fit.
Several multi-wavelength studies have been performed for Vela, using the radio, X-ray, and -ray data, in order to find constraints on and . We graphically summarise the best-fit and , with errors, from this and other works in Figure 5. We notice that the best fits generally prefer a large or or both. It is encouraging that many of the best-fit solutions lie near the inferred from the pulsar wind nebula (PWN) torus fitting , notably for the RVD -field. A significant fraction of fits furthermore lie near the diagonal, i.e., they prefer a small impact angle, most probably due to radio visibility constraints . For an isotropic distribution of pulsar viewing angles, one expects values to be distributed as between , i.e., large values are much more likely than small values, which seems to agree with the large best-fit values we obtain. There seems to be a reasonable correspondence between our results obtained for geometric models and those of other authors, but less so for the offset-PC dipole -field, and in particular for the SG -field case. The lone fit near may be explained by the fact that a very similar fit, but one with slightly worse , is found at . If we discard the non-optimal TPC / SG fits, we see that the optimal fits will cluster near the other fits at large and . Although our best fits for the offset-PC dipole -field are clustered, it seems that increasing leads to a marginal decrease in for the TPC model (light green) and opposite for SG (dark green), but not significantly. For our increased SG -field case (brown) we note that the fits now cluster inside the gray area above the fits for the static dipole and TPC, and offset-PC dipole for both the TPC and SG geometries.
We investigated the impact of different magnetospheric structures (i.e., static dipole, RVD, and a symmetric offset-PC dipole fields) on predicted -ray pulsar light curve characteristics. For the offset-PC dipole field we only considered the TPC (assuming uniform ) and SG (modulating the using the -field which is corrected for GR effects up to high altitudes) models. We concluded that the magnetospheric structure and emission geometry have an important effect on the predicted -ray pulsar light curves. However, the presence of an -field may have an even greater effect than small changes in the -field and emission geometries.
We fit our model light curves to the observed Fermi-measured Vela light curve for each -field and geometric model combination. We found that the RVD field and OG model combination fit the observed light curve the best for and an unscaled . As seen in Figure 4, for the RVD field an OG model is significantly preferred over the TPC model, given the characteristically low off-peak emission. For the other field and model combinations there was no significantly preferred model (per -field), since all the alternative models may provide an acceptable alternative fit to the data, within . The offset-PC dipole field for constant favoured smaller values of , and for variable larger values, but not significantly so (). When comparing all cases (i.e., all -fields), we noted that the offset-PC dipole field for variable was significantly rejected ().
Since we wanted to compare our model light curves to Fermi data we increased the usual low SG -field by a factor of 100
We found reasonable correspondence between our results obtained for geometric models and those of other independent studies. We noted that the optimal fits generally clustered near the other fits at large and . For our increased SG -field and offset-PC dipole combination, we noted that these fits now clustered at larger and .
There have been several indications that the SG -field may be larger than initially thought, as confirmed by this study. (i) Population synthesis studies found that the SG -ray luminosity may be too low, pointing to an increased -field and / or particle current through the gap, e.g., . (ii) If the -field is too low, one is not able to reproduce the observed spectral cutoffs of a few GeV (Section 2.3; ). (iii) A larger -field (increased by a factor of 100) led to statistically improved fits with respect to the light curves. (iv) The inferred best-fit and parameters for this -field clustered near the best fits of independent studies. (v) A larger SG -field also increased the particle energy gain rates leading to CRR being reached close to the stellar surface.
Independent multi-wavelength studies have considered many other pulsars, in addition to the Vela pulsar. For example, Ng & Romani [34, 35] used torus and jet fitting to constrain of a few X-ray pulsars, and obtained a consistent value of . Johnson et al.  and Pierbattista et al.  fitted the radio and -ray light curves of millisecond and younger pulsar populations respectively using standard geometric models. DeCesar et al.  constrained the and angles of a handful of pulsars using standard emission geometries coupled with the FF -field. Overall, there seems to be reasonable consistency between the best-fit geometries derived using the various models.
A number of studies have lastly considered signatures in the polarisation domain for different -field geometries, radiation mechanisms, and emission sites, e.g., [16, 9, 21]. This avenue may well prove very important in future to aid in differentiating between the various pulsar models, in addition to spectral and light curve measurements.
Acknowledgements.We thank Marco Pierbattista, Tyrel Johnson, Lucas Guillemot, and Bertie Seyffert for fruitful discussions. This work is based on the research supported wholly / in part by the National Research Foundation of South Africa (NRF; Grant Numbers 87613, 90822, 92860, 93278, and 99072). The Grantholder acknowledges that opinions, findings and conclusions or recommendations expressed in any publication generated by the NRF supported research is that of the author(s), and that the NRF accepts no liability whatsoever in this regard. A.K.H. acknowledges the support from the NASA Astrophysics Theory Program. C.V. and A.K.H. acknowledge support from the Fermi Guest Investigator Program.
- We therefore first scale the values using the optimal value obtained for a particular -field, and second we scale these using the overall optimal value irrespective of -field.
- This number is not unreasonable, especially in light of the observed spectral high-energy spectral cutoffs. Since pulsars have high local -field strengths and we expect that , such high -fields are realistic. A larger gap width is also likely, and this will further increase .
- A. A. Abdo, M. Ackermann, W. B. Atwood et al., Fermi Large Area Telescope Observations of the Vela Pulsar, ApJ 696 1084 (2009).
- A. A. Abdo, M. Ajello, A. Allafort et al., The Second Fermi Large Area Telescope Catalog of Gamma-Ray Pulsars, ApJS 208 17 (2013).
- W. B. Atwood, A. A. Abdo, M. Ackermann et al., The Large Area Telescope on the Fermi Gamma-Ray Space Telescope Mission, ApJ 697 1071 (2009).
- M. Barnard, C. Venter, and A. K. Harding, The Effect of an Offset Polar Cap Dipolar Magnetic Field on the Modeling of the Vela Pulsar’s Gamma-Ray Light Curves, ApJ 832 107 (2016).
- M. Breed, C. Venter, A. K. Harding, and T. J. Johnson, Implementation of an Offset-dipole Magnetic Field in a Pulsar Modelling Code, in proceedings of SAIP2013: the 58 Ann. Conf. of the SA Institute of Physics ed. R. Botha and T. Jili, 350 (2014).
- M. Breed, C. Venter, A. K. Harding, and T. J. Johnson, The Effect of Different Magnetospheric Structures on Predictions of Gamma-Ray Pulsar Light Curves, in proceeding of SAIP2012: the 57 Ann. Conf. of the SA Institute of Physics ed. J. Janse van Rensburg, 316 (2015).
- M. Breed, C. Venter, A. K. Harding, and T. J. Johnson, The Effect of an Offset-dipole Magnetic Field on the Vela Pulsar’s Gamma-Ray Light Curves, in proceedings of SAIP2014: the 59 Ann. Conf. of the SA Institute of Physics ed. C. Engelbrecht and S. Karataglidis, 311 (2015).
- B. Cerutti, A. A. Philippov, and A. Spitkovsky, Modelling High-Energy Pulsar Light Curves from First Principles, MNRAS 457 2401 (2016).
- B. Cerutti, J. Mortier, and A. A. Philippov, Polarized Synchrotron Emission from the Equatorial Current Sheet in Gamma-Ray Pulsars, MNRAS 463 L89 (2016).
- K. S. Cheng, C. Ho, and M. Ruderman, Energetic Radiation from Rapidly Spinning Pulsars. I Outer Magnetosphere Gaps. II VELA and Crab, ApJ 300 500 (1986).
- I. Contopoulos, D. Kazanas, and C. Fendt, The Axisymmetric Pulsar Magnetosphere, ApJ 511 351 (1999).
- J. K. Daugherty, and A. K. Harding, Electromagnetic Cascades in Pulsars, ApJ 252 337 (1982).
- M. E. DeCesar, Using Fermi Large Area Telescope Observations to Constrain the Emission and Field Geometries of Young Gamma-Ray Pulsars and to Guide Millisecond Pulsar Searches, PhD thesis, Univ. of Maryland, College Park (2013).
- A. J. Deutsch, The Electromagnetic Field of an Idealized Star in Rigid Rotation in Vacuo, AnAp 18 1 (1955).
- J. Dyks, and B. Rudak, Two-Pole Caustic Model for High-Energy Light Curves of Pulsars, ApJ 598 1201 (2003).
- J. Dyks, A. K. Harding, and B. Rudak, Relativistic Effects and Polarization in Three High-Energy Pulsar Models, ApJ 606 1125 (2004).
- P. Goldreich, and W. H. Julian, Pulsar Electrodynamics, ApJ 157 869 (1969).
- D. J. Griffiths, Introduction to Electrodynamics, ed.; San Francisco: Pearson Benjamin Cummings (1995).
- A. K. Harding, and A. G. Muslimov, Pulsar Pair Cascades in a Distorted Magnetic Dipole Field, ApJL 726 L10 (2011).
- A. K. Harding, and A. G. Muslimov, Pulsar Pair Cascades in Magnetic Fields with Offset Polar Caps, ApJ 743 181 (2011).
- A. K. Harding, and C. Kalapotharakos, in preparation, (2017).
- T. J. Johnson, C. Venter, A. K. Harding et al., Constraints on the Emission Geometries and Spin Evolution of Gamma-Ray Millisecond Pulsars, ApJS 213 6 (2014).
- S. Johnston, G. Hobbs, S. Vigeland et al., Evidence For Alignment of the Rotation and Velocity Vectors in Pulsars, MNRAS 364 1397 (2005).
- C. Kalapotharakos, and I. Contopoulos, Three-dimensional Numerical Simulations of the Pulsar Magnetosphere: Preliminary Results, A&A 496 495 (2009).
- C. Kalapotharakos, A. K. Harding, and D. Kazanas, Gamma-Ray Emission in Dissipative Pulsar Magnetospheres: From Theory to Fermi Observations, ApJ 793 97 (2014).
- C. Kalapotharakos, D. Kazanas, A. K. Harding, and I. Contopoulos, Toward a Realistic Pulsar Magnetosphere, ApJ 749 2 (2012).
- J. G. Li, Electromagnetic and Radiative Properties of Neutron Star Magnetospheres, PhD thesis, Princeton Univ., New Jersey (2014).
- J. Li, A. Spitkovsky, and A. Tchekhovskoy, Resistive Solutions for Pulsar Magnetospheres, ApJ 746 60 (2012).
- A. Lichnerowicz, Relativistic Hydrodynamics and Magnetohydrodynamics, ed.; New York: Benjamin, Inc. (1967).
- W. Lowrie, A Student’s Guide to Geophysical Equations, ed.; Cambridge University Press (2011).
- A. G. Muslimov, and A. K. Harding, Toward the Quasi-Steady State Electrodynamics of a Neutron Star, ApJ 485 735 (1997).
- A. G. Muslimov, and A. K. Harding, Extended Acceleration in Slot Gaps and Pulsar High-Energy Emission, ApJ 588 430 (2003).
- A. G. Muslimov, and A. K. Harding, High-Altitude Particle Acceleration and Radiation in Pulsar Slot Gaps, ApJ 606 1143 (2004).
- C.-Y. Ng, and R. W. Romani, Fitting Pulsar Wind Tori., ApJ 601 479 (2004).
- C.-Y. Ng, and R. W. Romani, Fitting Pulsar Wind Tori. II. Error Analysis and Applications, ApJ 673 411 (2008).
- J. Pétri, and G. Dubus, Implication of the Striped Pulsar Wind Model for Gamma-Ray Binaries, MNRAS 417 532 (2011).
- M. Pierbattista, A. K. Harding, I. A. Grenier et al., Light-curve Modelling Constraints on the Obliquities and Aspect Angles of the Young Fermi Pulsars, A&A 575 A3 (2015).
- R. W. Romani, and I.-A. Yadigaroglu, Gamma-Ray Pulsars: Emission Zones and Viewing Geometries, ApJ 438 314 (1995).
- A. Tchekhovskoy, A. Spitkovsky, and J. G. Li, Time-dependent 3D Magnetohydrodynamic Pulsar Magnetospheres: Oblique Rotators, MNRAS 435 L1 (2013).
- C. Venter, and O. C. de Jager, Accelerating High-energy Pulsar Radiation Codes, ApJ 725 1903 (2010).
- C. Venter, A. K. Harding, and L. Guillemot, Probing Millisecond Pulsar Emission Geometry Using Light Curves from the Fermi/Large Area Telescope, ApJ 707 800 (2009).
- K. P. Watters, R. W. Romani, P. Weltevrede, and S. Johnston, An Atlas for Interpreting Gamma-Ray Pulsar Light Curves, ApJ 695 1289 (2009).