DETERMINING NEUTRON STAR PROPERTIES BY FITTING OBLATE-STAR WAVEFORM MODELS TO X-RAY BURST OSCILLATIONS

2.1.1. Hot Spot

Burst oscillations are thought to be produced by X-ray emission from a region on the surface of the star that is hotter than the rest of the stellar surface and is rotating at or near the rotation frequency of the star. Such a hotter region could be produced either by heating of part of the stellar surface by thermonuclear burning or because a disturbance in the outer layers of the star, such as a global surface mode, has made a localized region hotter (see Watts 2012).

Oscillations with the relatively high amplitudes (≳5%–10%) required to derive significant constraints on M and ${R}_{\mathrm{eq}}$ are probably produced predominantly by a single hotter region (see Lamb et al. 2009a). The dominant role of a single, localized hotter region is also indicated by the absence in the waveform of substantial second or higher harmonics of the fundamental oscillation frequency. Hence, regardless of whether the localized hotter region is produced by heating of the stellar surface by nuclear burning or by a global surface mode, the waveform it produces can be modeled using a single hot spot. Lo et al. (2013) have shown that the waveform produced by a hot spot is relatively insensitive to elongation of the spot in the east–west or north–south directions, so modeling the spot as a circular area is expected to be adequate.

An important question is whether oscillations observed during the rise of bursts or during burst peaks and tails are likely to be more useful in deriving constraints on M and ${R}_{\mathrm{eq}}$. The nuclear-powered emission from the surface of the star is expected to be highly localized during the first fraction of a second of an X-ray burst. If adequate constraints could be obtained using data collected only when the emission comes from a hot spot with an angular radius ${\rm{\Delta }}{\theta }_{\mathrm{spot}}\lesssim 10^\circ$, waveforms could be modeled assuming the hot spot is a point, and it would be unnecessary to include ${\rm{\Delta }}{\theta }_{\mathrm{spot}}$ as a parameter in fitting waveform models to waveform data. If, on the other hand, data collected when the emission is larger must be used, ${\rm{\Delta }}{\theta }_{\mathrm{spot}}$ must be included as a parameter.

Lo et al. (2013) carried out a comprehensive analysis of the likely usefulness of data taken during burst rises versus burst peaks and tails (see their Section 2.2). They found that analyses of oscillations observed during burst tails may provide estimates of M and ${R}_{\mathrm{eq}}$ with uncertainties comparable to or possibly even smaller than the uncertainties provided by analyses of oscillations observed during burst rises. Our approach in this paper is based on their findings, which we therefore summarize here.

The precisions of M and ${R}_{\mathrm{eq}}$ estimates are most sensitive to the star's rotation rate and the inclinations of the hot spot and the observer, because these quantities strongly affect the special relativistic Doppler boost and aberration caused by the line-of-sight component of the surface rotational velocity, and therefore affect the harmonic content of the observed waveform. The uncertainties of M and ${R}_{\mathrm{eq}}$ estimates also depend on the amplitude of the oscillation and the total number of counts, because these determine the size of the statistical fluctuations in measurements of the waveform variation.

Using the results of their extensive parameter estimation study, Lo et al. (2013) showed that, other things being equal, the uncertainties in M and ${R}_{\mathrm{eq}}$ estimates obtained by fitting waveform models to waveform data scale approximately as ${{\mathcal{R}}}^{-1}$ (see their Equation (1), Section 4.2.1, and Table 3), where

$$\begin{eqnarray}&&{\mathcal{R}}\equiv {N}_{\mathrm{osc}}/\sqrt{{N}_{\mathrm{tot}}}=1.4\;{f}_{\mathrm{rms}}\sqrt{{N}_{\mathrm{tot}}}.\end{eqnarray} \tag{ 1 }$$

Here ${N}_{\mathrm{osc}}=1.4\;{f}_{\mathrm{rms}}{N}_{\mathrm{tot}}$ is the number of counts collected from the oscillating component of the waveform during the observation, ${N}_{\mathrm{tot}}$ is the total number of counts collected, and ${f}_{\mathrm{rms}}$ is the average fractional rms amplitude of the oscillation during the observation. ${N}_{\mathrm{osc}}$, ${N}_{\mathrm{tot}}$, and ${f}_{\mathrm{rms}}$ are all directly observable. ${{\mathcal{R}}}^{-1}$ is the fractional variation in the shape of the waveform produced by the counting noise and is therefore a useful figure of merit for evaluating and comparing data sets.

${N}_{\mathrm{osc}}$ is approximately equal to the integral of the semi-amplitude of the first harmonic (fundamental) component of the burst oscillation waveform over the duration of the data segment and ${f}_{\mathrm{rms}}$ is approximately equal to the rms amplitude ${f}_{\mathrm{rms}1}$ of this harmonic because the amplitudes of the higher harmonics in burst oscillation waveforms are much smaller than the amplitude of the first harmonic (see Watts 2012; Lo et al. 2013; and Table 3 below). ${N}_{\mathrm{tot}}={N}_{\mathrm{spot}}+{N}_{\mathrm{back}}$, where ${N}_{\mathrm{spot}}$ is the total number of counts detected from the hot spot and ${N}_{\mathrm{back}}$ is the total number of background counts, which we define as all counts not produced by photons from the hot spot. ${N}_{\mathrm{spot}}$ is not usually equal to ${N}_{\mathrm{osc}}$, because for many geometries the hot spot contributes an unmodulated component to the waveform, as well as an oscillating component, i.e., ${N}_{\mathrm{spot}}={N}_{\mathrm{const}}+{N}_{\mathrm{osc}}$. The value of ${N}_{\mathrm{back}}$ during bursts cannot currently be measured or predicted (see Section 2.1.2). Consequently, although ${N}_{\mathrm{tot}}$ is directly observable, its components ${N}_{\mathrm{spot}}$ and ${N}_{\mathrm{back}}$ can only be determined by fitting models to the observed waveform.

Lo et al. (2013) found that for systems with properties that allow interesting constraints to be derived on M and ${R}_{\mathrm{eq}}$, i.e., systems in which the stellar rotation rate is ≳300 Hz and the spot center and the observer's sightline are both within 30° of the star's rotational equator, waveform data with ${\mathcal{R}}\gtrsim 400$ can provide estimates of M and ${R}_{\mathrm{eq}}$ with uncertainties ≲10%.

A hot spot small enough to be treated as a point source (i.e., with ${\rm{\Delta }}{\theta }_{\mathrm{spot}}\lesssim 10^\circ$) spans ≲1% of the stellar surface and therefore cannot produce enough counts to obtain interestingly tight constraints on M and ${R}_{\mathrm{eq}}$, even if data are collected from many tens of bursts from a single star using a detector with a collecting area ∼10 m² and then combined. This is demonstrated by the existing observations of burst oscillations, which show that although the oscillations observed during the first fraction of a second of an X-ray burst typically have high fractional amplitudes, they are not bright enough and do not last long enough to provide the number of counts needed to derive tight constraints on M and ${R}_{\mathrm{eq}}$, for a practical number of burst observations from a single star.

For example, Lo et al. (2013) showed that observation of a bright burst from the bright X-ray burst source 4U 1636–536 using a detector with an effective area ∼10 m² could achieve an ${\mathcal{R}}$-value ∼33 during the first 1/16 s of the burst, when the emitting area is small. ${\mathcal{R}}$ scales as the square root of the number of counts, so achieving an ${\mathcal{R}}$-value ∼400, sufficient to provide estimates of M and ${R}_{\mathrm{eq}}$ with uncertainties ≲10%, would require combining data from the early rises of ∼150 such bursts, which is impractical. Although the oscillation amplitude diminishes as the burst rises, using data from the entire 1/4 s burst rise would yield a slightly larger ${\mathcal{R}}$-value ∼50, because the longer duration of the observation yields a much larger number of counts. Achieving an ${\mathcal{R}}$-value ∼400 would therefore require combining data from the entire rises of ∼65 such bursts.

Although the oscillations observed during burst tails usually have smaller fractional amplitudes, they last much longer than the first fraction of a second of burst rises. They also have larger emitting areas. Using data collected during a 2 s oscillation train observed in the tail of a burst from any one of the eight neutron stars that produce such oscillation trains would yield an ${\mathcal{R}}$-value ∼80 (see Lo et al. 2013, Section 2.2.1). Hence, combining data from ∼25 observations of such oscillation trains from a single star would yield ${\mathcal{R}}\sim 400$, sufficient to obtain ∼10% constraints on M and ${R}_{\mathrm{eq}}$ for systems that have favorable properties (again see Lo et al. 2013).

These results show that obtaining constraints on M and ${R}_{\mathrm{eq}}$ with uncertainties ≲10% will probably require using data from the peaks and/or tails of bursts, when the hot spot is not infinitesimal in extent.

Given that it will probably be necessary to include data taken when the hot spot is not infinitesimal, it is important to include the angular radius ${\rm{\Delta }}{\theta }_{\mathrm{spot}}$ of the spot as a parameter in waveform fits, because for a given observed waveform, the most probable value of ${\rm{\Delta }}{\theta }_{\mathrm{spot}}$ and the distance d to the star are related. Assuming point emission during the full rise, the peak of the burst or its tail would improperly remove this degeneracy, artificially reducing the uncertainties in M and ${R}_{\mathrm{eq}}$ estimates and possibly biasing them. As Lo et al. (2013) showed using a full Bayesian analysis, knowledge of the distance to the star can improve somewhat the precision of M and ${R}_{\mathrm{eq}}$ estimates by removing this degeneracy.

2.1.2. Backgrounds

2.1.3. Other Aspects of the Waveform Models

