Inclined Pulsar Magnetospheres in General Relativity: Polar Caps for the Dipole, Quadrudipole and Beyond

TL;DR

This paper develops a general analytical and numerical method to determine the shape and charge distribution of pulsar polar caps, incorporating general relativity effects for various magnetic field configurations, including dipole and quadrudipole.

Contribution

It introduces a novel combined analytical-numerical approach to model polar caps in relativistic pulsar magnetospheres with arbitrary magnetic fields.

Findings

Inclined dipole polar caps are significantly affected by general relativity.

Quadrudipole fields can produce thin annular polar caps.

The method enables realistic modeling of pulsar emission regions.

Abstract

In the canonical model of a pulsar, rotational energy is transmitted through the surrounding plasma via two electrical circuits, each connecting to the star over a small region known as a "polar cap." For a dipole-magnetized star, the polar caps coincide with the magnetic poles (hence the name), but in general, they can occur at any place and take any shape. In light of their crucial importance to most models of pulsar emission (from radio to X-ray to wind), we develop a general technique for determining polar cap properties. We consider a perfectly conducting star surrounded by a force-free magnetosphere and include the effects of general relativity. Using a combined numerical-analytical technique that leverages the rotation rate as a small parameter, we derive a general analytic formula for the polar cap shape and charge-current distribution as a function of the stellar mass, radius, rotation rate, moment of inertia, and magnetic field. We present results for dipole and quadrudipole fields (superposed dipole and quadrupole) inclined relative to the axis of rotation. The inclined dipole polar cap results are the first to include general relativity, and they confirm its essential role in the pulsar problem. The quadrudipole pulsar illustrates the phenomenon of thin annular polar caps. More generally, our method lays a foundation for detailed modeling of pulsar emission with realistic magnetic fields.

Equations

$${\mathcal{L}}_{\xi }g=0,\xi ={\partial }_t+\Omega {\partial }_{\phi },$$

$$u^{\mu }={\xi }^{\mu }/\sqrt{{\rm Y}},{\rm Y}=-g_{\mu \nu }{\xi }^{\mu }{\xi }^{\nu }.$$

$${\mathcal{L}}_{\xi }F=0.$$

$$F=\mathrm{d}{\varphi }_1\wedge \mathrm{d}{\varphi }_2.$$

$${\varphi }_1$$

$${\varphi }_2$$

$$F=\mathrm{d}{\psi }_1\wedge \mathrm{d}{\psi }_2,{\psi }_i={\psi }_i(r,\theta ,\phi -\Omega t).$$

$$\xi \cdot F=0,$$

$$\Lambda ({\psi }_1,{\psi }_2)=\sqrt{{\rm Y}}{J}_{\parallel },{\mathbf{J}}_c=J_{\parallel }{\mathbf{B}}_c.$$

$$d{s}^2$$

$$\displaystyle\quad\,+r^{2}\mathopen{}\mathclose{{}\left[\mathop{}\!\mathrm{d}\theta^{2}+\sin^{2}{\theta}\mathopen{}\mathclose{{}\left(\mathop{}\!\mathrm{d}\varphi-\Omega_{Z}\mathop{}\!\mathrm{d}t}\right)^{2}}\right],\quad r>R_{\star},$$

$${\Omega }_Z=\frac{2\hat{I}}{r^3}\Omega .$$

$$\displaystyle\Upsilon=\frac{r-2M}{r}-\mathopen{}\mathclose{{}\left[\mathopen{}\mathclose{{}\left(\Omega-\Omega_{Z}}\right)r\sin{\theta}}\right]^{2}.$$

$$R_L=\frac{1}{\Omega },$$

$$ϵ=\frac{R_{\star }}{R_L},\mathcal{C}=\frac{2M}{R_{\star }},\mathcal{I}=\frac{\hat{I}}{M{R}_{\star }^2},$$

$$F=\mathrm{d}{\psi }_1\wedge \mathrm{d}{\psi }_2,{\psi }_i={\psi }_i(r,\theta ,\phi -\Omega t;ϵ).$$

$$ϵ\to 0\text{at fixed }{R}_{\star }$$

$$ϵ\to 0\text{at fixed }{R}_L,$$

$$R_{\star }$$

$$R_L$$

$$M$$

$$\hat{I}$$

$$\Omega$$

$$r\ll R_L,t\ll 2\pi /\Omega ,$$

$$\displaystyle\quad r\gg R_{\star}\quad\mathopen{}\mathclose{{}\left(\textrm{and }\gg M,\hat{I}^{1/3}}\right).$$

$$\text{overlap region:}{R}_{\star }\ll r\ll R_L,t\ll 2\pi /\Omega ,$$

$$d{s}^2$$

$$\displaystyle\quad\,+r^{2}\mathopen{}\mathclose{{}\left[\mathop{}\!\mathrm{d}\theta^{2}+\sin^{2}{\theta}\mathop{}\!\mathrm{d}\varphi^{2}}\right]-2\Omega_{Z}\mathop{}\!\mathrm{d}t\mathop{}\!\mathrm{d}\varphi+\mathcal{O}_{\star}\mathopen{}\mathclose{{}\left(\epsilon^{2}}\right),$$

$${\psi }_1$$

$${\psi }_2$$

$$F$$

