A Split monopole

# Pulsar Magnetospheres: Beyond the Flat Spacetime Dipole

## Abstract

Most studies of the pulsar magnetosphere have assumed a pure magnetic dipole in flat spacetime. However, recent work suggests that the effects of general relativity are in fact of vital importance and that realistic pulsar magnetic fields will have a significant nondipolar component. We introduce a general analytical method for studying the axisymmetric force-free magnetosphere of a slowly rotating star of arbitrary magnetic field, mass, radius, and moment of inertia, including all the effects of general relativity. We confirm that spacelike current is generically present in the polar caps (suggesting a pair production region), irrespective of the stellar magnetic field. We show that general relativity introduces a correction to the formula for the dipolar component of the surface magnetic field inferred from spindown. Finally, we show that the location and shape of the polar caps can be modified dramatically by even modestly strong higher moments. This can affect emission processes occurring near the star and may help explain the modified beam characteristics of millisecond pulsars.

## I Introduction

Pulsars display a striking and diverse observational phenomenology, including pulsed emission across a variety of wavelengths and synchrotron emission in associated pulsar wind nebulae. The challenge to the pulsar theorist is to explain this array of observations from the rather simple theoretical beginnings of a rotating, magnetized, neutron star. After more than 40 years of effort, much has been understood, but the full solution remains elusive Arons (2012).

The presence of plasma outside the star is likely a key ingredient. Early on, it was realized that the rotation of the star in vacuum produces enormous “accelerating” electric fields (i.e., fields with component along the magnetic field) outside the star Goldreich and Julian (1969). Unless the work function is unnaturally large, these fields lift charged particles off the star and accelerate them to ultra-relativistic energies. The emitted curvature photons can then interact with the magnetic field (or other charged particles) to produce particle pairs Sturrock (1971), which act to screen the accelerating field. If this mechanism continuously produces charges to replace those that are accelerated away, one can hope to explain both the pulsed radio emission (from collective behavior in the bulk plasma) as well as the high energy particles in the wind. To fill in the details, however, one needs (at a minimum) the magnetic geometry, which must be determined self-consistently by solving for the global structure of the plasma magnetosphere.

Considerable progress has been made in constructing global models of pulsar magnetospheres under the assumption of abundant plasma supply. Since the energy of the magnetospheric plasma is much smaller than the energy in the magnetic field, the magnetohydrodynamic equations can be simplified to the low-inertia, or force-free, limit. For example, self-consistent force-free numerical solutions of axisymmetric Contopoulos et al. (1999); Gruzinov (2005); McKinney (2006); Timokhin (2006) and oblique Spitkovsky (2006); Kalapotharakos and Contopoulos (2009); Li et al. (2012); Kalapotharakos et al. (2012); Pétri (2012) pulsar magnetospheres were developed. These studies were extended to include the effects of plasma pressure and inertia Komissarov (2006); Tchekhovskoy et al. (2013), as well as general relativity Lehner et al. (2012); Palenzuela (2013); Ruiz et al. (2014); Pétri (2016), which are important in current layers and near the star, respectively. This has produced a canonical model for the dipole pulsar magnetosphere, with polar outflows along open field lines separated from a zone of closed field lines by thin return current layers, together with a strong current sheet beyond the light cylinder. Despite this success, the basic questions remain: is the plasma really there? And can it produce pulsed emission and relativistic particles for the wind?

Addressing these questions from first principles requires a full kinetic treatment, including the physics of pair formation and particle acceleration. Recent years have seen significant progress in kinetic plasma simulations, which have shed light on how and when plasma can be produced in a magnetosphere. Working with a dipolar field in flat spacetime, Refs. Chen and Beloborodov (2014) and Philippov et al. (2015a) found that, provided a source of particles at the stellar surface and allowing pair production close to the light cylinder, a nearly force-free magnetosphere develops in which pairs are produced in the thin current layers and the equatorial current sheet. Notably, however, there is no pair formation in the bulk of the polar cap in these models. Instead, the accelerating electric field is screened by a mildly relativistic charge-separated flow of particles lifted from the neutron star surface Shibata (1997); Beloborodov (2008); Timokhin and Arons (2013).

This slow outflow does not radiate the gamma photons needed to produce pairs, and is not promising for producing radio emission or pulsar wind particles. However, Ref. Philippov et al. (2015b) showed that the inclusion of general relativistic (GR) effects, in particular the dragging of inertial frames around the neutron star, leads to the ignition of a discharge at the polar cap of the aligned rotator. This phenomenon is associated with the current density becoming larger than the charge density (times the speed of light), which can no longer be supported by one sign of escaping charge and requires pair production. In relativistic language, the charge-current four-vector is spacelike. Since the configuration is force-free, this suggests a simple prescription for studying magnetospheric discharges Timokhin and Arons (2013): solve the equations of force-free electrodynamics and find the regions near the star where the current is spacelike.

Reference Philippov et al. (2015b) studied the aligned dipole in GR and found qualitative agreement between the spacelike-current regions of a force-free simulation and the pair-formation regions in the kinetic simulation. This supports the method, but much remains to be explored. In particular, there is little reason to believe that real pulsars have pure dipolar magnetic fields. Indeed, the magnetic field is determined by complicated internal processes, and to the extent that these are understood, the field is not dipolar Geppert and Viganò (2014); Gourgouliatos et al. (2016). In this paper, we introduce an analytical procedure to determine the near-field magnetosphere, including the globally determined current flow, of a slowly rotating star with an arbitrary axisymmetric magnetic field profile. We use the method of matched asymptotic expansions, solving the equations of force-free electrodynamics on curved spacetime to linear order in , the ratio of the stellar radius to the light cylinder radius.