In the OS approximation (Cadeau et al. 2007; Morsink et al. 2007), the spacetime exterior to the star is Schwarzschild but the stellar surface is oblate. Thus, in effect, the star is treated as an oblate shell of infinitesimal mass surrounding a concentrated sphere of gravitational mass M. In Schwarzschild coordinates, the line element anywhere outside the star is therefore

$$\begin{eqnarray}{{ds}}^{2} & = & -(1-2M/r){{dt}}^{2}+{{dr}}^{2}/(1-2M/r)\\ & & +{r}^{2}(d{\theta }^{2}+{\mathrm{sin}}^{2}\theta d{\phi }^{2}),\end{eqnarray} \tag{ 2 }$$

where t is the time as measured at infinity, r is the circumferential radius, and θ and ϕ are the standard polar and azimuthal angles. Here and henceforth we usually use geometrized units in which $G=c\equiv 1$. As in the S+D approximation discussed in Section 1, all special relativistic effects are treated exactly.

We assume that the rotating neutron star of interest is symmetric around its rotational axis. Morsink et al. (2007) found that for the families of stars they considered, if a star has an equatorial radius ${R}_{\mathrm{eq}}$ and an angular frequency Ω as seen at infinity, its radius as a function of colatitude θ is well described by (see their Equations (8)–(10) and their Table 1)

$$\begin{eqnarray}&&\displaystyle \frac{R(\theta )}{{R}_{\mathrm{eq}}}\approx 1+\underset{n=0}{\overset{2}{\sum }}{a}_{2n}(\zeta ,\epsilon ){P}_{2n}(\mathrm{cos}\theta ),\end{eqnarray} \tag{ 3 }$$

where ${P}_{2n}(\mathrm{cos}\theta )$ is the Legendre polynomial of order $2n$,

$$\begin{eqnarray}&&\zeta \equiv \displaystyle \frac{M}{{R}_{\mathrm{eq}}},\qquad \mathrm{and}\ \qquad \ \epsilon \equiv \displaystyle \frac{{{\rm{\Omega }}}^{2}{R}_{\mathrm{eq}}^{3}}{M}.\end{eqnarray} \tag{ 4 }$$

For neutron stars and hybrid quark stars,

$$\begin{eqnarray}{a}_{0} & = & -0.18\epsilon +0.23\zeta \epsilon -0.05{\epsilon }^{2},\\ {a}_{2} & = & -0.39\epsilon +0.29\zeta \epsilon +0.13{\epsilon }^{2},\mathrm{and}\\ {a}_{4} & = & 0.04\epsilon -0.15\zeta \epsilon +0.07{\epsilon }^{2}.\end{eqnarray} \tag{ 5 }$$

The area of an infinitesimal surface element at colatitude θ on the surface of an oblate star is (see Morsink et al. 2007, Equations (2) and (3))

$$\begin{eqnarray}&&{dS}(\theta )={R}^{2}(\theta )\mathrm{sin}\theta {\left[1+{f}^{2}(\theta )\right]}^{1/2}d\theta d\phi ,\end{eqnarray} \tag{ 6 }$$

where

$$\begin{eqnarray}&&f(\theta )\equiv \displaystyle \frac{{(1-2M/R)}^{-1/2}}{R}\displaystyle \frac{{dR}}{d\theta }.\end{eqnarray} \tag{ 7 }$$

The advantage of using the Schwarzschild spacetime rather than spacetimes that include the effect of the mass quadrupole of the star and frame dragging is that the spherical symmetry of the Schwarzschild spacetime guarantees that the path of any light ray in vacuum will lie in a plane. Thus, ray paths can be pre-computed and used in a lookup table, speeding up the computations enormously. We describe our table-lookup procedure later in this section.

When constructing OS synthetic waveforms and waveform models, we use the same ray-tracing codes described in Lo et al. (2013). Many of the equations and algorithms used in these codes were originally derived by Poutanen & Gierliński (2003).

We assume that the emission from the hot spot is azimuthally symmetric around the local surface normal, as seen in the surface comoving frame. The variation of the star's circumferential radius with colatitude creates an angle γ between this normal and the local outward radial vector given by (see Section 2.2 of Morsink et al. 2007)

$$\begin{eqnarray}&&\mathrm{cos}\gamma ={\left[1+{f}^{2}(\theta )\right]}^{-1/2}.\end{eqnarray} \tag{ 8 }$$

We assume further that the emission as seen in the surface comoving frame extends from the surface normal to the tangent to the surface, with a beaming function that is usually given by that expected for radiation from a uniform, semi-infinite, Thomson scattering atmosphere, which we approximate by (Lo et al. 2013)

$$\begin{eqnarray}&&{g}^{\prime }({\alpha }^{\prime })=0.42822+0.92236\mathrm{cos}{\alpha }^{\prime }-0.085751{\mathrm{cos}}^{2}{\alpha }^{\prime },\end{eqnarray} \tag{ 9 }$$

where ${\alpha }^{\prime }$ is the angle between the emitted ray and the surface normal, as seen in the surface comoving frame.

In addition to the minor modification to the element of surface area given by Equation (6), there are two more important differences between the OS approximation and the S+D approximation. The first is that the stellar radius varies with colatitude. As a result, if the equatorial circumferential radius is fixed, the maximum angular deflection of a photon leaving the star from a surface element not on the rotational equator is greater in the OS approximation than in the S+D approximation. The second difference is that in the S+D approximation, the maximum angle between the direction of an emitted photon and the local outward radial vector is $\pi /2$, assuming that the star is not more compact than the photon orbit $R=3M$ (for more compact stars, the maximum angle for an escaping photon is less than $\pi /2$; however, we do not consider such compact stars). In contrast, for an oblate star the tilt of the normal vector from the radial vector means that in the direction toward the closer rotational pole, the maximum angle is greater than $\pi /2$, whereas in the direction away from that pole, the maximum angle is less than $\pi /2$. More generally, if ψ is an azimuthal angle in the stellar surface that is 0 for the direction toward the closer pole and π for the direction away from it, then the angle χ between the local outward radial vector and the photon ray in direction ψ in the plane tangent to the stellar surface is given by $\mathrm{cos}\chi =-\mathrm{cos}\psi \mathrm{sin}\gamma$.

The two effects just described both tend to increase the minimum flux in the OS waveform relative to the S+D waveform for the same spot location and size, equatorial radius, gravitational mass, and observer colatitude, for observers in the same rotational hemisphere as the hot spot. To see why, consider a point on the surface somewhere between the pole and equator. The minimum in the waveform at a given energy will occur when the observer is on the opposite side of the star (modulo some effects related to relativistic aberration). In the OS approximation, the surface emission comes from a smaller radius and a ray can start in a direction with a smaller impact parameter to the center of the star than in the S+D approximation. Consequently, the gravitational light deflection is greater and causes the minimum flux to be greater than in the S+D approximation (see, e.g., Figures 3 and 4 of Cadeau et al. 2007 and Figures 3 and 4 of Morsink et al. 2007). In contrast, for observers in the opposite rotational hemisphere from the hot spot these two effects tend to decrease the minimum flux in the OS waveform relative to the S+D waveform (see, e.g., Figures 6 and 8 of Cadeau et al. 2007; but see their Figure 5 for a counterexample). The OS and S+D waveforms are very nearly the same if the hot spot is in the rotational equator (see, e.g., Figure 7 of Cadeau et al. 2007), because ${dR}/d\theta$ vanishes there.

Our algorithm computes the required ray paths in advance and stores them in a table that is then read by our code prior to calculating the required waveforms. In the Schwarzschild spacetime, the angular deflection of a ray to infinity depends on $R/M$, and not on R and M independently. We find that a table of ray paths for 440 evenly spaced values of $R/M$, from 3.6 to 8.0, and 1000 evenly spaced values of the angle between the ray and the outward radial direction, from 0 radians to 1.82 radians (recall that for oblate stars, the maximum angle from the outward radial direction can exceed $\pi /2\approx 1.57$) provides adequate accuracy: with this gridding, waveforms constructed using the ray table give a flux at any energy or phase that differs by no more than a few parts in 10⁵ from the corresponding directly traced waveform.

For a given initial radius and angle from the outward radial vector, we compute (1) the deflection angle to infinity (see Poutanen & Gierliński 2003 and Section A.1.1 of Lo et al. 2013), (2) the time delay of the ray relative to a radial ray (see Braje et al. 2000; Viironen & Poutanen 2004; and Section A.1.2 of Lo et al. 2013), and (3) the gravitational lensing factor for a cluster of rays near the fiducial ray.

Figure 1 compares an energy-integrated photon number flux waveform computed using our OS waveform code with the energy-integrated waveform for the same model parameters kindly computed by S. Morsink (2013, private communication), using the stellar shape given by the Stergioulas rotating neutron star code rns (Stergioulas 2003). In this figure only, the beaming pattern of the radiation from the hot spot (i.e., the angular distribution of the intensity from any point on the hot spot) was assumed to be isotropic in the surface comoving frame. The agreement between the two waveforms is excellent. Comparison of the two codes revealed that the difference between the two waveforms is due to the slight difference between the stellar shape computed using the rns code and the shape computed using the analytical method introduced by Morsink et al. (2007).

Zoom In Zoom Out Reset image size

Our statistical approach to estimating M and ${R}_{\mathrm{eq}}$ is based on Bayesian inference techniques, and follows closely the approach used in Lo et al. (2013). If we have a model with parameters ${\boldsymbol{y}}$, and if the prior probability distribution over those parameters is $q({\boldsymbol{y}}| I)$ (where I represents prior information), then Bayes' theorem states that the posterior probability distribution after analyzing data D is

$$\begin{eqnarray}&&P({\boldsymbol{y}}| D,I)\propto p(D| {\boldsymbol{y}},I)q({\boldsymbol{y}}| I).\end{eqnarray} \tag{ 10 }$$

