We consider massive STTs described by the action in the Einstein frame (G = c = 1) [70]. S [ g μ ν , ϕ ] = 1 16 π ∫ d 4 x − g R − 2 ∂ μ ϕ ∂ μ ϕ − V ( ϕ ) + L M ( A 2 ( ϕ ) g μ ν , χ ) , (1) with the metric g _μν , the curvature scalar R, the scalar field ϕ, the scalar potential V(ϕ), the nuclear matter action L _M , and the standard Brans-Dicke coupling function A ( ϕ ) = e − 1 3 ϕ . (2)

We choose two examples for the potential V(ϕ) [71]: V _I represents simply a mass term, and V _II is related to R ² gravity [28,63,72,73]. V I = 2 m ϕ 2 ϕ 2 , (3) V I I = 3 m ϕ 2 2 1 − e − 2 ϕ 3 2 . (4)

To see the relation with R ² gravity, we recall its action in the Jordan frame S [ g μ ν * ] = 1 16 π ∫ d 4 x − g * R * + a R * 2 + L M ( g μ ν * , χ ) , (5) where R* is the Ricci scalar associated with the metric g μ ν * . The parameter a is a positive constant with units of [length]², and controls the strength of the R ² deviation from General Relativity. This theory can be recast into a particular Brans-Dicke STT with Jordan frame action [63,65]. S [ g μ ν * , ψ ] = 1 16 π ∫ d 4 x − g * ψ R * − U ( ψ ) + L M ( g μ ν * , χ ) , (6) with scalar field ψ, and potential U(ψ) U ( ψ ) = R * ψ − f ( R * ) = 1 4 a ( ψ − 1 ) 2 , ψ = d f d R * . (7)

The transition to the Einstein frame follows with help of the relations g μ ν * = A 2 g μ ν ⇒ R * = A − 2 R − 6 A − 3 ∇ μ ∂ μ A , A − 2 = ψ = e 2 3 ϕ . (8) In the Einstein frame with the new metric g _μν and the scalar field ϕ we then obtain the action (1) with the potential V ( ϕ ) = 1 4 a 1 − e − 2 ϕ 3 2 . (9) The potentials in the two frames are related by U ( ψ ) = A − 4 V ( ϕ ) . (10) This action (1) now has an explicit kinetic term for the scalar field. The parameter a is related to the mass of the scalar, m ϕ = 1 6 a . (11) Note that, when the scalar field is weak, the potential is essentially quadratic, given by V ( ϕ ) ∼ ϕ 2 3 a . When a is set to zero, General Relativity is recovered with an infinitely massive scalar field, that is, hence, suppressed to zero. In the action in the Einstein frame (1) the matter and scalar field are non-minimally coupled. In contrast, in the Jordan frame (6) the scalar field is non-minimally coupled to the Ricci scalar, which makes the action highly non-linear. For further discussions on the Einstein and Jordan frames in STTs, see e.g., [72,73].

In this paper, we are working in the Einstein frame, and thus with the action (1). The field equation for the metric g _μν is then given by, G μ ν = T μ ν ( S ) + 8 π T μ ν ( M ) − 1 2 V ( ϕ ) g μ ν , (12) with Einstein tensor G μ ν = R μ ν − 1 2 R g μ ν . The energy-momentum tensor of the scalar field T ^(S) is given by T μ ν ( S ) = 2 ∂ μ ϕ ∂ ν ϕ − g μ ν ∂ σ ϕ ∂ σ ϕ , (13) and the energy-momentum tensor of the matter T ^(M), which is assumed to be a perfect fluid, is T μ ν ( M ) = ( ρ + p ) u μ u ν + p g μ ν , (14) where ρ is the energy density, p is the pressure, and u _μ is the 4-velocity of the fluid. We assume the existence of a barotropic equation of state relating the energy density and the pressure, determined by the properties of matter at high densities. Since this relation is typically calculated in the equivalent (physical) Jordan frame, we need to transform the energy density and the pressure from the Jordan frame to the Einstein frame. The relations with the pressure p ̂ and the density ρ ̂ defined in the Jordan frame are p = A 4 p ̂ , ρ = A 4 ρ ̂ , (15) and the equation of state is a relation of the form ρ ̂ = ρ ̂ ( p ̂ ) . The field equation for the scalar field in (1) is ∇ μ ∇ μ ϕ = − 4 π 1 A d A d ϕ T ( M ) + 1 4 d V d ϕ , (16) where T ^(M) is the trace of the energy-momentum tensor of the matter.

