Detection of Simultaneous Quasiperiodic Oscillation Triplets in 4U 1728-34 and Constraining the Neutron Star Mass and Moment of Inertia

We report simultaneous detection of twin kHz and ∼40 Hz quasiperiodic oscillations (QPOs) in the time-resolved analysis of the AstroSat/LAXPC observation of the neutron star low-mass X-ray binary, 4U 1728-34. The frequencies of the multiple sets of triplets are correlated with each other and are consistent with their identification as the orbital, periastron precession, and twice the nodal precession frequencies. The observed relations, along with the known spin of the neutron star, put constraints on the mass and the ratio of the moment of inertia to the mass of the neutron star to be M⊙*=1.92±0.01 and I45/M⊙*=1.07±0.01 under the simplistic assumption that the metric is a Kerr one. We crudely estimate that the mass and moment of inertia values obtained may differ by about 1% and 5%, respectively, if a self-consistent metric is invoked. Using the Tolman–Oppenheimer–Volkoff equations for computing the moment of inertia of a neutron star in slow rotation approximation, having different equations of state, we find that the predicted values of neutron star parameters favor stiffer equations of state. We expect more stringent constraints would be obtained using a more detailed treatment, where the equation of state–dependent metric is used to compute the expected frequencies rather than the Kerr metric used here. The results provide insight into both the nature of these QPOs and the neutron star interior.


