We present the novel algorithmically regularised integration method MSTAR for high accuracy ($|\Delta E/E| \gtrsim 10^{-14}$) integrations of N-body systems using minimum spanning tree coordinates. The two-fold parallelisation of the $\mathcal{O}(N_\mathrm{part}^2)$ force loops and the substep divisions of the extrapolation method allows for a parallel scaling up to $N_\mathrm{CPU} = 0.2 \times N_\mathrm{part}$. The efficient parallel scaling of MSTAR makes the accurate integration of much larger particle numbers possible compared to the traditional algorithmic regularisation chain (AR-CHAIN) methods, e.g. $N_\mathrm{part} = 5000$ particles on $400$ CPUs for $1$ Gyr in a few weeks of wall-clock time. We present applications of MSTAR on few particle systems, studying the Kozai mechanism and N-body systems like star clusters with up to $N_\mathrm{part} =10^4$ particles. Combined with a tree or a fast multipole based integrator the high performance of MSTAR removes a major computational bottleneck in simulations with regularised subsystems. It will enable the next generation galactic-scale simulations with up to $10^9$ stellar particles (e.g. $m_\star = 100 M_\odot$ for a $M_\star = 10^{11} M_\odot$ galaxy) including accurate collisional dynamics in the vicinity of nuclear supermassive black holes.

We study the orbital evolution and gravitational wave (GW) emission of supermassive black hole (SMBH) binaries formed in gas-free mergers of massive early-type galaxies using the hybrid tree-regularized N-body code KETJU. The evolution of the SMBHs and the surrounding galaxies is followed self-consistently from the large-scale merger down to the final few orbits before the black holes coalesce. Post-Newtonian corrections are included up to PN3.5-level for the binary dynamics, and the GW calculations include the corresponding corrections up to PN1.0-level. We analyze the significance of the stellar environment on the evolution of the binary and the emitted GW signal during the final GW emission dominated phase of the binary hardening and inspiral. Our simulations are compared to semi-analytic models that have often been used for making predictions for the stochastic GW background emitted by SMBHs. We find that the commonly used semi-analytic parameter values produce large differences in merger timescales and eccentricity evolution, but result in only $\sim 10\%$ differences in the GW spectrum emitted by a single binary at frequencies $f\gtrsim 10^{-1} \, \rm yr^{-1}$, which are accessible by current pulsar timing arrays. These differences are in part caused by the strong effects of the SMBH binaries on the surrounding stellar population, which are not included in the semi-analytic models.

Interpolation of data represented in curvilinear coordinates and possibly having some non-trivial, typically Riemannian or semi-Riemannian geometry is an ubiquitous task in all of physics. In this work we present a covariant generalization of the barycentric coordinates and the barycentric interpolation method for Riemannian and semi-Riemannian spaces of arbitrary dimension. We show that our new method preserves the linear accuracy property of barycentric interpolation in a coordinate-invariant sense. In addition, we show how the method can be used to interpolate constrained quantities so that the given constraint is automatically respected. We showcase the method with two astrophysics related examples situated in the curved Kerr spacetime. The first problem is interpolating a locally constant vector field, in which case curvature effects are expected to be maximally important. The second example is a General Relativistic Magnetohydrodynamics simulation of a turbulent accretion flow around a black hole, wherein high intrinsic variability is expected to be at least as important as curvature effects.

Ray-tracing is a central tool for constructing mock observations of compact object emission and for comparing physical emission models with observations. We present Arcmancer, a publicly available general ray-tracing and tensor algebra library, written in C++ and providing a Python interface. Arcmancer supports Riemannian and semi-Riemannian spaces of any dimension and metric, and has novel features such as support for multiple simultaneous coordinate charts, embedded geometric shapes, local coordinate systems and automatic parallel propagation. The Arcmancer interface is extensively documented and user-friendly. While these capabilities make the library well suited for a large variety of problems in numerical geometry, the main focus of this paper is in general relativistic polarized radiative transfer. The accuracy of the code is demonstrated in several code tests and in a comparison with GRTRANS, an existing ray-tracing code. We then use the library in several scenarios as a way to showcase the wide applicability of the code. We study a thin variable-geometry accretion disk model, and find that polarization carries information of the inner disk opening angle. Next, we study rotating neutron stars and determine that to obtain polarized light curves at better than $\sim1\%$ level of accuracy, the rotation needs to be taken into account both in the space-time metric as well as in the shape of the star. Finally, we investigate the observational signatures of an accreting black hole lensed by an orbiting black hole. We find that these systems exhibit a characteristic asymmetric twin-peak profile both in flux and polarization properties.