{'fraction_non_alphanumeric': 0.028581598149136855, 'fraction_numerical': 0.01185264281900694, 'mean_word_length': 4.770873985827257, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 2, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Dynamics of Fast Rotating Neutron Stars: An Approach in the Hilbert Gauge

We describe a set of time evolution equations and its numerical implementation 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. We derive the perturbation equations for the spacetime in the Hilbert gauge, while the hydrodynamical evolution is based on perturbations of the energy-momentum tensor. In our numerical implementation, we use Kreiss-Oliger dissipation in order to achieve a stable time evolution. Our code features high accuracy at comparably low computational expense and we are able to extract the frequencies of non-axisymmetric modes of compact objects with rotation rates up to the Kepler limit.


I. INTRODUCTION
The systematic study of non-axisymmetric neutron star oscillations began in the 1960s with the pioneering work of Thorne and collaborators [1,2] (and subsequent papers), in which they laid out the equations governing the perturbations of a perfect fluid. The numerical solution of these equations has proven highly challenging and it has taken nearly two decades before Lindblom & Detweiler found an advantageous formulation of the problem as an eigenvalue problem in 1983 that allowed them to numerically integrate the perturbation equations and determine the real and imaginary parts (the frequency and the damping time) of the acoustic modes of sufficiently realistic stellar models [3,4]. These results did not conclude the investigation of these modes; in particular Chandrasekhar & Ferrari turned to the perturbations of perfect fluids and shed further light on the structure of the equations and their solutions [5,6] (and subsequent papers).
At a similar time in the very early 1990s, after investigating a toy model [7], Kokkotas & Schutz applied techniques borrowed from studies of black hole quasi-normal modes (in their case the WKB approximation which they combined with numerical methods) to the solution in the exterior of the star, and proved that the dynamic spacetime of a neutron star exhibits its very own class of modes and christened them w-modes [8]. The large imaginary parts of the eigenvalues of these modes pose numerical challenges and this strategy of exploiting black hole science was successfully followed further. In particular, the application of the continued fraction method, the Wronskian method or integrating along anti-Stokes lines led to a considerably more accurate determination of w-modes (as well as fluid modes) of neutron stars and the discovery of the w II -modes [9,10]. These numerical solutions finally allowed for the first time to establish universal relations involving frequencies and damping times of os- * christian.krueger@tat.uni-tuebingen.de † kostas.kokkotas@uni-tuebingen.de cillation modes, enabling the solution of the inverse problem within the so-called asteroseismology, in which one constrains mass and/or radius of the neutron star, and eventually the nuclear equation of state, via the observation of oscillations [11][12][13][14].
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 Allen et al. in 1998 [15]. In contrast to the integration of ordinary differential equations, the evolution of hyperbolic equations in time is very sensitive to potential irregularities in the coefficients (which in realistically modelled neutron stars are, due to physical phase transitions in the matter, very difficult to avoid) and quickly develops numerical instabilities; Ruoff reformulated those equations by means of the ADM-formalism [16] and, crucially, introduced a non-uniform radial grid to overcome such instabilities when using realistic equations of state [17,18].
All the previously mentioned studies were concerned with non-rotating stars. 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 to consider even rotating stars as spherically symmetric [19]. In this so-called slow-rotation approximation, the perturbation equations gain considerably in complexity and have been written down in the Regge-Wheeler gauge [20] first by Kojima in 1992 [21]. Even though the problem remains one-dimensional, its solution is not straightforward as (amongst other technicalities) the outgoing-wave boundary condition at infinity is elusive. Notwithstanding, Andersson successfully applied this formalism and discovered that r-modes are prone to the so-called CFSinstability (named after their discoverers Chandrasekhar, Friedman, and Schutz) [22][23][24] at any rotation rate [25].
There has been continuing effort using the slow-rotation approximation concerning rotational modes [26,27], and also employing different gauges [28][29][30], 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 [31], 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), Yoshida & Eriguchi computed quadrupolar f -mode frequencies of rapidly rotating neutron stars and studied the associated CFS-instability in the late 90s [32,33]. Boutloukos & Nollert revived this approach and investigated the general properties of the spectrum of neutron stars regarding the acoustic and Coriolis-driven modes [34]. As a further step toward a more general relativistic treatment, Yoshida revisited the problem in the conformal flatness approximation [35]. However, with the mathematical difficulties of extending the eigenvalue formulation to include a dynamic spacetime, the focus shifted to solving this problem as a time evolution.
Despite their complexity, the perturbation equations for rapidly rotating relativistic stars have been written down by Priou already in 1992 [36], even though they were not approached numerically at that time. The most straightforward way of dodging the challenges arising from General Relativity is avoiding this theory entirely and treating the task in Newtonian theory. While this is obviously a crude approximation to the strong gravitational field of a neutron star, the numerical solution is less involved and may still result in qualitatively valuable results. Jones et al. completed this task formulated as a time evolution in 2002 and their code was later extended [37,38]. Returning to a general relativistic treatment, Gaertig & Kokkotas worked in the Cowling approximation and successfully extracted f -mode frequencies of arbitrarily fast rotating neutron stars by adding artificial viscosity (also known as Kreiss-Oliger dissipation) [39] to their evolution equations in order to stabilise their time evolutions [40,41].
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 obviously 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 threedimensional 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 [42,43] and those of (quasi-)radial modes in full General Relativity [44] had been reported. The non-linear codes kept evolving and were used to generate mode frequencies of f -modes in the conformal flatness approximation [45] or those of inertial modes in the Cowling approximation [46]. Not much later, the frequencies of non-axisymmetric modes in full General Relativity of rapidly rotating neutron stars obtained from fully non-linear simulations have finally been reported by Zink et al. in 2010 [47]. 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 this paper, we will present a perturbative approach to the computation of mode frequencies of rapidly rotating neutron stars accounting, for the first time, for a dynamic spacetime. We study perturbations on an axisymmetric background which allows us to reduce the dimensionality of the problem and keep the computational cost low. For this effort, we extend a code previously developed by us, which we used to investigate oscillations of differentially neutron stars in the Cowling approximation [48]. Except for the inclusion of a few spacetime terms, the hydrodynamical part of the problem remains largely unchanged; for the formulation of the evolution equations of the spacetime perturbations we choose the Hilbert gauge as this gauge will immediately lead to a numerically attractive, fully hyperbolic system. After performing several intrinsic code tests such as convergence tests concerning the grid resolution and the amount of artificial viscosity applied as well as examining the violation of the Hilbert gauge, we compute mode frequencies and compare them to published values from non-linear simulations. As this paper focuses on the technical details of the problem, we report the astrophysically relevant result in an accompanying paper [49].
Unless otherwise noted, we employ units in which c = G = M ⊙ = 1 throughout this paper.

