Radial modes of levitating atmospheres around Eddington-luminosity neutron stars

TL;DR

This paper analyzes the oscillation modes of a thin, radiation-supported atmosphere around a neutron star near the Eddington luminosity, deriving eigenfrequencies based on the atmosphere's properties.

Contribution

It provides a theoretical calculation of the eigenmodes and frequencies of levitating atmospheres in a Newtonian framework, highlighting the dependence on the adiabatic index.

Findings

Eigenfrequency ratio depends only on the adiabatic index, b3+1.

The atmosphere's oscillation frequencies are derived analytically.

The study enhances understanding of neutron star atmospheres near Eddington luminosity.

Abstract

We consider an optically thin radiation-supported levitating atmosphere suspended well above the stellar surface, as discussed recently in the Schwarzschild metric for a star of luminosity close to the Eddington value. Assuming the atmosphere to be geometrically thin and polytropic, we investigate the eigenmodes and calculate the frequencies of the oscillations of the atmosphere in Newtonian formalism. The ratio of the two lowest eigenfrequencies is $\sqrt{\gamma+1}$, i.e., it only depends on the adiabatic index.

Tables

Table 1. Table 1: Eigenmodes and eigenfunctions of lowest order modes.

$\nu$	Eigenfrequency squared, ${\sigma }^2$	Eigenfunction, $W$
0	0	1
1	1	$\left(2n-1\right)\eta$
2	$\left(2n+1\right)/n$	$\left(n-\frac{1}{2}\right)\left[{\eta }^2\left(2n+1\right)-1\right]$
3	$3\left(n+1\right)/n$	$\frac{\eta }{6}\left(4n^2-1\right)\left[\left(2n+3\right){\eta }^2-3\right]$
4	$2\left(2n+3\right)/n$	$\frac{1}{24}\left(4n^2-1\right)\left[\left(2n+3\right){\eta }^2\left(\left(2n+5\right){\eta }^2-6\right)+3\right]$

Equations

$$L(r)=L_{\infty }{\left(1-\frac{2r_{\mathrm{g}}}{r}\right)}^{-1}.$$

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

$$\frac{1}{\rho }\frac{dp}{dr}=-\frac{GM}{r^2\left(1-r_{\mathrm{g}}/r\right)}\left[1-\lambda {\left(1-\frac{2r_{\mathrm{g}}}{r}\right)}^{-1/2}\right].$$

$$R_{\mathrm{ECS}}\equiv \frac{2r_{\mathrm{g}}}{(1-{\lambda }^2)},$$

$$p=K{\rho }^{1+1/n},$$

$$(n+1)\frac{p}{\rho }+\Phi (r)=\mathrm{const}.$$

$$(n+1)\frac{p}{\rho }+\Phi (r)=(n+1)\frac{p_0}{{\rho }_0}+\Phi (r_0),$$

$$\frac{p}{\rho }$$

$$f(r)$$

$$\Phi (r)-\Phi (r_0)$$

$$f$$

$${\omega }_{r,0}^2=\frac{(1-{\lambda }^2{)}^4}{16{\lambda }^4}\frac{c^2}{r_{\mathrm{g}}^2}$$

$$\Delta r=2\beta r_0/{\omega }_{r,0}=4\sqrt{2n}\left(\frac{c_{\mathrm{s},0}}{c}\right)\left(\frac{{\lambda }^2}{1-{\lambda }^2}\right)$$

$$p$$

$$\rho$$

$$m=\int \mathrm{d}m={\int }_{z_0}^H\rho (z)\mathrm{d}z.$$

$$w(t)=\int g\mathrm{d}m={\int }_{z_0}^{H}g(z+\xi )\rho (z)\mathrm{d}z.$$

$$w(t)=w_0+{\omega }_{\perp }^2{\int }_{z_0}^H\xi (z,t)\rho (z)\mathrm{d}z.$$

$$m{\partial }^2{\xi }_{\mathrm{C}M}/\partial t^2=p-w.$$