We show that spacelike current generically appears in the current flow regions near the star, or “polar caps,” thereby extending the dipolar results of Ref. Philippov et al. (2015b) to generic fields. We show, however, that the location and shape of the polar caps changes significantly when higher multipoles are included. In particular, a polar cap can occupy a thin annular region of angular width displaced from the rotational pole instead of the classic region surrounding the pole. This could give rise to double-peaked emission and/or help explain the discrepant beam size scalings between ordinary and millisecond pulsars Kramer et al. (1998). Such shifted polar caps are also further favorable for pair production because of higher field line curvature and subsequent increase of the magnetic opacity Barnard and Arons (1982). In general, our analytical method provides the detailed near-field magnetosphere and may hence be used to construct models for near-zone pair creation and emission based on more complicated and realistic magnetic fields Arons (2012).

We also consider the spindown luminosity of the magnetosphere. We show that the luminosity depends only on the dipolar component of the surface magnetic field. We find that general relativity significantly decreases the luminosity for a fixed magnitude of the surface magnetic field, corresponding to a correction to the formula for the surface magnetic field inferred from measurements of the pulsar period and period derivative. This correction is distinct from that found recently in Refs. Ruiz et al. (2014); Pétri (2016), which occurs only for rapid spin. In these works, the comparison to flat spacetime was made at fixed asymptotic dipole moment. Defined this way, there are no corrections due to general relativity for slowly spinning pulsars. In general, simulations are limited to larger values of , whereas our analytical method explores the complimentary regime , which is valid for most pulsars.

The paper is organized as follows. In Sec. II, we describe our matched asymptotic expansion method and derive the spindown luminosity as well as general formulas for the polar cap regions and the near-zone current. We then apply the method to pulsars with dipole fields (Sec. III), both dipole and quadrupole fields (Sec. IV), and more general configurations suggested by recent work on neutron star magnetic fields (Sec. V). In Sec. VI, we summarize and discuss observational implications. We use Heaviside-Lorentz units and set .

## Ii General Method

We consider a slowly rotating star in general relativity. To linear order in the rotation velocity, the star is spherical and the metric outside is given by Hartle and Thorne (1968)

 ds2 =α2dt2+α−2dr2 (1) +r2[dθ2+sin2θ(dϕ−ΩZdt)2],r>R⋆,

where the “redshift factor” and the “frame-drag frequency” are

 α =√1−2Mr,ΩZ=2^Ir3Ω. (2)

Here, is the (areal) radius, is the mass, is the angular velocity, and is the moment of inertia (defined as the angular momentum over the angular velocity). The angular velocity defines the “light cylinder” on which a co-rotating observer would move with the speed of light. We define

 RL=1Ω, (3)

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

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

 ϵ=R⋆RL,C=2MR⋆,I=^IMR2⋆, (4)

corresponding to the surface rotation velocity, the stellar compactness, and the dimensionless moment of inertia, respectively. The compactness is subject to the theoretical bound Buchdahl (1959); realistic neutron star models generally have Haensel et al. (2007). Note that on the surface of the star, we have

 α=√1−C,ΩZ=CIΩ. (5)

The metric (1) is valid to linear order in . In the remainder of the paper, we will compute physical quantities to leading nontrivial order in .1

If the star is surrounded by a stationary, axisymmetric, force-free magnetosphere, then the electromagnetic field is characterized by the magnetic flux function , the field line rotation velocity , and the polar current as Gralla and Jacobson (2014)

 F=I(ψ)2πα2sinθdr∧dθ+dψ∧[dϕ−ΩF(ψ)dt]. (6)

We use the conventions of Refs. Tchekhovskoy et al. (2010); Gralla et al. (2015). For the field configuration to be smooth, the flux function and polar current must vanish quadratically in at both poles. Level sets of the flux function are called poloidal field lines, or field lines for short.

We assume that the star is perfectly conducting, which fixes the field line rotation to the stellar rotation,

 ΩF(ψ)=Ω. (7)

The force-free condition implies that and satisfy

 r2∂r(α2χ∂rψ) +sinθ∂θ(χsinθ∂θψ) (8) = −(r2πα)2I(ψ)I′(ψ),

where we introduced

 χ=1−(rsinθα)2(Ω−ΩZ)2. (9)

The star’s magnetic field provides a boundary condition

 ψ∣∣r=R⋆=ψ⋆(θ), (10)

where is the flux function on the stellar surface. Since the pulsar is isolated, the field lines should become radial at large , i.e.,

 ψ=A(θ)+O(r−1), (11)

for some to be determined by solving the equations. Together with the stream equation (8), this condition implies

 I=2πΩsinθ∂θψ,r→∞, (12)

which is equivalent to having electric and magnetic fields of asymptotically equal magnitude as .

This defines the problem: Given a choice of stellar parameters , one must solve Eq. (8) for and satisfying the boundary conditions (10) and (12). We will use the method of matched asymptotic expansions, solving separately in the near zone and the far zone , and matching the two in the overlap region . More formally, the near (respectively, far) expansions are fixing (respectively, ), and we match the large- behavior of the near expansion to the small- behavior of the far expansion.

### ii.1 Near-zone

In the near-zone , Eq. (8) becomes

 r2∂r(α2∂rψ)+sinθ∂θ(∂θψsinθ)=0. (13)

This is the vacuum magnetic stream equation in the Schwarzschild spacetime. We will label its solutions with “near.” The general solution vanishing at both poles and falling off at large is given in a multipole expansion as (see App. B)

 ψnear(r,θ)=∞∑ℓ=1CℓR>ℓ(r)Θℓ(θ). (14)

The normalization of the radial functions is chosen so that as .

To satisfy the boundary condition (10) on the star, we must have

 ψnear∣∣r=R⋆ =ψ⋆(θ). (15)

Provided the dipole moment is nonzero, the dipole term will always dominate at sufficiently large ,

 ψnear=μrsin2θ,r→∞. (16)

### ii.2 Far-zone

In the far-zone (where also ), Eq. (8) becomes

 r2∂r[(1−Ω2r2sin2θ)∂rψ] +sinθ∂θ[(1−Ω2r2sin2θ)sinθ∂θψ] (17) =−(r2π)2I(ψ)I′(ψ).

