Atmospheric oscillations provide simultaneous measurement of neutron star mass and radius

TL;DR

This paper introduces a novel method to measure neutron star mass and radius by analyzing oscillations of levitating atmospheres during X-ray bursts, providing a new tool for stellar parameter estimation.

Contribution

The paper derives relativistic oscillatory eigenmodes of neutron star atmospheres, including damping effects, and shows how frequency maxima can determine stellar mass and radius.

Findings

Identified characteristic maxima in damped oscillation frequencies.

Established that frequency ratios depend on the adiabatic index and luminosity.

Proposed a method to estimate neutron star parameters from observed oscillation frequencies.

Abstract

Neutron stars with near-Eddington observable luminosities were shown to harbor levitating atmospheres, suspended above their surface. We report a new method to simultaneously measure the mass and radius of a neutron star based on oscillations of such atmospheres. In this paper, we present an analytic derivation of a family of relativistic, oscillatory, spherically symmetric eigenmodes of the optically and geometrically thin levitating atmospheres, including the damping effects induced by the radiation drag. We discover characteristic maxima in the frequencies of the damped oscillations and show that using the frequency maxima, one can estimate mass and radius of the neutron star, given the observed frequency and the corresponding luminosity of the star during the X-ray burst. Thus, our model provides a new way to probe the stellar parameters. We also show that the ratio of any two undamped eigenfrequencies depends only on the adiabatic index of the atmosphere, while for the damped eigenfrequencies, this ratio varies with the luminosity. The damping coefficient is independent of the mode number of the oscillations. Signatures of these atmospheres' dynamics will be reflected in the source's X-ray light curves.

Tables

Table 1. Table 1 : Frequency maximum and luminosity are assumed to be normally distributed with the mean values as given in the first and second columns, respectively. Corresponding 1 ​ σ 1 𝜎 1\sigma errors, which are given in the third and fourth columns, are determined from the probability distributions of mass and radius.

${\overline{\nu }}_{\max }$ (Hz)	$\overline{L}$ ($\times {10}^{38}$ $\mathrm{erg}{s}^{-1}$)	$M$ ($M_{\odot }$)	$r_{\ast }$ (km)
300	1.8	${1.48}_{-0.32}^{+0.32}$	${19.3}_{-2.2}^{+2.0}$
300	2.0	${1.65}_{-0.35}^{+0.37}$	${20.4}_{-2.3}^{+2.1}$
600	1.8	${1.58}_{-0.32}^{+0.35}$	${14.1}_{-1.6}^{+1.5}$
600	2.0	${1.79}_{-0.40}^{+0.42}$	${15.0}_{-1.8}^{+1.7}$
900	1.8	${1.70}_{-0.37}^{+0.41}$	${11.9}_{-1.4}^{+1.3}$
900	2.0	${1.94}_{-0.44}^{+0.48}$	${12.6}_{-1.5}^{+1.4}$

Equations

$$L(r)=L_{\infty }{\left(1-\frac{2GM}{r{c}^2}\right)}^{-1}=L_{\ast }\left(1-\frac{2GM}{r_{\ast }{c}^2}\right){\left(1-\frac{2GM}{r{c}^2}\right)}^{-1},$$

$$L_{\mathrm{Edd}}{\left(1-\frac{2GM}{r_0{c}^2}\right)}^{-1/2}=L(r_0)=L_{\infty }{\left(1-\frac{2GM}{r_0{c}^2}\right)}^{-1},$$

$$\lambda =L_{\infty }/L_{\mathrm{Edd}}$$

$$\lambda ={\left(1-\frac{2GM}{r_0{c}^2}\right)}^{1/2}$$

$$r_0=\frac{2GM}{c^2(1-{\lambda }^2)}.$$

$$-\frac{GM}{r^2}\left[{\left(1-\frac{2GM}{r{c}^2}\right)}^{-1/2}-{\left(1-\frac{2GM}{r{c}^2}\right)}^{-1}{\left(1-\frac{2GM}{r_0{c}^2}\right)}^{1/2}\right]$$

$$\approx -{\left(\frac{GM}{r_0^{2}c}\right)}^2(1+z{)}^2\delta s\equiv -{\omega }^2\delta s.$$

$${\omega }_r=\frac{GM}{r_0^{2}c}=\frac{c^3(1-{\lambda }^2{)}^2}{4GM}.$$

$${\omega }_r=\frac{c^3}{GM}{\left(\frac{GM}{r_0{c}^2}\right)}^2,$$

$$d{s}^2=g_{ij}d{x}^{i}d{x}^j=-Bd{t}^2+B^{-1}d{r}^2+r^2(d{\theta }^2+{\sin }^2\theta d{\varphi }^2),$$

$$B\equiv 1-\frac{2M}{r}=-g_{tt}=g^{rr}.$$

$${\nabla }_{\mu }\left(\rho u^{\mu }\right)=0,$$

$${\nabla }_{\mu }{T}^{\mu \nu }=G^{\nu },$$

$$R^{tt}$$

$$R^{tr}$$

$$R^{rr}$$

$$\sin \alpha =\frac{r_{\ast }}{r}{\left(1-\frac{2M}{r}\right)}^{1/2}{\left(1-\frac{2M}{r_{\ast }}\right)}^{-1/2}.$$

$$u^{\mu }{\nabla }_{\mu }{u}^r+\left(g^{\mu r}+u^{\mu }{u}^r\right)\frac{{\nabla }_{\mu }p}{\rho }+\kappa u_{\mu }\left(R^{r\mu }+u^r{u}_{\nu }{R}^{\mu \nu }\right)=0.$$

$$\frac{1}{\rho }\frac{dp}{dr}=-\frac{M}{r^{2}B}\left(1-\frac{\lambda }{B^{1/2}}\right).$$

$$r_0=\frac{2M}{1-{\lambda }^2}.$$

$$\frac{1}{\rho }\frac{dp}{dr}=-\frac{{\mathcal{E}}^2}{2}\frac{d{\mathcal{U}}_{\mathrm{eff}}}{dr}.$$

$${\mathcal{U}}_{\mathrm{eff}}=-B^{-1}+\frac{2\lambda }{3}{B}^{-3/2}.$$

$$\frac{1}{\rho }\frac{dp}{dr}+\frac{{\mathcal{E}}_0^2}{2}\frac{d{\mathcal{U}}_{\mathrm{eff}}}{dr}=-\frac{1}{2}\left({\mathcal{E}}^2-{\mathcal{E}}_0^2\right)\frac{d{\mathcal{U}}_{\mathrm{eff}}}{dr},$$

$$(1+n)\frac{p}{\rho }+\frac{1}{2}{\mathcal{E}}_0^2{\mathcal{U}}_{\mathrm{eff}}+\psi =\text{const}.$$

$$\frac{p}{\rho }=\frac{p_0}{{\rho }_0}\left\{1-\frac{1}{n{c}_{\mathrm{s},0}^2}\left[\frac{{\mathcal{E}}_0^2}{2}\left({\mathcal{U}}_{\mathrm{eff}}-{\mathcal{U}}_{\mathrm{eff},0}\right)+\psi -{\psi }_0\right]\right\},$$

$${\frac{d{\mathcal{U}}_{\mathrm{eff}}}{dr}|}_{r_0}={\frac{d\psi }{dr}|}_{r_0}={\frac{d^2\psi }{d{r}^2}|}_{r_0}=0.$$

$$\frac{p}{\rho }\approx \frac{p_0}{{\rho }_0}\left[1-\frac{1}{n{c}_{\mathrm{s},0}^2}\frac{{\mathcal{E}}_0^2}{2}\frac{d^2{\mathcal{U}}_{\mathrm{eff}}}{d{r}^2}\frac{(r-r_0{)}^2}{2}\right]\equiv \frac{p_0}{{\rho }_0}f.$$

