Oscillation dynamics of scalarized neutron stars
OOscillation dynamics of scalarized neutron stars
Christian J. Krüger ∗ and Daniela D. Doneva
1, 2, † Theoretical Astrophysics, IAAT, University of Tübingen, 72076 Tübingen, Germany INRNE—Bulgarian Academy of Sciences, Sofia 1784, Bulgaria (Dated: February 24, 2021)
Abstract
Scalar-tensor theories are well studied extensions of general relativity that offer deviations whichare yet within observational boundaries. We present the time evolution equations governing theperturbations of a nonrotating scalarized neutron star, including a dynamic spacetime as well asscalar field within the framework of such scalar-tensor theories. We employ a theory that allowsfor a massive scalar field or a self-interaction term and we study the impact of those parameterson the non-axisymmetric f -mode. The time evolution approach allows for a comparatively simpleimplementation of the boundary conditions. We find that the f -mode frequency is no longer asimple function of the star’s average density when a scalar field is present. We also evaluate theaccuracy of different variants of the Cowling approximation commonly used in previous studies ofneutron star oscillation modes in alternative theories of gravity and demonstrate that it can giveus not only qualitatively correct results, but in some cases also good quantitative estimates of theoscillations frequencies. ∗ [email protected] † [email protected] a r X i v : . [ g r- q c ] F e b . INTRODUCTION Gravitational wave asteroseismology is a powerful tool to study the internal structure ofneutron stars through the observed gravitational wave signal emitted once various oscillationmodes are excited by some astrophysical process [1–3]. For this purpose a large number ofstudies were performed examining the different classes of oscillation modes both for static[2, 4–7] and rotating stars [8–13], including spacetime modes [14, 15]. While the effect ofdifferent equations of state, composition, temperature profile, differential rotation, etc. wasstudied in detail, little has been done in quantifying another source of uncertainty, that isthe underlying theory of gravity.Within the framework of alternative theories of gravity, perhaps the most widely studiedneutron star models were those in a class of scalar-tensor theories (STT) of gravity admittingthe so-called scalarization. The reason is first that STT are some of the most natural andunproblematic extensions of general relativity (GR). Second, scalarization is a nonlineareffect allowing to endow highly compact objects with scalar hair while leaving the weak-fieldregime equivalent to GR and thus fully in agreement with the observations. Scalarizationof neutron stars was first considered in the original work of Damour and Esposito-Farese[16], extended later to slow [17–19] and rapid rotation [20]. Scalarized stars with anisotropicpressure were studied in [21] while the effect of magnetic field was examined in [22, 23]. Themost stringent constraint to date on this class of theories comes from the observations ofpulsars in close binary systems and they limit the possible deviations from GR to a verysmall value [17, 24–26]. One possibility to evade these constraints is to consider a nonzeroscalar field mass that introduces a characteristic radius of the scalar field, associated with itsCompton wavelength, beyond which the scalar field drops to zero exponentially [27–30]. Inaddition, a self-interaction term in the scalar field potential can have a qualitatively similareffect [31]. Neutron star mergers in such theories were considered in [32] while the core-collapse was examined in [33–36]. An alternative way to evade the binary pulsar constraintsis to consider for example multiple scalar fields [37].The first study of neutron star oscillations in alternative theories of gravity was performedfor scalarized neutrons stars in STT—the polar fluid modes were examined in [38] whilethe spacetime w -modes were considered in [39]. Calculating the polar perturbations is,in general, a much more involved task which is why the (which we will later dub “full”)2owling approximation was employed in [38], i.e. assuming that the spacetime as well asthe scalar field are fixed and only the fluid perturbations are evolved. Even though thisseems like a crude approximation, similar approaches have proven to give qualitatively goodresults in general relativity [9, 40, 41] and that is why it is reasonable to adopt it as afirst approximation to study the leading order effects of STT. The results in [38] were latergeneralized to the case of rapid rotation [42] where also the effect of the scalar field onthe Chandrasekhar-Friedman-Schutz [43, 44] instability was examined in detail. The polarmodes in the Cowling approximation in f ( R ) gravity, which is mathematically equivalent toa particular class of STT [45], were considered in [46].The field developed further in the direction of calculating the axial modes in alternativetheories of gravity [47–50]. Torsional oscillations of scalarized neutron stars were consideredin [51]. Another major effort was the calculation of radial modes which was first approachedby keeping the spacetime metric fixed but allowing for the evolution of the fluid and thescalar field [52]. The full problem without approximation was addressed in [53] where notonly the stability of neutron stars in STT against small perturbations was proved but alsothe emergence of a new class of modes associated with the scalar field was demonstrated.These are the breathing modes that can be excited in processes such as core-collapse [54].The first study of non-radial neutron star polar modes with (cid:96) ≥ without approximation inalternative theories of gravity was performed in [55]. There, the fundamental f -mode andits overtones were calculated for compact objects in a special class of scalar-tensor theorywith a massive scalar field which is mathematically equivalent to R gravity.In the present paper, we will concentrate on studying the polar oscillations modes, andmore specifically the (cid:96) = 2 f -mode, of scalarized neutron stars in STT. We focus on thecases both with and without scalar field potential, even though only the former one can giveus large deviations from GR if one considers values of the parameters in agreement with theobservations. The oscillation modes are calculated by evolving the relevant perturbationequations in time. Even though this method is inferior in accuracy compared to solving theeigenvalue problem, the treatment of the boundary conditions is considerably simpler andallows for straightforward calculations both in the case of massive and massless scalar field.The paper is organized as follows. In Sec. II the formulation of the problem and the basicequations are discussed. The perturbations equations governing the neutron star oscillationsin STT are given in a separate appendix. The results are presented in Sec. III. The paper3nds with Conclusions.Unless otherwise noted, we work in units in which c = G = 1 . II. MATHEMATICAL FORMULATIONA. Background neutron star solutions in scalar-tensor theories
The action of scalar-tensor theories is given in the Einstein frame by S = 116 π (cid:90) d x √− g (cid:2) R − g µν ∇ µ ϕ ∇ ν ϕ − V ( ϕ ) (cid:3) + S matter (cid:0) A ( ϕ ) g µν , χ (cid:1) , (1)where R and ∇ µ are the Ricci scalar and the covariant derivative with respect to the Einsteinframe metric g µν . V ( ϕ ) is the scalar field potential and S matter is the action of the mattersources collectively denoted by χ . The Einstein frame and the Jordan (physical) frame arelinked via the conformal factor A ( ϕ ) , with the Jordan frame metric ˜ g µν = A ( ϕ ) g µν . Due tothe simplicity of the field equations in the Einstein frame we will adopt it in our calculationsand refer to the physical Jordan frame mainly for the final observable quantities that wecalculate. Unless otherwise specified, the quantities in the Jordan frame will be markedwith a tilde. A detailed description of the two frames and the transformations between bothespecially in the case of scalarized neutron stars can be found for example in [20]The field equations resulting from the action (1) are R µν − g µν R = 8 πT µν + 2 ∇ µ ϕ ∇ ν ϕ − g µν g αβ ∇ α ϕ ∇ β ϕ − V ( ϕ ) g µν , (2) ∇ µ ∇ µ ϕ = − πα ( ϕ ) T + ∂V ( ϕ ) ∂ϕ , (3)where we have defined the coupling function α ( ϕ ) := d ln A ( ϕ )d ϕ . By using the contractedBianchi identities, we find that the conservation law for the energy-momentum tensor in theEinstein frame takes the form ∇ µ T µν = α ( ϕ ) T ∇ ν ϕ, (4)where T is the trace of the energy-momentum tensor.We will work with a perfect fluid neutron star, for which the energy-momentum tensor(in the Einstein frame) is, as usual, given by T µν = ( (cid:15) + p ) u µ u ν + pg µν , (5)4here (cid:15) is the energy density, p the pressure, and u µ the four-velocity of the fluid.In order to calculate the background neutron star models and their oscillation spectrum,we have to employ an equation of state ˜ p = ˜ p (˜ (cid:15) ) which will be provided in the physicalJordan frame. For this purpose we need relations between the fluid quantities in the twoframes. Given that T µν = A ( ϕ ) ˜ T µν , one can easily show that for a perfect fluid (cid:15) = A ( ϕ )˜ (cid:15) , p = A ( ϕ )˜ p , and u µ = A − ( ϕ )˜ u µ . Since the energy density and pressure are the only twofluid quantities which we need to know in the Jordan frame, we will transform only thosetwo to the Jordan frame when necessary but will otherwise work with the Einstein frame.In the present study, we will consider non-rotating neutron stars, for which the lineelement in isotropic coordinates can be written as d s = − e ν d t + e ψ d r + e ψ r dΩ (6)where ν and ψ are the two metric potentials and Ω is the solid angle. As a result, the four-velocity of the fluid is given by u µ = ( e − ν , , , . The dimensionally reduced field equationsassuming the above form of the metric, can be found in [56].We restrict ourselves to the study of the dynamics of small perturbations around anequilibrium configuration for which we have to linearise Eqs. (2) to (4). We will use thesame formalism that has been presented in Ref. [13]; nonetheless, not only for completeness,but also since we use merely the “nonrotational subset” of the formalism presented there(however, extended for the presence of a scalar field ϕ ), we will repeat the basics for clarity.We will introduce time-dependent perturbations for which we will derive evolution equa-tions. First, we decompose the metric as g µν = g (0) µν + h µν , (7)where g (0) µν is the background metric and h µν its perturbation; we use the background metricto raise and lower the indices of the latter. As we will employ the Hilbert gauge, it will beadvantageous to work instead with the trace-reversed metric perturbation, defined by φ µν := h µν − g (0) µν h, (8)where h := h µµ is the trace of the metric perturbations. The Hilbert gauge, which is thegravitational equivalent to the well-known Lorenz gauge in electromagnetism, is specified by f µ := ∇ ν φ µν = 0 . (9)5e opt for this gauge as it provides us directly with wave equations for the spacetimeperturbations, which can, without further ado, be used for a time evolution (see Ref. [13]for a more in-depth discussion). As we consider perturbations on a spherically symmetricbackground, we will separate out the angular dependence of the perturbation variables bymeans of the spherical harmonics (and their vector and tensor counterparts).The precise definition of the perturbation variables for the spacetime, the fluid and thescalar field can be found in Appendix A while the resulting set of evolution equations isstated Appendix B. B. Fixing the theory
Let us discuss the particular class of STT we will concentrate on. The free functions thatdefine the theory are the conformal factor A ( ϕ ) (or equivalently the function α ( ϕ ) ) and thescalar field potential V ( ϕ ) . The simplest assumption we can make for the coupling α ( ϕ ) isthat it can be expanded in series with respect to ϕ . Keeping only terms up to linear orderin ϕ and assuming that the cosmological background value of the scalar field ϕ = 0 , wehave α ( ϕ ) = α + βϕ. (10)The case with α (cid:54) = 0 and β = 0 is equivalent to the famous Brans-Dicke theory thatis, though, severely constrained by the weak field observations [57]. In the present paperwe will focus on the Damour-Esposito-Farése (DEF) model [16] with α = 0 and β (cid:54) = 0 that is perturbatively equivalent to GR in the weak field regime but can lead to nonlineardevelopment of the scalar field for highly compact objects such as neutron stars, that is theso-called scalarization. Such scalarization can happen for negative β and the more negative β is, the stronger is the scalar field (for scalarization in special cases with large positive β werefer the reader to [58–60]). Binary pulsar observations set very tight constrains on the lowerlimit of β , more precisely β > − . [17, 24–26]. This criterion is dependent on the equationof state, but taking into account that scalarization may happen for roughly β < − . [17],it is evident that there is little room for deviations from GR.As discussed in the introduction, the scalar field mass term (and possibly a self-interactionone) in the potential can help us evade these constraints by introducing a finite range ofthe scalar field of the order of its Compton wavelength. In accordance with [29, 31], we will6onsider the following form of V ( ϕ ) V ( ϕ ) = 12 m ϕ ϕ + 14 λϕ , (11)where m ϕ is the scalar field mass and λ the self-interaction parameter. Below, we use unitsin which c = G = 1 and the dimensionless scalar field mass and self-interaction parameterare defined as m ϕ → m ϕ R and λ → λR , where R = 1 . is one half of the solargravitational radius. C. Numerical Implementation
The numerical implementation of the evolution equations follows what is described in theprevious paper [13] but is considerably less involved since we are dealing here with a purelyradial problem. We will briefly recall the relevant details; for further information, we referthe reader to [13].The numerical treatment of the evolution equations remains essentially the same: weuse finite differences to discretize spatial derivatives on the same radial grid and employa 3rd order Runge-Kutta scheme for time integration; we use Kreiss-Oliger dissipation todamp out numerical instabilities [61]. The boundary conditions of the fluid perturbationsat the surface of the star and those of the space-time perturbations at the outer boundaryof the numerical grid are identical to those in [13]. We treat the scalar field perturbation atthe outer boundary identical to the space-time perturbation (by applying an outgoing-waveboundary condition).The boundary conditions at the origin deserve an extra comment since we have sepa-rated out the spherical harmonics from our perturbation variables (see App. A for theirdefinitions) and, subsequently, our evolution equations (see App. B) now involve the ordernumber (cid:96) rather than the azimuthal parameter m . After expanding the evolution equationsaround the origin r = 0 in a Taylor expansion, we find that the perturbation variables {H , L , M , Q , Q , Q , δϕ } have to be zero at the origin, while the other four variables {K , P , Q , W } need to have a vanishing radial derivative there. Our problem is radial in the sense that our time evolution equations have no explicit angular dependenceas it is completely separated out by means of the spherical harmonics. This should not be confused with radial modes which have actually no angular dependence themselves; our study is indeed concerned withnonaxisymmetric modes. . 5 1 0 1 0 . 5 1 1 1 1 . 5 1 2 1 2 . 5 1 30 . 811 . 21 . 41 . 61 . 822 . 22 . 42 . 6 M /M (cid:1) R c [ k m ] m (cid:2) (cid:1) = 0 , (cid:3) (cid:1) = 0 G R (cid:2) = - 4 . 5 (cid:2) = - 5 . 0 (cid:2) = - 6 . 0 (cid:2) (cid:1) = - 6 , (cid:3) (cid:1) = 0 m (cid:2) (cid:1) = 0 . 0 0 5 m (cid:2) (cid:1) = 0 . 0 5 (cid:2) (cid:1) = - 6 , m (cid:2) (cid:1) = 0 (cid:3) = 1 0 (cid:3) = 1 (cid:3) = 0 . 1 Figure 1. The mass as a function of the radius for all of the sequences we employed.
III. RESULTS
The quantities that will be presented in this section will of course be in the physicalJordan frame. The mass definition in scalar-tensor theories of gravity is a subtle problemand it turns out that only the Einstein frame ADM mass has natural energy-like properties(see, e.g., the discussion in [20] and references therein), and that is why we will use it forthe plots. The physical circumferential radius can be obtained as ˜ R c = R EF A ( ϕ ) , (12)where R EF is the Einstein frame radius at which the pressure vanishes. The oscillationfrequencies, on the other hand, are the same in both frames [42]. A. Background models
Let us first discuss the sequences of background neutron star models that we have con-structed in order to study their dynamics. Our main aim is to explore in detail the behaviorof the oscillation frequency with the change of the different parameters of the theory. Thisis complicated by the fact that we have three independent parameters coming only from8he STT, that is the coupling parameter β , the scalar field field mass m ϕ , and the self-interaction parameter λ . There is one more free input in the problem that is the equation ofstate (EOS). We have decided to limit our studies to only one modern realistic EOS in orderto make the presentation of the results more tractable, the SLy EOS [62]; in fact, we use itspiecewise-polytropic approximation [63] for the ease of numerical implementation. As far asthe STT parameters are concerned we have considered a large range of parameters and thesequences used for the time evolution are plotted in Fig. 1. For each of the sequences thepart that is before the maximum of the mass is plotted with a solid line, while the unstableone is plotted with a dotted line. The same convention is used also in most of the figuresbelow.The coupling parameter β spans the interval from − . , which is marginally in agreementwith binary pulsar observations, until − . While β < − . is not allowed by the observationsin the massless case, this changes dramatically for even very small scalar field mass or smallvalues of λ . The qualitative effect of nonzero m ϕ and λ is quite similar – the scalar field issuppressed leading to smaller deviations from GR. While λ (cid:54) = 0 preserves the position of thebifurcation points (for a fixed β ), the range of central energy densities where scalarizationoccurs shrinks with increasing scalar field mass m ϕ . As it was demonstrated in [28, 56, 64], m ϕ ∼ − eV (or m ϕ ∼ − in our dimensionless units) will lead to a complete suppressionof the scalar dipole radiation in the binary pulsar observations, while the correspondingneutron star models are practically indistinguishable in their bulk properties, such as massand radius, from the corresponding model with the same β but m ϕ = 0 . That is why thecase of m ϕ = 0 and λ = 0 results in the maximum possible deviation from GR for a fixedvalue of β (even though it is rigorously speaking in contradiction to observations). In thecase of nonzero scalar field potential, the presented sequences are for m ϕ (cid:54) = 0 and λ = 0 ,and for m ϕ = 0 and λ (cid:54) = 0 in order to better distinguish between effects coming from thetwo different terms in Eq. (11). B. Oscillation spectrum
We perform a time evolution of the perturbation quantities and place an observer at somelocation inside the star. Our simulations usually cover around and the interior of thestar is well resolved by 240 grid points. As we are interested in the fluid modes, we use as9nitial data a perturbation that roughly resembles that of an f -mode eigenfunction of one ofthe fluid quantities. The power spectral density (henceforth PSD) of the time series takenby the observer then reveals the f -mode as well as a few of its overtones (the pressure modes p n ). ∂ t δ ϕ P S D o f ∂ t δ ϕ f p p fully dynamicscalar Cowling Figure 2. Time series (left panel) and PSD (right panel) of the model with (cid:15) c = 1 . × g/cm , β = − , m ϕ = 0 , and λ = 0 . . Shown is the time derivative of the scalar field perturbation. Theblack PSD corresponds to the fully dynamic case (shown in the left panel), whereas the red PSD isobtained while keeping the scalar field fixed; this approximation slightly increases the fluid modefrequencies. The f -mode (at .
034 kHz and .
249 kHz , respectively) and its first two overtones p and p are clearly visible. As a first characteristic example, we discuss the simulation of a scalarized neutron starmodel with massless scalar field and non-zero self-interaction parameter λ ; we show thetime series of the time derivative of the scalar field perturbation, ∂ t δϕ , in the left panel ofFig. 2. The PSD of that time series is shown in the right panel (in black) and the f -modeand its first two overtones are clearly distinguishable (their frequencies are independent ofwhich perturbation variable we analyse); as we started with initial data resembling an f -mode eigenfunction, higher overtones at even higher frequencies are only weakly excited andnot as clearly visible. In the same panel, we show in red the PSD of the time series of thesame stellar model, however, this time while keeping the scalar field fixed (“scalar Cowling”),i.e., δϕ = 0 , that will be discussed in detail in the Sec. III C. The numerical values of thefrequencies are shown in the top half of Tab. I.10 able I. The frequencies of the f -mode and its first two overtones of the two models shown in Figs. 2(top half) and 3 (bottom half); see those captions for details on the models. The first row containsthe frequencies in the fully dynamic case, whereas the second row contains the slightly increasedfrequencies due to keeping a degree of freedom fixed (see text for further detail). All frequencies aregiven in kHz. The naming convention of the different Cowling approximations will be introducedin Sec. III C. STT parameters evolution type f p p m ϕ = 0 , λ = 0 . (cf. Fig. 2) fully dynamic 2.034 6.867 10.040“scalar Cowling” 2.249 6.918 10.085 m ϕ = 0 . , λ = 0 (cf. Fig. 3) fully dynamic 2.122 6.571 9.644“classic Cowling” 2.552 7.300 10.050 In Fig. 3, we display another exemplary time series and its PSD (in black), this time fora scalarized neutron star model in which the scalar field has non-zero mass but instead theself-interaction parameter λ of the potential is zero. Again, we show in red for comparisonthe PSD of the same model, however, this time with the space-time perturbations heldconstant and allowing for evolution of δϕ (“classic Cowling”), that will be discussed in detailin the following Sec. III C. The precise frequencies of the visible peaks can be found in thebottom half of Tab. I.Our extracted frequencies are accurate to 1-2%; we have established this accuracy of ourcode in our previous study [13] and also by running a few of the present simulations atdifferent resolutions. Furthermore, our results are in very good agreement with the resultsby Blázquez-Salcedo et al. [55].Next, we will turn to the f -mode frequency and study how it is impacted by the differentscalar field parameters. We start with the simplest case of keeping the scalar field potentialzero, i.e., m ϕ = 0 = λ , and vary only the conformal factor via β ; as explained before,this will lead to the largest deviations from pure general relativity ( m ϕ > or λ > will suppress the scalar field, see also Fig. 1). As the f -mode is an acoustic mode andhence its frequency is closely linked to the size of the star, we can—at least qualitatively—predict the change in frequency from the mass-radius diagram in Fig. 1; for simplicity, wefocus on the most extreme case shown in red with upward triangles (and the precise results11 ∂ t δ ϕ P S D o f ∂ t δ ϕ f p p fully dynamicclassic Cowling Figure 3. Time series (left panel) and PSD (right panel) of the model with (cid:15) c = 1 . × g/cm , β = − , m ϕ = 0 . , and λ = 0 . Shown is the time derivative of the scalar field perturbation. Theblack PSD corresponds to the fully dynamic case (shown in the left panel), whereas the red PSD isobtained while keeping the space-time fixed; the Cowling approximation increases the fluid modefrequencies. The f -mode (at .
122 kHz and .
552 kHz , respectively) and its first two overtones p and p are clearly visible. The wider Fourier peaks and the apparently increased energy stored inthe overtones is an artifact of the Cowling approximation. are shown in the left panel in Fig. 4): starting from the low density range of the mass-radius curve, the scalarization first pushes the equilibrium configurations toward smallerradii once we cross the first bifurcation point, hence the f -mode will increase in frequencywhen compared to the non-scalarized models. As we continue along the mass-radius curvetoward compacter models (despite the large deviation from the GR curve, the average densitykeeps monotonically increasing), the scalar field flips its effect and increases the neutronstar’s radius (and decreases its average density), coming along with a decrease in frequency;this continues up until the upper bifurcation point. Note that this description is of purelyqualitative nature and the approximately linear dependence of the f -mode frequency on the(square root of) star’s average density that has been observed in GR no longer holds inSTT. The effect of the scalar field is nearly negligible for β = − . but becomes increasinglypronounced for smaller β (we considered models down to β = − ).Let us now investigate the impact of a non-vanishing scalar potential. For this, we keep β = − fixed and vary the other two parameters m ϕ and λ ; in order to separate the effects12 . 0 3 0 . 0 4 0 . 0 51 . 61 . 822 . 22 . 42 . 6 0 . 0 3 0 . 0 4 0 . 0 50 . 0 3 0 . 0 4 0 . 0 5 Frequency [kHz] m (cid:2) (cid:1) = 0 , (cid:3) (cid:1) = 0 G R (cid:2) = - 4 . 5 (cid:2) = - 5 . 0 (cid:2) = - 6 . 0 ~/ c RM ~/ c RM (cid:2) = - 6 , (cid:3) = 0 G R m (cid:2) = 0 m (cid:2) = 0 . 0 0 5 m (cid:2) = 0 . 0 5 ~/ c RM (cid:2) = - 6 , m (cid:2) (cid:1) = 0 G R (cid:3) (cid:1) = (cid:3) (cid:1) = (cid:3) (cid:1) = (cid:3) (cid:1) = Figure 4. The f -mode frequency as a function of the average density (cid:113) M/ ˜ R c . (Left panel) m ϕ = 0 , λ = 0 and β is varied. (Middle panel) β = − , m ϕ = 0 and λ is varied. (Right panel) β = − , λ = 0 and m ϕ is varied. of those two, we will vary only one at a time.First, we will focus on the middle panel of Fig. 4 where we increase λ . There, we showthe f -mode frequency as a function of the average density of the star and the red line withtriangles displays the largest deviation from general relativity. As we increase λ , the f -mode frequencies move closer to the frequencies as obtained in pure general relativity; thedifference to pure GR is only marginal for λ = 1 and for λ = 10 it is indistinguishable in thediagram. This is due to the fact that an increased value of λ suppresses the scalar field andhence the scalarized equilibrium configurations are closer to the purely general relativisticcase than if λ was zero.Next, we consider the impact of a massive scalar field by increasing m ϕ (and setting λ = 0 again), cf. the right panel of Fig. 4. We observe a similar effect to that of a non-zero λ : with increasing scalar field mass, the frequencies move closer to the general relativisticvalues. While the impact of the scalar field’s mass is rather small for m ϕ = 0 . , the f -mode frequency is nearly indistinguishable from the GR values for m ϕ = 0 . . However, themassive scalar field also impacts the bifurcation points between which scalarized equilibriumconfigurations different from the GR solution exist. The more massive the scalar field is,the smaller the parameter window for scalarized neutron stars becomes.13t is natural to expect that an additional class of modes will emerge in the spectrum whichis associated with the scalar field and this class was indeed identified in studies concerningradial oscillations [53, 55, 65]. It is a challenging task to extract these modes from a timeevolution, though, since they have a rather short damping time and even if they are present inthe signal, they are damped already after one or two oscillation cycles to such low amplitudesthat they are essentially swallowed by the fluid oscillations. Particular initial data arerequired that excite the scalar modes to sufficiently large amplitude while keeping the fluidoscillations low so that they can be distinguished for a longer period of time; nonetheless,the coupling to the fluid will inevitably set it in motion, leading to the same intricacy asjust described. Even though we had some indications that such modes are present in oursimulations, we were not able to perform a clear and unambiguous identification of thosemodes at this time. Since the main focus of our paper is the calculation of the fluid polarmodes, we leave their further investigation for future studies. C. Accuracy of the Cowling approximation
At the end, we will discuss in detail the accuracy of the Cowling approximation for thefundamental (cid:96) = 2 f -mode oscillations of the neutron stars in order to know how accuratequalitatively and quantitatively the previous results on the subject are [38, 46] and to beable to give a rough prediction whether the Cowling approximation is justified for otheralternative theories theories of gravity, where the corresponding perturbation equations canbe much more complicated.The Cowling approximation has been extensively used for studies in GR (see, e.g., [9, 11,40, 41] and references therein) and it can overestimate the f -mode frequencies up to .The qualitative behavior of the modes does not change, though, only the frequencies areshifted (increased) with respect to the true results, even in the case of rapid rotation [11–13].There is no unique definition of the Cowling approximation in alternative theories ofgravity. The most straightforward (and simplest from computational point of view) is toassume that both the spacetime and the scalar field perturbations are zero. We will callthis “full Cowling” and this is the most studied case because it leads not only to simplerequations, but also the computational domain is limited to the interior of the star. Thisapproach was undertaken for example in [38, 46, 51]. One can go one step further and14 . 0 3 0 . 0 4 0 . 0 5 0 . 0 61 . 61 . 822 . 22 . 42 . 62 . 83 0 . 0 3 0 . 0 4 0 . 0 5 0 . 0 6
Frequency [kHz] (cid:2) (cid:1) = - 6 , m (cid:2) (cid:1) = 0 , (cid:4) (cid:1) = 0 ~/ c RM (cid:2) (cid:1) = - 6 , m (cid:2) (cid:1) = 0 . 0 5 , (cid:4) (cid:1) = 0 ~/ c RM G R e x a c t G R C o w l i n gS T T , D y n a m i c a l (cid:3) (cid:1) : f u l l y d y n a m i c a l c l a s s i c a l C o w l i n gS T T , F i x e d (cid:3) (cid:1) : s c a l a r C o w l i n g f u l l C o w l i n g
Figure 5. Comparison of f -mode frequencies for several series of models with varying (cid:15) c , wherethe space-time, scalar field or both or none of those are held fixed. The three blue curves have adynamic space-time while the red curves represent a frozen space-time. The two solid lines are thepurely general relativistic case ( β = 0 ), whereas the dashed (scalar field dynamic) and dash-dotted(scalar field fixed) curves are obtained with β = − (and m ϕ = λ = 0) . allow for the evolution of the scalar field while keeping the metric fixed. This approach wasperformed in [52] and we will name it “classic Cowling”. The last option is to keep the scalarfield fixed while evolving the metric, which we will call “scalar Cowling”.The comparison between the PSD in the case without approximation and for “scalarCowling” is presented in the right panel of Fig. 2 for two representative scalarized neutronstar models while the numerical values of the extracted frequencies are shown in Tab. I.Quite similar to the Cowling approximation in GR, freezing the scalar field to its equilibriumvalue also pushes the fluid frequencies to higher values. The impact on the frequencies is ingeneral weak, though, with the strongest deviation reached for the f -mode. A comparisonbetween the PSD in the case without approximation and for “classic Cowling” is presented inthe right panel of Fig. 3. As before, and as well known from the Cowling approximation inpurely general relativistic studies, the frequencies of the acoustic modes are shifted to largervalues. In contrast to the “scalar Cowling” case before, where we kept the scalar field fixedbut evolved the spacetime, the impact of this “classic Cowling” approximation is stronger.To visualize the effect of the different Cowling variants defined above, we show the f -mode frequency for two representative combinations of scalar field parameters in Fig. 5;15he pure GR values are shown with solid lines for comparison. The left panel representsthe case with β = − and vanishing scalar field potential, while in the right one we haverelatively high scalar field mass reducing significantly the deviations from GR and shrinkingthe domain of existence of scalarized solutions ( β = − , m ϕ = 0 . and λ = 0 ). Lookingat the right panel first, it is obvious that the Cowling variant which induces the smallestabsolute deviation from the fully dynamic situation (shown with dashed lines) is the “scalarCowling” approximation (dotted lines), as expected. Both the “classic Cowling” and the“full Cowling” approximations (both shown in red) in which the spacetime is held fixed yielda qualitatively very similar picture compared to the cases with a dynamic spacetime (inblack). To a good approximation the “fixed spacetime” (red) sequences are merely moved tohigher frequencies with respect to the “dynamical spacetime” branches. The shift appearsto be quite similar for the whole range of average compactness and it is of the same orderas the shift between the Cowling (solid red) and full (solid black) GR sequences.Next, in the case of vanishing scalar field potential (left panel of Fig. 5), the deviationsfrom GR induced by the presence of a dynamic scalar field are qualitatively similar in bothcases of a dynamic and static spacetime. There, the parameter window in which scalarizationoccurs is split into two parts: one of lower average density in which the f -mode frequency isincreased with respect to GR, while it is decreased in the one of higher average density. Thepicture changes slightly when we compare the GR case to the scalar Cowling approximation;here, the f -mode frequency of the scalarized neutron stars is increased compared to GR forall average densities; nonetheless, the feature of a local maximum and a local minimum inthe frequency still occurs. Again, the impact of the spacetime being frozen can be veryroughly approximated by a constant frequency shift with respect to the fully dynamicalresults and also in this case the shift is quite similar to the difference between the Cowling(solid red) and full (solid black) GR sequences.Thus, we can provide a simple recipe to estimate the frequency of the f -mode in the fullydynamic case even without evolving the scalar field or the spacetime (which is coupled tothe scalar field) in time: first, calculate the spectrum in the “full Cowling” approximation,which should in general a comparatively feasible task; second, subtract from those valuesthe difference between purely general relativistic value and its Cowling counterpart—it iswell-known how to calculate this quantity. This simple calculation should give a useful (evenquantitatively!) estimate for the f -mode frequency in the fully dynamic case.16t is highly probable, that similar conclusions will hold also for other scalar-tensor types ofgravitational theories. Since the field equations can be extremely complicated for such theo-ries, having the opportunity to obtain qualitatively good results and even some quantitativeestimation by adopting some form of a Cowling approximation, is of great value. IV. CONCLUSIONS
In the present paper we study the oscillation modes of spontaneously scalarized neutronstars. We consider scalar-tensor theories with and without scalar field potential, includingboth a massive and a self-interacting term. For this purpose the perturbation equationsgoverning the neutron star oscillations are derived without approximations contrary to theprevious studies on this subject. A numerical code is developed for the solution of theseequations in the time domain. Even though this approach is typically somewhat inferiorin accuracy compared to the solution of the equations in the frequency domain, it stillleads to quantitatively accurate results, where the accuracy is of the order of 1-2 %. Moreimportantly, it is much more tractable and straightforward to handle.We have calculated the oscillation frequencies of a large number of models spanning awide range of the free parameters; these are the STT coupling constant β , the mass ofthe scalar field m ϕ and the self-interaction parameter λ . For all of these models we havecalculated the fundamental (cid:96) = 2 f -mode and its first two overtones. As expected fromthe behavior of the background equilibrium solutions, the decrease of β leads to an increaseof the deviations in the oscillations frequencies compared to pure GR. Nonzero m ϕ or λ suppress the scalar field and for large values they lead to an oscillation spectrum that ispractically indistinguishable from GR.It is well known that in GR the frequencies of the f -modes scale linearly with the (squareroot of) average density of the star. We demonstrate that this is not fulfilled in scalar-tensortheories of gravity where the deviation from a linear dependence can be significant. It mightbe possible to introduce a new type of universal asteroseismology relations with a propernormalization of the parameters which holds also for scalarized neutron stars; such a studyis currently in progress.The previous studies in the field were done mainly with the use of some kind of Cowlingapproximation, assuming that the spacetime metric and/or the scalar field are fixed. The17pplicability and accuracy of this approximation was evaluated only in pure general relativitywhere it was shown that it overestimates the f -mode frequency by up to 30%; despite thisdeviation, it preserves the qualitative behavior of the f -mode frequency well. In scalar-tensor theories, different variants of the Cowling approximation can be defined and we havedemonstrated that all of them produce to some good extent the same qualitative behaviourof the f -mode frequencies. We find that even some quantitative estimates for the truefrequencies (i.e., in the fully dynamic case) can be made based on the results in the classicCowling approximation: the difference between a dynamic and a frozen spacetime is to goodapproximation independent of whether we consider a purely general relativistic neutron staror a scalarized neutron star (with or without scalar potential). We hypothesize, that suchtype of approximations can be very useful to track the behavior of the neutron star oscillationspectrum and the related gravitational wave emission in other alternative theories of gravity,where the solution of the full perturbation equations is considerably more involved. ACKNOWLEDGMENTS
We would like to thank Stoytcho Yazadjiev for the continuous support and guidance aswell as many fruitful discussions in the preparation of this work. DD acknowledges financialsupport via an Emmy Noether Research Group funded by the German Research Foundation(DFG) under Grant no. DO 1771/1-1 and the partial support by the National Science Fundof Bulgaria under Contract No. K-06-N-38/12. DD is indebted to the Baden-WürttembergStiftung for the financial support of this research project by the Eliteprogramme for Post-docs. SY would like to thank the University of Tübingen for the financial support. Thepartial support by the Bulgarian NSF Grant KP-06-H28/7 and the Networking support12 by the COST Actions CA16104 and CA16214 are also gratefully acknowledged. CKacknowledges financial support through DFG research Grant No. 413873357.
Appendix A: Definitions of the Perturbation Variables
In the derivation of the evolution equations, we follow closely the approach taken in [13],to which we refer the reader for further details. However, in the present study, we limitourselves to non-rotating stars which simplifies the expressions significantly.18he spherically symmetric perturbation problem naturally decomposes into two indepen-dent subsets: polar and axial perturbations. As we are interested in polar perturbationsonly, three of the ten metric perturbations vanish. We define the other seven to be φ tt = − e ν H , (A1) φ tr = L , (A2) φ tθ = r M , (A3) φ rr = 2 e µ K , (A4) φ rθ = e µ r Q , (A5) φ θθ = 2 e µ r P , (A6) φ ˆ ϕ ˆ ϕ = 2 e ψ r sin θ W , (A7)i.e., we have introduced the 7 perturbations H , K , L , M , P , Q , W . Note that we use the hatin order to distinguish the azimuthal spherical coordinate ˆ ϕ from the scalar field ϕ . Forthe fluid quantities, we use the components of the “Cowling part” δT µνC of the perturbedenergy-momentum tensor δT µν , both of which are related via δT µν =: δT µνC − ph µν . (A8)The fluid perturbations then are Q := δT tt C , (A9) Q := δT tr C , (A10) Q := δT tθ C r, (A11) Q := δT ˆ ϕ ˆ ϕ C , (A12) Q := δT rr C r sin θ. (A13)By means of these definitions, it can easily be seen that Q = Q (in the non-rotating case).Furthermore, we will see later that Q is not independent, either.The perturbation of the scalar field will simply be δϕ . Summarizing, we will need toevolve eleven perturbation quantities in time: seven for the spacetime, three for the fluidand one for the scalar field.We exploit the spherical symmetry of our problem to remove the angular dependence ofour perturbation variables and express their angular dependence in suitable combinations of19he spherical harmonics Y lm ( θ, ˆ ϕ ) . We use the subscript “ ” to denote the newly defined per-turbation variables that will depend on t and r only. We have the following decompositionsinto spherical harmonics X = X r (cid:96) − Y lm for X ∈ {L , Q , Q , Q , δϕ } , (A14)and X = X r (cid:96) − ∂ θ Y lm for X ∈ {M , Q , Q } . (A15)The remaining four perturbation variables H , K , P , and W are defined via the linearcombinations H = 14 (cid:2) H − K − P + (cid:96) ( (cid:96) + 1) W (cid:3) r (cid:96) − Y lm , (A16) K = 14 (cid:2) −H + K − P + (cid:96) ( (cid:96) + 1) W (cid:3) r (cid:96) − Y lm , (A17) P = − (cid:104) ( H + K ) Y lm + (cid:2) θ∂ θ Y lm + (cid:96) ( (cid:96) + 1) Y lm (cid:3) W (cid:105) r (cid:96) − , (A18) W = − (cid:104) ( H + K ) Y lm − (cid:2) θ∂ θ Y lm + (cid:96) ( (cid:96) + 1) Y lm (cid:3) W (cid:105) r (cid:96) − . (A19)The factor of r (cid:96) − in these definitions corresponds to the natural behaviour of the per-turbation variables at the origin r = 0 (as discovered via a Taylor expansion) and ensuresregularity of the solution there. Appendix B: The Evolution Equations
In writing down the evolution equations, we use the comma notation to abbreviate partialderivatives of the background quantities, for example, ψ ,r := ∂ r ψ , and we use the shortcut κ := 4 πe ψ ( (cid:15) + p ) . 20he Hilbert condition (cf. Eq. (9)) yields three equations: e ψ ∂∂t H = − (cid:18) ν ,r + ψ ,r + (cid:96)r (cid:19) L + 2 (cid:96) ( (cid:96) + 1) r M + e ψ ∂∂t K − ∂∂r L + 2 e ψ ∂∂t P − e ψ (cid:96) ( (cid:96) + 1) ∂∂t W , (B1) e − ν ∂∂t L = − (cid:18) ν ,r + (cid:96) − r (cid:19) H + (cid:18) ν ,r + 2 ψ ,r + (cid:96) + 22 r (cid:19) K − (cid:18) ψ ,r + (cid:96)r (cid:19) P − (cid:96) ( (cid:96) + 1) r Q + (cid:96) ( (cid:96) + 1)2 r (cid:0) ψ ,r + (cid:96) (cid:1) W − ∂∂r H + 12 ∂∂r K − ∂∂r P + (cid:96) ( (cid:96) + 1)2 r ∂∂ W , (B2) e − ν ∂∂t M = − r H − r K + (cid:18) ν ,r + 3 ψ ,r + (cid:96) + 1 r (cid:19) Q − ( (cid:96) + 2)( (cid:96) − r W + ∂∂r Q . (B3)Since the Hilbert gauge is a purely general relativistic construct, the Hilbert conditions areunaltered by the presence of a scalar field. These three Hilbert conditions may be used tosimplify the evolution equations or to monitor the violation of the Hilbert gauge throughoutthe time evolution.The perturbed field equations (cf. Eq. (2)) provide us with seven wave equations for thespacetime perturbations: e ψ − ν ∂ ∂t L = − πe ψ +2 ν Q − (cid:20) ϕ ,r + 3 ν ,r + 2 ψ ,r ν ,r + ψ ,r + 1 r ( ν ,r + ψ ,r )( (cid:96) + 1) + 4 (cid:96)r (cid:21) L − (cid:18) ν ,r + ψ ,r − (cid:96) − r (cid:19) ∂∂r L + ∂ ∂r L + 2 (cid:96) ( (cid:96) + 1) r (cid:18) ψ ,r + 1 r (cid:19) M + 2 e ψ ν ,r ∂∂t K − e ψ ν ,r ∂∂t H + 4 e ψ ϕ ,r ∂∂t δϕ (B4) e ψ − ν ∂ ∂t M = − πe ψ +2 ν Q + 2 r (cid:18) ψ ,r + 1 r (cid:19) L + (cid:20) ϕ ,r − κ − ψ ,r + 1 r (cid:0) ν ,r ( (cid:96) − − ψ ,r ( (cid:96) + 3) (cid:1) − (cid:96) − r (cid:21) M − (cid:18) ν ,r + ψ ,r − (cid:96) − r (cid:19) ∂∂r M + 2 e ψ ν ,r ∂∂t Q + ∂ ∂r M , (B5) e ψ − ν ∂ ∂t Q = 2 r (cid:18) ψ ,r + 1 r (cid:19) K − r (cid:18) ψ ,r + 1 r (cid:19) P − (cid:20) ϕ ,r + 2 κ + ν ,r − ψ ,r ν ,r + 5 ψ ,r − r (cid:0) ν ,r (cid:96) + ψ ,r ( (cid:96) − (cid:1) + 4 (cid:96) + 2 r (cid:21) Q + (cid:18) ν ,r + ψ ,r + 2 (cid:96) − r (cid:19) ∂∂r Q + ∂ ∂r Q + 2( (cid:96) + (cid:96) − r (cid:18) ψ ,r + 1 r (cid:19) W e − ν ν ,r ∂∂t M + 4 r ϕ ,r δϕ , (B6) e ψ − ν ∂ ∂t H = − πe ψ +2 ν Q − πe ψ Q − (cid:20) κ + 2 ν ,r − r ( ν ,r + ψ ,r )( (cid:96) −
2) + 4 (cid:96) − r (cid:21) H + (cid:18) ν ,r + ψ ,r + 2 (cid:96) − r (cid:19) ∂∂r H + ∂ ∂r H + 2 (cid:20) ϕ ,r − κ + ν ,r − ψ ,r − r ψ ,r (cid:21) K + 4 ν ,r (cid:18) ψ ,r + 1 r (cid:19) P + 2 ν ,r (cid:96) ( (cid:96) + 1) (cid:18) ψ ,r + 1 r (cid:19) W − e − ν ν ,r ∂∂t L + 4 e ψ ∂V ( ϕ ) ∂ϕ δϕ , (B7) e ψ − ν ∂ ∂t K = 8 πe ψ +2 ν Q − πe ψ Q + 2 (cid:20) ϕ ,r + ν ,r + ψ ,r − r ψ ,r (cid:21) H − (cid:20) ϕ ,r + 2 ν ,r + 4 ψ ,r − r (cid:0) ν ,r ( (cid:96) −
2) + ψ ,r ( (cid:96) − (cid:1) + 4 (cid:96) + 2 r (cid:21) K + (cid:18) ν ,r + ψ ,r + 2 (cid:96) − r (cid:19) ∂∂r K + ∂ ∂r K + 4 (cid:96) ( (cid:96) + 1) r (cid:18) ψ ,r + 1 r (cid:19) Q + 4 (cid:20) ϕ ,r + ψ ,r − ψ ,r ν ,r + κ − r (cid:0) ν ,r − ψ ,r (cid:1) + 1 r (cid:21) P − (cid:96) ( (cid:96) + 1) r (cid:20) ϕ ,r + ψ ,r − ψ ,r ν ,r + κ − r (cid:0) ν ,r − ψ ,r (cid:1) + 1 r (cid:21) W + 4 e − ν ν ,r ∂∂t L + 4 (cid:20) e ψ ∂V ( ϕ ) ∂ϕ + 2 (cid:96) − r ϕ ,r (cid:21) δϕ + 8 ϕ ,r ∂∂r δϕ , (B8) e ψ − ν ∂ ∂t P = 8 πe ψ +2 ν Q − πe ψ Q + 2 (cid:20) κ − ψ ,r ν ,r − r ν ,r (cid:21) H + 2 (cid:20) ϕ ,r + ψ ,r − ψ ,r ν ,r + κ − r (cid:0) ν ,r − ψ ,r (cid:1) + 1 r (cid:21) K − (cid:20) ψ ,r − r (cid:0) ν ,r ( (cid:96) −
2) + ψ ,r ( (cid:96) − (cid:1) + 4 (cid:96)r (cid:21) P + (cid:18) ν ,r + ψ ,r + 2 (cid:96) − r (cid:19) ∂∂r P + ∂ ∂r P + 2 (cid:96) ( (cid:96) + 1) (cid:18) ψ ,r + 1 r (cid:19) W + 4 e ψ ∂V ( ϕ ) ∂ϕ δϕ , (B9) e ψ − ν ∂ ∂t W = 4 r (cid:18) ψ ,r + 1 r (cid:19) Q + 1 r (cid:20) ( ν ,r + ψ ,r )( (cid:96) − − (cid:96) − r (cid:21) W + (cid:18) ν ,r + ψ ,r + 2 (cid:96) − r (cid:19) ∂∂r W + ∂ ∂r W . (B10)As can be seen, the presence of a scalar field leads to only small modifications of the pertur-bation equations when compared to the purely general relativistic case: a few coefficientsgain an additional ϕ ,r term and most equations recruit a source term due to the scalar fieldperturbation δϕ . 22he conservation of energy-momentum, Eq. (4), results in three evolution equations forthe hydrodynamics: e ν ∂∂t Q = − e ν (cid:18) ν ,r + 3 ψ ,r + (cid:96)r (cid:19) Q − e ν ∂∂r Q + e ν (cid:96) ( (cid:96) + 1) r Q −
12 ( (cid:15) + p ) (cid:18) ∂∂t H + ∂∂t K + 2 ∂∂t P − (cid:96) ( (cid:96) + 1) W (cid:19) + α ( ϕ )( (cid:15) − p ) ∂∂t δϕ , (B11) e ψ ∂∂t Q = − e ν (cid:2) ν ,r + ϕ ,r α ( ϕ ) (cid:3) Q + e ψ (cid:18) ϕ ,r α ( ϕ ) − ν ,r − ψ ,r − (cid:96) − r (cid:19) Q − e ψ ∂∂r Q −
12 ( (cid:15) + p ) (cid:20) ϕ ,r α ( ϕ ) + 2 ν ,r + (cid:96) − r (cid:21) H −
12 ( (cid:15) + p ) ∂∂r H − e − ν ( (cid:15) + p ) ∂∂t L + (3 p − (cid:15) ) (cid:20) (cid:96) − r α ( ϕ ) + ϕ ,r ∂α ( ϕ ) ∂ϕ (cid:21) δϕ + α ( ϕ )(3 p − (cid:15) ) ∂∂r δϕ , (B12) e ψ ∂∂t Q = − e ψ r Q − r ( (cid:15) + p ) H − e − ν ( (cid:15) + p ) ∂∂t M + 1 r α ( ϕ )(3 p − (cid:15) ) δϕ . (B13)The presence of a dynamic scalar field is reflected in these three equations in a very similarway as it is in the wave equations for the spacetime perturbations.Finally, the scalar field equation (3) yields an evolution equation for the perturbation ofthe scalar field: e ψ − ν ∂ ∂t δϕ = − πe ψ +2 ν α ( ϕ ) Q + 12 πe ψ α ( ϕ ) Q − (cid:104) πe ψ ( (cid:15) + p ) α ( ϕ ) + ν r ϕ ,r (cid:105) H + (cid:20) πe ψ α ( ϕ )(3 p − (cid:15) ) − e ψ ∂V ( ϕ ) ∂ϕ + ( ν ,r + 2 ψ ,r ) ϕ ,r + 2 r ϕ ,r (cid:21) K − ϕ ,r (cid:18) ψ ,r + 1 r (cid:19) P + (cid:96) ( (cid:96) + 1) ϕ ,r (cid:18) ψ ,r + 1 r (cid:19) W + (cid:34) πe ψ (3 p − (cid:15) ) ∂α ( ϕ ) ∂ϕ − e ψ ∂ V ( ϕ ) ∂ϕ + 1 r ( ν ,r + ψ ,r )( (cid:96) − − (cid:96) − r (cid:35) δϕ + (cid:20) ν ,r + ψ ,r + 2 (cid:96) − r (cid:21) ∂∂r δϕ + ∂ ∂r δϕ (B14)All perturbation equations are written in the Einstein frame. For brevity, we have alsokept the fluid quantities (cid:15) and p (and their combination κ ) in this frame, even though weneed to keep in mind that their Jordan frame equivalent is the one that is physically relevant.In addition to the hyperbolic differential equations, we are equipped with one more al-gebraic equation that allows us to compute Q out of some other perturbation quantities:23he definition of the speed of sound, δ ˜ p = ˜ c s δ ˜ (cid:15) , can (in terms of our perturbation variables)be expressed as Q = e − ψ ˜ c s (cid:2) e ν Q + ( (cid:15) + p ) H (cid:3) + 4 e − ψ (cid:0) p − ˜ c s (cid:15) (cid:1) α ( ϕ ) δϕ . (B15)The scalar field perturbation δϕ enters this relation since we have to translate the definitionof the speed of sound from the Jordan frame (where it holds) to the Einstein frame (in whichwe perform our time evolution); this results in the relation δp = ˜ c s δ(cid:15) + 4( p − ˜ c s (cid:15) ) α ( ϕ ) δϕ . [1] N. Andersson and K. D. Kokkotas, Mon. Not. Roy. Astron. Soc. , 493 (1998), arXiv:gr-qc/9706010.[2] N. Andersson and K. D. Kokkotas, Mon. Not. Roy. Astron. Soc. , 1059 (1998), arXiv:gr-qc/9711088.[3] K. D. Kokkotas and B. G. Schmidt, Living Reviews in Relativity , 2 (1999), arXiv:gr-qc/9909058 [gr-qc].[4] O. Benhar, V. Ferrari, and L. Gualtieri, Phys. Rev. D , 124015 (2004), arXiv:astro-ph/0407529.[5] H. K. Lau, P. T. Leung, and L. M. Lin, Astrophys. J. , 1234 (2010), arXiv:0911.0131 [gr-qc].[6] J. L. Blázquez-Salcedo, L. M. González-Romero, and F. Navarro-Lérida, Phys. Rev. D ,044006 (2014), arXiv:1307.1063 [gr-qc].[7] C. Chirenti, G. H. de Souza, and W. Kastaun, Phys. Rev. D , 044034 (2015),arXiv:1501.02970 [gr-qc].[8] J. A. Font, H. Dimmelmeier, A. Gupta, and N. Stergioulas, Mon. Not. Roy. Astron. Soc. ,1463 (2001), arXiv:astro-ph/0012477.[9] E. Gaertig and K. D. Kokkotas, Phys. Rev. D , 064063 (2008), arXiv:0809.0629 [gr-qc].[10] C. Krüger, E. Gaertig, and K. D. Kokkotas, Phys. Rev. D , 084019 (2010), arXiv:0911.2764[astro-ph.SR].[11] D. D. Doneva, E. Gaertig, K. D. Kokkotas, and C. Krüger, Phys. Rev. D , 044052 (2013),arXiv:1305.7197 [astro-ph.SR].[12] C. J. Krüger and K. D. Kokkotas, Phys. Rev. Lett. , 111106 (2020), arXiv:1910.08370[gr-qc].
13] C. J. Krüger and K. D. Kokkotas, Phys. Rev. D , 064026 (2020), arXiv:2008.04127 [gr-qc].[14] K. D. Kokkotas and B. F. Schutz, Mon. Not. Roy. Astron. Soc. , 119 (1992).[15] N. Andersson, K. D. Kokkotas, and B. F. Schutz, Mon. Not. Roy. Astron. Soc. , 1230(1996), arXiv:gr-qc/9601015.[16] T. Damour and G. Esposito-Farese, Phys. Rev. Lett. , 2220 (1993).[17] T. Damour and G. Esposito-Farese, Phys. Rev. D , 1474 (1996), arXiv:gr-qc/9602056.[18] H. Sotani, Phys. Rev. D , 124036 (2012), arXiv:1211.6986 [astro-ph.HE].[19] P. Pani and E. Berti, Phys. Rev. D , 024025 (2014), arXiv:1405.4547 [gr-qc].[20] D. D. Doneva, S. S. Yazadjiev, N. Stergioulas, and K. D. Kokkotas, Phys. Rev. D , 084060(2013), arXiv:1309.0605 [gr-qc].[21] H. O. Silva, C. F. Macedo, E. Berti, and L. C. Crispino, Class. Quant. Grav. , 145008 (2015),arXiv:1411.6286 [gr-qc].[22] J. Soldateschi, N. Bucciantini, and L. Del Zanna, Astron. Astrophys. , A44 (2020),arXiv:2005.12758 [astro-ph.HE].[23] J. Soldateschi, N. Bucciantini, and L. Del Zanna, Astron. Astrophys. , A39 (2021),arXiv:2010.14833 [astro-ph.HE].[24] P. C. C. Freire, N. Wex, G. Esposito-Farese, J. P. W. Verbiest, M. Bailes, B. A. Jacoby,M. Kramer, I. H. Stairs, J. Antoniadis, and G. H. Janssen, Mon. Not. Roy. Astron. Soc. ,3328 (2012), arXiv:1205.1450 [astro-ph.GA].[25] J. Antoniadis et al. , Science , 6131 (2013), arXiv:1304.6875 [astro-ph.HE].[26] L. Shao, N. Sennett, A. Buonanno, M. Kramer, and N. Wex, Phys. Rev. X , 041025 (2017),arXiv:1704.07561 [gr-qc].[27] D. Popchev, Master’s thesis, University of Sofia (2015).[28] F. M. Ramazanoğlu and F. Pretorius, Phys. Rev. D , 064005 (2016), arXiv:1601.07475 [gr-qc].[29] D. D. Doneva and S. S. Yazadjiev, JCAP (2016), 019, arXiv:1607.03299 [gr-qc].[30] R. Rosca-Mead, C. J. Moore, U. Sperhake, M. Agathos, and D. Gerosa, Symmetry , 1384(2020), arXiv:2007.14429 [gr-qc].[31] K. V. Staykov, D. Popchev, D. D. Doneva, and S. S. Yazadjiev, Eur. Phys. J. C , 586 (2018),arXiv:1805.07818 [gr-qc].[32] L. Sagunski, J. Zhang, M. C. Johnson, L. Lehner, M. Sakellariadou, S. L. Liebling, C. Palen- uela, and D. Neilsen, Phys. Rev. D , 064016 (2018), arXiv:1709.06634 [gr-qc].[33] U. Sperhake, C. J. Moore, R. Rosca, M. Agathos, D. Gerosa, and C. D. Ott, Phys. Rev. Lett. , 201103 (2017), arXiv:1708.03651 [gr-qc].[34] P. C.-K. Cheong and T. G. F. Li, Phys. Rev. D , 024027 (2019), arXiv:1812.04835 [gr-qc].[35] R. Rosca-Mead, C. J. Moore, M. Agathos, and U. Sperhake, Class. Quant. Grav. , 134003(2019), arXiv:1903.09704 [gr-qc].[36] C.-Q. Geng, H.-J. Kuan, and L.-W. Luo, Eur. Phys. J. C , 780 (2020), arXiv:2005.11629[gr-qc].[37] D. D. Doneva and S. S. Yazadjiev, Phys. Rev. D , 104010 (2020), arXiv:2004.03956 [gr-qc].[38] H. Sotani and K. D. Kokkotas, Phys. Rev. D , 084026 (2004), arXiv:gr-qc/0409066.[39] H. Sotani and K. D. Kokkotas, Phys. Rev. D , 124038 (2005), arXiv:gr-qc/0506060.[40] L. Lindblom and R. J. Splinter, ApJ , 198 (1990).[41] H. Sotani and T. Takiwaki, Phys. Rev. D , 063025 (2020), arXiv:2009.05206 [astro-ph.HE].[42] S. S. Yazadjiev, D. D. Doneva, and K. D. Kokkotas, Phys. Rev. D , 064002 (2017),arXiv:1705.06984 [gr-qc].[43] S. Chandrasekhar, Phys. Rev. Lett. , 611 (1970).[44] J. L. Friedman and B. F. Schutz, Astrophys. J. , 281 (1978).[45] T. P. Sotiriou and V. Faraoni, Rev. Mod. Phys. , 451 (2010), arXiv:0805.1726 [gr-qc].[46] K. V. Staykov, D. D. Doneva, S. S. Yazadjiev, and K. D. Kokkotas, Phys. Rev. D , 043009(2015), arXiv:1503.04711 [gr-qc].[47] J. L. Blázquez-Salcedo, L. M. González-Romero, J. Kunz, S. Mojica, and F. Navarro-Lérida,Phys. Rev. D , 024052 (2016), arXiv:1511.03960 [gr-qc].[48] J. L. Blázquez-Salcedo and K. Eickhoff, Phys. Rev. D , 104002 (2018), arXiv:1803.01655[gr-qc].[49] J. L. Blázquez-Salcedo, D. D. Doneva, J. Kunz, K. V. Staykov, and S. S. Yazadjiev, Phys.Rev. D , 104047 (2018), arXiv:1804.04060 [gr-qc].[50] Z. Altaha Motahar, J. L. Blázquez-Salcedo, B. Kleihaus, and J. Kunz, Phys. Rev. D , 044032(2018), arXiv:1807.02598 [gr-qc].[51] H. O. Silva, H. Sotani, E. Berti, and M. Horbatsch, Phys. Rev. D , 124044 (2014),arXiv:1410.2511 [gr-qc].[52] H. Sotani, Phys. Rev. D , 064031 (2014), arXiv:1402.5699 [astro-ph.HE].
53] R. F. Mendes and N. Ortiz, Phys. Rev. Lett. , 201104 (2018), arXiv:1802.07847 [gr-qc].[54] D. Gerosa, U. Sperhake, and C. D. Ott, Class. Quant. Grav. , 135002 (2016),arXiv:1602.06952 [gr-qc].[55] J. L. Blázquez-Salcedo, F. Scen Khoo, and J. Kunz, EPL , 50002 (2020), arXiv:2001.09117[gr-qc].[56] S. S. Yazadjiev, D. D. Doneva, and D. Popchev, Phys. Rev. D , 084038 (2016),arXiv:1602.04766 [gr-qc].[57] C. M. Will, Living Rev. Rel. , 4 (2014), arXiv:1403.7377 [gr-qc].[58] R. F. Mendes, Phys. Rev. D , 064024 (2015), arXiv:1412.6789 [gr-qc].[59] R. F. Mendes and N. Ortiz, Phys. Rev. D , 124035 (2016), arXiv:1604.04175 [gr-qc].[60] R. F. Mendes and T. Ottoni, Phys. Rev. D , 124003 (2019), arXiv:1903.11638 [gr-qc].[61] H.-O. Kreiss and J. Oliger, Methods for the approximate solution of time dependent problems ,GARP publications series, Vol. 10 (World Meteorological Organization, Geneva, 1973).[62] F. Douchin and P. Haensel, A&A , 151 (2001).[63] J. S. Read, B. D. Lackey, B. J. Owen, and J. L. Friedman, Phys. Rev. D , 124032 (2009).[64] V. I. Danchev and D. D. Doneva, Phys. Rev. D , 024049 (2021), arXiv:2010.07392 [gr-qc].[65] D. D. Doneva, S. S. Yazadjiev, and K. D. Kokkotas, Phys. Rev. D , 044043 (2020),arXiv:2005.02750 [gr-qc]., 044043 (2020),arXiv:2005.02750 [gr-qc].