This is the stream equation in flat spacetime. In general, it must be solved numerically. We will label its solutions with “far.” To match the boundary condition (11), we must have asymptotically radial field lines

 ψfar=A(θ)+O(r−1), (18)

for some function determined by solving the equations. To match to the near solution (16) in the overlap region, the solution must satisfy

 ψfar=μrsin2θ,r→0. (19)

The solution to Eq. (17) satisfying Eqs. (18)-(19) corresponds to the small- limit of the classic dipole pulsar problem in flat spacetime. We use the force-free version of the HARM code Gammie et al. (2003); Tchekhovskoy et al. (2007) and set the inflow velocity into the current sheet to zero in order to minimize the numerical dissipation McKinney (2006). The simulation setup is described in Ref. Tchekhovskoy et al. (2016). Performing a sequence of simulations with decreasing , we find good convergence by . The resulting flux function is shown in Fig. 1. The configuration is symmetric about the equatorial plane and the last open field line(s) have the value Timokhin (2006)

 ψo=nμΩ,n≈1.23. (20)

There is no current flowing on the closed field lines,

 I(ψ)=0,ψ>ψo, (21)

while the current on the open field lines is given to an excellent approximation (see Fig. 2) by

 I(ψ)=±2πΩψ[2−ψψo−15(ψψo)3],ψ<ψo, (22)

where we take the sign (respectively, sign) for the northern (respectively, southern) hemisphere.

The first two terms in the brackets are the current for the split monopole [Eq. (61b)]; we see that the dipole provides a quartic correction. This correction describes the volume return current, which flows close to the edge of the polar cap. While Eq. (22) was obtained in the far expansion, it is independent of and hence holds everywhere.

### ii.3 Stellar multipole moments

To characterize the field on the star, it is useful to introduce stellar multipole moments by

 ψ⋆(θ)=R2⋆∞∑ℓ=1BℓΘℓ(θ). (23)

The factor of gives these moments units of magnetic field. The stellar moments are related to the moments appearing in Eq. (14) by

 Bℓ=CℓRℓ+2⋆Δℓ(C), (24)

where is a dimensionless function of and (and hence depends only on the compactness ). The precise form of any can be determined from the analysis of App. B. For example, the dipole and quadrupole correction factors are

 Δ1 =−32C3[C(C+2)+2log(1−C)], (25) Δ2 =−203C2[4−(4−3C)Δ1]. (26)

By construction, all the satisfy

 limC→0Δℓ=1. (27)

At the source and field moments are related as they are in flat spacetime, so we may view as the fractional correction due to general relativity.

### ii.4 Range of validity

We have assumed that the dipole moment is nonzero and included only that term in the overlap region [see Eqs. (16)-(19)]. This means that our approximation will only be reliable when the field is predominantly dipole in the overlap region . To be safe, the approximation is valid if for . Since is order unity, . Thus, the complete conditions for the validity of our approximation are

 ϵ≪1,BℓB1≪1ϵℓ−1,ℓ≥2. (28)

The latter condition is violated only when higher stellar moments dominate the dipole by rather extreme amounts. If this is the case, then the general method still applies, but one should instead assume the dominant moment in Eqs. (16) and (19).

### ii.5 Power output

The electromagnetic luminosity (energy loss rate) of the star is Beskin (2010)

 L=2Ω∫ψo0I(ψ)dψ. (29)

Plugging in Eqs. (22) and (20) yields

 L=4π4775n2μ2Ω4=0.95×(√4πμ)2Ω4. (30)

Note that the magnetic moment is in Gaussian units. Previous studies for in flat spacetime have found the power to be Contopoulos et al. (1999); Gruzinov (2005); McKinney (2006); Timokhin (2006), which is consistent with the small- result up to numerical accuracy.

In curved spacetime, is the magnetic moment inferred from the field in the overlap region (far from both the star and the light cylinder). We therefore prove that, when compared at fixed , there are no GR corrections at leading order in . The corrections found at fixed in Refs. Ruiz et al. (2014); Pétri (2016) will therefore disappear at smaller values of .

In terms of the stellar dipole moment , Eq. (30) becomes

 L=0.95×(√4πB1R3⋆Δ1(C))2Ω4. (31)

This formula is rather insensitive to the details of the stellar magnetic field, depending only on its dipole component . From an observation of pulsar spindown, we can therefore infer the dipole field . Noting that the angular momentum of the star is given by and using the general relation , to leading order in we have

 √4πB1=Δ1(C)2π√10.95IMR4⋆P˙P, (32)

where and are the period and period derivative, respectively. We can view the factor of as the fractional correction due to general relativity. At a typical compactness, we have , so GR makes a roughly correction.

These results pertain to force-free fields, showing that electromagnetic losses decrease with increasing compactness (at fixed surface magnetic field strength). Interestingly, the opposite is true for vacuum fields (for the oblique rotator): the losses increase with increasing compactness Rezzolla and Ahmedov (2004). The same reversal occurs for losses from stellar pulsations Abdikamalov et al. (2009).

### ii.6 Polar Caps

In the classic dipole pulsar, current flows off the star only in regions near the poles called the “polar caps.” When higher moments are allowed, there can be multiple regions of current outflow and they can be located anywhere on the star. We will show, however, that to leading order in , there are always exactly two regions, connecting to the northern and southern outflows in the asymptotic region. We will continue to refer to these regions as northern and southern polar caps, although they need not be located at the poles.

The regions of current outflow (open field lines) are described by in the far-zone. Respectively, the bounding field lines and intersect the star at angles and determined by [see Eq. (20)]

 ψ⋆(θ0) =0, (33a) ψ⋆(θc) =B1R2⋆nΔ1ϵ. (33b)

Since Eqs. (33) are to be solved to leading order in , all solutions for are parametrically near a solution for . The leading Taylor approximation for about the nearby is

 ψ⋆(θc)≈(θc−θ0)mm!ψ(m)(θ0), (34)

where is the first nonvanishing derivative of at . Thus, the leading-order solution for is

 θc=θ0+⎛⎝m!B1R2⋆ψ(m)⋆(θ0)nΔ1ϵ⎞⎠1/m. (35)

