A new class of quasinormal modes of neutron stars in scalar-tensor gravity

A new class of quasinormal modes of neutron stars in scalar-tensor gravity

Abstract

Detection of the characteristic spectrum of pulsating neutron stars can be a powerful tool not only to probe the nuclear equation of state, but also to test modifications to general relativity. However, the shift in the oscillation spectrum induced by modified theories of gravity is often small and degenerate with our ignorance of the equation of state. In this letter, we find that the coupling to additional degrees of freedom present in modified theories of gravity can give rise to new families of modes, with no counterpart in general relativity, which could be sufficiently well resolved in frequency space as to allow for a clear detection. We present a realization of this idea by performing a thorough study of radial oscillations of neutron stars in massless scalar-tensor theories of gravity. We anticipate astrophysical scenarios where the presence of this class of quasinormal modes could be probed with electromagnetic and gravitational wave measurements.

pacs:
04.50.Kd, 04.80.Cc

Introduction. The study of stellar pulsations has long provided a valuable tool to scrutinize the interior of stars Unno et al. (1989). With the detection of gravitational waves (GWs) by the Advanced LIGO and Virgo Collaborations et al. (LIGO Scientific Collaboration and Collaboration) (2016); LIGO Scientific Collaboration and VIRGO Collaboration (2016, 2017a); LIGO Scientific Collaboration and VIRGO Collaboration (2017b); LIGO Scientific Collaboration and VIRGO Collaboration (2017c), the study of relativistic pulsations of compact objects has finally met observations, giving rise to a wealth of possibilities to probe the workings of gravity in the strong field regime, from the geometry of black holes to the internal structure of relativistic stars Nollert (1999); Kokkotas and Schmidt (1999); Berti et al. (2018a, b). In order to build intuition into how modifications to general relativity (GR) can alter the emission of GWs, it is insightful to consider concrete frameworks in which Einstein’s theory can be generalized. A particular framework, subject of the present study, is that of scalar-tensor theories of gravity (STTs) Damour and Esposito-Farèse (1992); Fujii and Maeda (2003). These theories are extensions to GR which include, as an additional mediator of the gravitational interaction, a (possibly self-interacting) scalar field that couples non-minimally to the metric sector. Interestingly, for some choices of the non-minimal coupling function, STTs predict a radically different phenomenology for relativistic stars, due to the existence of a nonperturbative strong field effect known as spontaneous scalarization Damour and Esposito-Farèse (1993, 1996); Salgado et al. (1998); Harada (1998).

Relativistic stars that undergo the process of spontaneous scalarization have their internal structure altered with respect to GR. Consequently, their vibration modes are also modified. Stellar pulsations in STTs have been studied in a variety of contexts, often using strong approximation schemes Sotani and Kokkotas (2004, 2005); Silva et al. (2014); Sotani (2014); Yazadjiev et al. (2017). Here we examine the problem of radial oscillations of neutron stars (NSs) in STTs in its full generality.

In GR, the main interest in radial oscillations of relativistic stars lies in the information they provide about (in)stability against gravitational collapse Chandrasekhar (1964); Kokkotas and Ruoff (2001). However, in STTs the scalar sector of the gravitational field is dynamical even in spherical symmetry, and can carry energy—and information—away from the source, which adds to the relevance of studying spherical perturbations.

Our analysis of the full radial problem reveals the existence of entirely new families of quasinormal modes (QNM), with no counterpart in GR, which could not be captured by the Cowling approximation used previously Sotani (2014). Since these modes can be explicitly related to the scalar degree of freedom, we dub them scalar- or -modes. Indeed, one could expect the coupling to new degrees of freedom to introduce new families of modes, in addition to just shifting the original spectrum. This is the case of stellar oscillations in GR, where the coupling to metric degrees of freedom not only modifies the Newtonian modes but also introduces new spacetime- or -modes Kokkotas and Schutz (1986, 1992). Interestingly, the new family of scalar modes reported in this letter may be relatively long-lived, and lie in a lower frequency range (of kHz) more accessible to current GW detectors. Also, the -mode frequencies can be clearly distinguished from those of GR even if the coupling to the scalar field is small—i.e. for stars with small scalar charges—, in which case one might expect the GR fluid modes to only acquire a fine structure (easily mimicked by a change in the nuclear equation of state). Our findings suggest that the existence of new families of stellar QNMs could constitute a clear fingerprint of theories of modified gravity containing additional (scalar, vector…) degrees of freedom, and that looking for their signatures in astrophysical data could be a promising way to improve constraints on such theories.

