SYSTEMATIC UNCERTAINTIES IN THE SPECTROSCOPIC MEASUREMENTS OF NEUTRON-STAR MASSES AND RADII FROM THERMONUCLEAR X-RAY BURSTS. I. APPARENT RADII

We base our study on the X-ray burst catalog of Galloway et al. (2008a). We chose 12 out of the 48 sources in the catalog based on the following criteria:

• 1.  
  We considered sources that show at least two bursts with evidence for photospheric radius expansion, based on the definition of the latter used by Galloway et al. (2008a). This requirement arises from our ultimate aim, which is to measure both the mass and the radius of the neutron star in each system using a combination of spectroscopic phenomena (as in, e.g., Özel et al. 2009).
• 2.  
  We excluded dippers, accretion disk corona sources, or known high-inclination sources. This list includes EXO 0748−676, MXB 1659−298, 4U 1916−05, GRS 1747−312, 4U 1254−69, and 4U 1710−281, for which it was shown that geometric effects related to obscuration or reflection significantly affect the flux from the stellar surface that is measured by a distant observer (Galloway et al. 2008b).
• 3.  
  We did not consider the known millisecond pulsars SAX J1808.4−3658 and HETE J1900.1−2455 because the presence of pulsations in their persistent emission implies that their magnetic fields are dynamically important and, therefore, may affect the properties of X-ray bursts.
• 4.  
  We excluded the sources GRS 1741.9−2853 and 2E 1742.9−2929, as well as a small number of bursts from Aql X-1, 4U 1728−34, and 4U 1746−37, for which there is substantial evidence that their emission is affected by source confusion (Galloway et al. 2008a; Keek et al. 2010).
• 5.  
  For each source, we considered only bursts for which the flux in the persistent emission prior to the burst is at most 10% of the inferred Eddington flux for each source, i.e., γ ≡ F_per/F_Edd < 0.1, as calculated by Galloway et al. (2008a). Because we subtract the pre-burst persistent emission from the decay spectrum of each X-ray burst, this requirement reduces substantially the systematic uncertainties introduced by potential changes in the emission from the accretion flow during the burst.

Table 1 lists all the X-ray bursters that fulfill the above requirements, together with the number of bursts observed by RXTE for each source. This is the complete list of sources for which the masses and radii can be measured, in principle, using the spectroscopic method of Özel et al. (2009), with currently available data. Our analyses of the bursts from EXO 1745−248, 4U 1608−52, and 4U 1820−30 have been reported elsewhere (Özel et al. 2009; Güver et al. 2010a, 2010b) and will not be repeated here.

We analyzed the burst data for the 12 sources shown in Table 1 following the method outlined in Galloway et al. (2008a; see also Özel et al. 2009; Güver et al. 2010a, 2010b).

For each source, we extracted time-resolved 2.5–25.0 keV X-ray spectra using the seextrct ftool for the science event mode data and the saextrct ftool for the science array mode data from all of the RXTE/Proportional Counter Array (PCA) layers. Science mode observations provide high count-rate data with a nominal time resolution of 125 μs in 64 spectral channels over the whole energy range (2–60 keV) of the PCA detector. Following Galloway et al. (2008a), we extracted spectra integrated over 0.25, 0.5, 1, and 2 s time intervals, depending on the source count rate during the burst, so that the total number of counts in each spectrum is roughly constant. (In a few cases, data gaps during the observations result in X-ray spectra integrated over shorter exposure times.) We took the spectrum over a 16 s time interval prior to the onset of each burst as the spectrum of the persistent emission, which we subtracted from the burst spectra as background.

We generated separate response matrix files for each burst using the PCARSP version 11.7, HEASOFT release 6.7, and HEASARC's remote calibration database and took into account the offset pointing of the PCA during the creation of the response matrix files. This latest version of the PCA response matrix makes the instrument calibration self-consistent over the PCA lifetime and yields a normalization of the Crab pulsar that is within 1σ of the calibration measurement of Toor & Seward (1974) for that source. In Section 4, we discuss in some detail the effect of the uncertainties in the absolute flux calibration on the measurement of the apparent surface area of neutron stars. Finally, we corrected all of the X-ray spectra for PCA dead time following the method suggested by the RXTE team.²

