PSR J0030+0451 Mass and Radius from NICER Data and Implications for the Properties of Neutron Star Matter

Here, we briefly summarize the NICER X-ray Timing Instrument (XTI) data set that was used for the parameter estimation analyses presented here. More detailed descriptions of the observations, data reduction, and event folding procedures, and a thorough examination of the properties of the final event data can be found in Bogdanov et al. (2019a).

The NICER XTI data on PSR J0030+0451 were collected in a set of short observations (typically lasting a few hundred to ∼2000 s) over the period from 2017 July 24 to 2018 December 9. Due to the exceptional rotational stability of PSR J0030+0451 and the superb absolute timing accuracy of NICER, neither the large time span nor the discontinuities in the exposure had any adverse effects, such as smearing due to long-term pulse-phase drifts, on the pulse profile.

For this analysis, we used products from the calibration database (CALDB) version 20181105 and gain solution version optmv7. In addition to the standard pipeline processing event filtering that all NICER data undergo, the PSR J0030+0451 data set was subjected to screening criteria tailored to the inference analyses discussed later.

We included only data accumulated when the source subtended an angle greater than 80^◦ from the Sun, in order to minimize the contamination from optical loading in the lowest-energy detector channels. All events from the XTI detector with DET_ID 34 were excised because it frequently exhibits significantly elevated count rates relative to the other detectors (analogous to a "hot pixel" in charge-coupled device detectors).

In the final filtering stage, we constructed a time series with 16 s time bins using the event data in channels 25–800 and then discarded all the events in those time bins that had an average count rate greater than 3 counts s⁻¹. The long-term average count rate during the NICER observations of PSR J0030+0451 was ∼0.7 counts s⁻¹. This cut therefore excluded the event data in all the 16 s time bins that had count rates more than four times the long-term average count rate. The purpose of this cut was to improve the signal-to-noise ratio of the measured thermal emission from the pulsar by removing data contaminated by flares that were not filtered out by the standard pipeline processing.

The thermal emission pulse profile that was used to estimate the mass and radius of PSR J0030+0451 was constructed by averaging the data from the roughly 400 million rotation periods observed by NICER during more than 18 months of observations. We compared the profile measured during the first third of this observing time with the profile measured during the second third, and also compared the profile measured during the first half of this observing time with the profile measured during the second half. These comparisons showed no statistically significant differences between the compared profiles. We expect that any shorter-term flux variations were averaged out by using such a long data set.

To generate the final X-ray pulse profile as a function of pulsar rotational phase and energy, we used the "photons" plugin for the Tempo2 pulsar timing package (Hobbs et al. 2006) and the best publicly available radio timing solution, obtained by NANOGrav (Arzoumanian et al. 2018b), to assign pulse phases to each event. The number of the "pulse invariant" (PI) detector channel provides an estimate of the photon energy: the photon energy in eV is approximately equal to 10× the PI channel number in which the event occurs.

As discussed in Bogdanov et al. (2019a), the final X-ray pulse profile was cross-verified by comparing it with an independently generated profile in which pulse phases were assigned to events using the PINT pulsar timing software. The pulse phases assigned to individual events by the two procedures differed by ≲1 μs, and the differences between the two folded X-ray pulse profiles were negligible. The final event list used in the analysis reported in this Letter is available at https://doi.org/10.5281/zenodo.3524457.

To incorporate the performance of the NICER XTI in our model, we used version 1.02 of the redistribution matrix file (RMF), which gives the probability that a photon in a particular energy range is registered by a particular detector channel, and version 1.02 of the ancillary response file (ARF), which provides the effective area of the telescope as a function of energy, as measured on axis.

Spectral calibration of the NICER XTI has been carried out primarily using observations of the Crab pulsar and nebula. The response below channel 40 is more uncertain than the response at higher energies, because of current imperfections in the modeling of instrumental effects such as the detector trigger efficiency. Also, because at these low energies the Crab spectrum is strongly attenuated by the interstellar medium, only a relatively small number of counts are available to determine the spectrum at these energies. Fits to NICER observations of the Crab pulsar and nebula in the energy channels 40–299 used in this study (see Section 5.2) show systematic residuals, but these are typically ≲5% (Ludlam et al. 2018). They are likely due to imperfect modeling of the microphysics of the concentrator optics and astrophysical features, such as oxygen edges and emission lines.

The soft X-ray pulse waveform of PSR J0030+0451 is thought to be produced by hot regions on its rotating surface created by inward moving particles created in electron-positron pair cascades interacting with the outermost layers of the neutron star (Arons 1981; Harding & Muslimov 2001). Models of global pulsar magnetospheres (see, e.g., Kalapotharakos et al. 2014; Philippov et al. 2015; Brambilla et al. 2018) have delineated the current-closure patterns of pulsars, which show outward currents in the magnetic polar regions balanced by return currents over the outer parts of the polar caps that are connected to the current sheet. Some of the particles that make up the return currents (those that are accelerated in the current sheet) are responsible for the pulsar's γ-ray and nonthermal X-ray emission (Kalapotharakos et al. 2018; Philippov & Spitkovsky 2018). The polar cap pair cascades that take place in both outgoing and return current regions (Timokhin & Arons 2013) supply pairs for radio emission. Some of the particles produced in these pair cascades are accelerated downward and are expected to deposit energy deep in the neutron star's atmosphere (see Tsai 1974 and Bogdanov et al. 2007). This energy is expected to heat the atmosphere at large optical depths, producing thermal soft X-ray emission from the stellar surface (Harding & Muslimov 2001, 2002). The spectrum of this emission is expected to depend on the chemical composition of the star's atmosphere and, possibly, on the strength and orientation of the magnetic field at the locations of the hot spots (Potekhin 2014).

The models of the thermal emission from these hot spots that we have used in the analysis of the soft X-ray waveform of PSR J0030+0451 that we report here assume that the upper atmosphere of the pulsar is fully ionized pure hydrogen and that, for our purposes, the magnetic fields in the hot spots can be neglected. We now discuss these assumptions.

Composition of the upper atmosphere. Our model atmospheres were constructed assuming that the gas in the atmosphere is fully ionized, and hence that the dominant opacity is free–free absorption (Ho & Lai 2001). While opacity tables for partially ionized matter exist (Iglesias & Rogers 1996; Badnell et al. 2005; Colgan et al. 2016), they do not cover the full range of temperatures and energies needed for the present analysis. However, at the temperatures ∼10⁶ K that are relevant for the present study, the fraction of hydrogen that is neutral is very low; comparison of the outward specific intensity given by our fully ionized hydrogen atmosphere model differs by at most 1%–2% at 0.5–1 keV from that for an atmosphere with an effective temperature of 10⁶ K given by a model constructed using the Opacity Project (OP) opacity table (Badnell et al. 2005), which allows for partial ionization (see Bogdanov et al. 2019b).

Our assumption that the upper atmosphere is pure hydrogen is based primarily on the astrophysical abundance of hydrogen and calculations that show how the lightest element present in a neutron star atmosphere will float to the top within seconds. These calculations are extrapolations of earlier calculations by, for example, Alcock & Illarionov (1980). Given the astrophysical abundance of hydrogen and its tendency to float to the top, we expect it to completely dominate the composition of the upper atmosphere.

If the outer atmosphere of the pulsar were instead dominated by helium, the emission spectrum and beaming pattern in the energy range 0.5–1.0 keV would be expected to be very similar to those for a hydrogen atmosphere, but their normalizations would be slightly smaller than for hydrogen (see Bogdanov et al. 2019b). In particular, hydrogen and helium atmospheres with an effective temperature of 10⁶ K (the approximate inferred temperature of the heated regions on the surface of PSR J0030+0451) are predicted to produce specific intensities normal to the surface that differ from each other by ≲5% in the 0.5–1 keV energy range (again see Bogdanov et al. 2019b). Our current best fit of the zero-field NSX helium atmosphere (Ho & Lai 2001) to the NICER data on PSR J0030+0451 appears disfavored relative to our current best fit of the zero-field NSX hydrogen atmosphere (Ho & Lai 2001; see also Ho & Heinke 2009). Specifically, our current fit using the NSX helium model yields a best-fit mass of 2.7 M_⊙, with little posterior probability below 2.0 M_⊙. Although more work would be needed to completely rule out a helium atmosphere for PSR J0030+0451, this mass estimate is so large that we interpret it as disfavoring a helium atmosphere.

Thus, while there is at present no direct spectroscopic confirmation that the upper atmosphere of PSR J0030+0451 is pure hydrogen, the fact that the spectrum and beaming pattern predicted by hydrogen atmosphere models agree well with the spectrum and beaming of the soft X-ray emission from PSR J0030+0451 observed using NICER, and that the X-ray spectrum and beaming predicted by helium atmosphere models appear disfavored, supports our assumption that the upper atmosphere is hydrogen. Therefore, in computing model waveforms we used the local emission for a given surface gravity, effective temperature, angle from the surface normal, and photon energy that is predicted by the NSX code developed by Ho & Lai (2001; see also Ho & Heinke 2009) for a fully ionized hydrogen atmosphere.

