The effect of an offset polar cap dipolar magnetic field on the modeling   of the Vela pulsar's $\gamma$-ray light curves

Monica Barnard; Christo Venter; Alice K. Harding

arXiv:1701.02924·astro-ph.HE·January 12, 2017

The effect of an offset polar cap dipolar magnetic field on the modeling of the Vela pulsar's $\gamma$-ray light curves

Monica Barnard, Christo Venter, Alice K. Harding

PDF

TL;DR

This study models Vela pulsar gamma-ray light curves using various magnetic field configurations, revealing that an offset polar cap dipole field with outer gap emission best fits observations and influences particle acceleration and emission regions.

Contribution

It introduces the use of offset polar cap dipole magnetic fields in pulsar light curve modeling and compares their effects with static and retarded vacuum fields, including the impact of different emission geometries.

Findings

01

Retarded vacuum dipole with outer gap model provides the best fit.

02

Offset polar cap dipole field favors smaller epsilon for constant emissivity.

03

Enhanced slot gap E-field improves light curve fits and aligns with independent studies.

Abstract

We performed geometric pulsar light curve modeling using static, retarded vacuum, and offset polar cap (PC) dipole $B$ -fields (the latter is characterized by a parameter $ϵ$ ), in conjunction with standard two-pole caustic (TPC) and outer gap (OG) emission geometries. The offset-PC dipole $B$ -field mimics deviations from the static dipole (which corresponds to $ϵ = 0$ ). In addition to constant-emissivity geometric models, we also considered a slot gap (SG) $E$ -field associated with the offset-PC dipole $B$ -field and found that its inclusion leads to qualitatively different light curves. Solving the particle transport equation shows that the particle energy only becomes large enough to yield significant curvature radiation at large altitudes above the stellar surface, given this relatively low $E$ -field. Therefore, particles do not always attain the radiation-reaction limit.…

Tables1

Table 1. Table 1: Best-fit parameters for each B 𝐵 B -field and geometric model combination

Combinations			Our Best-fit Parameters						Other Multi-wavelength Fits
Model	$ϵ$	$χ^{2}$	$α$	$ζ$	A	$Δ ϕ_{L}$	$Δ ξ_{B}^{2}$	$Δ ξ_{all}^{2}$	$α$	$ζ$	Reference
		( $\times 10^{5}$ )	(^∘)	(^∘)					(^∘)	(^∘)
Static dipole B-field:
TPC	…	0.819	$73_{- 2}^{+ 3}$	$45_{- 4}^{+ 4}$	1.3	0.55	0.00	108.75
OG	…	0.891	$64_{- 3}^{+ 5}$	$86_{- 1}^{+ 1}$	1.3	0.05	8.44	126.75
RVD B-field:
TPC	…	3.278	$54_{- 5}^{+ 5}$	$67_{- 3}^{+ 5}$	0.5	0.05	723.50	723.50
OG	…	0.384	$78_{- 1}^{+ 1}$	$69_{- 1}^{+ 2}$	1.3	0.00	0.00	0.00
Offset-PC dipole $B$ -field for constant $ϵ_{ν}$ :
TPC	0.00	0.819	$73_{- 2}^{+ 3}$	$45_{- 4}^{+ 4}$	1.3	0.55	0.00	108.75
	0.03	0.834	$73_{- 2}^{+ 2}$	$43_{- 5}^{+ 4}$	1.3	0.55	1.76	112.50
	0.06	0.867	$73_{- 2}^{+ 2}$	$42_{- 5}^{+ 5}$	1.3	0.55	5.63	120.75
	0.09	0.882	$73_{- 2}^{+ 1}$	$41_{- 5}^{+ 3}$	1.3	0.55	7.39	124.50
	0.12	1.000	$74_{- 3}^{+ 1}$	$42_{- 6}^{+ 3}$	1.4	0.55	21.22	154.00
	0.15	0.948	$73_{- 2}^{+ 1}$	$39_{- 5}^{+ 3}$	1.4	0.55	15.12	141.00
	0.18	0.969	$73_{- 3}^{+ 2}$	$37_{- 4}^{+ 4}$	1.3	0.55	17.58	146.25
Offset-PC dipole $B$ -field for variable $ϵ_{ν}$ :
SG	0.00	1.587	$21_{- 3}^{+ 3}$	$71_{- 1}^{+ 1}$	0.5	0.85	40.52	300.75
	0.03	1.627	$73_{- 1}^{+ 1}$	$17_{- 3}^{+ 4}$	0.7	0.55	43.96	310.75
	0.06	1.525	$72_{- 1}^{+ 2}$	$14_{- 1}^{+ 5}$	0.5	0.60	35.18	285.25
	0.09	1.452	$73_{- 1}^{+ 1}$	$17_{- 1}^{+ 3}$	0.6	0.55	28.90	267.00
	0.12	1.437	$74_{- 1}^{+ 1}$	$27_{- 7}^{+ 1}$	0.8	0.55	27.61	263.25
	0.15	1.116	$76_{- 1}^{+ 3}$	$48_{- 11}^{+ 15}$	0.7	0.55	0.00	183.00
	0.18	1.119	$75_{- 1}^{+ 2}$	$40_{- 4}^{+ 6}$	0.5	0.55	0.26	183.75
Offset-PC dipole $B$ -field for variable $ϵ_{ν}$ ( $100 E_{∥}$ ):
SG	0.00	0.581	$75_{- 1}^{+ 3}$	$51_{- 5}^{+ 2}$	1.1	0.55	0.00	49.27
	0.03	0.634	$75_{- 2}^{+ 2}$	$49_{- 5}^{+ 5}$	1.1	0.55	8.73	62.48
	0.06	0.698	$75_{- 3}^{+ 3}$	$49_{- 6}^{+ 5}$	1.1	0.55	19.39	78.61
	0.09	0.774	$75_{- 3}^{+ 3}$	$50_{- 9}^{+ 5}$	1.1	0.55	31.90	97.54
	0.12	0.789	$77_{- 3}^{+ 2}$	$54_{- 8}^{+ 2}$	1.1	0.55	34.42	101.36
	0.15	0.845	$77_{- 4}^{+ 2}$	$55_{- 14}^{+ 1}$	0.9	0.55	43.62	115.28
	0.18	0.834	$78_{- 2}^{+ 1}$	$55_{- 5}^{+ 1}$	0.8	0.55	41.80	112.51
RVM									53	59.5	1
X-ray torus										63.6 ${}_{- 0.05}^{+ 0.07}$	2
RVD and TPC									62–68	64	3
RVD and OG									75	64	3
RVD and Symmetric SG									44 ${}_{- 1}^{+ 4}$	54 ${}_{- 5}^{+ 1}$	4
RVD and Asymmetric SG									65 ${}_{- 2}^{+ 1}$	65.5 ${}_{- 1}^{+ 2}$	4
RVD and OG									88 ${}_{- 3}^{+ 2}$	66.5 ${}_{- 1}^{+ 1}$	4
FF and Symmetric SG									15 ${}_{- 1}^{+ 1}$	68.5 ${}_{- 1}^{+ 1}$	4
FF and Asymmetric SG									55 ${}_{- 20}^{+ 10}$	54.5 ${}_{- 14}^{+ 4}$	4
FF and OG									80 ${}_{- 1}^{+ 1}$	53 ${}_{- 1}^{+ 1}$	4
RVD and PC									3 ${}_{- 3}^{+ 2}$	4 ${}_{- 2}^{+ 2}$	5
RVD and SG									45 ${}_{- 2}^{+ 2}$	69 ${}_{- 2}^{+ 2}$	5
RVD and OG									71 ${}_{- 2}^{+ 2}$	83 ${}_{- 2}^{+ 2}$	5
RVD and OPC									56 ${}_{- 2}^{+ 2}$	77 ${}_{- 2}^{+ 2}$	5

Equations54

B_{ret}

B_{ret}

μ (t) = μ (sin α cos Ω t \hat{x} + sin α sin Ω t \hat{y} + cos α \hat{z}),

μ (t) = μ (sin α cos Ω t \hat{x} + sin α sin Ω t \hat{y} + cos α \hat{z}),