If is even, then there are two solutions, corresponding to the ambiguity in the power of . This associates up to two values of with each root .

Among this potential panoply of solutions to Eqs. (33), only the values of closest to the poles actually connect to the far-zone outflow.2 In particular, the northern polar cap corresponds to for the smallest allowed and its associated , while the southern polar cap corresponds to for the largest allowed and its associated . To summarize, the northern/southern polar cap is found by finding all zeros of , computing associated values from Eq. (35), and choosing the pairs containing the smallest/largest values of .

Notice that the width of the polar cap region is set by the order of the first nonvanishing derivative at . When a polar cap lies at a pole or , the order is always at least , since all physical flux functions vanish quadratically on the axis. Generically, it will be exactly . On the other hand, if the polar cap is shifted, there is no reason to expect the first derivative to vanish, and generically we will have . Thus, apart from finely tuned situations,

 polar cap width∼{√ϵ% if on rotational axis,ϵif displaced. (36)

Note that while displaced polar caps are narrower, the total surface area is comparable since in both cases.

We have defined the polar cap as the angular region of current flow on the stellar surface. More generally, we could consider the angular region of the current flow at an arbitrary radius. The analysis of the section still applies using the flux function at that radius. As the radius increases and the dipole term starts to dominate, the polar cap will eventually move to the pole and assume the corresponding scaling, as indeed it must in order to match to the far-zone solution (Fig. 1). The transition region between the scalings can be handled by the methods discussed in App. C for the quadrudipole pulsar.

### ii.7 Near-zone Current

We can compute the charge-current from Eqs. (1), (6) and (7). To leading order in the near zone (i.e., fixing ), we find

 Jt =2(Ω−ΩZ)r(r−2M)[(r−3M)∂rψ+cotθ∂θψ], (37a) Jr =−∂θI(ψ)2πr2sinθ=−∂θI(ψ)2π√−g, (37b) Jθ =∂rI(ψ)2πr2sinθ=∂rI(ψ)2π√−g, (37c) Jϕ =0, (37d)

where a near-zone solution (14) should be used for . Note that the current is proportional to by the boundary condition (12) and in any case by the explicit formula (22). Eqs. (37) are then seen to be , noting that and is fixed in the near expansion. GR effects for the charge density were also discussed in Refs. Beskin (1990); Muslimov and Tsygan (1992).

The norm of the current is nonzero only at quadratic order , where we obtain

 J2 =−4(Ω−ΩZ)2r3(r−2M)[(r−3M)∂rψ+cotθ∂θψ]2 (38) +rr−2M[∂θI(ψ)2πr2sinθ]2+r2[∂rI(ψ)2πr2sinθ]2.

The current is only nonzero in two parametrically small regions whose footprints on the star are called the northern and southern polar caps (see Sec. II.6). Thus, in each region there is a small parameter to make further simplifications to the expressions (37)-(38) for the current.

## Iii Dipole Pulsar

Suppose that the field on the star is pure dipole,

 ψ⋆(θ)=B1R2⋆sin2θ. (39)

Then the near-zone solution is

 ψnear(r,θ)=μR>1(r)sin2θ, (40)

where and is given in Eq. (75a). The polar caps are found by solving Eqs. (33) using Eq. (39),

 sin2θ0=0,sin2θc=nΔ1ϵ, (41)

 θN∈(0,θ∗),θS∈(π−θ∗,π), (42)

where we defined

 θ∗=√nΔ1ϵ. (43)

The entire magnetosphere is reflection-symmetric, so we may focus on the northern cap for simplicity. The current norm on the surface may be obtained from Eq. (38) by plugging in Eq. (22) and Eq. (40). On the north cap, noting that and keeping to leading order in , the result is

 J2∣∣r=R⋆=(B1R⋆)2f(I,C,θθ∗)ϵ2, (44)

with

 Missing or unrecognized delimiter for \left (45)

The current norm is plotted for some representative values of and in Fig. 3. On the pole , we have

 f∣∣pole=16IC(2−IC1−C)=16α2(2−ΩZΩ)ΩZΩ, (46)

where for the second equality we have used Eq. (5). Since and for any reasonable star, we see that is always positive on the pole. This means that the current is always spacelike in a region near the pole. In light of the appearance of in Eq. (46), we can ascribe the spacelike current to the effects of frame-dragging Philippov et al. (2015b).

Ignoring the effects of gravity corresponds to taking . In this case, on the pole and it is straightforward to prove that everywhere else, as seen in the plot. This explains why work in flat spacetime missed the presence of spacelike current.

Concurrently with our work, Ref. Belyaev and Parfrey (2016) gave an independent analytical treatment of the current norm in the polar cap of a slowly rotating dipole pulsar in general relativity. Their method uses the monopole current (61b) along with the dipole flux function (40). By physical reasoning and comparison with numerical simulations, the authors recognized that this approach should only be reliable near the pole. Our result (22) for the dipole current quantifies the associated error, showing that their is valid to second order in . The error for is visible in their Fig. 1 (analogous to our Fig. 3), which lacks the upturn seen in Fig. 3. Our work also goes beyond that of Ref. Belyaev and Parfrey (2016) by systematizing the approximation method in terms of matched asymptotic expansions, considering higher multipole moments, and addressing the spindown luminosity.

We may also compare our results for the stellar current flow with previous force-free numerical simulations Philippov et al. (2015b); Belyaev and Parfrey (2016). These simulations were performed at relatively larger values of , so there is no a priori reason to expect agreement with our perturbative calculation. Nevertheless, there is full qualitative agreement of all features: spacelike current near the pole, transitioning to a timelike outflow, and finally, a timelike volume return current near the edge of the cap. This suggests that numerical simulations can be used to model the regime relevant to most pulsars, and conversely that our results may still be applicable for fast rotation.

We now consider a pulsar with both a dipole and a quadrupole moment,

 ψ⋆(θ) =R2⋆(B1sin2θ+B2cosθsin2θ) (47) =B1R2⋆(1+qcosθ)sin2θ, (48)

where we defined the ratio between the moments,

 q=B2B1. (49)

This configuration was studied previously in Ref. Barnard and Arons (1982), without constructing a self-consistent global model. Our approximation is valid when and . The near-zone flux function is

 ψnear=C1R>1(r)sin2θ+C2R>2(r)cosθsin2θ, (50)

where and . Expressions for and are given in Eqs. (75a) and (75b).

Since is equivalent to reflection about the equatorial plane, we may assume without loss of generality. The polar caps are found from Eqs. (33) using Eq. (39),

 sin2θ0(1+qcosθ0) =0, (51a) sin2θc(1+qcosθc) =nΔ1ϵ, (51b)

according to the method described in Sec. II.6.

For , the situation is similar to the dipolar case. The only solutions for are the poles and , and the associated polar caps are

 q<1:θN∈(0,θ+),θS∈(θ−,π), (52)

where

 θ+=√nϵ(1+q)Δ1,θ−=π−√nϵ(1−q)Δ1. (53)

On the other hand, for there is a new zero of the flux function located at

 θ0=π−arctan√q2−1, (54)

and the polar caps are

 q>1:θN∈(0,θ+),θS∈(θ−,θ0), (55)

where now

 θ−=θ0−q2(q2−1)3/2nΔ1ϵ. (56)

Thus, the northern polar cap is similar, but the southern cap is a narrow annulus ( instead of ) located away from the pole. The and results both break down for . This regime seems too highly tuned to be important in nature, but for completeness we discuss it in App. C. All three regimes can be handled simultaneously by working with exact solutions of Eqs. (51), which we also provide in App. C.

The current norm may be obtained from Eqs. (38) and (47). We give the full expression in Eq. (79). We can expand in in each cap of each regime of separately to produce analogs of Eq. (44), but the expressions are not transparent. Instead, we evaluate on the edge of the polar caps, where , and consider various limits.

First, for in the flat limit , we have

 J2N∣∣flat =0, (57) J2S∣∣flat =(2B1R⋆)2(q2+3Δ1)2(1−1q2)ϵ2q2. (58)

Here, is the value at and , while is evaluated at and . The northern current is null on the pole, while the southern current is always positive. This shows that the current can be spacelike even when gravity is ignored. Thus, while frame-dragging does produce spacelike current, in general it is not the only effect involved. The full expression for on the southern cap is given in Eq. (82).

When , the edges of the caps are the poles and

 J2N/S=(B1R⋆)2[16IC(2−IC1−C)(1±q2)]ϵ2. (59)

Here, is the value at and and corresponds to the plus sign, while is evaluated at and and corresponds to the minus sign. In particular, we see that the dipole formula (46) is corrected by a factor of of .

While we have focused on the fields at the stellar surface, the analysis of the angular extent of the current flow regions generalizes to arbitrary radius by replacing with the corresponding moment ratio computed at that radius. This defines a function that decreases roughly as . Near the star where , the width of the southern outflow scales as , while further away it scales as . The scaling in the transition regime is discussed in App. C, which also gives a master expression valid for all .

Fig. 4 illustrates our results for the quadrudipole pulsar with . We have also confirmed the overall field topology by running a numerical simulation in flat spacetime.

## V More General Models

We now discuss more general magnetic field configurations, including two that have arisen in recent research on neutron star magnetic field structure. We first consider the superposition of a dipole and octupole,

 ψ⋆(θ)=B1R2⋆[1+^q(1−5cos2θ)]sin2θ. (60)

It is straightforward to apply the method, but the formulas are complicated and we do not give details. We discuss the case , which was was found to be an attractor state in Hall-MHD studies of the magnetic field evolution in the neutron star crust Gourgouliatos and Cumming (2014). In this case, the polar caps are both at the poles, where there are regions of spacelike current. That is, the Hall attractor pulsar is qualitatively similar to the dipole pulsar.

On the other hand, if , then the polar caps are shifted away from the poles. By including a quadrupole moment as well, we could arrange to place the caps anywhere. Including even higher moments introduces a more complicated field structure in the closed zones, but does not change the overall topology. This is the limit of qualitatively new features accessible at small , so while it would be tempting to dramatize the power of the method by including higher moments, we see no reason beyond amusement to pursue the hexadecaquadrudipoctupole.

Another recent study suggests that the first several moments should be all roughly the same order of magnitude Geppert and Viganò (2014). In this case, the position of the polar caps is essentially random. However, the power output is still fixed by the dipole moment .

## Vi Discussion

In this paper, we developed an analytical technique to study the aligned force-free magnetosphere in the limit of slow stellar rotation with a generic magnetic field on the stellar surface, including all effects of general relativity. Previous numerical models were limited to pulsars with fast rotation (e.g., ) and dipolar surface fields. Our analytical results cover the regime of interest for most pulsars and extend to more general surface fields. Our only restriction on the magnetic field is that the dipole moment dominates far from the star (but still well inside the light cylinder). This assumption holds unless higher moments dominate the dipole by extreme amounts on the stellar surface.

We showed that the spindown luminosity is completely determined by the dipolar component of the surface magnetic field. This means that pulsar magnetic field strengths typically quoted from spindown measurements correspond only to the dipolar component. If pulsar magnetic fields contain many moments in equipartition Geppert and Viganò (2014), then the quoted value could differ from the true surface field strength by an order of magnitude or more. We also showed that ignoring GR effects leads to a roughly 60% error in the estimated surface field strength. This occurs because the radial dependence of the dipolar magnetic field near the star is modified from the perfect scaling of flat spacetime. The corrected formula for the dipolar component inferred from spindown is given in Eq. (32), where is one-half the value on the magnetic pole of a pure dipole. This formula will receive an additional order-unity correction when inclination is taken into account.

Our results show that aligned GR force-free magnetospheres generically have a region of spacelike current in the polar caps, which favors the existence of discharges Timokhin and Arons (2013). This is generally consistent with the first-principles kinetic simulation of Ref. Philippov et al. (2015b). One difference is that these authors also found a lower limit for the value of stellar compactness at which pair production can occur, whereas the force-free solution has spacelike current at any compactness. The discrepancy was attributed to lower values of the bulk current in the polar cap, a feature ultimately traced to the existence of a large vacuum region near the equatorial current sheet, where the volume return current flows in the force-free solution. If this gap is filled with plasma, then the polar cap current increases and reaches its force-free value Philippov et al. (2015b). The question of populating this region with pair plasma requires three-dimensional kinetic simulations that include pair creation due to gamma–gamma collisions, a computational challenge for the future.

We find that the presence of higher multipoles can significantly modify the geometry of the polar caps (defined as the angular region of current flow near the star). In particular, one or both polar caps can be shifted away from the rotational pole, occupying a thin annulus at some finite angle . The width of the annulus scales with the rotation rate as , in contrast to the scaling of an ordinary polar cap. Should these features persist in the case of an inclined rotator, they could potentially be observed if emission occurs close enough to the star for multipolar fields to be important. In ordinary pulsars, the radio emission is believed to be produced at distances , making the features unobservable. However, in millisecond pulsars the emission likely originates within a few stellar radii, where multipolar fields can potentially survive Kramer et al. (1998). Thus, we expect that observational signatures of shifted polar caps may be present for millisecond pulsars, but not ordinary pulsars.

One potential signature is double-peaked emission. If pulsed emission from an inclined rotator mimics the shape of its annular polar cap, then the profile would be strongly double-peaked as each edge of the hollow beam sweeps through the telescope. Although the details of the generation of pulsed emission are poorly understood, one might, in general, expect a more complicated beam shape when emission occurs near displaced polar caps, and hence more complicated beam shapes for millisecond pulsars. There is general consensus that millisecond pulsars do indeed display more complicated beam shapes than ordinary pulsars (at least marginally Kramer et al. (1998)), but this is difficult to quantify. The next generation of radio telescopes (e.g., SKA Karastergiou et al. (2015)) will significantly increase the sample of the observed millisecond pulsars and enable detailed studies of their emission geometry.

Another potential signature is a narrower pulse profile. Shifted polar caps are parametrically narrower ( instead of ) than caps at the poles. If only one of the two potential peaks is actually observed (for example, due to asymmetric synchrotron absorption of radio waves in the magnetosphere Beskin and Philippov (2012)), then the main signature of shifted polar caps is a modified beam width scaling with rotation rate. Tantalizingly, such a beam size scaling discrepancy has already been observed: millisecond pulsars are found to have systematically smaller beam widths compared to the usual scaling of normal radio pulsars Kramer et al. (1998).

A third potential signature is the phase-lag between gamma-ray and radio pulses. Recent studies suggest that current sheets in the magnetosphere are natural sources of gamma rays Bai and Spitkovsky (2010); Uzdensky and Spitkovsky (2014); Cerutti et al. (2016). Since the position of the current sheet is essentially unchanged, we expect gamma-ray light curves of multipolar magnetospheres to be identical to the dipolar case. However, the radio beam shape may be significantly modified, as discussed above. This could have important consequences for observations of radio to gamma-phase lags in millisecond pulsars Abdo et al. (2010).

Finally, shifted polar caps may have consequences for the pulsar wind by significantly enhancing the pair injection rate. On a shifted polar cap, the curvature radius of the magnetic field is generally much shorter compared to a dipolar field close to the pole, rendering curvature radiation more efficient Barnard and Arons (1982). We have shown that spacelike current generically appears, and in particular covers a large portion of the polar cap in the quadrudipole model that we study in detail (Fig. 4). These two facts favor efficient particle acceleration and emission of energetic curvature photons to produce secondary pairs.

## Acknowledgements

We thank Benoît Cerutti, Konstantinos Gourgouliatos and Anatoly Spitkovsky for fruitful discussions and Alexander Tchekhovskoy for providing the simulation setup. This research was supported by the NASA Earth and Space Science Fellowship Program (grant NNX15AT50H to A.P.), as well as NSF grants 1205550 to Harvard University and 1506027 to the University of Arizona. The simulations presented in this paper used computational resources supported by the PICSciE-OIT High Performance Computing Center and Visualization Laboratory.

## Appendix A Split monopole

The split monopole configuration can be solved analytically by perturbing the exact solution in the Schwarzschild spacetime Lyutikov (2011); Brennan et al. (2013). This is the analog of the procedure used by Blandford & Znajek Blandford and Znajek (1977) to derive the split monopole solution for a slowly-spinning black hole. Since the Kerr metric and the rotating star metric agree to linear order in , the only difference is the boundary condition at the compact object. In the black hole case, one finds , whereas here, we have . This yields

 ψ(r,θ) =ψo(1−|cosθ|), (61a) I(ψ) =±2πΩψ(2−ψψo)=±2πΩsin2θ, (61b)

where we take the sign (respectively, sign) for the northern (respectively, southern) hemisphere.

The last open field line lies on the equator,

 ψo=ψ(r,θ)∣∣θ=π/2=B0R2⋆. (62)

The power is

 L=8π3ψ2oΩ2. (63)

The four-current is given by

 Jt =2ψo(rα)2(Ω−ΩZ)|cosθ|, (64a) Jr =−2ψor2Ω|cosθ|, (64b) Jθ =0, (64c) Jϕ =0. (64d)

This agrees with the near-zone current expression Eq. (37) using Eq.(61), but Eqs. (64) are valid everywhere. The norm of the current is

 J2 =(2B0αR2⋆r2)2(2Ω−ΩZ)ΩZcosθ2 (65a) =B20R2⋆h(I,C,θ)ϵ2, (65b)

with

 h=4IC(2−IC1−C)cosθ2. (66)

We see that the current is spacelike everywhere for a reasonable star. When evaluated on the pole, this agrees with the dipole results in Eqs. (44) and (46), noting that is the value of the magnetic field on the pole.

## Appendix B Solutions of the near-zone equation

The solutions of Eq. (13) are presented in Beskin (2010); Gralla et al. (2016). The angular harmonics are proportional to the Gegenbauer polynomials ,

 _ℓ(θ)= _2F_1[ℓ2,-ℓ+12;12;cos^2θ]=-Γ(-ℓ2)Γ(ℓ+12)2πsin^2θ G^(3/2)_ℓ-1(cosθ) &ℓ odd, _2F_1[-ℓ2,ℓ+12;32;cos^2θ]cosθ=-(-1)^ℓ/2πΓ(ℓ2)4 Γ(ℓ+32)sin^2θ G^(3/2)_ℓ-1(cosθ) &ℓ even.

The first few harmonics are

 Θ1(θ) =sin2θ, (68) Θ2(θ) =cosθsin2θ, (69) Θ3(θ) =(1−5cos2θ)sin2θ, (70) Θ4(θ) =(1−73cos2θ)cosθsin2θ. (71)

We denote the radial harmonics regular at infinity by . These are given by

 Missing or unrecognized delimiter for \left (72)

where the are polynomials recursively defined for by

 P1(x) =x2+x2,P2(x)=4x4−x3−x26, (73a) Pℓ(x) =(2ℓ−1)[ℓ(ℓ−1)(2x−1)−1]xPℓ−1(x)−ℓ2(ℓ−2)x2Pℓ−2(x)(ℓ+1)(ℓ−1)2, (73b)

The normalization of the radial harmonics is chosen so that as , we have

 R>ℓ(r)=1rℓ+O(1rℓ+1). (74)

The first few harmonics are

 R>1(r) =−32r[3−4α2+α4+4logα(1−α2)3], (75a) R>2(r) =−103r2[17−9α2−9α4+α6+12(1+3α2)logα(1−α2)5], (75b) R>3(r) =−354r3[43+80α2−108α4−16α6+α8+24(1+8α2+6α4)logα(1−α2)7], (75c)

where was introduced in Eq. (2).

The exact solution of Eq. (51) is

 cosθ±=13q[1±i√322/31+3q2X+1∓i√324/3X−1], (76)

where we introduced

 X3 ≡Y+√Y2−4(1+3q2)3, (77a) Y ≡2+9(3nΔ1ϵ−2)q2. (77b)

For small , we can approximate this as

 θ+=√nϵ(1+q)Δ1,θ−=⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩π−√nϵ(1−q)Δ1q<1,π−√q−1+√2nϵΔ1+(q−1)2,q∼1,θ0−q2(q2−1)3/2nϵΔ1,q>1. (78)

The first and third lines correspond to Eqs. (53) and (56); now, we add an expression for the transition region. Equation (78) is convenient for seeing the scaling with in the various regimes, but for studying the problem globally it is generally easier to work with the single expression (76). We can view Eq. (76) as one convenient way to smoothly transition between the regimes by adding in higher-order terms in .

In terms of dimensionless coordinates and , the quadrudipole current norm is given everywhere by

 J2 Extra open brace or missing close brace (79)

where

 g Missing \left or extra \right (80a) Σ =Λ1+(Q−1)Λ2=Λ1(1+qcosθ),Q=1+Λ1Λ2qcosθ, (80b) A Extra open brace or missing close brace (80c) B =92C−3δ+3C(C−δ)Λ1−10(Q−1)[9−10δC+2(3−5δC)(C−δ)Λ1]. (80d)

Here, are dimensionless radial harmonics scaled to reduce to at the star,

 Λ1(δ)=−32C3[C(C+2δ)+2δ2log(1−Cδ)],Λ2(δ)=−203C2[4+(3C−4δ)Λ1(δ)],Λ1,2(1)=Δ1,2. (81)

For , when evaluated at (corresponding to ), the result is

 J2|θ=θ0 =(B1R⋆)2(2qϵ)21−C(1−1q2)2⎧⎨⎩1+(20Δ1Δ2C2)2[U2(1−C)(q2−1)−(Vq)2(1−IC1−C)2]⎫⎬⎭, (82a) U =43Δ1(Δ1+2)(1−C)−4,V=Δ21(1−C)C−4Δ1(1−C)2+2(2−3C). (82b)

### Footnotes

1. In perturbation theory, leading nontrivial order (or “leading order”) means the first order at which a given quantity is nonzero.
2. In the overlap region viewed from the near expansion, the field is dipolar but is parametrically small and hence, the last open field line is actually parametrically close to the pole. This means that the near-zone outflow is parametrically close to a field line that is asymptotically polar as in the near-zone, which corresponds to the solution to Eqs. (33) for which is closest to the pole. For the northern flow, , while for the southern flow, .

### References

1. P. Goldreich and W. H. Julian, ApJ 157, 869 (1969).
2. P. A. Sturrock, ApJ 164, 529 (1971).
3. I. Contopoulos, D. Kazanas,  and C. Fendt, ApJ 511, 351 (1999)arXiv:astro-ph/9903049 .
4. A. Gruzinov, Phys. Rev. Lett. 94, 021101 (2005)arXiv:astro-ph/0407279 .
5. J. C. McKinney, MNRAS 368, L30 (2006)arXiv:astro-ph/0601411 .
6. A. N. Timokhin, MNRAS 368, 1055 (2006)arXiv:astro-ph/0511817 .
7. A. Spitkovsky, ApJ 648, L51 (2006)arXiv:astro-ph/0603147 .
8. C. Kalapotharakos and I. Contopoulos, A&A 496, 495 (2009)arXiv:0811.2863 [astro-ph] .
9. J. Li, A. Spitkovsky,  and A. Tchekhovskoy, ApJ 746, 60 (2012)arXiv:1107.0979 [astro-ph.HE] .
10. C. Kalapotharakos, D. Kazanas, A. Harding, et al.ApJ 749, 2 (2012)arXiv:1108.2138 [astro-ph.SR] .
11. S. S. Komissarov, MNRAS 367, 19 (2006)arXiv:astro-ph/0510310 .
12. A. Tchekhovskoy, A. Spitkovsky,  and J. G. Li, MNRAS 435, L1 (2013)arXiv:1211.2803 [astro-ph.HE] .
13. L. Lehner, C. Palenzuela, S. L. Liebling, C. Thompson,  and C. Hanna, Phys. Rev. D 86, 104035 (2012)arXiv:1112.2622 [astro-ph.HE] .
14. C. Palenzuela, MNRAS 431, 1853 (2013)arXiv:1212.0130 [astro-ph.HE] .
15. M. Ruiz, V. Paschalidis,  and S. L. Shapiro, Phys. Rev. D 89, 084045 (2014)arXiv:1402.5412 [astro-ph.HE] .
16. A. Y. Chen and A. M. Beloborodov, ApJ 795, L22 (2014)arXiv:1406.7834 [astro-ph.HE] .
17. A. A. Philippov, A. Spitkovsky,  and B. Cerutti, ApJ 801, L19 (2015a)arXiv:1412.0673 [astro-ph.HE] .
18. S. Shibata, MNRAS 287, 262 (1997).
19. A. M. Beloborodov, ApJ 683, L41 (2008)arXiv:0710.0920 [astro-ph] .
20. A. N. Timokhin and J. Arons, MNRAS 429, 20 (2013)arXiv:1206.5819 [astro-ph.HE] .
21. A. A. Philippov, B. Cerutti, A. Tchekhovskoy,  and A. Spitkovsky, ApJ 815, L19 (2015b)arXiv:1510.01734 [astro-ph.HE] .
22. U. Geppert and D. Viganò, MNRAS 444, 3198 (2014)arXiv:1408.3833 [astro-ph.SR] .
23. K. N. Gourgouliatos, T. S. Wood,  and R. Hollerbach, Proc. Natl. Acad. Sci. USA 113, 3944 (2016)arXiv:1604.01399 [astro-ph.SR] .
24. M. Kramer, K. M. Xilouris, D. R. Lorimer, O. Doroshenko, A. Jessner, R. Wielebinski, A. Wolszczan,  and F. Camilo, ApJ 501, 270 (1998)arXiv:astro-ph/9801177 .
25. J. J. Barnard and J. Arons, ApJ 254, 713 (1982).
26. J. B. Hartle and K. S. Thorne, ApJ 153, 807 (1968).
27. H. A. Buchdahl, Phys. Rev. 116, 1027 (1959).
28. P. Haensel, A. Y. Potekhin,  and D. G. Yakovlev, eds., Neutron Stars 1 : Equation of State and Structure, Astrophysics and Space Science Library, Vol. 326 (2007).
29. S. E. Gralla and T. Jacobson, MNRAS 445, 2500 (2014)arXiv:1401.6159 [astro-ph.HE] .
30. A. Tchekhovskoy, R. Narayan,  and J. C. McKinney, ApJ 711, 50 (2010)arXiv:0911.2228 [astro-ph.HE] .
31. S. E. Gralla, A. Lupsasca,  and M. J. Rodriguez, Phys. Rev. D 92, 044053 (2015)arXiv:1504.02112 [gr-qc] .
32. C. F. Gammie, J. C. McKinney,  and G. Tóth, ApJ 589, 444 (2003)arXiv:astro-ph/0301509 .
33. A. Tchekhovskoy, J. C. McKinney,  and R. Narayan, MNRAS 379, 469 (2007)arXiv:0704.2608 [astro-ph] .
34. A. Tchekhovskoy, A. Philippov,  and A. Spitkovsky, MNRAS 457, 3384 (2016)arXiv:1503.01467 [astro-ph.HE] .
35. V. S. Beskin, MHD Flows in Compact Astrophysical Objects: Accretion, Winds and Jets, Astronomy and Astrophysics Library (Springer-Verlag, Berlin Heidelberg, 2010).
36. L. Rezzolla and B. J. Ahmedov, MNRAS 352, 1161 (2004)gr-qc/0406018 .
37. E. B. Abdikamalov, B. J. Ahmedov,  and J. C. Miller, MNRAS 395, 443 (2009)arXiv:0809.1191 [gr-qc] .
38. V. S. Beskin, Sov. Astron. Lett. 16, 286 (1990).
39. A. G. Muslimov and A. I. Tsygan, MNRAS 255, 61 (1992).
40. M. A. Belyaev and K. Parfrey, ArXiv e-prints  (2016), arXiv:1604.05670 [astro-ph.HE] .
41. K. N. Gourgouliatos and A. Cumming, MNRAS 438, 1618 (2014)arXiv:1311.7004 [astro-ph.SR] .
42. A. Karastergiou, S. Johnston, A. Karastergiou, S. Johnston, N. Andersson, R. Breton, P. Brook, C. Gwinn, N. Lewandowska, E. Keane, M. Kramer, J. P. Macquart, M. Serylak, R. Shannon, B. Stappers, J. van Leeuwen, J. Verbiest, P. Weltevrede,  and G. Wright, PoS(AASKA14)038  (2015), arXiv:1501.00126 [astro-ph.HE] .
43. V. S. Beskin and A. A. Philippov, MNRAS 425, 814 (2012)arXiv:1107.3775 [astro-ph.HE] .
44. X.-N. Bai and A. Spitkovsky, ApJ 715, 1282 (2010)arXiv:0910.5741 [astro-ph.HE] .
45. D. A. Uzdensky and A. Spitkovsky, ApJ 780, 3 (2014)arXiv:1210.3346 [astro-ph.HE] .
46. B. Cerutti, A. A. Philippov,  and A. Spitkovsky, MNRAS 457, 2401 (2016)arXiv:1511.01785 [astro-ph.HE] .
47. A. A. Abdo, M. Ackermann, M. Ajello, W. B. Atwood, M. Axelsson, L. Baldini, J. Ballet, G. Barbiellini, M. G. Baring, D. Bastieri,  and et al., ApJS 187, 460 (2010)