For static and spherically symmetric neutron stars we consider the following Ansatz for the metric d s 2 = g μ ν ( 0 ) d x μ d x ν = − e 2 ν ( r ) d t 2 + e 2 λ ( r ) d r 2 + r 2 ( d θ 2 + sin 2 ⁡ θ d φ 2 ) . (17) The scalar field, energy density and pressure are simply given by ϕ = ϕ ₀(r), ρ ̂ = ρ ̂ 0 ( r ) and p ̂ = p ̂ 0 ( r ) , respectively, and the four-velocity of the static fluid is given by u ⁽⁰⁾ = −e ^ν dt.

Inside the star, the equations for the static functions are 1 r 2 d d r r ( 1 − e − 2 λ ) = 8 π A 0 4 ρ ̂ 0 + e − 2 λ d ϕ 0 d r 2 + 1 2 V 0 , (18) 2 r e − 2 λ d ν d r − 1 r 2 ( 1 − e − 2 λ ) = 8 π A 0 4 p ̂ 0 + e − 2 λ d ϕ 0 d r 2 − 1 2 V 0 , (19) d p ̂ 0 d r = − ( ρ ̂ 0 + p ̂ 0 ) d ν d r + 1 A 0 d A 0 d ϕ 0 d ϕ 0 d r , (20) d 2 ϕ 0 d r 2 + d ν d r − d λ d r + 2 r d ϕ 0 d r = 4 π 1 A 0 d A 0 d ϕ 0 A 0 4 ( ρ ̂ 0 − 3 p ̂ 0 ) e 2 λ + 1 4 d V 0 d ϕ 0 e 2 λ , (21) where A ₀ = A(ϕ ₀) and V ₀ = V(ϕ ₀).

The second order differential Eq. 21 for ϕ is obtained from the scalar field Eq. 16 at the static level, while the others come from the Einstein Eq. 12. The system of equations has to be complemented with an equation of state ρ ̂ 0 = ρ ̂ 0 ( p 0 ̂ ) . Note that the static Einstein frame quantities p 0 = A ( ϕ 0 ) 4 p ̂ 0 and ρ 0 = A ( ϕ 0 ) 4 ρ ̂ 0 , will appear in some formulas below (to simplify notation).

Solutions describing neutron stars have to satisfy certain boundary conditions. At the center of the star, the configuration has to be regular. An expansion yields λ = 4 3 π A c 4 ρ ̂ c + 1 12 V c r 2 + o ( r 4 ) , (22) ν = ν c + 2 3 π A c 4 ( 3 p ̂ c + ρ ̂ c ) − 1 12 V c r 2 + o ( r 4 ) , (23) p ̂ = p ̂ c + ( p ̂ c + ρ ̂ c ) 2 π A c 3 3 3 d A 0 d ϕ 0 | c ( 3 p ̂ c − ρ ̂ c ) + 3 72 d V 0 d ϕ 0 | c + V c 12 − 2 π 3 A c 4 ( 3 p ̂ c + ρ ̂ c ) r 2 + o ( r 4 ) , (24) ϕ 0 = ϕ c + 2 π 3 A c 3 d A 0 d ϕ 0 c ( 3 p ̂ c − ρ ̂ c ) + 1 24 d V 0 d ϕ 0 c r 2 + o ( r 4 ) , (25) where ν _c = ν(0), p ̂ c = p ̂ ( 0 ) , ϕ _c = ϕ ₀(0), V _c = V ₀(ϕ _c ), A _c = A ₀(ϕ _c ), etc. This expansion has only three free parameters ν _c , p ̂ c and ϕ _c , the others being determined by the couplings and the equation of state. However, a global solution that is asymptotically flat will have only one free parameter at the center (typically chosen to be the central pressure p c ̂ in the numerical calculations).

The border of the neutron star is defined by the point r = r _s , where the pressure vanishes, p ̂ 0 ( r s ) = 0 . The physical (Jordan frame) radius of the star is given by R _s = A(ϕ ₀(r _s ))r _s , see [68] for more details. Outside the neutron star there is no fluid, thus p ̂ 0 = ρ ̂ 0 = 0 in Eqs 18–21. We are interested in asymptotically flat solutions with vanishing scalar field at infinity. However, for a ≠ 0, the scalar field outside the star does not vanish. Nonetheless, sufficiently far from the star, the scalar field decays exponentially with ϕ ∼ 1 r e − m ϕ r , and the background metric is essentially given by the Schwarzschild solution, with e ^2ν = e ^−2λ ∼ 1 − 2M/r, where M is the total mass of the neutron star.

