© The Authors 2026

1 Introduction

The direct detection of gravitational waves (GWs), which started with the binary black hole merger GW150914 (Abbott et al. 2016), has opened a new window onto the most extreme astrophysical phenomena. Subsequently, the first GW signal from a binary neutron star (NS) merger, GW170817 (Abbott et al. 2017b), was observed together with its electromagnetic counterpart (Abbott et al. 2017a,c), which provided unprecedented insight into the dynamics of this system (Abbott et al. 2018). These detections were made possible by the current-generation GW detectors: the Laser Interferometer Gravitational-wave Observatory (LIGO) (Aasi et al. 2015), the Virgo GW detector (Acernese et al. 2015), and the Kamioka Gravitational-wave detector (KAGRA) (Akutsu et al. 2021).

From the first numerical simulations of the merger of two NSs (Shibata & Uryū 2000; Shibata et al. 2005), it is thought that in the majority of cases the post-merger remnant of a binary NS is a hypermassive neutron star (HMNS; Baumgarte et al. 2000), which is supported by its differential rotation for a few tenths of seconds. After that, the rotation slows down and becomes rigid, making this object collapse to a black hole (Metzger 2020). During this period, the HMNS is expected to emit GWs in the kilohertz frequency band, in a very similar manner to excited differentially rotating isolated NSs (Stergioulas et al. 2011). Thus, the study of the properties of differentially rotating NSs (see, e.g., recent works by Weih et al. 2018; Muhammed et al. 2024; Szewczyk et al. 2025) is very relevant to our understanding of binary NS mergers (Bauswein et al. 2016).

Note that current GW detectors lack sufficient sensitivity in the kilohertz band to detect the GWs from HMNS oscillations (Abbott et al. 2020). However, next-generation detectors such as the Einstein Telescope (Punturo et al. 2010), Cosmic Explorer (Reitze et al. 2019), and Neutron Star Extreme Matter Observatory (Ackley et al. 2020) will have a significantly increased sensitivity in this frequency range, where future detections could be a probe of the NS internal structure and its equation of state (EoS). In order to maximize the scientific output from these signals by the time they arrive, it is thus necessary to have a solid understanding of how the EoS, rotation, and other physical parameters shape the emitted GW.

To this end, the ROXAS code (Relativistic Oscillations of non-aXisymmetric neutron stArS; Servignat & Novak 2025; Servignat et al. 2025) was recently developed to evolve isolated NSs. It is a full nonlinear general-relativistic hydrodynamics code, primarily aimed at obtaining oscillation modes from these objects (see, e.g., Baiotti et al. 2009, for a comparison with linear approaches). It uses pseudospectral methods and a formulation based on primitive variables (Servignat et al. 2023) together with the extended conformal flatness (xCFC) formulation developed in Cordero-Carrión et al. (2009). The conformal flatness condition (CFC; Isenberg 2008; Wilson et al. 1996) has been proved to be a convenient approximation of full General Relativity (GR), which allows us to simplify the formalism while keeping a high degree of accuracy (Iosif & Stergioulas 2014). These approaches reduce the computational cost with respect to other codes that rely on conserved schemes (Banyuls et al. 1997; Cipolletta et al. 2020, 2021; Dimmelmeier et al. 2005, 2006; Thierfelder et al. 2011; Pakmor et al. 2016; Lioutas et al. 2024; Kidder et al. 2017; Rosswog & Diener 2021; Yamamoto et al. 2008), obtaining a lightweight code that can be run on office computers and is thus ideally suited for parametric studies. In this paper, we present an update of ROXAS that allows us to dynamically evolve perturbed, differentially rotating NSs. We describe the main changes to the formalism and code, and present a series of axisymmetric and nonaxisymmetric simulations of differentially rotating stars known as the B sequence (Stergioulas et al. 2004; Dimmelmeier et al. 2006; Krüger et al. 2010; Iosif & Stergioulas 2021), which use a polytropic EoS. We extract their oscillation modes and, when possible, compare their frequencies to those available in the literature for code validation, while the nonaxisymmetric ones in CFC are reported for the first time.

The article is structured as follows. In Sect. 2, we review the formalism used by ROXAS and stress the difference with differential rotation. In Sect. 3, we present the changes implemented in the code and detail the parameters used in our simulations. The results are discussed in Sect. 4, and finally conclusions with some final remarks are given in Sect. 5.

Unless otherwise stated, we use a geometrized unit system with G = c = 1 and a four-metric signature (−, +, +, +). Greek indices denote the four-dimensional space-time coordinates, while Latin ones represent the three-dimensional space. We also use Einstein’s summation convention over repeated indices.

