A new method to constrain neutron star structure from quasi-periodic oscillations

TL;DR

This paper introduces a novel method to determine neutron star properties and constrain their equation of state by fitting observed quasi-periodic oscillation frequencies with theoretical models based on relativistic geodesic motion.

Contribution

The paper presents a new approach that uses geodesic frequency calculations and Bayesian inference to extract neutron star parameters from QPO data, improving constraints on dense matter physics.

Findings

Successfully applied to three neutron star systems.

Derived constraints on neutron star masses and radii.

Method is adaptable to other QPO models.

Abstract

We develop a new method to measure neutron star parameters and derive constraints on the equation of state of dense matter by fitting the frequencies of simultaneous Quasi Periodic Oscillation modes observed in the X-ray flux of accreting neutron stars in low mass X-ray binaries. To this aim we calculate the fundamental frequencies of geodesic motion around rotating neutron stars based on an accurate general-relativistic approximation for their external spacetime. Once the fundamental frequencies are related to the observed frequencies through a QPO model, they can be fit to the data to obtain estimates of the three parameters describing the spacetime, namely the neutron star mass, angular momentum and quadrupole moment. From these parameters we derive information on the neutron star structure and equation of state. We present a proof of principle of our method applied to pairs of kHz QPO frequencies observed from three systems (4U1608-52, 4U0614+09 and 4U1728-34). We identify the kHz QPOs with the azimuthal and the periastron precession frequencies of matter orbiting the neutron star, and via our Bayesian inference technique we derive constraints on the neutrons stars' masses and radii. This method is applicable to other geodesic-frequency-based QPO models.

Tables

Table 1. Table 1 Numerical coefficients of the empirical relation ( 7 ), which provides the NS circumferential radius as a function of the dimensionless spin parameter and of the stellar quadrupole. The best fit values for the exponents 𝒩 1 , 2 subscript 𝒩 1 2 {\cal N}_{1,2} are given by 𝒩 1 = 4.12566 subscript 𝒩 1 4.12566 {\cal N}_{1}=4.12566 and 𝒩 2 = 0.996284 subscript 𝒩 2 0.996284 {\cal N}_{2}=0.996284 (Pappas, 2015 ) .

	$\ell =0$	$\ell =1$	$\ell =2$
$A_{\ell }$	0.00927584	-0.0252801	0.0497335
$B_{\ell }$	-0.358824	3.15892	-5.30171
$C_{\ell }$	2.94923	-3.20369	6.02522

Table 2. Table 2 Median and 90% values for the parameters of the sources analysed.

source		$M/M_{\odot }$		$j$		$\alpha$
4U1608-52		${2.07}_{-0.15}^{+0.25}$		${0.10}_{-0.10}^{+0.16}$		${2.56}_{-1.06}^{+0.74}$
4U0614+09		${2.10}_{-0.27}^{+0.45}$		${0.20}_{-0.20}^{+0.24}$		${2.22}_{-0.68}^{+1.04}$
4U1728-34		${2.11}_{-0.34}^{+0.47}$		${0.27}_{-0.27}^{+0.24}$		${2.00}_{-0.46}^{+1.18}$

Equations

$$d{s}^2=f^{-1}\left[e^{2\zeta }\left(d{\rho }^2+d{z}^2\right)+{\rho }^{2}d{\phi }^2\right]-f{\left(dt-\omega d\phi \right)}^2,$$

$$M_2=-\alpha j^2{M}^3,S_3=-\beta j^3{M}^4,M_4=\gamma j^4{M}^5,$$

$$y_1$$

$$y_2$$

$$\mathcal{P}(\vec{\theta }|\vec{O})\propto {\mathcal{P}}_0(\vec{\theta })\mathcal{L}(\vec{O}|\vec{\theta }),$$

$${\chi }^2=\sum _{i=1}^{N_{\text{obs}}}\left(\frac{{\Delta }_{\varphi }^2}{{\sigma }_{{\nu }_{\varphi }}^2}+\frac{{\Delta }_{\mathrm{per}}^2}{{\sigma }_{{\nu }_{\mathrm{per}}}^2}\right),$$

$$R_{\text{NS}}/M=\sum _{i=0}^3[{\mathcal{B}}_i{j}^i+{\mathcal{A}}_i{j}^i{\alpha }^{{\mathcal{N}}_1/2}+{\mathcal{C}}_i{j}^i{\alpha }^{{\mathcal{N}}_2/2}],$$

$$\mathcal{P}(M,\chi ,\alpha )=\int \mathcal{P}(\vec{\theta }|\vec{O})d{r}_1\dots d{r}_n.$$

$$f(\rho ,z)$$

$$\omega (\rho ,z)$$

$$\zeta (\rho ,z)$$

$$A(\rho ,z)$$

