X-ray Lightcurves from Realistic Polar Cap Models: Inclined Pulsar Magnetospheres and Multipole Fields

TL;DR

This paper models realistic thermal X-ray lightcurves from pulsar polar caps considering complex magnetic field geometries, including quadrupole components, to better understand observed pulsar emission features.

Contribution

It introduces advanced polar cap temperature distributions incorporating multipole magnetic fields, improving upon simple models and explaining complex lightcurve features.

Findings

Realistic models show 5-10% differences in lightcurves for millisecond pulsars.

Quadrupole magnetic components create asymmetries and complex hotspot shapes.

Models can qualitatively reproduce features of PSR J0437-4715 lightcurve.

Abstract

Thermal X-ray emission from rotation-powered pulsars is believed to originate from localized "hotspots" on the stellar surface occurring where large-scale currents from the magnetosphere return to heat the atmosphere. Lightcurve modeling has primarily been limited to simple models, such as circular antipodal emitting regions with constant temperature. We calculate more realistic temperature distributions within the polar caps, taking advantage of recent advances in magnetospheric theory, and we consider their effect on the predicted lightcurves. The emitting regions are non-circular even for a pure dipole magnetic field, and the inclusion of an aligned magnetic quadrupole moment introduces a north-south asymmetry. As the aligned quadrupole moment is increased, one hotspot grows in size before becoming a thin ring surrounding the star. For the pure dipole case, moving to the more realistic model changes the lightcurves by $5-10\%$ for millisecond pulsars, helping to quantify the systematic uncertainty present in current dipolar models. Including the quadrupole gives considerable freedom in generating more complex lightcurves. We explore whether these simple dipole+quadrupole models can account for the qualitative features of the lightcurve of PSR J0437$-$4715.

Equations

$$\begin{array}{rl}d{s}^2& =-f(r)d{t}^2+f(r{)}^{-1}d{r}^2+\\ & r^2\left[d{\theta }^2+{\sin }^2\theta {\left(d\varphi -{\Omega }_Z(r)dt\right)}^2\right],r>R,\end{array}$$

$$\cos {\theta }^{\prime }$$

$$\tan {\varphi }^{\prime }$$

$${\psi }_{\star }({\theta }^{\prime })=R^2\sum _{\ell =1}^{\infty }{B}_{\ell }{\Theta }_{\ell }({\theta }^{\prime }),$$

$$\displaystyle\begin{split}J^{\hat{t}}&=\frac{\Omega-\Omega_{Z}}{\sqrt{f}r^{2}}\bigg{[}\partial_{\theta}\alpha\,\partial_{\theta}(\partial_{\phi}\beta)-\partial_{\theta}\beta\,\partial_{\theta}(\partial_{\phi}\alpha)-\partial_{\phi}\alpha\\ &\left(\frac{\partial_{\theta}(\sin{\theta}\partial_{\theta}\beta)}{\sin{\theta}}\right)+\partial_{\phi}\beta\left(f\partial_{r}(r^{2}\partial_{r}\alpha)+\frac{\partial_{\theta}(\sin{\theta}\partial_{\theta}\alpha)}{\sin{\theta}}\right)\!\bigg{]}\end{split}$$

$$J^{\hat{r}}$$

$$J^{\hat{\theta }}$$

$$J^{\hat{\varphi }}$$

$$\alpha (r,\theta ,\varphi )$$

$$\beta (r,\theta ,\varphi )$$

$$\displaystyle\Lambda(\alpha,\beta)=\begin{cases}\,\mp 2\Omega\bigg{[}J_{0}\left(2\arcsin{\sqrt{\alpha/\alpha_{0}}}\right)\cos{\zeta}\\ \quad\mp J_{1}\left(2\arcsin{\sqrt{\alpha/\alpha_{0}}}\right)\cos{\beta}\sin{\zeta}\bigg{]},&\alpha<\alpha_{0}\\ 0,&\alpha>\alpha_{0}\end{cases},$$

$${\alpha }_0$$

$$\mu$$

$$P_{\mathrm{in}}=(\bar{\gamma }-1)m_{\mathrm{e}}\chi \frac{|\vec{J}|}{e},$$

$$T={\left[\frac{m_e\chi (\bar{\gamma }-1)}{e\sigma }\right]}^{1/4}|\vec{J}{|}^{1/4}.$$