To analyze the spectra, we used the Interactive Spectral Interpretation System, version 1.4.9-55 (Houck & Denicola 2000). For each fit, we included a systematic error of 0.5% as suggested by the RXTE calibration team.³

We fit each spectrum with a blackbody function using the bbodyrad model (as defined in XSPEC; Arnaud 1996) and multiplied it by the tbabs model (Wilms et al. 2000) that takes into account the interstellar extinction, assuming interstellar medium (ISM) abundances. The model of the X-ray spectrum for each burst, therefore, depends on only three parameters: the color temperature of the blackbody, T_c, the normalization of the blackbody spectrum, A, and the equivalent hydrogen column density, N_H, of the interstellar extinction.

Allowing the hydrogen column density N_H to be a free parameter in the fitting procedure leads to correlated uncertainties between the amount of interstellar extinction and the temperature of the blackbody. Moreover, for every 10²² cm⁻² overestimation (underestimation) of the column density, the inferred flux of the thermal emission is systematically larger (smaller) by ≈10⁻⁹ erg cm⁻² s⁻¹ (Galloway et al. 2008a). Reducing the uncertainties related to the interstellar extinction for each burster is, therefore, very important in controlling systematic effects.

A recent analysis of high-resolution grating spectra from a number of X-ray binaries (Miller et al. 2009) showed that the individual photoelectric absorption edges observed in the X-ray spectra do not show significant variations with source luminosity or spectral state. This result strongly suggests that the neutral hydrogen column density is dominated by absorption in the ISM and does not change on short timescales. In order to reduce this systematic uncertainty and given the fact that there is no evidence of variable neutral absorption for each system, we fixed the hydrogen column density for each source to a constant value that we obtained in one of the following ways: (1) For a number of sources, the equivalent hydrogen column density was inferred independently using high-resolution spectrographs. (2) In cases where only an optical extinction or reddening measurement exists, we used the relation given by Güver & Özel (2009) to convert it to the equivalent hydrogen column density. (3) Finally, for three sources (4U 1702−429, KS 1731−260, SAX J1750.8−2900) there are no independent hydrogen column density measurements published in the literature. In these cases, we fit the X-ray spectra obtained during all the X-ray bursts of each source allowing the N_H value to vary in a wide range. We then found the resulting mean value of N_H and used this as a constant in our second set of fits to all the X-ray bursts. In Table 1, we show the adopted values for the hydrogen column density together with the measurement uncertainties, the method with which the values were estimated, and the appropriate references. Future observations of these sources with X-ray grating spectrometers on board Chandra and XMM-Newton can help decrease the uncertainty arising from the lack of knowledge of the properties of the interstellar matter toward these sources.

Our goal in this article is to study potential systematic uncertainties in the measurement of the apparent area of each neutron star during the cooling tails of X-ray bursts. Hereafter, we adopt the following working definition of the cooling tail. It is the time interval during which the inferred flux is lower than the peak flux of the burst or the touchdown flux for photospheric radius-expansion bursts. For the purpose of this definition, we use as a touchdown point the first moment at which the blackbody temperature reaches its highest value and the inferred apparent radius is lowest. In order to control the count-rate statistics, we also set a lower limit on the thermal flux during each cooling tail of 5 × 10⁻⁹ erg s⁻¹ cm⁻² (or 5 × 10⁻¹⁰ for the exceptionally faint sources 4U 1746−37 and 4U 0513−401).

