Contents

1 Introduction

Apparent horizons (AHs) are crucial for characterizing black holes and excising their interiors in numerical relativity (NR) simulations [6]. Connecting observational signatures to black hole physical parameters necessitates robust methods for identifying and tracking the BHs themselves within NR simulations. This is primarily accomplished using apparent horizon (AH) finding algorithms. An AH corresponds to the outermost marginally outer trapped surface (MOTS) of a BH, defined as a surface on which the expansion of outgoing null geodesics vanishes [2, 3, 4].

The BHaHAHA library [1] introduces the first open-source, infrastructure-agnostic library for AH finding in NR. It implements a hyperbolic flow-based approach, recasting the elliptic PDE for a MOTS as a damped nonlinear wave equation. This method preserves the full nonlinearity of the expansion equation \(\Theta =0\), avoiding the need for linearization or high-quality initial guesses, and is thus robust in practice.

The ET_BHaHAHA arrangement hosts the ET_BHaHAHA thorn, which serves as the Einstein Toolkit (ETK) [8] interface to the standalone BHaHAHA C library. Key responsibilities of ET_BHaHAHA include:

• Interpolating spacetime metric data (components of the 3-metric \(\gamma _{ij}\) and extrinsic curvature \(K_{ij}\)) from the host NR code’s grid (typically Cartesian and managed by Carpet [11, 12]) onto a spherical grid suitable for BHaHAHA.

• Managing persistent data for each horizon, such as previous horizon shapes, centroid locations, and search radii. This data is crucial for efficient tracking and for making good initial guesses for subsequent horizon finds.

• Orchestrating the horizon finding process on distributed memory architectures, handling MPI synchronization of persistent data. Due to ETK’s Fortran legacy, C structures cannot be directly synchronized via MPI. ET_BHaHAHA circumvents this by packing scalar horizon data into global CCTK_REAL variables, indexed by an enumeration, and then using Cactus reduction operations for synchronization. Horizon shape data (arrays) are also synchronized via reductions.

• Providing user-configurable parameters to control the behavior of the BHaHAHA solver and the ET_BHaHAHA interface.

• Outputting diagnostic information, similar in style to AHFinderDirect, such as horizon area, irreducible mass, centroid coordinates, and spin estimates derived from proper circumferences.

The BHaHAHA library itself incorporates several features to enhance performance and robustness, such as a multigrid-inspired refinement strategy, an over-relaxation technique, and OpenMP parallelization for the core computations. These features are leveraged by ET_BHaHAHA. An experimental ”BBH mode” is also available for automatically triggering common horizon searches during binary black hole merger simulations.

The hyperbolic relaxation system, central to BHaHAHA, evolves the horizon shape \(h(\theta , \phi )\) and an auxiliary ”velocity” \(v(\theta , \phi )\) in pseudo-time \(\tau \):

where \(\eta \) is a damping coefficient and \(c\) is the relaxation wave speed (set to 1 in BHaHAHA). The system relaxes towards a steady state where \(\Theta = 0\), the condition for a MOTS.

This thorn is inspired by AHFinderDirect [4, 6] and aims to provide a similarly robust and efficient tool for the numerical relativity community.

2 What ET_BHaHAHA Needs

The ET_BHaHAHA thorn requires a 3+1 decomposed spacetime, providing components of the 3-metric \(\gamma _{ij}\) and extrinsic curvature \(K_{ij}\). These are typically provided by any ADM evolution thorn that populates the corresponding ADMBase grid functions.

ET_BHaHAHA inherits from and depends on several other Einstein Toolkit thorns and arrangements for its operation, summarized in Table 1.

The metric data is expected to be provided on a Cartesian grid, potentially with adaptive mesh refinement (AMR) as managed by Carpet. ET_BHaHAHA then interpolates this data onto its internal spherical grid.

3 Obtaining and Compiling ET_BHaHAHA

The ET_BHaHAHA thorn is part of the ET_BHaHAHA arrangement for the Einstein Toolkit. The core apparent horizon finding capabilities are provided by the BHaHAHA library, which is generated using NRPy+ [9, 10] and is included within the thorn’s source tree (src/BHaHAHA/). The thorn ET_BHaHAHA itself is written in C.

To compile ET_BHaHAHA, ensure that the ET_BHaHAHA arrangement is part of your Cactus checkout. The compilation process is standard for Einstein Toolkit thorns and uses the main Cactus make system. The README file within the thorn’s main directory provides instructions for regenerating the BHaHAHA library if needed. The BHaHAHA library components and the ETK interface code are compiled using gcc by default, as specified in the Makefiles (src/BHaHAHA/Makefile and src/make.code.defn). These files also indicate support for OpenMP (CFLAGS include -fopenmp -std=gnu99 -O2 -march=native), which is enabled if detected by the compiler. No special external libraries beyond a C compiler supporting C99 (GNU extensions are used, e.g. __typeof__) and OpenMP are strictly required by ET_BHaHAHA itself, as its core dependencies (like GSL for table utilities or ADMBase functionality) are handled by the Einstein Toolkit’s build system when those thorns are active.

4 ET_BHaHAHA Parameters

The behavior of the ET_BHaHAHA thorn is controlled by several parameters defined in its param.ccl file. These parameters can be broadly categorized. Parameters are referred to by their param.ccl names.

4.1 Overall Thorn Parameters

max_num_horizons

(Type: CCTK_INT, Default: 3)
Specifies the maximum number of apparent horizons that the thorn will attempt to find and track. This determines the size of various internal arrays. This parameter must be set to be at least as large as the number of horizons you intend to find. In BBH mode, it must be 3. The hardcoded limit is 100 (BHAHAHA_MAX_HORIZONS_SUPPORTED). Range: 0:100

find_every

