New Constraints on the Nuclear Equation of State from the Thermal Emission of Neutron Stars in Quiescent Low-mass X-Ray Binaries

The surface emission from NSs in qLMXBs is now routinely used to obtain measurements of ${M}_{\mathrm{NS}}$ and ${R}_{\mathrm{NS}}$. While during outbursts, the accreted matter dominates the X-ray emission, the thermal emission from the surface of the NS becomes visible in quiescence. The source of this thermal emission is internal and originates from the heat deposited by nuclear reactions in the crust during accretion episodes (e.g., Haensel & Zdunik 2008). As this emission, reprocessed by the NS outer layers, is observed in the X-rays and modeled with realistic atmosphere models (Zavlin et al. 1996; Heinke et al. 2006; Ho & Heinke 2009; Haakonsen et al. 2012), one can measure the redshifted temperature and the size of the emission area. In this way, the X-ray spectra of qLMXBs provide a measurement of ${R}_{\infty }$, defined as

$$\begin{eqnarray}&&{R}_{\infty }={R}_{\mathrm{NS}}\ \left(1+z\right)={R}_{\mathrm{NS}}\ {\left(1-\displaystyle \frac{2{{GM}}_{\mathrm{NS}}}{{R}_{\mathrm{NS}}{c}^{2}}\right)}^{-1/2}.\end{eqnarray} \tag{ 1 }$$

This requires knowing the distance to the source, and qLMXBs located in GCs have provided ${R}_{\mathrm{NS}}$ measurements, since their distances can be independently and rather precisely measured (see Section 2.2).

The qLMXBs inside GCs also present another advantage of exhibiting a remarkable flux stability at all timescales (Heinke et al. 2006, 2014; Guillot et al. 2011; Servillat et al. 2012). While LMXBs in the field of the Galaxy often exhibit flux variability attributed to changes in the nonthermal and/or thermal components (e.g., Rutledge et al. 2002; Campana et al. 2004), which complicate the spectral modeling, the spectra of known qLMXBs located in GCs are purely thermal, without signs of nonthermal emission (e.g., Guillot et al. 2013). Overall, this reinforces the scenario in which we are observing the uncontaminated thermal cooling of NSs.

Another advantage of NSs in qLMXBs over other subgroups of NSs for the purpose of radius measurements is the relatively straightforward modeling of their emergent spectra. While the atmospheric composition of isolated NSs may be uncertain (e.g., Burwitz et al. 2003; Ho & Heinke 2009), the atmosphere of NSs in LMXBs consists of a single-composition layer of a fully ionized light element. Since the accreted matter settles gravitationally within 10–100 s (Alcock & Illarionov 1980; Bildsten et al. 1992), the outermost layer of a transiently accreting NS is thought to be composed of the lightest accreted element, usually hydrogen (H). Moreover, the magnetic fields of these old sources are thought to be weak, as supported by the fact that their presumed descendants, MSPs (Alpar et al. 1982; Bhattacharya & van den Heuvel 1991; Tauris & van den Heuvel 2006), have inferred dipole fields B ∼ 10⁸–${10}^{9}\,{\rm{G}}$, compared to 10¹¹–10¹² G for the younger, "classical" pulsars, which have not undergone accretion. Such low B fields do not affect the emergent spectrum, and it can therefore be assumed that the NS atmosphere is nonmagnetic. For these reasons, H-atmosphere models, and in some cases helium- (He-) atmosphere models (see below), have been used to fit the spectra of the NS in qLMXBs and extract measurements of ${M}_{\mathrm{NS}}$ and ${R}_{\mathrm{NS}}$.

It is generally accepted that the atmosphere of an NS in a qLMXB is composed of pure H, since the atmospheric composition would be that of the lightest element present in the companion star. Unless the companion is completely devoid of H, the matter transferred onto the NS will contain some H, and therefore the material present in the outermost layer will be H. Diffusive burning of H into He may happen in the hot photosphere, but this is expected to happen on timescales of 10³–10⁴ yr (Chang & Bildsten 2004; Chang et al. 2010), whereas the atmosphere (of thickness $\sim 1\,\mathrm{cm}$ and mass ${M}_{\mathrm{atm}}\sim {10}^{-20}\,{M}_{\odot }$; Bogdanov et al. 2016) can rapidly be replenished by H matter from the stellar companion, even at very low accretion rates of $\sim {10}^{-13}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$. More importantly, observational evidence demonstrated the presence of H in the qLMXB systems in 47 Tux X-5 (Edmonds et al. 2002) and 47 Tuc X-7 (Bogdanov et al. 2016) and in the GC ω Cen (Haggard et al. 2004). Searches for Hα emission at the position of the qLMXB in NGC 6397 were unsuccessful, only placing upper limits on the equivalent width of the spectral line and thus on the accretion rate (Heinke et al. 2014). It was therefore argued that this qLMXB was devoid of H, and the authors advocated for an He atmosphere instead. This conclusion was supported by the low ${R}_{\mathrm{NS}}$ found from the spectral analyses with an H atmosphere (∼8 km in the earlier work of Guillot et al. 2011, 2013), while an He atmosphere resulted in an ${R}_{\mathrm{NS}}$ value compatible with that of other NSs (Heinke et al. 2014). However, the stellar companion was only detected in the R band, and limits on its photometric colors made it compatible with both the white dwarf sequence and the main sequence of the host GC. However, as discussed below, the proper modeling of pileup (an instrumental effect; Davis 2001) in the Chandra X-ray Observatory spectra of this qLMXB is sufficient to yield radii in the range⁷ 10–11 km.

With all of these considerations in mind, qLMXBs located inside GCs are ideal objects that provide a well-understood scenario to measure the radii of NSs. As mentioned above, obtaining constraints on the EoS from qLMXBs requires combining them into a statistical analysis. Here we analyze the spectra of the qLMXB in the GCs M13 (NGC 6205), M28 (NGC 6266), M30 (NGC 7099), NGC 6304, NGC 6397, ω Cen (NGC 5139), and 47 Tuc (NGC 104) X-7. We excluded 47 Tuc X-5 because of its eclipses, flux variability, and variable line-of-sight absorption, which make the spectral modeling rather uncertain (Bogdanov et al. 2016). Some information about these sources is detailed in Table 1.

Table 1.  Observational Information on the Seven qLMXB Sources Considered in Our Analysis

Globular	R.A.a	Decl.a	XMM Exp.	Chandra Exp.	S/N	Groupb	Distances	Distances	Referencesc
Cluster Host	(J2000)	(J2000)	Time (ks)	Time (ks)			Dist. 1 (kpc)	Dist. 2 (kpc)
47 Tuc (X-7)	00:24:03.53	−72:04:52.2	0	181	122	A, A'	4.53 ± 0.08	4.50 ± 0.06	(1)
M28	18:24:32.84	−24:52:08.4	0	327	113	A, A'	5.5 ± 0.3	5.50 ± 0.13	(2), (3)
NGC 6397	17:40:41.50	−53:40:04.6	0	340	82	A, A'	2.51 ± 0.07	2.30 ± 0.05	(4)
ω Cen	13:26:19.78	−47:29:10.9	36	291	49	B, B'	4.59 ± 0.08	5.20 ± 0.09	(5)
M13	16:41:43.75	+36:27:57.7	29	55	36	B, A'	7.1 ± 0.62	7.10 ± 0.10	(6)
M30	21:40:22.16	−23:10:45.9	0	49	32	B, B'	8.2 ± 0.62	8.10 ± 0.12	(6)
NGC 6304	17:14:32.96	−29:27:48.1	0	97	28	B, B'	6.22 ± 0.26	5.90 ± 0.14	(7)

Notes. All distance uncertainties are given at the 1σ confidence level.

^aCoordinates of the qLMXB in each of the GCs. ^bGroups A and B denote the sources with a higher S/N ($\gt 60$) and lower S/N ($\lt 60$), respectively. Groups A' and B' denote the sources for which we obtain a peaked and flat posterior distribution of the NS mass, respectively (see Section 5 for more details). ^cGaia Collaboration et al. (2018c), from which the distances were obtained from the individual X, Y, Z coordinate values, as given in their Table C.3, using ${r}_{\mathrm{GC},\odot }=\sqrt{{X}^{2}+{Y}^{2}+{Z}^{2}}$.References: (1) Bogdanov et al. (2016); (2) Harris (2010); with uncertainties estimated in (3) Servillat et al. (2012); (4) Heinke et al. (2014); (5) Watkins et al. (2013); (6) O'Malley et al. (2017); (7) Recio-Blanco et al. (2005).

