First frequency-domain phenomenological model of the multipole asymmetry in gravitational-wave signals from binary-black-hole coalescence

Gravitational-wave signals from binaries that contain spinning black holes in general include an asymmetry between the $+m$ and $-m$ multipoles that is not included in most signal models used in LIGO-Virgo-KAGRA (LVK) analysis to date. This asymmetry manifests itself in out-of-plane recoil of the final black hole and its inclusion in signal models is necessary both to measure this recoil, but also to accurately measure the full spin information of each black hole. We present the first model of the anti-symmetric contribution to the dominant co-precessing-frame signal multipole throughout inspiral, merger and ringdown. We model the anti-symmetric contribution in the frequency domain, and take advantage of the approximations that the anti-symmetric amplitude can be modelled as a ratio of the (already modelled) symmetric amplitude, and analytic relationships between the symmetric and anti-symmetric phases during the inspiral and ringdown. The model is tuned to single-spin numerical-relativity simulations up to mass-ratio 8 and spin magnitudes of 0.8, and has been implemented in a recent phenomenological model for use in the fourth LVK observing run. However, the procedure described here can be easily applied to other time- or frequency-domain models.

A non-eccentric BBH is characterised by the black-hole masses m 1 and m 2 , and each black hole's angular momentum, S i , which are usually represented in geometric units as the dimensionless vectors χ i = S i /m 2 i .The dominant effect of the spins on the GW signal is due to the spin components aligned with the binary's orbital angular momentum, which affect the rate of inspiral, and can therefore be measured through their effect on the signal's phase.The remaining (in-plane) spin components have little effect on the inspiral rate.They instead induce orbital and spin precession, which lead to modulations in the signal's amplitude and phase [30].In most cases this is a weaker contribution to the signal and more difficult to measure, in turn making it difficult to measure the full spin information of the binary.Spin misalignment also leads to an asymmetry in the power emitted above and below the orbital plane, and can lead to large out-of-plane recoils of the final black hole [31].Most signals to date have been too weak to observe precession and recoil (with the notable exception of several analyses of the signal GW200129 065458 signal -which we refer to as GW200129 in the rest of the text) [3,32,33], but more signals with measurable spin misalignment are expected as detector sensitivities improve.
Most current generic-binary models separately consist of a model of the signal in a non-inertial frame that tracks the precession of the orbital plane (a "co-precessing" frame), and a model of the time-or frequency-dependent precession angles.If the signal is decomposed into spin-weighted spherical harmonics, the dominant contributions in the co-precessing frame are the (ℓ = 2, |m| = 2) multipoles.As discussed in more detail in Sec.II, current Phenom and EOB models assume the symmetry h CP 22 = h CP * 2,−2 in the co-precessing frame.Precessing binaries also include an anti-symmetric contribution.There have been indications for some time that neglecting this contribution could lead to measurement biases [34,35], and more recently explicit examples of such biases have been found [36].One model that does include the anti-symmetric contribution is the NR surrogate model NRSur7dq4 [37], and this likely plays an important role in being able to accurately infer the primary spin in GW200129 [32,33], demonstrating the need to include the anti-symmetric contribution in Phenom and EOB models.
In this paper we present a simple method to model the antisymmetric contribution to the (ℓ = 2, |m| = 2) co-precessingframe multipoles, taking into account the phenomenology of how the anti-symmetric contribution depends on the in-plane spin direction and relates to the symmetric contribution.Note that all of the examples shown in this paper are constructed using either numerical relativity waveforms or post-Newtonian estimates and that the model introduced here is versatile in that it can be integrated into any frequency-domain approximant.
To motivate our focus on only the anti-symmetric contribution to the dominant multipoles, Fig. 1 shows the frequencydomain amplitude of the co-precessing-frame multipoles for a signal with mass-ratio q = m 1 /m 2 = 2, spin on the larger black hole of χ = 0.8, and spin misalignment with the orbital angular momentum of θ LS = 90 • , i.e., the spin initially lies entirely in the orbital plane.(This is case CF 38 in Ref. [38].)We see that the anti-symmetric (2, 2) amplitude is of comparable strength to the symmetric (3,3); since the (3, 3) extends to higher frequencies, it will contribute more power than the anti-symmetric (2, 2) at high masses, and comparable power arXiv:2310.16980v1[gr-qc] 25 Oct 2023 in low-mass systems.The next-strongest anti-symmetric contribution is to the (3,3), and we see that this is significantly weaker than the symmetric (4,4).This suggests that any model that includes symmetric contributions up to ℓ ≤ 4 need only include the dominant (2, 2) anti-symmetric contribution.If we wish to accurately model the signal to the level of the symmetric (5,5) contribution, then we must also include the anti-symmetric (3,3).Current Phenom models include symmetric multipoles up to ℓ = 4, and so we limit our attention to only the dominant anti-symmetric contribution.(Note that the anti-symmetric (3,3) multipole is also weaker than the symmetric (2,1) and (3,2) in this configuration.) We find that we can model the (2, 2) anti-symmetric contribution using numerical relativity (NR) simulations that cover only the reduced parameter space of the binary's mass ratio, the larger black hole's spin magnitude, and its misalignment angle; to a first approximation we do not need to sample the in-plane spin direction, which can be treated analytically.A more complete model that includes sub-dominant in-planespin-direction effects, and two-spin effects, is left to future work.FIG.1: Frequency-domain amplitude of co-precessing-frame symmetric and anti-symmetric contributions for a (q = 2, χ = 0.8, θ LS = 90 • ) configuration.The figure shows the symmetric contributions to (2, 2), (3,3), (4,4) and (5,5), and the anti-symmetric contributions to (2, 2) and (3,3).
In Sec.II we explain the motivation behind our modelling approach.We describe the preparation of the NR data that we used to calibrate our model in Sec.III.In Sec.IV we construct a model of the ratio of the anti-symmetric and symmetric amplitudes, and in Sec.V, we present our method to construct the anti-symmetric phase from the symmetric phase and the precession angle, α.We discuss the accuracy of our prescription in Sec.VI.