2 Theoretical formalism

The formalism used here is based on the one in Servignat & Novak (2025) and Servignat et al. (2023), with some modifications to account for a differentially rotating star. The full details are provided in these references, and here we only summarize the main aspects and the key differences introduced in this paper.

2.1 Review of formalism

We considered a space-time described by the 3+1 formalism, so that the metric takes the form $$g_{\mu v}d{x}^{\mu }d{x}^{\nu }=-N^{2}d{t}^2+{\gamma }_{ij}(d{x}^i+{\beta }^{i}dt)(d{x}^j+{\beta }^{j}dt),$$(1)

where N is the lapse, βⁱ the shift, and γ_ij the induced spatial metric. We assumed CFC, so that $${\gamma }_{ij}={\mathrm{\Psi }}^4{f}_{ij},$$(2)

where Ψ is called the conformal factor and f_ij is a flat metric. We denoted the covariant derivative associated with γ_ij as D_i and, as described by Bonazzola et al. (2004), we worked in maximal slicing and the Dirac gauge, which is automatically satisfied in CFC. The flat metric f_ij is expressed in spherical coordinates (r, θ, φ), as in Servignat & Novak (2025).

On the hydrodynamics side, we considered a simulated NS to be a perfect fluid, whose energy-momentum tensor is given by $$T^{\mu \nu }=(e+p)u^{\mu }{u}^{\nu }+p{g}^{\mu \nu },$$(3)

where u^μ is the time-like unitary four-velocity, e is the energy density in the fluid frame, and p is the pressure. The fluid is composed of baryons with a number density n_B, at zero temperature, and in β equilibrium. The fluid thermodynamics can then be described by a barotropic EoS, p = p(n_B), which we assume to be polytropic throughout this work, that is, $$p=\kappa n_B^{\gamma },$$(4)

where γ = 2 and κ = 0.02689 in LORENE units¹, defined by c = r₀ = m₀ = 1, with r₀ = 10 km and m₀ ≈ 1.66 × 10²⁹ kg. We also defined the log-enthalpy H, $$H=\ln \left(\frac{e+p}{m_B{n}_B}\right),$$(5)

where m_B ≈ 1.66 × 10⁻²⁷ kg is the nucleon mass. Additionally, we considered the Lorentz factor with respect to the Eulerian observer Γ, the Eulerian velocity of the fluid Uⁱ, and its coordinate velocity vⁱ, which satisfy the equations $$\mathrm{\Gamma }={(1-U_j{U}^j)}^{-1/2},$$(6) $$u^{\alpha }=\frac{\mathrm{\Gamma }}{N}(1,v^i),U^i=\frac{1}{N}(v^i+{\beta }^i).$$(7)

Within this framework, the basic evolution equations are (Servignat et al. 2023; Servignat & Novak 2025) $$\begin{array}{rl}{\mathrm{\partial }}_{t}H=& -v^i{D}_{i}H-c_s^2\frac{{\mathrm{\Gamma }}^{2}N}{{\mathrm{\Gamma }}^2-c_s^2({\mathrm{\Gamma }}^2-1)}\\ & \times [K_{ij}{U}^i{U}^j+D_i{U}^i-\frac{U^i}{{\mathrm{\Gamma }}^2}{D}_{i}H],\end{array}$$(8) $$\begin{array}{rl}{\mathrm{\partial }}_t{U}_i=& -v^j{D}_j{U}_i-U_j{D}_i{v}^j+N{U}_j{D}_i{U}^j-\frac{N}{{\mathrm{\Gamma }}^2}{D}_i(H+\ln N)\\ & +U_i{U}^j{D}_{j}N+\frac{c_s^{2}N{U}_i}{{\mathrm{\Gamma }}^2-c_s^2({\mathrm{\Gamma }}^2-1)}{D}_j{U}^j\\ & +U_i\frac{{\mathrm{\Gamma }}^2(c_s^2-1)}{{\mathrm{\Gamma }}^2-c_s^2({\mathrm{\Gamma }}^2-1)}N{U}^l{U}^j{K}_{lj}\\ & +\frac{(1-c_s^2)N}{{\mathrm{\Gamma }}^2-c_s^2({\mathrm{\Gamma }}^2-1)}{U}_i{U}^j{D}_{j}H.\end{array}$$(9)

In a similar way to Servignat & Novak (2025), we decomposed each field into the equilibrium value, denoted with a subscript “eq”, plus the deviation of the total quantity from this equilibrium, denoted with a bar, $$\begin{array}{rl}& U_i=U_{i,\mathrm{e}\mathrm{q}}+{\overline{U}}_i,H=H_{\mathrm{e}\mathrm{q}}+\overline{H},\\ & {\beta }^i={\beta }_{\mathrm{e}\mathrm{q}}^i+{\overline{\beta }}^i,N=N_{\mathrm{e}\mathrm{q}}+\overline{N}.\end{array}$$(10)

