Stamatis Vretinaris

This thesis studied numerical modeling of compact binary systems across regimes in which the relevant physical and numerical difficulties are qualitatively different. Comparable-mass binaries require the nonlinear dynamics of spacetime and, for neutron stars, matter. Extreme-mass-ratio systems require a controlled treatment of an unresolved small body and its self-field. Across both regimes, the reliability of the calculation depends on preserving continuum structure at the discrete or statistical level: geometric quantities must retain their meaning, evolution systems must remain stable, singular or unresolved sources must be represented consistently, and inference problems must be constrained by the available physics.

The work presented here consists of four connected contributions. The first constructs high-order summation-by-parts discretizations for curvilinear coordinates with coordinate singularities. The second proposes a finite-resolution formulation for scalar self-force calculations. The third concerns the quasi-local geometry of binary black-hole mergers. The fourth develops a robust inference strategy for binary-neutron-star post-merger gravitational waves. Together, these results form a thesis on compact-binary modeling rather than on a single computational technique.

High-order SBP operators for curvilinear singularities
This thesis begins with the fundamental requirement of numerical relativity: that the discretized evolution equations should preserve the mathematical structure of the underlying continuum physics. Coordinates adapted to symmetry, such as spherical coordinates, are natural in many relativistic calculations. They also contain coordinate singularities. The factors 1/rp that appear at the origin are not physical singularities, but they can destabilize a finite-difference method if the discrete operators do not reproduce the correct integration-by-parts structure.

The work on summation-by-parts operators constructed high-order finite-difference schemes for first-order hyperbolic systems in such coordinates. The construction starts from Cartesian SBP operators but modifies the norm and derivative structure to account for the curvilinear volume measure and the distinction between scalar and vector-valued quantities. Non-diagonal vector-valued norm operators are allowed. Accuracy conditions near the origin determine the closure coefficients, and the covariant divergence operator is derived from the SBP identity itself.

Explicit operators were constructed up to sixth order. In spherical scalar-wave evolutions, the schemes conserve the discrete energy to floating-point accuracy under reflecting boundary conditions. A Gaussian wave packet can pass through the origin and re-emerge cleanly, with fourth- and sixth-order convergence observed for the corresponding operators. These tests demonstrate that the origin can be included on a collocated radial grid without sacrificing the energy stability expected from an SBP discretization.

This result supplies numerical infrastructure for relativistic evolution problems. It provides a systematic alternative to excising the origin or using staggered grids solely to avoid the coordinate singularity. Because the construction is based on weighted norms, compatible gradient and divergence operators, and the SBP identity, it is not restricted in principle to scalar waves. The method can be extended to arbitrary tensor rank and spacetime dimension, although its integration into full numerical-relativity systems remains future work.

Finite-resolution scalar self-force modeling
The scalar self-force work addressed the extreme-mass-ratio side of the thesis. In such systems, a small compact object moves in the field of a much larger body and is influenced by its own field. A literal point-particle description introduces singular fields, which must then be regularized. In self-consistent evolutions, acceleration-dependent singular-field terms can lead to implicit or higher-derivative structures in the equations of motion. This is undesirable for an initial-value evolution scheme.

The approach developed here starts from a different premise: a numerical calculation has finite resolution. The relevant hierarchy is one in which the large scale is resolved, the small object is not, and the grid scale lies between them. In this regime the small object can be modeled as finite-sized below the resolution scale. The particle field is treated as an additional basis function, and the evolved field is decomposed into this particle field plus a regular remainder field. Only the sum has physical meaning.

The formulation was developed for a scalar particle coupled to a scalar field in flat spacetime. The equations were written explicitly in the lab frame, and the state vector was chosen before deriving the equations of motion. This leads to a self-consistent initial-value system that includes the radiation produced by particle acceleration without introducing infinities or requiring higher time derivatives of the self-acceleration as independent data.

Two discretizations were implemented: a finite-volume method and a spectral method. In one-dimensional tests, both showed excellent agreement with analytic solutions. The fully self-interacting case also agreed well with analytic approximations. The calculations emphasized that effective sources generated by particle acceleration must be projected properly into the discrete function space; this projection can be performed analytically or numerically, with the latter being attractive when the source expressions become complicated.

The present work is a model problem, but its conceptual implication is broader. It suggests that effective-source modeling can be formulated directly at finite resolution. In higher dimensions the effective sources will have stronger local structure, including directional discontinuities. This motivates discretizations adapted to the particle neighborhood, such as finite-element methods with local sub-element structure. Window functions can be used to localize the particle field and reduce computational cost. Extending the construction to electrodynamics and ultimately to general relativity would require additional degrees of freedom, including spin for compact objects, but the central idea remains the same: the unresolved object should be represented in a way that is consistent with the numerical scale at which the field is evolved.

Black-hole horizon mergers and cusp formation
The study of cusp formation in merging black-hole horizons approached the binary black-hole merger through marginally outer trapped surfaces. In this quasi-local picture, the merger can be understood as the merger of the two individual MOTSs associated with the progenitor black holes with the inner MOTS associated with the remnant. For the head-on collision of two nonspinning black holes, the inner MOTS develops a cusp at the touching time.

