A new code for computing differentially rotating neutron stars

TL;DR

The paper develops and validates new initial-data solvers in the FUKA framework for constructing differentially rotating neutron stars in full general relativity, implementing both Quasi-Isotropic Coordinates (QIC) and eXtended Conformal Thin Sandwich (XCTS) formulations with the Komatsu–Eriguchi–Hachisu (KEH) rotation law. It extends EOS support to 3D tabulated tables via GRHayLEOS and demonstrates exponential spectral convergence, cross-solver consistency, and agreement with the RNS code across a range of differential-rotation profiles, axis ratios, and equations of state. It also investigates how initial-data resolution affects dynamical evolutions using the Einstein Toolkit and IllinoisGRMHD, showing that evolution resolution dominates error growth and confirming robust bar-mode dynamics for polytropic and tabulated EOSs. The work lays groundwork for future expansions to include magnetic fields, non-isentropic flows, and binary systems, offering a publicly available, high-accuracy toolkit for modeling differentially rotating relativistic stars with realistic microphysics.

Abstract

We present new initial data codes for constructing stationary, axisymmetric equilibrium models of differentially rotating neutron stars in full general relativity within the Frankfurt University/KADATH (FUKA) suite of initial data codes. FUKA leverages the KADATH spectral library to solve the Einstein equations under the assumption of an isentropic fluid without magnetic fields while incorporating GRHayLEOS to support 3D tabulated equations of state in \textit{stellar collapse} format. The two solvers explored in this work include one using quasi-isotropic coordinates (QIC) in Spherical coordinates while the other solves the eXtended Conformal Thin Sandwich (XCTS) decomposition in Cartesian coordinates, enabling the construction of equilibrium configurations with high accuracy and efficiency. In this work we adopt the Komatsu-Eriguchi-Hachisu differential rotation law, however, the code is designed to be extensible to other rotation laws, allowing for exploration of physically relevant sequences and critical rotation thresholds. Furthermore, we perform convergence tests demonstrating the exponential accuracy of the spectral approach, we validate QIC and XCTS solutions against models well-studied in the literature, and we also compare FUKA solutions against the well-known RNS code. Finally, we explore the impact that initial data resolution has on dynamical simulations and recover the convergence order of the evolution scheme, the dominate source of error in this study. The new FUKA codes and results presented here lay the foundation for future extensions to more general configurations, including magnetic fields, removal of isentropic assumptions, and binary systems, and have been made publicly available to support community efforts in modeling differentially rotating relativistic stars.

Figures (10)

• Figure 1: Self-convergence of the XCTS (left) and QIC (right) solvers for U13 and SB6 configurations, respectively. Exponential convergence is observed up to $\tilde{N} \sim 22$ at which point the error is dominated by spectral noise and solutions begin to deteriorate. We include $\varepsilon_{\rm STOP}$, the stopping threshold for the Newton-Raphson solver in KADATH, to highlight the expected precision floor.
• Figure 2: Comparison between diagnostic quantities computed using QIC and XCTS solutions for a range of radial collocation points $R$. We find excellent agreement of $M_K$ and $M_{ADM}$ as the radial resolution is increased, however, fluid quantities such as $M_b$ and $\rho_c$ only show convergence up to $10^{-3}$.
• Figure 3: Comparison of solutions between the FUKA QIC solver (circles) and the RNS code (solid lines). Here we explore a range of differential rotation profiles $\hat{\mathcal{A}} = \left\lbrace 0.5 \,, 1.0 \,, 2.0 \right\rbrace$ (left, middle, right panes, respectively) and axis ratios $r_p / r_e = \left\lbrace 0.4 \,, 0.6 \,, 0.8 \,, 1.0 \right\rbrace$ (red, orange, pink, and blue lines, respectively) over a range of central densities. No solutions were found with FUKA for $\hat{\mathcal{A}} = 2.0$ and $r_p / r_e = 0.4$, but the RNS solution is included for completeness and to highlight the change in behavior in this region of the parameter space.
• Figure 4: Convergence analysis of the Hamiltonian constraint violations along the $x$-axis at $t = 41 M$ where we find the expected second-order convergence as predicted by the hydro evolution code, IllinoisGRMHD.
• Figure 5: Top: Evolution of the $m = 2$ Fourier mode of the stellar interior for the XCTS solution (left) and QIC solution (right) for the U13 configuration. Both are consistent with behavior reported in Ref. Baiotti2006. Bottom: Convergence analysis of the $m = 2$ Fourier mode of the stellar interior during dynamical evolution for the U13 configuration to assess the impact of initial data resolution $R = \left\lbrace11 \,, 13 \,, 17\right\rbrace$ at evolution resolution $\Delta_H x = 0.32$ and $\Delta_M x = 0.4$. For the QIC results, we have aligned the datasets at the onset of the bar-mode instability to focus on the long-term behavior of each configuration.