II. MATHEMATICAL FORMULATION
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 energymomentum 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 isotropic coordinates takes the form 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 equation of state (henceforth 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 [50][51][52].

A. Introducing Perturbations
We will now introduce time-dependent perturbations 1 for which we will derive evolution equations. First, we decompose the metric as 1 Note that we assume a harmonic azimuthal dependence for all perturbations, i.e.,Q(t, r, θ, ϕ) = Q(t, r, θ)e imϕ for any perturbation variable Q. This results in the reduction of the dimensionality of the problem and a substantial decrease of computational expense.
where g (0) µν is the background metric and h µν its perturbation; we use the background metric to raise and lower the indices of the latter. 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 Regge-Wheeler gauge [20], and hence to reduce the number of perturbation equations. In this study, however, we take 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 electro-magnetism, 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 wave-like 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 [15,18,28,53] or the ones without any gauge choice [36] where the field equations split into subsets of hyperbolic and elliptic equations which have to be solved simultaneously. The fully hyperbolic character of the perturbation equations in the Hilbert gauge make this gauge particularly convenient for the numerical implementation in a time evolution.
Having opted for the Hilbert gauge, we are dealing with ten functions describing the spacetime perturbations. In choosing our perturbed line element, we are inspired by the one suggested by Priou [36], which has also been used by Stergioulas & Friedman in the study of zero-frequency f -modes in rapidly rotating stars [54]; we apply it to the trace-reversed metric perturbations and introduce a few small modifications (mostly factors of r and sin θ) which turn out crucial in avoiding numerical instabilities. We use calligraphic letters to denote the ten spacetime perturbations A, B, H, K, L, M, P, Q, W, and Y and write the perturbed metric such that the individual com-ponents take the form φ tt = −2e 2ν H + 2e 2ψ ωr sin θ e 2ν Y + ωr sin θW , (10) φ tr = L + e 2ψ ωr sin θA, (11) φ tθ = rM + e 2ψ ωr 2 sin θB, (12) φ tϕ = −e 2ψ r sin θ e 2ν Y + 2ωr sin θW , With this particular choice for the metric components, we can spell out the four constraint equations emerging from the Hilbert gauge, Eq. (8); we show them in Appendix B.

B. Fluid Perturbations
Our description for the fluid perturbations is inspired by the very convenient formulation that has been developed and successfully used in [48,[55][56][57] for studies within the Cowling approximation. The system of equations describing the hydrodynamics is essentially the same as in the cited literature, except that the equations gain numerous terms accounting for the spacetime perturbations. Nonetheless, we will give a brief overview of their derivation here.
The fundamental idea of the approach is to use the perturbations of the energy-momentum tensor as evolution variables, rather than the typical fluid primitives such as pressure or fluid velocity; this results in a significant shortening of the fluid perturbation equations. In order to avoid confusion regarding the definitions of the perturbation variables in comparison to previous studies, we retain the shortcuts defined in the Cowling approximation and write the perturbed energy-momentum tensor as where δT µν C is the perturbation of the energy-momentum tensor obtained in the Cowling approximation and the minus sign in front of the last term is to account for the fact that h µν = −g µα g νβ h αβ . We will use the six shortcuts Note that we have, in fact, slightly modified the definitions of Q 2 , Q 4 , and Q 5 when compared to previous studies using this formalism. By doing so, we remove the degeneracy of the corresponding base vectors at the rotation axis and considerably improve the numerical stability of the time evolution; technically, these modifications merely shift about a few factors of r and sin θ in the perturbation equations but do not substantially alter them, which is why we decide to accept the minor deviation in definition from the previous studies without introducing differently named variables. The fluid primitives, i.e. the perturbations of energy density, pressure and the fluid velocities, can be reconstructed by inverting the definitions of the variables Q 1 , Q 2 , Q 3 , and Q 4 . Further, it can be shown that both Q 5 and Q 6 are not independent variables but are rather a linear combination of the other perturbation variables, We show the full relations for those two variables in Appendix C.

C. Perturbation Equations
We specified 14 perturbation variables that we need to evolve in time: 10 for the spacetime and 4 for the fluid. The evolution equations follow in a very straightforward manner from the perturbed Einstein equations and the perturbed law for the conservation of energymomentum We show the former set of equations, governing the dynamics of the spacetime, in Appendix D and the latter set, governing the dynamics of the fluid, in Appendix E. One special case of the perturbation equations, which we will consider for testing purposes later in more detail, cf. Section IV, is the non-rotating limit. Setting the star's rotation rate to zero, Ω = 0, immediately implies one of the metric potentials to disappear, ω = 0, and that the polar derivative of the all background quantities vanishes, ∂ θ = 0. Furthermore, in non-rotating stars, the oscillation modes of different azimuthal index m are degenerate and we may set m = 0 without loss of generality. This results not only in a dramatic shortening of the evolution equations, but the set of the 14 evolution equations also decomposes into two fully decoupled sets of equations: one set for polar perturbations described by the variables H, K, L, M, P, Q, W, Q 1 , Q 3 , and Q 4 ,and one set for axial perturbations described by the variables A, B, Y, and Q 2 . Beside being an obvious test case for our code, the non-rotating limit has been studied to great detail using time evolution codes as well as eigenvalue codes (see, e.g., [8,15,18,[58][59][60]) and will not offer any new scientific insight to us.

D. Boundary Conditions
The perturbation equations are complemented by boundary conditions which describe the behaviour of the perturbations on the boundaries of the numerical domain. As we will explain in Section III, our numerical domain is bounded by the star's rotation axis and the equatorial plane in the polar direction; in the radial direction, we have to apply boundary conditions at the star's surface for the fluid perturbations as well as at the outer edge of the numerical grid for the spacetime perturbations. As our description of the fluid is largely the same as in [48], most of the boundary conditions for the fluid can be found there but due to slight modifications in the definitions and for the sake of completeness, we will repeat them here.
• Having the equatorial plane, characterised by θ = π/2, as a boundary of our numerical domain allows us to select one of two disjoint sets of oscillations modes that behave differently under reflection about the equatorial plane. Let R be the corresponding reflection operator, which in spherical coordinates is given by the mapping (r, θ, ϕ) → (r, π − θ, ϕ); it is easy to show that for any spherical harmonic, Y lm , the relation RY lm = (−) l−m Y lm holds. In the non-rotating limit, all perturbation variables can be expanded into spherical, vector and tensor harmonics [20] and, according to their behaviour under the operator R, they can be categorised into the two sets R 1 := {A, H, K, L, P, W, Y, Q 1 , Q 2 , Q 3 } and R 2 := {B, M, Q, Q 4 }. Depending on whether l − m is even or odd, 2 the perturbation variables of one of those two sets are either zero at the equatorial plane or their θ-derivative vanishes there; we show a summary in Table I.
• To derive the boundary conditions at the rotation axis, which is given by θ = 0, we will again utilise the angular behaviour of the perturbation variables as specified by the spherical harmonics and their tensorial equivalents. In this instance, we have to evaluate the limits of the perturbations as θ approaches 0 for each azimuthal index m; however, we will restrict our study to the astrophysically relevant cases of m ≥ 2. We find that for m ≥ 3, all perturbations vanish along the rotation axis, while for m = 2 they separate into the set of perturbations Θ 1 := {B, P, W}, whose θ-derivative vanishes at the rotation axis, and the set of perturbations Θ 2 := {A, H, K, L, M, Q, Y, Q 1 , Q 2 , Q 3 , Q 4 }, which become zero along the rotation axis. We show a summary of the boundary conditions at the rotation axis in Table II.
• The origin is mapped onto a boundary line of our numerical grid. Physically, however, it remains part of the rotation axis and the equatorial plane. When a perturbation variable is shown to vanish on either of the neighbouring boundary lines, we set it to zero at the origin as well; otherwise, we use the same value as on the next grid point.
• At the stellar surface, we need to apply boundary conditions to the fluid perturbations. The variables Q 3 and Q 4 vanish by virtue of their definition because pressure and energy density vanish there. Q 1 and Q 2 merely need to exhibit "numerically smooth" behaviour there and we extrapolate them linearly in radial direction from the interior of the star. For the case of rotating and hence oblate stars, the stellar surface traverses our grid. As an approximation to the true stellar surface, we apply the boundary conditions for the fluid variables at the very first grid point (counting from the origin) that is located outside the star; see also Figure 1.
• Last, we consider the outer edge of the numerical grid which is typically located about a hundred or so star radii away from the origin and can be chosen arbitrarily; there is no physical boundary or mathematical singularity at this location. We employ a Sommerfeld boundary condition at this edge in order to remove the energy carried by gravitational waves from the grid. While it accounts correctly only for spherical wave fronts but not for waves traveling in arbitrary directions, this condition offers a highly attractive ratio of simplicity of implementation when compared to the amount of artificially reflected energy and has also been used in 3-dimensional simulations [61,62]; furthermore, at the distances mentioned above, a spherical wave front will be a very good approximation to the outgoing radiation for our relevant scenarios. This is confirmed in our simulations, in which we observe only marginal reflection that does not spoil the later stages of our time evolutions. Table I. Boundary condition of the perturbation variables at the equatorial plane, θ = π/2. The two sets R1 and R2 are defined in the text and the stated boundary condition applies to each member of those sets.
We need to comment further on the boundary conditions in polar direction, i.e., those applied at the equatorial plane and at the rotation axis; we have derived those from the angular behaviour of the spherical harmonics Table II. Boundary conditions of the perturbation variables along the rotation axis, θ = 0. The two sets Θ1 and Θ2 are defined in the text and the stated boundary condition applies to each member of those sets.

Θ1
Θ2 which requires the background configuration to be spherically symmetric. While it has been shown that those boundary conditions are valid also in the slow-rotation approximation [21,63], it is not immediately clear that they also hold in rapidly rotating stars. The equatorial plane remains a plane of symmetry also when the star is oblate. As none of the two hemispheres is distinguished, the oscillation modes remain eigenfunctions of our previously defined reflection operator R. Further, it has been argued within the slow-rotation formalism that all oscillation modes, even though gaining contributions from other modes, keep their overall parity [26,64]; hence, their eigenvalue with respect to R remains unchanged. Further, despite the simplification of the slow-rotation approximation, the relevant couplings between the different mode classes is already present and it remains the same in rapidly rotating stars; thus, the argument can be extended to the case of rapid rotation without further ado, leaving the boundary conditions unaltered.
In order to understand the boundary conditions at the rotation axis without relying on the star's spherical symmetry, we study the behaviour of the perturbation variables under rotation about the rotation axis, (r, θ, ϕ) → (r, θ, ϕ ′ ). All of our evolution variables describe scalar, vectorial or tensorial (of rank 2) perturbations; in order for them to be well-defined at the rotation axis, they may have a non-zero value only if their azimuthal dependence resembles that which is specified by the azimuthal parameter m. Hence, for m ≥ 3 all perturbations have to vanish along the rotation axis; when m = 2, all scalar and vectorial perturbations have to vanish there, too, and only the tensorial components (that have rank 2) of the perturbed metric, i.e. B, P, and W, may be non-zero. Instead, it is their polar derivative which has to vanish there, i.e. ∂ θ = 0, to ensure smoothness of the solution across the poles.

III. NUMERICAL METHODS
After decomposing the perturbed quantities with respect to ϕ, the resulting evolution equations form a twodimensional problem in the spherical coordinates r and θ. We evolve the perturbations on a grid with radial coordinate s ∈ [0, 1] and angular coordinate t ∈ [0, 1]; the mapping to the spherical coordinates r and θ is given by where r e is the equatorial coordinate radius of the star and a and b are two parameters with which we can adjust the radial grid. This particular radial grid function, r(s), has a varying grid spacing across the entire radial domain and allows us to avoid the introduction of different grids for different parts of the radial domain and grid interpolation between those different grids. In particular, the grid spacing is comparably coarse near the centre of the star, becomes much finer in the region around the star's surface, and then becomes coarser again in the exterior of the star where it approaches a near uniform grid spacing at the outer edge; we show a sketch of our computational grid in the vicinity of a fast rotating star in Figure 1. Furthermore, it comes with two parameters a and b with which we can steer the behaviour of the grid; in particular, we choose a and b such that the outer edge of the grid is roughly a hundred star radii away from the star and that the physical grid spacing in radial direction is limited by half the star's radius in that region. We find that this grid offers a reasonable compromise between having a grid extending sufficiently into the far field, while keeping the number of grid points and hence the computational expense low.
The grid coordinates are not fitted to the shape of the star, which means that the surface of rotating and thus oblate stars lies in between the grid points of our computational domain. We approximate the surface of the star to be located at those grid points for which the energy density first vanishes when moving outward from the centre along lines of constant polar angle; this is where we apply the surface boundary condition for the fluid variables.
We use central differences of second order to discretise the spatial derivatives. For the time evolution, we use the method of lines and apply the third-order Runge-Kutta method as time integrator. Similar to previous studies within the same perturbative framework [41,48], this code is prone to numerical instabilities. As the hydrodynamical equations are essentially identical to those in [48], we also encounter the same instability near the origin of the star; we apply second-order Kreiss-Oliger dissipation [39] to cure those instabilities. The spacetime perturbations do not suffer from similar instabilities near the origin as potential instabilities emerging from truncation errors will be radiated away; however, we observe some unstable behaviour in the spacetime perturbations in the far field where the grid becomes coarser. These instabilities stem from the fact that our grid has difficulties resolving the gravitational waves in the wave zone. We could increase the grid resolution to overcome this issue, however, we find that applying a very small amount (typically 10 to 100 times smaller than that for the fluid) of Kreiss-Oliger dissipation to the spacetime grid eliminates those instabilities, too. Furthermore, even an increased grid resolution cannot entirely prevent this issue since gravitational waves of sufficiently high frequency will confront us with the same difficulties.
The typical grid resolution that we employ for our simulations is 3120 × 50 points, where the radial grid extends to 100R e . The interior of the star is then covered by 181 × 50 grid points. For this resolution, we vary the dissipation coefficients for the fluid between ≈ 1 × 10 −7 − 5 × 10 −6 , even though we distinguish between dissipation coefficients for the radial and for the angular direction (in some cases a larger dissipation in angular direction is required); from our experience, the amount of dissipation required to achieve a stable numerical simulation depends on the background model and varies with central energy density and rotation rate. The spacetime perturbations behave in a more benign manner and we leave the dissipation coefficient of 10 −8 largely untouched. In general, finding those dissipation coefficients that are "optimal" in the sense that the least amount of artificial viscosity is closest to the physical reality (ignoring the fact that we are dealing with highly idealised neutron star models) is within certain bounds a matter of trial and error.
In order to extract the mode frequencies, we observe the temporal evolution of one of our perturbations variables at one freely chosen point on our grid and then cal-culate the discrete Fourier Transform of this time series. A considerable improvement of the obtained frequencies can be achieved by locating the roots of the derivative of the power spectral density, which we approximate with second order central finite differences.

IV. CODE TESTS
We have described the mathematical formulation of the problem as well as the numerical methods that we have employed to solve the problem. It is crucial for a newly developed code to be thoroughly tested. We will start by showing results from several convergence tests for various modes of non-rotating and rotating configurations; along with those tests, we will list for future reference the f -mode frequencies of the well known BU sequence which has served as a popular test case in various studies [41,45,46,65]. Next, we will turn to our treatment of the Hilbert gauge and how its violation evolves throughout our simulations. Last, we will reproduce some mode frequencies and compare them to previously published data.
A. Convergence Test

Non-Rotating Neutron Star Model
Our very first convergence test is for the f -mode of the regularly employed BU0 model, which is a simple polytropic model with no spin, characterised by the polytropic constants N = 1, K = 100 and a central energy density of ǫ c = 0.891 · 10 15 g/cm 3 ; these parameters result in a neutron star model with typical values of M = 1.400M ⊙ and R e = 14.16 km. Even though this model has previously been studied extensively, we performed a detailed convergence test in order to understand the general behaviour of our code, not impacted by potential numerical challenges stemming from the oblateness of rapidly rotating models or tabulated EoSs.
We perform a series of simulations using several different grid resolutions; starting from our "standard" grid resolution of 3120 × 50 points, we decrease or increase the number of grid points in uniform steps. Since we have chosen a highly non-uniform radial grid (cf. Section III), we show some of its characteristics in Table III; it is apparent that the grid spacing or the number of points inside the star is not proportional to the number of grid points in radial direction, but this does not pose a problem as long as the map between them is strictly monotonic. We emphasise that we show physical distances relating to the evolution grid of the BU0 model in Table III for intuition only, however, they are not universal and depend on the particular neutron star model, in particular its oblateness; they scale with the star's equatorial radius, R e , and hence vary only moderately. Table III. We show some characteristic values of our radial grid for different grid resolutions. It is in the nature of our grid function that the part of the grid that covers the star depends on the overall number of grid points (nstar is the number of grid points inside the star along the equatorial plane). In particular, the radial grid spacing at the origin, ∆r(0), is comparably large, while the radial spacing at the surface, ∆r(Re), is very fine, which allows a good resolution of the star's surface. For all listed grids, the outer edge of the grid is located at 100Re and the radial spacing at this edge is 5950 m. The values shown here are indicative and depend mildly on the oblateness of the star. When comparing obtained frequencies for different grid resolutions, we have to take into account the necessity for artificial viscosity which will impact the physical observables to a small extent. For the BU0 model, we find that the lower bound for the dissipation coefficients is somewhere around 10 −7 ; we run the simulations for different grid resolutions as well as for three different fluid dissipation coefficients and show the resulting frequencies of the f -mode and the p 1 -mode for l = 2 in Tables IV and  V; further, we show the frequencies of both modes and linear fits in Figure 2. We want to make a note on the accuracy of the listed frequencies: In addition to locating the roots of the derivative of the power spectral density, we also analyse smaller chunks of the time series. By varying size and position of the time series window, we obtain marginally different values (due to spectral leakage) over which we then average; this procedure results in a precision of well below 1 Hz. For this method to work, it is crucial that we ensure the surface of the star to match up with an angular grid line (the location at which we apply the surface boundary conditions); this requires a meticulous fine-tuning of the grid parameters and is in this way, unfortunately, not possible for rotating stars as their surface traverses the grid in between the grid points. Without this fine-tuning of the grid, the obtained frequencies would appear to randomly deviate by several Hz from the listed ones and not result in a meaningful convergence graph.
When performing the convergence test for both the f -mode and the p 1 -mode, we observe that the frequencies obtained using the 2790 × 45 grid deviate somewhat from an otherwise pretty linear convergence behaviour displayed for all other grid resolutions. At this time, it is not clear to us why this happens and we exclude the results from those simulations from the convergence test (but nonetheless show them in the graph, cf. Figure 2). Further, it is obvious that the impact of the artificial dissipation on the extracted oscillation frequencies is minuscule and can safely be ignored. Overall, the convergence Table IV. Frequencies of the 2 f -mode of the non-rotating BU0 model for different grid resolutions (the grid resolution labelled "∞" denotes extrapolation to the continuum-limit) and different dissipation coefficients for the fluid. The dissipation coefficient for the spacetime is fixed at 10 −8 . The simulation duration is 15 ms. All frequencies are given in kHz.
Grid points Dissipation coefficient 5 · 10 −7 7 · 10 −7 10 · 10 −7 ∞ test demonstrates convergence of our code and by extrapolation, we yield the frequencies 1.586 kHz for the fmode and 3.723 kHz for the p 1 -mode. Those frequencies deviate only by 0.5% from the frequencies obtained by a (naturally much more precise) eigenvalue code: they are 1.5780 kHz and 3.7068 kHz, respectively. As both convergence tests are consistent by themselves, the remaining deviation must be a systematic error, stemming from numerical inaccuracies: the most prominent candidate for error is the fact that we are using two two-dimensional codes (the rns-code for the background configuration and our time evolution code) to tackle a problem which, due to its spherical symmetry, is inherently one-dimensional; while this is mathematically not a problem, we are limited by the computational expense to somewhat lower accuracy. Furthermore, we will inevitably introduce some small interpolation error when transferring the background quantities from the grid of the rns-code to our computational grid. Normalized frequency BU0 f-mode BU0 p 1 -mode BU1 f-mode (m=-2) BU1 f-mode (m=2) BU8 f-mode (m=2) Figure 2. Graphical visualisation of the convergence tests for the five selected oscillation modes discussed in the text. The frequencies of the modes are normalised by their respective limiting value; the straight lines are linear fits to the data. The results clearly show linear convergence of the code and even for fairly coarse grids, the frequencies are recovered with only a few percent error. The two data points depicted by triangles are removed from the linear fits as outliers. For actual data refer to Tables IV, V, and VI.

Rotating Neutron Star Model
We continue with convergence tests for rotating configurations. Having begun with the non-rotating model BU0 of the BU sequence in the previous section, we now turn to its rotating members. Along this sequence, the central energy density is kept fixed at ǫ c = 0.891 · 10 15 g/cm 3 while the axis ratio r p /r e is lowered in uniform steps of 0.05 until the Kepler limit is reached [65]; this happens at an axis ratio of r p /r e = 0.576 and we call this model BUK.
We perform convergence tests for the f -mode frequency of the model BU1, which is the slowest rotating member of the BU sequence, and the model BU8, which is very close to the mass-shedding limit. As explained before, it is not possible for models of rotating stars to have their surface perfectly aligned with the grid lines (cf. also Figure 1 for the computational grid employed for the BU8 model); nonetheless, the impact of the non-alignment of the surface with the grid on the mode frequencies seems to even out over the angular direction and we observe the frequencies converging with increasing grid resolution. We list the obtained frequencies in Table VI and show them also in Figure 2. We keep the dissipation coefficient for the spacetime fixed at 10 −8 ; for the fluid, the dissipation coefficients is 10 −6 . As in the convergence test for the BU0 model, we observe only marginal dependence of the frequencies on the artificial dissipation and hence do not show those results separately.
Our convergence tests prove that our code allows us to determine the frequencies of oscillation modes to an accuracy of around 1 − 2% or better when sticking to the  [41,48]), the frequencies in full General Relativity are, as expected, lower: for this particular sequence, the frequencies of both branches experience a shift of roughly 200 − 300 Hz with the exception of the fastest rotating models BU8 and beyond; for those, the shift in frequency for the counter-rotating mode is closer to ≈ 100 Hz. This results in a slightly enhanced parameter window for the CFS-instability. Table VII. Frequencies of the counter-and co-rotating f -mode of the BU sequence as obtained from simulations on a grid with 3120 × 50 points (we have extended the sequence by the two models BU8a and BUK). We show the axis ratio, rp/re, for completeness. The simulation time is usually between 15− 30 ms; to achieve the displayed accuracy of the 2 f2-mode of the BU7 model, we evolved this particular model for 300 ms (and obtained 11.3 Hz). All frequencies are given in kHz.
Model rp/re σ( 2 f2) σ( 2 f−2) BU0 In the derivation of the evolution equations for the spacetime perturbations, we have made use of the Hilbert conditions, cf. Eq. (8), in order to simplify the lengthy equations. Nonetheless, the four gauge conditions, f µ = 0, supplement the ten field equations and should be fulfilled at all times, even though the latter are sufficient to evolve the system in time (we are ignoring the fluid perturbations for the moment as they are not affected by the gauge choice). The apparent overdetermination of the system is remedied by the fact that the field equations preserve the satisfaction of the gauge conditions if they are satisfied on the initial time slice. That this holds true has been shown for vacuum spacetimes by Barack & Lousto [66]. In the presence of neutron star matter, the proof is a little more involved but yields a comparable result. Since the divergence of the Einstein tensor vanishes, the perturbation thereof must vanish, too, δ (∇ ν G µν ) = 0. By virtue of the contracted Bianchi identities, it then follows that this condition is equivalent to the homogeneous wave equation for the divergence of φ µν . Hence, if we pick initial data such that the Hilbert conditions are fulfilled on the initial time slice and that their first time derivative vanishes, (∂ t f µ )| t=0 = 0, they will remain satisfied throughout the evolution. Clearly, the previous statement holds true only in the continuum-limit and the gauge conditions will inevitably be violated in our simulations due to various truncation and round-off errors in the numerical implementation. There is the possibility that such gauge violations grow out of control over the course of the time evolution and render it void of physical meaning. However, we do not observe such destructive behaviour in our simulations; in contrast, the gauge violation remains bounded without further ado. Furthermore, we observe that the violation decreases with increasing grid resolution. We show the evolution of the Hilbert gauge violations in Figure 3 during the simulation of a polar perturbation of a nonrotating star. More precisely, the graph shows the maximum norm of the three non-zero gauge conditions across the time slices; the three gauge violations clearly remain bounded by 3 · 10 −3 and 10 −4 , respectively. 3 The apparently severely amplified gauge violation during the first half of a millisecond of the evolution (after it started at machine precision in the initial data) can be justified by noting that we show the absolute violation rather than the one relative to the perturbation amplitude; as our initial Gaussian pulse travels toward the origin of our spherical coordinate system at the beginning of the simulation, the Gaussian pulse grows in amplitude and so does the absolute violation of the Hilbert conditions. As the graph clearly shows, this comparatively large violation is only temporary and once the energy is radiated away, it reduces to numerically expected levels.  We now turn to the Hilbert conditions concerning rotating configurations. The mathematical complexity of the problem escalates as rotation sets in, and then the Hilbert conditions (and their first time derivatives) form an involved and strongly coupled set of partial differential equations, for which we were not able to construct a straightforward appropriate solution. However, as we are currently concerned with the dynamics of fluid perturbations which are not affected by the gauge choice, a potential violation of the Hilbert conditions is irrelevant for the observables of interest. This is confirmed by the fact that, first, we observe the identical spectrum for nonrotating stars, irrespective of whether or not we choose to satisfy the gauge conditions on the initial time slice, and second, we find excellent agreement between our and previously published results (cf. Section IV C).
A numerical solution satisfying the Hilbert gauge (to numerically acceptable precision) is clearly more desirable and even indispensable when investigating w-modes or the gravitational echo of rotating neutron stars; we will, therefore, implement consistent initial data also for rotating configurations in the future, as well as implementing constraint damping mechanism [67] in order to improve the long-term accuracy and reliability of our code (even though we are already able to evolve for sev-eral hundreds of milliseconds with no problem). However, at the current stage, we are not concerned by the arbitrariness in the initial data with regards to the accuracy of the observables in our simulation as the system continues to oscillate at its natural frequencies.

C. Comparison to Published Values
While the present work is the first perturbative and comprehensive study of fluid modes of fast rotating neutron stars accounting for a dynamic spacetime, some sporadic values have been published during the past years.
The largest set of f -mode frequencies has been calculated by Zink et al. [47] (henceforth ZKSS) who investigated the CFS-instability of the f -mode in rapidly rotating polytropic stars using non-linear simulations. They constructed two previously used sequences of neutron star models with two different polytropic indices, Γ = 2 (sequence S) and Γ = 2.5 (sequence C), and computed the frequencies of the l = |m| = 2 and l = |m| = 3 f -mode for those models.
We constructed the same nine neutron star models as shown in Tables II and III in ZKSS, 4 and juxtapose the fmode frequencies obtained by our perturbative code and the published values in Tables VIII and IX. A graphical representation of the comparison is shown in Figure 4. The frequencies show excellent agreement with a deviation of less than 2 % in the majority of cases; the relative difference between the two codes is naturally considerably larger for the potentially unstable f -mode where longer evolution times are necessary to achieve the same accuracy for the f -mode frequency.   Table VIII. Comparison of the frequencies of the l = |m| = 2 f -mode obtained with the present code to the frequencies published by Zink et al. [47] for their sequences S and C. σ( 2 f2) is the frequency of the counter-rotating f -mode, σ( 2 f−2) is the frequency of the co-rotating mode. All frequencies are given in kHz. For completeness, we also report the f -mode frequen-cies of the non-rotating models as obtained by our eigenvalue code. The model S0 yields f -mode frequencies of σ( 2 f ) = 2.197 kHz and σ( 3 f ) = 2.811 kHz, respectively, while for the model C0 we report σ( 2 f ) = 2.520 kHz and σ( 3 f ) = 3.291 kHz. A comparison to the frequencies listed in Tables VIII and IX shows that all three codes produce values that are in strong agreement with each other. Beside the work by ZKSS, we are aware of only one further data point at this time. Chaurasia et al. [68] report the f -mode frequency extracted from full numerical relativity simulations of highly eccentric spinning binary neutron stars. Owing to the absence of results for fmode frequencies in full General Relativity at their time of publishing, they resorted to an intricate estimate for the f -mode frequency based on perturbation theory by starting from values and universal relations obtained in the Cowling approximation and applying rough corrections for effects of a dynamic spacetime. They ultimately arrive at a perturbation theory estimate for the frequency of the co-rotating f -mode of 2.20 kHz, which they claim to be ∼ 5% larger than their observed frequency; thus, the f -mode frequency observed in their simulation should be ∼ 2.09 kHz. For comparison, our proposed fitting formula, cf. Eq. (6) in [49], yields an estimate of 2.156 kHz.
Following their analysis, we constructed a neutron star model using the tabulated EoS SLy that has a baryon mass of M 0 = 1.504 M ⊙ and an angular rotation rate of Ω = 2π · 191 Hz; according to the rns-code, we find that this model has a central energy density of ǫ c = 0.962 · 10 15 g/cm 3 and is rotating at Ω/Ω K = 0.166 of its Kepler limit (where Ω K is the angular velocity of a neutron star at the mass-shedding limit with the same central energy density). For our test neutron star, our evolution code yields the frequencies 1.655 kHz and 2.105 kHz for the counter-and co-rotating branch of the f -mode, respectively. At a relative difference of less than 1%, the results by Chaurasia et al. [68] and our perturbative time evolution code are in excellent agreement.

V. SUMMARY AND OUTLOOK
We present 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, while our treatment of the fluid is-with the exception for the correctional terms due to the spacetime perturbations-borrowed from previous studies on fast rotating neutrons stars.
For the numerical implementation, we have constructed a particular radial grid with a varying physical resolution; it resolves the surface of oblate neutron stars well, but becomes coarser with increasing radius, allow-ing it to extend far into the wave zone with a considerably less than proportional increase in computational expense. Like in previous studies, our hydrodynamical evolution is spoiled by numerical instabilities; we apply Kreiss-Oliger dissipation to damp out those spurious oscillations. We choose 3120 × 50 grid points as our "standard resolution" and on such a grid, the dissipation coefficients required for a stable evolution are of the order of 10 −8 for the spacetime perturbations and ≈ 10 −7 − 10 −6 for the fluid perturbations. We sketch a proof that the Hilbert gauge (mathematically) remains satisfied throughout the evolution if it is satisfied in the initial data; in our simulations, we observe only very small numerical violations of the Hilbert gauge that remain bounded if we pick such initial data (even though we are not required to do so when concerned with fluid modes). We are planning to implement constraint damping mechanism for improved accuracy and also to construct initial data consistent with the Hilbert gauge for rotating configurations.
We perform convergence tests to study the accuracy of our code. With increasing grid resolution, we observe a clear convergence of the obtained frequencies toward a limit for both non-rotating and rotating background models. The convergence tests show 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. For non-rotating models, we can compare our obtained frequencies also with results from high-accuracy eigenvalue codes; the agreement between the results from our time evolution code and the values from the eigenvalue code is excellent. Furthermore, we reproduced previously published results from non-linear codes for several neutron star models at varying rotation rate with our perturbative code; again, the agreement between our values and the published ones is excellent for all considered models, irrespective of the rotation rate.
Our formulation of the problem allows to 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 CPU-hours; this enables us to study broad ranges of parameters and various EoSs. In the future, we expect to reduce the computational cost even further by deriving a simplified set of perturbation equations.
This concludes the presentation of the perturbation equations, their numerical implementation and accuracy tests of our code. Equipped with this code, we are able to produce a plethora of results concerning various oscillation modes of neutron stars; the most fundamental results concerning the fundamental mode of neutron stars are published in an accompanying article [49].
Beside the investigation of mode frequencies, we will also study the emission of gravitational waves associated with the various modes, in which the quadrupole formula seems the most suitable candidate as extraction method to complement our code.

ACKNOWLEDGMENTS
We would like to thank the anonymous referee for constructive comments that helped us improve the manuscript. This work was supported by DFG research grant 413873357. A part of the computations were performed on Trillian, a Cray XE6m-200 supercomputer at UNH supported by the NSF MRI program under grant PHY-1229408.

Appendix A: Conventions and Abbreviations
We introduce a few notational conventions and abbreviations that we use throughout all further appendices for the sake of brevity and clarity. denotes the set of all ten spacetime perturbations.
Furthermore, Q 6 is not an independent variable, too, but can be written as a linear combination of Q 1 , Q 2 , and some spacetime perturbations. To see this, we need to employ the definition of the speed of sound of the fluid, and express the perturbations of pressure and energy density in terms of the Q i 's and the spacetime perturbations.
After some algebraic manipulations, we arrive at Compared to the corresponding expression in the Cowling approximation (see, e.g., Eq. (11) in [48]), this expression has gained some correctional terms from the spacetime perturbations (and accounts for the slightly modified definition of Q 2 ) but is otherwise unaltered.

Appendix D: Equations Governing the Spacetime Dynamics
In this appendix, we present the wave equations governing the dynamics of the ten spacetime perturbations. In order to simplify the expressions, we make liberal use of the background equations, i.e. those equations that emerge from the unperturbed Einstein equations, cf. Eq. (1), as well as the Hilbert conditions, cf. Eqs. (B1) to (B4). The ten wave equations all take the same form and for any spacetime perturbation, say, X its corresponding wave equation can be written as where we denote with the symbol [X | ] the coefficient of the function appearing in the wave equation for X , and Σ 1 is the set containing all ten spacetime perturbations, their first time derivatives, their first and second radial and polar derivatives, and the fluid perturbations, i.e.