If Poisson noise is the only source of fluctuations in the data, the likelihood of the "observed" data, given a particular set ${\boldsymbol{y}}$ of values for the model parameters, is

$$\begin{eqnarray}&&{\mathcal{L}}\equiv p(D| {\boldsymbol{y}},I)=\displaystyle \underset{i}{\prod }\displaystyle \frac{{m}_{i}{({\boldsymbol{y}})}^{{d}_{i}}}{{d}_{i}!}{e}^{-{m}_{i}({\boldsymbol{y}})},\end{eqnarray} \tag{ 11 }$$

where the product is over all the oscillation phase and energy bins, d_i is the measured number of counts in the ith bin, and ${m}_{i}({\boldsymbol{y}})$ is the number of counts in the ith bin predicted by the model for the trial set ${\boldsymbol{y}}$ of parameter values. If the number of counts in a given bin exceeds a few tens, then the Poisson likelihood may be replaced by a Gaussian.

Unlike frequentist statistics such as ${\chi }^{2}$, the value of ${\mathcal{L}}$ itself has no implication for whether the fit is good in an absolute sense. Instead, comparisons are made between different sets of parameter values and depend only on the ratio between different posterior probabilities, which means that they depend on the product of the ratio of the prior probabilities and the ratio of the likelihoods. The common factor ${\prod }_{i}(1/{d}_{i}!)$ therefore cancels out. In our analysis we adopt flat priors for all of our main parameters, within the range of values that we search for each parameter. In real situations, it is possible to have additional information about some of the parameters, and if so, that information should be included in the analysis via the prior. For example, Lo et al. (2013) found that knowledge of $M/{R}_{\mathrm{eq}}$ via an identified atomic line in the spectrum from the stellar surface, or of the observer's inclination angle, can tighten the constraints on M and ${R}_{\mathrm{eq}}$ considerably, whereas knowledge of the distance improves the constraints only modestly.

We adopt the common procedure of working with the log likelihood, which after removal of the $\sum \mathrm{log}(1/{d}_{i}!)$ term is

$$\begin{eqnarray}&&\mathrm{log}{\mathcal{L}}=\underset{i}{\sum }{d}_{i}\mathrm{log}{m}_{i}({\boldsymbol{y}})-\underset{i}{\sum }{m}_{i}({\boldsymbol{y}}).\end{eqnarray} \tag{ 12 }$$

It is the ratio of the likelihoods and thus the difference between the log likelihoods that matters. When we quote uncertainties at a certain level of confidence, we use Wilks' theorem (Wilks 1938), which states that ${\rm{\Delta }}{\chi }^{2}\approx -2{\rm{\Delta }}\mathrm{log}{\mathcal{L}}$ for a reasonably large total number of counts, and ${\chi }^{2}$ tables.

In our problem, we are primarily interested in M and ${R}_{\mathrm{eq}}$. If we designate the other parameters (called nuisance parameters in this context) by ${\boldsymbol{z}}$, then the correct way to find the final posterior probability distribution for only M and ${R}_{\mathrm{eq}}$ is to marginalize over the other parameters, i.e.,

$$\begin{eqnarray}&&P(M,{R}_{\mathrm{eq}})\propto \int d{\boldsymbol{z}}p(D| M,{R}_{\mathrm{eq}},{\boldsymbol{z}},I)q(M,{R}_{\mathrm{eq}},{\boldsymbol{z}}| I).\end{eqnarray} \tag{ 13 }$$

However, our waveform model has a large number of nuisance parameters. As explained in Section 2.1.1 and Lo et al. (2013), one must include in the waveform model an oscillation-phase-independent background component with an arbitrary magnitude and energy spectrum, which is specified by the number of phase-independent counts in each energy channel. The number of parameters in the background model is therefore equal to the number of energy channels. In addition, for a given set of candidate parameters one must consider an arbitrary shift in the phase of the model waveform relative to the waveform data, adding another parameter. It is not practical to marginalize fully over so many parameters. Consequently, we instead maximize the likelihood over the parameters in the background model and the start time of the model waveform, using a bisection method.

We have tested whether maximizing the likelihood of these nuisance parameters gives results comparable to marginalizing them, by fitting a Gaussian to the likelihood distributions found during the bisection procedure, analytically integrating over the Gaussian, and comparing the results with those obtained by simply using the maximum likelihood values of these parameters. We compared the log likelihood differences between these methods for combinations of $(M,{R}_{\mathrm{eq}},{\theta }_{\mathrm{obs}},{\theta }_{\mathrm{spot}},{\rm{\Delta }}{\theta }_{\mathrm{spot}})$ ranging from the combination that maximized $\mathrm{log}{\mathcal{L}}$ to combinations that gave $\mathrm{log}{\mathcal{L}}$ values 20 less than the maximum. For five parameters, Wilks' theorem (Wilks 1938) suggests that ${\rm{\Delta }}\mathrm{log}{\mathcal{L}}=2.94$ is approximately equivalent to 1σ. We found that for the background model, the standard deviation in the difference of $\mathrm{log}{\mathcal{L}}$ between the maximization and marginalization procedures was only 0.007. That is, even though the value of $\mathrm{log}{\mathcal{L}}$ for a given parameter combination in the marginalization procedure has an offset from the value of $\mathrm{log}{\mathcal{L}}$ for the same parameter combination in the maximization procedure, the offset is almost exactly constant from one parameter combination to the next. Thus, the differences between log likelihoods are preserved and hence maximization of the likelihood over the background parameters is functionally equivalent to marginalization. When we performed a similar comparison of maximization to marginalization for the shift in the start time of the waveform, we found that the standard deviation of the difference in $\mathrm{log}{\mathcal{L}}$ between the two methods is 0.3. This is therefore a $0.1\sigma$ shift, which is too small to affect any of our results.

A challenge faced in Lo et al. (2013) was that if the waveform data are informative, i.e., if they tightly constrain M and ${R}_{\mathrm{eq}}$, a grid search over the values of the parameters in the waveform model requires a grid so fine that a truly blind search would require a prohibitive number of waveform computations, even if large computational resources are available.

In this section, we describe a new waveform data processing procedure that is based on the approach used by Lo et al. (2013) but performs blind searches far more efficiently. In this new procedure we first determine the volume of the waveform parameter space of interest by performing a Markov chain Monte Carlo (MCMC) search using ray-tracing. We then construct and interpolate in a table of waveforms to further localize the volume of interest and compute bounding ellipsoids that encompass the points in this reduced volume that have interestingly high likelihoods, finally sampling the volumes within these ellipsoids by performing a Monte Carlo integration using direct ray-tracing. The details are as follows:

• 1.  
  We explore ($M,{R}_{\mathrm{eq}},{\theta }_{\mathrm{obs}},{\theta }_{\mathrm{spot}},{\rm{\Delta }}{\theta }_{\mathrm{spot}}$) space using an MCMC search with ray-tracing (see Section 3.2 of Lo et al. 2013 for our MCMC approach and, e.g., von Toussaint 2011 for a discussion of the Metropolis algorithm for MCMC sampling). In general, one would also need to explore a range of values of the surface comoving temperature T and the distance d to the star, but in the approach we use here we either assume that the redshifted temperature $T{(1-2M/{R}_{\mathrm{eq}})}^{1/2}$ is known and maximize the likelihood over d, or assume that d is known and maximize the likelihood over T (see Section 3.3.3 of Lo et al. 2013 for an explanation and justification of this approach).
• 2.  
  We determine the maximum log likelihood ($\mathrm{log}{{\mathcal{L}}}_{\mathrm{max}}$) of the data given the waveform model, in the volume of the waveform parameter space explored in the previous step. We then select the $(M,{R}_{\mathrm{eq}},{\theta }_{\mathrm{obs}},{\theta }_{\mathrm{spot}},{\rm{\Delta }}{\theta }_{\mathrm{spot}})$ combinations in this volume that have log likelihoods within 30 of the maximum and determine the range of each of these five parameters that spans the selected combinations. Finally, we generate a table of 10⁵ template waveforms, one for each point on the five-dimensional grid of 10 evenly spaced values that span the previously determined range of each parameter.
• 3.  
  We perform a second MCMC calculation using waveforms computed by interpolating in the table of template waveforms generated in the previous step. In this calculation, we pick an $(M,{R}_{\mathrm{eq}})$ pair and marginalize by performing a Monte Carlo integration with 1000 points spanning $({\theta }_{\mathrm{obs}},{\theta }_{\mathrm{spot}},{\rm{\Delta }}{\theta }_{\mathrm{spot}})$. To do this, we first compute for the selected $(M,{R}_{\mathrm{eq}})$ the log likelihood on an evenly spaced $10\times 10\times 10$ grid in the three angular variables, over the full range of each variable determined in step 1. Using the points on this grid with log likelihoods within 40 of the maximum found in the previous step, we determine, for each value of ${\theta }_{\mathrm{obs}}$ on our grid, the minimum ellipse that contains all $({\theta }_{\mathrm{spot}},{\rm{\Delta }}{\theta }_{\mathrm{spot}})$ pairs with $\mathrm{log}{\mathcal{L}}\gt \mathrm{log}{{\mathcal{L}}}_{\mathrm{max}}-40$ (see the Appendix for a description of the algorithm we used to find a minimum bounding ellipse). If, for our chosen $(M,{R}_{\mathrm{eq}})$ pair, any of these ellipses are nonvanishing, we perform a 1000-point Monte Carlo integration over the $({\theta }_{\mathrm{obs}},{\theta }_{\mathrm{spot}},{\rm{\Delta }}{\theta }_{\mathrm{spot}})$ volume identified by the nonvanishing bounding ellipses. The use of bounding ellipses typically reduces the volume over which the integral must be performed by a factor ∼30, which allows us to sample the relevant volume much more densely, reducing the fractional error in the integral by a factor $\sim {30}^{1/2}\sim 5$ for a given number of evaluations.