Download table as:  ASCIITypeset image

In this paper, we work with a set of distances obtained from a heterogeneous set of methods, including dynamical (Watkins et al. 2013) and photometric (other references in Table 1) distance measurements. In most cases, these are recent measurements or measurements discussed in previous qLMXB analyses, which we used for convenient comparison (e.g., Bogdanov et al. 2016). These distances used and their uncertainties are listed in Table 1 as Dist. 1.

To evaluate the impact of the choice of distances, we also considered distances from a more uniform set of measurements. The determination of accurate astrometric distances to large samples of GCs has now become a tangible reality, thanks to the exquisite data provided by the European Space Agency's Gaia space mission (Gaia Collaboration et al. 2016). Within the framework of Gaia's Data Release 2 (DR2; Gaia Collaboration et al. 2018b), trigonometric parallaxes have already become available for large numbers of stars belonging to dozens of GCs. Still, as discussed in detail by Pancino et al. (2017), and more recently also emphasized by Gaia Collaboration et al. (2018a), systematic uncertainties still preclude the determination of reliable distances based on the available Gaia data for such crowded fields as Galactic GCs—even though, by the end of the mission, GC distances that are accurate to within the 1% level can be expected (Pancino et al. 2017). Confronting the Gaia DR2 data with distances from the literature, as independently compiled in the Harris (1996, 2010) catalog, a relatively small systematic offset, at the level of 0.029 mas, was found (Gaia Collaboration et al. 2018c), in the sense that parallaxes derived by Gaia are smaller than those implied by the distances given in Harris (2010). In any case, at this stage, the Gaia Collaboration et al. (2018c, 2018a) is using the latter distances, as opposed to those implied by the Gaia parallaxes, in its analyses of the Hertzsprung–Russell diagram and GC orbits.

Using the Harris (2010) distances, the Gaia Collaboration et al. (2018c) rederived the X, Y, Z coordinates of the GCs with respect to the Sun, given the improved positional information obtained by the Gaia mission. For our uniform set of distance measurements to the seven GCs studied here (Dist. 2), we used the distances calculated from the X, Y, Z coordinates in the Gaia Collaboration et al. (2018c). We note that these distances are in most cases consistent with those of Dist. 1, albeit with smaller uncertainties. The most significant difference between the two sets is for the GC ω Cen, although it has been noted that the dynamical measurement for this cluster (Watkins et al. 2013) may suffer from systematics. Finally, we note that Chen et al. (2018) reported a distance to 47 Tuc of $4.45\pm 0.01\pm 0.12\,\mathrm{kpc}$ (statistical and systematic uncertainties) obtained from a careful treatment of the Gaia DR2 parallaxes. This result is fully consistent with the values used in our sets Dist. 1 and Dist. 2. Using these two sets allows us to study the impact of the distance choices on the analyses of X-ray spectra of thermally emitting NSs.

The processing of the XMM-Newton and Chandra data sets is performed with the XMMSAS v15.0 (Gabriel et al. 2004) and CIAO v4.8 (Fruscione et al. 2006), respectively, following their respective standard procedures. The spectra are created from flare-filtered event files by extracting counts in circular regions. Background spectra are chosen from circular regions near the qLMXB, on the same CCD chip, and devoid of other sources. Finally, we grouped energy channels to ensure a minimum of 20 counts per bin. A detailed description of the data preparation is available in Guillot et al. (2013), and here we follow similar data reduction recipes.

The analysis of the qLMXB spectra is performed with PyXSPEC, the Python interface to the fitting package XSPEC (Arnaud 1996). This allows us to employ an MCMC approach to sample the parameter space, as described in Section 4. The spectral model used is the NS H-atmosphere model nsatmos (Heinke et al. 2006), modulated by absorption of soft X-rays by the interstellar medium (ISM). For the Galactic absorption, we used the recent model tbabs (Wilms et al. 2000). We also add a power-law component to account for a possible excess of counts above 2 keV that may originate from nonthermal emission. The exponent of this power law is fixed to 1.5, and we fit for the normalization. As will be shown below, the contribution of this power-law component is consistent with being null for all qLMXBs.

A pileup component is also added for all Chandra spectra, even those qLMXBs with low count rates inducing a pileup fraction ≲1%. As was pointed out by Bogdanov et al. (2016), uncorrected pileup, even at a low pileup fraction of ∼1%, can significantly bias the radius measurement. Specifically, for NGC 6397, the low ${R}_{\mathrm{NS}}$ obtained with H-atmosphere models was a consequence of the unmodeled pileup of photons in the X-ray spectra.

In summary, for each NS qLMXB in our sample, the spectral parameters of the model are

• 1.  
  the parameter α in the pileup model,
• 2.  
  the column density of neutral hydrogen N_H from the tbabs model,
• 3.  
  the NS surface temperature kT_eff in the nsatmos model,
• 4.  
  the NS mass in the nsatmos model,
• 5.  
  the NS radius in the nsatmos model,
• 6.  
  the NS distance (set as a prior; see Table 1) in the nsatmos model, and
• 7.  
  the power-law normalization (model powerlaw with fixed Γ = 1.5).

In addition, multiplicative constants are used to account for absolute flux cross-calibration uncertainties between different detectors (XMM-pn, XMM-MOS, and Chandra). Therefore, for sources with spectra obtained with multiple detectors, multiplicative constants are added to the spectral model, as is commonly done.⁸ In this work, all NSs are assumed to be described by the same EoS. Therefore, their masses and radii will be tied together by the parameterized EoS described in the following section.

Our EoS spans from the outer crust of NSs down to their dense core. We consider the HP94 model for the outer crust, which represents it as a Coulomb lattice of spherical nuclei immersed in a gas of electrons (Haensel & Pichon 1994). In this model, the nuclear masses are the experimental ones when available, supplemented by a theoretical mass formula (Möller & Nix 1992) for the more exotic nuclei. The inner crust starts when the energy density reaches $3.285\times {10}^{11}$ g cm⁻³, and we consider the tabulated SLY EoS (Douchin & Haensel 2001) obtained from a compressible liquid drop model based on the Skyrme interaction SLy4 (Chabanat et al. 1998). A test of the sensitivity on the crust EoS can be performed by replacing the SLY EoS with another one, such as the FPS one. These two tabulated EoSs can be downloaded.⁹

For numerical reasons, the transition between the crust and the core is guided within and log ρ–log P cubic spline matching the values and derivatives at both boundaries. The two boundaries are taken to be ${n}_{\mathrm{sat}}/10$ for the lower bound and ${n}_{\mathrm{sat}}$ for the upper bound. The sensitivity of this procedure to the choice of the boundaries is found to be small. Its impact on the total NS radius is less than 100 m, which is much smaller than current measurement uncertainties (Margueron et al. 2018a, 2018b).

In this work, we considered that the NS interior is made only of purely nucleonic matter whose properties are obtained from the extrapolation of the known saturation properties of nuclear matter. These properties are encoded in the so-called empirical parameters of nuclear matter, which are defined as being the coefficients of the series expansion in terms of the density parameter $x=(n-{n}_{\mathrm{sat}})/(3{n}_{\mathrm{sat}})$ of the energy per particle in symmetric matter,

$$\begin{eqnarray}&&{e}_{\mathrm{sat}}={E}_{\mathrm{sat}}+\displaystyle \frac{1}{2}{K}_{\mathrm{sat}}{x}^{2}+\displaystyle \frac{1}{3!}{Q}_{\mathrm{sat}}{x}^{3}+\displaystyle \frac{1}{4!}{Z}_{\mathrm{sat}}{x}^{4}+\,...,\end{eqnarray} \tag{ 2 }$$

and of the symmetry energy per particle,