$${\partial }^2{\xi }_{\mathrm{C}M}/\partial t^2+(\gamma +1){\omega }_{\perp }^2{\xi }_{\mathrm{C}M}=0,$$

$$\frac{\partial }{\partial t}\delta v_r+\frac{\partial }{\partial r}\left(\frac{\delta p}{\rho }\right)$$

$$\frac{\partial }{\partial t}(\delta \rho )+\frac{1}{r^2}\frac{\partial }{\partial r}\left(r^2\rho \delta v_r\right)$$

$$\delta v_r$$

$${\omega }^{2}W$$

$$\frac{{\partial }^{2}W}{\partial r^2}+\left(\frac{n}{f}\frac{\partial f}{\partial r}+\frac{2}{r}\right)\frac{\partial W}{\partial r}+\frac{{\omega }^2}{c_s^2}W$$

$$(1-{\eta }^2)\frac{{\partial }^{2}W}{\partial {\eta }^2}-2n\eta \frac{\partial W}{\partial \eta }+2n{\sigma }^{2}W$$

$$(\Delta p{)}_{f=0}={\left(\delta p+\xi \frac{\partial p}{\partial r}\right)}_{f=0}=0,$$

$${\left[{\rho }_0{f}^n{\omega }_{r,0}\left(W\sigma -\frac{\eta }{\sigma }\frac{\partial W}{\partial \eta }\right)\right]}_{f=0}=0.$$

$${\omega }_{\nu }=\sqrt{\frac{\nu (\nu +2n-1)}{2n}}{\omega }_{r,0}.$$

$$\frac{{\omega }_b}{{\omega }_f}=\sqrt{2+\frac{1}{n}}=\sqrt{1+\gamma }.$$

Full text

††pubyear: 2016††pagerange: Radial modes of levitating atmospheres around Eddington-luminosity neutron stars–References

1 Introduction

Recent observations have mapped the existence of neutron stars with highly super-Eddington luminosities in ultra-luminous X-ray source (ULX) systems NGC 7793 P13, NGC 5907, NuSTAR J09551+6940.8 (Bachetti et al., 2014; Israel et al., 2016, 2017). It has long been known that neutron stars in X-ray binaries, A0538-66, SMC X-1 and GRO J1744–28 are unusually luminous (Skinner et al., 1982; Sazonov et al., 1997; Coe et al., 1981). Also, studies of X-ray bursts report that, during very energetic Type I X-ray bursts, the X-ray luminosity is close to the Eddington limit (Lewin et al., 1997). For neutron stars with such high luminosities, it is possible that close to the star’s surface, radiative forces become dominant over gravity. Under such conditions a “levitating atmosphere” may be formed, the properties of which were derived in the spherically symmetric case by Wielgus et al. (2015, 2016), and had previously been discussed for the approximately plane- parallel case in the context of near Eddington luminosity accretion discs by Fukue (1996), who termed them ‘photon floaters’. An additional argument of the observational relevance for these levitating atmospheres has been afforded in a new study, in which they have been shown to affect the propagation of light rays from the central source, thereby affecting the observed appearance of the central object (Rogers, 2017).

X-ray bursters exhibit high frequency (300 – 600 Hz) ‘burst oscillations’ (Strohmayer, 2001), which may be related to the Eddington-luminosity of the sources. In bright X-ray bursts, neutron star often reaches near Eddington-limit luminosities during which the dominant radiation pressure drives the photospheric layers off the neutron star; this phenomenon is termed as photospheric radius expansion (PRE). Oscillations in the case of PRE bursts are observed during the rise and decay phase of the burst (Muno et al., 2002). However, bursts without PRE also show evidence for oscillations (van Straaten et al., 2001). During the decay phase of the burst, frequency is observed to increase with time (Strohmayer & Bildsten, 2006). While the source of these oscillations is still open for debate, here we try to understand if the radial oscillations of the levitating atmospheres could be relevant to such phenomena.

