NEUTRON STAR MASS–RADIUS CONSTRAINTS USING EVOLUTIONARY OPTIMIZATION

We computed pulse profiles for a set of nine test models with the parameter values listed in Table 1. The properties that were kept the same for all test models are the spot's temperature, the phase offset ϕ, the observer's inclination angle i, and the distance to the star d. The spot emission model is a 2 keV blackbody (in the frame comoving with the neutron star's surface) with a limb-darkening atmosphere, approximated by the Hopf function (Chandrasekhar 1960), appropriate for a Thompson-scattered atmosphere. Since for real data we would independently have the spot temperature at infinity instead of in the star's comoving frame, this parameter should be allowed to vary within a narrow range (explored in Section 3.7), where the appropriate range would be determined by the temperature at infinity. The model assumes that any emission from the surface of the neutron star outside the spot is negligible. The observer's inclination angle i = 60° and the phase offset ϕ = 0 for all cases. The computed pulse profiles are normalized to an average flux of 1, so that the distance d to the star does not affect the pulse profiles.

Table 1.  Summary of Test Models

Name	Description	Shape	νspin	M	R	i	θ	ϕ	M/R	$\sin i\sin \theta$	$\cos i\cos \theta$	Amp	β
			(Hz)	(${\text{}}{M}_{\odot }$)	(km)	(deg)	(deg)
A	Fiducial	Sphere	600	1.6	12	60	20	0	0.197	0.296	0.470	0.373	0.057
AO	Oblate	Oblate	600	1.6	12	60	20	0	0.197	0.296	0.470	0.373	0.057
θ37	θ = 37°	Sphere	600	1.6	12	60	37	0	0.197	0.521	0.399	0.720	0.101
θ37O	θ = 37°	Oblate	600	1.6	12	60	37	0	0.197	0.521	0.399	0.720	0.101
θ60	θ = 60°	Sphere	600	1.6	12	60	60	0	0.197	0.750	0.250	1.305	0.145
(M/R)hi	High M/R	Sphere	600	1.68	11	60	20.8	0	0.226	0.308	0.467	0.350	0.057
(M/R)lo	Low M/R	Sphere	600	1.35	13	60	19.8	0	0.153	0.293	0.471	0.423	0.057
βhi(M/R)hi	High β, high M/R	Sphere	600	1.68	11	60	64	0	0.226	0.778	0.219	1.235	0.145
ν400	Low νspin	Sphere	400	1.60	12	60	20	0	0.197	0.296	0.470	0.373	0.038

A machine-readable version of the table is available.

Download table as:  DataTypeset image

The parameters that were changed for different test models are M, R, θ, the star's shape (spherical or oblate), and ν_spin. For oblate neutron star models, R is defined to be the radius of the neutron star at the spot. The formalism for the oblate model is detailed in Morsink et al. (2007) and does not require the addition of any extra parameters.

The parameters for the fiducial model, model A, were chosen to be representative of the masses and radii of accreting millisecond neutron stars. Our fiducial mass M = 1.6 ${\text{}}{M}_{\odot }$ is larger than the mass M = 1.4 ${\text{}}{M}_{\odot }$ typically measured for slowly rotating radio pulsars, since we expect that the neutron star has been spun up by and gained mass from accretion. The radius of 12 km is consistent with other radius estimations (Leahy 2004; Steiner et al. 2010). Rapid rotation with ν_spin ≥ 550 Hz is seen in burst oscillations for at least seven neutron stars that exhibit Type I X-ray bursts (Watts 2012). We have chosen ν_spin = 600 Hz for the fiducial model since it is representative of these rapid rotators. The angles were chosen to provide a pulse amplitude (Amp = 0.373) comparable to the largest pulse amplitudes seen, in order to reduce the parameter degeneracy. In particular, the neutron star 4U 1636–536 has rms pulse amplitudes as high as 0.25 after subtracting the pre-burst emission (Strohmayer et al. 1998; Galloway et al. 2008), where Amp = $\sqrt{2}\,\mathrm{rms}$. The full set of theoretical parameters is given in the row labeled "A" in Table 1. The resulting pulse profiles in two energy bands are shown in Figure 1 with black and red solid curves labeled "True." Our fiducial model is similar to the model described as "low inclination" by Lo et al. (2013), except that they consider a slightly smaller R and a slower ν_spin. As will be shown in Section 3.5, small changes in R do not qualitatively change the results of this paper. Our choice of a faster ν_spin improves the accuracy of the fits, but the faster ν_spin is also more appropriate for the neutron stars that we hope to apply this method to.

Zoom In Zoom Out Reset image size

Theoretical pulse profiles were calculated for the nine different test models. Each theoretical pulse profile was converted to a set of 20 synthetic observations by adding noise from a Poisson distribution to the pulse profile. The standard case (a low-noise model) assumes 25,000 photon counts per phase bin and no background count rate, as if the background were negligible in comparison with the hot spot emission. By not including a background, we are underestimating the Poisson fluctuations and thus underestimating the error bars. This is done so we can test the suitability of evolutionary optimization on the best-possible quality of synthetic data. For a variation on the model θ₆₀, a higher noise level was used: θ₆₀C₆₂₅₀ assumes 6250 photon counts per phase bin with no background. In total, there are 10 sets of 20 simulated observed pulse shapes.

Each of the 200 simulated pulse profiles was then fit using the Ferret algorithm (described in the next subsection), which searches for different values of the free parameters in order to minimize the ${\chi }^{2}$ fit statistic. Five parameters were allowed to vary within physical ranges: 1.0 ${\text{}}{M}_{\odot }$ ≤ M ≤ 2.5 ${\text{}}{M}_{\odot }$, 6.0 km ≤ R ≤ 16.0 km, 0° ≤ (i, θ) ≤ 90°, and ϕ defined cyclically from 0 to 2π. The star's shape (spherical or oblate) was fixed to be the same in the fitting procedure as in the theoretical waveform. The temperature is kept fixed at 2 keV since it has been shown by Lo et al. (2013) that observations in multiple energy bands allow for a good determination of the gravitationally redshifted spot temperature. We do not allow the spot size to vary, as explained in the previous subsection. For applications to real data, it would not be difficult to add variations in spot size (and shape) and temperature and to use unnormalized fluxes to allow for a distance measurement. However, the addition of the extra free parameters for the number of synthetic waveforms considered in this study would be impractical. We carried out preliminary trials (summarized in Section 3.7), allowing these parameters to vary, but did not find much change in our final results.

Since there are two energy bands with 32 phase bins each and five free parameters, there are 59 degrees of freedom (dof). The fitting yielded best-fit model parameters and also allowed computation of confidence regions in the M–R plane, from which we can assess the effect of various parameters on the uncertainties in parameter determination.

Evolutionary optimization algorithms provide a useful but heuristic approach to search through large parameter spaces, to optimize the fit of a model to data. Such algorithms use principles inspired by biology to evolve a population of candidate parameter sets, over many generations, toward an optimal solution. The heuristic nature of the search aims to sample the parameter space efficiently, but not completely, since exhaustive search is not practical for problems involving many parameters. Therefore, it cannot be guaranteed that the true, globally optimal solution has been found, although this would also be true for any other optimization algorithm, besides an exhaustive search. Evolutionary algorithms have been well studied and found to be useful in many fields of science and engineering.

The oldest, and most commonly known, type of evolutionary optimizer is the genetic algorithm (Holland 1975; Goldberg 1989, 2002). Classic genetic algorithms encode their search parameters on a (usually) binary string (genotype), which is decoded into a model (phenotype). A population of candidate parameter sets is normally initialized as a set of random bit strings, which are expressed as a model, and in turn evaluated by a fitness function. This information is used to probabilistically select good parameter sets (individuals) to propagate to the next generation, in analogy to the "survival of the fittest" principle of Darwinian evolution. Parameter sets undergo random bitwise mutations to help perturb solutions into previously unexplored regions of parameter space. Information is shared between individuals by means of a crossover operator that cuts a pair of bit strings at a random position and recombines them into a new "offspring" configuration, in analogy with sexual reproduction. This process of evaluation, selection, mutation, and selection is performed iteratively, over many generations, until a convergence criterion terminates the search. The genetic algorithm can be viewed as a directed stochastic search, which makes use of random noise to evade local minima, while a mostly deterministic selection operator pushes solutions toward an optimal solution.

In this paper, we used versions 5.3–5.5 of Ferret from the Qubist Global Optimization Toolbox for MATLAB (Fiege 2010; Rogers & Fiege 2011), a commercially available software package, to find the best fits to our synthetic data. This is an alternative to the Markov chain Monte Carlo approach in Lo et al. (2013) to fit pulse profiles. Ferret's development began in 2002, as a variant of the multi-objective genetic algorithm, which made use of a real-valued parameter encoding and real-valued mutation and selection operators, rather than the traditional binary encoding discussed above. Numerous other features have been added to Ferret since then, which go beyond the usual genetic algorithm paradigm. Most notably, the code contains a unique linkage-learning algorithm, which detects a certain type of nonlinearity (linkage) between parameters, with the goal of dividing a large parameter space into several smaller and nearly independent subspaces, thus greatly simplifying the problem. This capability is discussed at length in the software user's manual (Fiege 2010). Ferret also contains an algorithm that allows several of its most important control parameters, including the typical strength of mutation and crossover events, to be automatically adapted and optimized during a run. This auto-adaptation capability is similar in spirit to the self-adaptation used in Evolution Strategies (ES) codes, although ES codes do not typically employ a crossover operator. Ferret does not adhere to the strict definition of a genetic algorithm, due to these enhancements and others, but the code still remains closer to the genetic algorithm paradigm than to others within the family of evolutionary optimizers.

As an example of Ferret's inner workings, consider the parameters i and θ. Ferret selects sets of values from the allowed region of parameter space (0 ≤ i, θ ≤ 90°) for the population of a generation and fits pulse profile models. It then keeps the best-fitting parameter sets and selects new sets to explore more of the parameter space. Within a few generations, Ferret discovers the degeneracy between i and θ, since their values can be swapped with minimal impact on the ${\chi }^{2}$ fit. A more detailed explanation of Ferret in relation to pulse profile modeling can be found in Chapter 4 of Stevens (2013). Evolutionary optimization algorithms have also been used in gravitational lensing (Rogers & Fiege 2011, 2012), medical physics (Fiege et al. 2011), star formation (Franzmann 2014), and X-ray spectral fitting (Rogers et al. 2015) applications.

Our problem is quite easy for Ferret, which consistently finds the global minimum χ² fit statistic within a few generations. We note that the true global minimum is known, since our results are based on artificial data tests. After finding the χ² minimum, the algorithm was allowed to run for approximately 100 generations, as Ferret accumulated solutions within the confidence region. No specific convergence criterion was implemented; rather, we terminated the run manually when it was evident from the software's graphical user interface that no further improvements were being made to the minimum χ², and enough solutions had been accumulated within the confidence region to make contours maps of acceptable quality. Ferret was executed in MATLAB R2011a on 12-core AMD Linux servers (dual socket Opteron 2439 SE, with 32Gb RAM, running Red Hat version 4.1.2), using Ferret's built-in parallel computing features. Each fit took approximately 8–20 hr.

Twenty realizations of the Poisson noise were made to simulate observations of the fiducial model A. The resulting flux values with error bars for the first Poisson realization of the data are shown as black and red points labeled "A1" on Figure 1. The data points with Poisson errors were used as input to Ferret, which searched for the best fit to the A1 data set by minimizing ${\chi }^{2}$. In the case of the A1 data set, the pulse profiles that best fit the data are displayed with blue dashed and dotted curves in Figure 1. The parameters corresponding to the best fit of the A1 data set (${\chi }^{2}$ = 58.2) are shown in the row labeled "1" in Table 3. For reference we have included the derived quantities M/R, $\sin i\sin \theta$, $\cos i\cos \theta$, Amp, and β in this table.

Table 3.  Summary of Best Fits for Model A

Model	${\chi }^{2}$	M	R	i	θ	ϕ	M/R	$\sin i\sin \theta$	$\cos i\cos \theta$	Amp	β
	59 dof	(${\text{}}{M}_{\odot }$)	(km)	(deg)	(deg)						×10−2
True		1.600	12.00	60.0	20.0	0.0000	0.1969	0.296	0.470	0.373	5.741
1	58.2	1.575	11.97	48.7	25.0	0.0071	0.1943	0.318	0.598	0.347	6.120
2	60.8	1.517	12.39	19.9	57.3	−0.0011	0.1808	0.286	0.508	0.362	5.577
3	58.6	1.344	9.39	39.4	34.4	0.0021	0.2114	0.358	0.638	0.357	5.565
4	71.7	1.532	10.70	21.7	59.8	0.0000	0.2114	0.319	0.467	0.383	5.658
5	63.9	1.666	11.46	71.8	17.6	−0.0062	0.2147	0.287	0.297	0.426	5.471
6	56.7	1.445	10.26	36.7	36.0	0.0064	0.2080	0.351	0.649	0.350	5.932
7	64.6	1.469	10.81	37.2	34.4	0.0074	0.2007	0.342	0.657	0.345	6.007
8	53.9	1.694	12.39	68.7	17.3	−0.0047	0.2019	0.277	0.347	0.404	5.584
9	42.3	1.417	11.76	42.8	27.5	0.0044	0.1779	0.313	0.651	0.338	5.773
10	67.0	1.510	10.32	21.9	60.6	−0.0007	0.2161	0.324	0.456	0.388	5.583
11	57.1	1.581	11.29	21.7	58.1	0.0026	0.2068	0.314	0.491	0.372	5.828
12	68.9	1.433	12.08	29.1	39.9	0.0062	0.1752	0.311	0.671	0.331	5.869
13	58.4	1.458	10.88	35.6	35.8	0.0075	0.1979	0.341	0.659	0.345	5.994
14	63.7	1.702	13.86	15.9	67.6	−0.0057	0.1813	0.254	0.366	0.390	5.541
15	58.7	1.447	9.59	38.1	36.9	0.0067	0.2228	0.371	0.629	0.360	6.006
16	45.5	1.746	13.31	18.3	63.1	0.0015	0.1937	0.280	0.429	0.375	5.983
17	51.8	1.435	12.96	30.7	36.3	0.0086	0.1635	0.302	0.693	0.323	6.009
18	53.1	1.488	10.57	38.4	34.2	0.0076	0.2079	0.349	0.648	0.348	6.076
19	72.0	1.591	12.22	19.9	59.0	0.0003	0.1923	0.292	0.485	0.367	5.727
20	58.1	1.438	9.68	38.5	36.1	0.0060	0.2195	0.367	0.632	0.359	5.962
Ave	59.3	1.524	11.39	34.7	41.8	0.0028	0.1989	0.318	0.549	0.363	5.813
Std	7.6	0.107	1.23	15.0	15.2	0.0046	0.0159	0.032	0.120	0.025	0.203

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

The best-fit model (star) and contours of constant ${\chi }^{2}={\chi }_{\min }^{2}+{\rm{\Delta }}{\chi }^{2}$ (black curves) for the A1 data set are shown in the M–R plane in Figure 1. The contours correspond to values of ${\rm{\Delta }}{\chi }^{2}$ = 2.3, 6.2, and 11.8, corresponding to 1, 2, and 3${\sigma }_{\mathrm{degen}}$ confidence levels, respectively, in mass–radius space for two free parameters, M and R. This plot uses the profile likelihood (Murphy & Van Der Vaart 2000) to eliminate the nuisance parameters i, θ, and ϕ; a grid is defined in M and R, and for each (M, R) grid point, Ferret finds the lowest ${\chi }^{2}$ within each grid cell, allowing for any values of i, θ, and ϕ. The 1${\sigma }_{\mathrm{degen}}$ confidence region (calculated by finding the maximum and minimum values on the 1${\sigma }_{\mathrm{degen}}$ contour and dividing by 2) gives an approximate 1${\sigma }_{\mathrm{degen}}$ error of 1.9 km for the radius and a 0.18 ${\text{}}{M}_{\odot }$ error for the mass. Note that since the contours are not ellipses, these are only approximate 1${\sigma }_{\mathrm{degen}}$ limits.

Ferret was used to fit each of the 20 Poisson realizations of model A. The independent best-fit results for the 20 different realizations are shown in Table 3. The average and standard deviation for each parameter are displayed at the bottom of Table 3 and also appear in the second line of Table 2. For the individual angles i and θ, it can be seen from Table 3 that the determinations of these angles are very poor. This is due to the already well-known degeneracy (for spherical stars; Poutanen & Gierliński 2003) that occurs since the equations for light bending and the Doppler effect only depend on the combinations $\sin i\sin \theta$ and $\cos i\cos \theta$ and not on their individual values. The $i-\theta$ degeneracy can be seen in Figure 2, where the normalized low-energy band pulse profiles for model A and another model with i and θ swapped (both labeled "Sphere") are indistinguishable. For many of the Poisson realizations, the best-fit values for i and θ shown in Table 3 are swapped from their true values, and as a result, the average and standard deviations for the individual angles are really not meaningful, except to illustrate that it is the trigonometric combinations of the angles that can be reliably determined. However, since there are independent methods for constraining i and θ through optical (Wang et al. 2013) or gamma-ray (Venter et al. 2012) observations, it is still useful to discuss these two angles separately.

Zoom In Zoom Out Reset image size

To ensure that a suitable number of independent realizations of the data were made, a histogram of the best-fit radii is plotted in Figure 3. The resulting histogram is approximately Gaussian (70% are within 1 standard deviation, 95% within 2 standard deviations, 100% within 3 standard deviations), indicating that 20 trials are sufficient to illustrate general trends. A scatter plot, also shown in Figure 3, shows the values of mass and radius (red crosses) for the 20 different random realizations of the data. The large black cross indicates the average and standard deviation values for the mass and radius fits.

Zoom In Zoom Out Reset image size

As previously stated, parameter degeneracy plays a major role in all aspects of uncertainty in our results, even the standard deviation. By testing 20 different realizations of the data for each model (a relatively small number), we still find that some models converge on the true parameter values while others do not. The models with stronger parameter degeneracy have larger standard deviations on the parameters. This provides an initial assessment of which types of models and pulse profile shapes have stronger degeneracies. Even with very small error bars on the pulse profile, there would still be relatively large standard deviations over many models due to the inherent degeneracies.

In this section, we examine the effect of changing the hot spot's colatitude θ on the accuracy and precision of the pulse profile fits. By changing θ while keeping M, R, and i constant, we alter the projected velocity β and the approximate pulse amplitude Amp (see Table 1). Models A (with θ = 20°), θ₃₇, and θ₆₀ have increasing values of θ, β, and Amp. Due to the i–θ degeneracy, this is also equivalent to keeping θ fixed and varying i.

The θ₃₇1 pulse profile (i.e., for the first Poisson realization of the θ₃₇ model) is shown in the left panel of Figure 4, and the constant ${\chi }^{2}$ contours are shown in the right panel. This case is one of the "outlier" pulse shapes (of the 20 Poisson realizations) with a very large best-fit R that only includes the true R at the 3${\sigma }_{\mathrm{degen}}$ level. The 1${\sigma }_{\mathrm{degen}}$ error regions for the θ₃₇1 model (from the right panel of Figure 4) are smaller than the A1 model, close to 1.1 km for R and 0.11 ${\text{}}{M}_{\odot }$ for M. The 1${\sigma }_{\mathrm{degen}}$ regions for this particular Poisson realization are somewhat larger than the standard deviation 1σ computed for the ensemble of θ₃₇ models.

Zoom In Zoom Out Reset image size

Similarly, the θ₆₀1 pulse profile (first Poisson realization of the θ₆₀ model) is shown in the left panel of Figure 5, and the constant ${\chi }^{2}$ contours are shown in the right panel. The 1${\sigma }_{\mathrm{degen}}$ error regions for the θ₆₀1 model are 0.7 km for R and 0.06 ${\text{}}{M}_{\odot }$ for M, larger than the ensemble standard deviation 1σ by a factor of about 6.

Zoom In Zoom Out Reset image size

We find a strong trend as θ increases: the average M and R fit values are closer to the true values, and the standard deviation is smaller. For example, σ in M decreases from 0.11 to 0.08 to 0.01 M_⊙ for models A, θ₃₇, and θ₆₀. This improvement in fitting accuracy is due to models θ₃₇ and θ₆₀ having a larger Amp and β than model A. This is consistent with the general trend seen by Lo et al. (2013) in individual model fits.

During an X-ray burst, θ is expected to change. This set of models roughly approximates this process, which has been explored in detail by Mahmoodifar & Strohmayer (2016). If a series of different burst oscillation pulse shapes are found for a neutron star during the same burst, future work could include a simultaneous multi-epoch fit, allowing θ (and spot size ρ) to be dependent on epoch.

Rotation alters the shape of a neutron star, making it an oblate spheroid. Although the change in shape is small for stars spinning at rates seen in accreting systems, the alteration in the star's shape changes the positions on the neutron star's surface for which photons can reach the observer (Morsink et al. 2007), leading to large changes in the pulse profile. As an example, a 1.6 M_⊙ neutron star with an interior given by the APR equation of state (Akmal et al. 1998) has an equatorial radius of 11.7 km and a ratio of polar to equatorial radii of 0.93. This oblate geometry has a larger effect on the pulse profile than other effects due to rotation, such as frame dragging, as has been discussed in detail in, e.g., Morsink et al. (2007).

In this section, we investigate the effect of the oblate shape on the accuracy and precision of pulse profile fitting and reproducing the input parameters. To do this, we construct oblate versions of two models, A and θ₃₇. The oblate versions, designated with the letter "O," are constructed so that the radius of the star at the location of the spot is the same as for the corresponding spherical model. This means that the values of the star's compactness M/R and the projected velocity β are the same for the spherical and oblate models. As a result of this definition for the radius, the star's equatorial radius is larger than that listed in Table 1. The general trend is for the oblate neutron star's pulse profile to have a smaller pulse amplitude than the spherical neutron star with the same parameters (when i and θ are in the same hemisphere), since the visibility condition makes it easier to see the far side of the neutron star when the star is oblate (Cadeau et al. 2007; Morsink et al. 2007). In cases where the spot is visible for all phases (as in both models A and θ₃₇), consider the emission when the spot and the observer have the same azimuthal angle. At this moment, the light from the spherical star is emitted close to the normal to the surface. For the oblate star, the light is emitted in the same direction in space, but it is at an angle farther from the surface's normal, due to the tilt of the surface. Since the intensity of light is proportional to the cosine of the angle between the original direction of emission and the normal to the surface, the spherical star's spot appears brighter at this phase. The opposite is true when the spot is on the opposite side of the star from the observer. In this case the light is emitted close to the tangent to the surface for the spherical star and closer to the normal for the oblate star, making the spot appear dimmer for the spherical star at this phase. The overall effect is to create a less modulated pulse profile for the oblate star.

In model AO, the neutron star's equatorial radius is 12.7 km, while the radius at the spot's latitude is 12 km. As a comparison with the spherical model, the low-energy band pulse profiles for both AO and A are plotted as black dot-dashed and solid curves, respectively, in Figure 2. The high-energy band (not shown for clarity) has a similar decrease in modulation. In model θ₃₇O the equatorial radius is 12.5 km and has a lower pulse amplitude than model θ₃₇. Since the pulse profiles and M–R confidence regions for models AO and θ₃₇O look quite similar to those for models A and θ₃₇, they are not shown here. The mean and standard deviation of the best-fit parameters for the 20 realizations for each oblate model are given in Table 2. For most parameters, the σ and accuracies are smaller for the oblate case than for the fiducial spherical case (we also found that the 1${\sigma }_{\mathrm{degen}}$ regions were smaller by about a half). A similar improvement in the precision of the results was also noted by Miller & Lamb (2015); however, since they only considered one Poisson realization of the data, they could not rule out the possibility that this was due to a statistical fluctuation. In our case, since we are comparing a sample of Poisson realizations for each model, the increase in precision and accuracy is most likely due to the properties of the oblate model pulse profiles.

The improvement in the accuracy and precision in the determination of most of the parameters for oblate models is most likely due to the partial lifting of the degeneracy between i and θ. For an oblate star, the direction in which the normal to the surface points depends on the shape of the star, which is a function of θ, but is independent of i. This introduces small differences between the normalized pulse profiles for models with i and θ switched, while for spherical stars, the normalized pulse profiles are the same when the angles are switched. This is illustrated in Figure 2, where the two oblate models representing the swapped inclination and spot angles, shown with black dot-dashed and red dotted curves, are clearly different from each other. However, since they are still fairly similar to each other, there is still a partial degeneracy when the angles are swapped.

There is a similar but much smaller lifting of degeneracy for large hot spots on spherical stars. A large spot extends over a range of θ values, while the observer's inclination i is just one fixed value, so swapping the two angles if the spot is large will not yield the same pulse shape. However, the magnitude of the effect is much smaller than the magnitude of the change that occurs when the angles are swapped on an oblate star.

We also tested the effect of using the wrong shape model by using the oblate data corresponding to the 20 AO models as input and fitting them with pulse shapes for spherical stars. In this case, ${\chi }^{2}$ increased by a small amount, but the best-fit values for M and R became very inaccurate. For the set of fits to the Poisson realizations, the average fit returned ${\chi }^{2}$ = 64.7, M = 1.91 ${\text{}}{M}_{\odot }$, and R = 12.3 km (which should be compared with the values in the fourth row of Table 2). The average M and R fits are accurate to 19.4% and 2.5%, respectively. Similar results were found by Miller & Lamb (2015). In practice, this is not a problem, since we know that the stars are actually oblate, so the correct shape can be used. The spherical shape model is computationally somewhat cheaper; hence, it is used for the other models tested in this paper. Furthermore, the oblate shape model used is still an approximation. It would be worthwhile, in future research, to compare the slightly different shape models that have been used in this paper with other models used by other groups (e.g., Bauböck et al. 2013; Miller & Lamb 2015).

The average photon count rate per phase bin of 25,000 is the expected best-case scenario, assuming typical burst flux, large detector effective area, and long observation time (Feroci et al. 2012), and similar to the count rate used in other pulse profile analyses (Lo et al. 2013; Psaltis et al. 2014). The pulse profile is noisier when fewer counts are detected. We tested the effect of statistical noise on model θ₆₀ by generating different Poisson realizations, assuming a reduced count rate, of the same theoretical pulse profile. Comparing model θ₆₀ with the reduced count-rate model, θ₆₀C₆₂₅₀, allows us to explicitly compare the degeneracy-related uncertainty versus error. The average counts per bin are 25,000 and 6250 for models θ₆₀ and θ₆₀C₆₂₅₀, respectively.

The pulse profile for model θ₆₀C₆₂₅₀ looks very much like that for model θ₆₀, but with larger errors, so it is not shown. The ${\chi }^{2}$ contours for model θ₆₀C₆₂₅₀1 are shown in Figure 6. The sizes of the 1${\sigma }_{\mathrm{degen}}$ contours in Figures 5 and 6 are very similar, suggesting that the confidence contours for one particular mock observation are mainly dominated by the parameter degeneracy.

Zoom In Zoom Out Reset image size

From Table 2, the standard deviations of the parameters for the low count-rate model are ≃2–3 times as large as those of model θ₆₀. The factor of 2 is expected purely on statistics (from the factor of 4 reduction in counts). The actual degradation is somewhat worse, likely because partial degeneracy between the parameters makes the parameter determination degrade more rapidly than the errors. Model θ₆₀C₆₂₅₀ still gives well-determined parameters even at the lower count-rate value. The accuracies in M and R determination are 2% and 4%, respectively, for the low count-rate case.

Changes in the compactness ratio (M/R) are expected to impact the quality of the fits, since the compactness affects the pulse amplitude Amp (evident in Equation (1)). Increasing the compactness decreases Amp, which might be expected to decrease the precision and/or accuracy of the fits. To test this effect, we generated models with differing values for the compactness ratio but the same values of the projected spot velocity β.

We created two models similar to model A with a larger and smaller compactness ratio, (M/R)_hi and (M/R)_lo, and generated a set of pulse shapes with Poisson noise. The results of the fits can be seen in Table 2. All three models with the same value of β have similar accuracy and precision for most of the parameters. For instance, the accuracy of determining the mass ranges from 4% to 5%, while the precision ranges from 6% to 7%, so the change in compactness does not greatly affect the fits. However, while there is little change in the precision of the radius measurement (9%–10%), in the case of the low compactness model the accuracy was improved.

As a further test on the effect of compactness, model β_hi(M/R)_hi was created to compare with the high spot colatitude model θ₆₀. The parameters were the same except that θ was adjusted to give the same β as model θ₆₀. This also necessitated that Amp was somewhat lower than model θ₆₀ (see Equations (1) and (2)). Previously for model θ₆₀, the accuracy and precision for both M and R were better than 1%; by fitting model β_hi(M/R)_hi, we found that increasing the compactness degrades both the accuracy and precision in M and R, but only so that they range from 1% to 3%.

From these results, given a value of β, it appears that the accuracy and precision to which the other parameters can be determined do not strongly depend on the compactness ratio of the star. Thus, our fit results would not be significantly different if we had chosen a different M and R for the fiducial model A.

The last parameter we adjusted to be different from the fiducial model was ν_spin. We compared model A (ν_spin = 600 Hz) with a 400 Hz model, labeled ν₄₀₀. Model ν₄₀₀ has all parameters the same as model A except ν_spin, and hence a different projected spot velocity (β) (see Table 1). This model has parameter values very similar to the "low-inclination" model presented by Lo et al. (2013).

The fits results are summarized in Table 2. The precision (standard deviations) and accuracy (difference of mean and true values) are worse for model ν₄₀₀, as expected. For M, R, M/R, and Amp, the precisions are about twice as large for the 400 Hz model as for the standard 600 Hz model. However, β has the same absolute precision for the 400 Hz model as for the 600 Hz model.

Along with the sets of trials summarized in the previous subsections, we tested some of the model assumptions used in the fits. These tests include the assumptions we have made about the spectral emissivity, spot size, presence of a constant background flux, and the star's shape. Errors are not reported for the tests of this subsection since fits were only carried out on the A1 synthetic data set, the first Poisson realization of the model A theoretical curve.

We first tested the effect of using the wrong atmosphere model in the fits. To do this, we used the synthetic data corresponding to the A1 model (see Figure 1) as input to Ferret. The synthetic data were generated using the Hopf limb-darkening function, but Ferret was run assuming a perfectly isotropic blackbody spectrum. The algorithm was unable to find a good fit after 300 generations (about double the number of generations normally required), and the best fit had ${\chi }^{2}$ = 316.6 for 59 dof. The best-fit M and R were 1.21 ${\text{}}{M}_{\odot }$ and 11.3 km, which are quite inaccurate (24.4% for M, 5.8% for R). As was seen by Lo et al. (2013), one does not get a good fit using an incorrect atmosphere model.

We next tested fitting a variable spot size model to a pulse profile that was generated with an infinitesimally small spot. The infinitesimal size of the spot is a necessary assumption due to the large number of trials carried out in this work. In order to model a larger spot, the spot has to be cut into a number of segments, and separate computations of the deflection angles must be computed, increasing the run time linearly with the number of segments. The pulse shape for a large spot tends to have a lower pulse amplitude than a small spot centered at the same latitude, due to the effect of averaging over many latitudes. As a test of this assumption, we used the A1 data set as input to the Ferret algorithm, but added the angular radius ρ of the spot as a free parameter that was allowed to vary between 1° and 60°. The resulting best fit has ${\chi }^{2}$ = 59 (for 58 dof) and converged on an angular radius of ρ = 1°, the smallest angle allowed. The best-fit M and R for this case are 1.51 ${\text{}}{M}_{\odot }$ and 11.6 km. While these values are less accurate (5.6% and 3.3%, respectively) than the best-fit value shown in the first row of Table 3 corresponding to a fit with an infinitesimal spot, they are within the 1σ limits computed for the set of model A fits. Since using multiple spot segments leads to a significantly longer computation time, it was not feasible to use variable spot sizes in this paper. However, in future work when real data are being fitted, it will be necessary to allow the spot size to vary and to compute the observed flux from multiple spot segments.

We also tested the effect of allowing the temperature of the spot to be a free parameter in the fitting program. The theoretical model was constructed with a spot temperature of 2 keV (as measured in the comoving frame at the star's surface). The Ferret program was used to fit the data with the addition of a local temperature parameter that was allowed to vary between 1 and 3 keV. The best fit returned ${\chi }^{2}$ = 57.9 (for 58 dof), M = 1.67 ${\text{}}{M}_{\odot }$, R = 11.8 km, and a temperature of 2.1 keV. The accuracies in M and R are 4.8% and 1.7%, respectively. Although these M and R values are within the 1σ limits for this data set, the addition of the extra parameter does degrade the accuracy. This is most likely due to the degeneracy introduced since we measure the redshifted temperature with the light curves. For the case of real data, it would be important to allow the temperature to vary in the fits, since the temperature of the spot in the comoving frame is unknown. Furthermore, making use of multiple energy bands would improve the accuracy of the temperature measured at infinity.

The effect of an unknown background flux has been investigated in detail by Lo et al. (2013), and for comparison we test it for one simple case. The synthetic data set A1 was constructed without a background count rate. We used this as input and allowed Ferret to add two new parameters corresponding to a constant (in time) background flux in each energy band. These background values were allowed to vary between 0 and 1 (recall that our pulse profiles are normalized to 1). The resulting fit had ${\chi }^{2}$ = 57.7 (for 57 dof), and the best-fit M and R were 1.52 ${\text{}}{M}_{\odot }$ and 11.7 km. M and R were accurate to 5.0% and 2.5%, respectively. It is possible to add a more realistic background model by adding emission at a lower temperature from the rest of the star (Lo et al. 2013; Elshamouty et al. 2016) or by adding light scattered from the disk (Morsink & Leahy 2011), but at this point it is not obvious what the most realistic model would be. Real data are also likely to contain non-negligible emission from the surface of the neutron star outside the spot region and from the accretion disk.

The effect of these tests on the best-fit values of M and R are summarized in Table 4. These tests highlight the importance of having as realistic an atmosphere model as possible, since fitting with the wrong model drastically affected the quality and accuracy of the fits. Thus, this method provides a sensitive test for atmosphere models. Adding a free parameter for the temperature, unknown background flux, or spot size did not significantly detract from the accuracy in M and R, and our code was able to replicate the additional parameter quite well. Ideally, for real data we would allow these parameters to also vary in the fits, even if it is not feasible to do so in the present study.

Table 4.  Summary of Tests on A1 Data

Model	${\chi }^{2}$/dof	M	R	i	θ	ϕ	M/R	$\sin i\sin \theta$	$\cos i\cos \theta$	Amp	β
		(${\text{}}{M}_{\odot }$)	(km)	(deg)	(deg)	×10−3					×10−2
True		1.600	12.00	60.0	20.0	0.0000	0.1969	0.296	0.470	0.373	5.741
A1 best fit	58.2/59	1.575	11.97	48.7	25.0	0.0071	0.1943	0.318	0.598	0.347	6.120
Wrong atmosphere	316.6/59	1.210	11.29	40.2	39.2	0.0104	0.1583	0.408	0.592	0.495	7.004
Spot size varies	59.0/58	1.508	11.57	33.3	37.0	0.0090	0.1925	0.331	0.667	0.338	6.137
Temperature varies	57.9/58	1.668	11.81	23.1	55.1	0.0057	0.2086	0.322	0.526	0.364	6.266
Background	57.7/57	1.522	11.75	45.8	26.7	0.0066	0.1913	0.321	0.623	0.344	6.045

Download table as:  ASCIITypeset image