New constraints on the nuclear equation of state from the thermal emission of neutron stars in quiescent low-mass X-ray binaries

TL;DR

This study uses thermal emission data from neutron stars in quiescent low-mass X-ray binaries to constrain the nuclear equation of state, providing new empirical parameter measurements and reconciling previous radius estimates with nuclear physics.

Contribution

It introduces a novel analysis combining MCMC with empirical EOS parametrization, measuring key nuclear matter parameters for the first time from astrophysical data.

Findings

Neutron star radii are around 12 km, consistent with nuclear physics.

Measured symmetry energy slope $L_{sym}$ as approximately 37 MeV.

First empirical constraints on $K_{sym}$ and $Q_{sat}$ from neutron star observations.

Abstract

This paper presents a new analysis of the thermal emission from the neutron star surface to constrain the dense matter equation of state. It is based on the use of a Markov-Chain Monte Carlo algorithm combined with an empirical parametrization of the equation of state, as well as the consistent treatment of seven neutron star quiescent low-mass X-ray binaries in globular clusters with well-measured distances. Previous analyses have indicated that the thermal emission of these neutron stars tends to prefer low neutron star radii, questioning basic knowledge from nuclear physics. We show that it is possible to reconcile the thermal emission analyses with nuclear physics knowledge, with or without including a prior on the slope of the symmetry energy $L_{\rm sym}$. We obtain radii of the order of about 12~km without worsening the fit statistic. With an empirical parametrization of the equation of state, we obtain the following values for the slope of the symmetry energy, its curvature $K_{\rm sym}$, and the isoscalar skewness parameter $Q_{\rm sat}$: $L_{\rm sym}=37.2^{+9.2}_{-8.9}$ MeV, $K_{\rm sym}=-85^{+82}_{-70}$ MeV, and $Q_{\rm sat}=318^{+673}_{-366}$ MeV. For the first time, we measure the values of the empirical parameters $K_{\rm sym}$ and $Q_{\rm sat}$. These values are only weakly impacted by our assumptions, such as the distances or the number of free empirical parameters, provided they are taken within a reasonable range. We also study the weak sensitivity of our results to the set of sources analyzed, and we identify a group of sources that dominates the constraints. The resulting masses and radii obtained are also discussed in the context of the independent constraints from GW 170817 and its electromagnetic counterpart, AT 2017gfo.

Figures

Tables

Table 1. Table 1: Observational information on the 7 qLMXB sources considered in our analysis.

Globular	R.A.aaCoordinates of the qLMXB in each of the GC.	Decl.aaCoordinates of the qLMXB in each of the GC.	XMM Exp.	Chandra Exp.	S/N	GroupbbThe groups A and B denote the sources with a high S/N ($>60$) and lower S/N (), respectively. The 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).	Distances	Distances [8]
Cluster host	(J2000)	(J2000)	time (ks)	time (ks)			Dist #1 (kpc)	Dist #2 (kpc)
47Tuc (X-7)	00:24:03.53	–72:04:52.2	0	181	122	A,A’	$4.53\pm 0.08$ [1]	$4.50\pm 0.06$
M28	18:24:32.84	–24:52:08.4	0	327	113	A,A’	$5.5\pm 0.3$ [2,3]	$5.50\pm 0.13$
NGC 6397	17:40:41.50	–53:40:04.6	0	340	82	A,A’	$2.51\pm 0.07$ [4]	$2.30\pm 0.05$
$\omega$ Cen	13:26:19.78	–47:29:10.9	36	291	49	B,B’	$4.59\pm 0.08$ [5]	$5.20\pm 0.09$
M13	16:41:43.75	+36:27:57.7	29	55	36	B,A’	$7.1\pm 0.62$ [6]	$7.10\pm 0.10$
M30	21:40:22.16	–23:10:45.9	0	49	32	B,B’	$8.2\pm 0.62$ [6]	$8.10\pm 0.12$
NGC 6304	17:14:32.96	–29:27:48.1	0	97	28	B,B’	$6.22\pm 0.26$ [7]	$5.90\pm 0.14$

Table 2. Table 2: Standard values and domain of variation of the empirical parameters considered in this analysis; taken from Margueron et al. ( 2018a ) . See Section  3.1 for the description of the parameters.

Emp. param.	$E_{\mathrm{sat}}$	$E_{\mathrm{sym}}$	$n_{\mathrm{sat}}$	$L_{\mathrm{sym}}$	$K_{\mathrm{sat}}$	$K_{\mathrm{sym}}$	$Q_{\mathrm{sat}}$	$Q_{\mathrm{sym}}$	$Z_{\mathrm{sat}}$	$Z_{\mathrm{sym}}$	$m^{\ast }$	$\Delta m^{\ast }$
	(MeV)	(MeV)	(fm-3)	(MeV)	(MeV)	(MeV)	(MeV)	(MeV)	(MeV)	(MeV)	($m_N$)	($m_N$)
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	–	–	–	–	–

Table 3. Table 3: Distribution of empirical parameters L sym subscript 𝐿 sym L_{\rm{\rm sym}} , K sym subscript 𝐾 sym K_{\rm{\rm sym}} and Q sat subscript 𝑄 sat Q_{\rm{\rm sat}} for various cases. Group A contains high S/N sources (peaked masses): NGC6397, 47-Tuc, M28. Group B contains low S/N sources (flat masses): ω 𝜔 \omega  Cen , NGC6304, M13 and M30. Group A ′ contains sources with peaked masses: NGC6397, 47-Tuc, M28 and M13. Group B ′ contains sources with almost-flat masses: ω 𝜔 \omega  Cen , NGC6304 and M30. See text for more details.

Framework	Sources	Distances	prior	$L_{\mathrm{sym}}$	$K_{\mathrm{sym}}$	$Q_{\mathrm{sat}}$	$R_{1.45}$	${\chi }_{\nu }^2$	nb. of	d.o.f.
			$L_{\mathrm{sym}}$	(MeV)	(MeV)	(MeV)	(km)		param.
1	all	Dist #2	yes	37.2${}_{-8.9}{}^{}{}_{}{}^{+9.2}$	-85${}_{-70}{}^{}{}_{}{}^{+82}$	318${}_{-366}{}^{}{}_{}{}^{+673}$	$12.35\pm 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\pm 0.34$	1.07	49	1126
3	all	Dist #1	yes	38.6${}_{-8.7}{}^{}{}_{}{}^{+9.2}$	-95${}_{-36}{}^{}{}_{}{}^{+80}$	$300$	$12.25\pm 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\pm 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\pm 0.41$	1.08	43	700
6	all/NGC6397	Dist #1	yes	42.6${}_{-9.5}{}^{}{}_{}{}^{+9.9}$	-77${}_{-96}{}^{}{}_{}{}^{+129}$	623${}_{-544}{}^{}{}_{}{}^{+757}$	$12.58\pm 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\pm 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\pm 0.29$	1.04	21	874
9	A′	Dist #2	yes	37.5${}_{-8.9}{}^{}{}_{}{}^{+9.0}$	-88${}_{-70}{}^{}{}_{}{}^{+76}$	263${}_{-361}{}^{}{}_{}{}^{+764}$	$12.22\pm 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\pm 0.43$	1.19	28	255
11	B′	Dist #2	yes	50.3${}_{-9.6}{}^{}{}_{}{}^{+9.8}$	-1${}_{-143}{}^{}{}_{}{}^{+134}$	881${}_{-705}{}^{}{}_{}{}^{+671}$	$12.98\pm 0.40$	1.18	23	178