In this paper, we consider a radial perturbation of a spherically symmetric, optically thin levitating atmosphere. The goal is to analytically investigate the eigenmodes of the oscillations of the atmosphere. Abarca & Kluźniak (2016) have found the fundamental radial mode of oscillations of an optically thin levitating atmosphere. They have shown that including the radiation drag stabilizes the atmosphere against the radial perturbations that could excite the fundamental radial mode. Here, we would like to investigate all possible modes of radial oscillations for this given system. Assuming the atmosphere to be geometrically thin, we neglect the relativistic corrections and perform linear perturbation analysis in the Newtonian approximation. We then determine all the possible radial eigenmode solutions of the system. Our calculations are performed by perturbing the levitating atmospheres of Wielgus et al. (2015), which have the form of spherical shells, but the calculations should be equally valid for the plane parallel case of the photon floaters of Fukue (1996), and have some application as well to the eigenmodes of other accretion structures with a reflection symmetry, such as the slender tori of Blaes et al. (2006).

Optically thin levitating atmospheres are especially likely to survive in the pulsating ULXs, which can be understood as slowly rotating neutron stars, a special case of extremely luminous X-ray pulsars (Bachetti et al., 2014), in which the super-Eddington accretion flow near the poles is confined by the strong stellar magnetic field (Basko & Sunyaev, 1976; Mushtukov et al., 2017).

A study of optically thick atmospheres would have been more directly relevant to X-ray burst oscillations, however we are hampered at present by the lack of a suitable analytic solution for the background equilibrium solution. The optically thick super-Eddington atmosphere solution (Wielgus et al., 2016) requires a proper treatment of the radiative transfer that has only been achieved numerically.

2 Levitating atmospheres

Abramowicz et al. (1990) have studied the fully general relativistic radial motion of test particles in combined gravitational and radiation fields around a spherical, compact star radiating isotropically. They have shown that under these circumstances, a spherical shell of matter may form around the neutron star whose location varies with luminosity from the surface of the star to infinity. This surface is defined by the radius where the gravitational forces balance the radiative forces due to stellar luminosity. Unlike in the Newtonian theory, where both the gravity and radiative flux fall off as $1/r^{2}$, in general relativity, the radiative flux has a stronger radial dependence than the gravitational acceleration. Therefore, radiative forces may dominate over gravitational pull close to the surface of the compact star, while the gravitational forces take over at larger radii. A similar situation may obtain above nearly Eddington-luminosity accretion disc, as discussed in several papers (Fukue, 1996; Kitabatake & Fukue, 2003; Itanishi & Fukue, 2017).

For a star radiating at super-Eddington luminosities, gravitational forces will balance the radiation forces on a sphere at a particular radius from the star. It has been referred to as the ‘Eddington capture sphere’ (ECS) (Stahl et al., 2012). A test particle moving with a reasonable velocity would bind to this sphere without falling on to the stellar surface also if initially in non-radial motion as the Poynting-Robertson-radiation drag removes the angular momentum of the particle reducing it to zero (Bini et al., 2009; Oh et al., 2010; Stahl et al., 2013). The captured particles levitate above the stellar surface, as every point on ECS corresponds to an equilibrium state that is stable in radial direction (Abramowicz et al., 1990) and neutral in directions tangent to the sphere (Stahl et al., 2012). These properties of ECS result in an atmosphere in hydrostatic equilibrium, which is detached from the star and levitating above its surface (Wielgus, 2016). In the case of an optically thin atmosphere, analytical solutions show that both the pressure and density fall steeply on either sides of the ECS, while the peak value is attained at the ECS (Wielgus et al., 2015). The term “levitation” was first introduced by Inogamov & Sunyaev (1999) in the context of the luminous equatorial accretion belt (boundary layer) formed at the interface of the neutron star surface and the accretion disc; the effective gravity in the boundary layer is further reduced by the rapid azimuthal motion of the fluid, so the balance considered by these authors was between gravity and three forces: pressure gradient, centrifugal and radiation forces. Our simplified treatment neglects rotation.

3 Hydrostatic Equilibrium Configuration

In the Schwarzschild metric around a neutron star of mass $M$, a local observer at a certain radius $r$ measures the redshifted luminosity of the star, $L(r)$ as

[TABLE]