B_{OPCs}^{'}

B_{OPCs}^{'}

Δ r_{PC} ≃ R θ_{PC} [1 - θ_{PC}^{ϵ}],

Δ r_{PC} ≃ R θ_{PC} [1 - θ_{PC}^{ϵ}],

W^{'} (r^{'}, θ^{'})

W^{'} (r^{'}, θ^{'})

B_{EOD}^{'} (r^{'}, θ^{'})

B_{EOD}^{'} (r^{'}, θ^{'})

B_{OPCs}^{'} (r^{'}, θ^{'}, ϕ^{'})

B_{OPCs}^{'} (r^{'}, θ^{'}, ϕ^{'})

E_{∥, low}

E_{∥, low}

E_{∥, high}

E_{∥, high}

E_{∥, SG} ≃ E_{∥, low} exp [- (η - 1) / (η_{c} - 1)] + E_{∥, high} .

E_{∥, SG} ≃ E_{∥, low} exp [- (η - 1) / (η_{c} - 1)] + E_{∥, high} .

χ^{2}

χ^{2}

Y_{m}^{'} (ϕ_{L, i}) = \frac{Y _{m} ( ϕ _{L, i} )}{( k _{1} + ϵ _{0} )} A (k_{2} - BG) + BG \approx \frac{Y _{m} ( ϕ _{L, i} )}{k _{1}} k_{2},

Y_{m}^{'} (ϕ_{L, i}) = \frac{Y _{m} ( ϕ _{L, i} )}{( k _{1} + ϵ _{0} )} A (k_{2} - BG) + BG \approx \frac{Y _{m} ( ϕ _{L, i} )}{k _{1}} k_{2},

χ_{opt}^{2}

χ_{opt}^{2}

ξ^{2}

ξ^{2}

ξ_{opt}^{2}

ξ_{opt}^{2}

Δ ξ^{2}

Δ ξ^{2}

\overset{γ}{˙} = \overset{γ}{˙}_{gain} + \overset{γ}{˙}_{loss} = \frac{e E _{∥, new}}{m _{e} c} - \frac{2 e ^{2} γ _{e}^{4}}{3 ρ _{curv}^{2} m _{e} c} = \frac{1}{m _{e} c ^{2}} [ce E_{∥, new} - \frac{2 c e ^{2} γ _{e}^{4}}{3 ρ _{curv}^{2}}],

\overset{γ}{˙} = \overset{γ}{˙}_{gain} + \overset{γ}{˙}_{loss} = \frac{e E _{∥, new}}{m _{e} c} - \frac{2 e ^{2} γ _{e}^{4}}{3 ρ _{curv}^{2} m _{e} c} = \frac{1}{m _{e} c ^{2}} [ce E_{∥, new} - \frac{2 c e ^{2} γ _{e}^{4}}{3 ρ _{curv}^{2}}],

E_{CR} \sim 4 E_{∥, 4}^{3/4} ρ_{curv, 8}^{1/2} GeV,

E_{CR} \sim 4 E_{∥, 4}^{3/4} ρ_{curv, 8}^{1/2} GeV,

B^{'} (r^{'}, θ^{'}, ϕ^{'}) = B_{r}^{'} (r^{'}, θ^{'}, ϕ^{'}) \hat{r}^{'} + B_{θ}^{'} (r^{'}, θ^{'}, ϕ^{'}) \hat{θ}^{'} + B_{ϕ}^{'} (r^{'}, θ^{'}, ϕ^{'}) \hat{ϕ}^{'} .

B^{'} (r^{'}, θ^{'}, ϕ^{'}) = B_{r}^{'} (r^{'}, θ^{'}, ϕ^{'}) \hat{r}^{'} + B_{θ}^{'} (r^{'}, θ^{'}, ϕ^{'}) \hat{θ}^{'} + B_{ϕ}^{'} (r^{'}, θ^{'}, ϕ^{'}) \hat{ϕ}^{'} .

\begin{array}[]{ccrcrcr}{\hat{\rm{\mathbf{r}}}}^{\prime}&=&(\sin\theta^{\prime}\cos\phi^{\prime}){\hat{\rm{\mathbf{x}}}}^{\prime}&+&(\sin\theta^{\prime}\sin\phi^{\prime}){\hat{\rm{\mathbf{y}}}}^{\prime}&+&(\cos\theta^{\prime}){\hat{\rm{\mathbf{z}}}}^{\prime},\\ {\hat{\rm{\bm{\theta}}}}^{\prime}&=&(\cos\theta^{\prime}\cos\phi^{\prime}){\hat{\rm{\mathbf{x}}}}^{\prime}&+&(\cos\theta^{\prime}\sin\phi^{\prime}){\hat{\rm{\mathbf{y}}}}^{\prime}&-&(\sin\theta^{\prime}){\hat{\rm{\mathbf{z}}}}^{\prime},\\ {\hat{\rm{\bm{\phi}}}}^{\prime}&=&-(\sin\phi^{\prime}){\hat{\rm{\mathbf{x}}}}^{\prime}&+&(\cos\phi^{\prime}){\hat{\rm{\mathbf{y}}}}^{\prime}.&\end{array}

\begin{array}[]{ccrcrcr}{\hat{\rm{\mathbf{r}}}}^{\prime}&=&(\sin\theta^{\prime}\cos\phi^{\prime}){\hat{\rm{\mathbf{x}}}}^{\prime}&+&(\sin\theta^{\prime}\sin\phi^{\prime}){\hat{\rm{\mathbf{y}}}}^{\prime}&+&(\cos\theta^{\prime}){\hat{\rm{\mathbf{z}}}}^{\prime},\\ {\hat{\rm{\bm{\theta}}}}^{\prime}&=&(\cos\theta^{\prime}\cos\phi^{\prime}){\hat{\rm{\mathbf{x}}}}^{\prime}&+&(\cos\theta^{\prime}\sin\phi^{\prime}){\hat{\rm{\mathbf{y}}}}^{\prime}&-&(\sin\theta^{\prime}){\hat{\rm{\mathbf{z}}}}^{\prime},\\ {\hat{\rm{\bm{\phi}}}}^{\prime}&=&-(\sin\phi^{\prime}){\hat{\rm{\mathbf{x}}}}^{\prime}&+&(\cos\phi^{\prime}){\hat{\rm{\mathbf{y}}}}^{\prime}.&\end{array}

B^{'} (r^{'}, θ^{'}, ϕ^{'}) = B_{x}^{'} (r^{'}, θ^{'}, ϕ^{'}) \hat{x}^{'} + B_{y}^{'} (r^{'}, θ^{'}, ϕ^{'}) \hat{y}^{'} + B_{z}^{'} (r^{'}, θ^{'}, ϕ^{'}) \hat{z}^{'} .

B^{'} (r^{'}, θ^{'}, ϕ^{'}) = B_{x}^{'} (r^{'}, θ^{'}, ϕ^{'}) \hat{x}^{'} + B_{y}^{'} (r^{'}, θ^{'}, ϕ^{'}) \hat{y}^{'} + B_{z}^{'} (r^{'}, θ^{'}, ϕ^{'}) \hat{z}^{'} .

\begin{array}[]{ccrcrcr}B_{x}^{\prime}(r^{\prime},\theta^{\prime},\phi^{\prime})&=&B_{r}^{\prime}\sin\theta^{\prime}\cos\phi^{\prime}&+&B_{\theta}^{\prime}\cos\theta^{\prime}\cos\phi^{\prime}&-&B_{\phi}^{\prime}\sin\phi^{\prime},\\ B_{y}^{\prime}(r^{\prime},\theta^{\prime},\phi^{\prime})&=&B_{r}^{\prime}\sin\theta^{\prime}\sin\phi^{\prime}&+&B_{\theta}^{\prime}\cos\theta^{\prime}\sin\phi^{\prime}&+&B_{\phi}^{\prime}\cos\phi^{\prime},\\ B_{z}^{\prime}(r^{\prime},\theta^{\prime},\phi^{\prime})&=&B_{r}^{\prime}\cos\theta^{\prime}&-&B_{\theta}^{\prime}\sin\theta^{\prime}.&\end{array}

\begin{array}[]{ccrcrcr}B_{x}^{\prime}(r^{\prime},\theta^{\prime},\phi^{\prime})&=&B_{r}^{\prime}\sin\theta^{\prime}\cos\phi^{\prime}&+&B_{\theta}^{\prime}\cos\theta^{\prime}\cos\phi^{\prime}&-&B_{\phi}^{\prime}\sin\phi^{\prime},\\ B_{y}^{\prime}(r^{\prime},\theta^{\prime},\phi^{\prime})&=&B_{r}^{\prime}\sin\theta^{\prime}\sin\phi^{\prime}&+&B_{\theta}^{\prime}\cos\theta^{\prime}\sin\phi^{\prime}&+&B_{\phi}^{\prime}\cos\phi^{\prime},\\ B_{z}^{\prime}(r^{\prime},\theta^{\prime},\phi^{\prime})&=&B_{r}^{\prime}\cos\theta^{\prime}&-&B_{\theta}^{\prime}\sin\theta^{\prime}.&\end{array}

\begin{array}[]{ccl}x^{\prime}&=&r^{\prime}\sin\theta^{\prime}\cos\phi^{\prime},\\ y^{\prime}&=&r^{\prime}\sin\theta^{\prime}\sin\phi^{\prime},\\ z^{\prime}&=&r^{\prime}\cos\theta^{\prime}.\end{array}

\begin{array}[]{ccl}x^{\prime}&=&r^{\prime}\sin\theta^{\prime}\cos\phi^{\prime},\\ y^{\prime}&=&r^{\prime}\sin\theta^{\prime}\sin\phi^{\prime},\\ z^{\prime}&=&r^{\prime}\cos\theta^{\prime}.\end{array}

B^{'} (x^{'}, y^{'}, z^{'}) = B_{x}^{'} (x^{'}, y^{'}, z^{'}) \hat{x}^{'} + B_{y}^{'} (x^{'}, y^{'}, z^{'}) \hat{y}^{'} + B_{z}^{'} (x^{'}, y^{'}, z^{'}) \hat{z}^{'} .

B^{'} (x^{'}, y^{'}, z^{'}) = B_{x}^{'} (x^{'}, y^{'}, z^{'}) \hat{x}^{'} + B_{y}^{'} (x^{'}, y^{'}, z^{'}) \hat{y}^{'} + B_{z}^{'} (x^{'}, y^{'}, z^{'}) \hat{z}^{'} .

\left[\begin{array}[]{c}B_{x}\\ B_{y}\\ B_{z}\end{array}\right]=\left[\begin{array}[]{ccc}\cos\alpha&0&\sin\alpha\\ 0&1&0\\ -\sin\alpha&0&\cos\alpha\end{array}\right]\left[\begin{array}[]{c}B_{x}^{\prime}\\ B_{y}^{\prime}\\ B_{z}^{\prime}\end{array}\right].

\left[\begin{array}[]{c}B_{x}\\ B_{y}\\ B_{z}\end{array}\right]=\left[\begin{array}[]{ccc}\cos\alpha&0&\sin\alpha\\ 0&1&0\\ -\sin\alpha&0&\cos\alpha\end{array}\right]\left[\begin{array}[]{c}B_{x}^{\prime}\\ B_{y}^{\prime}\\ B_{z}^{\prime}\end{array}\right].

\begin{array}[]{ccrcrcr}{\rm{\mathbf{B}}}(x^{\prime},y^{\prime},z^{\prime})&=&B_{x}(x^{\prime},y^{\prime},z^{\prime}){\hat{\rm{\mathbf{x}}}}^{\prime}&+&B_{y}(x^{\prime},y^{\prime},z^{\prime}){\hat{\rm{\mathbf{y}}}}^{\prime}&+&B_{z}(x^{\prime},y^{\prime},z^{\prime}){\hat{\rm{\mathbf{z}}}}^{\prime}.\end{array}

\begin{array}[]{ccrcrcr}{\rm{\mathbf{B}}}(x^{\prime},y^{\prime},z^{\prime})&=&B_{x}(x^{\prime},y^{\prime},z^{\prime}){\hat{\rm{\mathbf{x}}}}^{\prime}&+&B_{y}(x^{\prime},y^{\prime},z^{\prime}){\hat{\rm{\mathbf{y}}}}^{\prime}&+&B_{z}(x^{\prime},y^{\prime},z^{\prime}){\hat{\rm{\mathbf{z}}}}^{\prime}.\end{array}

B (x, y, z) = B_{x} (x, y, z) \hat{x} + B_{y} (x, y, z) \hat{y} + B_{z} (x, y, z) \hat{z} .

B (x, y, z) = B_{x} (x, y, z) \hat{x} + B_{y} (x, y, z) \hat{y} + B_{z} (x, y, z) \hat{z} .

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Full text

The effect of an offset polar cap dipolar magnetic field on the modeling of the Vela pulsar’s $\gamma$ -ray light curves

M. Barnard11affiliation: Centre for Space Research, North-West University, Potchefstroom Campus, Private Bag X6001, Potchefstroom 2520, South Africa 33affiliation: Formerly Monica Breed and C. Venter11affiliation: Centre for Space Research, North-West University, Potchefstroom Campus, Private Bag X6001, Potchefstroom 2520, South Africa

A. K. Harding22affiliation: Astrophysics Science Division, NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA

Abstract

We performed geometric pulsar light curve modeling using static, retarded vacuum, and offset polar cap (PC) dipole $B$ -fields (the latter is characterized by a parameter $\epsilon$ ), in conjunction with standard two-pole caustic (TPC) and outer gap (OG) emission geometries. The offset-PC dipole $B$ -field mimics deviations from the static dipole (which corresponds to $\epsilon=0$ ). In addition to constant-emissivity geometric models, we also considered a slot gap (SG) $E$ -field associated with the offset-PC dipole $B$ -field and found that its inclusion leads to qualitatively different light curves. Solving the particle transport equation shows that the particle energy only becomes large enough to yield significant curvature radiation at large altitudes above the stellar surface, given this relatively low $E$ -field. Therefore, particles do not always attain the radiation-reaction limit. Our overall optimal light curve fit is for the retarded vacuum dipole field and OG model, at an inclination angle $\alpha=78{{}_{-1}^{+1}}^{\circ}$ and observer angle $\zeta=69{{}_{-1}^{+2}}^{\circ}$ . For this $B$ -field, the TPC model is statistically disfavored compared to the OG model. For the static dipole field, neither model is significantly preferred. We found that smaller values of $\epsilon$ are favored for the offset-PC dipole field when assuming constant emissivity, and larger $\epsilon$ values favored for variable emissivity, but not significantly so. When multiplying the SG $E$ -field by a factor of 100, we found improved light curve fits, with $\alpha$ and $\zeta$ being closer to best fits from independent studies, as well as curvature radiation reaction at lower altitudes.

gamma rays: stars — pulsars: individual (PSR J0835 $-$ 4510) — stars: magnetic field — stars: neutron

1 INTRODUCTION

The field of $\gamma$ -ray pulsars has been revolutionized by the launch of the Fermi Large Area Telescope (LAT; Atwood et al., 2009). Over the past eight years, Fermi has detected over 200 $\gamma$ -ray pulsars and has furthermore measured their light curves and spectral characteristics in unprecedented detail. Fermi’s Second Pulsar Catalog (2PC; Abdo et al., 2013) describes the properties of some 117 of these pulsars in the energy range 100 MeV $-$ 100 GeV. This catalog includes the Vela pulsar (Abdo et al., 2009), one of the brightest persistent sources in the GeV sky. Recently, H.E.S.S. detected pulsed emission from the Vela pulsar in the 20 $-$ 120 GeV range (A. Abramowski et al. 2016, in preparation). This followed the detection of the Crab pulsar by VERITAS and MAGIC in the very-high-energy ( $>$ 100 GeV) regime (and now possibly up to 1 TeV; Aleksić et al., 2011, 2012; Aliu et al., 2011). In this paper, we will focus on the GeV band light curves of the Vela pulsar.

Despite the major advances made after nearly 50 years since the discovery of the first pulsar (Hewish et al., 1968), many questions still remain regarding the electrodynamical character of the pulsar magnetosphere, including details of the particle acceleration and pair production, current closure, and radiation of a complex multi-wavelength spectrum. Physical emission models such as the slot gap (SG; Muslimov & Harding, 2003) and outer gap (OG; Romani & Yadigaroglu, 1995) fall short of explaining these global magnetospheric characteristics. More recent developments include the global magnetospheric properties. One example is the force-free (FF) inside and dissipative outside (FIDO) model (Kalapotharakos & Contopoulos, 2009; Kalapotharakos et al., 2014) that assumes FF electrodynamical conditions (infinite plasma conductivity, $\sigma_{\rm c}\rightarrow\infty$ ) inside the light cylinder and dissipative conditions (finite $\sigma_{\rm c}$ ) outside. The wind models of, e.g., Pétri & Dubus (2011) provide an alternative picture where dissipation takes place outside the light cylinder.

Although much progress has been made using the physical models, geometric light curve modeling still presents a crucial avenue for probing the pulsar magnetosphere in the context of traditional pulsar models, as these emission geometries may be used to constrain the pulsar geometry (inclination angle $\alpha$ and the observer viewing angle $\zeta$ with respect to the spin axis $\bm{\Omega}$ ), as well as the $\gamma$ -ray emission region’s location and extent. This may provide vital insight into the boundary conditions and help constrain the accelerator geometry of next-generation full radiation models. Geometric light curve modeling has been performed by, e.g., Dyks et al. (2004a), Venter et al. (2009), Watters et al. (2009), Johnson et al. (2014), and Pierbattista et al. (2015) using standard pulsar emission geometries, including a two-pole caustic (TPC, of which the SG is its physical representation; Dyks & Rudak, 2003), OG, and pair-starved polar cap (PC) (Harding et al., 2005) geometry.

A notable conclusion from the 2PC was that the spectra and light curves of both the millisecond pulsar (MSP) and young pulsar populations show remarkable similarities, pointing to a common radiation mechanism and emission geometry (tied to the $B$ -field structure). The assumed $B$ -field structure is essential for predicting the light curves seen by the observer using geometric models, since photons are expected to be emitted tangentially to the local $B$ -field lines in the corotating pulsar frame (Daugherty & Harding, 1982). Even a small difference in the magnetospheric structure will therefore have an impact on the light curve predictions. For all of the above geometric models, the most commonly employed $B$ -field has been the retarded vacuum dipole (RVD) solution first obtained by Deutsch (1955). However, other solutions also exist. One example is the static dipole (non-rotating) field, a special case of the RVD (rotating) field (Dyks & Harding, 2004). Bai & Spitkovsky (2010b) furthermore modeled high-energy (HE) light curves in the context of OG and TPC models using an FF $B$ -field geometry (assuming a plasma-filled magnetosphere), proposing a separatrix layer model close to the last open field line (tangent to the light cylinder at radius $R_{\rm LC}=c/\Omega$ where the corotation speed equals the speed of light $c$ , with $\Omega$ the angular speed), which extends from the stellar surface up to and beyond the light cylinder. In addition, the annular gap model of Du et al. (2010), which assumes a static dipole field, has been successful in reproducing the main characteristics of the $\gamma$ -ray light curves of three MSPs. This model does not, however, attempt to replicate the non-zero phase offsets between the $\gamma$ -ray and radio profiles.

In this paper, we investigate the impact of different magnetospheric structures (i.e., the static dipole and RVD) on the predicted $\gamma$ -ray pulsar light curves. Additionally, we incorporate an offset-PC dipole $B$ -field solution (Harding & Muslimov, 2011a, b) into our geometric modeling code. This analytic heuristic model mimics deviations from the static dipole such as what occur in complex solutions, e.g., dissipative (e.g., Kalapotharakos et al., 2012; Li et al., 2012; Tchekhovskoy et al., 2013; Li, 2014) and FF (e.g., Contopoulos et al., 1999) fields. These complex $B$ -fields usually only have numerical solutions, which are limited by the resolution of the spatial grid. Hence, it is simpler to investigate the main effects of these structures using analytical approximations such as the offset-PC dipole solution. In combination with the different $B$ -field solutions mentioned above, we assume standard TPC and OG emission geometries.

Geometric models assume constant emissivity $\epsilon_{\nu}$ in the rotational frame. We have also incorporated an SG $E$ -field associated with the offset-PC dipole $B$ -field (making this latter case an emission model), which allows us to calculate the $\epsilon_{\nu}$ in the acceleration region in the corotating frame from first principles. We have only considered the TPC (assuming uniform $\epsilon_{\nu}$ ) and SG (assuming variable $\epsilon_{\nu}$ as modulated by the $E$ -field) models for the offset-PC dipole $B$ -field, since we do not have $E$ -field expressions available for the OG model within the context of an offset-PC dipole $B$ -field. The fact that we have an $E$ -field solution enables us to solve the particle transport equation on each $B$ -field line, yielding the particle energy (Lorentz factor $\gamma_{\rm e}$ ) as a function of position. We can then use this factor to test whether the particle reaches the curvature radiation reaction (CRR) limit, i.e., where acceleration balances curvature radiation losses.

We have implemented a chi-squared ( $\chi^{2}$ ) method to search the multivariate solution space for optimal model parameters when we compare our predicted model light curves with Fermi LAT data for the Vela pulsar. In this way, we are able to determine which $B$ -field and geometric model combination yields the best light curve solution, how the different light curve predictions compare with each other, and which pulsar geometry ( $\alpha$ , $\zeta$ ) is optimal (Breed et al., 2014, 2015a, 2015b).

In Section 2 we describe the geometric pulsar models and $B$ -field structures we considered. Section 3 details the implementation of an offset-PC dipole $B$ -field in our code. We also discuss the implementation of the associated SG $E$ -field and the matching of the low-altitude and high-altitude solutions using a matching parameter (scaled radius) $\eta_{\rm c}$ . We briefly describe the implementation of a $\chi^{2}$ method in order to find the best-fit ( $\alpha$ , $\zeta$ ) for the different model combinations. In Section 4, we present our solution of the transport equation for the offset-PC dipole $B$ -field, as well as our light curve predictions and best-fit ( $\alpha$ , $\zeta$ ) contours for the Vela pulsar. In Section 5 we investigate the effects of lowering the minimum photon energy as well as multiplying the $E$ -field by a factor of 100 when constructing model light curves, before we compare, in Section 6, our results to previous multi-wavelength studies from other works. Our conclusions follow in Section 7. Given the length of this paper, we summarize our main conclusions here:

The SG $E$ -field is relatively low. We thus find that the particle energy only becomes large enough to yield significant curvature radiation at large altitudes above the stellar surface, so that particles do not always attain the radiation-reaction limit. 2. 2.

Our overall optimal light curve fit is for the RVD field and OG model; for this $B$ -field, the TPC model is statistically disfavored. 3. 3.

For the static dipole field, neither OG nor TPC model is significantly preferred. 4. 4.

We found that a smaller PC offset is favored for the offset-PC dipole field when assuming constant emissivity, and a larger PC offset favored for variable emissivity, but not significantly so. 5. 5.

When multiplying the SG $E$ -field by a factor of 100, we found improved light curve fits, with $\alpha$ and $\zeta$ being closer to best fits from independent studies, as well as CRR at lower altitudes.

2 GEOMETRIC MODEL DESCRIPTION

2.1 TPC and OG Geometries

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_{\rm max}{\simeq}0.8{R_{\rm LC}}$ (Muslimov & Harding, 2004). This was later extended to $R_{\rm max}{\simeq}1.2{R_{\rm 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, b) 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_{\rm 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_{\parallel}\neq 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.

2.2 $B$ -field Structures

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, b). 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}\rightarrow{0}$ ) and FF ( $\sigma_{\rm c}\rightarrow\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 $\gamma$ -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., , ; , ; , ; , ; , ). The general expression for this RVD field in the rotation frame (where $\hat{{\mathbf{z}}}||\bm{\Omega}$ ) is given by