The 205.53 Hz rotation frequency of PSR J0030+0451 (Arzoumanian et al. 2018a) is small enough that the only effects of the rotation that need to be taken into account in computing pulse waveforms are the special relativistic Doppler shift and aberration caused by the motion of the emitting surface and the rotationally induced oblateness of the star; the effect of frame-dragging and the effect on the exterior spacetime of the rotationally induced stellar mass quadrupole are negligible (see Morsink et al. 2007 and Bogdanov et al. 2019b). We therefore used the so-called oblate-star Schwarzschild-spacetime (OS) approximation (Cadeau et al. 2007; Morsink et al. 2007; AlGendy & Morsink 2014) when tracing light rays from the stellar surface to the observer. Computations of pulse waveforms like those used here were verified by members of the NICER Lightcurve Working Group, as described in Bogdanov et al. (2019b).

Strength and structure of the stellar magnetic field. We can estimate an approximate lower bound on the strength of the surface magnetic field of PSR J0030+0451 by assuming that the field is a centered dipole and estimating the strength of the dipole moment using the best braking models currently available (see, e.g., Contopoulos & Spitkovsky 2006). The period P of PSR J0030+0451 is 0.00487 seconds and its period derivative $\dot{P}$ is 1.02 × 10⁻²⁰ (Arzoumanian et al. 2018a). If we use Equation (12) of Contopoulos & Spitkovsky (2006) and, for the sake of argument, choose a magnetic field inclination θ of π/2 and the field configuration described by α = 0, then the inferred strength of the magnetic field at the stellar surface is B ∼ 4 × 10⁸ G. In principle, this represents an approximate lower bound on the strength of the total magnetic field in the atmospheres of the hot spots. For comparison, the electron cyclotron energy in a magnetic field of this strength is ${\hslash }{\omega }_{c}\approx 4\,\mathrm{eV}$, which is more than 50 times smaller than the lowest photon energies that can be observed using NICER and ∼100 times smaller than the lowest energies that we consider in our fits.

We note, however, that the three hot spots in the "southern" rotational hemisphere (we denote the rotational hemisphere that includes the line of sight to the observer as the "northern" rotational hemisphere) that can explain the pulse waveform of PSR J0030+0451 observed using NICER suggest that the magnetic field of PSR J0030+0451 is at least an offset dipole, and may have higher moments that dominate the field strength near the stellar surface. A complex field would be consistent with the structures of the magnetic fields of other millisecond rotation-powered pulsars inferred from their radio and γ-ray emission (Ruderman 1991; Ruderman et al. 1998), and the inferred magnetic field configurations of the accretion-powered millisecond pulsars that are thought to be the progenitors of the millisecond rotation-powered pulsars (Lamb et al. 2009a, 2009b).

If the magnetic field of PSR J0030+0451 does have significant higher multipole moments, they may not be negligible at the light cylinder, which is closer to the stellar surface in millisecond pulsars than in pulsars with lower spin rates. Moreover, in some spin-down models (again see Contopoulos & Spitkovsky 2006), the spin-down rate depends on the magnetic flux enclosed by the open-field-line region, which is defined by the boundary of the closed magnetosphere and may not extend all the way out to the light cylinder. The way this flux evolves with spin rate can be affected by the presence of multipolar field components, complicating the relation between the strength of the surface magnetic field and the spin-down rate.

Given this complexity, it is difficult to infer the strength and configuration of the magnetic field on the surface of PSR J0030+0451. In current braking models, the braking torque scales with the strength of the dipole component of the magnetic field at some characteristic distance from the star. For a given strength of the dipole component at this distance, the strength of the magnetic field at the stellar surface scales as (R_e/d)³, where R_e is the circumferential radius of the star and 2d is the distance through the star from one magnetic pole to the other. Thus, if the large-scale magnetic field is essentially dipolar and the corresponding magnetic poles are a significant distance south of the rotational equator, the field strength at the surface of the star could be 5–10× greater than it would be for a centered dipole field with the same strength a large distance from the stellar surface. In this case, the indicated strength of the magnetic fields at the hot spots could be as large as B ∼ 4 × 10⁹ G.

Effect of the magnetic field on the soft X-ray emission from the stellar surface. The model atmospheres that we have used to compute the waveforms that we fitted to the NICER waveform data assume that the effects of the stellar magnetic field on the structure and radiative properties of the atmosphere are negligible. This approximation is valid for magnetic field strengths that are less than ${B}_{c}={e}^{3}{m}_{e}^{2}\,c/{{\hslash }}^{3}=2\times {10}^{9}$ G. The upper bound on the strength of the magnetic field that we estimated above by assuming that it is dipolar and using the best currently available braking models is an order of magnitude smaller than B_c, but our estimates of the possible field strengths at the hot spots are a factor of two larger than B_c.

If the magnetic field exceeds B_c, its effect is non-perturbative and dominant (see Lai 2001): the radiation field becomes polarized, the beaming pattern becomes strongly anisotropic, and spectral features associated with the electron cyclotron resonance and bound species become important (Potekhin 2014). The absence of any obvious spectral feature in the spectra observed using XMM-Newton (see, e.g., Bogdanov & Grindlay 2009) and in the ∼0.3–3 keV spectra obtained using NICER rules out B = (3–30) × 10¹⁰ G, assuming a redshift of 20%. In addition, while an electron cyclotron line produced by a magnetic field B ≫ 10¹¹ G could be hidden at E > 3 keV, such fields would produce pencil or fan beams that are unlikely to be able to produce an energy-dependent pulse profile that can match the profile observed using NICER.

The changes in the spectrum and the beaming pattern that occur for magnetic fields B ≳ 10¹⁰ G are somewhat larger than the changes in the spectrum and beaming pattern that occur when going from a purely H model atmosphere to a purely He model atmosphere without any magnetic field. While spectra that assume such strong fields might yield statistically acceptable fits to the NICER observations of PSR J0030+0451, it is likely, though not certain, that such fits would, as in the case of the He atmosphere models, yield mass and radius estimates that are so extreme that we would consider these models disfavored.

If the total strengths of the magnetic fields in the atmospheres of the hot spots of PSR J0030+0451 are less than B_c, they are expected to have little effect on the spectrum and shape of its soft X-ray waveform. The agreement of our energy-resolved waveform models with the energy-resolved waveform data obtained using NICER are consistent with this interpretation.

In modeling the soft X-ray pulse waveform of PSR J0030+0451, we started with the simplest possible description of the heated region(s) on the stellar surface and then, as necessary, considered more complicated descriptions, guided at each step by the differences between the best-fit simpler pulse waveform models and the waveform observed using NICER. Along the way, we verified our parameter estimation algorithms by fitting appropriate models to several synthetic waveforms that were generated using the models of the heated regions on the stellar surface that emerged from our fits to the observed waveform, and then comparing the values of the parameters obtained from these fits with the values of these same parameters that were used to generate the synthetic waveforms.

Visual examination of the PSR J0030+0451 waveform observed using NICER shows evidence for two peaks, and waveform models constructed using a single hot spot fail to reproduce the qualitative features of the observed waveform. Consequently, we began our modeling of the observed waveform using models that had two uniform circular hot spots with possibly different locations, areas, and temperatures, but these models did not adequately describe the observed waveform. In particular, the (relatively large) spot sizes that were required to reproduce the weak higher harmonic content of the observed waveform appeared to conflict with the (relatively small) spot areas that were required to reproduce the observed X-ray flux.

We therefore considered models with three and four different uniform-temperature circular hot spots, and finally models with two and three different uniform-temperature oval hot spots. In all cases we allowed the spots to overlap in any way that improved the fit; our procedure for treating overlapping spots is described below. This algorithm allowed the heated regions in the model to evolve toward a variety of different complex shapes as the fitting process proceeded, including multiple separate hot spots, multiple-temperature spots in which each spot has hotter and cooler regions, elongated or more complicated configurations produced by two or more partially or completely overlapping or adjacent circular or oval hot spots, crescent-shaped hot regions produced by cold oval spots partially or completely overlapping hot oval spots, and so on. These explorations found no evidence for different temperatures within the same hot spot and no evidence for more than three heated regions on the stellar surface.

We found that a model with two different, non-overlapping, uniform oval spots is preferred over any models with two overlapping spots, adequately describes the observed waveform, and gives a much better fit to the data than a model with two circular spots. The larger east–west extent of the best-fitting oval spots allows them to reproduce the observed weak higher harmonic content of the waveform, while their smaller north–south extent allows them to have the relatively small areas required to reproduce the observed X-ray flux.

