FUKA Suite: Relativistic Initial Data Framework

• FUKA Suite is a high-order spectral-initial-data framework that models equilibrium states for compact objects like neutron stars and black holes.
• It utilizes multidomain pseudospectral decomposition and Newton–Raphson solvers with automatic differentiation to solve general-relativistic constraint equations.
• Validation studies demonstrate exponential convergence and robust performance across both QIC and XCTS solver modules.

The Frankfurt University/KADATH (FUKA) suite is a publicly available, high-order spectral-initial-data framework designed for the construction of equilibrium models and initial data for a broad spectrum of relativistic astrophysical systems. Targeting compact objects such as isolated and binary neutron stars, binary black holes, and black hole–neutron star binaries, FUKA provides accurate solutions of the general-relativistic constraint equations via state-of-the-art spectral methods. It leverages the KADATH spectral library as its computational backbone, incorporating multidomain pseudospectral decomposition, Newton–Raphson nonlinear solvers with automatic differentiation, and sophisticated boundary/domain management. The suite facilitates precise modeling of rotating and differentially rotating neutron stars, including full general-relativistic configurations with 3D tabulated equations of state, and serves as a frontend for generating initial data compatible with major evolution codes in numerical relativity (Tootle et al., 8 Jan 2026, 0909.1228).

1. Architecture and Computational Foundations

FUKA is architected as a modular, layered system in which high-level astrophysical problem setup is mapped onto multidomain spectral PDE solvers provided by KADATH. KADATH handles field and PDE definitions—expressed via a LaTeX-like interface—across multiple geometric subdomains (e.g., spherical, bispherical, cylindrical-critical), each equipped with optimal spectral bases (Chebyshev or Legendre polynomials radially, trigonometric or associated Legendre functions angularly, and Fourier modes azimuthally). Domains can be adapted so that interfaces precisely coincide with physical surfaces, such as the stellar surface, yielding spectral regularity and eliminating Gibbs phenomena at surface discontinuities.

The solver workflow proceeds via global assembly of the discretized nonlinear system $\mathbf{f}(\mathbf{u})=0$ in coefficient space, where $\mathbf{u}$ collects all spectral coefficients of the unknowns. Nonlinear systems are iteratively solved using a Newton–Raphson scheme, with Jacobians computed by automatic differentiation using dual numbers. Convergence is tracked by monitoring $\|\mathbf{f}\|_\infty$ and relative updates $\|\Delta\mathbf{u}\|/\|\mathbf{u}\|$, with default residual tolerance $\varepsilon_\text{STOP}=10^{-8}$ (0909.1228). Parallelization for large systems employs MPI and ScaLAPACK for distributed direct linear algebra.

2. Solver Modules: QIC and XCTS Approaches

FUKA implements two complementary solver modules, each adapted to distinct coordinate and decomposition strategies for neutron star equilibrium:

A. Quasi-Isotropic Coordinate (QIC) Solver:

This solver operates in axisymmetric spherical coordinates, $(t, r, \theta, \phi)$, with the metric expressed as

$$g_{\mu\nu}dx^\mu dx^\nu = -\alpha^2 dt^2 + A^2(dr^2 + r^2 d\theta^2) + B^2 r^2 \sin^2\theta (d\phi + \beta^\phi dt)^2$$

with four scalar fields $(\nu=\ln\alpha,\, \omega=-\beta^\phi,\, A,\, B)$ as unknowns. Einstein’s equations under stationarity and axisymmetry reduce to a coupled set of Laplacian-like PDEs, with fluid profiles governed by the Komatsu–Eriguchi–Hachisu (KEH) differential rotation law:

$$\Omega(r,\theta)=\frac{\Omega_c}{1+\hat{A}^2 r^2 \sin^2\theta}$$

and the corresponding hydrostationary first integral. The domain comprises a polar nucleus, adapted polar shells fit to the stellar surface, and a compactified exterior shell ($r\to\infty$).

B. eXtended Conformal Thin Sandwich (XCTS) Solver:

This solver formulates the initial data problem in the conformally flat, maximal-slicing XCTS framework, using Cartesian spectral expansions mapped onto spherical-like domains. The elliptic system is solved for the conformal factor $\Psi$, shift vector $\beta^i$, and conformal lapse $\alpha\Psi$, under the conditions $\gamma_{ij} = \Psi^4 \delta_{ij}$ and $K = 0$. Fluid sources and the hydrostationary first integral mirror those in the QIC module, but the underlying fields are 3D. A hybrid initialization chain—starting with a converged QIC solution mapped to Cartesian XCTS collocation—enables rapid high-precision convergence.

3. Spectral Methods and Domain Decomposition

KADATH implements tensor-product spectral expansions, with bases adapted to each coordinate type. Radially, Chebyshev or Legendre polynomials are used, while the angular directions employ appropriate trigonometric, Legendre, or Fourier expansions. For 3D fields,

$$f(r^*,\theta,\varphi)=\sum_{k=0}^{N_r}\sum_{\ell,m} a_{k\ell m}\,\Phi_k(r^*)\,\Theta_{\ell m}(\theta)\,\Phi_m(\varphi)$$

