CONSTRAINING PULSAR MAGNETOSPHERE GEOMETRY WITH γ-RAY LIGHT CURVES

To start, one must choose a field line structure. In the earliest light curve modeling (Morini 1983; Chiang & Romani 1994) static dipole magnetic field lines were assumed. Emission was assumed beamed along these field lines, which were taken to occupy the rotating frame, and the resulting aberrations and time delay formed the observed pulse. While this "static" (Stat) model had the right topology (closed zone plus flaring open zone), and already predicted (for outer magnetosphere emission) the generic double γ-ray pulse lagged from the low altitude radio pulsations, it clearly is not a physically consistent picture.

The next step, introduced by Romani & Yadigaroglu (1995) and used widely by other modelers, is to use the retarded dipole instead of the static case. This is commonly referred to as the Deutsch (1955) field approximation. However, at large radius from the star, this is effectively just a rotating point dipole which is most simply expressed by Kaburaki (1980) and can be re-written as

$$\begin{eqnarray} &&{\vec{B}} = -[ {\vec{m}}(t_r) + r {\dot{ \vec{m}}}(t_r) +r^2 {\ddot{ \vec{m}}}(t_r) ]/r^3 \nonumber\\ &&\quad + {\hat{r}} ({\hat{r}}\cdot [ 3 {\vec{m}}(t_r) + 3r {\dot{ \vec{m}}}(t_r) +r^2 {\ddot{ \vec{m}}}(t_r) ]/r^3)\nonumber\\ &&{\vec{m}}(t_r) = m ({\rm sin}\alpha \, {\rm cos} \omega t_r\, {\hat{x}} + {\rm sin}\alpha \,{\rm sin} \omega t_r \,{\hat{y}} + {\rm cos} \alpha \,{\hat{z}}) \end{eqnarray} \tag{ 4 }$$

with m being the dipole magnetic moment evaluated at the retarded time t_r = t − r/c. This is the field structure used in Watters et al. (2009) and so we refer to it as the "atlas" (Atl) field. It has the virtue of having the intuitively expected "sweepback" as one approaches the light cylinder (Figure 1).

Zoom In Zoom Out Reset image size

This field is stationary in the non-inertial frame rotating with respect to the lab frame, where particles are assumed to follow the field lines. In Watters et al. (2009) and earlier computations by a number of groups, the particle velocity in the lab frame was incorrectly computed as a co-rotation induced boost to a relativistic particle following the rotating field lines. Bai & Spitkovsky (2009a) correctly point out that the particle velocities in the two frames are instead connected by a simple coordinate transformation. These authors assumed that the earlier computations (incorrectly) used the field (Equation (4)) in the rotating frame. The boost error of Watters et al. (2009) is conceptually different, but mathematically equivalent to this assumption. In practice, the coordinate transformation allows one to compute the velocity β'_∥ along the magnetic field line for the particle forced to co-rotate with the star at a highly relativistic lab frame velocity $|\beta _0|\longrightarrow 1$

$$\begin{equation} {\vec{\beta }_0} = {\beta ^\prime _\parallel } {\hat{B}} + {\vec{\Omega }}\times {\vec{r}}/c \end{equation} \tag{ 5 }$$

which gives

$$\begin{equation} \beta ^\prime _\parallel =-{\hat{B}} \cdot ({\vec{\Omega }}\times {\vec{r}}/c) +\lbrace [{\hat{B}} \cdot ({\vec{\Omega }}\times {\vec{r}}/c)]^2 - ({\vec{\Omega }}\times {\vec{r}}/c)^2 +1\rbrace ^{1/2} \end{equation} \tag{ 6 }$$

so that substituting into Equation (5) gives the direction of the particle as viewed in the lab frame. Note that β'_∥ can become very small for particles traveling nearly along ${\vec{\Omega }}\times {\vec{r}}$ for r_⊥ ∼ r_LC, i.e., nearly tangent to the light cylinder in the "forward" direction. In this case, small effects of particle inertia or field line perturbations will "break open" such field lines and the particles should leave the light cylinder and contribute little to the pulse emission. Although we do not follow these physical effects, such field lines do not affect our light curves, since we taper the emissivity (see Section 2.4 below).

In any case, to have the plasma static in the rotating frame requires the presence of the co-rotation charge density that produces the lab frame field