Framework. We focus on a class of massless STTs described by the action (we use unless specified)

 S=116π∫d4x√−g(R−2∇μϕ∇μϕ)+Sm[Ξm;a(ϕ)2gμν], (1)

where represents matter fields with action , that couple to the conformally rescaled (Jordan-frame) metric . In the following, we consider two representative coupling models,

 Model 1 (M1): a(ϕ)=[cosh(√3βϕ)]1/(3β), (2) Model 2 (M2): a(ϕ)=eβϕ2/2, (3)

where . M1 is an analytical approximation to the coupling function of a “standard” massless scalar field non-minimally coupled to gravity (see Sec. IIB in Ref. Mendes and Ortiz (2016)), whereas M2 is the simplest model leading to spontaneous scalarization Damour and Esposito-Farèse (1993). Both models agree up to the cubic term in an expansion around the “cosmological” value 1, and in this particular case have the same post-Newtonian expansion as GR Damour and Esposito-Farèse (1996).

Variation of Eq. (1) yields the field equations

 Gμν−2∇μϕ∇νϕ+gμν∇σϕ∇σϕ=8πa(ϕ)2~Tμν, (4) ∇μ∇μϕ=−4πa(ϕ)4α(ϕ)~T, (5)

where , , and is the (Jordan-frame) stress-energy tensor of matter fields, which is covariantly conserved (). We model NSs by a perfect fluid stress-energy tensor , with a barotropic equation of state (EoS) relating the pressure and the rest-mass density . The total energy density is obtained from the first law of thermodynamics, . Conservation of baryon number further implies , where denotes the four-velocity of fluid elements. For our purposes, it is enough to adopt a polytropic EoS, , consisting of two phases, a stiff core () and a soft crust (). We choose , and the transition density to be , where g/cm 2.

Background. The unperturbed equilibrium star is assumed to be static and spherically symmetric, with line element of the form

 ds2=−e2νdt2+e2λdr2+r2(dϑ2+sin2ϑdφ2), (6)

where , . The explicit form of the field equations (4) and (5) in the case of hydrostatic equilibrium can be found, e.g. in Ref. Damour and Esposito-Farèse (1993). Given a central rest-mass density value, , the background equations can be numerically integrated subject to regularity conditions at the origin and asymptotic conditions and . The stellar models thus constructed are characterized by their ADM mass and scalar charge , such that asymptotically. The stellar radius is determined by .

Families of equilibrium solutions for Models 1 and 2 are presented in Fig. 1 for and 3. For low enough central densities, the only equilibrium solutions have , in which case the fluid and metric functions recover their GR values. These solutions are depicted in solid gray in Fig. 1. Nonetheless, above a certain critical central density, the scalar field can be suddenly activated: the trivial solutions become unstable under scalar perturbations Harada (1997); Mendes (2015), and branches of “scalarized” configurations develop—i.e., solutions with a non-trivial scalar field profile. For , scalarized solutions in both M1 and M2 have the same qualitative behaviour, but they are radically different when Mendes and Ortiz (2016).

Perturbations. We describe radial perturbations of the spacetime in a gauge where the metric retains the form (6), by setting and , where the subscript refers to background quantities. The perturbed scalar field can be written as . The Jordan-frame metric becomes , with

 ~hμν=a2(0)hμν+2g(0)μνa2(0)α(0)δϕ, (7)

where and .

The radially perturbed configuration is characterized by a Lagrangian displacement vector field , in terms of which we can write the Eulerian change in the fluid four-velocity, . Under the assumption that the perturbed fluid has a negligible temperature and the same EoS as the unperturbed configuration, the pressure and energy density perturbations, and , can be related to the rest-mass density perturbation through and , where is the adiabatic exponent.