To assess whether a model with two different uniform-temperature oval spots is sufficient, we also considered configurations with three different uniform-temperature oval spots. We found that a model with three different, non-overlapping, uniform oval spots is preferred over any models in which the spots overlap, describes the observed waveform adequately, and gives a fit to the data that is slightly, but not significantly, better than the best-fit model with two non-overlapping, uniform oval spots. The two main spots in the best-fit model with three spots are very similar to the spots in the two-spot model, while the third spot has a much higher temperature but a much smaller area than either of the other two spots, is located close to the rotational pole on the far side of the star from the observer, and makes only a very small contribution to the waveform.

The pulse waveforms and the mass and radius estimates given by the best-fit two- and three-spot models are statistically indistinguishable. We therefore concluded that both models adequately describe the observed waveform, that there is no evidence that a more complicated model is required, and that the mass and radius estimates given by these models are reliable. We slightly prefer the radius and mass estimates given by the model with three oval spots, only because the evidence for this model is slightly—although not significantly—greater than the evidence for the model with two oval spots.

In Section 4, we show and discuss the results of several of the analyses of synthetic waveform data that we performed to gain confidence in the accuracy and reliability of our parameter estimation procedure. Specifically, we show the results that we obtained by fitting a model of the waveform produced by two circular spots to synthetic waveform data generated using this model, and by fitting a model of the waveform produced by two oval spots to synthetic waveform data generated using this more complicated model. In Section 5, we describe in detail our systematic approach to modeling the pulse waveform of PSR J0030+0451 and the evidence that led us to focus on the waveform model with three oval spots that we used to estimate the radius and mass of PSR J0030+0451. But first we detail the various waveform models that we constructed, the modeling procedure that led us to consider these models, and the algorithms that we used to estimate the best-fit values of the parameters in these models. We also describe the procedures that we used to evaluate and compare these models, the motivation for our selection of the NICER data we analyzed, and our assessment of the effects of the errors in the NICER response.

In the analyses that we present here, we used pulse waveform models with the primary parameters listed and defined in Table 1. There are 12 primary parameters in the model that has two circular spots (this model assumes f₁ = 1 = f₂ and has no parameters for a third spot), 14 primary parameters in the model that has two oval spots, and 19 primary parameters in the model that has three oval spots.

Table 1.  Primary Parameters of the Pulse Waveform Models Considered in this Work

Parameter	Definition	Assumed Prior
GM/(c2Re)	Stellar compactness	0.125–0.3125
M	Gravitational mass	1.0–2.4 M⊙
θc1	Spot 1 center	0.1–3.1 rad
Δθ1	Spot 1 radius	0–3 rad
f1	Spot 1 oval ratio	0.1–20
kTeff,1	Spot 1 eff. temp.	0.011–0.5 keV
Δϕ2	Spot 2 longitude diff.	0–1 cycles
θc2	Spot 2 center	0.1–3.1 rad
Δθ2	Spot 2 radius	0–3 rad
f2	Spot 2 oval ratio	0.1–20
kTeff,2	Spot 2 eff. temp.	0.011–0.5 keV
Δϕ3	Spot 3 longitude diff.	0–1 cycles
θc3	Spot 3 center	0.1–3.1 rad
Δθ3	Spot 3 radius	0–3 rad
f3	Spot 3 oval ratio	0.1–20
kTeff,3	Spot 3 eff. temp.	0.011–0.5 keV
θobs	Observer inclination	0.1-π/2 rad
NH	H column density	0–2.5 × 1020 cm−2
D	Distance	$\exp [-{(D-0.325{\rm{kpc}})}^{2}/2{(0.017{\rm{kpc}})}^{2}]$

Note. We assumed the above priors, which are flat in the indicated range except for the distance D, and a rotation frequency of 205.53 Hz as seen by a distant observer (a) when we analyzed the synthetic PSR J0030+0451 waveform generated using two oval hot spots, which had parameter values (see Table 4) chosen so that it would mimic the waveform observed by NICER, and (b) when we analyzed the actual PSR J0030+0451 waveform observed by NICER. When we analyzed the synthetic waveform that was generated using two circular hot spots, we used the priors listed above for most of the parameters, but different priors for four parameters. This is because the waveform was constructed to mimic the waveform of PSR J0437−4715 observed by NICER, and in generating it we used values for these four parameters that are different from those listed in Table 2. In generating this synthetic waveform, we assumed a rotation frequency of 173.6 Hz and flat priors in the intervals listed for the following four parameters: M: 1.0–1.9M_⊙, θ_obs: 0.6–0.9 rad, N_H: 0–2.5 × 10²⁰ cm², D: 0.13–0.18 kpc. See the text for further details about the assumed prior probability distributions for each parameter.

Download table as:  ASCIITypeset image

For the reasons discussed previously, all the waveform models that we present in this report assume that the effective temperature of the emission from a given hot spot does not vary across the spot. The effective temperature listed in Table 1 is the effective temperature that would be measured by a local observer on the stellar surface.

When constructing an oval spot centered at a given colatitude and having a given latitudinal extent, we determined the longitudinal extent it would have if it were instead circular and then multiplied this extent by the colatitude-independent factor f₁ for spot 1, f₂ for spot 2, and f₃ for spot 3. The possible values of f₁, f₂, and f₃ are bounded from above by the requirement that no spot have a longitudinal extent greater than 2π.

We defined the northern hemisphere as the hemisphere containing the sightline to the observer. We defined longitude 0 as the longitude of the center of spot 1; after generating a trial waveform, we rotated its phase to give the best fit to the waveform data (see below).

As noted earlier, when fitting model waveforms to the observed waveform, we allowed hot spots to overlap one another. This approach allowed the fitting process to create more complicated heated regions, such as multiple isolated spots, spots with hotter cores and cooler annuli, spots with cooler cores and hotter annuli, noncircular heated regions having two different temperatures, and so on, if these configurations were favored by the data. For the reasons described in Section 3.2, we discuss here only our modeling of waveforms with at least two hot spots and at most three.

In the pixels on the stellar surface where spots overlap, our fitting algorithm chooses the temperature of the lowest-numbered spot. For example, if spot 2 and spot 3 overlap, all the pixels in the overlap region are assumed to emit with the effective temperature currently assigned to spot 2. If spot 1, spot 2, and spot 3 all overlap, all pixels in the overlap region are assumed to emit with the effective temperature currently assigned to spot 1.

In addition to the primary parameters listed in Table 1, we introduced two types of ancillary parameters. The first set of ancillary parameters describes the counts in each energy channel that do not come from any of the hot spots. We refer to these as background counts. These background counts allow for any unmodulated emission from the star that is not produced by any of the hot spots, any other unmodulated emission from the pulsar system, any emission within the field of view that does not come from the pulsar system, and the instrumental backgrounds. Another ancillary parameter describes the overall starting phase of the pulse waveform. In our fitting procedure, we chose to marginalize over these ancillary parameters, rather than fitting for them explicitly. Separately marginalizing the likelihood over the ancillary parameters is fully justified when, as here, the values of the primary parameters are uncorrelated with the values of the ancillary parameters.

We estimated the number of background counts in each energy channel and marginalized the likelihood with respect to the number of background counts in each channel as follows. We first generated a trial waveform by choosing a trial set of values for all the parameters in the waveform model and a trial starting phase for the waveform. Then, for each energy channel, we independently fit a Gaussian to the curve of the likelihood versus the number of background counts assumed present in that channel. By integrating over each Gaussian, we were able to estimate the number of background counts in each energy channel and marginalize the likelihood over the number of background counts in each energy channel, for each trial waveform. This procedure added one parameter to the waveform model for each energy channel that we considered. As we discuss in detail below, for our final fits to the NICER pulse waveform data we chose to use NICER energy channels 40 through 299. This approach therefore added 260 parameters to the waveform model.

We performed a similar procedure to determine the best-fit overall phase of the waveform and marginalize the likelihood over this phase. Namely, for each trial offset of the overall phase (defined as the difference between the longitude of the center of the first spot relative to the longitude of the line of sight to the observer), we computed the likelihood over all energy channels. We then fit a Gaussian to the likelihood as a function of the overall phase, which allowed us to estimate the overall phase. We then marginalized the likelihood with respect to the overall phase by integrating the phase offset over the Gaussian curve, assuming a prior for the overall phase that was uniform from 0 to 1 cycles and zero outside this range. This added one more parameter to the waveform model.

Our procedure for estimating the number of background counts in each NICER energy channel assumes that we have no independent knowledge of the number or spectrum of the background counts. In reality, we have some information about astrophysical and sky backgrounds and the backgrounds contributed by solar system, magnetospheric, terrestrial, and International Space Station sources. We have developed codes that can take this information into account, but many of these backgrounds are time-variable and are often not accurately known. The values of the waveform parameters that we derived when we included models of the properties of these backgrounds differ little from the values we found when we used the procedure just described. We have therefore chosen to adopt the conservative approach of using this procedure, which does not rely on any models of the various backgrounds.