The mass–radius relation for most of the background neutron star configurations investigated in the following is exhibited in Figure 1, where the mass M is given in solar masses M _⊙ and the physical (Jordan frame) radius R _s in km. The chosen potential is V _II (R ² gravity), and the scalar field mass m _ϕ assumes several values in the physically interesting mass range [74,75]. The color coding is purple, cyan, red, orange and green for m _ϕ = 1.08, 0.343, 0.108, 0.0343 and 0.0108 neV, respectively. Also shown are the general relativistic limit (black), and the case of a massless scalar field (blue). For the SLy equation of state the maximum mass of static neutron stars in General Relativity is slightly above two solar masses. In the STTs studied, the value of the maximum mass is increased the more, the smaller the scalar field mass is, reaching about 2.25 solar masses. At the same time, the neutron star radii become larger except for rather small neutron star masses.

In this section we present the quasinormal mode formalism, focusing on polar perturbations. In the Einstein frame, we perturb the background metric, the scalar field and the fluid in the following way: g μ ν = g μ ν ( 0 ) ( r ) + ϵ h μ ν ( t , r , θ , φ ) , (26) ϕ = ϕ 0 ( r ) + ϵ δ ϕ ( t , r , θ , φ ) , (27) ρ = ρ 0 ( r ) + ϵ δ ρ ( t , r , θ , φ ) , (28) p = p 0 ( r ) + ϵ δ p ( t , r , θ , φ ) , (29) u μ = u μ ( 0 ) ( r ) + ϵ δ u μ ( t , r , θ , φ ) , (30) where ϵ ≪ 1 is the perturbation parameter. The zeroth order of the configurations describes a static and spherically symmetric background, discussed in the previous section, while the perturbations are in general dependent on time, radial coordinate and angular directions.

Assuming that the perturbation functions can be expanded as a product of radial, temporal and angular components, they can be further separated into classes of axial and polar perturbations, depending on the transformation of the angular component under parity [2,76–89]. For polar perturbations, the spherical harmonics transform as Y _lm (θ, φ) → Y _lm (π −θ, π + φ) = (−1) ^l Y _lm (θ, φ). For a study of the current theory in the axial case we refer the reader to [68].

We follow the standard gauge convention previously used for neutron stars in General Relativity [78,83]. The corresponding Ansatz for the polar perturbations of the metric is h μ ν ( polar ) = ∑ l , m ∫ r l e 2 ν H 0 Y l m − i ω r l + 1 H 1 Y l m 0 0 − i ω r l + 1 H 1 Y l m r l e 2 λ H 2 Y l m 0 0 0 0 r l + 2 K Y l m 0 0 0 0 r l + 2 ⁡ sin 2 ⁡ θ K Y l m e − i ω t d ω , (31) in the order of (t, r, θ, φ) in the rows and columns of the matrix. The functions H ₀, H ₁, H ₂, K only depend on the radial coordinate r, the integer multipole numbers l, m, and the complex wave frequency ω, where ω = ω _R + iω _I for ω R , ω I ∈ R . The functions Y _lm are the standard spherical harmonics. The scalar field can be decomposed as δ ϕ = ∑ l , m ∫ r l ϕ 1 Y l m e − i ω t d ω , (32) where again the function ϕ ₁ depends on r, l, m and ω.

Concerning the perturbation of the fluid inside the star, the scalar quantities, such as the energy density ρ and pressure p, can be decomposed similarly to the scalar field δ ρ = ∑ l , m ∫ r l E 1 Y l m e − i ω t d ω , δ p = ∑ l , m ∫ r l Π 1 Y l m e − i ω t d ω , (33) while the perturbation of the 4-velocity is given by δ u μ = ∑ l , m ∫ 1 2 r l e ν H 0 Y l m r l i ω e − ν e λ W / r − r H 1 Y l m − i ω r l e − ν V ∂ θ Y l m − i ω r l e − ν V ∂ ϕ Y l m e − i ω t d ω , (34) The functions V and W depend on r, l, m, ω.