Standard manipulations of the perturbed field equations [see Supplemental Material (SM)] lead to master equations for the Lagrangian displacement and the scalar field perturbation , in terms of which the remaining perturbation variables, , , and , can be determined. They read:

 e3λ(0)(~ϵ(0)+~p(0))¨ξ−⎡⎣Γ1~p(0)a4(0)r2eλ(0)+3ν(0)(e−ν(0)a4(0)r2ξ)′⎤⎦′ +Aξξ+Aδϕδϕ+Aδϕ′δϕ′=0, (8)
 e2λ(0)−2ν(0)δ¨ϕ−δϕ′′+Bδϕ′δϕ′+Bδϕδϕ+Bξ′ξ′+Bξξ=0, (9)

where primes and dots indicate radial and temporal differentiation, respectively, and the - and -coefficients depend only on background quantities. They are displayed in the SM, together with explicit expressions for .

We solve the coupled second-order system (A new class of quasinormal modes of neutron stars in scalar-tensor gravity)-(9) in the frequency domain with the ansätze and , . We require and to be everywhere regular, and to be purely outgoing at spatial infinity. In general, solutions satisfying these conditions exist only for a discrete set of complex frequencies , which we aim to determine. As a consistency check, we also solve the system (A new class of quasinormal modes of neutron stars in scalar-tensor gravity)-(9) in the time domain, and perform a Fourier analysis to unveil the vibration modes of the solution. Details on boundary conditions and numerical techniques are deferred to the SM.

Radial quasinormal modes of NSs in STTs. Since gravity in GR is not dissipative in spherical symmetry, radial fluid oscillations are confined to the stellar interior, and are described by a set of normal modes. In STTs, stellar solutions with present the same spectrum of fluid oscillations as in GR: in this case, Eqs. (A new class of quasinormal modes of neutron stars in scalar-tensor gravity) and (9) decouple, and Eq. (A new class of quasinormal modes of neutron stars in scalar-tensor gravity) simply describes perturbations in GR. Moreover, these GR-like stellar models display additional families of modes associated with the scalar sector. For scalarized solutions, fluid and scalar field perturbations are coupled, and the spectrum changes accordingly. Fluid modes have their frequency shifted and gain an imaginary part, while scalar modes dictate the stability of the star. We elaborate on these ideas below for the cases of and .

Figure 2a displays the frequency of the first lowest-frequency radial modes for scalarized solutions in Models 1 and 2 with , as a function of stellar compactness . The fundamental () and first overtone () frequencies for general relativistic stars are displayed in solid gray, and have values kHz and kHz, respectively. These modes are shifted in the presence of spontaneous scalarization, in a way analogous to that described by Sotani Sotani (2014) within the Cowling approximation—see SM for a comparison. Note that in M1, the and mode frequencies are less shifted than in M2, which might be traced to the fact that the scalar charge of scalarized solutions in M1 is typically smaller than in M2 4. The inverse of the damping time is shown in Fig. 2b as a function of the compactness . For all modes, it is of the order of the dynamical timescale (ms) of the star.

Besides the QNMs that correspond to deformations of GR fluid modes, Fig. 2 reveals the presence of a new family of pulsation modes, which have no analogue in GR. We call them scalar- or -modes. For , they have smaller frequencies (more than 1.6 kHz lower than the mode in GR), and a slightly lower damping time than the shifted GR modes. Intuitively, the lower frequency of -modes may come from the fact that scalar field oscillations are not confined to the star, thus their wavelength is typically larger than that of fluid modes.

Figure 2c presents a more complete picture of the fundamental -mode. For stellar models with , the fundamental -mode is shown in solid gray; while it is shown in colors for scalarized solutions. Before the critical compactness for spontaneous scalarization (), scalar perturbations are decoupled from fluid oscillations, and the fundamental -mode frequency goes to zero as is approached. For solutions with , the fundamental -mode becomes unstable—this is the same instability originaly discussed in Ref. Harada (1997). As for scalarized solutions, the frequency and damping time of the fundamental -mode were already visible in Figs. 2a and 2b. Remarkably, it is the -mode that determines (in)stability to gravitational collapse for scalarized solutions: for M2, we see that its frequency goes to zero at the compactness corresponding to the turning point in the mass-density diagram 1a, becoming unstable () thereafter. For M1, all scalarized solutions are stable. For , all equilibrium solutions again have a constant scalar field profile, and the scalar modes decouple from fluid oscillations.