Parameter estimation and model comparison. To estimate the values of the parameters in our pulse waveform models and compare the evidence for different models using synthetic or NICER data, we used standard Bayesian methods. For a given model ${ \mathcal M }$ with parameters ${\boldsymbol{\theta }}$ that have a normalized prior probability density $q({\boldsymbol{\theta }})$, the normalized posterior probability density is given by Bayes' Theorem:

$$\begin{eqnarray}&&P({ \mathcal M },{\boldsymbol{\theta }}| \mathrm{data})=\displaystyle \frac{{ \mathcal L }(\mathrm{data}| { \mathcal M },{\boldsymbol{\theta }})q({\boldsymbol{\theta }})}{P({ \mathcal M }| \mathrm{data})},\end{eqnarray} \tag{ 1 }$$

where ${ \mathcal L }(\mathrm{data}| { \mathcal M },{\boldsymbol{\theta }})$ is the likelihood of the data given the model with parameter values ${\boldsymbol{\theta }}$, and $P({ \mathcal M }| \mathrm{data})\,\equiv {\int }_{{\boldsymbol{\theta }}}{ \mathcal L }(\mathrm{data}| { \mathcal M },{\boldsymbol{\theta }})q({\boldsymbol{\theta }})d{\boldsymbol{\theta }}$ is the evidence for model ${ \mathcal M }$.

Given a model ${ \mathcal M }$, we can compute the joint posterior probability density distribution of the values of any subset of the parameters ${\boldsymbol{\theta }}$ by marginalizing the other parameters, i.e., by integrating $P({ \mathcal M },{\boldsymbol{\theta }}\,| \,\mathrm{data})$ over the possible values of all parameters except those in the subset. For example, the 1D probability density distribution of the value of any single parameter can be computed by integrating $P({ \mathcal M },{\boldsymbol{\theta }}\,| \,\mathrm{data})$ over the possible values of all the other parameters.

The relative probability of two different models is given by their odds ratio, which is the ratio of their evidences (their Bayes factors) multiplied by the ratio of the prior probabilities of the two models.

Sampling methods. We used three statistical sampling methods to estimate the values of the parameters in our various waveform models and compute the evidence for these models: the publicly available nested sampler MultiNest (MN; Feroz et al. 2009), the parallel-tempering Markov chain Monte Carlo sampler (PT-emcee) that is included in the publicly available emcee package version 2.2.1 (Foreman-Mackey et al. 2013), and a hybrid sampling algorithm (see below) that we call MN+PT-emcee.

The MN and PT-emcee sampling algorithms are complementary: nested samplers are optimized for computing the Bayesian evidence and are therefore well suited for computing the odds ratios of different models, whereas Markov Chain Monte Carlo algorithms are optimized for parameter estimation. We carried out many tests of MN, PT-emcee, and MN+PT-emcee in the context of our waveform-fitting problem, using both synthetic and NICER waveform data sets.

For relatively simple waveform models and data sets, such as those with two uniform circular hot spots, we were able to adequately sample the full parameter space using either MN or PT-emcee in a reasonable amount of clock time. For these simpler models and data sets, the posterior probability density distributions given by these two algorithms appeared well converged and were in excellent agreement.

For more complex models and data sets, such as our fits of models with two or three oval spots to synthetic data or to the NICER data on PSR J0030+0451, we found that, at least for the values of its sampling parameters that we explored, the MN algorithm abandoned regions of parameter space that contained solutions with log likelihoods comparable to the highest found in previous searches, did not give reproducible results, and—when we used it to analyze synthetic data—produced 1σ and 2σ credible regions that too often did not contain the values of the parameters that had been assumed in generating the synthetic data, even after we ran the code for a very long time.

In contrast, even with poor initial seeds, the PT-emcee algorithm appeared able to sample well the parameter space of complicated models and data sets, gave reproducible results, and produced 1σ and 2σ credible regions that contained the values of the parameters that had been assumed in generating the synthetic data about as often as one would expect, i.e., about 68% of the time for 1σ credible regions and about 95% of the time for 2σ credible regions. We therefore adopted the PT-emcee algorithm for analyzing the more complex models and data sets that we considered. However, with poor initial seeds, PT-emcee took large amounts of clock time to complete each parameter estimation.

In order to speed up parameter estimation using PT-emcee, we sought to choose high likelihood points as initial guesses. Having the MN nested sampling algorithm readily available, we typically used it to provide an initial guess. We tested this hybrid approach, which we refer to as MN+PT-emcee, extensively, and found that it was able to find correct solutions in much less clock time than was required if PT-emcee was not given a high likelihood initial guess. Consequently, we used MN+PT-emcee to estimate the values of the parameters in the models with two and three oval spots, and used the resulting samples to estimate the evidence for these models, given the data.

In our initial MultiNest surveys, we used 1000 active points (essentially the number of parameter combinations being actively updated), an efficiency of 0.01 (the target fraction of new parameter combinations that update the list of active points), and a tolerance of 0.1; for precise definitions of these parameters; see Feroz et al. (2009). Based on our experience sampling the synthetic and NICER waveform data, we limited the number of modes to five for models with two spots and 10 for models with three spots. In our PT-emcee analyses, we used 400 walkers when fitting models with two spots and 1000 walkers when fitting models with three spots, and the default temperature ladder.

Completeness of models. Although our parameter estimations and model comparisons are fully Bayesian, we also performed χ² analyses to determine whether our models adequately describe the NICER data. In computing χ², we use the expression for model variance originally advocated by Pearson (1900), namely,

$$\begin{eqnarray}&&{\chi }^{2}=\displaystyle \sum _{i}\displaystyle \frac{{\left({m}_{i}-{d}_{i}\right)}^{2}}{{m}_{i}},\end{eqnarray} \tag{ 2 }$$

where m_i is the number of counts in bin i that are expected in the model, and d_i is the observed number of counts in bin i, rather than the common data variance form in which the m_i in the denominator is replaced by d_i. The model variance form has the advantage that if the counts d_i are drawn from a Poisson distribution with expected values m_i, then the expected value of χ² is equal to the number of bins for any values of m_i. This is not true for the data variance form. We note that for our purposes the computation of χ² is a one-way test: if the value of χ² for a given number of degrees of freedom is large enough to be highly improbable then we are warned that our model of the spots and/or the instrument is likely to be incomplete, but if the value of χ² is reasonably probable we cannot conclude that the model is correct.

The method that we used to make the final parameter estimates we report here takes into account the substantially larger uncertainties in the NICER response below 0.4 keV (channel 40) and the evidence we found that including data from these channels might bias our results on PSR J0030+045.

As we discuss in Section 5.2, we found that including data from channels below channel 40 biases waveform-based estimates of the distance to PSR J0030+0451 toward values much larger than the independently known distance. The source of this bias is currently unknown. Since its origin is unknown, it could also bias waveform-based estimates of other parameters, including the radius and mass of the star. Given this risk, we judged it important to address this issue.

It might be possible to gain a better understanding of the origin of this distance bias if we also had reliable, independent estimates of the values of other stellar and waveform parameters, but for PSR J0030+0451, we only have a reliable independent estimate of the distance. This provides too little information to identify the source(s) of this bias or to devise sophisticated corrections. We therefore considered only simple ways to address this bias. We also wanted to be careful not to do anything that might bias other, even more important model parameters.

As was mentioned in Section 2.3, the NICER response below channel 40 is more uncertain than the response above channel 40. It was therefore natural to explore the effect of ignoring the data in channels below channel 40. As we show in detail in Section 5.2, discarding the data in these low channels eliminates the bias in our waveform-based estimates of the distance to PSR J0030+045. It is our hope that discarding these data also reduces any potential biases in our estimates of the other parameters in the waveform model. Because this data selection procedure is so simple and does not touch the data that we do use, we think it is less likely to bias other waveform parameters than a more complex procedure.

In analyzing the NICER data on PSR J0030+0451, we found no detectable modulation in the energy channels above channel 299, and we therefore did not include data from these channels in our analyses of the NICER data on PSR J0030+0451.

One would expect that discarding the data in channels 25–39 would reduce the precisions of the parameter estimates. We investigated this possibility by comparing the precisions of the parameter estimates obtained by analyzing the data in channels 40–299 produced by a synthetic waveform that is similar to the PSR J0030+0451 waveform observed by NICER with the precisions obtained by analyzing the data in channels 25–299 produced by the same synthetic waveforms and found that discarding the data in channels 25–39 did not significantly reduce the precisions of parameter estimates (see Section 4.2). We also compared the posterior probability distributions for the radius and mass of PSR J0030+0451 obtained using the NICER data in channels 40–299 with the probability distributions obtained using the data in channels 25–299 and found no evidence that the precisions of the parameter estimates were reduced by discarding the data in channels 25–39 (see Section 5.2).