$$\begin{eqnarray}\begin{array}{rcl}{e}_{\mathrm{sym}} & = & {E}_{\mathrm{sym}}+{L}_{\mathrm{sym}}x+\displaystyle \frac{1}{2}{K}_{\mathrm{sym}}{x}^{2}\\ & & +\ \displaystyle \frac{1}{3!}{Q}_{\mathrm{sym}}{x}^{3}+\displaystyle \frac{1}{4!}{Z}_{\mathrm{sym}}{x}^{4}+\,...,\end{array}\end{eqnarray} \tag{ 3 }$$

where the symmetry energy ${e}_{\mathrm{sym}}$ is defined as the isospin polarization energy

$$\begin{eqnarray}&&{e}_{\mathrm{sym}}=\displaystyle \frac{1}{2}\displaystyle \frac{{\partial }^{2}e}{\partial {\delta }^{2}},\end{eqnarray} \tag{ 4 }$$

$\delta =({n}_{n}-{n}_{p})/({n}_{n}+{n}_{p})$ is the isospin asymmetry parameter, and $e(n,\delta )$ is the nuclear energy per particle.

Here ${E}_{\mathrm{sat}}$ and ${E}_{\mathrm{sym}}$ are the saturation and symmetry energy at the saturation density ${n}_{\mathrm{sat}}$, and L_sym is the slope of the symmetry energy; since the saturation is an equilibrium point, there is no slope of the energy per particle in symmetric matter. Here ${K}_{\mathrm{sat}/\mathrm{sym}}$ stands for the curvature, ${Q}_{\mathrm{sat}/\mathrm{sym}}$ for the skewness, and ${Z}_{\mathrm{sat}/\mathrm{sym}}$ for the kurtosis of the energy per particle in symmetric matter and the symmetry energy, respectively. The values of these empirical parameters are determined from experimental measurements with different accuracies. Reviews of their experimental determination can be found in Margueron et al. (2018a, 2018b, references therein).

We consider the meta-modeling ELFc (for empirical local density functional) proposed in Margueron et al. (2018a, 2018b), which is based on the decomposition of the nuclear energy per particle in terms of a kinetic term t and a potential term v, as

$$\begin{eqnarray}&&e(n,\delta )=t(n,\delta )+v(n,\delta ).\end{eqnarray} \tag{ 5 }$$

The kinetic energy is defined as that of the Fermi gas plus medium corrections to the bare mass (encoded in the parameters ${\kappa }_{\mathrm{sat}/\mathrm{sym}}$),

$$\begin{eqnarray}&&t\left(n,\delta \right)=\displaystyle \frac{{t}_{\mathrm{sat}}}{2}{\left(\displaystyle \frac{n}{{n}_{\mathrm{sat}}}\right)}^{2/3}\left[\left(1+{\kappa }_{\mathrm{sat}}\displaystyle \frac{n}{{n}_{\mathrm{sat}}}\right){f}_{1}\left(\delta \right)+{\kappa }_{\mathrm{sym}}\displaystyle \frac{n}{{n}_{\mathrm{sat}}}{f}_{2}(\delta )\right],\end{eqnarray} \tag{ 6 }$$

where ${t}_{\mathrm{sat}}=3{{\hslash }}^{2}/(10m){\left(3{\pi }^{2}/2\right)}^{2/3}{n}_{\mathrm{sat}}^{2/3}$, m is the nucleon mass, and the functions ${f}_{1/2}$ are defined as

$$\begin{eqnarray}&&{f}_{1}\left(\delta \right)={\left(1+\delta \right)}^{5/3}+{\left(1-\delta \right)}^{5/3},\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{f}_{2}\left(\delta \right)=\left[{\left(1+\delta \right)}^{5/3}-{\left(1-\delta \right)}^{5/3}\right]\delta .\end{eqnarray} \tag{ 8 }$$

The potential term is expressed as

$$\begin{eqnarray}&&v\left(n,\delta \right)=\displaystyle \sum _{\alpha =0}^{N}\left({v}_{\alpha }^{\mathrm{is}}+{\delta }^{2}{v}_{\alpha }^{\mathrm{iv}}\right)\displaystyle \frac{{x}^{\alpha }}{\alpha !}u(x),\end{eqnarray} \tag{ 9 }$$

where the function u(x) takes into account the corrections due to the truncation N at low density, as

$$\begin{eqnarray}&&u\left(x\right)=1-{\left(-3x\right)}^{N+1-\alpha }\exp \left(-{bn}/{n}_{\mathrm{sat}}\right).\end{eqnarray} \tag{ 10 }$$

Fixing $b=10\mathrm{ln}2\approx 6.93$, as in Margueron et al. (2018a, 2018b), implies that the function u converges quickly to 1 as the density increases from zero. It ensures that $v(n,\delta )\to 0$ for $n\to 0$ for any order N. The larger the N, the smaller the correction u(x). The parameters ${v}_{\alpha }^{\mathrm{is}/\mathrm{iv}}$ entering into the series expansion of the potential term have a one-to-one relation with the empirical parameters.

The ability of this meta-modeling to reproduce existing EoSs increases as the order N increases. For N = 4, the meta-modeling can very accurately (at the few percent accuracy, in the worst case) reproduce the binding energy, pressure, and sound velocity of a large number of existing EoSs up to $4{n}_{\mathrm{sat}}$, as shown in Margueron et al. (2018a, 2018b).

In the present work, we use the flexibility of the meta-modeling to sample the parameter space of the empirical parameters using an MCMC approach. The range of variation for each of the empirical parameters considered in this analysis is given in Table 2. We fix the value of the lowest-order empirical parameters at saturation density to be ${E}_{\mathrm{sat}}\,=-15.8\,\mathrm{MeV}$, ${E}_{\mathrm{sym}}=32\,\mathrm{MeV}$, ${n}_{\mathrm{sat}}=0.155$ fm⁻³, and ${K}_{\mathrm{sat}}\,=230\,\mathrm{MeV}$. The parameters ${\kappa }_{\mathrm{sat}/\mathrm{sym}}$ are adjusted so that the Landau mass in symmetric matter is ${m}^{* }/m=0.75$ and the splitting between the neutron and proton Landau masses $({m}_{n}^{* }-{m}_{p}^{* })/m$ in neutron matter is 0.1 (see Table 2). The ${M}_{\mathrm{NS}}\mbox{--}{R}_{\mathrm{NS}}$ relation is known to be mostly influenced by the empirical parameters L_sym, K_sym, and Q_sat, since the EoS in the density range going from ${n}_{\mathrm{sat}}$ to approximately $3{n}_{\mathrm{sat}}$ depends most strongly on them (Margueron et al. 2018a, 2018b). Here L_sym and K_sym (respectively, Q_sat) control the density dependence of the symmetry energy (respectively, the energy per particle in symmetric matter) above saturation density. The higher-order empirical parameters are poorly known, but they impact the EoS at higher densities. They could, in principle, be deduced from ${M}_{\mathrm{NS}}$ and ${R}_{\mathrm{NS}}$ measurements for high-mass NSs.

Table 2.  Standard Values and Domain of Variation of the Empirical Parameters Considered in This Analysis; Taken from Margueron et al. (2018a, 2018b)

	${E}_{\mathrm{sat}}$	${E}_{\mathrm{sym}}$	${n}_{\mathrm{sat}}$	Lsym	${K}_{\mathrm{sat}}$	Ksym	Qsat	${Q}_{\mathrm{sym}}$	${Z}_{\mathrm{sat}}$	${Z}_{\mathrm{sym}}$	m*	Δm*
	(MeV)	(MeV)	(fm−3)	(MeV)	(MeV)	(MeV)	(MeV)	(MeV)	(MeV)	(MeV)	(mN)	(mN)
Standard	−15.8	32.0	0.155	60	230	−100	300	0	−500	−500	0.75	0.1
Variation	⋯	⋯	⋯	20–90	⋯	−400–200	−1300–1900	⋯	⋯	⋯	⋯	⋯

Note. See Section 3.1 for the description of the parameters.

Download table as:  ASCIITypeset image

A priori, we do not know which region of NS masses will be reached by our analysis. Anticipating our results, however, we find that the NS masses do not exceed 1.5–1.6 $\,{M}_{\odot }$, which implies that the central densities of these NSs are not very large, and the meta-model can reasonably be applied.