$$B(\rho ,z)$$

$$H(\rho ,z)$$

$$G(\rho ,z)$$

$$F(\rho ,z)$$

$$f(\rho ,z)$$

$$\zeta (\rho ,z)=\frac{1}{2}\log \left(\frac{r^2-M^2+r_{-}{r}_{+}}{2r_{-}{r}_{+}}\right)+\frac{j^2{M}^4{\rho }^2\left[(1-3\alpha ){\rho }^2+4(3\alpha -2)z^2\right]}{4r^8},$$

Full text

1. Introduction

Neutron stars (NSs), the densest stable stellar objects in the Universe, provide key information on the properties of cold, supra-nuclear density matter, strong gravitational fields and a variety of astrophysical processes that take place in or around them. Observations of NSs extend over the entire electromagnetic spectrum, from radio frequencies to gamma rays. The recent detection of gravitational waves (GWs) from two coalescing NSs has opened another, entirely different, observational window. Combinations of measurements of NS parameters such as mass, radius, moment of inertia, angular frequency, tidal deformability or quadrupole moments, which have already been obtained (or are within reach of present and future instrumentation), hold the potential to constrain the equation of state (EoS) of dense matter and the structure of NSs to high precision and accuracy (see e.g. Lattimer & Prakash 2007; Hinderer et al. 2018). To this aim different methods exploiting different diagnostics have been pursued: for instance X-ray based methods have been devised to infer the mass and radius of NSs from the evolution of radius expansion Type I X-ray bursts, broadened Fe K$\alpha$line profiles, fastest rotation frequencies, periodic modulations resulting from NS rotation (for a review see Watts et al. 2016; Lattimer 2019). NICER X-ray observations of a rotation-powered NS have recently provided $\sim 8-12$% precise mass and radius measurements, based on the latter method (Riley et al., 2019a; Miller et al., 2019a). Constraints on the EoS of NSs have recently been obtained also from the GW signal of GW170817 (Abbott et al., 2018). Present limitations in the determination of the EoS of NS arise from data paucity and/or quality (insufficient sensitivity and signal to noise ratio especially), modelling uncertainties and control of systematics.

The fast quasi-periodic oscillations (QPOs) in the X-ray flux of NS low mass X-ray binaries are among the observables that can yield measurements of NS mass and radius. These signals appear as narrow features in the power spectra of the light curves of accreting NS and black holes. While their interpretation is still debated, QPOs are believed to be produced in the inner regions of the accretion disc surrounding the compact object. They often display different modes, some which are excited at the same time, and undergo correlated frequency changes (for a review see van der Klis 2006). Much of what is currently known about QPOs derives from observations with the Rossi X-Ray Timing Explorer (RXTE, Swank et al. 1994). RXTE detected a large number of low frequency QPOs (LF QPOs, observed below several tens of Hz) and enabled the discovery of high frequency QPOs (HF QPOs), often observed in pairs, with frequencies up to over a kHz in NS systems (the upper and lower kHz QPOs).

Following the realisation that kHz QPO frequencies are close to dynamical frequencies around compact objects (van der Klis et al., 1996; Strohmayer et al., 1996), the potential of QPOs for diagnosing matter motion in strong gravitational fields became apparent. This stimulated the development of models involving the innermost regions of accretion flows where general relativistic departures from Newtonian gravity are large (see e.g.Belloni & Stella 2014). In a generic model that has been widely adopted since, the upper kHz QPO signal is directly associated to the azimuthal frequency of matter orbiting in the inner disk region, as Syunyaev (1973) and others had suggested decades earlier). More complex local models aiming at interpreting also other QPO modes, the low frequency QPOs and the lower kHz QPO in particular, were soon proposed, building on the idea that QPOs are excited at a specific radius in the disk. Among these are the Relativistic Precession Model (RPM Stella & Vietri 1998, 1999) and the epicyclic resonance model (ERM Abramowicz & Kluzniak 2001; Kluzniak & Abramowicz 2001, 2002; Fragile et al. 2016). The frequencies of the signals predicted for different QPO modes consist of combinations of the fundamental frequencies $\nu_{\phi}$, $\nu_{\theta}$ and $\nu_{r}$ of quasi-circular geodesics in the strong field regime. Also QPO models involving global oscillations of the inner disk, such as g- and p-modes (Nowak & Wagoner, 1992, 1991; Nowak et al., 1997; Wagoner et al., 2001; Kato, 2004; kato, 2012; Kato, 2012) and corrugation and warping modes (Markovic & Lamb, 1998; Armitage & Natarajan, 1999; Kato, 2001) predict frequencies related to the frequencies of geodetic motion.