A remaining question is whether the fits of our waveform models to the NICER waveform data on PSR J0030+0451 are affected by the unavoidable systematic errors present in our model of the NICER response in the energy channels 40–299 that we used. As we described in Section 2.3, the NICER response in these channels was calibrated by comparing the results of NICER observations of the Crab Nebula, which is traditionally assumed to have a power-law spectrum, with the power-law spectrum of the nebula determined by observations made using other instruments. When NICER observations of the Crab Nebula spectrum are compared with the spectrum obtained by folding the power law determined by other instruments through the nominal NICER response, the residuals in the energy channels we used in this study (channels 40–299) are ≲5% (Ludlam et al. 2018). As mentioned in Section 2.3, these residuals are likely due to imperfect modeling of the microphysics of the concentrator optics and astrophysical features, such as oxygen edges and emission lines.

We explored the likely effect of deviations of the nominal NICER response from the true NICER response curve by multiplying the nominal effective area curve of NICER by the ratio between the count rate as a function of energy given by the assumed spectrum of the Crab Nebula and the count rate as a function of energy measured by NICER (Ludlam et al. 2018). We then used this modified effective area curve to analyze both synthetic waveform data constructed to mimic the NICER data on PSR J0030+0451, and the actual NICER waveform data on PSR J0030+0451.

We found that modifying the nominal effective area in this way had a negligible effect on the quality of our fits to synthetic waveform data and to the actual waveform data on PSR J0030+0451. For example, the probability of the χ² obtained by comparing the energy-resolved waveform given by the best-fit model with two oval spots to the NICER data in energy channels assigned to 32 phase bins changed from 0.105 to 0.112 when we used the modified effective area curve instead of the nominal effective area curve. This indicates that the current systematic errors in the nominal effective area of NICER in channels 40–299 have only a minor effect on our fits.

Given these results, we chose to ignore all events in the channels below channel 40 and use the current nominal NICER effective area curve for channels 40–299, both when analyzing the actual NICER data on PSR J0030+0451 and when analyzing synthetic NICER data.

The first analysis of synthetic data we present here shows the results we obtained by fitting a 12-parameter waveform model with two possibly different and overlapping uniform circular spots to a synthetic waveform generated assuming two different, non-overlapping uniform circular spots, using our MN+PT-emcee sampling algorithm. Table 2 lists the values of the waveform parameters that we assumed when we generated the synthetic waveform data. We note that the temperature we assumed for one of the two spots when we generated this synthetic waveform is substantially less than the temperatures of either of the two main spots we found when we analyzed the NICER data on PSR J0030+0451. We show our results for this synthetic waveform because it was the synthetic waveform that was selected—in advance of the analyses of the PSR J0030+0451 waveform—to verify the parameter estimation procedures that would be used by the NICER team to analyze pulse waveforms. The results obtained by others on the team when they analyzed this waveform are consistent with the results we obtained when we analyzed it (Bogdanov et al. 2019b).

Table 2 shows the results that we obtained by fitting the waveform model with two possibly different and overlapping uniform circular spots to synthetic waveform data generated assuming two different, non-overlapping uniform circular spots. The table lists the median value of each model parameter computed using its 1D posterior probability distribution and the boundaries of the ±1σ and ±2σ credible intervals for each parameter.

We found that when this waveform model was fit to the synthetic waveform data in NICER energy channels 25–299, the values of 10 of the 12 parameters assumed in generating the synthetic waveform are within the relevant ±1σ credible intervals for these parameters while the values of the remaining two parameters (θ_c1 and N_H) assumed in generating the synthetic waveform are within their ±2σ credible intervals. We found that when this model was instead fit only to the synthetic data in NICER energy channels 40–299, the assumed values of seven of the 12 parameters were within the relevant ±1σ credible intervals, the assumed values of four others were within the relevant ±2σ credible intervals, and the assumed value of the other parameter was within the relevant ±3σ credible interval. Thus, for both energy channel choices, the 1D posterior probability distributions derived for the waveform parameters are consistent with the values of the parameters assumed when the synthetic waveform was generated. These results are summarized in Table 3.

Figure 1 shows that the joint probability density distribution for M and R_e that was obtained by fitting this model to the synthetic waveform data in NICER energy channels 25–299 is consistent with the values of the stellar mass and radius assumed when the synthetic waveform was generated.

Zoom In Zoom Out Reset image size

The panels in the top row of Figure 2 show the 1D probability distributions of M, R_e, and D obtained by fitting this model to the synthetic waveform data in NICER energy channels 25–299, whereas the panels in the bottom row show the results obtained by fitting to the data in channels 40–299. This figure illustrates two noteworthy aspects of these distributions.

Zoom In Zoom Out Reset image size

First, these 1D probability distributions are consistent with each other and with the values of M, R_e, and D that were assumed in generating the synthetic waveform. The widths of the distributions for M and R_e obtained using only the data in channels 40–299 underestimate the true uncertainties in the estimates of these parameters, because the temperature assumed for the lower-temperature spot when generating this synthetic waveform was sufficiently low that discarding the data in energy channels 25–39 narrowed these distributions. Our tests (see, e.g., Figure 4 below) have shown that—in contrast to the reduced widths of the M, R_e, and D distributions obtained from this synthetic waveform when the data in channels 25–39 are discarded—discarding the data in these channels has no significant effect on the distributions of M, R_e, and D obtained by analyzing the NICER data on PSR J0030+0451, which has higher-temperature hot spots.

Second, although the 1D probability distributions for D obtained by fitting this model to data selected using both these energy-channel criteria are consistent with the value of D assumed in generating the synthetic waveform data (the cumulative probability near the assumed value of D is substantial for both channel selections and fits with wider priors show that the probability density eventually decreases with increasing distance), these results show that pulse waveform analysis does not determine D very precisely. These analyses demonstrate the importance of using the precise and accurate independent measurement of D that is available for PSR J0030+0451 (see below) when analyzing the pulse waveform of this pulsar.

The second analysis of synthetic data we present here shows the results we obtained by fitting a 14-parameter waveform model with two possibly different and overlapping uniform oval spots to a synthetic waveform generated assuming two different, non-overlapping uniform oval spots, using our MN+PT-emcee sampling algorithm. Table 4 lists the values of the waveform parameters that we assumed when we generated the synthetic waveform data.

The results that we obtained by fitting this waveform model to the synthetic waveform data in NICER energy channels 25–299 and 40–299 are shown in Table 4, which lists the median value of each model parameter computed using its 1D posterior probability distribution and the boundaries of the ±1σ and ±2σ credible intervals for each parameter.

We found that when the model waveform was fit to the synthetic waveform data in NICER energy channels 25–299, the values of 10 of the 14 parameters assumed in generating the synthetic waveform are within the ±1σ credible intervals for these parameters, the values assumed for three others are within the relevant ±2σ credible intervals, and the value assumed for the remaining parameter is within the relevant ±3σ credible interval. When this model waveform was fit only to the synthetic data in NICER energy channels 40–299, the assumed values of 12 of the 14 parameters fall within their ±1σ credible intervals and the assumed values of the other two parameters fall within their ±2σ credible intervals. The waveform parameter values that were assumed when generating the synthetic waveform data are all consistent with the 1D posterior probability densities obtained by analyzing the data using either energy cut. These results are summarized in Table 5.

Figure 3 shows that the joint probability density distribution for M and R_e obtained by fitting this model to the synthetic waveform data in NICER energy channels 25–299, and by fitting it to the synthetic data in channels 40–299, are both consistent with the values of the stellar mass and radius assumed when the synthetic data was generated. Figure 3 also demonstrates that discarding the data in NICER energy channels 25–39 does not significantly degrade the precision of the M and R_e estimates.

Zoom In Zoom Out Reset image size

The panels in the top row of Figure 4 show the 1D probability distributions of M, R_e, and D obtained by fitting this model to the synthetic waveform using the data in NICER energy channels 25–299, whereas the panels in the bottom row show the results obtained using only the data in energy channels 40–299. We note several aspects of these distributions.

Zoom In Zoom Out Reset image size

First, these 1D probability distributions are consistent with each other and with the values of M, R_e, and D that were assumed in generating the synthetic waveform. Second, these distributions are not significantly affected by using only the data in energy channels 40 and above. Finally, although the 1D probability distributions for D obtained by fitting this model to these data are consistent with the value of D assumed in generating the synthetic waveform (the cumulative probability near the assumed value of D is substantial for both energy cuts and fits with wider priors show that the probability density eventually decreases with increasing distance), waveform analysis again does not determine D very precisely. This further emphasizes the value of using the precise and accurate independent measurement of D (see below) when analyzing the NICER waveform data on PSR J0030+0451.