$$\frac{-g_{tt}(r_0){\mathcal{E}}_0^2}{2g_{rr}(r_0)}{\frac{d^2{\mathcal{U}}_{\mathrm{eff}}}{d{r}^2}|}_{r_0}=\frac{{\lambda }^6}{2}{\frac{d^2{\mathcal{U}}_{\mathrm{eff}}}{d{r}^2}|}_{r_0}={\left(\frac{M}{r_0^2}\right)}^2,$$

$${\omega }_r^2$$

$$\beta$$

$$x$$

$$\eta$$

$$\frac{p}{\rho }\approx \frac{p_0}{{\rho }_0}\left(1-\frac{{\omega }_r^2}{{\beta }^2}{x}^2\right)=\frac{p_0}{{\rho }_0}(1-{\eta }^2)=\frac{p_0}{{\rho }_0}f.$$

Full text

††pubyear: 2018††pagerange: Atmospheric oscillations provide simultaneous measurement of neutron star mass and radius–Atmospheric oscillations provide simultaneous measurement of neutron star mass and radius

1 Introduction

Neutron stars are the most compact non-singular objects observed in the Universe. The theoretical densities of these objects range from about $10^{4}~{}\rm{g}~{}\rm{cm}^{-3}$ at the surface of the crust (Chamel & Haensel, 2008) to above nuclear density of around $10^{14}~{}\rm{g}~{}\rm{cm}^{-3}$ in the core (Rhoades & Ruffini, 1974). The physics governing the equation of state at such large densities becomes uncertain, and a wide range of equations of state has been proposed, each producing different mass–radius relations, for neutron stars (Arnett & Bowers, 1977; Cooperstein, 1988; Cook et al., 1994; Lattimer, 2012). Accurate measurements of both the mass and radius of neutron stars are necessary for constraining the equation of state, and thus the physics of highly dense material. Precise measurements of the mass are possible for pulsars in binary systems (Hulse & Taylor, 1975; Özel & Freire, 2016). In semidetached low-mass X-ray binaries (LMXBs), stellar masses can only be roughly estimated by the usual methods of determining the orbital parameters (Charles, 2011).

Measurements of radius are more difficult and they are usually performed for X-ray sources that are not radio pulsars. Most methods involve X-ray spectroscopy, for example during thermonuclear bursts, when the redshift may be inferred by comparing the flux at maximum radius expansion of the burst with that at “touchdown” (Ebisuzaki, 1987; van Paradijs & Lewin, 1987; Damen et al., 1990; Nättilä et al., 2016), or in the quiescent state of LMXBs (Brown et al., 1998; Marino et al., 2018). The determined radius usually depends on the mass – in this sense the spectral method yields simultaneous measurements of both the mass and the radius, although the confidence contours give rather large (banana-shaped) allowed areas in the mass–radius plot (e.g., Bogdanov et al., 2016). To date, sufficiently accurate results have not been obtained, or at least in enough systems, to effectively measure the equation of state, although constraints have been imposed.

We have found a new method for an accurate simultaneous measurement of both the mass and the radius based on the expected timing properties of accreting neutron stars at near-Eddington luminosities, if the distance is known. Potentially, if a certain type of atmospheric oscillation is detected over a wide range of X-ray fluxes in an X-ray burst source, the distance to the source will also be directly measurable from the timing properties.

Super-Eddington luminous neutron stars are observed in various astrophysical phenomena and it is understood that such high luminosities are powered by the accretion of matter on to the stellar surface or surface nuclear burning. Recent studies have reported highly luminous neutron stars to be such accretors in pulsing ultraluminous X-ray sources NGC 7793 P13, NGC 5907, NuSTAR J09551+6940.8, NGC 300 ULX1 (Bachetti et al., 2014; Israel et al., 2017a, b; Carpano et al., 2018). Super-Eddington luminosity has long been known to occur during the outbursts in a few transient X-ray binaries e.g., A0538–66, SMC X-1 and GRO J1744–28 (Skinner et al., 1982; Sazonov et al., 1997; Coe et al., 1981), and during the X-ray bursts.

Type-I X-ray bursts, are powered by thermonuclear burning of accreted material, during which neutron stars often reach Eddington luminosities (Lewin et al., 1984; Tawara et al., 1984; Lewin et al., 1993). Spectral-timing analysis of these X-ray bursts shows periodic intensity variations that are termed as burst oscillations (Strohmayer & Bildsten, 2006). Such oscillations have been observed to occur during the rise, peak and/or decay phases of the X-ray bursts with typical frequency range within 245–620 Hz (Strohmayer et al., 1996; Watts, 2012). Burst oscillations during the rise phase are usually attributed to the rotational modulation of the hot spot corresponding to flame propagation on the stellar surface. Observations of burst oscillations from nuclear powered pulsars that emit persistent accretion-powered pulsations have shown that the burst oscillation is correlated with the spin of the neutron star in SAX J108.4–3658 (Chakrabarty et al., 2003) and 4U 1636–53 (Strohmayer & Markwardt, 2002). However, the origin of burst oscillations that occur during the decay phase of the X-ray bursts is not yet well understood. Here, we explore another possibility of understanding these frequencies along with a new method to measure the mass and radius of accreting neutron stars in the context of the radial oscillations of “levitating atmospheres".

For compact stars with Eddington luminosities, radiation pressure becomes quite significant close to the stellar surface and thus has a strong influence on the dynamics of the ambient material. There is an effect unique to general relativity when examining the forces on the fluid outside a highly luminous star. If we take the fluid to be optically thin, then there exists exactly one radius, $r_{0}$, where the gravitational acceleration is balanced by the radiative force of a spherically symmetric star (Abramowicz et al., 1990). The force balance can be described by an effective potential well centred at $r_{0}$. Later studies on the motion of test particles in the combined field of gravitation and radiation have shown that the radiation drag has significant effects on the particle dynamics (Bini et al., 2009; Sok Oh et al., 2011; Stahl et al., 2012). A particle orbiting a highly luminous star has its angular momentum removed by the radiation field due to the Poynting–Robertson effect (or radiation drag effect), causing it to come to rest at this special radius $r_{0}$. For this reason, this radius is also deemed as the “Eddington capture sphere", or ECS (Wielgus et al., 2012; Stahl et al., 2013). These results generalize to a shell of optically thin fluid, in which case a stable levitating atmosphere suspended around $r_{0}$ is formed (Wielgus et al., 2015), which is completely supported by the radiation pressure, and not in contact with the star at all. It has been proposed that oscillations of these radiation supported atmospheres could be good candidates to explain the oscillation frequencies observed during Type I X-ray bursts. Abarca & Kluźniak (2016) examined the lowest order, incompressible, radial mode of such atmospheres, and found it to be overdamped by radiation drag. Encouraged by the possibility of higher order modes, Bollimpalli & Kluźniak (2017), performed a full analysis in a Newtonian framework without radiation drag, finding a family of eigenmodes with promising eigenfrequencies. For values typical for a near-Eddington neutron star, the natural frequency is right around 600 Hz. What remains to be done is a similar analysis in a general relativistic framework with a proper treatment of radiation drag, in quest for underdamped oscillations. In this paper, we perform such an analysis.