Introduction
The power density spectra (PDS) obtained from the X-ray light curves of X-ray binaries (XRBs) often reveal narrow features that are called quasiperiodic oscillations (QPOs; van der Klis 1989Klis , 2000;;Ingram & Motta 2019).QPOs, with frequencies ranging from a few mHz to around 1.2 kHz, have been observed in different XRBs, which are broadly classified into mHz QPOs, low-frequency (LF) QPOs (0.1-100 Hz), and high-frequency (HF) QPOs above 300 Hz.The HF QPOs provide a useful probe into the inner accretion flows around the compact objects.The HF QPOs as well as the LF QPOs are observed in both black hole (BH) and neutron star (NS) lowmass X-ray binaries (LMXBs).Understanding the HF QPOs is expected to provide insights into the behavior of matter under a strong gravitational field.For NS systems, the HF QPOs further provide the potential to constrain both mass and radius and, hence, information about the equation of state (EOS) of the NS.The HF QPOs in NS-LMXBs are often observed as pairs, which are called lower and upper kHz QPOs (van der Klis et al. 1997).
NS-LMXBs are subdivided into Z sources and atoll sources on the basis of their evolution on the color-color diagram (Hasinger & van der Klis 1989).The first kHz QPOs at ∼800 and ∼1100 Hz were discovered in a Z source Sco-X1 from the observations made by the Rossi X-ray Timing Explorer in 1996 (Berger et al. 1996;van der Klis et al. 1996).Apart from this discovery, van der Klis et al. (1996) also reported correlations between the frequencies of the kHz and LF QPOs.Similar kHz QPO features have also been reported for atoll sources, and among these atoll sources, the kHz QPOs were first discovered in 4U 1728-34 (Strohmayer et al. 1996).In 4U 1728-34, Strohmayer et al. (1996) reported three pairs of kHz QPOs in which the lower kHz QPO frequency increased from 637.5 ± 3.6 Hz to 716.0 ± 1.1 Hz, whereas the upper kHz QPO frequency increased from 988 ± 5.9 Hz to 1058.3 ± 12.1 Hz.In all three PDS, the difference between the lower and upper kHz QPOs frequencies was nearly equal to the burst oscillation (BO; ∼363 Hz) observed in the same observation, which is identified as the spin frequency of the NS.This led to the beat frequency model where the upper kHz QPO was identified as a Keplerian frequency, while the lower one was the beat between the orbital frequency and spin frequency of the NS (Alpar & Shaham 1985;Lamb et al. 1985;Wijnands et al. 2003).However, more observations were in contradiction since the frequency difference turned out to be sometimes significantly smaller than the NS spin frequency (Strohmayer et al. 1996;Miller et al. 1998;Psaltis et al. 1998;Salvo et al. 2001).Psaltis et al. (1998) showed that the upper and lower kHz QPO frequencies are correlated as a power-law for Sco-X1, and a similar correlation was obtained when data from nine other sources were included.Moreover, correlations between the QPO frequencies and characteristic frequencies of the broadband noise were discovered for both NS and BH systems (Psaltis et al. 1999).Thus, it was established that the Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
frequencies of different QPOs and the broadband noise were probably related to an underlying common process.
These correlations led to the relativistic precession model (RPM), based on the motion of particles around a spinning BH (Stella & Vietri 1997;Stella et al. 1999).In general relativity, the azimuthal frequency of a particle moving in a circular orbit around a Kerr BH of mass M and specific angular momentum ( ˜) = a J M (where J is the angular momentum), in an equatorial plane, is given by (Bardeen et al. 1972) If such a particle's motion is perturbed along the polar angle θ and radial distance r, then it would start oscillating in the directions of θ and r.The frequency of these small oscillations within the plane of the orbit (ν r ) and perpendicular to the plane (ν θ ) is as follows (Okazaki et al. 1987;Kato 1990): where G and c have been taken to be unity and the positive and negative signs are for prograde and retrograde orbits, respectively.These three coordinate frequencies give rise to two kinds of relativistic precessions: periastron and nodal precessions, and the corresponding precession frequencies are defined as ν per = ν f − ν r and ν nod = ν f − ν θ , respectively.According to the RPM, ν f , ν per , and ν nod are the frequencies of upper kHz, lower kHz, and LF QPOs, respectively.An important aspect of the model is that the mass and spin parameter (a = Jc/GM 2 ) of the compact object can, in principle, be estimated using the correlations between the observed frequencies.
There have been several indications that the RPM is applicable to the QPO frequencies observed in NS systems (Ford & van der Klis 1998;Stella et al. 1999).Stella et al. (1999) showed that correlations reported by Psaltis et al. (1999), where a large number of QPO frequencies from different systems were used, broadly follow the predictions of the model.The model has also been extended to BH systems.For two BH systems, the mass and spin parameter of the BHs have been estimated using observation of a QPO triplet for each case (Motta et al. 2014a(Motta et al. , 2022)).In the absence of simultaneous observation of a QPO triplet, the correlation between the frequencies of two QPOs, consisting of an LF QPO and an HF QPO (or a broad-noise), has been used to estimate the BH spin parameter, using the BH mass estimated through optical measurements (Motta et al. 2014b;Bhargava et al. 2021).Detection of several simultaneous QPO triplets for a single source would allow us to verify the RPM and estimate the mass and spin parameter of the NS.If the spin is known independently, this will allow for an estimate of the star's moment of inertia and subsequently constrain the NS's EOS.
In this work, we report the detection of three simultaneous QPOs in the power spectra of 13 time segments as observed by AstroSat/LAXPC for the NS source 4U 1728-34.The data are unique because all 13 segments contain QPO triplets corresponding to a single observation, showing a significant variation in the QPO frequencies.We use the data to verify the RPM and constrain the mass and moment of inertia of the NS, which will eventually help to constrain the EOS of a rotating NS.We also numerically compute the moment of inertia of NS spinning at 300 Hz for ten EOSs.We take the EOS tables that are freely available and explore some by using them in our stellar models.The stellar model is constructed by solving the Tolman-Oppenheimer-Volkoff (TOV) equation with a numerical integration algorithm.This gives the stellar model of a single star with a given core density and given EOS.The process is repeated for all the EOSs discussed in this work.To get the plot of the moment of inertia versus mass, we also take into account the constraints from the tidal deformability parameter and tidal love number.