(Type: CCTK_INT, Default: 1, Steerable: ALWAYS)
This integer parameter controls the execution frequency of the main horizon finding routine BHaHAHA_find_horizons. The finder will run when the current Cactus iteration cctk_iteration is an integer multiple of find_every. If set to 0, ET_BHaHAHA will be disabled globally. Individual horizon search cadences can override this via find_every_individual. Range: 0:*

bah_verbosity_level

(Type: CCTK_INT, Default: 1, Steerable: ALWAYS)
Controls the amount of diagnostic information printed to standard output during the horizon search.

• 0: Essential output only (e.g., major errors, final timing).

• 1: Physical summary of found horizons (default). This includes basic information about each horizon found, like its area and centroid, and BBH mode logging.

• 2: Full algorithmic details. This includes verbose information about the convergence process for each multigrid level, iteration counts, internal timings, MPI sync messages, and possibly debugging information.

Range: 0:2

max_Ntheta

(Type: CCTK_INT, Default: 32)
The maximum number of grid points used to sample the \(\theta \) (polar) direction on ET_BHaHAHA’s internal spherical surfaces. This corresponds to the highest resolution in the multigrid hierarchy if multigrid is enabled. It must match the last entry in bah_Ntheta_array_multigrid. Range: 2:128 (even values only)

max_Nphi

(Type: CCTK_INT, Default: 64)
The maximum number of grid points used to sample the \(\phi \) (azimuthal) direction on ET_BHaHAHA’s internal spherical surfaces. This corresponds to the highest resolution in the multigrid hierarchy. It must match the last entry in bah_Nphi_array_multigrid. Range: 2:256 (even values only)

bah_num_resolutions_multigrid

(Type: CCTK_INT, Default: 3)
Specifies the number of multigrid levels used in the horizon search. The actual resolutions for each level are defined in bah_Ntheta_array_multigrid and bah_Nphi_array_multigrid. Range: 1:16 (practical limit related to array sizes of bah_Ntheta_array_multigrid)

bah_Ntheta_array_multigrid[16]

(Type: CCTK_INT ARRAY, Default: -100 for all, Steerable: ALWAYS)
An array defining the number of \(\theta \) points for each multigrid level, starting from the coarsest. The number of entries actually used is determined by bah_num_resolutions_multigrid. Example: [8, 16, 32]. The last entry must equal max_Ntheta. Poison value -100 indicates an unset entry. Range (for each element used): 2: * (even values only, up to max_Ntheta)

bah_Nphi_array_multigrid[16]

(Type: CCTK_INT ARRAY, Default: -100 for all, Steerable: ALWAYS)
An array defining the number of \(\phi \) points for each multigrid level, starting from the coarsest. The number of entries used is determined by bah_num_resolutions_multigrid. Example: [16, 32, 64]. The last entry must equal max_Nphi. Poison value -100 indicates an unset entry. Range (for each element used): 2: * (even values only, up to max_Nphi)

bah_enable_eta_varying_alg_for_precision_common_horizon

(Type: CCTK_INT, Default: 0)
Enables a specialized algorithm where the damping parameter \(\eta \) varies to achieve high precision, particularly for common horizon finding. For general-purpose use, this should be disabled (0). Range: 0:1

sync_data_across_ranks_after_each_search__warning_slow

(Type: CCTK_INT, Default: 0)
If set to 1, all persistent horizon data (centers, radii, shapes) will be synchronized across all MPI ranks after each horizon search sequence using MPI reductions. This allows for restarting a simulation with a different number of MPI ranks while maintaining consistent horizon history. However, this can significantly slow down the simulation (potentially by a factor of \(>\)2x for ET_BHaHAHA) as all ranks must wait for the slowest horizon search. If 0 (default), only minimal synchronization required for the thorn’s logic (e.g., for BBH mode) is performed at the beginning of a find, and rank 0 might have stale data for horizons owned by other ranks for checkpointing purposes. BHaHAHA is generally robust enough to recover from a lost horizon if this is disabled. Range: 0:1

interpolator_name

(Type: CCTK_STRING, Default: "Hermite polynomial interpolation")
Specifies the name of the Cactus interpolator to be used for interpolating metric data from the simulation grid to ET_BHaHAHA’s internal spherical grid. ”Hermite polynomial interpolation” is highly recommended as it often leads to faster overall convergence in BHaHAHA for comparable accuracy. ”Lagrange polynomial interpolation” is an alternative. Range: Any valid Cactus interpolator name.

interpolator_pars

(Type: CCTK_STRING, Default: "order=3")
A string containing parameters for the chosen interpolator, passed to Util_TableCreateFromString(). For ”Hermite polynomial interpolation”, ”order=3” is recommended. For ”Lagrange polynomial interpolation”, ”order=4” is recommended. Range: Any string accepted by the selected interpolator’s parameter table creation.

4.2 Per-Horizon Parameters

These parameters are arrays, indexed from 0 to max_num_horizons-1.

find_every_individual[101]

(Type: CCTK_INT ARRAY, Default: -1 for all, Steerable: ALWAYS)
Allows overriding the global find_every for individual horizon i. If set to -1, uses global find_every. If \(\ge 1\), specifies the search cadence for this horizon. If 0, finding for this horizon is disabled. Range (for each element): -1, 0, 1:*

bah_max_search_radius[101]

(Type: CCTK_REAL ARRAY, Default: 1.5 for all, Steerable: ALWAYS)
The maximum coordinate radius of the spherical shell where ET_BHaHAHA will search for horizon i, relative to its current guessed center. This defines the outer boundary of the interpolation region for metric data. This parameter persists throughout the simulation for each horizon. Range (for each element): 0:* (must be positive)

bah_Nr_interp_max[101]