As geodesics are determined by the spacetime metric, the above QPO models as well as models that extend them in different ways (see e.g. Torok et al. 2014; Stuchlík & Kološ 2016 and reference therein) hold the potential to probe the spacetime itself and measure key parameters of the compact object. In applications to black holes, the Kerr spacetime is generally used. The varying frequencies of two (and in one case three) QPO modes measured from a few black hole transients were successfully fit to the frequencies predicted by the RPM, and black hole mass and spin inferred (Motta et al., 2014a, b). Similarly the $3:2$ frequency ratio in some black hole QPOs were exploited in the ERM to infer black hole parameters (e.g. Török et al. 2005). The spacetime around rotating NSs is more complex, as it is characterised by mass, angular momentum, and higher order multipole moments. Analytical approximations to the fundamental frequencies of geodesic motion around NSs were used in early applications of the RPM (e.g. Stella & Vietri 1998). Other studies adopted instead the numerically calculated spacetimes around rotating NS models, for selected EoSs (Morsink & Stella, 1999; Stella et al., 1999).

In this paper we introduce a general method to compute geodesic motion in the strong-field spacetimes generated by NSs, which builds on the finding that higher order multipole moments of NSs are related, to a good approximation, by some“three-hair relations” (Pappas & Apostolatos, 2014; Stein et al., 2014; Yagi et al., 2014; Yagi & Yunes, 2017). In this framework, fundamental frequencies of geodesic motion can be computed for any stationary and axisymmetric spacetime once the mass, angular momentum and quadrupole moment are specified (for some additional applications, see Maselli et al., 2015a, b). These frequencies can then be used within QPO models to fit the observed QPO frequencies, test the assumptions of any given model (such as, e.g., the association of certain predicted frequencies with a given set of QPOs), and ultimately measure NS parameters that can be translated into constraints to the EoS. An explicit form of the space-time metric in terms of mass, angular momentum, and higher order multipole moments has been derived by Pappas (2017)111We note that our metric is a different, more accurate approach than the second order Hartle-Thorne spacetime which neglects multipole moments higher than the quadrupole, and is not valid for larger rotation rates.. Such a metric enables to reproduce the orbital features of rotating NSs computed in fully relativistic numerical simulations up to the innermost stable circular orbit, with better than 1% accuracy (Pappas, 2017). In this work we adopt the above metric, which enables to extend previous work on NS QPOs using analytic spacetimes (Pappas, 2012; Pappas & Apostolatos, 2013; Pappas, 2015; Tsang & Pappas, 2016; Pappas, 2017).

Our aim is to measure/constrain the NS parameters by modelling the QPOs in the X-ray light curve of accreting NS binaries, making use of an accurate description of the NS surrounding spacetime. In this paper we present a first application of a new method to do so that differs from similar methods previously proposed in the use of (i) an accurate metric to describe the spacetime of rapidly rotating NS, (ii) a Bayesian inference technique, which yields better constraints as the number of QPOs considered increases as it can take into account the information of all the QPOs consiered at the same time. We use QPOs measured with RXTE from NS LMXBs with known rotation period (from oscillations observed during type-I X-ray bursts, see, e.g., Watts, 2012). We considered the sample of van Doesburgh & van der Klis (2017), and selected three sources with a large number of QPO triplets (i.e. two kHz QPOs and a LF QPO observed simultaneously), namely 4U1608-52, 4U0614+09, and 4U1728-34.

Here we consider only the twin kHz QPOs; in a follow-up paper (Maselli et al., in prep) we will report on the results obtained by considering QPO triplets. We follow the prescriptions of, e.g., the RPM, and we identify the upper and lower kHz QPOs with the azimuthal frequency $\nu_{\phi}$, and the periastron precession frequency, $\nu_{\textnormal{per}}=\nu_{\phi}-\nu_{r}$ (Stella & Vietri 1999)222According to the RPM, a third frequency (historically known as horizontal branch oscillations, see van der Klis, 1995) is associated with the nodal precession frequency, $\nu_{\textnormal{nod}}=\nu_{\phi}-\nu_{\theta}$ (Stella & Vietri 1998; Stella et al 1999).. The two QPO signals (three, when a third QPO is considered) are assumed to be generated at the same orbital radius. Correlated QPO frequency variations are thus the result of variations in the radius at which the QPOs are emitted (see also Motta et al., 2014a).

In this paper we shall use geometric units $G=c=1$. The mass will be expressed either in kilometres, or in solar masses $M_{\odot}=1.4768$ km.

2. The metric for a rotating neutron star

The metric around a rotating NS is expressed as a stationary and axisymmetric spacetime; it can be parametrised in terms of the first five relativistic multipole moments (Geroch, 1970a, b; Hansen, 1974; Fodor et al., 1989), i.e., the mass $M$, the angular momentum $J$, the mass quadrupole $M_{2}$, the spin octupole $S_{3}$, and the mass hexadecapole $M_{4}$ (Pappas, 2017). The line element for such a spacetime can be written as (Papapetrou, 1953),