The paper is outlined as follows. In Section 2, we give a pedagogical introduction to the radius of “force balance” and to the fundamental frequencies. We shall then discuss in Section 3, the necessary equilibrium conditions for levitating atmospheres and derive a family of relativistic eigenmodes and eigenfrequencies of the radial oscillations. Damping of oscillations owing to radiation drag effect is investigated and an analytical expression for the damping coefficient is also found. In Section 4, we demonstrate how a new method to measure the mass and radius of neutron stars emerges from a unique maximum frequency treated as a function of stellar luminosity, mass and radius, as well as the oscillation mode. In Section 5, we show how observing changes in the frequency over a wide range of X-ray fluxes would allow a determination of the distance to the source as well. We conclude the paper by a discussion (Section 6) and a summary of the results (Section 7).

2 Fundamentals of levitating atmosphere oscillations

Radial oscillations of levitating atmospheres have been previously investigated for special cases (Abarca & Kluźniak, 2016; Bollimpalli & Kluźniak, 2017), and the current paper aims to present a general study of these oscillations, including radiation drag in the relativistic regime. Before turning to rigorous calculations (Section 3 and following), we present a heuristic derivation of the fundamental mode frequency of the levitating atmosphere.

The usual Eddington luminosity derivation assumes that both the radiative flux and the force of gravity are inversely proportional to the square of the radius, $r$, of a spherical surface concentric with the star. In general relativity this is no longer the case – because it is redshifted, the radiative flux decays with distance more rapidly than the gravitational acceleration – and (in the spherically symmetric case) there will be only one particular radius, $r_{0}$, at which the radiative force on the electron of a stationary hydrogen atom can balance the gravitational acceleration of the proton. At lower radii, radiation force dominates over gravity and so (in Newtonian language) the net force is directed away from the star, towards the ECS at $r_{0}$. At larger radii, beyond $r_{0}$, gravitational attraction dominates the radiation force and so the net force is directed towards the star, i.e., again towards the ECS. A stable, optically thin fluid shell may therefore be positioned at $r\approx r_{0}$. In hydrostatic equilibrium the net force (“effective gravity”) is balanced by the fluid pressure gradient – thus at $r_{0}$, where the effective gravity vanishes, the pressure is at its maximum.

If the radial extent of such a fluid shell is much smaller than $r_{0}$ (“thin atmosphere limit”), its fundamental mode of radial oscillations corresponds to uniform radial displacement by a distance $\delta s(r,t)=\delta s(t)=(1+z)\delta r(t)$ – the redshift factor ($1+z$) accounting for the difference between proper and coordinate distance – with the shell suffering a restoring force per unit mass of $-\omega^{2}\delta s(t)$, resulting in harmonic radial oscillations at the (locally observed) frequency $\omega=\sqrt{({\rm d}\,f_{g}/{\rm d}\,r)(\delta r/\delta s)}$, with $f_{g}(r)$ being the value of the net force per unit mass at $r$, and the derivative evaluated at $r_{0}$. 111The higher radial modes are more difficult to derive, as pressure variations play a role in determining the oscillatory motion. For example, the next higher mode is a breathing mode in which the oscillatory motion can be described as expansion away and contraction towards the ECS, with the (only) node at the ECS (Bollimpalli & Kluźniak, 2017).

To derive $r_{0}$ and $\omega$, we note that for an observer measuring the stellar luminosity, both the frequency of the photons and the rate of the arrival of photons are gravitationally redshifted, introducing two redshift factors, i.e., $(1+z)^{-2}$. Therefore, for a static, spherically symmetric field around a star with mass $M$ and radius $r_{*}$, a local observer at radius $r$ measures the stellar luminosity $L(r)$ to be decreasing with the radial coordinate as

[TABLE]

where $L_{\infty}$ is the luminosity measured by a distant observer $(r~{}\rightarrow~{}\infty)$, $L_{*}$ would be measured at the stellar surface, and $1+z=\left[1-{2GM}/\left({rc^{2}}\right)\right]^{-1/2}$ is the redshift factor.

Now, the gravitational acceleration of a static observer at radius $r$ is given by $GMr^{-2}(1+z)$. Therefore, at the radius $r_{0}$, where the gravitational acceleration and the radiative force of the star balance each other at local luminosity $L(r_{0})$, one has

[TABLE]

where $L_{\rm Edd}=4{\mathrm{\pi}}cGM/\kappa$ is the usual expression for the Eddington luminosity, with the opacity given by $\kappa$. Defining

[TABLE]

for our convenience, we immediately obtain

[TABLE]

from eq. (2), and the radius of equilibrium follows as (Abramowicz et al., 1990; Stahl et al., 2013)

[TABLE]

For a fixed value of $\lambda$, the radius scales directly with mass, as is usual in general relativity.

The radiative force per unit mass is $\kappa L(r)/(4{\mathrm{\pi}}cr^{2})=\lambda(GM/r^{2})[1-2GM/(rc^{2})]^{-1}$, where we used the first equality in eq. (1) and the definition of $\lambda$ and $L_{\mathrm{\rm Edd}}$. The locally observed eigenfrequency of the fundamental mode, $\omega$, and its counterpart observed at infinity, $\omega_{r}=\omega/(1+z)$, can now easily be obtained by expanding to first order in $r-r_{0}$ the expression for the effective gravity (difference between the radiative force per unit mass and the acceleration of gravity)

[TABLE]

Finally (Abarca & Kluźniak, 2016),

[TABLE]

Note the $1/M$ general relativistic (GR) scaling of the frequency:

[TABLE]

the expressions in parentheses being constant for fixed $\lambda$ in eqs. (5) – (7). As the observable frequency depends on the observed flux and the neutron star mass, it is then possible to determine the stellar mass from the two observables, $\omega_{r}$ and $L_{\infty}$. In the following sections, we will see that when radiation drag is included in the calculation, the frequency has a maximum corresponding to a location of the levitating atmosphere that is quite close to the stellar surface, $r_{0}\approx r_{*}$. As $r_{0}$ is also determined from the two observables, this allows a determination of the stellar radius as well. More precise statements can be found in the remainder of the paper.

3 Levitating atmospheres: Hydrostatic Equilibrium and Radial Oscillations

In this section, we present the full general relativistic calculations for radial perturbations of an optically thin levitating atmosphere. We adopt the natural system of units in which $G=1=c$. We assume a static spherically symmetric space–time for the background solution, thereby the line element for a Schwarzschild metric with signature ($-$,+,+,+) is given by

[TABLE]

where we have introduced

[TABLE]

We assume a perfect fluid with stress energy tensor $T^{\mu\nu}=(\epsilon+p+\rho)u^{\mu}u^{\nu}+pg^{\mu\nu}$, where $u^{\mu}$ is the four-velocity, $\epsilon$ is the internal energy, and $\rho$, $p$ are rest mass density and pressure of the flow, respectively. For the relevant gas temperatures, $T<10^{10}\,K$, it is safe to assume that $\rho\gg p+\epsilon$. We keep this simplifying assumption throughout the work.

We start with the governing equations of fluid dynamics; the continuity equation reads

[TABLE]

where $\nabla_{\mu}$ is the covariant derivative. The conservation of energy–momentum is given by

[TABLE]

where $G^{\nu}$ is the radiation four-force density (Mihalas & Mihalas, 1984). For electron-scattering-dominated regions, $G^{\nu}=\kappa\rho F^{\nu}$, where $\kappa$ is the scattering opacity and $F^{\nu}$ is the radiation flux. Using the projection tensor $h^{\nu}_{\ \mu}=\delta^{\nu}_{\ \mu}+u^{\nu}u_{\mu}$, this quantity can be determined through the radiation stress–energy tensor, $R^{\mu\nu}$, as $F^{\nu}=h^{\nu}_{\ \lambda}R^{\mu\lambda}u_{\mu}$. Under the optically thin limit, all components of $R^{\mu\nu}$ for an isotropically radiating star were first presented in Abramowicz et al. (1990) and we adopt the same here. The relevant components of the radiation stress–energy tensor are