II. ASYMMETRY BACKGROUND
An aligned-spin binary is invariant under reflection across the orbital plane.If we choose a coordinate system where the orbital plane is the x-y plane and perform a decomposition of the gravitational-wave signal into spin-weighted spherical harmonics, then this symmetry arises in the signal multipoles as This relationship is useful when constructing a model of the multipoles of an aligned-spin binary: we need only explicitly model the positive-m multipoles, and the negative-m multipoles follow from Eq. ( 1).When the spins are not aligned with the orbital angular momentum, both the spins and the orbital plane may precess [30,39].In systems with mis-aligned spins Eq. ( 1) no longer holds, even if there is no orbital precession.The simplest illustration of this is the "superkick" configuration [31,40,41]: here the black holes are of equal mass and of equal spin, and the spins lie in the orbital plane but point in opposite directions.The symmetry of this system implies that the orbital plane does not precess, and although the spins will precess in the orbital plane, they both precess at the same rate and so remain oppositely directed to each other, so that the vector sum of the two spins is zero at all times.Although the direction of the orbital plane remains fixed, emission of linear momentum perpendicular to the orbital plane causes the entire system to move up and down.This linear momentum emission and "bobbing" [42] manifests itself in the gravitationalwave signal as an asymmetry in the positive-and negative-m multipoles, i.e., a violation of Eq. ( 1) [31].
The symmetry of Eq. ( 1) will remain broken regardless of any rotations performed on the multipoles [43].This point becomes important when constructing signal models, where we regularly make use of a "co-precessing frame".In this frame the signal during the inspiral can be approximated as that of a non-precessing binary [44] and so many current waveform models are split into a model for aligned-spin binaries and a model for the precession dynamics, and the precession dynamics are then used to "twist up" a non-precessing-binary waveform to produce the complete precessing-binary waveform [13,15,16,19,24,25,28,45,46].However, since the non-precessing-binary waveform respects the symmetry in Eq. ( 1), the model cannot reproduce the asymmetry that should be present in the true precessing-binary signal.
Several studies have considered the impact of neglecting these multipole asymmetries.Ref. [35] compares the multipoles from precessing-binary waveforms with those from nominally equivalent non-precessing binaries, to test a number of assumptions that go into the construction of many commonly used waveform models, including neglecting the multipole asymmetry.Ref. [34] compares NR waveforms from configurations with different in-plane spin directions and magnitudes, and argues that neglecting the multipole asymmetry may lead to parameter biases even at moderate SNRs, and that including the multipole asymmetry will be necessary to clearly measure in-plane spins and identify precessing systems.Finally, Ref. [36] uses the surrogate model NRSur7dq4 to identify the level of bias in binary measurement examples, and confirms that neglecting the multipole asymmetry leads to biases in measuring in-plane spins (but the masses and effective aligned spin χ eff are less affected).They also confirm the importance of the multipole asymmetry in precession measurements, showing that it had a significant impact on measurement of the properties of the signal GW200129 [32,33].
In the next section we will summarise the leading-order PN contribution to the asymmetry, which provides some insight into the phenomenology of the multipole asymmetry, and also motivate our modelling procedure.Although the multipole asymmetry has been known for some time, and indeed is included in the standard PN expressions that we use here, and is also discussed in detail in Ref. [43], we are not aware of any prior treatment that discusses the amplitude and phasing of the anti-symmetric (2,2) contribution in relation to the symmetric contribution, or notes the simple dependence of the relative phase between different in-plane spin directions, which is a key feature of the asymmetry that we exploit in constructing our model.