• 4.  
  We take all the $(M,{R}_{\mathrm{eq}})$ pairs from the previous step that yielded nonvanishing marginalized posterior probabilities and use ray-tracing rather than waveform interpolation to compute the log likelihood in a full, uniformly spaced $10\times 10\times 10$ grid over $({\theta }_{\mathrm{obs}},{\theta }_{\mathrm{spot}},{\rm{\Delta }}{\theta }_{\mathrm{spot}})$. Then, as in the previous step, we construct bounding ellipses in $({\theta }_{\mathrm{spot}},{\rm{\Delta }}{\theta }_{\mathrm{spot}})$ at each grid value of ${\theta }_{\mathrm{obs}}$ for each $(M,{R}_{\mathrm{eq}})$ pair and perform a Monte Carlo integration with 1000 points per $(M,{R}_{\mathrm{eq}})$ combination, using direct ray-tracing rather than waveform interpolation. This yields the marginalized posterior probability at points in the $({R}_{\mathrm{eq}},M)$ plane that, because of the MCMC sampling procedure, are concentrated around the maximum posterior probability.
• 5.  
  We normalize the marginalized posterior probability to its maximum and then determine and output the 1σ, 2σ, and 3σ contours in the $(M,{R}_{\mathrm{eq}})$ plane.

2.5.1. Lo et al. (2013)

2.5.2. Cadeau et al. (2007)

All the synthetic observed waveforms we analyze here were generated using the OS approximation, assuming a stellar gravitational mass M of $1.6\;{M}_{\odot }$, a circular hot spot with an angular radius ${\rm{\Delta }}{\theta }_{\mathrm{spot}}$ of 25°, and a distance d to the neutron star of 10 kpc. Table 1 lists the rotational frequency ${\nu }_{\mathrm{rot}}$ of the hot spot as seen at infinity, the stellar equatorial radius ${R}_{\mathrm{eq}}$, the hot-spot inclination ${\theta }_{\mathrm{spot}}$, and the observer inclination ${\theta }_{\mathrm{obs}}$ used to generate each synthetic observed waveform. We assumed ${\nu }_{\mathrm{rot}}$ equal to either 300 or 600 Hz, ${R}_{\mathrm{eq}}$ equal to either 11.8 or 15 km, and ${\theta }_{\mathrm{spot}}$ and ${\theta }_{\mathrm{obs}}$ both equal to either 60° or 90°. We independently generate the synthetic waveform for each case that we analyze in Figures 2 and 3. In case 4 we used the same spot data as in case 3(a), but we independently generated new background data. In cases 5(a) and (b), we used the same waveform realization but two different phase binnings.

All the synthetic observed waveforms except one were generated assuming a uniform hot spot that emits radiation with a blackbody spectrum having a temperature of 2 keV as measured in the surface comoving frame. In case 3(d), we consider a hot spot with a temperature that varies with latitude (see below). Our assumption that the emission from the hot spot has a blackbody spectrum and normalization is formally inconsistent with our assumption that the beaming pattern of the radiation from the hot spot is that of an electron-scattering atmosphere (see Lo et al. 2013); when real data are analyzed, one should use emission spectra and normalizations computed using appropriate model atmospheres.

We assume that the background does not vary at frequencies commensurate with the spot rotational frequency. We include a placeholder background in the synthetic observed waveform by assuming uniform emission from the entire surface of a star with the same mass and radius as the star used in modeling the hot-spot waveform, with a spectrum having the shape of a Planck spectrum as seen in a frame comoving with the surface of the star, which for this purpose is assumed to be rotating with the same frequency as the hot spot. In all cases except one, we assume that the background emission has a temperature of 1.5 keV. In case 4, we assume that the temperature of the background emission is the same as the temperature of the emission from the hot spot, i.e., 2 keV, in order to explore whether this weakens the derived constraints on M and ${R}_{\mathrm{eq}}$. We normalize the background component to produce the desired number of expected background counts.

In all cases except one, we Poisson-sampled the synthetic waveforms to generate count data in 16 equally spaced phase bins and 30 equally spaced energy channels. Each synthetic waveform therefore consists of the number of counts in each of $16\times 30=480$ phase-energy bins. In case 5(a), we first Poisson-sampled the synthetic waveform in 32 bins. Thus the synthetic waveform in case 5(a) consists of the number of counts in each of 32 × 30 = 960 phase-energy bins. We then grouped the data in this sampled waveform into 16 bins (case 5(b)), to test whether using 16 phase bins provides adequate phase resolution. For all the synthetic waveforms except the one analyzed in case 3(d), the centroids of the energy channels are spaced 0.3 keV apart and run from 3.65 to 12.35 keV; for the waveform analyzed in case 3(d), the centroids are spaced 0.3 keV apart and run from 1.85 to 10.55 keV.

In all the synthetic observed waveforms, the expected number of counts from the spot is 10⁶, whereas the expected number of counts from the background is $9\times {10}^{6}$ (the actual numbers vary because of the Poisson sampling). As noted in Section 2.1, these numbers produce a realistic modulation amplitude and a total number of counts comparable to the number that could be obtained by the accepted NICER mission and the proposed LOFT and AXTAR missions, by combining data from many bursts from a given star.

The standard S+D and OS waveform models that we fit to the synthetic observed waveform data both assume that the temperature of the hot spot is uniform as seen in the surface comoving frame. In fitting these models to the synthetic waveform data, we assume that the rotational frequency of the star is the same as the rotational frequency ${\nu }_{\mathrm{rot}}$ of the hot spot and that ${\nu }_{\mathrm{rot}}$ is known. Both our standard waveform models have seven primary adjustable parameters (M, ${R}_{\mathrm{eq}}$, ${\theta }_{\mathrm{spot}}$, ${\rm{\Delta }}{\theta }_{\mathrm{spot}}$, ${{kT}}_{\mathrm{spot}}$, ${\theta }_{\mathrm{obs}}$, and d) and 31 ancillary adjustable parameters (the phase-independent background in each of the 30 energy channels and the overall time shift). Each waveform model therefore has 38 adjustable parameters, and hence there are $480-38$ = 442 (or, in case 5(a), $960-38$ = 922) dof. We assume uniform priors over the allowed ranges of all the parameters in our model waveforms. We perform a blind search over M, from $1.2$ to $2.2\;{M}_{\odot }$; over ${R}_{\mathrm{eq}}/M$, from 4 to 8; and over ${\theta }_{\mathrm{spot}}$, ${\rm{\Delta }}{\theta }_{\mathrm{spot}}$, and ${\theta }_{\mathrm{obs}}$, from 0.1 to $\pi /2$ radians. In principle, ${\theta }_{\mathrm{obs}}$ could range from 0 to π radians. However, we find that allowing this larger range does not change the confidence regions for M and ${R}_{\mathrm{eq}}$. Restricting ${\theta }_{\mathrm{obs}}$ to $\leqslant \pi /2$ radians allows us to sample the angular range of interest with a higher density of points.

In all cases except 3(d), we assume that ${{kT}}_{{\infty }}$, the radiation temperature measured at infinity, is ${{kT}}_{\mathrm{spot}}{(1-2M{/R}_{\mathrm{eq}})}^{1/2}$. This assumption is justified because the energy of the spectral peak will in practice be very precisely measured. For all cases except 3(d), we maximize the likelihood over the distance d for a given combination of the other parameters as described in Section 2.3, rather than marginalizing over d. For a more detailed explanation and justification of this approach, see Section 3.3.3 of Lo et al. (2013). In case 3(d), the synthetic waveform data were generated with a surface comoving temperature that varies with latitude over the hot spot. Hence, in analyzing these waveform data we assume that we know d and maximize the likelihood over ${{kT}}_{\mathrm{spot}}$, rather than marginalizing over ${{kT}}_{\mathrm{spot}}$. When real data are analyzed, it will be necessary to marginalize the likelihood over both ${{kT}}_{\mathrm{spot}}$ and d.

In order to save computational time, we determine the magnitude and spectrum of the background and the time shift of the model waveform relative to the synthetic waveform by maximizing the likelihood, as described in Section 2.3, rather than by marginalizing these parameters.

Figure 2 shows the constraints on M and ${R}_{\mathrm{eq}}$ obtained by fitting our standard S+D waveform model to synthetic observed waveform data generated using the OS approximation. Figures 2(a)–(c) show the results obtained for two values of ${\nu }_{\mathrm{rot}}$ and ${R}_{\mathrm{eq}}$ when the colatitude of the hot spot and the inclination of the observer are both 60°, whereas Figure 2(d) shows a typical example of the results obtained when the colatitude of the hot spot and the inclination of the observer are both 90°. In all these cases, fitting the data using our standard S+D model does not significantly bias the estimates of M but does significantly bias the estimates of ${R}_{\mathrm{eq}}$ in cases 2(b)–(d). The ${\chi }^{2}$ values for these fits (see the last column of Table 1) indicate that they are all statistically good. We now discuss each of these cases in turn.

Zoom In Zoom Out Reset image size