[TABLE]

where $r_{*}$ is the radius of the star, and $\alpha$ is the viewing angle, which is the maximum polar angle up to which an observer at radius $r$ can see the photons from the stellar surface. For simplicity, we assumed that the stellar radius is beyond the photon orbit, $r_{*}>3M$, which gives222Refer to Abramowicz et al. (1990) for details.

[TABLE]

The dependence of $\sin\alpha$ on $r_{0}$ is shown in Figure 1 for several stellar radii.

The radial component of the relativistic Euler equation for an optically thin gas subject to the radiation field becomes

[TABLE]

We have assumed a uniformly radiating star, which may be a good description of a neutron star undergoing an X-ray burst. For an accreting and rotating neutron star, a rotating, radiating and levitating belt in the boundary layer has been shown to form due to the balance of gravitation, centrifugal and radiation forces; as the accretion rate approaches the near-Eddington luminosities, this belt extends from the equatorial plane up to the pole, allowing the entire surface of the neutron star to radiate (Inogamov & Sunyaev, 1999). Thus, the spherically symmetric envelope may be a useful approximation even for a rotating neutron star accreting through a disc.

3.1 Geometrically thin background solution

Initially, the atmosphere is in hydrostatic equilibrium with a normalized four-velocity, $u^{t}$, given by $(B^{-1/2},0,0,0)$. Following equation (16), the hydrostatic equilibrium condition for a spherically symmetric and optically thin levitating atmosphere around a uniformly radiating star is (Wielgus et al., 2015)

[TABLE]

Note that the location of the pressure maximum from the above equation, at $\lambda=B_{0}^{1/2}$, corresponds to the location of the ECS, $r_{0}$, see also equation (5),

[TABLE]