Note that we did not assume that the latter quantities needed to be small.

The next step was to substitute these expressions into Eqs. (8) and (9), and use the identities satisfied at equilibrium to simplify the resulting equations. This is what we call “well-balanced formulation”, and where the main differences between rigid and differential rotation arise. This is not strictly speaking a well-balanced method as defined, for example, by Dumbser et al. (2024), but it is inspired by it.

2.2 Differential rotation

An axisymmetric star at equilibrium with an angular velocity, Ω(r, θ), so that $v_{\text{eq}}^i=\mathrm{\Omega }(r,\theta )r\sin \theta {\delta }_{\phi }^i$, satisfies the equation (Bonazzola et al. 1993; Iosif & Stergioulas 2021) $$H_{\mathrm{e}\mathrm{q}}+\ln N_{\mathrm{e}\mathrm{q}}-\ln {\mathrm{\Gamma }}_{\mathrm{e}\mathrm{q}}+{\int }_{{\mathrm{\Omega }}_p}^{\mathrm{\Omega }}F\left({\mathrm{\Omega }}^{\mathrm{\prime }}\right)d{\mathrm{\Omega }}^{\mathrm{\prime }}=\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{s}\mathrm{t}.,$$(11)

where $$F=u^t{u}_{\phi }=\frac{{\mathrm{\Gamma }}_{\mathrm{e}\mathrm{q}}^2}{N_{\mathrm{e}\mathrm{q}}}{U}_{\phi ,\mathrm{e}\mathrm{q}}r\sin \theta$$(12)

and Ω_p is the angular velocity at the star’s pole.

The introduction of differential rotation thus involves the definition of a rotation law that relates these two variables, F(Ω). Historically, due to its simplicity, the most standard choice has been the Komatsu-Eriguchi-Hachisu (KEH) profile (Komatsu et al. 1989), $$F(\mathrm{\Omega })=A^2({\mathrm{\Omega }}_c-\mathrm{\Omega }),$$(13)

where Ω_c is the central angular velocity and A is a constant that determines the length scale over which the angular velocity, Ω, changes. This can easily be seen in the Newtonian limit, where u^t ≈ 1, g_φφ ≈ (r sin θ)² and u^φ ≈ Ω, which yields the simple rotation law $$\mathrm{\Omega }(r,\theta {)}_{\mathrm{N}\mathrm{e}\mathrm{w}\mathrm{t}}=\frac{{\mathrm{\Omega }}_c}{1+(r\sin \theta /A{)}^2}.$$(14)

In particular, A → ∞ corresponds to the rigidly rotating limit. The KEH profile produces a maximum frequency at the center of the star, which monotonically decreases with the radius, as can be seen in Fig. 3.

More recently, other differential rotation laws have been proposed that more realistically describe the situation in the post-merger remnant, such as those in Uryū et al. (2017). However, these profiles introduce additional computational challenges, as discussed in Sect. 3. In addition, for code validation purposes, it is useful to consider rotation profiles that have been more broadly studied in previous references. For these reasons, in the present paper, our numerical simulations only deal with KEH profiles, leaving more realistic and complex rotation laws for future work. In any case, the formalism described in this section remains general, with no particular assumption regarding the shape of F(Ω).

2.3 The well-balanced formulation with differential rotation

One of the main changes to obtaining hydrodynamics equations with differential rotation in the well-balanced formulation is the integral term that appears in Eq. (11). The other relevant identity that impacts the dynamical equations is $$v_{\mathrm{e}\mathrm{q}}^j{D}_j{U}_{i,\mathrm{e}\mathrm{q}}+U_{j,\mathrm{e}\mathrm{q}}{D}_i{v}_{\mathrm{e}\mathrm{q}}^j=U_{\phi ,\mathrm{e}\mathrm{q}}r\sin \theta D_i\mathrm{\Omega },$$(15)