(Type: CCTK_INT ARRAY, Default: 48 for all, Steerable: ALWAYS)
Maximum number of radial points on ET_BHaHAHA’s internal 3D spherical grid for interpolating metric components for horizon i. The actual number of points used will be determined adaptively by bah_radial_grid_cell_centered_set_up based on the current search radii, but will not exceed this value. Recommended range is 48-64. Increasing this value mainly increases the cost of the initial interpolation from the ETK grid. Range (for each element): 4:*

bah_M_scale[101]

(Type: CCTK_REAL ARRAY, Default: 1.0 for all, Steerable: ALWAYS)
Characteristic mass scale for horizon i, in code units. Used to make tolerances and damping dimensionless. Use puncture bare masses or expected irreducible masses. Range (for each element): 0:* (must be positive)

bah_cfl_factor[101]

(Type: CCTK_REAL ARRAY, Default: 0.98 for all, Steerable: ALWAYS)
CFL factor for the hyperbolic relaxation solver for horizon i. This controls the pseudo-timestep. A value of 0.98 is robust. Range (for each element): 0:1.0 (must be positive, \(\le 1.05\) technically for SSPRK3 but \(\le 1.0\) is safer)

bah_max_iterations[101]

(Type: CCTK_INT ARRAY, Default: 10000 for all, Steerable: ALWAYS)
Maximum hyperbolic relaxation steps before the horizon search for horizon i is terminated, even if convergence criteria are not met. For typical tracking, 1000-10000 is sufficient. Range (for each element): 0:*

bah_Theta_L2_times_M_tolerance[101]

(Type: CCTK_REAL ARRAY, Default: 1e-2 for all, Steerable: ALWAYS)
Convergence criterion based on the area-weighted \(L_2\) norm of \(\Theta \times M_{scale}\) for horizon i. A value of \(10^{-2}\) is recommended to balance speed and accuracy for typical horizon tracking, often deferring to the \(L_{\infty }\) norm. Convergence requires both L2 and Linf criteria to be met. Range (for each element): 0:*

bah_Theta_Linf_times_M_tolerance[101]

(Type: CCTK_REAL ARRAY, Default: 1e-5 for all, Steerable: ALWAYS)
Convergence criterion based on the \(L_{\infty }\) norm (maximum absolute value) of \(\Theta \times M_{scale}\) for horizon i. A value of \(10^{-5}\) is recommended as a typical target for accuracy. Range (for each element): 0:*

bah_eta_damping_times_M[101]

(Type: CCTK_REAL ARRAY, Default: 7.0 for all, Steerable: ALWAYS)
The dimensionless damping parameter \(\eta \times M_{scale}\) for the hyperbolic relaxation for horizon i. A value of 7.0 is highly recommended. Range (for each element): 0:*

bah_KO_strength[101]

(Type: CCTK_REAL ARRAY, Default: 0.0 for all, Steerable: ALWAYS)
Strength of the Kreiss-Oliger dissipation applied during the hyperbolic relaxation for horizon i. This usually hinders convergence or accuracy for apparent horizon finding, so it is disabled (0.0) by default. Range (for each element): 0:*

bah_initial_grid_x_center[101]

(Type: CCTK_REAL ARRAY, Default: 0.0 for all, Steerable: ALWAYS)
The initial x-coordinate for the center of ET_BHaHAHA’s spherical search grid for horizon i. This is used at \(t=0\) or if a horizon is lost and needs a reset. Subsequently, the center is extrapolated. Range (for each element): *:*

bah_initial_grid_y_center[101]

(Type: CCTK_REAL ARRAY, Default: 0.0 for all, Steerable: ALWAYS)
The initial y-coordinate for the center of ET_BHaHAHA’s spherical search grid for horizon i. Range (for each element): *:*

bah_initial_grid_z_center[101]

(Type: CCTK_REAL ARRAY, Default: 0.0 for all, Steerable: ALWAYS)
The initial z-coordinate for the center of ET_BHaHAHA’s spherical search grid for horizon i. Range (for each element): *:*

4.3 BBH Mode Parameters

These parameters control the behavior when ET_BHaHAHA is run in its experimental Binary Black Hole (BBH) mode.

bah_BBH_mode_enable

(Type: CCTK_INT, Default: 0)
Enable (1) or disable (0) the BBH mode. When enabled, ET_BHaHAHA will manage the activation/deactivation of searches for two individual inspiralling black holes and a common horizon. Requires max_num_horizons = 3. Range: 0:1

bah_BBH_mode_inspiral_BH_idxs[2]

(Type: CCTK_INT ARRAY, Default: [0,0])
An array of two integers specifying the zero-offset horizon indices for the two inspiralling black holes. Example: [0,1]. These indices must be unique and different from bah_BBH_mode_common_horizon_idx. Range (for each element): 0:2 (Indices must be distinct and within \(0 \dots \text {\texttt {max\_num\_horizons}}-1\); typically \(0,1,2\) if \(\text {\texttt {max\_num\_horizons}}=3\))

bah_BBH_mode_common_horizon_idx

(Type: CCTK_INT, Default: 2)
The zero-offset horizon index for the common (merged) horizon. Example: 2. This index must be unique and different from those in bah_BBH_mode_inspiral_BH_idxs. Range: 0:2 (Index must be distinct and within \(0 \dots \text {\texttt {max\_num\_horizons}}-1\))

5 Provided Variables

The ET_BHaHAHA thorn declares several variables in its interface.ccl to store persistent horizon data. These variables are primarily used internally for tracking horizon evolution across time steps and for MPI synchronization, but advanced users or other thorns might find them useful. All variables are public. The index for horizon-specific variables (e.g., x_center_m1[h]) is 0-based in C, ranging from 0 to max_num_horizons-1.

bah_prev_horizons (Group)

