Tutorial: AsymptoticBondiData

The scri.asymptotic_bondi_data.AsymptoticBondiData() class (or ABD for short) is a convenient way to store and work with the standard set of asymptotic waveform quantities: the five Weyl scalars \((\Psi_4, \Psi_3, \Psi_2, \Psi_1, \Psi_0)\) and the Newman-Penrose shear \(\sigma\). These waveform quantities can be expanded in a spin-weighted spherical harmonic (SWSH) basis, and the coefficients of each mode will be stored in the ABD object.

The ABD class makes it simple to compute BMS charges and fluxes, and it serves as a vehicle for elaborate calculations invovling asymptotic data. See the class documentation for more details on all it has to offer. This tutorial is by no means an exhaustive overview of the class.

Creating an AsymptoticBondiData Object

To create an ABD object, we first need to define the time series that all the waveform quantities share, the highest value of \(\ell\) in the SWSH expansion, and how we will treat values of \(\ell\) larger than ell_max during operations:

>>> abd = scri.asymptotic_bondi_data.AsymptoticBondiData(
...     time = np.linspace(0,10,11),
...     ell_max = 8,
...     multiplication_truncator = max,
... )

When multiplying two fields, the SWSH expansion of the product will require a higher ell_max. Setting multiplication_truncator=sum will choose the ell_max of the product of two fields to be the sum of ell_max for each field. Here we defined each waveform quantity in abd to have ell_max=8, which means that the product \(\sigma\Psi_4\) will be represented by ell_max = 16. This means that the total modes stored for ell_max = 16 is 289 modes, as opposed to 81 modes for ell_max = 8. Although it is more correct to keep all these modes, given the precision of numerical relativity (NR) waveforms it is often sufficient to set multiplication_truncator=max as we have done here. This will instead choose the ell_max of the product of the two fields to be the max of the values of ell_max for each field. In this case, the product \(\sigma\Psi_4\) will be represented by ell_max=8.

Now that abd is set up, we can proceed to add some asymptotic data to it. The data should be a numpy.ndarray of shape (n_times, (ell_max + 1)**2) , where n_times is the size of the time array passed into the definition of abd above. This array must contain the coefficients of the SWSH modes starting from \(\ell=0\), even these modes are all zeros. The \((\ell,m)\) modes must also be in the following order:

(0,0), (1,-1), (1,0), (1,1), (2,-2), (2,1), (2,0),…

If we have an array psi4_mode_data of coefficients for \(\Psi_4\), then we can load this into abd by:

>>> abd.psi4 = psi4_mode_data

You can load the data for abd.sigma, abd.psi3, abd.psi2, abd.psi1, and abd.psi0 in the same way. There are two crucial points to note for numerical relativists. First, the ABD class uses the Moreschi-Boyle convention, which is not the standard convention used in NR. See Appendices B and C of arXiv:2010.15200 for more information. Second, the Newman-Penrose shear \(\sigma\) is not the gravitational-wave strain \(h\). In the Moreschi-Boyle convention, we have the relation \(h = 2\bar{\sigma}\), but this relation is not gauranteed in every convention and only holds for asymptotic quantities.

To make the ABD interface more nicely with NR waveforms, there is a tool for directly creating an ABD object from the HDF5 files of the SXS waveform catalog or any other catalog using the same file format:

>>> abd = scri.SpEC.create_abd_from_h5(
...     "SXS",
...     h    = "/path/to/rhOverM_file.h5",
...     Psi4 = "/path/to/rMPsi4_file.h5",
...     Psi3 = "/path/to/r2Psi3_file.h5",
...     Psi2 = "/path/to/r3Psi2OverM_file.h5",
...     Psi1 = "/path/to/r4Psi1OverM2_file.h5",
...     Psi0 = "/path/to/r5Psi0OverM3_file.h5",
... )

The first argument specifies three different ways that the HDF5 files can be internally structured. We currently support SXS for extrapolated asymptotic files in the the SXS catalog, CCE for asymptotic files produced by SpECTRE CCE, and RPXMB for compressed waveform files. See the documentation for scri.SpEC.file_io.create_abd_from_h5() for more information about these formats. This loader will convert the data into the Moreschi-Boyle convention and set \(\sigma = \bar{h}\) for abd.sigma.

Calculations with ModesTimeSeries

With a fully loaded ABD in hand, let’s do some computations! The time array can be accessed by:

>>> abd.t   # or abd.u

and the data for any individual quantity can be accessed by:

>>> abd.sigma
>>> abd.psi4
>>> abd.psi3
>>> abd.psi2
>>> abd.psi1
>>> abd.psi0

Individual modes of abd.psi4 (for example) can be accessed by the abd.psi4.index function. Alternatively, you can use the spherical_functions.LM_index function from the spherical_functions module. This can be aliased to lm for convenience, as done below. The third argument of LM_index is ell_min, but we always have ell_min=0 for ABD.

>>> # Get the (2,1) mode of Psi4
>>> l, m = 2, 1
>>> abd.psi4[:, abd.psi4.index(l,m)]

>>> # Alternatively:
>>> from spherical_functions import LM_index as lm
>>> abd.psi4[:, lm(l,m,0)]

The data for each quantity is stored as a scri.modes_time_series.ModesTimeSeries():

>>> type(abd.sigma)
<class 'scri.modes_time_series.ModesTimeSeries'>

There are many built-in functions that can be performed with a ModesTimeSeries, and multiple operations can be composed easily. Multiple operations will be performed in order from left to right.

>>> # take a derivative or two
>>> psi4_dot  = abd.psi4.dot
>>> psi4_ddot = abd.psi4.ddot

>>> # Integrate once or twice
>>> psi4_int  = abd.psi4.int
>>> psi4_iint = abd.psi4.iint

>>> # Apply the eth or ethbar operator (using the GHP definition,
>>> # which is the one that is natural to the Moreschi-Boyle convention)
>>> eth_psi4    = abd.psi4.eth_GHP
>>> ethbar_psi4 = abd.psi4.ethbar_GHP

>>> # If you really want to use the Newman-Penrose eth operators then you can do:
>>> eth_psi4    = abd.psi4.eth
>>> ethbar_psi4 = abd.psi4.ethbar

>>> # Get fancy and combine them together
>>> abd.sigma.dot.eth_GHP.eth_GHP

These operations will respect and update the spin-weight accordingly:

>>> # Check the spin-weight
>>> abd.sigma.s
2
>>> abd.sigma.bar.s
-2
>>> abd.sigma.dot.eth_GHP.eth_GHP.s
4

We can add ABD quantities together, but only in a way that makes sense. An error will be thrown for adding quantities of different spin weights:

>>> # This will throw an error
>>> abd.psi4 + abd.psi3

>>> # This will work!
>>> abd.psi4.eth_GHP + abd.psi3

There are two ways to perform multiplication with ModesTimeSeries quantities. Here we do not mean multiplying the mode coefficients together, but properly multiplying the fields and return the mode coefficients of the product. The more straightforward way to perform a multiplication is:

>>> sigma_psi4 = abd.sigma * abd.psi4

This will combine the mode coefficients with Wigner-3j symbols to compute the resulting mode coefficients of the product. The benefits of this approach are that it is straightforward to code and it does not suffer from aliasing effects. The downside is that it take a long time to run, so make sure you store the result as a variable if you will be using it more than once!

The second way to perform a multiplication is:

>>> sigma_psi4 = abd.sigma.grid_multiply(abd.psi4)

This will convert abd.sigma and abd.psi from a spectral representation to a physical-space representation by evaluating the fields on a grid of points on \(\mathscr{I}^+\). The mutliplication is performed pointwise and then transformed back to a spectral representation. This approach is much faster. However, if the value of ell_max is not high enough then aliasing effects might arise. See the documentation on scri.modes_time_series.ModesTimeSeries.grid_multiply() for options to adjust ell_max during the operation.

This just scratches the surface of all you can do with the ABD class. See the documentation on scri.asymptotic_bondi_data.AsymptoticBondiData() and scri.modes_time_series.ModesTimeSeries() to explore more functions.

Caveat

When taking the real part, imaginary part, or the complex conjugate, be very careful to know whether you are acting on the quantity as a field or just the modes of the quantity. For example:

>>> # This takes the real part of the quantity Psi2, and then
>>> # returns the modes of Re(Psi2). The mode weights are still complex!
>>> abd.psi2.real

>>> # This returns the real part of the mode weights of Psi2.
>>> # This is an array of real numbers!
>>> abd.psi2.ndarray.real

>>> # This returns the modes of the complex conjugate of sigma.
>>> # This is usually what you want.
>>> abd.sigma.bar