In Sections 3 and 4, we discuss in detail our approach of quantifying the systematic uncertainties in the measurements of the apparent radii during the cooling tails of thermonuclear bursts, using the sources KS 1731−260, 4U 1728−34, and 4U 1636−536 as case studies. For another analysis of the cooling tails of X-ray bursts from 4U 1636−536, see Zhang et al. (2011). In Section 5, we repeat this procedure systematically for seven additional sources from Table 1. Finally, the cooling tails of the bursts observed from 4U 0513−40 and Aql X-1 show irregular behavior, and we report our analysis of them in Appendixes C and D.

In order to quantify the degree to which any of these effects influence the spectra of X-ray bursters, we fit all cooling-tail spectra from each source with a blackbody model, allowing for a level of systematic uncertainty σ_sys, determined in the following way. We reduced the degree of freedom in this procedure by setting σ_sys in each spectral bin to be a constant fraction ξ of the Poisson error σ_formal and fixed the value of ξ to a constant for each source. We then added this systematic uncertainty to the formal Poisson error of each measurement in quadrature, i.e., σ² = σ²_formal + σ_sys², and defined a statistic X² such that

$$\begin{equation} X^2\big(\sigma _{\rm sys}^2\big)= \frac{1}{1+\xi ^2} \chi ^2\;. \end{equation} \tag{ 2 }$$

Here, χ² is the standard statistic calculated for each fit. The posterior distribution of the new statistic X² may not be the same as that of the χ² statistic and will depend on the nature of the systematic uncertainties. In order to explore the degree and nature of systematic uncertainties, we fit for each source the X² distribution with the χ² distribution expected for the number of degrees of freedom used in the fits, with ξ as a free parameter. Note that the expected distribution that we use is formally correct if the count-rate data in the individual spectral bins have uncorrelated errors.

The expected χ² and the observed X² distributions for the optimal ξ value for KS 1731−260, 4U 1728−34, and 4U 1636−536 are shown in Figure 1. The ξ value required to make the observed X² distributions for KS 1731−260 consistent with the expected ones is 0.55 (see Table 2). This value corresponds to systematic uncertainties that are very small. As an example, we consider a typical 0.25 s integration for KS 1731−260, when its spectrum is characterized by a color temperature of 2.5 keV. In this case, the average number of counts in each spectral bin in the 3–10 keV energy range is 110, and the corresponding Poisson uncertainty is ≃ 9.5%. Multiplying this by the inferred value ξ = 0.52 for this source leads to the conclusion that the systematic uncertainties required to render the observed spectrum consistent with a blackbody are ≃ 5%. For the case of 4U 1728−34, which is brighter than KS 1731−260, the formal uncertainties are ≃ 6% and the systematic uncertainties amount to ≃ 3%.

Zoom In Zoom Out Reset image size

Figure 1 shows that the resulting X² distributions can be well approximated by the χ² distribution but have weak tails extending to higher values. These tails most likely arise from the inconsistency of only a small number of spectra with blackbody functions even though we cannot rule out the possibility that the X² statistic follows a different posterior distribution than χ².

Figure 1 allows us also to identify the X-ray spectra that are statistically inconsistent with blackbodies. For each source, there is a maximum value of X² per degree of freedom beyond which the distribution of X² values deviates from the theoretical expectation. For the case of KS 1731−260, this limiting value of the reduced X² is 1.7, for 4U 1728−34 it is 1.7, and for 4U 1636−536 it is 1.9 (see Table 2). The spectra with high reduced X² values often occur at the late stages of the cooling tails, where the subtraction of the persistent emission is the most problematic. However, unacceptable X² values may also be found in other, seemingly random places of the cooling tails. Nevertheless, the fraction of spectra that is inconsistent with blackbodies is ≲ 3%, ≲ 6%, and ≲ 2% for KS 1731−260, 4U 1728−34, and 4U 1636−536, respectively. Hereafter, we consider only the spectra that we regard to be statistically acceptable, given the values of the reduced X².