[TABLE]

where ($\rho,z$) are the Weyl-Papapetrou coordinates and the metric components $f,\;\omega,$ and $\zeta$, shown in Appendix A, are functions of the multipole moments and of the coordinates $(\rho,z)$333This metric is an approximate vacuum solution of Einstein’s field equations, that is accurate up to $M_{4}$ in the moments and up to sixth order in $M/\sqrt{\rho^{2}+z^{2}}$..

To adjust the spacetime to the stellar structure of the central object, the right set of multipole moments must be specified. Recent work (Pappas & Apostolatos, 2014; Yagi et al., 2014) has shown that for NSs the first few relativistic multipole moments can be expressed as,

[TABLE]

where $M$ is the mass and $j=J/M^{2}$ is the spin parameter, with $J$ being the angular momentum of the star. For NSs the coefficients $\alpha$, $\beta$, and $\gamma$ can be much larger than $1$, in contrast to Kerr black holes (see Doneva & Pappas (2018) for a review). Furthermore it has been shown that for realistic EoSs, based on microphysical calculations, the higher order NS multipole moments (higher than $M_{2}$) can be expressed in terms of the quadrupole, angular momentum and mass (Pappas & Apostolatos (2014); Stein et al. (2014); Yagi et al. (2014); we refer the reader to Pappas (2017) for a detailed discussion). The spin octupole and the mass hexadecapole of a NS are related to the quadrupole by the relations 444Eqs. (3), (4) are accurate up to $5\%$ for all state-of-the-art hadronic NS EoSs (Yagi et al., 2014; Pappas, 2017). In the case of quark stars, similar relations hold with slightly different values of the coefficients.

[TABLE]

where $y_{1}=\sqrt[3]{-\bar{S}_{3}}=\sqrt[3]{\beta}$, $y_{2}=\sqrt[4]{\bar{M}_{4}}=\sqrt[4]{\gamma}$, $x=\sqrt{-\bar{M}_{2}}=\sqrt{\alpha}$, and $\bar{M}_{n}=\frac{M_{n}}{j^{n}M^{n+1}}$, $\bar{S}_{n}=\frac{S_{n}}{j^{n}M^{n+1}}$ are the (dimensionless) reduced moments. For NSs $\alpha$ varies in the range between $1.5$ and $10$ for masses between $1M_{\odot}$ and up to the maximum mass (which is different for each EoS), where smaller values correspond to larger masses.

Therefore, the description of the spacetime and of the various geodesic frequencies will only depend on three parameters: the mass $M$ (that we express in units of kilometers in $G=c=1$ units), the dimensionless spin parameter $j$, and the dimensionless reduced quadrupole $\alpha\equiv-\bar{M}_{2}$. This is especially relevant for our analysis, as it reduces the number of parameters to be constrained.

By relating QPO frequencies to geodesic frequencies through the adoption of a given model (the RPM in our present case), the mass $M$, reduced spin $j$ and reduced quadrupole moment $\alpha$ can be measured as independent parameters determining the characteristics of the spacetime. These three parameters, if measured precisely, provide constraints on the NS EoS. This top-down approach involves no a priori assumption on the EoS and contrasts with other studies in which geodesic frequencies are compared to QPO frequencies based on the spacetime calculated for individual EoS and specific values of $M$ and $j$ (see e.g. Morsink & Stella 1999; Torok et al. 2016)

The spacetime considered in this work is equivalent to that introduced in (Pappas, 2017) and valid for all rotation rates (see Appendix A for technical details). Our metric is different and more accurate than slowly-rotating approaches, as the second order Hartle-Thorne spacetime (see Urbancova et al. (2019) for recent and detailed analyses using this approach). It extends previous work on NS QPOs based on analytic spacetimes (Pappas, 2012; Pappas & Apostolatos, 2013; Pappas, 2015; Tsang & Pappas, 2016; Pappas, 2017). The accuracy of the epicyclic frequencies used in this paper has been tested against fully numerical solutions showing an agreement better than 99% down to the innermost stable circular orbit (Pappas, 2017), which outperforms the approaches described above.

3. Numerical analysis

In order to test our method with the QPO frequencies from the sources in our sample we adopted a Bayesian approach. For a given set of $n$ observations $\vec{O}$ we wish to determine the posterior probability distribution of the system’s parameters $\vec{\theta}=(r_{i=1,\ldots n},M,j,\alpha)$, i.e.

[TABLE]

where ${\cal P}_{0}(\vec{\theta})$ represents the prior information on the parameters. ${\cal L}(\vec{O}|\vec{\theta})$ is the likelihood function, which we assume proportional to a chi-square variable, ${\cal L}\propto e^{-\frac{1}{2}\chi^{2}}$, given by:

[TABLE]

where $N_{\textnormal{obs}}\geq N_{\textnormal{min}}$, $\Delta_{k}\equiv\nu_{k}^{\rm obs}-\nu_{k}(\vec{\theta})$, and $\nu_{j}$ can be either the azimuthal or the periastron precession frequency.

Since we use here only pairs of kHz QPO frequencies (”doublets”), each source provides $2N_{\textnormal{obs}}$ frequencies, which are used to determine $3+N_{\textnormal{obs}}$ unknown parameters, i.e. the NS parameters $(M,j,\alpha)$ and the circumferential radii $r_{i}$ where each is produced. Thus, we need at least $N_{\textnormal{min}}=3$ doublets to characterise each source. As the errors of the QPO observed frequencies are in general asymmetric (see Table 3), i.e. $\nu_{j}=\nu_{j}^{\textnormal{obs}}\pm\sigma^{(\pm)}_{\nu_{j}}$, with $\sigma^{(+)}_{\nu_{j}}\neq\sigma^{(-)}_{\nu_{j}}$, for the sake of simplicity we compute the chi-square functions of Eq. (6) by using their $\sigma_{\nu_{j}}=[\sigma^{(+)}_{\nu_{j}}+\sigma^{(-)}_{\nu_{j}}]/2$. This has a negligible effect on inferred values of the source parameters. Moreover, given the accuracy of our metric, the relative difference between the actual frequencies and those computed using (1) is subdominant with respect to the observational errors, and therefore it will be neglected in our analysis (Pappas, 2017).

We sample the posterior distribution (5) using a Markov Chain Monte Carlo (MCMC) approach, based on the Metropolis-Hastings algorithm (Gilks et al., 1996). The random jump within the parameter space is chosen according to a multivariate Gaussian distribution, whose covariance matrix is continuously updated through a Gaussian adaptation scheme (Müller & Sbalzarini, 2010), which increases the mixing of the chains and boosts the convergence to the target distribution. For each set of data, we run four independent chains of $2\times 10^{6}$ samples, generally discarding the first $10\%$ of the simulation as a burn in. The convergence of the processes is then assessed by a standard Rubin test.

We consider flat prior distributions for all the parameters within the following ranges: $M\in[0.7,3]M_{\odot}$, $j\in[0,0.7]$, $\alpha\in[1.5,15]$, $r_{i}\in[R_{\textnormal{NS}},15M]$. Note that the prior on $\chi$ allows only co-rotating orbital motion. We also set the prior on $r_{i}$ such that it is always larger than the stellar equatorial circumferential radius $R_{\textnormal{NS}}$. The latter can be expressed with very good accuracy (a few percent) as a function of $M,j$, and $\alpha$:

[TABLE]

were the numerical coefficients $({\cal B}_{i},{\cal A}_{i},{\cal C}_{i},{\cal N}_{1},{\cal N}_{2})$ are listed in Table 1.

Once the ${\cal P}(\vec{\theta}|\vec{O})$ is sampled by the MCMC, we derive the probability distribution of the source parameters, by marginalising the joint posterior distribution over the emission radii:

[TABLE]

In the following sections we apply this analysis to the kHz QPO of 4U1608-52, 4U0614+09, 4U1728-34, and compute ${\cal P}(M,\chi,\alpha)$ using various combinations of the doublets for each source, (see list in Table 3). The latter step is a crucial to our analysis as it provides a self-consistency test for the applicability of the RPM in conjunction with the geodesic frequencies calculated in our method. If the assumptions of the method are correct, different sets of doublets for a given source must necessarily provide probability distributions for masses, spins and quadrupole moments which are consistent with each other to within the uncertainties. On the contrary, an inconsistency of the different probability distributions would signal problems with the adopted geodesic and/or QPO model.

4. Results

Since we analysed kHz QPOs doublets, i.e. the QPOs that in the RPM correspond to the periastron and the azimuthal frequencies, the MCMC requires for each source at least three sets of QPO frequencies, each containing the QPO doublet relevant for the model, i.e. $(\nu_{\phi},\nu_{\textnormal{per}})_{i=1,2,3}$, to be solved in terms of $(M,\alpha,j)$ and three emission radii $(r_{1},r_{2},r_{3})$. Therefore, for each source, we randomly selected sub-sets of three doublets from the available data. Figure 1 shows the posterior probabilities obtained for different sets of doublets, while the associated numerical values are reported in Table 2. As a representative case, we describe in detail the results we obtained for 4U1608-52.