>>> # This returns the complex conjguate of the modes of sigma.
>>> # This is usually NOT what you want.
>>> np.conjugate(abd.sigma.ndarray)

BMS Transformations

Spacetime translations, supertranslations, frame rotations, and boosts can all be performed with the scri.asymptotic_bondi_data.AsymptoticBondiData.transform() function. See the documentation of the function for details. All the quantities in the ABD object will be transformed together. The transformation is not performed in place, so it will return a new ABD with the transformed data:

>>> abd_prime = abd.transform(
...     space_translation=[-1., 4., 0.2],
...     boost_velocity=[0., 0., 1e-2],
... )

These components of a BMS transformation can also all be stored in the scri.bms_transformations.BMSTransformation() class, which allows for things like reording the components of the BMS transformation, inverting BMS transformations, and composing BMS transformations. For more, see Tutorial: BMSTransformation/LorentzTransformation.

BMS Frames

All waveforms at future null infinity (and all waveforms more generally) are functions of coordinates. Therefore, there are certain “frames” which may be more useful than others, like that of a rest frame. For waveforms at future null infinity, the number of coordinate freedoms, i.e., the symmetries, that they exhibit is infinite and is summarized by a group known as the BMS group. This controls the types of frames that one may map waveforms to. Because GR is covariant, there is no preferred frame. However, for performing analysis on waveforms or building waveform models, it turns out that there are certain frames that are more useful than others. In particular, within GR one can extend the notion of a rest frame to something called a “superrest frame” (see arXiv:2405.08868 or arXiv:2208.04356 for more details), which typically yields waveforms that are easier to understand/analyze. Effectively, mapping to this frame amounts to mapping the system to be in the center-of-mass frame, with no instananeous memory, and its angular velocity in the z-direction. For example, for a remnant black hole, this corresponds to making the coordinates match those of the usual Kerr metric and is therefore incredibly useful (and necessary) for fitting QNMs to NR waveforms. There is another choice of frame which we call the Post-Newtonian BMS frame (or PNBMS frame). In this frame, the waveform will agree with the Post-Newtonian description of the system at early times.

The function scri.asymptotic_bondi_data.map_to_superrest_frame.map_to_superrest_frame() can map to one of these frames. In particular, it takes as input:

• t_0, the time at which to map to the superrest frame;

• target_PsiM_input, the target Moreschi supermomentum; this should be None to map to the superrest frame, but to map to the PN BMS frame one should input the PN Moreschi supermomentum (see arXiv:2208.04356).

• target_strain_input, the target strain; this should be None to map to the superrest frame, but to map to the PN BMS frame one should input the PN strain (see arXiv:2208.04356).

• padding_time, the time window about t_0 to be used when finding the BMS transformation to the superrest frame.

Loading CCE data and adjusting the BMS frame

For processing the output of SpECTRE CCE, one may use the function scri.SpEC.file_io.create_abd_from_h5. This function takes as input the path to SpECTRE CCE’s output file (via the option file_name) and creates an abd object from said file. It also can perform a number of other important post-processing steps, such as:

• time translate the time array of the waveforms by the radius of the worldtube; this ensures that the CCE waveforms are more closely aligned (in time) with extrapolated waveform. This is performed via the radius option.

• scale out the total Christoudoulou mass of the system from each waveform. This is performed via the ch_mass option.

• interpolate to a coarser time array, such as the time array of the worldtube. This is performed via the t_interpolate option.

• map to the superrest BMS frame at some time. This is performed via the t_0_superrest and padding_time options. E.g., to make reasonable-looking waveforms, one should map to the superrest frame at some time after junk radiation; t_0_superrest is the time at which to map to this frame, and padding_time is the window around t_0_superrest that is used when computing this BMS transformation. t_0_superrest - padding_time should be after junk radiation. padding_time should be a few hundred \(h=2\overline{\sigma}\) (where \(\sigma\) is the shear), e.g., two orbital periods. The function used to do this is abd.map_to_superrest_frame (see the “BMS Frames” section).

We recommend including all of these post-processing steps when processing SpECTRE CCE output.

To obtain the strain \(h\) from the abd object, one can use the function scri.asymptotic_bondi_data.map_to_superrest_frame.MT_to_WM() via h = MT_to_WM(2.0*abd.sigma.bar). This is because the strain \(h\) is related to the shear \(\sigma\) via \(h=2\overline{\sigma}\).