A. Inspiral
To gain insight into the phenomenology of the multipole asymmetry during the inspiral, we consider the leading-order post-Newtonian contributions to a binary where only one black hole is spinning and the spin lies entirely in the orbital plane.The binary consists of two black holes with masses m 1 and m 2 and the dimensionless spin on the primary is χ = S 1 /m 2 1 , where S 1 is the magnitude of the black hole's angular momentum.We use the post-Newtonian expressions from Ref. [47], where in this single-spin case the symmetric and anti-symmetric spin contributions are χ s = χ a = χ/2.The in-plane spin components incline the total angular momentum J with respect to the normal to the orbital plane (and direction of the Newtonian orbital angular momentum, L) by an angle ι, and the azimuthal precession angle of L around J is α; this is also the azimuthal angle of the total in-plane spin.As such, if we choose the instantaneous orbital plane to coincide with the x-y plane, then the entirely in-plane spin can be written as χ = χ(cos(α), sin(α), 0).
We start with the multipole h 22 as given in Eq. (B1) in Ref. [47].Requiring symmetry due to exchanging black holes (see the discussion prior to Eq. (4.15) in the same paper), leads to the relation h ℓm (Φ) = (−1) ℓ+m h * ℓ−m (Φ + π), where Φ is the orbital phase.We can enter the instantaneous orbital plane by setting ι = α = 0, and from our choice of spin we can then substitute χ ax = χ sx = χ cos(α)/2 and χ ay = χ sy = χ sin(α)/2; see sec.VI.B of Ref. [48] for details.We then have an approximation to the symmetric and anti-symmetric contributions to the signal in the co-precessing frame, h where the overall amplitude is A = √ 64π/5Mηv 2 /D L , M is the total mass, η = m 1 m 2 /M 2 is the symmetric mass ratio, δ = (m 1 − m 2 )/M, v is the relative speed of the two black holes, and D L is the luminosity distance to the source.Note that in Ref. [47] the symbol Φ denotes the orbital phase in the instantaneous orbital plane, but here we use it to denote the total orbital phase that enters in the waveform.
We may immediately note several important features from Eqs. ( 2) and (3).The spin does not enter the amplitude of the symmetric contribution at this order.The anti-symmetric contribution enters at O(v 2 ) lower than the symmetric contribution.We see that we may also consider the anti-symmetric amplitude |h CP,a 22 | as a simple rescaling of the symmetric amplitude, |h CP,s 22 |.
The in-plane spin direction, α, does not enter into the amplitude of the anti-symmetric contribution, but does modify the phase.The physical interpretation is that the phase of the anti-symmetric contribution depends on the direction of the in-plane spin relative to the separation vector of the two black holes.This will vary with the orbital phase, Φ, but also the (slower) precession rotation of the spin, given by α.This is also consistent with the observation in studies of out-of-plane recoil, that the recoil amplitude depends sinusoidally on the initial direction of the in-plane spin [31].
Finally, we note a key observation for the model that we will produce: if we modify the initial in-plane spin direction by ∆α, this will induce a simple overall phase shift in the anti-symmetric contribution, h CP,a 22 .This suggests that, given a set of single-spin numerical-relativity (NR) waveforms that cover the parameter space of mass ratio, aligned-spin magnitude and in-plane spin magnitude, we will have enough information to build a model of the anti-symmetric contribution to single-spin waveforms without the need to also sample multiple initial in-plane spin directions.We have just such a set of waveforms to hand, as used to construct the first NR-tuned full inpiral-merger-ringdown model of (the symmetric contribution to) precessing-binary waveforms [48], and discussed in Ref. [38].

B. Merger and ringdown
Before proceeding to construct a model, we consider the phenomenology of the anti-symmetric contribution through merger and ringdown, and inspect which inspiral features hold for the entire waveform.
One of the main features of the anti-symmetric contribution that we see in the leading-order inspiral single-spin expressions (2) and ( 3) is that an in-plane rotation of the spin by an angle ∆α results in a corresponding shift in the antisymmetric (2,2) phase by ∆α.This is evident from Fig. 2; the two configurations considered here correspond to the superkick configuration described earlier.⃗ S ⊥ 1 denotes the initial in-plane spin vector of the primary and ⃗ r 12 is the initial separation vector pointing from the primary to the secondary.It is clear that the asymmetry phase for a configuration with ⃗ S ⊥ 1 ⊥ ⃗ r 12 can be easily produced by applying a phase shift of π/2 to the anti-symmetric waveform of a configuration with ⃗ S ⊥ 1 ∥ ⃗ r 12 .We note that the simple phase relationship does not appear to hold as well through merger and ringdown, but the deviation is small enough that this could be due to numerical error, and requires more detailed study in future.
A second key feature of the anti-symmetric contribution is that its frequency is roughly half that of the symmetric contribution (plus a small correction from the spin-precession rate α).Fig. 3 shows the time-domain GW frequency of the symmetric and anti-symmetric contributions for a configuration where only the larger black hole is spinning, with the spin (χ = 0.7) entirely in-plane.We see that during the inspiral the anti-symmetric frequency is approximately half that of the symmetric, as we expect.
During ringdown the two frequencies are equal.This is consistent with our expectations from perturbation theory.In the ringdown, where perturbation theory results are applicable, the (ℓ = 2, m = ±2) multipoles will have the same (constant) frequency, and will decay exponentially at the  same rate.We therefore expect both the symmetric and antisymmetric combinations of h 22 and h 2,−2 to have the same frequency, and for the ratio of the symmetric and anti-symmetric amplitudes to be constant throughout the ringdown.
The third property we took from the inspiral expressions (2) and (3) was that the anti-symmetric amplitude can be considered as a rescaling of the symmetric amplitude.Fig. 4 illustrates that this holds for the entire waveform.It shows the frequency-domain ampiltude of the symmetric and antisymmetric (2,2) contributions for a configuration with (q = 1, χ = 0.4, θ LS = 60 • ), case CF 7 in Ref. [38].We see in par-ticular that the turnover to ringdown occurs at the same frequency (the ringdown frequency is at f M ∼ 0.09 for this configuration).We also see that, as discussed above, the symmetric and anti-symmetric contributions decay at the same rate.We will now exploit these features to construct a model for the anti-symmetric amplitude and phase.