Figure 2(a) shows the constraints on M and ${R}_{\mathrm{eq}}$ obtained by fitting our standard S+D waveform model to the synthetic waveform produced by a hot spot at ${\theta }_{\mathrm{spot}}=60^\circ$, rotating at ${\nu }_{\mathrm{rot}}=300$ Hz, on the surface of a star with a moderate radius (${R}_{\mathrm{eq}}=11.8$ km). As noted above, in fitting our waveform models to synthetic waveform data we assume that the rotational frequency of the star is the same as the rotational frequency of the hot spot; consequently, in this case we assume that the stellar rotational frequency is 300 Hz. The oblateness of the stellar surface, and hence the deviation of the synthetic observed waveform from the waveform in the S+D model, scales as the square of the stellar rotational frequency. We therefore expect, and find, that for the moderate rotational frequency of this case, the bias in the estimated value of ${R}_{\mathrm{eq}}$ is not significant: the edge of the 1σ confidence region touches the values of M and ${R}_{\mathrm{eq}}$ that were used in generating the synthetic observed waveform, which are indicated by the black square. For the moderate spot colatitude, spot rotational frequency, and stellar radius of this case, the amplitudes of the overtones of ${\nu }_{\mathrm{rot}}$ in the synthetic observed waveform are very low, and the constraints on M and ${R}_{\mathrm{eq}}$ are therefore expected to be weak (see, e.g., Lo et al. 2013). This expectation is confirmed by Figure 2(a): the 1σ contour is large, and the 2σ and 3σ contours are even larger, extending toward high M and ${R}_{\mathrm{eq}}$ and intersecting the lower boundary of our search domain at $M={R}_{\mathrm{eq}}/8$.

Figure 2(b) displays the constraints on M and ${R}_{\mathrm{eq}}$ obtained by fitting our standard S+D waveform model to the synthetic waveform produced by a hot spot that is again at ${\theta }_{\mathrm{spot}}=60^\circ$ on a star with ${R}_{\mathrm{eq}}=11.8$ km, but now rotating at ${\nu }_{\mathrm{rot}}=600$ Hz. Even though the rotational frequency is twice that in the case featured in Figure 2(a), M and ${R}_{\mathrm{eq}}$ are still poorly constrained, because the hot spot is not near the rotational equator. The oblateness of the ${R}_{\mathrm{eq}}=11.8$ km star in Figure 2(b) should be $\approx 4$ times larger than the oblateness of the 11.8 km star featured in Figure 2(a) (if the star featured in Figure 2(a) were spun up to 600 Hz, the increase in its rotational distention would cause its equatorial radius to be slightly larger than the 11.8 km radius assumed in Figure 2(b)). Figure 2(b) shows that the $\approx 4$ times larger oblateness of this star is sufficient to introduce a significant bias in the values of M and ${R}_{\mathrm{eq}}$ estimated by fitting our standard S+D waveform model, which assumes the star is spherical. The values of M and ${R}_{\mathrm{eq}}$ that were used to generate the synthetic observed waveform are well outside the 3σ contour, partly because the contours in this case are modestly smaller than in the case featured in Figure 2(a), due to the star's higher rotational frequency.

Figure 2(c) shows the constraints on M and ${R}_{\mathrm{eq}}$ obtained by fitting our standard S+D waveform model to the synthetic waveform produced by a hot spot that is again at ${\theta }_{\mathrm{spot}}=60^\circ$ and rotating at ${\nu }_{\mathrm{rot}}=600$ Hz, but now on the surface of a star with ${R}_{\mathrm{eq}}=15$ km. The contours are smaller than those in Figure 2(b) because ${R}_{\mathrm{eq}}$, and hence the surface rotational velocity, is larger, but still extend to large M and ${R}_{\mathrm{eq}}$, again intersecting the lower boundary of our search domain at $M={R}_{\mathrm{eq}}/8$. The estimated values of M and ${R}_{\mathrm{eq}}$ are again significantly biased.

Figure 2(d) illustrates the constraints on M and ${R}_{\mathrm{eq}}$ typically obtained by fitting our standard S+D waveform model to OS synthetic waveforms produced by a hot spot near the rotational equator, when it is viewed by an observer at a high inclination relative to the rotational axis. The confidence regions shown in this figure were obtained by fitting the synthetic OS waveform data produced by a hot spot with ${\nu }_{\mathrm{rot}}=600$ Hz, centered on the rotational equator of a star with ${R}_{\mathrm{eq}}=11.8$ km, and viewed by an observer who is in the plane defined by the star's rotational equator. The constraints are much tighter than in the previous figures because of the high surface rotational velocity and the large projection of the velocity along the line-of-sight to the observer. The fractional biases in the estimates of M and ${R}_{\mathrm{eq}}$ are smaller than in the cases discussed previously but are still statistically significant, because the constraints on M and ${R}_{\mathrm{eq}}$ are much tighter. If the spot were infinitesimal in extent, the biases would be zero, because the stellar oblateness has no effect on the waveform produced by a point source if ${dR}/d\theta =0$ at the source, and ${dR}/d\theta$ is zero on the rotational equator where this spot is located. However, the synthetic observed waveform used in this case was generated using a hot spot with an angular radius of 25°, leading to a small but significant bias. Because this fit, like those discussed previously, is formally statistically acceptable, the quality of the fit does not by itself provide an indication that the M and ${R}_{\mathrm{eq}}$ estimates are biased.

These results show that fitting our standard S+D waveform model to OS synthetic waveform data tends to produce estimates of the star's equatorial radius that are significantly larger than its true radius, but produces estimates of the star's mass that do not differ from the true mass by a statistically significant amount. The oscillations produced by a hot spot near the star's rotational equator have the relatively strong overtones needed to provide tight constraints on M and ${R}_{\mathrm{eq}}$. As Figure 2(d) shows, the estimates of M and ${R}_{\mathrm{eq}}$ derived from such waveforms are also less susceptible to bias caused by the oblateness of the star.

We conclude that if a neutron star has a large radius or a rotational frequency ≳300 Hz, one should fit OS waveform models, rather than S+D waveform models, to the waveform data.

Figure 3 shows examples of the constraints on M and ${R}_{\mathrm{eq}}$ obtained by fitting our standard OS waveform model to synthetic waveform data generated using the OS approximation. Figures 3(a) and (b) show the results obtained by fitting this model to the synthetic waveforms produced by a hot spot on the star's rotational equator, observed at an inclination of 90°, for rotation rates of 600 and 300 Hz, whereas Figures 3(c) and (d) show the results obtained by fitting this model to the waveforms produced by a hot spot at a colatitude of 60°, observed at an inclination of 60°, for a rotation rate of 600 Hz (the waveform in case 3(d) was generated assuming a temperature variation across the hot spot, as discussed below). The ${\chi }^{2}$ values for these fits (see the last column of Table 1) indicate that they are all statistically good. The approximate 1σ uncertainties in M and ${R}_{\mathrm{eq}}$ for each of these fits are listed in Table 2. Comparison of this table with Table 2 of Lo et al. (2013) shows that the constraints on M and ${R}_{\mathrm{eq}}$ obtained by fitting our standard OS waveform model to OS waveform data are similar to the constraints obtained there by fitting our standard S+D waveform model to S+D waveform data. We now discuss these results in more detail.

Zoom In Zoom Out Reset image size

Table 2.  Constraints on M and ${R}_{\mathrm{eq}}$ Obtained Using Our Standard OS Waveform Model^a

Figure	${\nu }_{\mathrm{rot}}$ b	${\theta }_{\mathrm{spot}}$ c	${\theta }_{\mathrm{obs}}$ d	${R}_{\mathrm{eq}}$ e	${\rm{\Delta }}{R}_{\mathrm{eq}}$ f	${\rm{\Delta }}M$ g	$\delta {R}_{\mathrm{eq},1}$ h,j	$\delta {M}_{1}$ i,j
(case)	(Hz)	(deg)	(deg)	(km)	(km)	$({M}_{\odot })$	(%)	(%)
3(a)	600	90	90	11.8	11.57–12.27	1.57–1.66	2.9	2.8
3(b)	300	90	90	11.8	11.14–12.97	1.44–1.61	7.6	5.6
3(c)	600	60	60	15	13.83–15.83	1.58–1.80	6.7	6.5
3(d)	600	60	60	11.8	9.99–11.49	1.41–1.62	7.0	6.9
4	600	90	90	11.8	11.35–12.19	1.49–1.59	3.6	3.2
5(a)	600	90	90	11.8	11.63–12.03	1.56–1.64	1.7	2.5
5(b)	600	90	90	11.8	11.68–12.24	1.56–1.64	2.3	2.5

Notes.

^aThe 1σ M and ${R}_{\mathrm{eq}}$ uncertainties listed here somewhat understate the actual constraints on $M({R}_{\mathrm{eq}})$, because $M/{R}_{\mathrm{eq}}$ is usually better constrained than M or ${R}_{\mathrm{eq}}$ considered separately (see, e.g., Figure 3(a)). All the synthetic waveform data were generated assuming $M=1.6\;{M}_{\odot }$. The data analyzed in case 3(d) were generated assuming a temperature variation across the hot spot (see text). ^bRotational frequency of the hot spot as seen at infinity assumed in generating the synthetic observed waveform data. ^cInclination (colatitude) of the hot-spot center assumed in generating the synthetic observed waveform data. ^dInclination of the observer assumed in generating the synthetic observed waveform data. ^eEquatorial radius assumed in generating the synthetic observed waveform data. ^fRange of the 1σ contour projected onto the ${R}_{\mathrm{eq}}$ axis. ^gRange of the 1σ contour projected onto the M axis. ^hApproximate 1σ fractional uncertainty in ${R}_{\mathrm{eq}}$ computed by dividing one-half the 1σ range of ${R}_{\mathrm{eq}}$ by its central value. ⁱApproximate 1σ fractional uncertainty in M computed by dividing one-half the 1σ range of M by its central value. ^jThe definitions of the 1σ uncertainties in M and ${R}_{\mathrm{eq}}$ used here differ from those used in Lo et al. (2013), where they were estimated by projecting the full extent of the 1σ contour onto the M and R_eq axes, because in that work some of the 1σ confidence regions were highly asymmetric, with best-fit values of M and R_eq far from the center.