Observation and Data Analysis
Large Area X-ray Proportional Counter (LAXPC) is an instrument on board AstroSat satellite launched by the Indian Space Research Organization on 2015 September 28.LAXPC consists of three coaligned proportional counter units named LAXPC10, LAXPC20, and LAXPC30 with a large effective area of about 6000 cm 2 , making it a perfect instrument to study even the faint X-ray sources.Each LAXPC detector records the arrival time of each photon with a time resolution of 10 μs.In event analysis mode, each LAXPC instrument works in the energy range of 3-80 keV with an energy resolution of about 15% at 30 keV (Yadav et al. 2016;Antia et al. 2017).
In the course of its performance verification phase, AstroSat/ LAXPC observed the source 4U 1728-34 covering 20 orbits during 2016 March 7-8 (observation ID: T01_041T01_ 9000000362). 6For one of these orbits (orbit number 2398), Chauhan et al. (2017) reported a kHz QPO at ∼800 Hz, which showed evolution with time and BO at ∼363 Hz, and performed spectral analysis of the persistent emission.In this work, we analyze the full data set.We used the LAXPCsoftware7 to compute the background-subtracted light curves and power spectra.All three LAXPC units were used for the analysis.There are four thermonuclear bursts detected during this observations that were removed for this analysis.We defer a study of these bursts to a future work.Figure 1 shows the 3-30 keV background-subtracted light curve of persistent emission for the whole observation, while Figure 2 shows the corresponding PDS.The software subtracts the expected dead-time corrected Poisson noise from the PDS.However, there was a small residual power at high frequencies, indicating that there is still some residual Poisson noise contribution.Hence, we included a constant represented by a power-law with zero index to take into account this excess contribution.Six Lorentzian components were required to fit the time-averaged power spectrum and their parameters, and the fit statistics are listed in Table 1.At the lowest frequencies, there is a broad feature that is modeled as a Lorentzian with centroid 10 Hz (fixed) and width ∼33 Hz.A prominent LF QPO is seen at ∼41 Hz, along with a faint broad component at ∼160 Hz.In the kHz region, complex shapes are seen, which are modeled by three relatively narrow Lorentzians.The complex shape, especially near 800 Hz, is due to varying QPO frequencies, and hence, the data have to be analyzed in segments to detect such variations.
To investigate the variation of the centroid frequencies of QPOs with time, we divided the observation into sixteen segments, as shown in Figure 1.The segments were chosen such that, for the generated PDS, the HF features could be represented by either one or two Lorenztian(s) for the upper and lower kHz QPOs.Each individual PDS was fitted with a similar model used for the time-averaged one, except that the width of the 10 Hz Lorentzian was not constrained and hence was fixed at 33.36 Hz, corresponding to the time-averaged value.The broadband noise at ∼160 Hz generally had a low Qfactor (∼1.4) and did not appear in some of the segments.In this work, we concentrate on three distinct QPOs, namely the LF QPO and the lower and upper kHz QPOs.Representative PDS of segments 6 and 7 and the fitting are shown in Figures 3  and 4, and the best-fit parameters for these three QPOs are listed in Table 2.For three of the segments (1, 2, and 5), there was no detectable upper kHz QPO, and hence, there were 13 segments for which a triplet of QPOs was detected.