(Type: CCTK_REAL ARRAY, Size: max_num_horizons*max_Ntheta*max_Nphi, Distribution: CONSTANT)
This group holds the shapes \(h(\theta , \phi )\) of previously found horizons. It is a single flat 1D array, but logically contains three sets of surfaces for each of the max_num_horizons. For a given horizon index hidx (0 to max_num_horizons-1), the shape data for \(h(\theta \_j, \phi \_k)\) is accessed with an offset: hidx * max_Ntheta * max_Nphi + itheta * max_Nphi + iphi (assuming \(\phi \) varies fastest, though the actual C indexing macro IDX2 uses itheta + NUM_THETA * iphi).

• prev_horizon_m1: Surface shape \(h(\theta , \phi )\) of the most recently found horizon (time \(t-\delta t\)).

• prev_horizon_m2: Surface shape \(h(\theta , \phi )\) of the second most recently found horizon (time \(t-2\delta t\)).

• prev_horizon_m3: Surface shape \(h(\theta , \phi )\) of the third most recently found horizon (time \(t-3\delta t\)).

These are used for extrapolating an initial guess for the current time step. All ranks maintain a copy, synchronized by BHaHAHA_sync_persistent_across_ranks.

bah_coord_info (Group)

(Type: CCTK_REAL SCALAR, Distribution: CONSTANT, Effectively arrays of size max_num_horizons)
Stores scalar information from previous horizon finds for each horizon. Each member is a Cactus scalar grid function that behaves as an array of size max_num_horizons.

• t_m1[h], t_m2[h], t_m3[h]: Simulation times of the last three successful finds for horizon h. Initialized to -1.0 if not found.

• r_min_m1[h], r_max_m1[h], ... r_min_m3[h], r_max_m3[h]: Min/max coordinate radii (centroid-based) for horizon h from the last three finds.

• x_center_m1[h], y_center_m1[h], z_center_m1[h], ... x_center_m3[h], y_center_m3[h], z_center_m3[h]: Cartesian coordinates of the horizon’s centroid for horizon h from the last three finds. These are global Cactus coordinates.

Synchronized by BHaHAHA_sync_persistent_across_ranks.

bah_guess (Group)

(Type: CCTK_INT SCALAR, Distribution: CONSTANT, Effectively an array of size max_num_horizons)

• use_fixed_radius_guess_on_full_sphere[h]: Flag (0 or 1) for horizon h. If 1, ET_BHaHAHA uses a simple spherical initial guess. If 0, it uses an extrapolated shape. Updated based on solver success.

Synchronized by BHaHAHA_sync_persistent_across_ranks.

bah_is_horizon_active (Group)

(Type: CCTK_INT SCALAR, Distribution: CONSTANT, Effectively an array of size max_num_horizons)

• bah_horizon_active[h]: Flag (0 or 1) for horizon h. If 1, a search is attempted. If 0, the search is skipped. Primarily managed by BBH mode.

Synchronized by BHaHAHA_sync_persistent_across_ranks.

6 Scheduling and MPI Synchronization

The ET_BHaHAHA thorn schedules its main driver function and manages storage for its interface.ccl variables as defined in schedule.ccl.

6.1 Storage

The following grid function groups are declared for storage:

• bah_prev_horizons

• bah_coord_info

• bah_guess

• bah_is_horizon_active

Memory for these variables is allocated by the Cactus framework.

6.2 Scheduled Functions

BHaHAHA_find_horizons

Scheduled at CCTK_ANALYSIS with options: global-early.
This is the main C language driver function for ET_BHaHAHA. It is responsible for:

1. Checking if a horizon search is due based on find_every and find_every_individual.

2. Performing initial setup at cctk_iteration == 0, including initializing persistent variables and BBH mode flags.

3. Synchronizing persistent data across MPI ranks using BHaHAHA_sync_persistent_across_ranks (see Section 6.3 below).

4. Populating the BHaHAHA parameter structure (bhahaha_params_and_data_struct) for each active horizon, including reading persistent data from Cactus GFs (via read_persistent_data_store_to_bhahaha_params_and_data).

5. Applying robustness improvements (e.g., adjusting CFL, multigrid resolutions for first finds).

6. Extrapolating initial guesses for horizon centers and search radii using bah_xyz_center_r_minmax.

7. Executing BBH mode logic for activating/deactivating common and individual horizon searches.

8. For each active horizon:

   • Setting up the radial interpolation grid using bah_radial_grid_cell_centered_set_up.

   • Calling BHaHAHA_interpolate_metric_data to interpolate ADM metric components from the ETK grid to the BHaHAHA spherical grid. This is done only on the MPI rank owning the horizon.

   • Invoking the core BHaHAHA solver (bah_find_horizon) on the owning MPI rank.

   • Processing results: if successful, updating persistent Cactus GFs using write_persistent_from_bhahaha_params_and_data and outputting diagnostics using bah_diagnostics_file_output. If failed, setting flags to ensure a robust guess next time.

9. Optionally, performing a final synchronization of all persistent data if sync_data_across_ranks_after_each_search__warning_slow is enabled.

The global-early option is chosen to be consistent with AHFinderDirect, aiming to find horizons early in the analysis phase.

6.3 MPI Synchronization Details

Persistent horizon data (previous shapes, centers, radii, active flags) needs to be consistent across all MPI ranks for robust tracking and BBH mode operation. ET_BHaHAHA employs a custom MPI synchronization strategy implemented in BHaHAHA_sync_persistent_across_ranks.c:

• Data Ownership: Each horizon is ”owned” by a specific MPI rank (typically \(h \pmod {\text {num\_ranks}}\)). This rank is responsible for performing the horizon find and updating its primary data.