Example usage of this function could be:

>>> import scri
>>> import matplotlib.pyplot as plt
>>> abd = scri.create_abd_from_h5(
      file_name="CharacteristicExtractVolume_R0292.h5",
      file_format="spectrecce_v1",
      # ch_mass=1.0,  # Optional; helpful if known
      # t_interpolate=t_worldtube,  # Optional; for some specified values of `t_worldtube`
      t_0_superrest=1600,
      padding_time=200
    )
>>> h = abd.h
>>> plt.plot(h.t, h.data[:, h.index(2,2)])
>>> plt.show()

For more on the scri.WaveformModes() class, i.e., what \(h\) is in the above code, see https://github.com/moble/scri/blob/main/docs/tutorial_waveformmodes.rst.

The CCE output is sometimes very densely sampled in time, which causes the BMS transformation to the superrest frame or PNBMS frame to take a long time (more than 10 minutes). An effective way around this is to down sample the CCE data, compute the BMS transformation to the appropriate frame with the down sampled data, and then apply to the full data.

Fixing the superrest frame

The superrest frame of the CCE waveforms can be fixed using the following two functions.

def load_waveform(filename: str, dt: float = 1.0) -> scri.AsymptoticBondiData:
    abd = scri.SpEC.file_io.create_abd_from_h5("SpECTRECCE_v1", file_name=filename)
    cce_times = np.arange(abd.time[0], abd.time[-1], step=dt)
    return scri.SpEC.file_io.create_abd_from_h5("SpECTRECCE_v1", file_name=filename, t_interpolate=cce_times)

def load_waveform_and_fix_to_superrest(filename: str, dt: float = 1.0,
                                       t0_supperrest: float = 1800.0,
                                       superrest_padding: float = 200.0,
                                       superrest_dt: float = 5.0) -> scri.AsymptoticBondiData:
    abd = load_waveform(filename, dt) # Load full data
    # Interpolate to some lower-resolution in time ABD. This is to speed up the BMS optimization
    # which can take an hour or two if you don't down-sample in time.
    abd_for_bms = abd.interpolate(np.arange(t0_supperrest - 1.1 * superrest_padding,
                                            t0_supperrest + 1.1 * superrest_padding,
                                            step=superrest_dt))
    # Compute the BMS transformation to the superrest frame on the low-resolution ABD.
    _, BMS, _ = abd_for_bms.map_to_superrest_frame(t_0=t0_supperrest, padding_time=superrest_padding)
    # Transform the full data with the BMS transformation from the low-resolution computation. There is
    # minimal to no drawback from this.
    return abd.transform(supertranslation=BMS.supertranslation,
                         frame_rotation=BMS.frame_rotation.components,
                         boost_velocity=BMS.boost_velocity,)

Fixing the PNBMS frame using PN CoM charge

The PNBMS frame fixing of CCE waveforms requires target_PsiM_input (the PN Moreschi supermomentum), and target_strain_input (the PN strain). These post-Newtonian waveform modes can be generated using the PostNewtonian.jl module. We can generate these waveform modes and provide them as input arguments to the map_to_superrest_frame function for fixing the PNBMS frame. The generation of these waveform modes and fixing the frame can be performed using the following functions.

from sxs.julia import PNWaveform
import numpy as np
import sxs
import scri
from scri.pn.boosted_comcharge import PN_charges, analytical_CoM_func