We illustrate here the impact of the empirical parameters L_sym, K_sym, and Q_sat on the ${M}_{\mathrm{NS}}\mbox{--}{R}_{\mathrm{NS}}$ relation. Since the rotation of the sources studied here is unknown, we consider nonrotating NS models, whose ${M}_{\mathrm{NS}}\mbox{--}{R}_{\mathrm{NS}}$ relation for a given EoS is obtained by solving the well-known Tolman–Oppenheimer–Volkoff (TOV) equations (Oppenheimer & Volkoff 1939; Tolman 1939). Only if the frequency is larger than 300 Hz (period <3 ms) could the rotational effects bias the ${R}_{\mathrm{NS}}$ measurements (Morsink et al. 2007). For an NS with a spin frequency of 600 Hz, its nonrotating radius would be underestimated by 2%–5%, depending on the NS size (Bauböck et al. 2013).

The EoS selection criteria include those that satisfy the requirements of causality and positiveness of the symmetry energy, as well as being compatible with a maximum mass above $1.9\,{M}_{\odot }$. This mass limit corresponds approximately to the 2σ lower limits of the measurements for PSR J1614–2230, $1.908\pm 0.016\,{M}_{\odot }$ (Demorest et al. 2010; Fonseca et al. 2016; Arzoumanian et al. 2018), and PSR J0348+0432 (Antoniadis et al. 2013), $2.01\pm 0.04\,{M}_{\odot }$.

Figure 1 shows the effect of varying the empirical parameters L_sym, K_sym, and Q_sat on the ${M}_{\mathrm{NS}}\mbox{--}{R}_{\mathrm{NS}}$ relation. Specifically, an increase of any of these three parameters shifts the high-mass part of the ${M}_{\mathrm{NS}}\mbox{--}{R}_{\mathrm{NS}}$ relation to larger radii. For clarity, only two parameters are varied in each of the top and bottom panels, the third parameter being kept fixed (${Q}_{\mathrm{sat}}=300\,\mathrm{MeV}$ in the top panel and ${K}_{\mathrm{sym}}=-85\,\mathrm{MeV}$ in the bottom panel). There are four groups of curves corresponding to the same value of L_sym and coinciding for very low mass NSs (${M}_{\mathrm{NS}}\lt 0.2\,{M}_{\odot }$). As ${M}_{\mathrm{NS}}$ increases, the central density increases as well, since we consider only the stable branch, and the different values for K_sym change the ${M}_{\mathrm{NS}}\mbox{--}{R}_{\mathrm{NS}}$ curves associated with the different EoSs. Overall, varying L_sym, K_sym, and Q_sat over the whole range allowed by nuclear physics, together with the requirement of supporting a $1.9\,{M}_{\odot }$ NS, yields radii between 11.5 and 14.2 km at 1.4 $\,{M}_{\odot }$. The effect of varying the parameter Q_sat is most noticeable for ${M}_{\mathrm{NS}}$ above 1.0–1.2 $\,{M}_{\odot }$. Being of higher order in the density expansion, Q_sat influences the EoS at a high density only or, equivalently, at a high ${M}_{\mathrm{NS}}$ only. Depending on the value of Q_sat, the EoS can be stiffer at high density, as reflected in the curves that go straight up, or softer at high density, letting the ${M}_{\mathrm{NS}}\mbox{--}{R}_{\mathrm{NS}}$ curve populate the low-${R}_{\mathrm{NS}}$ space at high ${M}_{\mathrm{NS}}$. There is, however, a limitation in the radii that can be explored based on the nucleonic EoS. As suggested in Margueron et al. (2018a, 2018b), low-mass NSs with ${R}_{\mathrm{NS}}\lt 11$ km (at $\sim 1.4\,{M}_{\odot }$) cannot be described by nucleonic EoS models that respect causality and must support a 1.9 $\,{M}_{\odot }$ NS.

Zoom In Zoom Out Reset image size

While there are various EoSs that pass through a point in the ${M}_{\mathrm{NS}}\mbox{--}{R}_{\mathrm{NS}}$ diagram, their paths are different. The degeneracy between different EoSs thus requires the knowledge of a set of ${M}_{\mathrm{NS}}\mbox{--}{R}_{\mathrm{NS}}$ points as distant as possible from each other. In conclusion of this analysis, the empirical parameters L_sym, K_sym, and Q_sat allow the exploration of a wide domain of ${M}_{\mathrm{NS}}$ and ${R}_{\mathrm{NS}}$ with various paths. Therefore, it may be possible to constrain the values of these parameters by confronting them with the observational data from the thermal X-ray emission of NSs.

For all of the cases considered here, the priors on the parameters are chosen so as to minimize any a priori assumption on the parameter distributions. All astrophysical parameters (except the distances to the sources) are sampled with uniform distributions within the boundaries allowed by the spectral model (defined in XSPEC).

The distances D are strongly coupled to the NS radii and effective surface temperatures. Letting this parameter explore a uniform prior would increase the uncertainties in the analysis enormously. For the two sets of distances presented in Table 1, we limit the qLMXB distances to Gaussian priors given by the central values and 1σ uncertainties listed. To do so, we add to the likelihood ${\chi }^{2}$ a penalty for each source i, proportional to the difference between the MCMC sampled distance ${D}_{\mathrm{mcmc},i}$ and the actual measured data ${D}_{\mathrm{data},i}$ (from Table 1), taking into account their standard deviations ${\sigma }_{i}$. The distance penalty reads ${\chi }_{D}^{2}={\sum }_{i=0}^{N}{\chi }_{D,i}^{2}$, where ${\chi }_{D,i}^{2}$ for each source is given by

$$\begin{eqnarray*}&&{\chi }_{D,i}^{2}=\displaystyle \frac{{\left({D}_{\mathrm{mcmc},i}-{D}_{\mathrm{data},i}\right)}^{2}}{{\left({\sigma }_{i}\right)}^{2}}.\end{eqnarray*}$$

The MCMC approach permits efficient sampling of our parameter space with high dimensionality: 49 parameters in total, including three nuclear physics EoS parameters, six astrophysical parameters per qLMXB (those listed in Section 2.3, except for the radii, which are obtained given the sampled EoS parameters and NS masses, after solving the TOV equations), and four multiplicative normalization constants (for the cross-calibration between the XMM-pn, XMM-MOS, and Chandra for the qLMXBs in M13 and ω Cen).

We use the python emcee package (Foreman-Mackey et al. 2013) with the stretch-move algorithm (Goodman & Weare 2010), which we applied as follows (see also the flowchart in Figure 2).

• 1.  
  Step 0. A large number of chains, or "walkers," are initialized, each one corresponding to a random point in the multidimensional parameter space defined by the set of parameters described above. We use 426 walkers (a multiple of the number of CPU cores available for our study).
• 2.  
  Step 1. We solve the TOV equations for each walker, providing 426 ${M}_{\mathrm{NS}}\mbox{--}{R}_{\mathrm{NS}}$ relations at each iteration.
• 3.  
  Step 2. For each walker, the sampled masses of the seven NSs are associated to seven calculated radii according to the ${M}_{\mathrm{NS}}\mbox{--}{R}_{\mathrm{NS}}$ relation. Using those ${M}_{\mathrm{NS}}$ and ${R}_{\mathrm{NS}}$ and the other sampled astrophysical parameters, we calculate the global ${\chi }^{2}$ between the emission models (NS atmosphere) and the data for the seven NSs.
• 4.  
  Step 3. Given the calculated probability (its likelihood multiplied by the distance Gaussian priors mentioned above), the evolution of the walkers in the parameter space is decided according to the stretch-move algorithm. To determine the new position, each walker is randomly paired with another and will move along the line joining the two current points in the parameter space. The amount of "stretch" is determined by the scale parameter (the only adjustable parameter in this algorithm), which has been chosen to be a = 2.0 as prescribed in Goodman & Weare (2010). The new position is accepted or rejected depending on its probability. For more details about the stretch-move algorithm, see Goodman & Weare (2010) and Foreman-Mackey et al. (2013).