Quantities with subscript [math] are evaluated at the ECS, e.g., $B_{0}\equiv B(r_{0})=1-2M/r_{0}$. Analytical solutions of an optically and geometrically thin polytropic atmosphere show that the density and pressure have Gaussian-like profiles (Wielgus et al., 2015) with the maximum located at $r_{0}$. Here, we present such solutions in a form usable for the calculation of normal modes. The procedure is quite similar to the one employed in the task of finding normal modes of oscillating tori (see, e.g. Blaes et al. (2006). The main difference is that in the currently considered case the radiation force, rather than the centrifugal force, balances gravity in the background solution. The simplification is that the current problem has a spherical symmetry; hence, we obtain ordinary (and not partial) differential equations. Denoting $\mathcal{E}=-u_{t}=B^{1/2}$, we define an effective potential $\mathcal{U}_{\rm eff}$ such that

[TABLE]

Using equation (17) and integrating, we find

[TABLE]

Casting the equation of hydrostatic equilibrium as

[TABLE]

allows us to express the right-hand side of the equation as a gradient of a scalar function, $\psi$. For a polytropic fluid with an adiabatic index $n$, $p\propto\rho^{(1+1/n)}$, we may integrate equation (21) to find the Bernoulli equation in the following form:

[TABLE]

After evaluating the constant at the ECS and some reordering, this integral condition can be cast as

[TABLE]

with the speed of sound (squared) $c_{s}^{2}=dp/d\rho$.

So far we did not explicitly assume that the background solution should be geometrically thin. We make this simplification now, by approximating the right-hand side of equation (23) with the second-order Taylor expansion around $r_{0}$, noticing that

[TABLE]

The approximation holds for atmospheres that are sufficiently geometrically thin to neglect the gradient of the potential, i.e., the radial variation of space–time curvature and the radiation field, in comparison with the gradient of pressure. Hence, it is acceptable for atmospheres with thickness much less than $r_{0}$. The following approximate formula is obtained:

[TABLE]

Seeing how

[TABLE]

this result can be cast in a more compact form using the following substitutions:

[TABLE]

with a simple interpretation: $\omega_{r}$ corresponds to the fundamental oscillation angular frequency, already familiar from equation (6), parameter $\beta$ controls the geometric thickness of the atmosphere, $x$ and $\eta$ correspond to a convenient scaling of the radial coordinate. Thus, equation (25) for the radial profile of pressure may be rewritten simply as

[TABLE]

3.2 Perturbations of the geometrically thin solution

We assume coherent, spherically symmetric perturbations of physical variable $X$ in the form of

[TABLE]

and linearly perturb the relevant hydrodynamic equations. The first thing to notice is that $\delta u^{t}=0$ is necessary to fulfil the perturbed four-velocity norm, as $u_{r}=0$ in

[TABLE]

The linearly perturbed version of the continuity equation takes the following form

[TABLE]

The perturbed radial Euler equation yields

[TABLE]

where $\delta\mathcal{D}=\kappa\left(u_{t}^{2}R^{tt}+g_{rr}R^{rr}\right)\delta u^{r}$ relates to the radiation drag acting against the radial motion of the particle and can be rewritten as

[TABLE]

cf. equations (12)–(16). Although $\chi$ has a radial dependence through the cosine and sine functions of the viewing angle, in the geometrically thin limit $r\sim r_{0}$, hence $\chi$ will be taken to be constant for given stellar parameters and stellar luminosity.

Following Ipser & Lindblom (1992) and Abramowicz et al. (2006), we define a new variable,

[TABLE]

From the above definition, and equation (35), it follows that in the geometrically thin limit

[TABLE]

Substituting $W$ in equation (34) and combining with equation (38) gives a single second-order ordinary differential equation. In the geometrically thin atmosphere limit ($\beta\rightarrow 0$), this equation simplifies to

[TABLE]

with $\omega_{r}^{2}$ given by equation (27). Note that equation (39) can be cast into a convenient dimensionless form

[TABLE]

with

[TABLE]

and with the specific viewing angle given by

[TABLE]

As we will see, $\tau=1/(\omega_{r}\gamma)$ is the decay constant (damping time scale) of the oscillatory solutions, so it is important to compare the value of the damping coefficient $\gamma$ to the values $\sigma$ of the (dimensionless) oscillation frequencies for the undamped solutions.

It is easy to show that outside the photon orbit, i.e., for $r_{0}\geq 3M$, the viewing angle decreases monotonically with $\lambda$ (or with the ECS radius $r_{0}$), and so does the damping coefficient $\gamma$. The maximum value of $\gamma=4/3$ is attained at the stellar radius $r_{0}=r_{*}$, and $\gamma\rightarrow 1$ as $r_{0}\rightarrow\infty$, $\alpha\rightarrow 0$, $\lambda\rightarrow 1$. Thus, the radiative drag is always present in the problem at hand, with $4/3\geq\gamma>1$. None the less, it is instructive to consider a simplified problem in which radiative drag is artificially neglected.

3.3 Undamped oscillations

Let us now consider the artificial problem of undamped oscillations. This is obtained by blithely neglecting the damping term in equation (40), i.e., putting $\gamma=0$, which leads to the Gegenbauer equation,

[TABLE]

The corresponding solutions represent the undamped oscillations with the relativistic eigenfrequencies observed at infinity given by333Which can also be identified directly from equation (55) with the damping coefficient set to zero.

[TABLE]

where

[TABLE]

The relativistic frequencies differ from the Newtonian ones (Bollimpalli & Kluźniak, 2017) by a redshift related factor of $\lambda^{2}$, while the relativistic eigenmodes are similar to the Newtonian eigenmodes. Nevertheless $\sigma_{k}$, the ratio of normal mode frequencies to the fundamental mode frequency, remains the same for both the relativistic and non-relativistic cases. The ratio of frequencies of the second (“breathing") mode and the fundamental mode is $\sigma_{2}=\sqrt{2+1/n}$.

The eigenfrequencies grow with the mode number; the first two are $\sigma_{1}=1<\gamma$ and $\sigma_{2}^{2}={2+1/n}>2>\gamma^{2}$. Thus, in the full problem, which includes radiative damping, the radial mode ($k=1$) is overdamped (Abarca & Kluźniak, 2016), but all the higher modes ($k\geq 2$) have an oscillatory character (are underdamped). This is also clear when we plot the frequencies, $\nu=\omega_{k}(r_{0})/(2\pi)$, of the first ten modes of undamped oscillations (left to right), as a function of Eddington parameter or atmosphere location (Figure 2).

The frequencies are plotted in Hz to facilitate comparison with observations. The values on the left vertical axis correspond to $M=1.4\,{\rm M}_{\odot}$. For a different mass, $M$, the frequency at a given $r_{0}$ and a given mode can be obtained with the scaling factor $M/(1.4{\rm M}_{\odot})$ shown on the left vertical axis. Frequencies computed for $2.1{\rm M}_{\odot}$ are shown with similar scaling on the right vertical axis. The radius in units of the stellar mass has to be assumed a priori to compute these frequencies, we took $r_{*}=5\,M$. For a given $\lambda$ or $r_{0}$, higher number modes have larger frequencies. The further the atmosphere is located from the stellar radius, the lower the frequencies get as the magnitude of the restoring forces decreases with the radius. This radial dependence is fully accounted for by the functional form of $\omega_{r}$, given in equation (27). We note that there is a range of luminosities for which the frequencies of these oscillations fall in the range of 300–600 Hz, typical for observed frequencies of X-ray burst oscillations.

3.4 Damped oscillations

We now turn to the full problem including the damping effects of the ever-present radiation drag in radiation-supported atmospheres.

It is easiest to solve equation (40) by considering the real ($\sigma_{\rm R}$) and imaginary ($\sigma_{\rm I}$) parts of $\sigma$ separately,

[TABLE]

and the corresponding relation $\omega=\omega_{\rm R}+i\omega_{\rm I}$, with $\omega_{\rm R}=\omega_{r}\sigma_{\rm R}$ and $\omega_{\rm I}=\omega_{r}\sigma_{\rm I}$. Substituting complex $\sigma$ into equation (40) results in a differential equation with real and imaginary parts that are separable and can be solved simultaneously. The real part is an eigenvalue problem in the form of a Gegenbauer differential equation,

[TABLE]

while the imaginary part yields

[TABLE]

A particularly simple solution corresponds to the overdamped case of purely imaginary $\sigma$. The eigenvalues are given by, $-\sigma_{\rm I}^{2}-2\gamma\sigma_{\rm I}=\sigma_{k}^{2}$, giving,

[TABLE]

In fact, as follows from the discussion above, this solution only occurs for $k=1$, for which $\sigma_{k}=1$, corresponding to the imaginary frequency of the overdamped solution of

[TABLE]

with

[TABLE]

The time dependence of this decaying solution is $\exp(-i\omega\,t)=\exp(\omega_{\rm\,OD}\,t)$.

In the general case, equation (48) can be solved with a proper boundary condition of pressure vanishing on the atmosphere boundary. The corresponding eigenmodes are represented by the Gegenbauer polynomials $C^{\alpha}_{k}(\eta)$ with $\alpha=n-\frac{1}{2}$ and $k=0,1,2,3,...$ which can be calculated with a recursive formula (Koornwinder et al., 2010)

[TABLE]

They form a complete set of modes, that is, any spherically symmetric oscillation of the levitating atmospheres is necessarily a combination of the eigenmodes described by the equation (53). For example, the eigenfunctions of the second mode ($k=2$) and the third mode ($k=3$), take the form of $\left(n-\dfrac{1}{2}\right)[\eta^{2}(2n+1)-1]$ and $\dfrac{\eta}{6}(4n^{2}-1)[(2n+3)\eta^{2}-3]$, respectively. 444As a reminder, $n$ is the adiabatic index and $\eta$ is a scaled radial coordinate, (see eq. 30). The above imply a non-vanishing eigenfunction W, and for $k>1$ also $\sigma_{\rm R}\neq 0$, so equation (49) simply gives the damping rate

[TABLE]

Using this, the corresponding eigenfrequencies of the oscillations are computed from the Gegenbauer relation,

[TABLE]

where $k=2,3,4...$. Note that the time dependence of the solution is $\exp(-i\omega_{\rm R}\,t-\omega_{r}\gamma\,t)$, so $1/(\omega_{r}\gamma)$ is the damping time. The damping coefficient $\gamma$ has been defined in equation (42) and the fundamental frequency $\omega_{r}$ in equation (27).

The frequency $\omega_{\rm R}^{2}$ is a product of the fundamental frequency $\omega_{r}$ squared, which is a rapidly and monotonically decreasing function of the ECS radius $r_{0}$, and of a monotonically increasing function of the same variable, $\sigma_{k}^{2}-\gamma^{2}(r_{0})$. The product has a maximum, close to the stellar surface, even though the damping coefficient $\gamma$ of equation (42) is a slowly varying function of the viewing angle, and $\sigma_{k}$ is a constant. This is because close to the stellar surface, the angle $\alpha$ itself varies rapidly with $r_{0}$. Indeed,

[TABLE]

and the second factor dominates the radial derivative of $\omega_{\rm R}^{2}$, as $\cos\alpha\rightarrow 0$ at $r_{0}\rightarrow r_{*}$, where $\alpha\rightarrow{\mathrm{\pi}}/2$, while the other terms are regular everywhere:

[TABLE]

[TABLE]

Thus, close to the stellar surface $\omega_{\rm R}$ is a growing function of $r_{0}$, while at larger radii, over most of its range, $\omega_{\rm R}$ decreases with distance to the star. Therefore, for any mode higher than the fundamental ($k>1$), a maximum of the atmospheric oscillation frequency occurs close to the stellar surface. The fundamental mode ($k=1$) is overdamped, and hence the maximum occurs in the inverse decay constant $|\omega_{\rm\,OD}|$, as can be seen in Figure 3.

Figure 3 presents the frequencies $\omega_{\rm R}(r_{0})/(2\pi$) of the initial few lowest modes of the underdamped oscillations in an increasing order of mode number from left to right. One may immediately notice from the figure that the frequencies of the modes have now decreased slightly due to damping. For the fundamental mode, we plot the magnitude, $|\omega_{\rm\,OD}|/(2\pi)$, shown by the thin dashed line in the plot. The overdamped solutions are not of much astrophysical interest, so in further discussions we shall only consider the second and higher modes, as they allow for oscillations.

Observed oscillations in the X-ray bursts occur mainly during the rise and decay phase of the outburst (Muno et al., 2001). Interestingly, even with the damping, the observed range of 300–600 Hz still falls within the frequency range of these mode oscillations. Oscillations during the outburst decay phase are found to have increasing frequency with time (Strohmayer & Bildsten, 2006). Eq. (55) and Figure 3 show that the frequencies increase with decreasing luminosity (as long as the frequency maximum is not reached), which is in accord with the observed oscillations, since luminosity decreases with time during the decay phase of the outbursts, as the name suggests. However, the magnitude of the frequency increase in the optically thin model exceeds the observed changes.

In the following section, we present a method of simultaneous mass and radius measurements from the frequency maximum of the damped oscillations.

4 Characteristic frequency maximum: Mass and radius determination

While the frequency of the oscillations in general decreases with luminosity, a distinctive feature to be noted in Figure 3 is the maximum in the frequencies of underdamped oscillations close to the stellar surface. For given stellar parameters, $M$, $r_{*}$, we mark the radius at which the atmosphere can oscillate with maximum frequency, $\nu_{\rm max}$, as $r_{\rm max}$. The maximum for the second mode is clearly pronounced. For higher modes, since $r_{\rm max}$ is quite close to the stellar radius, it would be difficult to observe the decrease of frequency to the left of the maximum as the luminosity decreases. The presence of the frequency maximum can be identified as due to the steep radial decrease of damping close to the star which is mostly attributed to a rapid decrease in viewing angle, as remarked in Section 3.4, in the discussion of eq. (56).

Figure 4 shows for the first few underdamped modes the coordinate distance between the location of the atmosphere that corresponds to the maximum frequency and the stellar surface. Lower modes are comparatively farther away from the stellar surface, while the distance also increases with the stellar radius.

For a particular mode, this maximum of the frequency, $\nu_{\rm max}$, depends on two parameters ($r_{*},~{}M$). In Figure 5, we plot the frequencies of the second mode for $1.4\,{\rm M}_{\odot}$ and $2.1\,{\rm M}_{\odot}$ stellar mass and various stellar radii. The following can be noted from the figure:

a) For a given stellar mass, $\nu_{\rm max}$ decreases with increasing $r_{*}/M$.

b) For a given $r_{*}/M$, $\nu_{\rm max}$ is inversely proportional to the stellar mass, as expected for GR scaling.