$$T(\theta ,\varphi )\equiv \{\begin{array}{ll}{T}_0|\vec{J}/J_0{|}^{1/4},& |\rho |<|\vec{J}|(\text{spacelike})\\ 0,& |\rho |>|\vec{J}|(\text{timelike})\end{array},$$

$$I(\theta ,\varphi ;ϵ,\Theta )=\frac{{ϵ}^3}{\exp [ϵ/T(\theta ,\varphi )]-1}\mathcal{B}(\Theta ).$$

$${\Theta }_1(\theta )$$

$$R_1^{>}(r)$$

$$\alpha =B_1{R}^2\left(\frac{R_1^{>}(r)}{R_1^{>}(R)}\right){\sin }^2{\theta }^{\prime },\beta ={\varphi }^{\prime }.$$

$$\begin{array}{rl}\alpha (r)& =B_1{R}^2\left(\frac{R_1^{>}(r)}{R_1^{>}(R)}\right){\sin }^2{\theta }^{\prime }+\\ & B_2{R}^2\left(\frac{R_2^{>}(r)}{R_2^{>}(R)}\right)\cos {\theta }^{\prime }{\sin }^2{\theta }^{\prime },\end{array}$$

$$R_1^{>}(r)$$

$$R_2^{>}(r)$$

Full text

pulsars: general, pulsars: individual (PSR J0437$-$4715), stars: neutron, X-rays: stars, magnetic fields

1 Introduction

Since their discovery in the 1960s, pulsars have offered a unique window onto the exotic physics of neutron stars (for a recent review, see Özel & Freire, 2016). Radio pulsars that are gravitationally bound to a companion in a binary orbit lead to precise measurements of the neutron-star masses and to global constraints on the properties of ultra-dense baryonic matter (Demorest et al., 2010; Antoniadis et al., 2013, 2016; Fonseca et al., 2016). Observations of pulsar glitches and of the rate of thermal cooling of young pulsars offer additional information about their moments of inertia (Link et al., 1999) and their internal compositions (Page et al., 2006), respectively. Spectroscopic analyses of the thermal emission from bursting and quiescent neutron stars have provided direct measurements of the radii of a handful of sources (Özel et al., 2009; Güver et al., 2010; Özel et al., 2010; Steiner et al., 2010; Guillot et al., 2013; Guillot & Rutledge, 2014; Heinke et al., 2014; Nättilä et al., 2016; Özel et al., 2016; Bogdanov et al., 2016). More recently, the detection of gravitational waves from the neutron-star merger event GW170817 placed constraints on the tidal deformability (Abbott et al., 2017; The LIGO Scientific Collaboration et al., 2018) and, consequently, on the radii of the two neutron stars (Raithel et al., 2018; De et al., 2018; Zhao & Lattimer, 2018; The LIGO Scientific Collaboration et al., 2018; Radice et al., 2018; Coughlin et al., 2018).

A complementary approach to measuring neutron-star radii involves modeling the lightcurves of the surface thermal emission from rotation-powered millisecond pulsars. This emission is thought to arise from localized “hotspots” on the stellar surface that come in and out of view as the star rotates, causing the observed brightness to oscillate periodically in time. The self-lensing of this emission in the gravitational field of a neutron star of mass $M$ and radius $R$ suppresses the amplitude of pulsations by factors that depend primarily on the compactness of the neutron star $GM/Rc^{2}$, where $G$ and $c$ are the gravitational constant and the speed of light, respectively (Pechenick et al., 1983). As a result, fitting the detailed properties of X-ray lightcurves observed from millisecond pulsars of known mass can lead, in principle, to a measurement of their radii (see, e.g., Zavlin et al., 1995; Zavlin & Pavlov, 1998; Bogdanov et al., 2007; Bogdanov, 2013; Psaltis et al., 2014; Özel & Freire, 2016; Watts et al., 2016; Bogdanov et al., 2019).

The Neutron-star Interior Composition ExploreR (NICER) is now measuring X-ray lightcurves from millisecond pulsars with sensitivity that is poised to make such measurements possible (Gendreau et al., 2016). However, to exploit the full capabilities of these observations, a new generation of theoretical modeling is required to extract physical parameters from these lightcurves and thereby advance efforts to constrain the neutron-star equation of state.

In the canonical pulsar model, the neutron star surface is heated by relativistic electrons and positrons arising from the pair production that supports the global magnetospheric circuit (Arons, 1981, 1983). The current flows along “open” field lines (i.e., those field lines that cross the light cylinder), whose intersection with the stellar surface defines the polar caps.111Here, the light cylinder $R_{\rm L}\equiv c/\Omega$ is defined as the distance at which an object co-rotating with a star with angular spin frequency $\Omega$ would move at the speed of light. Note that we will retain the term “polar caps” even when these regions are neither polar nor cap-shaped. For a pure dipole field configuration, the polar caps are circles centered on the magnetic poles (at least at slow rotation), but more generally a variety of shapes are allowed. Furthermore, the bombardment (and thus the resulting emission) occurs over only a sub-region of the polar caps determined by the detailed properties of the current (Timokhin & Arons, 2013). In this framework, the configuration and multipolar structure of the magnetic field of the neutron star determines the properties of the emitting regions.

Past efforts to model the thermal surface emission from pulsars have usually used circular hotspots that coincide with the magnetic poles under the implicit assumption of a purely dipolar magnetic field (see, e.g., Zavlin & Pavlov, 1998). This canonical choice was driven by the relative (numerical) simplicity of simulating the dipolar configuration as well as the lack of an obvious, physically motivated alternative. However, there is little reason to believe that pulsar fields are purely dipolar. Indeed it is natural to suspect that the complexity of pulsar emission, and its variation across sources (Lorimer, 2008; Rankin et al., 2017), are driven by similar complexity and variation in the pulsar magnetic field. Furthermore, for X-ray pulse profile modeling itself, the lightcurves and spectra observed from the prototypical source PSR J0437$-$4715 have required multiple concentric regions of emission with different temperatures, as well as two hotspots of emission that are not antipodal (Bogdanov, 2013). On the theoretical side, recent force-free and particle-in-cell (PIC) simulations (Spitkovsky, 2006; Bai & Spitkovsky, 2010; Chen & Beloborodov, 2014; Philippov et al., 2015a, b; Brambilla et al., 2018) and analytic calculations (Gralla et al., 2016; Belyaev & Parfrey, 2016; Gralla et al., 2017) have revealed a similar degree of complexity in polar cap properties.

In this paper, we model the shapes and thermal properties of polar caps on neutron stars with dipolar and “quadrudipolar” (i.e., dipole+quadrupole) magnetic fields using earlier semi-analytic general-relativistic calculations of force-free currents in their magnetospheres (Gralla et al., 2017). We combine these calculations with basic models of the resulting surface emission and with ray-tracing calculations in the external spacetimes of the neutron stars (see, e.g., Bauböck et al., 2012, 2013; Psaltis & Özel, 2014) to predict the expected lightcurves from these polar-cap models. Our goal is to explore new pulse morphologies and properties that result from more realistic polar-cap models and understand the magnitude of potential biases that may be introduced by the simplified configurations commonly used.

We find that realistic polar-cap models naturally give rise to regions of surface emission that have different temperatures and to lightcurves with primary and secondary peaks that are not offset by 180 degrees, as also inferred from observations. Furthermore, the realistic polar-cap models introduce complexity in the resulting lightcurves at the $5-10$% level. These effects will mask and alter similar complexity in the lightcurves introduced by rotational Doppler effects that depend primarily on the neutron-star radii (Psaltis et al., 2014). As a result, not accounting for the thermal structure of the pulsar polar caps can lead to imprecise model selection and may bias the radii measurements obtained with this approach.

2 Method

In this section, we introduce the various ingredients of our method, beginning with a set of neutron star parameters and arriving at a simulated lightcurve. Throughout the paper, we consider a neutron star of mass $M$ and radius $R$, rotating with angular frequency $\Omega$. We denote polar coordinates around the spin axis of the star by the colatitude $\theta$ and azimuth $\phi$, and the inclination of the observer’s line of sight with respect to the spin axis by $\theta_{o}$. We consider the case where the magnetic field of the neutron star is axisymmetric with respect to a “magnetic axis” that co-rotates with the star at an inclination angle $\zeta$ relative to the axis of rotation. We also set $G=c=1$, unless explicitly stated otherwise.

Our procedure is composed of three steps:

From a choice of stellar parameters and magnetic field configuration, we compute the four-current density $J^{\mu}(\theta,\phi)$ near the star using the model of a slowly rotating, perfectly conducting star with a force-free magnetosphere (Gralla et al., 2017). 2. 2.

Given the current $J^{\mu}(\theta,\phi)$ near the star, we determine the surface temperature $T(\theta,\phi)$ via different phenomenological prescriptions for Joule heating. From this temperature, we then determine the specific intensity of the emerging radiation using a blackbody spectrum and a simple beaming function which we define later on. 3. 3.

Finally, we simulate the predicted lightcurve using general-relativistic ray-tracing of null geodesics from the stellar surface to an observer at infinity (Bauböck et al., 2012, 2013; Psaltis & Özel, 2014).

Steps 1 and 3 rest on well-established underlying physics and contribute little systematic uncertainty to the results. By contrast, the phenomenological prescriptions in step 2 are motivated but by no means definitive. We therefore pay careful attention to the dependence of the lightcurves on the choice of prescription. To compare with observations would require a further step 4 of taking into account foregrounds and instrumental response, which is beyond the scope of the current work.

2.1 From Magnetization to Current

Our first task, given the parameters of the star and its magnetic field configuration, is to determine the current flow at the surface. This flow is sensitive to the global structure of the magnetospheric circuit (extending past the light cylinder), which can in general be determined only with expensive numerical simulations. However, in the slow-rotation limit, there is a decoupling of near and far regions that enables the full general-relativistic, non-dipolar problem to be solved analytically with input from a single simulation of a dipole pulsar in flat spacetime (Gralla et al., 2017).

We follow Gralla et al. and keep terms to leading order in the rotation parameter $\Omega R/c$. For NICER’s main target, and targets of similar interest, this parameter is estimated to be $\Omega R/c\sim 0.05$. By keeping to leading order in rotation one can therefore expect relative corrections of $(0.05)^{2}$ from higher-order terms, which can safely be ignored. We briefly review this approach and present the formulae we will need for lightcurve modeling.

To leading order in rotation, the spacetime geometry outside of a general-relativistic star is given by the linearized Kerr metric,

[TABLE]

where $f(r)=1-2M/r$ is the square of the redshift factor, $\Omega_{Z}=2I\Omega/r^{3}$ is the frame-dragging frequency, and $I$ is the moment of inertia (angular momentum over angular velocity). Note that the surface of the star is spherically symmetric at this order, although the spacetime is only axisymmetric.

Assuming a perfectly conducting star, the force-free magnetosphere is completely determined by specifying the radial component of the magnetic field on the stellar surface, $B^{r}_{\star}(\theta,\phi,t)$. Let primed angles $(\theta^{\prime},\phi^{\prime})$ denote polar coordinates with respect to the magnetic axis, related to the unprimed coordinates by a rotation,222We work at an instant of time and make the arbitrary choice to align the angular velocity with the $z$-axis.

[TABLE]

Since the field is axisymmetric, we can express $B^{r}_{\star}$ in terms of the magnetic flux in the co-rotating frame, $\psi_{\star}(\theta^{\prime})=\int B^{r}_{\star}\sqrt{-g}d\theta^{\prime}$. We may further decompose this flux into a set of multipole moments $B_{\ell}$,

[TABLE]

where $\Theta_{\ell}$ are the eigenfunctions for magnetic flux (Gralla et al., 2016). To leading order in $\Omega R$, the magnetosphere is completely determined by the choice of parameters $\{M,R,I,\Omega,\zeta,B_{\ell}\}$ defined above.

Using this framework, Gralla et al. (2017) have shown that the four-current density $J^{\mu}$ near the star has the following orthonormal frame components:

[TABLE]

where all the terms are to be evaluated at $r=R$. The scalars $\alpha$ and $\beta$ are “Euler potentials” for the magnetic field near the star, given explicitly by

[TABLE]

where $R_{\ell}^{>}(r)$ are radial eigenfunctions (Gralla et al. (2016) and Eq. 18 below). Finally, $\Lambda(\alpha,\beta)$ is a conserved quantity on field lines related to the current on each line. Gralla et al. (2017) provide an analytic expression for $\Lambda(\alpha,\beta)$, determined by fitting numerical data, in terms of Bessel functions $J_{n}$,333$\mp$ refers to the north/south hemisphere. The internal $\mp$ appearing in front of $J_{1}$ was missing from the original formula in Gralla et al. (2017), and has been corrected here.

[TABLE]

with the boundary of the polar cap delineated by

[TABLE]

Here $\mu$ is the magnitude of the dipole moment at the light cylinder, which Eq. (9) relates to the dipole moment $B_{1}$ at the star by compactness-dependent factors given below (see Eq. 15). The polar caps are defined by $0<\alpha<\alpha_{0}$. Their physical size and shape follows from the magnetic field, the magnetic inclination, the rotation rate, the moment of inertia, and the compactness via Eq. (5) and Eq. (8). The area of the polar cap scales as $(\Omega R)^{2}$ when the other parameters are held fixed. In this paper, we focus on the effects of magnetic field and inclination, setting the compactness and moment of inertia to fiducial values given in § 2.3.

These expressions allow the current $J^{\mu}$ on the star to be computed analytically from the parameters $\{M,R,I,\Omega,\zeta,B_{\ell}\}$.

2.2 From Current to Emerging Radiation

The mechanism most commonly invoked to explain the thermal surface emission from millisecond pulsars is the bombardment of the stellar atmosphere by relativistic electrons and positrons, produced in so-called “gaps” (regions of non-zero $\mathbf{E}\cdot\mathbf{B}$) just outside the star and accelerated by the fields (Arons, 1981, 1983). While current flows over the entirety of the polar cap (by definition), pair production (leading to bombardment) occurs mainly in regions where the four-current in the corresponding force-free solution is spacelike ($J_{\mu}J^{\mu}>0$), which requires the presence of both signs of charge (Timokhin & Arons, 2013; Philippov et al., 2015a). This spacelike region defines the shape of the hotspot in our model.444Pair production is also expected to occur in the thin current sheet on the last open field lines, and possibly in nearby regions of volume return current (Timokhin & Arons, 2013; Philippov & Spitkovsky, 2018). This could give rise to further X-ray emission from near the edges of the polar cap, an effect we neglect at present.

The detailed plasma physics governing the deposition of energy by magnetospheric currents in the surface layers of the neutron star is explored in Bauböck et al. (2019). As electrons and positrons bombard the surface of the star, their kinetic energy is deposited onto a thin layer of the stellar atmosphere. Assuming that a fixed fraction $\chi$ of the total current $|\vec{J}|$ is carried by relativistic particles with typical Lorentz factor $\bar{\gamma}$ traveling towards the stellar surface, we write the rate at which this energy is deposited in the surface layers as

[TABLE]

where $e$ and $m_{\rm e}$ are the electron charge and mass, respectively. The rate at which this energy is deposited must equal the power radiated, which can be written in terms of the effective temperature using the Stefan-Boltzmann law as $P_{\rm{out}}=\sigma T^{4}$. Setting $P_{\rm{in}}=P_{\rm{out}}$ gives

[TABLE]

Using this scaling and neglecting, for simplicity, the details of the atmospheric structure leads to the simple prescription for the effective temperature on the polar cap

[TABLE]

where $|\vec{J}|=\sqrt{(J^{\hat{r}})^{2}+(J^{\hat{\theta}})^{2}+(J^{\hat{\phi}})^{2}}$ is the three-current density, $\rho=J^{\hat{t}}$ is the charge density, $J_{0}$ is the maximum value of $|\vec{J}|$ for a neutron star with an aligned dipolar magnetic field, and $T_{0}$ is the corresponding maximum effective temperature for the aligned dipole.

In order to convert the effective temperature to a model for the emerging specific intensity of the radiation field, we need, in principle, a detailed calculation of the thermal structure of the neutron-star atmosphere (as done, e.g., in Bauböck et al., 2019). This is beyond the scope of this paper. Instead, we will use a simplified approximation of blackbody emission at the effective temperature, with beaming that is potentially non-isotropic,

[TABLE]

In this last equation, we denote the photon energy by $\epsilon$, the angle to surface normal by $\Theta$, and the beaming function by $\mathcal{B}(\Theta)$. We will mainly consider the case of isotropic emission ($\mathcal{B}=1$) except in §5, where we introduce a simple model for limb-darkening, $\mathcal{B}=1+h\cos^{2}\Theta$, with $h$ a free parameter. Equation (13) represents a compromise between realism and simplicity, given current understanding of the underlying micro-physical processes. Ultimately we expect that a more accurate prescription will be determined by simulations; in the meantime, we estimate the importance of the precise prescription by comparing to other reasonable models (§3 below).

2.3 From Emerging Radiation to Lightcurves

Using equation (13) as a boundary condition, we then use the numerical algorithm of Psaltis & Özel (2014) to trace the light rays along the null geodesics of the spacetime and map the surface of the simulated compact object onto the image plane of an observer at infinity. We keep only terms that are up to first order in the neutron-star spin frequency, treat the star as spherical (i.e., neglecting any mass quadrupole terms; see, e.g., Bauböck et al., 2013), and set the moment of inertia to $(2/5)MR^{2}$ for simplicity.

The simulations presented in this paper have eight parameters: the mass $M$, radius $R$, and spin frequency $f_{\rm s}$ of the star (or equivalently, the angular frequency $\Omega=2\pi f_{\rm s}$); the inclination angle $\theta_{o}$ between the spin axis and the line of sight of the observer; the inclination angle $\zeta$ of the magnetic axis from the spin axis; the beaming factor $h$; the temperature constant $T_{0}$, which relates current density to temperature; and the magnetic quadrupole-to-dipole ratio $q$. When comparing to observations of a known millisecond pulsar, the spin frequency is known a priori.

In the following sections we present analysis of lightcurves from pulsars with non-negligible dipolar and quadrupolar components of the magnetic fields. The procedure to include higher-order terms is straightforward in principle. We take as our default choice of stellar parameters $M=1.5M_{\odot},R=11$ km, and $f_{\rm s}=300~{}$Hz, which corresponds to a dimensionless compactness $M/R\simeq 0.2$ and surface velocity $\Omega R\simeq 0.07c$. Unless stated otherwise, all plots are shown for these fiducial values.

3 Pulsars with Pure Dipolar

Magnetic Fields

The simplest non-trivial magnetic field configuration is a pure dipole, where only the dipole moment $B_{1}$ in Eq. (3) is non-zero. The eigenfunctions for magnetic flux reduce to555Note that the expression for $R^{>}_{1}$ present in Gralla et al. (2017) contains a typographical error.

[TABLE]

where we remind the reader that $f=(1-2M/r)$. Equation (5) for the Euler potentials becomes

[TABLE]

Substituting equation (16) into equation (4) gives the current flow on the star, from which equation (12) provides the temperature distribution. The polar caps are antipodally symmetric.

Figure 1 shows the shape and temperature distribution on the polar caps for six different magnetic inclination angles. For our canonical parameters, the half-opening angle of the polar cap lies in the range $14-16\degree$ depending on the inclination. For small inclinations, the emitting region (i.e., the region with spacelike currents) is smaller than the polar cap but centered on the magnetic axis. As the inclination increases, the emitting region becomes larger and its center is displaced from the magnetic axis (i.e., the center of the cap) until, at a critical angle ($\zeta\simeq 70\degree$ for the fiducial parameters), a second lobe forms.

In the interest of simplicity and clarity of comparison to previous dipole models, we calculate lightcurves from neutron stars with dipolar magnetic fields by considering only isotropic surface emission, i.e., setting the beaming factor $h=0$ so that $\mathcal{B}=1$. In order to understand our results below, we first consider the lightcurve due to a single polar cap; the full lightcurve is then just the sum of the contribution from each polar cap separately. Figure 2 shows the lightcurve from a single polar cap on a neutron star with magnetic inclination $\zeta=75\degree$, spinning at three different frequencies, and observed at an inclination angle of $\theta_{o}=90\degree$. As has been shown earlier, the combination of Doppler and time-delay effects introduces a skewness to the lightcurve, making the rising (left) edge of the curve steeper and the waning (right) edge shallower (Braje et al., 2000).

When two of these skewed lightcurves are added together, $180\degree$ out of phase, to generate the lightcurve of a neutron star with two polar caps in the pure dipole field configuration, they give rise to the asymmetric troughs shown in Figure 3. The most striking features are the increased complexity in the lightcurves arising from such simple geometric polar caps and the fact that their time-reversal symmetry is broken.

In order to explore the dependence of the simulated lightcurves on the shape and temperature profile across the polar caps, we compare the results from three different prescriptions:

1. (A)

constant temperature over the circular polar cap (this is the geometry that has been employed in most previous analyses); 2. (B)

constant temperature over just the spacelike-current region; 3. (C)

the phenomenological model described in §2, with a temperature $T\propto|\vec{J}|^{1/4}$ within the spacelike-current region (Eq. 12).

Figure 4 compares the normalized lightcurves for two neutron stars, one with moderate spin (300 Hz) and one with high spin (600 Hz), calculated with these three different phenomenological prescriptions. The detailed shape and temperature profiles of the polar caps introduce complexity to the predicted lightcurves at amplitudes that increase with the neutron-star spin frequency and can reach levels as large as $\sim 5-10$%. A large part of this complexity is caused by the range of temperatures that are present on the more realistic polar caps. As the star rotates, the relative projection of different parts of the polar cap that have different temperatures evolves with time, causing the observed spectrum to also change with time. The net result is a pulse amplitude from realistic polar caps that is larger than that of a polar cap with constant temperature (cf. Bauböck et al., 2015).

This spectral evolution with pulse phase also results in a pulse amplitude that depends on photon energy. Figure 5 shows the RMS fractional pulse amplitude666The RMS fractional amplitude is defined as the square root of the time-average of the squared deviation from the mean flux. as a function of photon energy for the configurations shown in Figure 4. The weak dependence on photon energy seen in the pulse amplitudes of configurations with constant temperatures across the polar caps is caused entirely by rotational Doppler effects (Psaltis & Özel, 2014). The larger amplitudes and stronger dependence on photon energy seen in the simulations with realistic polar caps is a direct outcome of the presence of regions with different temperatures across each polar cap.

The increased amplitudes and strong photon energy dependence of the pulse profiles generated with realistic polar caps competes with and may mask the consequences of the rotational Doppler effects, which are qualitatively similar. If we were to attempt to fit the pulse profiles shown in Figures 4 and 5 using models with circular, constant temperature polar caps, we would be forced to increase artificially the rotational effects in order to account for the large amplitudes. For a neutron star of known spin frequency and observer orientation, the only way to achieve the latter is to increase artificially the stellar radius. As a result, modeling pulse profiles with circular, constant temperature emitting regions introduces the risk of biasing the inferences of the neutron-star radii.

4 Pulsars with Quadrudipolar Magnetic Fields

We now consider the effect on the simulated pulse profiles of adding a quadrupole moment to the magnetic field. This amounts to including the next term in the multipole expansion of the Euler potential,

[TABLE]

where the radial functions are given by

[TABLE]

We refer to this family of dipole-plus-quadrupole fields as “quadrudipoles” and parameterize them by the ratio $q\equiv B_{2}/B_{1}$, with $q=0$ corresponding to the pure dipole.

The addition of a magnetic quadrupole moment changes both the current density and the size and geometry of the polar caps, as shown in Figure 6. For $q<1$, the effect of the magnetic quadrupole moment is to spread out the field lines in the southern hemisphere and contract those in the north.777We take north to point in the direction of the spin vector according to the right-hand rule. As a result, the northern cap shrinks and grows hotter with increasing quadrupole, while the southern cap expands and cools. When $q>1$, the southern cap is no longer circular but becomes an annular strip encircling the star. When the magnetic field is substantially inclined from the spin axis, a large temperature gradient emerges along this ring (see bottom right insert in Figure 6), since only part of this region has spacelike four-current.

An important consequence of the quadrudipole model is that the northern and southern emitting regions have different effective temperatures. We can quantify this difference by calculating the ratio of the average temperatures of the two emitting regions as a function of the magnetic quadrupole moment $q$ and the magnetic inclination $\zeta$. Figure 7 shows that this temperature ratio depends primarily on $q$, initially increasing with $q$ for $q\leq 1$. This is a direct result of the fact that the total current flow onto each polar cap is the same. As the field lines spread out with increasing $q$, the surface area of the southern cap grows and the current density (which determines the temperature) decreases. A quadrudipole field with $q=1$ marks the critical point at which the southern emitting region transitions from being a spot to a ring. Beyond this point, the area of the southern region decreases with increasing $q$, causing its temperature to rise.

The quadrudipole configuration clearly breaks antipodal symmetry because the emitting regions in general have dramatically different shapes. The one symmetry that is maintained is reflection symmetry across the $\Omega\!-\!\mu$ plane, i.e. the plane determined by the spin ($\vec{\Omega}$) and magnetic dipole ($\vec{\mu}$) vectors. This reflection symmetry guarantees that the phases of greatest flux emitted from each hotspot are still $180\degree$ apart. However, when the full ray-tracing calculation is carried out — including Doppler shift, time-delays, and lensing — the phases of peak flux observed are no longer evenly spaced. The quadrudipole configuration makes possible additional features in the lightcurve that are not seen with the antipodal configuration of the pure dipole.

Figure 8 shows the simulated lightcurve from a pulsar with a quadrudipole magnetic field ($q=1.5$) for a particular choice of parameters. The interesting feature of this pulse profile is the appearance of an ‘interpulse’, which is offset in phase from the main peak by less than $180\degree$. This interpulse, which we might have naively identified with an off-axis hotspot, in fact arises because the emission regions from the northern spot and the southern ring are not antipodally symmetric, and their contributions to the pulse profile are affected differently when subject to the Doppler shift and time delays. Overall, we find temperature ratios as large as $\sim 2.0$ between the two emitting regions (Figure 7) and pulse profiles with offset interpulses (Figure 8), potentially providing a natural explanation to the multi-temperature blackbody models and non-antipodal emission geometries that are required to fit the observed spectra of millisecond pulsars (see, e.g., Bogdanov, 2013).

5 Application to X-ray Observation

of PSR J0437$-$4715

PSR J0437$-$4715 is the nearest known rotation-powered millisecond pulsar and the best candidate for measuring neutron-star properties via pulse profile modeling (Arzoumanian et al., 2014). It has a spin frequency $f_{\rm s}=174\rm{Hz}$ and orbits a white dwarf companion with an orbital period of 5.7 days. Pulsar timing measurements yield a pulsar mass of $M=1.44\pm 0.07M_{\odot}$ (Reardon et al., 2016).

The X-ray properties of this pulsar have posed significant challenges to modeling (Bogdanov, 2013). As discussed earlier, fitting the X-ray spectra of the surface emission requires multi-temperature blackbody models. The lightcurve itself exhibits an interpulse that is not evenly-spaced between the main pulses, but is rather shifted forward by roughly 20 degrees. In order to model this lightcurve, Bogdanov (2013) employed an “offset dipole” model, in which the magnetic field is dipolar but the dipole is not located at the geometric center of the star. This gives rise to emitting regions that are not antipodal. Nevertheless, the shapes of the hotspots are assumed to be circular. Employing this model to fit the observations requires four parameters in addition to the ones required for the model we explore here: angular offsets $\Delta\theta$ and $\Delta\phi$ that specify the location of one of the hotspots with respect to the other, the ratio of the two hotspot radii, and, independently, the ratio of their temperatures.

The quadrudipole model offers the potential of reproducing the key qualitative features of PSR J0437$-$4715 in a more satisfactory way by introducing the single parameter $q$ (the ratio of the aligned magnetic quadrupole to dipole moments) while retaining axisymmetry. Figure 9 shows a model lightcurve that is suggestive of a potential fit to PSR J0437$-$4715. The underlying magnetic field structure is quadrudipolar with $q=0.5$.

For the purposes of this plot, we have assumed a magnetic inclination of $\zeta=30\degree$ and an observer inclination of $\theta_{o}=42\degree$, which is the value inferred for this pulsar from the inclination of its binary orbit (Reardon et al., 2016). In order to generate a large pulse amplitude, we have also used non-isotropic beaming, as expected for the atmospheres of pulsars that are bombarded by magnetospheric charges (Bauböck et al., 2019), and chose $h=0.4$. In the simulated pulse profile, the main peak corresponds to the northern emitting region, which is seen almost head-on at its peak intensity. The interpulse is due to the southern emitting region, seen at a glancing angle.

This quadrudipole model offers several advantages over the offset dipole. First, we only need to add a single new parameter to the classic dipole model, without breaking axisymmetry in the magnetic field. Second, a natural consequence of the quadrudipole model is that the hotspot with the smaller surface area will also have the higher temperature. This is consistent with the relationship found in the multi-temperature blackbody fits of the data from PSR J0437$-$4715. Finally, the appearance of the offset interpulse is robust enough that we are also able to set the observer angle to $\theta_{o}=42\degree$, as is inferred for this source.

This quadrudipole model offers several advantages over the offset dipole. First, we only need to add a single new parameter to the classic dipole model, without breaking axisymmetry in the magnetic field. Second, a natural consequence of the quadrudipole model is that the hotspot with the smaller surface area will also have the higher temperature. This is consistent with the relationship found in the multi-temperature blackbody fits of the data from PSR J0437$-$4715. Finally, the appearance of the offset interpulse is robust enough that we are also able to set the observer angle to $\theta_{o}=42\degree$, as is inferred for this source.

Performing a detailed fit of our model to the data would require taking into account the effects of interstellar absorption, the response of the detector, as well as the pulsed magnetospheric emission that contaminates the measurement of the surface emission (Guillot et al., 2016). We expect to report on these issues in the future.

6 Conclusions

We conclude by summarizing our key results and discussing some limitations of the present work as well as promising future directions. We have shown that the shape and temperature profile across the pulsar hotspots plays an important role in lightcurve modeling that cannot be ignored, and that significant complexity in lightcurve features is possible simply with the addition of a quadrudipolar magnetic field.

This work is based on an analytic model of the pulsar magnetosphere that is guided by state-of-the-art numerical simulations. We can expect refinements of the magnetospheric structure and currents as numerical simulations improve. For example, volume return currents, which we ignore in our present model, may also play a role in atmospheric heating (Timokhin & Arons, 2013; Philippov & Spitkovsky, 2018). However, we expect that the basic qualitative results presented here will remain unchanged and that this approach will aid in measuring neutron star parameters and ultimately help tighten constraints on the neutron-star equation of state.

Acknowledgments

We thank Sasha Philippov and Deepto Chakrabarty for helpful discussions and valuable comments on the manuscript. We gratefully acknowledge support from NASA grant NNX16AC56G and NSF grant PHY-1752809.