In Figure 3, we present the analogue of Fig. 2 for M1 with (our analysis shows that all scalarized solutions for M2 with are unstable; see SM). Remarkably, despite the fact that scalarized solutions in M1 with typically have very small scalar charges 5, extremely abrupt modifications to the spectrum can be seen in Fig. 3a. Figure 3b reveals that the typical damping is now much longer than the stellar dynamical timescale, reaching for the -mode (in agreement with simulations performed in Ref. Mendes and Ortiz (2016)). In Fig. 3c it is clear how, at the critical compactness for spontaneous scalarization, the scalar mode becomes unstable along the sequence of trivial solutions, whereas the coupling to fluid oscillations stabilizes this mode for scalarized configurations lying before the turning point in the mass diagram 1b.

Discussion. Our findings suggest that even when equilibrium stellar models in modified theories of gravity are close to GR—especially taking into account the uncertainties in the nuclear EoS—their oscillation spectrum may be radically different. In particular, new families of modes are expected to arise due to coupling to the additional degrees of freedom.

In the case of spherical oscillations in STTs, we find that the new -modes have lower frequencies and relatively long damping times compared to GR. They could be excited in a variety of astrophysical processes. For instance, Soft Gamma Repeaters—presumably magnetars Duncan and Thompson (1992); Thompson and Duncan (1995)—manifest as giant flares followed by X-ray tails with decaying times of the order of hundreds of seconds. Such tails exhibit quasi-periodic oscillations with frequencies in the range Hz Watts and Strohmayer (2007), which could in principle excite -modes, with possible imprints in the electromagnetic emission.

-modes could also be excited in binary neutron star systems. In the inspiral phase, they could become resonant with the orbital motion, draining energy from it, and thus altering the GW phase evolution Gold et al. (2012); Hinderer et al. (2016). Due to the lower -mode frequency, this effect could occur in STTs for larger binary separations than in GR, thus in a frequency band more suitable to current GW detectors. In the post-merger phase, quasi-radial scalar oscillations (i.e., the modified radial -modes in rotating systems) could be excited if a hypermassive neutron star forms Baumgarte et al. (2000). These could have non-spherical components that would emit GWs directly detectable by Advanced LIGO/Virgo.

Other astrophysical scenarios in which -modes could be excited include gravitational collapse, and phase transitions in the core of NSs.

Important ramifications of our work include studies on nonradial oscillations of spherical stars and quasi-radial modes of rotating systems, extensions to other modified theories of gravity, and detectability analyses in various astrophysical scenarios.

Acknowledgements.
We are glad to thank William East, Eric Poisson, Huan Yang, Luis Lehner, Ryan Westernacher-Schneider, Job Feldbrugge, Stephen Green, and Jonah Miller for insightful discussions and suggestions. This work was started while RM was benefiting from a CITA National Fellowship at the University of Guelph. It was supported in part by the Natural Sciences and Engineering Research Council of Canada. Research at the Perimeter Institute is supported by the Government of Canada through Industry Canada and by the Province of Ontario through the Ministry of Research and Innovation.

Appendix A Supplemental Material

a.1 Other perturbed quantities and coefficients of master equations

We describe how the perturbed metric functions and , as well as the perturbation to the rest-mass density , relate to the Lagrangian displacement and scalar field perturbations, which obey Eqs. (A new class of quasinormal modes of neutron stars in scalar-tensor gravity) and (9). (See e.g. Ref. Friedman and Stergioulas (2013) for background on perturbation theory of relativistic fluids.)

The perturbed rest-mass density obeys

 δ~ρ= −~ρ(0)δλ−3~ρ(0)α(0)δϕ−~ρ(0)ξ′−ξ~ρ′(0) −ξ~ρ(0)(2r+3α(0)Ψ(0)+λ′(0)), (10)

where . This equation follows from the perturbed equation for conservation of baryon mass, . The -component of the perturbed field equations can be solved for , resulting in

 δλ=rΨ(0)δϕ−4πre2λ(0)a4(0)(~ϵ(0)+~p(0))ξ, (11)