Our algorithm for developing a successful model of the heated regions on the surface of PSR J0030+0451 was to proceed systematically from the simplest possible model to more complicated models, guided at each step by the differences between the waveform given by the best-fit simpler models and the waveform observed using NICER, and ending when adding additional hot spots or complexity fails to increase significantly the evidence for the model. The sequence of models that we considered in our quantitative modeling is presented in Table 6. We now discuss this sequence.

As mentioned in Section 3.2, visual examination of the PSR J0030+0451 waveform shows evidence for two peaks. As a result, waveform models constructed using a single hot spot fail to reproduce the qualitative features of the observed waveform. We therefore began our quantitative modeling of the observed waveform using models with two uniform circular hot spots, allowing the spots to have different locations, areas, and temperatures. We also allowed the two circular spots to overlap partially or completely. The sampling algorithm therefore explored configurations in which the two circular spots (1) are disjoint, with the same or different temperatures; (2) partially overlap, potentially producing a somewhat oval-shaped uniform-temperature spot if the two component spots have the same temperature, or a spot with a complicated boundary and a complex temperature structure, if the two component spots have different temperatures; or (3) fully overlap, producing (3a) a single, uniform-temperature circular spot, if one of the spots completely covers the other, (3b) a spot with a centered or off-center core-annulus structure in which the core is hotter or colder than the annulus, if the two spots have different temperatures, or (3c) a circular spot partially surrounded by a hotter or colder crescent-shaped region. None of these models produced waveforms that adequately describe the observed waveform. Among other problems, the (relatively large) spot sizes required to reproduce the weak higher harmonic content of the observed waveform appeared to conflict with the (relatively small) spot areas required to reproduce the observed X-ray flux. We therefore considered more complicated heated regions.

In the second step of our exploration of waveform models, we considered models with three and four uniform circular hot spots. In all cases we allowed the spots to overlap in any way that improved the fit. This algorithm allowed the heated regions in the model to evolve toward a variety of different complex shapes as the fitting process proceeded (see Section 3.2). These explorations found no significant evidence for configurations with more than three circular hot spots on the stellar surface, for configurations in which three hot spots were near one another, or for overlapping or adjacent hot spots with significantly different temperatures, such as one would expect if there were spots on the stellar surface that had regions with different temperatures or a significant temperature gradient. We again found that the larger spot sizes required to reproduce the weak higher harmonic content of the waveform appeared to conflict with the smaller sizes required to reproduce the observed X-ray flux. This conflict led us to consider uniform-temperature oval hot spots as the next step in our sequence of trial models.

The first oval hot-spot model that we considered examined heated regions that could be modeled using two uniform-temperature oval spots. We again allowed the two spots to have different properties and to overlap one another, partially or completely. The sampling algorithm could have collapsed this model to a model with two uniform-temperature circular spots, if that was preferred by the data, but it was not. The sampling algorithm therefore explored configurations in which the two oval spots (1) are disjoint, with the same or different temperatures; (2) partially overlap, potentially producing a uniform-temperature spot with a complicated boundary, if both spots have the same temperature, or a spot with a complicated boundary and a complex temperature structure, if the two spots have different temperatures; or (3) fully overlap, producing (3a) a single, uniform-temperature oval spot, if one of the spots completely covers the other, (3b) a possibly oval spot with a potentially oval centered or off-center core-annulus structure in which the core is hotter or colder than the annulus, if the two spots have different temperatures, or (3c) a possibly oval spot partially surrounded by a hotter or colder crescent-shaped region.

Our analysis showed that a configuration with two non-overlapping, uniform oval spots with almost the same temperatures is preferred over any models with two overlapping spots, adequately describes the observed waveform, and is strongly favored over the best-fit model with two circular spots. Indeed, the log evidence for the best-fit model with two oval spots was 11.6 larger than the log evidence for the best-fit model with two circular spots. The fit with two uniform oval spots favors a configuration in which both spots are elongated in the east–west direction and both have almost the same temperature. The larger east–west extent of the best-fitting oval spots allows them to successfully reproduce the observed weak higher harmonic content of the waveform while the smaller north–south extent allows them to have the relatively small areas required to reproduce the observed X-ray flux.

To assess whether a model with two uniform-temperature oval spots is sufficient, we also considered configurations with three, possibly different, uniform-temperature oval spots. We found that a configuration with three different non-overlapping oval spots also describes the observed waveform adequately. This model is again strongly favored over the best-fit model with two uniform circular spots: the log evidence for the best-fit model with three oval spots was 13.3 greater than the log evidence for the best-fit model with two circular spots. We again found no evidence for more than one temperature in any of the heated regions. Two of the three spots in the best-fitting three-spot model are almost identical to the two spots in the best-fitting two-spot model. Like the two spots in the best-fitting two-spot model, the two larger spots in the best-fitting three-spot model are elongated in the east–west direction and have almost the same temperatures. The third spot in the best-fitting three-spot model has a much higher temperature but a much smaller area than either of the other two spots, is located close to the rotational pole on the far side of the star from the observer, and makes only a very small contribution to the waveform. The waveform produced by the best-fitting model with three oval spots is very similar to the waveform produced by the model with two oval spots.

The best-fitting models with two and three oval spots both provide acceptable fits to the phase-channel and bolometric waveform data in the χ² sense, and are not statistically distinguishable from each other. Most importantly, the best-fit values of the stellar radius and mass given by the model with three oval spots are almost identical to, and are statistically indistinguishable from, the best-fit values of the stellar radius and mass given by the model with two oval spots. The evidence for the best-fit model with three oval spots is slightly, but not significantly, greater than the evidence for the model with two oval spots.

Based on these results, we conclude that the models with two and three oval spots both provide adequate descriptions of the observed waveform, that there is no evidence that a more complicated model is required to describe the observed waveform, and that the mass and radius estimates inferred using these two models are reliable. We slightly prefer, and therefore report here, the radius and mass estimates given by the model with three oval spots, only because the evidence for this model is slightly—although not significantly—greater than the evidence for the model with two oval spots.

These considerations and results are summarized in Table 6. We now provide the details.

The lowest-energy NICER data that was included in the standard pipeline data set that we analyzed was in energy channel 25 (see Section 2.2). As noted previously (see Section 3.6), we found no detectable modulation of the X-ray flux in NICER energy channels 300 and above. We therefore considered in our analyses only the data in channels 25–299.

The calibration errors in the data from energy channels 40–299 are relatively small. For example, the residuals between the Crab Nebula spectrum determined from NICER measurements and the Crab Nebula spectrum determined from measurements made using other instruments are ≲5% for channels 40–299 (see Sections 2.3 and 3.6; see also Ludlam et al. 2018). In contrast, there are still substantial and poorly understood errors in the calibration of energy channels 25–39. We were therefore concerned that including data from channels 25–39 in our analysis would bias our estimates of the mass and radius of PSR J0030+0451.

The distance to PSR J0030+0451 has been accurately measured to be 0.325 ± 0.009 kpc, independent of any waveform analyses (see Arzoumanian et al. 2018a). We therefore decided to investigate the effect of including or not including the data in channels 25–39 on estimates of the distance made using the PSR J0030+0451 waveform. We found that if we used the NICER waveform data in channels 25–299, thereby including the data in channels 25–39, and assumed a wide, flat prior for the distance, all our analyses gave posterior probability distributions for the distance that peaked at distances substantially larger than the independently measured distance of 0.325 kpc and excluded 0.325 kpc with high confidence. This is illustrated by panel (c) in Figure 5, which shows the 1D posterior distributions for M, R_e, and D obtained by fitting the waveform model with two possibly different and overlapping uniform oval hot spots to the NICER data in channels 25–299. In contrast, when we used only the well-calibrated data in channels 40–299, we found posterior probability distributions for the distance that naturally peak close to the independently measured distance. This is illustrated by panel (f) in Figure 5, which shows the 1D posterior distributions for M, R_e, and D obtained by fitting the waveform model with two possibly different and overlapping uniform oval hot spots to the NICER data in channels 40–299.

Zoom In Zoom Out Reset image size

These results suggested to us that including the poorly calibrated data in channels 25–39 in the analysis was likely to be biasing the distance estimate and could also be biasing the estimates of M and R_e, about which we have no independent information. We therefore thought it was important to address the effects of the apparent calibration errors in channels 25–39, but we wanted to do this in a way that would not itself bias our results. Consequently, we decided that, rather than attempting to correct for the calibration errors in channels 25–39, which are poorly understood at present, the safest approach would be to discard the data in all these channels. This would prevent the calibration errors in these channels from biasing our M and R_e estimates, but it could in principle increase the uncertainties in our estimates of these parameters because the amount of data we used would be reduced.