• Packing/Unpacking and Reduction: Scalar data for each horizon (from groups bah_coord_info, bah_guess, bah_is_horizon_active) are packed by the owner rank into a local CCTK_REAL buffer. Non-owner ranks contribute zeros. A global sum reduction (CCTK_ReduceLocArrayToArray1D) broadcasts the owner’s data to all ranks. The BHaHAHA_PackedScalarFields enum (defined in ET_BHaHAHA.h) ensures consistent indexing into this buffer. Horizon shape arrays (from bah_prev_horizons) are also synchronized segment by segment using sum reductions.

This synchronization occurs at the beginning of BHaHAHA_find_horizons (if running in parallel) and optionally at the end if sync_data_across_ranks_after_each_search__warning_slow is enabled.

7 Physics Background and Numerical Methods

The ET_BHaHAHA thorn interfaces with the BHaHAHA library, which implements a specific set of methods for finding apparent horizons. The core physics and numerics are detailed in [1].

7.1 Apparent Horizons and MOTS

An apparent horizon is defined as the outermost marginally outer trapped surface (MOTS). A MOTS is a 2-surface where the expansion \(\Theta \) of outgoing null geodesics vanishes. In the ADM 3+1 formalism [5], where the line element is \begin {equation} ds^2 = -\alpha ^2 dt^2 + \gamma _{ij} (dx^i + \beta ^i dt) (dx^j + \beta ^j dt), \label {eq:adm_line_element} \end {equation} with \(\alpha \) the lapse, \(\beta ^i\) the shift, and \(\gamma _{ij}\) the spatial 3-metric, the expansion for a surface \(F(x^k)=r-h(\theta ,\phi )=0\) is [1]: \begin {equation} \Theta = \frac {1}{\lambda } (\partial _i F \partial _j \gamma ^{ij} + \gamma ^{ij} \partial _i \partial _j F) - \frac {1}{\lambda ^3} \partial _k \lambda \gamma ^{kj} \partial _j F + s^k \partial _k (\ln \sqrt {\gamma }) - K + s^i s^j K_{ij} = 0, \label {eq:MOTS_expansion_detail} \end {equation} where \(s\_i = \partial _i F / \lambda \) is the unit normal, \(\lambda = \sqrt {\gamma ^{kl}\partial _k F \partial _l F}\), \(\gamma = \det (\gamma _{ij})\), \(K_{ij}\) is the extrinsic curvature, and \(K = \gamma ^{ij}K_{ij}\).

7.2 Hyperbolic Relaxation

BHaHAHA recasts the elliptic PDE (2) into a hyperbolic system by introducing a pseudo-time \(\tau \) evolution for the surface shape \(h(\theta , \phi )\) and an auxiliary ”velocity” field \(v(\theta , \phi )\):

where \(\eta \) is a damping coefficient (units of \(1/M\), related to bah_eta_damping_times_M), and \(c\) is the relaxation wave speed (dimensionless, set to 1 in BHaHAHA). As \(\tau \rightarrow \infty \), the system is driven towards a steady state where \(\partial _\tau v = 0 \implies \Theta = 0\).

7.3 Numerical Implementation Details within BHaHAHA

The BHaHAHA library’s numerical implementation, leveraged by ET_BHaHAHA, involves several key steps:

7.3.1 Metric Data Handling and Interpolation

• Input Data and Transformation: ET_BHaHAHA interpolates the ADM 3-metric \(\gamma _{ij}\) and extrinsic curvature \(K_{ij}\) (Cartesian basis) from the host Cactus grid onto a 3D cell-centered spherical grid used by BHaHAHA (via BHaHAHA_interpolate_metric_data.c). Internally, BHaHAHA (in bah_numgrid__external_input_set_up.c) converts these to conformal BSSN-like variables in a spherical basis using a reference-metric formalism [13, 14, 15] to analytically handle coordinate singularities at the poles.

• Precomputation of Metric Derivatives: First spatial derivatives of these BSSN-like variables (both angular and radial) are precomputed on the 3D spherical BHaHAHA ”interpolation source” grid (interp_src_gfs in bah_numgrid__interp_src_set_up.c) using sixth-order finite differencing. This is done once per multigrid level. Angular derivatives are computed in bah_hDD_dD_and_W_dD_in_interp_src_grid_interior.c, and radial derivatives at boundaries in bah_apply_bcs_r_maxmin_partial_r_hDD_upwinding.c.

• Interpolation to Trial Surface: During the pseudo-time evolution, metric quantities and their derivatives on the current trial surface \(r=h(\tau , \theta , \phi )\) are obtained via 1D Lagrange interpolation along radial spokes from the precomputed 3D BHaHAHA grid data using bah_interpolation_1d_radial_spokes_on_3d_src_grid.c. This is a crucial performance optimization.

7.3.2 Initial Data for Hyperbolic Relaxation

The initial guess \(h(0, \theta , \phi )\) for the hyperbolic relaxation is determined by ET_BHaHAHA (within BHaHAHA_find_horizons.c before calling the library’s bah_initial_data.c):

• If no prior horizon information is available (e.g., first timestep, or after a find failure indicated by use_fixed_radius_guess_on_full_sphere=1), a simple spherical guess \(h = 0.8 \times R\_{\text {search\_max}}\) is used, with metric data interpolated onto a full spherical volume.

• Otherwise, \(h(0, \theta , \phi )\), the grid center \((x\_c, y\_c, z\_c)\), and search radii \((r\_{\min }, r\_{\max })\) are extrapolated (up to quadratic order via bah_xyz_center_r_minmax.c) from up to three previous successful finds. This allows the metric interpolation to be restricted to a potentially much smaller spherical shell.

The initial auxiliary field is set to \(v(0, \theta , \phi ) = \eta h(0, \theta , \phi )\) by bah_initial_data.c to start near the \(\partial _\tau h = 0\) condition.

7.3.3 Pseudo-Time Evolution