• 5.  
  Step 4. Steps 1–3 are repeated numerous times until the walkers have converged in the region of the parameter space resulting in the highest likelihood or minimum ${\chi }^{2}$ value.
• 6.  
  Step 5. When the MCMC loop stops, the statistical posterior distributions are calculated and marginalized to create the outputs.

Zoom In Zoom Out Reset image size

Before running the code on the full data set and with the most general meta-modeling in Section 5, we first test it considering the constant-${R}_{\mathrm{NS}}$ toy model. In addition to its simplicity, this test is interesting because it allows us to compare with results that have already been reported in the literature.

We first consider the constant-${R}_{\mathrm{NS}}$ model (Guillot et al. 2013), which assumes that all NSs have the same radius, i.e., that the EoS is represented in ${M}_{\mathrm{NS}}\mbox{--}{R}_{\mathrm{NS}}$ space by a vertical line in which ${R}_{\mathrm{NS}}$ is independent of ${M}_{\mathrm{NS}}$ (which remains as a free parameter). This is a simple toy-model approximation motivated by the observations that most nucleonic EoSs (the ones consistent with 2 $\,{M}_{\odot }$) have a rather weak dependence on ${M}_{\mathrm{NS}}$, between 1 and 2 $\,{M}_{\odot }$. The purpose of this toy model is mainly to test our code and MCMC approach.

After running the analysis described in Section 4.1, we obtain the results shown in Figure 3, i.e., the ${R}_{\mathrm{NS}}$ posterior distributions (the same for all seven NSs) considering the distance sets Dist. 1 and Dist. 2, marginalized over the other parameters of the model. For both distance sets, the radius distributions are ${R}_{\mathrm{NS}}={11.09}_{-0.36}^{+0.38}\,\mathrm{km}$ (Dist. 1, for ${\chi }_{\nu }^{2}=1.06$) and ${R}_{\mathrm{NS}}={11.04}_{-0.35}^{+0.39}\,\mathrm{km}$ (Dist. 2, for ${\chi }_{\nu }^{2}=1.07$). These values are consistent with the recent results of Guillot (2016) but at odds with older results (e.g., Guillot et al. 2013; Guillot & Rutledge 2014). The differences are likely due to the inclusion of new sources (47 Tuc X-7) and data, the use of recent distance measurements, the improvement of the analysis (e.g., the new absorption model tbabs), and the inclusion of the pileup correction model for all sources (including those with an $\sim 1 \%$ pileup fraction; see Section 2.1). Overall, we find a radius distribution that is easier to reconcile with the nuclear physics models of Section 3, i.e., having nonnegligible probabilities for an NS radius larger than about 11 km.

Zoom In Zoom Out Reset image size

There is, however, also a large fraction of the posterior probability distribution that is located below 11 km, in conflict with nuclear physics expectations (e.g., Margueron et al. 2018a, 2018b; Tews et al. 2018b), as well as our illustrative Figure 1. For instance, one could deduce from these figures that ${R}_{\mathrm{NS}}\lesssim 11\,\mathrm{km}$ requires ${L}_{\mathrm{sym}}\lesssim 20\,\mathrm{MeV}$, which contradicts nuclear physics expectations (Lattimer & Lim 2013).

In the following section, we address the question of the compatibility between the thermal emission modeling and the nuclear EoS by including the meta-model directly in the global spectral data analysis. In this way, we show that there is no inconsistency between the observational data and the nuclear EoS, and we extract an estimation for the nuclear EoS parameters.

We remind the reader that the main features of our work are that we (i) simultaneously fit seven NS qLMXB sources, (ii) impose the same EoS on all of these sources, and (iii) treat the EoS and the astrophysical model parameters equally.

Only a few nuclear EoS parameters are taken as free parameters. We recall that the nuclear meta-modeling is governed by a set of empirical parameters (see Section 3). Some of these empirical parameters can be well constrained by nuclear experiments (see the discussion in Margueron et al. 2018a, 2018b), and they are kept fixed in the present analysis: ${n}_{\mathrm{sat}}$, ${E}_{\mathrm{sat}}$, ${E}_{\mathrm{sym}}$, and ${K}_{\mathrm{sat}}$ (values in Table 2). The more influential and less known parameters, L_sym, K_sym, and Q_sat, are fitted in our analysis. The values of K_sym and Q_sat are currently unknown, while there exist constraints on L_sym from nuclear experiments and theoretical predictions. These constraints indicate that L_sym has a value around 50 MeV with an uncertainty of about ±10 MeV (Lattimer & Lim 2013). We therefore incorporate this knowledge from nuclear physics by considering a Gaussian prior on L_sym centered at 50 MeV with a width of 10 MeV. Since K_sym and Q_sat are unknown, we consider a uniform distribution in the wide ranges listed in Table 2. The higher-order empirical parameters, ${Q}_{\mathrm{sym}}$ and ${Z}_{\mathrm{sat}/\mathrm{sym}}$, are not known. However, since they influence only the high-density part of the EoS, they will not be tightly constrained by the present analysis. Therefore, they can be fixed to the values listed in Table 2 (Margueron et al. 2018a, 2018b).

The MCMC routine (described in Section 4) was run considering the seven qLMXB sources mentioned above. We have considered the chains that converged to the global minimum, excluding a few percent (1%–5%) stuck at higher ${\chi }^{2}$ values (typically for reduced ${\chi }^{2}$ above 10). We tested the presence of these "stuck chains" with repeated iterations of the exact same MCMC run. In each case, the minimum best-fit ${\chi }^{2}$ is always found to be the same, and a small fraction of chains remain in the high-${\chi }^{2}$ parts of the parameter space. After 150,000 iterations, the reduced ${\chi }^{2}$ distribution is centered around 1.10 ± 0.02 for 1126 d.o.f., and the best fit corresponds to ${\chi }^{2}=1.08$, giving a null hypothesis probability of 3.1%.

The marginalized posterior probabilities for the empirical EoS parameters L_sym, K_sym, and Q_sat are shown in Figure 4. We observed from the marginalized distributions that L_sym peaks at lower values than the one imposed by the prior ($50\pm 10\,\mathrm{MeV}$) but remains consistent with it: ${L}_{\mathrm{sym}}={37.2}_{-8.9}^{+9.2}\,\mathrm{MeV}$. This somewhat reflects the tension driving the fit toward low radii at low masses (see Figure 1 and related discussion). The empirical parameter ${K}_{\mathrm{sym}}=-{85}_{-70}^{+82}\,\mathrm{MeV}$ is rather well constrained compared to the uniform prior, showing that this parameter is important for our data set. Notice that it is also remarkably compatible with the one $-100\pm 100\,\mathrm{MeV}$ extracted from analysis of the chiral effective field theory (EFT) calculations (Margueron et al. 2018a, 2018b). Finally, the empirical parameter ${Q}_{\mathrm{sat}}={318}_{-366}^{+673}\,\mathrm{MeV}$ is less constrained, but there is a preference for the lower values of the uniform prior distribution. The values of the empirical parameters for this run are reported in the first row of Table 3. We point out that, despite the rather large uncertainties on the empirical parameters K_sym and Q_sat, this is the first time that these parameters are extracted from data.

Zoom In Zoom Out Reset image size

Table 3.  Distribution of Empirical Parameters ${L}_{\mathrm{sym}}$, ${K}_{\mathrm{sym}}$, and ${Q}_{\mathrm{sat}}$ for Various Cases

Framework	Sources	Distances	Prior	Lsym	Ksym	Qsat	${R}_{1.45}$	${\chi }_{\nu }^{2}$	No. of	d.o.f.
			Lsym	(MeV)	(MeV)	(MeV)	(km)		Param.