$$\begin{equation} {\vec{E}} = ({\vec{\Omega }} \times {\vec{r}})\times {\vec{B}}/c. \end{equation} \tag{ 7 }$$

This is assumed present in the "Atlas" and similar computations, and is similarly assumed here. We refer to the revised computation as the "pseudo force-free" (PFF) case, since it contains co-rotation enforcing charges, but no currents. These PFF computations, for both TPC and OG-type models, give light curves slightly distorted from the "Atlas"-type computations and appear to match well to the sample antenna patterns and light curves of Bai & Spitkovsky (2009a). The principal effect, as noted by these authors, is that for small w < 0.02 the low altitude emission that dominates the second pulse of the TPC model has decreased phase folding, making the Jacobian non-singular. There is a pulse peak but without a true caustic there is no sharp pulse. This also weakens somewhat the first TPC peak for such high ${\dot{E}}$ pulsars. For pulsars with larger w and for the higher altitude emission of the OG model there are departures from the "Atlas"-style computation, but the effects are more subtle. In the Appendix, we present arrays of example light curves for the TPC and OG models. These may be compared with, and considered as updates to, the similar light-curve panels in Watters et al. (2009). We also update the figures describing the pulse width and flux correction factors and their variation with viewing geometry, for comparison with that paper.

The gap-accelerated radiating charges alone produce some currents in the open zone, and we expect the gap-closing pair front to produce densities comparable to the co-rotation value for at least some of the open field lines. This current will perturb the magnetic field of Equations (4).

In the spirit of attempting analytic amendments to the vacuum model, we applied the approximate perturbation field computed by Muslimov & Harding (2009) for a pair-starved current flow in the open zone. This is expressed in magnetic coordinates (r_B, θ_B, ϕ_B), where the perturbation amplitude ' determines B' = ' [2m Ω/(r² c)] with m being the magnitude of the star's dipole moment

$$\begin{eqnarray*} B^\prime _{r_B} &=& B^\prime \, s \, (1-s^2)^{-1}, \nonumber\\ B^\prime _{\theta _B} &=& B^\prime \bigg [(1-s^2)^{-1} {{ \partial s}\over {\partial \theta _B}}+ \lbrace {\rm sin}\theta _B(1-s^2)^{3/2}\rbrace ^{-1} {{ \partial s}\over {\partial \phi _B}}\bigg ], \nonumber\\ B^\prime _{\phi _B} &=& B^\prime \bigg [-(1-s^2)^{-3/2} {{ \partial s}\over {\partial \theta _B}}+ \lbrace {\rm sin}\theta _B (1-s^2)\rbrace ^{-1} {{ \partial s}\over {\partial \phi _B}}\bigg ], \ \ \, (8) \end{eqnarray*} \tag{ 8 }$$

and s = cos θ_rot = cos α cos θ_B + sin α sin θ_B cos ϕ_B. This perturbation field is defined with respect to the static vacuum dipole. We take ' positive for outward-flowing electrons (i.e., negative current) above the magnetic pole passing closest to the Earth line of sight. This pole is generally inferred to produce the low-altitude radio pulsar emission. If this pole has the opposite current (i.e., opposite magnetic polarity) we have ' < 0. Of course an electrodynamically self-consistent model will have a particular sign of ' for, e.g., the positive magnetic pole, but our PFF model includes no current and so is insensitive to the sign of B; only this perturbation current breaks the degeneracy.

To approximate the retarded solution, we matched this perturbation to our magnetosphere, mapping B(r_B,  0,  ϕ_B) to the magnetic axis of the swept-back dipole. The perturbation field becomes singular along the rotation axis where $s \longrightarrow \pm 1$. By exponentially tapering this singularity ($\propto 1-e^{(|s|-1)/\sigma _s}$, with σ_s = 0.25) we were able to integrate the field lines to determine the last closed field line surface (see below). Unfortunately, for lines near ϕ_B = 0,  π the residual effects of this singularity still dominate, making the cap boundary at these azimuths unstable (Figure 2).

Zoom In Zoom Out Reset image size

Accordingly we proceeded with a simpler toy perturbation field computed by integrating a constant line current passing through the star along the magnetic axis (the current extends to cylindrical radius r_⊥ = rsinθ = 1.2R_LC = 1.2c/Ω to minimize end effects). This perturbation field