with the right-hand side vanishing for rigid rotation. Both these equations produce additional terms in Eq. (9) with respect to the rigidly rotating case, leading to the final evolution equation for U_i: $$\begin{array}{rl}{\mathrm{\partial }}_t{\overline{U}}_i=& -[v_{\mathrm{e}\mathrm{q}}^j{D}_j{\overline{U}}_i+{\overline{v}}^j{D}_j{U}_{i,\mathrm{e}\mathrm{q}}+{\overline{v}}^j{D}_j{\overline{U}}_i\\ & +U_{j,\mathrm{e}\mathrm{q}}{D}_i{\overline{v}}^j+{\overline{U}}_j{D}_i{v}_{\mathrm{e}\mathrm{q}}^j+{\overline{U}}_j{D}_i{\overline{v}}^j\\ & -N(U_{j,\mathrm{e}\mathrm{q}}{D}_i{\overline{U}}^j+{\overline{U}}_j{D}_i{U}_{\mathrm{e}\mathrm{q}}^j+{\overline{U}}_j{D}_i{\overline{U}}^j)]\\ & +N{\overline{U}}_j(2U_{\mathrm{e}\mathrm{q}}^j+{\overline{U}}^j)D_i(H_{\mathrm{e}\mathrm{q}}+\ln N_{\mathrm{e}\mathrm{q}})\\ & -\frac{N}{{\mathrm{\Gamma }}^2}{D}_i(\overline{H}+\ln (1+\frac{\overline{N}}{N_{\mathrm{e}\mathrm{q}}}))+U_i{U}^j{D}_{j}N\\ & +\frac{c_s^{2}N{U}_i}{{\mathrm{\Gamma }}^2-c_s^2({\mathrm{\Gamma }}^2-1)}(6U_{\mathrm{e}\mathrm{q}}^{\phi }{\overline{D}}_{\phi }\ln \mathrm{\Psi }+D_j{\overline{U}}^j)\\ & +U_i\frac{{\mathrm{\Gamma }}^2(c_s^2-1)}{{\mathrm{\Gamma }}^2-c_s^2({\mathrm{\Gamma }}^2-1)}N{U}^l{U}^j{K}_{lj}\\ & +\frac{(1-c_s^2)N}{{\mathrm{\Gamma }}^2-c_s^2({\mathrm{\Gamma }}^2-1)}{U}_i{U}^j{D}_{j}H\\ & +\frac{\overline{N}}{N_{\mathrm{e}\mathrm{q}}}{U}_{\phi ,\mathrm{e}\mathrm{q}}r\sin \theta D_i\mathrm{\Omega },\end{array}$$(16)

where the last term is a combination of the two described contributions from differential rotation. Two additional terms appear that were not present in Servignat & Novak (2025), corresponding to higher-order contributions that do not contribute to the evolution in a relevant way in most cases. However, it is safer to account for them, especially for high velocities or if large deviations from equilibrium are involved.

Regarding the equation for the enthalpy, Eq. (8) can be slightly simplified by using the same identity used for rigid rotation, $v_{\text{eq}}^i{D}_i{H}_{\text{eq\hspace{0.17em}}}=0$, leading to the same evolution equation for H: $$\begin{array}{rl}{\mathrm{\partial }}_t\overline{H}& =-{\overline{v}}^i{D}_i{H}_{\mathrm{e}\mathrm{q}}-v_{\mathrm{e}\mathrm{q}}^i{D}_i\overline{H}-{\overline{v}}^i{D}_i\overline{H}\\ -c_s^2& \frac{{\mathrm{\Gamma }}^{2}N}{{\mathrm{\Gamma }}^2-c_s^2({\mathrm{\Gamma }}^2-1)}[K_{ij}{U}^i{U}^j-K+D_i{U}^i-\frac{U^i}{{\mathrm{\Gamma }}^2}{D}_{i}H].\end{array}$$(17)

3 Numerical simulations

The numerical simulations performed in this work were performed with a modified version of the code ROXAS (Servignat & Novak 2025)². The changes made to the code – mainly motivated by the introduction of differential rotation – are summarized below.

3.1 Updates to the code

ROXAS is a pseudospectral code for the dynamical evolution of perturbed, isolated rotating NSs. It makes use of the LORENE (Langage Objet pour la RElativité NumériquE; Gourgoulhon et al. 2016) infrastructure and, in particular, it relies on it to produce the equilibrium star that is perturbed and evolved. Therefore, differential rotation needs to be implemented in both LORENE and ROXAS.

Although tools for computing the equilibrium configuration of a differentially rotating star have been previously implemented in LORENE (Saijo & Gourgoulhon 2006), these solvers were based on the full GR formalism. We implemented new solvers under the CFC assumption, which reduces the number of metric variables and simplifies the system of Einstein equations to be solved.

Despite the fact that we focus on KEH rotation profiles throughout this text, the goal is to be able to use more complex profiles in the future. Profiles such as those in Uryū et al. (2017) need to be specified through a noninjective function Ω(F), which cannot be inverted as F(Ω). While KEH profiles are generally defined through the F(Ω) in Eq. (13), this relation can be inverted and also expressed as an Ω(F), $$\mathrm{\Omega }(F)={\mathrm{\Omega }}_c-\frac{F}{A^2}.$$(18)