C. Structure of the model
Based on the observations in the previous section we construct an approximate model of the anti-symmetric contribution to the (2,2) multipole in the co-precessing-frame as follows.We start with a model for the symmetric contribution, which provides us with the symmetric amplitude, A s ( f ), and symmetric phase, ϕ s ( f ).In the examples we consider here the symmetric contribution is calculated from NR waveforms.In a full model, we would start with the symmetric contribution from an already existing model.An explicit example is the multi-mode precessing-binary model described in Ref. [49], but the anti-symmetric model described in this paper can be applied to any existing frequency-domain precessing-binary model that separately provides A s ( f ), and symmetric phase, ϕ s ( f ), and the precession angle α( f ), as we describe below.
To construct the anti-symmetric amplitude, A a ( f ), we model the ratio κ( f ) = A a ( f )/A s ( f ).In the inspiral we use a post-Newtonian estimate of the amplitude ratio, and find that we can model the amplitude ratio accurately through to the ringdown by adding only one higher-order term, which we fit to our NR data.In the ringdown we treat the amplitude ratio as a constant, as motivated in the previous section.The amplitude model is presented in Sec.IV To construct the anti-symmetric phase, during the inspiral we combine the symmetric phase, ϕ s ( f ), and the precession angle α( f ), as prescribed by Eq. (3), i.e., ϕ a ( f ) ∼ ϕ s ( f )/2 + α( f ).In the merger-ringdown the anti-symmetric phase will behave as ϕ a ( f ) ∼ ϕ s ( f ).We apply an overall time and phase shift to smoothly connect these two functional forms at some transition frequency.We note that we find that it is possible to produce an accurate model of the antisymmetric phase ϕ a ( f ) using the same transition function for any binary, i.e., the parameters of the transition function do not need to be fit across the binary parameter space.The phase model is presented in Sec.V

III. NUMERICAL RELATIVITY DATA
The model is tuned to 80 single-spin NR waveforms generated using the BAM code [50,51].They cover mass ratios q = {1, 2, 4, 8}, spin magnitudes on the larger black holes of χ 1 = {0.2,0.4, 0.6, 0.8}, and angles of spin misalignment with the orbital angular momentum of θ LS = {30 • , 60 • , 90 • , 120 • , 150 • }.In labelling the configurations, the cases are ordered according to the mass ratios, then the spin magnitudes, then the misalignment angles; for example, case 57 corresponds to (q = 4, χ 1 = 0.8, θ LS = 60 • ).This is the same indexing as in Ref. [38], which provides full details on the catalog of simulations.To motivate the model and test our modelling assumptions we have also used families of simulations that consider variations in the initial in-plane-spin direction, based on those in Ref. [34].One notable addition to this set were simulations of superkick configurations where the black holes were given an additional out-of-plane momentum, to remove a secular drift in the centre-of-mass.
Several processing steps are performed to prepare the NR data for modelling, using the tools in Ref. [52,53].We start with the waveform data for the ℓ = 2 multipoles of the Weyl scalar, ψ 4,2m , in an inertial frame.We apply a Hann window to remove "junk" radiation from the beginning of the waveforms (a non-physical artefact of the initial data), and to remove numerical noise in the post-ringdown waveform.Furthermore, to ensure that the frequency-domain step-size is sufficiently small, the time-domain data are padded with zeros to the right.More details are given in Ref. [38].
We then transform to a frame that tracks Ĵ(t); as described in Ref. [48], this retains the approximation that the direction of Ĵ(t) is constant throughout the waveform.Modelling deviations from this approximation are left to future work, and are discussed further in Ref. [48], along with the procedure to track Ĵ(t) and perform this transformation.At the level of approximation and accuracy of the anti-symmetric model presented here, we do not expect this approximation to have any appreciable effect on the final model.We then transform the ψ 4,2±2 multipoles to a co-precessing frame, the "quadrupole-aligned" (QA) [54] or "optimal emission direction" [55,56] frame.In this frame the multipoles are significantly simplified, with the (l = 2, |m| = 2) multipoles having the strongest amplitude and the precession-induced modulations minimised.
The (2, ±2) multipoles of time-domain co-precessing-frame Ψ 4 are now transformed to the frequency domain, To compute the strain, we note that Ψ 4 (t) = − ḧ(t), so in the frequency domain we may write, where ω = 2π f .The anti-symmetric and symmetric components of the waveform in the QA frame are computed from The symmetric and anti-symmetric strains are complex quantities that can be written as and we can easily calculate their amplitude, A NR , and phase, ϕ NR , as their absolute value and argument, respectively.We denote the ratio between the anti-symmetric and symmetric amplitudes as