The posterior distributions obtained by using different datasets are in good agreement with each other, and the different parameter distributions are all remarkably consistent (even though in some cases shifts between the peaks of the different marginal distributions are apparent). However the quadrupole moment of the star remains essentially unconstrained, its probability distributions being almost flat between $\alpha\simeq 1.5$ and $\alpha\simeq 4$. This is to be expected for two reasons: firstly, the quadrupole moment gives a sub-leading contribution to the spacetime metric relative to mass and spin; and secondly, the effect of the quadrupole moment is largest on the $\nu_{\textnormal{nod}}$ frequency - which we do not consider in this work - while $\nu_{\phi}$ and $\nu_{\textnormal{per}}$ are only weakly affected by its variation. However, it should be noted that adopting the approach described here with a larger set of QPO doublets would allow the quadrupole moment to be better constrained. The box plot shown in Fig. 2 supports the use of our approach together with the RPM model: it is clearly seen that the median of the distributions are all consistent with each other, and the interquartile ranges overlap with good accuracy.

Motivated by the above results, we performed the same Bayesian analysis by progressively increasing the number of doublets. We found that in all cases the posterior distributions are consistent and yield acceptable parameters. We show the results of these steps in the triangle plot in Fig. 3. The diagonal and off-diagonal panels show the marginalised and the 2D joint distributions of $M$, $j$, and $\alpha$, respectively. The mass is the parameter that we determine with the highest precision: at $90\%$ confidence level555Following (Abbott et al., 2019) for parameters bounded by the prior, as the quadrupole, we quote the $90\%$ one-sided upper limit, while for all the other parameters we quote the $90\%$ highest-posterior density interval. we find $M\in[1.92,2.32]M_{\odot}$, with a median of $M=2.07M_{\odot}$.

Our analysis also provides a poorly constrained spin parameter with median $j=0.1$ within the $90\%$ interval $j\in[0,0.26]$. The 2D distributions of Fig. 3 show some degree of positive correlation between of $j$ and $M$. The inclusion of all datasets does not modify significantly the constraints on $\alpha$, for which we still obtain a flat posterior distribution within $\alpha\in[1.50,3.30]$ at $90\%$ (see the bottom row of Fig. 3). In fact the values of the quadrupole sampled by the Monte-Carlo simulation are essentially degenerate with respect to $M$ and $j$. As noted above, this result is somewhat expected as the azimuthal and periastron precession frequencies depend weakly on $\alpha$.

The MCMC also leads to constraints on the orbital radius associated with each doublet in the analysis. We find values $r_{i}/M\lesssim 7$ in all cases. Figure 4 shows the the $90\%$ confidence intervals for the 7 values. As noted in applications of simple geodesic models (e.g. Miller et al. 1998) upper limits on the QPO radius derived the from the highest observed QPO frequencies provide an upper bound on the stellar radius. Requiring that the oscillations is generated within the accretion disc at orbital distances larger than $R_{\textnormal{NS}}$ yields the limit $R_{\textnormal{NS}}\lesssim 6.4M$ at 90% confidence level from the highest frequency Doublet of 4U1608-52 (see $r_{6}$ in Fig. 4).

The analysis of the kHz QPOs of 4U0614+09 and 4U1728-34 yields similar results to that of the 4U1608-52 kHz QPOs. Figure 5 shows that also for these systems there exists a general agreement of the inferences from different subsets of doublets, with 4U1728-34 displaying the largest spread in the posterior distributions of the parameters when the full set of doublet is taken into account. 4U0614+09 and 4U1728-34 are also characterised by a mass distribution which peaks around $\sim 2M_{\odot}$. The posteriors of $M$ are shown in Fig. 6; different groupings of doublets appear to give consistent results for each system. and converge to a similar distribution. These are fairly large mass values compared to those typically derived for, e.g., isolated pulsars (Özel et al., 2012), but we note that other studies recovered similarly large NS masses (e.g., Stella & Vietri (1999); Torok et al. (2014)). The quadrupole moments inferred for these two systems would suggest more compact objects (i.e, with larger values of $M/R$) with faster rotation rates with respect to 4U1608-52, although we stress that the inferred spin parameters are not very well constrained. The median of all parameters for each sources, together with $90\%$ uncertainties, are reported in Table 2.

The values of $(M,j,\alpha)$ provide an estimate for the NS equatorial radius through (7). We used the values we determined from the analysis of the full set of doublets of 4U1608-52, 4U0614+09, and 4U1728-34 to build the joint 2D distribution of ${\cal P}(M,R_{\textnormal{NS}})$. The contour plot of Fig. 7 shows the $90\%$ region of the mass-radius posteriors of the three sources. Fig. 7 shows also the expected mass-radius relations for three hadronic EoS (APR (Akmal et al., 1998), SLy4 (Douchin & Haensel, 2001), and UU (Wiringa et al., 1988)), are all consistent with current constraints from X-ray observations (Ozel & Freire, 2016; Riley et al., 2019b; Miller et al., 2019b; Raaijmakers et al., 2019) as well as gravitational wave data (Abbott et al., 2018). They are also consistent with the most recent results inferred by the NICER experiment from observations of a rotation-powered pulsar (Riley et al., 2019b; Miller et al., 2019b; Raaijmakers et al., 2019).