while from the -component, together with Eq. (11), one gets

 δν′= −4πe4λ(0)a4(0)[1+8πr2a4(0)~p(0)](~ϵ(0)+~p(0))ξ +e2λ(0)[Ψ(0)+8πr~p(0)a4(0)(2α(0)+rΨ(0))]δϕ +rΨ(0)δϕ′+4πre2λ(0)a4(0)δ~p. (12)

Equation (A new class of quasinormal modes of neutron stars in scalar-tensor gravity) is obtained from the -component of the perturbed equations of motion, , together with Eqs. (11) and (A.1), whereas Eq. (9) stems from the perturbed scalar field equation. The - and -coefficients in Eqs. (A new class of quasinormal modes of neutron stars in scalar-tensor gravity) and (9) are given in terms of the background quantities as follows:

 Aξ = −14r2(~p(0)+~ϵ(0))eλ(0)+2ν(0){e4λ(0)+6e2λ(0)−7+8rα(0)Ψ(0)(e2λ(0)+1+r2Ψ2(0))+12r2Ψ2(0)α2(0) + 16πr2e2λ(0)a4(0)[4πe2λ(0)~p2(0)r2a4(0)+~p(0)(e2λ(0)+1+3α2(0)+2rΨ(0)α(0))−~ϵ(0)(α(0)+rΨ(0))2] + Ψ2(0)r2(2e2λ(0)+6+r2Ψ2(0)−4α′(0))}+Γ12reλ(0)+2ν(0){Ψ(0)~p(0)[rΨ(0)(3e2λ(0)+7−r2Ψ2(0)+2α′(0)) − α(0)(e2λ(0)+7−9r2Ψ2(0))−8rΨ(0)α2(0)]−Ψ(0)~ϵ(0)(α(0)−rΨ(0))(e2λ(0)−1+2rα(0)Ψ(0)+r2Ψ2(0)) − 8πre2λ(0)~p(0)a4(0)[~p(0)α(0)(3α(0)−7rΨ(0))−~ϵ(0)(α(0)−2rΨ(0))(α(0)+rΨ(0))]} + ~p(0)Ψ(0)Γ′1eλ(0)+2ν(0)(α(0)−rΨ(0)), Aδϕ = 12r(~p(0)+~ϵ(0))eλ(0)+2ν(0)[8πr2~p(0)e2λ(0)a4(0)(α(0)+rΨ(0))−α(0)(3e2λ(0)−3+5r2Ψ2(0))−6rΨ(0)α2(0) + rΨ(0)(e2λ(0)+1−r2Ψ2(0)+2α′(0))]+Γ12reλ(0)+2ν(0)[−8πr2~p(0)~ϵ(0)e2λ(0)a4(0)(α(0)+rΨ(0)) + ~ϵ(0)(e2λ(0)−1)(3α(0)+rΨ(0))+r~ϵ(0)Ψ(0)(5rα(0)Ψ(0)+6α2(0)+r2Ψ2(0))+2r~p(0)Ψ(0)(e2λ(0)−3α′(0))] − ~p(0)Γ′1eλ(0)+2ν(0)(3α(0)+rΨ(0)), Aδϕ′ = eλ(0)+2ν(0)[(~p(0)+~ϵ(0))(α(0)+rΨ(0))−Γ1~p(0)(3α(0)+rΨ(0))],
 Bδϕ′ = −1r(e2λ(0)+1)+4πre2λ(0)a4(0)(~ϵ(0)−~p(0)), Bδϕ = 4πa4(0)~p(0)e2λ(0){(3α(0)+rΨ(0))[α(0)(3Γ1−5)+rΨ(0)(Γ1−3)]−3α′(0)} + 4πa4(0)~ϵ(0)e2λ(0)[(α(0)+rΨ(0))2+α′(0)]−2e2λ(0)Ψ2(0), Bξ′ = 4πa4(0)e2λ(0){~p(0)[(3Γ1−1)α(0)+(Γ1−1)rΨ(0)]−~ϵ(0)(α(0)+rΨ(0))}, Bξ = −2πe2λ(0)a4(0)Γ1r~p(0){−(~p(0)+~ϵ(0))2(α(0)+rΨ(0))(8πr2a4(0)e2λ(0)~p(0)+e2λ(0)−1+2rα(0)Ψ(0)+r2Ψ2(0)) + × [α(0)(8πr2a4(0)e2λ(0)(2~p(0)−~ϵ(0))−e2λ(0)−1)−6rα(0)Ψ(0)(α(0)+rΨ(0)) + rΨ(0)(8πr2a4(0)e2λ(0)(~p(0)−~ϵ(0))+2e2λ(0)−2−r2Ψ2(0))]}.

