Possible π/2 phase shift missing in PhenomPv2

[SSL32081]

It appears that the PhenomPv2 waveforms are flipped along the x-axis as compared to LALSim, and even PhenomD.
Since this model consists of only the dominant modes, this issue could be due to either a π/2 phase shift or a flip in sign, which will require further investigation.

Here is a figure to illustrate the issue:

Blue is the waveform computed with Ripple, and the orange dashed one is from LALSim.
On the right is the results from PhenomD, and PhenomPv2 is on the left.
Both columns were injected with the same parameters (non-precessing), and it is clear that the orange ones are consistent in both columns whereas the blue one is flipped.
The legend also shows the (frequency-domain) match; in the PhenomPv2 case, it is a -1, a perfect anti-phase.

I am also attaching the MWE to produce this result.

import jax
import jax.numpy as jnp
jax.config.update("jax_enable_x64", True)

from jimgw.single_event.utils import Mc_eta_to_m1_m2
from jimgw.single_event.waveform import RippleIMRPhenomPv2, RippleIMRPhenomD

import numpy as np
from scipy.interpolate import interp1d
from scipy.integrate import simpson

import lalsimulation as lalsim
from lalsimulation import IMRPhenomD, IMRPhenomPv2
from lal import MSUN_SI, PC_SI

MPC_SI = 1e6 * PC_SI

import matplotlib.pyplot as plt


def waveform_fd_to_td(waveform_fd):
    '''
    Convert frequency domain waveform to time domain
    '''
    return np.fft.irfft(waveform_fd) * sampling_frequency

freqs, raw_psd = np.loadtxt('test_jim_psd/GW150914_psd_H1.dat', unpack=True)
interp_psd = interp1d(freqs, raw_psd, 
                      fill_value=(raw_psd[0], raw_psd[-1]), 
                      bounds_error=False)(conjugate_frequencies)

def noise_weighted_inner_product(frequencies, waveform_1, waveform_2):
    '''
    Convert frequency domain waveform to time domain
    '''

    numerator = np.conj(waveform_1) * waveform_2 
    return 2 * simpson(
        numerator / interp_psd,
        x=frequencies,
    )

def FD_match(frequencies, waveform_1, waveform_2):
    norm1 = noise_weighted_inner_product(frequencies, waveform_1, waveform_1).real * 2
    norm2 = noise_weighted_inner_product(frequencies, waveform_2, waveform_2).real * 2
    inner_product = noise_weighted_inner_product(frequencies, waveform_1, waveform_2).real * 2
    return inner_product / np.sqrt(norm1 * norm2)


# Prepare conjugate time and frequency arrays
sampling_frequency = 2048.
duration = 8.
delta_t = 1 / sampling_frequency
N_samples = int(jnp.round(duration * sampling_frequency))
N_freqs = int(jnp.round(N_samples / 2) + 1)
time_array = jnp.linspace(0.0, duration - delta_t, N_samples)
conjugate_frequencies = jnp.linspace(0.0, sampling_frequency / 2, N_freqs)

f_min = 10.
f_max = sampling_frequency / 2
reference_frequency = 50.0  # Hz
frequency_mask = (conjugate_frequencies >= f_min) & (conjugate_frequencies < f_max)
masked_frequencies = conjugate_frequencies[frequency_mask]

# Injection parameters
injection_parameters = {
    'M_c': 30.0,
    'eta': 0.20,
    's1_x': 0.0,
    's1_y': 0.0,
    's1_z': 0.1,
    's2_x': 0.0,
    's2_y': 0.0,
    's2_z': 0.2,
    'd_L': 1.0,
    'luminosity_distance': 1.0,
    'phase_c': 0.0,
    'phase': 0.0,
    'iota': 1.3,
}

m1, m2 = Mc_eta_to_m1_m2(
    injection_parameters['M_c'],
    injection_parameters['eta'] 
)
injection_parameters['mass_1_SI'] = float(m1) * MSUN_SI
injection_parameters['mass_2_SI'] = float(m2) * MSUN_SI
injection_parameters['luminosity_distance_SI'] = \
    injection_parameters['luminosity_distance'] * MPC_SI

# Prepare Jim/Ripple waveforms
PhenomPv2 = RippleIMRPhenomPv2(f_ref=reference_frequency)
PhenomD = RippleIMRPhenomD(f_ref=reference_frequency)

ripple_PhenomD_hpc = PhenomD(masked_frequencies, injection_parameters)
ripple_PhenomPv2_hpc = PhenomPv2(masked_frequencies, injection_parameters)

ripple_D_hp, ripple_D_hc = np.zeros((2, N_freqs), dtype=np.complex128)
ripple_D_hp[frequency_mask] = ripple_PhenomD_hpc['p']
ripple_D_hc[frequency_mask] = ripple_PhenomD_hpc['c']
ripple_Pv2_hp, ripple_Pv2_hc = np.zeros((2, N_freqs), dtype=np.complex128)
ripple_Pv2_hp[frequency_mask] = ripple_PhenomPv2_hpc['p']
ripple_Pv2_hc[frequency_mask] = ripple_PhenomPv2_hpc['c']