The hyperbolic system (Eqs. ??-??) is evolved using a third-order Strong Stability Preserving Runge-Kutta (SSPRK3) scheme [16] by bah_MoL_step_forward_in_time.c. The pseudo-timestep \(\Delta \tau \) is determined by a CFL condition (via bah_cfl_limited_timestep_based_on_h_equals_r.c) based on the wave speed \(c=1\) and the proper grid spacings on the evolving surface \(h\). Angular derivatives of \(h\) (via bah_rhs_eval.c) needed for \(\Theta \) are recomputed at each SSPRK3 substep.

7.3.4 Convergence and Stopping Conditions

The pseudo-time evolution stops when the dimensionless expansion \(M\_{\text {scale}}\Theta \) satisfies user-defined tolerances for both its \(L_\infty \) and area-weighted \(L_2\) norms (controlled by parameters bah_Theta_Linf_times_M_tolerance and bah_Theta_L2_times_M_tolerance), or when bah_max_iterations is reached. \(M\_{\text {scale}}\) is a user-provided characteristic mass for the horizon.

7.3.5 Convergence Acceleration Techniques

• Multigrid-Inspired Refinement: The system is first solved on a coarse angular grid. The solution is then prolonged (interpolated) to serve as the initial guess for a finer grid, and so on, up to the target resolution specified by bah_Ntheta_array_multigrid and bah_Nphi_array_multigrid.

• Over-Relaxation: Periodically, an over-relaxation step (bah_over_relaxation.c) attempts to accelerate convergence by extrapolating the horizon shape based on its recent evolution: \(h\_{\text {trial}}(\alpha ) = h\_c + \alpha (h\_c - h\_p)\), where \(h\_c\) is the current shape, \(h\_p\) is a prior shape, and \(\alpha > 1\) is an extrapolation factor. The \(h\_{\text {trial}}\) that minimizes \(\|M\_{\text {scale}}\Theta \|_{L_\infty }\) is adopted if it offers significant improvement.

7.4 Binary Black Hole (BBH) Mode in ET_BHaHAHA

When the parameter bah_BBH_mode_enable = 1, ET_BHaHAHA activates specialized logic for binary black hole merger simulations, typically requiring max_num_horizons = 3.

1. Initial Phase: Horizon searches are active for two individual inspiralling BHs (whose indices are given by bah_BBH_mode_inspiral_BH_idxs). The search for a common horizon (index bah_BBH_mode_common_horizon_idx) is initially inactive.

2. Common Horizon Activation: The thorn monitors the separation distance between the centroids of the two inspiral BHs and their last known maximum radii. If the condition \((\text {distance}(BH\_1, BH\_2) + r\_{\text {max},1} + r\_{\text {max},2}) \le (2.0 \times \text {\texttt {bah\_max\_search\_radius}}[\text {common\_idx}])\) is met, and both individual BHs have been successfully found previously, the common horizon search is activated. The initial guess for the common horizon’s center is set to the center-of-mass of the two inspiral BHs (weighted by their respective bah_M_scale values). A full sphere guess is forced for its first few attempts.

3. Deactivation of Individual BHs: Once the common horizon is consistently found (i.e., use_fixed_radius_guess_on_full_sphere[common_idx] == 0) and remains active, searches for the individual inspiral BHs are deactivated.

To improve robustness for initial common horizon finds, if the common horizon has not been found successfully for three consecutive attempts, ET_BHaHAHA adjusts parameters for the common horizon search: the CFL factor is reduced by 10%, a fixed-radius guess is enforced, and if at least two multigrid levels are used, the coarsest multigrid resolution is temporarily increased to match that of the second level, with the maximum iterations doubled.

7.5 Diagnostics Computation

Upon convergence, various physical diagnostics are computed by the BHaHAHA library from the found MOTS \(h(\theta , \phi )\), including:

• Area \(A\) and irreducible mass \(M\_{\text {irr}} = \sqrt {A / (16\pi )}\).

• Coordinate centroid \(x\_c^i\) (area-weighted).

• Minimum, maximum, and mean coordinate radii relative to the centroid.

• Proper circumferences (\(C\_{xy}, C\_{xz}, C\_{yz}\)) in coordinate planes through the centroid.

• Dimensionless spin estimates \(|a^i/M|\) from ratios of these circumferences, by inverting Kerr metric relations [17].

These are handled by functions such as bah_diagnostics_area_centroid_and_Theta_norms.c, bah_diagnostics_min_max_mean_radii_wrt_centroid.c, and bah_diagnostics_proper_circumferences.c.

8 Output Data

The ET_BHaHAHA thorn, through the BHaHAHA library function bah_diagnostics_file_output.c, produces two main types of ASCII output files per found horizon, stored in the directory specified by the Cactus parameter IO::out_dir. These files are designed to be gnuplot-friendly.

8.1 Horizon Diagnostics File

• Filename Convention: [IO::out_dir]/BHaHAHA_diagnostics.ah[H].gp, where [H] is the 1-based horizon number (from bhahaha_params_and_data.which_horizon).

• Content: This file contains a time series of diagnostic quantities for the specified horizon. The file is appended to at each successful find. A header is written if the file is new.

