# ROXAS: a new pseudospectral non-linear code for general relativistic oscillations of fast rotating isolated neutron stars

Gaël Servignat, Jérôme Novak

► To cite this version:

Gaël Servignat, Jérôme Novak. ROXAS: a new pseudospectral non-linear code for general relativistic oscillations of fast rotating isolated neutron stars. Class.Quant.Grav., 2025, 42 (9), pp.095015. 10.1088/1361-6382/adcd1e. hal-04758037

HAL Id: hal-04758037

https://hal.science/hal-04758037v1

Submitted on 8 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
PAPER • OPEN ACCESS

ROXAS: a new pseudospectral non-linear code for
general relativistic oscillations of fast rotating
isolated neutron stars

To cite this article: Gaël Servignat and Jérôme Novak 2025 *Class. Quantum Grav.* **42** 095015

View the [article online](#) for updates and enhancements.

You may also like

- [Black Hole Waterfall: a unitary phenomenological model for black hole evaporation with Page curve](#)
  P M Alsing

- [Supermassive binary black holes and current status of their multimessenger observations](#)
  Youjun Lu

- [An asymptotic characterisation of the Kerr spacetime](#)
  R Sansom and J A Valiente Kroon
# ROXAS: a new pseudospectral non-linear code for general relativistic oscillations of fast rotating isolated neutron stars

Gaël Servignat¹,³,* and Jérôme Novak²,³

¹ Université Paris Cité, CNRS, Astroparticule et Cosmologie, F-75013 Paris, France

² Université de Strasbourg, CNRS, Observatoire astronomique de Strasbourg, UMR 7550, F-67000 Strasbourg, France

³ Laboratoire Univers et Théories, Observatoire de Paris, Université PSL, CNRS, Université Paris-Cité, 92190 Meudon, France

E-mail: servignat@apc.in2p3.fr

Received 22 October 2024; revised 11 March 2025

Accepted for publication 15 April 2025

Published 2 May 2025

<figure><img src="image_2.png" /></figure>

CrossMark

## Abstract

Next-generation gravitational wave detectors are expected to increase their sensitivity in the kHz band where binary neutron star (NS) remnants are expected to emit. In this context, robust predictions of oscillation modes of the post-merger object are desirable. To that end, we present here our code Relativistic Oscillations of non-axisymmetric neutron star (ROXAS) that is aimed at simulating isolated rotating NSs. It is based on our previously published formalism relying on primitive variables, along with the extended conformal-flatness condition formulation of Einstein equations. The equations are written under a well-balanced formulation. The code uses pseudospectral methods for both the metric and the hydrodynamics. Since standard test beds were already presented in a previous paper, we here mainly focus on presenting the numerical ingredients of ROXAS, as well as on the frequency extraction for axisymmetric and non-axisymmetric modes. We compare our results with semi-analytic perturbative methods. Spherically symmetric modes are also recovered. The code currently supports fast, rigidly rotating isolated NSs, described with a cold equation of state. We impose symmetries about the equatorial plane and by rotation of $\pi$ around the vertical z axis. These symmetries include the dominant quadrupolar $\ell = |m| = 2$ mode. A notable feature of the code is its ability to

<figure><img src="image_3.png" /></figure>

<footnote>* Author to whom any correspondence should be addressed.</footnote>
follow the surface of the star at every step. It is very lightweight and can be run
on office computers.

Keywords: hydrodynamics, primitive variables, general relativity,
3+1 formalism, pseudospectral methods, neutron stars, oscillation modes

# 1. Introduction

The LIGO-Virgo-Kagra collaboration has been observing compact binary mergers since almost a decade [1]. Three observation campaigns have led the collaboration to release a catalogue of approximately 90 events [2], and the fourth campaign is currently ongoing. Successors of current detectors, including Einstein Telescope [3–5] and Cosmic Explorer [6, 7] will have an increased sensitivity in the kHz band where the post-merger remnant of binary neutron stars (NSs) is expected to emit. Depending on the masses of the binary components, this object is thought to be a metastable hypermassive NS (HMNS), supported by its rotation, therefore surviving for a few tenths of seconds before the rotation slows down and the HMNS collapses to a black hole, see [8] for a review. Before collapsing, the HMNS emits a large quantity of gravitational waves that, if detected, would help obtain information about the interior of the NS, and therefore on the equation of state (EoS). Having robust predictions of these oscillation modes is therefore a scientific goal for the collaboration. To that end, we developed a new code, Relativistic Oscillations of non-aXisymmetric neutron stArS (ROXAS) [9], dedicated to the simulation of an isolated NS with the aim to provide oscillation mode frequencies of a HMNS. It is based for the hydrodynamics on the formalism developed in [10] that uses primitive variables (velocity and thermodynamic quantities) to write evolution equations of the enthalpy and velocity, and on the extended conformal-flatness condition formulation (xCFC) formulation of Einstein equations developed in [11]. As most existing codes [12–21] rely on conserved schemes [22], a recovery step is needed to recover the primitive variables at each time step. Being very time-consuming, and a possible source of code failure, configurations that do not exhibit shocks might rely on non-conserved formulations. The present code combines the advantages of spectral methods and primitive variables to be as lightweight as possible. A spherically symmetric version of the code has been tried and tested in [10], so that we focus on the extension of the code to two and three spatial dimensions, and test it through the extraction of axisymmetric and non-axisymmetric oscillation mode frequencies. Spherically symmetric and non-spherically symmetric oscillations were studied in the non-rotating case [23, 24], where both works are based on perturbative approaches, and the latter was later compared to non-linear codes [25].

The paper is organized as follows: in section 2, we introduce the fluid variable notations and thermodynamic properties, write the evolution equations and then reformulate them for numerical implementation. We also introduce the quadrupole formula and some of its extensions that account for general relativistic effects. In section 3, we describe the numerical ingredients of ROXAS: three-dimensional grids, numerical variables, regularization techniques, boundary conditions (BCs), domain matching, interpolation techniques and the evolution algorithm; in section 4 we show the results of the frequency extraction that we compare to previously published values. We summarize and conclude in section 5.

Unless otherwise stated, we use the set of geometrical units where $G = c = 1$ and a 4-metric $g_{\mu\nu}$ of signature $(-, +, +, +)$, of which the associated covariant derivative (connection) is denoted $\nabla_\mu$. We use Einstein's convention of summation over repeated indices. Greek indices are running from 0 to 3, whereas Latin lowercase ones are running from 1 to 3. Latin uppercase indices are used to denote various thermodynamic species present in the fluid. Unit vectors
in $\mathbb{R}^3$ are denoted with a hat, e.g. the orthonormal spherical triad is denoted $(\hat{e}_r, \hat{e}_\theta, \hat{e}_\varphi) = \begin{pmatrix} \frac{\partial}{\partial r}, \frac{1}{r} \frac{\partial}{\partial \theta}, \frac{1}{r \sin \theta} \frac{\partial}{\partial \varphi} \end{pmatrix}$.

## 2. Theoretical background

In this section we give the definitions and notations that we use to describe NS evolution. We give the thermodynamic assumptions, define the useful hydrodynamic quantities, and write the 3+1 decomposition of spacetime that we use in ROXAS. Although mostly repeated from [10], the first three subsections are provided for completeness. We then provide the set of equations for the hydrodynamic and metric sectors. Finally, the quadrupole formula, used to extract gravitational waveforms in the code, as well as some of its extension taking General Relativity into account, are presented.

### 2.1. Thermodynamics

We consider a perfect fluid carried by the timelike unitary 4-velocity $u^\mu$ that carries all species and that satisfies $u^\mu u_\mu = -1$. The energy-momentum tensor of a perfect fluid is:

$$T^{\mu\nu} = (e+p)u^{\mu}u^{\nu} + pg^{\mu\nu}, \qquad (1)$$

where $e$ is the energy density in the fluid frame and $p$ the pressure. The fluid is assumed to be composed of baryons at zero temperature and in $\beta$ equilibrium. Out-of-$\beta$-equilibrium stars can be described by adding appropriate evolution equations [10, 26]. The thermodynamics of the fluid can therefore be described with a barotropic EoS, i.e. that takes a single thermodynamic parameter as an argument. The baryons, which have a mass $m_B$, are described by their number density $n_B$. We define the log-enthalpy $H$:

$$H = \ln \left( \frac{e+p}{m_B n_B} \right), \qquad (2)$$

and the sound speed $c_s$ which, in the barotropic case, is given by:

$$c_s^2 = \frac{dp}{de}. \qquad (3)$$

### 2.2. 3+1 decomposition

We consider an asymptotically flat spacetime. The spacetime metric is given in the standard 3+1 form:

$$g_{\mu\nu}dx^\mu dx^\nu := -N^2dt^2 + \gamma_{ij}(dx^i + \beta^i dt)(dx^j + \beta^j dt), \qquad (4)$$

where $N$ is the lapse function, $\beta^i$ the shift 3-vector and $\gamma_{ij}$ the induced 3-metric on spacelike 3-hypersurfaces, associated with its covariant derivative $D_i$. The extrinsic curvature tensor is defined by

$$K_{ij} := -\frac{1}{2} \mathcal{L}_n \gamma_{ij}, \qquad (5)$$

where $\mathcal{L}_n$ is the Lie derivative along the vector $n^\mu$, a unitary, future-oriented, timelike 4-vector normal to a hypersurface $\Sigma_t$ of constant time $t$:

$$n_\mu := N \nabla_\mu t. \qquad (6)$$
$K := K_i^i$ denotes the trace of the extrinsic curvature tensor. Finally, we consider the following hydrodynamic quantities: the Lorentz factor with respect to the Eulerian observer will be denoted $\Gamma$, the Eulerian velocity $U^\mu$ and the coordinate velocity $v^i$. These quantities are related according to:

$$ \Gamma := -n_{\mu}u^{\mu} = (1 - U_{i}U^{i})^{-1/2}, \qquad (7) $$

and

$$ u^{\alpha} = \Gamma(n^{\alpha} + U^{\alpha}) = \frac{\Gamma}{N}(1, v^{i}), \quad U^{i} = \frac{1}{N}(v^{i} + \beta^{i}). \qquad (8) $$

## 2.3. Conformal decomposition and choice of coordinates

In order to fix the coordinates, we must choose a gauge and a foliation. Following [27], we use the conformal decomposition of spacelike hypersurfaces. Let $f_{ij}$ be a flat metric, $\bar{D}_i$ its associated covariant derivative, $\gamma = \det(\gamma_{ij})$ and $f = \det(f_{ij})$. We require that $\gamma_{ij}$ is $f_{ij}$ at spatial infinity. By introducing the so-called conformal factor

$$ \Psi = \left( \frac{\gamma}{f} \right)^{1/12}, \qquad (9) $$

the conformal metric is defined as

$$ \tilde{\gamma}_{ij} = \Psi^{-4} \gamma_{ij}. \qquad (10) $$

We choose the maximal slicing as a foliation:

$$ K = 0, \qquad (11) $$

and the Dirac gauge:

$$ \bar{D}_i \tilde{\gamma}^{ij} = 0. \qquad (12) $$

## 2.4. Hydrodynamics sector

The hydrodynamics is based on the formalism of [10], which uses only primitive variables. With all the hydrodynamic and thermodynamic quantities at hand, the evolution equations of the log-enthalpy $H$ and the Eulerian velocity $U_i$ are:

$$ \partial_t H = -v^j D_i H - c_s^2 \frac{\Gamma^2 N}{\Gamma^2 - c_s^2 (\Gamma^2 - 1)} \left[ K_{ij} U^i U^j + D_i U^i - \frac{U^i}{\Gamma^2} D_i H \right] \qquad (13) $$

$$ \begin{aligned} \partial_t U_i = & -v^j D_j U_i + U_j D_i \beta^j - D_i N + U_i U^j D_j N \\ & + \frac{c_s^2 N U_i}{\Gamma^2 - c_s^2 (\Gamma^2 - 1)} D_j U^j + U_i \frac{\Gamma^2 (c_s^2 - 1)}{\Gamma^2 - c_s^2 (\Gamma^2 - 1)} N U^i U^j K_{ij} \\ & - \frac{N}{\Gamma^2} \left( D_i H - \frac{\Gamma^2 (1 - c_s^2)}{\Gamma^2 - c_s^2 (\Gamma^2 - 1)} U_i U^j D_j H \right). \end{aligned} \qquad (14) $$

For a more robust numerical evolution, these equations will be rewritten in order to mathematically enforce equilibrium conditions. The first integral of motion of a general relativistic rotating self-gravitating object is, for an axisymmetric stationary spacetime [28]:

$$ D_i (H_{\mathrm{eq}} + \ln N_{\mathrm{eq}} - \ln \Gamma_{\mathrm{eq}}) = 0. \qquad (15) $$
Here, the subscript eq indicates that it describes an equilibrium state. This expression is implicitly present in (14). It is possible to make it appear explicitly to remove the sum of gradients. In order to do so, we write that the hydrodynamic, thermodynamic and metric quantities are the sum of their equilibrium value and the deviation to the equilibrium state (which is not necessarily small):

$$U_i = U_{i,eq} + \bar{U}_i, \qquad (16)$$

$$H = H_{eq} + \bar{H}, \qquad (17)$$

$$\beta^i = \beta_{eq}^i + \bar{\beta}^i, \qquad (18)$$

$$N = N_{eq} + \bar{N}. \qquad (19)$$

Note that since the conformal factor $\Psi$ does not appear explicitly in the hydrodynamic equations, we do not need to write $\Psi = \Psi_{eq} + \bar{\Psi}$. In particular, the coordinate velocity $v^i$ is transformed into

$$v^i = v_{eq}^i + \bar{v}^i, \qquad (20)$$

where

$$v_{eq}^i = N_{eq}U_{eq}^i - \beta_{eq}^i \qquad (21)$$

$$\bar{v}^i = -\bar{\beta}^i + N_{eq}\bar{U}^i + \bar{N}U_{eq}^i + \bar{N}\bar{U}^i. \qquad (22)$$

The equilibrium quantities are static: $\partial_t U_{i,eq} = \partial_t H_{eq} = 0$. Two more equilibrium conditions can be enforced, based on the assumption of rigid rotation. In this particular case, the coordinate velocity components are, in a spherical triad: $v^r = v^\theta = 0$, $v^\varphi = \Omega r \sin\theta$, with $\Omega$ a rotation rate in rad · s⁻¹ (if f is the rotation frequency in Hz, $\Omega = 2\pi f$). The conditions are:

$$D_j U_{eq}^j = 0, \qquad (23)$$

$$v_{eq}^j D_j U_{i,eq} + U_{j,eq} D_i v_{eq}^j = 0. \qquad (24)$$

The condition (24) can be shown analytically when assuming that the conformal metric is flat; we will explain in the next subsection that we place the present study within this assumption. Then, we write that $\beta^i = N U^i - v^i$, use the identity $D_i \Gamma = \Gamma^3 U_j D_i U^j$ and inject the conditions (15), (23) and (24) to obtain the evolution equation of $\bar{U}_i$:

$$\begin{aligned} \partial_t \bar{U}_i = & -[v_{eq}^j D_j \bar{U}_i + \bar{v}^j D_j U_{i,eq} + \bar{v}^j D_j \bar{U}_i + U_{j,eq} D_i \bar{v}^j + \bar{U}_j D_i v_{eq}^j + \bar{U}_j D_i \bar{v}^j \\ & -N(U_{j,eq} D_i \bar{U}^j + \bar{U}_j D_i U_{eq}^j + \bar{U}_j D_i \bar{U}^j)] - \frac{N}{\Gamma^2} D_i \left(\bar{H} + \ln\left(1 + \frac{\bar{N}}{N_{eq}}\right)\right) \\ & + U_i U^j D_j N + \frac{N}{\Gamma^2} \frac{\Gamma^2 (1 - c_s^2)}{\Gamma^2 - c_s^2 (\Gamma^2 - 1)} U_i U^j D_j H \\ & + \frac{c_s^2 N U_i}{\Gamma^2 - c_s^2 (\Gamma^2 - 1)} D_j \bar{U}^j + U_i \frac{\Gamma^2 (c_s^2 - 1)}{\Gamma^2 - c_s^2 (\Gamma^2 - 1)} N U^i U^j K_{lj}. \end{aligned} \qquad (25)$$

The evolution equation of $\bar{H}$ is:

$$\begin{aligned} \partial_t \bar{H} = & -\bar{v}^j D_i H_{eq} - v_{eq}^i D_i \bar{H} - \bar{v}^j D_i \bar{H} \\ & - c_s^2 \frac{\Gamma^2 N}{\Gamma^2 - c_s^2 (\Gamma^2 - 1)} \left[ K_{ij} U^i U^j + D_i \bar{U}^i - \frac{U^i}{\Gamma^2} D_i H \right]. \end{aligned} \qquad (26)$$
This way to write the evolution equation is referred to in the literature as a *well-balanced* formulation. Since we evolve only the ‘perturbation’, this formulation helps exploit all digits of the numerical representation of double precision floating-points, and allows as shown above to enforce equilibrium conditions exactly. Dumbser et al [29] provide a short review on well-balanced formulations. The evolution equations (25) and (26) are the ones that are implemented in the code.

## 2.5. Metric sector

The metric quantities $N$, $\Psi$, $\beta^i$ and $\gamma_{ij}$ all enter (explicitly or implicitly) the hydrodynamics of a perfect fluid. Moreover, the spacetime responds to changes of its matter content, meaning that evolving $H$ and $U_i$ should have an effect on the metric. The Einstein equations are therefore needed to compute the metric along the hydrodynamics. We use the fully constrained formalism of [27] and assume the CFC, which amounts to set $\tilde{\gamma}^{ij} = f^{ij}$. We use the extended formulation of [11]; the Einstein equations reduce to a hierarchical set of two linear vector and two non-linear scalar Poisson-like partial differential equations (PDEs). As described in [10], the code comprises two distinct numerical representations of the conformal factor $\Psi$. The first one, denoted $\Psi_{ev}$, is obtained from the following evolution equation:

$$ \partial_t \ln \Psi_{ev} = \beta^i \bar{D}_i \ln \Psi_c + \frac{1}{6} \bar{D}_i \beta^i. \qquad (27) $$

$\Psi_c$ is the second representation of the conformal factor, and is obtained from the numerical solution of the Hamiltonian constraint of the Einstein equations. Since the CFC kills all radiative degrees of freedom of the spacetime, the code does not naturally produce gravitational waves. Gravitational waves can still be extracted from the simulations thanks to Einstein's quadrupole formula.

## 2.6. Gravitational waves extraction

To make the dimension of the formulas clearer, we bring back $c$ and $G$ in this subsection. The quadrupole formula, as derived by Einstein [30], relates the gravitation quadrupole radial field $h_{ij}^{TT}$ to the variations of the quadrupole mass tensor $Q_{ij}$, and is given by:

$$ h_{ij}^{TT}(\mathbf{x}, t) = \frac{2G}{c^4 r} P_{ij}^{kl}(\mathbf{n}) \ddot{Q}_{kl} \left(t - \frac{r}{c}\right), \qquad (28) $$

where $\mathbf{n} = \mathbf{x}/r$ is the direction of propagation, $P_{ij}^{kl}$ a projection operator on the direction of propagation, and the dots correspond to the derivative with respect to time. The quadrupole formula is defined in a purely Newtonian context, and $\rho = m_B n_B$ is the rest-mass density of the matter. It has one flaw from the numerical point of view: it involves a second derivative with respect to time. Numerical derivatives tend to be expensive to evaluate and to degrade the accuracy of the representation, especially time derivatives (in contrast with spatial derivatives which, in this context, are computed with spectral methods). To take advantage of the spectral representation, we follow the formalism of [31] to replace the first spatial derivative with a well-behaved integral:

$$ \frac{d}{dt} Q_{ij} = \int_{\mathbb{R}^3} \rho^* \left( x_i v_j + x_j v_i - \frac{2}{3} \delta_{ij} x_k v^k \right) d^3 x, \qquad (29) $$

where $\rho^* = \Psi^6 \Gamma \rho$ is the coordinate rest-mass density, $(x_i)_{i=1,2,3}$ a coordinate system (e.g. $x_1 = x$, $x_2 = y$, $x_3 = z$ for the Cartesian coordinates), and $v_i$ is the coordinate velocity. Equation (6.5) of [31] further expresses the second derivative of $Q_{ij}$ as an integral of a Newtonian-like
potential (i.e. that is a solution of the Poisson equation $\Delta\Phi = 4\pi G\rho^*$), but we found that this expression (the so-called stress formula) would introduce an offset in the waveforms. It is likely that the higher-order terms that were neglected to write the formula contribute significantly in our very compact configurations. For the purpose of the present study we found that a first order backwards scheme to compute the second time derivative was sufficient to produce clean waveforms and leave the improvements of the formalism for future investigations. We then use the plus/cross polarization decomposition of $h_{ij}^{TT}$ (see e.g. [32]):

$$h_{ij}^{TT}(\mathbf{x}, t) = \frac{1}{r} (A_+ \hat{\mathbf{e}}_+ + A_\times \hat{\mathbf{e}}_\times), \qquad (30)$$

where

$$\hat{\mathbf{e}}_+ = \hat{\mathbf{e}}_\theta \otimes \hat{\mathbf{e}}_\theta - \hat{\mathbf{e}}_\varphi \otimes \hat{\mathbf{e}}_\varphi, \qquad (31)$$

$$\hat{\mathbf{e}}_\times = \hat{\mathbf{e}}_\theta \otimes \hat{\mathbf{e}}_\varphi + \hat{\mathbf{e}}_\varphi \otimes \hat{\mathbf{e}}_\theta, \qquad (32)$$

with $\otimes$ the tensor product. If the observer is in the equatorial plane with respect to the source, and denoting $\ddot{I}_{ij} = \frac{2G}{c^4}\ddot{Q}_{ij}$, the expressions of $A_+$ and $A_\times$ are:

$$A_+^e = \ddot{I}_{zz} - \ddot{I}_{yy}, \qquad (33)$$

$$A_\times^e = -2\ddot{I}_{yz}. \qquad (34)$$

The final formula is:

$$A_+^e = \frac{G}{c^4} \frac{d}{dt} \left\{ \int_{\mathbb{R}^3} \rho^* (zv_z - yv_y) d^3 x \right\}, \qquad (35)$$

$$A_\times^e = -\frac{2G}{c^4} \frac{d}{dt} \left\{ \int_{\mathbb{R}^3} \rho^* (zv_y + yv_z) d^3 x \right\}. \qquad (36)$$

### 3. Numerical representation

The numerical implementation of ROXAS is based upon LORENE [33], a dedicated C[4]++ library that comes with built-in three-dimensional Poisson solvers, based on Chebyshev pseudospectral methods [34–37], and using spherical coordinates. Pseudospectral methods are used to compute the spatial derivatives in the sources of the evolution and structure equations. The ideas described in section 4 of [10] are mostly the same here, up to one major difference: the present three-dimensional code comprises two distinct grids. The first one, referred to as the *metric grid* is used to solve the Einstein equations, and the second one, referred to as the *hydro grid*, is where the primitive variables are evolved. Since the hydro grid follows the star’s surface that is not spherically symmetric, and since the three-dimensional Poisson solvers rely on the sphericity of the numerical domains, the metric grid and the hydro grid need to be distinct. We here give details on each of those two grids.

#### 3.1. Metric grid

The metric grid consists in a spherical nucleus, an arbitrary number of spherical shells and a compactified external domain (CED). If $N$ is the total number of domains, and $(r_i)_{i=1,...,N-1}$ a set of $N-1$ real numbers such that

$$0 < r_1 < \dots < r_{N-1} < +\infty, \qquad (37)$$
then the nucleus spans $[0, r_1] \times S$, the ith shell spans $[r_i, r_{i+1}] \times S$ and the CED spans $[r_{N-1}, +\infty] \times S$, where $S = \{(\theta, \varphi), \theta \in [0, \pi], \varphi \in [0, 2\pi)\}$. The physical variables $(r, \theta, \varphi)$ are mapped to the numerical coordinates $(\xi, \theta', \varphi')$ such that we always have $\theta = \theta'$ and $\varphi = \varphi'$, while $\xi \in [0, 1]$ in the nucleus and $\xi \in [-1, 1]$ in the shells and the CED. The radial variable is changed according to

$$ 
\left\{ 
\begin{aligned}
r &= r_1\xi, && (\text{nucleus}) \\
r &= \frac{r_{i+1} - r_i}{2}\xi + \frac{r_{i+1} + r_i}{2}, && (\text{ith shell}) \\
r &= \frac{2r_{N-1}}{1-\xi}. && (\text{CED})
\end{aligned}
\right. 
\qquad (38) 
$$

This changes the radial differential operator:

$$ 
\left\{ 
\begin{aligned}
\frac{\partial}{\partial r} &= \frac{1}{r_1} \frac{\partial}{\partial \xi}, && (\text{nucleus}) \\
\frac{\partial}{\partial r} &= \frac{2}{r_{i+1} - r_i} \frac{\partial}{\partial \xi}, && (\text{ith shell}) \\
\frac{\partial}{\partial r} &= \frac{(\xi - 1)^2}{2r_{N-1}} \frac{\partial}{\partial \xi}. && (\text{CED})
\end{aligned}
\right. 
\qquad (39) 
$$

The metric grid is used for several purposes. The main one is to solve for the metric. The metric quantities $\beta^i$, $N$, $\Psi_c$ and $\Psi_{ev}$ are defined on the metric grid. The source of equation (27) is computed on the metric grid so that $\Psi_{ev}$ is evolved there. The metric grid also serves to store the equilibrium values of the initial profiles $H_{eq}$, $U_{i,eq}$, $N_{eq}$ and $\beta_{eq}^i$. The angular directions are set by the number of grid points in the $\theta$ and $\varphi$ directions, respectively denoted $N_\theta$ and $N_\varphi$. We consider the spatial symmetry that leaves scalar fields unchanged under the transformation

$$ (\theta, \varphi) \rightarrow (\pi - \theta, \varphi + \pi). \qquad (40) $$

Under this assumption, the angular variables only need to span $\theta \in [0, \pi/2]$ and $\varphi \in [0, \pi)$, and the angular collocation points are chosen at

$$ \forall \ell \in \{0, \dots, N_{\theta} - 1\}, \theta_{\ell} = \frac{\ell}{N_{\theta} - 1} \frac{\pi}{2}, \qquad (41) $$

$$ \forall m \in \{0, \dots, N_{\varphi} - 1\}, \varphi_{m} = \frac{m}{N_{\varphi}} \pi. \qquad (42) $$

Since in this setting the angular variables span smaller intervals, fewer points are needed to achieve the same resolution as in the non-symmetric case. We emphasize here that imposing this symmetry reduces the angular basis to spherical harmonics with only even quantum numbers $\ell$ and $m$. Since the dominant quadrupolar modes correspond to $\ell = |m| = 2$, this restriction allows to study this particular dominant mode. These conditions also do not refrain from describing the whole space as the basis functions can be evaluated everywhere. All domains share the same angular directions.

## 3.2. Hydro grid

The hydro grid consists in a spherical nucleus and an arbitrary number of shells. The last shell's outer boundary is deformed so that it is adapted to the surface of the star. Since most simulations were run with exactly one shell we here write the mapping between the physical and numerical variables for two domains in the hydro grid. The case where there are no shells is only used for spherically symmetric simulations and the code described in [10] is to be preferred in this case. Let $r_{nuc}$ be the radius of the nucleus, and $R_S(\theta)$ the equilibrium coordinate
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Cross-section view of the rotating BU4 model (see table 1) in the (x, z) plane. The nucleus is spherical and the shell is adapted to the surface of the star. Each dashed line represents one of the interior angular directions of the grid.</figcaption></figure>

radius profile of the rotating star. At equilibrium, the nucleus spans $[0, r_{nuc}] \times S$ and the shell spans $[r_{nuc}, R_S(\theta)] \times S$. This is illustrated in figure 1. An important feature of the grid is that the angular directions of the hydro grid are chosen to be the same as the ones of the metric grid. This helps increase the efficiency of interpolation procedures that we will describe in a further subsection. As the star evolves and oscillates, the outer boundary of the shell co-moves with the surface of the star. In order to do so, we impose that the value of $r$ on the outer boundary of the shell, denoted $R(t, \theta, \varphi)$, is evolved with the following PDE:

$$ \frac{\partial R}{\partial t} = v_r(R(t, \theta, \varphi), t) - \frac{v_\theta(R(t, \theta, \varphi), t)}{R} \frac{\partial R}{\partial \theta} - \frac{v_\varphi(R(t, \theta, \varphi), t)}{R \sin\theta} \frac{\partial R}{\partial \varphi}. \quad (43) $$

The derivation of this equation is detailed in appendix A, and the results of numerical experiments are shown in appendix B: we perform simulations of rotating ellipsoids and check for spectral convergence of the space scheme as well as convergence of the time integration scheme.
In the nucleus we map the physical coordinates $(t, r, \theta, \varphi)$ to the numerical coordinates $(t', \xi, \theta', \varphi')$ by:

$$ \left\{ \begin{array}{l} t = t', \\ r = r_{\text{nuc}}\xi, \\ \theta = \theta', \\ \varphi = \varphi'. \end{array} \right. \qquad (44) $$

The differential operators are then changed to:

$$ \left\{ \begin{array}{l} \displaystyle \frac{\partial}{\partial t} = \frac{\partial}{\partial t'}, \\ \displaystyle \frac{\partial}{\partial r} = \frac{1}{r_{\text{nuc}}} \frac{\partial}{\partial \xi}, \\ \displaystyle \frac{\partial}{\partial \theta} = \frac{\partial}{\partial \theta'}, \\ \displaystyle \frac{\partial}{\partial \varphi} = \frac{\partial}{\partial \varphi'}. \end{array} \right. \qquad (45) $$

In the shell we set:

$$ \left\{ \begin{array}{l} t = t', \\ r = \frac{R(t', \theta', \varphi') - r_{\text{nuc}}}{2}\xi + \frac{R(t', \theta', \varphi') + r_{\text{nuc}}}{2}, \\ \theta = \theta', \\ \varphi = \varphi'. \end{array} \right. \qquad (46) $$

Then the differential operators become:

$$ \left\{ \begin{array}{l} \displaystyle \frac{\partial}{\partial t} = \frac{\partial}{\partial t'} - \frac{2v_g}{R(t', \theta', \varphi') - r_{\text{nuc}}} \frac{\partial}{\partial \xi}, \\ \displaystyle \frac{\partial}{\partial r} = \frac{2}{R(t', \theta', \varphi') - r_{\text{nuc}}} \frac{\partial}{\partial \xi}, \\ \displaystyle \frac{\partial}{\partial \theta} = \frac{\partial}{\partial \theta'} - \frac{\xi + 1}{R(t', \theta', \varphi') - r_{\text{nuc}}} \frac{\partial R}{\partial \theta'} \frac{\partial}{\partial \xi}, \\ \displaystyle \frac{\partial}{\partial \varphi} = \frac{\partial}{\partial \varphi'} - \frac{\xi + 1}{R(t', \theta', \varphi') - r_{\text{nuc}}} \frac{\partial R}{\partial \varphi'} \frac{\partial}{\partial \xi}. \end{array} \right. \qquad (47) $$

We have denoted with $v_g$ the grid advection velocity that arises from using a moving grid (see [10] for some more insight on grid advection). Its expression is

$$ v_g = \frac{1}{2}(\xi + 1) \frac{\partial R}{\partial t'}, \qquad (48) $$

and can therefore be computed thanks to the equation (43).

### 3.3. Scalar fields regularity

The simulations can exhibit a diverging behavior if two regularity conditions on scalar fields are not met. They are linked to the singularity of spherical coordinates on the origin and on the vertical axis. These regularization processes can be thought of as BC consistency on the rectangle $(\theta, \varphi) \in [0, \pi/2] \times [0, \pi]$ and we explain hereafter why they are needed. An alternative to our approach is, for instance, the use of stereographic and gnomic maps to cure angular singularities related to spherical coordinates, as described in [38].
**3.3.1. Regularization at the origin.** The regularization scheme described in this paragraph was introduced by Nicoules in his PhD thesis [39]. The spherical coordinates are degenerate at the origin, and therefore we expect that this particular physical point must receive a particular attention. An analytical field must not be multivalued at the origin; numerically, the point $r = 0$ is represented by as many grid point as there are angular directions. For example, if we choose $N_\theta = 5$ and $N_\varphi = 4$, there are 20 grid points that represent the same $r = 0$ physical point, and those grid points have a priori no reason to remain identically valued during a simulation. By changing how the division by $\xi$ is done in the nucleus, the divergences linked to multivalued fields at the origin can be cured. By default, the division by $r$ is done on the regular part of the field, i.e. if we want to divide a given function $f(r)$ by $r$, the actual division that is done is $(f(r) - f(0))/r$. As in the nucleus, the radial variable $r$ is changed to the numerical variable $\xi$ by $r = r_{nuc}\xi$, this corresponds to a division by $\xi$, which is done in the coefficient space, thanks to the recurrence relation between the Chebyshev polynomials. To show the new procedure, we assume that $f$ is an even scalar field (if it is odd $f(0) = 0$ by construction and there is no particular treatment):

$$f(\xi) = \sum_{k=0}^{N_r-1} a_i T_{2i}(\xi), \qquad (49)$$

with $T_{2i}$ the even Chebyshev polynomials and $N_r$ the number of collocation points in the radial direction. Then, the expression $f(r) - f(0)$ can be interpreted in the coefficient space as $f(\xi) - f(0)T_0(\xi)$, which amounts to change the first coefficient of the spectral representation of $f$. This operation can be performed on the last coefficient by writing $f(\xi) \leftarrow f(\xi) - f(0)T_{2N_r-1}(\xi)/T_{2N_r-1}(0)$, which looks on the coefficient like:

$$a_{N_r-1} \leftarrow a_{N_r-1} - (-1)^{N_r-1} \sum_{k=0}^{N_r-1} (-1)^k a_k. \qquad (50)$$

This procedure has been implemented as the default one in LORENE.

**3.3.2. Regularization on the vertical axis.** Once again, analytical fields must not be multi-valued on the vertical axis; a grid point on the vertical axis has $N_\varphi$ representatives. If a regular scalar field $f$ is expanded on the spherical harmonics basis like

$$f(r, \theta, \varphi) = \sum_{\ell=0}^{\ell_{\max}} \sum_{m=-\ell}^{\ell} f_{\ell m}(r) Y_{\ell m}(\theta, \varphi), \qquad (51)$$

then for $m \neq 0$, $f_{\ell m} Y_{\ell m}$ must behave as $\sin^m(\theta)$ when $\theta \to 0$. Because $f_{\ell m}$ only depends upon the radial variable and $Y_{\ell m} \propto \sin^m(\theta)$, the condition is automatically fulfilled when the fields are represented on the spherical harmonics basis. Since derivation routines in LORENE are written on the Fourier basis of sines and cosines, the physical fields are represented on this basis by default. This means that the regularization procedure consists in simply performing a basis change on scalar fields from Fourier to spherical harmonics, and back to Fourier again. The transformation costs $O(N \log N)$ in each direction so if $N = \max(N_r, N_\theta, N_\varphi)$, this amounts to a procedure that costs $O((N \log N)^3)$.

### 3.4. Filters

As we have seen in section 3.2, the angular differential operators involve an additional Jacobian term which depends on the angular variables. Since it is non-linear in $R$, it may provoke aliasing
issues as the multiplication is performed in the collocation space. To cure the issues associated with the Jacobian, we introduce exponential filters [40, 41] that we apply on the spectral representation ($f_{i\ell m}$) of a scalar field $f$ expanded on Chebyshev polynomials and spherical harmonics (which is simply considering in equation (51) the expansion of $f_{\ell m}$ on Chebyshev polynomials). Given three positive integers ($p_r, p_\theta, p_\varphi$), we apply the following operation:

$$ \mathcal{F}(f_{i\ell m}) = f_{i\ell m} e^{-\alpha[i/(N_r-1)]^{pr}} e^{-\alpha(\ell/\ell_{\max})^{p\theta}} e^{-\alpha(m/N_\varphi)^{p\varphi}}, \quad (52) $$

where $\alpha = \ln(10^{-16})$ so that the last coefficient is always set below the machine accuracy. It is possible to choose not to filter in a given direction; typically the enthalpy and velocity profiles are not filtered in the radial direction. The radial filter is used only on the sources of the Einstein equations. This procedure can be directly applied on scalar fields such as the enthalpy profile, but to filter vector fields that are by default represented with the spherical triad, we must first rotate the triad to get the Cartesian components of the vector field, filter those components and rotate the triad back to spherical.

### 3.5. BCs

#### 3.5.1. Angular BCs.
Since the angular variables $(\theta, \varphi)$ span $[0, \pi/2] \times [0, \pi)$, the numerical domain has borders where $\theta = 0, \pi/2$ and $\varphi = 0, \pi$. These numerical borders are of course not physical border, and the angular BCs are taken care of with the Galerkin method that consists in choosing a basis that fulfills the desired BC so that it is automatically fulfilled by the numerical representation of the field.

#### 3.5.2. External BCs.
Here we describe the BCs imposed on the outermost boundaries of the numerical domains. On the metric grid, we impose the conditions that $\Psi \to 1, N \to 1$ and $\beta^i \to 0$ for $r \to +\infty$ (corresponding to $\xi \to 1$ in the CED) as the spacetime is asymptotically flat. On the hydro grid, we impose the impermeable BC that requires that no matter crosses the numerical domain outermost boundary. It is already accounted for by choosing the comoving grid with equation (43) (see also appendix A).

#### 3.5.3. Domain matching.
Now the communication between the domains remains to be addressed. On the metric grid, the metric fields $\Psi_c, N$ and $\beta^i$ are matched directly by the Poisson solvers which require the continuity of the functions and their first radial derivative. Equation (27) is an advection equation with characteristic speed $-\beta^r$ and is matched in an upwind fashion: given two neighboring domains $l_z$ and $l_z + 1$ and an angular direction, if $\beta^r < 0$ the value of the last grid point of the domain $l_z$ is copied to the first grid point of the domain $l_z + 1$, and vice-versa if $\beta^r > 0$. On the hydro grid we match the enthalpy profile with the tau method. The deformed domain along with equation (47) implies that the angular coefficients of $R(\theta, \varphi)$ are mixed with those of $\partial_\xi H$ so the spectral representation of $H$ cannot be used directly for this purpose. Instead, we perform the inverse Fourier transform on the angular coefficients of $H$ so that the resulting object is a collection of radial functions $\tilde{H}_{\ell m}(r) = H(r, \theta_\ell, \varphi_m)$, where $\theta_\ell$ and $\varphi_m$ are defined in equation (41). Then, every $\tilde{H}_{\ell m}$ is matched according to the procedure described in appendix C. Assuming a subsonic flow, there is one ingoing characteristic at the outer boundary of the nucleus and one at the inner boundary of the shell. We therefore impose the continuity of all $\tilde{H}_{\ell m}$ and $\partial_\ell \tilde{H}_{\ell m}$.

Even though it is not motivated mathematically, we found that the equilibrium term $\tilde{H} + \ln(1 + \tilde{N}/N_{eq})$ (see equation (14)) also needed to be matched when computing the sources of the hydrodynamic equations. We employ the procedure of appendix C again.
### 3.6. Evolution algorithm

We presented all the technical components of the ROXAS code and now turn to sum them up in a description of the successive steps to perform the simulation. Since there are two distinct grids for the metric and the hydrodynamics, a number of interpolation steps are involved to compute the sources of the Einstein equations and of the hydrodynamic equations, that we detail hereafter.

* Computation of the source terms of the metric equations (see equations (30)–(33) of [11]).

(i) We first interpolate $\Psi_{ev}$ to the hydro grid to define $\gamma_{ij} = \Psi_{ev}^4 f_{ij}$ there, which allows to compute the Lorentz factor $\Gamma = (1 - U_i U^i)^{-1/2} = (1 - \gamma_{ij} U^i U^j)^{-1/2}$, and then to compute the matter terms $E = \Gamma^2(e+p) - p$, $S = 3p + (\Gamma^2 - 1)(e+p)$ and $p_j = \Gamma^2(e+p)U_j$.

(ii) Then we interpolate $E$ and $p_j$ to the metric grid and form the quantities $E^* = E\Psi_{ev}^6$ and $p_j^* = p_j\Psi_{ev}^6$.

(iii) The first two metric equations can be solved. $\Psi_c$ is now known, and we compute $S^* = S\Psi_c^6$ where $S$ has been computed on the hydro grid and then interpolated to the metric grid. The last two metric equations can be solved.

* Computation of the source terms of the hydrodynamic equations.

(i) We interpolate $\Psi_c$ to the hydro grid to define $\gamma_{ij} = \Psi_c^4 f_{ij}$.

(ii) We interpolate the lapse and the shift on the hydro grid, but this time we separate the equilibrium and the bar quantities. $N_{eq}$ and $\beta_{eq}^i$ are simply interpolated from the metric to the hydro grid, while to compute $\bar{N}$ and $\bar{\beta}^i$ we interpolate $N - N_{eq}$ and $\bar{\beta}^i - \bar{\beta}_{eq}^i$: these differences are computed on the metric grid before interpolation. The total lapse and shift on the hydro grid are formed with $N_{eq} + \bar{N}$ and $\beta_{eq}^i + \bar{\beta}^i$ where we use the interpolated values previously computed.

Then, the algorithm to perform a time step is the following:

(i) At a given time $t$, $H$, $U_i$ and $\Psi_{ev}$ are known thanks to their evolution equations (14), (13) and (27). $\Gamma$ can be computed from $U_i$, and along with the EoS, $H$ allows computing $p$ and $e$.

(ii) Following the interpolation procedure we compute the sources of the Einstein equations and solve them with the relaxation iterative method. The output is the value of $\Psi_c$ consistent with the Hamiltonian constraint, which replaces the value of $\Psi_{ev}$ computed with equation (27) as soon as it is available, as well as the lapse $N$ and the radial component of the shift vector $\beta^r$. A number of simulations will be run using the Cowling approximation which amounts to fix the spacetime to its initial equilibrium value; in this case the current step is skipped.

(iii) We use the knowledge of the radius of the star $R$, as well as $H$, $U_i$, $N$, $\beta^i$ and $\Psi_c$ together with the EoS and the interpolation procedure to compute the sources of equations (14), (13), (43) and (27).

(iv) $H$, $U_i$, $R$ and $\Psi_{ev}$ are evolved with an explicit time integration scheme.

(v) The hydro grid of the next time step is computed.

(vi) The BCs are imposed as described in section 3.5.

The evolution is performed with an Adams-Bashforth scheme of order 3 (AB3 scheme), which uses the knowledge of the current and previous two time steps to compute the next
one. We therefore initialize the first two time steps with a Runge–Kutta scheme of order 4 (RK4 scheme). The reason why we choose the Adams–Bashforth method lies in this multistep nature. All the data relative to the previous time steps are already computed and can be stored, whereas the exclusive use of an RK4 scheme demands to compute several sub-steps to compute a single time step. In each sub-step we have to compute for the metric again, which significantly increases the computation time even though the RK4 scheme allows for greater time step lengths than the AB3 scheme. All initial data are computed using the initial data solvers of LORENE [42, 43]. Since we compute the equilibrium state of rigidly rotating stars, the initial equilibrium state is an axisymmetric (therefore two-dimensional) configuration. To trigger the oscillations, the star is perturbed with the profile given by:

$$ \delta H = \varepsilon \left( x^2 - y^2 - \frac{2}{3} r^2 \right) \left( 1 - \frac{r^2}{R_S(\theta)^2} \right). \qquad (53) $$

It is a mixture between $\ell = m = 0$ and $\ell = m = 2$ modes. The radial dependence of the perturbation profile is arbitrary, but we choose it to be a low order polynomial so that it is easily represented by the spectral decomposition, as well as fulfilling the initial condition $H = 0$ on the border of the star. The parameter $\varepsilon$ allows to control the strength of the perturbation. Note that, when this perturbation is added, initial data depend on all three coordinates $(r, \theta, \varphi)$ and are thus fully three-dimensional.

## 4. Tests

To study oscillation modes of rotating stars, we use the widely studied sequence of uniformly rotating polytropic stars BU [15, 44–48]. They consist in a fixed polytropic EoS described with $\kappa = 100$ and $\gamma = 2$, so that the pressure $p$ and baryon number density $n_B$ are linked by $p = \kappa n_B^{\gamma}$, and a fixed central enthalpy $H_c = 0.227823853$. The members of the sequence are numbered from 0, the larger the number the greater the rotation rate. We here study the BU sequence starting at BU0, which is non-rotating, to BU7 which spins at the rate of 831 Hz, close to the Keplerian frequency of the star (see table 1). In particular, [15] gives values for the frequencies of axisymmetric modes in CFC, [45, 46] give values for the frequencies of axisymmetric and non-axisymmetric modes using the Cowling approximation and [48] gives values for the frequencies of non-axisymmetric modes in full GR. We give the macroscopic properties of the sequence BU in table 1. Even though the mode frequencies are present in the waveforms computed with equations (35) and (36), we extract the frequencies from the appropriate coefficient of the decomposition of the radius on spherical harmonics, i.e. quasi-radial modes are extracted from the $\ell = m = 0$ coefficient, axisymmetric modes from the $\ell = 2, m = 0$ coefficient and non-axisymmetric modes from the $\ell = |m| = 2$ coefficients. The grid settings are the following:

* The hydro grid is composed of two domains: a spherical nucleus which is fixed and set to span [0, 0.7$r_p$] where $r_p$ is the coordinate radius at $\theta = 0$ i.e. the North Pole, and a shell that describes the rest of the star. We set $N_r = 17$ in both domains, and $N_\theta$ is set as stated in table 1 (except for BU0 where we set $N_\theta = 9$). The 2D simulations are run with $N_\varphi = 1$ and all 3D simulations are run with $N_\varphi = 8$.

* In the metric grid, we set a nucleus and a shell that respectively span [0, 8 km] and [8 km, $r_p$]. Next we choose an arbitrary number of shells so that they span intervals of the form $[R_S(\theta_i), R_S(\theta_j)]$, $i < j$. Beyond $r_e$ we set a final shell and the CED. In the first two domains
<table><caption>**Table 1.** Properties of the BU sequence as computed with the initial data solvers of LORENE [42, 43]. All stars are built with a central enthalpy of $H_c = 0.227823853$. $N_\theta$ is the number of points in the $\theta$ direction needed to describe the surface of the star. We give the gravitational mass $M$ and the circumferential equatorial radius $R$, as well as the ratio between the coordinate radii at the North Pole $r_p$ and at the equator $r_e$. $\Omega_K = 2\pi \times 853 \text{ rad} \cdot \text{s}^{-1}$ is the Kepler frequency of the star. The ratio $\Omega/\Omega_K$ is given in %.</caption><thead><tr><th>Model</th><th>f(Hz)</th><th>Ω/Ω<sub>K</sub></th><th>M[M<sub>☉</sub>]</th><th>R(km)</th><th>r<sub>e</sub>(km)</th><th>r<sub>p</sub>/r<sub>e</sub></th><th>N<sub>θ</sub></th></tr></thead><tbody><tr><td>BU0</td><td>0</td><td>0</td><td>1.400</td><td>14.15</td><td>11.99</td><td>1</td><td>1</td></tr><tr><td>BU1</td><td>347</td><td>41</td><td>1.431</td><td>14.51</td><td>12.30</td><td>0.950</td><td>9</td></tr><tr><td>BU2</td><td>487</td><td>57</td><td>1.465</td><td>14.92</td><td>12.64</td><td>0.900</td><td>9</td></tr><tr><td>BU3</td><td>590</td><td>69</td><td>1.502</td><td>15.37</td><td>13.03</td><td>0.850</td><td>17</td></tr><tr><td>BU4</td><td>673</td><td>79</td><td>1.542</td><td>15.90</td><td>13.49</td><td>0.800</td><td>17</td></tr><tr><td>BU5</td><td>740</td><td>87</td><td>1.585</td><td>16.51</td><td>14.03</td><td>0.750</td><td>25</td></tr><tr><td>BU6</td><td>793</td><td>93</td><td>1.627</td><td>17.26</td><td>14.70</td><td>0.700</td><td>25</td></tr><tr><td>BU7</td><td>831</td><td>97</td><td>1.665</td><td>18.16</td><td>15.55</td><td>0.649</td><td>33</td></tr></tbody></table>

we set $N_r = 25$; in the shells spanning $[r_p, r_e]$ we set $N_r = 25$ or 33 depending on the size of the domain; we finally set $N_r = 25$ in the last shell and $N_r = 17$ in the CED. The angular directions of the metric grid are the same as the hydro grid and therefore the values of $N_\theta$ and $N_\phi$ are identical.

We once again employ exponential angular filters in the hydro grid. The value of the order for the $\theta$ filter is decreased (i.e. its strength is increased) as $N_\theta$ increases; we use an order 10 filter for $N_\theta = 9$ and 17, an order 8 filter for $N_\theta = 25$ and an order 6 filter for $N_\theta = 33$. When $N_\phi > 1$, we always set $N_\phi = 8$ and in that case we use an order 12 filter. We filter the enthalpy and Eulerian velocity profiles, as well as the sources of the Einstein equations, as described in section 3.4. For the non-rotating star BU0 we set only angular filters of order 20 and 22 for $\theta$ and $\phi$ respectively. For the rotating stars we set an order 4 radial filter, an order 6 $\theta$ filter and an order 16 $\phi$ filter. The simulations are performed with a time step $\Delta t = 1.67 \times 10^{-7}$ [s], except for the simulations of the BU7 model where we use $\Delta t = 3.34 \times 10^{-8}$ [s]. We give the extracted frequencies of the modes in tables 2 and 3, that we can compare to tabulated frequencies of the literature, reported in table 4. We extract the fundamental quasi-radial $F$ mode, its first overtone $H_1$, the fundamental $\ell = 2$, $m = 0$ $2f_0$ mode, and the counter- and corotating $\ell = |m| = 2$ $2f_2$ and $2f_{-2}$ modes. We compute the relative differences between our computed values and the tabulated ones and report them in table 5. Some relative differences are of the order of 3%–4% and some surpass 10%, but the latter corresponds to some of the lowest frequencies (~100–200 Hz) for which we do not expect an accurate estimation with 25 ms simulations. Otherwise, most frequencies are recovered at the level of 1% or less. We therefore consider that the frequency extraction is successful. This success also justifies the CFC approach to extract oscillation frequencies, compared to full GR evolutions that are typically much more computationally expensive. The lower frequencies could be re-extracted with longer simulation using future improved versions of ROXAS. We also note that the frequencies of the $F$ and $H_1$ mode for the non-rotating model BU0 produced with the three-dimensional code (although extracted from two-dimensional simulations) are consistent with the frequencies produced with the early spherically symmetric version of the code [10]. We plot on figures 2 and 3 the GW signal from the simulation performed with the BU4 model, along with its Fourier spectrum, and highlight the peaks with their corresponding mode.
<figure><img src="image_5.png" /><figcaption>**Figure 2.** $A_+^e$ GW waveform for the BU4 model. The Fourier spectrum is shown on figure 3. The value of $A_+^e$ corresponds, for a source at 10 Mpc, to a strain $h_+ \sim 10^{-24}$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Fourier spectrum of the signal on figure 2, with the mode peaks highlighted.</figcaption></figure>
<table><caption>**Table 2.** Oscillation frequencies of the BU sequences, solving for the metric. The frequencies are given in kHz. The *F* mode is the fundamental quasi-radial mode, the *H*₁ mode is its first overtone, the ²*f*₀ mode is the fundamental *ℓ* = 2, *m* = 0 mode and the ²*f*₂ and ²*f*₂⁻ modes are the counter- and corotating *ℓ* = |*m*| = 2 modes. All frequencies are obtained with the perturbation (53) using *ε* = 10⁻². The displayed frequencies for the *F*, *H*₁ and ²*f*₀ mode are obtained with 2D simulations, but we checked that the 3D simulations gave consistent frequencies for these particular modes. δ*f*₂D,₃D = 1/tₘₐₓ are the discretization steps in the Fourier space, for 2D and 3D simulations. It is dictated by the length of the simulation, tₘₐₓ being the final simulation time. The ²*f*₂ mode of the BU7 star has not been extracted as its expected frequency is approximately 11 Hz i.e. a period of 91 ms, which the simulations do not resolve since they are shorter.</caption><thead><tr><th>Model</th><th><i>F</i></th><th><i>H</i><sub>1</sub></th><th><sup>2</sup><i>f</i><sub>0</sub></th><th>δ<i>f</i><sub>2D</sub> (Hz)</th><th><sup>2</sup><i>f</i><sub>2</sub></th><th><sup>2</sup><i>f</i><sub>-2</sub></th><th>δ<i>f</i><sub>3D</sub> (Hz)</th></tr></thead><tbody><tr><td>BU0</td><td>1.446</td><td>3.956</td><td>1.596</td><td>40.0</td><td>1.565</td><td>1.565</td><td>40.0</td></tr><tr><td>BU1</td><td>1.423</td><td>3.883</td><td>1.609</td><td>40.0</td><td>1.080</td><td>1.953</td><td>40.0</td></tr><tr><td>BU2</td><td>1.396</td><td>3.884</td><td>1.642</td><td>40.0</td><td>0.840</td><td>2.052</td><td>40.0</td></tr><tr><td>BU3</td><td>1.360</td><td>3.915</td><td>1.680</td><td>40.0</td><td>0.640</td><td>2.123</td><td>40.0</td></tr><tr><td>BU4</td><td>1.322</td><td>3.934</td><td>1.715</td><td>40.0</td><td>0.451</td><td>2.153</td><td>40.0</td></tr><tr><td>BU5</td><td>1.282</td><td>3.958</td><td>1.723</td><td>40.0</td><td>0.280</td><td>2.158</td><td>40.0</td></tr><tr><td>BU6</td><td>1.258</td><td>3.998</td><td>1.744</td><td>66.67</td><td>0.130</td><td>2.143</td><td>43.5</td></tr><tr><td>BU7</td><td>1.202</td><td>4.001</td><td>1.751</td><td>66.67</td><td>×</td><td>2.105</td><td>52.6</td></tr></tbody></table>

<table><caption>**Table 3.** Same as table 2 but in the Cowling approximation.</caption><thead><tr><th>Model</th><th><i>F</i></th><th><i>H</i><sub>1</sub></th><th><sup>2</sup><i>f</i><sub>0</sub></th><th>δ<i>f</i><sub>2D</sub> (Hz)</th><th><sup>2</sup><i>f</i><sub>2</sub></th><th><sup>2</sup><i>f</i><sub>-2</sub></th><th>δ<i>f</i><sub>3D</sub> (Hz)</th></tr></thead><tbody><tr><td>BU0</td><td>2.683</td><td>4.561</td><td>1.881</td><td>40.0</td><td>1.881</td><td>1.881</td><td>40.0</td></tr><tr><td>BU1</td><td>2.641</td><td>4.475</td><td>1.885</td><td>40.0</td><td>1.367</td><td>2.259</td><td>40.0</td></tr><tr><td>BU2</td><td>2.602</td><td>4.411</td><td>1.889</td><td>40.0</td><td>1.121</td><td>2.364</td><td>40.0</td></tr><tr><td>BU3</td><td>2.561</td><td>4.403</td><td>1.883</td><td>40.0</td><td>0.918</td><td>2.443</td><td>40.0</td></tr><tr><td>BU4</td><td>2.522</td><td>4.402</td><td>1.877</td><td>40.0</td><td>0.721</td><td>2.473</td><td>40.0</td></tr><tr><td>BU5</td><td>2.484</td><td>4.403</td><td>1.837</td><td>40.0</td><td>0.520</td><td>2.480</td><td>40.0</td></tr><tr><td>BU6</td><td>2.461</td><td>4.401</td><td>1.762</td><td>40.0</td><td>0.348</td><td>2.477</td><td>40.0</td></tr><tr><td>BU7</td><td>2.439</td><td>4.393</td><td>1.681</td><td>40.0</td><td>0.163</td><td>2.396</td><td>40.0</td></tr></tbody></table>

Furthermore, we performed the simulation of a rotating star described with the barotropic SLy4 EoS [49]. To do so, we used the pseudo-polytropic fit described in [26] and for which the coefficients therein. We choose a rotation frequency of *f* = 191 Hz with almost the same central enthalpy as the BU sequence: *H*<sub>c</sub> = 0.22656, so that it has a baryon mass of *M*<sub>b</sub> = 1.504*M<sub>☉</sub>*, a gravitational mass of *M* = 1.361*M<sub>☉</sub>* and Ω/Ω<sub>K</sub> ≈ 16%. Since its equilibrium coordinate polar radius is *r*<sub>p</sub> = 9.34 km, we set the nucleus of the metric grid up to 7 km instead of 8 km, and set *N*<sub>r</sub> = 33 in all metric domains (except the CED where we keep *N*<sub>r</sub> = 17). The angular filter for the hydro grid are set at order 10 for the *θ* direction and 12 for the *φ* direction, and the filters for the metric grid are set to order 4 (radial), order 6 (*θ*) and order 16 (*φ*). We use equation (53) with *ε* = 10⁻⁶ to perturb the star, and let it evolve for 25 ms. We extract the frequencies of the ²*f*₂ and ²*f*₂⁻ modes for which we find frequencies at 1.675 kHz and 2.185 kHz. We compare with [48] where they study the same star, and they find 1.655 kHz and 2.105 kHz respectively, which constitutes a relative difference of 1.2% and 3.7%. Once again the frequency extraction is satisfactory and this is an important step for the validation of the ROXAS code with realistic EoSs. We also provide the frequencies of the other modes: the *F* and ²*f*₀ are respectively found at frequencies of 3.085 kHz and 1.964 kHz. We finally provide
<table><caption>**Table 4.** Frequencies extracted from table II of [45] for axisymmetric modes with the Cowling approximation, table 5 of [15] for axisymmetric modes in CFC, digitized with engauge-digitizer [50] from figure 1 of [46] for non-axisymmetric modes with the Cowling approximation and table VII of [48] for non-axisymmetric modes in full GR.</caption><thead><tr><th rowspan="2">Model</th><th colspan="2">F</th><th colspan="2">H<sub>1</sub></th><th colspan="2"><sup>2</sup>f<sub>0</sub></th><th colspan="2"><sup>2</sup>f<sub>2</sub></th><th colspan="2"><sup>2</sup>f<sub>-2</sub></th></tr><tr><th>Cowling</th><th>CFC</th><th>Cowling</th><th>CFC</th><th>Cowling</th><th>CFC</th><th>Cowling</th><th>Full GR</th><th>Cowling</th><th>Full GR</th></tr></thead><tbody><tr><td>BU0</td><td>2.679</td><td>1.458</td><td>4.561</td><td>3.971</td><td>1.890</td><td>1.586</td><td>1.881</td><td>1.578</td><td>1.881</td><td>1.578</td></tr><tr><td>BU1</td><td>2.638</td><td>1.413</td><td>4.466</td><td>3.915</td><td>1.890</td><td>1.611</td><td>1.374</td><td>1.079</td><td>2.257</td><td>1.942</td></tr><tr><td>BU2</td><td>2.605</td><td>1.380</td><td>4.435</td><td>3.907</td><td>1.906</td><td>1.635</td><td>1.122</td><td>0.836</td><td>2.362</td><td>2.054</td></tr><tr><td>BU3</td><td>2.570</td><td>1.343</td><td>4.409</td><td>3.921</td><td>1.895</td><td>1.669</td><td>0.909</td><td>0.636</td><td>2.428</td><td>2.118</td></tr><tr><td>BU4</td><td>2.539</td><td>1.304</td><td>4.410</td><td>3.950</td><td>1.875</td><td>1.698</td><td>0.709</td><td>0.456</td><td>2.459</td><td>2.155</td></tr><tr><td>BU5</td><td>2.500</td><td>1.281</td><td>4.400</td><td>3.964</td><td>1.844</td><td>1.714</td><td>0.534</td><td>0.293</td><td>2.469</td><td>2.127</td></tr><tr><td>BU6</td><td>2.484</td><td>1.219</td><td>4.392</td><td>4.010</td><td>1.794</td><td>1.729</td><td>0.355</td><td>0.146</td><td>2.467</td><td>2.182</td></tr><tr><td>BU7</td><td>2.456</td><td>1.207</td><td>4.394</td><td>4.018</td><td>1.703</td><td>1.720</td><td>0.181</td><td>0.011</td><td>2.441</td><td>2.188</td></tr></tbody></table>
<table><caption>**Table 5.** Relative difference in absolute value between the frequencies of tables 2 and 3 computed with our code and the tabulated frequencies reported in table 4, in %.</caption><thead><tr><th rowspan="2">Model</th><th colspan="2">F</th><th colspan="2">H<sub>1</sub></th><th colspan="2">2<sup>f<sub>0</sub></sup></th><th colspan="2">2<sup>f<sub>2</sub></sup></th><th colspan="2">2<sup>f<sub>-2</sub></sup></th></tr><tr><th>Cowling</th><th>CFC</th><th>Cowling</th><th>CFC</th><th>Cowling</th><th>CFC</th><th>Cowling</th><th>Full GR</th><th>Cowling</th><th>Full GR</th></tr></thead><tbody><tr><td>BU0</td><td>0.15</td><td>0.83</td><td>0.00</td><td>0.38</td><td>0.48</td><td>0.63</td><td>0.00</td><td>0.83</td><td>0.00</td><td>0.83</td></tr><tr><td>BU1</td><td>0.11</td><td>0.70</td><td>0.20</td><td>0.82</td><td>0.27</td><td>0.12</td><td>0.51</td><td>0.09</td><td>0.09</td><td>0.56</td></tr><tr><td>BU2</td><td>0.12</td><td>1.15</td><td>0.54</td><td>0.59</td><td>0.90</td><td>0.43</td><td>0.98</td><td>0.48</td><td>0.09</td><td>0.10</td></tr><tr><td>BU3</td><td>0.35</td><td>1.25</td><td>0.14</td><td>0.15</td><td>0.63</td><td>0.65</td><td>0.09</td><td>0.63</td><td>0.61</td><td>0.24</td></tr><tr><td>BU4</td><td>0.67</td><td>1.36</td><td>0.20</td><td>0.40</td><td>0.64</td><td>0.99</td><td>1.66</td><td>1.11</td><td>0.57</td><td>0.09</td></tr><tr><td>BU5</td><td>0.64</td><td>0.08</td><td>0.07</td><td>0.15</td><td>0.38</td><td>0.52</td><td>2.69</td><td>4.64</td><td>0.44</td><td>1.44</td></tr><tr><td>BU6</td><td>0.93</td><td>3.10</td><td>0.20</td><td>0.30</td><td>1.82</td><td>0.86</td><td>2.01</td><td>12.3</td><td>0.40</td><td>1.82</td></tr><tr><td>BU7</td><td>0.70</td><td>0.42</td><td>0.02</td><td>0.57</td><td>1.31</td><td>1.77</td><td>11.0</td><td>×</td><td>1.88</td><td>3.94</td></tr></tbody></table>
the values of the modes with simulations performed using the Cowling approximation: the frequency extraction gives 4.415kHz, 2.400kHz, 2.141kHz and 2.640kHz for the $F, ^2f_0, ^2f_2$ and $^2f_{-2}$ modes respectively.

**5. Summary and conclusion**

We presented the main numerical ingredients that compose the new code ROXAS, a general relativistic spectral code dedicated to NS oscillations. Using LORENE as its foundations, we use the xCFC formulation of Einstein equations to compute the metric, and the hydrodynamics sector is based upon the previously published formulation that uses only primitive variables, which is already tried and tested in spherical symmetry. To be able to perform simulations of rotating stars perturbed with non-axisymmetric profiles we implemented and developed a number of techniques which allowed to perform simulations with spherical coordinates, including the origin in the grid: regularization schemes, a multidomain approach, filters, an adapted mapping that allows to follow the surface of the star throughout the simulations, two distinct grids that respectively support metric computation and hydrodynamics, along with interpolation schemes, and the use of a well-balanced scheme. Every tool or technique has been thoroughly tested and therefore allowed to perform simulations of polytropic stars and stars described with an analytical representation of a realistic EoS. The standard test beds for NS general relativistic evolution codes such as the migration test and the collapse to a black hole were previously presented in [10], and we here presented the frequency extraction for axisymmetric and non-axisymmetric simulations. The results were compared to previously published results and were found to be very accurate. In particular, we emphasize that although the simulations are performed under the CFC approximation, the reference frequencies for non-axisymmetric simulations were computed in full General Relativity, and the frequency extraction is very satisfactory. This means that the CFC approximation seems to be very good for frequency extraction; using it represents a significant gain of computational resources compared to a possible implementation of General Relativity. We also emphasize that the imposed symmetries still allow to extract the dominant quadrupolar mode $ℓ = |m| = 2$. All in all, ROXAS fulfills its goal to be a lightweight code that runs on office computers and provides the possibility to perform parametric studies of NS oscillations using a wide panel of EoS. Future versions will also support EoS with more than one argument by taking out-of-β-equilibrium and temperature effects into account, relax the spatial symmetries, go beyond CFC and support differential rotation profiles, for example those of [51] that are relevant for simulations of a post-coalescence object like a HMNS.

**Data availability statement**

The data that support the findings of this study are openly available at the following URL/DOI: https://doi.org/10.5281/zenodo.14849547.

**Acknowledgments**

The authors gratefully acknowledge the Agence Nationale de la Recherche (ANR) under Contract ANR-22-CE31-0001-01; the Italian Instituto Nazionale de Fisica Nucleare (INFN), the French Centre National de la Recherche Scientifique (CNRS) and the Netherlands
Organization for Scientific Research for the construction and operation of the Virgo detector and the creation and support of the EGO consortium. GS thanks all the ROC team of LUTH for creating a work environment conducive to the success of his PhD. The authors acknowledge the precious help of Santiago Jaraba for correcting issues in the implementation of the quadrupole formula.

## Appendix A. Derivation of the impermeable BC

In this appendix we show the derivation of the expression of the coordinates of a normal vector defined on an arbitrarily deformed sphere. We then apply this to write the impermeable BC with a given coordinate velocity field ($v^r, v^\theta, v^\varphi$). To find out how to write the condition, we first assume that the deformed boundary is mapped from $S^2$, the unit sphere of $\mathbb{R}^3$ by:

$$ \mathbf{R}: \mathbf{x} \in S^2 \mapsto R(\theta, \varphi) \mathbf{x}. \qquad (A.1) $$

Then we compute the (rescaled) angular gradient of the radius:

$$ \mathbf{g} = \frac{\nabla_{\theta\varphi} R}{R}, \qquad (A.2) $$

where

$$ \nabla_{\theta\varphi} = \begin{pmatrix} 0 \\ \frac{\partial}{\partial\theta} \\ \frac{1}{\sin\theta} \frac{\partial}{\partial\varphi} \end{pmatrix}. \qquad (A.3) $$

Then we compute the tangent vector:

$$ \mathbf{h} = \mathbf{g} - (\mathbf{g} \cdot \mathbf{x}) \mathbf{x}. \qquad (A.4) $$

Finally, the (unnormalised) normal vector to the deformed sphere at $\mathbf{R}(\mathbf{x})$ is:

$$ \mathbf{n} = \mathbf{x} - \mathbf{h}. \qquad (A.5) $$

Its coordinates in the orthonormal spherical triad are:

$$ n_r = x_r - h_r = 1 - (g_r - g_r) = 1, \qquad (A.6) $$

$$ n_{\theta} = x_{\theta} - h_{\theta} = -g_{\theta} = -\frac{1}{R} \frac{\partial R}{\partial \theta}, \qquad (A.7) $$

$$ n_{\varphi} = x_{\varphi} - h_{\varphi} = -g_{\varphi} = -\frac{1}{R \sin\theta} \frac{\partial R}{\partial \varphi}. \qquad (A.8) $$

Therefore, with $v_{\text{border}} = \frac{\partial R}{\partial t} \hat{e}_r$, the impermeable BC is expressed as:

$$ v_r(R(t), t) - \frac{v_{\theta}(R(t), t)}{R} \frac{\partial R}{\partial \theta} - \frac{v_{\varphi}(R(t), t)}{R \sin\theta} \frac{\partial R}{\partial \varphi} = \frac{\partial R}{\partial t}. \qquad (A.9) $$

This expression was found by [52] when looking for non-axisymmetric equilibrium configurations of self-gravitating objects. We challenged the expression numerically, the results can be found in appendix B.
# Appendix B. Numerical test of the impermeable BC

In this appendix we provide numerical evidence that equation (A.9) corresponds to the correct impermeable BC for a general velocity field inside a self-gravitating object. To do so, we perform the simulation of a rotating ellipsoid and show a comparison with the theoretical rotating configuration, which can be derived analytically in some simple cases. We start with some definitions, write equation (A.9) for a rotating ellipsoid around a given axis, derive the analytical formula with which we compare the code, and show some numerical results.

## B.1. Introduction

**B.1.1. Definition** An ellipsoid is a surface of $\mathbb{R}^3$ that can be obtained with an affine transformation applied on a sphere. Its general definition in Cartesian coordinates is:

$$ \frac{x^2}{a^2} + \frac{y^2}{b^2} + \frac{z^2}{c^2} = 1, \qquad (B.1) $$

where a, b, c are the length of the three principal axes.

**B.1.2. Parametrization** Using polar spherical coordinates $(r, \theta, \varphi)$, the surface of the ellipsoid can be parametrized as:

$$ R(\theta, \varphi) = \left( \frac{\cos^2 \varphi \sin^2 \theta}{a^2} + \frac{\sin^2 \varphi \sin^2 \theta}{b^2} + \frac{\cos^2 \theta}{c^2} \right)^{-1/2}, \qquad (B.2) $$

or, alternatively:

$$ R(\theta, \varphi) = \frac{abc}{\sqrt{c^2 \sin^2 \theta (b^2 \cos^2 \varphi + a^2 \sin^2 \varphi) + a^2 b^2 \cos^2 \theta}}. \qquad (B.3) $$

## B.2. Rotation field

Now that we have given the parametrization of an ellipsoid, let us turn to the velocity field that describes a rotation around a given axis. We define it as

$$ \mathbf{v} = R(\theta, \varphi) \mathbf{\Omega}, \qquad (B.4) $$

where $\mathbf{\Omega}$ is the rotation vector, of which the norm $\mathbf{\Omega} = |\mathbf{\Omega}|$ is the angular rotation rate and direction is defined with respect to some point of the ellipsoid $\mathbf{M} = R\hat{\mathbf{e}}_r$ and the chosen rotation axis. It means that $\mathbf{v} = \mathbf{\Omega}\mathbf{M} \times \hat{\mathbf{e}}_r$, where $\hat{\mathbf{e}}_r$ is a generic axis in $\mathbb{R}^3$. This yields three possible rotation laws associated to the three Cartesian axes (of which the coordinates are given in the spherical orthonormal basis):

$$ \mathbf{\Omega}_x / \mathbf{\Omega} = \hat{\mathbf{e}}_x \times \hat{\mathbf{e}}_r = \begin{pmatrix} 0 \\ -\sin\varphi \\ -\cos\theta\cos\varphi \end{pmatrix}, \qquad (B.5) $$

$$ \mathbf{\Omega}_y / \mathbf{\Omega} = \hat{\mathbf{e}}_y \times \hat{\mathbf{e}}_r = \begin{pmatrix} 0 \\ \cos\varphi \\ -\cos\theta\sin\varphi \end{pmatrix}, \qquad (B.6) $$

$$ \mathbf{\Omega}_z / \mathbf{\Omega} = \hat{\mathbf{e}}_z \times \hat{\mathbf{e}}_r = \begin{pmatrix} 0 \\ 0 \\ \sin\theta \end{pmatrix}, \qquad (B.7) $$
<figure><img src="image_7.png" /><figcaption>Figure B1. Convergence plots for a rotating ellipsoid. The quantity plotted is defined by equation (B.12). On the left, we fixed $t_{max} = 100$, $\Delta t = 10^{-3}$ and $N_{\varphi} = 32$ and varied $N_{\theta}$. On the right, we fixed $t_{max} = 10$, $N_{\theta} = 65$ and $N_{\varphi} = 64$ and varied $\Delta t$.</figcaption></figure>

and therefore three PDEs that describe the evolution of a rotating ellipsoid:

$$ \frac{1}{\Omega} \frac{\partial R}{\partial t} = \sin \varphi \frac{\partial R}{\partial \theta} + \frac{\cos \varphi}{\tan \theta} \frac{\partial R}{\partial \varphi} \quad (\text{rotation around x axis}), \quad (B.8) $$

$$ \frac{1}{\Omega} \frac{\partial R}{\partial t} = -\cos \varphi \frac{\partial R}{\partial \theta} + \frac{\sin \varphi}{\tan \theta} \frac{\partial R}{\partial \varphi} \quad (\text{rotation around y axis}), \quad (B.9) $$

$$ \frac{1}{\Omega} \frac{\partial R}{\partial t} = - \frac{\partial R}{\partial \varphi} \quad (\text{rotation around z axis}). \quad (B.10) $$

For a given constant angular rotation rate $\Omega$, the quantity $\alpha$, defined as the angle between the ellipsoid and a given axis and assumed to be initially zero, varies as

$$ \alpha(t) = \Omega t. \quad (B.11) $$

With this convention, the ellipsoid rotates counterclockwise with respect to the chosen axis.

### B.3. Numerical results

For all simulations we set $\Omega = 1$. The angular derivatives are computed thanks to the spectral representation of $R(\theta, \varphi)$, and the evolution is performed thanks to an RK4 scheme. Rotation around the z axis is straightforward, however when the ellipsoid rotates around the x or the y axis, it crosses the z axis on which the angular coordinates are singular. To cure this issue, we perform the regularization process described in section 3.3.2. We compute the rotation around the y axis for an ellipsoid with the parameters $a = 2$, $b = 1.5$, $c = 1$. On figure B1 we show the results of a convergence test. We measure the error on a simulation by computing

$$ \delta R/R = \max_{t \in [0, t_{\text{max}}], \theta \in [0, \pi], \varphi \in [0, 2\pi]} \left| \frac{R_{\text{numerical}}(t, \theta, \varphi) - R_{\text{analytical}}(t, \theta, \varphi)}{R_{\text{analytical}}(t, \theta, \varphi)} \right|, \quad (B.12) $$

where $R_{\text{numerical}}$ is the ellipsoid configuration computed with the code and $R_{\text{analytical}}$ the one computed by writing the general expression of a tilted ellipsoid around a given axis, and assuming the rotation law (B.11).
By fixing $\Delta t$ and $N_\phi$ to vary $N_\theta$ we show the exponential convergence of the spatial scheme, which is here shown by the presence of a straight line in the semilog plot before reaching the plateau at $N_\theta = 33$, and then fixing $N_\theta$ and $N_\phi$ to vary $\Delta t$ we show the order of the time scheme: since we use the RK4 scheme we expect a $O(\Delta t^4)$ behaviour; the log-log plot shows a straight line with a slope of 4, which validates the convergence test for the time integration scheme.

## Appendix C. Scalar field matching

In this appendix we describe the so-called tau method to match a scalar field $f$ between a nucleus and a single shell by imposing the continuity of $f$ and $\partial_{\nu}f$. The spectral expansion of $f$ is:

$$f_{\text{nuc}}(\xi) = \sum_{k=0}^{N_{r,\text{nuc}}-1} f_{\text{nuc},k} T_{2k}(\xi) \qquad (C.1)$$

in the nucleus, and:

$$f_{\text{sh}}(\xi) = \sum_{k=0}^{N_{r,\text{sh}}-1} f_{\text{sh},k} T_k(\xi) \qquad (C.2)$$

in the shell. Writing the continuity of the field and of its first derivative allows writing a linear system with $f_{n,N_{r,\text{nuc}}-1}$ and $f_{s,N_{r,\text{sh}}-1}$ as unknowns by writing the following conditions:

$$f_{\text{nuc}}(1) = f_{\text{sh}}(-1), \qquad (C.3)$$

$$\frac{1}{\alpha_{\text{nuc}}} f'_{\text{nuc}}(1) = \frac{1}{\alpha_{\text{sh}}} f'_{\text{sh}}(-1), \qquad (C.4)$$

and making use of the following properties of the Chebyshev polynomials [41]:

$$\forall k \in \mathbb{N}, T_k(\pm 1) = (\pm 1)^k, T'_k(\pm 1) = (\pm 1)^{k+1} k^2, \qquad (C.5)$$

which allow to rewrite the junction conditions:

$$\sum_{k=0}^{N_{r,\text{nuc}}-1} f_{\text{nuc},k} = \sum_{k=0}^{N_{r,\text{sh}}-1} (-1)^k f_{\text{sh},k}, \qquad (C.6)$$

$$\frac{4}{\alpha_{\text{nuc}}} \sum_{k=0}^{N_{r,\text{nuc}}-1} k^2 f_{\text{nuc},k} = \frac{1}{\alpha_{\text{sh}}} \sum_{k=0}^{N_{r,\text{sh}}-1} (-1)^{k+1} k^2 f_{\text{sh},k}. \qquad (C.7)$$

We can rewrite this in a matrix form $AX = B$ with:

$$A = \begin{pmatrix} 1 & -1 \\ \frac{4(N_{r,\text{nuc}} - 1)^2}{\alpha_{\text{nuc}}} & \frac{(N_{r,\text{sh}} - 1)^2}{\alpha_{\text{sh}}} \end{pmatrix}, \qquad (C.8)$$

$$X = \begin{pmatrix} f_{n,N_{r,\text{nuc}}-1} \\ f_{s,N_{r,\text{sh}}-1} \end{pmatrix}, \qquad (C.9)$$

$$B = \begin{pmatrix} \sum_{k=0}^{N_{r,\text{nuc}}-2} f_{\text{nuc},k} - \sum_{k=0}^{N_{r,\text{sh}}-2} (-1)^k f_{\text{sh},k} \\ \frac{4}{\alpha_{\text{nuc}}} \sum_{k=0}^{N_{r,\text{nuc}}-2} k^2 f_{\text{nuc},k} - \frac{1}{\alpha_{\text{sh}}} \sum_{k=0}^{N_{r,\text{sh}}-2} (-1)^{k+1} k^2 f_{\text{sh},k} \end{pmatrix}, \qquad (C.10)$$
then solve for X by computing the inverse of A, and finally replace the last coefficients of $f_{nuc}$ and $f_{sh}$ by the first and second coefficients of the vector X. We make use of the explicit inverse of a 2×2 matrix to write:

$$A^{-1} = \left( \begin{array}{cc} \frac{\alpha_{nuc} (N_{r,sh} - 1)^2}{\alpha_{nuc} (N_{r,sh} - 1)^2 + 4\alpha_{sh} (N_{r,nuc} - 1)^2} & -\frac{\alpha_{nuc}\alpha_{sh}}{\alpha_{nuc} (N_{r,sh} - 1)^2 + 4\alpha_{sh} (N_{r,nuc} - 1)^2} \\ \frac{4}{\alpha_{nuc} (N_{r,sh} - 1)^2 + 4\alpha_{sh} (N_{r,nuc} - 1)^2} & -\frac{\alpha_{nuc}\alpha_{sh}}{\alpha_{nuc} (N_{r,sh} - 1)^2 + 4\alpha_{sh} (N_{r,nuc} - 1)^2} \end{array} \right) .(C.11)$$

## Appendix D. Relative precision on the frequency extraction

<table><caption>**Table D1.** This table gives the total time (rounded) needed to capture a hundred periods of every BU0 to BU7 model. The values of the $^2f_2$ modes are taken from [48]. The period is the inverse of the frequency and it is then multiplied by 100 to give the total simulation length that captures as many periods.</caption><thead><tr><th>Model</th><th><sup>2</sup>f<sub>2</sub> (kHz)</th><th>Period (ms)</th><th>Needed simulation length [ms]</th></tr></thead><tbody><tr><td>BU0</td><td>1.578</td><td>6.337 &times; 10<sup>-1</sup></td><td>63</td></tr><tr><td>BU1</td><td>1.079</td><td>9.268 &times; 10<sup>-1</sup></td><td>93</td></tr><tr><td>BU2</td><td>0.836</td><td>1.196</td><td>120</td></tr><tr><td>BU3</td><td>0.636</td><td>1.572</td><td>157</td></tr><tr><td>BU4</td><td>0.456</td><td>2.193</td><td>219</td></tr><tr><td>BU5</td><td>0.293</td><td>3.413</td><td>341</td></tr><tr><td>BU6</td><td>0.146</td><td>6.849</td><td>685</td></tr><tr><td>BU7</td><td>0.011</td><td>91</td><td>9100</td></tr></tbody></table>

In this appendix we provide a table that shows, for every BU model that we evolved, the theoretical simulation length necessary to guarantee a relative difference of 1% on the extraction of the $^2f_2$ mode (i.e. obtaining at least the first two digits), as given by the Fourier sampling that only depends on the simulation length: $\Delta f = 1/t_{max}$. We use the tabulated values of [48]. This time corresponds to a hundred periods captured by the simulation.

We see on the values of table 5 that the code does not need to capture a hundred periods to achieve 1% accuracy in the extraction thanks to the interpolation of the Fourier spectrum that we use to extract the frequencies from the signal. However, table D1 gives a good order of magnitude of the simulation length needed to achieve good accuracy.

## ORCID iDs

Gaël Servignat https://orcid.org/0000-0003-0057-922X

Jérôme Novak https://orcid.org/0000-0002-6029-4712

## References

[1] Abbott B P et al (LIGO Scientific Collaboration and Virgo Collaboration) 2016 GW150914: the advanced LIGO detectors in the era of first discoveries Phys. Rev. Lett. 116 131103

[2] Abbott R et al 2023 GWTC-3: compact binary coalescences observed by ligo and virgo during the second part of the third observing run Phys. Rev. X 13 041039
[3] Punturo M et al 2010 The Einstein telescope: a third-generation gravitational wave observatory
*Class. Quantum Grav.* **27** 194002

[4] Maggiore M et al 2020 Science case for the Einstein telescope *J. Cosmol. Astropart. Phys.* JCAP03(2020)050

[5] Branchesi M et al 2023 Science with the Einstein telescope: a comparison of different designs *J. Cosmol. Astropart. Phys.* JCAP07(2023)068

[6] Reitze D et al 2019 Cosmic explorer: the U.S. contribution to gravitational-wave astronomy beyond LIGO vol 51 p 35 (arXiv:1907.04833)

[7] Evans M et al 2021 A horizon study for cosmic explorer: science, observatories, and community (arXiv:2109.09882)

[8] Metzger B D 2020 Kilonovae *Living Rev. Relativ.* **23** 1

[9] Servignat G, Novak J and Jaraba S 2025 ROXAS (1.0) Zenodo (https://doi.org/10.5281/zenodo.14849547)

[10] Servignat G, Novak J and Cordero-Carrión I 2023 A new formulation of general-relativistic hydrodynamic equations using primitive variables *Class. Quantum Grav.* **40** 105002

[11] Cordero-Carrión I, Cerdá-Durán P, Dimmelmeier H, Luis Jaramillo J, Novak J and Gourgoulhon E 2009 Improved constrained scheme for the Einstein equations: an approach to the uniqueness issue *Phys. Rev. D* **79** 024017

[12] Cipolletta F, Kalinani J V, Giangrandi E, Giacomazzo B, Ciolfi R, Sala L and Giudici B 2021 Spritz: general relativistic magnetohydrodynamics with neutrinos *Class. Quantum Grav.* **38** 085021

[13] Cipolletta F, Kalinani J V, Giacomazzo B and Ciolfi R 2022 Zenodo (https://doi.org/10.5281/zenodo.4350072)

[14] Dimmelmeier H, Novak J, Font J A, Ibáñez J M and Müller E 2005 Combining spectral and shock-capturing methods: a new numerical approach for 3D relativistic core collapse simulations *Phys. Rev. D* **71** 064023

[15] Dimmelmeier H, Stergioulas N and Font J A 2006 Non-linear axisymmetric pulsations of rotating relativistic stars in the conformal flatness approximation *Mon. Not. R. Astron. Soc.* **368** 1609–30

[16] Thierfelder M, Bernuzzi S and Brügmann B 2011 Numerical relativity simulations of binary neutron stars *Phys. Rev. D* **84** 044012

[17] Pakmor R, Springel V, Bauer A, Mocz P, Munoz D J, Ohlmann S T, Schaal K and Zhu C 2016 Improving the convergence properties of the moving-mesh code AREPO *Mon. Not. R. Astron. Soc.* **455** 1134–43

[18] Lioutas G, Bauswein A, Soultanis T, Pakmor R, Springel V and Röpke F K 2024 General relativistic moving-mesh hydrodynamic simulations with arepo and applications to neutron star mergers *Mon. Not. R. Astron. Soc.* **528** 1906–29

[19] Kidder L E et al 2017 SpECTRE: a task-based discontinuous Galerkin code for relativistic astrophysics *J. Comput. Phys.* **335** 84–114

[20] Rosswog S and Diener P 2021 SPHINCS_bssn: a general relativistic smooth particle hydrodynamics code for dynamical spacetimes *Class. Quantum Grav.* **38** 115002

[21] Yamamoto T, Shibata M and Taniguchi K 2008 Simulating coalescing compact binaries by a new code (SACRA) *Phys. Rev. D* **78** 064054

[22] Banyuls F, Font J A, Ibanez J M, Marti J M and Miralles J A 1997 Numerical 3 + 1 General relativistic hydrodynamics: a local characteristic approach *Astrophys. J.* **476** 221–31

[23] Kokkotas K D and Ruoff J 2001 Radial oscillations of relativistic stars *Astron. Astrophys.* **366** 565–72

[24] Bernuzzi S and Nagar A 2008 Gravitational waves from pulsations of neutron stars described by realistic equations of state *Phys. Rev. D* **78** 024024

[25] Baiotti L, Bernuzzi S, Corvino G, De Pietri R and Nagar A 2009 Gravitational-wave extraction from neutron-star oscillations: comparing linear and nonlinear techniques *Phys. Rev. D* **79** 024002

[26] Servignat G, Davis P J, Novak J, Oertel M and Pons J A 2024 One- and two-argument equation of state parametrizations with continuous sound speed for neutron star simulations *Phys. Rev. D* **109** 103022

[27] Bonazzola S, Gourgoulhon E, Grandclément P and Novak J 2004 Constrained scheme for the Einstein equations based on the Dirac gauge and spherical coordinates *Phys. Rev. D* **70** 104007

[28] Bonazzola S, Gourgoulhon E, Salgado M and Marck J A 1993 Axisymmetric rotating relativistic bodies: a new numerical approach for ‘exact’ solutions *Astron. Astrophys.* **278** 421–43
[29] Dumbser M, Zanotti O, Gaburro E and Peshkov I 2023 A well-balanced discontinuous Galerkin method for the first-order Z4 formulation of the Einstein–Euler system *J. Comput. Phys.* **504** 112875

[30] Einstein A 1918 Über Gravitationswellen *Sitzungsberichte der Königlich Preussischen Akademie der Wissenschaften* 154–67

[31] Blanchet L, Damour T and Schaefer G 1990 Post-Newtonian hydrodynamics and post-Newtonian gravitational wave generation for numerical relativity *Mon. Not. R. Astron. Soc.* **242** 289–305

[32] Misner C W, Thorne K S and Archibald Wheeler J 1973 *Gravitation* (W. H. Freeman)

[33] Gourgoulhon E, Grandclément P, Marck J-A, Novak J and Taniguchi K 2016 LORENE: Spectral methods differential equations solver (Astrophysics Source Code Library)

[34] Grandclément P, Bonazzola S, Gourgoulhon E and Marck J A 2001 A multidomain spectral method for scalar and vectorial Poisson equations with noncompact sources *J. Comput. Phys.* **170** 231–60

[35] Gottlieb D and Orszag S A 1989 *Numerical Analysis of Spectral Methods: Theory and Applications (CBMS-NSF Regional Conference Series in Applied Mathematics. Society for Industrial and Applied Mathematics)*

[36] Boyd J P 2001 *Chebyshev and Fourier Spectral Methods* 2nd ed., rev edn (Dover Publications)

[37] Grandclément P and Novak J 2009 Spectral methods for numerical relativity *Living Rev. Relativ.* **12** 1

[38] Gómez R, Barreto W and Frittelli S 2007 Framework for large-scale relativistic simulations in the characteristic approach *Phys. Rev. D* **76** 124029

[39] Nicoules J 2023 Numerical evolutions in general relativity: application to the schwarzschild black hole and teukolsky gravitational waves *PhD Thesis* Université Paris Cité

[40] Canuto C, Hussaini M Y, Quarteroni A and Zang T A 2006 *Spectral Methods* (Springer)

[41] Hesthaven J S, Gottlieb S and Gottlieb D 2007 *Spectral Methods for Time-Dependent Problems* (Cambridge University Press)

[42] Gourgoulhon E, Haensel P, Livine R, Paluch E, Bonazzola S and Marck J A 1999 Fast rotation of strange stars *Astron. Astrophys.* **349** 851–62

[43] Lin L-M and Novak J 2006 Rotating star initial data for a constrained scheme in numerical relativity *Class. Quantum Grav.* **23** 4545–61

[44] Stergioulas N, Apostolatos T A and Font J A 2004 Non-linear pulsations in differentially rotating neutron stars: mass-shedding-induced damping and splitting of the fundamental mode *Mon. Not. R. Astron. Soc.* **352** 1089–101

[45] Gaertig E and Kokkotas K D 2008 Oscillations of rapidly rotating relativistic stars *Phys. Rev. D* **78** 064063

[46] Krüger C, Gaertig E and Kokkotas K D 2010 Oscillations and instabilities of fast and differentially rotating relativistic stars *Phys. Rev. D* **81** 084019

[47] Zink B, Korobkin O, Schnetter E and Stergioulas N 2010 Frequency band of the f-mode Chandrasekhar-Friedman-Schutz instability *Phys. Rev. D* **81** 084055

[48] Krüger C J and Kokkotas K D 2020 Dynamics of fast rotating neutron stars: an approach in the Hilbert gauge *Phys. Rev. D* **102** 064026

[49] Gulminelli F and Raduta A R 2015 Unified treatment of subsaturation stellar matter at zero and finite temperature *Phys. Rev. C* **92** 055803

[50] Mitchell M, Muftakhidinov B, Winchen T, Wilms A and Van Schaik B, Badshah400, Mo-Gul, The Gitter Badger, Zbigniew Jedrzejewski-Szmek, Kensington, and Kylesower 2020 Zenodo markumitchell/engauge-digitizer (https://doi.org/10.5281/zenodo.3941227)

[51] Iosif P and Stergioulas N 2021 Equilibrium sequences of differentially rotating stars with post-merger-like rotational profiles *Mon. Not. R. Astron. Soc.* **503** 850–66

[52] Uryu K and Eriguchi Y 1996 Existence of non-axisymmetric polytropes sustained by internal motions *Mon. Not. R. Astron. Soc.* **282** 653–64