Before continuing, we note that although we have given the perturbations in the Einstein frame, they can be alternatively defined in the Jordan frame. Perturbations in the Jordan frame for the scalar and the metric (δψ and h μ ν * , respectively) can be expressed as a combination of the perturbations in the Einstein frame, δ ψ = 2 3 e 2 3 ϕ 0 δ ϕ , h μ ν * = e − 2 3 ϕ 0 h μ ν − 2 3 g μ ν ( 0 ) δ ϕ . (35) Concerning the energy density and pressure, the barotropic equation of state in the Jordan frame implies a relation between the perturbations δp, δρ and δϕ δ ρ = d ρ ̂ 0 d p ̂ 0 δ p + 4 A 0 3 d A 0 d ϕ 0 ρ ̂ 0 − p ̂ 0 d ρ ̂ 0 d p ̂ 0 δ ϕ . (36)

Note that the complex wave frequency ω is the same in the Jordan frame and in the Einstein frame.

Employing the Ansatz shown in the previous section for the perturbations on the field Eq. 12 and Eq. 16, leads to a system of ordinary differential equations in r, that is characterized by the eigenvalue ω, and the multipole number l, but that is independent of m because of the spherical symmetry.

The modified Einstein Eq. 12 result in six ordinary differential equations d 2 K d r 2 = − 8 π e 2 λ E 1 + 1 2 r 2 2 d ϕ 0 d r 2 r 2 − 4 r d λ d r + ( l e 2 λ + 2 ) ( l + 1 ) H 2 + 2 d ϕ 0 d r d ϕ 1 d r + 1 6 r 3 e 2 λ r a e − 4 3 ϕ 0 − e − 2 3 ϕ 0 − 12 l d ϕ 0 d r ϕ 1 + 1 2 r 2 2 l r d λ d r + e 2 λ ( l ( l + 1 ) − 2 ) − 2 l ( l + 2 ) K + r d λ d r − 2 l − 3 r d K d r + 1 r d H 2 d r + 1 8 r 2 64 π ρ r 2 e 2 λ + 8 d ϕ 0 d r 2 r 2 + e 2 λ r 2 a e − 2 3 ϕ 0 − 1 2 − 16 r d λ d r + 8 ( 1 − e 2 λ ) H 0 , (37) d K d r = 1 r H 2 + d ν d r r − l − 1 r K − 2 d ϕ 0 d r ϕ 1 − 8 π e λ p + ρ r W + 1 8 r 4 l ( l + 1 ) − 16 d λ d r r e − 2 λ + 8 ( e − 2 λ − 1 ) H 1 + 1 8 r 64 π ρ r 2 − r 2 a e − 2 3 ϕ 0 − 1 2 + 8 d ϕ 0 d r 2 r 2 e − 2 λ H 1 , (38) d H 1 d r = e 2 λ r H 2 + d λ d r − d ν d r − l + 1 r H 1 + e 2 λ r K − 16 π r e 2 λ p + ρ V , (39) d H 0 d r = − 8 π e 2 λ r Π 1 − 2 r ω 2 e − 2 ν H 1 + l 2 r ( l + 1 ) e 2 λ − 2 H 0 − 2 r d ϕ 0 d r d ϕ 1 d r − 1 8 r 64 π e 2 λ p r 2 − r 2 a e 2 λ e − 2 3 ϕ 0 − 1 2 + 8 e 2 λ H 2 + 1 2 r 2 e 2 ( λ − ν ) ω 2 r 2 − e 2 λ ( l − 1 ) ( l + 2 ) + 2 l r d ν d r + 1 K + 1 2 3 r a e 2 λ e − 2 3 ϕ 0 − e − 4 3 ϕ 0 − 2 l d ϕ 0 d r ϕ 1 + r d ν d r + 1 d K d r , (40) d H 0 d r = − d ν d r + 1 r H 2 − r ω 2 e − 2 ν H 1 − d ν d r + l − 1 r H 0 + l r K + 4 d ϕ 0 d r ϕ 1 + d K d r , (41) d 2 K d r 2 − d 2 H 0 d r 2 = 16 π e 2 λ Π 1 − 4 d ϕ 0 d r d ϕ 1 d r + d ν d r + 1 r d H 2 d r + 2 r ω 2 e − 2 ν d H 1 d r − 2 e − 2 ν ω 2 r d λ d r − l − 2 H 1 − 1 2 r 2 l 2 r d λ d r − 4 r d ν d r + l e 2 λ + e 2 λ − 2 l H 0 + 1 2 r 2 − 2 e − 2 ν + 2 λ ω 2 r 2 + 4 d ϕ 0 d r 2 r 2 − 4 r 2 d λ d r d ν d r + 4 d ν d r 2 r 2 + 4 r 2 d 2 ν d r 2 H 2 + 1 2 r 2 + l 2 e 2 λ + 2 l r d ν d r − 4 r d λ d r + l e 2 λ + 4 r d ν d r + 2 l H 2 + 1 4 r 2 64 π e 2 λ p r 2 − 4 e − 2 ν + 2 λ ω 2 r 2 + 8 r 2 d λ d r d ν d r − 8 d ϕ 0 d r 2 r 2 − 8 d ν d r 2 r 2 K + 1 4 r 2 4 r ( l + 2 ) d λ d r − d ν d r − 8 r 2 d 2 ν d r 2 − 4 l ( l + 1 ) − r 2 a e 2 λ 1 − e − 2 3 ϕ 0 2 K + 1 3 r 3 r a e 2 λ e − 4 3 ϕ 0 − e − 2 3 ϕ 0 − 12 l d ϕ 0 d r ϕ 1 + 1 r ( 2 l + 1 ) − d λ d r + 2 d ν d r d H 0 d r + d λ d r − d ν d r − 2 r ( l + 1 ) d K d r . (42)