The Lagrangian displacement is only defined in the stellar interior (). Outside the star (), the only nonvanishing coefficients in Eq. (9) are

 Bδϕ′=−1+e2λ(0)r,Bδϕ=− 2e2λ(0)Ψ2(0).

a.2 Unstable modes of scalarized solutions with β>0

Here we analyze the stability of scalarized solutions in STTs with . For theories with , and in particular for the model (our Model 2) popularized by the works of Damour and Esposito-Farèse, stability of scalarized solutions was investigated by Harada Harada (1998) through a turning point technique. In Model 2 with , however, turning point techniques do not seem to be directly applicable. As a by-product of our work, here we examine instability to collapse for scalarized solutions in M1 and M2 with (but our results seem generic for any ).

To interpret our results it is useful to recall our thermodynamic expectations. It is well known that in GR a sequence of stellar models is secularly unstable on one side of an extremum in a mass diagram Sorkin (1982). For spherical models, a turning point also signals the onset of dynamical instability, when a radial mode changes its frequency from real to imaginary. It is easy to see how dynamical and secular stabilities are related in this case. The point of marginal stability is characterized by the existence of a zero-frequency mode, which is just a time-independent perturbation linking two nearby equilibrium solutions. From the first law of thermodynamics, a perturbation that keeps the star in equilibrium satisfies where is the injection energy and is the baryon mass. Thus, for a zero frequency perturbation involving no change in baryon mass we have , which is the necessary condition for an extremum in a mass diagram. The same argument holds in STTs, as discussed by Harada Harada (1998). The only difference is that now the control space of spherically symmetric equilibrium models is not only the baryon number, but also the asymptotic value of the scalar field, , so that the first law reads , where is the scalar charge. Thus for a zero frequency perturbation that involves no change in or .

Our results conform to these thermodynamic considerations. Figure 4a shows the inverse of the instability timescale (such that ), as a function of the normalized central rest-mass density for M1 with . We see that the branch of scalarized solutions with displays unstable modes for central densities , which is consistent with the turning point in Fig. 1b. The same applies for the family of scalarized solutions with . As can be seen back in Fig. 1b, these solutions exist for stellar models with central densities , and closer inspection reveals that the turning point in the mass diagram occurs at . In agreement with thermodynamic expectations, in Fig. 4a we see that, in the case, unstable modes only exist for . Therefore, there is only a small range of densities for which scalarized solutions with have no unstable radial modes. At around a third branch of scalarized solutions develop. From Fig. 4a we see that they all are unstable. This is also the case for families of solutions with higher values of . In all cases, the instability timescale roughly decreases as as the central density increases.

With respect to Model 2, in Ref. Mendes and Ortiz (2016) we found no evidence of the existence of stable scalarized solutions with , so we conjectured that all such scalarized solutions were unstable. This conjecture is confirmed by the perturbative analysis in this work. In Fig. 4b we show for scalarized stars in M2 with , for the three lowest- branches. We see that as soon as scalarized solutions appear, they quickly develop unstable radial modes. As the density approaches from the right, the instability timescale goes to zero—the pattern repeats for higher families. This is the critical density at which the trace of the energy-momentum tensor vanishes in the stelar center; below this density there are no scalarized solutions Mendes and Ortiz (2016).

In the main text, we interpreted the lowest frequency oscillation mode of scalarized stars as a deformation of a pure scalar mode of configurations. Correspondingly, overtones of scalarized stars were interpreted as deformations of pure fluid modes of those trivial solutions. Here we give further support for this interpretation. As an example, we focus on scalarized solutions in Model 1 with . The lower the scalar charge—i.e. the closer the compactness to the merging points with the branch of GR solutions, which occur at and —, the lowest-frequency mode is expected to be predominantly a scalar perturbation, while the overtones are expected to be mainly fluid perturbations. This expectation is verified in Fig. 5, where we compare the radial eigenfunctions of the Lagrangian displacement and the scalar field perturbation at the surface of the star (), for a sequence of scalarized solutions. For the lowest frequency mode, we see that the ratio goes to zero as the transition to GR solutions is approached, while it diverges for the overtone modes.