Comparison of the results obtained by analyzing synthetic waveform data in NICER energy channels 25–299 and 40–299 that were generated assuming two different, non-overlapping uniform circular spots (see Section 4.1) shows that using only the data in the reduced energy range does not significantly bias the estimates of M and R_e or decrease the precisions of the estimates (compare panels (a) and (d) and panels (b) and (e) of Figure 2)). (We note that for this synthetic waveform, the temperature assumed for one of the hot spots was substantially lower than any of the temperatures inferred for the hot spots of PSR J0030+0451. Consequently, the cut at channel 40 discarded much more of the spectrum of this waveform than it does for the spectrum of PSR J0030+0451 observed by NICER.)

Similarly, comparison of the results obtained by analyzing synthetic waveform data in NICER energy channels 25–299 and 40–299 that were generated assuming two different, non-overlapping uniform oval spots (see Section 4.2) shows that using only the data in the reduced energy range does not significantly bias the estimates of M and R_e or decrease the precisions of the estimates (compare panels (a) and (d) and panels (b) and (e) of Figure 4).

Finally, comparison of the results obtained by analyzing the NICER data on the waveform of PSR J0030+0451 in energy channels 25–299 and 40–299 shows that using only the data in the reduced energy range does not significantly decrease the precisions of the estimates of M and R_e (compare panels (a) and (d) and panels (b) and (e) of Figure 5).

Taking into account all these results, we decided to use only the data in energy channels 40–299 to estimate the values of the parameters in our models of the PSR J0030+0451 pulse waveform. Consequently, in the remainder of this Letter we present only results obtained using the data in these energy channels. Given that the independently measured distance of 0.325 ± 0.009 kpc is very precise, we imposed a Gaussian prior on the distance that uses this measurement and an uncertainty that takes into account the uncertainty in this measurement and the global (in energy) uncertainty in the effective area of the NICER instrument (see below).

In this section, we present and discuss the results of our analysis of the NICER data on PSR J0030+0451, focusing on our estimates of its mass and radius. For the reasons described in Section 5.2, we present only the results we obtained by fitting the NICER data in energy channels 40–299. When estimating the values of all the parameters except the distance, we assumed prior probability distributions constant within the intervals listed in Table 1 and zero outside these intervals. When estimating the distance, we assumed a Gaussian prior. We chose the width of the Gaussian by combining linearly, rather than quadratically, the 0.009 kpc uncertainty in the measured distance to PSR J0030+0451 (see Arzoumanian et al. 2018a) and the estimated uncertainty in the overall normalization of the NICER effective area, converted to an uncertainty in the distance. We assumed a 5% uncertainty in the overall normalization of the NICER effective area (see Section 3.6). This translates into a distance uncertainty of 2.5%. The distance to PSR J0030+0451 has been independently measured to be 0.325 kpc (Arzoumanian et al. 2018a). Consequently, a 1σ uncertainty of 2.5% in the distance is ≈0.008 kpc, yielding a total 1σ width for the Gaussian of 0.009 kpc + 0.008 kpc =0.017 kpc.

We considered two waveform models, one with two possibly different and overlapping uniform-temperature oval spots, and one with three possibly different and overlapping uniform-temperature oval spots. The model with two oval spots has 14 primary parameters. These are listed in Table 7, along with their median values and the ±1σ and ±2σ boundaries of their 1D credible regions, computed using the 1D posterior probability density distribution we found for each parameter. The model with three oval spots has 19 primary parameters. These are listed in Table 8, along with their median values and the ±1σ and ±2σ boundaries of their 1D credible regions, computed using the 1D posterior probability density distribution we found for each parameter.²⁰ Both fits favored models in which the spots do not overlap.

Table 9 compares the median values of R_e and M and the boundaries of the ±1σ credible regions computed from the posterior probability density distributions obtained by fitting the two models described here.

Figure 6 compares the 1D posterior probability density distributions for M, R_e, and D. Note that the waveform data modifies the posterior probability distributions for the distance only slightly compared to the Gaussian prior that we assumed based on the independently measured distance to PSR J0030+0451 and the uncertainty in the NICER effective area (see Section 5.2), showing that the distance estimates obtained by fitting both models to the waveform data are consistent with the measured distance.

Zoom In Zoom Out Reset image size

Figure 7 compares the joint posterior probability density distributions for R_e and M that were obtained by fitting the two models to the NICER waveform data.

Zoom In Zoom Out Reset image size

Figure 8 compares the bolometric waveforms given by the best-fit waveform models with two and three oval spots. The waveform components are very similar for the two models, and the full waveforms are almost identical.

Zoom In Zoom Out Reset image size

All these comparisons show that the models with two and three uniform-temperature oval spots give mutually consistent credible regions for all the parameters that they have in common, most notably for the mass and radius of the pulsar. The log evidence for the best fit of the model with three oval spots is higher than the log evidence for the best fit of the model with two oval spots, but only by 1.7, which is less than the estimated uncertainty in the log evidence. Thus, both models are comparably good. This consistency suggests that, at least for the types of waveform models considered here, no further complexity is needed to adequately describe the observed waveform. The 68% credible intervals for M and R_e given by the fit of the model with three oval spots to the NICER data from PSR J0030+0451 are 1.30–1.59 M_⊙ and 11.96–14.26 km.

Our results for mass and radius of PSR J0030+0451 may be compared with the recent results of Steiner et al. (2018), who examined the evidence from neutron stars that are transient X-ray sources and conclude that the radius of a 1.4 M_⊙ neutron star is most likely between 10 and 14 km, which is consistent with our estimates of the mass and radius of PSR J0030+0451.

Our results may also be compared with those of Nättilä et al. (2017), who estimated that the mass of the (neutron star) X-ray burst source 4U 1702−429 is between 1.6 and 2.2 M_⊙ and that its circumferential radius is between 12.0 and 12.8 km (with 68% credibility). Although this estimate of the mass of 4U 1702−429 is higher than our estimate of the mass of PSR J0030+0451, for most neutron star EoS currently under consideration, the radius of a neutron star depends only weakly on its mass in this mass range, in which case this estimate of the radius of 4U 1702−429 is consistent with our estimate of the radius of PSR J0030+0451.

The full posterior probability density distributions for the parameters in the fit of the model with two oval spots to the NICER waveform data are shown in Figure 15 in Appendix A. The full posterior probability density distributions for the parameters in the fit of the model with three oval spots are shown in Figure 16 in Appendix B. Among other things, these "corner" plots show that there is a correlation between estimates of the stellar radius and the distance, but that this correlation is fairly weak. Thus, any errors in the data that bias the estimate of the distance may not significantly bias the estimate of the radius. There are indications that our sampling of the posterior distribution of the stretching parameter for the second oval spot did not explore a wide enough range, but this probably did not significantly affect our results, because there is no correlation between the value of this parameter and the values of the mass, radius, or distance parameters.

We performed χ² analyses (see Equation (2)) to determine whether the models we have used to estimate the radius and mass of PSR J0030+0451 adequately describe the NICER data. Specifically, we compared our best-fit energy-resolved model waveforms with the energy-resolved waveform data (i.e., the pulse-phase–energy-channel data) collected by NICER, and compared the energy-integrated (bolometric) waveforms given by our best-fit energy-resolved waveform models with the bolometric waveforms observed by NICER, using the values of χ² given by these comparisons. (In principle, we could also compare the pulse-phase-integrated spectra given by our models with the pulse-phase-integrated spectrum observed by NICER by computing the relevant values of χ², but our procedure for modeling the observed unmodulated counts—see Section 3.4—guarantees a nearly perfect description of the photon energy spectrum, so this comparison would be uninformative.)

Assigning the NICER data to 32 phase bins and 260 energy channels, the best-fit energy-resolved waveform model with two oval spots gives a χ² of 8204.68 when compared with this data set, which consists of the number of counts in each of 32 × 260 = 8320 phase-energy bins. As discussed in Section 3.3, the model waveform with two oval spots has 14 primary parameters, 260 ancillary parameters related to the non-star emission, and one additional parameter that describes its overall phase, yielding a total of 8320-260-14-1 = 8045 degrees of freedom. The resulting χ²/degrees of freedom (dof) is therefore 8204.68/8045. If this model is correct, the probability of finding a value of χ²/dof this large or larger is 0.104. Thus, according to the χ² test this model provides an acceptable description of this data set.

Again assigning the NICER data to 32 phase bins and 260 energy channels, the best-fit energy-resolved waveform model with three oval spots gives a χ² of 8188.99 when compared with the NICER data binned in this way. This model has 19 primary parameters, yielding a total of 8040 degrees of freedom. The resulting χ²/dof is therefore 8188.99/8040. If this model is correct, the probability of finding a χ²/dof this large or larger by chance is 0.120, so the fit of this model to this data set is also acceptable, according to the χ² test. In Figure 9, we show the value of χ in each of the 8320 phase-energy bins, for this fit. No patterns in the values of χ are evident as a function of phase or energy, which is what one would expect for a good fit.

Zoom In Zoom Out Reset image size