$$\begin{equation} \qquad \qquad {\vec{B} }^\prime ({\vec{r}}) \propto \int _{B(r_b,0,\phi _B)} {\rm d}{\vec{I}} \times {\vec{r}} \end{equation} \tag{ 9 }$$

is treated as static in the rotating frame and is summed with the field of Equation (4). The effect of this simple perturbation field is dominated by the twisting of field lines along the magnetic axis, as can be seen by the shift of the "notch" foot-points of the last closed field line cap boundaries in Figure 2. This rotation is also present in the Muslimov & Harding (2009) field, but is dominated by the singularity effects except for nearly aligned rotators α ≈ 0. For ease of comparison, we scale the perturbation amplitude for fields (Equations (8) and (9)) to the underlying dipole field (Equation (4)) so that

$$$ \epsilon = \langle B^{\prime }(0.5r_{\rm LC},\theta,\phi) \rangle / \langle | B(0.5r_{\rm LC},\theta,\phi) | \rangle $$$

where we average over a sphere at r = 0.5r_LC. Again the sign of is determined by B · Ω.

None of these fields is fully realistic. In particular, the current-induced B field perturbations used above are not complete consistent field/current systems. Nevertheless, given that the vacuum field geometry has already provided encouraging successes in approximating the observed light curves, it is useful to explore how current systems can perturb the modeled beam patterns and pulse shapes. As an example, Figure 3 shows the "antenna patterns" in the hemisphere containing emission from above the null charge surface (the OG case below). The unperturbed beam pattern (middle) suffers opposite distortions for α < π/2 and α>π/2.

Zoom In Zoom Out Reset image size

For each field structure, we assume the appropriate co-rotation charge density so that the field pattern is stationary in the rotating frame. We define the open zone by identifying field lines tangent to the light cylinder (actually at slightly smaller radius for numerical expediency). These field lines are traced back to the surface where they define the open-/closed-zone boundary. In Figure 1, we show the foot-point locations, equatorial plane field lines (for α = 90°), and the last closed field line structure for an inclined α = 80° rotator. Our radiating surface is always in the open zone toward the magnetic axis from this boundary. This plot shows the static ("Stat") and point retarded dipole ("Atl" & "PFF") field structures. The active field lines for the OG are lighter (plotted in blue).

With the emitting field lines defined, we must choose the extent of the radiating surface. For the TPC model, we follow the original definition (Dyks & Rudak 2003; Dyks et al. 2004), taking the emission surface to be the full set of last-closed field lines, but stopping emission once the radial distance is r > r_LC or once the distance from the rotation axis exceeds r_⊥ > 0.75r_LC. This cut-off is required to avoid the higher altitude pulse components (as used in the OG picture), otherwise the modeled light curves generally have three or four pulse components. As in Watters et al. (2009) we make this model more physical by placing the emission surface in the open zone, separated from the last closed surface by the thickness w appropriate to the pulsar ${\dot{E}}$ (Equation (1)).

In the OG model radiation is produced on field lines crossing the null-charge surface (${\vec{\Omega }} \cdot {\vec{B}}(R_{\rm NC}) = 0)$, with emission starting near this crossing point and extending toward the light cylinder. Again, we assume a gap thickness w. Field lines that do not cross the null charge surface before r_LC do not radiate.

The OG emission should be dominated by a section of the flared cone of field lines above a given pole. For these simple vacuum-based models, some field lines just inside the last closed surface extend for a very long distance before exiting the light cylinder or crossing the null charge surface. For example, some field lines arc over the star and cross ${\vec{\Omega }} \cdot {\vec{B}}=0$ for the first time far from the star with a path length (in units of R_LC) s ≫ 1. These field lines define a disjoint, "high altitude" second gap surface, that does not in most cases appear to be active (although the HF pulse components of the Crab might represent such emission; Moffet & Hankins 1996). Physically, we expect that such lines arcing near R_LC will be opened in a real magnetosphere with currents and particle inertia. Also, gaps on field lines starting at large r_NC will likely be inactive if the soft emission needed to close the gap arises from other field lines at much smaller r_NC. Accordingly, as in Romani & Yadigaroglu (1995), we apply a path-length s cutoff in the gap surface emissivity. Here the functional form used is

