Fast Rotating Neutron Stars: Oscillations and Instabilities

In this review article, we present the main results from our most recent research concerning the oscillations of fast rotating neutron stars. We derive a set of time evolution equations for the investigation of non-axisymmetric oscillations of rapidly rotating compact objects in full general relativity, taking into account the contribution of a dynamic spacetime. Using our code, which features high accuracy at comparably low computational expense, we are able to extract the frequencies of non-axisymmetric modes of compact objects with rotation rates up to the Kepler limit. We propose various universal relations combining bulk properties of isolated neutron stars as well as of binary systems before and after merger; these relations are independent of the true equation of state and may serve as a valuable tool for gravitational wave asteroseismology. We also present an introductory example using a Bayesian analysis.


I. INTRODUCTION
Oscillations and instabilities of neutron stars were always considered among the promising sources for gravitational waves. The systematic study of non-axisymmetric neutron star oscillations began in the 1960s with the pioneering works of Thorne and collaborators (cf. [1,2] (and subsequent papers)), in which they laid out the equations governing the perturbations of compact stars in general relativity. The numerical solution of these equations has proven highly challenging and it has taken nearly two decades before Lindblom and Detweiler found an advantageous formulation of the eigenvalue problem that allowed them to determine the complex frequencies of the acoustic modes of sufficiently realistic stellar models ( [3,4]). These results did not conclude the investigation of these modes; in particular [5,6] turned to the perturbations of relativistic stars and studied their oscillations as a scattering problem.
In the mid-1980s, inspired by a toy model [7] suggested that the dynamic spacetime of a neutron star exhibits its very own class of modes and christened them w-modes ( [8]). The short damping times of these modes (comparable to the ones of black-holes) pose numerical challenges, but the application of the continued fraction method ( [9]) and numerical integration along anti-Stokes lines ( [10]) improved considerably the accuracy in the calculation of both fluid and spacetime modes and led to the discovery of the new class of w-modes, the so-called w II -modes ( [11]). This advancement in computational accuracy allowed the study of many neutron star equations of state leading to the discovery of universal relations involving frequencies and damping times of oscillation modes. These universal relations provided ways for solving the inverse problem within the so-called gravitational wave asteroseismology, allowing to set constrains on the mass and/or radius of the neutron star, and eventually the nuclear equation of state, via the observation of oscillations (see, e.g., [12][13][14][15]).
With increasing computational power during the 1990s, the first attempts were made to increase the dimensionality of the hitherto (due to the restriction to spherical symmetry) purely radial problem and include time as a second dimension. The first successful time evolutions of perturbations of relativistic neutron stars were reported by [16]. [17,18] reformulated those evolution equations by means of the ADM-formalism ( [19]) and introduced a non-uniform radial grid to ensure stable numerical evolution when using realistic equations of state .

A. Asteroseismology
The different oscillation patterns of a neutron star are characterized by their restoring force, e.g., p(pressure)-modes, g(gravity)-modes, i(Coriolis)-modes, s(shear)-modes or w(wave)modes. The f -mode is the fundamental mode of the p-mode sequence and it is the oscillation mode most likely to be excited in violent processes such as neutron star formation by supernova core collapse ( [20,21]), the pre-merger interaction of neutron stars (see, e.g., [22][23][24][25][26][27][28][29][30]), the early-post-merger oscillations of the final object ( [31][32][33][34][35][36][37][38][39][40][41]). In the case that the merging neutron stars are of relatively small mass and the post-merger object is a fast spinning neutron star, the unstable f -mode oscillations can lead to its spin-down ( [42]). The f -mode is associated with major density variations and thus can potentially be an emitter of copious amounts of gravitational radiation. The emission of gravitational waves is the primary reason for the mode's rapid damping, at least for newly born neutron stars.
The efforts to associate the patterns of oscillations with the bulk parameters of the stars, e.g., their mass, radius or equation of state (henceforth EoS) was initiated in the mid-1990s and continued for almost two decades, advancing the field of gravitational wave asteroseismology ( [12-15, 43, 44]). To date, very robust empirical relations have been derived for non-rotating neutron stars, connecting observables such as frequency, damping time, or moment of inertia I to the bulk properties; for example, relations of the form [12,44]) could provide the average density or the moment of inertia of the star if the f -mode frequency σ 0 is known.
In the era of gravitational wave astronomy, the various oscillation patterns (traced already in non-linear numerical simulations, e.g., [45][46][47][48][49][50][51]), if observed, can provide a wealth of information about the emitting sources, and their effects can leave their imprints both in the gravitational but also in the electromagnetic spectrum. Moreover, recent studies relate the f -mode frequencies to the Love numbers ( [52][53][54][55]) and even to the postmerger short gamma-ray bursts ( [56]).