where $L_{\infty}$ is the luminosity of the star measured at infinity and $r_{\mathrm{g}}=GM/c^{2}$ is the gravitational radius. For a compact star of radius $R_{*}$ and luminosity $L_{*}$, we have $L_{\infty}=L_{*}\left(1-\dfrac{2r_{\mathrm{g}}}{R_{*}}\right)$. For convenience, we shall use the following parametrized luminosity throughout our calculations,

[TABLE]

where $L_{\mathrm{Edd}}$ is the Eddington luminosity.

We now consider a static atmosphere levitating above the surface of the neutron star radiating at super-Eddington luminosities. The equation of hydrostatic equilibrium for an optically thin atmosphere under relativistic formalism has been derived previously (Wielgus et al., 2015):

[TABLE]

The first term on the right-hand side comes from the gravitational force, while the second term is related to the radiation force in an optically thin fluid. In deriving the above Eq. (3), it has been assumed that the atmospheric temperatures are non-relativistic, i.e., of the order of few keV. At such temperatures, one can neglect the internal energy and pressure of the fluid in comparison to the rest mass energy density.

From the above equation, it follows that the gravitational force and the force due to radiation balance each other at radius

[TABLE]

and the pressure gradient at this radius becomes zero, corresponding to the maximum pressure at $R_{\mathrm{ECS}}$.

We consider barotropic fluid, $p=p(\rho)$, and in this case we assume it to be a polytrope. Therefore

[TABLE]

where $K$ is the polytropic constant and $n$ is the polytropic index.

We may cast Eq. (3) in terms of an effective potential $\Phi(r)$, in such a way that upon its integration using Eq. (5), it results in the Bernoulli equation,

[TABLE]

We may evaluate this constant on the equilibrium surface, $r_{0}\equiv R_{\mathrm{ECS}}$, to get

[TABLE]

where subscript zero denotes the quantities evaluated at the ECS, $r=r_{0}$. Following simple algebraic steps, one could write

[TABLE]

where we define the function $f(r)$ as

[TABLE]

and $c_{\mathrm{s}}^{2}=\dfrac{(n+1)}{n}\dfrac{p}{\rho}$ is the square of the sound speed.

We now evaluate this function $f$ in the vicinity of $r_{0}$ by the second-order Taylor expansion of the effective potential, which is valid for thin atmospheres. We have

[TABLE]

where, we have introduced a new coordinate $x=(r-r_{0})/r_{0}$, centred at $r_{0}$. This will immediately give us

[TABLE]

where we define $\beta^{2}=2nc_{\mathrm{s,0}}^{2}/r_{0}^{2}$ and

[TABLE]

for our convenience. We shall later show that $\omega_{r,0}$ is indeed the frequency of the fundamental mode. The boundaries of the atmosphere are defined at $f=0$ as both the pressure and density vanish at this point. Following this, the thickness of the atmosphere would be

[TABLE]

In Fig. 1, we explore the dependence of the atmosphere thickness, $\Delta r$ on $\beta$ and $\lambda$. Each solid and dash-dot line, representing a particular value of $\lambda$, shows that $\Delta r$ increases with $\beta$. Specifically, solid lines represent the regime where the second order approximation may hold good, i.e., $x\ll 1$ (we simply chose $x=0.005$ for the plot). Therefore, for a star with a given luminosity, $\lambda$ remains constant, and the extent to which the atmosphere extends is determined by $\beta$ and in the slender limit $\beta\rightarrow 0$.

As $\beta$ implicitly depends on temperature at $r_{0}$, in Fig. 1 we considered the temperature range within $[10^{5},10^{8}]~{}\mathrm{K}$ and we take the mean molecular weight of the fluid, $\mu$ to be $1/2$. The dashed line in Fig. 1 is plotted to give a glimpse of how $\Delta r$ varies with $\lambda$. As shown in the figure, the slenderness of the atmosphere increases as $\lambda$ decreases. Here the dashed line is plotted for the temperature value, $10^{8}~{}\mathrm{K}$, but the trend would remain the same for other values of temperature.