$$\displaystyle\quad\,-\Omega t\mathopen{}\mathclose{{}\left[\mathop{}\!\mathrm{d}\alpha\wedge\mathop{}\!\mathrm{d}\mathopen{}\mathclose{{}\left(\mathop{}\!\partial_{\phi}\beta}\right)+\mathop{}\!\mathrm{d}\mathopen{}\mathclose{{}\left(\mathop{}\!\partial_{\phi}\alpha}\right)\wedge\mathop{}\!\mathrm{d}\beta}\right]$$

$$\displaystyle\quad\,+\epsilon\mathopen{}\mathclose{{}\left(\mathop{}\!\mathrm{d}\alpha\wedge\mathop{}\!\mathrm{d}\tilde{\beta}+\mathop{}\!\mathrm{d}\tilde{\alpha}\wedge\mathop{}\!\mathrm{d}\beta}\right)+\mathcal{O}_{\star}\mathopen{}\mathclose{{}\left(\epsilon^{2}}\right).$$

$$F^{(0)}=\mathrm{d}\alpha \wedge \mathrm{d}\beta .$$

$$J={\mathcal{O}}_{\star }(ϵ).$$

$$\alpha$$

Full text

I Introduction

Fifty years after the basic elements of pulsar theory were established Goldreich and Julian (1969), self-consistent modeling remains challenging. To simplify the problem, the community has mainly focused on a canonical choice of magnetic field: the pure dipole. This endeavor has come to fruition in the last decade, as the dipole pulsar has been self-consistently modeled using increasingly detailed descriptions of the surrounding plasma, from force-free electrodynamics Contopoulos et al. (1999); Gruzinov (2005); McKinney (2006); Timokhin (2006); Spitkovsky (2006); Kalapotharakos and Contopoulos (2009); Li et al. (2012); Kalapotharakos et al. (2012); Pétri (2012); Lehner et al. (2012); Palenzuela (2013); Ruiz et al. (2014); Pétri (2016) to magnetohydrodynamics Komissarov (2006); Tchekhovskoy et al. (2013) to kinetic theory Philippov et al. (2014); Chen and Beloborodov (2014); Belyaev (2015); Cerutti et al. (2015); Philippov et al. (2015a, b).

In contrast, alternative magnetic field configurations have remained relatively unexplored. This is due mainly to the high cost of numerical simulations and the lack of an obvious alternative field configuration to choose. However, there is little reason to believe in a pure dipole field, and the wide variation in emission properties among the known pulsars seems to be naturally accounted for by a correspondingly wide variation in the stellar magnetic field. It is therefore of interest to develop efficient techniques for exploring the effects of more realistic magnetic fields.

Building on our recent work in the axisymmetric case Gralla et al. (2016a) (hereafter Paper I; see also Ref. Belyaev and Parfrey (2016)), in this paper, we introduce a general method for determining the near-field charge and current flow (i.e., the pulsar polar caps) for a given magnetic field geometry on a given general relativistic star. The key observation underpinning our analysis is that the pulsar problem contains a small parameter: the ratio of the stellar radius $R_{\star}$ to the light cylinder radius $R_{L}=c/\Omega$ (where $\Omega$ denotes the rotation rate and $c$ is the speed of light). This parameter $\epsilon=R_{\star}/R_{L}$ is proportional to the surface rotation velocity and ranges from $10^{-4}$ to $10^{-1}$ for rotation-powered pulsars.

The small value of $\epsilon$ makes numerical work challenging, since both scales $R_{\star}$ and $R_{L}$ must be resolved. Indeed, most numerical simulations are run at large values $\epsilon\sim 1/5$ to reduce the dynamic range. In contrast, analytic methods can shine at small values $\epsilon\ll 1$. To resolve both scales, one can use the method of matched asymptotic expansions Gralla et al. (2016a). This involves making separate far $(r\gg R_{\star})$ and near $(r\ll R_{L})$ expansions, and matching in the overlap region $R_{\star}\ll r\ll R_{L}$ of shared validity.

We apply this method to a force-free magnetosphere surrounding a perfectly conducting general relativistic star Hartle and Thorne (1968). In the far region, the equations reduce to the force-free magnetosphere of a rotating point dipole in flat spacetime, which we solve numerically. In the near region, the equations become those of a static vacuum magnetic field in the Schwarzschild spacetime (representing the intrinsic magnetic field of the star), whose general solution is known in closed form. Our near-far matching is expedited by a conserved quantity on magnetic field lines that generalizes the familiar field-aligned current (see App. A).

The main quantities of interest are the leading-order charge and current near the star. The charge is essentially induced by rotation, and is easily computed from the stellar magnetic field and rotation rate. For the current, we fit for the conserved quantity in the far zone and then “paint” the associated current onto each field line, following it down to the surface of the star. This last step amounts to finding Euler potentials for the stellar magnetic field, and hence is highly non-trivial in general. It is, however, trivial to find Euler potentials for fields that are axisymmetric about some axis (not necessarily aligned with the rotation axis). We thereby provide definite analytical formulae for the polar cap structure of an inclined pulsar with an arbitrary axisymmetric magnetic field. For intrinsically nonaxisymmetric fields, the Euler potentials must in general be found numerically. However, this task is still considerably simpler than running a complete self-consistent simulation.

While the magnetization is the most important, all of the stellar parameters affect the polar cap shapes and properties. In particular, the compactness and moment of inertia control the redshift and frame-drag effects (respectively) of general relativity, and may be used to quantify the importance of these effects. We emphasize, however, that in our approach, general relativistic effects are not added in piecemeal or “by hand”; we simply begin with force-free electrodynamics in the spacetime of a rotating, conducting star and compute self-consistently to leading order in the rotation rate $\epsilon$. In this limit, the magnetosphere is completely described by the stellar mass, radius, moment of inertia, rotation rate, and magnetization, all of which may be independently chosen.