Figure 3 (left panels) shows the dependence of the emerging flux on color temperature for all the spectra in the cooling tails of KS 1731−260, 4U 1728−34, and 4U 1636−536 that we consider to be statistically acceptable (see Section 3). We chose these three sources to use as detailed examples because of the high number spectra obtained for each and the fact that they span the range of behavior in cooling tracks that we will discuss below. If the whole neutron star is emitting as a single-temperature blackbody and the color correction factor is independent of color temperature, then F_cool should scale as T⁴_c. Our aim here is to investigate the systematic uncertainties on the measurement of the apparent surface area in each source at different flux levels. We, therefore, divided the data into a number of flux bins and plotted in the same figure (right panels) the distribution of the blackbody normalization values for some representative bins. The blackbody normalization for each spectrum is defined formally as A ≡ F_cool/σ_SBT⁴_c, although in practice this is one of the two measured parameters and the flux is derived from the above definition. According to Equation (3), the blackbody normalization is equal to A = f⁻⁴_cR_app²/D², and we report it in units of (km/10 kpc)². If we do not correct for the dependence of the color correction factor on the flux, we expect the normalization A to show a dependence on color temperature that is the mirror symmetric of that shown in Figure 2.

Zoom In Zoom Out Reset image size

The flux–temperature diagrams of KS 1731−260, 4U 1728−34, and 4U 1636−536 share a number of similarities but are also distinguished by a number of differences. In all three cases, the vast majority of data points lie along a very well defined correlation. This reproducibility of the cooling curves of tens of X-ray bursts per source, combined with the lack of large-amplitude flux oscillations during cooling tails of bursts, provides the strongest argument that the thermal emission engulfs the entire neutron-star surface with very small temperature variations at different latitudes and longitudes.

In 4U 1728−34, a deviation of the data points from the F_cool ∼ T⁴_c correlation is evident at high fluxes and color temperatures (T_c ⩾ 2.5 keV), which may be due to the evolution of the color correction factor at high Eddington fluxes (see the discussion in Section 4 and Figure 2). The same deviation is not evident in KS 1731−260 or 4U 1636−536, for which the highest temperatures encountered in the cooling tails were ≲ 2.5 keV.

Finally, in all three sources, a number of outliers exist at the lowest flux levels, with normalizations that deviate from the above correlation. In KS 1731−260 and 4U 1636−536, the outliers correspond to higher normalization values, whereas in 4U 1728−34 they correspond to lower normalization values with respect to the majority of the data points.

Any combination of the effects discussed earlier in this and in the previous section may be responsible for the outliers. Non-uniform cooling of the neutron-star surface will lead to a reduced inferred value for the apparent surface area. Reflection of the surface X-ray emission off a geometrically thin accretion disk will cause an increase in the inferred value for the apparent surface area. Finally, Comptonization of the surface emission in a corona may have either effect, depending on the Compton temperature of the radiation.

Our goal in this work is not to understand the origin of the outliers, but rather to ensure that their presence does not introduce any biases in the measurements of the apparent surface areas inferred from the vast majority of spectra. Indeed, including the outliers in a formal fit of the flux–temperature correlation will cause a systematic increase of the apparent surface area with decreasing flux in KS 1731−260 and in 4U 1636−536 and a systematic decrease in 4U 1728−34 (cf. Damen et al. 1990; Bhattacharyya et al. 2010). This issue was largely avoided in our earlier work (Özel et al. 2009; Güver et al. 2010a, 2010b) by considering only the relatively early time intervals in the cooling tail of each burst, which correspond only to the brightest flux bins. In order to go beyond this limitation here, we consider all spectral data for each burst and employ a Bayesian Gaussian mixture algorithm, which is a standard procedure for outlier detection in robust statistics (e.g., Titterington et al. 1985; McLachlan & Peel 2000; Huber & Ronchetti 2009).

Our working hypothesis that the main peak of normalizations in each flux interval corresponds to the signal and the remainder are outliers is shaped by two aspects of the observations. First, at high fluxes, each histogram can be described by a single Gaussian with no evidence or room for a second distribution of what we would call outliers. At lower fluxes, when the histograms can be decomposed into two Gaussians, the distribution with the highest peak has properties that smoothly connect to those of the single Gaussians at higher fluxes. Second, the Gaussians that we call our signal always contain the majority of the data points compared with the distribution of what we call the outliers. We take these as our criteria for defining our signal.