[TABLE]

with the magnetic moment

[TABLE]

where $\dot{\bm{\mu}}(t)$ and $\ddot{\bm{\mu}}(t)$ are the first and second time derivatives of the magnetic moment, and ( $\hat{\bf{x}}$ , $\hat{\bf{y}}$ , $\hat{\bf{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_{\rm LC}\rightarrow 0$ , the retarded field simplifies to the non-aligned static dipole $B$ -field ( $\alpha\neq 0$ ). For the static dipole the field lines are symmetric about the $\bm{\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 ${\bm{\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_{\rm 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, b) 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 to 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{{\mathbf{z}}}^{\prime}\parallel{\bm{\mu}}$ ) is as follows (Harding & Muslimov, 2011b)

[TABLE]

where $\mu^{\prime}=B_{0}R^{3}$ is the magnetic moment, $R$ the stellar radius, $B_{0}$ 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^{\prime}_{0}=0$ or $\phi^{\prime}_{0}=\pi$ the offset is in the $x^{\prime}$ direction (i.e., along the $x^{\prime}$ -axis). If $\phi^{\prime}_{0}=\pi/2$ or $\phi^{\prime}_{0}=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 $\epsilon$ (related to the magnitude of the shift of the PC from the magnetic axis) and $\phi_{0}^{\prime}$ (we choose $\phi^{\prime}_{0}=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

[TABLE]

where $\theta_{\rm PC}=(\Omega{R}/c)^{1/2}$ is the standard half-angle of the PC. This effective shift of the PCs is a fraction of $\theta_{\rm 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_{\rm 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

[TABLE]

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=-\nabla W$ :

[TABLE]

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

[TABLE]

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, b) derived these perturbed components of the distorted magnetic field while satisfying the solenoidality condition $\nabla\cdot B=0$ .

2.3 Geometric Modeling Code

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_{\rm 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_{\rm ovc}$ and $l_{\rm ovc}$ ). After the footpoints of the field lines on a ( $r_{\rm ovc}$ , $l_{\rm 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 $\zeta$ , i.e., a sky map is formed by plotting the bin contents (divided by the solid angle subtended by each bin) for a given $\alpha$ , and it is therefore a projection of the radiation beam. To simulate light curves, one chooses a sky map corresponding to a fixed $\alpha$ , then fixes $\zeta$ 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_{\rm 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_{\rm LC}$ ). This implies a revised aberration formula, which we have implemented in our code.

3 IMPLEMENTATION OF AN OFFSET-PC DIPOLE $B$ -FIELD AND ASSOCIATED SG $E$ -FIELD AND A SEARCH FOR THE BEST-FIT LIGHT CURVE

3.1 Transformation of a $B$ -field from the Magnetic to the Rotational Frame

Since the offset-PC dipole field is specified in the magnetic frame ( $\hat{\rm{\mathbf{z}}}^{\prime}\parallel\bm{\mu}^{\prime}$ ; Harding & Muslimov, 2011b), it is necessary to transform this solution to the (corotating) rotational frame ( $\hat{\rm{\mathbf{z}}}\parallel{\bm{\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.

3.2 Finding the PC Rim and Extending the Range of $\epsilon$

The object is to find the polar angle $\theta_{*}$ at each azimuthal angle $\phi$ at the footpoints of the last open $B$ -field lines, lying within a bracket $\theta_{\rm min}<\theta_{*}<\theta_{\rm 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$ ( $\epsilon\lesssim{0.05-0.1}$ , depending on $\alpha$ ). We improve the range of $\epsilon$ by varying the colatitude parameters $\theta_{\rm min}$ and $\theta_{\rm 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_{\rm LC}$ ). We obtain a progressively larger range of $\epsilon$ upon decreasing $\theta_{\rm min}$ and increasing $\theta_{\rm max}$ . We find a maximum $\epsilon=0.18$ valid for the full range of $\alpha$ . 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_{\rm min}$ and found (by trial and error) that the following linear equation $\theta_{\rm min}=[(-31/18)\epsilon+0.6]\theta_{\rm PC}$ , for a fixed $\theta_{\rm max}=2.0$ , resulted in $\theta_{\rm min}$ that yielded maximum values for $\epsilon$ .

We illustrate the PC shape for a few cases of $\alpha$ and $\epsilon$ in Figure 1. We plot the PC rims ( $r_{\rm ovc}=1$ ) in the $x^{\prime}-y^{\prime}$ plane (in the magnetic frame, in units of $r_{\rm PC}$ ), assuming that the $\bm{\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 $\alpha$ in the range $10^{\circ}-90^{\circ}$ with increments of $10^{\circ}$ . For each $\alpha$ we plot the PC shape for $\epsilon$ values of [math] (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 $\epsilon$ is increased. As $\alpha$ and $\epsilon$ are increased the PC shape changes considerably. For larger $\epsilon$ the PC offset is larger along the $-x^{\prime}$ -axis (in the direction of “unfavorably curved” $B$ -field lines). Also, as $\alpha$ increases for each $\epsilon$ 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 $\epsilon$ increases.

3.3 Incorporating an SG $E$ -field

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)

[TABLE]

with $\mathcal{E}_{0}=(\Omega{R}/c)^{2}(B_{r}/B)B_{0}$ , $B_{r}$ the radial $B$ -field component, $\nu_{\rm SG}\equiv({1/4})\Delta\xi_{\rm SG}^{2}$ , and $\Delta\xi_{\rm SG}$ the colatitudinal gap width in units of dimensionless colatitude $\xi=\theta/\theta_{\rm PC}$ . Also, $x=r/R_{\rm LC}$ is the normalized radial distance in units of $R_{\rm LC}$ . Here, $\kappa\approx{0.15I_{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 cm2, $I$ the moment of inertia, $R_{6}=R/10^{6}$ cm, $\eta=r/R$ the dimensionless radial coordinate in units of $R$ , $e_{\rm 1A}=1+{a}(\eta^{3}-1)/3$ , $e_{\rm 2A}=(1+3a)\eta^{(1+a)/2}-2a$ and $e_{\rm 3A}=[({5-3a})/{\eta^{(5-a)/2}}]+2a$ . The magnetic azimuthal angle $\phi_{\rm PC}$ is defined for usage with the $E$ -field, being $\pi$ out of phase with $\phi^{\prime}$ (one chooses the negative $x$ -axis toward $\bm{\Omega}$ to coincide with $\phi_{\rm 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_{\ast}$ is the dimensionless colatitude labeling the gap field lines (defined such that $\xi_{\ast}=0$ corresponds to the field line in the middle of the gap and $\xi_{\ast}=1$ at the boundaries; Muslimov & Harding, 2003).

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

[TABLE]

with $f(\eta)\sim{1+0.75y+0.6y^{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.16y^{2}-0.5(\kappa/\epsilon_{\rm g}^{3})y^{3}(1-0.25y-0.21y^{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}=2GM/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_{\rm LC}=R_{\rm 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_{\rm LC}$ (i.e., over the entire length of the gap) is constructed as follows (see Equation (59) of Muslimov & Harding, 2004)

[TABLE]

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

3.4 Determining the Matching Parameter $\eta_{\rm c}$

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_{\rm 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 $\eta$ .

We solve for the matching parameter in the following way. First, we calculate $E_{\parallel,\rm low}$ , which is independent of $\eta_{\rm c}$ , along the $B$ -field. If $-E_{\parallel,\rm low}<0$ for all $\eta$ , it will never intersect with $E_{\parallel,\rm high}$ and we set $\eta_{\rm c}=1.1$ , thereby basically using $E_{\parallel,\rm SG}\approx E_{\parallel,\rm high}$ . Second, we step through $\eta_{\rm c}$ (in the range $1.1-5.1$ ), calculating $E_{\parallel,\rm SG}$ and $E_{\parallel,\rm high}$ as well as the ratio $S_{i}=S(\eta_{i})=E_{\parallel,\rm SG}(\eta_{i})/E_{\parallel,\rm low}(\eta_{i})$ for $i=1,..,N$ at different radii $\eta_{i}$ . If $S_{i}>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,\rm low}>-E_{\parallel,\rm high}$ (i.e., we basically fit $E_{\parallel,\rm SG}$ to $E_{\parallel,\rm low}$ when $-E_{\parallel,\rm low}>-E_{\parallel,\rm 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_{\rm cut}>\eta_{\rm LC}$ (i.e., $E_{\parallel,\rm low}$ and $E_{\parallel,\rm 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,\rm low}$ does not decrease as rapidly (e.g., as in Figure 2(b)) we find reasonable solutions. We note that $E_{\parallel,\rm SG}$ (referred to as $E_{\parallel,\rm old}$ in Figure 2) produces a bump when $-E_{\parallel,\rm low}$ decreases more rapidly. To circumvent this problem we test whether $-E_{\parallel,\rm SG}<-E_{\parallel,\rm high}$ and in this case we use the intersection radius $\eta_{\rm cut}$ of $E_{\parallel,\rm low}$ and $E_{\parallel,\rm high}$ , rather than $\eta_{\rm c}$ , to match our solutions (calling this new solution $E_{\parallel,\rm new}$ ; see Figure 2(c)). We lastly observe that for $\phi_{\rm PC}=\pi$ (on “unfavorably curved” field lines) for larger $\alpha$ , the $-E_{\parallel,\rm low}$ field changes sign resulting in a small $\eta_{\rm c}=\eta_{\rm cut}=1.7$ value (Figure 2(d)).

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}>\eta_{\rm LC}$ .) Since the $E$ -field solutions have an $x^{a}=x^{\epsilon\cos(\phi^{\prime}-\phi^{\prime}_{0})}=x^{-\epsilon\cos\phi_{\rm PC}}$ factor dependence, a larger (non-zero) offset results in different matching contours vs. 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_{\rm 1A}$ (with an $\epsilon$ dependence) being initially larger at low $\eta$ for $\phi_{\rm PC}=0$ than for $\phi_{\rm PC}=\pi$ ( $x^{a}$ dominates), but rapidly decreasing with $\eta$ ( $e_{\rm 1A}$ dominates), leading to a lower value of $\eta_{\rm c}$ for $\phi_{\rm 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 $\eta$ and eventually become lower than the second. In Figure 3 one can see that there is no $\phi_{\rm PC}$ -dependence for $\alpha=0$ , which is not the case for Figure 4. For a slightly larger $\alpha$ the second terms in Equations (7) and (8) start to contribute to the radiation. This is due to the $\sin\alpha$ term with an $\epsilon$ 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>20^{\circ}$ (Figure 4) the second term $\sim\cos\phi_{\rm PC}$ starts to dominate and thus we find solutions of $\eta_{\rm c}\sim{5.1}$ for $\phi_{\rm PC}\simeq 0$ and systematically smaller solutions for $\phi_{\rm PC}\simeq\pi$ as $\alpha$ increases and the second term $\propto\sin\alpha$ becomes increasingly important (in both cases of $\epsilon$ ). At $\alpha=90^{\circ}$ we obtain the same solution as in Figure 3 where the second term dominates (for this case $-E_{\parallel,\rm SG}<0$ for all $\eta$ , since the Goldreich $-$ Julian charge density $\rho_{\rm GJ}$ becomes positive). We note that the $\eta_{\rm c}$ distribution reflects two symmetries (one about $\phi_{\rm PC}=\pi$ and one about $\xi=0.975$ , i.e., $\xi_{\ast}=0$ , given our gap boundaries): that of the $\cos\phi_{\rm PC}$ term and that of the $(1-\xi_{\ast}^{2})$ term in the $E_{\parallel}$ solutions. After solving for $\eta_{\rm c}$ , we can solve the particle transport equation along each $B$ -field line (see Section 4.1).

3.5 Chi-squared Fitting Method

We apply a standard $\chi^{2}$ statistical fitting technique to assist us in objectively finding the pulsar geometry ( $\alpha$ , $\zeta$ ) which best describes the observed $\gamma$ -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

[TABLE]

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,...,N_{{\rm bins}}$ , with $N_{\rm 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_{{\rm 2}}$ by using the following expression:

[TABLE]

with $k_{1}={\rm max}(Y_{{\rm m}}(\phi_{{\rm L},i}))$ , $k_{2}={\rm 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 $\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 $\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 $\Delta\phi_{\rm L}$ ). Importantly, we note that we have not changed the relative position (the radio-to- $\gamma$ phase lag $\delta$ ), since this is a crucial model prediction. The radio and $\gamma$ -ray emission regions are tied to the same underlying $B$ -field structure, and $\delta$ therefore reflects important physical conditions (or model assumptions) such as a difference in emission heights of the radio and $\gamma$ -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^{\circ}$ resolution), $0.5<A<1.5$ with $0.1$ resolution, and $0<\Delta\phi_{\rm L}<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 $\alpha$ and $\zeta$ . Once we have iterated over the entire parameter space ( $\alpha$ , $\zeta$ , $A$ , $\Delta\phi_{\rm L}$ ), we obtain a global minimum value for $\chi^{2}$ (also called the optimal $\chi^{2}$ ):

[TABLE]

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_{{\rm opt}}^{2}/N_{\rm dof}\gg 1$ . We therefore need to rescale (to compensate for the uncertainty in $\sigma_{{\rm m,}i}$ ) the $\chi^{2}$ values by $\chi_{{\rm opt}}^{2}$ and multiply by the number of degrees of freedom $N_{\rm dof}$ (the difference between $N_{\rm bins}$ and number of free parameters). The scaled $\chi^{2}$ is presented by (Pierbattista et al., 2015):

[TABLE]

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

[TABLE]

with $\xi_{\rm opt}^{2}/N_{\rm dof}=\xi_{\rm opt,\nu}^{2}=1$ the reduced $\xi_{\rm 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_{{\rm opt}}^{2}$ and the $\xi^{2}$ of the alternative models:

[TABLE]

The confidence intervals can be estimated by reading the $\Delta\xi^{2}$ (i.e., $\Delta\xi_{1\sigma,\mu_{{\rm dof}}}^{2}$ , $\Delta\xi_{2\sigma,\mu_{{\rm dof}}}^{2}$ , and $\Delta\xi_{3\sigma,\mu_{{\rm dof}}}^{2}$ ) values from a standard $\chi^{2}$ table for the specified confidence interval at $\mu_{{\rm dof}}=2$ (corresponding to the two-dimensional ( $\alpha,\zeta$ ) grid Lampton et al., 1976). Using these values for $\Delta\xi^{2}$ and $\xi_{{\rm opt}}^{2}=N_{{\rm dof}}$ , we can determine $\xi^{2}=\xi^{2}_{\rm opt}+\Delta\xi^{2}=N_{\rm dof}+\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 ${\rm log_{10}\xi^{2}}$ on an ( $\alpha$ , $\zeta$ ) grid, with a minimum value of ${\rm log_{10}}\xi_{{\rm opt}}^{2}={\rm log_{10}}(N_{{\rm dof}})=1.98$ (corresponding to the best-fit solution by construction, i.e., after rescaling, with $N_{\rm dof}=100-4=96$ in our study). The best-fit solution is therefore positioned at $\xi_{{\rm opt}}^{2}=96$ and enclosed by the confidence contours with values of $\xi_{1\sigma,\mu_{{\rm dof}}}^{2}=96+2.30$ , $\xi_{2\sigma,\mu_{{\rm dof}}}^{2}=96+6.17$ , and $\xi_{3\sigma,\mu_{{\rm dof}}}^{2}=96+11.8$ (see Equation (15); Press et al., 1992). We determine errors on $\alpha$ and $\zeta$ for the best-fit solution of each $B$ -field and model combination using the $3\sigma$ interval connected contours. We choose errors of $1^{\circ}$ for cases when the errors were smaller than one degree (given a model resolution of $1^{\circ}$ ). 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.

4 RESULTS

4.1 Solving the Particle Transport Equation

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

[TABLE]

with $\dot{\gamma}_{\rm gain}$ the gain (acceleration) rate, $\dot{\gamma}_{\rm 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 $\log_{10}$ of $-E_{\parallel,{\rm high}}$ (solid cyan line), $-E_{\parallel,{\rm low}}$ (solid blue line), the general $-E_{\parallel,{\rm SG}}$ field (using $\eta_{\rm c}$ as the matching parameter; $-E_{\parallel,{\rm 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_{\rm cut}$ ; $-E_{\parallel,{\rm new}}$ , dashed red line), $\dot{\gamma}_{\rm gain}$ (solid yellow line), $\dot{\gamma}_{\rm 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_{\rm 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,{\rm low}}\propto{x^{a}}[1+a(\eta^{3}-1)/3]$ , where $a=-\epsilon\cos\phi_{\rm PC}$ . The sign of $x^{a}$ stays the same but at small values of $\eta$ , $x^{a}\approx{3}$ for $\phi_{\rm PC}=0$ (top panel) and $x^{a}\approx{1/3}$ for $\phi_{\rm PC}=\pi$ (bottom panel). This explains why $-E_{\parallel,{\rm low}}$ is higher for $\phi_{\rm PC}=0$ than for $\phi_{\rm PC}=\pi$ (for $\epsilon\neq 0$ ) at low $\eta$ . In the case of $\phi_{\rm PC}=\pi/2$ (middle panel), $a=0$ and the values of $E_{\parallel,{\rm low}}$ are very nearly the same for both $\epsilon=0$ and $\epsilon\neq 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_{\rm PC}=0$ ). The term $[1+a(\eta^{3}-1)/3]$ changes sign beyond some $\eta$ for $\phi_{\rm PC}=0$ , explaining the behavior of $-E_{\parallel,{\rm 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_{\rm PC}=0$ or $\phi_{\rm PC}=2\pi$ . This is contrary to the case of $\phi_{\rm PC}=\pi$ (“unfavorably curved” $B$ -field lines) where the sign stays the same for all $\eta$ . Correspondingly, $\eta_{\rm c}\approx 5$ around $\phi_{\rm PC}\neq 0$ in Figure 4 for $\alpha\approx 0^{\circ}$ . For $-E_{\parallel,{\rm high}}\propto{x^{a}[1+2\eta/\eta_{\rm LC}]}$ , we note that for $\phi_{\rm PC}=(0,\pi/2,\pi)$ , $x^{a}$ changes from the following low to high values: $(3\rightarrow 1,1\rightarrow 1,1/3\rightarrow 1)$ (for $\epsilon\neq 0$ , otherwise $x^{a}=1$ ) while $[1+2\eta/\eta_{\rm LC}]$ changes from $1\rightarrow 3$ (the number to the left of $\rightarrow$ is valid for $\eta\ll 1$ and the number to the right of $\rightarrow$ is valid for $\eta\gg 1$ ). The combined effect of these two terms is such that $-E_{\parallel,{\rm high}}\propto(3\rightarrow 3,1\rightarrow 3,1/3\rightarrow 3)$ for $\phi_{\rm PC}=(0,\pi/2,\pi)$ and $\epsilon\neq 0$ , and $-E_{\parallel,{\rm high}}\propto\rightarrow 3$ for $\epsilon=0$ , for all $\phi_{\rm PC}$ . Therefore the high-altitude $E$ -field is again enhanced in the direction of the “favorably curved” $B$ -field lines at low values of $\eta$ (and suppressed for “unfavorably curved” $B$ -field lines), coinciding at high $\eta$ for the different $\phi_{\rm PC}$ . It follows that $\dot{\gamma}_{\rm gain}$ , $\gamma_{\rm e}$ , and $\dot{\gamma}_{\rm loss}$ are higher for $\phi_{\rm PC}=0$ , and lower for $\phi_{\rm PC}=\pi$ ( $\epsilon\neq 0$ ). We note that $E_{\parallel,{\rm high}}$ is not so strongly dependent on $\eta$ as $E_{\parallel,{\rm low}}$ . There is no effect in changing $\phi_{\rm PC}$ in the case of $\alpha=0^{\circ}$ and $\epsilon=0$ , since the first term does not depend on $\phi_{\rm 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,{\rm high}}$ displays the same behavior at low $\eta$ as previously: for $\phi_{\rm PC}=0$ , $-E_{\parallel,{\rm high}}^{\epsilon\neq 0}>-E_{\parallel,{\rm high}}^{\epsilon=0}$ ; these are nearly equal for $\phi_{\rm PC}=\pi/2$ , and $-E_{\parallel,{\rm high}}^{\epsilon\neq 0}<-E_{\parallel,{\rm high}}^{\epsilon=0}$ for $\phi_{\rm PC}=\pi$ . For $\phi_{\rm PC}=\pi/2$ , the first term of $-E_{\parallel,{\rm high}}$ dominates the second, and for $\phi_{\rm PC}=\pi$ , the second term of $-E_{\parallel,{\rm high}}$ is always negative, but the positive first term dominates and therefore $-E_{\parallel,{\rm high}}$ does not change sign as $\eta$ increases. A similar behavior is also seen for $-E_{\parallel,{\rm low}}$ (boosted for non-zero $\epsilon$ and $\phi_{\rm PC}=0$ ). For $\alpha=45^{\circ}$ (second column of Figure 5), the second term $\propto\sin\alpha$ now contributes, stopping $-E_{\parallel,{\rm low}}$ from changing sign along $\eta$ for $\phi_{\rm PC}=0$ (vs. the case of $\alpha=0$ , first column, first row). The second term of $-E_{\parallel,{\rm low}}\sim x^{a}\cos\phi_{\rm PC}$ is comparable to the first at low $\eta$ , but quickly dominates as $\eta$ increases for $\phi_{\rm PC}=0$ . The second term of $-E_{\parallel,{\rm low}}$ remains positive so that we find $\eta_{\rm c}=5.1$ in this case (compare Figures 3 and 4). For $\phi_{\rm PC}=\pi/2$ we note that $-E_{\parallel,{\rm low}}^{\epsilon\neq 0}$ becomes negative with $\eta$ . For $\phi_{\rm PC}=\pi$ , the second term of $-E_{\parallel,{\rm low}}$ is negative, forcing this field to change sign; this change takes slightly longer to occur when $\epsilon\neq 0$ . The fact that $-E_{\parallel,{\rm high}}$ is positive leads to a “recovery” of the total $E$ -field, so that it becomes positive again at larger $\eta$ . 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}$ -field444As mentioned above, the sign of $\rho_{\rm GJ}$ changes near $\alpha\approx 90^{\circ}$ and $\phi_{\rm PC}\approx\pi$ , and we ignored such $B$ -field lines when calculating the emission from the pulsar. This is also the reason why our plots of $\gamma_{\rm e}$ vs. $s/R$ only go up to $\alpha=85^{\circ}$ , since we only consider electrons to be emitters of HE $\gamma$ -rays in our model. for $\alpha=85^{\circ}$ . For $\phi_{\rm PC}=0$ , the behaviors of $-E_{\parallel}$ , $\gamma_{\rm e}$ , and $\dot{\gamma}$ are very similar to the first and second columns in Figure 5. For $\phi_{\rm PC}=\pi/2$ (i.e., $a=0$ , as in the case for $\epsilon=0$ ) the first and second terms are positive for all $\eta$ and comparable in magnitude for $-E_{\parallel,{\rm low}}$ so that the sign change in the case of $\epsilon\neq 0$ may be ascribed to the structure of the offset-PC dipole $B$ -field ( $\phi_{\rm PC}$ increases along the field line since $B_{\phi}\neq 0$ so that $\cos\phi_{\rm PC}$ becomes negative). For $\phi_{\rm PC}=\pi$ , $-E_{\parallel,{\rm low}}^{\epsilon\neq 0}$ is not smaller than $-E_{\parallel,{\rm low}}^{\epsilon=0}$ at low $\eta$ . In this case the second term of $-E_{\parallel,{\rm low}}$ is always negative, and the sum of the two terms are very close for both $\epsilon\neq 0$ and $\epsilon=0$ so that $-E_{\parallel,{\rm low}}^{\epsilon\neq 0}\approx-E_{\parallel,{\rm low}}^{\epsilon=0}$ . Again, we see that $-E_{\parallel,{\rm low}}$ becomes negative, but the positive $-E_{\parallel,{\rm 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.7R_{\rm LC}}/R$ for $\phi_{\rm PC}=0$ and $\alpha=0^{\circ}$ , and beyond $\eta\approx{R_{\rm LC}}$ for $\phi_{\rm PC}=0$ and $\alpha=45^{\circ}$ . The notable exception occurs at large $\alpha$ 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_{\rm 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,{\rm 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 $>100$ 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).

4.2 Phase plots and Light Curves

We next perform simulations using a geometric modeling code (Section 2.3) which has the following free parameters: $\alpha$ , $\zeta$ , and $\epsilon$ (in case of the offset-PC dipole field). We fix the scaled co-latitude of the innermost ring of the gap ( $r_{\rm ovc}^{\rm min}=0.95$ ), PC rim ( $r_{\rm ovc}^{\rm max}=1.00$ ), and gap width $w=r_{\rm ovc}^{\rm max}-r_{\rm ovc}^{\rm min}=0.05$ . We choose $\alpha\in[0^{\circ},90^{\circ}]$ with a $1^{\circ}$ 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 $\alpha$ 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 $\zeta$ (the radiation pattern is a mirror image about $\zeta=90^{\circ}$ , also including a phase shift of $\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.5^{\circ},179.5^{\circ}]$ , both with a resolution of $2^{\circ}$ . Unique emission characteristics are visible in the light curves depending on the choice of $\alpha$ , $\zeta$ , $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 $\alpha$ the caustics extend over a larger range in $\zeta$ , 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 $\zeta$ below $45^{\circ}$ 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 $\epsilon$ results in qualitatively different phase plots and light curves, e.g., modulation at small $\alpha$ . Also, the caustics occupy a slightly larger region of phase space and seem more pronounced for larger $\epsilon$ and $\alpha$ values, accompanied by the same change in the position of the PCs as in Figure 6. The light curve shapes are also slightly different.

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 $\gamma$ -rays; see Sections 3.3). The caustic structure and resulting light curves are qualitatively different for various $\alpha$ compared to the constant $\epsilon_{\nu}$ case. The caustics appear smaller and less pronounced for larger $\alpha$ values (since $E_{\parallel}$ becomes lower as $\alpha$ increases), and extend over a smaller range in $\zeta$ . 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 $\alpha$ and $\zeta\approx(0^{\circ},180^{\circ})$ . This reflects the boosted $E_{\parallel}$ -field on the “favorably curved” $B$ -field lines (with $E_{\parallel}\propto{x^{a}}\cos\alpha$ , with $a=-\epsilon{\cos\phi_{\rm PC}}$ and $\phi_{\rm 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 $\alpha$ and $\zeta$ increase, more peaks become visible, with emission still visible from both poles as seen for larger $\alpha$ and $\zeta$ values, e.g., $\alpha=75^{\circ}$ and $\zeta=75^{\circ}$ .

If we compare Figure 6 with Figure 8 (also Figure 7 with Figure 9), we notice that when we take $E_{\parallel}$ 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 $>100$ 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 $\epsilon$ , 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 $\epsilon$ 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 $\epsilon$ . For $\epsilon=0$ the light curve has a single peak and as $\epsilon$ increases, the peak becomes slightly narrower. Also, for larger $\epsilon$ values, the caustic structure becomes slightly broader and more pronounced. Figures 10(c) and 10(d) represent the offset-PC dipole field and SG model with variable $\epsilon_{{\rm\nu}}$ . When we compare the phase plots of Figures 10(a) and 10(c), the caustics are dimmer and smaller due to the low SG $E$ -field, and the light curves display less off-peak emission. As $\epsilon$ increases, some small features become more pronounced.

4.3 $\chi^{2}(\alpha,\zeta)$ Contours and Best-fit Light Curves

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 ${\rm log_{10}\xi^{2}}$ (color bar) as well as the optimal ( $\alpha$ , $\zeta$ ) combination. We determine errors on $\alpha$ and $\zeta$ 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 $\alpha$ and $\zeta$ which surrounds the $3\sigma$ contour. We illustrate this in Figure 11(a), with the white lines indicating the bounding box $[\alpha_{{\rm min}},\alpha_{{\rm max}}]$ and $[\zeta_{{\rm min}},\zeta_{{\rm max}}]$ (see enlargement in bottom left corner of panel (a)). We choose errors of $1^{\circ}$ for cases when the errors were smaller than $1^{\circ}$ (given our chosen resolution of $1^{\circ}$ ). 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 $\Delta\phi_{\rm L}=0$ by a white star. The curved region ranging from $\zeta=70^{\circ}$ downward to $\zeta=10^{\circ}$ , over the entire range of $\alpha$ , is caused by the caustic structure as seen in the phase plots (i.e., there is no emission visible for low values of $\alpha$ and $\zeta$ – the turquoise bottom-left region).

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 ${\rm 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 $\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 $\alpha$ and $\zeta$ 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 ${\rm 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 $\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.

5 FURTHER INVESTIGATION

5.1 Light Curves in a Different Waveband

Since the SG $E$ -field (see Section 3.3) is low, CRR is reached in most cases but only at high $\eta$ and small $\alpha$ . 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_{\rm 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)

[TABLE]

with $\rho_{\rm curv,8}\sim\rho_{\rm curv}/{10^{8}}$ cm the curvature radius of the $B$ -field line and $E_{\parallel,4}\sim E_{\parallel}/{10^{4}}$ statvolt cm*-1*, the $E$ -field parallel to the $B$ -field. For our given SG $E$ -field with a magnitude of $E_{\parallel}\sim 10^{2}$ statvolt cm*-1* the estimated cutoff is $E_{\rm CR}\sim 90$ MeV. This leads to pulsar emission being emitted in the hard X-ray waveband, and cannot be compared via $\chi^{2}$ to Fermi ( $>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_{\rm min}>1$ 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 $\alpha$ , e.g., at $\alpha=90^{\circ}$ the visibility is enhanced. The peaks are also wider at low $\alpha$ . Sometimes extra emission features appear, leading to small changes in the light curve shapes.

5.2 Effect of Increasing the $E$ -field

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_{\rm CR}\sim 4$ GeV which is in the energy range of Fermi ( $>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 $\alpha$ . 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 $\alpha$ and $\zeta$ values. This leads to small changes in the light curve shapes. At smaller $\alpha$ values, the emission around the PC forms a circular pattern that becomes smaller as $\alpha$ 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 $\alpha$ the background becomes feature-rich, but not at significant intensities, however.

When we boost the low $E$ -field, we find that the CRR limit is in fact reached almost immediately at lower $\eta$ for certain parameter combinations of $\alpha$ , and $\phi_{\rm 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^{+3}_{-1}$ , $\zeta=51^{+2}_{-5}$ , $A=1.1$ , and $\Delta\phi_{\rm L}=0.55$ . In Figure 15 we show our significance contour ${\rm 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).

6 COMPARISON OF BEST-FIT PARAMETERS FOR DIFFERENT MODELS

We next follow the same approach as Pierbattista et al. (2015) to compare the various optimal solutions of the different models, in two ways: (i) per $B$ -field and model combination and (ii) overall (for all $B$ -field and model combinations). We determine the difference between the scaled555We therefore first scale the $\chi^{2}$ values using the optimal value obtained for a particular $B$ -field, and second we scale these using the overall optimal value irrespective of $B$ -field. $\chi^{2}$ of the optimal model, $\xi_{\rm opt}^{2}$ , and the other models ( $\xi^{2}$ ) using Equation (15), substituting $N_{\rm dof}=96$ , as summarized in Table The effect of an offset polar cap dipolar magnetic field on the modeling of the Vela pulsar’s $\gamma$ -ray light curves. The best-fit parameters for each $B$ -field and geometric model combination, including the case for $100E_{\parallel}$ , are summarized in Table The effect of an offset polar cap dipolar magnetic field on the modeling of the Vela pulsar’s $\gamma$ -ray light curves. The table includes the different model combinations, the optimal unscaled $\chi^{2}$ value for each combination, the best-fit free parameters with $3\sigma$ errors on $\alpha$ and $\zeta$ , and the comparison between models per $B$ -field ( $\Delta\xi^{2}_{\rm B}$ ) and overall ( $\Delta\xi^{2}_{\rm all}$ , with $\Delta\xi^{2}=0$ representing the best-fit solution for each $B$ -field or the overall optimal fit; Pierbattista et al., 2015). We also include several multi-wavelength independent fits (all for the Vela pulsar).

In Figure 16 we label the different $B$ -field structures assumed in the various models as well as the overall comparison along the $x$ -axis, and plot $\Delta\xi^{2}_{\rm B}$ and $\Delta\xi^{2}_{\rm all}$ on the $y$ -axis. We represent the TPC geometry with a circle, the OG with a square, and for the offset-PC dipole field we represent the various $\epsilon$ values for constant $\epsilon_{\nu}$ by different colored stars, for variable $\epsilon_{\nu}$ by different colored left pointing triangles, and for the case of $100E_{\parallel}$ by different colored upright triangles, as indicated in the legend. The dashed horizontal lines indicate the confidence levels we obtained by calculating the expected $\Delta\xi^{2}$ values using an online $\chi^{2}$ statistical calculator666http://easycalculation.com/statistics/critical-value-for-chi-square.php for $N_{\rm dof}=96$ degrees of freedom,777We note that Pierbattista et al. (2015) assumed that $\Delta\xi^{2}$ follows a $\chi^{2}$ distribution with $N_{\rm dof}$ degrees of freedom. We will follow this approximation here, assuming that the best-fit model provides a good fit to the observed light curves. The degrees of freedom may in reality slightly differ, however, and the matter is complicated by the fact that we want to statically compare non-nested models. A Monte Carlo approach would be preferable to find these significance levels. However, our main conclusions will not change for slight changes in these levels (which may be different for each $B$ -field and model combination), and so we do not pursue this matter any further. i.e., using $p$ -values of $p_{1\sigma}=1-0.682$ , $p_{2\sigma}=1-0.954$ , and $p_{3\sigma}=1-0.9973$ . We found critical values of $\Delta\xi^{2}=102.06$ ( $1\sigma$ ), $120.60$ ( $2\sigma$ ), and $139.05$ ( $3\sigma$ ) respectively. These confidence levels are used as indicators of when to reject or accept an alternative fit compared to the optimum fit. The last column represents fits for all models, irrespective of $B$ -field.

For the static dipole field the TPC model gives the optimum fit and the OG model lies within $1\sigma$ , implying that the OG geometry may provide an acceptable alternative fit to the data in this case. For the RVD field the TPC model is significantly rejected beyond the $3\sigma$ level (not shown on plot), and the OG model is preferred. We show three cases for the offset-PC dipole field, including the TPC model assuming constant $\epsilon_{\nu}$ , the SG model assuming variable $\epsilon_{\nu}$ , and the latter with an $E_{\parallel}$ -field increased by a factor of $100$ . The optimal fits for the offset-PC dipole field and TPC model reveal that a smaller offset ( $\epsilon$ ) is generally preferred for constant $\epsilon_{\nu}$ , while a larger offset is preferred for variable $\epsilon_{\nu}$ (but not significantly), with all alternative fits falling within $1\sigma$ . However, when we increase $E_{\parallel}$ , a smaller offset is again preferred for the SG and variable $\epsilon_{\nu}$ case. When we compare all model and $B$ -field combinations with the overall best fit (i.e., rescaling the $\chi^{2}$ values of all combinations using the optimal fit involving the RVD $B$ -field and OG model), we notice that the static dipole and TPC model falls within $2\sigma$ , whereas the static OG model lies within $3\sigma$ . We also note that the usual offset-PC dipole $B$ -field and TPC model combination (for all $\epsilon$ values) is above $1\sigma$ (with some fits $<2\sigma$ ), but the offset-PC dipole $B$ -field and SG model combination (for all $\epsilon$ values) is significantly rejected ( $>3\sigma$ ). However, the case of the offset-PC dipole field and a higher SG $E_{\parallel}$ for all $\epsilon$ values leads to a recovery, since all the fits fall within $2\sigma$ and delivers an overall optimal fit for $\epsilon=0$ , second only to the RVD and OG model fit.

Several multi-wavelength studies have been performed for Vela, using the radio, X-ray, and $\gamma$ -ray data, in order to find constraints on $\alpha$ and $\zeta$ . We only fit the $\gamma$ -ray light curve, because we did not want to bias our results by using a geometric radio emission model (DeCesar, 2013). However, Johnston et al. (2005) determined the radio polarization position angle from polarization data by applying a rotating vector model (RVM) fit to the datafinding best-fit values of $\alpha=53^{\circ}$ and $\zeta=59.5^{\circ}$ , with an impact angle of $\beta=\zeta-\alpha=6.5^{\circ}$ . Ng & Romani (2008) applied a torus-fitting technique (Ng & Romani, 2004) to fit the Chandra data in order to constrain the Vela X pulsar wind nebula (PWN) geometry, deriving a value of $\zeta=63.6{{}_{-0.05}^{+0.07}}^{\circ}$ represented by the dashed black line in Figure 17. Watters et al. (2009) modeled light curves using the PC, TPC, and OG geometries in conjunction with an RVD field, thereby constraining the geometrical parameters $\alpha$ , $\zeta$ , and also finding small $\beta$ in the case of the PC model. They found a good fit for their TPC model at $\alpha=62^{\circ}-68^{\circ}$ and $\zeta=64^{\circ}$ , and for the OG geometry at $\alpha=75{}^{\circ}$ and $\zeta=64^{\circ}$ . We find that our best-fit values for the RVD field, for both the TPC and OG models, are in good agreement with those found by Watters et al. (2009). DeCesar (2013) followed a similar approach to ours, but for the RVD and FF $B$ -fields combined with emission geometries such as the SG (symmetric and asymmetric cases) and OG. They have different free model parameters including $\alpha$ , $\zeta$ , $w$ (gap width), and $R_{\rm max}$ (maximum emission radius), and determined errors on their best fits using the $3\sigma$ confidence intervals. They found best-fit solutions for the RVD and OG model at $\alpha=88{{}_{-3}^{+2}}^{\circ}$ and $\zeta=66.5{{}_{-1}^{+1}}^{\circ}$ , which is within $10^{\circ}$ or less compared to our best-fit solution. Their overall best fit was for the FF $B$ -field and OG geometry, with $\alpha={80^{+1}_{-1}}^{\circ}$ and $\zeta={53^{+1}_{-1}}^{\circ}$ . Pierbattista et al. (2015) found a best-fit solution for Vela using the RVD field and OG model combination at $\alpha={71_{-2}^{+2}}^{\circ}$ and $\zeta=83{{}_{-2}^{+2}}^{\circ}$ , with $\zeta$ exceeding the best-fit solution we found by nearly $15^{\circ}$ . However, they fit both the $\gamma$ -ray and radio light curves, which may explain this discrepancy. We summarize all these multi-wavelength fits and more in Table The effect of an offset polar cap dipolar magnetic field on the modeling of the Vela pulsar’s $\gamma$ -ray light curves.

We graphically summarize the best-fit $\alpha$ and $\zeta$ , with errors, from this and other works in Figure 17. We notice that the best fits generally prefer a large $\alpha$ or $\zeta$ or both. It is encouraging that many of the best-fit solutions lie near the $\zeta$ inferred from the PWN torus fitting (Ng & Romani, 2008), notably for the RVD $B$ -field. A significant fraction of fits furthermore lie near the $\alpha-\zeta$ diagonal, i.e., they prefer a small impact angle, probably due to radio visibility constraints (Johnson et al., 2014). For an isotropic distribution of pulsar viewing angles, one expects $\zeta$ values to be distributed as $\sin(\zeta)$ between $\zeta=[0^{\circ},90^{\circ}]$ , i.e., large $\zeta$ values are much more likely than small $\zeta$ values, which seems to correspond to the large best-fit $\zeta$ values we obtain. There seems to be a reasonable correspondence between our results obtained for geometric models and those of other authors, but less so for the offset-PC dipole $B$ -field, and in particular for the SG $E$ -field case. The lone fit near $(20^{\circ},70^{\circ})$ may be explained by the fact that a very similar fit, but one with slightly worse $\chi^{2}$ , is found at $(50^{\circ},80^{\circ})$ . If we discard the non-optimal TPC / SG fits, we see that the optimal fits will cluster near the other fits at large $\alpha$ and $\zeta$ . Although our best fits for the offset-PC dipole $B$ -field are clustered, it seems that increasing $\epsilon$ leads to a marginal decrease in $\zeta$ for the TPC model (light green) and opposite for SG (dark green), but not significantly (see Table The effect of an offset polar cap dipolar magnetic field on the modeling of the Vela pulsar’s $\gamma$ -ray light curves). For our increased SG $E$ -field case (brown) we note that the fits now cluster inside the gray area above the fits for the static dipole and TPC, and offset-PC dipole for both the TPC and SG geometries.

7 CONCLUSIONS

We investigated the impact of different magnetospheric structures on predicted $\gamma$ -ray pulsar light curve characteristics. We extended our code, which already included the static dipole and RVD $B$ -fields, by implementing an additional $B$ -field, i.e., the symmetric offset-PC dipole field (Harding & Muslimov, 2011a, b) characterized by an offset $\epsilon$ of the magnetic PCs. We also included the full accelerating SG $E$ -field corrected for GR effects up to high altitudes. For the offset-PC dipole field we only considered the TPC (assuming uniform $\epsilon_{\nu}$ ) and SG (modulating the $\epsilon_{\nu}$ using the $E$ -field) models, since we do not have $E$ -field expressions available for the OG model for this particular $B$ -field. We matched the low-altitude and high-altitude solutions of the SG $E_{\parallel}$ by determining the matching parameter $\eta_{\rm c}(P,\dot{P},\alpha,\epsilon,\xi,\phi_{\rm PC})$ on each field line in multivariate space. Once we calculated the general $E$ -field we could solve the particle transport equation. This yielded the particle energy $\gamma_{\rm e}(\eta)$ , necessary for determining the CR $\epsilon_{\nu}$ and to test whether the CRR limit is attained. For the case of a variable $\epsilon_{\nu}$ , we found that the CRR limit is reached for many parameter combinations (of $\alpha,\epsilon$ and $\phi_{\rm PC}$ ; see Figure 5), albeit only at large $\eta$ . A notable exception occurred at large $\alpha$ where the first term of each $E$ -field expression (e.g., Equations (7) and (8)) became lower and the second term played a larger role, leading to smaller gain rates and therefore smaller Lorentz factors $\gamma_{\rm e}$ .

We concluded that the magnetospheric structure and emission geometry have an important effect on the predicted $\gamma$ -ray pulsar light curves. However, the presence of an $E$ -field may have an even greater effect than small changes in the $B$ -field and emission geometries: When we included an SG $E$ -field, thereby modulating $\epsilon_{\nu}$ , the resulting phase plots and light curves became qualitatively different compared to the geometric case.

We fit our model light curves to the observed Fermi-measured Vela light curve for each $B$ -field and geometric model combination. We found that the RVD field and OG model combination fit the observed light curve the best for $(\alpha,\zeta,A,\Delta\phi_{\rm L})=({78_{-1}^{+1}}^{\circ},{69_{-1}^{+2}}^{\circ},1.3,0.00)$ and an unscaled $\chi^{2}=3.84\times{10}^{4}$ . As seen in Figure 16, for the RVD field an OG model is significantly preferred over the TPC model, given the characteristically low off-peak emission. For the other field and model combinations there was no significantly preferred model (per $B$ -field), since all the alternative models may provide an acceptable alternative fit to the data, within $1\sigma$ . The offset-PC dipole field for constant $\epsilon_{\nu}$ favored smaller values of $\epsilon$ , and for variable $\epsilon_{\nu}$ larger $\epsilon$ values, but not significantly so ( $<1\sigma$ ). When comparing all cases (i.e., all $B$ -fields), we noted that the offset-PC dipole field for variable $\epsilon_{\nu}$ was significantly rejected ( $>3\sigma$ ).

We further investigated the effect which the SG $E_{\parallel}$ had on our predicted light curves in two ways. First, we lowered the minimum photon energy from $E_{\rm min}=100$ MeV to $E_{\rm min}=1$ MeV, leading to emission in the hard X-ray waveband. We noted new caustic structures and emission features on the resulting phase plots and light curves that were absent when $E_{\rm min}>100$ MeV. Since we wanted to compare our model light curves to Fermi data we increased the usual low SG $E$ -field by a factor of 100 (with a spectral cutoff $E_{\rm CR}\sim 4$ GeV). When solving the particle transport equation, we noticed that the CRR limit is now reached in most cases at lower $\eta$ . The increased $E$ -field also had a great impact on the phase plots, e.g., extended caustic structures and new emission features as well as different light curve shapes emerged. We also compared the best-fit light curves for the offset-PC dipole $B$ -field and $100E_{\parallel}$ combination for each $\epsilon$ (Figure 16) and noted that a smaller $\epsilon$ was again preferred (although not significantly; $<1\sigma$ ). However, when we compared this case to the other $B$ -field and model combinations, we found statistically better $\chi^{2}$ fits for all $\epsilon$ values with an optimal fit at $\alpha={75^{+3}_{-1}}^{\circ}$ , $\zeta={51^{+2}_{-5}}^{\circ}$ , $A=1.1$ , and $\Delta\phi_{\rm L}=0.55$ for $\epsilon=0$ , being second in quality only to the RVD and OG model fit.

We graphically compared the best-fit $\alpha$ and $\zeta$ , with errors, from this and other works in Figure 17. We noted that many of the best-fit solutions cluster inside the gray area at larger $\alpha$ and $\zeta$ . Some fits lie near the $\alpha-\zeta$ diagonal (possibly due to radio visibility constraints in some cases) as well as near the $\zeta$ inferred from the PWN torus fitting (Ng & Romani, 2008), notably for the RVD $B$ -field. There was reasonable correspondence between our results obtained for geometric models and those of other independent studies. When we discarded the non-optimal TPC / SG fits, we saw that the optimal fits clustered near the other fits at large $\alpha$ and $\zeta$ . For our increased SG $E$ -field and offset-PC dipole combination (brown) we noted that these fits clustered at larger $\alpha$ and $\zeta$ .

There have been several indications that the SG $E$ -field may be larger than initially thought. For example, (i) population synthesis studies found that the SG $\gamma$ -ray luminosity may be too low, pointing to an increased $E$ -field and / or particle current through the gap (e.g., Pierbattista et al., 2015). Furthermore, if the $E$ -field is too low, one is not able to reproduce the (ii) observed spectral cutoffs of a few GeV (Section 5.1; Abdo et al., 2013). We found additional indications for an enhanced SG $E$ -field. A larger $E$ -field (increased by a factor of 100) led to (iii) statistically improved $\chi^{2}$ fits with respect to the light curves. Moreover, the inferred best-fit $\alpha$ and $\zeta$ parameters for this $E$ -field (iv) clustered near the best fits of independent studies. We additionally observed that a larger SG $E$ -field also (v) increased the particle energy gain rates and therefore yielded a larger particle energy $\gamma_{\rm e}$ (giving CR that is visible in the Fermi band) as well as leading to a CRR regime already close to the stellar surface. These evidences may point to a reconsideration of the boundary conditions assumed by Muslimov & Harding (2004) which suppressed the $E_{\parallel}$ at high altitudes. They assumed equipotentiality of the SG boundaries as well as the steady state drift of charged particles across the SG $B$ -field lines, implying $E_{\perp}\approx 0$ at high-altitudes, with the flux of charges remaining constant up to high altitudes. One possible way to bring self-consistency may be implementation of the newly developed FIDO model that includes global magnetospheric properties and calculates the $B$ -field and $E$ -field self-consistently.

We envision several future projects that may emanate from this study. One could continue to extend the range of $\epsilon$ for which our code finds the PC rim, since more complex field solutions, e.g., the dissipative and FF field structures, may be associated with larger PC offsets. However, the offset-PC dipole solutions have limited applicability to outer magnetosphere emission since they use the static dipole frame and do not model the field line sweep back. Therefore, it would be preferable to investigate the $B$ -fields and $E$ -fields of the dissipative models and solve the transport equation to test if the particles reach the CRR limit. The effect of these new fields on the phase plots and light curves can also be studied. There is also potential for multi-wavelength studies, such as light curve modeling in the other energy bands, e.g., combining radio and $\gamma$ -ray light curves (see Seyffert et al., 2010, 2012; Pierbattista et al., 2015). One could furthermore model energy-dependent light curves, such as those available for Vela and other bright pulsars using Fermi data (e.g., Abdo et al., 2009). Lastly, model phase-resolved spectra can be constructed which is an important test of the $E_{\parallel}$ -field magnitude and spatial dependence.

We thank Marco Pierbattista, Tyrel Johnson, Lucas Guillemot, and Bertie Seyffert for fruitful discussions. This work is based on the research supported wholly / in part by the National Research Foundation of South Africa (NRF; Grant Numbers 87613, 90822, 92860, 93278, and 99072). The Grantholder acknowledges that opinions, findings and conclusions or recommendations expressed in any publication generated by the NRF supported research is that of the author(s), and that the NRF accepts no liability whatsoever in this regard. A.K.H. acknowledges the support from the NASA Astrophysics Theory Program. C.V. and A.K.H. acknowledge support from the Fermi Guest Investigator Program.

APPENDIX

TRANSFORMATION OF THE $B$ -FIELD FROM THE MAGNETIC TO THE ROTATIONAL FRAME

Consider a general $B$ -field specified in the magnetic frame (indicated by primed coordinates), in terms of spherical coordinates

[TABLE]

The spherical unit vectors may be expressed in terms of Cartesian unit vectors (e.g., Griffiths, 1995):

[TABLE]

To transform the $B$ -vector components from a spherical base unit vector set to a Cartesian one, substitute Equation (19) into (18) and rearrange according to unit basis vectors and compare with

[TABLE]

This yields

[TABLE]

To transform the spherical coordinates in Equation (21) to Cartesian coordinates in the magnetic frame, use

[TABLE]

We then obtain

[TABLE]

Lastly, we transform the $B$ -field components and coordinates from the magnetic frame to the rotational frame using a rotation of axes. We rotate the Cartesian frame through an angle $-\alpha$ (the angle between the $\bm{\Omega}$ and $\bm{\mu}$ axes), thereby transforming the $B$ -field from the magnetic to the rotational frame (indicated by the unprimed coordinates), letting ${\rm{\mathbf{y}}}\parallel{\rm{\mathbf{y}}}^{\prime}$ :

[TABLE]

This yields

[TABLE]

Lastly, we transform the the Cartesian coordinates of the position vector from the magnetic to the rotational frame using a similar rotation matrix as Equation (24), and substitute it into Equation (25) to obtain

[TABLE]

Bibliography65

The reference list from the paper itself. Each links out to its DOI / PubMed record.

1Abdo et al. (2009) Abdo, A. A., Ackermann, M., Atwood, W. B., et al. 2009, Ap J, 696, 1084
2Abdo et al. (2013) Abdo, A. A., Ajello, M., Allafort, A., et al. 2013, Ap JS, 208, 17
3Aleksić et al. (2011) Aleksić, J., Alvarez, E. A., Antonelli, L. A., et al. 2011, Ap J, 742, 43
4Aleksić et al. (2012) Aleksić, J., Alvarez, E. A., Antonelli, L. A., et al. 2012, A&A, 540, A 69
5Aliu et al. (2011) Aliu, E., Arlen, T., Aune, T., et al. 2011, Science, 334, 69
6Arendt & Eilek (1998) Arendt, Jr., P. N., & Eilek, J. A. 1998, ar Xiv:astro-ph/9801257
7Arons (1983) Arons, J. 1983, Ap J, 266, 215
8Atwood et al. (2009) Atwood, W. B., Abdo, A. A., Ackermann, M., et al. 2009, Ap J, 697, 1071

TL;DR

Contribution

Findings

Abstract

Peer Reviews

Videos

The effect of an offset polar cap dipolar magnetic field on the modeling of the Vela pulsar’s γ\gammaγ-ray light curves

Abstract

1 INTRODUCTION

2 GEOMETRIC MODEL DESCRIPTION

2.1 TPC and OG Geometries

2.2 BBB-field Structures

2.3 Geometric Modeling Code

3 IMPLEMENTATION OF AN OFFSET-PC DIPOLE BBB-FIELD AND ASSOCIATED SG EEE-FIELD AND A SEARCH FOR THE BEST-FIT LIGHT CURVE

3.1 Transformation of a BBB-field from the Magnetic to the Rotational Frame

3.2 Finding the PC Rim and Extending the Range of ϵ\epsilonϵ

3.3 Incorporating an SG EEE-field

3.4 Determining the Matching Parameter ηc\eta_{\rm c}ηc​

3.5 Chi-squared Fitting Method

4 RESULTS

4.1 Solving the Particle Transport Equation

4.2 Phase plots and Light Curves

4.3 χ2(α,ζ)\chi^{2}(\alpha,\zeta)χ2(α,ζ) Contours and Best-fit Light Curves

5 FURTHER INVESTIGATION

5.1 Light Curves in a Different Waveband

5.2 Effect of Increasing the EEE-field

6 COMPARISON OF BEST-FIT PARAMETERS FOR DIFFERENT MODELS

7 CONCLUSIONS

APPENDIX

The effect of an offset polar cap dipolar magnetic field on the modeling of the Vela pulsar’s $\gamma$ -ray light curves

2.2 $B$ -field Structures

3 IMPLEMENTATION OF AN OFFSET-PC DIPOLE $B$ -FIELD AND ASSOCIATED SG $E$ -FIELD AND A SEARCH FOR THE BEST-FIT LIGHT CURVE

3.1 Transformation of a $B$ -field from the Magnetic to the Rotational Frame

3.2 Finding the PC Rim and Extending the Range of $\epsilon$

3.3 Incorporating an SG $E$ -field

3.4 Determining the Matching Parameter $\eta_{\rm c}$

4.3 $\chi^{2}(\alpha,\zeta)$ Contours and Best-fit Light Curves

5.2 Effect of Increasing the $E$ -field