c) Irrespective of the values of $r_{*}/M$ and $M$, $r_{\rm max}$ is always located close to the stellar surface.

d) There is a degeneracy in the frequency maximum value that occurs for different combinations of $r_{*}/M$, $r_{\rm max}/M$ and $M$. For example, two different stars with masses $1.4\,{\rm M}_{\odot}$ and $2.1\,{\rm M}_{\odot}$ and radii $5.945\,M$ and $4.835\,M$ exhibit the same frequency maximum value of $700$ Hz when their respective locations of the atmosphere are at $6.1\,M$ and $4.95\,M$. These frequencies are plotted as green dash-dotted lines in Figure 5. As we shall discuss below, the stellar luminosity breaks this degeneracy and allows us to obtain the stellar parameters unique to the corresponding frequency maximum.

e) However, for a given $r_{*}/M$, the frequency maximum occurs at the same radius (in units of $M$), independent of the stellar mass. This follows immediately from equation (55), which is a function of $r_{*}/M$ and $r_{0}/M$ alone, once the $1/M$ dependence of $\omega_{r}$ is taken out. So, the condition for the maximum, i.e., the first derivative of the frequency with respect to $r_{0}/M$ being set to zero, gives a relation just between the two radii $r_{*}/M$ and $r_{\rm max}/M$.

Given the frequency maximum value and the corresponding luminosity, we can determine the mass and radius of the neutron star. Indeed, we already know that for any single mode exhibiting a frequency maximum, we can now find the value of the dimensionless stellar radius, $r_{*}/M$, if we know $r_{\rm max}/M$. We require one last relation to close the system of equations that enable us to solve for the two parameters ($r_{*},\,M$). Equation (18) fulfils this requirement by providing a relation between $M$ and $r_{\max}$ through $\lambda$ or stellar luminosity.

To derive the stellar parameters, the set of nonlinear relations are solved semi-analytically using the following protocol:

i) We derive $\lambda$ as a function of $M$ and $L_{\infty}$ from equation (3), assuming the Eddington limit for ionized hydrogen. An atmosphere with a different composition would have a different Eddington limit and the $x$-axis of the plots in Figure 6(a) and 6(b) would then vary accordingly.

ii) With the direct dependence on $\lambda$, $r_{\rm max}/M$ is then determined as a function of $M$ and $L_{\infty}$, using equation (5) with $r_{0}=r_{\rm max}$.

iii) For the assumed mode number (say $k=2$ or $k=3$), we set the first derivative of frequency with respect to $r_{0}$ to zero to derive $r_{*}/M$ as a function of $r_{\rm max}/M$, and therefore as a function of $M$ and $L_{\infty}$.

iv) The final equation for the frequency maximum is thus solely a function of mass and luminosity, so that given any two quantities, the other can be determined. Thus, given that we know the observed frequency maximum and luminosity, we can determine the mass of the star.

v) Knowing $M$ and $r_{*}/M$ we determine the stellar radius.

We admit a source of ambiguity in the measurements of radius and mass due to the mode number. The mode number is a discrete parameter, and so assuming different mode numbers introduces large discrete changes in the measured mass and radius. We expect the higher modes to be less likely to be excited than the breathing mode ($k=2$), given their more complicated eigenfunctions. However, in any case, we believe it would in principle be possible to determine the mode number by observing the frequency variation with varying flux (luminosity) – the higher the mode number, the more rapid the variation, as is clear from Figure 3.

In Figures 6(a) and 6(b), we show a grid of stellar mass and radius contours for a given range of luminosity and frequency maximum computed as described above for the second and third modes, respectively. The coloured contours (sloping upwards with luminosity) are for stellar radius in increasing steps of 0.5 km from top to bottom. The nearly vertical black contours represent mass in range of $1\,{\rm M}_{\odot}$ to $4\,{\rm M}_{\odot}$ in increasing steps of $0.05\,{\rm M}_{\odot}$ from left to right. At a given luminosity, $\nu_{\rm max}$ of a star increases with its compactness. Since we assumed that the stellar surface lies above the photon orbit, we omit solutions with $r_{*}<3M$. Hence, for combinations of higher luminosity and higher frequency maximum, we see no contours in the figures.

As we can note from the grid in the plots, each pair of stellar luminosity and frequency maximum thus gives a unique measure of the radius and mass of the neutron star. For example, if the observed frequency maximum is around 600 Hz and the corresponding luminosity of the star observed at this frequency is around $1.9\times 10^{38}~{}\mathrm{erg~{}s}^{-1}$, then assuming that the frequency corresponds to the second mode ($k=2$), the determined stellar radius is $14.5$ km and mass is $1.65~{}{\rm M}_{\odot}$. If we assume the third mode ($k=3$), then the determined stellar radius increases to $18.5$ km and the mass estimate decreases to $1.6~{}{\rm M}_{\odot}$.

Inferring the mass and radius of the neutron star through this method thus requires only a single observation of two quantities, luminosity and frequency. However, this method is subject to certain systematic uncertainties, of which the most important is the uncertainty in the distance to the source, which translates directly into uncertainty in the luminosity, thereby leading to errors in the estimated mass and radius. We discuss these in the next subsection.

4.1 Error budget

The determination of the neutron star mass and radius, described in this section, would be quite accurate if the distance to the source were known quite precisely, as in the 47 Tuc globular cluster, with quoted distance $d=4.53\pm 0.06\,$kpc (Woodley et al., 2012). Unfortunately, the distance errors to other globular clusters are less constrained, with uncertainties sometimes approaching $10\%$ (Guillot et al., 2013; Steiner et al., 2018). Typical errors on the distances to the neutron star LMXBs in the field are even larger, falling between $10-20\%$ (Galloway et al., 2008). Such large uncertainties in the distance imply that the typical uncertainties in luminosity can be around $20\%$ or higher.

Frequencies of oscillations, once detected, are typically much more precise. Planned, sensitive X-ray instruments e.g. STROBE-X (Ray et al., 2019) and eXTP (Zhang et al., 2016), should be able to measure the frequencies to exquisite precision, given by the ratio of the oscillation period to the time-scale for frequency change (as the luminosity varies), which should be less than one part in a hundred. To be conservative, in the remainder of this section we assume a frequency error of $3\%$. All the analysis throughout the rest of this paper considers only the second mode, i.e., $k=2$.