Based on Eq. (8) and the equation of state, Eq. (5), it is convenient to scale the pressure and density in equilibrium as follows:

[TABLE]

which clearly ensures that the pressure and density drop to zero at the boundaries while retaining their maximum value at the centre. The density gradient is always in the direction of the effective gravity (taking into account radiation forces as well), and so is Rayleigh-Taylor stable.

In passing, we note an interesting property of the numerical solution for the optically thick atmosphere, which may make our analysis an interesting limiting case of future (perhaps numerical) calculations of the oscillations of the optically thick levitating atmosphere. The optically thin edges of the optically thick atmosphere coincide with this analytically derived optically thin solution, while the radial density has a universal profile (independent of the mass of the atmospheric shell) in the optically thick regions, satisfying over the bulk of the shell the condition of density decreasing with radial distance from the source, i.e., increasing in the direction of the effective force of gravity (Wielgus, 2016).

4 Elementary discussion of radial oscillations

To develop intuitions relevant to fluid motion in a gravitational field of magnitude increasing with height we present an elementary discussion in a Newtonian formalism of the lowest frequency radial modes of a thin levitating atmosphere in the plane-parallel approximation.

Choose the $z$ axis in the vertical direction, pointing either towards or away from the neutron star, and consider a fluid atmosphere in a gravitational field directed towards the $z=0$ plane: $\vec{g}(z)=-\hat{z}g(z)$, with $g(z)=\omega^{2}_{\perp}z$. In hydrostatic equilibrium, the surface of maximum pressure coincides with the $z=0$ plane, and the atmosphere is symmetric with respect to reflection in that plane.

Let us consider the lowest-frequency modes of the atmosphere, one with a symmetric displacement function $z\mapsto z+\xi(z,t)$, $\xi(-z,t)=\xi(z,t)$, and no nodes (the “fundamental"), and one with antisymmetric displacements, $\xi(-z,t)=-\xi(z,t)$, and one node (the “first overtone"). As we will see, the fundamental mode corresponds to rigid motion, with the plane of non-varying maximum pressure moving harmonically up and down with respect to the $z=0$ plane. The first overtone is a breathing-like motion characterized by reflection symmetry in the $z=0$ plane. In general, as the atmosphere oscillates, its pressure is a function of position and time $p(z,t)$.

To derive the eigenfrequencies, consider — without any loss of generality — the part of the atmosphere that is initially above the horizontal plane at $z=z_{0}>0$ plane. For a polytropic atmosphere the density goes to zero at the upper boundary initially at, say, $z=H>0$ (and another one at $z=-H$). Thus, we are focusing on atmospheric slab of matter confined by the upper boundary (plane) at $H+\xi(H,t)$ and an imagined lower boundary at $z_{0}+\xi(z_{0},t)$. The conserved mass per unit area in this slab is

[TABLE]

The slab suffers two external forces, a body force owing to gravity, and a surface force at the lower boundary at $z_{0}+\xi$. The latter force per unit area is just the pressure $p=p_{0}+\Delta p(t)$, with $p_{0}$ the equilibrium pressure at $z=z_{0}$, and $\Delta p$ the Lagrangian change in pressure: $\Delta p(t)=p(z_{0}+\xi_{0},t)-p_{0}$, with $\xi_{0}\equiv\xi(z_{0},t)$. The body force is just the weight of the slab, per unit area

[TABLE]

We can now specify to the case at hand, $w_{0}=m\omega^{2}_{\perp}z_{\mathrm{C}M}$ with $z_{\mathrm{C}M}=(\int z\,\mathrm{d}m)/m$,

[TABLE]

The equation of motion for the displacement of the centre of mass, $\xi_{\mathrm{C}M}(t)\equiv\xi(z_{\mathrm{C}M},t)$, is

[TABLE]