Table 4. Table 4: Distributions of all the 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 subscript 𝐿 sym L_{\rm sym} , variation over the 3 empirical parameters L sym subscript 𝐿 sym L_{\rm{\rm sym}} , K sym subscript 𝐾 sym K_{\rm{\rm sym}} and Q sat subscript 𝑄 sat Q_{\rm{\rm sat}} ).

Parameter	Source	Dist #1	Dist #2
	M13	0.15 – 0.96	0.15 – 0.95
	$\omega$ Cen	0.23 – 0.97	0.20 – 0.97
	47Tuc	0.15 – 0.65	0.16 – 0.67
pile-up $\alpha$	M28	0.38 – 0.54	0.39 – 0.55
	M30	0.27 – 0.97	0.24 – 0.97
	NGC6304	0.20 – 0.97	0.20 – 0.97
	NGC6397	0.33 – 0.81	0.38 – 0.85
	M13	1.56 – 4.44	2.47 – 5.43
	$\omega$ Cen	17.19 – 21.01	15.76 – 19.32
	47Tuc	2.86 – 4.64	2.93 – 4.76
$N_H\left({\mathrm{\hspace{0.17em}10}}^{20}{\mathrm{cm}}^{-2}\right)$	M28	35.17 – 37.98	35.34 – 37.88
	M30	3.01 – 6.37	3.18 – 6.47
	NGC6304	44.96 – 56.61	45.98 – 58.13
	NGC6397	16.60 – 17.73	17.38 – 19.44
	M13	79.06 – 90.42	76.04 – 86.57
	$\omega$ Cen	70.48 – 80.96	73.62 – 86.05
	47Tuc	104.73 – 112.71	104.54 – 112.79
$k{T}_{\mathrm{eff}}\left(\mathrm{eV}\right)$	M28	109.01 – 118.76	108.49 – 116.35
	M30	86.69 – 99.99	86.51 – 99.17
	NGC6304	91.93 – 108.39	90.26 – 105.36
	NGC6397	62.87 – 68.88	61.16 – 65.88
	M13	0.77 – 1.95	0.74 – 1.95
	$\omega$ Cen	0.80 – 2.01	0.85 – 2.06
	47Tuc	0.66 – 1.44	0.67 – 1.47
$M_{\mathrm{NS}}\left(M_{\odot }\right)$	M28	0.70 – 1.51	0.68 – 1.43
	M30	0.78 – 2.00	0.80 – 1.99
	NGC6304	0.87 – 2.07	0.85 – 2.06
	NGC6397	0.72 – 1.62	0.70 – 1.47
	M13	7.72 – 8.50	7.07 – 7.19
	$\omega$ Cen	4.50 – 4.65	5.18 – 5.25
	47Tuc	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
	NGC6304	6.15 – 6.43	5.86 – 6.01
	NGC6397	2.49 – 2.56	2.29 – 2.34
	M13	4.55 – 15.57	5.64 – 16.48
	$\omega$ Cen	1.21 – 5.20	0.89 – 4.87
	47Tuc	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
	NGC6304	3.87 – 15.78	4.42 – 15.66
	NGC6397	2.83 – 8.19	2.89 – 8.30
$C_1$	M13	0.89 – 1.09	0.89 – 1.09
	$\omega$ Cen	0.99 – 1.16	0.99 – 1.18
$C_2$	M13	0.79 – 0.96	0.78 – 0.99
	$\omega$ Cen	0.96 – 1.12	0.96 – 1.12

Equations

$$\mbox{$R_{\infty}$}=\mbox{$R_{\rm NS}$}\left(1+z\right)=\mbox{$R_{\rm NS}$}\left(1-\frac{2G\mbox{$M_{\rm NS}$}}{\mbox{$R_{\rm NS}$}c^{2}}\right)^{-1/2}.$$

$$e_{\mathrm{sat}}=E_{\mathrm{sat}}+\frac{1}{2}{K}_{\mathrm{sat}}{x}^2+\frac{1}{3!}{Q}_{\mathrm{sat}}{x}^3+\frac{1}{4!}{Z}_{\mathrm{sat}}{x}^4+...,$$

$$e_{\mathrm{sym}}=E_{\mathrm{sym}}+L_{\mathrm{sym}}x+\frac{1}{2}{K}_{\mathrm{sym}}{x}^2+\frac{1}{3!}{Q}_{\mathrm{sym}}{x}^3+\frac{1}{4!}{Z}_{\mathrm{sym}}{x}^4+...,$$

$$e_{\mathrm{sym}}=\frac{1}{2}\frac{{\partial }^{2}e}{\partial {\delta }^2},$$

$$e(n,\delta )=t(n,\delta )+v(n,\delta ).$$

$$t(n,\delta)=\frac{t_{\rm sat}}{2}\left(\frac{n}{n_{\rm sat}}\right)^{2/3}\Big{[}\left(1+\kappa_{\rm sat}\frac{n}{n_{\rm sat}}\right)f_{1}(\delta)+\kappa_{\rm sym}\frac{n}{n_{\rm sat}}f_{2}(\delta)\Big{]},$$

$$f_1(\delta )$$

$$f_2(\delta )$$

$$v(n,\delta )=\sum _{\alpha =0}^N\left(v_{\alpha }^{\mathrm{is}}+{\delta }^2{v}_{\alpha }^{\mathrm{iv}}\right)\frac{x^{\alpha }}{\alpha !}u(x),$$

$$u(x)=1-(-3x{)}^{N+1-\alpha }\exp (-bn/n_{\mathrm{sat}}).$$

$${\chi }_{D,i}^2=\frac{{\left(D_{\mathrm{mcmc},\mathrm{i}}-D_{\mathrm{data},\mathrm{i}}\right)}^2}{{\left({\sigma }_i\right)}^2}.$$

Full text

dense matter — equation of state — stars: neutron

††software: emcee (Foreman-Mackey et al., 2013), corner (Foreman-Mackey, 2016), HEAsoft (Nasa High Energy Astrophysics Science Archive Research Center (2014), Heasarc), Xspec (and PyXspec, Arnaud 1996), XMMSAS (Gabriel et al., 2004), and CIAO (Fruscione et al., 2006). ††thanks: sebastien.guillot@irap.omp.eu

1 Introduction

