THE EFFECT OF AN OFFSET POLAR CAP DIPOLAR MAGNETIC FIELD ON THE MODELING OF THE VELA PULSAR'S γ-RAY LIGHT CURVES

The geometric TPC pulsar model was first introduced by Dyks & Rudak (2003). Muslimov & Harding (2003) revived the physical SG model of Arons (1983), including general relativistic (GR) corrections, and argued that the SG model may be considered a physical representation of the TPC model. This gap geometry has a large radial extent, spanning from the neutron star (NS) surface along the last closed field line up to the light cylinder. The original definition stated that the maximum radial extent reached ${R}_{\max }\simeq 0.8{R}_{\mathrm{LC}}$ (Muslimov & Harding 2004). This was later extended to ${R}_{\max }\simeq 1.2{R}_{\mathrm{LC}}$ for improved fits (e.g., Venter et al. 2009, 2012). Typical transverse gap extents of 0%–5% of the PC angle have been used (Venter et al. 2009; Watters et al. 2009).

The OG model was introduced by Cheng et al. (1986a, 1986b) and elaborated by Romani & Yadigaroglu (1995). They proposed that when the primary current passes through the neutral sheet or null-charge surface (NCS; with a radius of ${R}_{\mathrm{NCS}}$, i.e., the geometric surface across which the charge density changes sign) the negative charges above this sheet will escape beyond the light cylinder. A vacuum gap region is then formed (in which the E-field parallel to the local B-field, ${E}_{}\ne 0$). Analogously, the geometric OG model has a radial extent spanning from the NCS to the light cylinder. We follow Venter et al. (2009) and Johnson et al. (2014) who considered a one-layer model with a transverse extent along the inner edge of the gap.

