Precessing waveforms

Full precessing waveforms from the aligned-spin mlgw_bns surrogate.

The surrogate of Model is aligned-spin: it predicts the frequency-domain multipoles \(\tilde{h}_{\ell m}(f)\) of a binary whose spins are parallel to the orbital angular momentum. A precessing binary, though, is well described — this is the observation underlying every “twisted” approximant, TEOBResumS included — by those same multipoles, taken to live in the co-precessing frame that follows \(\hat{L}(t)\), and then rotated into the inertial frame by three time-dependent Euler angles:

\[h^{\rm inertial}_{\ell m} = \sum_{n} D^{\ell}_{m n}(\alpha, \beta, \gamma) \, h^{\rm co-prec}_{\ell n} \,.\]

This module puts the two halves together: the Euler angles come from the PN spin-precession dynamics of twist_waveform (a reproduction of what TEOBResumS itself integrates), the co-precessing multipoles from coprecessing_modes_dict(), and the rotation from twist_modes().

Twisting in the frequency domain

The rotation above is a statement about the multipoles at a given time, whereas the surrogate lives in the frequency domain. In the stationary-phase approximation each multipole’s Fourier transform at a frequency \(f\) is dominated by the single time \(t\) at which the multipole’s instantaneous frequency equals \(f\); since the Euler angles vary on the (much longer) precession timescale, the twist carries over to the frequency domain by evaluating the angles at that stationary time. This is the standard frequency-domain twisting-up procedure, as used by e.g. IMRPhenomXPHM (arXiv:2004.06503).

Two consequences shape the implementation:

  • the stationary time differs from multipole to multipole. The co-precessing multipole \(n\) has phase \(\simeq n \Phi_{\rm orb}\), so at frequency \(f\) it is stationary where \(M \Omega_{\rm orb} = 2 \pi f M / n\). The angles are therefore looked up separately for each \(n\), against the orbital frequency that the PN integration provides alongside them.

  • with the \(\tilde{h}(f) = \int h(t) e^{2 \pi i f t} \mathrm{d}t\) convention used throughout mlgw_bns, the \(n > 0\) co-precessing multipoles have their support at \(f > 0\) and the \(n < 0\) ones at \(f < 0\). The two halves of the sum over \(n\) must then be twisted separately: the \(n > 0\) half builds \(\tilde{h}^{\rm inertial}_{\ell m}(f)\), the \(n < 0\) half builds \(\tilde{h}^{\rm inertial}_{\ell m}(-f)\). Both are needed, because precession breaks the equatorial symmetry that would otherwise let \(h_\times\) be recovered from \(h_+\) alone.

The whole construction reduces, identically, to predict() when the in-plane spin components vanish; check_aligned_spin_limit() asserts exactly that.

The orbital phase at the reference frequency