With this change, the integral term in the equilibrium equation, Eq. (11), needs to be computed as $${\int }_{F_p}^F{F}^{\mathrm{\prime }}\frac{d\mathrm{\Omega }}{dF}\left(F^{\mathrm{\prime }}\right)d{F}^{\mathrm{\prime }},$$(19)

but otherwise the formalism remains the same. As using Ω(F) is therefore a more generic framework, we chose this way to define the rotation profile for our implementation, in contrast to the original implementation of differential rotation in LORENE.

Note that this step alone is insufficient to deal with the profiles in Uryū et al. (2017), since other challenges remain. For instance, these profiles involve a higher number of free parameters whose values are usually determined through other variables with more clear physical meanings, as can be seen in Iosif & Stergioulas (2021). This implies that these parameters change during the iterative process to find the star’s equilibrium configuration, possibly leading to the algorithm failing or to nonphysical configurations if not treated carefully, especially in the most rapidly rotating cases. Another issue is that these profiles show a local maximum in the angular velocity at a nonzero radius, contrary to KEH profiles (see Fig. 3), which requires more resolution around this point in the dynamical evolution and usually leads to less stable simulations.

Once a NS equilibrium configuration is found, the resulting fields are stored for later usage with ROXAS, including Ω (for use in Eq. (16), for example) instead of F, which can be recovered from other variables. ROXAS was also modified to read equilibrium configurations of differentially rotating stars, as well as to account for the additional terms described in Sect. 2 for a consistent dynamical evolution. Before evolving the star, a perturbation is added to either the Eulerian velocity, U_i, or to the log-enthalpy, H. In the latter case, the perturbation can be chosen to be $$\delta H=\epsilon (x^2\pm y^2-\frac{2}{3}{r}^2)(1-\frac{r^2}{R_S(\theta {)}^2}),$$(20)

where R_S(θ) is the radius of the equilibrium star in the direction of θ. In this equation, the minus sign (−) corresponds to a nonaxisymmetric perturbation that mixes l = m = 0 and l = m = 2 modes, as defined and used in Servignat & Novak (2025), while the plus sign (+) yields axial symmetry, which is useful for saving simulation time if only an axisymmetric study is required.

Finally, ROXAS is now able to export one-dimensional profiles and two-dimensional slices of the hydrodynamic and metric grids during the simulation, where the output rate can be specified independently. They can easily be plotted afterwards, and allow us to have a better understanding of how the simulation internally behaves, and to spot potential issues. Examples with the initial profile and slice of the enthalpy perturbation of a nonaxisymmetric simulation are shown in Figs. 1 and 2, respectively.

Fig. 1 Radial profile in different directions of the initial enthalpy perturbation of the nonaxisymmetric B4 simulations.

Fig. 2 Slice in the xz plane of the initial enthalpy perturbation of the nonaxisymmetric B4 simulations. The dashed green lines correspond to the boundaries between the different domains of the metric grid, as described in Sect. 3.2.

3.2 Simulation parameters

In order to test the new implemented changes, we checked the consistency with other numerical simulations of oscillations in differentially rotating stars in the literature. In particular, we used the B sequence described in Dimmelmeier et al. (2006) and compared the resulting frequencies with those provided in this reference.

For the B sequence, the central value for the log-enthalpy is the same as in the BU sequence used in Servignat & Novak (2025), which in LORENE units is H_c = 0.2280. For the rotation profile in Eq. (13), we set A = r_e, with r_e being the star radius at the equator. The B sequence is defined by the star’s flattening (ratio between the star’s radius at the pole, r_p, and r_e). However, LORENE does not take this as an input parameter, and requires the central angular velocity Ω_c. Therefore, we took these values from Dimmelmeier et al. (2006) and checked that the final configurations found by the algorithm yielded similar values of r_p/r_e, masses, and ratios between kinetic and potential energies, T/W. A summary of these parameters and their relative errors with respect to those in Dimmelmeier et al. (2006) is provided in Table 1. We can see that the agreement is excellent, even for the fastest rotating stars. The largest errors correspond to T/W for B1 and B2, due to the smallness of these quantities compared to the precision of the reference values.

With respect to the dynamical evolution with ROXAS, the simulation properties are very similar to the ones in Servignat & Novak (2025). Two types of grids are used, with the following parameters:

• The hydro grid consists of two domains: the spherical nucleus, with a maximum radius of 0.7 r_p, and a shell that covers the rest of the star. In both domains, we set N_r = N_θ = 17, with N_i, i = r, θ, φ, the number of grid points for each coordinate. In addition, N_φ is set to 8 for the nonaxisymmetric simulations and to 1 for the axisymmetric ones. In Servignat & Novak (2025), N_θ is increased for the most demanding (i.e., more rapidly rotating) simulations, but we found them to be more stable by keeping the same value;

• The metric grid consists of several spherical domains: a nucleus that extends up to 8 km (7 km for B9, which has r_p ≈ 7.5 km), a first shell from this radius up to r_p, four more shells up to r_e, and then one more up to 20 km followed by the compactified external domain, in which a change in variable u = 1/r is made to take into account boundary conditions at r → ∞ (Bonazzola et al. 1998). The four shells between r_p and r_e are chosen so that their edges match R_S(θ) for θ uniformly spaced between 0 and π/2.

4 Results

We carried out an initial test of the code to check its ability to keep stable equilibrium initial data, with the evolution of unperturbed differentially rotating NS. Some angular velocity profiles corresponding to an unperturbed B4 are thus shown in Fig. 3, both initially and after a dynamical evolution in CFC of 25 ms. The fact that these profiles remain constant after this time period (many rotation periods) shows that the formalism, the dynamical code, and the equilibrium star solvers work together consistently and in a stable way.

Then, for each of the differentially rotating NSs in equilibrium reported in Table 1, we performed four simulations: two in the Cowling approximation and another two with dynamical space-time under CFC. For each pair, one of the simulations is evolved in axisymmetry after adding an axisymmetric perturbation of the form (20) (+ sign) to the equilibrium star, while the other one is evolved in 3D³ after adding a nonaxisymmetric perturbation of the form (20) (− sign). In all cases, the amplitude is chosen to be ϵ = 10⁻³ in code units. The GWs were extracted following the procedure in Servignat & Novak (2025), based on Einstein’s quadrupole formula; an example waveform is provided in Fig. 6.

Fig. 3 Angular velocity profiles in different directions for an unperturbed B4 simulation. The solid lines correspond to the initial profiles, while the superimposed points are plotted after 25 ms of dynamical evolution in CFC.

4.1 Cowling approximation

The Cowling approximation assumes that the metric stays static and equal to the values at equilibrium. This approximation has been widely used for the simulation of relativistic stars (see, e.g., McDermott et al. 1983; Villain & Bonazzola 2002; Sotani et al. 2023; Montefusco et al. 2025) due to the reduction in the computational cost of simulations. While this is not particularly relevant for ROXAS, running some simulations in the Cowling approximation allows us to perform more extensive code validation against existing literature. In particular, axisymmetric modes can be compared to those in Stergioulas et al. (2004), which used an axisymmetric code for dynamical evolution in GR under the Cowling approximation. On the other hand, nonaxisymmetric modes can be compared to Krüger et al. (2010), where these frequencies were obtained with a perturbative approach, once again limited to the Cowling approximation.

When possible, these simulations were run for 25 ms; however, for the most demanding configurations, the nonaxisymmetric simulations are only stable up to shorter times, which are reported in Table 2. In these cases, at late stages, the fields develop unphysical, large, high-order coefficients in the spectral representation, leading to an eventual numerical crash, without preceding signatures of physical instabilities such as shock formation. Possible causes include insufficient resolution or aliasing effects; we will investigate these further and attempt to mitigate them in future works.

An example frequency spectrum for the nonaxisymmetric B4 simulation is provided in Fig. 4, where we can distinguish the different oscillation modes present in the simulation. In particular, a secondary fundamental mode, F_II, appears around 2.6 kHz, as was reported for the first time in Stergioulas et al. (2004). As F_II and ²f₋₂ have similar values, they generate a single joint peak in Fig. 4. However, this secondary fundamental mode can be identified independently in the absence of an ²f₋₂, such as in axisymmetric simulations or in spectra corresponding to l = m = 0 modes, as in Fig. 5 presented in the next subsection.

The frequencies for the full B sequence are presented in Table 2, where the precise values are extracted from the Fourier transform of the decomposition of the star’s radius, R_S, into spherical harmonics. These time series include fewer oscillation modes per multipole order, leading to more isolated and cleaner peaks and, therefore, to a more precise frequency determination.

The axisymmetric frequencies we extracted are the fundamental modes F and F_II, their first overtone H₁, the l = 2, m = 0 mode ²f₀, and its first overtone ²p₁. The values shown in Table 2 correspond to those obtained from the axisymmetric simulations, given their better stability and lower number of modes, which generate cleaner peaks. However, as these frequencies are also present in the nonaxisymmetric simulations, we extracted them from these simulations as well, and obtained the same values up to maximum relative errors of 0.9%. On the other hand, the extracted nonaxisymmetric frequencies are the l = |m| = 2 modes, denoted by ²f₂ and ²f₋₂, and appear only in the nonaxisymmetric simulations.