The B-field is one of the basic assumptions of the geometric models (others include the gap region's location and the ${\epsilon }_{\nu }$ profile in the gap). Several B-field structures have been studied in this context, including the static dipole (Griffiths 1995), the RVD (a rotating vacuum magnetosphere which can in principle accelerate particles but do not contain any charges or currents; Deutsch 1955), the FF (filled with charges and currents, but unable to accelerate particles, since the accelerating E-field is screened everywhere; Contopoulos et al. 1999), and the offset-PC dipole (which analytically mimics deviations from the static dipole near the stellar surface; Harding & Muslimov 2011a, 2011b). A more realistic pulsar magnetosphere, i.e., a dissipative solution (Kalapotharakos et al. 2012; Li et al. 2012; Tchekhovskoy et al. 2013; Li 2014), would be one that is intermediate between the RVD and the FF fields. The dissipative B-field is characterized by the plasma conductivity ${\sigma }_{{\rm{c}}}$ (e.g., Lichnerowicz 1967), which can be chosen in order to alternate between the vacuum (${\sigma }_{{\rm{c}}}\to 0$) and FF (${\sigma }_{{\rm{c}}}\to \infty$) cases (see Li et al. 2012). We studied the effect of different magnetospheric structures (static dipole, RVD, and offset-PC dipole, further discussed below) and emission geometries (TPC and OG) on pulsar visibility and γ-ray pulse shape, particularly for the case of the Vela pulsar.

The solution for a B-field surrounding a rotating NS in vacuum was first derived by Deutsch (1955; see, e.g., Yadigaroglu 1997; Arendt & Eilek 1998; Jackson 1999; Cheng et al. 2000; Dyks & Harding 2004). The general expression for this RVD field in the rotation frame (where $\hat{{\boldsymbol{z}}}| | {\boldsymbol{\Omega }}$) is given by

$$\begin{eqnarray}{{\boldsymbol{B}}}_{\mathrm{ret}} & = & -\left[\displaystyle \frac{{\boldsymbol{\mu }}(t)}{{r}^{3}}+\displaystyle \frac{\dot{{\boldsymbol{\mu }}}(t)}{{{cr}}^{2}}+\displaystyle \frac{\ddot{{\boldsymbol{\mu }}}(t)}{{c}^{2}r}\right]\\ & & +\hat{{\boldsymbol{r}}}\hat{{\boldsymbol{r}}}\,\cdot \,\left[3\displaystyle \frac{{\boldsymbol{\mu }}(t)}{{r}^{3}}+3\displaystyle \frac{\dot{{\boldsymbol{\mu }}}(t)}{{{cr}}^{2}}+\displaystyle \frac{\ddot{{\boldsymbol{\mu }}}(t)}{{c}^{2}r}\right],\end{eqnarray} \tag{ 1 }$$

with the magnetic moment

$$\begin{eqnarray}&&{\boldsymbol{\mu }}(t)=\mu (\sin \alpha \cos {\rm{\Omega }}t\hat{{\boldsymbol{x}}}+\sin \alpha \sin {\rm{\Omega }}t\hat{{\boldsymbol{y}}}+\cos \alpha \hat{{\boldsymbol{z}}}),\end{eqnarray} \tag{ 2 }$$

where $\dot{{\boldsymbol{\mu }}}(t)$ and $\ddot{{\boldsymbol{\mu }}}(t)$ are the first and second time derivatives of the magnetic moment, and ($\hat{{\boldsymbol{x}}}$, $\hat{{\boldsymbol{y}}}$, $\hat{{\boldsymbol{z}}}$) are the Cartesian unit vectors. See Dyks & Harding (2004) for the RVD solution in spherical or Cartesian coordinates in the laboratory frame. In the limit of $r/{R}_{\mathrm{LC}}\to 0$, the retarded field simplifies to the non-aligned static dipole B-field ($\alpha \ne 0$). For the static dipole the field lines are symmetric about the ${\boldsymbol{\mu }}$-axis, whereas the RVD is distorted due to sweepback of the field lines as the NS rotates. This has implications for the definition of the PC (see Section 2.3).

The offset-PC dipole is a heuristic model of a non-dipolar magnetic structure where the PCs are offset from the ${\boldsymbol{\mu }}$-axis. The B-field lines are therefore azimuthally asymmetric compared to those of a pure static dipole field. This leads to field lines having a smaller curvature radius ${\rho }_{\mathrm{curv}}$ over half of the PC (in the direction of the PC offset) compared to those of the other half. Such small distortions in the B-field structure are due to retardation and asymmetric currents, thereby shifting the PCs by small amounts in different directions.

Harding & Muslimov (2011a, 2011b) considered two cases, i.e., symmetric and asymmetric PC offsets. The symmetric case involves an offset of both PCs in the same direction and applies to NSs with some interior current distortions that produce multipolar components near the stellar surface. The asymmetric case is associated with asymmetric PC offsets in opposite directions and applies to PC offsets due to retardation and/or currents of the global magnetosphere. These non-dipolar B-field geometries are motivated by observations of thermal X-ray emission, e.g., pulse profiles from MSPs such as PSR J0437−4715 (Bogdanov et al. 2007) and PSR J0030+0451 (Bogdanov & Grindlay 2009), with the B-fields of NSs in low-mass X-ray binaries being even more distorted (Lamb et al. 2009). Magnetic fields such as the FF solution (characterized by different PC currents than those assumed in space-charge limited flow models; Contopoulos et al. 1999; Timokhin 2006) will undergo larger sweepback of field lines near the light cylinder, and consequently display a larger offset of the PC toward the trailing side (opposite the rotation direction) than in the RVD field (which has offset PCs due to rotation alone). In what follows, we decided to study the effect of the simpler symmetric case (which does not mimic field line sweepback of FF, RVD, or dissipative magnetospheres) on predicted light curves. In the future, one can also include the more complex asymmetric case.

The general expression for a symmetric offset-PC dipole B-field in spherical coordinates $({r}^{\prime },{\theta }^{\prime },{\phi }^{\prime })$ in the magnetic frame (indicated by the primed coordinates, where ${\hat{{\boldsymbol{z}}}}^{\prime }\parallel {\boldsymbol{\mu }}$) is as follows (Harding & Muslimov 2011b)

$$\begin{eqnarray}{{\boldsymbol{B}}}_{\mathrm{OPCs}}^{\prime } & \approx & \displaystyle \frac{{\mu }^{\prime }}{{r}^{\prime 3}}\left[\cos {\theta }^{\prime }{\hat{{\boldsymbol{r}}}}^{\prime }+\displaystyle \frac{1}{2}(1+a)\sin {\theta }^{\prime }{\hat{{\boldsymbol{\theta }}}}^{\prime }\right.\\ & & -\epsilon \sin {\theta }^{\prime }\cos {\theta }^{\prime }\sin ({\phi }^{\prime }-{\phi }_{0}^{\prime }){\hat{{\boldsymbol{\phi }}}}^{\prime }],\end{eqnarray} \tag{ 3 }$$

where ${\mu }^{\prime }={B}_{0}{R}^{3}$ is the magnetic moment, R the stellar radius, B₀ the surface B-field strength at the magnetic pole, ${\phi }_{0}^{\prime }$ the magnetic azimuthal angle defining the plane in which the offset occurs, and $a=\epsilon \cos ({\phi }^{\prime }-{\phi }_{0}^{\prime })$ characterizes the offset direction in the ${x}^{\prime }$−${z}^{\prime }$ plane. If ${\phi }_{0}^{\prime }=0$ or ${\phi }_{0}^{\prime }=\pi$ the offset is in the ${x}^{\prime }$ direction (i.e., along the ${x}^{\prime }$-axis). If ${\phi }_{0}^{\prime }=\pi /2$ or ${\phi }_{0}^{\prime }=3\pi /2$ the offset is in the ${y}^{\prime }$ direction. The B-field lines are distorted in all directions. This distortion depends on parameters (related to the magnitude of the shift of the PC from the magnetic axis) and ${\phi }_{0}^{\prime }$ (we choose ${\phi }_{0}^{\prime }=0$ in what follows). If we set $\epsilon =0$ the symmetric case reduces to a symmetric static dipole.

The distance by which the PCs are shifted on the NS surface is given by

$$\begin{eqnarray}&&{\rm{\Delta }}{r}_{\mathrm{PC}}\simeq R{\theta }_{\mathrm{PC}}[1-{\theta }_{\mathrm{PC}}^{\epsilon }],\end{eqnarray} \tag{ 4 }$$

where ${\theta }_{\mathrm{PC}}={({\rm{\Omega }}R/c)}^{1/2}$ is the standard half-angle of the PC. This effective shift of the PCs is a fraction of ${\theta }_{\mathrm{PC}}$; therefore it is a larger fraction of R for pulsars with shorter periods (Harding & Muslimov 2011a). Harding & Muslimov (2011b) found that for the RVD solution, $\epsilon =0.03-0.1$, where offsets as large as 0.1 are associated with MSPs with large ${\theta }_{\mathrm{PC}}$. However, $\epsilon =0.09-0.2$ is expected for FF fields, with the larger offset values related to MSPs (Bai & Spitkovsky 2010b).

The difference between our offset-PC field and a dipole field that is offset with respect to the stellar center can be most clearly seen by performing a multipolar expansion of these respective fields. Lowrie (2011) gives the scalar potential W for an equatorially offset dipole (EOD) field as

$$\begin{eqnarray*}&&{W}^{\prime }({r}^{\prime },{\theta }^{\prime })\\ &&\quad =\left[\displaystyle \frac{{\mu }^{\prime }}{{r}^{\prime 2}}\cos {\theta }^{\prime }+\displaystyle \frac{{\mu }^{\prime }d}{{r}^{\prime 3}}\sin {\theta }^{\prime }\cos {\theta }^{\prime }+\displaystyle \frac{{\mu }^{\prime }d}{2{r}^{\prime 3}}{\sin }^{2}{\theta }^{\prime }\right],\end{eqnarray*}$$

with d being the offset parameter, and with the first few leading terms in $d/{r}^{\prime }$ listed above. From this potential, we may construct the magnetic field using $B=-{\rm{\nabla }}W$:

$$\begin{eqnarray}&&{{\boldsymbol{B}}}_{\mathrm{EOD}}^{\prime }({r}^{\prime },{\theta }^{\prime })={{\boldsymbol{B}}}_{\mathrm{dip}}^{\prime }({r}^{\prime },{\theta }^{\prime })\\ &&\qquad +\left[\displaystyle \frac{3{\mu }^{\prime }d}{{r}^{\prime 4}}\sin {\theta }^{\prime }\cos {\theta }^{\prime }+\displaystyle \frac{3{\mu }^{\prime }d}{2{r}^{\prime 4}}{\sin }^{2}{\theta }^{\prime }\right]\hat{{{\boldsymbol{r}}}^{\prime }}\\ &&\qquad -\left[\displaystyle \frac{{\mu }^{\prime }d}{{r}^{\prime 4}}{\cos }^{2}{\theta }^{\prime }+\displaystyle \frac{{\mu }^{\prime }d}{{r}^{\prime 4}}{\sin }^{2}{\theta }^{\prime }-\displaystyle \frac{{\mu }^{\prime }d}{{r}^{\prime 4}}\sin {\theta }^{\prime }\cos {\theta }^{\prime }\right]\hat{{{\boldsymbol{\theta }}}^{\prime }},\\ &&\quad ={{\boldsymbol{B}}}_{\mathrm{dip}}^{\prime }({r}^{\prime },{\theta }^{\prime })+{\text{}}O\left(\displaystyle \frac{1}{{r}^{\prime 4}}\right).\end{eqnarray} \tag{ 5 }$$

This means that an offset dipolar field may be expressed (to lowest order) as the sum of a centered dipole and two quadropolar components. Conversely, our offset-PC field may be written as

$$\begin{eqnarray}&&{{\boldsymbol{B}}}_{\mathrm{OPCs}}^{\prime }({r}^{\prime },{\theta }^{\prime },{\phi }^{\prime })\approx {{\boldsymbol{B}}}_{\mathrm{dip}}^{\prime }({r}^{\prime },{\theta }^{\prime })+{\text{}}O\left(\displaystyle \frac{\epsilon }{{r}^{\prime 3}}\right).\end{eqnarray} \tag{ 6 }$$

Therefore, we can see that the EOD consists of a centered dipole plus quadropolar and other higher-order components (Equation (5)), while our offset-PC model (Equation (6)) consists of a centered dipole plus terms of order $a/{r}^{\prime 3}$ or $\epsilon /{r}^{\prime 3}$. Since $a\sim 0.2$ and $\epsilon \sim 0.2$, the latter terms present perturbations (e.g., poloidal and toroidal effects) to the centered dipole. Harding & Muslimov (2011a, 2011b) derived these perturbed components of the distorted magnetic field while satisfying the solenoidality condition ${\rm{\nabla }}\cdot B=0$.

We performed geometric light curve modeling using the code first developed by Dyks et al. (2004a). We extended this code by implementing an offset-PC dipole B-field (for the symmetric case), as well as the SG ${E}_{\parallel }$-field corrected for GR effects (see Section 3). We solve for the PC rim as explained in Section 3.2. The shape of the PC rim depends on the B-field structure at ${R}_{\mathrm{LC}}$. Once the PC rim has been determined, the PC is divided into self-similar (interior) rings. These rings are calculated by using open-volume coordinates (${r}_{\mathrm{ovc}}$ and ${l}_{\mathrm{ovc}}$). After the footpoints of the field lines on a (${r}_{\mathrm{ovc}}$, ${l}_{\mathrm{ovc}}$) grid have been determined, particles are followed along these lines in the corotating frame and emission from them is collected in bins of pulse phase ${\phi }_{{\rm{L}}}$ and ζ, i.e., a sky map is formed by plotting the bin contents (divided by the solid angle subtended by each bin) for a given α, and it is therefore a projection of the radiation beam. To simulate light curves, one chooses a sky map corresponding to a fixed α, then fixes ζ and plots the intensity per solid angle.

The code takes into account the structure/geometry of the B-field (since the photons are emitted tangentially to the local field line), aberration of the photon emission direction (due to rotation, to first order in $r/{R}_{\mathrm{LC}}$), and time-of-flight delays (due to distinct emission radii) to obtain the caustic emission beam (Morini 1983; Dyks et al. 2004b). However, Bai & Spitkovsky (2010a) pointed out that previous studies assumed the RVD field to be valid in the instantaneously corotating frame, but actually it is valid in the laboratory frame (implying corrections that are of second order in $r/{R}_{\mathrm{LC}}$). This implies a revised aberration formula, which we have implemented in our code.

Since the offset-PC dipole field is specified in the magnetic frame (${\hat{{\boldsymbol{z}}}}^{\prime }\parallel {{\boldsymbol{\mu }}}^{\prime };$ Harding & Muslimov 2011b), it is necessary to transform this solution to the (corotating) rotational frame ($\hat{{\boldsymbol{z}}}\parallel {\boldsymbol{\Omega }}$). In order to do so, we first need to rotate the Cartesian coordinate axes specified in the rotational frame through an angle $+\alpha$ to move to the magnetic frame, and then perform transformations between the Cartesian and spherical coordinates in the magnetic frame. Only then can we transform the B-vector components from a spherical base unit vector set to a Cartesian one. We lastly perform a rotation of the Cartesian B-vector components through $-\alpha$, to move from the magnetic to the rotational frame. See the Appendix for a systematic discussion.

The object is to find the polar angle ${\theta }_{* }$ at each azimuthal angle ϕ at the footpoints of the last open B-field lines, lying within a bracket ${\theta }_{\min }\lt {\theta }_{* }\lt {\theta }_{\max }$, such that the field line is tangent to the light cylinder. The PC rim is thus defined. The magnetic structure at the light cylinder therefore determines the PC shape (Dyks & Harding 2004; Dyks et al. 2004a).

After initial implementation of the offset-PC dipole field in the geometric code, we discovered that we could solve for the PC rim in a similar manner as for the RVD B-field, but only for small values of the offset parameter ($\epsilon \lesssim 0.05-0.1$, depending on α). We improve the range of by varying the colatitude parameters ${\theta }_{\min }$ and ${\theta }_{\max }$, which delimit a bracket ("solution space") in colatitude thought to contain the footpoint of last open field line (tangent to the light cylinder ${R}_{\mathrm{LC}}$). We obtain a progressively larger range of upon decreasing ${\theta }_{\min }$ and increasing ${\theta }_{\max }$. We find a maximum $\epsilon =0.18$ valid for the full range of α. Choosing a maximal solution bracket in colatitude would in principle work, but the code would take much longer to find the PC rim compared to when a smaller bracket (that does contain the correct solution) is used. Therefore, we generalize the search for optimal ${\theta }_{\min }$ and found (by trial and error) that the following linear equation ${\theta }_{\min }=[(-31/18)\epsilon +0.6]{\theta }_{\mathrm{PC}}$, for a fixed ${\theta }_{\max }=2.0$, resulted in ${\theta }_{\min }$ that yielded maximum values for .

We illustrate the PC shape for a few cases of α and in Figure 1. We plot the PC rims (${r}_{\mathrm{ovc}}=1$) in the ${x}^{\prime }$−${y}^{\prime }$ plane (in the magnetic frame, in units of ${r}_{\mathrm{PC}}$), assuming that the ${\boldsymbol{\mu }}$-axis is located perpendicularly to the page at $({x}^{\prime },{y}^{\prime })=(0,0)$ and that ${\phi }^{\prime }$ is measured counterclockwise from the positive ${x}^{\prime }$-axis. Each PC is for a different value of α in the range 10°–90° with increments of 10°. For each α we plot the PC shape for values of 0 (green solid circle), 0.09 (red dashed circle), and 0.18 (blue dashed–dotted circle). The horizontal line at ${x}^{\prime }=0$ (black dotted line) serves as a reference line to show the magnitude and direction of offset as is increased. As α and are increased the PC shape changes considerably. For larger the PC offset is larger along the $-{x}^{\prime }$-axis (in the direction of "unfavorably curved" B-field lines). Also, as α increases for each the PC shape along the ${x}^{\prime }$-axis becomes narrower and irregular, e.g., compare the cases of $\epsilon =0$ and $\epsilon =0.18$ for $\alpha =90^\circ$. We note that the PC also becomes slightly narrower along the ${y}^{\prime }$-axis as increases.

Zoom In Zoom Out Reset image size

It is important to take the accelerating E-field into account when such expressions are available, since this will modulate the emissivity ${\epsilon }_{\nu }(s)$ (as a function of arclength s along the B-field line) in the gap as opposed to geometric models where we assume constant ${\epsilon }_{\nu }$ per unit length in the corotating frame. For the SG case we implement the full E-field in the rotational frame corrected for GR effects (e.g., Muslimov & Harding 2003, 2004). This solution consists of a low-altitude and high-altitude limit which we have to match on each B-field line. The low-altitude solution is given by A. K. Harding (2016, private communication)

$$\begin{eqnarray}{E}_{\parallel ,\mathrm{low}} & \approx & -\,3{{ \mathcal E }}_{0}{\nu }_{\mathrm{SG}}{x}^{a}\\ & & \times \left\{\displaystyle \frac{\kappa }{{\eta }^{4}}{e}_{1{\rm{A}}}\cos \alpha +\displaystyle \frac{1}{4}\displaystyle \frac{{\theta }_{\mathrm{PC}}^{1+a}}{\eta }[{e}_{2{\rm{A}}}\cos {\phi }_{\mathrm{PC}}\right.\\ & & \left.+\displaystyle \frac{1}{4}\epsilon \kappa {e}_{3{\rm{A}}}(2\cos {\phi }_{0}^{\prime }-\cos (2{\phi }_{\mathrm{PC}}-{\phi }_{0}^{\prime }))\right]\\ & & \times \sin \alpha (1-{\xi }_{* }^{2}),\end{eqnarray} \tag{ 7 }$$

with ${{ \mathcal E }}_{0}={({\rm{\Omega }}R/c)}^{2}({B}_{r}/B){B}_{0}$, B_r the radial B-field component, ${\nu }_{\mathrm{SG}}\equiv (1/4){\rm{\Delta }}{\xi }_{\mathrm{SG}}^{2}$, and ${\rm{\Delta }}{\xi }_{\mathrm{SG}}$ the colatitudinal gap width in units of dimensionless colatitude $\xi =\theta /{\theta }_{\mathrm{PC}}$. Also, $x=r/{R}_{\mathrm{LC}}$ is the normalized radial distance in units of ${R}_{\mathrm{LC}}$. Here, $\kappa \approx 0.15{I}_{45}/{R}_{6}^{3}$ is a GR compactness parameter characterizing the frame-dragging effect near the stellar surface (Muslimov & Harding 1997), ${I}_{45}=I/{10}^{45}$ g cm², I the moment of inertia, ${R}_{6}=R/{10}^{6}\,\mathrm{cm}$, $\eta =r/R$ the dimensionless radial coordinate in units of R, ${e}_{1{\rm{A}}}=1+a({\eta }^{3}-1)/3$, ${e}_{2{\rm{A}}}=(1+3a){\eta }^{(1+a)/2}\ -2a$, and ${e}_{3{\rm{A}}}=[(5-3a)/{\eta }^{(5-a)/2}]+2a$. The magnetic azimuthal angle ${\phi }_{\mathrm{PC}}$ is defined for usage with the E-field, being π out of phase with ${\phi }^{\prime }$ (one chooses the negative x-axis toward ${\boldsymbol{\Omega }}$ to coincide with ${\phi }_{\mathrm{PC}}=0$, labeling the "favorably curved" B-field lines). We define ${\phi }^{\prime }=\arctan ({y}^{\prime }/{x}^{\prime })$ the magnetic azimuthal angle used when transforming the B-field (Section 3.1). Lastly, ${\xi }_{* }$ is the dimensionless colatitude labeling the gap field lines (defined such that ${\xi }_{* }=0$ corresponds to the field line in the middle of the gap and ${\xi }_{* }=1$ at the boundaries; Muslimov & Harding 2003).

We approximate the high-altitude SG E-field by (Muslimov & Harding 2004)

$$\begin{eqnarray}{E}_{\parallel ,\mathrm{high}} & \approx & -\,\displaystyle \frac{3}{8}{\left(\displaystyle \frac{{\rm{\Omega }}R}{c}\right)}^{3}\displaystyle \frac{{B}_{0}}{f(1)}{\nu }_{\mathrm{SG}}{x}^{a}\\ & & \times \left\{\left[1+\displaystyle \frac{1}{3}\kappa \left(5-\displaystyle \frac{8}{{\eta }_{{\rm{c}}}^{3}}\right)+2\displaystyle \frac{\eta }{{\eta }_{\mathrm{LC}}}\right]\cos \alpha \right.\\ & & \left.+\displaystyle \frac{3}{2}{\theta }_{\mathrm{PC}}H(1)\sin \alpha \cos {\phi }_{\mathrm{PC}}\right\}(1-{\xi }_{* }^{2}),\end{eqnarray} \tag{ 8 }$$

with $f(\eta )\sim 1+0.75y+0.6{y}^{2}$, a GR correction factor of order 1 for the dipole component of the magnetic flux through the magnetic hemisphere of radius r in a Schwarzchild metric. The function $H(\eta )\sim 1-0.25y-0.16{y}^{2}-0.5$ $(\kappa /{\epsilon }_{{\rm{g}}}^{3})$ ${y}^{3}(1-0.25y-0.21{y}^{2})$ is also a GR correction factor of order 1, with $y={\epsilon }_{{\rm{g}}}/\eta$, ${\epsilon }_{{\rm{g}}}={r}_{{\rm{g}}}/R$, and ${r}_{{\rm{g}}}=2{GM}/{c}^{2}$ the gravitational or Schwarzchild radius of the NS (with G the gravitational constant and M the stellar mass). The factors $f(\eta )$ and $H(\eta )$ account for the static part of the curved spacetime metric and have a value of 1 in flat space (Muslimov & Harding 1997). The critical scaled radius ${\eta }_{{\rm{c}}}={r}_{{\rm{c}}}/R$ is where the high-altitude and low-altitude E-field solutions are matched, with ${r}_{{\rm{c}}}$ the critical radius and ${\eta }_{\mathrm{LC}}={R}_{\mathrm{LC}}/R$. This high-altitude solution (excluding the factor x^a) is actually valid for the SG model assuming a static (GR-corrected, non-offset) dipole field. We therefore scale the E-field by a factor x^a to generalize this expression for the offset-PC dipole field. The general E-field valid from R to ${R}_{\mathrm{LC}}$ (i.e., over the entire length of the gap) is constructed as follows (see Equation (59) of Muslimov & Harding 2004)

$$\begin{eqnarray}&&{E}_{\parallel ,\mathrm{SG}}\simeq {E}_{\parallel ,\mathrm{low}}\exp [-(\eta -1)/({\eta }_{{\rm{c}}}-1)]+{E}_{\parallel ,\mathrm{high}}.\end{eqnarray} \tag{ 9 }$$

A more detailed discussion of the electrodynamics in the SG geometry may be found in Muslimov & Harding (2003, 2004). In the next section, we solve for ${\eta }_{{\rm{c}}}(P,\dot{P},\alpha ,\epsilon ,\xi ,{\phi }_{\mathrm{PC}})$, where P is the period and $\dot{P}$ its time derivative.

At first, we matched the low-altitude and high-altitude E-field solutions by setting ${\eta }_{{\rm{c}}}=1.4$ for simplicity (Breed et al. 2014). However, we realized that ${\eta }_{{\rm{c}}}$ may strongly vary for the different parameters. Thus, we had to solve ${\eta }_{{\rm{c}}}(P,\dot{P},\alpha ,\epsilon ,\xi ,{\phi }_{\mathrm{PC}})$ on each B-field line. In what follows we consider electrons to be the radiating particles, and our discussion will therefore generally deal with the negative of the E-field. Since particle orbits approximately coincide with the B-field lines in the corotating frame, it is important to consider the behavior of the E-field as a function of arclength s rather than η.

We solve for the matching parameter in the following way. First, we calculate ${E}_{\parallel ,\mathrm{low}}$, which is independent of ${\eta }_{{\rm{c}}}$, along the B-field. If $-{E}_{\parallel ,\mathrm{low}}\lt 0$ for all η, it will never intersect with ${E}_{\parallel ,\mathrm{high}}$ and we set ${\eta }_{{\rm{c}}}=1.1$, thereby basically using ${E}_{\parallel ,\mathrm{SG}}\approx {E}_{\parallel ,\mathrm{high}}$. Second, we step through ${\eta }_{{\rm{c}}}$ (in the range $1.1\mbox{--}5.1$), calculating ${E}_{\parallel ,\mathrm{SG}}$ and ${E}_{\parallel ,\mathrm{high}}$ as well as the ratio ${S}_{i}=S({\eta }_{i})={E}_{\parallel ,\mathrm{SG}}({\eta }_{i})/{E}_{\parallel ,\mathrm{low}}({\eta }_{i})$ for $i=1,..,N$ at different radii ${\eta }_{i}$. If ${S}_{i}\gt 1$ we use $1/{S}_{i}$. We next calculate a test statistic $T{({\eta }_{{\rm{c}}})={\sum }_{i}^{N}({S}_{i}-1)}^{2}/N$ using only E-field values where $-{E}_{\parallel ,\mathrm{low}}\gt -{E}_{\parallel ,\mathrm{high}}$ (i.e., we basically fit ${E}_{\parallel ,\mathrm{SG}}$ to ${E}_{\parallel ,\mathrm{low}}$ when $-{E}_{\parallel ,\mathrm{low}}\gt -{E}_{\parallel ,\mathrm{high}}$). We then minimize T to find the optimal ${\eta }_{{\rm{c}}}$ (similar to what was done in Figure 2 of Venter et al. 2009). In Figure 2(a), the intersection radius ${\eta }_{\mathrm{cut}}\gt {\eta }_{\mathrm{LC}}$ (i.e., ${E}_{\parallel ,\mathrm{low}}$ and ${E}_{\parallel ,\mathrm{high}}$ do not intersect within the light cylinder) and therefore we impose the restriction that the solution of ${\eta }_{{\rm{c}}}$ should lie at or below 5.1. When $-{E}_{\parallel ,\mathrm{low}}$ does not decrease as rapidly (e.g., as in Figure 2(b)) we find reasonable solutions. We note that ${E}_{\parallel ,\mathrm{SG}}$ (referred to as ${E}_{\parallel ,\mathrm{old}}$ in Figure 2) produces a bump when $-{E}_{\parallel ,\mathrm{low}}$ decreases more rapidly. To circumvent this problem we test whether $-{E}_{\parallel ,\mathrm{SG}}\lt -{E}_{\parallel ,\mathrm{high}}$ and in this case we use the intersection radius ${\eta }_{\mathrm{cut}}$ of ${E}_{\parallel ,\mathrm{low}}$ and ${E}_{\parallel ,\mathrm{high}}$, rather than ${\eta }_{{\rm{c}}}$, to match our solutions (calling this new solution ${E}_{\parallel ,\mathrm{new}};$ see Figure 2(c)). We lastly observe that for ${\phi }_{\mathrm{PC}}=\pi$ (on "unfavorably curved" field lines) for larger α, the $-{E}_{\parallel ,\mathrm{low}}$ field changes sign resulting in a small ${\eta }_{{\rm{c}}}={\eta }_{\mathrm{cut}}=1.7$ value (Figure 2(d)).

Zoom In Zoom Out Reset image size

We present ${\eta }_{{\rm{c}}}$ contours in Figure 3 for an offset parameter $\epsilon =0$ and in Figure 4 for $\epsilon =0.18$. (For plotting purposes, we set ${\eta }_{{\rm{c}}}=11$ when ${\eta }_{{\rm{c}}}\gt {\eta }_{\mathrm{LC}}$.) Since the E-field solutions have an ${x}^{a}={x}^{\epsilon \cos ({\phi }^{\prime }-{\phi }_{0}^{\prime })}={x}^{-\epsilon \cos {\phi }_{\mathrm{PC}}}$ factor dependence, a larger (non-zero) offset results in different matching contours versus the case for $\epsilon =0$. In the case of $\alpha =0$, the first term $\propto \cos \alpha$ is the only contribution to the E-field, with the factor ${x}^{a}{e}_{1{\rm{A}}}$ (with an dependence) being initially larger at low η for ${\phi }_{\mathrm{PC}}=0$ than for ${\phi }_{\mathrm{PC}}=\pi$ (x^a dominates), but rapidly decreasing with η (${e}_{1{\rm{A}}}$ dominates), leading to a lower value of ${\eta }_{{\rm{c}}}$ for ${\phi }_{\mathrm{PC}}=0$ (compare the first panel of Figure 4 to the first column of Figure 5). One should therefore note that the magnitude of one instance of the E-field with low ${\eta }_{{\rm{c}}}$ may initially be higher than another instance with high ${\eta }_{{\rm{c}}}$, but the first will decrease rapidly with η and eventually become lower than the second. In Figure 3 one can see that there is no ${\phi }_{\mathrm{PC}}$-dependence for $\alpha =0$, which is not the case for Figure 4. For a slightly larger α the second terms in Equations (7) and (8) start to contribute to the radiation. This is due to the $\sin \alpha$ term with an dependence that delivers an extra contribution (Figure 4) which is zero in the case for $\epsilon =0$ (Figure 3). At $\alpha =20^\circ$ the effects of the first and second terms seem to balance each other and therefore we find the same solution of ${\eta }_{{\rm{c}}}=5.1$ everywhere except on the SG model boundary (at $\xi \in [0.95,1.0]$) where ${\eta }_{{\rm{c}}}=1.1$, just as in the case of $\epsilon =0$ and $\alpha =0^\circ$. For values of $\alpha \gt 20^\circ$ (Figure 4) the second term $\sim \cos {\phi }_{\mathrm{PC}}$ starts to dominate and thus we find solutions of ${\eta }_{{\rm{c}}}\sim 5.1$ for ${\phi }_{\mathrm{PC}}\simeq 0$ and systematically smaller solutions for ${\phi }_{\mathrm{PC}}\simeq \pi$ as α increases and the second term $\propto \,\sin \alpha$ becomes increasingly important (in both cases of ). At $\alpha =90^\circ$ we obtain the same solution as in Figure 3 where the second term dominates (for this case $-{E}_{\parallel ,\mathrm{SG}}\lt 0$ for all η, since the Goldreich–Julian charge density ${\rho }_{\mathrm{GJ}}$ becomes positive). We note that the ${\eta }_{{\rm{c}}}$ distribution reflects two symmetries (one about ${\phi }_{\mathrm{PC}}=\pi$ and one about $\xi =0.975$, i.e., ${\xi }_{* }=0$, given our gap boundaries): that of the $\cos {\phi }_{\mathrm{PC}}$ term and that of the $(1-{\xi }_{* }^{2})$ term in the ${E}_{}$ solutions. After solving for ${\eta }_{{\rm{c}}}$, we can solve the particle transport equation along each B-field line (see Section 4.1).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We apply a standard ${\chi }^{2}$ statistical fitting technique to assist us in objectively finding the pulsar geometry (α, ζ) which best describes the observed γ-ray light curve of the Vela pulsar. We use this ${\chi }^{2}$ method to determine the best-fit parameters for each of our B-field and geometric model combinations (spanning a large parameter space). The general expression is given by

$$\begin{eqnarray}&&{\chi }^{2}=\displaystyle \sum _{i=1}^{{N}_{\mathrm{bins}}}\displaystyle \frac{{({Y}_{{\rm{d}},i}-{Y}_{{\rm{m}},i})}^{2}}{{\sigma }_{{\rm{m}},i}^{2}}\approx \displaystyle \sum _{i=1}^{{N}_{\mathrm{bins}}}\displaystyle \frac{{({Y}_{{\rm{d}},i}-{Y}_{{\rm{m}},i})}^{2}}{{Y}_{{\rm{d}},i}},\end{eqnarray} \tag{ 10 }$$

where ${Y}_{{\rm{d}},i}({\phi }_{{\rm{L}},i})$ and ${Y}_{{\rm{m}},i}({\phi }_{{\rm{L}},i})$ are the number of counts of the observed and modeled light curves (relative units at phase ${\phi }_{{\rm{L}},i}$), and ${\sigma }_{{\rm{m}},i}({\phi }_{{\rm{L}},i})$ the uncertainty of the model light curves in each phase bin $i=1,\ldots ,{N}_{\mathrm{bins}}$, with ${N}_{\mathrm{bins}}$ the number of bins. Since we do not know the uncertainty of the model, we approximate the model error by the data error, assuming ${\sigma }_{{\rm{m}},i}^{2}({\phi }_{{\rm{L}},i})\approx {Y}_{{\rm{d}},i}({\phi }_{{\rm{L}},i})$ for Poisson statistics. Since we use geometric models, with an uncertainty in the absolute intensity, we assume that the shape of the light curve is correct. The data possess a background that is also uncertain. Furthermore, Fermi has a certain response function that influences the intrinsic shape of the light curve, which reflects the sum of counts from many pulsar rotations. Given all these uncertainties, we incorporate a free amplitude parameter A to allow more freedom in terms of finding the best fit of the model light curves to the data. We normalize the model light curve to range from 0 to the maximum number of observed counts ${k}_{2}$ by using the following expression:

$$\begin{eqnarray}{Y}_{{\rm{m}}}^{\prime }({\phi }_{{\rm{L}},i}) & = & \displaystyle \frac{{Y}_{{\rm{m}}}({\phi }_{{\rm{L}},i})}{({k}_{1}+{\epsilon }_{0})}A({k}_{2}-\mathrm{BG})\\ & & +\mathrm{BG}\approx \displaystyle \frac{{Y}_{{\rm{m}}}({\phi }_{{\rm{L}},i})}{{k}_{1}}{k}_{2},\end{eqnarray} \tag{ 11 }$$

with ${k}_{1}=\max ({Y}_{{\rm{m}}}({\phi }_{{\rm{L}},i}))$, ${k}_{2}=\max ({Y}_{{\rm{d}}}({\phi }_{{\rm{L}},i}))$, ${\epsilon }_{0}$ a small value added to ensure that we do not divide by zero, A a free normalization parameter, and BG the background level of ${Y}_{{\rm{d}}}({\phi }_{{\rm{L}},i})$. We treat the data as being cyclic so we need to ensure that the model light curve is cyclic as well. The model light curve has to be re-binned in order to have the same number of bins in ${\phi }_{{\rm{L}}}$ as the data (Abdo et al. 2013). We use a Gaussian Kernel Density Estimator function to rebin and smooth the model light curve (Parzen 1962). Furthermore, we also introduce the free parameter ${\rm{\Delta }}{\phi }_{{\rm{L}}}$ which represents an arbitrary phase shift of the model light curve so as to align the model and data peaks. We choose the phase shift ${\rm{\Delta }}{\phi }_{{\rm{L}}}$ as a free parameter due to the uncertainty in the definition of ${\phi }_{{\rm{L}}}=0$ (see, e.g., Johnson et al. 2014 who also used A and ${\rm{\Delta }}{\phi }_{{\rm{L}}}$). Importantly, we note that we have not changed the relative position (the radio-to-γ phase lag δ), since this is a crucial model prediction. The radio and γ-ray emission regions are tied to the same underlying B-field structure, and δ therefore reflects important physical conditions (or model assumptions) such as a difference in emission heights of the radio and γ-ray beams.

After preparation of the model light curve, we search for the best-fit solution for each of our B-field and gap combinations over a parameter space of $\alpha \in [0^\circ ,90^\circ ]$, $\zeta \in [0^\circ ,90^\circ ]$ (both with 1° resolution), $0.5\lt A\lt 1.5$ with 0.1 resolution, and $0\lt {\rm{\Delta }}{\phi }_{{\rm{L}}}\lt 1$ with 0.05 resolution. For a chosen B-field and model geometry we iterate over each set of parameters and search for a local minimum ${\chi }^{2}$ value at a particular α and ζ. Once we have iterated over the entire parameter space (α, ζ, A, ${\rm{\Delta }}{\phi }_{{\rm{L}}}$), we obtain a global minimum value for ${\chi }^{2}$ (also called the optimal ${\chi }^{2}$):

$$\begin{eqnarray}{\chi }_{\mathrm{opt}}^{2} & \approx & \displaystyle \sum _{i=1}^{{N}_{\mathrm{bins}}}\displaystyle \frac{{({Y}_{{\rm{d}},i}-{Y}_{\mathrm{opt},i})}^{2}}{{Y}_{{\rm{d}},i}}.\end{eqnarray} \tag{ 12 }$$

If faint pulsars are modeled, Poisson statistics will be sufficient to describe the observations. For the bright Vela, however, we assume Gaussian statistics which yields small errors, since the emission characteristics are more significant than those of faint pulsars. However, these small errors on the data yield large values for the reduced optimal ${\chi }^{2}$ value ${\chi }_{\mathrm{opt}}^{2}/{N}_{\mathrm{dof}}\gg 1$. We therefore need to rescale (to compensate for the uncertainty in ${\sigma }_{{\rm{m}},i}$) the ${\chi }^{2}$ values by ${\chi }_{\mathrm{opt}}^{2}$ and multiply by the number of degrees of freedom ${N}_{\mathrm{dof}}$ (the difference between ${N}_{\mathrm{bins}}$ and number of free parameters). The scaled ${\chi }^{2}$ is presented by (Pierbattista et al. 2015):

$$\begin{eqnarray}&&{\xi }^{2}={N}_{\mathrm{dof}}\displaystyle \frac{{\chi }^{2}}{{\chi }_{\mathrm{opt}}^{2}}.\end{eqnarray} \tag{ 13 }$$

From Equation (13) the ${\xi }^{2}$ for the optimal model are as follows:

$$\begin{eqnarray}&&{\xi }_{\mathrm{opt}}^{2}={N}_{\mathrm{dof}}\displaystyle \frac{{\chi }_{\mathrm{opt}}^{2}}{{\chi }_{\mathrm{opt}}^{2}}={N}_{\mathrm{dof}},\end{eqnarray} \tag{ 14 }$$

with ${\xi }_{\mathrm{opt}}^{2}/{N}_{\mathrm{dof}}={\xi }_{\mathrm{opt},\nu }^{2}=1$ the reduced ${\xi }_{\mathrm{opt}}^{2}$.

If one wishes to compare the optimal model to alternative models, e.g., in our case a B-field combined with several geometric models, confidence contours for 68% ($1\sigma$), $95.4 \%$ ($2\sigma$), and $99.73 \%$ ($3\sigma$) can be constructed by estimating the difference in the ${\xi }_{\mathrm{opt}}^{2}$ and the ${\xi }^{2}$ of the alternative models:

$$\begin{eqnarray}&&{\rm{\Delta }}{\xi }^{2}={\xi }^{2}-{\xi }_{\mathrm{opt}}^{2}={N}_{\mathrm{dof}}({\chi }^{2}/{\chi }_{\mathrm{opt}}^{2}-1).\end{eqnarray} \tag{ 15 }$$

The confidence intervals can be estimated by reading the ${\rm{\Delta }}{\xi }^{2}$ (i.e., ${\rm{\Delta }}{\xi }_{1\sigma ,{\mu }_{\mathrm{dof}}}^{2}$, ${\rm{\Delta }}{\xi }_{2\sigma ,{\mu }_{\mathrm{dof}}}^{2}$, and ${\rm{\Delta }}{\xi }_{3\sigma ,{\mu }_{\mathrm{dof}}}^{2}$) values from a standard ${\chi }^{2}$ table for the specified confidence interval at ${\mu }_{\mathrm{dof}}=2$ (corresponding to the two-dimensional ($\alpha ,\zeta$) grid; Lampton et al. 1976). Using these values for ${\rm{\Delta }}{\xi }^{2}$ and ${\xi }_{\mathrm{opt}}^{2}={N}_{\mathrm{dof}}$, we can determine ${\xi }^{2}={\xi }_{\mathrm{opt}}^{2}+{\rm{\Delta }}{\xi }^{2}={N}_{\mathrm{dof}}+{\rm{\Delta }}{\xi }^{2}$ (i.e., ${\xi }_{1\sigma }^{2}$, ${\xi }_{2\sigma }^{2}$, and ${\xi }_{3\sigma }^{2}$) from Equation (15), which is the value at which we plot each confidence contour. To enhance the contrast of the colors on the filled ${\chi }^{2}$ contours, we plot ${\mathrm{log}}_{10}{\xi }^{2}$ on an (α, ζ) grid, with a minimum value of ${\mathrm{log}}_{10}{\xi }_{\mathrm{opt}}^{2}={\mathrm{log}}_{10}({N}_{\mathrm{dof}})=1.98$ (corresponding to the best-fit solution by construction, i.e., after rescaling, with ${N}_{\mathrm{dof}}=100-4=96$ in our study). The best-fit solution is therefore positioned at ${\xi }_{\mathrm{opt}}^{2}=96$ and enclosed by the confidence contours with values of ${\xi }_{1\sigma ,{\mu }_{\mathrm{dof}}}^{2}=96+2.30$, ${\xi }_{2\sigma ,{\mu }_{\mathrm{dof}}}^{2}=96+6.17$, and ${\xi }_{3\sigma ,{\mu }_{\mathrm{dof}}}^{2}=96+11.8$ (see Equation (15); Press et al. 1992). We determine errors on α and ζ for the best-fit solution of each B-field and model combination using the $3\sigma$ interval connected contours. We choose errors of 1° for cases when the errors were smaller than one degree (given a model resolution of 1°). See Section 4.3. For the TPC and RVD model combination, we encountered poor qualitative and statistical fits using the ${\chi }^{2}$ method, thus an alternative solution had to be selected even though the ${\chi }^{2}$ value was larger.

Once we solved ${\eta }_{{\rm{c}}}$ (see Section 3.4), we could calculate the general E-field (${E}_{\parallel ,\mathrm{new}}$) in order to solve the particle transport equation (in the corotating frame) to obtain the particle energy ${\gamma }_{{\rm{e}}}(\eta ,\phi ,{\xi }_{* })$, necessary for determining the CR emissivity (Breed et al. 2014)

$$\begin{eqnarray}\dot{\gamma } & = & {\dot{\gamma }}_{\mathrm{gain}}+{\dot{\gamma }}_{\mathrm{loss}}=\displaystyle \frac{{{eE}}_{\parallel ,\mathrm{new}}}{{m}_{{\rm{e}}}c}-\displaystyle \frac{2{e}^{2}{\gamma }_{{\rm{e}}}^{4}}{3{\rho }_{\mathrm{curv}}^{2}{m}_{{\rm{e}}}c}\\ & = & \displaystyle \frac{1}{{m}_{{\rm{e}}}{c}^{2}}\left[{{ceE}}_{\parallel ,\mathrm{new}}-\displaystyle \frac{2{{ce}}^{2}{\gamma }_{{\rm{e}}}^{4}}{3{\rho }_{\mathrm{curv}}^{2}}\right],\end{eqnarray} \tag{ 16 }$$

with ${\dot{\gamma }}_{\mathrm{gain}}$ the gain (acceleration) rate, ${\dot{\gamma }}_{\mathrm{loss}}$ the loss rate, e the electron charge, ${m}_{{\rm{e}}}$ the electron mass, and ${m}_{{\rm{e}}}{c}^{2}$ the rest-mass energy; CRR (taking only CR losses into account) occurs when the energy gain balances the losses and $\dot{\gamma }=0$.

In Figure 5 we plot the ${\mathrm{log}}_{10}$ of $-{E}_{\parallel ,\mathrm{high}}$ (solid cyan line), $-{E}_{\parallel ,\mathrm{low}}$ (solid blue line), the general $-{E}_{\parallel ,\mathrm{SG}}$ field (using ${\eta }_{{\rm{c}}}$ as the matching parameter; $-{E}_{\parallel ,\mathrm{old}}$, solid green line) and a corrected E-field in cases where a bump was formed using the standard matching procedure (see Section 3.4, i.e., setting ${\eta }_{{\rm{c}}}={\eta }_{\mathrm{cut}};$ $-{E}_{\parallel ,\mathrm{new}}$, dashed red line), ${\dot{\gamma }}_{\mathrm{gain}}$ (solid yellow line), ${\dot{\gamma }}_{\mathrm{loss}}$ (solid magenta line), and ${\gamma }_{{\rm{e}}}$ (solid black line) as a function of arclength s/R along the B-field line. For each case we show $\epsilon =0$ (thick lines) and $\epsilon =0.18$ (thin lines) on the same plot. We note that the values for ${\phi }_{\mathrm{PC}}$ appearing on the figure are actually values on the stellar surface; these may change along the B-field lines. In the first column we set $\alpha =0^\circ$ so that only the first term $\propto \cos \alpha$ in Equations (7) and (8) contributes. In this case, $-{E}_{\parallel ,\mathrm{low}}\propto {x}^{a}[1+a({\eta }^{3}-1)/3]$, where $a=-\epsilon \cos {\phi }_{\mathrm{PC}}$. The sign of x^a stays the same but at small values of η, ${x}^{a}\approx 3$ for ${\phi }_{\mathrm{PC}}=0$ (top panel) and ${x}^{a}\approx 1/3$ for ${\phi }_{\mathrm{PC}}=\pi$ (bottom panel). This explains why $-{E}_{\parallel ,\mathrm{low}}$ is higher for ${\phi }_{\mathrm{PC}}=0$ than for ${\phi }_{\mathrm{PC}}=\pi$ (for $\epsilon \ne 0$) at low η. In the case of ${\phi }_{\mathrm{PC}}=\pi /2$ (middle panel), a = 0 and the values of ${E}_{\parallel ,\mathrm{low}}$ are very nearly the same for both $\epsilon =0$ and $\epsilon \ne 0$, the only difference stemming from the B-field structure that enters into Equation (7) through ${{ \mathcal E }}_{0}$. The low-altitude E-field is therefore enhanced in the direction of the "favorably curved" B-field lines (${\phi }_{\mathrm{PC}}=0$). The term $[1+a({\eta }^{3}-1)/3]$ changes sign beyond some η for ${\phi }_{\mathrm{PC}}=0$, explaining the behavior of $-{E}_{\parallel ,\mathrm{low}}$ beyond $s\approx 2R$. This effect is also noticeable in Figure 4 for $\alpha \approx 0^\circ$, where smaller values of the matching altitude ${\eta }_{{\rm{c}}}$ are found for ${\phi }_{\mathrm{PC}}=0$ or ${\phi }_{\mathrm{PC}}=2\pi$. This is contrary to the case of ${\phi }_{\mathrm{PC}}=\pi$ ("unfavorably curved" B-field lines) where the sign stays the same for all η. Correspondingly, ${\eta }_{{\rm{c}}}\approx 5$ around ${\phi }_{\mathrm{PC}}\ne 0$ in Figure 4 for $\alpha \approx 0^\circ$. For $-{E}_{\parallel ,\mathrm{high}}\propto {x}^{a}[1+2\eta /{\eta }_{\mathrm{LC}}]$, we note that for ${\phi }_{\mathrm{PC}}=(0,\pi /2,\pi )$, x^a changes from the following low to high values: $(3\to 1,1\to 1,1/3\to 1)$ (for $\epsilon \ne 0$, otherwise x^a = 1) while $[1+2\eta /{\eta }_{\mathrm{LC}}]$ changes from $1\to 3$ (the number to the left of $\to$ is valid for $\eta \ll 1$ and the number to the right of $\to$ is valid for $\eta \gg 1$). The combined effect of these two terms is such that $-{E}_{\parallel ,\mathrm{high}}\propto (3\to 3,1\to 3,1/3\to 3)$ for ${\phi }_{\mathrm{PC}}\ =(0,\pi /2,\pi )$ and $\epsilon \ne 0$, and $-{E}_{\parallel ,\mathrm{high}}\propto \to 3$ for $\epsilon =0$, for all ${\phi }_{\mathrm{PC}}$. Therefore the high-altitude E-field is again enhanced in the direction of the "favorably curved" B-field lines at low values of η (and suppressed for "unfavorably curved" B-field lines), coinciding at high η for the different ${\phi }_{\mathrm{PC}}$. It follows that ${\dot{\gamma }}_{\mathrm{gain}}$, ${\gamma }_{{\rm{e}}}$, and ${\dot{\gamma }}_{\mathrm{loss}}$ are higher for ${\phi }_{\mathrm{PC}}=0$, and lower for ${\phi }_{\mathrm{PC}}=\pi$ ($\epsilon \ne 0$). We note that ${E}_{\parallel ,\mathrm{high}}$ is not so strongly dependent on η as ${E}_{\parallel ,\mathrm{low}}$. There is no effect in changing ${\phi }_{\mathrm{PC}}$ in the case of $\alpha =0^\circ$ and $\epsilon =0$, since the first term does not depend on ${\phi }_{\mathrm{PC}}$.

In the second and third columns of Figure 5 we set $\alpha =45^\circ$ and $\alpha =85^\circ$, respectively. We note that $-{E}_{\parallel ,\mathrm{high}}$ displays the same behavior at low η as previously: for ${\phi }_{\mathrm{PC}}=0$, $-{E}_{\parallel ,\mathrm{high}}^{\epsilon \ne 0}\gt -{E}_{\parallel ,\mathrm{high}}^{\epsilon =0};$ these are nearly equal for ${\phi }_{\mathrm{PC}}=\pi /2$, and $-{E}_{\parallel ,\mathrm{high}}^{\epsilon \ne 0}\lt -{E}_{\parallel ,\mathrm{high}}^{\epsilon =0}$ for ${\phi }_{\mathrm{PC}}=\pi$. For ${\phi }_{\mathrm{PC}}=\pi /2$, the first term of $-{E}_{\parallel ,\mathrm{high}}$ dominates the second, and for ${\phi }_{\mathrm{PC}}=\pi$, the second term of $-{E}_{\parallel ,\mathrm{high}}$ is always negative, but the positive first term dominates and therefore $-{E}_{\parallel ,\mathrm{high}}$ does not change sign as η increases. A similar behavior is also seen for $-{E}_{\parallel ,\mathrm{low}}$ (boosted for non-zero and ${\phi }_{\mathrm{PC}}=0$). For $\alpha =45^\circ$ (second column of Figure 5), the second term $\propto \,\sin \alpha$ now contributes, stopping $-{E}_{\parallel ,\mathrm{low}}$ from changing sign along η for ${\phi }_{\mathrm{PC}}=0$ (versus the case of $\alpha =0$, first column, first row). The second term of $-{E}_{\parallel ,\mathrm{low}}\sim {x}^{a}\cos {\phi }_{\mathrm{PC}}$ is comparable to the first at low η, but quickly dominates as η increases for ${\phi }_{\mathrm{PC}}=0$. The second term of $-{E}_{\parallel ,\mathrm{low}}$ remains positive so that we find ${\eta }_{{\rm{c}}}=5.1$ in this case (compare Figures 3 and 4). For ${\phi }_{\mathrm{PC}}=\pi /2$ we note that $-{E}_{\parallel ,\mathrm{low}}^{\epsilon \ne 0}$ becomes negative with η. For ${\phi }_{\mathrm{PC}}=\pi$, the second term of $-{E}_{\parallel ,\mathrm{low}}$ is negative, forcing this field to change sign; this change takes slightly longer to occur when $\epsilon \ne 0$. The fact that $-{E}_{\parallel ,\mathrm{high}}$ is positive leads to a "recovery" of the total E-field, so that it becomes positive again at larger η. The effect of matching the E-field is seen in the evolution of ${\gamma }_{{\rm{e}}}(s)$ since ${\gamma }_{{\rm{e}}}$ is determined by ${E}_{\parallel }$.

In the third column of Figure 5, we indicate the ${E}_{\parallel }$-field⁴ for $\alpha =85^\circ$. For ${\phi }_{\mathrm{PC}}=0$, the behaviors of $-{E}_{}$, ${\gamma }_{{\rm{e}}}$, and $\dot{\gamma }$ are very similar to the first and second columns in Figure 5. For ${\phi }_{\mathrm{PC}}=\pi /2$ (i.e., a = 0, as in the case for $\epsilon =0$) the first and second terms are positive for all η and comparable in magnitude for $-{E}_{\parallel ,\mathrm{low}}$ so that the sign change in the case of $\epsilon \ne 0$ may be ascribed to the structure of the offset-PC dipole B-field (${\phi }_{\mathrm{PC}}$ increases along the field line since ${B}_{\phi }\ne 0$ so that $\cos {\phi }_{\mathrm{PC}}$ becomes negative). For ${\phi }_{\mathrm{PC}}=\pi$, $-{E}_{\parallel ,\mathrm{low}}^{\epsilon \ne 0}$ is not smaller than $-{E}_{\parallel ,\mathrm{low}}^{\epsilon =0}$ at low η. In this case the second term of $-{E}_{\parallel ,\mathrm{low}}$ is always negative, and the sum of the two terms are very close for both $\epsilon \ne 0$ and $\epsilon =0$ so that $-{E}_{\parallel ,\mathrm{low}}^{\epsilon \ne 0}\approx -{E}_{\parallel ,\mathrm{low}}^{\epsilon =0}$. Again, we see that $-{E}_{\parallel ,\mathrm{low}}$ becomes negative, but the positive $-{E}_{\parallel ,\mathrm{high}}$ leads to a recovery.

We lastly notice that the CRR limit is in fact reached in some cases, but only at high altitudes (the yellow and magenta lines reach the same value): e.g., beyond $\eta \approx 0.7{R}_{\mathrm{LC}}/R$ for ${\phi }_{\mathrm{PC}}=0$ and $\alpha =0^\circ$, and beyond $\eta \approx {R}_{\mathrm{LC}}$ for ${\phi }_{\mathrm{PC}}=0$ and $\alpha =45^\circ$. The notable exception occurs at large α where the first term of the E-field becomes lower and the second term plays a larger role, leading to smaller gain rates and therefore smaller Lorentz factors. We note the importance of actually solving ${\eta }_{{\rm{c}}}(P,\dot{P},\alpha ,\epsilon ,\xi ,{\phi }_{\mathrm{PC}})$ on each B-field line. Previously we set ${\eta }_{{\rm{c}}}=1.4$ for all cases and found that the particles did not attain the CRR limit (Breed et al. 2014). Only when we allowed larger values of ${\eta }_{{\rm{c}}}$ was $-{E}_{\parallel ,\mathrm{low}}$ boosted and we found particles reaching the CRR limit in many more cases. The relatively low SG E-field leads to small caustics on the phase plots constructed for photon energies $\gt 100\,\mathrm{MeV}$ (see Section 4.2). We therefore anticipate that a higher E-field should lead to CRR being reached at lower altitudes, as well as to extended caustic structures on these phase plots, resulting in qualitatively different light curve shapes (see Section 5.2).

We next perform simulations using a geometric modeling code (Section 2.3) which has the following free parameters: α, ζ, and (in case of the offset-PC dipole field). We fix the scaled co-latitude of the innermost ring of the gap (${r}_{\mathrm{ovc}}^{\min }=0.95$), PC rim (${r}_{\mathrm{ovc}}^{\max }=1.00$), and gap width $w={r}_{\mathrm{ovc}}^{\max }-{r}_{\mathrm{ovc}}^{\min }=0.05$. We choose $\alpha \in [0^\circ ,90^\circ ]$ with a 1° resolution, since the simulations show symmetry in both the northern and southern hemispheres of the pulsar (as per assumption). This implies that the emission signature for α is the same as to that of $\pi -\alpha$, except that the phase is shifted by half a rotation, i.e., ${\phi }_{{\rm{L}}}=0.5$. This model symmetry is also visible in ζ (the radiation pattern is a mirror image about $\zeta =90^\circ$, also including a phase shift of ${\rm{\Delta }}{\phi }_{{\rm{L}}}=0.5$). For the offset-PC dipole we choose for the TPC (assuming uniform ${\epsilon }_{\nu }$) and SG (assuming variable ${\epsilon }_{\nu }$) models a range of $\epsilon \in [0.00,0.18]$ with a resolution of 0.03. Our sky maps use ranges of ${\phi }_{{\rm{L}}}\in [-0.5,0.5]$ and $\zeta \in [0\buildrel{\circ}\over{.} 5,179\buildrel{\circ}\over{.} 5]$, both with a resolution of 2°. Unique emission characteristics are visible in the light curves depending on the choice of α, ζ, B-field structure, and emission geometry (Dyks et al. 2004a; Seyffert et al. 2015).

In Figure 6 we present the phase plots and light curves we obtained for the offset-PC dipole B-field and TPC model combination, for $\epsilon =0$ (equivalent to the static dipole solution). For larger values of α the caustics extend over a larger range in ζ, with the emission forming a "closed loop," which is a feature of the static dipole B-field at $\alpha =90^\circ$. The TPC model is visible at nearly all angle combinations, since some emission occurs below the NCS for this model, in contrast to the OG model. However, for $\alpha =90^\circ$ and ζ below 45° 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 6 with the phase plots and light curves for the static dipole case and the TPC model, we notice that they are identical at all angle combinations. This important test case implies that we successfully transformed our offset-PC dipole B-field, as discussed in the Appendix. In Figure 7 we chose an offset parameter $\epsilon =0.18$ assuming constant ${\epsilon }_{\nu }$. If we compare Figures 6 and 7, we notice that this larger 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, accompanied by the same change in the position of the PCs as in Figure 6. The light curve shapes are also slightly different.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Figure 8 we present phase plots and light curves for the offset-PC dipole B-field and $\epsilon =0$, obtaining a variable ${\epsilon }_{\nu }(s)$ due to using an SG E-field solution (with CR the dominating process for emitting γ-rays; see Sections 3.3). The caustic structure and resulting light curves are qualitatively different for various α compared to the constant ${\epsilon }_{\nu }$ case. The caustics appear smaller and less pronounced for larger α values (since ${E}_{}$ becomes lower as α increases), and extend over a smaller range in ζ. In Figure 9 we chose $\epsilon =0.18$, finding a variable ${\epsilon }_{\nu }$. If we compare Figure 9 with Figure 8 we note a new emission structure close to the PCs for small values of α and $\zeta \approx (0^\circ ,180^\circ )$. This reflects the boosted ${E}_{\parallel }$-field on the "favorably curved" B-field lines (with ${E}_{}\propto {x}^{a}\cos \alpha$, with $a=-\epsilon \cos {\phi }_{\mathrm{PC}}$ and ${\phi }_{\mathrm{PC}}=0;$ see Figure 5). In Figure 9 there is also more phase space filled than in Figure 8. The light curves generally display only one broad peak with less off-peak emission compared to Figure 6. As α and ζ increase, more peaks become visible, with emission still visible from both poles as seen for larger α and ζ values, e.g., $\alpha =75^\circ$ and $\zeta =75^\circ$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

If we compare Figure 6 with Figure 8 (also Figure 7 with Figure 9), we notice that when we take E_∥ into account the phase plots and light curves change considerably. For example, for $\alpha =90^\circ$ in the constant ${\epsilon }_{\nu }$ 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 ${\epsilon }_{\nu }$ case. Therefore, we see that both the B-field and E-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 $\gt 100\,\mathrm{MeV}$. Given the relatively low E-field there are only a few photons with energies exceeding 100 MeV.

To further illustrate the effect of changing , we present phase plots for $\alpha =70^\circ$ in Figure 10(a) and their corresponding light curves at $\zeta =50^\circ$ associated with the particular phase plot in Figure 10(b), using an offset-PC dipole field and TPC model, with ranging from 0.00 to 0.18 with intervals of 0.03, and assuming constant ${\epsilon }_{\nu }$. The caustic structure is slightly different for different values of . For $\epsilon =0$ the light curve has a single peak and as increases, the peak becomes slightly narrower. Also, for larger values, the caustic structure becomes slightly broader and more pronounced. Figures 10(c) and (d) represent the offset-PC dipole field and SG model with variable ${\epsilon }_{\nu }$. When we compare the phase plots of Figures 10(a) and (c), the caustics are dimmer and smaller due to the low SG E-field, and the light curves display less off-peak emission. As increases, some small features become more pronounced.

Zoom In Zoom Out Reset image size

In this section we present our best-fit solutions of the simulated light curves using the Vela data from Fermi. We plot some example contours of ${\mathrm{log}}_{10}{\xi }^{2}$ (color bar) as well as the optimal (α, ζ) combination. We determine errors on α and ζ for the optimal solution of each B-field and gap model combination using a bounding box delimited by a minimum and maximum value in both α and ζ which surrounds the $3\sigma$ contour. We illustrate this in Figure 11(a), with the white lines indicating the bounding box $[{\alpha }_{\min },{\alpha }_{\max }]$ and $[{\zeta }_{\min },{\zeta }_{\max }]$ (see enlargement in bottom left corner of panel (a)). We choose errors of 1° for cases when the errors were smaller than 1° (given our chosen resolution of 1°). In Figure 11(a) we indicate by a white star our overall best statistical fit for an OG model using an RVD field at $\alpha ={78}_{-1}^{+{1}^{^\circ }}$, $\zeta ={69}_{-1}^{+{2}^{^\circ }}$, A = 1.3, and ${\rm{\Delta }}{\phi }_{{\rm{L}}}=0$. The curved region ranging from $\zeta =70^\circ$ downward to $\zeta =10^\circ$, over the entire range of α, is caused by the caustic structure as seen in the phase plots (i.e., there is no emission visible for low values of α and ζ—the turquoise bottom-left region).

Zoom In Zoom Out Reset image size

The corresponding light curve fit of the model (solid red line) for the best-fit geometry to the Vela data (blue histogram) is also shown (Figure 11(b)). The observed light curve represents weighted counts per bin as a function of normalized phase ${\phi }_{{\rm{L}}}=[0,1]$ (Abdo et al. 2013). The model light curve yields a qualitatively good fit to the Vela data, exhibiting distinct qualitative features including the three peaks at the same phases, with roughly the same width, as seen in the Vela data. The OG model fits the data qualitatively better than the TPC model since the OG model displays less off-pulse emission, as seen in the phase plots and light curves in Section 4.2.

In Figure 11(c) we present our significance contour ${\mathrm{log}}_{10}{\xi }^{2}$ and in Figure 11(d) the corresponding best-fit light curve for a TPC model assuming an offset-PC dipole field, with $\epsilon =0.00$ and a constant ${\epsilon }_{\nu }$. We find an optimal solution at $\alpha ={73}_{-2}^{+{3}^{^\circ }}$, $\zeta ={45}_{-4}^{+{4}^{^\circ }}$, A = 1.3, and ${\rm{\Delta }}{\phi }_{{\rm{L}}}=0.55$. Disconnected confidence intervals are visible in this case, with the $3\sigma$ errors (using only the small connected confidence contour) on α and ζ yielding larger values than in Figure 11(a). The best-fit model light curve yields a less satisfactory fit to the Vela data, although the model exhibits one peak coinciding with the second peak in the data. The model also displays a low level of off-peak emission similar to the data.

In Figure 11(e) we present our significance contour ${\mathrm{log}}_{10}{\xi }^{2}$ and in Figure 11(f) the corresponding best-fit light curve for an SG model using an offset-PC dipole field, with $\epsilon =0.15$ and a variable ${\epsilon }_{\nu }$. For this combination, we find a best-fit solution at $\alpha ={76}_{-1}^{+{3}^{^\circ }}$, $\zeta ={48}_{-11}^{+{15}^{^\circ }}$, A = 0.7, and ${\rm{\Delta }}{\phi }_{{\rm{L}}}=0.55$. The model light curve yields a reasonable fit to the Vela data, but the peaks are low (constrained by the low level of off-peak emission, i.e., the ${\chi }^{2}$ prefers a small value for A), with the first peak being very broad and a small bump preceding the second peak when compared to the data.

Since the SG E-field (see Section 3.3) is low, CRR is reached in most cases but only at high η and small α. This low E-field also causes the phase plots to display small caustics which result in "missing structure." Therefore, we investigate the effect on the light curves of the offset-PC dipole B-field and SG model combination when we lower the minimum photon energy ${E}_{\min }$ from 100 to 1 MeV, above which we construct phase plots. In the CRR limit we can determine the CR cutoff of the CR photon spectrum as follows, using the formula of Venter & de Jager (2010)

$$\begin{eqnarray}&&{E}_{\mathrm{CR}}\sim 4{E}_{\parallel ,4}^{3/4}{\rho }_{\mathrm{curv},8}^{1/2}\quad \mathrm{GeV},\end{eqnarray} \tag{ 17 }$$

with ${\rho }_{\mathrm{curv},8}\sim {\rho }_{\mathrm{curv}}/{10}^{8}$ cm the curvature radius of the B-field line and ${E}_{\parallel ,4}\sim {E}_{}/{10}^{4}$ statvolt cm⁻¹ the E-field parallel to the B-field. For our given SG E-field with a magnitude of ${E}_{}\sim {10}^{2}$ statvolt cm⁻¹, the estimated cutoff is ${E}_{\mathrm{CR}}\sim 90\,\mathrm{MeV}$. This leads to pulsar emission being emitted in the hard X-ray waveband, and cannot be compared via ${\chi }^{2}$ to Fermi ($\gt 100$ MeV) data for the Vela pulsar. As an illustration, we present the phase plots and light curves in Figure 12 for $\epsilon =0.18$ and ${E}_{\min }\gt 1\,\mathrm{MeV}$. If we compare Figure 12 with Figure 9 we notice that a larger region of phase space is filled by caustics, especially at larger α, e.g., at $\alpha =90^\circ$ the visibility is enhanced. The peaks are also wider at low α. Sometimes extra emission features appear, leading to small changes in the light curve shapes.

Zoom In Zoom Out Reset image size

Additionally, we investigate what the effect is on the light curves when we increase the E-field. As a test we multiply the E-field by a factor of 100. Using Equation (17) we estimate a cutoff energy ${E}_{\mathrm{CR}}\sim 4\,\mathrm{GeV}$ which is in the energy range of Fermi ($\gt 100$ MeV). We present the phase plots and light curves for this larger E-field in Figure 13 for the offset-PC dipole and SG model solution with $\epsilon =0.18$. If we compare Figure 13 with Figure 9 we notice that more phase space is filled by caustics, especially at larger α. At $\alpha =90^\circ$ 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 E-field is boosted, leading to an increase in bridge emission as well as higher signal-to-noise ratio. At low α the background becomes feature-rich, but not at significant intensities, however.

Zoom In Zoom Out Reset image size

When we boost the low E-field, we find that the CRR limit is in fact reached almost immediately at lower η for certain parameter combinations of α and ${\phi }_{\mathrm{PC}}$, as shown in Figure 14. We also obtained a better ${\chi }^{2}$ best-fit solution for this larger E-field compared to the usual one, for $\epsilon =0.00$ at $\alpha ={75}_{-1}^{+3}$, $\zeta ={51}_{-5}^{+2}$, A = 1.1, and ${\rm{\Delta }}{\phi }_{{\rm{L}}}=0.55$. In Figure 15 we show our significance contour ${\mathrm{log}}_{10}{\xi }^{2}$ on the left and the corresponding best-fit light curve on the right. This offset-PC dipole B-field and SG model (using the increased E-field) combination therefore provides an overall optimal fit, second only to the RVD and OG model combination (see Section 6).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size