We apply our method to the inclined dipole as well as to the inclined “quadrudipole,” a superposition of dipole and quadrupole fields Barnard and Arons (1982); Gralla et al. (2016a). These are the first results on polar cap properties that include both general relativity and non-zero inclination. Our findings in the dipole case indicate the presence of pair production, and hence confirm the self-consistency of the dipole model. The quadrudipole case illustrates the phenomenon of thin annular polar caps, confirming that the shape seen in the aligned case Gralla et al. (2016a) persists for non-zero inclination. For axisymmetric (but possibly inclined) magnetic fields, these properties exhaust the allowed shapes: the polar caps will be either circular or annular. However, any shape is allowed for intrinsically nonaxisymmetric fields, a case we expect to consider in future work.

The charge-current distribution on the polar cap plays a direct role in several types of pulsar emission. First, thermal X-ray emission likely arises from heating of the polar cap by pulsar return current Arons (2012). Thus, our formulae can be used to model X-ray emission as a function of magnetic field geometry, and hence infer the magnetic geometry from observations (at least in principle). Our results can also be used to help check the robustness of upcoming neutron star radius measurements by NICER NIC (icer) to variation in intrinsic magnetic field. Second, the regions in which electron-positron pair production occurs are determined by bulk current flow Timokhin and Arons (2013). In most models, the pair production is directly related to radio emission and the formation of the pulsar wind. Given such a model, our results can predict the radio emission signature and pair loading of the pulsar wind as a function of magnetic field geometry. The gamma-ray emission is likely related to the magnetospheric current sheet Bai and Spitkovsky (2010); Cerutti et al. (2016), whose properties are largely insensitive to the magnetic field. As such, variations between gamma-ray and other types of emission can be a useful probe of magnetic field geometry.

The derivation of our results relies on the spacetime approach to force-free electrodynamics Gralla and Jacobson (2014), making extensive use of differential forms and focusing on invariant properties. Those readers who are unfamiliar with this approach, or otherwise uninterested in the derivation, may skip directly to Sec. III, where we state the assumptions and provide a detailed prescription for applying the method. Finally, in Sec. IV, we present results for the dipole and quadrupole pulsars. Our metric has signature $(-,+,+,+)$ and we use Heaviside-Lorentz units with $G=c=1$.

II Method

We consider the exterior metric of a rotating body, by which we mean a metric with a Killing field $\xi$ that takes the form $\mathop{}\!\partial_{t}+\Omega\mathop{}\!\partial_{\varphi}$ in cylindrical coordinates far from the star, with $\Omega$ a constant. We can always choose coordinates such that $\xi$ takes this form everywhere, in which case our assumptions become

[TABLE]

where $\mathcal{L}$ denotes the Lie derivative. We do not assume at this stage that the spacetime is separately stationary and axisymmetric. For example, a nonaxisymmetric rotating body (such as a neutron star with a “mountain” on it) would still possess the symmetry (1) when gravitational radiation is neglected. Wherever $\xi$ is timelike, we can define co-rotating observers with four-velocity

[TABLE]

We call $\Upsilon$ the redshift factor of the orbit $u^{\mu}$. For typical fluid stars, $\xi$ is timelike sufficiently close to the axis of rotation. The boundary of the timelike region, where $\xi$ becomes null, is called the light cylinder. It is the radius at which an observer co-rotating with the star would be moving at the speed of light.

Now suppose that the star is magnetized. If the intrinsic magnetization does not change with time (i.e., the only changes are due to rotation), then the electromagnetic field will share the symmetry of the metric. Working with the Maxwell two-form $F$, we therefore assume that

[TABLE]

If the electromagnetic field is degenerate ($F\wedge F=0$ or, equivalently, $\bm{E}\cdot\bm{B}=0$), then we may always introduce scalar potentials $\phi_{1}$ and $\phi_{2}$ such that Carter (1979); Uchida (1997a); Gralla and Jacobson (2014)

[TABLE]

These potentials are the relativistic generalization of Euler potentials. When the field is magnetically dominated ($F^{2}>0$ or, equivalently, $\bm{B}^{2}>\bm{E}^{2}$), the two-surfaces of constant $\phi_{1}$ and $\phi_{2}$ represent worldsheets of magnetic field lines, or “field sheets” for short Gralla and Jacobson (2014).

The symmetry (3) implies that $\phi_{1}$ and $\phi_{2}$ can be chosen such that Uchida (1997b); Gralla and Jacobson (2014)

[TABLE]

where $\kappa$ is a constant. The electric field measured by co-rotating observers is $F\cdot u=F\cdot\xi/\sqrt{\Upsilon}=\kappa\mathop{}\!\mathrm{d}\psi_{1}/\sqrt{\Upsilon}$. If the star is perfectly conducting, then this must vanish on the surface of the star, thereby fixing $\kappa=0$ (provided $F\neq 0$ on the star). Therefore, we may write

[TABLE]

In particular, the magnetosphere everywhere satisfies

[TABLE]

which means that the Killing field is tangent to the field sheets.111The notation $v\cdot\omega$ indicates the contraction of the vector $v$ into the first index of a differential form $\omega$. In this case, \mathopen{}\mathclose{{}\left(\xi\cdot F}\right)_{\nu}=F_{\mu\nu}\xi^{\mu}. In flat spacetime, Eq. (7) corresponds to the formula $\bm{E}+\bm{V}\times\bm{B}=0$ (see App. A.2 for details).