1	All	Dist. 2	Yes	${37.2}_{-8.9}^{+9.2}$	−85${}_{-70}^{+82}$	${318}_{-366}^{+673}$	12.35 ± 0.37	1.08	49	1126
2	All	Dist. 1	Yes	${38.3}_{-8.9}^{+9.1}$	−91${}_{-71}^{+85}$	${353}_{-484}^{+696}$	12.42 ± 0.34	1.07	49	1126
3	All	Dist. 1	Yes	${38.6}_{-8.7}^{+9.2}$	−95${}_{-36}^{+80}$	300	12.25 ± 0.30	1.07	48	1127
4	All	Dist. 1	No	${27.2}_{-5.3}^{+10.9}$	−59${}_{-74}^{+103}$	${408}_{-430}^{+735}$	12.37 ± 0.30	1.07	49	1126
5	All/47 Tuc	Dist. 1	Yes	${43.4}_{-9.3}^{+9.7}$	−66${}_{-102}^{+137}$	${622}_{-560}^{+763}$	12.57 ± 0.41	1.08	43	700
6	All/NGC 6397	Dist. 1	Yes	${42.6}_{-9.5}^{+9.9}$	−77${}_{-96}^{+129}$	${623}_{-544}^{+757}$	12.58 ± 0.40	1.09	43	961
7	All/M28	Dist. 1	Yes	${42.5}_{-9.5}^{+9.5}$	−80${}_{-91}^{+124}$	${597}_{-510}^{+717}$	12.46 ± 0.37	1.07	43	846
8	A	Dist. 2	Yes	${38.6}_{-8.9}^{+9.4}$	−91${}_{-76}^{+81}$	${343}_{-431}^{+805}$	12.18 ± 0.29	1.04	21	874
9	${A}^{{\prime} }$	Dist. 2	Yes	${37.5}_{-8.9}^{+9.0}$	−88${}_{-70}^{+76}$	${263}_{-361}^{+764}$	12.22 ± 0.32	1.06	29	945
10	B	Dist. 2	Yes	${49.12}_{-10.0}^{+10.0}$	−6.66${}_{-138}^{+137}$	${804}_{-675}^{+709}$	12.88 ± 0.43	1.19	28	255
11	${B}^{{\prime} }$	Dist. 2	Yes	${50.3}_{-9.6}^{+9.8}$	−1${}_{-143}^{+134}$	${881}_{-705}^{+671}$	12.98 ± 0.40	1.18	23	178

Note. Group A contains high-S/N sources (peaked masses): NGC 6397, 47 Tuc, and M28. Group B contains low-S/N sources (flat masses): ω Cen, NGC 6304, M13, and M30. Group A' contains sources with peaked masses: NGC 6397, 47 Tuc, M28, and M13. Group B' contains sources with almost-flat masses: ω Cen, NGC 6304, and M30. See text for more details.

Download table as:  ASCIITypeset image

The correlations among empirical parameters are also visible in Figure 4. There is a weak anticorrelation between L_sym and K_sym and a stronger anticorrelation between K_sym and Q_sat. These correlations reflect the causality and stability requirements, implying, for instance, that a large value for K_sym shall be compensated by a small value of L_sym or Q_sat to limit the upper bound for the sound velocity, and vice versa for the lower bound. The anticorrelation between L_sym and K_sym was already found in Margueron et al. (2018a, 2018b), but the empirical parameters L_sym/K_sym and Q_sat were found to be correlated for a stiff EoS (if the direct Urca process occurs for ${M}_{\mathrm{NS}}\lt 2\,{M}_{\odot }$), while for a soft EoS (no direct URCA possible for ${M}_{\mathrm{NS}}\lt 2\,{M}_{\odot }$), no correlations were found. The anticorrelation between K_sym and Q_sat is therefore a new feature coming from the fit to the thermal X-ray emission.

The ${M}_{\mathrm{NS}}\mbox{--}{R}_{\mathrm{NS}}$ posterior probability distributions corresponding to the MCMC runs with the distances of set Dist. 1 (top panel) and set Dist. 2 (bottom panel) are displayed in Figure 5. It is reassuring to notice that the ${M}_{\mathrm{NS}}\mbox{--}{R}_{\mathrm{NS}}$ posterior probability distribution is almost insensitive to the set of distances considered, as was also observed for the constant-${R}_{\mathrm{NS}}$ test runs (Figure 3). The global features of the probability distribution are the same for the two distance sets: the radius that we obtain is between ∼11.5 and $13.0\,\mathrm{km}$ for a 1.4 $\,{M}_{\odot }$ NS, with a preference for low masses, although the 90% credible intervals are compatible with 2 $\,{M}_{\odot }$.

Zoom In Zoom Out Reset image size

In Figure 6 (top panels), we present the 90% credible interval ${M}_{\mathrm{NS}}\mbox{--}{R}_{\mathrm{NS}}$ posterior distributions of individual sources, obtained with the distance sets Dist. 2 (panel (a)) and Dist. 1 (panel (b)). Most sources have credible intervals that reach $\sim 1.9\,{M}_{\odot }$ or higher. In a few cases, the 90% credible intervals only reach masses around 1.4–1.5 $\,{M}_{\odot }$, which appear to favor low-mass NSs. However, this is compatible with current distributions of MSP masses (Antoniadis et al. 2016; Özel & Freire 2016) that descend from LMXBs; we note that the lowest known NS mass is 1.174 ± 0.004 $\,{M}_{\odot }$ (for PSR J0453+1559; Martinez et al. 2015). Therefore, at the moment, there are no discrepancies with our current knowledge of NS formation mechanisms and their expected masses.

Zoom In Zoom Out Reset image size

We note that the masses of the individual sources are not constrained as well as the radius (Figure 6). This is inherent to the method used in the present work, in which the measurable physical quantity is ${R}_{\infty }$. The constraints on the radius emerge from the combination of (1) the general shape in ${M}_{\mathrm{NS}}\mbox{--}{R}_{\mathrm{NS}}$ space of most EOS models in our nucleonic parameterization, and (2) the shape of the quantity ${R}_{\infty }$ extracted from qLMXB spectra (see previous works, e.g., Guillot et al. 2011; Bogdanov et al. 2016; Shaw et al. 2018 for the ${M}_{\mathrm{NS}}\mbox{--}{R}_{\mathrm{NS}}$ constraints from single qLMXBs, for which a significant portion of the ${M}_{\mathrm{NS}}\mbox{--}{R}_{\mathrm{NS}}$ contours appear at constant ${R}_{\mathrm{NS}}$). Therefore, unless the mass of the NS is measured independently, the degeneracy between ${M}_{\mathrm{NS}}$ and ${R}_{\mathrm{NS}}$ from observations of qLMXB can only be minimized by the implementation of an EoS parameterization as done in this work. Other events involving NSs enable measurements of ${M}_{\mathrm{NS}}$ and ${R}_{\mathrm{NS}}$ independently. For example, type I X-ray bursts with photospheric radius expansion bursts provide two observables (the Eddington flux and the cooling tail normalization, which both depend on ${M}_{\mathrm{NS}}$ and ${R}_{\mathrm{NS}}$; see Özel & Freire 2016 for a review). The GW signals of two merging NSs, on the other hand, provide measurements of the merging masses and tidal deformability, which can be used to derive the radius (e.g., Abbott et al. 2018; De et al. 2018).

The 50%, 90%, and 99% contours resulting from the present analysis can also be compared with the constraints from AT 2017gfo (Bauswein et al. 2017) and GW170817 (Abbott et al. 2018; Annala et al. 2018; Tews et al. 2018b). There is a good agreement between these different constraints. Nonetheless, the width of the distribution obtained from the present work appears narrower than those from analyses of GW170817, indicative of more restrictive constraints on the NS radius, but this could also be due to the fact that we do not consider phase transitions in the meta-modeling of our work.

The most likely EoS properties can also be deduced from our MCMC analyses, since the EoSs are described by empirical parameters. In Figure 7, we show the boundaries of the relation between the total pressure P and the energy density ρ resulting from our analysis, considering the nucleon and lepton contributions in β-equilibrium and for the two distance sets Dist. 1 and Dist. 2. As noted before, the two distance sets do not significantly affect the most likely EoSs defined by those boundaries. Our predictions for the EoS are contrasted with a prediction based on different constraints and labeled as MM in Figure 7. The meta-model MM is constrained by quantum Monte Carlo predictions in low-density nuclear matter up to n_sat and based on two- and three-nucleon forces from the chiral EFT Hamiltonians given in Tews et al. (2018a). The extrapolation beyond n_sat is controlled by causality and stability requirements, as well as positiveness of the symmetry energy and maximum observed NS masses predicted in Tews et al. (2018b). There is a good overlap between the EoS deduced from our analysis and the one from the MM analysis. It is, however, interesting to note that the intersection between the bands generated from our analysis and the MM prediction could potentially further reduce the possibilities for the EoS.