Download table as:  ASCIITypeset image

Figure 3(a) illustrates the constraints on M and ${R}_{\mathrm{eq}}$ typically obtained when our standard OS waveform model is fit to OS synthetic waveforms produced by a hot spot near the rotational equator, viewed by an observer at a high inclination relative to the rotational axis. The confidence regions shown in this figure were obtained by analyzing the OS waveform produced by a hot spot rotating at ${\nu }_{\mathrm{rot}}=600$ Hz, centered on the rotational equator of a star with ${R}_{\mathrm{eq}}=11.8$ km, when seen by an observer who is in the plane defined by the star's rotational equator. The constraints in Figure 3(a) are tighter than in Figure 2(d) and the estimates of M and ${R}_{\mathrm{eq}}$ are not significantly biased. The constraints in Figure 3(a) (1σ uncertainties of 2.9% in ${R}_{\mathrm{eq}}$ and 2.8% in M; see Table 2) are much tighter than in Figure 3(b), which shows results for a much lower rotational frequency (300 Hz), and much tighter than in Figure 3(c), which shows results for a hot spot that is not near the rotational equator (${\theta }_{\mathrm{obs}}=60^\circ$).

It is useful to compare Figure 3(a) with Figure 2(e) of Lo et al. (2013). In both cases, the synthetic waveform corresponds to a hot spot rotating at ${\nu }_{\mathrm{rot}}=600$ Hz, centered on the rotational equator of a star with ${R}_{\mathrm{eq}}=11.8$ km, and seen by an observer in the plane defined by the star's rotational equator. The 1σ confidence region in Figure 2(a) is comparable in size to the 1σ confidence region in Figure 2(e) of Lo et al. (2013), but the 3σ region in Figure 3(a) is much smaller than the 3σ region in Figure 2(e) of Lo et al. (2013).⁵ There are at least four possible explanations for this difference: (1) a statistical fluctuation in the synthetic waveform data of Lo et al. (2013) made that waveform realization less constraining than it would typically be, (2) a statistical fluctuation in the synthetic waveform data used here made this waveform realization more constraining than it would typically be, (3) our new analysis pipeline does a better job of representing the true constraints, or (4) something about OS waveforms actually yields better 3σ (but not 1σ) constraints in this situation. It is not clear without further exploration which, if any, of these explanations is responsible for this difference.