If the magnetosphere is also force-free ($J\cdot F=0$), then as explained in App. A, Eq. (7) implies the existence of a quantity $\Lambda$ that is constant on each field sheet. In terms of the three-current $\bm{J}_{c}$ and magnetic field $\bm{B}_{c}$ measured by the co-rotating observers, this conserved quantity may be written as

[TABLE]

These formulae hold only inside the light cylinder; outside the light cylinder, $\Lambda$ must be computed in a different manner.

The conserved quantity $\Lambda$ has an amusing cultural history. The flat spacetime version was discovered by Mestel (1973) and subsequently used in Beskin et al. (1983), before being rediscovered by Uchida (1998) and re-rediscovered by Gruzinov (2005). In App. A, we reveal the geometric interpretation and find the full generalization of $\Lambda$: we show that such a conserved quantity always exists in any spacetime $g$ and field configuration $F$ sharing a symmetry $\xi$ that satisfies $\xi\cdot F=0$.

II.1 Metric

We now specialize to the case where the stellar surface rotation velocity is much less than the speed of light ($\epsilon\ll 1$). To leading order in rotation, the metric outside of a relativistic star is Hartle and Thorne (1968)

[TABLE]

where the “frame-drag frequency” $\Omega_{Z}$ is

[TABLE]

Here, $R_{\star}$ is the (areal) radius, $M$ is the mass, $\Omega$ is the angular velocity, and $\hat{I}$ is the moment of inertia (defined as the angular momentum over the angular velocity). The norm of the co-rotation Killing field is

[TABLE]

The light cylinder is defined to be the locus where $\Upsilon$ vanishes. We will work with a characteristic light cylinder radius defined by

[TABLE]

which agrees with the actual radius in the slow-rotation limit.

We can characterize the problem by an overall scale and three dimensionless parameters,

[TABLE]

corresponding to the surface rotation velocity, the stellar compactness, and the dimensionless moment of inertia, respectively. The compactness has a theoretical upper limit $\mathcal{C}<8/9$ Buchdahl (1959), with realistic neutron star models having $\mathcal{C}\sim 1/2$ Haensel et al. (2007). For clarity of discussion, we will regard the metric (9) as exact, although it of course receives higher-order corrections in $\epsilon$. These corrections do not influence the leading-order calculations performed below.

II.2 Field strength and expansions

We wish to solve the equations of force-free electrodynamics perturbatively in $\epsilon$. Formally, we may imagine having solved the problem at finite $\epsilon$. This defines a family of solutions

[TABLE]

Perturbation theory consists of Taylor-expanding in $\epsilon$, but we have a choice of what to hold fixed. We define two expansions,

[TABLE]

with the other dimensionless parameters $\mathcal{C}$ and $\mathcal{I}$ fixed in both. We introduce order symbols $\mathcal{O}_{\star}$ and $\mathcal{O}_{L}$ to represent the scalings in each limit. As a trivial example, we have

[TABLE]

Other dimensionful parameters scale as

[TABLE]

(For example, $M=\mathcal{C}R_{\star}/2=\epsilon\mathcal{C}R_{L}/2$.) Since the coordinates are held fixed in both expansions, the scalings of the dimensionful parameters define coordinate regions corresponding to each expansion. Since the metric and field strength depend on $r/M$, $\hat{I}/r^{3}$, and $\Omega t$, the regimes of validity are222Note that $r\gg M$ implies $R\gg R_{\star}$ because general relativity does not allow stars to become too compact (i.e., we have $R_{\star}\lesssim M$). This condition is also implied by the metric (9).

[TABLE]

(The arbitrary choice of working near $t=0$ is inherited from the choice of regarding the Euler potentials as functions of $\varphi-\Omega t$.) Both expansions are valid in the

[TABLE]

which exists because $\epsilon=R_{\star}/R_{L}\ll 1$. This observation allows the large-$r$ behavior of the near expansion to be matched to the small-$r$, small-$t$ behavior of the far expansion.

II.3 Near zone

In the near expansion, the geometry reduces to the Schwarzschild metric plus corrections due to rotation,

[TABLE]

where we remind the reader that $\Omega_{Z}=\mathcal{O}_{\star}(\epsilon)$. Recalling that $\Omega=\epsilon/R_{\star}=\mathcal{O}_{\star}(\epsilon)$, the near expansion of the Euler potentials takes the form

[TABLE]

where $\alpha,\beta,\tilde{\alpha},\tilde{\beta}$ are functions of the spatial coordinates only. Then the field strength is written as

[TABLE]

where we remind the reader that $\Omega=\epsilon/R_{\star}=\mathcal{O}_{\star}(\epsilon)$. The first term in Eq. (II.3) represents the magnetic field of the star when it is not rotating. We denote it by

[TABLE]

The remaining terms in (II.3) are $\mathcal{O}_{\star}(\epsilon)$. The second term represents the leading electric field induced by the rotation. The third term (the second line) contains some time-dependence due to the rotation of the star. The final term (the third line) represents corrections to the magnetic field of the star.

We will make the basic assumption that the currents are due only to rotation,

[TABLE]

This describes an isolated pulsar, where current flows only because of unipolar induction due to the rotating conductor. In the aligned case, Eq. (25) can be proven from the assumption of asymptotically radial field lines (e.g., Eq. (12) of Paper I), and it is supported by numerical simulations done in isolation, including those described in App. B.