For a uniform displacement, $\xi(z,t)=\xi_{\mathrm{C}M}(t)$, i.e., for a rigid motion, there are no nodes, and we obtain the “fundamental" frequency $\omega_{1}=\omega_{\perp}$ from the simple harmonic equation $\partial^{2}\xi_{\mathrm{C}M}/\partial t^{2}=-\omega^{2}_{\perp}\xi_{\mathrm{C}M}$, following from $\partial\xi/\partial z=0$ that implies $p=p_{0}=w_{0}$, with $p-w=-m\omega^{2}_{\perp}\xi_{\mathrm{C}M}$ by Eqs. (15) and (17).

For the breathing mode, with the (easy to derive) ansatz $\xi=z\xi_{\mathrm{C}M}/z_{\mathrm{C}M}$, we have $\Delta p=-\gamma p_{0}\xi_{\mathrm{C}M}/z_{\mathrm{C}M}=-\gamma m\omega^{2}_{\perp}\xi_{\mathrm{C}M}$, and $w(t)=w_{0}+m\omega^{2}_{\perp}\xi_{\mathrm{C}M}(t)$. Thus

[TABLE]

and the eigenfrequency of the breathing mode is $\omega_{2}=\sqrt{1+\gamma}\omega_{1}$. The same ratio of $\sqrt{1+\gamma}$ can be obtained for the breathing and vertical mode eigenfrequencies of a slender torus orbiting a black hole (Blaes et al., 2006), if the Newtonian limit is taken (in which the epicyclic frequencies are equal to one another, $\omega_{r}=\omega_{\perp}$). In general, just like the fundamental mode just considered, the vertical mode of the slender torus corresponds to rigid motion at the vertical epicyclic frequency, $\omega_{1}=\omega_{\perp}$ (Kluźniak, 2005), and the breathing mode of the levitating atmosphere is a direct analogue of the breathing mode of a slender torus.

5 Radial oscillations in Newtonian framework

We subject the static atmosphere to small radial perturbations of the form $\delta X(r,t)=\delta X_{*}(r)e^{i\omega t}$, where $\omega$ is the frequency with which perturbations propagate. Because of these small radial perturbations, fluid attains certain radial velocity given by $\delta v_{r}$. We then perform the linear perturbation analysis and derive the wave equation that governs the propagation of these perturbations. To do so, we assume that atmosphere is geometrically thin, i.e, $x\ll 1$. Relativistic corrections to energy density are negligible on such small scales to which these atmospheres extend. Therefore, Newtonian theory would suffice to evaluate $\omega$. In this work, we ignore the radiation drag, but note that Eq. (3) is valid in any case as the equilibrium atmosphere is static and the initial velocities being zero do not contribute to the pressure balance.

We start with the equations of motion for a radially moving fluid in Newtonian framework. Linearized Euler and continuity equations are given as

[TABLE]

We define a scalar potential $W=\delta p/(\omega\rho)$, which identifies with the fluid enthalpy, and simplifies our problem to great extent (Papaloizou & Pringle, 1984). Considering the time derivative of the perturbation quantities would give us the following relations

[TABLE]

Combining both of these equations, Eq. (21) and Eq. (22) and using Eq. (14), we obtain a second-order differential equation in $W$,

[TABLE]

In general it is difficult to solve the above equation analytically but in the slender limit $(\beta\rightarrow 0)$ this equation simplifies to

[TABLE]

where we have made use of the notations $\eta=x\omega_{r,0}/\beta$, with $f=1-\eta^{2}$, and $\sigma=\omega/\omega_{r,0}$. The eigenvalue problem given by Eq. (24) has to be solved with appropriate boundary conditions. In this case, it is required that pressure and density vanish at the boundaries. Since we consider small, linear perturbations, this condition translates to vanishing of the Lagrangian perturbation in pressure at the boundaries in the form

[TABLE]

where $\xi$ is the Lagrangian displacement. After rewriting the Lagrangian displacements in terms of $W$, the surface boundary condition becomes,

[TABLE]

The requirement that $W$ and its first derivative be finite at $f=0$ is a sufficient condition for this boundary condition to hold. It is also possible that the singular eigenfunctions may appear on the surface boundaries. But such solutions must be rejected as they are not compatible with Eq. (24). Therefore, the eigenfunctions are simply given by the Gegenbauer polynomials $C^{\mu}_{\nu}(\eta)$ with order $\mu=n-1/2$ and degree $\nu=0,1,2,3...$, where $\nu(\nu+2n-1)=2n\sigma^{2}$ (Abramowitz & Stegun, 1972). Hence,