In addition, it imposes the constraint H 2 = H 0 . (43) The scalar field Eq. 16 results in d 2 ϕ 1 d r 2 = − 4 3 π e 2 λ ( E 1 − 3 Π 1 ) + l 2 r d ϕ 0 d r H 0 + d ϕ 0 d r e − 2 ν ω 2 r H 1 − l r d ϕ 0 d r K + 1 6 r 2 − e 2 λ r 2 6 e − 2 ν ω 2 + e − 2 3 ϕ 0 − 2 e − 4 3 ϕ 0 1 a + 6 l r d λ d r − d ν d r + 6 l ( 1 + l ) ( e 2 λ − 1 ) ϕ 1 + 1 2 r − 2 r d ϕ 0 d r d λ d r + 2 r d ϕ 0 d r d ν d r + 2 r d 2 ϕ 0 d r 2 + d ϕ 0 d r l + 4 d ϕ 0 d r H 2 + d λ d r − d ν d r − 2 r ( l + 1 ) d ϕ 1 d r + 1 2 d ϕ 0 d r d H 2 d r + d H 0 d r − 2 d K d r . (44)

The barotropic condition on the equation of state (Eq. 36) becomes a relation between the energy density perturbation, pressure perturbation and scalar field perturbation, E 1 = d ρ ̂ 0 d p ̂ 0 Π 1 − 4 3 e − 4 3 ϕ 0 ρ ̂ 0 − p ̂ 0 d ρ ̂ 0 d p ̂ 0 ϕ 1 . (45)

By tedious algebraic manipulations the nine Eqs 37–45 can be simplified. For this purpose it is convenient to define the following function of the perturbations [90], X = ω 2 p ̂ 0 + ρ ̂ 0 e − ν V − 1 r d p ̂ 0 d r e ν − λ W + 1 2 p ̂ 0 + ρ ̂ 0 e ν H 0 . (46) The resulting minimal system of differential equations is given by a set of six first order differential equations for the functions Ψ = ( K , H 1 , W , X , ϕ 1 , d ϕ 1 d r ) , which takes the form d d r Ψ + σ Ψ = 0 , (47) where σ is a matrix that depends in a complicated way on the static functions ν , λ , ϕ 0 , p ̂ 0 , ρ ̂ 0 , and also on the eigenvalue ω and the multipole number l.

Note that inside the star, the perturbation is described by the functions (K, H ₁), which parametrize the metric perturbation, (W, X) which parametrize the fluid perturbation and ( ϕ 1 , d ϕ 1 d r ) which is the scalar field perturbation. Outside the star, the system simplifies, since p ̂ = ρ ̂ = 0 . For instance, W = X = 0 when there is no fluid, and the system reduces to a system of four first order differential equations for (K, H ₁) (metric) and ( ϕ 1 , d ϕ 1 d r ) (scalar field).

Let us note that in STTs all perturbation equations are coupled with each other, whereas in the general relativistic limit, the system decouples on one hand into the metric and fluid perturbations, and on the other hand into the scalar field perturbation, which is simply governed by the minimally coupled scalar test field equation in the general relativistic limit.