Equation (25) means that $F^{(0)}$ is a vacuum (no charge or current) Maxwell solution in the Schwarzschild metric. Since $\alpha$ and $\beta$ are independent of time, this solution is also stationary and purely magnetic. The stationary, magnetic, vacuum solutions in Schwarzschild spacetime are known in closed form Anderson and Cohen (1970). Provided the dipole moment $\mu$ is non-zero, at sufficiently large radius the field always becomes dipolar, and without loss of generality, we may take the dipole moment to instantaneously point along the $x$ axis at an inclination $\iota$ relative to the rotation axis ($z$ axis). This means that the Euler potentials may always be chosen such that

[TABLE]

where $\theta^{\prime}$ and $\varphi^{\prime}$ are polar coordinates about the dipole axis, related to $\theta$ and $\varphi$ by

[TABLE]

Note that in Eqs. (26), the regime $r\to\infty$ does not refer to asymptotic infinity because we have already taken the near limit. Instead, $r\to\infty$ refers to the overlap region where we will match to the far zone.

Our choice of a dipole field (26) in the overlap region corresponds to stars whose dipole component of field dominates before the light cylinder is reached, though the dipole can still be sub-dominant at the stellar surface. For a star with a different moment dominating in the overlap region, Eq. (26) should be suitably modified to include that moment. The rest of the analysis presented in this paper (including a new numerical simulation and fit for $\Lambda$) can then be repeated for such a case.

II.4 Far zone

In the far expansion, the metric becomes flat,

[TABLE]

The Euler potentials do not simplify, and we simply name the leading piece $\chi_{i}\equiv\psi_{i}(r,\theta,\varphi-\Omega t;0)$, so that

[TABLE]

Then, the field strength is

[TABLE]

The small-$r$, small-$t$ behavior of the far expansion must match the large-$r$ behavior of the near expansion. That is, Eqs. (26) require the boundary condition

[TABLE]

These equations hold at $t=0$, an arbitrary choice of time at which we chose to align the dipole moment with the $x$ axis [see Eqs. (26) and (27)]. By the co-rotation symmetry $\xi=\mathop{}\!\partial_{t}+\Omega\mathop{}\!\partial_{\varphi}$ of the system, similar equations hold at other times, except that $\theta^{\prime}$ and $\varphi^{\prime}$ must be aligned with the rotating dipole axis. Thus, the small-$r$ boundary condition in the far zone is a rotating point dipole.

The large-$r$ boundary condition is that the magnetosphere is isolated. In practice, it is handled numerically by making the simulation box large enough that edge effects cannot affect the solution before it reaches steady state.

We see that the leading-order field in the far expansion is the force-free magnetosphere of an isolated, rotating, conducting, inclined point dipole in flat spacetime. This field can be determined numerically by considering a sequence of increasingly smaller rotating conducting stars. We can then fit for the conserved quantity $\Lambda$ as a function of $\chi_{1}$ and $\chi_{2}$. In the numerical simulation, this is most conveniently done near or on the star, which by definition is in the overlap region where $\chi_{1}=\alpha$ and $\chi_{2}=\beta$ of the matched asymptotic expansion [see Eqs. (26) and (31)].

In finding a suitable fit, our goal is to capture the qualitative features while maintaining simplicity. We find that an excellent fit (see App. B) is given by

[TABLE]

where $J_{0}$ and $J_{1}$ are Bessel functions of the first kind and the “last open field line” $\alpha_{o}$ is given by

[TABLE]

The upper/lower sign in Eq. (II.4) corresponds to the northern/southern flow.

Our motivations for this fit are the following. In the limit $\epsilon\to 0$, the polar cap occupies a vanishingly small portion of the sphere, and may therefore be approximated as a disk. The natural fitting functions are then the Bessel harmonics $J_{n}(\rho)e^{\pm in\gamma}$ (for $n\geq 0$ with $(\rho,\gamma)$ polar coordinates in the disk). These harmonics form representations of the Euclidean group $\mathsf{E}_{2}$ and are the disk analogs of spherical harmonics. The angular coordinate $\gamma$ must equal $\phi^{\prime}=\beta$ by axisymmetry. The radial coordinate $\rho$ should be proportional to $\theta^{\prime}$, but the proportionality constant is free since we have not yet fixed the size of the disk. Noting that $\alpha\propto\sin^{2}{\theta^{\prime}}$ because the field is dipolar in this region, we have $\theta^{\prime}=\arcsin{\sqrt{\alpha/\alpha_{o}}}$ for some constant $\alpha_{o}$. We find that $\rho=2\theta^{\prime}$ works well, resulting in the argument $2\arcsin\sqrt{\alpha/\alpha_{o}}$ for the Bessel functions. The $\arcsin$ function is undefined for $\alpha>\alpha_{o}$, so $\alpha_{o}$ naturally delineates the boundary of the polar cap. We fit the simple functional form (33) to the size of the polar cap. Remarkably, an excellent fit for the current $\Lambda$ is now obtained by simply taking a linear combination of the Bessel harmonics with prefactors $\cos{\iota}$ and $\sin{\iota}$, as shown in Eq. (II.4).

Interestingly, an analogous expression involving Bessel functions provides a qualitative description Gralla and Jacobson (2015) of the force-free jets produced by a spinning black hole embedded in a misaligned magnetic field Palenzuela et al. (2010).