To quantify the mass and radius errors, we do a simple $\chi^{2}$-error analysis with chosen uncertainties in luminosities and frequencies. We assume that the luminosity and frequency maximum measurements are normally distributed, with different mean values for different cases. We consider the systematic uncertainties on the luminosity and frequency to be $20$ per cent and $3$ per cent respectively, which are taken to be the standard deviations of the corresponding distributions. We consider six cases in total with the mean luminosities ($\bar{L}$) to be $1.8\times 10^{38}$ and $2.0\times 10^{38}~{}{\rm erg~{}s}^{-1}$, and the mean frequency maxima ($\bar{\nu}_{\rm max}$) to be 300, 600, and 900 Hz. Figure 7 shows the joint confidence regions for mass and radius for the six different cases with $\bar{L}$ and $\bar{\nu}_{\rm max}$ given in the title of the each plot. The confidence regions are determined based on the $\Delta\chi^{2}$ values calculated with the chosen uncertainties. Overplotted black contours denote the $1\sigma$, $2\sigma$, and $3\sigma$ errors that correspond to the $68.27\%$, $95.45\%$, and $99.73\%$ confidence levels, respectively.

A projection of the $1\sigma$ contour ($\Delta\chi^{2}=1$) on to the corresponding axes yields upper limits to the $1\sigma$ errors on mass and radius, but the actual $1\sigma$ errors on the parameters are always smaller. As we cannot assert that these parameters have a Gaussian distribution, we turn to a Monte Carlo analysis to generate the probability distributions for mass and radius and to accurately determine the errors on their values.

We once again consider the six cases with the assumption that the frequency and luminosity are normally distributed with the mean values as given in the first two columns of Table 1 and with the same $3\%$ and $20\%$ uncertainties as before. For each such case, we perform a Monte Carlo simulation to generate a large ensemble of $10^{8}$ pairs of luminosity and frequency maximum values that are normally distributed with the chosen $\bar{\nu}_{\rm max}$ and $\bar{L}$. For each such pair of luminosity and frequency maximum, we compute the corresponding pair of mass and radius to generate the posterior distributions for the mass and radius. In Figure 8 we show the distributions for mass and radius for one such case ($\bar{\nu}_{\max}=600$ Hz and $\bar{L}=2.0\times 10^{38}~{}{\rm erg~{}s}^{-1}$). The red solid line represents the median (50th percentile), while the red dashed lines denote $1\sigma$ error bounds (16th and 84th percentile). The central plot shows the contours of joint confidence regions at the levels 0.68 and 0.95. The $1\sigma$ errors on mass and radius determined through this method are given in the third and fourth columns of Table 1. The typical errors on mass and radius from this table are slightly larger than $20$ per cent and $10$ per cent, respectively.

Therefore, we can say that with respective errors of $20\%$ and $3\%$ in luminosity and frequency, we can estimate the parameters of the neutron star with an accuracy of $20\%$ for mass and $10\%$ for radius. It is easy to see that the errors on mass and radius will be significantly lowered if the errors on luminosity can be decreased. Indeed, the uncertainty in the determined values of mass and radius are of the order of the uncertainty in the luminosity.

Interestingly, if the oscillation frequency is observed over a wide range of fluxes, the neutron star mass and radius could be determined without prior knowledge of the distance to the source. In the following section, we discuss how to simultaneously determine the mass, radius, and distance to the source, and show that the errors on $M$, and $r_{*}$ would be decreased by an order of magnitude relative to the errors discussed above.

5 Mass, radius, and distance determination

One way to circumvent the problem with distance uncertainties is to consider the ratio of luminosities that essentially translates into flux ratio, which is a directly measurable quantity like frequency. The two green dashed-dotted lines in Figure 5 that share the same frequency maximum for different masses and radii show visibly different frequency variation with $\lambda$ ($\propto~{}L/M$), but the mass, and therefore $\lambda$ as well, is not yet known, so in practice we cannot use this plot. In fact, let us assume that we do not know the distance to the source, so that we do not know the luminosity. None the less, if we consider plotting the frequency variation with the flux ratio, which does not depend on the distance, we can distinguish the curves for different masses and radii. Suppose that we have a set of observed data points for frequency variation with flux. In principle, as per our model, there exists a unique pair of mass and radius, within error, that can fit all the observed data points. Although, in reality, if we do not have a sufficient number of data points, or data spanning a sufficiently wide range over frequencies, then a wide range of masses and radii can potentially fit the observed data points, thus making it difficult to stringently constrain the mass and radius (unless we already know the luminosity, in which case we can proceed as in Section 4).

To elaborate on this, we consider an example of frequency variation with luminosity for the parameters $M=1.5~{}M_{\odot}$ and $r_{*}=13.72$ km. The curve has $\nu_{\rm max}~{}=$ 600 Hz and we chose five data points on this frequency curve corresponding to 550, 500, 450, 400 and 350 Hz with $1\%$ uncertainty. We can use these five data points as a proxy for the real observation data set and check to within what accuracy we can recalculate the values of mass and radius by finding the best fits to these points. Since we do not know the distance, these points can only be represented on the frequency-flux plane if we scale the corresponding luminosities of each point with the luminosity of any one particular point. This is the same as scaling the fluxes of the observed data points to the flux of one data point. In this example, we chose this reference point to be the one with frequency 600 Hz.

First, let us suppose that we have a precise measurement for the frequency maximum. In such a case, we can try to fit through all the points with the frequency curves for various mass and radius pairs that have the same $\nu_{\rm max}~{}=~{}600$ Hz, as shown in Figure 9(a). The $x$-axis for each curve is scaled with its corresponding flux (${\rm F}_{\nu_{\rm max}}$) at $\nu_{\rm max}$. The lower we go in the frequency data compared to $\nu_{\rm max}$, the better we can constrain the mass and radius measurements. An inspection of the plot close to the error bars located at $400$ Hz and $350$ Hz shows that we can already limit the errors on mass to a much better accuracy than $\pm~{}0.1~{}{\rm M}_{\odot}$. This can be better visualized in Figure 9(b), which shows fits to the error bounds for each data point. This way the errors on mass and radius determined from the fits to the bounds of the 350 Hz error bar are a factor of $6$ less compared to the errors from the fits to the 550 Hz error bar. So given a precise measurement for $\nu_{\rm max}$, and a wide range of frequency measurements, the mass and radius can be estimated with relatively minor errors of $<2\%$ ($\pm 0.027\,{\rm M}_{\odot}$) and $<1\%$ ($\pm 0.121$ km) respectively.

Another potential source of error is from a measurement of the frequency maximum with given uncertainty just like the remaining data points. We fit three curves to the centre, lower end, and upper end of the error bars as shown in Figure 10. The respective masses are labelled in the plot, and we see all the curves passing through the “data” points in Figs. 9 and 10 yield values of $M$ and $r_{*}$ within an error of less than 2%.

Thus, we conclude that the overall error in the neutron star mass and radius determination by this method is of the same order as the assumed error in the measured frequencies ($~{}1$%), which would be an excellent result.

5.1 A note on distance estimates

With the above derived constraints on mass and radius for a known $\nu_{\rm max}$, we can determine the luminosity and thereby the distance to the source with a high accuracy. In Figure 11, we show a grid of $\nu_{\rm max}$ (nearly diagonal lines in black) and luminosity contours (coloured downward sloping lines) plotted for a wide range of masses and radii. The red crosses lying close to the intersection of $\nu_{\rm max}=600$ Hz and $L=1.70\times 10^{38}\,{\rm erg~{}s}^{-1}$ contours mark the mass and radius errors estimated with $1\%$ uncertainty in $\nu_{\rm max}$, as just described, and with these, the luminosity can be determined with an accuracy of $1\%$. This implies that the distance to the source can also be constrained to within $1\%$ errors.