At the center of the star we impose regularity of the perturbations. This means that at the center of the star the perturbation functions satisfy the constraints l ( l + 1 ) H 1 ( 0 ) − 2 l K ( 0 ) − 16 π A c 4 ( p ̂ c + ρ ̂ c ) W ( 0 ) = 0 , l Π 1 ( 0 ) + l 3 A c 4 ( 3 p ̂ c − ρ ̂ c ) ϕ 1 ( 0 ) + A c 4 ( p ̂ c + ρ ̂ c ) l 2 K ( 0 ) − ω 2 e − 2 ν c W ( 0 ) = 0 , l e − ν c p ̂ c + ρ ̂ c X ( 0 ) − l 2 K ( 0 ) = ( 8 π l 9 A c 4 ( 3 p ̂ c + 2 ρ ̂ c ) − e − 2 ν c ω 2 − l 72 a ( e − 4 3 ϕ c − 4 e − 2 3 ϕ c + 3 ) ) W ( 0 ) , d ϕ 1 d r ( 0 ) = 0 , l V ( 0 ) + W ( 0 ) = 0 , H 2 ( 0 ) − K ( 0 ) = 0 , H 0 ( 0 ) − K ( 0 ) = 0 . (48)

At infinity the massive scalar field is exponentially suppressed. Therefore the scalar perturbation is effectively asymptotically decoupled from the metric perturbations. Consequently, sufficiently far from the star, the oscillation can be described by the standard Zerilli function Z _g (r) given by H 1 ( r ) = 3 M ( n r + M ) − n r 2 r l + 1 ( 2 M − r ) ( n r + 3 M ) Z g ( r ) + r r l ( r − 2 M ) F g ( r ) , (49) K ( r ) = n r 2 ( n + 1 ) + 3 M ( n r + 2 M ) r l + 2 ( n r + 3 M ) Z g ( r ) + 1 r l F g ( r ) , (50) where n = l(l + 1)/2, and F _g (r) is some supplementary function. Then the two first order equations for (K, H ₁) can be rewritten as a more standard Schrödinger-like equation for Z _g d 2 Z g d y 2 = ( G g ( r ) − ω 2 ) Z g (51) with tortoise coordinate d y d r = e λ − ν and the effective potential G _g (r) for space-time perturbations G g ( r ) = 2 ( r − 2 M ) r 4 ( n r + 3 M ) 2 n 2 ( n + 1 ) r 3 + 3 M n 2 r 2 + 9 M 2 n r + 9 M 3 . (52)

Note that asymptotically, the potential goes to zero, G _g (r → ∞) → 0. Therefore the asymptotic behaviour of the perturbation is given by the combination of an ingoing and an outgoing solution of the form Z g ∼ A i n e − i ω y + A o u t e i ω y . (53)

We now consider the scalar field. Since this component decouples exponentially from the metric perturbation sufficiently far from the star, the scalar field equation can also be reformulated as a Schrödinger-like equation. Defining ϕ 1 ( r ) = r l Z s ( r ) , (54) the equation then becomes d 2 Z s d y 2 = ( G s ( r ) − ω 2 ) Z s (55) with the (same) tortoise coordinate y, and G _s (r) denotes the effective potential for scalar perturbations G s ( r ) = 1 − 2 M r 1 r 2 l ( l + 1 ) + 2 M r + 1 6 a . (56) Note that in this case, the potential goes asymptotically to G s ( r → ∞ ) → 1 6 a = m ϕ 2 . Therefore the scalar perturbation is asymptotically also given by a combination of outgoing and ingoing waves, but now of the form Z s ∼ A i n e − i Ω y + A o u t e i Ω y , (57) where Ω satisfies the following dispersion relation for the scalar field perturbations Ω 2 = ω 2 − m ϕ 2 . (58)

The appearance of Ω is related to the fact that the scalar field perturbations cannot propagate at the speed of light, like the space-time perturbations, since for finite values of a, the scalar field is massive. Only in the limit of a → ∞ the scalar field becomes massless and we recover Ω = ω.

We now briefly summarize the numerical implementation of the quasinormal mode calculations. The very first step is the calculation of the background configurations. To this end we solve the static equations using Colsys [91]. The solutions are obtained by employing a compactified coordinate x = r r s + r , which allows to impose the physical boundary conditions exactly at the center of the star, at its surface r _s and at infinity. The input parameters are the central pressure p ̂ c and the scalar field mass m _ϕ . The system has to be complemented with an equation of state. In this paper we shall focus on one particular choice for the equation of state, the SLy EOS [92], which is a representative model that captures the basic features of realistic neutron star models. We implement this equation of state by using a piece-wise polytropic approximation [93].