The simple fit (II.4) is adequate for the applications we have in mind, and hence we use it for the remainder of the paper. However, we emphasize that improvements to the simulation and fitting could surely produce a more accurate expression for $\Lambda(\alpha,\beta)$. It is straightforward to repeat the analysis of the rest of the paper for any such expression.

II.5 Near-zone charge-current

We may now determine the near-zone four-current to leading order $\mathcal{O}_{\star}(\epsilon)$. We can obtain three components from the conservation law \bm{J}=\mathopen{}\mathclose{{}\left(\Lambda/\sqrt{\Upsilon}}\right)\bm{B} [see Eq. (8)]. Since $\Lambda=\mathcal{O}_{\star}(\epsilon)$ and $\bm{J}=\mathcal{O}_{\star}(\epsilon)$, we need only the $\mathcal{O}_{\star}(1)$ piece of $\bm{B}$ (i.e., $F^{(0)}=\mathop{}\!\mathrm{d}\alpha\wedge\mathop{}\!\mathrm{d}\beta$) to determine the leading, $\mathcal{O}_{\star}(\epsilon)$ current. A straightforward calculation yields

[TABLE]

The fourth component can be obtained by taking the divergence $J^{\nu}=\nabla_{\mu}F^{\mu\nu}$ directly from Eq. (II.3),

[TABLE]

Notice that the result involves only $\alpha$ and $\beta$; the higher-order corrections $\tilde{\alpha}$ and $\tilde{\beta}$ do not appear.333In flat spacetime, this corresponds to the statement that we may compute $\rho_{e}=\bm{\nabla}\cdot\bm{E}$ from knowledge of $\bm{B}$ only, because $\bm{E}=-\bm{V}\times\bm{B}$ with $\bm{V}=\rho\hat{\bm{\phi}}$. This is essential to the method, since the higher-order corrections are unknown. Similar direct computations of the spatial components would result in expressions featuring $\tilde{\alpha}$ and $\tilde{\beta}$. These expressions could be set equal to the current found from the conservation of $\Lambda$ [Eq. (34)] to provide equations for the corrections $\tilde{\alpha}$ and $\tilde{\beta}$.

According to Eqs. (34) and (35), the charge-current is only non-zero where the conserved quantity $\Lambda(\alpha,\beta)$ is non-zero. There is a subtlety, however, in that $\Lambda$ is a conserved quantity along field sheets (equivalently, field lines), but $\alpha$ and $\beta$ may not uniquely label the field lines. The charge-current (II.4)-(33) is determined in the overlap region, and hence charge and current flow only on the field lines that reach the overlap region (which are also the field lines that open up to asymptotic infinity). Thus, in case there are multiple field lines with $\alpha<\alpha_{o}$ in the near region, one is to use Eqs. (34) and (35) only for the field lines that reach the overlap region where the labeling becomes unique. On all other field lines, the current flow is zero. This ensures that there are only two polar caps; see Paper I for further discussion of the axisymmetric case.

III Statement of method

We now summarize the assumptions and describe the method we have derived. We have assumed that:

The star is spherical (mass $M$ and areal radius $R_{\star}$) and rigidly rotating with angular velocity $\bm{\Omega}=\Omega\hat{\bm{z}}$ such that $\Omega R_{\star}\ll c$. 2. 2.

The stellar surface is perfectly conducting and the magnetosphere is force-free. The star is isolated. 3. 3.

The star is magnetized in such a way that when it is not rotating, the magnetic field is dipolar at large distances. Furthermore, the angular velocity is small enough that the dipole component dominates before the light cylinder radius $r\sim R_{L}=c/\Omega$ is reached.

Under these assumptions, we provide analytic formulae for the near-field charge and current associated with a given choice of stellar parameters (13) and magnetization. The magnetization is described by a static, vacuum, asymptotically dipolar magnetic field solution in the Schwarzschild spacetime. This solution must be expressed in terms of Euler potentials $\alpha(t,r,\phi)$ and $\beta(t,r,\phi)$ in Schwarzschild coordinates. (That is, the Schwarzschild coordinate components of $F$ are given by $F_{\mu\nu}=\mathop{}\!\partial_{\mu}\alpha\mathop{}\!\partial_{\nu}\beta-\mathop{}\!\partial_{\nu}\alpha\mathop{}\!\partial_{\mu}\beta$.) For axisymmetric fields (not necessarily aligned with the rotation axis), one can accomplish this by picking an axisymmetric flux function $\psi(r,\theta)$ (e.g., App. B2 of Ref. Gralla et al. (2016b)) and setting

[TABLE]

where $(\theta^{\prime},\varphi^{\prime})$ denote spherical coordinates rotated to the desired inclination angle $\iota$ [see Eqs. (27)].

The Euler potentials $\alpha$ and $\beta$ label magnetic field lines. The labeling is not necessarily one-to-one, but it becomes so at sufficiently large $r$ as the field becomes dipolar. Current flows on the field lines that have $\alpha<\alpha_{o}$ [see Eq. (33)] and enter the region where the labeling is unique.444Near the star, there can be additional field lines labeled by $\alpha<\alpha_{o}$; these field lines do not connect to the far zone, and hence no current flows on them. The portion of space occupied by these field lines has non-zero charge and current given by555Here, we present orthonormal-frame components of static Schwarzschild observers, using the notation and definitions of Ch. II of Ref. Thorne et al. (1986). This agrees with the measurements of co-rotating observers to leading order in the rotation $\Omega R_{\star}$.