B. Rotation
Nearly all the previously mentioned studies were concerned with non-rotating stars. In nature, neutron stars will always rotate and their rotation rate may reach extreme values.
The inclusion of rotation proves difficult since the extreme rotation rates that neutron stars may (and do) reach do not allow to neglect the star's oblateness which removes the spherical symmetry from the system; this in turn makes the mathematical formulation much more involved. As a first approximation, rotation was treated perturbatively, too, which allowed considering even rotating stars as spherically symmetric ( [57]). In this so-called slow-rotation approximation, the perturbation equations gain considerably in complexity and have been written down in a gauge introduced by [58] first by [59]. Even though the problem remains one-dimensional, its solution is not straightforward as (among other technicalities) the outgoing-wave boundary condition at infinity is elusive. Notwithstanding, [60] successfully applied this formalism and discovered that r-modes are prone to the so-called CFS-instability, named after their discoverers [61] and [62,63], at any rotation rate. There has been continuing effort using the slow-rotation approximation concerning rotational modes ( [64,65]), and also employing different gauges ( [66][67][68]), but with the pressing need for frequencies of rapidly rotating neutron stars, the interest slowly faded.
Even though the slow-rotation approximation has proven fruitful in the understanding of neutron star physics, it is no longer applicable when considering neutron stars at arbitrary rotation rates, which is essential for nascent neutron stars or post-merger configurations in the immediate aftermath of a binary merger. Without the spherical symmetry of the problem, one has to account for at least two spatial dimensions which complicates the equations further and amplifies the computational expense; furthermore, it remains elusive how to formulate the outgoing-wave boundary condition at infinity for the spacetime perturbations in two spatial dimensions which essentially removes the possibility to formulate a corresponding eigenvalue problem. This issue can be circumvented by adopting the Cowling approximation ( [69]), in which the spacetime is considered static, also leading to a considerable simplification of the perturbation equations. Ignoring the impact of a dynamic spacetime (which is most severe for the quadrupolar f -mode), [70,71] computed quadrupolar f -mode frequencies of rapidly rotating neutron stars and studied the associated CFS-instability in the late 1990s. [72] revived this approach and investigated the general properties of the spectrum of neutron stars regarding the acoustic and Coriolis-driven modes. As a further step toward a more general relativistic treatment, [73] revisited the problem in the conformal flatness approximation. However, with the mathematical difficulties of extending the eigenvalue formulation to include a dynamic spacetime, the focus shifted to study the oscillation spectra by evolving the perturbation equations in time.
Despite their complexity, the perturbation equations for rapidly rotating relativistic stars have been written down by [74], even though they were not approached numerically at that time. [75,76] worked in the Cowling approximation and successfully extracted f -and g-mode frequencies of arbitrarily uniformly and fast rotating neutron stars and even of differentially rotating ones [77] by adding artificial viscosity (also known as [78] dissipation) to their evolution equations in order to stabilise their time evolutions.
During the first decade of the new millennium, substantial advances were made in the time evolution of the unperturbed, non-linear Einstein equations, mostly driven by the aim to simulate compact binary mergers but also applicable to isolated neutron stars. These systems have hardly any symmetries that can be exploited to reduce the complexity of the problem, requiring to carry out the time evolutions on a three-dimensional grid. The upside of which is that essentially no constraints have to be placed on the rotational profile when simulating the dynamics of a neutron star. Such codes have been seen as a promising new approach to the calculation of mode frequencies of rapidly (and differentially) rotating neutron stars and already at the beginning of the decade, the frequencies of axisymmetric modes in the Cowling approximation ( [79,80]) and those of (quasi-)radial modes in full general relativity ( [81]) had been reported. The non-linear codes kept evolving and were used to generate mode frequencies of f -modes in the conformal flatness approximation ( [82]) or those of inertial modes in the Cowling approximation ( [83]). [48]. Not much later, the frequencies of non-axisymmetric modes in full general relativity of non-rotating polytropic neutron stars ( [48]) and soon those of rapidly rotating polytropic neutron stars ( [49]) were obtained from fully non-linear simulations. Even though successful, this approach to computing the frequencies of non-axisymmetric modes, however, has not been followed closely, which is also due to the computational expense associated with such numerical simulations and the accompanying limited accuracy.
In fact, from the point of view of gravitational wave detectability of oscillation modes, the most relevant scenarios are likely to involve rapidly rotating stars. Unfortunately, the aforementioned empirical relations cannot be trivially extended to rotating stars. Rotation splits the oscillation spectra in a similar fashion as the Zeeman splitting of the spectral lines due to the presence of magnetic fields. In rotating stars, the splitting leads to perturbations propagating in the direction of rotation (so-called co-rotating modes) and perturbations traveling in the opposite direction (counter -rotating modes). The oscillation frequency as observed by an observer at infinity will either increase or decrease depending on the propagation direction of the waves; for slow rotation there will be a shift of the form where m is the angular harmonic index, κ a mode and stellar model-dependent constant, and Ω the angular rotation rate of the star. If the spin of the star exceeds a critical value, which depends on, e.g., the EoS and its mass-i.e., when the pattern velocity σ/m of the backward moving mode becomes smaller than the star's rotation rate Ω/2π-then the star becomes unstable to the emission of gravitational radiation; this is the aforementioned CFS instability. This instability is generic (independent of the degree of rotation) for the r-modes ( [60,84]) while it can be excited only for relatively high spin values (Ω 0.8Ω K , with Ω K the Kepler velocity) for the quadrupolar f -modes. An extensive discussion can be found in [85] and [86].
This review is based on the highlights of four recent articles published by the authors, which are [55,[87][88][89].
Throughout this article, we employ units in which c = G = M = 1.

A. Background configuration
We are going to work with the Einstein equations along with the law for the conservation of energy-momentum, where G µν is the Einstein tensor and T µν is the energy-momentum tensor.
We restrict ourselves to the study of the dynamics of small perturbations around an equilibrium configuration which allows us to linearise equations (1). We assume an axisymmetric, stationary background configuration for which the metric written in quasi-isotropic coordinates takes the form = −e 2ν dt 2 + e 2ψ r 2 sin 2 θ(dϕ − ωdt) 2 + e 2µ (dr 2 + r 2 dθ 2 ).
Here, ν, ψ, µ, and ω are the four unknown metric potentials, depending only on r and θ.
We model the neutron star to be a perfect fluid without viscosity for which the corresponding energy-momentum tensor takes the form where is the energy density, p is the pressure, and u µ the 4-velocity of the fluid. The only two non-vanishing components of the 4-velocity are linked via the star's angular rotation rate, u ϕ = Ωu t , and by means of the normalisation of the 4-velocity they are given by After specifying an EoS, which may be a polytropic or a tabulated one, linking energy density and pressure to each other, we generate uniformly rotating equilibrium configurations using the rns-code ( [90][91][92]).
We considered sequences of neutron stars along which we keep either the central energy density constant or the baryon mass fixed (the latter are also known as evolutionary sequences) with rotation rates up to their respective mass-shedding limit. Our neutron star models were based on polytropic and realistic EoSs. We considered polytropic models with three

B. Perturbation Equations
As usual in perturbative studies, we decompose the metric as where g (0) µν is the background metric and h µν its perturbation. As we will work in the Hilbert gauge, it will be advantageous to work instead with the trace-reversed metric perturbation, defined by where h := h µ µ is the trace of the metric perturbations. The metric perturbations are not unique but possess gauge freedom which can be utilised in different ways. Often, the gauge freedom is used to eliminate some of the spacetime perturbations, e.g., by using the well-known gauge by [58], and hence to reduce the number of perturbation equations. In our studies, however, we followed a different approach (which we will reason below) and opt for the Hilbert gauge, which is the gravitational equivalent to the well-known Lorenz gauge in electromagnetism, specified by In the Hilbert gauge, the perturbed Einstein tensor takes the form where R α µ β ν , R αβ , and R are the background Riemann tensor, Ricci tensor and scalar curvature, respectively. The advantage of the Hilbert gauge is that the evolution equations for the metric perturbations will take the form of ten coupled wavelike equations (note that in the above expression, the d'Alembert operator, defined with respect to the background metric, is the only differential operator acting on the metric perturbations) while the mixing of temporal and spatial derivatives is avoided. This is in contrast to other common gauge choices and the perturbed law for the conservation of energy-momentum The evolution equations themselves are quite lengthy and not very enlightening and we refer the reader to [88] for details in the derivation and their implementation.

A. Universal Relations for single neutron stars
As shown in [88], our code produces results in excellent agreement with previously published values ( [27,49]) and our convergence tests demonstrate an accuracy of the obtained frequencies of 1 − 2%. In this article, we will provide some highlighted results in order to demonstrate the existence of asteroseismological relations of various types and we lay out the way that one can make use of these relations in analyzing gravitational wave signals.
More specifically, we will show different universal relations providing accurate estimates for the f -mode frequency given some bulk parameters of the star and vice versa. First, we observe a universal behavior of the f -mode frequency σ i as observed in the inertial frame as a function of the star's angular spinning frequency Ω along sequences of fixed central energy density models when we normalize both frequencies with the f -mode frequency σ 0 of the corresponding non-rotating star. Figure 1 displays this behavior for more than 230 different neutron star models of each the co-rotating (i.e., stable) and counter-rotating (i.e., potentially unstable) branches of the f -mode for seven EoSs and various central energy densities (with corresponding central rest mass densities ρ c ∈ [2.2, 7.3]ρ 0 , where ρ 0 = 2.7 × 10 14 g/cm 3 is the nuclear saturation density); we model the universal behavior using the quadratic function The results of a least squares fit are a u 1 = −0.193 and a u 2 = −0.0294 for the potentially unstable branch and a s 1 = 0.220 and a s 2 = −0.0170 for the stable branch of the f -mode. The quadratic fit accounts well for the increasing oblateness of the star with its rotation; however, close to the Kepler limit, deviations from this simple model become visible. As this deviation is most pronounced for the less realistic polytropic EoSs, we do not take them into account for the quadratic fits. The root mean square of the residuals is 0.024 for the counter-rotating branch and 0.048 for the co-rotating branch.
We point out that our model predicts that the unstable branch of the quadrupolar f -mode becomes susceptible to the CFS instability once the angular rotation rate of the star exceeds Ω ≈ (3.4 ± 0.1) σ 0 (when considering sequences of constant central energy density); note that the given uncertainty is a bound, not a confidence interval. This finding regarding the critical value complements the well-known threshold of T /|W | ≈ 0.08 ± 0.01 in terms of the ratio of rotational to gravitational potential energy ( [96,97]), which is confirmed in our simulations and is in contrast to the widely used Newtonian result of T /|W | ≈ 0.14.
The stable branch of the f -mode can be fitted more accurately when switching to the comoving frame and considering sequences of constant baryon mass. The frequency σ c observed in the comoving frame is related to the frequency observed in the inertial frame via σ c = σ i + mΩ/2π. We show our results for more than 120 different neutron star models using four realistic EoSs in Figure 2. We fit our results to the quadratic function note that we use the Kepler velocity Ω K to normalize the star's rotation rate in this formula. gravitational mass and moment of inertia, inspired by [44]. We will be guided by the model employed in the Cowling approximation by [100] which reproduces the f -mode frequency of a 1 Similar relations were presented in [98,99] but in the Cowling approximation. particular neutron star from its rotation rate, gravitational mass, and effective compactness.
We propose the fitting formulâ whereσ i :=M σ i / kHz andΩ :=M Ω/ kHz; note that we set d 2 = 0 as it turns out that this coefficient would be afflicted with a large uncertainty. Using around 100 models based on polytropic as well as around 400 models based on realistic EoSs, the resulting coefficients from a least-squares fit for the counter-rotating branch of the f -mode are of magnitude with those published by [100] in the Cowling approximation; comparing the special case of no rotation,Ω = 0, our fitting formula yields roughly 20 % lower frequencies in our fully general relativistic setup, which is in accordance with expectations.
The lines of constantΩ in Figure 3 may give the impression that the CFS instability operates more easily in stars with low (effective) compactness, seemingly in contrast to the finding that post-Newtonian effects tend to enhance this instability [101]. This paradox can be resolved by noting that relativistic effects mainly shift the f -mode frequency to lower values while the inclination of the lines of constantΩ is largely unaltered (cf. Figure 3 in only the mass of the postmerger object will be needed to constrain its parameters by using only the asteroseismological relations (see [98,100,102]). This will be an independent yet complementary constraint in the estimation of the radius in addition to those based on the Love numbers (see [103][104][105][106]).
As a graphical illustration of the behavior of the f -mode frequency across the entire parameter space of stable equilibrium models, we present in Figure 4 the frequency of the counter-rotating branch as obtained from the time evolutions exemplary for EoS H4; the graph will be qualitatively similar for other EoSs. We constructed several hundred neutron star models across the entire M -R e plane along so-called evolutionary sequences; those are sequences of differently fast rotating neutron stars that share the same baryon mass. The two-dimensional plane of equilibrium models has four distinct boundaries: the static limit at the lower left along which the non-rotating models are located; second, the mass-shedding limit bounds the equatorial radius to the right; next, in the top left of the graph, the limit of stability with respect to quasi-radial perturbations connects the (non-rotating) maximum The graph shows that for this particular EoS each neutron star model may become CFSunstable if it is sufficiently spun up. For heavier neutron stars, this happens considerably below the mass-shedding limit and sufficiently heavy supramassive models (approximately M 2.25M ) will inevitably be CFS-unstable (with respect to the quadrupolar f -mode); i.e. those stars will be stabilized merely by viscous mechanisms counteracting the CFS-instability.

B. Universal Relations involving long-lived remnants from BNS mergers
Inspired by the works on universal relations for single neutron stars ( [13,15,110]), the last five years have given rise to universal relations for BNS: they relate the pre-merger neutron stars to the early post-merger remnant, and have been developed using numerical relativity simulations ( [51,111,112]).
These works have primarily focused on relating the tidal deformability of the pre-merger stars, which impact the dynamics of the pre-merger gravitational waves at leading order through the reduced tidal deformability or binary tidal deformabilityΛ ( [113,114]), to various stellar parameters of the early remnant.
Recently, [115] investigated empirical relations for BNS mergers based on the extensive CoRe data set of numerical relativity gravitational wave simulations ( [116]). Covering a wide range of mass ratios, they find an extensive set of universal relations involving the various peak frequencies of the post-merger gravitational wave signal, involving, e.g., the chirp mass and characteristic radius of a 1.6 M neutron star. In particular, they also find universal relations between the binary tidal deformability and the primary f -mode frequency of the post-merger signal (as in [112]), however, this time involving the chirp mass of the BNS.
Other established universal relations connect the post-merger peak frequency to the tidal coupling constant ( [41,51]).
The universal relation between the binary tidal deformability of the BNS and the stable, co-rotating f -mode frequency σ s of the early, differentially rotating remnant proposed in [112], led us in [55] to derive a similar relation for a potentially long-lived, uniformly rotating remnant: the relation takes the form whereσ s = M M σ s kHz is the normalized co-rotating f -mode frequency, and q = M 1 M 2 ≤ 1 the gravitational mass ratio of the pre-merger stars. For rapidly rotating, long-lived remnants (with rotation frequencyΩ ≥ 800 Hz), this relation achieves an average relative error of 1.3 %.
We also derive a relation for the potentially unstable, counter-rotating f -mode frequency of the long-lived remnant, presenting the possibility of predicting the onset of the earlier mentioned CFS-instability.
Combining these results with the universal relation, Eq. (11), for fast rotating neutron stars between the stable, co-rotating f -mode frequency and the effective compactness η = M 3 /I 45 (as defined in the previous section), we also derived a combined relation of the form that relates the pre-merger binary tidal deformability of the BNS with the effective compactness of the long-lived remnant. For rapidly rotating remnants, this relation achieves an average relative error of 2.4 %.
Finally, by directly relating these quantities without going via the f -mode, we obtain a universal relation of the form This relation achieves improved accuracy, reaching an average relative error of ∼ 1.5 % for remnants with any rotation frequency. We show the quadratic fit in Fig. 5 for the symmetric case of q = 1.
We finally also consider a direct relation between the binary tidal deformability and the compactness C = M/R of the long-lived remnant. Such a relation would allow the direct estimation of the remnant's radius R using independent estimates of its gravitational mass.
We propose a relation of the formM which, however, only achieves an accuracy an order of magnitude worse than for the effective compactness relation, reaching an average relative error of ∼ 8.8 %. A graphic representation of this fit along with the data is shown in Figure (6).
The results presented in [55] represent a first step towards finding universal relations between the pre-merger neutron stars and the potential long-lived remnant of a BNS merger  using perturbative calculations. Our approach can be freely extended to, e.g., hot EoSs, phase transitions, as well as differential rotation for the remnant to, to cover, e.g., earlier parts of the post-merger phase.

IV. APPROACHING THE INVERSE PROBLEM
Having laid out the procedure to extract fluid oscillation frequencies from time evolution in a general relativistic framework and having developed various universal relations between bulk properties of neutron stars and their f -mode frequencies, we now turn to one of the possible applications. The field of gravitational wave asteroseismology attempts to invert the above outlined procedure and aims to constrain bulk properties of neutron stars such as mass and radius given one or more oscillations frequencies. The above derived universal relations are a valuable tool for this as they do not depend on the underlying and hitherto poorly understood nuclear equation of state. However, even without universal relations at hand, the inverse problem can be approached.
We will summarise an approach of the inverse problem using Bayesian methods as laid out by [89]. Here, the fundamental idea is tested on the f -mode frequency, but it may easily be extended to other fluid modes or combined with their damping times. The use of the Bayesian framework is advantageous (compared to, e.g., an analytical inversion of the non-linear universal relation) as it directly allows to incorporate error bars of the measurement of the f -mode into the calculation. Besides the Bayesian approach and universal relations, there are in principle also semi-classical techniques, e.g., WKB theory, which can be used to address the inverse spectrum problem in simplified cases; see [117]¸for using axial w-modes of spherically symmetric neutron stars and [118] using a similar approach for ultra compact stars.
In a first step, we pick an EoS and a particular (rotating) neutron star model for which we then calculate two f -mode frequencies. We assume those to have a relative error of 3 % This analysis required the assumption of a particular EoS and so far we used the same EoS to reconstruct the model that we also used to generate the underlying neutron star model.
If we use a different EoS to perform the Bayesian analysis, we will obviously reconstruct a different model; this is shown in orange in Figure 7. The reconstructed star in this case is considerably smaller and lighter. Note that the apparent cut of the MPA1 EoS posteriors (in particular for the mass) is not an indication that this EoS has not been used for the injection, since a similar behavior can also be found for the H4 EoS, if the injected parameters are closer to the edge of the H4 EoS neutron star parameter space. One needs further information, such as the precise mass, the tidal deformability, or perhaps a damping time, in order to rule out this EoS.
Instead of making an educated guess for the "correct" EoS and reconstructing mass and radius as just described, we may also use the universal relation proposed in Eq. (12). We will not need to assume any EoS but the universal relation will provide us-given two f -mode frequencies and a good estimate of the star's mass-with an estimate for the stars's rotation rate Ω and its effective compactness η.
We again use the Bayesian method to invert the universal relation and we apply it to the same model as in the previous analysis; the H4 model with M = 1.8 M and R = 15 km.
In Figure 8 we show the results of the analysis. On the x-axis, we show the star's effective compactness and rotation rate, both normalised to the actual value corresponding to the equilibrium configuration. The two graphs then show the posterior distributions. We ran the same analysis twice, once assuming a 30 % relative error on the prior mass (solid lines) and once assuming a 10 % relative error (dashed lines). As expected, with a less informed prior on the mass, the posterior distribution for the desired quantities have a considerably larger variance.
The blue and orange curves correspond to the previously describe method where we assume a particular EoS. We have described this above and here we show the posteriors for η and Ω. The green curves show the posterior distribution employing the universal relation.
It is evident that the H4 EoS method (blue lines) and universal relation (UR) method (green lines) yield very similar results for the rotation rate Ω, while assuming the MPA1 EoS (orange lines) indicates a value that is larger than the correct one. Note that both observations hold independent of the specific prior knowledge of M assumed here (30 % and 10 %).
The situation for the effective compactness η is qualitatively different. First, the prior knowledge of M plays a big role for the UR method, but is less important for the EoS methods. For those we find that the correctly assumed H4 EoS is almost independent of uncertainties in M , while the posterior distribution obtained by the MPA1 EoS is shifted.
Note that the rather different scaling behavior of the UR method is in agreement with the findings described above.
Finally, while the posteriors of Ω are very smooth, one observes small "bumps" for the H4 EoS, e.g., at η/η 0 ≈ 1.02. We have verified that this likely is an artifact from the finite resolution and particular range of the used H4 f -mode data that is available to us. This directly sets the scale of how precise our currently implemented EoS data can be used to resolve the underlying parameters, which is of order percent level.

V. SUMMARY AND OUTLOOK
We report the first extraction of frequencies of the l = |m| = 2 f -mode of general relativistic, rapidly rotating neutron stars without the commonly used slow-rotation or Cowling approximation to an extent that allows us to generalize the findings into universal relations. This concludes a long-standing open problem, building upon the effort from numerous studies throughout the past five decades.
For this, we have derived a set of time evolution equations governing the perturbations of the fluid of a rapidly rotating neutron star as well as its surrounding spacetime, derived in a perturbative framework in full general relativity. We have opted for the Hilbert gauge in order to arrive at a set of fully hyperbolic equations for the spacetime perturbations whose implementation does not pose major obstacles.
Convergence tests to study the accuracy of our code reveal that in a typical simulation, the obtained frequency usually deviates less than 1 % − 2 % from the limiting value or has an accuracy of about 10 Hz. As we model axisymmetric configurations only, we can reduce the problem to two spatial dimensions which drastically lowers the computational expense of our numerical time evolutions in comparison to non-linear codes that perform three-dimensional simulations. The evolution of the perturbations of a neutron star for 15 ms on a grid with 3120 × 50 points, which facilitates a decent frequency resolution, requires only a few dozen of CPU hours; this enables us to study broad ranges of parameters and various EoSs. We expect to further reduce the computational cost even further by deriving a simplified set of perturbation equations.
As a result, we provide different universal relations for the frequencies of the l = |m| = 2 f -modes of uniformly rotating neutron stars at zero temperature which are independent of the EoS; the proposed formulae are calibrated to several hundred neutron star models that are constructed using both polytropic and realistic EoSs and are scattered across the entire parameter space of equilibrium solutions. Furthermore, it is possible to link the pre-merger binary tidal deformability to the effective compactness of the late post-merger remnant in a universal manner, too. Such universal relations will be an essential piece in the asteroseismological toolkit once the third-generation GW observatories will be able to pick up the ring-down and fluid ringing signal following the merger of a binary neutron star system; they allow to solve the inverse problem, leading to significantly tighter constraints for mass and radius of the postmerger object. For this task, it is elementary to have a smorgasbord of universal relations at hand, which allows to make a practical choice, depending on which observables are available or in which of the star's properties one is interested. We will extend the present list of such universal relations in future articles, utilizing different combinations of bulk properties of the star; while we may obviously be (and already have been) inspired by previously published fitting formulae that were derived using different approximative frameworks, we need to be open-minded about models involving novel combinations of observables.
We also report the discovery of an accurate estimate for the onset of the CFS-instability when the f -mode frequency of the non-rotating member of the family is known and verified the corresponding critical value of T /|W |.
A natural extension of the present study will be a more comprehensive investigation of the spectrum of neutron stars (i.e. higher multipole f -modes, low p-modes and g-modes as well as w-modes) which may be excited in different astrophysical processes. Furthermore, we are going to extend our code to account for differentially rotating neutron stars and hot EoSs that are particularly relevant for nascent neutron stars or postmerger configurations in the immediate aftermath of a binary merger, both of which will have a considerable impact on the vibration frequencies or the onset of the CFS-instability (and via two further scaling parameters also on the universal relations) during a very short but dynamic interval of their lives.
Furthermore, a precise knowledge of the spectrum of compact objects is invaluable as also isolated neutron stars as well as those in inspiraling binary systems may possess high spin rates and various oscillation modes, which may be excited, e.g., via glitches or tidal coupling, may impact their electromagnetic emission or the dynamic of the whole system.