Results
As mentioned earlier, following the RPM, ν f , ν per , and ν nod are identified with the frequencies of the upper kHz, lower kHz, and LF QPOs, respectively.It is convenient to express ν per and ν nod in terms of ν f by eliminating the radius r from Equations (1)-(3).Introducing the constants G and c into the expressions and noting that the specific angular momentum can be written as 2 s , where ν s is the spin frequency, one can write r as a function of ν f : where M and I/M are mass and ratio of moment of inertia to the mass, respectively.The expressions for ν per and ν nod to have M and I/M as free parameters can be written as While fitting the observed QPO frequencies with the RPM, we have taken I in units of 10 45 g cm 2 (I 45 ) and M in solar mass (  M*).Moreover, we have considered the motion of the particle in a prograde orbit.These expressions can then be compared with the observed correlations between the frequencies of the lower kHz and LF QPOs with respect to the upper kHz QPO frequency.* is significantly larger than that expected for NSs based on theoretical studies on their EOS.One possibility is that the LF QPO frequency may be equal to twice the nodal precession frequency.Indeed, some collective observational data points for a few atolltype sources (Morsink & Stella 1999;Stella et al. 1999) suggest that this is the case, and Stella et al. (1999) have argued that the geometry of the tilted inner accretion disk might produce a stronger signal at the second harmonic of the nodal precession frequency.Thus, we also considered the situation where the LF QPO frequency is twice the nodal precession frequency, which  resulted in 1.07 0.01 45 * (equivalently a = 0.145 ± 0.001).Figure 7 shows the plot of the LF QPO, lower kHz QPO, and upper kHz QPO frequencies as a function of the LF QPO frequency.Solid lines represent orbital, periastron precession, and twice the nodal precession frequencies corresponding to the best-fit values of  = M 1.92 * and a = 0.145.The vertical dashed line indicates the maximum QPO frequencies at r ISCO = 5.52 r g .The figure may be compared to the ones shown in the literature for BH systems (Motta et al. 2014a(Motta et al. , 2022)), where the theoretical curves were plotted against a single triplet QPO observation.
While using Equations (1)-(3) we have assumed the metric to be a Kerr one, while the actual metric around a spinning NS will be different and will depend on its EOS (Morsink & Stella 1999).Since such a metric has to be generally obtained numerically, we defer a detailed analysis to a later work.Here, we attempt to estimate the effect of the Kerr metric assumption by ignoring the second-order ã2 terms in the equation and fitting the correlations to obtain

64.19/63
Note.The errors quoted here are within 1σ confidence level.
Thus, considering these systematic effects, we estimate the approximate values to be  =  M 1.92 0.03 * and  = I M 45 *  1.07 0.05.The systematics estimated by ignoring the second order of spin is a crude estimation.In fact, the metric around an NS for a given EOS would depend on a number of parameters such as stellar oblateness, higher-order mass moments, magnetic field, etc.Hence, there may exist a significant degree of uncertainty associated with these systematics.

Constraining EOS
In this section, we make an approximate estimate of the moment of inertia of an NS, given an EOS, i.e., pressure as a function of density P(ρ).In particular, we numerically solve TOV equations for an EOS to obtain pressure as a function of radius P(r).We note that the TOV equations assume that the underlying metric is a Schwarzschild one.The moment of inertia is estimated numerically using Equation (11) mentioned in Ravenhall & Pethick (1994), given by Here it is assumed that the NS is slowly rotating corresponding to a spin of ∼300 Hz.
Figure 8 shows the results of the calculations with solid lines representing the variation of .The gravitational waveform for the event GW170817 leads to constraints on the dimensionless tidal deformability number, Λ (Abbott et al. 2018).The cyan (blue) colored points in Figure 8 represent the 50% (90%) confidence levels obtained from these Λ constraints, which we have obtained by solving the relevant equations (Thorne & Campolattaro 1967;Hinderer 2008;Hinderer et al. 2010;Chatziioannou 2020).The constraints on M and I/M obtained from this work are marked in Figure 8 with a red point with solid and dotted bars representing the statistical and systematic errors, respectively.Figure 9 shows a zoomed-in version of Figure 8, which reveals the location of the parameter space favored by the RPM fitting to the AstroSat data.

Summary and Discussion
To study kHz QPOs, we have thoroughly analyzed 20 orbits data of 4U 1728-34 obtained from the observation made by AstroSat/LAXPC during 2016 March 7-8.After removing thermonuclear bursts, we extracted 3-30 keV PDS of persistent emission for the whole observation (Figure 1).For PDS extraction, we used all layers of three LAXPCs and kept the frequency bin at 2 Hz.We have fitted the PDS for the entire data set with six Lorentzians and one power law with zero spectral index to account for some Poisson noise residuals at higher frequencies.The time-averaged power spectrum clearly shows the LF QPO at ∼41 Hz and the twin kHz QPOs at ∼800 and ∼1100 Hz, along with a broad feature at ∼160 Hz (Figure 2 and Table 1).A complex feature near 800 Hz showing a double-peaked lower kHz QPO was detected in the PDS.This complex feature disappeared once the data were divided into segments, which was indicative of significant variation in the lower kHz QPO for the whole observation.
To see the time evolution of QPOs, we divided the entire data set into 16 segments.For each segment, we computed 3-30 keV PDS with a 2 Hz frequency bin.The QPO frequencies were found to vary with segment number, manifesting a dynamic nature of QPOs.The LF QPO and lower and upper kHz QPO frequencies varied from ∼35 Hz to 45 Hz, ∼674 Hz to 853 Hz, and ∼1036 Hz to 1136 Hz, respectively (Table 2).For 3 out of 16 segments, upper kHz QPO was not detected.Thus, we could get 13 QPO triplets, making it the first-ever detection of multiple sets of QPO triplets from a single observation. .We used the RPM to fit the observed QPO frequencies although the metric around a spinning NS deviates from that of the Kerr metric and