[TABLE]

where

[TABLE]

with $J^{t}$ given in Eq. (35) and $\Lambda$ given in Eq. (II.4). These formulae express the charge and current at leading order in $\epsilon$ near the star ($r\ll R_{L}$). When evaluated on the star, they provide the polar cap structure.

IV Results: dipole and quadrudipole

We now present results for two specific choices of stellar magnetization. First, we consider the canonical case of a pure dipole. The flux function $\psi$ is given (e.g.) in Eq. (40) of Paper I. Following Eq. (36), the Euler potentials are given by

[TABLE]

(Although it is not apparent, this solution does approach $\mu\sin^{2}{\theta^{\prime}}/r$ at large $r$, as required.) In Fig. 1, we plot the norm of the charge-current, $J^{2}=J_{\mu}J^{\mu}=\bm{J}^{2}-\rho_{e}^{2}$, for a variety of parameters. The main result is that for realistic compactness, there are regions of spacelike current ($J^{2}>0$) at all inclinations. Since spacelike current ensures the pair production necessary to sustain a high-multiplicity magnetosphere (Timokhin and Arons, 2013; Chen and Beloborodov, 2014; Philippov et al., 2015a, b), our results support the basic consistency of the dipole force-free model. We also show flat spacetime (zero-compactness) results to compare with previous work and illustrate the role of general relativity. These are the first results for the polar cap structure of the general relativistic force-free inclined dipole pulsar.

As a more complicated example, we consider a quadrudipole field (a superposition of dipole and quadrupole fields). We choose the moment ratio $q=3$ in the notation of Sec. IV of Paper I. The flux function $\psi(\theta)$ (and hence the Euler potential $\alpha=\psi(\theta^{\prime})$) is given by Eq. (50) of Paper I. Instead of reproducing the formula, we include a plot of the level sets of $\alpha$ (i.e., the poloidal field lines) in Fig. 2.

The current norm $J^{2}$ is plotted in Fig. 3 for a variety of inclinations. Notice that the southern polar cap takes a thin annular shape, a typical feature of non-dipolar fields. We refer the reader to Paper I for a discussion of the potential importance of this feature, including the possibility that it accounts for the modified beam characteristics of millisecond pulsars.

We now discuss an interesting feature of the orthogonal case $\iota=90^{\circ}$. Notice that when the polar cap is circular (dipole or northern cap of quadrudipole), the current is entirely spacelike, while for the annular cap there are both timelike and spacelike currents. This can be understood from the fact that the charge density scales like \rho_{e}\sim\mathopen{}\mathclose{{}\left(\bm{\Omega}-\bm{\Omega}_{Z}}\right)\cdot\bm{B}. In the orthogonal case, the circular polar caps lie in a region where the rotation and magnetic field become orthogonal, making the charge density anomalously small. Thus, the four-current is entirely spacelike. On the other hand, the annular polar cap is shifted away from this region, and both kinds of currents are present.

Finally, we discuss the presence of volume return current. By “return current” we mean a region near the star satisfying $\mathbf{J}\cdot\mathbf{n}/\rho_{e}<0$, where $\mathbf{n}$ is the co-rotating normal vector to the star and the charge and current are those measured by co-rotating observers. This corresponds to inward flow of plasma when only a single sign of charge is present. Noting that $\Lambda\sim\bm{J}\cdot\bm{B}$, the regions of volume return current are those with $\pm\Lambda/\rho_{e}<0$, where $+$ and $-$ correspond to the northern ($\mathbf{B}\cdot\mathbf{n}>0$) and southern ($\mathbf{B}\cdot\mathbf{n}<0$) caps, respectively.

We describe the qualitative features rather than presenting plots. For non-orthogonal inclinations, the charge density does not change sign on the polar caps. We find that there is a small region of volume return current at the outer edge (the edge further from the magnetic pole) in all the non-orthogonal cases we consider. In the orthogonal case, the charge density changes sign, but we find that there is no volume return current for both the dipolar and quadrudipolar magnetic geometries.

Acknowledgements

We thank Feryal Özel, Dimitrios Psaltis, Anatoly Spitkovsky, and Alexander Tchekhovskoy for helpful conversations. This work was supported by NSF grants 1205550 to Harvard University and 1506027 to the University of Arizona, by the NASA Earth and Space Science Fellowship Program (grant NNX15AT50H to A.P.) and by the Porter Ogden Jacobus Fellowship, awarded to A.P. by the graduate school of Princeton University. The simulations presented in this paper used computational resources supported by the PICSciE-OIT High Performance Computing Center and Visualization Laboratory, and by the NASA/Ames HEC Program (SMD-16-6663, SMD-16-7816).

Appendix A Conserved quantity

Force-free electrodynamics in an arbitrary curved spacetime may be expressed in terms of three scalar potentials $(\phi_{1},\phi_{2},\lambda)$ as (see footnote 6 of Ref. Gralla and Jacobson (2014))

[TABLE]

Here, $\phi_{1}$ and $\phi_{2}$ are relativistic Euler potentials. There is no standard name for $\lambda$, but it acts as a stream function for the charge-current as measured in units of magnetic field strength. In particular, an observer with four-velocity $u^{\mu}$ measures current density $\bm{J}$ and magnetic field $\bm{B}$ related by

[TABLE]

From studies of stationary force-free magnetic fields, we are used to this field-aligned current $J_{\parallel}$ being conserved along field lines. To generalize this statement, we suppose that the metric and field configuration have a symmetry tangent to the field sheets,