The relative errors of our results with respect to Stergioulas et al. (2004) and Krüger et al. (2010) are provided in Table 2. Most of them are similar or below 2%. The larger values of δ²f₂ from B6 on can be explained by the smallness of these frequencies compared to our frequency resolution. Similarly, the 3.16% error in F_II for B2 is due to the fact that the peaks for F and F_II are completely blended in our B2 spectrum, which leads us to report them as equal.

Fig. 4 Gravitational waves spectrum of the nonaxisymmetric B4 simulation in the Cowling approximation.

Fig. 5 Fourier spectrum of the l = m = 0 RS coefficient for the nonaxisymmetric B4 simulation in the Cowling approximation (left) and a dynamical space-time (right).

4.2 Full dynamics in CFC

Having checked that ROXAS produces consistent results with other articles within the Cowling approximation, we were able to relax this assumption and work with the Einstein-CFC system solved dynamically (see Servignat & Novak (2025) for details), using the full potential of the code. In this case, the waveform spectra show a single fundamental mode, unlike the double peak seen in the Cowling approximation. We can see this more clearly by comparing the spectra of the l = m = 0 coefficient of R_S, from which the fundamental frequencies and their overtones are extracted, as can be seen in Fig. 5. In this figure, we can see the two fundamental-mode peaks around 2.6 kHz in the Cowling simulation, while there is a single fundamental-mode peak around 1.3 kHz in the CFC simulation. This confirms that the secondary fundamental mode is an artifact of the Cowling approximation in differential rotation, as was suspected in Stergioulas et al. (2004). While previous works such as Dimmelmeier et al. (2006) studied axisymmetric modes without the Cowling approximation, this is the first time this difference is explicitly addressed by comparing both types of simulation.

We show an example waveform in Fig. 6, with its spectrum in Fig. 7, that corresponds to the nonaxisymmetric B4 simulation in CFC. In Table 3, we report the frequencies for the same axisymmetric and nonaxisymmetric modes as those described in the previous section, this time for the B sequence in a dynamical space-time. The only exception is the secondary fundamental mode, which is absent in this case. Once again, the axisymmetric modes in Table 3 correspond to the axisymmetric simulations, while we use the nonaxisymmetric ones to check the consistency of these modes and to extract the nonaxisymmetric frequencies. The maximum relative error between the axisymmetric modes of both types of simulation is 6.5%, which corresponds to H₁ for B9. This is an outlier caused by a higher-frequency peak (≈4.45 kHz), which merges with the one corresponding to H₁ in the nonaxisymmetric simulation. Both of them can be resolved in the axisymmetric simulation, as would likely be the case with a longer nonaxisymmetric simulation. Other than that, the largest relative error is 2.4%.

The nonaxisymmetric modes are reported here for the first time, to our knowledge. The ²f₂ frequency becomes negative from B7 onward, indicating that this mode rotates in the opposite direction with respect to the previous cases, in the star’s co-rotating reference frame. For low rotation frequencies, the (2, 2) mode corresponds to a prograde rotation, while the (2, −2) is retrograde. The regime where the (2, 2) mode is also retrograde is called the CFS (Chandrasekhar-Friedman-Schutz) instability (Chandrasekhar 1970; Friedman & Schutz 1978; Friedman & Schutz 1978), as reported in Krüger et al. (2010) within the Cowling approximation and in Zink et al. (2010) and Krüger & Kokkotas (2020) for rigid rotation. With our approach, both rotation directions can be distinguished with a Fourier transform of the complex l = m = 2 R_S coefficient. We can see the different behaviors in Fig. 8.

On the other hand, the axisymmetric modes can be compared to the values quoted in Dimmelmeier et al. (2006), in which axisymmetric, dynamical simulations were performed in CFC. The relative errors between this work and ours are also included in Table 3.

We can see that the agreement of H₁, ²f₀, and ²p₁ is of the order of 1% or lower in most of the cases. For the H₁ corresponding to B9, a similar issue to ours might have also arisen in Dimmelmeier et al. (2006), as their reported frequency is closer to the misleading 4.45 kHz found in our simulations, which leads to a similar relative error of 6.5%. The fundamental frequency, F, however, shows larger discrepancies, especially in the fastest rotating cases. Taking into account the uncertainties purely associated with the frequency resolution limited by the total simulation time (see Fig. 9), there is only a clear tension for the B7 and B9 simulations. However, the fundamental frequencies in Dimmelmeier et al. (2006) show a faster decaying trend than ours. The origin of these discrepancies is unclear at this point, as the codes, resolutions, and simulation times of both works are different, which requires a more detailed investigation. Note that the CFC approximation, which is also used in Dimmelmeier et al. (2006), is not the cause. While a comparison with frequencies obtained directly from GR is not possible here, it was done for the rigidly rotating case in Servignat & Novak (2025) with excellent agreement, showing the reliability of CFC for NS oscillation-mode studies.