5. Discussion and Conclusions

Based on an analytic description of the spacetime around a NS in terms of three independent parameters (mass, spin, and quadrupole moment), we developed a novel method to calculate accurately the frequencies of geodesic motion in the closest vicinity of the star, without resorting to a specific EoS, but instead adopting priors to our treatment that only allow specific parameter ranges for an unspecific EoS. Once these frequencies are related to the QPO frequencies observed in low mass X-ray binaries through a QPO model, they can be fit to the data in order to obtain estimates of the NS mass, spin and quadrupole moment, from which, in turn, an estimate of the NS radius and mass can be derived. This provides constraints to the EoS, and in principle can inform models of supranuclear density matter in a complementary way to other X-ray based techniques.

We presented a proof of principle of the method in application to the observed pairs of kHz QPO frequencies (doublets) in three systems (4U1608-52, 4U0614+09 and 4U1728-34) by applying the QPO frequency identification proposed, for instance in the relativistic precession model (Stella & Vietri, 1999). Through a Bayesian analysis we obtained mass estimates around $\sim 2.05-2.16$ M*⊙*, close to (but not exceeding) the maximum of the observed distribution of NS masses (Lattimer, 2019).

Our results appear to favour stiff EoSs, i.e. large NS masses and relatively large NS radii, as indicated by the marginalised distributions (see Fig. 7). Our method yields relatively low values of the spin parameter $j$, all strictly smaller than 0.3. One may ask, how is it possible to produce any constrains on the parameter $\alpha$, since the geodesic frequencies used are such weak functions of that parameter. This is due to the fact that the spin $j$ and $\alpha$ also determine both the equatorial radius of the star and the ISCO. This means that the radius where the orbital motion occurs is limited by either the ISCO or the surface of the star, which limits the possible range of frequencies. These limitations in turn result in constraints on the parameters $j$ and $\alpha$.

Once the values of $\alpha,j$ and $M$ are known, we can derive the probability distribution of the neutron star rotation rate $f$. To this aim we use the semi-analytic fit derived in Pappas (2015) (see eq. B1 and the discussion in Pappas 2015). Due to the poor constraints on the quadrupole moment $\alpha$, the resulting bounds on $f$ are rather loose. We find that $f$ is constrained at 90% confidence interval within $\sim[0,582]$Hz ($f\sim[0,739]$Hz at 95% confidence level) and $\sim[20,1000]$Hz ($f\sim[18,1100]$Hz at 95% confidence level) for 4U1608-52 and 4U0614+09. $f$ is unconstrained for 4U1728-34, for which we have a nearly flat distribution between 0 and 1000 Hz. These values are consistent (marginally, in the case of 4U1608-52) with previous estimates of the NS spin frequencies for the above binary systems, i.e. $f\approx 620$Hz, $f\approx 415$Hz and $f\approx 363$Hz, respectively (Watts, 2012). We stress that being this work a proof of concept, we decided to adopt priors with only weak limitations in order to assess the goodness of our approach. If additional information, chiefly the NS rotation period, were included as a prior in the analysis, significanlty tighter bounds would result. A precise determination of the NS mass, moment of inertia and quadrupole moment that can be obtained in this way would provide an unprecedented three-parameter constraint on the EoS, for which suitable 3D mappings of the EoS (as opposed to the familiar 2D mapping involving mass and radius, or mass and moment of inertia, Lattimer & Prakash 2007) would be needed. Moreover, such measurements could be used together with current and future constraints form gravitational wave sources, to infer multi-messenger bounds on the stellar structure (Fasano et al., 2019).

The above mass and spin estimates are compatible with those from the early applications of the RPM, as the relevant geodesic frequencies depend only weakly on the spin and quadrupole moment (see, e.g., Stella & Vietri 1998, but also du Buisson et al. 2019). Instead, the EoS-independent bounds on the radius and mass that we derived represent original results of our new method. We note that the $\sim 2$ M*⊙* region in the NS mass-radius diagram is virtually unconstrained at present by observations in the radio and X-ray bands, while the limits on the tidal deformability from the GW170817 event, once translated into mass-radius bounds, are compatible with our results (Abbott et al., 2018). We remark that a precise measurement of the mass quadrupole moment can be used to directly constrain the EoS, without an explicit determination of the radius (see Pappas, 2012; Pappas & Apostolatos, 2014). The limits on the quadrupole that we obtained in this work are very loose due to the low number of QPOs considered (up to 8 for a given source) and no stringent conclusion can be drawn. Instead, significantly tighter bounds are obtained by considering a larger number of doublets (Motta et al. in prep).