$$\begin{equation} F(s) \propto e^{-[(s-2s_{\rm NC,min}-1)/\sigma _s]^2}, \qquad \sigma _s=0.1 \end{equation} \tag{ 10 }$$

with s_NC,min being the lowest null charge crossing for any active field line and all distances in units of R_LC. A detailed surface emissivity weighting awaits a physically realistic spectral emission model.

While this work was being prepared for publication, a new discussion of light curve formation in fully force-free models has been presented by Bai & Spitkovsky (2009b). These authors note the limitations with pure vacuum modeling mentioned above. They also argue that in fully force-free magnetospheres if the radiation is associated with the return current sheet, neither TPC nor OG pictures produce realistic light curves, and propose a new "Annular Gap" (AG) geometry that has emission extending beyond the light cylinder. The foot-points of the AG field lines are the same as those used here for the TPC geometry (for a given w) but in Bai & Spitkovsky (2009b) the co-rotating pulsed emission is taken to extend well beyond the light cylinder. Thus, this model covers our OG zone plus additional emission at both high and low altitudes. As we see below, the vacuum dipole field (possibly with current-induced B perturbations) does a rather good job of matching the observed pulsations—it will be interesting to see if the force-free numerical fields using the AG or OG geometries can match the data comparably well.

We identify the last closed field line surface (those field lines tangential to the speed of light cylinder) and trace these field lines to the PC surface at r_*. We follow the field line structure to <10⁻⁴r_LC, allowing good mapping even for pulsar periods of ⩾1 s. A single magnetosphere model can serve a range of pulsar periods, since the OG and TPC model features arise from distances measured in fractions of r_LC. The pulsar period only affects R_*/R_LC and hence the inner boundary of the emission zone.

Models are computed for all magnetic inclinations α. Photons are emitted along the particle path in the lab frame (moving tangent to the local field line in the co-rotating frame) and then phased, including light travel time across the magnetosphere, to add to the skymap antenna pattern, as would be seen by distant observers. We assume uniform emissivity along all active field lines (except for the taper at large pathlength s as noted above). In practice, we trace the trajectories for ∼10⁴ photons per field line, which produces adequately smooth light curves for the present computation. Models are computed for a range of gap widths from w = 0.01 to w = 0.3, to accommodate a range of pulsar ${\dot{E}}$.

Our next task is to take this family of light curves and compare with the observed pulse profiles. We take our γ-ray data from the full band E > 0.1 GeV pulse profiles published in Abdo et al. (2010). These use the first 6 months of Fermi data and give good statistics for the relatively bright pulsars considered here. In this Fermi pulsar catalog, a background level is estimated from the surrounding sky. Interestingly, some pulsars have significant steady (DC) emission. Matching this emission is an important test of the model—accordingly we fit to the total emission, pulsed and un-pulsed, after subtracting the background, as estimated in Abdo et al. (2010). In particular, we do not fit the un-pulsed flux as a separate degree of freedom—the excess above the LAT pulsar catalog background is determined by the magnetosphere model.

We take the measured pulsar ${\dot{E}}$ and compute w according to Equation (1), then use the closest value from the model grid. For this grid we take all α (in one degree intervals) and fit the model light curve for each α and ζ. To compare with the data we bin a given model curve into the phase bins (25 or 50/cycle) defining the light curve for the pulsar in Abdo et al. (2010). The total model counts are normalized to the background-subtracted data counts and we then search for the best-fit model to the combined (pulsar plus background emission) light curve.

A variety of fitting statistics work well. Simple χ² fitting generally finds the best solutions, but is controlled by the more numerous off-pulse and "bridge" emission phase bins and so occasionally provides poor discrimination against models with unacceptably weak peaks. Choosing the model with the smallest maximum model-data difference for any bin Max[Abs(Mod_i − D_i)] does an excellent job of matching the peaks, but less well on the bridge flux. We find that an exponentially tapered weighting of the largest model-data differences is robust and produces sensible results. Thus we minimize

$$$ \chi _n = {{({\rm{Mod}_i-D_i})^2} \over D_i} \, {e}^{{-i}(|{{\rm{Mod}_i}-{D_i}}|)/n} $$$

with i(|Mod_i − D_i|) being the index of the phase bin differences, sorted large to small. Small n ≪ 1 just uses the biggest difference bin (and thus focuses on the peaks), large n goes over to classical χ² weighting. The method is only weakly sensitive to n and an e-folding weighting n ≈ 3 worked well.