Zoom In Zoom Out Reset image size

In Figure 8, we show both the squared sound speed in units of the speed of light, ${\left({v}_{s}/c\right)}^{2}$, and the symmetry energy, ${e}_{\mathrm{sym}}$, as functions of the energy density ρ. The colored bands in Figure 8 are the same as in Figure 7. Here too, one can remark that there is little impact of the distance sets on ${\left({v}_{s}/c\right)}^{2}$ and ${e}_{\mathrm{sym}}$, as for the general EoS shown in Figure 7. Similar to Figure 7, there is also a good overlap of our predictions with the MM ones.

Zoom In Zoom Out Reset image size

The posterior ranges at 98% confidence for the qLMXB emission model parameters are given in Table 4 for the two distances considered here, Dist. 1 and Dist. 2. First, we note that all parameters resulting from the Dist. 1 run are consistent with those of Dist. 2. The small differences observed between the results of these two runs are not significant—only the seven distance posterior distributions differ, since they are driven by the priors imposed. The NS temperatures and masses are consistent with previously reported values (Guillot et al. 2013; Guillot & Rutledge 2014; Heinke et al. 2014; Bogdanov et al. 2016). Interestingly, none of the NSs studied have masses going over $\sim 2.1\,{M}_{\odot }$ at 98% confidence. The best-fit absorption values N_H are also consistent with the expected values in the direction of the host GCs (see, e.g., neutral H maps; Dickey & Lockman 1990; Kalberla et al. 2005). Finally, we note that the obtained power-law normalizations ${N}_{\mathrm{pl}}$ are consistent with zero. Although this might not be readily obvious from the quantile ranges provided in Table 4, the seven posterior distributions (not shown in the paper) do indeed have nonzero probabilities for ${N}_{\mathrm{pl}}=0$. This lends further evidence for the absence of nonthermal emission in these objects. We have nonetheless considered the possible existence of nonthermal emission in our analyses by including a power-law component in the spectral model.

Table 4.  Distributions of All Model Parameters with Quantiles Corresponding to the 98% Credible Interval, Except the Empirical Parameters Given in Table 3 for the Reference Calculation (Distances of Dist. 2, Prior on L_sym, and Variation over the Three Empirical Parameters ${L}_{\mathrm{sym}}$, ${K}_{\mathrm{sym}}$, and ${Q}_{\mathrm{sat}}$)

Parameter	Source	Dist. 1	Dist. 2
	M13	0.15–0.96	0.15–0.95
	ω Cen	0.23–0.97	0.20–0.97
	47 Tuc	0.15–0.65	0.16–0.67
Pileup α	M28	0.38–0.54	0.39–0.55
	M30	0.27–0.97	0.24–0.97
	NGC 6304	0.20–0.97	0.20–0.97
	NGC 6397	0.33–0.81	0.38–0.85
	M13	1.56–4.44	2.47–5.43
	ω Cen	17.19–21.01	15.76–19.32
	47 Tuc	2.86–4.64	2.93–4.76
${N}_{{\rm{H}}}\ \left({10}^{20}\,{\mathrm{cm}}^{-2}\right)$	M28	35.17–37.98	35.34–37.88
	M30	3.01–6.37	3.18–6.47
	NGC 6304	44.96–56.61	45.98–58.13
	NGC 6397	16.60–17.73	17.38–19.44
	M13	79.06–90.42	76.04–86.57
	ω Cen	70.48–80.96	73.62–86.05
	47 Tuc	104.73–112.71	104.54–112.79
${{kT}}_{\mathrm{eff}}\ \left(\mathrm{eV}\right)$	M28	109.01–118.76	108.49–116.35
	M30	86.69–99.99	86.51–99.17
	NGC 6304	91.93–108.39	90.26–105.36
	NGC 6397	62.87–68.88	61.16–65.88
	M13	0.77–1.95	0.74–1.95
	ω Cen	0.80–2.01	0.85–2.06
	47 Tuc	0.66–1.44	0.67–1.47
${M}_{\mathrm{NS}}\ ({M}_{\odot })$	M28	0.70–1.51	0.68–1.43
	M30	0.78–2.00	0.80–1.99
	NGC 6304	0.87–2.07	0.85–2.06
	NGC 6397	0.72–1.62	0.70–1.47
	M13	7.72–8.50	7.07–7.19
	ω Cen	4.50–4.65	5.18–5.25
	47 Tuc	4.50–4.65	4.47–4.59
$D\ \left(\mathrm{kpc}\right)$	M28	5.48–5.89	5.467–5.65
	M30	8.02–8.77	8.06–8.21
	NGC 6304	6.15–6.43	5.86–6.01
	NGC 6397	2.49–2.56	2.29–2.34
	M13	4.55–15.57	5.64–16.48
	ω Cen	1.21–5.20	0.89–4.87
	47 Tuc	1.15–11.50	1.17–12.44
${N}_{\mathrm{pl}}$	M28	2.30–10.63	2.23–10.32
	M30	4.01–19.99	4.16–19.54
	NGC 6304	3.87–15.78	4.42–15.66
	NGC 6397	2.83–8.19	2.89–8.30
C1	M13	0.89–1.09	0.89–1.09
	ω Cen	0.99–1.16	0.99–1.18
C2	M13	0.79–0.96	0.78–0.99
	ω Cen	0.96–1.12	0.96–1.12

Note. Here C₁ and C₂ are the multiplicative coefficients that account for absolute flux cross-calibration uncertainties between the XMM-pn, XMM-MOS, and Chandra detectors (see Section 2.3).

Download table as:  ASCIITypeset image

This section presents a sensitivity analysis of our results in which modifications of the main framework are tested, such as reducing the number of empirical parameters to vary, changing the set of distances considered (notice that this was already largely explored in the previous subsection), or reducing the number of qLMXB sources considered.

We report in Table 3 the global results of the sensitivity analysis, where the impact of the changes is given for a few parameters: the EoS empirical parameters L_sym, K_sym, and Q_sat; the radius ${R}_{1.45}$ for a $1.45\,{M}_{\odot }$ NS; and the best ${\chi }_{\nu }^{2}$. We also give the number of fitting parameters and d.o.f. for each run. The rows of Table 3 represent the various frameworks considered. The first two rows represent the framework already explored, considering all seven qLMXB sources, the two sets of distances Dist. 1 and Dist. 2, the prior on L_sym, and the variation of the three empirical parameters. They are considered hereafter as our reference results, around which we slightly perturb the framework to extract the sensitivity of this reference to small corrections.

In the first approach for the sensitivity analysis, we modify the distance sets. For framework 3, we reduce the number of free EoS parameters by fixing ${Q}_{\mathrm{sat}}=300\,\mathrm{MeV}$. The impact on the centroid and width for L_sym, K_sym, and ${R}_{1.45}$ is marginal. The minimum ${\chi }_{\nu }^{2}$ changes minimally. For framework 4, we replace the Gaussian prior on L_sym with a uniform prior ranging from 20 to 120 MeV. The change in the minimum ${\chi }_{\nu }^{2}$ is marginal, indicating that the fit statistic is not affected by adding/removing the prior on L_sym. Furthermore, removing the prior on L_sym somewhat reduces the posterior value, which should have the effect of decreasing the overall radius (as indicated by Figure 5). However, we observe that ${R}_{1.45}$ remains unchanged (marginal decrease; see Table 3) because lower values of L_sym are compensated for by larger values of K_sym and Q_sat. This emerges from the anticorrelation between L_sym and K_sym, as discussed in Section 5.2 and reported in Margueron et al. (2018a, 2018b). This observed compensation between the empirical parameters originates from the fact that in our meta-model, ${e}_{\mathrm{sym}}$ is constrained to positive values up to a threshold central density that produces 1.9 $\,{M}_{\odot }$ NSs. Therefore, if L_sym becomes too low, the other parameters (mostly K_sym) readjust to satisfy the condition on ${e}_{\mathrm{sym}}$. In the spectral analyses of this work, this translates into a stable average radius (Table 3), as required by the observational data.