The method we presented is amenable to further, more extensive applications which exploit different QPO datasets and/or alternative QPO models. In a forthcoming study we will fit QPO triplets, each consisting of a low-frequency QPO simultaneous to the two kHz QPOs, where the former is associated with nodal precession frequency, in addition to the azimuthal and periastron precession frequencies that were used in the present study. This is expected to yield higher precision estimates of the parameters governing the NS spacetime, since the nodal precession frequency depend quite strongly on NS spin and quadrupole moment. Our method is also applicable to other models relating QPOs to geodesic frequencies. This includes models in which QPO frequency identifications are different than those of the RPM, such as the epicyclic resonance models and global disk oscillation models in their various versions (e.g. Torok et al. (2014); Stuchlík & Kološ (2016) and refs therein). In all cases a key requirement is that the set of independent QPO frequencies predicted by a model and fit to the datais sufficiently large that the three parameters describing NS spacetime can be derived.

The reliability of results from our new method, like that of inferences based on QPOs in general, depend crucially on the correctness of the association of QPO signals to geodesic frequencies. Extensive, high sensitivity and large signal to noise QPO measurements to be obtained with next generation, large area X-ray missions, such as eXTP (De Rosa et al., 2019), Athena (Barcons et al., 2015) and STROBE-X (Ray et al., 2019) will allow the detection and precise characterisation of a significantly higher number of QPOs from many sources. In conjunction with detailed applications of competing models and advanced methods like our own, they may provide the key to resolving long-standing ambiguities in QPO interpretation and placing stringent constraints on NS structure and EoS.

Acknowledgments.— We thank Tiziano Abdelsalhin for having carefully read this manuscript. AM, GP, and PP acknowledge financial support provided under the European Union’s H2020 ERC, Starting Grant agreement no. DarkGRA–757480. AM acknowledges partial financial support provided under the European Union’s H2020 ERC Consolidator Grant “Matter and strong-field gravity: New frontiers in Einstein’s theory” grant agreement no. MaGRaTh–646597 during the development of this project. LS acknowledges financial contributions from ASI-INAF agreements 2017-14-H.O and I/037/12/0 and from “iPeska” research grant (P.I. Andrea Possenti) funded under the INAF call PRIN-SKA/CTA (resolution 70/2016). This project has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No 690904. SEM acknowledges the Violette and Samuel Glasstone Research Fellowship programme, and the Oxford Centre for Astrophysical Surveys, which is funded through generous support from the Hintze Family Charitable Foundation. The authors would like to acknowledge networking support by the COST Action CA16104 and support from the Amaldi Research Center funded by the MIUR program “Dipartimento di Eccellenza” (CUP: B81I18001170001).

Appendix A Metric components of rotating neutron stars

In this appendix we present the explicit form of the metric functions introduced in Sec. 2 for the line element given in Eq. (1), which describes the spacetime around a rotating NS. We refer the reader to Pappas (2017) for further details. The components of the metric are functions of the coordinates $(\rho,z)$, and of the multipole moments, and are given by the following expressions:

[TABLE]

where

[TABLE]

The spacetime defined above can be given in an even more convenient form so as to have the right Schwarzschild limit when the rotation goes to zero, i.e., $j\rightarrow 0$. To this aim, we resum the expansions of $f(\rho,z)$ and $\zeta(\rho,z)$, using the variable $r=\sqrt{\rho^{2}+z^{2}}$, such that when the rotation vanishes the metric coincides with the exact Schwarzschild solution in its standard form. With this procedure, we obtain:

[TABLE]

and

[TABLE]

where $r_{\pm}=\sqrt{(M\pm z)^{2}+\rho^{2}}$ and we have used the definitions (2) for the moments. When $j\rightarrow 0$, from Eq. (A2) it follows that $\omega\rightarrow 0$ and the functions $f$ and $\zeta$ take their Schwarzschild form. This metric, as the previous one, is accurate up to $M_{4}$ in the moments and up to order $\mathcal{O}(M^{6}/r^{6})$ with respect to the vacuum field equations.

It is worth remarking that while the spacetime is given in Weyl-Papapetrou coordinates, which differ from the usual Schwarzschild-like or quasi-isotropic ones, the various geodesic frequencies are coordinate-independent quantities, while the relevant radii on the equatorial plane can be expressed in terms of the circumferential radius which is also a geometric and coordinate-independent quantity.

Appendix B QPO frequencies

In this appendix we give the QPO frequencies and corresponding uncertainties for the three LMXB systems considered in this paper. Frequencies are all from RXTE/PCA observations and were taken from van Doesburgh & van der Klis (2017).