The algorithm, which we describe below in detail, allows us also to measure the degree of intrinsic variation in the apparent surface area for each source that is consistent with the width of the flux–temperature correlation. Even though we compare spectra in relatively narrow flux bins, the different observing modes and the different number of PCUs used in each observation result in a wide range of count rates and, hence, in a wide range of formal errors. This necessitates the use of the Bayesian approach we describe below as opposed to, e.g., performing χ² fits of the distributions over blackbody normalization.

Consider first the situation in which there are no systematic uncertainties and the inferred blackbody normalization from every spectrum reflects the apparent surface area of the entire neutron star. If this were the case, the intrinsic distribution of blackbody normalizations would be a delta function, while the observed distribution would be broader, with a width equal to the average formal uncertainties of the measurements.

In a more realistic situation, there is an intrinsic range of the inferred blackbody normalizations, caused by a combination of the various effects discussed earlier. We assume that, in each flux bin, the intrinsic distribution of normalization values is the Gaussian:

$$\begin{equation} P_{\rm int}(A| A_{\rm int}, \sigma _{\rm int}) =\frac{1}{\sqrt{2\pi \sigma _{\rm int}^2}} \exp \left[-\frac{(A-A_{\rm int}^2)}{2\sigma _{\rm int}^2}\right]\;, \end{equation} \tag{ 4 }$$

where A_int is its mean and σ_int is its standard deviation. The observed distribution of blackbody normalizations will be again broader than the intrinsic distribution because of the formal uncertainties in each measurement. Moreover, at low flux levels, a small number of outliers exist, which skews and introduces tails to the observed distribution. Our goal here is to quantify the most probable value of the blackbody normalization, i.e., the parameter A_int in the distribution (Equation (4)), as well as the intrinsic range of values in each flux bin, i.e., the parameter σ_int, while excluding the outliers.

When the inspection of an observed distribution requires an outlier detection scheme (as, e.g., in the lower right panels of Figure 3), we model the distribution of outliers by a similar Gaussian, i.e.,

$$\begin{equation} P_{\rm out}(A| A_{\rm out}, \sigma _{\rm out}) =\frac{1}{\sqrt{2\pi \sigma _{\rm out}^2}} \exp \left[-\frac{(A-A_{\rm out}^2)}{2\sigma _{\rm out}^2}\right]\;, \end{equation} \tag{ 5 }$$

with a mean A_out and a standard deviation σ_out. Our model distribution function is, therefore, in general the Gaussian mixture:

$$\begin{eqnarray} P_{\rm model}(A| A_{\rm int}, \sigma _{\rm int}, A_{\rm out}, \sigma _{\rm out},\eta)&=& P_{\rm int}(A| A_{\rm int}, \sigma _{\rm int})\nonumber\\ &&+\eta P_{\rm out}(A| A_{\rm out}, \sigma _{\rm out}),\quad\quad \end{eqnarray} \tag{ 6 }$$

where η measures the relative fraction of outliers in our sample.

We are looking for the parameters in these distributions that maximize the posterior probability P(A_int, σ_int, A_out, σ_out|data), which measures the likelihood that a particular set of parameters is consistent with the data. Using Bayes' theorem, this probability distribution is

$$\begin{eqnarray} &&P(A_{\rm int}, \sigma _{\rm int}, A_{\rm out}, \sigma _{\rm out},\eta \vert {\rm data})=C P({\rm data}\vert A_{\rm int}, \sigma _{\rm int}, A_{\rm out},\nonumber\\ &&\quad\quad \sigma _{\rm out},\eta) P(A_{\rm int})P(\sigma _{\rm int}) P(A_{\rm out}) P(\sigma _{\rm out}) P(\eta)\;, \end{eqnarray} \tag{ 7 }$$