Determining the equation of state of dense matter (EoS) – the relation between the pressure $P$ and energy-density $\rho$ beyond the nuclear saturation energy-density \rho_{\rm sat}\sim 2.4\mbox{\times 10^{14}}\mbox{,{\rm g\mbox{$\,{\rm cm^{-3}}$}}} – is an important goal of fundamental physics and astrophysics, with far-reaching implications. Observations of neutron stars (NSs) offer extraordinary tools to investigate dense matter properties, which are complementary to experimental studies (e.g., Lattimer & Prakash, 2007; Kramer, 2008; Lattimer & Prakash, 2010; Lattimer & Lim, 2013; Hebeler et al., 2013). For instance, macroscopic properties of NSs, such as masses, radii, moments of inertia or tidal deformabilities, provide constraints on dense matter at energy-densities beyond $\rho_{\rm sat}$ (e.g., Lattimer et al., 1990; Lattimer & Schutz, 2005; Flanagan & Hinderer, 2008; Lattimer & Prakash, 2010; Abbott et al., 2018).

A variety of methods exist to constrain the EoS from NSs. Besides electromagnetic observations described below, the recent observation of the gravitational wave signal from a NS-NS merger and its electromagnetic counterpart has been analyzed to better constrain the stiffness of matter inside NSs. Specifically, the signal GW 170817, detected by the LIGO and Virgo gravitational wave (GW) detectors on 2017 August 17th, resulted in constraints on the tidal deformability of the NSs from the quadrupole moment in the space-time surrounding the NS merger (Abbott et al., 2017). Following the discovery of GW 170817, several articles proposed constraints on the EoS and the radius of these NSs using information from the GW signal and the simultaneous GRB 170817, including its afterglow AT 2017gfo (e.g., Bauswein et al. 2017; Radice et al. 2018; Annala et al. 2018; Raithel et al. 2018; Tews et al. 2018b; De et al. 2018, with the most recent one from Abbott et al. 2018). The conclusions of these papers are consistent, although the real quantitative information extracted from this first ever detection may not yet compete with nuclear physics knowledge (Tews et al., 2018b). The future of this detection method is however promising and will certainly constrain present EoS models.

All other methods to constrain the EoS make use of electromagnetic observations of NSs. More generally, they rely on mass $M_{\rm NS}$ and radius $R_{\rm NS}$ measurements (or other related properties). For example, the modelling of the pulse profile of millisecond pulsars (MSP) can provide measurements of $M_{\rm NS}$ and $R_{\rm NS}$ (e.g., Bogdanov et al., 2007; Bogdanov, 2013). The currently operating Neutron Star Interior Composition ExploreR (NICER) is routinely observing MSPs with this aim (Gendreau et al., 2016; Gendreau & Arzoumanian, 2017). Measuring the moment of inertia of pulsars using radio timing observations of pulsars in binary systems via spin-orbit coupling effects is also being envisaged to constrain the EoS (e.g., Lattimer & Schutz, 2005; Kramer, 2008). Finally, the thermal emission from NSs provides a promising technique to obtain $M_{\rm NS}$ and $R_{\rm NS}$ (see Miller, 2013; Heinke, 2013; Özel & Freire, 2016, for recent reviews). While this could, in principle, be achieved with all cooling NSs, some of them may be affected by systematic uncertainties that may alter the measurements. For example, the spectral modeling of X-ray dim isolated NSs may be complicated by uncertainties about their atmospheres, their magnetic field B\sim\mbox{10^{11-12}}\mbox{,G}, and the presence of X-ray pulsations indicating a non-uniform surface emission (e.g., Pons et al., 2002); which may require phase-resolved spectroscopy (Hambaryan et al., 2017). Similarly, central compact objects are likely affected by the same effects, although not all CCOs show pulsations (e.g., Klochkov et al., 2015).

The cooling tails of Type-I bursts from NSs in X-ray binaries have also been used for EoS constraints (e.g., Suleimanov et al., 2011; Nättilä et al., 2016). Furthermore, when these bursts reach the Eddington flux, the peak flux provides an additional observable with which to break the degeneracy between $M_{\rm NS}$ and $R_{\rm NS}$ (e.g. Özel et al., 2010; Güver & Özel, 2013). However, difficulties may arise from the spectral modeling with a Planck function and the use of a color correction from theoretical atmosphere models (see Güver et al., 2012b, a; Kajava et al., 2014; Nättilä et al., 2016; Özel et al., 2016, for discussions). To remedy these issues, recent work fitted such atmosphere models to each spectrum during the cooling tail of the NS 4U 1702–429 to obtain $M_{\rm NS}$ and $R_{\rm NS}$ measurements (Nättilä et al., 2017), instead of relying on color corrections.

All methods using thermally emitting NSs require precise knowledge of the source distances. For this reason, quiescent low-mass X-ray binaries (qLMXBs) located inside globular clusters (GCs) have provided reliable constraints on the EoS. The distances to GCs can be measured independently with uncertainties of $\sim$5–10% (Harris, 1996, 2010); compared to the $\sim$ 30–50% uncertainties of LMXBs in the field of the Galaxy. Furthermore, qLMXBs present other advantages that we describe in Section 2. While initially these sources where analyzed individually to constrain the EoS (e.g., Heinke et al., 2006; Webb & Barret, 2007; Guillot et al., 2011), it has become clear in recent years that statistical analyses combining multiple qLMXBs would provide more useful constraints on dense matter (Guillot et al., 2013; Guillot & Rutledge, 2014; Guillot, 2016; Lattimer & Steiner, 2014; Özel et al., 2016; Bogdanov et al., 2016; Steiner et al., 2018).

This article presents one such analysis in which the spectra of a sample of qLMXBs are simultaneously analyzed to constrain the EoS. Because the red-shifted radius, measured from the modelling of the observed spectrum, depends on both the gravitational mass and the physical radius, a simultaneous analysis of several qLMXB sources can help break degeneracies between these two properties of NSs, assuming these objects are governed by the same $M_{\rm NS}$–$R_{\rm NS}$ relation, i.e., the same EoS. This can in turn be used to infer the properties of the dense NSs matter. For practical reasons, such method requires parameterizing the EoS, i.e., representing it as a function of some parameters either in $M_{\rm NS}$–$R_{\rm NS}$ space, or in $P$–$\rho$ space. Previous work used analytical parameterizations, such as a toy-model constant-$R_{\rm NS}$ (Guillot et al., 2013; Guillot & Rutledge, 2014; Guillot, 2016) or piecewise polytrope representations111A sequence of connected power laws, $P=k_{i}\rho^{\gamma_{i}}$, where $i$ typically runs up to 3 or 5 (e.g., Read et al., 2009; Raithel et al., 2016). (Lattimer & Steiner, 2014; Özel et al., 2016; Bogdanov et al., 2016; Steiner et al., 2018).

In this work, we employ a representation of the EoS based on nuclear physics empirical parameters. The model is presented in Margueron et al. (2018a, b) and offers the possibility to easily incorporate nuclear physics knowledge. In Section 2, we summarize the characteristics of qLMXBs and present the reasons that make them ideal sources for EoS constraints. We also describe the data reduction and spectral extraction of our qLMXBs sample, as well as the surface emission model of these NSs. Section 3 summarizes various aspects of the EoS meta-model of Margueron et al. (2018a, b) that we used to fit our spectral data of qLMXBs. Section 4 presents the Markov-Chain Monte Carlo (MCMC) approach used to find the best fit EoS model to the qLMXBs spectra and Section 5 presents the results, and compares them with previous constraints on the EoS. Finally, the conclusions in Section 6 summarize this work.