Our reference phase ϕ = 0 occurs when the dipole axis, spin axis, and Earth line of sight lie all in the same plane. It is generally assumed that radio emission arises at relatively low altitudes (but see Karastergiou & Johnston 2007), and that the radio pulse occupies a fraction of the open zone field lines at low altitudes. In Abdo et al. (2010) phase zero is taken to be the amplitude peak of the main radio pulse. However, the radio pulse profile often exhibits "patchy" illumination of the emission zone (Lyne & Manchester 1988). Thus, the radio peak may not mark the phase of the closest approach of the magnetic axis. Additional geometrical information comes from the linear polarization sweep of the radio pulse, with the maximum sweep rate associated with the line-of-sight passage past the magnetic axis. This maximum often occurs well away from the intensity peak. Finally, if the radio emission arises well above the star surface, relativistic effects can perturb the phases of the closest approach (ideally pulse intensity maximum) and sweep rate. Blaskiewicz et al. (1991) argue that the approximate phase shifts are

$$$ \phi _{\rm Int} \approx -r/r_{\rm LC} \qquad \qquad \phi _{\rm Pol} \approx +3r/r_{\rm LC}. $$$

Thus, polarization information can give additional constraints on the absolute phase, but altitude (plus mode changing and other sweep perturbations) can lead to substantial additional uncertainty. Accordingly, we have the option to fit light curves with the model phase free. Generally, we allow ϕ to shift by as much as ±0.1 of a rotation, and retain the best-fit model. When no radio pulse is known we allow full δϕ = ±0.5 freedom in the pulse phase. The phase for the best-fit model is recorded in the (α, ζ) map, as well as the model light curve and properties of the pulsed emission, for use in further analysis. Recalling that the radio phase ϕ = 0 in Abdo et al. (2010) is set at the pulse intensity peak, we expect that our fit, which estimates ϕ₀ for the true magnetic axis, will tend to give small δϕ>0 for emission from finite pulse heights. However, we expect that δϕ will be less than ϕ_Pol. In practice, patchy pulses and polarization sweep uncertainties make δϕ of any modest amplitude plausible.

In considering the values of the fit statistic χ_n, we must remind the reader that these are not complete models. Detailed radiation physics will certainly change the light-curve shape from the uniform weighting assumed here. Also, the lack of self-consistent magnetospheric currents must, at some level, cause light curve shape errors. This, coupled with the excellent signal-to-noise of the Fermi light curves, guarantees that these will not be statistically "good" fits. Monte Carlo experiments with the χ_n fit statistic show that for Poisson realizations of a 50 bin light curve, the typical value is 〈χ₃〉 = 7.5 while 99% and 99.9% values are 15.5 and 19. Thus fit differences of Δχ₃ ≈ 8 are statistically significant at the normal 2.5σ level. This gives a guide to interpreting the relative goodness of fits to the various models.

The model fits in (α,  ζ) space are especially useful when we have external constraints on these angles. Some care is needed, however, in applying angle constraints from other wavebands. For Vela and a number of other young pulsars with bright X-ray pulsar wind nebulae (PWNe) we have been able to place relatively model-independent constraints on viewing angle by fitting for the inclination $\tilde{\zeta }$ of the relativistic Doppler-boosted PWN torus (Ng & Romani 2008). This method does not, however, distinguish the sign of the spin axis. Therefore, the X-ray torus fits allow two possible viewing angles between the positive spin axis and the Earth line of sight $\zeta =\tilde{\zeta }$ and $\zeta ^\prime =180^\circ -\tilde{\zeta }$. In fact, in Ng & Romani (2008) the torus position angle was restricted to be 0° ⩽ Ψ < 180°, so the angles quoted in this paper may represent either ζ or ζ'.

In the context of the rotating vector model (Radhakrishnan & Cooke 1969) radio polarization data can also constrain the viewing angles. In most cases, the small range of phase illuminated by the radio pulse allows only an estimate of the magnetic impact parameter

$$\begin{equation} \quad \beta = \zeta - \alpha \approx {\rm sin^{-1}} [{\rm sin} \alpha /(d\Psi /d\phi)_{\rm max}], \end{equation} \tag{ 11 }$$