• Columns:

  1. Current simulation iteration (cctk_iteration)

  2. Current simulation time (cctk_time)

  3. x-coordinate of the horizon centroid (global grid coordinates)

  4. y-coordinate of the horizon centroid

  5. z-coordinate of the horizon centroid

  6. Minimum coordinate radius of the horizon surface, relative to its centroid

  7. Maximum coordinate radius of the horizon surface, relative to its centroid

  8. Mean coordinate radius of the horizon surface, relative to its centroid (area-weighted)

  9. Proper circumference in the xy-plane (at centroid’s z-height)

  10. Proper circumference in the xz-plane (at centroid’s y-height)

  11. Proper circumference in the yz-plane (at centroid’s x-height)

  12. Apparent horizon area

  13. Irreducible mass \(M\_{\text {irr}} = \sqrt {A / (16\pi )}\)

  14. \(M\_{\text {scale}} \times ||\Theta ||_{L_\infty }\) (max norm of expansion)

  15. \(M\_{\text {scale}} \times ||\Theta ||_{L_2}\) (L2 norm of expansion)

  16. Spin x-component estimate from \(C\_{xy}/C\_{yz}\) ratio (or NaN if not applicable)

  17. Spin x-component estimate from \(C\_{xz}/C\_{yz}\) ratio (or NaN)

  18. Spin y-component estimate from \(C\_{yz}/C\_{xz}\) ratio (or NaN)

  19. Spin y-component estimate from \(C\_{xy}/C\_{xz}\) ratio (or NaN)

  20. Spin z-component estimate from \(C\_{xz}/C\_{xy}\) ratio (or NaN)

  21. Spin z-component estimate from \(C\_{yz}/C\_{xy}\) ratio (or NaN)

8.2 Horizon Surface Shape File

• Filename Convention: [IO::out_dir]/h.t[ITER].ah[H].gp, where [ITER] is the 7-digit zero-padded Cactus iteration number, and [H] is the 1-based horizon number.

• Content: This file contains the Cartesian coordinates \((x,y,z)\) of points on the found apparent horizon surface. The coordinates are relative to the global grid origin (i.e., the input center for BHaHAHA plus the found shape \(h(\theta ,\phi )\) transformed to Cartesian). Data is organized for easy plotting with gnuplot’s splot command using set pm3d map. Theta rings are separated by blank lines.

• Columns:

  1. x-coordinate

  2. y-coordinate

  3. z-coordinate

9 Monitoring ET_BHaHAHA’s Status

The verbosity of ET_BHaHAHA is controlled by the bah_verbosity_level parameter:

• bah_verbosity_level = 0: Only essential messages are printed, primarily error messages via CCTK_VError or CCTK_VWarn, and a summary message upon successful completion of all horizon finds for an iteration (e.g., ET_BHaHAHA total elapsed time...).

• bah_verbosity_level = 1 (Default): Prints a summary for each horizon found by the BHaHAHA library after its final multigrid level. This includes information about final expansion norms, area, irreducible mass, centroid location (global), mean/min/max radii, and estimated spin components. Also includes timing information for the BHaHAHA search for each horizon. For BBH mode, it logs activation/deactivation of horizon searches.

• bah_verbosity_level = 2: Provides the most detailed output. This includes everything from level 1, plus:

  • Detailed convergence information from the BHaHAHA library for each multigrid resolution level (iteration number, norms of \(\Theta \), etc., printed every output_diagnostics_every_nn internal iterations of BHaHAHA).

  • Timings for various sub-stages within BHaHAHA_find_horizons and BHaHAHA_interpolate_metric_data, such as time spent in MPI synchronization or metric interpolation for each horizon.

  • Information printed by read_persistent_data_store_to_bhahaha_params_and_data when DEBUG_MODE is enabled in ET_BHaHAHA C code.

  • Pre-interpolation state information: guess coordinates, active flags, etc.

  This level is primarily intended for debugging or detailed performance analysis.

Error messages from BHaHAHA are converted to descriptive strings by bah_error_message() and reported via CCTK_VWarn or CCTK_VInfo, often indicating failure to find a horizon and the reason.

10 Visualization

The output files from ET_BHaHAHA are formatted for easy visualization, primarily with gnuplot.

• Diagnostic Time Series: Files named BHaHAHA_diagnostics.ah[H].gp contain time series data. For example, to plot the irreducible mass of horizon 1 over time:

           gnuplot> plot "BHaHAHA_diagnostics.ah1.gp" using 2:13 with lines title "M_{irr} H1"
           
  

  Refer to Section 8.1 for column definitions.

• Horizon Shapes: Files named h.t[ITER].ah[H].gp contain the (x,y,z) coordinates of the horizon surface at a specific iteration. These can be plotted as a 3D surface in gnuplot:

           gnuplot> splot "h.t0000000.ah1.gp" with lines title "AH1 at t=0"
           # For a color map view
           gnuplot> set view map
           gnuplot> splot "h.t0000000.ah1.gp" with pm3d title "AH1 at t=0"
           
  

  Other 3D visualization tools like VisIt or ParaView can also be used by importing the ASCII data.

11 Accuracy and Testing

The BHaHAHA library and its ET_BHaHAHA interface have undergone several tests to assess accuracy and robustness, as detailed in [1].

• Comparison with AHFinderDirect: For dynamical binary black hole spacetimes (e.g., GW150914-like inspiral and merger), ET_BHaHAHA with default tolerances (\(L_\infty (\Theta \cdot M\_{\text {scale}}) \sim 10^{-5}\)) achieves results (horizon area, centroid trajectories) consistent with AHFinderDirect. The dominant source of discrepancy at these tolerance levels is typically the error from interpolating metric data from the host NR code’s grid to the horizon finder’s internal grid (around \(10^{-5}\) for typical ETK setups). ET_BHaHAHA tends to be significantly faster (approx. 2.1x) than AHFinderDirect when achieving comparable (interpolation-limited) accuracy in such dynamical scenarios.

• Scale Invariance: BHaHAHA is designed with scale-invariant convergence criteria (tolerances are relative to bah_M_scale). Tests on Brill-Lindquist two black hole initial data, where the total mass was varied over 16 orders of magnitude while keeping dimensionless quantities fixed, showed that BHaHAHA maintained consistent accuracy in normalized horizon areas (\(A/M^2\)). In contrast, AHFinderDirect (with absolute tolerances) showed deviations at very small or very large mass scales.