The technique may have astrophysical applications. Type-I X-ray bursts with photospheric radius expansion have been treated as standard candles to evaluate distances to globular clusters containing these X-ray bursters with uncertainties on the level of $30\%$ (Kuulkers et al., 2003). Measurements of LMXBs in quiescence seem to be more precise, with errors less than $10\%$ (Guillot et al., 2013; Steiner et al., 2018). With our model, if it is successfully applied to X-ray bursters, we could place much more stringent constraints on the distances to these globular clusters. For clusters such as 47 Tuc, where the distance has been estimated accurately through other methods (Woodley et al., 2012), one would get a valuable cross-check on the validity of our approach.

6 Further discussion

We have derived the modes of oscillation of an atmosphere suspended above the neutron star surface by the force of radiation. Although the fundamental mode is critically damped due to radiation drag, the higher oscillatory modes of optically thin atmospheres may be detected by future X-ray missions. We have shown that with the knowledge of the frequency and the corresponding near-Eddington luminosity from the X-ray observations, the mass and radius of the neutron star can be derived. It would be interesting to check the applicability of this method with the currently available X-ray data, to compare with the already available estimates of the neutron star mass and radius.

Recent works on constraining the mass and radius of neutron stars rely on X-ray spectroscopy of the thermal emission from the hot neutron stars residing in quiescent LMXBs. An analysis of LXMB sources in the Galactic plane gave mass and radius constraints on two sources with respective errors at $90\%$ confidence level between $~{}30\%$ and $~{}50\%$ (Marino et al., 2018), and comparable, or even higher, errors for several globular cluster sources (Steiner et al., 2018). These error estimates assume a precisely known distance, the errors further increase when accounting for distance uncertainties. For quiescent LMXBs X7 in the globular cluster 47 Tuc, the distance to which is known accurately, Bogdanov et al. (2016) obtains a $7\%$ error in the radius, assuming a particular value for the neutron star mass. Our method, which is described in Section 4, is potentially much more accurate (and yields both, the mass and radius), but for uncertain distances it yields errors as large as those of the distance determination. However, we note that the slope (in the mass–radius diagram) of the confidence level contours differs between the methods. The spectral technique employed to various LMXB sources in globular clusters yields mass and radius constraints with probabilities shifting towards larger radii for lower masses in most of the sources (Steiner et al., 2018). This picture when overlapped with Figure 7 of this paper in which the probabilities shift towards larger radii for larger masses can further narrow down the limits on mass and radius. Thus in an overall comparison, burst oscillation study can be a very useful technique to constrain the mass and radius with a high accuracy.

As we showed in Section 5, if the oscillation is detected over a sufficiently wide range of X-ray fluxes, the distance to the source may be directly determined by comparing observations with the theoretical frequency–flux curves. In such a case, the mass and radius, as well as the distance, could be determined to high accuracy with no reference to the spectral method.

Apart from the oscillations and their relevance to the observed frequencies, the levitating atmosphere may have its own observational signatures. A study by Rogers (2017) suggests that such levitating atmospheres can deflect the light rays coming from the central compact source and this can significantly affect the appearance of the central object in the observations.

Similar atmospheres are found in the corona of accretion discs and described as photon floaters, formed at the critical height where the radiation force from the disc is balanced by the gravitational force of the central object (Fukue, 1996). It would be interesting to see if similar oscillations persist in the disc corona. Under the plane-parallel approximation, the eigenfunctions of these atmospheric oscillations may remain the same as for the spherically symmetric geometry presented here, while the frequencies of such oscillations would be different.

So-called burst oscillation frequencies are detected during the Type I X-ray bursts. If these are identified with the actual oscillations of an optically thin atmosphere, the methods described in this paper can be directly applied to the X-ray burst frequencies in order to determine the neutron star mass and radius. However, the radiating surface in X-ray bursters is optically thick and throughout this work we have assumed an optically thin fluid. The construction of optically thick atmospheres was given by Wielgus et al. (2016) and the procedure involved numerical methods. So an analytical analysis of the oscillation modes in such a case is non-viable and their analysis is beyond the scope of this paper. Clearly, they are of great astrophysical interest and could possibly be subject to less radiation drag since radiation transport inside the optically thick atmosphere is due to diffusion. It is hoped that future work will treat the oscillations of the optically thick levitating atmospheres, especially in the context of radius expansion X-ray bursts. Regardless of this, future instruments with high time resolution and good sensitivity to higher energy photons, such as STROBE-X (Ray et al., 2019) and eXTP (Zhang et al., 2016), may detect optically thin atmospheres, and their oscillations, outside Eddington luminosity neutron stars.

7 Conclusion

Neutron stars with super-Eddington luminosity have a strong radiation field close to the stellar surface that dominates over gravity. The strong radiation force pushes the surrounding matter away from the stellar surface to a certain critical radius $r_{0}$, beyond which gravity prevails, thereby forming an atmospheric shell at $r_{0}$ that levitates above the stellar surface. These levitating atmospheres are thus supported by the strong radiation flux from the stellar surface. The same flux is a source of strong radiation drag that eventually damps any oscillations of the atmosphere.

In this paper, we analytically study in general relativity the oscillations of levitating atmospheres that are optically and geometrically thin, including the radiative terms in the perturbation equation. Radiation drag induces damping of oscillations in all regimes and we find that the damping coefficient is independent of the mode number. Radiation drag prevents oscillations only for the first mode, while higher modes are underdamped. The frequency range observed for the burst oscillations in the decay phase of the X-ray bursts is 300–600 Hz, which lies within the obtained frequencies of the damped oscillations. The frequency of these oscillations increasing with decreasing luminosity is in qualitative agreement with the observations for the oscillations in the decay phase of the X-ray burst where luminosity decreases with time.

The frequency of these oscillations exhibits a characteristic maximum that is more pronounced for the lower order modes. Noting that the frequency is only dependent on the stellar luminosity, mass and radius, we compute the mass and radius as a function of the maximum frequency, the luminosity. However, the accuracy of this method is limited by the uncertainties in distance to the source, typically of $\sim 10\%$, which lead to a comparative relative error on mass and radius estimates.

An alternative method to determine the stellar parameters is also proposed, which does not rely on prior knowledge of the distance. This method focuses on the variation of the two directly observable quantities with respect to each other, i.e., the frequency variation with the flux. We show that with this method, the stellar parameters could be estimated as accurately as the frequencies are measured, i.e., to $\sim 1\%$ should such oscillations be observed over a sufficiently wide range of frequencies. The new constraints on the mass and radius when combined with the knowledge of the frequency maximum allow us to also directly infer the luminosity of the source with an accuracy of $<1\%$. Therefore, the distance can also be constrained to within a $1\%$ error, by comparing the inferred luminosity with the measured flux. This model thus potentially establishes a new way to determine, with exceptionally high accuracy, the distance to the neutron star, and its mass and radius.

Acknowledgements

The authors would like to thank the anonymous referee for the interesting and helpful comments that helped improve the quality of the paper. This research was partly supported by the Polish NCN grants UMO-2013/08/A/ST9/00795, 2013/10/M/ST9/00729 and in part by the National Science Foundation under Grant No. NSF PHY11-25915. MW acknowledges support by the Black Hole Initiative at Harvard University, through the grant from the John Templeton Foundation.