def generate_h_and_Psi_M_from_abd(abd, params, t0_pnbms, pnbms_padding):
    """This function computes the target strain and target PsiM waveform
    required for PNBMS frame fixing.

    Parameters
    ----------
    abd : AsymptoticBondiData
    AsymptoticBondiData object from which the PN waveforms will be computed.
    params: tuple
    List of intrinsic parameters corresponding to the simulation. Its
    content should follow the form:
    params = M1, M2, chi1, chi2
    where M1, M2 are floats,
    chi1 : list of floats (n,3)
    chi2 : list of floats (n,3)
    t0_pnbms : float
    When to map to the PNBMS frame.
    pnbms_padding : float
    Amount by which to pad around t0 to speed up computations.
    This also determines the range over which certain BMS charges will be
    computed.

    Returns
    -------
    h_pn_scri, Psi_M_scri : scri.WaveformModes object
    """
    # First we create a window for frame-fixing using t0 and padding time.
    # The window is padded by an additional 200M on both sides of t0.
    t_left = t0_pnbms - (pnbms_padding + 200)
    t_right = t0_pnbms + (pnbms_padding + 200)

    left_idx = np.abs(abd.t - t_left).argmin()
    right_idx = np.abs(abd.t - t_right).argmin()

    h = abd.h[left_idx : right_idx]
    t = h.t - h.t[0]

    N = h.copy()
    N.dataType = scri.hdot
    N.data = h.data_dot
    ω = np.linalg.norm(N.angular_velocity(), axis=1)[0]

    # The convention for choosing the phase from h_{21} mode and the factor of
    # pi/2 is described in <https://arxiv.org/abs/2603.24661>. This choice
    # does not work for q=1.
    h_21 = h.data[:, h.index(2, 1)]
    θ = -1 * np.unwrap(np.angle(-h_21))[0] + np.pi/2

    R_i = np.array([np.cos(θ / 2), 0, 0, np.sin(θ / 2)]).T

    M1, M2, chi1, chi2 = params

    h_PN = PNWaveform(M1, M2, chi1, chi2, ω, modes_function=sxs.julia.h_bang,
        R_i=R_i, saveat=t)
    Psi_M_PN = PNWaveform(M1, M2, chi1, chi2, ω,
        modes_function=sxs.julia.Ψ_M_bang, R_i=R_i, saveat=t)

    # The PNWaveform function needs the kwarg times to start from t=0. Thus, we
    # translate the time array of PN waveform modes to match with the inspiral
    # window constructed from the abd object at the end.
    h_PN.t += t0_pnbms - (pnbms_padding + 200)
    Psi_M_PN.t += t0_pnbms - (pnbms_padding + 200)

    # The sxs.WaveformModes objects should be converted to scri.WaveformModes
    # before passing to map_to_superrest_frame.
    h_pn_scri = scri.WaveformModes.from_sxs(h_PN)
    Psi_M_scri = scri.WaveformModes.from_sxs(Psi_M_PN)

    h_pn_scri.data *= -1

    return h_pn_scri, Psi_M_scri

def load_waveform_and_fix_to_pnbms(
    filename: str,
    params,
    dt: float = 1.0,
    t0_pnbms: float = 1800.0,
    pnbms_padding: float = 200.0,
    ) -> scri.AsymptoticBondiData:
    abd = load_waveform(filename, dt)
    # params should take the same signature as it does in the previous
    # function.
    M1, M2, chi1, chi2 = params
    M = M1 + M2
    ν = M1 * M2 / M**2

    # We use the parameters of the simulation, and the window information
    # (t0_pnbms and pnbms_padding) to generate PN strain and Moreschi
    # supermomentum modes data during that window.
    h_PN, Psi_M_PN = generate_h_and_Psi_M_from_abd(abd, params, t0_pnbms, pnbms_padding)

    # Compute the BMS transformation to the PNBMS frame and the transformed
    # abd object.
    abd_prime, BMS, _ = abd.map_to_superrest_frame(
        t_0=t0_pnbms,
        padding_time=pnbms_padding,
        target_PsiM_input=Psi_M_PN,
        target_strain_input=h_PN,
        Gfun=analytical_CoM_func,
        Gparams0=np.zeros(8),
        Gargsfun=PN_charges(M, ν),
    )

    return abd_prime

This is the recent method for fixing the PNBMS frame described in Khairnar et al. https://arxiv.org/abs/2603.24661. The function map_to_superrest_frame takes 3 additional keyword arguments - Gfun, Gparams0 and Gargsfun. Gfun is the analtyical function used for fitting the CoM charge, Gparams0 is the initial guess for the fitting parameters, and Gargsfun are additional arguments passed to Gfun. See the documentation of scri.asymptotic_bondi_data.map_to_superrest_frame.com_transformation_to_map_to_superrest_frame() in scri.asymptotic_bondi_data.map_to_superrest_frame() for additional details on the signature of Gfun and Gargsfun.

Note that we have used \(h_{21}\) mode to define the initial phase for the PN waveform. This choice won’t work for symmetric mass ratio systems because of the degeneracy in the choice of phase angle by \(\pi\). Therefore, we will have to carefully choose the initial phase to break this degeneracy for such cases.

The parameter dt can be used to choose the time sampling for the waveform ABD object you get back, with default of 1M. Reducing the time sampling for the waveform ABD is also useful if the CCE data was written too frequently to work with.