Figure 1 .
Figure1.3-30 keV background-subtracted light curve of persistent emission after removing thermonuclear bursts.The light curve has been extracted from all layers of LAXPC10, LAXPC20, and LAXPC30 with a time bin of 1 s.The vertical red lines indicate the segments used to extract the PDS.S1, S2, S3, . . .S16 stand for segments 1, 2, 3, . . .16, respectively.

Figure 2 .
Figure 2. 3-30 keV PDS of persistent emission in the frequency range from 10 to 1500 Hz.It shows three clear QPOs at ∼40, ∼800, and ∼1100 Hz, along with broadband noise at ∼160 Hz.The Lower kHz QPO shows double features at ∼723 and ∼806 Hz.

Figure 3 .
Figure 3. 3-30 keV PDS of segment 6 in the frequency range from 10 to 1500 Hz.

Figure 4 .
Figure 4. 3-30 keV PDS of segment 7 in the frequency range from 10 to 1500 Hz.

Figure 5 .
Figure 5. Correlation between LF QPO and upper kHz QPO fitted with the RPM.The best-fit parameter values obtained are a = 0.269 ± 0.01 (equivalently  =  I M 2.21 0.02 45 * ) and 2.12 (fixed) with c = 0.57 red 2 in the parameters obtained, which is <1% for  M* and about 5% for  I M 45 * , may reflect the deviation caused by assuming a Kerr metric.

Figure 6 .
Figure 6.Correlation between the lower and upper kHz QPOs.The best-fit parameter values obtained are a = 0.269 (fixed; equivalently  = I M 2.21 45 * ) and 2.12 ± 0.01 with c = 1.06 red 2

Figure 7 .
Figure 7.The LF QPO and the lower and upper kHz QPOs are plotted as a function of LF QPO.The solid lines represent characteristic frequencies of the RPM, corresponding to  = M 1.92 * and  = I M 1.07 45 *.The orbital frequency, the periastron precession frequency, and twice the nodal precession frequency are represented by solid black, red, and magenta lines, respectively.The vertical dashed line depicts the maximum QPO frequencies at the innermost stable circular orbit.

Figure 8 .
Figure 8.In this plot, the solid lines represent the variation of moment of inertia (I) with the mass of the NS (M) for different equations of state.The cyan-colored dots represent the GW170817 constraints for the 50% confidence interval, and blue dots represent the GW170817 constraints for the 90% confidence interval.The red dot represents the observed values of I and M from the QPO correlation with the RPM, and the red error bar (solid/dotted) denotes the statistical/systematic error in measurement.The equations of state considered here are as follows: ALF2 (Alford et al. 2005), AP3 (Akmal et al. 1998), AP4 (Akmal et al. 1998), ENG (Engvik et al. 1995), MPA1 (Müther et al. 1987), MS1 (Mueller & Serot 1996), MS1B (Mueller & Serot 1996), SLY (Douchin & Haensel 2001), WFF1 (Wiringa et al. 1988), and WFF2 (Wiringa et al. 1988).

Figure 9 .
Figure 9.The above plot is a zoomed-in version of Figure 8. Red cross in the plot indicates the location of the parameter space constrained by this work.

Table 1
The Best-fit Parameters of 3-30 keV PDS of Persistent Emission in 4U 1728-34 for Entire Observation

Table 2
The Best-fit Parameters of the Lower and Upper kHz QPOs Observed in 4U 1728-34 for Each Segment