Polar Quasinormal Modes of Neutron Stars in Massive Scalar-Tensor Theories

We study polar quasinormal modes of relativistic stars in scalar-tensor theories, where we include a massive gravitational scalar field and employ the standard Brans-Dicke coupling function. For the potential of the scalar field we consider a simple mass term as well as a potential associated with $R^2$ gravity. The presence of the scalar field makes the spectrum of quasinormal modes much richer than the spectrum in General Relativity. We here investigate radial modes ($l=0$) and quadrupole modes ($l=2$). The general relativistic $l=0$ normal modes turn into quasinormal modes in scalar-tensor theories, that are able to propagate outside of the stars. In addition to the pressure-led modes new scalar-led $\phi$-modes arise. We analyze the dependence of the quasinormal mode frequencies and decay times on the scalar field mass.


INTRODUCTION
Following the first direct detections of gravitational waves, mostly emitted from merging black holes (Abbott et al., 2016b(Abbott et al., ,a, 2017b(Abbott et al., ,a, 2019a, there have been numerous further detections. Detections of merging neutron stars are much rarer events (Abbott et al., , 2020, including, in particular, the observation of the electromagnetic counterpart, which has opened the age of multi-messenger gravitational wave astronomy (Coulter et al., 2017;Abbott et al., 2017cAbbott et al., ,d, 2019b. With this new channel of observation, it is possible to directly test the strong gravity regime. Therefore, the study of the properties of the gravitational waves emitted by astrophysical sources has become ever more important. Gravitational waves from merging compact objects possess three phases, inspiral, merger and ringdown. The resonant frequencies that dominate the ringdown can be studied using quasinormal modes (QNMs) (Kokkotas and Schmidt, 1999;Nollert, 1999;Berti et al., 2009;Konoplya and Zhidenko, 2011). Whereas currently the precision of the detected signals is mostly considered too low to extract information from the ringdown, detectors are expected to achieve the necessary sensitivity in the near future (Berti et al., 2018a;Barack et al., 2019). However, already now first attempts have been made to extract not only the dominant QNM but also subdominant QNMs from the ringdown phase (Giesler et al., 2019;Bhagwat et al., 2020;Jiménez Forteza et al., 2020;Capano et al., 2021). Clearly, QNMs represent crucial observables.
In the case of the QNMs of neutron stars, one additional difficulty is that their properties including their QNM spectrum depend on the equation of state (EOS), which describes the effective relation between energy density and pressure of the matter that composes the interior of the star. So far, the EOS is not well understood, and numerous models have been proposed (see e.g., (Haensel et al., 2007)). Nonetheless, current observations allow to impose many constraints on the EOS (Lattimer and Steiner, 2014;Antoniadis et al., 2013;Özel and Freire, 2016;Most et al., 2018) and since the spectrum of neutron stars is much richer than the spectrum of black holes, even constraints on the theory of gravity are possible (Berti et al., 2018a(Berti et al., ,b, 2015. Moreover, universal relations between properly scaled global parameters of the neutron stars provide almost EOS independent relations, which can strongly enhance the analysis of neutron stars (Yagi and Yunes, 2017;Doneva and Pappas, 2018).
The spectrum of QNMs of neutron stars has already been well studied in General Relativity Kokkotas, 1996, 1998;Kokkotas et al., 2001;Benhar et al., 2004;Blázquez-Salcedo et al., 2013Mena-Fernández and González-Romero, 2019;Völkel and Kokkotas, 2019). When the background spacetime is static and spherically symmetric, the perturbations decouple into two independent channels, the axial (odd-parity) channel and the polar (even-parity) channel. Axial perturbations couple only to spacetime oscillations, and possess the so called rapidly damped w-modes (Kokkotas and Schutz, 1992). Polar perturbations couple also to the matter of the star, thus their spectrum is much richer. In the simplest case one finds pressure-driven modes (the fundamental f-mode and the excited p-modes) as well as spacetime modes (w-modes).
A neutron star may possess, however, still further modes. In General Relativity these are a set of radial normal modes. These represent undamped modes, that turn into unstable modes when for a given equation of state the maximum of the neutron star mass is passed (Chandrasekhar, 1964b,a;Bardeen et al., 1966;Meltzer and Thorne, 1966;Chanmugam, 1977;Glass and Lindblom, 1983;Vaeth and Chanmugam, 1992;Datta et al., 1998). The usual nomenclature for these modes is the fundamental F-mode and the excited H-modes. In General Relativity these modes are irrelevant for the ringdown, though, since radial perturbations cannot propagate outside the neutron stars.
Because of their extreme compactness neutron stars represent excellent laboratories to test General Relativity and alternative theories of gravity (see e.g., (Will, 2006;Faraoni and Capozziello, 2011;Berti et al., 2015;Saridakis et al., 2021)). For instance, universal relations in alternative theories of gravity can differ appreciably from those of General Relativity (Yagi and Yunes, 2017;Doneva and Pappas, 2018). Moreover, the spectrum of neutron stars is typically much richer in alternative theories of gravity, since additional degrees of freedom are present. Therefore constraints on alternative theories of gravity are possible (Berti et al., 2018a(Berti et al., ,b, 2015.
We now focus on QNMs in alternative theories of gravity with an additional scalar degree of freedom (see e.g. ). The QNMs of neutron stars in scalar-tensor theories STTs, were first considered in (Sotani and Kokkotas, 2004). Here the polar f-and p-modes were obtained by making use of the Cowling approximation, which simplifies the calculations considerably, since the perturbations of the spacetime and the scalar field are frozen (see also (Yazadjiev et al., 2017) for the inclusion of rotation). The gravitational axial w-modes were studied first in (Sotani and Kokkotas, 2005). As noted above, there is no coupling to the fluid or to the scalar field in the axial case, making these studies much simpler than the polar ones. For the axial QNMs the universal relations were investigated (Altaha Motahar et al., 2018), also for a massive scalar field with self-interaction (Altaha . When one ventures beyond the Cowling approximation (see (Blázquez-Salcedo et al., 2020;Krüger and Doneva, 2021) for its effects on the fundamental quadrupole mode), qualitatively new types of polar modes arise, since scalar radiation can be produced. The detection of scalar radiation would represent a most important discovery, while the non-detection of scalar radiation would allow to put constraints on the theory. So far, the presence of scalar modes in STTs has only been studied in an exploratory way. The sector of radial neutron star oscillations has been explored first in (Mendes and Ortiz, 2018), where scalar radiation from spontaneously scalarized neutron stars was addressed. Recently, some scalar modes and also quadrupole modes have been obtained in massive Brans-Dicke theories (Blázquez-Salcedo et al., 2020).
Such massive Brans-Dicke theories are closely related to f (R) theories, since the latter can be reformulated in terms of STTs (Sotiriou and Faraoni, 2010;De Felice and Tsujikawa, 2010;Capozziello and De Laurentis, 2011). Among these theories, in particular, R 2 gravity has received much attention in recent years. R 2 gravity is based on the Lagrangian f (R) = R+aR 2 , with coupling parameter a. When reformulated and studied in the Einstein frame, a scalar field potential results, where the mass for the scalar field depends on the coupling parameter a, m ϕ ∼ a −1/2 Staykov et al., 2014;Yazadjiev et al., 2015;Astashenok et al., 2017). For a → 0 the general relativistic limit is obtained, in contrast, for a → ∞ a particular Brans-Dicke theory arises.
Neutron star properties have been studied in R 2 gravity in (Orellana et al., 2013;Staykov et al., 2014;Yazadjiev et al., 2014Yazadjiev et al., , 2015Astashenok et al., 2017). Axial QNMs have been obtained in Blázquez-Salcedo et al., 2019). In R 2 gravity the frequencies typically deviate significantly from the values in General Relativity, but the damping times differ appreciably only for small coupling constant a. Polar modes have been studied in the Cowling approximation in R 2 theory (Staykov et al., 2015)). Only recently, the full polar QNMs have been addressed without making use of the Cowling approximation, where, in particular, ultra long lived modes were shown to exist in the radial sector (Blázquez-Salcedo et al., 2020).
The present paper is devoted to a more detailed study of the QNMs of the full polar perturbations of realistic neutron stars in massive STTs, extending the previous analyses (Staykov et al., 2015;Blázquez-Salcedo et al., 2020). The paper is organized as follows. In section 2 we revisit the static neutron stars to fix the notation. In section 3 we provide the full polar perturbations for neutron stars, and derive the equations and boundary conditions describing the oscillations. In section 4 we make use of this formalism in order to calculate the polar modes. We here focus on the l = 2 fundamental mode and the l = 0 radial modes, analyzing their dependence on the parameters of the STTs, and the total mass of the configurations. We end the paper with our conclusions and an outlook.

Theoretical framework
We consider massive STTs described by the action in the Einstein frame (G = c = 1) (Wagoner, 1970) 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 We choose two examples for the potential V (φ) (Faraoni, 2009): V I represents simply a mass term, and V II is related to R 2 gravity (Faraoni and Gunzig, 1999;Yazadjiev et al., 2014Yazadjiev et al., , 2015Bhattacharya and Majhi, 2017) To see the relation with R 2 gravity, we recall its action in the Jordan frame where R * is the Ricci scalar associated with the metric g * µν . The parameter a is a positive constant with units of [length] 2 , and controls the strength of the R 2 deviation from General Relativity. This theory can be recast into a particular Brans-Dicke STT with Jordan frame action ,  S with scalar field ψ, and potential U(ψ) The transition to the Einstein frame follows with help of the relations In the Einstein frame with the new metric g µν and the scalar field φ we then obtain the action (1) with the potential The potentials in the two frames are related by This action (1) now has an explicit kinetic term for the scalar field. The parameter a is related to the mass of the scalar, Note that, when the scalar field is weak, the potential is essentially quadratic, given by V (φ) ∼ φ 2 3a . 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., (Faraoni and Gunzig, 1999;Bhattacharya and Majhi, 2017).
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, with Einstein tensor G µν = R µν − 1 2 Rg µν . The energy-momentum tensor of the scalar field T (S) is given by T and the energy-momentum tensor of the matter T (M ) , which is assumed to be a perfect fluid, is 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 pressurep and the densityρ defined in the Jordan frame are and the equation of state is a relation of the formρ =ρ(p). The field equation for the scalar field in (1) is where T (M ) is the trace of the energy-momentum tensor of the matter.
This is a provisional file, not the final typeset article

Static neutron stars with scalar hair
For static and spherically symmetric neutron stars we consider the following Ansatz for the metric The scalar field, energy density and pressure are simply given by φ = φ 0 (r),ρ =ρ 0 (r) andp =p 0 (r), respectively, and the four-velocity of the static fluid is given by u (0) = −e ν dt.
Inside the star, the equations for the static functions are The second order differential equation (21) for φ is obtained from the scalar field equation (16) at the static level, while the others come from the Einstein equations (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 ) 4p 0 and ρ 0 = A(φ 0 ) 4ρ 0 , will appear in some formulas below (to simplify notation).

Frontiers
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 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 pressurep 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(φ 0 (r s ))r s , see  for more details. Outside the neutron star there is no fluid, thusp 0 =ρ 0 = 0 in equations (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 2 gravity), and the scalar field mass m φ assumes several values in the physically interesting mass range (Naf and Jetzer, 2010;Brito et al., 2017). 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.

Setup
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: 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 (Regge and Wheeler, 1957;Zerilli, 1970;Thorne and Campolattaro, 1967;Price and Thorne, 1969;Thorne, 1969;Campolattaro and Thorne, 1970;Thorne, 1980;Detweiler and Lindblom, 1985;Chandrasekhar and Ferrari, 1991;Chandrasekhar et al., 1991b,a;Ipser and Price, 1991;Kojima, 1992;Fernandez-Jambrina and Gonzalez-Romero, 2003). For polar perturbations, the spherical harmonics transform as . For a study of the current theory in the axial case we refer the reader to .
The Ansatz for the polar perturbations of the metric is in the order of (t, r, θ, ϕ) in the rows and columns of the matrix. The functions H 0 , H 1 , H 2 , 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 where again the function φ 1 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 while the perturbation of the 4-velocity is given by 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, Concerning the energy density and pressure, the barotropic equation of state in the Jordan frame implies a relation between the perturbations δp, δρ and δφ

Equations of the polar perturbations
Employing the Ansatz shown in the previous section for the perturbations on the field equations (12) and (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 equation (12) result in six ordinary differential equations In addition, it imposes the constraint The scalar field equation (16) results in The barotropic condition on the equation of state (36) becomes a relation between the energy density perturbation, pressure perturbation and scalar field perturbation, By tedious algebraic manipulations the nine equations (37)-(45) can be simplified. For this purpose it is convenient to define the following function of the perturbations (Lindblom and Detweiler, 1983), 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 dr ), which takes the form d dr 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 1 ), which parametrize the metric perturbation, (W, X) which parametrize the fluid perturbation and (φ 1 , dφ 1 dr ) which is the scalar field perturbation. Outside the star, the system simplifies, sincep =ρ = 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 1 ) (metric) and (φ 1 , dφ 1 dr ) (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.

Boundary conditions and asymptotic behaviour
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 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 K(r) = nr 2 (n + 1) + 3M(nr + 2M) r l+2 (nr + 3M) where n = l(l + 1)/2, and F g (r) is some supplementary function. Then the two first order equations for (K, H 1 ) can be rewritten as a more standard Schrödinger-like equation for Z g with tortoise coordinate dy dr = e λ−ν and the effective potential G g (r) for space-time perturbations G g (r) = 2(r − 2M) r 4 (nr + 3M) 2 n 2 (n + 1)r 3 + 3Mn 2 r 2 + 9M 2 nr + 9M 3 .
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 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 the equation then becomes with the (same) tortoise coordinate y, and G s (r) denotes the effective potential for scalar perturbations Note that in this case, the potential goes asymptotically to G s (r → ∞) → 1 6a = m 2 φ . Therefore the scalar perturbation is asymptotically also given by a combination of outgoing and ingoing waves, but now of the form where Ω satisfies the following dispersion relation for the scalar field perturbations 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 Ω = ω.

Outline of the numerical method
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 (Ascher et al., 1979). The solutions are obtained by employing a compactified coordinate x = r rs+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 pressurep 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 (Douchin and Haensel, 2001), 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 (Read et al., 2009).
Once a particular background star with mass M is calculated, it is used to calculate the coefficients of the matrix σ of equation (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 equation (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 φ 0 (r i ) 10 −4 . 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 equations (51)  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 φ .

RESULTS
Because of the nature of the polar perturbations, the spectrum of polar perturbations is much richer than the spectrum of axial perturbations, which involve only spacetime perturbations. Polar modes on the contrary feature different families of modes. In General Relativity static and spherically symmetric neutron stars not only possess spacetime modes (w-modes), they also possess modes related to the fluid perturbation. For l ≥ 2, the spectrum is dominated by the fundamental mode (f-mode), a nodeless fluctuation driven by pressure oscillations inside the star, which typically possesses the lowest frequency. But there are also excited pressure modes (p-modes) with higher values of the frequency.
When considering radial perturbations in General Relativity, stars are seen to possess a family of normal modes with ω 2 ∈ R, that become unstable beyond the maximum mass of the equation of state. These modes are well known in the literature (Chandrasekhar, 1964a,a;Bardeen et al., 1966;Meltzer and Thorne, 1966;Chanmugam, 1977;Glass and Lindblom, 1983;Vaeth and Chanmugam, 1992;Datta et al., 1998). However, in General Relativity radial perturbations cannot propagate gravitational radiation outside the neutron star.
In STTs there is another degree of freedom, the scalar field φ. Therefore similarly to the scalar modes of hairy black holes (Blázquez-Salcedo et al., 2016b, 2017Blázquez-Salcedo et al., 2019), there are additional modes for scalarized neutron stars. The l = 0 modes can propagate outside the neutron stars in STTs, which makes them relevant for the study of gravitational waves. In addition, of course, dipole (l = 1) modes arise. These and the l = 0 modes are supported by the scalar field, but they are coupled via the field equations with oscillations of the metric and the neutron star fluid.

Quadrupole modes
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 The frequency is normalized in terms of the corresponding frequency for a star of given mass in General Relativity. This value is obtained asymptotically as the mass of the scalar field is increased.
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 (Naf and Jetzer, 2010;Brito et al., 2017).
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 . 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 φ as well as for General Relativity. 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.
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.

Radial modes
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 pressureled mode, the F-mode, and its excitations, the H 1 -mode, the H 2 -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 (Kokkotas and Schutz, 1986)): the presence of the gravitational scalar field now allows for these l = 0 modes to propagate outside the star.
We exhibit in Fig. 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.
On the stable neutron star branches, in contrast, the very small positive values of ω I indicate the inverse damping time. While numerical inaccuracy does not allow us to precisely extract these values, the calculations indicate that the damping time is on the order of ∼ 10 5 years or larger. This means that these F-modes are ultra long lived (Blázquez-Salcedo et al., 2020). In fact, these modes exist up to m φ → ∞, where they turn into the normal modes of the respective general relativistic stars. Interestingly, long lived scalar radiation was also detected in core collapse processes in massive STTs (Sperhake et al., 2017;Rosca-Mead et al., 2020b;Rosca-Mead, 2020;Rosca-Mead et al., 2020a).
In Figure 6 we show the frequency ω R in kHz (left) and the inverse damping time ω I = 1/τ in 1/ms (right) versus the Compton wavelength L φ = 1/m φ in km for the F-mode and the excited H-modes, H 1 -H 3 , for a fixed value of the central density on the stable neutron star branch.
We note, that also for the excited modes the spectrum is qualitatively similar to the GR case. But again the fundamental difference is that the normal modes of General Relativity turn into quasinormal modes in the STTs, that are allowed to propagate outside the star, since the pressure-led modes induce also oscillations in the scalar field, which itself is coupled to the metric functions in the field equations. Moreover, the excited H-modes are ultra long lived, as well. As an aside we note, that in STTs with spontaneously scalarized neutron stars, analogous ultra long lived quasinormal modes have not been observed (Mendes and Ortiz, 2018). 10 100 For comparison Figure 6 (left) also shows the frequency which simply corresponds to the inverse of the Compton wavelength of the scalar field, and the size of the neutron star where R s represents the radius of the family of stars. The figure shows, that for small values of the Compton wavelength L φ the frequencies ω R of the F-mode and the excited H-modes assume almost constant values, that increase with increasing excitation.
The (almost) m φ -independence ends for the fundamental F-mode when L φ reaches the size of the star, and, analogously, for the excited H-modes, when appropriate fractions of the size of the star are reached. For the F-mode this happens when the Compton wavelength reaches L φ = 2R s with m φ = 0.052 neV, and for the first three H-modes for m φ = 0.072, 0.16, 0.19 neV, respectively. For larger values of the Compton wavelength L φ the frequencies decrease according to ω R = 1/2πL φ , thus the frequencies simply follow the scale given by the mass of the scalar field, as depicted by the orange curve in the left figure.
Figure 6 (right) shows that the inverse decay times 1/τ of the excited modes are indeed also very small. When estimating numerically the damping distances for these modes one finds that they should be equal or larger than ∼ 10 5 ly. This corresponds to the size of a large galaxy like the Milky Way, and possibly even larger.
As discussed above, in addition to the pressure-led l = 0 modes there are also the scalar-led l = 0 modes, the φ-modes. Figure 6 also shows the fundamental l = 0 φ-mode. Its frequency ω R follows always the scalar mass m φ , as seen by the overlap of the green dots (φ-mode) and the orange curve (left). However, its imaginary part ω I = 1/τ increases rapidly with increasing Compton wavelength L φ , showing that these modes possess much shorter lifetimes.
In Figure 7 we exhibit the frequency ω R (kHz) (left) and the damping time τ (ms) (right) versus the total mass of the neutron star M (M ⊙ ) for the l = 0 φ-mode. As already seen for the l = 2 φ-mode, the frequency ω R changes only very little with the neutron star mass. The damping time shows so little dependence only in the limiting case of General Relativity, and for the massless scalar field.
As a comment, let us note here that scalar-led l = 0 modes appear also in simulations of oscillating and collapsing neutron stars in chamaleon theories, in addition to the fluid modes (Dima et al., 2021).

Comparison of potentials
We finally briefly address the effect of the two different potentials, V I and V II (see Eqs. (4)), recalling that V I represents a STT with a mass term only, whereas V II is the potential derived from R 2 gravity. In Figure 8 we exhibit a comparison of the l = 2 f-mode for both potentials and scalar mass m φ = 0.108 neV. The insets highlight the differences.
In Figures 9 and 10 we exhibit a comparison of the radial modes. Figure 9 shows the l = 0 F-mode versus the neutron star mass for several values of the scalar field mass. Figure 10 shows the l = 0 F-mode and the l = 0 φ-mode versus the Compton wavelength of the scalar field for fixed central density.
From the figures we conclude, that for the range of masses m φ of the scalar field that are of interest (Naf and Jetzer, 2010;Brito et al., 2017), i.e., 10 −9 m φ 10 −13 eV, there is very little deviation between the quasinormal modes in these two theories. Deviations are most noticeable for the most compact configurations, close to the maximum mass, and in particular in the damping times.  Figure 10. Frequency ω R in kHz (left) and and inverse damping time ω I = 1/τ in 1/ms (right) versus the Compton wavelength L φ = 1/m φ in km for the l = 0 F-mode and l = 0 φ-mode at fixed central densitŷ ρ 0 . The colors represent the two potentials V I and V II .

CONCLUSIONS
In this paper we have studied polar quasinormal modes of neutron stars in two STTs, one possessing a simple scalar field mass term and one corresponding to R 2 gravity. We have presented the set of equations and the boundary conditions necessary to obtain these modes. The scalar field mass m φ leads to a dispersion relation between the frequency of the spacetime oscillations and the scalar field oscillations as infinity is asymptotically approached.
We have analyzed the quasinormal modes for the SLy equation of state, a realistic equation of state that yields for static neutron stars in General Relativity a maximum mass slightly above two solar masses. In the STTs studied, the value of the neutron star maximum mass increases with decreasing scalar field mass, reaching about 2.25 solar masses. For the scalar field we have covered, in particular, the physically interesting mass range (Naf and Jetzer, 2010;Brito et al., 2017).
We have investigated the l = 2 modes, starting with the fundamental f-mode, and studied the effect of the mass of the gravitational scalar field of the STTs. Moreover, we have studied the l = 2 φ-modes, which represent a set of additional quadrupole radiation modes, present only in STTs. While they depend strongly on the scalar field mass, they depend only weakly on the neutron star mass.
In General Relativity, the quadrupole quasinormal modes represent the lowest multipole modes, that can propagate gravitational radiation. In STTs there is, of course, additional gravitational radiation, corresponding to l < 2 modes. Here we have analyzed the l = 0 quasinormal modes present in these STTs. In fact, the radial perturbations lead to a rich spectrum of modes.
We have shown that the pressure normal modes of General Relativity, the fundamental F-mode and its excitations, turn into propagating quasinormal modes in these STTs. For the observationally relevant range of the scalar mass these modes are ultra long lived. Thus the damping distance of the corresponding gravitational radiation corresponds to large distances of 10 5 or more light years.
When considering the dependence of these F-modes and the excited H-modes on the Compton wavelength of the scalar field, one finds almost constant values for their frequencies as long as the Compton wavelength is smaller than the size of the star (F-mode) or appropriate fractions of it (H-modes). Beyond these points, the frequencies decay with the inverse of the Compton wavelength.
Besides these pressure-led modes, there are also the l = 0 φ-modes in these STTs. The frequencies of the φ-modes show again very little dependence on the neutron star mass as seen already for the l = 2 φ-modes. In contrast to the pressure-led modes their damping times are on the order of ms.
The next steps will include the study of the l = 1 dipole modes, and the inclusion of rotation, at least at a perturbative level. Moreover, a whole set of realistic equations of state will be employed with the aim to derive universal relations for the polar modes, analogously to the axial case . Also other closely-related theories as, for example, STTs with spontaneous scalarization will be studied.
Numerical simulations show that the ringdown after a merger is dominated by three modes (Bauswein and Stergioulas, 2015;Bernuzzi et al., 2015) that are interpreted as containing the fundamental f-mode and a mixture of the f-mode and a quasiradial F-mode of the remnant star (Takami et al., 2015). The mixing occurs because of the decrease in symmetry when going from spherical to axial symmetry.
Whereas merger simulations have been mostly performed in General Relativity so far, recently the merger of neutron stars has also been considered in STTs (Sagunski et al., 2018). Making use of an effective model to estimate the merger of neutron stars in R 2 gravity, two modes have been found to dominate the ringdown (Sagunski et al., 2018). On the other hand, the analysis of the merger of neutron stars in chameleon theory has shown that the l = 1 dipole modes tend to be suppressed by the screening mechanism . It will be interesting to perform full merger simulations in STTs to extract the effects of the scalar field on the ringdown spectrum and its dependence on the scalar field mass.