[TABLE]

with $\omega_{r,0}^{2}=(1-\lambda^{2})^{4}c^{6}/(16G^{2}M^{2}\lambda^{4})$ by Eq. (12).

6 Results and Discussions

When the atmospheric fluid element in equilibrium is radially displaced, the restoring forces due the gravity and pressure variation act upon the fluid trying to restore it back to equilibrium. But owing to inertia, the displacement overshoots and thus the atmosphere is set into harmonic oscillations. Table LABEL:table1 summarizes the eigenfunctions and eigenfrequencies of the first few eigenmodes of this system. Degree $\nu=0$ corresponds to a trivial solution with eigenfunction $1$, i.e., vanishing displacement by Eq. (21).

The lowest order mode, $\nu=1$, is the fundamental mode given by the frequency $\omega_{\mathrm{f}}\equiv\omega_{1}=\omega_{r,0}$. In this mode the fluid oscillates around the ECS radially back and forth while maintaining the Lagrangian pressure profile in the atmosphere constant. Thus it is an incompressible mode. The displaced fluid is set into oscillations by the restoring force due to gravity.

For the next lowest order mode, $\nu=2$, pressure restoring forces become important. We term this as the breathing mode since the atmosphere expands and contracts symmetrically with respect to the ECS, with a frequency $\omega_{\mathrm{b}}\equiv\omega_{2}=\sqrt{(2n+1)/n}\,\omega_{r,0}$.

In Fig. 2, we plot the radial profiles of the velocity perturbation. We consider $n=3/2$ and $\lambda=0.85$. We can scale the perturbation by an arbitrary value so that assumptions made to pursue linear perturbation analysis hold well. The domain of the plot corresponds to the extent of the equilibrium atmosphere. The vertical dotted line marks the location of the ECS, which for the given parameters is $\approx 7.207~{}r_{\mathrm{g}}$. The horizontal dashed line is plotted to identify the nodes at which velocity perturbations vanish. The fluid in the fundamental mode oscillates with uniform velocity depicted by the horizontal blue line. This means that at a given time, all the fluid particles in the atmosphere are moving with same speed in same direction. As the system evolves in this mode, the fluid oscillates back and forth through the horizontal dashed line. Velocity perturbation for the breathing mode holds a linear relation with radius as shown by the green line with a node at the ECS. While expanding and contracting, fluid in the atmosphere either recedes away from the ECS or approaches ECS, while the centre stays still marking a node at the ECS. As the system evolves, the green curve pivots around its node at the ECS. As we move further to the higher modes, the number of nodes simply increases.

In the last figure, Fig. 3, the top and bottom panel show the radial profiles of the Eulerian perturbation in pressure and density respectively. The parameters considered here are same as before. The vertical dotted line carries the same meaning as discussed for the previous figure. The horizontal dashed line is plotted to identify the nodes at which the perturbations in density and pressure vanish. Note that the pressure and density perturbations vanish at the boundaries too as required by the boundary condition, but they are not to be counted as nodes. The fundamental mode plotted by the dark blue curve is an incompressible mode. But to see this from an Eulerian perturbation plot, we need the information of velocity perturbations. Since the velocity perturbation in the fundamental mode is uniform (see Fig. 3), its divergence vanishes. This implies from the continuity equation that density or the pressure in a given volume remains the same. For the breathing mode, the pressure drop at the centre is accompanied by an increase in pressure close to the boundaries as shown by the green curve with two nodes in the plot. Since we considered a polytropic fluid, the density perturbation profiles look analogous to the pressure perturbation profiles. Note that the number of nodes in pressure perturbations profile is one greater than the number of nodes in the velocity perturbations for a given mode.