where the maximum rate of the polarization position angle sweep Ψ(ϕ) occurs at ϕ_pol = 0, near the closest approach to the magnetic axis. Here the sign of the sweep is meaningful, determining whether the line of sight is closer to or farther from the positive rotation axis than the observed magnetic pole (at inclination α). Thus, when we have an independent PWN estimate of $\tilde{\zeta }$, we have two possible α values—e.g., for a negative sweep we get β < 0 and α = ζ − β or α = ζ' − β. In this paper, we tabulate the possible values of ζ and ζ', but save space by plotting only in the positive rotation hemisphere. Thus there are two values possible for α from Equation (10), but when α or ζ>90°, we reflect back into the positive hemisphere for plotting purposes.

Occasionally, when the radio pulse is very broad or when the pulse profile presents an interpulse, the radio polarization can make meaningful estimates of both α and ζ, from fits to the full polarization sweep

$$$ {\rm tan}(\Psi + \Psi _0) = {{{\rm sin}\alpha \, {\rm sin} (\phi -\phi _0)} \over {{\rm sin}\zeta \,{\rm cos}\alpha - {\rm cos}\zeta \,{\rm sin}\alpha \,{\rm cos}(\phi -\phi _0)} }, $$$

where the polarization has the absolute position angle Ψ₀ at ϕ_Pol = ϕ₀. Keith et al. (2009) have recently presented several examples. When the radio illumination is insufficient for a full solution, we can still break the degeneracy of the fits using the X-ray $\tilde{\zeta }$ constraints. Finally, if the magnetic inclinations are, at least initially, isotropic we expect the inclination choice giving α closest to 90° to be preferred on statistical grounds.

In Figure 4, the two possible rotating vector model (RVM; Radhakrishnan & Cooke 1969) α solutions for Vela (Johnston et al. 2006) consistent with the PWN data are shown by the green circles. Note that the OG and TPC models have minima close to the solution at larger α (bold circle). In general, with constrained α, ζ values one can make a discrimination between the models. Note that for the OG model the position of the minimum is closest to the preferred angle for the relatively physical "PFF" case. The values of the fit statistic are substantially worse for TPC at all locations near these preferred angles (Table 1). Figure 5 shows the PFF light curves at the RVM angles, and at the global fit minimum, for the TPC and OG scenarios.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 1. PFF Model Angles and Efficiencies

Name	log(${\dot{E}}$)	w	da	α/ζ	α/ζ'	α/ζ/	χ/δϕ	fΩ	α/ζ/	χ/δϕ	fΩ	F> 0.1GeVa	ηOG
	(erg s-1)			$\longleftarrow$ External $\longrightarrow$		$\longleftarrow$      TPC      $\longrightarrow$			$\longleftarrow$      OG      $\longrightarrow$
Vela	36.84	0.012	0.287+0.019−0.017	72/64	122/116[58/64]	72/64/0	201/+0.02	1.06	72/64/0	120/−0.03	1.03	879.4 ± 5.4	0.013
				Best	in reg				71/64/+0.2	88/−0.02	0.98
				Best	Global				82/71/−0.2	34/−0.03	0.86
J1952 + 3252	36.57	0.016	2.0 ± 0.5	b	b	71/84/0	19/+0.05	1.10	66/78/0	14/+0.05	0.84	13.4 ± 0.9	0.015
				Best	in reg	86/66/−0.2	18/+0.05	1.10	71/71/−0.2	14/+0.05	0.70
J1709 − 4429	36.53	0.017	1.4 − 3.6	34/53	108/127[72/53]	31/49/0	147/+0.10	1.30	36/56/0	30/−0.04	0.89	124 ± 2.6	0.24
				Best	Global				47/53/+0.2	19/−0.01	0.76
Geminga	34.52	0.175	0.25+0.12−0.06	⋅⋅⋅	⋅⋅⋅	50/26/0	258/+0.41	0.93	21/89/0	241/+0.87	0.13	338.1 ± 3.5	0.05
J1057 − 5226	34.48	0.183	0.71 ± 0.2		105/111[75/69]	75/69/0	414/+0.12	0.91	75/69/0	106/+0.11	0.82	27.2 ± 1.0	0.45
				Best	Global	41/08/−0.2	37/+0.03	0.38	7/87/−0.2	5/+0.11	0.09