2 Thermal emission from quiescent low-mass X-ray binaries

In this section, we detail our present understanding of qLMXB thermal emission in GCs as well as host GC distance measurements. We also give details on our X-ray spectral data analysis and spectral model.

2.1 Low-mass X-ray binaries in quiescence

The surface emission from NSs in qLMXBs is now routinely used to obtain measurements of $M_{\rm NS}$ and $R_{\rm 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 red-shifted 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:

[TABLE]

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

The qLMXBs inside GCs also present the other advantage of exhibiting a remarkable flux stability at all timescales (Heinke et al., 2006; Guillot et al., 2011; Servillat et al., 2012; Heinke et al., 2014). While LMXBs in the field of the Galaxy often exhibit flux variability, attributed to changes in the non-thermal 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 non-thermal 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 sub-groups 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$\,{\rm sec}$ (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 is thought to be weak, as supported by the fact that their presumed descendants, millisecond pulsars (Alpar et al., 1982; Bhattacharya & van den Heuvel, 1991; Tauris & van den Heuvel, 2006), have inferred dipole fields $B\sim 10^{8}$–10^{9}\mbox{,G}, compared to $10^{11}$–10^{12}$$\,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 non-magnetic. 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_{\rm NS}$ and $R_{\rm NS}$.

It is generally accepted that the atmosphere of a 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^{3}$–$10^{4}$ yrs (Chang & Bildsten, 2004; Chang et al., 2010), whereas the atmosphere (of thickness \sim 1\mbox{,{\rm cm}} and mass M_{\rm atm}\sim 10^{-20}\mbox{,{\rm 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}\mbox{\rm,\mbox{$\,{\rm M}_{\odot}$}\mbox{$\,{\rm 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 $\omega$ Cen (Haggard et al., 2004). Searches for H$\alpha$ 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 a He atmosphere instead. This conclusion was supported by the low $R_{\rm NS}$ found from the spectral analyses with a H atmosphere (\sim 8\hbox{,{\rm km}}, in the earlier work of Guillot et al. 2011, 2013), while a He atmosphere resulted in a $R_{\rm 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 pile-up (an instrumental effect, Davis 2001) in the Chandra X-ray Observatory spectra of this qLMXB is sufficient to yield radii in the range222It was demonstrated that pile-up effects, even at the 1%-level, can significantly shift the peak of the thermal spectrum to higher energies, and therefore result in underestimated radii (Bogdanov et al., 2016). 10–11$\,{\rm km}$.

With all 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, $\omega$ Cen (NGC 5139), and 47 Tuc (NGC 104) X-7. We excluded 47 Tuc X-5 because of its eclipses, its flux variability, and variable line-of-sight absorption, which make the spectral modelling rather uncertain (Bogdanov et al., 2016). Some information about these sources is detailed in Table 1.

2.2 On the distances of globular clusters

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 have now become a tangible reality, thanks to the exquisite data provided by the European Space Agency’s (ESA’s) Gaia space mission (Gaia Collaboration et al., 2016). Within the framework of Gaia’s Data Release 2 (DR2; Gaia Collaboration et al., 2018a), 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. (2018c), 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., 2018b), 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. (2018b, c) 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. (2018b) 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. (2018b). 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 $\omega$ 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\mbox{,{\rm 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.

2.3 X-ray spectral data analysis and spectral model

The processing of the XMM-Newton and Chandra data sets is performed with the XMMSASv15.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. For the Galactic absorption, we used the recent model tbabs (Wilms et al., 2000). We also add a power-law component to account for possible excess of counts above 2$\,{\rm keV}$ that may originate from non-thermal 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 pile-up component is also added for all Chandra spectra, even those qLMXBs with low count rates inducing a pile-up fraction $\lesssim 1\%$. As was pointed out by Bogdanov et al. (2016), uncorrected pile-up, even at low pile-up fraction $\sim 1\%$, can significantly bias the radius measurement. Specifically, for NGC 6397, the low $R_{\rm NS}$ obtained with H atmosphere models was a consequence of the unmodelled pile-up of photons in the X-ray spectra.

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

• •

the parameter $\alpha$ in the pileup model,

• •

the column density of neutral hydrogen $N_{\rm H}$, from the tbabs model,

• •

the NS surface temperature $kT_{\rm eff}$ in the nsatmos model,

• •

the NS mass in the nsatmos model,

• •

the NS radius in the nsatmos model,

• •

the NS distance (set as a prior; see Table 1) in the nsatmos model,

• •

the power-law normalization (model powerlaw with fixed $\Gamma=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 commonly done333In those case, the constant for XMM-pn is fixed to unity, while the ones for the XMM-MOS and Chandra spectra, $C_{1}$ and $C_{2}$ respectively, are fitted parameters.. 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.

3 The dense matter equation of state

For the present analysis, the dense matter EoS is provided by a meta-modeling described in Margueron et al. (2018a, b), instead of the toy-model constant-$R_{\rm NS}$ representation of the EoS (Guillot et al., 2013; Guillot & Rutledge, 2014; Guillot, 2016), or instead of the polytropes (Steiner et al., 2013; Özel et al., 2016; Steiner et al., 2018) used in previous works. The meta-modeling employed here is able to accurately reproduce existing nucleonic EoSs and smoothly interpolate between them. It is based on a Taylor expansion in the baryon density $n=n_{n}+n_{p}$, where $n_{n}$ and $n_{p}$ are the neutron and proton densities, around the nuclear saturation density $n_{\mathrm{sat}}\approx 0.16$ fm*-3*. Note that the nuclear saturation density is expressed as baryon number per unit volume and it coincides with the energy-density $\rho_{\mathrm{sat}}$ introduced previously. Such an approach is realistic up to 3–4 $n_{\rm sat}$, where one could expect the onset of new degrees of freedom (hyperons, quarks, pion condensation, etc). This meta-model may therefore break down for high-mass NSs (at around or above 2$\,{\rm M}_{\odot}$). Fortunately, these high masses seem not to be favored in the present analysis and for the present sources. For completeness, we briefly describe our modeling for the crust and the core of the NSs in this section.

3.1 Equation of state for cold catalyzed neutron stars

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*-3*, 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 by another one, such as the FPS one. These two tabulated EoS can be downloaded from the following website444http://www.ioffe.ru/astro/NSG/NSEOS/.

For numerical reasons, the transition between the crust and the core is guided within and log$\rho$–log$P$ cubic spline matching the values and derivatives at both boundaries. The two boundaries are taken to be $n_{\rm sat}/10$ for the lower bound and $n_{\rm sat}$ for the upper one. 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., 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_{\rm sat})/(3n_{\rm sat})$ of the energy per particle in symmetric matter,

[TABLE]

and of the symmetry energy per particle

[TABLE]

where the symmetry energy $e_{\rm sym}$ is defined as the isospin polarization energy

[TABLE]

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

$E_{\rm sat}$ and $E_{\rm sym}$ are the saturation and symmetry energy at the saturation density $n_{\rm sat}$. $L_{\rm sym}$ is the slope of the symmetry energy, and since the saturation is an equilibrium point, there is no slope of the energy per particle in symmetric matter. $K_{\rm sat/sym}$ stands for the curvature, $Q_{\rm sat/sym}$ for the skewness, and $Z_{\rm sat/sym}$ for the kurtosis of the energy per particle in symmetric matter and of the symmetry energy, respectively. The values of these empirical parameters are determined from experimental measures, with different accuracies. Reviews of their experimental determination can be found in Margueron et al. (2018a) and in references therein.

We consider the meta-modeling ELFc proposed in Margueron et al. (2018a), 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

[TABLE]

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

[TABLE]

where $t_{\rm sat}=3\hbar^{2}/(10m)(3\pi^{2}/2)^{2/3}n_{\rm sat}^{2/3}$, $m$ is the nucleon mass, and the functions $f_{1/2}$ are defined as,

[TABLE]

The potential term is expressed as,

[TABLE]

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

[TABLE]

Fixing $b=10\ln 2\approx 6.93$, as in Margueron et al. (2018a), implies that the function $u$ converges quickly to $1$ as the density increases from 0. It ensures that $v(n,\delta)\rightarrow 0$ for $n\rightarrow 0$ for any order $N$. The larger $N$, the smaller the correction $u(x)$. The parameters $v_{\alpha}^{\rm is/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 EoS increases as the order $N$ increases. For $N=4$, the meta-modeling can very accurately (at the % accuracy, in the worst case) reproduce binding energy, pressure, and sound velocity of a large number of existing EoS up to $4n_{\rm sat}$, as shown in Margueron et al. (2018a).

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_{\rm sat}=-15.8$ MeV, $E_{\rm sym}=32$ MeV, $n_{\rm sat}=0.155$ fm*-3* and $K_{\rm sat}=230$ MeV. The parameters $\kappa_{\rm sat/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_{\rm NS}$–$R_{\rm NS}$ relation is known to be mostly influenced by the empirical parameters $L_{\rm sym}$, $K_{\rm sym}$, and $Q_{\rm sat}$, since the EoS in the density range going from $n_{\rm sat}$ to approximately $3n_{\rm sat}$ depends most strongly on them (Margueron et al., 2018b). $L_{\rm sym}$ and $K_{\rm sym}$ (respectively, $Q_{\rm 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_{\rm NS}$ and $R_{\rm NS}$ measurements for high-mass NSs.

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 $\,{\rm M}_{\odot}$, which implies that the central densities of these NSs are not very large, and the meta-model can reasonably be applied.

3.2 The effect of the empirical parameters on the $M_{\rm NS}$–$R_{\rm NS}$ relation

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

The EoS selection criteria include those which satisfy the requirements of causality and positiveness of the symmetry energy, as well as being compatible with a maximum mass above 1.9\mbox{,{\rm M}{\odot}}. This mass limit corresponds approximately to the $2\sigma$ lower limits of the measurements for PSR J1614–2230, 1.908\pm 0.016\mbox{,{\rm 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\mbox{,{\rm M}_{\odot}}.

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

While there are various EoSs which pass through a point in the $M_{\rm NS}$–$R_{\rm NS}$ diagram, their paths are different. The degeneracy between different EoSs thus requires the knowledge of a set of $M_{\rm NS}$–$R_{\rm NS}$ points, as distant as possible from each other. In conclusion of this analysis, the empirical parameters $L_{\rm sym}$, $K_{\rm sym}$, and $Q_{\rm sat}$ allow the exploration of a wide domain of $M_{\rm NS}$ and $R_{\rm NS}$ with various paths. Therefore, it may be possible to constrain the values of these parameters by confronting them to the observational data from the thermal X-ray emission of NSs.

4 Confronting the equation of state with the data

In this section, we detail the methodology of our analysis: we employ an MCMC approach with the stretch-move algorithm (Goodman & Weare, 2010) to consistently analyze the seven qLMXB sources, and the nuclear matter EoS meta-modeling is included directly in the analysis. The result is that the astrophysical (NS properties) and nuclear physics (EoS) parameters are adjusted together, without over-constraining one or the other. We can solve the so-called inverse problem and obtain constraints on the EoS properties directly from the data analysis. This is the first time that the thermal emission from NSs is analyzed in this manner.

4.1 MCMC approach with the stretch-move algorithm

For all 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\sigma$ 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_{\rm{mcmc,i}}$ and the actual measured data $D_{\rm{data,i}}$ (from Table 1), taking into account their standard deviations $\sigma_{i}$. The distance penalty reads $\chi^{2}_{D}=\sum_{i=0}^{N}\chi^{2}_{D,i}$, where the $\chi^{2}_{D,i}$ for each source are given by:

[TABLE]

The MCMC approach permits efficient sampling of our parameter space with high dimensionality: 49 parameters in total, including 3 nuclear physics EoS parameters, plus 6 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), plus 4 multiplicative normalization constants (for the cross-calibration between the XMM-pn, XMM-MOS and Chandra, for the qLMXBs in M13 and $\omega$ 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 flow-chart in Figure 2):

• •

Step 0: a large number of chains or ”walkers” are initialized, each one corresponding to a random point in the multi-dimensional 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).

• •

Step 1: we solve the TOV equations for each walker, providing 426 $M_{\rm NS}$–$R_{\rm NS}$ relations at each iteration.

• •

Step 2: for each walker, the sampled masses of the seven NSs are associated to seven calculated radii according to the $M_{\rm NS}$–$R_{\rm NS}$ relation. Using those $M_{\rm NS}$ and $R_{\rm NS}$ and the other astrophysical parameters, we calculate the global $\chi^{2}$ between the emission models (NS atmosphere) and the data for the seven NSs.

• •

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), that have been chosen to $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); Foreman-Mackey et al. (2013).

• •

Step 4: Steps 1 to 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.

• •

Step 5: When the MCMC loop stops, the statistical posterior distributions are calculated and marginalized to create the outputs.

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_{\rm NS}$ toy-model. In addition to its simplicity, this test is interesting since it allows us to compare with results that have already been reported in the literature.

4.2 Tests using a constant radius toy model

We first consider the constant-$R_{\rm NS}$ model (Guillot et al., 2013), which assumes that all NSs have the same radius, i.e., that the EoS is represented in $M_{\rm NS}$–$R_{\rm NS}$ space by a vertical line in which $R_{\rm NS}$ is independent of $M_{\rm NS}$ (which remain as free parameters). This is a simple toy-model approximation motivated by the observations that most nucleonic EoSs (the ones consistent with 2$\,{\rm M}_{\odot}$) have a rather weak dependence on $M_{\rm NS}$, between 1$\,{\rm M}_{\odot}$ and 2$\,{\rm 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_{\rm 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 \mbox{R_{\rm NS}}=11.09^{+0.38}_{-0.36}\hbox{,{\rm km}} (Dist #1, for \mbox{\chi^{2}{\nu}}=1.06) and \mbox{R{\rm NS}}=11.04^{+0.39}_{-0.35}\hbox{,{\rm km}} (Dist #2, for \mbox{\chi^{2}_{\nu}}=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 (47Tuc X-7) and new data, the use of recent distance measurements, the improvement of the analysis (e.g., the new absorption model tbabs), and the inclusion of the pile-up correction model for all sources (including those with a $\sim 1\%$ pile-up 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 non-negligible probabilities for a NS radius larger than about 11$\,{\rm km}$.

There is however also a large fraction of the posterior probability distribution which is located below 11$\,{\rm km}$, in conflict with nuclear physics expectations (e.g., Margueron et al., 2018b; Tews et al., 2018b) as well as our illustrative Figure 1. For instance, one could deduce from these figures that \mbox{R_{\rm NS}}\lesssim 11\hbox{,{\rm km}} requires L_{\rm sym}\lesssim 20\mbox{,{\rm 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.

5 Results

In this section, the main results of our novel approach are presented and discussed.

5.1 Framework

We remind the reader that the main features of our work are that i) we fit simultaneously seven NS qLMXB sources, ii) we impose the same EoS to all these sources, and iii) we 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), and they are kept fixed in the present analysis: $n_{\rm sat}$, $E_{\rm sat}$, $E_{\rm sym}$, and $K_{\rm sat}$ (values in Table 2). The more influential and less known parameters, $L_{\rm sym}$, $K_{\rm sym}$ and $Q_{\rm sat}$, are fitted in our analysis. The values of $K_{\rm sym}$ and $Q_{\rm sat}$ are currently unknown, while there exist constraints on $L_{\rm sym}$ from nuclear experiments and nuclear theoretical predictions. These constraints indicate that $L_{\rm sym}$ has a value around 50 MeV with an uncertainty of about $\pm$10 MeV (Lattimer & Lim, 2013). We therefore incorporate this knowledge from nuclear physics by considering a Gaussian prior on $L_{\rm sym}$ centered at 50 MeV with a width of 10 MeV. Since $K_{\rm sym}$ and $Q_{\rm sat}$ are unknown, we consider a uniform distribution in the wide ranges listed in Table 2. The higher-order empirical parameters, $Q_{\rm sym}$ and $Z_{\rm sat/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., 2018b).

5.2 Main results

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\pm 0.02$, for 1126 degrees of freedom, and the best fit corresponds to \mbox{\chi^{2}}=1.08, giving a null hypothesis probability of 3.1%.

The marginalized posterior probabilities for the empirical EoS parameters $L_{\rm sym}$, $K_{\rm sym}$ and $Q_{\rm sat}$ are shown in Figure 4. We observed from the marginalized distributions that $L_{\rm sym}$ peaks at lower values than the one imposed by the prior (50\pm 10\mbox{,{\rm MeV}}), but remains consistent with it: L_{\rm sym}=37.2^{+9.2}_{-8.9}\mbox{,{\rm MeV}}. This somewhat reflects the tension driving the fit towards low radii at low masses (see Figure 1 and related discussion). The empirical parameter K_{\rm sym}=-85^{+82}_{-70}\mbox{,{\rm 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\mbox{,{\rm MeV}} extracted from analysis of chiral effective field theory (EFT) calculations (Margueron et al., 2018a). Finally, the empirical parameter Q_{\rm sat}=318^{+673}_{-366}\mbox{,{\rm 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_{\rm sym}$ and $Q_{\rm sat}$, this is the first time that these parameters are extracted from data.

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

The $M_{\rm NS}$–$R_{\rm NS}$ posterior probability distributions corresponding to the MCMC runs with the distances of set Dist #1 (upper panel) and set Dist #2 (lower panel) are displayed in Figure 5. It is reassuring to notice that the $M_{\rm NS}$–$R_{\rm NS}$ posterior probability distribution is almost insensitive to the set of distances considered, as was also observed for the constant-$R_{\rm 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 $\sim 11.5$ and 13.0\hbox{,{\rm km}} for a 1.4$\,{\rm M}_{\odot}$ NS, with a preference for low masses, although the 90% credible intervals are compatible with 2$\,{\rm M}_{\odot}$.

In Figure 6 (top panels), we present the 90% credible interval $M_{\rm NS}$–$R_{\rm 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\mbox{,{\rm M}_{\odot}} or higher. In a few cases, the 90% credible intervals only reach masses around 1.4–1.5$\,{\rm M}_{\odot}$, which appear to favor low-mass NSs. However, this is compatible with current distribution of MSP masses (Antoniadis et al., 2016; Özel & Freire, 2016) which descend from LMXBs; we note that the lowest known NS mass is 1.174$\pm$0.004$\,{\rm 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.

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_{\rm NS}$–$R_{\rm 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_{\rm NS}$–$R_{\rm NS}$ constraints from single qLMXBs, for which a significant portion of the $M_{\rm NS}$–$R_{\rm NS}$ contours appear at constant $R_{\rm NS}$.). Therefore, unless the mass of the NS is measured independently, the degeneracy between $M_{\rm NS}$ and $R_{\rm 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_{\rm NS}$ and $R_{\rm 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_{\rm NS}$ and $R_{\rm NS}$, see Özel & Freire 2016 for a review). The GW signals of two merging NSs, on the other hand, provides measurements of the merging masses and of the 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 AT2017gfo (Bauswein et al., 2017) and from GW170817 (Annala et al., 2018; Abbott 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 GW 170817, 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 $\rho$ resulting from our analysis, considering the nucleon and lepton contributions in $\beta$-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 for the EoS based on different constraints and labelled 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 2 and 3-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 from the MM prediction could potentially further reduce the possibilities for the EoS.

In Figure 8, we show both the squared sound speed in units of the speed of light, $(v_{s}/c)^{2}$, and the symmetry energy, $e_{\rm sym}$, as functions of the energy density $\rho$. The coloured 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 $(v_{s}/c)^{2}$ and $e_{\rm sym}$, as for the general EoS shown in Figure 7. Similarly to Figure 7, there is also a good overlap of our predictions with the MM ones.

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 distances 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\mbox{,{\rm M}_{\odot}} at 98% confidence. The best-fit absorption values $N_{\rm 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 power-law normalizations $N_{\rm pl}$ obtained 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 non-zero probabilities for $N_{\rm pl}=0$. This lends further evidence for the absence of non-thermal emission in these objects. We have nonetheless considered the possible existence of non-thermal emission in our analyses by including a power-law component in the spectral model.

5.3 Sensitivity analysis

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 sub-section), 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_{\rm sym}$, $K_{\rm sym}$ and $Q_{\rm sat}$, the radius $R_{1.45}$ for a 1.45\mbox{,{\rm M}_{\odot}} NS, and the best $\chi^{2}_{\nu}$. We also give the number of fitting parameters and the number of degrees of freedom (d.o.f.) for each run. The rows of Table 3 represent the various frameworks considered. The two first rows represent the framework already explored, considering all the seven qLMXB sources, the two sets of distances Dist #1 and Dist #2, the prior on $L_{\rm 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 parameter by fixing Q_{\rm sat}=300\mbox{,{\rm MeV}}. The impact on the centroid and width for $L_{\rm sym}$, $K_{\rm sym}$ and $R_{1.45}$ is marginal. The minimum $\chi^{2}_{\nu}$ changes minimally. For the Framework #4, we replace the Gaussian prior on $L_{\rm sym}$ by a uniform prior ranging from 20 to 120 MeV. The change in the minimum $\chi^{2}_{\nu}$ in marginal, indicating that the fit statistic is not affected by adding/removing the prior on $L_{\rm sym}$. Furthermore, removing the prior on $L_{\rm 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_{\rm sym}$ are compensated by larger values of $K_{\rm sym}$ and $Q_{\rm sat}$. This emerges from the anti-correlation between $L_{\rm sym}$ and $K_{\rm sym}$, as discussed in Section 5.2, and reported in Margueron et al. (2018b). This observed compensation between the empirical parameters originates from the fact that in our meta-model, $e_{\rm sym}$ is constrained to positive values up to a threshold central density that produces 1.9$\,{\rm M}_{\odot}$ NSs. Therefore, if $L_{\rm sym}$ becomes too low, the other parameters (mostly $K_{\rm sym}$) re-adjust to satisfy the condition on $e_{\rm 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 analyse the sensitivity of the result to the modification of the qLMXB source set by removing a single qLMXB. In the Framework #4, #6, and #7, we removed 47Tuc X-7, NGC 6397, or M28, respectively. Since these sources have the largest signal-to-noise ratio (S/N), their removal allows to check to which extent they contribute to drive the results. While marginal, there are indeed some systematic effects: $L_{\rm sym}$ is increased by about 5–6$\,{\rm MeV}$, $K_{\rm sym}$ by 10–20$\,{\rm MeV}$, the value for $Q_{\rm sat}$ is almost doubled, and the radius $R_{1.45}$ is increased by up to 0.23$\,{\rm km}$. These systematic corrections remain inside the original uncertainty estimated for the reference results (Framework #1 and #2). These different approaches are summarized in Figure 6 (bottom) which shows the 90% credible interval of the $M_{\rm NS}$–$R_{\rm NS}$ posterior distributions of Frameworks #1 to #7, and demonstrating that they broadly overlap in the 12–13$\,{\rm 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}$). The groups $A$ and $B$ are defined with respect to the S/N ($A$ for S/N $>$ 60, $B$ otherwise, see Table 1). The 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 the groups $A$ and $A^{\prime}$, as well as $B$ and $B^{\prime}$ are almost identical. The groups $A$ and $A^{\prime}$ prefer the lower values for $L_{\rm sym}$, $K_{\rm sym}$ and $Q_{\rm sat}$ comparable to the reference results. They favor lower radii $R_{1.45}\approx 12.2\pm 0.3$ km. By contrast, the groups $B$ and $B^{\prime}$ tend to increase the values for $L_{\rm sym}$, $K_{\rm sym}$, $Q_{\rm 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_{\rm sym}$ and for the radius $R_{1.45}$. We also note that, for the groups $B$ and $B^{\prime}$, the $L_{\rm sym}$ values are essentially identical to the prior given on that parameter (L_{\rm sym}=50\pm 10\mbox{,{\rm MeV}}); implying that these two groups have little weight in the constraints on $L_{\rm 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_{\rm sym}$. These are the qLMXBs in $\omega$ Cen, M13, M30, and NGC 6304. An improvement in the analysis of the qLMXB thermal emission will require more statistics especially for these sources.

5.4 Comparison with previous work

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}$ meausurements, e.g., Heinke et al. 2006; Guillot et al. 2011, to more consistent approaches. In a first attempt to consistently analyse several qLMXB sources combined, a constant radius EoS model was proposed, inspired by the qualitative behaviour of most of the nuclear EoSs (Guillot et al., 2013; Guillot, 2016).

Because these early results did not consider a full treatment of the pile-up instrumental effects in the Chandra data (which are significant even at low pile-up 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 pile-up and which contain similar inputs as in our analysis. Recently, Steiner et al. (2018) found that the radius of a 1.4\mbox{,{\rm M}_{\odot}} NS is most likely between 10.4 to 13.7$\,{\rm km}$ at 68% confidence level, considering all cases tested in that work. Assuming a pure H atmosphere for all objects, they found $R_{\rm NS}$ in the range 11.2–12.3$\,{\rm 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_{\rm sym}$, Steiner et al. (2018) found 38.94–58.09$\,{\rm 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) determine a $M_{\rm NS}$–$R_{\rm 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$\,{\rm km}$, for masses ranging from 1 to 2$\,{\rm M}_{\odot}$, which is a smaller estimation than ours. In a more recent analysis, Bogdanov et al. (2016) included the same twelve 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$\,{\rm 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 bursts data, which we have not, (3) they used polytropes to parameterize the EoS. Another main difference to 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_{\rm NS}$ approximation prefer radii around \sim 11.1\pm 0.4\hbox{,{\rm km}}, consistent with the estimates in Özel et al. (2016) and Steiner et al. (2018), while including nuclear EoS and a prior on the empirical parameter $L_{\rm sym}$, the radius can increase up to $\sim$12.0–12.5$\,{\rm 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\hbox{,{\rm km}} at 68% credibility, for a mass more difficult to constrain, in the range 1.4–2.2$\,{\rm M}_{\odot}$.

Observations of millisecond pulsars 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>10.4\hbox{,{\rm 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 (priv. communication, and Riley et al., 2019, in press; Miller et al., 2019, submitted to ApJ), 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 \mbox{R_{\rm NS}}=13.1\mbox{{}^{+0.9}_{-0.7}}\hbox{,{\rm km}}, 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 GW 170817 allowed to get an estimation of the radii of the two stars as well as constraints on their EoS through the tidal deformability parameter $\Lambda$ (Abbott et al., 2018). Further analyses of the GW and electromagnetic signals lead to the constraints on the radii drawn in Figure 5. A good agreement with our analysis is also to be noticed.

6 Conclusions

We have used a collection of X-ray spectra coming from seven qLMXBs and have analyzed their surface thermal emission assuming a NS H-atmosphere, and assuming a flexible meta-modeling for the nuclear EoS which has been implemented directly in the fit. For the first time, the emission model and the EoS parameters have been treated on equal footing, avoiding overconstraints which were potentially present in previous analyses. In all our analyses, the instrumental phenomenon of pile-up and the absorption of X-rays in the ISM have been taken into account using the new tbabs absorption model, as well as a power-law component accounting for non-thermal emission. We modeled the surface thermal emission using the NSATMOS model, which requires the mass and the radius of the sources as inputs, so that we can implement the $M_{\rm NS}$–$R_{\rm NS}$ relation, obtained from the EoS parameterization, directly in the spectral modeling. Because of the degeneracy between the radius of a source and its distance to the observatory in the thermal photon flux (Rutledge et al., 1999), we have investigated the sensitivity of all our results to the distances of the sources. We have used two sets of distance measurements and showed that their differences have a rather small impact on the EoS parameter estimation.

The MCMC method based on the stretch-move algorithm has been used to sample the whole parameter space (49 dimensions in our reference runs), and we found the best set of parameters reproducing the observational data. The method employed here has first been tested on the constant-$R_{\rm NS}$ approximation (Guillot et al., 2013), giving \mbox{R_{\rm NS}}=11.1\pm 0.4\hbox{,{\rm km}}, consistent with recent analyses (Guillot, 2016).

When applied with the meta-modeling for the nuclear EoS (Margueron et al., 2018a), our MCMC method permitted obtaining, for the first time, some constraints on the most determinant parameters: L_{\rm sym}=27.2\mbox{{}^{+10.9}{-5.3}}\mbox{,{\rm MeV}}, K_{\rm sym}=-59\mbox{{}^{+103}{-74}}\mbox{,{\rm MeV}} and Q_{\rm sat}=408\mbox{{}^{+735}{-430}}\mbox{,{\rm MeV}}. When considering current knowledge of nuclear physics as input (prior) for the value of $L_{\rm sym}$ (Lattimer & Lim, 2013), we find slightly better constrained parameters, as expected: L_{\rm sym}=37.2\mbox{{}^{+9.2}{-8.9}}\mbox{,{\rm MeV}}, K_{\rm sym}=-85\mbox{{}^{+82}{-70}}\mbox{,{\rm MeV}} and Q_{\rm sat}=318\mbox{{}^{+673}{-366}}\mbox{,{\rm MeV}}. We stress that the values of $K_{\rm sym}$ and $Q_{\rm sat}$ we reported are the first estimations for these empirical parameters extracted from observational data. These quantities are not yet accessible in nuclear physics experiments and are therefore poorly constrained (Margueron et al., 2018a), since their effects are mainly situated far from saturation density, such as in NS matter. We also obtained an anti-correlation between $K_{\rm sym}$ and $Q_{\rm sat}$, induced by the causality and stability requirements. The distributions of these empirical parameters are not affected by the choice in the distance set.

As a product of our analyses, we also provide the average radius (at 1.45$\,{\rm M}_{\odot}$) of the statistically preferred EoS. When the prior on $L_{\rm sym}$ is included ,We find R_{1.45}=12.42\pm 0.34\hbox{,{\rm km}} for the set of distances Dist #1 and R_{1.45}=12.35\pm 0.37\hbox{,{\rm km}} for the set of distances Dist #2. These resulting radius distributions are narrower than the range of radii allowed by the meta-model used (see Figure 5), i.e., than the prior range imposed by our choice of nuclear physics input. One can note that the radius obtained here is at the upper bound of previous analyses, (e.g., Steiner et al., 2018), resulting from the fact that we took into account the nuclear physics knowledge through the prior on $L_{\rm sym}$. Adding that prior did not degrade the fit statistics. Furthermore, the average radius was constant under this change, implying that it is required by the data, and not driven by the $L_{\rm sym}$ prior. The only nuclear physics input in our model is the well accepted condition of an EOS respecting causality and reaching at least 1.9$\,{\rm M}_{\odot}$. Leaving $L_{\rm sym}$ free results in a posterior distributions at values significantly lower than the prior, but it is compensated by an adjustment of the other two parameters, $K_{\rm sym}$ and $Q_{\rm sat}$, that keeps the radius essentially constant, while supporting a 1.9$\,{\rm M}_{\odot}$ NS. We further note that there are major differences between the meta-model and the constant-radius assumption. The latter does not require to satisfy causality and does not impose a condition on the maximum mass of NS. These conditions together naturally make the radius at 1.4$\,{\rm M}_{\odot}$ larger for the meta-model than in the constant-radius toy-model.

While previous analyses invoked the need for He atmosphere model555As noted in Section 2.1, applying He atmosphere models to NS emission spectra produces NS radii larger by $\sim$ 30–50% (Servillat et al., 2012; Catuneanu et al., 2013; Heinke et al., 2014; Steiner et al., 2018). to reconcile the otherwise small radii obtained from qLMXB spectra (e.g., Guillot et al., 2013; Guillot & Rutledge, 2014), we demonstrated that the use of our meta-model produces radii in the 12–13 km range, with or without prior on $L_{\rm sym}$.

We have also investigated the impact of the selection of the sources on the results and found that we can separate the sources in two ways: according to the S/N (group A and B presented in Table 1), or according to the posterior distribution of the mass (groups A’ and B’). When using only the sources with a high S/N, or a peaked posterior mass distribution, we found slightly smaller radii $R_{1.45}=12.2\pm 0.3$ km compared to our reference results. On the other hand, selecting the sources with lower S/N or a flat mass distribution increased the radius up to $R_{1.45}=12.9\pm 0.4$ km. These results therefore advocate for improving the statistics for the sources in $\omega$ Cen, NGC 6304, M30 and M13.

In the future, we foresee two possibilities:

The mass and radius predictions presented in this work are consistent with those obtained by other means (e.g., pulse waveform modelling with NICER, constraints obtained with future GW signals from NS-NS mergers), which would bring support to the nuclear physics assumptions we made; 2. 2.

Future analyses prefer low-radius NS, which would suggests tension with these nuclear physics assumptions. Such tension would open up the possibility to learn about dense matter, to eventually reject some of the assumptions or advocate for the presence of phase transitions in dense nuclear matter, which goes beyond our present model.

We plan to improve the nuclear EoS modeling by implementing strong first-order phase transitions and by calibrating its parameters on the data in the same way as we have done here for the empirical parameters. We believe that this will shed light on the need for first order phase transitions to reproduce the thermal spectrum of qLMXBs. The model selection could also include constraints from other observables, such as those expected from NICER as well as the wealth of new results expected from the LIGO-Virgo collaboration.

Some limitations due to flux calibrations will always remain for methods that rely on broad band X-ray spectroscopy, such as that in the present work, since F_{\rm X}\propto(\mbox{R_{\infty}}/D)^{2}. In addition to the multiplicative constants accounting for flux cross-calibration uncertainties between the instrument, we have also included 3% systematic uncertainties to each spectral bin, as was done in previous works Guillot et al. (2013); Guillot & Rutledge (2014); Bogdanov et al. (2016), to include flux calibrations uncertainties into our final results. We note, however, that at the moment, other sources of uncertainties (e.g., on the distances to the sources) likely dominate over flux calibration uncertainties.

The authors thank the anonymous referee for their useful comments that improved the discussion in this article. We acknowledge the support of ECOS-CONICYT collaboration grant C16U01. The authors are grateful to the LABEX Lyon Institute of Origins (ANR-10-LABX-0066) of the Université de Lyon for its financial support within the program ”Investissements d’Avenir” (ANR-11-IDEX-0007) of the French government, operated by the National Research Agency (ANR). NB and JM were partially supported by the IN2P3 Master Project MAC. The authors also thank the ”NewCompStar” COST Action MP1304 and PHAROS COST Action MP16214 for the conferences where this project was born. SG and NAW acknowledge the support of the French Centre National d’Études Spatiales (CNES), and of the FONDECYT Postdoctoral Project 3150428 in the early phases of this work. Additional support for MC is provided by the Chilean Ministry for Economy, Development, and Tourism’s Millennium Science Initiative through grant IC 120009, awarded to the Millennium Institute of Astrophysics (MAS). The work of MC and AR is funded by the Center for Astronomy and Associated Technologies (CATA; CONICYT project Basal AFB-170002). AR also acknowledges support from FONDECYT grant #1171421.