From Table LABEL:table1, the frequency of oscillations in the fundamental mode is simply $\omega_{r,0}$. This is independent of $n$ because the fundamental mode oscillations are caused purely by gravity. Remaining modes have the contribution from the fluid’s pressure. The exact analytical expressions for eigenfrequencies allows us to calculate the ratio of the frequencies of the breathing mode, $\omega_{\mathrm{b}}$ to the fundamental mode, $\omega_{\mathrm{f}}$,

[TABLE]

For $n=3/2$, this ratio is $\sqrt{8/3}\approx 1.633$. Interestingly, the ratio of frequencies of any two modes is independent of $\lambda$, the stellar luminosity.

For a neutron star with mass $1.4M_{\odot}$ and $\lambda=0.8$, we calculate the Newtonian frequency of oscillations for the $\nu=1$ fundamental mode to be $1.171~{}\mathrm{kHz}$. Our result matches with the fundamental mode oscillation frequency earlier derived for full general relativistic equations (Abarca & Kluźniak, 2016), if we consider the following relativistic effects. First, the invariant line element differs from coordinate displacement $\mathrm{d}r$, inducing a correction factor, $1/\sqrt{g_{rr}}\approx\lambda$. Secondly, in general relativity, proper time and coordinate time are different and the frequency measured by an observer at infinity would be redshifted by a factor $\sqrt{g_{tt}}\approx\lambda$. Owing to these two factors, if we multiply the Newtonian frequency by $\lambda^{2}$ we get the relativistic frequency $\tilde{\omega}=\dfrac{(1-\lambda^{2})^{2}c}{4r_{\mathrm{g}}}$. For the above considered parameters, the relativistic frequency for the fundamental mode is then $\approx 750~{}\mathrm{Hz}$, for $\lambda=0.8$, which is the same as obtained in Abarca & Kluźniak (2016). For $0.8<\lambda$ the fundamental frequency decreases in value with increasing luminosity, but similarly accounting for the relativistic corrections, the eigenfrequencies of the higher modes for $0.8<\lambda<0.98$ may fall in the observed range of 300–600 Hz for X-ray burst quasi-periodic oscillations (QPOs). Abarca & Kluźniak (2016) have also shown that for optically thin solutions, radiation drag can efficiently damp the fundamental mode oscillation. Here we find that the higher modes have larger frequencies than the damping rate found by those authors, so it is likely that the oscillations of the higher modes can be underdamped. We leave this for our future work, where we plan to study the full general relativistic motion of particles including radiation drag.

Also, note that $\omega_{r,0}$ increases with decreasing luminosity111For an impulsive change in luminosity the optically thin atmosphere may be ejected (Kluźniak, 2013; Mishra & Kluźniak, 2014).. So for any given mode, assuming that the polytropic index of the atmosphere remains the same, frequency of any of the eigenmodes would increase with time if the luminosity of the source decreases with time. This is similar to what is usually observed in the decay phase of the X-ray bursts (Strohmayer, 2001; Strohmayer & Bildsten, 2006).

7 Conclusions

In this paper, we provide a complete Newtonian analytic solution for the radial oscillation eigenmodes of a thin levitating atmospheric shell centred at the ECS around a luminous non-rotating stellar source. We derive the complete set of eigenfunctions for the modes, and the eigenfrequencies for all the mode oscillations222The atmospheric oscillations described here are very different from the oscillations of the ECS described in (Wielgus et al., 2012), which are due to accretion feedback and stellar luminosity variability. are given by an exact analytical expression, Eq. (27). We provide the ratio of the frequencies of breathing mode to the fundamental mode, $\sqrt{1+\gamma}$, which entirely depends on the polytropic index of the atmosphere. For $\gamma=5/3$, we find this ratio to be $\sqrt{8/3}\approx 1.633$, and for $\gamma=4/3$ it is $\sqrt{7/3}\approx 1.53$.

Acknowledgements

We thank Maciek Wielgus and David Abarca for useful discussions throughout the work. This research was supported by the Polish NCN grants 2013/08/A/ST9/00795 and 2013/10/M/ST9/00729 and also in part by the National Science Foundation under Grant No. NSF PHY11-25915.