The main result is that the time evolution of the horizon mass and mass multipole moments is smooth except at this touching time. At the cusp, the curvature contains a distributional contribution, and the quasi-local quantities exhibit discontinuous jumps. The jump in the mass follows from the additivity of the area at the merger. The jumps in the higher mass multipole moments are controlled by the Gauss invariant, which determines the strength of the distributional curvature at the cusp. Although these jumps are not the entire change between the initial progenitor horizons and the final remnant, they can be among the dominant contributions.

This result gives a source-side description of part of the merger dynamics. It shows that horizon multipoles contain information that is localized at the merger event and that can be lost rapidly in the subsequent smooth evolution of the common horizon. Accurate numerical resolution of the MOTSs near the touching time is therefore essential. The analysis also suggests a broader program: to connect horizon multipoles, angular momentum, recoil, and gravitational-wave mode amplitudes in generic black-hole mergers. The present work treated a symmetric head-on configuration. Extending the analysis to spinning, precessing, and non-axisymmetric binaries is the natural next step.

Post-merger inference for binary neutron stars
The binary-neutron-star post-merger work addressed a different aspect of comparable-mass dynamics: the extraction of physical parameters from the gravitational-wave signal of the remnant. The post-merger spectrum contains a dominant frequency fpeak and subdominant structures associated with the dynamics of the differentially rotating remnant. These features are potentially valuable probes of dence nuclear matter, but robustly inferring them is difficult because the signal is short, high-frequency, and multi-component.

The method developed in this thesis combines an analytic damped-oscillator model with empirical relations and a classification of post-merger spectra. Inspiral-measurable quantities, in particular the chirp mass Mchirp and the effective tidal deformability Λ˜, are used to construct informed priors for the main post-merger frequencies. The spectral classification determines which components should be emphasized and which priors should be broad. This structure reduces the multimodality that can arise when broad uninformative priors allow different oscillators to exchange roles.

The method was validated on fifteen numerical-relativity waveforms of binary neutron-star mergers forming differentially rotating remnants. The set includes thirteen waveforms from the CoRe catalog and two additional waveforms, mostly equal-mass systems and two cases with mass ratio q = 1.5. The models span neutron-star radii of approximately 11.5 km to 16 km and masses from approximately 1.2 M⊙ to 1.7 M⊙, covering the relevant post-merger spectral types considered in the analysis.

The resulting fitting factors have median values in the range 0.935–0.980. At post-merger network signal-to-noise ratio 50, the typical one-sigma uncertainty in the dominant frequency fpeak is about 0.22%. For the dominant secondary peak it is about 1.4%, while for weaker secondary or post-peak frequencies it is about 6.2%. The informed-prior strategy gives converged posterior distributions in cases where less constrained analyses can produce ambiguous or multimodal results. Replacing dynesty with pocoMC yields nearly an order-of-magnitude speedup in wall time, making it possible to explore more cases and noise realizations.

The broader significance of this result is that post-merger inference should not be treated as a blind fit of a flexible model. The inspiral and post-merger phases are physically connected, and the inference method should use that connection. Future work should test the robustness of the approach against additional physical effects that are absent or simplified in the waveform set considered here, including magnetic fields, neutrino transport, turbulence modeling, bulk viscosity, possible dark-matter imprints, and deviations from general relativity. It will also be important to compare analytic oscillator models with reduced-order or machine-learning representations of the post-merger spectrum.

Synthesis and outlook
The four parts of the thesis are connected by a common principle: numerical models should preserve the structures that make the continuum problem well defined. In the SBP study, the relevant structure is the energy estimate; preserving it discretely yields stable high-order evolution through a coordinate singularity. In the self-force study, the relevant structure is the scale hierarchy; the small object is unresolved and should be modeled at finite resolution rather than as a singularity that the discretization cannot represent. In the horizon study, the relevant structure is geometric and topological; the cusp carries a distributional curvature contribution that controls jumps in quasi-local quantities. In the post-merger study, the relevant structure is physical information shared between inspiral and post-merger observables; using that information in the priors turns an ambiguous inference problem into a robust one.

The comparable-mass direction opened by this work is to connect source-side diagnostics and detector-side observables more directly. For black-hole mergers, this means extending the cusp and multipole analysis to generic spinning binaries and investigating whether horizon multipoles can be quantitatively related to waveform modes and recoil. For binary neutron stars, it means applying robust post-merger inference to waveforms with increasingly realistic microphysics and to detector networks expected in the third-generation era.

The extreme-mass-ratio direction is to move from scalar model problems to gravitational self-force calculations in curved spacetime. The finite-resolution viewpoint developed here should be tested in settings where long-time phase accuracy can be measured. It must also be combined with numerical methods capable of handling localized effective sources in multiple dimensions.

The numerical-methods direction is to incorporate the high-order curvilinear SBP operators into more complete evolution systems. This includes coupling to nonlinear equations, matter sources, gauge dynamics, and realistic boundary conditions. If successful, such operators would allow simulations to use symmetry-adapted coordinates while retaining the stability guarantees of an energy method.

Taken together, the results of this thesis advance compact-binary modeling across the mass-ratio spectrum. They do so by combining quasi-local geometry, gravitational-wave inference, finite-resolution source modeling, and stable high-order discretization. The unifying lesson is that robust numerical relativity is not only a matter of increasing resolution or computational cost. It is a matter of formulating the discrete problem so that the essential physics and mathematics of the continuum problem remain visible in the computation.