a.4 Comparison to the Cowling approximation

In GR, the Cowling approximation consists in keeping the metric fixed and only allowing for fluid perturbations (see for instance Ref. Friedman and Stergioulas (2013)). It becomes increasingly accurate for fluid perturbations that are confined to a small region inside the star or have a short wavelength Lindblom and Splinter (1990); Yoshida and Kojima (1997); Yoshida and Eriguchi (2001). This is the case, for instance, of torsional oscillations confined to the stellar crust and high multipole QNMs. The rationale behind it can be easily understood in the Newtonian setting: If is the characteristic wavelength of the mode, then the Laplacian of the gravitational potential is roughly , and Poisson’s equation implies that , which goes to zero as .

In STTs, the presence of a scalar field as an additional mediator of the gravitational interaction introduces an extra dilemma of whether or not to keep the scalar field fixed in the Cowling approximation scheme. Indeed, previous works made different choices, either keeping both the metric and the scalar field fixed (as in Ref. Sotani and Kokkotas (2004); Yazadjiev et al. (2017)), or allowing for scalar field perturbations while keeping the Jordan-frame metric fixed (as in Refs. Silva et al. (2014); Sotani (2014)). Here we compare our results with those obtained with the latter procedure, which was employed in Ref. Sotani (2014) in the study of radial oscillations of spherically symmetric scalarized stars. In this case, one sets , and Eq. (7) implies that . Manipulations of the perturbed field equations allow all perturbation variables to be written in terms of the Lagrangian displacement, which follows the master equation (3.15) of Ref. Sotani (2014). (Strictly speaking, this equation is only valid for M2, but it suffices to substitute the combination by and the remaining ’s by to obtain the general expression.) This version of the Cowling approximation allows scalar waves to be sourced by fluid perturbations (see Eq. (3.12) of Ref. Sotani (2014)), although they do not cause damping to the fluid motion, which is described by a set of normal modes.

Figure 6 shows a comparison between the mode frequencies obtained with our full analysis and with the Cowling approximation for M2 with . We see that the Cowling approximation overestimates the mode frequencies, especially for the fundamental fluid mode, and does not properly capture information about instability to collapse. More importantly, it does not reveal the presence of the additional families of modes reported in this work. Therefore, although the Cowling approximation can yield relatively accurate results, especially for higher overtones and higher angular momentum modes, it may fail to capture key aspects of stellar QNMs in modified theories of gravity.

a.5 Frequency domain approach

Here we describe the frequency domain techniques that we implement in order to solve the perturbation system (A new class of quasinormal modes of neutron stars in scalar-tensor gravity)-(9), with the oscillatory ansätze

 ξ(t,r)=ξ(r)eiωt,δϕ(t,r)=δϕ(r)eiωt,ω∈C. (13)

Boundary conditions

Appropriate boundary conditions for the solutions are the following.

• The radial eigenfunctions and are required to be analytic around . Equations (A new class of quasinormal modes of neutron stars in scalar-tensor gravity) and (9) imply that, in a neighbourhood of ,

 ξ(r)=N∑i=0a2i+1r2i+1,δϕ(r)=N∑i=0b2ir2i,

for some , where are free constants. The remaining coefficients are determined recursively in terms of and the background quantities.

• At the perturbed stellar surface, the pressure vanishes, implying that its Lagrangian variation is zero, . But since

 Δ~p= −Γ1~p(0){rΨ(0)δϕ+ξ′ + ξ[3α(0)Ψ(0)−ν′(0)+rΨ2(0)+2/r]},