Fig. 6 Gravitational wave strain extracted from the nonaxisymmetric B4 simulation in CFC.

Fig. 7 Gravitational wave spectrum of the nonaxisymmetric B4 simulation in CFC.

Fig. 8 Fourier spectrum of the complex l = m = 2 RS coefficient for the nonaxisymmetric B4 (left) and B8 (right) simulations in CFC.

Fig. 9 Fundamental frequencies obtained for the B sequence simulations and comparison with those in Dimmelmeier et al. (2006). The error bars are the inverse of the simulation times, which correspond to the frequency resolutions.

5 Conclusions

In this work, we present an update of the publicly available code ROXAS, which can now simulate differentially rotating stars, together with several additional improvements. We ran a series of both nonaxisymmetric and axisymmetric simulations, first in the Cowling approximation and then with a dynamical space-time in CFC. We verified their consistency with different references and obtained excellent overall agreement for both axisymmetric and nonaxisymmetric mode frequencies. In addition, we show that, as conjectured in Stergioulas et al. (2004), the secondary fundamental peak appearing in simulations within the Cowling approximation does not appear in dynamical space-times under the same initial conditions.

Furthermore, the nonaxisymmetric frequencies ²f₂ and ²f₋₂ in CFC are reported here for the first time, illustrating the ability of ROXAS to explore a wider class of configurations and generate novel results. We also emphasize that these simulations are performed with a lightweight code that runs efficiently on standard computers within short computation times, even though it is not parallel. For radial (one-dimensional) simulations, ROXAS typically performs about five times faster (Servignat et al. 2023) than the reference core-collapse and NS oscillation code CoCoNuT (Dimmelmeier et al. 2005) used in its sequential version. Assuming that this improvement scales equally in all dimensions and no paralellization, this effect should be enhanced for axisymmetric (giving roughly a factor of 5² = 25) and nonaxisymmetric (a factor of 5³ = 125) simulations. Nevertheless, these are order-of-magnitude estimates and a direct comparison through benchmarking of multi-dimensional simulations remains to be done in order to properly assess the difference in performance.

Together with Servignat & Novak (2025), this work represents an initial step toward a versatile tool for studying perturbed, rotating NSs with moderate computational resources. Building on this, future works will include differential rotation profiles, which are more motivated by binary NS simulations, such as those in Uryū et al. (2017), as well as more realistic EoSs, which will require further extensions of the present framework.

Acknowledgements

The authors acknowledge support from the Agence Nationale de la Recherche (ANR) under contract ANR-22-CE31-0001-01.

References

All Tables

All Figures

	Fig. 1 Radial profile in different directions of the initial enthalpy perturbation of the nonaxisymmetric B4 simulations.
In the text

	Fig. 2 Slice in the xz plane of the initial enthalpy perturbation of the nonaxisymmetric B4 simulations. The dashed green lines correspond to the boundaries between the different domains of the metric grid, as described in Sect. 3.2.
In the text

	Fig. 3 Angular velocity profiles in different directions for an unperturbed B4 simulation. The solid lines correspond to the initial profiles, while the superimposed points are plotted after 25 ms of dynamical evolution in CFC.
In the text

	Fig. 4 Gravitational waves spectrum of the nonaxisymmetric B4 simulation in the Cowling approximation.
In the text

	Fig. 5 Fourier spectrum of the l = m = 0 RS coefficient for the nonaxisymmetric B4 simulation in the Cowling approximation (left) and a dynamical space-time (right).
In the text

	Fig. 6 Gravitational wave strain extracted from the nonaxisymmetric B4 simulation in CFC.
In the text

	Fig. 7 Gravitational wave spectrum of the nonaxisymmetric B4 simulation in CFC.
In the text

	Fig. 8 Fourier spectrum of the complex l = m = 2 RS coefficient for the nonaxisymmetric B4 (left) and B8 (right) simulations in CFC.
In the text

	Fig. 9 Fundamental frequencies obtained for the B sequence simulations and comparison with those in Dimmelmeier et al. (2006). The error bars are the inverse of the simulation times, which correspond to the frequency resolutions.
In the text