With precession the orbital phase is no longer a choice of observer: it sets the angle between the orbital separation and the in-plane spins. As in LALSuite (f_ref, phiRef, LIGO-T1500606), NRSur7dq4, SEOBNR and TEOBResumS (at the start of its integration), it is fixed at the frequency where the spins are given, reference_frequency_hz; the aligned-spin surrogate instead puts it to zero at the merger. The stationary-phase transform \(X_{\ell m}(f) = \Psi_{\ell m} - f \Psi'_{\ell m}\) of a frequency-domain multipole’s phase does not change under time shifts, and at \(f = m F\) equals

\[X_{\ell m}(m F) = m \phi_{\rm orb} + \delta_{\ell m} + \pi / 4 \,,\]

with \(\phi_{\rm orb}\) the orbital phase when the orbital frequency is \(F\), and \(\delta_{\ell m}\) the multipole’s phase at zero orbital phase. PrecessingModel.reference_orbital_phase() reads \(\phi_{\rm orb}\) off the surrogate’s own \((2, 2)\) at the reference frequency (modulo \(\pi\), the branch from the \((2, 1)\) or \((3, 3)\)), taking for \(\delta_{\ell m}\) its leading-order value (LEADING_ORDER_MODE_PHASES), and the multipoles are rotated so that it equals PrecessingParametersWithExtrinsic.reference_phase. Like NRSur7dq4’s, this orbital phase is defined by the waveform: TEOBResumS’ dynamical one differs from it by its multipoles’ higher-order phase corrections, ~5e-3 m rad at 10 Hz for a binary neutron star. The spins, on the other hand, are placed on the PN orbital-frequency track the Euler angles are integrated (and looked up) along, as in TEOBResumS; the two read the same instant slightly differently (TEOBResumS’ EOB and PN spin dynamics start at 9.506 and 9.5 Hz), which a comparison with it has to account for (PrecessingModel.teob_reference_phase()).

class PrecessingModel(model: Model)[source]

A precessing waveform model built on an aligned-spin Model.

Parameters:

model (Model) – The aligned-spin surrogate providing the co-precessing multipoles.

Examples

>>> from mlgw_bns.model import Model
>>> precessing = PrecessingModel(Model.default_for_testing())
>>> hp, hc = precessing.predict(frequencies, params)
coprecessing_amplitudes_and_phases(frequencies: ndarray, params: PrecessingParametersWithExtrinsic, source: str = 'surrogate') → Dict[Tuple[int, int], Tuple[ndarray, ndarray]][source]

The co-precessing multipoles, before the twist.

Those of coprecessing_amplitudes_and_phases() for params.aligned(), with the orbital phase at the reference frequency set to params.reference_phase if params.reference_frequency_hz is given (see the module docstring). predict_modes_dict() twists these.

Parameters:
Returns:

Mapping (l, m) -> (amplitude, phase).

Return type:

dict[tuple[int, int], tuple[np.ndarray, np.ndarray]]

euler_angles(params: PrecessingParametersWithExtrinsic, initial_frequency_hz: float, omega_dot: Callable[[float], float] | None = None, anchor_to_reference_phase: bool = False) → EulerAngles[source]

The precession angles for a source, from the PN dynamics.

See euler_angles(), of which this is the bound version.

Parameters:

  • params (PrecessingParametersWithExtrinsic) – Source parameters.

  • initial_frequency_hz (float) – Lowest frequency, in Hz, at which the waveform is wanted.

  • omega_dot (callable, optional) – Orbital-frequency rate to march the precession along; see euler_angles().

    Default: None

  • anchor_to_reference_phase (bool, optional) – If True (and omega_dot is not given), build omega_dot from this model’s own aligned-spin \((2, 2)\) phase for params – i.e. integrate the precession along the surrogate’s EOB-accurate orbital-frequency track rather than the PN flux. This is the surrogate analogue of TEOBResumS’ SPIN_FLX_EOB hand-off.

    Default: False

Returns:

The angles, tabulated against \(M \Omega_{\rm orb}\).

Return type:

EulerAngles

jax_predict(modes: Sequence | None = None, n_steps: int | None = None) → Callable[source]

A pure JAX function computing predict() for a batch of binaries.

See mlgw_bns.batched_precession.precessing_waveform() for its signature. Wrap it in jax.jit(); every new input shape compiles again. Requires the jax extra.

Parameters:

  • modes (sequence of (l, m), optional) – Co-precessing multipoles to twist; defaults to all of the model’s.

    Default: None

  • n_steps (int, optional) – Steps of each leg of the precession integration; defaults to mlgw_bns.batched_precession.N_STEPS.

    Default: None

jax_predict_modes(modes: Sequence | None = None, n_steps: int | None = None) → Callable[source]

A pure JAX function giving each co-precessing multipole and its twist, for a batch of binaries.

See mlgw_bns.batched_precession.precessing_mode_components(); the arguments are those of jax_predict().

property largest_mode_m: int

The largest \(m\) among the surrogate’s modes.

mode_components(frequencies: ndarray, params: PrecessingParametersWithExtrinsic, source: str = 'surrogate', angles: EulerAngles | None = None) → Dict[Tuple[int, int], Tuple[ndarray, ndarray, ndarray]][source]

Each co-precessing multipole and its twist coefficients.

The numpy counterpart of jax_predict_modes(): predict_numpy() is sum(c_plus * coprecessing), sum(c_cross * coprecessing) over the multipoles (see twist_coefficients()).

Parameters:
Returns:

Mapping (l, m) -> (coprecessing, c_plus, c_cross).

Return type:

dict[tuple[int, int], tuple[np.ndarray, np.ndarray, np.ndarray]]

predict(frequencies: ndarray, params: PrecessingParametersWithExtrinsic) → Tuple[ndarray, ndarray][source]

Predict the two polarizations of a precessing waveform.

Evaluates jax_predict() for one binary (compiled on the first call, and again for every new length of frequencies), so it requires the jax extra. For many binaries, call the function jax_predict() returns on a batch of them instead.

Parameters:

  • frequencies (np.ndarray) – Frequencies at which to evaluate the waveform, in Hz.

  • params (PrecessingParametersWithExtrinsic) – Source parameters; reference_frequency_hz must be set.

Returns:

The complex polarizations (h_plus, h_cross), in the same convention as predict().

Return type:

tuple[np.ndarray, np.ndarray]

predict_modes_dict(frequencies: ndarray, params: PrecessingParametersWithExtrinsic, source: str = 'surrogate', angles: EulerAngles | None = None) → Tuple[Dict[Tuple[int, int], ndarray], Dict[Tuple[int, int], ndarray]][source]

Inertial-frame multipoles of the precessing waveform.

Parameters:

  • frequencies (np.ndarray) – Frequencies at which to evaluate the waveform, in Hz.

  • params (PrecessingParametersWithExtrinsic) – Source parameters.

  • source (str, optional) – Where the co-precessing multipoles come from: "surrogate" (default), "post_newtonian" or "eob". See coprecessing_modes_dict().

    Default: 'surrogate'

  • angles (EulerAngles, optional) – Precomputed Euler angles, to save re-integrating the PN dynamics when comparing two sources of co-precessing multipoles for the same binary. Computed from params if None.

    Default: None

Returns:

(modes_positive_f, modes_negative_f); see twist_modes_frequency_domain().

Return type:

tuple[dict, dict]

predict_numpy(frequencies: ndarray, params: PrecessingParametersWithExtrinsic, source: str = 'surrogate', angles: EulerAngles | None = None) → Tuple[ndarray, ndarray][source]

The numpy reference implementation of predict().

Painfully slow (~20 s a waveform: the precession is integrated with scipy from the bottom of the band at every call), and only meant for testing predict() and the twist: unlike it, it can twist TEOBResumS’ own co-precessing multipoles (source="eob"), take precomputed Euler angles, and accept reference_frequency_hz = None.

Parameters:
Returns:

The complex polarizations (h_plus, h_cross), in the same convention as predict().

Return type:

tuple[np.ndarray, np.ndarray]

reference_orbital_phase(params: PrecessingParametersWithExtrinsic, source: str = 'surrogate') → float[source]

Orbital phase of the bare co-precessing multipoles at the reference frequency.

The multipoles as coprecessing_modes_dict() returns them for params.aligned(), at the time their orbital frequency is half of params.reference_frequency_hz. From the \((2, 2)\), modulo \(\pi\); the branch is the one closer to the estimate of the \((2, 1)\) (else the \((3, 3)\)), and irrelevant if the model has no odd-\(m\) multipole. See the module docstring.

Parameters:
Returns:

The orbital phase, in radians, in \((-\pi, \pi]\).

Return type:

float

teob_reference_phase(params: PrecessingParametersWithExtrinsic, start_momega: float, start_orbital_phase: float, start_phase_22: float) → float[source]

The reference_phase that reproduces a TEOBResumS precessing run.

TEOBResumS imposes the spins, and puts its orbital phase to zero, at the first sample of its integration, which two of its clocks read differently: its PN spin dynamics, along whose orbital frequency the Euler angles are read, starts at 0.95 initial_frequency (in the frequency domain); its EOB dynamics, which the waveform follows, at a slightly different frequency (~9.506 Hz for 9.5). params must give the spins at the former, reference_frequency_hz = 0.95 * initial_frequency. The orbital phase is then carried there from the latter along this model’s own phase (a ~0.4 rad advance in those 6 mHz for a binary neutron star). At the start of the integration it is not zero in the sense of the module docstring, but the offset of TEOBResumS’ \((2, 2)\) from its leading-order phase (~5e-3 rad, the multipole’s higher-order phase corrections).

Parameters:

  • params (PrecessingParametersWithExtrinsic) – The binary, with reference_frequency_hz where TEOBResumS imposes the spins.

  • start_momega (float) – \(M \Omega\) of TEOBResumS’ EOB dynamics at its first sample (MOmega[0] of the dynamics EOBRunPy returns with arg_out).

  • start_orbital_phase (float) – Its orbital phase there (Phi[0], zero).

  • start_phase_22 (float) – The phase of its time-domain \((2, 2)\) multipole there, in its own convention \(h_{22} = A e^{-i \phi}\).

Returns:

The orbital phase at params.reference_frequency_hz, in \((-\pi, \pi]\).

Return type:

float

class PrecessingParametersWithExtrinsic(mass_ratio: float, lambda_1: float, lambda_2: float, chi_1: Sequence[float], chi_2: Sequence[float], distance_mpc: float, inclination: float, total_mass: float, azimuth: float = 0.0, reference_phase: float = 0.0, merger_time: float = 0.0, reference_frequency_hz: float | None = None)[source]

Parameters of a single precessing waveform.

The aligned-spin counterpart is ParametersWithExtrinsic, of which this is the natural generalisation: the two scalar spins are promoted to full three-vectors, and an azimuthal angle is added to the inclination, since for a precessing source the two are no longer degenerate with a global phase.

The frame is the usual one at the reference frequency: \(\hat{z}\) along the orbital angular momentum \(\hat{L}_0\) there, and the spin vectors are the ones at that frequency.

Parameters:

  • mass_ratio (float) – Mass ratio \(q = m_1 / m_2 \geq 1\).

  • lambda_1 (float) – Tidal polarizabilities of the larger and smaller star.

  • lambda_2 (float) – Tidal polarizabilities of the larger and smaller star.

  • chi_1 (Sequence[float]) – Dimensionless spin vectors \(\vec{\chi}_i = \vec{S}_i / m_i^2\) of the two stars, as (x, y, z) triples.

  • chi_2 (Sequence[float]) – Dimensionless spin vectors \(\vec{\chi}_i = \vec{S}_i / m_i^2\) of the two stars, as (x, y, z) triples.

  • distance_mpc (float) – Distance to the source, in Megaparsecs.

  • inclination (float) – Polar angle \(\iota\) of the line of sight, in radians, measured from \(\hat{L}_0\).

  • azimuth (float, optional) – Azimuthal angle \(\varphi\) of the line of sight, in radians. Defaults to 0.

    Default: 0.0

  • total_mass (float) – Total mass of the binary, in solar masses.

  • reference_phase (float, optional) – Orbital phase at reference_frequency_hz, in radians (see the module docstring), which for a precessing binary is the angle between the orbital separation and the in-plane spins; at the merger, as for the aligned-spin model, if reference_frequency_hz is None. Defaults to 0.

    Default: 0.0

  • merger_time (float, optional) – Time of the merger, in seconds. Defaults to 0.

    Default: 0.0

  • reference_frequency_hz (float, optional) – The \((2, 2)\) gravitational-wave frequency, in Hz, at which chi_1, chi_2, \(\hat{L}_0 = \hat{z}\) and reference_phase are given — TEOBResumS’ initial_frequency, LALSuite’s f_ref. Precession rotates both the in-plane spins and \(\hat{L}\) between any two frequencies, so the same vectors given at two different frequencies are two different binaries. None (default) means wherever the PN spin-precession integration starts, which is well below the requested band (see euler_angles()).

    Default: None

aligned() → ParametersWithExtrinsic[source]

The aligned-spin parameters describing the co-precessing frame.

Only the spin components along \(\hat{L}_0\) survive; these are the ones the aligned-spin surrogate is trained on. The inclination is set to zero, since the co-precessing multipoles returned by coprecessing_modes_dict() carry no sky projection — the projection happens after the twist. Its coalescence phase is reference_phase only if reference_frequency_hz is None: otherwise PrecessingModel sets the orbital phase at the reference frequency itself.

Returns:

Parameters to hand to the aligned-spin Model.

Return type:

ParametersWithExtrinsic

property chi_1_vector: ndarray

Spin of the larger star as a length-3 array.

property chi_2_vector: ndarray

Spin of the smaller star as a length-3 array.

property eta: float

Symmetric mass ratio \(\eta = q / (1+q)^2\).

classmethod from_lvk(mass_1: float, mass_2: float, theta_jn: float, phi_jl: float, tilt_1: float, tilt_2: float, phi_12: float, a_1: float, a_2: float, phase: float, reference_frequency_hz: float, lambda_1: float, lambda_2: float, distance_mpc: float, merger_time: float = 0.0) → PrecessingParametersWithExtrinsic[source]

A binary given by the LVK spin angles, at reference_frequency_hz.

Masses in solar masses, the heavier first; see mlgw_bns.spin_conversion.lvk_to_precessing().

class EulerAngles(momega: ndarray, alpha: ndarray, beta: ndarray, gamma: ndarray)[source]

The precession Euler angles, as functions of the orbital frequency.

integrate_pn_spin_precession() returns \(\alpha, \beta, \gamma\) and \(M \Omega_{\rm orb}\) sampled along the PN inspiral; since \(M \Omega_{\rm orb}\) increases monotonically it can be used as the independent variable, which is what the stationary-phase lookup needs.

momega

Orbital frequency \(M \Omega_{\rm orb}\), increasing.

Type:

np.ndarray

alpha, beta, gamma

The three Euler angles at those frequencies, in radians.

Type:

np.ndarray

at_momega(momega: ndarray) → Tuple[ndarray, ndarray, ndarray][source]

Interpolate the three angles at the given orbital frequencies.

Outside the integrated range the angles are held at their boundary values: below, the PN integration simply started later than the requested frequency; above, it stopped at merger, past which the twist is in any case only a crude extrapolation.

Parameters:

momega (np.ndarray) – Orbital frequencies \(M \Omega_{\rm orb}\) at which the angles are wanted.

Returns:

(alpha, beta, gamma), each of the shape of momega.

Return type:

tuple[np.ndarray, np.ndarray, np.ndarray]

Notes

\(\alpha\) and \(\gamma\) are unwrapped before interpolating: they wind through many turns over an inspiral, and interpolating the wrapped values would smear each \(2\pi\) jump across a whole interval.

A monotone cubic Hermite interpolant (PCHIP) is used in place of the plain linear one used before, closer in spirit to TEOBResumS’ own gsl_interp_cspline lookup of these angles (twist_hlm_FD/twist_hlm_TD in TEOBResumSWaveform.c): the PN integration returns only a sparse set of accepted ODE steps (see integrate_pn_spin_precession()), and linear interpolation between them is a coarser approximation of the (smooth) angle trajectory. PCHIP is shape-preserving (no overshoot between knots, whatever their spacing, which spans orders of magnitude here), and it still reproduces the ODE solution exactly at every returned node.

euler_angles(params: PrecessingParametersWithExtrinsic, initial_frequency_hz: float, largest_mode_m: int = 4, omega_dot: Callable[[float], float] | None = None) → EulerAngles[source]

Integrate the PN spin-precession dynamics for a given source.

Thin wrapper around integrate_pn_spin_precession() that takes care of the units and of choosing a low enough starting frequency.

Parameters:

  • params (PrecessingParametersWithExtrinsic) – Source parameters; only the intrinsic ones matter here.

  • initial_frequency_hz (float) – Lowest frequency, in Hz, at which the waveform will be evaluated.

  • largest_mode_m (int, optional) – Largest \(|m|\) among the multipoles that will be twisted. The multipole \(m\) reaches the frequency initial_frequency_hz when the orbit is only at \(2/m\) of the corresponding \((2,2)\) frequency, so the integration has to start that much earlier. Defaults to 4.

    Default: 4

  • omega_dot (callable, optional) – M \Omega_{\rm orb} \mapsto \mathrm{d}(M\Omega_{\rm orb})/ \mathrm{d}t in geometric units. When given, the precession is integrated against orbital frequency, marched along this rate instead of the PN flux – reproducing the SPIN_FLX_EOB hand-off TEOBResumS does once the spin dynamics reaches the EOB band. Build it from an accurate \((2, 2)\) phase with eob_orbital_frequency_rate(). When None (default), the integration marches in time with the 3.5PN flux, exactly as TEOBResumS’ SPIN_FLX_PN branch.

    Default: None

Returns:

The angles, tabulated against \(M \Omega_{\rm orb}\).

Return type:

EulerAngles

backward_gamma(gamma, momega, in_plane_l, reference_momega, xp=<module 'numpy' from '/home/docs/checkouts/readthedocs.org/user_builds/mlgw-bns/envs/latest/lib/python3.13/site-packages/numpy/__init__.py'>)[source]

\(\gamma\), continued through the reference point.

At the reference point \(\hat{L} = \hat{z}\), so the in-plane part of \(\hat{L}\) passes through zero there and changes direction: \(\alpha = \arctan(L_y / L_x)\) jumps by \(\pi\), while \(\gamma\), integrated, does not. The rotation \(R(\alpha + \pi, \beta, \gamma) = R(\alpha, -\beta, \gamma + \pi)\) then jumps by \(\pi\) about \(\hat{L}\), flipping the sign of every odd-\(m\) co-precessing multipole below the reference frequency. TEOBResumS’ backward integration does the same, invisibly: its multipoles start at the reference point. Adding \(\pi\) to \(\gamma\) where \(\alpha\) is flipped (below the reference, wherever \(\hat{L} \neq \hat{z}\)) makes the rotation continuous; with no in-plane spin nothing changes.

LEADING_ORDER_MODE_PHASES: Dict[Tuple[int, int], float]

Phase of a co-precessing multipole at zero orbital phase, to leading (Newtonian) order: \(\delta_{\ell m}\) in \(X_{\ell m}(m F) = m \phi_{\rm orb} + \delta_{\ell m} + \pi / 4\) (see the module docstring). TEOBResumS’ multipoles at the first sample of its integration have these plus ~5e-3 m rad at 10 Hz (a binary neutron star), the same for every binary.

twist_modes_frequency_domain(coprecessing_modes: ~typing.Dict[~typing.Tuple[int, int], ~numpy.ndarray], frequencies: ~numpy.ndarray, angles: ~mlgw_bns.precessing_model.EulerAngles, mass_sum_seconds: float, xp=<module 'numpy' from '/home/docs/checkouts/readthedocs.org/user_builds/mlgw-bns/envs/latest/lib/python3.13/site-packages/numpy/__init__.py'>) → Tuple[Dict[Tuple[int, int], ndarray], Dict[Tuple[int, int], ndarray]][source]

Twist frequency-domain co-precessing multipoles into the inertial frame.

Applies the Wigner-\(D\) rotation of twist_modes() frequency by frequency, with the Euler angles evaluated at each co-precessing multipole’s own stationary-phase point (see the module docstring for why).

Parameters:

  • coprecessing_modes (dict[tuple[int, int], np.ndarray]) – Co-precessing multipoles \(\tilde{h}^{\rm co-prec}_{\ell n}(f)\) for \(n > 0\), as returned by coprecessing_modes_dict().

  • frequencies (np.ndarray) – The (positive) frequency grid, in Hz.

  • angles (EulerAngles) – Precession angles, tabulated against orbital frequency: anything with an at_momega method, such as TabulatedAngles.

  • mass_sum_seconds (float) – Total mass of the binary, in seconds, used to convert the frequencies to geometric units.

  • xp (module, optional) – Array namespace, numpy (default) or jax.numpy; with the latter, frequencies and mass_sum_seconds may carry a batch axis, (N, k) and (N, 1).

    Default: <module 'numpy' from '/home/docs/checkouts/readthedocs.org/user_builds/mlgw-bns/envs/latest/lib/python3.13/site-packages/numpy/__init__.py'>

Returns:

(modes_positive_f, modes_negative_f): the inertial-frame multipoles \(\tilde{h}^{\rm inertial}_{\ell m}(f)\) and \(\tilde{h}^{\rm inertial}_{\ell m}(-f)\), both keyed by \((\ell, m)\) and both spanning all \(-\ell \leq m \leq \ell\) for every \(\ell\) present in the input.

Return type:

tuple[dict, dict]

polarizations_from_inertial_modes(modes_positive_f: ~typing.Dict[~typing.Tuple[int, int], ~numpy.ndarray], modes_negative_f: ~typing.Dict[~typing.Tuple[int, int], ~numpy.ndarray], inclination: float, azimuth: float, xp=<module 'numpy' from '/home/docs/checkouts/readthedocs.org/user_builds/mlgw-bns/envs/latest/lib/python3.13/site-packages/numpy/__init__.py'>) → Tuple[ndarray, ndarray][source]

Project the inertial-frame multipoles onto the observer’s sky.

On the positive-frequency grid the two independent combinations of the polarizations are

\[\begin{split}\tilde{h}_+ + i \tilde{h}_\times &= \sum_{\ell m} \tilde{h}_{\ell m}(f)\; {}_{-2}Y_{\ell m}(\iota, \varphi) \\ \tilde{h}_+ - i \tilde{h}_\times &= \sum_{\ell m} \tilde{h}_{\ell m}^*(-f)\; {}_{-2}Y_{\ell m}^*(\iota, \varphi)\end{split}\]

the second following from the reality of the time-domain strain. For an aligned-spin source the two are related by the equatorial symmetry of the multipoles and only the first is needed; under precession that symmetry is broken and both are required.

Parameters:

  • modes_positive_f (dict[tuple[int, int], np.ndarray]) – Inertial-frame multipoles at \(+f\) and \(-f\), as returned by twist_modes_frequency_domain().

  • modes_negative_f (dict[tuple[int, int], np.ndarray]) – Inertial-frame multipoles at \(+f\) and \(-f\), as returned by twist_modes_frequency_domain().

  • inclination (float) – Polar angle \(\iota\) of the line of sight, in radians.

  • azimuth (float) – Azimuthal angle \(\varphi\) of the line of sight, in radians.

  • xp (module, optional) – Array namespace, numpy (default) or jax.numpy; the angles may then be arrays broadcasting against the multipoles.

    Default: <module 'numpy' from '/home/docs/checkouts/readthedocs.org/user_builds/mlgw-bns/envs/latest/lib/python3.13/site-packages/numpy/__init__.py'>

Returns:

The complex polarizations (h_plus, h_cross), in the same convention as predict().

Return type:

tuple[np.ndarray, np.ndarray]

twist_coefficients(keys: ~typing.Sequence[~typing.Tuple[int, int]], frequencies: ~numpy.ndarray, angles: ~mlgw_bns.precessing_model.EulerAngles, mass_sum_seconds: float, inclination: float, azimuth: float, xp=<module 'numpy' from '/home/docs/checkouts/readthedocs.org/user_builds/mlgw-bns/envs/latest/lib/python3.13/site-packages/numpy/__init__.py'>) → Dict[Tuple[int, int], Tuple[ndarray, ndarray]][source]

The twist and sky projection of each co-precessing multipole.

The polarizations are linear in the co-precessing multipoles:

\[\tilde{h}_{+, \times}(f) = \sum_{\ell m} c^{+, \times}_{\ell m}(f) \, \tilde{h}^{\rm co}_{\ell m}(f) \,.\]

twist_modes_frequency_domain() builds the inertial multipoles at \(-f\) from the complex conjugates of the co-precessing ones, but polarizations_from_inertial_modes() conjugates those again, so the coefficients are those of a unit multipole twisted and projected. They vary only on the precession time scale — what mode-by-mode relative binning interpolates (Leslie, Dai and Pratten, arXiv:2109.09872, eq. 3).

Parameters:

  • keys (sequence of (l, m)) – Co-precessing multipoles, \(m > 0\).

  • frequencies (np.ndarray) – The (positive) frequency grid, in Hz.

  • angles (EulerAngles) – Precession angles; anything with an at_momega method.

  • mass_sum_seconds (float) – Total mass of the binary, in seconds.

  • inclination (float) – Line of sight, in radians; see PrecessingParametersWithExtrinsic.

  • azimuth (float) – Line of sight, in radians; see PrecessingParametersWithExtrinsic.

  • xp (module, optional) – Array namespace, numpy (default) or jax.numpy.

    Default: <module 'numpy' from '/home/docs/checkouts/readthedocs.org/user_builds/mlgw-bns/envs/latest/lib/python3.13/site-packages/numpy/__init__.py'>

Returns:

Mapping (l, m) -> (c_plus, c_cross), on frequencies.

Return type:

dict[tuple[int, int], tuple[np.ndarray, np.ndarray]]

newtonian_time_to_merger(eta: float, momega: float) → float[source]

Leading-order time left to merger, in units of the total mass.

\[\frac{t_{\rm c}}{M} = \frac{5}{256\, \eta\, v^8} \,, \qquad v = (M \Omega_{\rm orb})^{1/3}\]

Used only to bound the PN spin-precession integration: it runs until the merger fit stops it, and a cutoff comfortably past the Newtonian estimate keeps the integrator from being the thing that stops it. For a binary neutron star entering the band at a few tens of Hz this is of order \(10^8 M\), several orders of magnitude beyond the default cutoff of integrate_pn_spin_precession().

Parameters:

  • eta (float) – Symmetric mass ratio.

  • momega (float) – Orbital frequency \(M \Omega_{\rm orb}\).

Returns:

Time to merger, in units of the total mass.

Return type:

float

check_aligned_spin_limit(model: Model, frequencies: ndarray, params: PrecessingParametersWithExtrinsic, rtol: float = 1e-08) → None[source]

Assert that the twist is the identity for aligned spins.

With no in-plane spin the orbital angular momentum never tilts, so \(\beta \equiv 0\) and the Wigner-\(D\) matrices collapse to the identity: the precessing pipeline must then return exactly what predict() returns. This exercises every convention in the chain at once — the sign of the rotation, the split between the \(\pm f\) halves of the twist, and the normalisation of the multipoles — so it is the check worth running after touching any of them.

Parameters:

  • model (Model) – The aligned-spin surrogate.

  • frequencies (np.ndarray) – Frequencies at which to compare, in Hz.

  • params (PrecessingParametersWithExtrinsic) – Source parameters. Its in-plane spin components are ignored: the comparison is made for the aligned-spin source with the same \(\chi_z\).

  • rtol (float, optional) – Relative tolerance of the comparison. Defaults to 1e-8.

    Default: 1e-08

Raises:

AssertionError – If the two waveforms differ by more than rtol.

Batched, in JAX

Batched precessing waveforms, in JAX.

PrecessingModel.predict_numpy for N binaries in one call, as a pure JAX function: it can be wrapped in jax.jit(), vmap-ed and differentiated. PrecessingModel.jax_predict builds one. Every stage is the numpy one, run on jax.numpy:

  • the co-precessing multipoles come from the batched surrogate, Model.jax_modes_amp_phase;

  • their orbital phase is set at the reference frequency with orbital_phase_from_transforms(), from the stationary-phase transform \(X = \Psi - f \Psi'\) of the surrogate’s own phase and its analytic derivative (return_tf) instead of the cubic window fit of stationary_phase_transform(), which it matches to ~1e-6 rad;

  • the twist and the projection are twist_coefficients(), which applies twist_modes_frequency_domain() and polarizations_from_inertial_modes() to a unit multipole: the polarizations are linear in the co-precessing multipoles, and precessing_mode_components() returns them and their coefficients separately, as mode-by-mode relative binning needs.

The one stage that is not shared is the integration of the PN spin-precession equations (the same right-hand side, mlgw_bns.twist_waveform._pn_precession_derivatives()). The numpy path integrates them in time with SciPy’s adaptive DOP853, 3–12 s per binary; here they are integrated with DOP853’s own coefficients on a fixed grid (integrate_angles()), against the orbital frequency, which is the same trajectory: \(\mathrm{d}y/\mathrm{d}\Omega = (\mathrm{d}y / \mathrm{d}t) / \dot\Omega\), with the full 3.5PN + tidal \(\dot\Omega\) (positive throughout the BNS parameter space). The grid is uniform in

\[x = \ln \Omega - \kappa \, \Omega_{\rm lo} / \Omega ,\]

which is uniform in \(1/\Omega\) at low frequency, where \(\hat{L}\) winds around \(\vec{J}\) hundreds of times (the precession phase grows as \(\Omega^{-1}\)), and in \(\ln \Omega\) near the merger; KAPPA sets the balance. The state between the nodes is a cubic Hermite interpolant with the exact derivatives of the right-hand side, and the Euler angles are computed from the interpolated \(\hat{L}\) (TabulatedAngles).

precessing_waveform(model: Model, modes: Sequence | None = None, n_steps: int = 2048) → Callable[source]

A JAX function reproducing PrecessingModel.predict_numpy, batched.

Returns predict(intrinsic, frequencies, total_mass, distance_mpc, inclination, azimuth, reference_phase, reference_frequency_hz, merger_time=0.0) -> (h_plus, h_cross), where intrinsic has shape (N, 9), rows \([q, \Lambda_1, \Lambda_2, \vec{\chi}_1, \vec{\chi}_2]\), frequencies (k,) or (N, k) (positive and increasing, in Hz) and the rest are scalars or of shape (N,); the polarizations have shape (N, k). The parameters are those of PrecessingParametersWithExtrinsic; the reference frequency is required. The Euler angles are integrated from the lowest frequency of each row, as the numpy path does. The sum over the modes of precessing_mode_components().

Parameters:

  • model (Model) – The aligned-spin surrogate.

  • modes (sequence of (l, m), optional) – Co-precessing multipoles to twist; defaults to all of model.modes.

    Default: None

  • n_steps (int, optional) – Steps of each leg of the precession integration (integrate_angles()).

    Default: 2048

precessing_mode_components(model: Model, modes: Sequence | None = None, n_steps: int = 2048) → Callable[source]

A JAX function giving each co-precessing multipole and its twist, batched.

The precessing polarizations are linear in the co-precessing multipoles,

\[\tilde{h}_{+, \times}(f) = \sum_{\ell m} c^{+, \times}_{\ell m}(f) \, \tilde{h}^{\rm co}_{\ell m}(f) \,,\]

where the coefficients \(c_{\ell m}\) hold the twist (the Euler angles at the multipole’s own stationary-phase point) and the projection on the sky (twist_coefficients()). The co-precessing multipoles carry the fast orbital phase, the coefficients only the slow precession: this is the split that mode-by-mode relative binning needs (Leslie, Dai and Pratten, arXiv:2109.09872, eqs. 1–4).

Returns predict(..., angles=None) -> (coprecessing, c_plus, c_cross), with the arguments of precessing_waveform() and three complex arrays of shape (N, n_modes, k), the modes in the order of modes. angles, from precession_angles(), are the precession angles to use; by default they are integrated from the lowest frequency of each row. coprecessing already has its orbital phase set at the reference frequency and the merger-time shift; precessing_waveform() is (sum(c_plus * coprecessing), sum(c_cross * coprecessing)) over the modes.

Parameters:

  • model (Model) – The aligned-spin surrogate.

  • modes (sequence of (l, m), optional) – Co-precessing multipoles; defaults to all of model.modes.

    Default: None

  • n_steps (int, optional) – Steps of each leg of the precession integration (integrate_angles()).

    Default: 2048

precession_angles(model: Model, n_steps: int = 2048) → Callable[source]

A JAX function integrating the precession of a batch of binaries, for precessing_mode_components() (its angles argument).

Returns angles(intrinsic, total_mass, reference_frequency_hz, start_frequency_hz) -> tuple, with intrinsic of shape (N, 9) as for precessing_waveform() and the rest scalars or of shape (N,): the fields of each binary’s TabulatedAngles, stacked along a first axis of length N, integrated from the lowest frequency any multipole of model needs at start_frequency_hz. The angles do not depend on the frequencies the waveform is then evaluated at, so a long frequency array can be evaluated in pieces with one integration.

Parameters:

  • model (Model) – The aligned-spin surrogate.

  • n_steps (int, optional) – Steps of each leg of the integration (integrate_angles()).

    Default: 2048

batch_arguments(params: Sequence['PrecessingParametersWithExtrinsic'], frequencies) → tuple[source]

The arguments of a precessing_waveform() function for the binaries params, as numpy arrays.

Parameters:
integrate_angles(mass_ratio, lambda_1, lambda_2, chi_1, chi_2, reference_frequency_22, initial_frequency_22, n_steps: int = 2048, final_frequency_22=None, full_state: bool = False)[source]

The PN spin-precession dynamics of one binary, on a fixed grid.

What euler_angles() integrates (the spins given at reference_frequency_22, integrated down to initial_frequency_22 and up to TEOBResumS’ stopping frequency \(1.1\,\Omega_{\rm mrg}\)), in JAX; see the module docstring.

Parameters:

  • mass_ratio (scalar) – \(q \geq 1\) and the tidal polarizabilities.

  • lambda_1 (scalar) – \(q \geq 1\) and the tidal polarizabilities.

  • lambda_2 (scalar) – \(q \geq 1\) and the tidal polarizabilities.

  • chi_1 (array of shape (3,)) – Dimensionless spin vectors at the reference frequency.

  • chi_2 (array of shape (3,)) – Dimensionless spin vectors at the reference frequency.

  • reference_frequency_22 (scalar) – \((2, 2)\) GW frequencies in geometric units (\(M f\)): where the spins are given, and the lowest needed.

  • initial_frequency_22 (scalar) – \((2, 2)\) GW frequencies in geometric units (\(M f\)): where the spins are given, and the lowest needed.

  • n_steps (int, optional) – Steps of each of the two legs.

    Default: 2048

  • final_frequency_22 (scalar, optional) – \((2, 2)\) frequency, geometric units, to integrate up to instead of TEOBResumS’ stopping frequency.

    Default: None

  • full_state (bool, optional) – Also return the whole state \([\vec{S}_A, \vec{S}_B, \hat{L}, \alpha - \gamma, \Omega]\) and its derivative at the nodes.

    Default: False

Returns:

Or (angles, state, derivative) with full_state, the last two of shape (2 n_steps + 1, 11).

Return type:

TabulatedAngles

class TabulatedAngles(x: Any, state: Any, derivative: Any, kappa_omega: Any, omega_lo: Any, omega_hi: Any, alpha_reference: Any)[source]

The Euler angles of one binary, from integrate_angles().

The JAX counterpart of EulerAngles, with the same at_momega(). The state [\hat{L}, \alpha - \gamma] is tabulated at 2 n + 1 nodes uniform in \(x\) on either side of the reference, x[n], with its derivatives.

at_momega(momega)[source]

\((\alpha, \beta, \gamma)\) at the orbital frequencies momega, held at the end values outside the integrated range, as EulerAngles.at_momega.

From the LVK spin angles

Precessing binaries from the LVK spin angles.

Parameter estimation codes (LALInference, bilby, …) describe the orientation of a precessing binary with the angles of LALSuite’s SimInspiralTransformPrecessingNewInitialConditions, all at the reference frequency \(f_{\rm ref}\):

  • theta_jn, the angle between the total angular momentum \(\vec{J}\) and the line of sight \(\hat{N}\);

  • phi_jl, the azimuth of \(\hat{L}\) on its cone about \(\vec{J}\), measured from the plane of \(\vec{J}\) and \(\hat{N}\);

  • tilt_1, tilt_2, the angles between the spins and \(\hat{L}\), phi_12 the azimuth of the second spin about \(\hat{L}\) measured from the first, and a_1, a_2 their magnitudes;

  • phase, the orbital phase.

lvk_to_precessing() is a port of that function, vectorised and generic over numpy and jax.numpy: it gives the spins in the frame of PrecessingParametersWithExtrinsic, \(\hat{z} = \hat{L}_0\), and the line of sight in it. As in LALSuite, \(\vec{J} = \vec{L} + \vec{S}_1 + \vec{S}_2\) with the 1PN orbital angular momentum.

The orbital phase enters only as reference_phase: the spins are those LALSuite gives for phiRef = 0, with the line of sight at azimuth \(\pi / 2\), and the orbital separation is rotated by phase with respect to them, which is the same binary as LALSuite’s (the spins and the line of sight rotated by -phase about \(\hat{L}_0\), the separation along \(\hat{x}\)). So phase only rotates the co-precessing multipole \(m\) by \(e^{i m \, {\rm phase}}\), and the precession angles do not depend on it.

lvk_to_precessing(theta_jn, phi_jl, tilt_1, tilt_2, phi_12, a_1, a_2, mass_1, mass_2, reference_frequency_hz, phase, xp=<module 'numpy' from '/home/docs/checkouts/readthedocs.org/user_builds/mlgw-bns/envs/latest/lib/python3.13/site-packages/numpy/__init__.py'>) → PrecessingOrientation[source]

The orientation of a precessing binary given its LVK spin angles.

Masses in solar masses, the heavier first (a_1, tilt_1 are its spin’s); all angles in radians, at reference_frequency_hz. See the module docstring.

Returns:

inclination, azimuth, reference_phase, chi_1 and chi_2 of PrecessingParametersWithExtrinsic.

Return type:

PrecessingOrientation

class PrecessingOrientation(inclination: object, azimuth: object, reference_phase: object, chi_1: object, chi_2: object)[source]

The output of lvk_to_precessing(); arrays of the input shape, the spins with a trailing axis of length 3.

lal_precessing_spins(theta_jn, phi_jl, tilt_1, tilt_2, phi_12, a_1, a_2, mass_1, mass_2, reference_frequency_hz, phase, xp=<module 'numpy' from '/home/docs/checkouts/readthedocs.org/user_builds/mlgw-bns/envs/latest/lib/python3.13/site-packages/numpy/__init__.py'>)[source]

SimInspiralTransformPrecessingNewInitialConditions, vectorised.

Masses in solar masses, the heavier first. Returns (inclination, s1x, s1y, s1z, s2x, s2y, s2z): the angle between \(\hat{L}_0\) and \(\hat{N}\), and the dimensionless spins in LALSuite’s frame, \(\hat{z} = \hat{L}_0\), \(\hat{x}\) along the orbital separation, and the line of sight at azimuth \(\pi / 2 - {\rm phase}\).

orbital_angular_momentum(total_mass, eta, reference_frequency_hz, xp=<module 'numpy' from '/home/docs/checkouts/readthedocs.org/user_builds/mlgw-bns/envs/latest/lib/python3.13/site-packages/numpy/__init__.py'>)[source]

\(|\vec{L}|\) at 1PN, in units of the solar mass squared, as in SimInspiralTransformPrecessingNewInitialConditions: \(\eta M^2 / v \, (1 + v^2 (3/2 + \eta / 6))\), \(v = (\pi M f_{\rm ref})^{1/3}\).

The precession dynamics

Standalone reproduction of the TEOBResumS waveform “twisting” procedure.

TEOBResumS generates the precessing waveform by first evolving the dynamics in a co-precessing frame (where only the m=+-2 multipoles are large) and then rotating (“twisting”) the co-precessing multipoles h_{l,n}(t) = A_{l,n}(t) * exp(-i*phi_{l,n}(t)) (n = 1 .. l) into the fixed inertial frame using three time-dependent Euler angles alpha(t), beta(t), gamma(t):

h^inertial_{l,m}(t) = sum_n D^l_{m,n}(alpha,beta,gamma) h^coprec_{l,n}(t)

with the D-matrix built out of Wigner-d functions d^l_{m,n}(beta) and the usual m=<0 symmetry of the co-precessing multipoles, H_{l,-n} = (-1)^l H_{l,n}^*.

This module reimplements, line for line, the algebra of:

  • rmap_twist() in C/src/TEOBResumSUtils.c

  • twist_hlm_TD() in C/src/TEOBResumSWaveform.c

  • compute_hpc() in C/src/TEOBResumSWaveform.c (bonus: h+, hx)

as well as the PN spin-precession dynamics used to produce alpha(t), beta(t), gamma(t) in the EOBPars->project_spins == 0 (default), the EOBPars->spin_flx == SPIN_FLX_PN (default) branch:

  • eob_spin_dyn_rhs_PN() in C/src/TEOBResumSDynamics.c (SPIN_FLX_PN case)

  • eob_spin_dyn() “ “ (initial data)

  • eob_spin_dyn_alpha/beta() “ “

  • alpha_initial_condition() “ “

  • eob_mrg_momg() in C/src/TEOBResumSFits.c (stop frequency)

(The Wigner d-function and the spin-weighted spherical harmonics of C/src/SpecialFuns.c live in mlgw_bns.special_func, which is shared with the rest of the package.)

Convention (matches the C code throughout):

h_{l,m} = A_{l,m} * exp(-i * phi_{l,m}), phi_{l,m} increasing in time.

For the surrogate-facing, frequency-domain use of all this — taking the aligned-spin mlgw_bns multipoles as the co-precessing ones and twisting them into a full precessing waveform — see mlgw_bns.precessing_model.

integrate_pn_spin_precession(nu, chi1vec, chi2vec, f0, t_max=1000000.0, max_step=None, refinement=32, rtol=1e-09, atol=1e-11, independent_variable='time', omega_dot=None, f_start=None, lambdas=(0.0, 0.0))[source]

Integrate the closed PN spin-precession ODE system that TEOBResumS solves in eob_spin_dyn()/eob_spin_dyn_integrate() for EOBPars->use_spins == MODE_SPINS_GENERIC, EOBPars->project_spins == 0 (the default) and EOBPars->spin_flx == SPIN_FLX_PN (the default).

independent_variable selects what the system is marched against:

  • "time" (default) reproduces TEOBResumS’ SPIN_FLX_PN branch – state [SA, SB, Lh, gamma, Momega], Momega advanced by the 3.5PN energy-balance flux. The Euler angles are accurate as functions of time, but the orbital frequency they end up labelled with drifts from the true one through the late inspiral, because the PN flux runs fast there. In the C code this is why, once the spin dynamics reaches the EOB band, TEOBResumS stops using the PN flux and drives Momega from the EOB dynamics instead (the SPIN_FLX_EOB hand-off in eob_spin_dyn_integrate). 't' in the returned dict is the integration time.

  • "orbital_frequency" marches against Momega itself – dy/dMomega = (dy/dt) / omega_dot. This is what reproduces the SPIN_FLX_EOB hand-off: pass omega_dot as a callable Momega -> dMomega/dt taken from an accurate (EOB / surrogate) (2,2) phase and the angles come out both integrated along and labelled by the true orbital frequency. With omega_dot=None the 2PN-truncated PN flux domg_leading is used instead, which on its own barely differs from the "time" result. 't' is reconstructed by quadrature of dt/dMomega alongside the angles.

omega_dotcallable or None

Only used when independent_variable == "orbital_frequency". omega_dot(Momega) -> dMomega/dt in geometric units; supplies the orbital-frequency evolution the precession is marched along, in place of the PN flux. Must be positive and defined over [pi*f0, omg_stop].

Note on alpha, beta: in the C code these are formally ODE state variables, but their r.h.s. is identically zero (see the “alpha and beta are not actually evolved” comment in eob_spin_dyn_rhs_PN); after every accepted forward step they are instead overwritten algebraically from the new Lhat: alpha = atan2(Lhat_y, Lhat_x), beta = acos(Lhat_z). That is mathematically identical to treating them as derived quantities of Lhat(t) rather than as independent state, which is what is done here; only SA, SB, Lhat, gamma and Momega are integrated.

Parameters:

  • nu (symmetric mass ratio (sets MA = m1/M >= MB = m2/M via) – MA = nu_to_X1(nu))

  • chi1vec (length-3 sequences) – Dimensionless spin vectors of body 1 (mass fraction MA) and body 2 (mass fraction MB).

  • chi2vec (length-3 sequences) – Dimensionless spin vectors of body 1 (mass fraction MA) and body 2 (mass fraction MB).

  • f0 (reference GW (2,2) frequency, M*f0 in geometric units (so the) – reference M*Omega_orb = pi*f0, matching EOBPars->initial_frequency): where chi1vec, chi2vec are given and Lhat = z.

  • lambdas ((float, float), optional) – Quadrupolar tidal polarizabilities of bodies 1 and 2, entering the "time" branch’s dOmega/dt through the LO 5PN tidal term (see tidal_flux_coefficient()), as in TEOBResumS for a BNS. Leaving it out on a BNS makes Momega(t) run slow through the late inspiral: the angles, labelled by Momega, drift from TEOBResumS’ by ~0.1 rad by 400 Hz.

    Default: (0.0, 0.0)

  • f_start (float or None, optional) – GW (2,2) frequency, geometric units, to start the returned angles from. When below f0 the system is also integrated backwards from the reference state down to it, as TEOBResumS does (eob_spin_dyn_integrate_backwards) when its EOB dynamics starts below initial_frequency; the spins keep their meaning at f0. None (default) or >= f0: start at f0.

    Default: None

  • t_max (safety cutoff on the integration time (geometric units, M=1)., optional) –

    Default: 1000000.0

  • refinement (int, optional) – Number of sub-intervals each accepted step of the integrator is split into, using its dense output, before the angles are returned. The integrator is accurate but takes long steps — a whole inspiral can be covered in a couple of hundred of them — and the returned arrays are meant to be interpolated, so the grid has to resolve the precession, not merely carry it. Without this the angles are perfectly accurate at the returned points and visibly wrong between them. Costs only polynomial evaluations, no extra right-hand sides.

    Default: 32

  • rtol (float, optional) – Tolerances passed to the integrator.

    Default: 1e-09

  • atol (float, optional) – Tolerances passed to the integrator.

    Default: 1e-11

  • independent_variable ({"time", "orbital_frequency"}, optional) – What the system is marched against; see the discussion above.

    Default: 'time'

  • omega_dot (callable or None, optional) – Orbital-frequency rate for the "orbital_frequency" branch; see above.

    Default: None

Returns:

  • dict with keys ‘t’, ‘alpha’, ‘beta’, ‘gamma’, ‘Momega’, ‘SA’, ‘SB’, ‘Lh’

  • (ndarrays; SA/SB/Lh have shape (N,3)). In the "orbital_frequency"

  • branch 'Momega' is the (exact) integration grid and 't' its

  • quadrature; in the "time" branch it is the reverse.

twist_modes(hlm_coprec, alpha, beta, gamma, lm_inertial, include_positive_n=True, include_negative_n=True, xp=<module 'numpy' from '/home/docs/checkouts/readthedocs.org/user_builds/mlgw-bns/envs/latest/lib/python3.13/site-packages/numpy/__init__.py'>)[source]

Twist co-precessing multipoles into the inertial frame.

Parameters:

  • hlm_coprec (dict[(l, n)] -> complex ndarray) – Co-precessing multipoles for n = 1 .. l (m>0 only; the m<0 co-precessing multipoles are reconstructed from the standard symmetry H_{l,-n} = (-1)^l H_{l,n}^*, exactly as the C code does). h = A*exp(-i*phi) convention. All arrays share the same grid.

  • alpha (ndarray) – Euler angles (radians) on that same grid.

  • beta (ndarray) – Euler angles (radians) on that same grid.

  • gamma (ndarray) – Euler angles (radians) on that same grid.

  • lm_inertial (iterable of (l, m), m > 0) – Inertial-frame multipoles to produce.

  • include_positive_n (bool, optional) – Whether to include, in the sum over the co-precessing index n, the n>0 multipoles and their n<0 partners respectively. Both default to True, which is what the time-domain twist of the C code does. Splitting the two halves apart is what lets the same routine be used in the frequency domain, where the n>0 and n<0 co-precessing multipoles have their support at opposite signs of the frequency and so must be twisted separately — see twist_modes_frequency_domain().

    Default: True

  • include_negative_n (bool, optional) – Whether to include, in the sum over the co-precessing index n, the n>0 multipoles and their n<0 partners respectively. Both default to True, which is what the time-domain twist of the C code does. Splitting the two halves apart is what lets the same routine be used in the frequency domain, where the n>0 and n<0 co-precessing multipoles have their support at opposite signs of the frequency and so must be twisted separately — see twist_modes_frequency_domain().

    Default: True

  • xp (module, optional) – Array namespace, numpy (default) or jax.numpy; the arrays may then have any (common) shape.

    Default: <module 'numpy' from '/home/docs/checkouts/readthedocs.org/user_builds/mlgw-bns/envs/latest/lib/python3.13/site-packages/numpy/__init__.py'>

Returns:

hTlm, hTlm_neg, hTl0 – Twisted inertial multipoles for m>0, m<0 (keyed by the negative m) and m=0 respectively, h = A*exp(-i*phi) convention.

Return type:

dict[(l, m)] -> complex ndarray

compute_hpc(hTlm, hTlm_neg, hTl0, phi, iota)[source]

Combine twisted inertial multipoles into h+, hx for a given azimuth phi and inclination iota, mirroring compute_hpc() in the “generic spins” (precessing) branch: no +m/-m symmetry is assumed, m<0 and m=0 multipoles are used directly.

phi is the azimuth of the line of sight, entering the harmonics as exp(i m phi). TEOBResumS parametrises the same freedom through its coalescence_angle parameter, with phi = pi/2 - coalescence_angle, so its default of 0 corresponds to phi = pi/2. With that identification this function reproduces the C code’s h+ and hx to machine precision; see visualization/validate_twist_against_teob.py.

Returns h = h+ - i*hx as a complex ndarray.

Special functions

Shared special functions and small numerical helpers.

The Wigner \(d\) and spin-weighted spherical harmonic implementations here are the single source of truth for the whole package: model uses them to project the co-precessing modes onto the observer’s sky, higher_order_modes uses them for the \(s=-2\) specialisation, and twist_waveform uses them to build the Wigner-\(D\) matrices of the precessing twist. They mirror wigner_d_function() and spinsphericalharm() in TEOBResumS’ C/src/SpecialFuns.c.

Both accept either a scalar angle or a numpy array of angles, and run on jax.numpy arrays as well when passed xp=jax.numpy.

wigner_d_function(l: int, m: int, s: int, i: float | ndarray, xp=<module 'numpy' from '/home/docs/checkouts/readthedocs.org/user_builds/mlgw-bns/envs/latest/lib/python3.13/site-packages/numpy/__init__.py'>) → float | ndarray[source]

Wigner \(d\)-function \(d^{\ell}_{m,s}(\iota)\).

As given in Eq. (II.8) of arXiv:0709.0093.

Parameters:

  • l (int) – Upper index \(\ell\).

  • m (int) – First lower index.

  • s (int) – Second lower index.

  • i (float or np.ndarray) – Argument of the Wigner \(d\)-function, in radians. May be an array, in which case an array of the same shape is returned.

  • xp (module, optional) – Array namespace, numpy (default) or jax.numpy.

    Default: <module 'numpy' from '/home/docs/checkouts/readthedocs.org/user_builds/mlgw-bns/envs/latest/lib/python3.13/site-packages/numpy/__init__.py'>

Returns:

The value of \(d^{\ell}_{m,s}(\iota)\).

Return type:

float or np.ndarray

spinsphericalharm(s: int, l: int, m: int, phi: float | ndarray, i: float | ndarray, xp=<module 'numpy' from '/home/docs/checkouts/readthedocs.org/user_builds/mlgw-bns/envs/latest/lib/python3.13/site-packages/numpy/__init__.py'>) → Tuple[float | ndarray, float | ndarray][source]

Spin-weighted spherical harmonic \({}_{s}Y_{\ell m}(\varphi, \iota)\).

As given in Eq. (II.7) of arXiv:0709.0093.

Parameters:

  • s (int) – Spin weight.

  • l (int) – Multipolar index \(\ell\).

  • m (int) – Multipolar index \(m\).

  • phi (float or np.ndarray) – Azimuthal angle, in radians.

  • i (float or np.ndarray) – Polar angle, in radians.

  • xp (module, optional) – Array namespace, numpy (default) or jax.numpy.

    Default: <module 'numpy' from '/home/docs/checkouts/readthedocs.org/user_builds/mlgw-bns/envs/latest/lib/python3.13/site-packages/numpy/__init__.py'>

Returns:

(rY, iY): the real and imaginary parts of \({}_{s}Y_{\ell m}\).

Return type:

tuple

Raises:

ValueError – If (l, m) are not a valid pair of indices.