Notes. ^aFrom Abdo et al. (2010), flux units 10⁻¹¹erg cm² s^-1. ^bdΨ/dϕ|_max ≈ +35 deg^-1, Weisberg et al. (1999).

Download table as:  ASCIITypeset image

The best-fit models in Figure 4 are several σ off of the PWN/RVM angles. This suggests either that there are systematic errors in these angle estimates or that the true magnetosphere structure is slightly different than the PFF estimate above. The measurement of such offsets for a number of pulsars should enable us to map the required geometry differences and test physical models for their origin. As an example, here we illustrate the perturbing effect of magnetospheric currents. Figure 6 shows the χ plane for fits including a magnetic axis current (Equation (9)). The current shifts the fit minimum and distorts the light curves. For a best fit in the RVM-determined region positive currents are required for the OG model. The best models within ∼1σ of the PWN+RVM angles also prefer a small positive current. For negative the best-fit region shifts further from the externally determined angles, although the match to the pulse shape is quite good, reaching a global minimum (χ = 34) at (α,  ζ) = (82°, 71°) for = −0.2. Since the global minimum does not shift precisely to the PWN/RVM angles, we infer that (unsurprisingly) the true field geometry is only approximated by this simple model. The TPC model fits are not improved by including finite ; the best minimum anywhere near the external PWN/RVM angles has χ ≈ 200.

Zoom In Zoom Out Reset image size

To summarize, for Vela with the external constraint, PFF is clearly preferred and the best-fit models are relatively close to the known α, ζ. The rotation axis is inferred to be oriented out of the plane of the sky, we are viewing at ζ ≈ 65° and the magnetic inclination is large at α ≈ 75°. Positive-current field perturbations are preferred. The γ-ray pulse fits imply that ϕ₀ lies earlier than the radio peak (and earlier than the maximum sweep rate). This would suggest that Vela has patchy radio emission dominated by the "trailing" side of the radio zone.

The next most energetic pulsar in our test set is PSR J1952+3253. As for Vela, the three field structures give quite similar χ-plane maps. In Figure 7, we show the χ maps for the PFF model with three values of B'. This pulsar has two narrow peaks separated by Δ = 0.49 in phase, so we find that the TPC model can produce reasonable solutions, as indicated by the relatively dark regions in the lower panels. This object is presently interacting with the shell of its supernova remnant, CTB80, and so produces a PWN bow shock rather than a torus, precluding a $\tilde{\zeta }$ measurement. However from the max sweep rate of the radio polarization from Weisberg et al. (1999), we can infer a constraint on the magnetic impact parameter (two green diagonal bands). The sweep maximum appears to lie toward the tail of the broad radio pulse, consistent with our fit δϕ ≈ 0.05; again the radio pulse peak would represent patchy emission, in this case from the leading edge of the emission region at relatively low altitudes.

Zoom In Zoom Out Reset image size

Both models prefer relatively large α and ζ, and each can produce acceptable light curves within the region allowed by the radio data. If improved radio measurements or other angle constraints can be obtained, it may be possible to select between these two models for this object and also measure the altitude of the radio emission zone.

Next and only slightly less energetic is PSR B1706 − 44 = PSR J1709 − 4429 (Figure 8). This object, in contrast to PSR J1952+3252, has a narrow double pulse with Δ = 0.25. Here OG models clearly prefer the small α RVM solution. While the best OG fits have χ ⩽ 30, the TPC models in the PWN-determined ζ strip are poor with χ>90. Polarization sweep data are limited for this pulsar. The published plots in Johnston et al. (2006) give α ≈ 34°, which is best matched with negative currents ⩽ 0. These fits however, imply a surprisingly large offset for the magnetic pole of δϕ = −0.04 = 14°, which would put the magnetic axis at the leading edge of the radio pulse. Intriguingly, the best fit of all parameters is found with an OG PFF = 0.05 model at ζ = 53°, α = 47° (χ = 19). This very good light curve (Figure 9) has the magnetic axis near the radio peak with δϕ = −0.01, but would have a polarization sweep rate ∼3× larger that seen in present radio data. Again, improved radio polarization data might help shed light on the preferred inclination angle.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The next pulsar, PSR B1055 − 52 = PSR J1057 − 5226, is appreciably older. It shows a pulse that is nearly a square wave in the full energy band, although the energy-dependent light curves of Abdo et al. (2010) indicate that it is likely a tight double of separation Δ ≈ 0.2 with strong bridge emission. The pulsar is too old to show a bright PWN, so no ζ is available from the X-rays. Luckily in the radio it has both a very broad pulse and interpulse and radio polarization fitting can strongly constrain both α and ζ. We adopt the values from Weltevrede & Wright (2009) and note the phase of the maximum polarization sweep is clearly at ϕ = 0.08 (later) than plotted in Abdo et al. (2010). Indeed, we find that the model fits prefer positive δϕ, although the value is slightly larger than expected. The peaks on the model light curves are also stronger than in the data for these parameters; the match can be significantly improved for models with finite δw (not detailed here). As noted by Watters et al. (2009), no viable TPC models lie anywhere near the RVM-determined angles (Figure 10). For OG there is a shallow fit minimum centered on the RVM angles. Superposed on this are χ stripes, caused by the relatively coarse (25 bins) LAT light curve. With improved LAT data the fits should become more discriminating.