and , this condition is automatically satisfied as long as the term within brackets is regular at the stellar surface. Therefore, it is enough to impose regularity of and at . Dividing Eq. (A new class of quasinormal modes of neutron stars in scalar-tensor gravity) by , and taking into account our choice of EoS in the stellar crust, we see that the terms which are potentially diverging in the limit are those proportional to . Therefore, we require the coefficient of the term to vanish at , which amounts to requiring that the quantity

 f(r):=γ2ξ′[e2λ(0)−1+rΨ(0)(2α(0)+rΨ(0))] +2rδϕ′(α(0)+rΨ(0))+ξ2r{−4e2λ(0)−2ν(0)r2ω2−5γ2 +6(γ2−1)(e2λ(0)+2r2α2(0)Ψ2(0))−e4λ(0)(γ2+1)+7 +4rα(0)Ψ(0)[(1+e2λ(0))(γ2−2)+2r2(γ2−1)Ψ2(0)] +r2Ψ2(0)[4α′(0)−6−2e2λ(0)+4γ2+r2(γ2−1)Ψ2(0)]} +δϕ{(γ2−1)α(0)(3e2λ(0)−3+5r2Ψ2(0)) +rΨ(0)[eλ(0)(γ2+1)+2α′(0)+(γ2−1)(r2Ψ2(0)−1)] +6r(γ2−1)α2(0)Ψ(0)} (14)

vanishes at the stellar surface, i.e. .

• We impose a purely outgoing boundary condition for the scalar field perturbation at spatial infinity,

 limr→∞δϕ(t,r)∝eiω(t−r). (15)

To write this more precisely, we first notice that setting , Eq. (9) for can be cast in the alternative form

 d2Zdr2∗+[ω2+e2ν(0)(2Ψ2(0)+e−2λ(0)−1r2)]Z=0, (16)

where the tortoise coordinate is defined through . Equation (16) admits two linearly independent solutions, , with asymptotic expansions

 Z±(r→∞)=e±iωr∗N∑i=0c±iri, (17)

where and are specifiable constants, in terms of which the coefficients , , can be determined by solving Eq. (16) order by order. The purely outgoing boundary condition singles out as the physical solution.

In summary, our task is to solve the coupled system (A new class of quasinormal modes of neutron stars in scalar-tensor gravity)-(9) subject to a set of four boundary conditions: regularity at [specifically, and ], regularity at [specifically, with defined in Eq. (A.5.1)], and an outgoing boundary condition for the scalar perturbation at spatial infinity. However, since the perturbation equations (A new class of quasinormal modes of neutron stars in scalar-tensor gravity)-(9) are homogeneous, there is an additional overall normalization freedom to the perturbed quantities. This means that the boundary value problem posed above is overdetermined. Therefore, solutions exist only for a discrete set of values of . Below we detail the numerical techniques that we employ to search for these frequencies.

Solution to the inner problem

In the domain , the perturbation equations (A new class of quasinormal modes of neutron stars in scalar-tensor gravity) and (9) with the ansatz (13) are to be solved subject to the boundary conditions , , and . To fix the normalization freedom we require .

Our strategy relies on the fact that for a homogeneous system of the form

 dx(t)dt=M(t)x(t), (18)

where is a vector function and is a matrix function, every solution can be expressed as

 x(t)=c1x1(t)+⋯+cnxn(t),

where are constants and is a set of linearly independent solutions of the system (see, e.g. Ref. Lindblom and Detweiler (1983)). It is straightforward to rewrite the perturbation equations (A new class of quasinormal modes of neutron stars in scalar-tensor gravity) and (9) for the state vector in the form (18). At , we select two linearly independent vectors satisfying the boundary conditions and [say, and ] and numerically integrate the differential equations from to (where is the transition point of the two polytropic phases, ) with such initial conditions. The integration yields two solutions, and defined in the domain and satisfying the correct boundary conditions at . The unique solution that satisfies all the boundary conditions must be a linear combination of the form for , and some constants , . We repeat this procedure for . Namely, we select three independent vectors that satisfy the boundary condition at [say, , , and , with determined in each case by solving ] and numerically integrate the system from to . This yields three solutions, , , and , defined in the domain , which satisfy the correct boundary condition at . Again, the unique solution that satisfies all the boundary conditions must be a linear combination of the form for and some constants , . Finally, we require that satisfies the proper (dis)continuity condition at , namely,

 0 =(0,Δ,0,0