When the bolometric waveforms predicted by the models with two oval spots and three oval spots that best fit the NICER phase-channel data with 32 phases are compared with the 32-phase bolometric waveform constructed using the NICER data, the values of χ²/dof are considerably larger. When the bolometric waveform given by the model with two oval spots is compared with the 32-phase bolometric waveform constructed using the NICER data, the value of χ² is 40.6. There are 32-14-1-1 = 16 dof in this 32-phase bolometric data. The χ²/dof is therefore 40.6/16, which has a probability 6.3 × 10⁻⁴ of occurring by chance. This probability is low enough that it indicates that this model provides a description of this data set that is at least incomplete. When the bolometric waveform given by the model with three oval spots is compared with the 32-phase bolometric waveform constructed using the NICER data, the value of χ² is 41.7. There are 32-19-1-1 = 11 dof in this 32-phase bolometric data. The χ²/dof is therefore 41.7/11, which has a probability 2 × 10⁻⁵ of occurring by chance. This probability again indicates that this model provides a description of this data set that is at least incomplete.

Noting these low probabilities, we performed exploratory fits specifically designed to minimize the value of χ² obtained when the 32-phase bolometric waveforms predicted by these two energy-resolved waveform models are compared with the 32-phase bolometric waveform constructed from the NICER data. The probabilities of the models given by these fits were not significantly larger than the values given by the best fits of these models to the full phase-channel data.

Next, we constructed 64-phase bolometric waveforms using the NICER data and compared these with the 64-phase bolometric waveforms predicted by the best fit of the 32-phase energy-resolved waveform model with two oval spots to the 32-phase energy-resolved NICER waveform data, re-fitting the phase shift and the number of phase-independent background counts in each energy channel. The resulting minimum value of χ²/dof is 58.64/48, which has a chance probability of 0.140, indicating that this model provides an acceptable description of the 64-phase bolometric waveform data. When the predictions of this model for the 64-phase energy-resolved waveform are compared with 64-phase energy-resolved NICER data, the resulting χ²/dof is 16363.9/16365, which has a probability of 0.501, much higher than the probability we found when we divided the data into only 32 phase bins, and indicating that this energy-resolved waveform model also provides an acceptable description of the energy-resolved NICER waveform data with 64 phase bins.

Finally, we compared the 64-phase bolometric waveforms constructed using the NICER data with the 64-phase bolometric waveforms predicted by the energy-resolved waveform model with three oval spots that best fits the 32-phase energy-resolved NICER waveform data (again re-fitting the phase shift and the number of phase-independent counts in each energy channel). The results are shown in Figure 10. The minimum value of the bolometric χ²/dof is 59.6/43, which has a probability of 0.0473, indicating that this 64-phase model also provides an acceptable description of the 64-phase bolometric waveform data. When the predictions of this best-fit model for the energy-resolved waveform with 64 phase bins are compared with the 64-bin energy-resolved NICER waveform data, the resulting χ²/dof is 16347.5/16360, which has a probability of 0.526, again much higher than the probability that we found when we divided the data into only 32 phase bins, and indicating that this energy-resolved waveform model provides an acceptable description of the energy-resolved NICER waveform data with 64 phase bins.

Zoom In Zoom Out Reset image size

These results indicate that our best-fit models with two and three oval spots provide good descriptions of the NICER waveform data at high phase resolutions, and that the radius and mass estimates inferred from them are therefore credible. Why the bolometric waveforms given by the models that best fit the 32- and 64-phase energy-resolved NICER waveform data differ from the 32-phase bolometric waveforms constructed from the NICER data but agree with 64-phase bolometric waveforms is unclear, but is most likely due to temporal fluctuations in the NICER data that are not yet understood. This question deserves further study.

5.5.1. Locations of the Hot Spots

Figure 11 compares the best-fit locations and shapes of the hot spots obtained by fitting the models with two and three oval spots to the NICER waveform data. Both hot spots in the model with two oval spots and all three hot spots in the model with three oval hot spots are located in the southern rotational hemisphere, well away from the sightline to the observer (recall our convention that the sightline to the observer defines the northern hemisphere). Fits that force any of the spots to be in the same hemisphere as the sightline to the observer are strongly disfavored. As we noted earlier, when we discussed the phase-channel χ² values for the fits of both models to the NICER waveform data, these spot locations provide formally excellent fits.

Zoom In Zoom Out Reset image size

The fundamental reason spot locations in the southern hemisphere are favored appears to be that the modulation fraction of the waveform observed by NICER is high and its harmonic content is substantial. Both waveform properties favor spot locations in the southern hemisphere, because when the spots are visible for only a small fraction of a rotational cycle, both the modulation amplitude of the waveform and the strengths of the higher harmonics of the spin frequency are larger. Spots located in the northern hemisphere are visible for most of each rotational cycle and therefore produce waveforms with smaller modulation amplitudes and weaker high harmonics. In contrast, radio observations as well as joint fitting of radio and γ-ray profiles suggest geometries in which one of the two radio pulses comes from the northern hemisphere (Johnson et al. 2014). However, these studies assume a star-centered dipolar field, whereas the NICER results suggest that the magnetic field of PSR J0030+0451 is not a centered dipole (see also Lockhart et al. 2019). Although the locations of the soft, thermal X-ray emission and the nonthermal radio emission need not be the same, this difference should be investigated further.

5.5.2. Colatitude of the Observer

5.5.3. Spectrum of the Emission from the Hot Spots

An independent estimate of the total X-ray emission from PSR J0030+0451 can be made using observations made with XMM-Newton (Bogdanov & Grindlay 2009). Although the XMM-Newton EPIC MOS1/2 observations have substantially fewer source counts, they have a much lower background in the point-source spatial extraction region compared to the non-imaging NICER data, which means that to first order all of the counts detected by XMM-Newton come from the star rather than from unassociated sources. If the phase-integrated data from XMM-Newton are fit using a two-temperature non-magnetic hydrogen atmosphere model (as in Bogdanov & Grindlay 2009) and the predictions of the model are folded through the NICER response matrix, then we obtain an estimate for the total number of NICER counts from the star and the spectral shape that we expect to see in NICER. This predicted spectral flux can be compared with the total spectral flux expected from all the hot spots in our best-fit models of the pulse waveform emission. Figure 12 shows this comparison for our models with two oval spots and three oval spots. In both cases, the combined spectral flux expected from the hot spots at low photon energies falls short of the spectral flux predicted by the XMM-Newton observations.

Zoom In Zoom Out Reset image size

A first thought would be that the missing emission might be unpulsed thermal emission from a substantial fraction of the stellar surface. There are, however, several serious difficulties with such an interpretation. First, the emission would have to be almost exactly axisymmetric around the stellar rotation axis in order to avoid generating detectable pulsed emission. This forces one to consider emission patterns that are highly tuned: filled circular emitting regions around one or both rotation poles, annular emitting regions centered around one or both rotation axes, or some combination of these would be required to avoid producing detectable flux modulation. Second, if the emission is thermal, the total area of the axisymmetric emitting region(s) would have to be a very small fraction of the stellar surface.

We illustrate these difficulties by an example. In order to make up the observed deficit in the flux at low energies and not produce detectable modulation, a circular spot centered on the north rotational pole (the pole nearer the observer), would have to have an angular radius of just 0.075 rad, assuming thermal emission at the best-fit effective temperature of kT_eff ≈ 0.075 keV. Such a spot would also have to be very nearly circular and centered on the north rotational pole: a deviation of more than ∼0.01 rad would cause a flux modulation that would be inconsistent with the observed waveform. Because of its smaller projected area, a circular spot centered on the south rotational pole would have to have an angular radius ∼0.7 rad, about 10 times larger than a spot at the north rotational pole, but would have to be even more precisely circular and centered than a spot at the north rotational pole: a deviation of more than ∼0.005 rad would cause a flux modulation inconsistent with the observed waveform.

Other thermal emission geometries with the same projected area are theoretically possible (e.g., one or more very thin axisymmetric annuli), but these seem to us even more fine-tuned. The spectrum could also be produced by processes in which the peak in the spectrum required to make up the deficit does not reflect the temperature of the emitting plasma, but is instead created by a steep decrease of the optical depth of the emitting region with increasing photon energy. High-harmonic cyclotron emission in a region where the optical depth to cyclotron absorption decreases with increasing photon energy is one such possibility (Psaltis & Lamb 1999), but usually requires a population of higher-energy electrons, in which case the optically thick (Rayleigh–Jeans) portion of the spectrum would have to have a very small area, in order to avoid producing too much emission.

It does appear worthwhile to explore other possibilities. For example, the spot pattern is clearly not that of a centered dipole. If this means that the total magnetic field at the surface of the star has significant multipolar components, and is therefore much stronger than the nominal dipole component of ∼2 × 10⁸ G, then perhaps cyclotron emission could contribute to the missing emission.