Zoom In Zoom Out Reset image size

For the TPC model there is a large region with modest χ values at small α and ζ, with the best solution listed in Table 1. However these models are not only inconsistent with the radio data, but all produce single pulses with one broad, approximately Gaussian, component and are unlikely to represent the true solution. If one allows any viewing angle, exceptionally good fits can in fact be found for OG with χ slightly below the statistical minimum value. Despite an excellent match to the observed γ-ray pulse, this is also unlikely to represent the true solution. This underlines the fact that sensible solutions must be compatible with external angle constraints and, when these constraints are applied with care, one may use the multi-wavelength data to rule out otherwise viable pulse models.

Finally, we discuss Geminga, the archetype of the old γ-selected pulsars, and the only one known prior to the launch of Fermi. This object does show an X-ray PWN, but like many old objects it is dominated by bow shock structure (Pavlov et al. 2006) and so it will be very difficult to extract torus/jet constraints on ζ. However we do know that despite very sensitive radio searches, no convincing detection of pulsed radio emission has been found, implying that |β| = |ζ − α| is large. It has also been suggested that the thermal surface pulsations are most consistent with a nearly aligned rotator (small α; Pavlov et al. 2006).

Here, with no information from the radio phase we must allow the ϕ to be freely fit. The results are interesting. For all B field structures, and in particular the PFF structure, the OG model only fits satisfactorily in a small region at small α, large ζ. This means that one is viewing a nearly aligned rotator from near the spin equator, i.e., very far from the magnetic axis. This is in good agreement with the expectations of Romani & Yadigaroglu (1995) and guarantees that the pulsar should be radio faint. In Figure 11, we show surfaces computed for w = 0.18 and three trial values. TPC shows best solutions over a range of angles and . However, these have Δχ>17 larger than the best OG solution and show only a single peak with bridge emission (Figure 12). Thus while only few OG models are suitable, the light curve shape fits better and the viewing angles for Geminga are very well constrained. It would be very interesting to measure ζ and/or α from data at other wavelengths.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Table 1, we collect these model fit results, giving the values for the externally constrained fits for the unperturbed = 0 field geometry for both models. In addition, we quote best fits close to the RVM angles and the best global fits, when appropriate. For Geminga, with no strong radio constraint we only list the global best fits for each model.

Now that we have estimated the corrections to the all-sky flux f_Ω, it is worth checking how well the pulsars agree with the simple luminosity law Equation (1). A major challenge to this check is the distance imprecision; of the five pulsars considered here only Vela and Geminga have parallax distance estimates. The other distances (derived from dispersion measure, DM) are subject to appreciable uncertainty (Figure 13).

Zoom In Zoom Out Reset image size

Geminga is in reasonable agreement with the rule only if it lies at the upper end of the present distance range. PSR J1057 − 5226 is nominally somewhat more luminous than expected, but can be accommodated if the DM distance is a slight overestimate. PSR J1709 − 4429, on the other hand, is much brighter than expected unless the true distance is ∼0.7 kpc or the flux correction factor is as small as f_Ω ∼ 0.07 (or some combination of the two factors). Such changes would be quite surprising. It would be very interesting to obtain a parallax distance estimate for PSR J1709 − 4429 to eliminate that source of uncertainty.