Domain decomposition splits the physical domain into nucleus (for $r=0$), adapted shells (stellar interior, exterior), and an outer compactified shell. Each unknown is spectrally expanded with boundary regularity and matching enforced via the $\tau$-method and domain-matching conditions. Spectral collocation points are chosen (e.g., Gauss–Lobatto), and PDEs are discretized via a weighted residual formalism, combining $\tau$- and Galerkin approaches for strict regularity.

Automatic differentiation propagates through all field operations, providing full Jacobian assembly for the Newton–Raphson solver. For typical neutron star initial data, exponential convergence is observed with effective spectral resolutions $\tilde{N}\sim 20-25$, beyond which roundoff and physical gradients limit accuracy (Tootle et al., 8 Jan 2026).

4. Equations of State, Rotation Laws, and Fluid Modeling

FUKA supports a wide class of equations of state (EOS):

• Piecewise-polytropic EOS,
• Cold tabulated EOS (LORENE format),
• 3D tabulated stellar-collapse EOS via GRHayLEOS interface.

In its current implementation, FUKA restricts to cold, beta-equilibrium slices ($ds=0$), with future development planned for arbitrary ($\rho, T, Y_e$) thermodynamic dependence. Rotation is governed by the KEH law, with extensibility for alternative differential and multiparameter rotation profiles. Major fluid unknowns—enthalpy $h$, Lorentz factor $W$, three-velocity $U$—and stress-energy projections are computed consistently for both QIC and XCTS modules, with identical hydrostationary integrals.

5. Validation, Self-Convergence, and Cross-Verification

Extensive convergence and validation studies demonstrate the robustness and accuracy of FUKA:

• Self-convergence: Both QIC and XCTS solvers exhibit exponential error decay in global diagnostics (ADM mass, Komar mass, baryonic mass, central density) up to $\tilde N \simeq 22$, with convergence floors comparable to Newton tolerance.
• QIC–XCTS agreement: For reference configurations (SB6, U13 of Baiotti et al.), ADM and Komar masses match at $\sim10^{-8}$, while baryonic mass and central density differences are $\sim10^{-3}$.
• QIC–RNS code comparison: For sequences scanning differential rotation ($\hat{A}=\{0.5,1,2\}$) and axis ratio ($r_p/r_e=\{0.4,0.6,0.8,1\}$), FUKA's $M_b(\rho_c)$ curves agree with those of the RNS code at the percent level, with minor discrepancies at the most extreme rotation.

Dynamic evolutions using generated initial data confirm that, at fixed spectral (initial data) resolution $R \gtrsim 13$, the dominant error source in subsequent evolutions is the finite-difference evolution scheme, not the initial data resolution.

6. Extensibility, Community Infrastructure, and Impact

FUKA’s modular architecture, built atop KADATH, admits rapid extension to:

• New rotation laws (e.g., Uryu et al., multi-parameter families, Camêlio et al. models),
• Non-isentropic, finite-temperature fluids (by solving for $T$ in future releases),
• Magnetic fields (by incorporating Maxwell stress tensors in the XCTS module),
• Non-axisymmetric, fully 3D equilibrium solutions (including triaxial or one-armed instabilities),
• Binary initial data, including binary neutron star, binary black hole, and black hole–neutron star systems,
• Microphysical neutrino and transport coupling via further GRHayLEOS/GRHayL development.

The suite is maintained with comprehensive documentation and is integrated with major evolution codes (Einstein Toolkit, SpECTRE, BAM, SACRA, SPHINCS_BSSN). FUKA's open-source availability supports a broad community engaged in the modeling and simulation of relativistic stars and compact-object mergers, establishing it as a versatile and high-accuracy platform for equilibrium and initial-data studies (Tootle et al., 8 Jan 2026).

7. KADATH Library: Technical Integration and Numerical Methods

KADATH is a modular library applying multidomain spectral methods for generalized PDE systems relevant to theoretical physics (0909.1228). Its architectural layers comprise:

• Geometry/domain layer: Abstract “Space” and “Domain” classes (spherical, bispherical, and cylindrical-critical geometries with domain-matching logic).
• Spectral basis and collocation: Supports Chebyshev, Legendre, and Fourier bases, with interpolation ($I$) and differentiation ($D$) matrices, regularity via Galerkin constraints.
• System assembly: Encapsulates unknown fields, constants, and LaTeX-like equation input for global system construction; PDEs discretized using the weighted-residual method.
• Solver: Newton–Raphson with automatic differentiation, integrated distributed linear algebra (LAPACK, ScaLAPACK), and MPI-based parallelization.

KADATH’s application examples corroborate accuracy and performance for various theoretical physics domains (gauge theory, gravitational collapse, rotating black holes, and binaries) with exponential convergence and scalable solve times (see performance table in (0909.1228)).

Integration of KADATH into FUKA is achieved via C++ API, with high-level FUKA problem setup translated to KADATH field and equation definitions, spectral builds, and output formats. This composite framework underpins FUKA’s breadth and computational efficiency, establishing it as a reference platform for relativistic initial data modeling.