Once a particular background star with mass M is calculated, it is used to calculate the coefficients of the matrix σ of Eq. 47. Then we proceed to calculate the quasinormal modes for a particular configuration and a particular multipole number l. To do so, we solve the perturbation equations for a fixed value of ω in three different steps.

First we obtain two independent solutions of Eq. 47 inside the star, satisfying the conditions (48), and imposing X(r _s ) = 0. Second, these two solutions are continued outside the star, by requiring continuity of the metric and scalar field perturbations. These solutions are obtained up to some point r _i > r _s . We choose this r _i so that ϕ ₀(r _i ) ≲ 10⁻⁴. In the third step, we calculate the phases of the metric and scalar perturbations in the background of a Schwarzschild solution with mass M. We impose purely outgoing wave solutions at infinity for both phases. To do so, we make use of the exterior complex scaling method on Eq. 51 and Eq. 55. For more details we refer the reader to [35,36,51,55,56,68,94,95].

Then we check if the phases obtained for the asymptotic behaviour of the perturbation, which satisfy the outgoing wave behaviour, match with the full perturbative solution obtained in a region 0 < r < r _i . If they match, then ω is the eigenvalue of a quasinormal mode. If not, we repeat the process for different values of ω until matching is achieved. In this way, we investigate the spectrum for numerous configurations for several values of the scalar field mass m _ϕ . The real part of the eigenvalue, ω _R , determines the frequency of oscillation. The inverse of the imaginary part determines the characteristic time of the perturbation, τ = 1/ω _I (the damping time for stable perturbations, the instability timescale for unstable ones).

We start our discussion of the results with a detailed analysis of the properties of the quadrupole f-mode. The l = 2 fundamental mode is probably the most interesting mode as regards to astrophysical scenarios, since simulations in General Relativity show, that it tends to dominate the ringdown spectrum after a merger.

In Figure 2 we show the l = 2 fundamental mode for several values of the scalar field mass m _ϕ . In the left panel we show the frequency as a function of the total mass (in units of solar masses), and in the right panel the damping time (in seconds) as a function of the total mass. Different colors symbolize different values of the parameter a, and thus the scalar field mass m _ϕ , with purple, cyan, red, orange and green for m _ϕ = 1.08, 0.343, 0.108, 0.0343 and 0.0108 neV, respectively. In blue we show the massless case, and in black the general relativistic values for comparison. Note that for 1.08 neV the values are already very close to those of General Relativity, while below m _ϕ = 1 peV the values do not deviate significantly from the blue curve. The range of masses considered, 0.0108 ≤ m _ϕ ≤ 1.08 neV, is compatible with current observations and constraints on the mass of a hypothetical ultra light boson [74,75].

Overall we observe that the frequency and the damping time do not change drastically when going to these STTs, finding typical variations within 10% of the general relativistic value. This behaviour is reminiscent of the previous observations for the axial modes [68]. Figure 2 (left) shows that a decrease of the mass of the scalar field leads to a decrease of the frequency, except for values of the stellar mass close to one solar mass, where the frequency rises slightly from the general relativistic value. Regarding the damping time, Figure 2 (right) shows that the overall value of τ decreases the more, the lighter the scalar field is.

To better understand the dependence of the fundamental mode on the mass m _ϕ of the scalar field, we now fix the mass M of the star, while we vary m _ϕ . In Figure 3 we show the l = 2 fundamental mode as a function of the scalar field mass m _ϕ . Each curve corresponds to a family of stars with fixed value of the mass, with orange, blue, pink and green for M = 2, 1.8, 1.5, and 1.2M _⊙, respectively. On the left we show the frequency and on the right the damping time, both normalized to the respective value in General Relativity. The figure clearly shows that the larger the scalar mass, the closer the frequency comes to the general relativistic value. In fact, the frequency and the damping time deviate only significantly when the scalar field mass is below the neV. The shortest damping time occurs for massless scalar fields, and the largest deviation occurs for not very massive neutron stars with M = 1.2M _⊙. We note that for M = 1.2M _⊙, the frequency rises slightly as the scalar mass is decreased. However, in general for sufficiently massive neutron stars, the maximum deviation of the frequency appears also for massless scalar fields, with a reduction of the frequency of around 10%.