IV. MODEL OF THE AMPLITUDE RATIIO
We model the amplitude of the anti-symmetric contribution A a ( f ) as a ratio of the symmetric contribution A s ( f Our first step in constructing the ratio model is to compute the ratio in the framework of PN theory as a PN expansion in terms of v/c, where v is the relative velocity of the two black holes and c is the speed of light, and we choose geometric units where c = 1.We again restrict our analysis to singlespin binaries. To compute the PN ratio, we follow a procedure similar to the illustrative calculation in Sec.II.We obtain from Ref. [47] the complex 1.5PN expressions of the (ℓ = 2, |m| = 2) multipoles, h PN 2±2 , for spinning, precessing binaries with generic inclination angle ι moving on nearly circular orbits.The strains of the ℓ = |m| = 2 multipoles can then be transformed to a co-precessing frame that follows the instantaneous orbital plane.To achieve this, we choose a simple co-precessing frame where we set to zero the inclination angle of the orbital angular momentum relative to the total angular momentum, ι = 0, and also set the precession angle, α = 0. We then use an approximate reduction to a single-spin system [48], where χ s = (χ 1 + χ 2 )/2, χ a = (χ 1 − χ 2 )/2 are the symmetric and anti-symmetric spins defined in Ref. [47] and θ LS is the inclination of the spin from the orbital angular momentum vector.
We then find that the symmetric and anti-symmetric amplitudes are, The ratio of the two amplitudes is then given by, where η = m 1 m 2 /M 2 is the symmetric mass ratio and δ = (m 1 − m 2 )/M > 0 is a fractional mass difference.The expression of the PN ratio of the anti-symmetric amplitude over the symmetric amplitude depends on the symmetric mass ratio, η, the spin magnitude, χ, and the angle θ LS of the system.Consequently, Eq. ( 15) can be used to compute the PN ratio of any configuration as a function of frequency since v = (π f M) 1/3 .We cannot expect the PN estimate of the amplitude ratio, Eq. ( 15), to be accurate through merger and ringdown.To capture this behaviour we investigated the addition of unknown higher-order terms and fit their coefficients to the NR data.The simplest approach is to add only one additional term, for example, where b is fit to the NR data.To do this it was necessary to choose an appropriate frequency range over which to perform the fit.The NR data are noisy in the early inspiral and in the post-ringdown frequencies, so we first identified a consistent definition of a frequency range that could be used for all 80 NR simulations, based on the ringdown frequency, f RD of each NR configuration.The frequency range that we used in the final fits was given by M f min = (M f RD − 0.01) /5 and M f max = M f RD − 0.002.We investigated fits to the amplitude ratio of the form Eq. ( 16) with n = 3, 5, 7. Fig. 5 shows the results for one configuration, and illustrates that n = 5 provides the most accurate fit.In fact, we find that n = 5 consistently provides the most accurate fit across all 80 NR simulations.0.00 0.02 0.04 0.06 0.08 0.10 0.12 Mf FIG. 5: Amplitude ratio for the (q = 1, χ = 0.4, θ LS = 60 o ) configuration, for the NR data, the PN ratio κ PN , and fits to higher-order corrections as in Eq. ( 16), with n = 3, 5, 7. The vertical dashed line indicates the ringdown frequency.
It is also clear from Fig. 5 that our fit is not accurate beyond the ringdown frequency.Beyond this point the amplitude ratio appears to plateau.This is consistent with our expectations from perturbation theory, as discussed in Sec.II B, which pre-dicts that the amplitude ratio will be constant throughout ringdown.To include this feature in our model, we fix κ( f ) to the value κ( f RD ) at frequencies f > f RD .To avoid a sharp transition we use a moving average algorithm such that, We use a symmetric window of an equal number of points (k = 40) on either side of the frequency f to calculate the moving average.Here N is the length of the frequency series and the algorithm updates κ( f ) for 40 ≤ n ≤ N − 40.
We fit the coefficient b in Eq. ( 16) (with n = 5) to each of the 80 single-spin NR waveforms from the NR catalogue in Ref. [38].Fig. 6 shows the resulting fit for the amplitude ratio for the (q = 1, χ = 0.4, θ LS = 60 • ) configuration.We see that the final fit, including the levelling off of the amplitude ratio through ringdown, agrees well with the NR data.The lower panel of the figure shows the resulting estimate for the anti-symmetric amplitude.We see that this agrees well with the NR data up to the point where the NR amplitude becomes dominated by noise.
The fit cannot reproduce the NR amplitude ratio in all cases.In many cases the NR amplitude ratio oscillates during the inspiral.Fig. 7 shows an extreme example of this feature, from the (q = 1, χ = 0.8, θ LS = 30 • ) configuration.It is not clear what causes these oscillations.Oscillations in the coprecessing-frame amplitude ratio during the inspiral can be due to the choice of co-precessing frame, as we will discuss later.However, if that were the cause of the oscillations in the NR data, we would expect there to be some correlation with the degree of precession in the configuration.Instead we do not find any relationship between the amplitude of the oscillations, which in some cases are negligible (as in Fig. 6), and in others (like Fig. 7) lead to significant variations in the amplitude ratio.We do find, though, that our model passes through a reasonable estimate of the mean of the oscillations.The largest impact is on the constant value of the amplitude that the model settles on for the ringdown regime; in the example in Fig. 7 the model's estimate of the anti-symmetric amplitude during the ringdown is roughly 20% below the NR value.
Finally, in some cases we also found that the amplitude ratio in the NR data did not level off during the ringdown.Fig. 8 is an example of this.We have not been able to determine the reason for this.As noted previously, we expect that since the (2, 2) and (2, −2) multipoles decay at the same rate, that the ratio between their amplitudes would remain constant during the ringdown.It is possible that this effect is obscured by numerical noise.Regardless of the cause, and in the absence of any compelling explanation for alternative behaviour, we have chosen to impose the behaviour that we expect from perturbation theory.
The coefficient b is shown as a function of symmetric mass ratio and misalignment angle in Fig. 9; the values for all four spin magnitudes are shown together.This plot illustrates the general trend of b across the parameter space.We find that there is no general trend with respect to the spin magnitude, and this is illustrated more clearly in Fig. 10, which shows b as 0.00 0.02 0.04 0.06 0.08 0.10 0.12 Mf 0.00 0.02 0.04 0.06 0.08 0.10 0.12 Mf A a (f) a function of symmetric mass ratio, for each spin magnitude, for the subset of cases with θ LS = 90 • .The PN amplitude ratio already includes a linear dependence on the spin magnitude, and given the uncertainty in our fits, we do not attempt to identify a higher-order spin dependence.We motivate this point further in Sec.VI (Fig. 14).We then include all simulations into a two-dimensional parameter-space fit of the form, where we find b 0 = 18.0387, b 1 = 15.4509,b 2 = 55.1140 and b 3 = −203.6290.From Eq. ( 18), we notice that b does not go to zero when θ LS is 0 • or 180 • .However, the presence of the sin θ LS term at the numerator of the ansatz of the ratio model in Eq. ( 15) ensures that the multipole asymmetry goes to zero at these points.The amplitude as predicted by this fit is shown on each of our figures and labelled as "final model".Two caveats to this approach are worth noting.One is the choice of co-precessing frame.Previous work has shown that the symmetric (2,2) contribution takes a simple form in the QA co-precessing frame; indeed, the amplitude evolution can be approximated by the (monotonic) amplitude of an equivalent aligned-spin configuration.This is not necessarily the case for the anti-symmetric contribution, and this is 0.00 0.02 0.04 0.06 0.08 0.10 0.12 Mf 0.00 0.02 0.04 0.06 0.08 0.10 0.12 Mf A a (f) one possible cause of the oscillations that we see (although, as noted previously, it does not show a clear correlation with the strength of precession).Conversely, we found that the amplitude evolution of the anti-symmetric (2,2) multipole was monotonic for PN waveforms if we choose ι = α = 0 in their construction (which is equivalent to choosing a co-precessing frame that tracks the Newtonian orbital angular momentum, i.e., the normal to the orbital plane), but if we consider the PN waveforms in the QA frame then the anti-symmetric (2,2) amplitude shows strong modulations.This illustrates that the anti-symmetric component can depend strongly on the choice of co-precessing frame, and although we do not expect this to be a significant issue at the level of accuracy or approximation of the current model, it may require better understanding in future refinements of anti-symmetric models.
The second point is that our model is constructed based on the phenomenology of single-spin binaries.If the model is to be used for generic binaries, one much choose a method to treat two-spin configurations.One option, which is employed in the PhenomX implementation [48,49], is as follows.
We can describe the spin of an equivalent single-spin system using Eqs.( 16) and ( 17) in Ref. [48], but diverge slightly 0.00 0.02 0.04 0.06 0.08 0.10 0.12 Mf 0.00 0.02 0.04 0.06 0.08 0.10 0.12 Mf A a (f) in the definition of the in-plane spin as follows: where χ p is the effective precession spin as defined in Ref. [57].The anti-symmetric amplitude for an equal-mass binary with both spins equal in magnitude, entirely in-plane and pointing in the same direction, must drop to zero.The symmetric in-plane spin combination of Eq. ( 15) in Ref. [48] cannot account for this effect in superkick configuration.Therefore, we can instead use an anti-symmetric in-plane spin combination, to appropriately map two-spin to single-spin systems for generating the anti-symmetric waveform.

V. PHASE MODEL
GW phases for chirp signals in the frequency domain are typically quite featureless and therefore not conducive to modeling.Following a suite of models for the symmetric phase [11,12,17], we focus first on the anti-symmetric phase derivative.Fig. 11 demonstrates that the anti-symmetric phase derivative (i.e., frequency) behaviour in the time-domain as discussed in Sec.II B is preserved in the frequency domain as well.Remarkably, we find that it is possible to construct a map of the symmetric phase derivative to the anti-symmetric phase derivative that is independent of the binary's parameters.Therefore, we do not need to produce any parametric fits for the map over the intrinsic parameter space of BH binaries, which makes this model extremely simple.
Our model of the anti-symmetric phase derivative is defined The specific configuration shown here is for a single-spin binary with (q = 1, χ = 0.2, θ LS = 30 • ).
by the piecewise function, where the phase derivative in the intermediate region is given by, As is evident from Fig. 11, the functional form of region I transitions to the functional form of region II in the intermediate region.To ensure smooth transition in the intermediate region an obvious choice would be a tanh windowing function.
Noting that a tanh function can be computationally inefficient during model evaluation, we instead use the Taylor expansion of the tanh function up to second order with appropriate normalization to construct the phase derivative functional form of Eq. ( 22).The phase derivative ansatz was calibrated to NR simulations by treating the transition frequency, f T , the width of the transition window, f w , and the phase coefficients (A 0 and B 0 ) as free parameters.The window parameters were not particularly sensitive to the tuning and the minor variations could be attributed to the noise in the anti-symmetric phase derivatives; f T = 0.85 f RD and f w = 0.005 was found to be an optimal choice across the entire single-spin parameter space.
Once the parameters of the window function were fixed, we investigated the impact of fixing the phase coefficients.A best fit to data for B 0 = 0 was found to be consistent across the parameter space; A 0 on the other hand showed some variation, but no specific trend.Furthermore, choosing the alge-braic mean of A 0 for the set of 80 simulations did not significantly impact the quality of the fit.This shows (1) the fitting algorithm tries to find the best A 0 for continuity of the phase derivative at f T , and (2) the variation in optimal A 0 across the parameter space is more likely due to noisy data and not a function of intrinsic parameters of the binary.
Applying a shift to the phase derivative is equivalent to an overall time shift of the waveform.We exploit this freedom by fixing the symmetric phase derivative minima to be 0 at f RD .This imposes, Figure 12 shows that A 0 obtained from the fitting algorithm and from NR data using Eq. ( 23) reasonably well for a majority of the cases (cf.y-axis on Fig. 11).The anti-symmetric model will be used with a phenomenological symmetric waveform model, so an A 0 derived from the symmetric waveform model makes the phase construction self-consistent and robust.Furthermore, Fig. 12 highlights that the gain in accuracy by making a model to capture the near-stochastic behaviour of A 0 may be outweighed by errors introduced by over-modelling.As such, our model of the anti-symmetric phase is NR-informed but, noting that the data for the antisymmetric waveform is often close to numerical noise, we prioritised our understanding of the physics rather than model optimization.The phase of the anti-symmetric waveform is obtained by integrating the two pieces, where the integration constant ϕ B0 is determined by continuity at f T Finally, the phase of the asymmetry is modulated by the inplane spin direction α; therefore, the initial phase, ϕ A0 , is the value of α at a reference frequency i.e., ϕ A0 = α( f ref ).

VI. MODEL ACCURACY
A standard measure of waveform accuracy used extensively in the literature is the match of the waveform model with NR waveforms, defined as, where h( f ) = h + ( f ) − ih × ( f ) is a complex frequency sequence constructed from the two polarizations of the waveform.As such, calculating matches of just the anti-symmetric waveform is not physically meaningful.Furthermore, a true measure of performance of precessing waveforms in data analysis can only be obtained by calculating matches of the full waveform in the inertial frame, making considerations for precession as well as extrinsic parameters.Therefore, matches of just the anti-symmetric waveform in the co-precessing frame provide some indication of the accuracy of one ingredient in the full waveform, but do not indicate the overall accuracy of the corresponding precessing waveform.However, an inner product like that in Eq. ( 26) is a useful measure of agreement between two complex frequency series.Since a match-like calculation for the anti-symmetric contribution in the co-precessing frame cannot be interpreted in terms of either signal detection efficiency or parameter measurement accuracy, there is no reason to include the detector sensitivity, and so we will use a simpler inner product of the form, where f 1 is the starting frequency of the NR waveform in geometric units, M f 1 = 0.02, and M f 2 = 0.15, after which point the amplitude of the NR waveform is below the noise floor of the data.We consider normalised waveforms, ĥ = h/ √ ⟨h|h⟩, so that the maximum value of the inner product is one.We used the standard implementation of this inner product in pycbc [58], a python software package for GW data analysis, for our match computations.We then consider the "mismatch" between the two waveforms, In Fig. 13 we show the mismatches of the anti-symmetric waveform constructed from our model with the 80 NR waveforms that were used to calibrate the model.To determine the accuracy of the individual components, we also computed matches of the amplitude (phase) model complemented by phase (amplitude) constructed from NR data, with the full anti-symmetric waveform constructed from NR data.The overall accuracy of both the models are comparable.We note that although the model was verified using the same waveforms as used for modelling, since the NR tuning was relatively simple -i.e., a single co-efficient in Eq. ( 16) fit to the four-parameter ansatz Eq. ( 18) across a two-dimensional parameter space -using a much smaller subset of waveforms would have produced a model of similar accuracy, and the simplicity of this model and the single-spin parameter space obviates any concerns about over-fitting or unexplored corners of parameter space.
To investigate the quality of the surface fit in Eq. ( 18), we also computed mismatches for the amplitude model using fit coefficients b of Eq. ( 16).As is evident from Fig. 14, for most cases the performance is unchanged and for the handful of cases where the mismatch changes, the difference is not very significant.This further illustrates that capturing the nonlinear dependence on spin magnitude is unlikely to make significant improvement to the amplitude model.In addition to the argument made for using Eq. ( 23) to calculate A 0 from the symmetric waveform and precession angle, α, we calculated mismatches for the different choices of A 0 in the phase model -i.e., A f it 0 and A equation 0 in Fig. 12 -to confirm that there was no reduction in model accuracy.
Note that the anti-symmetric waveform model is downstream from the symmetric waveform model as well as the precession angle models.Therefore, enhancement in performance of the overall model due to the addition of an asymmetry model must always be discussed in the context of the underlying symmetric, precessing waveform model.This is beyond the scope of present work and will be discussed in the context of the current generation frequency-domain precessing-binary PhenomXPNR model [49].

VII. CONCLUSION
We have presented a method to model the anti-symmetric contribution to the (ℓ = 2, |m| = 2) multipoles in the coprecessing frame.This is a general approach that can be applied to any frequency-domain model that provides the symmetric contribution to the (2,2) amplitude, phase and precession angle, α.We expect that the main insights of this model can be used to easily generalise the procedure to the timedomain and be useful in including multipole asymmetry in current generation time-domain models [45,46].We summarise the key insights here.
For the amplitude, we observe that the anti-symmetric amplitude can be easily modeled as a rescaling of the symmetric amplitude.As shown in Fig. 4, both amplitudes have the same basic structure, in particular the same ringdown frequency and decay rate.In the inspiral our model of the amplitude ratio is based on a PN expression for single-spin systems, and in the late inspiral and merger a higher-order correction to the PN expression is calibrated to 80 NR simulations of single-spin binaries [38].In the ringdown we make use of the prediction from perturbation theory that the ratio of the symmetric and anti-symmetric amplitudes will be constant.The amplitude model ratio is presented in Eqs. ( 15), ( 16) and the fit to the higher-order contribution given in Eq. (18).
For the phase, we use the facts that in the inspiral the frequency of the anti-symmetric contribution equals the orbital frequency plus the spin precession frequency, and that during the ringdown the symmetric and anti-symmetric frequencies are the same.We are able to construct a mapping from the symmetric phase and precession angle, ϕ s and α, to the antisymmetric phase, ϕ a , motivated by the 80 single-spin simula- tions, but without the need of any explicit tuning.The model of the phase is given in Eqs. ( 23), ( 24) and ( 25).An additional crucial observation is that a rotation of the initial in-plane spin direction of ∆α introduces a corresponding shift of ∆α into the anti-symmetric phase.It is this observation that allows us to model the dependence on in-plane-spin direction, even though we do not have a set of NR simulations that span this subspace.
The procedure presented here is not in itself a signal model.As already noted, one must also provide the symmetric amplitude, phase and precession angle.In addition, having constructed a model for the anti-symmetric contribution to the (2,2) multipoles in the co-precessing frame, one must then "twist them up" to produce the multipoles in the inertial frame.This has been done for the recent extension of the multimode frequency-domain precessing-binary IMRPhenomXPNR model [49].This is the first phenomenological full inspiral-mergerringdown model of the anti-symmetric multipole contributions, and there are many directions for improvement and issues to be resolved.The first limitation of the model is that it is based only on single-spin binaries.We expect that generic two-spin systems can in most cases be modelled to a good approximation by single-spin systems, but given that the antisymmetric contribution is itself a weak effect, it is unclear how well this approximation can be used.It would be useful to study the applicability of the single-spin approximation for the anti-symmetric contribution, and, indeed, to extend the model to two-spin systems.The anti-symmetric model is also limited to the (ℓ = 2, |m| = 2) multipoles.We argue in Fig. 1 that the (2,2) contribution will be sufficient for most applications, since the anti-symmetric amplitude is far weaker than the symmetric, but for high signal-to-noise-ratio (SNR) systems, higher-order anti-symmetric multipoles will ultimately be required.
When decomposing the signal into a co-precessing frame and corresponding time-or frequency-dependent angles, one will find (depending on the co-precessing frame chosen) oscillations in the amplitude and/or phase of the anti-symmetric contribution; recall that the co-precessing-frame maps exactly to a non-precessing aligned-spin waveform only for the leading-order quadrupole terms [44,48,49].In this work we have removed oscillations in the inspiral amplitude by simply setting the inclination angle ι and precession angle α to zero in the PN expressions we have used for the signal multipoles [47]; oscillations remain in the co-precessing-frame signals for the NR waveforms, but our model consists of only a monotonic fit through these oscillations.A better understanding of these oscillations, or at least a model that captures them, is a necessary next step in modelling the multipole asymmetries.
We have assumed that the only affect of the initial in-plane spin direction is to introduce an overall offset in the antisymmetric phase.This approximation will not be exact during the merger-ringdown.In particular, in the ringdown we have made the approximation that the amplitude ratio will be independent of the initial in-plane spin direction.However, the amplitude ratio will be determined by the relative phase of the

2 FIG. 3 :
FIG.3:The symmetric and anti-symmetric frequencies obtained from NR data of a (q = 2, χ = 0.7, θ LS = 90 • ) configuration.During inspiral the anti-symmetric frequency, ω a , is about half of the symmetric frequency, ω s (or nearly equal to the orbital frequency), and close to merger quickly catches up with the symmetric frequency.As expected from perturbation theory, ω s = ω a during ringdown.

FIG. 11 :
FIG.11:The anti-symmetric phase derivative (thick black line), the symmetric phase derivative (grey dashed line) and the combination of half of the symmetric phase derivative with the derivative of the precession angle α (thin grey line).The specific configuration shown here is for a single-spin binary with (q = 1, χ = 0.2, θ LS = 30 • ).

FIG. 12 :
FIG. 12: Comparison of phase coefficient A 0 , obtained from fitting the phase ansatz to NR data (light purple triangles) across the parameter space and from NR data using Eq.(23) (blue circles), for cases with χ = 0.4 and 0.8.The algebraic mean of the set of coefficients are shown by horizontal lines in corresponding colors.The case index is as described in Sec.III.

FIG. 13 :
FIG.13: Mismatches of the anti-symmetric waveform model in the co-precessing frame with NR data.The black triangles show the mismatches for the combined amplitude and phase model while the magenta squares (cyan circles) show the mismatches for just the amplitude (phase) model with the phase (amplitude) constructed from NR data.The dashed magenta and cyan lines show the average mismatch for the amplitude and phase model, respectively; on average they are of comparable accuracy.The black solid line shows the average mismatch for the overall model.The x-axis denotes the case index of the NR simulation as usual i.e., five different θ LS for each spin magnitude shown in figure, for q = 1,2,4 and 8.
⃗ r 12 in dashed black and ⃗ S ⊥ 1 ⊥ ⃗ r 12 in solid grey, shows a constant phase offset of π/2.Inset: Dashed blue line shows that the anti-symmetric waveform for ⃗ S ⊥1 ⊥ ⃗ r 12 can be constructed by just applying a π/2 phase shift to the ⃗ S ⊥ 1 ∥ ⃗ r 12 waveform even in the strong-field regime close to merger (t merger = 1784M).