# Prepare LALSimulation waveforms
parameters = [injection_parameters[key] for key in(
    'mass_1_SI', 'mass_2_SI', 's1_x', 's1_y', 's1_z',
    's2_x', 's2_y', 's2_z', 'luminosity_distance_SI', 'iota', 'phase',
)]
lal_PhenomD_hp, lal_PhenomD_hc = lalsim.SimInspiralChooseFDWaveform(
            *parameters,
            0, 0.0, 0.0,
            1/duration, f_min, f_max,
            reference_frequency, None, IMRPhenomD,
        )
lal_PhenomPv2_hp, lal_PhenomPv2_hc = lalsim.SimInspiralChooseFDWaveform(
            *parameters,
            0, 0.0, 0.0,
            1/duration, f_min, f_max,
            reference_frequency, None, IMRPhenomPv2,
        )

print('PhenomD, hp:', FD_match(conjugate_frequencies, ripple_D_hp, lal_PhenomD_hp.data.data))
print('PhenomD, hc:', FD_match(conjugate_frequencies, ripple_D_hc, lal_PhenomD_hc.data.data))
print('PhenomPv2, hp:', FD_match(conjugate_frequencies, ripple_Pv2_hp, lal_PhenomPv2_hp.data.data))
print('PhenomPv2, hc:', FD_match(conjugate_frequencies, ripple_Pv2_hc, lal_PhenomPv2_hc.data.data))

fig, Axes = plt.subplots(2, 2, figsize=(9, 3*2), sharex=True, sharey='row',
                         constrained_layout=True)

for ax, (jim_wfm, lal_wfm) in zip(
    Axes.T.flatten(), 
    [
        (ripple_D_hp, lal_PhenomD_hp.data.data),
        (ripple_D_hc, lal_PhenomD_hc.data.data),
        (ripple_Pv2_hp, lal_PhenomPv2_hp.data.data),
        (ripple_Pv2_hc, lal_PhenomPv2_hc.data.data),
    ]
):
    match = FD_match(conjugate_frequencies, jim_wfm, lal_wfm)
    jim_td_wfm = waveform_fd_to_td(jim_wfm)
    lal_td_wfm = waveform_fd_to_td(lal_wfm)
    ax.plot(time_array, np.roll(jim_td_wfm, -100), label='Ripple')
    ax.plot(time_array, np.roll(lal_td_wfm, -100), ls='--', label=fr'${{\cal M}} = {match:.3f}$')
    ax.set_xlabel(r'$t \,/\,{\rm s}$')
    ax.grid(False)
    ax.tick_params(axis='both', which='both',
                   top=True, bottom=True, left=True, right=True)
    ax.legend(loc='upper left')

Axes[0, 0].set_title('PhenomD')
Axes[0, 1].set_title('PhenomPv2')
Axes[0, 0].set_ylabel(r'$h_+(t)$')
Axes[1, 0].set_ylabel(r'$h_\times(t)$')
ax.set_xlim(7.5, 8)

plt.savefig('compare_ripple_lal_Pv2.png', dpi=300)

[kazewong]

My memory is vague on this one so I need to double check later but I don't think lalsim takes phase of coalescence as an input to their waveform generator https://lscsoft.docs.ligo.org/lalsuite/lalsimulation/group__lalsimulation__inspiral.html. If I remember correctly this is handled down stream in constructing this likelihood. @SSL32081 Can you double check if that is the case?

[kazewong]

At some point I think this discussion sprung up in the past, perhaps it is useful to look around this commit to see what's going on b7e6b41

[SSL32081]

My memory is vague on this one so I need to double check later but I don't think lalsim takes phase of coalescence as an input to their waveform generator https://lscsoft.docs.ligo.org/lalsuite/lalsimulation/group__lalsimulation__inspiral.html. If I remember correctly this is handled down stream in constructing this likelihood. @SSL32081 Can you double check if that is the case?

I think it is the peak time shift that is being handled in likelihood. But the phase is used in lalsimulation, at least when I feed a phase with the π/2 shift to SimInspiralChooseFDWaveform, I get the perfect match again.

I am guessing this was not noticed because the match was mostly computed with an absolute value?

[SSL32081]

I believe I have located the issue, adding the first few lines to PhenomPv2_utils to catch for aligned-spin BBHs solves the issue:

    no_inplane_J0 = jnp.abs(J0x_sf) < MAX_TOL_ATAN and jnp.abs(J0y_sf) < MAX_TOL_ATAN
    phiJ_sf = jax.where(
        no_inplane_J0,
        jnp.pi/2. - phiRef,
        jnp.arctan2(J0y_sf, J0x_sf)
    )

It is there in lalsimulation, but not found in the current ripple implementation.

Since this is only an issue with aligned-spin binaries, I tried to check the precessing cases as well. And as expected, even without the fix, the matches are perfect.

If you want, I can make a quick PR to fix this.