[TABLE]

The current three-form $J=\mathop{}\!\mathrm{d}\star F$ will also respect the symmetry, so

[TABLE]

This means that $\xi\cdot\mathop{}\!\mathrm{d}\lambda$ is a function only of $\phi_{1}$ and $\phi_{2}$,

[TABLE]

If $\xi$ is timelike, then we may define an associated observer $\mathcal{O}$ with four-velocity $u^{\mu}$ and redshift factor $1/\sqrt{\Upsilon}$,

[TABLE]

It is then straightforward to show that the magnetic field $\bm{B}_{\mathcal{O}}$ and three-current $\bm{J}_{\mathcal{O}}$ seen by these observers satisfy

[TABLE]

Thus, in regions where the Killing field $\xi$ is timelike, we may interpret the conserved quantity $\Lambda$ as a redshifted field-aligned current. (Note that there is no component of field perpendicular to the current, as guaranteed by force-free electrodynamics.) Note also that by Eq. (42), we may always choose the potentials to satisfy Uchida (1997b); Gralla and Jacobson (2014)

[TABLE]

We next compare to previous work, showing how this conserved quantity generalizes those previously known in various different contexts.

A.1 Force-free magnetic fields

If spacetime is flat and $\xi=\mathop{}\!\partial_{t}$ is the time-translation Killing field, then our assumptions (42) correspond to

[TABLE]

a setup generally known as a force-free magnetic field. In this case, $\Upsilon=1$ identically, so the conserved quantity $\Lambda$ is just the field-aligned current. By Eq. (47), we may always choose $\phi_{1}$ and $\phi_{2}$ to be independent of $t$, making them the usual Euler potentials, which label the field lines as $\bm{B}=\bm{\nabla}\phi_{1}\times\bm{\nabla}\phi_{2}$. Thus, we recover the standard result that $J_{\parallel}$ is constant on field lines.

A.2 Rotating configuration in flat spacetime

If $\xi=\mathop{}\!\partial_{t}+\Omega\mathop{}\!\partial_{\phi}$ in flat spacetime, then the assumptions (42) correspond to

[TABLE]

where $\bm{V}=\omega\rho\hat{\bm{\phi}}$. Inside the light cylinder, $\rho<\Omega^{-1}$, we may define co-rotating observers with a four-velocity and redshift

[TABLE]

The current density and magnetic field according to these observers are

[TABLE]

where $\bm{J}$ and $\rho_{e}$ are the current and charge densities in the fixed frame. Thus, by Eq. (46), we may write

[TABLE]

the form given by Uchida (1998). The left side of this equation is sometimes called the current in a rotating frame. We will avoid this terminology since we reserve “frame” for the physical measurements of some class of observers. An alternative form is given by Gruzinov (2006),

[TABLE]

Equations (52) and (53) are equivalent under the assumptions (49) (and the force-free equations).

A.3 Axisymmetric solutions

In a general axisymmetric (circular) spacetime, one may characterize axisymmetric solutions with non-zero poloidal field by the flux function $\psi(r,\theta)$, the polar current $I(\psi)$, and field line velocity $\Omega(\psi)$. We adopt the conventions of Paper I, where the flux function is the magnetic flux through a loop of revolution divided by $2\pi$, and the polar current is minus the current through the loop of revolution. Then, the conserved quantity $\Lambda$ is

[TABLE]

For example, the aligned dipole pulsar has $I(\psi)$ given to an excellent approximation by Gralla et al. (2016a)

[TABLE]

The $+$ refers to field lines that asymptote to the northern hemisphere, while the $-$ refers to field lines that asymptote to the southern hemisphere. The expression holds for $\psi<\psi_{o}$, with $\psi_{o}$ given by

[TABLE]

Thus, in this case, the conserved quantity is given by

[TABLE]

where the upper/lower sign refers to the northern/southern flow.

Appendix B Simulations and fit

We carried out a number of time-dependent 3D simulations of oblique pulsar magnetospheres in flat spacetime in the force-free approximation with the code by Spitkovsky (2006) in a Cartesian grid. We consider a perfectly conducting star of radius $R_{\star}$ with a magnetic dipole field of dipole moment $\mu$ that makes an angle $\iota$ with the rotational axis. The star rotates with angular velocity $\Omega$. We performed simulations for a range of inclination angles \iota=\mathopen{}\mathclose{{}\left\{0^{\circ},30^{\circ},60^{\circ},90^{\circ}}\right\} and rotation \epsilon=\mathopen{}\mathclose{{}\left\{0.4,0.2,0.1,0.067}\right\}. The stellar radius is resolved by 40 computational cells in all simulations. We find that the conserved quantity is essentially converged, in that $\Lambda(\alpha,\beta)/\epsilon$ changes very little between $\epsilon=0.1$ and $\epsilon=0.067$. The rescaled polar cap size $\alpha_{0}/\epsilon$ shows a small increase with decreasing $\epsilon$, tending toward the reported value (33). The non-uniformity of the asymptotic magnetic field Tchekhovskoy et al. (2016), which is linked to the functional form of $\Lambda(\alpha,\beta)$, does not depend on $\epsilon$. The increase of the polar cap size (or the value of the open magnetic flux) with decreasing $\epsilon$ leads to an increase in the spin-down power Philippov et al. (2014). We perform the fit for $\epsilon=0.067$. Figure 4 shows the results.