where P(data|A_int, σ_int, A_out, σ_out, η) is the posterior probability that the particular set of normalization data have been measured for a given set of parameters and the last five terms are the prior probabilities for the model parameters; the constant C ensures that the probability density is normalized. We consider a flat prior probability distribution for the central value of each Gaussian in an interval that extends beyond the observed range of normalization values for each source. We also consider a flat prior probability distribution for the standard deviation for each Gaussian from 10⁻³ (km/10 kpc)² (which is practically zero) to a value equal to the width of the range of normalizations. Finally, we take a flat prior distribution of η between zero and one. Our results are extremely weakly dependent on the ranges of parameter values we considered.

Given a set of N observed normalization values A_i and their uncertainties σ_i, we calculate the posterior probability of the data given a set of model parameters as

$$\begin{eqnarray} &&P({\rm data}\vert A_{\rm int}, \sigma _{\rm int}, A_{\rm out} \sigma _{\rm out},\eta)=C _2 \prod _{i=1}^N\int dA \frac{1}{\sqrt{2\pi \sigma _i^2}}\nonumber\\ &&\quad\quad\quad\times\exp \left[-\frac{(A-A_{\rm i}^2)}{2\sigma _{\rm i}^2}\right] [P_{\rm int}(A| A_{\rm int}, \sigma _{\rm int})\nonumber\\ &&\quad\quad\quad +\eta P_{\rm out}(A| A_{\rm out}, \sigma _{\rm out})]. \end{eqnarray} \tag{ 8 }$$

Here, A_i is the value of each measurement and σ_i is its corresponding uncertainty, which we assume to be normally distributed; C₂ is a normalization constant. Using this posterior probability distribution, we identify the most probable values of the model parameters, as well as their uncertainties, following standard procedures.

Figure 4 shows, as an example, the application of the Gaussian mixture algorithm to the distribution of blackbody normalizations obtained from fitting the spectra of 4U 1728−34, when the burst flux was in the range (1–2) × 10⁻⁸ erg cm⁻² s⁻¹. The histogram shows the distribution of blackbody normalizations measured in this flux interval. The parameters of the Gaussian mixture that maximize the posterior probability distribution calculated using Equation (7) correspond to an intrinsic distribution with a mean and standard deviation of A_int = 126.1 (km/10 kpc)² and σ_int = 13.4 (km/10 kpc)², respectively, and a distribution of outliers with a mean of A_out = 62.3 (km/10 kpc)², a standard deviation of σ_out = 6.1 (km/10 kpc)², and a relative normalization of η = 0.2. In order to compare the Gaussian mixture with the observed histogram, we first convolve it with a third Gaussian distribution with a standard deviation of σ_formal = 16.5 (km/10 kpc)² that is equal to the average formal errors of all measurements in this flux interval. The result is shown as a red solid line in Figure 4, which demonstrates that our assumption of a Gaussian distribution for the majority of the data and for the outliers is consistent with the observations.

Zoom In Zoom Out Reset image size

We can appreciate the importance of excluding the outliers from the sample using the Bayesian Gaussian mixture algorithm in the following way. If we repeat the above procedure for the same flux interval in 4U 1728−34 but do not allow for the presence of outliers, the parameters of the intrinsic distribution will be A_int = 118.6 (km/10 kpc)² and σ_out = 26.4 (km/10 kpc)², which are significantly different from the results quoted above. Moreover, had we assumed that there are no systematic errors but that the underlying distribution was a delta function, the most probable value would have been A_int = 117.1^+0.2_−0.8 (km/10 kpc)². The formal errors in such a measurement would have been substantially smaller than the systematic uncertainties.

In the following, we show in detail the application of this Gaussian mixture algorithm to all the flux intervals for the three sources.