Figure 3(b) shows the constraints on M and ${R}_{\mathrm{eq}}$ obtained by analyzing a synthetic waveform produced by the same hot spot and observer geometry as in Figure 3(a) (the hot spot is on the rotational equator of a star with ${R}_{\mathrm{eq}}=11.8$ km and is seen by an observer who is in the plane defined by the star's rotational equator), but for a rotation rate of 300 Hz, half the rotation rate assumed in Figure 3(a). The slower rotation rate decreases the harmonic content of the waveform, which increases the sizes of the confidence regions, but the 1σ uncertainties (5.6% in M and 7.6% in ${R}_{\mathrm{eq}}$; see Table 2) are still interestingly small. This figure shows that fitting the waveforms of stars that rotate at a moderate rate can provide interesting constraints on M and ${R}_{\mathrm{eq}}$, provided that the hot spot is near the rotational equator and the observer is at a high inclination.

Figure 3(c) displays the constraints on M and ${R}_{\mathrm{eq}}$ obtained by analyzing a synthetic waveform produced by a hot spot on a star with a larger radius (${R}_{\mathrm{eq}}=15$ km) again rotating at 600 Hz, but with the spot at a colatitude of 60° and the observer at an inclination of 60°. The constraints are much less precise than for the case shown in Figure 3(a) and less precise overall than for the case shown in Figure 3(b), even though the star has a larger radius, because the hot spot is not near the rotational equator and the observer's inclination is not close to 90°. Even so, the 1σ uncertainties (6.5% in M and 6.7% in ${R}_{\mathrm{eq}}$; see Table 2) are still interestingly small.

Figure 3(d) shows a case in which the properties of the hot spot assumed in our standard OS waveform model are different from the properties used to generate the synthetic waveform data to which the model was fit. In this case, the fitted model assumes that the surface temperature is uniform across the hot spot, as seen in the frame comoving with the surface, whereas the synthetic waveform data were generated assuming that the temperature varies by 25% with colatitude, from 2 keV at the center of the spot to 1.5 keV at the top and bottom edges of the spot. This case tests whether such a variation, which is physically plausible, produces a significant bias in the estimated values of M and ${R}_{\mathrm{eq}}$. As Figure 3(d) shows, there is a slight bias, but it is not statistically significant: the best-fit values of M and ${R}_{\mathrm{eq}}$ differ by only ∼1.5σ from the values used to generate the synthetic waveform data.

Figure 4 explores the effects on the mass and radius constraints if there is less contrast between the energy spectrum of the emission from the hot spot and the spectrum of the phase-independent background. For this analysis, we took the counts from the hot spot in each phase-energy bin generated for the analysis shown in Figure 3(a), but generated a new background having $\approx 9\times {10}^{6}$ counts and a spectrum identical to that of the hot spot (i.e., a Planck spectrum with a surface comoving temperature of 2.0 keV instead of the value of 1.5 keV used in generating our other synthetic waveforms). Hence, in this case the emission from the spot is simply extra emission on top of the background, rather than extra emission with a different spectrum. The constraints obtained are only marginally worse than in case 3(a). Given the randomness inherent in the generation of the synthetic data and in the MCMC analysis of the data, this difference may not be significant. This result shows that even if the phase-independent background has the same spectrum as the emission from the hot spot, one can obtain good constraints on M and ${R}_{\mathrm{eq}}$.

Zoom In Zoom Out Reset image size

Figure 5 explores whether 16 phase bins are adequate for the cases considered here. We first generated a synthetic observed waveform using the OS approximation and 32 phase bins, for a hot spot centered on the rotational equator of a star rotating at 600 Hz and an observer in the plane of the star's rotational equator. We then regrouped these same data into 16 phase bins. Figure 5 shows that the constraints on M obtained by analyzing the data binned in these two different ways are essentially identical, whereas the constraint on ${R}_{\mathrm{eq}}$ appears slightly tighter when 32 phase bins are used. The randomness in the generation of each synthetic waveform and in the MCMC analysis of a given synthetic waveform produces variations in the derived constraints; for example, Table 2 indicates that the 1σ constraints derived in case 5(b) are slightly tighter than in case 3(a), even though both used waveform data with 16 phase bins. Thus, the apparent slight improvement in the constraint on ${R}_{\mathrm{eq}}$ when 32 phase bins are used may not be significant.

Zoom In Zoom Out Reset image size

Understanding the uncertainties $\delta M$ and $\delta {R}_{\mathrm{eq}}$ in M and ${R}_{\mathrm{eq}}$ estimates requires understanding how the properties of the observed waveforms constrain M and ${R}_{\mathrm{eq}}$. This is explained in detail by Lo et al. (2013; see also Psaltis et al. 2014). The asymmetry and harmonic content of the waveform constrain the component of the velocity of the emitting region in the observer's direction, primarily via special relativistic Doppler boosts and aberration. Because the rotational frequency of the emitting region is accurately known from the oscillation frequency, knowing the line-of-sight velocity of the emitting region constrains the stellar radius. The amplitude of the waveform constrains the colatitude of the hot spot, the observer's inclination, and the compactness ($M/R$) of the star, the last primarily via general relativistic light-bending effects.

Although the harmonic content of the waveform encodes information about the rotational velocity of the emitting element—and hence the cylindrical radius of the element (see Lo et al. 2013, Section 2.2.2)—other unrelated aspects of the system also affect the harmonic structure of the waveform. In particular, the semi-amplitude C₂ of the second harmonic in the waveform is not uniquely related to the rotational velocity of the emitting element. Consequently, there is in general no simple way to extract information about ${R}_{\mathrm{eq}}$ from the harmonic content of the waveform; model fitting is required.

Effects other than the rotational velocity that contribute to C₂ include (1) anisotropic beaming of the radiation from each emitting element and (2) occultation of part or all of the hot spot by the star, as it rotates (see Poutanen & Beloborodov 2006, which provides a useful guide to waveform properties, even though the results reported there assume the hot spot is infinitesimal in extent and were derived using an approximate analytic expression for the light deflection). To lowest order, the line-of-sight linear velocity of the emitting gas ${v}_{\mathrm{los}}$ contributes a second harmonic with semi-amplitude ${C}_{2}\propto ({v}_{\mathrm{los}}/c){C}_{1}$, where c is the speed of light, C₁ is the semi-amplitude of the first harmonic (the fundamental), and the coefficient of proportionality depends on the spectrum of the emission (see Poutanen & Beloborodov 2006, Equation (66)). The linear velocity produced by the rotation of the emitting gas therefore contributes a second harmonic with semi-amplitude

$$\begin{eqnarray}&&{C}_{2}\propto (2\pi {\nu }_{\mathrm{rot}}{R}_{\mathrm{eq}}/c)\mathrm{sin}\theta \mathrm{sin}{\theta }_{\mathrm{obs}},\end{eqnarray} \tag{ 14 }$$

where θ is the colatitude of the surface element from which the radiation is emitted. Because this contribution to C₂ is ${\mathcal{O}}({v}_{\mathrm{los}}/c){C}_{1}$ and ${v}_{\mathrm{los}}/c\ll 1$, it is generally $\ll {C}_{1}$. In contrast, the beaming of emission from the atmosphere of the hot spot (see Section 2.1.3) contributes a second harmonic with semi-amplitude ${C}_{2}\approx {C}_{1}$. In particular, emission with the anisotropic beaming pattern $I({\alpha }^{\prime })={I}_{0}(1+h\mathrm{cos}{\alpha }^{\prime })$ produces a second harmonic with semi-amplitude ${C}_{2}\approx h\mathrm{sin}\theta \mathrm{sin}{\theta }_{\mathrm{obs}}{C}_{1}$ (see Poutanen & Beloborodov 2006, Equation (50)); for a burst atmosphere, $h\approx 0.92$ (see Equation (9)). Because the component of the second harmonic contributed by the anisotropy of the radiation from the hot spot is independent of the star's rotational frequency, it is ≈C₁ even if the star is rotating slowly. Occultation of part or all of the hot spot by the star can also contribute a second harmonic component with a semi-amplitude ${C}_{2}\approx {C}_{1}$. Thus, in burst oscillations the second harmonics generated by other effects can be comparable to or larger than the second harmonic produced by the rotational velocity of the star.

In order to obtain useful constraints on M and ${R}_{\mathrm{eq}}$, the velocity of the hot spot surface must make a significant contribution to the harmonic content of the waveform. Equation (14) shows that this contribution is greater if the spot and the observer's sightline are closer to the star's rotational equator. Then the relativistic Doppler shift and aberration are greater, the waveform depends more sensitively on the radius of the star, and the constraints on M and ${R}_{\mathrm{eq}}$ are tighter.

The uncertainties in estimates of M and ${R}_{\mathrm{eq}}$ also depend on the fractional amplitude of the oscillation and the total number of counts, including the number of background counts, which we define as all counts not produced by photons from the hot spot. As explained in Section 2.1.1, if all the other properties of the system that affect the waveform are kept fixed, the uncertainties in estimates of M and ${R}_{\mathrm{eq}}$ obtained by fitting waveform models to waveform data scale approximately as ${{\mathcal{R}}}^{-1}$, where ${\mathcal{R}}=1.4{f}_{\mathrm{rms}}\sqrt{{N}_{\mathrm{tot}}}$. Here ${f}_{\mathrm{rms}}$ is the average fractional rms amplitude of the oscillation during the observation and ${N}_{\mathrm{tot}}$ is the total number of detected counts.

The dependence of the waveform harmonic content on the rotational frequency of the star, the inclinations of the spot and the observer, and ${\mathcal{R}}$ are illustrated by the results presented in Table 3 (compare Table 4 of Lo et al. 2013). The amplitudes of the higher harmonics are substantially smaller in case 3(b) than in case 3(a) because the lower rotational frequency in the former case produces a line-of-sight velocity a factor of two smaller than in case 3(a). Comparing case 3(c) with case 3(a) demonstrates the sensitivity of the harmonic content of the waveform to the inclinations of the spot and the observer. Even though the line-of-sight velocity in case 3(c) is 95% of that in case 3(a), because the larger radius of the star in case 3(c) almost compensates for the higher inclinations of the spot and observer, the higher inclinations reduce the amplitudes of all the harmonics in case 3(c) compared to case 3(a). The ∼20% reduction in the total rms amplitude in case 3(c) reduces ${\mathcal{R}}$ by a comparable amount. The lower amplitudes of the higher harmonics in cases 3(b) and (c) increase the uncertainties $\delta M$ and $\delta {R}_{\mathrm{eq}}$ in these cases (see the final columns of Table 3).

Table 3.  Waveform Harmonic Content, ${\mathcal{R}}$ Values, and Uncertainties in M and ${R}_{\mathrm{eq}}$ ^a

Figure	${\nu }_{\mathrm{rot}}$	${\theta }_{\mathrm{spot}}$	${\theta }_{\mathrm{obs}}$	${f}_{\mathrm{rms}1}$ b	${f}_{\mathrm{rms}2}$ c	${f}_{\mathrm{rms}3}$ d	${f}_{\mathrm{rms}}$ e	${\mathcal{R}}$ f	$\delta {R}_{\mathrm{eq},1}$	$\delta {M}_{1}$
(case)	(Hz)	(deg)	(deg)	(%)	(%)	(%)	(%)		(%)	(%)
3(a)	600	90	90	10.0	3.9	1.4	11	479	2.9	2.8
3(b)	300	90	90	9.5	2.7	0.6	10	440	7.6	5.6
3(c)	600	60	60	8.1	2.6	0.7	9	377	6.7	6.5
4	600	90	90	10.0	3.8	1.3	11	479	3.6	3.2
5(b)	600	90	90	9.9	3.8	1.3	11	471	2.3	2.5

Notes.

^aSee the notes to Table 2 for the definitions of ${\nu }_{\mathrm{rot}}$, ${\theta }_{\mathrm{spot}}$, ${\theta }_{\mathrm{obs}}$, $\delta {R}_{\mathrm{eq},1}$, and $\delta {M}_{1}$. As noted there, the uncertainties listed here somewhat understate the actual constraints on $M({R}_{\mathrm{eq}})$, because $M/{R}_{\mathrm{eq}}$ is usually better constrained than M or ${R}_{\mathrm{eq}}$ considered separately. All the synthetic waveform data were generated assuming $M=1.6\;{M}_{\odot }$. ^bFractional rms amplitude of the first harmonic (fundamental) component of the synthetic waveform. ^cFractional rms amplitude of the second harmonic component of the synthetic waveform. ^dFractional rms amplitude of the third harmonic component of the synthetic waveform. ^eFractional rms amplitude of the total variation of the synthetic waveform. ^fSynthetic observed waveform figure of merit ${\mathcal{R}}\equiv 1.4{f}_{\mathrm{rms}}\sqrt{{N}_{\mathrm{tot}}}$; see Equation (1).

Download table as:  ASCIITypeset image

The dependence of $\delta M$ and $\delta {R}_{\mathrm{eq}}$ on the rotational frequency of the star and the inclinations of the spot and the observer are also illustrated by the trends in the confidence regions for M and ${R}_{\mathrm{eq}}$ in Figures 3(a)–(c) (see also Section 4.2.3 of Lo et al. 2013). As noted above, the lower rotational frequency in case 3(b) produces a smaller line-of-sight velocity, and hence a larger confidence region for M and ${R}_{\mathrm{eq}}$ than in case 3(a). The smaller inclinations of the spot and the observer in case 3(c) produce a larger confidence region for M and ${R}_{\mathrm{eq}}$ than in case 3(a), even though the line-of-sight velocity in case 3(c) is almost the same as in case 3(a).

These results and those of Lo et al. (2013) show that $\delta M$ and $\delta {R}_{\mathrm{eq}}$ are most sensitive to the colatitude of the hot spot and the inclination of the observer. They are sensitive, but less so, to the radius and the rotational frequency of the star. They depend more weakly on the background.

Three of the cases featured in Table 3 (cases 3(a), 4, and 5(b)) analyze different realizations of synthetic waveforms generated using identical values of the OS waveform parameters. Comparing these cases therefore provides insight into the sizes of the statistical and sampling errors in our results. The differences between the amplitudes of the harmonics in these three waveforms reflect the different shapes of the waveforms produced by fluctuations in the number of counts in each energy and phase bin; they indicate that the variations in the waveforms caused by these fluctuations are ≲2%. The fractional differences in the 1σ uncertainties in M and ${R}_{\mathrm{eq}}$ are much larger and are probably caused by the limitations of our sampling of the parameters of the model waveform during the fitting process; these differences indicate that the sampling errors are ∼20% in $\delta {R}_{\mathrm{eq}}$ and ∼10% in $\delta M$.

As a specific example, suppose that (1) ${\nu }_{\mathrm{rot}}=600$ Hz; (2)${\theta }_{\mathrm{spot}}=90^\circ$ and ${\theta }_{\mathrm{obs}}=90^\circ$; (3) the number of counts from the hot spot is equal to the number of background counts, producing a fractional rms amplitude of ∼54%; (4) the total number of detected counts is $\sim 2\times {10}^{6}$; and (5) the values of all the parameters in the waveform model except M and ${R}_{\mathrm{eq}}$ are known independently of the waveform-fitting procedure. The uncertainties in the estimates of M and ${R}_{\mathrm{eq}}$ would then be ∼2% in M and ∼1% in ${R}_{\mathrm{eq}}$ (see Section 4.2.1 and Table 2 of Lo et al. 2013).

Realistically, in addition to M and ${R}_{\mathrm{eq}}$, the values of some or all of the other parameters in the waveform model will have to be determined as part of the waveform-fitting procedure. Determining these additional parameters as part of the waveform-fitting process produces larger uncertainties in the estimates of M and ${R}_{\mathrm{eq}}$. The reason is that the effects on the waveform of changing different parameters in the waveform model are often very similar (see Section 4.2.2 of Lo et al. 2013 for a detailed discussion of these degeneracies). These degeneracies with respect to changes in the values of waveform model parameters are an inherent property of any physical model based on a rotating hot spot and cannot be removed by improving the model.

The number of background counts in observed burst oscillation waveforms is expected to be much greater than the number of counts collected from the hot spot (see Section 2.1.1 above and Section 2.2.1 of Lo et al. 2013). A large background increases the effects of the parameter degeneracies, because it increases the statistical fluctuations in the observed waveform. As a result, a wider range of model waveforms will adequately fit the waveform data. If the number of background counts is much greater than the number of oscillating counts, and the number of oscillating counts and the geometry remain unchanged, the uncertainties in M and ${R}_{\mathrm{eq}}$ increase with the number of background counts as ${{\mathcal{R}}}^{-1}\propto \sqrt{{N}_{\mathrm{back}}}$ (see Section 2.1.1).

The uncertainties in M and ${R}_{\mathrm{eq}}$ are much more sensitive to the inclinations of the hot spot and observer than to the background. This can be seen by comparing cases 2(a), (c), and (e) of Lo et al. (2013), which have the same inclinations but very different backgrounds. The relatively modest effect of an unknown background on the uncertainties in M and ${R}_{\mathrm{eq}}$ can be seen by comparing cases 5(a) and (b) of Lo et al. (2013), which assume the background is known exactly, with their corresponding cases 2(c) and (d), which assume the background is unknown and must be determined as part of the waveform-fitting procedure.

As an example of a realistic situation, suppose that (1)${\nu }_{\mathrm{rot}}=600$ Hz; (2) ${\theta }_{\mathrm{spot}}=90^\circ$ and ${\theta }_{\mathrm{obs}}=90^\circ$; (3) the number of background counts is 9 times the number of counts from the hot spot, producing a fractional rms amplitude of 11%; (4) the total number of detected counts is 10⁷; and (5) the values of all the parameters in the waveform model, including M and ${R}_{\mathrm{eq}}$, must be determined as part of the waveform-fitting procedure. Then the uncertainties in the derived estimates of M and ${R}_{\mathrm{eq}}$ would be ∼3% (see Table 2).

Independent knowledge of some of the system parameters can reduce or eliminate degeneracies, reducing the uncertainties in M and ${R}_{\mathrm{eq}}$. For example, accurate a priori knowledge of the observer's inclination can significantly improve the constraints, if the spot and observer inclinations are high; a priori knowledge of the distance to the system can also help (see Lo et al. 2013, Section 4.2.5).

Psaltis et al. (2014) proposed a simple formula for estimating the fractional uncertainty in estimates of ${R}_{\mathrm{eq}}$ obtained by fitting waveform models to waveform data. This formula is

$$\begin{eqnarray}&&\displaystyle \frac{\delta {R}_{\mathrm{eq}}}{{R}_{\mathrm{eq}}}={\left[\left(\displaystyle \frac{4\pi {\nu }_{\mathrm{rot}}{R}_{\mathrm{eq}}}{c}\right)\mathrm{sin}{\theta }_{\mathrm{spot}}\mathrm{sin}{\theta }_{\mathrm{obs}}\right]}^{-1}\left(\displaystyle \frac{\sqrt{S+B}}{{c}_{1}\;S}\right),\end{eqnarray} \tag{ 15 }$$

where S is the number of counts from the hot spot, B is the number of background counts, and c₁ is the fractional rms amplitude of the fundamental (first harmonic) in the waveform, defined in terms of S; the other symbols have their previous meanings. This formula is based on the following assumptions and approximations: (1) a single-energy or bolometric analysis is adequate; (2) the hot spot is infinitesimal in extent; (3) the distortion of the waveform produced by the surface rotational velocity is the dominant source of a nonzero second harmonic in the waveform; (4) the surface rotational velocity contributes a second harmonic with amplitude ${C}_{2}=2({v}_{\mathrm{los}}/c)\;{C}_{1}$; (5) ${\nu }_{\mathrm{rot}}$, ${R}_{\mathrm{eq}}$, $\mathrm{sin}{\theta }_{\mathrm{spot}}$, and $\mathrm{sin}{\theta }_{\mathrm{obs}}$ are all known independently of the waveform-fitting process.

Several of the assumptions and approximations on which Equation (15) is based can be questioned for X-ray burst oscillations: (1) the use of energy-resolved waveforms provides better results than single-energy or bolometric waveforms (see, e.g., Lo et al. 2013; Psaltis et al. 2014); (2) it will probably be necessary to use data from the peaks and/or tails of bursts, when the hot spot is not infinitesimal in extent (see Section 2.1.1); (3) the distortion of the waveform produced by the surface rotational velocity is not the dominant source of a nonzero second harmonic for the radiation beaming pattern expected for X-ray burst oscillations (see above and Section 2.1.3); (4) the assumed coefficient of proportionality between C₂ and C₁ is not the one derived by Poutanen & Beloborodov (2006), which depends on the spectrum of the emission (see their Equation (66)); (5) in most cases, some or all of the parameters ${R}_{\mathrm{eq}}$, $\mathrm{sin}{\theta }_{\mathrm{spot}}$, and $\mathrm{sin}{\theta }_{\mathrm{obs}}$ will not be known a priori, and will have to be determined by the waveform-fitting process. Assuming that the values of these parameters are known a priori completely eliminates the degeneracies among them, and therefore greatly overestimates the precision with which ${R}_{\mathrm{eq}}$ can realistically be determined using a given data set. (Although the fractional rms amplitude c₁ of the first harmonic component of the waveform, the number of counts S from the hot spot, and the number of background counts B are not directly observable (see Section 2.1.2) and would have to be determined by waveform fitting, the total number of counts ${N}_{\mathrm{tot}}=S+B$ and the number of counts ${c}_{1}S$ in the first harmonic component of the waveform are directly observable.)

How much do the simplifying assumptions made in deriving Equation (15) affect the computed uncertainty in ${R}_{\mathrm{eq}}$? To investigate this, we assume that the values of all the system parameters needed to evaluate Equation (15) are somehow known without doing any waveform fitting and evaluate Equation (15) using the values of these parameters that we used to generate the synthetic waveforms in our cases 3(a)–(c). We then compare the values of $\delta {R}_{\mathrm{eq}}$ given by Equation (15) with the 1σ uncertainties in ${R}_{\mathrm{eq}}$ we obtained by fitting our standard OS waveform model to the corresponding synthetic waveform data, using our Bayesian approach. In case 3(a), ${\nu }_{\mathrm{rot}}=600$ Hz, ${R}_{\mathrm{eq}}=11.8$ km, ${\theta }_{\mathrm{obs}}={\theta }_{\mathrm{spot}}=90^\circ$, $S={10}^{6}$, $B=9\times {10}^{6}$, and c₁ = 1.077 (recall that the fractional rms amplitude can exceed unity). Inserting these values into Equation (15) yields $\delta {R}_{\mathrm{eq}}/{R}_{\mathrm{eq}}=0.0099$, whereas our Bayesian statistical analysis yields $\delta {R}_{\mathrm{eq}}/{R}_{\mathrm{eq}}=0.029$. In case 3(b), ${\nu }_{\mathrm{rot}}=300$ Hz, ${R}_{\mathrm{eq}}=11.8$ km, ${\theta }_{\mathrm{obs}}={\theta }_{\mathrm{spot}}=90^\circ$, $S={10}^{6}$, $B=9\times {10}^{6}$, and c₁ = 0.990. Inserting these values into Equation (15) yields $\delta {R}_{\mathrm{eq}}/{R}_{\mathrm{eq}}=0.022$, whereas our Bayesian analysis yields $\delta {R}_{\mathrm{eq}}/{R}_{\mathrm{eq}}=0.076$. Finally, in case 3(c), ${\nu }_{\mathrm{rot}}=600$ Hz, ${R}_{\mathrm{eq}}=15$ km, ${\theta }_{\mathrm{obs}}={\theta }_{\mathrm{spot}}=60^\circ$, $S={10}^{6}$, $B=9\times {10}^{6}$, and c₁ = 0.856. Inserting these values into Equation (15) yields $\delta {R}_{\mathrm{eq}}/{R}_{\mathrm{eq}}=0.013$, whereas our Bayesian analysis yields $\delta {R}_{\mathrm{eq}}/{R}_{\mathrm{eq}}=0.067$. Thus, in these cases the uncertainties in ${R}_{\mathrm{eq}}$ obtained by fitting our standard OS waveform model to synthetic waveform data using our Bayesian statistical approach are ∼3–5 times larger than the uncertainties given by Equation (15). In reality, the values of all the system parameters needed to evaluate Equation (15) usually will not be known.

Equations (14) and (15) suggest that the uncertainty in ${R}_{\mathrm{eq}}$ should be smaller for larger values of ${\nu }_{\mathrm{rot}}$, ${R}_{\mathrm{eq}}$, $\mathrm{sin}{\theta }_{\mathrm{spot}}$, and $\mathrm{sin}{\theta }_{\mathrm{obs}}$, other things being equal. Our waveform-fitting results (see Table 2) are qualitatively consistent with this behavior but, as discussed above, reveal important quantitative differences from the scaling implied by these equations. This is not surprising, because Equations (14) and (15) do not include some important effects, such as occultation of the hot spot by the star. This can occur when the spot is near the rotational equator, and if it does, it will reduce further the uncertainties in M and ${R}_{\mathrm{eq}}$. Thus, there is no assurance that the actual uncertainty in ${R}_{\mathrm{eq}}$ will scale with ${\nu }_{\mathrm{rot}}$, ${R}_{\mathrm{eq}}$, $\mathrm{sin}{\theta }_{\mathrm{spot}}$, and $\mathrm{sin}{\theta }_{\mathrm{obs}}$ as simply as suggested by Equations (14) and (15). Consequently, choosing targets for observation and planning observational campaigns should be done using uncertainties in M and ${R}_{\mathrm{eq}}$ obtained by fitting waveform models to appropriate synthetic waveform data.

Tables 2 and 3 show the results of our Bayesian analysis of the uncertainties in M and ${R}_{\mathrm{eq}}$ that can be obtained by fitting our standard OS waveform model to burst oscillation waveform data, represented here by synthetic waveform data generated using the OS approximation. Figures 3, 4 and 5 shows the 1σ, 2σ, and 3σ confidence contours in the M–${R}_{\mathrm{eq}}$ plane for these cases. We find that M and ${R}_{\mathrm{eq}}$ can both be estimated with 1σ uncertainties ≲7% if (1) the star's rotation rate is ≳600 Hz, (2) the hot spot is located at a colatitude ≳60°, (3) the star is observed at a rotational inclination ≳60°, (4) the oscillations have a fractional rms modulation ≳10%, and (5) ≳10⁷ total counts are collected from the star.