• Precision Tests (Three-BH Critical Radius): For the demanding Thornburg three-black-hole ”critical radius” test [4], BHaHAHA (when run within its native BlackHoles@Home environment with high-order finite differencing/interpolation and strict tolerances, e.g., \(M_{scale}||\Theta ||_{L_\infty } < 3 \times 10^{-8}\)) reproduces \(R\_{\text {crit}} \approx 1.1954995\) to at least 8 significant digits, matching results from AHFinderDirect and other specialized solvers [18]. This demonstrates the intrinsic accuracy of the underlying hyperbolic relaxation algorithm.

The default tolerance for ET_BHaHAHA (bah_Theta_Linf_times_M_tolerance = 1e-5) is typically sufficient for most production NR simulations, as accuracy is often limited by metric interpolation error rather than the horizon finder’s intrinsic precision.

12 Examples

Example parameter files for using ET_BHaHAHA can be found in the ET_BHaHAHA/par/ directory.

12.1 Single Shifted Kerr Black Hole

The parameter file shifted_kerr_schild.par sets up a single, spinning Kerr black hole using the ShiftedKerrSchild initial data thorn. Key ET_BHaHAHA parameters in this file:

ActiveThorns = "ET_BHaHAHA"
ET_BHaHAHA::max_num_horizons = 1
ET_BHaHAHA::bah_verbosity_level = 2
ET_BHaHAHA::bah_max_search_radius[0] = 1.5
# BHaHAHA::bah_M_scale[0] = 1.0 (Matches ShiftedKerrSchild::BH_mass)

12.2 Binary Black Hole Inspiral and Merger (BBH Mode)

The parameter file GW150914_baikalvacuum.par simulates a GW150914-like binary black hole system using TwoPunctures initial data and evolved with BaikalVacuum. This example demonstrates the BBH mode of ET_BHaHAHA. Key ET_BHaHAHA parameters for BBH mode:

ActiveThorns = "ET_BHaHAHA"
BHaHAHA::max_num_horizons = 3
BHaHAHA::bah_BBH_mode_enable = 1
BHaHAHA::bah_BBH_mode_common_horizon_idx = 2
BHaHAHA::bah_BBH_mode_inspiral_BH_idxs = [0, 1]
BHaHAHA::find_every = 44 # Search cadence
# Per-horizon parameters (indices 0, 1 for inspiral; 2 for common)
BHaHAHA::bah_M_scale = [0.5538, 0.4461, 1.0]
BHaHAHA::bah_initial_grid_x_center[0] =  4.461538
BHaHAHA::bah_initial_grid_x_center[1] = -5.538461
BHaHAHA::bah_max_search_radius = [0.6, 0.5, 1.4]

This setup tracks the two inspiralling black holes. When they get close enough, ET_BHaHAHA automatically starts searching for the common horizon.

12.3 Brill-Lindquist Distorted Horizons

Files like brill_lindquist_highly_distorted_horizon.par (and variants with different mass scales, e.g., -1e-5.par, -1e8.par) set up two Brill-Lindquist punctures with a small separation, resulting in a highly distorted common apparent horizon. These test robustness and scale invariance. BBH mode is disabled, so three horizons (common, and two small individual ones near punctures) are sought independently.

13 Acknowledgements

The ET_BHaHAHA thorn was developed by Zachariah B. Etienne. The underlying BHaHAHA library, also primarily developed by Zachariah B. Etienne, benefitted from input and contributions from Thiago Assumpção, Leonardo Rosa Werneck, and Samuel D. Tootle, as detailed in [1]. Special thanks to Jonathan Thornburg, whose AHFinderDirect thorn has been a cornerstone of apparent horizon finding in the numerical relativity community and served as a template and high standard for BHaHAHA and ET_BHaHAHA. Thanks also to Roland Haas for valuable input during the development of the BHaHAHA library.

This research made use of the resources of the High Performance Computing Center at Idaho National Laboratory, which is supported by the Office of Nuclear Energy of the U.S. Department of Energy and the Nuclear Science User Facilities under Contract No. DE-AC07-05ID14517. In addition, it made use of the Falcon supercomputer, operated by the Idaho C3+3 Collaboration. Z.B.E. gratefully acknowledges support from NSF awards PHY-2110352/2508377, PHY-2409654, OAC-2004311/2227105, OAC-2411068, and AST-2108072/2227080, as well as NASA awards ISFM-80NSSC18K0538, TCAN-80NSSC18K1488, and ATP-80NSSC22K1898. TA is thankful for support from NSF grants OAC-2229652 and AST-2108269. L.R.W. gratefully acknowledges support from NASA award LPS-80NSSC24K0360. Z.B.E. and ST received additional support from NASA award ATP-80NSSC22K1898 and the University of Idaho P3-R1 Initiative.

14 License

The ET_BHaHAHA thorn is licensed under the BSD 2-Clause License:

BSD 2-Clause License

Copyright (c) 2025, Zachariah Etienne
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

* Redistributions of source code must retain the above copyright notice, this
  list of conditions and the following disclaimer.

* Redistributions in binary form must reproduce the above copyright notice,
  this list of conditions and the following disclaimer in the documentation
  and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

References

15 Parameters

16 Interfaces

General

Implements:

bhahaha

Inherits:

admbase

boundary

grid

Grid Variables

16.0.1 PUBLIC GROUPS

17 Schedule

This section lists all the variables which are assigned storage by thorn EinsteinAnalysis/ET_BHaHAHA. Storage can either last for the duration of the run (Always means that if this thorn is activated storage will be assigned, Conditional means that if this thorn is activated storage will be assigned for the duration of the run if some condition is met), or can be turned on for the duration of a schedule function.

Storage

Scheduled Functions

CCTK_ANALYSIS

  bhahaha_find_horizons

  driver function for bhahaha