Figure 5 shows the result of the outlier detection procedure as applied to the various flux intervals of KS 1731−260, 4U 1728−35, and 4U 1636−536 (see also Table 3). Each dot represents the most likely centroid of the intrinsic distribution of blackbody normalizations, while each error bar represents its most likely width.

Zoom In Zoom Out Reset image size

The normalization of the blackbody in the case of 4U 1728−34 shows a strong dependence on flux when the source is emitting at near-Eddington rates. This is not seen in KS 1731−260 and 4U 1636−536, for which the color temperature never reached temperatures above 2.5 keV. In fact, in our observed sample, the strong dependence of the normalization on the flux is seen in all three sources with color temperatures in excess of 2.5 keV, namely, 4U 1728−34, 4U 1702−429, and 4U 1724−307 (see Figure 6 below), but is absent from the other sources. The evolution of the blackbody normalization at near-Eddington fluxes depends on the local gravitational acceleration (g ∼ M/R²) on the neutron-star surface and on its composition. Therefore, modeling the decline of the blackbody normalization at large fluxes would allow us, in principle, to measure the combination M/R² of the mass and radius of the neutron star, as has been attempted by Majczyna & Madej (2005) and Suleimanov et al. (2011).

At low burst fluxes, the observed normalizations show a small albeit statistically significant trend toward lower values. In the case of 4U 1728−34, where the potential decrease is the largest, the normalization changed from 133.7 ± 16.4 (km/10 kpc)² to 114.0 ± 15.6 (km/10 kpc)² as the flux declined from 6 × 10⁻⁸ erg s⁻¹ cm⁻² to 0.5 × 10⁻⁸ erg s⁻¹ cm⁻². This corresponds to a ≲ 15% reduction in the normalization. In other sources, the decline is even smaller. This weak dependence seen in the data argues against the models with solar composition, as shown in Figure 2.

The decline in the normalization can, in principle, be accounted for with a ≲ 4% increase in the color correction factor (see Equation (3)) toward low temperatures. As can be seen in Figure 2, current atmosphere models do not predict such an evolution. However, as we discussed in the beginning of Section 4, a number of effects related to the physics of bursts on the neutron-star surface (uneven cooling and burst oscillations), as well as the reduced sensitivity of the PCA at low energies, are capable of causing the observed trend in the normalization. Moreover, the ≲ 4% effect in the data we would aim to model is comparable to the uncertainty in inferring theoretically the color correction factor from the atmosphere models and is also comparable to the deviation of the observed and theoretical spectra from blackbodies. The latter concern can be remedied by fitting directly theoretical model spectra to the data, but reducing the theoretical uncertainties present in the models to less than roughly a few percent is significantly more challenging.

An alternative approach is to allow for a range of values for the color correction factor that span the spread of theoretical uncertainties, metallicities, and fluxes of the cooling tails. In earlier work (Güver et al. 2010a, 2010b) we allowed for a 7% range in the color correction factor (from 1.3 to 1.4), which is adequate to account for the evolution in the blackbody normalizations shown by the data. Following this approach, we need to identify the range of blackbody normalizations we obtained from the data at different flux intervals as a systematic uncertainty in the measurement. In order to achieve this, we applied the Bayesian Gaussian mixture algorithm to the combined data set for each source for all flux intervals that correspond to a color temperature ⩽2.5 keV. This way, we computed the most probable value for the blackbody normalization for each source in a wide flux range, as well as the systematic uncertainties in that measurement, which account for the potential decline of the normalization with decreasing flux.

In Figure 5, we identify the flux intervals we used for each source with a filled circle in the middle of each error bar. We also depict with a solid line the most probable value of the normalization in this wide flux range and with dashed lines the range of systematic uncertainties. Our results are summarized in Table 3. The blackbody normalizations for KS 1731−260, 4U 1728−34, and 4U 1636−536 are 96.0 ± 7.9 (km/10 kpc)², 134.4 ± 14.9 (km/10 kpc)², and 130.7 ± 20.9 (km/10 kpc)², respectively, with the uncertainties dominated entirely by systematics.