In a second approach, we analyze the sensitivity of the result to the modification of the qLMXB source set by removing a single qLMXB. In frameworks 4, 6, and 7, we removed 47 Tuc X-7, NGC 6397, and M28, respectively. Since these sources have the largest signal-to-noise ratio (S/N), their removal allows us to check to what extent they contribute to drive the results. While marginal, there are indeed some systematic effects: L_sym is increased by about 5–6 MeV, K_sym is increased by 10–20 MeV, the value for Q_sat is almost doubled, and the radius ${R}_{1.45}$ is increased by up to 0.23 km. These systematic corrections remain inside the original uncertainty estimated for the reference results (frameworks 1 and 2). These different approaches are summarized in Figure 6 (bottom), which shows the 90% credible interval of the ${M}_{\mathrm{NS}}\mbox{--}{R}_{\mathrm{NS}}$ posterior distributions of frameworks 1–7 and demonstrates that they broadly overlap in the 12–13 km radius range.

In a third approach for the sensitivity analysis, we split the qLMXB sources into the different groups (A, B) and (${A}^{{\prime} }$, $B^{\prime}$). Groups A and B are defined with respect to the S/N (A for S/N > 60, B otherwise; see Table 1). Groups $A^{\prime}$ and ${B}^{{\prime} }$ are defined with respect to the posterior mass distribution (${A}^{{\prime} }$ if the posterior mass distribution is well peaked, ${B}^{{\prime} }$ if it is almost flat). There is a nice correlation between the S/N ratio and the posterior mass distribution ($A=A^{\prime}$ and $B=B^{\prime}$), except for the source M13, which has a low S/N but a well-peaked mass distribution (see Table 1). As a consequence, the results for groups A and $A^{\prime}$, as well as B and $B^{\prime}$, are almost identical. Groups A and $A^{\prime}$ prefer lower values for L_sym, K_sym, and Q_sat comparable to the reference results. They favor lower radii ${R}_{1.45}\approx 12.2\pm 0.3$ km. By contrast, groups B and $B^{\prime}$ tend to increase the values for L_sym, K_sym, Q_sat, and ${R}_{1.45}$ to values that are still compatible with the uncertainty of the reference results, albeit with some tension. Naturally, the uncertainty on these values is also increased, especially for the parameter K_sym and radius ${R}_{1.45}$. We also note that, for groups B and $B^{\prime}$, the L_sym values are essentially identical to the prior given on that parameter (${L}_{\mathrm{sym}}=50\pm 10\,\mathrm{MeV}$), implying that these two groups have little weight in the constraints on L_sym.

As a conclusion of this sensitivity analysis, we can state that our reference results are only marginally impacted by small changes in the crucial input parameters, such as the distance set, the number of free EoS parameters, and the selection of qLMXB sources. In addition, we identified a group of qLMXB sources with low S/N (subsets B and $B^{\prime}$), which do not contribute significantly to the constraints on the empirical parameters, especially L_sym. These are the qLMXBs in ω Cen, M13, M30, and NGC 6304. An improvement in the analysis of the qLMXB thermal emission will require more statistics especially for these sources.

Since the seminal papers of Brown et al. (1998) and Rutledge et al. (2002), the thermal emission from qLMXBs has been analyzed by several authors in order to better constrain the properties of matter at high density. Over the years, atmosphere models have been improved (e.g., Heinke et al. 2006; Haakonsen et al. 2012) and the number of sources used in the analysis has increased (Guillot & Rutledge 2014; Bogdanov et al. 2016). The theoretical description of the EoS has also been improved, from the unconstrained case where masses and radii are considered independently of each other (i.e., directly extracted from ${R}_{\infty }$ measurements; e.g., Heinke et al. 2006; Guillot et al. 2011) to more consistent approaches. In a first attempt to consistently analyze several qLMXB sources combined, a constant-radius EoS model was proposed, inspired by the qualitative behavior of most nuclear EoSs (Guillot et al. 2013; Guillot 2016).

Because these early results did not consider a full treatment of the pileup instrumental effects in the Chandra data (which are significant even at low pileup fractions; Bogdanov et al. 2016), we only compared our results to the most recent ones in which qLMXBs are analyzed including the effects of the pileup and which contain similar inputs as in our analysis. Recently, Steiner et al. (2018) found that the radius of a $1.4\,{M}_{\odot }$ NS is most likely between 10.4 and 13.7 km at the 68% confidence level, considering all cases tested in that work. Assuming a pure H atmosphere for all objects, they found ${R}_{\mathrm{NS}}$ in the range 11.2–12.3 km, which is consistent with our results. In comparison, the interval of possible radii in the present work is narrower, since we disregarded the possible occurrence of a strong phase transition. For L_sym, Steiner et al. (2018) found 38.94–58.09 MeV, which is also consistent with our findings, while the uncertainty band is also larger in their case.

However, the main difference between our analysis and that of Steiner et al. (2018) is that we have implemented the EoS parameters in the fitting procedure, while Steiner et al. (2018) determined an ${M}_{\mathrm{NS}}\mbox{--}{R}_{\mathrm{NS}}$ posterior probability independent from the EoS and, in a second step, fit different EoS scenarios to this posterior result. It is reassuring to find that our results agree.

In another analysis, Özel et al. (2016) analyzed the thermal emission of the same sources as ours, except 47 Tuc X-7, in addition to data from six type I X-ray bursts. They found radii between 10.1 and 11.1 km for masses ranging from 1 to 2 $\,{M}_{\odot }$, which is a smaller estimation than ours. In a more recent analysis, Bogdanov et al. (2016) included the same 12 sources as Özel & Freire (2016), with the addition of 47 Tuc X-5 and X-7, and found radii ranging from 9.9 to 11.2 km. These two analyses favor a rather soft EoS, at odds with our results. There are still some differences between these analyses and ours, including (1) different values of the distances; (2) they included the X-ray burst data, which we did not; and (3) they used polytropes to parameterize the EoS. Another main difference from our work is that they deduced the radii of NSs from the marginalized posterior mass distributions (as in Steiner et al. 2018), while in our case, the radii are calculated consistently with the masses for each considered EoS. In our analysis, we have shown that without nuclear physics inputs, the constant-${R}_{\mathrm{NS}}$ approximation prefers radii around $\sim 11.1\pm 0.4\,\mathrm{km}$, consistent with the estimates in Özel et al. (2016) and Steiner et al. (2018), while including a nuclear EoS and a prior on the empirical parameter L_sym can increase the radius up to ∼12.0–12.5 km. This demonstrates the advantage of fitting the thermal emission model parameters together with the ones of the EoS.

Recently, Nättilä et al. (2017) performed the first direct atmosphere model spectral analysis of five hard-state type I X-ray burst cooling tails from the LMXB 4U 1702–429. They extracted a precise estimation of the radius, $12.4\pm 0.4\,\mathrm{km}$ at 68% credibility, for a mass more difficult to constrain, in the range 1.4–2.2 $\,{M}_{\odot }$.

Observations of MSPs also provided measurements of the NS radius and therefore constraints on the EoS. While the early analyses provided lower limits on the radius (e.g., $R\gt 10.4\,\mathrm{km}$ for PSR J0437−4715; Bogdanov 2013), the recent NS parameter estimation resulting from X-ray pulse profile analyses of NICER data resulted in better constrained radii (with ∼10% uncertainties; Miller et al. 2019; Riley et al. 2019), compatible with those reported here. A different analysis, which exploited the far-ultraviolet and soft X-ray emission of PSR J0437−4715 fitted to a low-temperature atmosphere model, resulted in ${R}_{{\rm{N}}{\rm{S}}}\,={13.6}_{-0.8}^{+0.9}\,{\rm{k}}{\rm{m}}$, compatible with our results here, although with some moderate tension (Gonzalez-Caniulef et al. 2019).

Finally, the recent observation of the NS–NS merger event GW170817 resulted in an estimation of the radii of the two stars, as well as constraints on their EoS through the tidal deformability parameter Λ (Abbott et al. 2018). Further analyses of the GW and electromagnetic signals lead to the constraints on the radii shown in Figure 5. A good agreement with our analysis should also be noticed.