The presence of the scalar field allows for an additional type of l = 2 mode. In General Relativity this mode would correspond to the mode of an independent minimally coupled scalar field in the background of the neutron star, and therefore not be of much interest. In STTs, however, the equations for the scalar field perturbation and the matter and metric perturbations are coupled for polar modes. Therefore this new type of mode is necessarily present, and is dubbed l = 2 ϕ-mode, distinguishing it from the previously discussed l = 2 f-mode.

In Figure 4 we show the l = 2 ϕ-mode as a function of the mass M of the neutron star for several values of the scalar field mass m _ϕ . For reference, we include the mode that results from a minimally coupled scalar perturbation on the General Relativity (GR) background (black curve). Interestingly, there is very little dependence of the frequency and the damping time on the neutron star mass. Moreover, for most of the relevant mass range of the scalar field in these STTs the frequencies of these ϕ-modes are lower than the frequencies of the corresponding f-modes. In contrast, the damping times of the ϕ-modes (of the order ms) are much shorter than the damping times of the corresponding f-modes. In a ringdown spectrum the ϕ-modes will therefore fast decay, leaving the f-modes to dominate the spectrum.

In Table 1 we show some values of the frequency and damping time of both quadrupole modes: the f-mode and the ϕ-mode. In particular, we show values for neutron stars of mass M = 2 M _⊙ for several values of the scalar mass m _ϕ . For comparison we also show the respective values for the case of General Relativity.

Let us also mention here that, although we have not made a systematic study of the l = 2 excited modes, they also exist in the STTs. Our numerical results indicate that the overtones for both the f-mode and the ϕ-mode always possess larger frequencies and shorter damping times than the corresponding ground states.

Gravitational monopole radiation does not exist in GR. Hence, the ringdown phase in astrophysically realistic scenarios is expected to be dominated by the l = 2 modes. In neutron stars the dominant mode would be the f-mode, discussed above. However, the additional scalar degree of freedom present in STTs implies that, in astrophysical scenarios, apart from the l = 2 f-mode, also the radial l = 0 modes could play an important role, when the neutron stars carry scalar hair.

The analysis of the l = 0 modes in STTs reveals indeed a rich spectrum, consisting of two families of modes, named again according to their general relativistic limits. First, there are the modes that reduce to the oscillations of the nuclear matter in the general relativistic limit. These are the fundamental pressure-led mode, the F-mode, and its excitations, the H₁-mode, the H₂-mode, etc. Second, there are the modes that reduce to oscillations of an independent minimally coupled scalar field in the general relativistic limit. These are the l = 0 ϕ-modes.

As discussed above, in General Relativity the equations for these two types of modes are completely decoupled. The physically interesting modes in General Relativity are the pressure-led modes. They represent normal modes, since l = 0 modes are confined to the interior of the stars here. The fundamental pressure-led mode, the F-mode, is of particular relevance, since it shows, that neutron stars become radially unstable beyond their maximum mass.

For the F-mode the eigenvalue ω is a positive real number, as long as the mass M of the neutron star increases with increasing central pressure, and the star is radially stable. Then ω becomes zero as the maximal neutron star mass is reached. Thus here a zero mode is encountered where the star is only marginally stable. Beyond the maximum mass of the star, ω becomes purely imaginary with ω _I negative, i.e., the star becomes radially unstable.

The ϕ-modes, on the other hand, could in principle propagate outside the star, in case some external scalar field were minimally coupled to General Relativity, and a quasinormal mode analysis does yield a spectrum of such damped radial modes. However, the motivation for the presence of such a field in the environment of a neutron star would be typically lacking.

In STTs the picture changes almost completely. Only the instability of the stars beyond the maximum mass is still revealed by a purely imaginary eigenvalue with negative ω _I . Most importantly, however, scalar radiation is a natural effect in STTs, since the scalar field is a relevant gravitational degree of freedom, that is intimately coupled with the tensor degrees of freedom. Consequently, the stable normal modes from General Relativity now turn into propagating and thus quasinormal modes (see e.g., the toy model studied in [98]): the presence of the gravitational scalar field now allows for these l = 0 modes to propagate outside the star.

We exhibit in Figure 5 the frequency ω _R in kHz (left) and the imaginary part ω _I in 1/ms (right) for the F-mode versus the total mass M in solar masses M _⊙ for the neutron stars. As before, the colors represent different values of the scalar field mass m _ϕ , and the general relativistic limit is shown in black. On the unstable neutron star branches beyond the maximum mass ω _I = 1/τ represents the inverse instability timescale.