{'fraction_non_alphanumeric': 0.0347065214199724, 'fraction_numerical': 0.013624218222391873, 'mean_word_length': 4.546002279758997, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
New perspectives on neutron star and black hole spectroscopy and dynamic tides

We elaborate on a powerful tidal interaction formalism where the multipole dynamics is kept generic and encoded in a linear response function. This response function is the gravitational counterpart of the atomic spectrum and can become of similar importance with the rise of gravitational wave astronomy. We find that the internal dynamics of nonrotating neutron stars admit a harmonic oscillator formulation yielding a simple interpretation of tides. A preliminary investigation of the black holes case is given. Our results fill the gap between Love numbers and dynamic tides.


I. INTRODUCTION
Analytic models for gravitational interaction of compact objects in General Relativity (GR) are plagued by potentially very complicated internal dynamics. Recent progress on such tidal interactions is mostly focused on nondynamical models [1][2][3], which in particular can not account for oscillation modes. This situation was already criticized and improved in [4,5]. However, it may be difficult to extend this approach to objects other than neutron stars (NS) and the internal dynamics is developed around the Newtonian limit (but the adopted GR corrections seem to be sufficient for most applications). Here we devise a substantially more powerful tidal interaction formalism based on an effective field theory (EFT) approach [6]. This approach was proposed in the context of black hole absorption [7] and consists in effectively replacing the extended object by a point particle comprising dynamic covariant multipolar degrees of freedom (DOF). In this article the dynamics of the multipoles is kept generic and encoded by a linear response function to external tidal fields.
Motion of extended bodies in General Relativity (GR) has been subject to question from the very beginning of the theory and gives the most important way to test gravity. Describing the dynamics of these objects is complicated and approximate methods have been developed, such as multipole expansion schemes along the lines of Mathisson, Papapetrou, and Dixon [8][9][10][11], between many others. However, the definition of covariant compactsource multipoles in GR according to Dixon is only useful for test bodies. The extension to self gravitating objects is not fully understood, though it is clear that Dixon's multipoles should be renormalized [12]. The adopted EFT approach implies a definition of covariant source multipoles of self-gravitating objects in GR. This definition is implicit until an explicit matching of the point-particle description to the actual extended object is worked out. This is the main purpose of the present work, but also highly nontrivial. (For instance, the determination of NS multipoles is already quite subtle for nonperturbed stationary spacetimes [13], see also [14].) This article substantially improves the situation for linear perturbations around a static (nonlinear) background such as NS, Black Holes (BH), White Dwarfs and likely even Boson Stars.
We illustrate our formalism with a simple Neutron Star (NS) model. The tidal constants (Love numbers and yet undetermined constants) are easily extracted from the response function. It turns out that, as long as linear perturbations are applicable, the internal dynamics of NS admits a formulation in terms of harmonic oscillator amplitudes [15][16][17] similar to the Newtonian case, making the multipolar DOF composite. This astonishing result leads to simple and intuitive interpretations of tidal interaction in GR. A full analysis of the black hole case is still in progress. Its outcome is hard to foresee and thus for sure will bear surprises.
Although numerical simulations capture the nonlinear aspects of tidal interactions, complementary analytic models stimulate invaluable (at least qualitative) interpretations of the physical processes at hand. An important aspect of our analytic dynamic tidal model is to naturally account for resonances between external tidal fields and oscillation modes of the NS in GR (see [16,[18][19][20][21][22][23][24][25]). Such resonances are of great importance. For instance, it was suggested recently that the oscillations excited by these resonances can be strong enough to shatter the NS crust, thus producing a weak short Gamma Ray Burst (GRB) [26] (more precisely, a weak precursor to the main flare of the GRB produced by the merger of the binary). Besides such spectacular effects, resonances can of course leave more subtle, but invaluable, imprints on the internal structure in the Gravitational Wave (GW) signal. Numeric relativity simulations suggest that oscillations excited by resonances can even be driven into the nonlinear regime and thus contribute significantly to GW [27].
The next revolution in GR will certainly arise from GW observatories like Advanced LIGO and VIRGO. These detectors will begin its operation soon and likely detect GW from binary NS mergers on a regular basis [28]. Such GW signals encode a tremendous amount of information on the internal structure of NS. This expectation is supported by recent numerical simulations, which reveal imprints of the equations of state [29,30] or the formation of a metastable hypermassive NS [31]. Simultaneous detection of GW and GRB can provide for the first time persuasive evidence for certain GRB scenarios [26,32].
The present article is a continuation of our work in [17] (on the Newtonian case) and we adopt notations and conventions therein.

II. EFFECTIVE ACTION
Our approach to account for the innumerable internal DOF of compact objects follows along the lines of thermodynamics. The obstacle is to identify state variables, which by definition describe the system on large scales (infrared, IR). In the case of gravitational interaction, these state variables reduce to the source multipole moments: At the same time they encode the IR field and the motion [11] of the object. An effective point-particle action along the lines of [6,7,33] is most natural to implement covariant multipoles as macroscopic variables, where m is the mass of the NS and E ab is the electric part of the Weyl tensor. For simplicity, we only discuss the covariant electric type quadrupole Q ab here, but inclusion of other multipoles is straightforward, see [33,Eq. (1)]. The indices a, b indicate the spatial components in a local Lorentz frame comoving with the NS. The worldline parameter τ is the proper time here. We consider linear perturbations of compact objects, so we expect a linear response of the quadrupole to the (quadrupolar) tidal field E ab , where the tilde denotes Fourier transformation from τ to ω, andF is the linear response function (or propagator). The main objective of the present article is to determinẽ F from a matching procedure. As explained in [17], from a Taylor-expansioñ the tidal constants µ 2 and µ ′ 2 emerge. The first parameter µ 2 is related to the dimensionless (relativistic, quadrupolar, 2nd-kind) tidal Love number k 2 = 3Gµ 2 /2R 5 , where R is the radius and G is the Newton constant, in agreement with definitions in [2]. Furthermore, µ ′ 2 parametrizes the tidal response beyond the adiabatic case. Though it was formally introduced in [34], it was not determined numerically yet. It obviously comes out as a byproduct within our approach. The constant λ is related to absorption [7], see also [35][36][37] and for a non-EFT treatment see, e.g., [38]. The time dependence of the mass parameter in the effective action due to absorption is discussed in [37]. The response function is analogous to the refractive index in optics, where imaginary parts also encode absorption. This analogy enlightens the matching procedure. Indeed, the phase shift between ingoing and outgoing waves encodes the real part of the response (Love number/refractive index) while the change in amplitude is due to absorption. But the nonlinear nature of GR makes the interpretation of phase shifts more subtle.
We should stress that besides encoding all quadrupolar tidal constants in a single function F , our approach can naturally accommodate the presence of oscillation modes that are obviously missed by a Taylor expansion (3). This possibility was not discussed in [7], where the focus is on absorption.
Generic extensions of the point-mass action were considered in [39] and the resulting EOM were related to Dixon's results. This can readily be applied to (1). Explicit expressions for the stress tensor in terms of Dirac delta distributions can be found in [40,41]. However, it should be noted that the relation between Dixon's covariant multipole moments and the covariant moments used in the action (1) is more of a formal nature when self-gravitating objects are considered.
It is straightforward to derive the Newtonian interaction potential for binaries belonging to the effective action. Even the first post-Newtonian (PN) correction for a generic quadrupole was already worked out [42] (though not from the effective action; see also [43] for the impact on GW). However, the dynamics of the quadrupole was essentially left open and only made explicit for the adiabatic case. The present work fills this gap by providing a dynamical quadrupole model. It should be stressed that even if the effective action is applied to Newtonian or PN approximations, the response functionF encodes strong field aspects of GR. This is the eminent advantage of the EFT approach. PN interaction potentials including tidal coefficients were derived in [34,42,44].

III. PERTURBED COMPACT OBJECTS
Without going into detail, we just mention here that perturbations of static compact objects can be determined from a system of coupled ordinary differential equations with the radial coordinate as the variable and the frequency ω entering as a parameter. (The specific case of nonrotating spherical symmetric NS perturbation goes back to [45,46], for reviews see [47][48][49].) In the exterior, the perturbation equations are given by the famous Zerilli [50] or Regge-Wheeler (RW) [51] equations for electric-or magnetic-parity type perturbations, respectively. These are the same equations that describe perturbations of Schwarzschild BH. The Zerilli equation can be cast into the (simpler) RW form [52], such that the discussion can be restricted to the latter. Moreover, the RW equation possesses analytic series solutions [53], see also [54,55], which are central for the present work. Our approach consists in solving the perturbation equations numerically in the interior and connecting to the analytic vacuum solutions by imposing appropriate boundary conditions at the surface.
As the RW equation is a second order homogeneous differential equation, its generic solutions can be represented by a linear combination of two independent solutions. In [53], the pair of analytic solutions X ν 0 and X −ν−1 0 represented by series of Gauss Hypergeometric Functions converges at the BH horizon but not at spatial infinity, while for X ν C and X −ν−1 C (series of Confluent Hypergeometric Functions) it is the other way around. Here ν is the renormalized angular momentum, which is fixed by requiring convergence of the analytic series solutions [53]. We review the needed elements of [53] in Appendix A 2. The solutions can be matched as where K ν is given by [53,Eq. (4.2)] or (A34).
In the absence of dissipation,F should be real so it is natural to work with manifestly real quantities. Before proceeding, we thus introduce normalization constants N ν such that the analytic solutions X ν N := N ν X ν C are real. We also require that the asymptotic amplitude is 1, i.e., N ν is uniquely defined by the requirement where r * = r + 2M log(r/2M − 1) is the tortoise coordinate, r is the Schwarzschild radial coordinate, and we identify M = Gm. It is straightforward to work out an explicit series representation for N ν from the formulas provided in [53,54]. The result is shown in Appendix A 2 c. Now the RW function X in the exterior can be decomposed as The main numeric result needed for our investigation is encoded in the amplitudes A 1 and A 2 . The introduction of ǫ 4 is suggested by an analysis of K ν and K −ν−1 for small ǫ.
We checked our implementation of the analytic solutions [53] against a direct integration method [56].

IV. EFT CALCULATION
The generic idea is to replace the compact object by an effective source encoded by the action (1), such that the RW function X at large distance (in the IR) is reproduced. This singular source can be expressed in terms of Dirac delta distributions. In the effective theory we therefore need to solve an inhomogeneous RW equation where l is the angular momentum quantum number. The source term S can be derived by projecting the stress tensor following from the action (1) onto tensor spherical harmonics, completely analogous to a point-mass source [57,58]. Explicit expressions for the quadrupole case l = 2 are supplied in Appendix A 3. The principle is the same for other values of l.
As the distributional source S should mimic the compact object, it must be located at r = 0. However, due to the (regular) singularity at r = 2M , the inhomogeneous RW equation (7) then does not seem to make sense. This problem is resolved by understanding the solutions as expanded in the post-Minkowskian expansion parameter M or ǫ = 2M ω. Expanding the solutions in ǫ is subtle due to various poles arising from Gamma Functions. If one keeps l generic and performs the limit l → 2 after the expansion, one ends up with a different set of solutions to the RW equation denoted X ν δ l and X −ν−1 δ l , where δ l = l − 2 represents the deviation from the quadrupole case l = 2. It holds Notice that the δ l -poles in the coefficients are canceled by poles contained in the solutions X ν δ l and X −ν−1 δ l . Despite these complication we work with X ν δ l and X −ν−1 δ l . The reason is that by keeping l generic one can easily identify the terms produced by the source S.
Still the solutions may be singular at r = 0, so regularization techniques are needed to handle the delta distributions contained in the source S. We choose a smooth ultraviolet "cutoff" in the form of a Riesz-kernel representation, see, e.g., [59]. Furthermore, as the calculation of the source S for generic l is not an easy one, we use a more ad hoc approach. We take S for l = 2 only and multiply the Riesz kernel by r −δ l to augment it with a fractional multipole character. Finally we represent δ(r) by where d = 3 is the number of spatial dimensions. µ 0 and c l are arbitrary parameters of unit inverse length formally introduced to make the expression dimensionally correct. The limit δ → 0 is understood to be taken in final expressions.
We are going to construct the solution to the inhomogeneous equation from the solution to the homogeneous one using the standard method of variation of parameters/constants. This allows us to reinterpret the analytic solutions X ν δ l and X −ν−1 δ l , which converge in the IR, as the ones belonging to the EFT. So the ansatz reads where C 1 and C 2 are yet unknown functions of r. If C 1 and C 2 were constants, then this would just be a solution of the homogeneous equation. It is an elementary result that where H 1 and H 2 are yet undetermined integration constants and W * is the Wronskian w.r.t. the tortoise coordinate, By virtue of Abel's identity, this Wronskian is actually a constant and can be evaluated at r = ∞ by analyzing the asymptotic behavior of the analytic solutions, see (A53). This procedure to solve the RW equation with a delta source has some similarities to the construction of the gravitational Green function, which was also obtained from analytic solutions very recently [60]. It can be interesting to further study this connection in the future.
Next one must constrain the integration constants H 1 and H 2 . The solution they represent is not allowed to correspond to further delta-type sources at r = 0. This argument leads to H 2 = 0, as the sourced terms are proportional to r −l−1 ∼ X −ν−1 δ l , while X ν δ l ∼ r l . If we would set l = 2 in the very beginning, then one must identify the constraint on H 1 and H 2 by iteratively solving the RW equation in ǫ. However, the coefficients in the numeric solution (6) scale differently by 4 orders in ǫ, so one must iterate to the same order in ǫ just for the leading order result. This is the reason for keeping l generic here, at the cost of introducing poles in δ l into the calculation. Dimensional regularization would solve this problem in a similar manner, see [61]. Interestingly in the static limit the dependence on the dimension can be absorbed into l [61].
It should be emphasized that the lower integration bound in (12) is r = 0, which is only possible if the integrand is understood as expanded in ǫ. Then the singularity of dr * /dr and the oscillatory behavior of the analytic solutions at the horizon are removed. Here nonconvergence of X ν δ l and X −ν−1 δ l at the horizon is actually not a flaw, but a feature. Furthermore, Eqs. (12) only start to depend on the upper integration bound r ′ at linear order in δ, so C 1 and C 2 are actually constant for ω → 0. Finally, the quadrupole components entering the source S must be determined according to (2). This requires to evaluate the tidal fieldẼ ab at r = 0, or where again the Riesz kernel provides the necessary regularization and the integrand is expanded in ǫ before evaluation. Note thatẼ ab (r, ω) follows from the RW master function (11) and the background. See Appendix A 3 d for its representation in RW gauge.

V. RESULTS
Comparing the numerically obtained solution (6) with the effective one (11) leads to two conditions. These are solved for the yet undetermined quantities H 1 andF in terms of A 1 and A 2 , completing the computation. At the end of the day, one arrives at Let us insist again that though this seems to be an expansion in ǫ, the numerical quantity A 2 /A 1 can still capture strong field effects from the interior. We define a renor-malizedF MS by dropping the poles in δ and δ l , analogous to minimal subtraction (MS) in dimensional regularization. We apply our method to the astrophysically most relevant case of NS. The used system of differential equations is derived in Appendix A 1. We use a simple polytropic equation of state (EOS) with index 1 for the nuclear matter. The results presented below are for a NS with m = 1.27m ⊙ and R = 8.85 km. The complex quasinormal mode frequencies for this specific NS model were reported in [47] and our numeric implementation reproduces them very well (except for the damping of the curvature modes).
An excellent fit forF MS for this NS model turns out to be provided we also fit the renormalization scale µ 0 . We numerically generated a set of 350 data points (with higher density near the poles) and the fit deviates from all of  (16) and (17). For the circles the logarithmic scale dependence was ignored (µ0 ∼ ω) and for the crosses only the leading order in (15) was taken into account. This shows that the ǫ 2 corrections in (15) are essential for a good fit, while the corrections from µ0 start to contribute only beyond the f-mode.
them by at most 2% (see Fig. 1). The optimal fit parameters are given by It is straightforward to infer the tidal constants defined by (3). It should be noted that all fit parameters are essentially independent of c l , which we varied from ω to µ 0 . It is remarkable that the quadrupole propagator can be approximated by a sum of response functions of harmonic oscillators, just like in the Newtonian case [17]. The relativistic case thus seems to admit an amplitude formulation [15,16] analogous to the Newtonian case. The a priori very complex internal dynamics of the NS is then approximated just by a set of harmonic oscillators, which are the more fundamental effective DOF composing the dynamical quadrupole. The constants q f R 3 /G and q p R 3 /G can be understood as GR versions of the overlap integrals. Resonances are quantitatively described by forced harmonic oscillators. An extension of this mechanical picture to nonlinear oscillators and/or mode coupling can offer a demonstrative phenomenological way to model even more realistic situations.
The frequency dependent Love number introduced recently in [62], where the formalism in [4,5] was used, should be related to our response function in the low frequency regime and a comparison is most interesting due to the very different setup (single object in perturbation theory vs. complete binary). It is further an interesting question whether the expansion (3) including the tidal coefficient µ ′ 2 introduced in [34] is enough to find agreement within the regime where a comparison is possible.
It is well known that the f-mode properties (frequency) basically just depend on the mean density of the NS [47].
To a good approximation this should also be true for the tidal response of the quadrupole, as we find here that it is largely dominated by the f-mode (at least for the adopted NS model). This aspect may be related to the universal relations discussed in [62][63][64]. Also the BH case can be readily investigated using (15). One can even work out analytic formulas for A 1 and A 2 from X BH ∝ X ν 0 + X −ν−1 0 [53], reading The nontrivial character of this analytic result becomes apparent once (A34), (A48), (A46), and (A47) are inserted. Because of the absorption due to the horizon one can not expect poles for ω ∈ R like in the NS case. This makes the analysis more complicated, as the whole complex plane must be considered. However, one can immediately obtain an expansion ofF in ǫ, This implies that the BH Love number µ 2 vanishes, in agreement with the findings in [2,3,61,65]. Unfortunately the unspecified parameter c l can substantially influence the next order tidal coefficient µ ′ 2 . This makes clear that for a rigorous investigation one should first redo the EFT calculation within a better regularization method like dimensional regularization.
But it should be stressed that the Love number is independent of c l . In fact, it is possible to obtain the leading order of (15) by setting δ l = 0 throughout the computation using a shortcut. Though l is not available to identify the constraint on H 1 and H 2 , a simple argument can be given at leading order. H 1 and H 2 must correspond to the linear combination of X ν N and X −ν−1 N that cancels the z −2 -term in the solution at orders ǫ 4 and ǫ 5 . Unfortunately at order ǫ 6 [corresponding to ǫ 2 in (15)] this approach breaks down due to z −2 log z contributions to the solution, which can only be interpreted by an iteration of the field equations. Still it is highly desirable to reach the next to leading order, as illustrated by Fig. 1.

VI. CONCLUSIONS
The validity of our results is supported in many ways. First of all, by comparing the adiabatic limit of the analytic solutions used here against the analytic zerofrequency solutions used in [2] one can show that the definition of µ 2 in [2] agrees with the definition through (3) and (15). Second, the term linear in ǫ in (19) was already derived in [7] and agrees with our findings. Finally, an intermediate result is the dependence of quadrupole componentsQ ab from (14) on the scale µ 0 , which in fact agrees with the beta function found in [33] using dimensional regularization (if it is assumed that c l is independent of µ 0 ). But unfortunately our current results depend on another parameter c l with unclear interpretation. This situation can be improved by applying dimensional regularization to the EFT calculation, which is also most useful for applications to post-Newtonian theory. But we expect our results to be good approximations for the dimensional regularized ones (the matching scale µ 0 will be slightly different). Still this calls for a clearer connection between our formalism and Ref. [33] (where also the background Schwarzschild metric is generated within the EFT). At the same time, higher multipoles (including magnetic/axial) should be treated as well.
However, these current problems with the regularization method play no role for the static limit. Thus no ambiguities for the definition of the Love numbers emerge. The predictions for the RW function from the EFT are "simply" matched to the numeric results, as in [61] for the black hole case in dimensional regularization. Consequently there is no need to interpret the definition of k 2 as relative to the BH case in the current approach. Furthermore, our computation is based on the exterior solution, so it is applicable to arbitrary (nonrotating) compact objects.
Another obvious next step is an application of our method to more realistic NS models. Besides realistic EOS, an investigation of the NS crust is most promising due to a possible connection to precursors flares in short GRB [26] by a shattering of the crust.
Further, realistic NS are rotating. Neutron star modes can become unstable in the rotating case (including the f-mode [66]), which can also give rise to violent astrophysical processes. An extension of the present method can be tried within a slow rotation approximation, see, e.g., [67,68] and references therein. The r-modes of rotating NS are of particular interest for resonances [16].
Analytic predictions for GW including tidal effects from an Effective One-Body (EOB) approach agree even quantitatively with numeric simulations [44,69]. Yet the tides are modeled by Love numbers and absorption [70] only. Inclusion of the dynamical multipole response function into the EOB formalism is expected to establish the impact of resonances between NS modes and orbital motion on GW in a reliable manner.
Our approach can probably also be evolved into a method for finding oscillation modes. In the conservative case, modes can be found by "just" integrating the perturbation equations for real frequencies. Estimates for the damping times can be obtained using the quadrupole formula (and basically correspond to the overlap integrals). However, if the modes are damped by, e.g., dissipative effects in the nuclear matter or mode coupling, then the poles of the response function should have a nonvanishing imaginary part. This also illustrates that our approach separates properties of the star from effects due to the surrounding spacetime, which is nontrivial as the background is nonlinear and does not admit superposition arguments.
The BH case is largely left unexplored for now. At the same time, the prospects are fascinating. If a mechanical oscillator model for BH can be formulated, then one can further elaborate on the thermodynamic analogy of these macroscopic DOF. This can lead to insights on macroscopic concepts like BH entropy and temperature from an EFT point of view. Analogies to the AdS/CFT correspondence discussed in [7,71] can probably be made more explicit, too. The highly damped modes and branch cuts (eventually introduced by nonanalytic terms like log ǫ) can be difficult to handle. Extension to the case of rotating BH should be almost straightforward, as the perturbation master equations [72], the appropriate analytic solution [73], and the effective action [35] are readily available. The Love number of rotating BH should even come out unambiguously if the current (improvable) regularization method is applied.
For an application to BH scattering a simplistic fit of the response function can be accurate enough. This can lead to interesting connections to, e.g., the scattering thresholds discussed in [74,75]. Finally, if the present method is applied to perturbations of massive scalar fields around rotating black holes, the floating orbits existing for extreme mass ratios [76] can possibly be constructed for comparable mass binaries using PN methods.
Appendix A: Formulas and Implementations

a. Preliminaries
A static spherically symmetric star configuration is described by a metric where dΩ is the line element on the unit sphere, and a perfect fluid with stress tensor is given by where ρ is the density, P is the pressure, g µν is the metric components, u µ is the four-velocity such that u µ u µ = −1, and Greek indices belong to the spacetime coordinate basis. In the NS case, the description is usually completed by a barotropic EOS relating P and ρ (neglecting temperature), according to for a given functionρ. In this article we considered a polytropic EOS defined by The dynamics of the compact object is then given by the Einstein equations and the conservation equation of the stress tensor The even parity metric perturbations around a spherically symmetric background are given in the Regge-Wheeler gauge by where the functions h ∈ {h 0 , h 1 , h 2 , k} depend on the coordinates x µ according to h = h(r) exp(iωt)Y lm (Ω), Y lm being the scalar spherical harmonics. The total metric then reads ds 2 = ds 2 0 + ds 2 1 , where the subscripts 0 and 1 denote the background and perturbation, respectively. The perturbations to the matter fields are given by Given an equation of state, the perturbation to the density is given in terms of the perturbation to the pressure: which defines the speed of sound c s . In the following, it will be useful to introduce the function U such that The solution for u µ 1 is given in terms of the function U and of the metric perturbation where dΩ i = (dθ, dϕ).
This transformation has the advantage to turn the vacuum equation (U = 0) to the simple Regge-Wheeler equation (in tortoise coordinates) where r * is the tortoise coordinate defined by Note that the Regge-Wheeler equation originally describes the odd sector of the metric perturbations. The system of equations describing the even sector is called the Zerilli equation. However, these sectors are isospectral and are actually equivalent, by virtue of the transformation (A14). Here we use the Regge-Wheeler equation since the form of it is of Heun's equation and is better suited for the construction of the analytic series solutions of [53].

c. Boundary Conditions and Series
Within the NS interior it is difficult to extend the numeric integration up to the boundary points r = 0 and r = R. This problem is solved by terminating the numeric integration very close to these points and use analytic series solutions to extend the numeric solutions to the boundaries. Furthermore, certain boundary conditions must be fulfilled, which are directly implemented into the series solutions here.
At the origin r = 0, the relevant boundary condition is simply the regularity of the perturbation master functions. This imposes two independent conditions, so we need only two (of four) integration constants to parametrize the solution around r = 0. The leading order series solutions read where U 0 and k 0 are the integration constants.
The boundary condition at the surface is given by the requirement that the Lagrangian (comoving) perturbation of the pressure or density vanishes for r = R, or explicitly However, the behavior of the perturbations near the surface r = R crucially depends on the EOS. In the following analysis, we restrict to the case that the EOS near the surface is a polytrope with index 1 ≤ n < ∞. Then the boundary condition is actually equivalent to regularity of the perturbation master functions. The boundary condition allows us to eliminate one of the four integration constants, so we are left with U (R), k(R), and h 1 (R). At the surface the functions k and h 1 must be continuous, which provides two further boundary conditions. We are therefore able to express k(R) and h 1 (R) in terms of the RW function X(R) and its derivative X ′ (R) (which describe the exterior perturbation). The boundary series finally reads where U R , X R , X ′ R are the functions U, X, X ′ evaluated at the radius R. If the RW equation is solved numerically using a direct integration method [56], then one must derive series solutions for the RW function X at r = ∞ in a similar manner (and for BH also at the horizon, where the physical boundary condition only permits an ingoing flux).
The relation between the solutions is given by [53, Eq. (4.1)] where [53, Eq. (4.2)] (A34) The value of r ∈ Z is in principle arbitrary, which can also be checked numerically. For definiteness, we chose r = 0.
We proceed along the lines of [53,55] provided that the continued fractions converge in the specified limit. The corresponding solution to the three-term recurrence relation is called the minimal solution in the specific limit and is guaranteed to exist. But the minimal solutions for n → ∞ and n → −∞ are not necessarily the same, e.g., in general one can fulfill either (A41) or (A42), but not both at the same time. However, requiring both (A41) and (A42) fixes the renormalized angular momentum ν. This is dictated by the convergence of the analytic solutions to the RW equation.

c. Normalization of Analytic Solutions and Wronskian
It is straightforward to determine the asymptotic behavior of X ν C as with the complex amplitudes (2ǫ) iǫ a ν n Γ(n − iǫ + ν − 1)Γ(n − iǫ + ν + 1) Γ(n + iǫ + ν + 1)Γ(n + iǫ + ν + 3) . (A47) Then we obtain for the normalization Notice that N ν is multivalued (bi-valued). Eventually the second root must be used in some frequency regimes in order to make results continuous. In the present investigation this is necessary at about ωR/2π > 0. 23. It follows that = 1 2 e −i(ωr * +αν ) + e i(ωr * +αν ) , with α ν : = cos(ωr * + α ν ), as envisaged. Because the RW equation has real coefficients, it is guaranteed that a real solution for α ν exists. It is also straightforward to obtain the Wronskian (13) from this analysis as The Wronskian (13) based on the solutions X ν δ l and X −ν−1 δ l follows by expanding (A53) in ǫ for generic l.

a. Inhomogeneous RW Equation
The homogeneous Zerilli and Regge-Wheeler equations describe vacuum perturbations. Considering additional matter fields will source the vacuum perturbation equations. If the right hand side of the perturbed Einstein equation