bilby-dev / bilby

A unified framework for stochastic sampling packages and gravitational-wave inference in Python. Note that we are currently transitioning from git.ligo.org/lscsoft/bilby, please bear with us!
https://bilby-dev.github.io/bilby/
MIT License
57 stars 63 forks source link

Phase Marginalization Inconsistency in Tutorial #807

Open bilby-bot opened 3 months ago

bilby-bot commented 3 months ago

In GitLab by @git.ligo:jacob.golomb on Jun 25, 2024, 22:28

@jake.summers noticed that this line and this line in two data tutorials both say that phase marginalization is formally invalid for IMRPhenomPv2, but one of these tutorials turns it on and the other turns it off.

Isn't this only invalid for waveforms with HMs? It seems like phase marginalization should be able to be enabled in both tutorials since they use IMRPhenomPv2.

bilby-bot commented 3 months ago

In GitLab by @git.ligo:michael.williams on Jun 26, 2024, 15:20

My understanding is that phase marginalization is not formally valid when using a precessing waveform, however, if the effect of precession is small then it should approximately hold.

bilby-bot commented 3 months ago

In GitLab by @git.ligo:jacob.golomb on Jun 26, 2024, 17:38

It was my understanding that the Bessel function in the phase marginalized likelihood comes from being able to analytically integrate over the likelihood when you have the phasing terms that go like $e^{im\phi_{\rm orb}}$ for $m = \pm 2$. I agree that precession should excite HMs, but I would think the assumptions made in the marginalized likelihood means that it is valid for any waveform that only includes $\mathcal{l} = |m| = 2$ modes (such as IMRPhenomPv2).

The reason I want to clarify this in the tutorial is because new users use the instructions in the tutorial for recommendations for what settings to use when running bilby.

bilby-bot commented 3 months ago

In GitLab by @git.ligo:colm.talbot on Jun 26, 2024, 22:04

Precession induces mode mixing (2|2| -> 2{|1|,0}) which breaks the approximation. I'm not sure if this is actually the case for Pv2, I would have to double check, the easiest way to check is to see if the waveform rotates like e^{2iphi}.

bilby-bot commented 3 months ago

In GitLab by @git.ligo:jacob.golomb on Jun 26, 2024, 23:19

I don't think the twisting-up procedure in Pv2 mixes modes between |m|=2 and |m|=/=2. I'm pretty sure it just sums over the |m|=2 modes into new |m|=2 modes.

bilby-bot commented 3 months ago

In GitLab by @git.ligo:jacob.golomb on Jun 26, 2024, 23:23

import bilby
import numpy as np
import matplotlib.pyplot as plt
mass_1 = 50
mass_2 = 30

injection_parameters = dict(
    mass_1=mass_1,
    mass_2=mass_2,
    a_1=0.2,
    a_2=0.1,
    cos_tilt_1=0.5, 
    cos_tilt_2=0.1,
    luminosity_distance=10.0,
    theta_jn=0.4,
    psi=2.659,
    phase=0,
    geocent_time=1126259642.413,
    ra=1.375,
    dec=-1.2108,
)

duration = 6
sampling_frequency = 2048
start_time = injection_parameters["geocent_time"] + 2 - duration

waveform_arguments = dict(
    waveform_approximant="IMRPhenomPv2",
    reference_frequency=20.0,
    minimum_frequency=20.0,
)

waveform_generator = bilby.gw.WaveformGenerator(
    duration=duration,
    sampling_frequency=sampling_frequency,
    frequency_domain_source_model=bilby.gw.source.lal_binary_black_hole,
    parameter_conversion=bilby.gw.conversion.convert_to_lal_binary_black_hole_parameters,
    waveform_arguments=waveform_arguments,
)

interferometer = bilby.gw.detector.get_empty_interferometer("H1")
interferometer.minimum_frequency = 20
interferometer.sampling_frequency=sampling_frequency
interferometer.duration=duration
interferometer.start_time=start_time

strain_phase0 = waveform_generator.frequency_domain_strain(injection_parameters)
detector_strain_phase0 = interferometer.get_detector_response(strain_phase0, parameters = injection_parameters)

new_injection_parameters = injection_parameters.copy()
new_injection_parameters.update({'phase': 1})
strain_phase1 = waveform_generator.frequency_domain_strain(new_injection_parameters)
detector_strain_phase1 = interferometer.get_detector_response(strain_phase1, parameters = new_injection_parameters)

rotated_waveform = detector_strain_phase0 * np.exp(2j*1) #m=2, phase=1
np.sum(np.abs(rotated_waveform-detector_strain_phase1))

This outputs 2.094249845908886e-32, so I think IMRPhenomPv2 obeys the expected transformation. (Note the injected configuration precesses)

bilby-bot commented 3 months ago

In GitLab by @git.ligo:colm.talbot on Jun 27, 2024, 00:08

It isn't very highly spinning, can you run this again with large mostly in plane spins.

bilby-bot commented 3 months ago

In GitLab by @git.ligo:colm.talbot on Jun 27, 2024, 00:08

We could also ask a waveforms person @git.ligo:patricia-schmidt or @git.ligo:geraint.pratten do you have an authoritative take?

bilby-bot commented 3 months ago

In GitLab by @git.ligo:jacob.golomb on Jun 27, 2024, 04:52

Changing the spin parameters to:

    a_1=0.9,
    a_2=0.9,
    cos_tilt_1=0.5, 
    cos_tilt_2=-0.8,

which corresponds to $\chi_p \approx 0.78$ in the injected configuration also leads to an output (sum of absolute strain differences) $\mathcal{O}(10^{-33}$. I tried with several similar high-precessing parameter combinations and it seems to consistently obey the transformation (very small difference between the explicitly generated and manually rotated waveform).

bilby-bot commented 3 months ago

In GitLab by @git.ligo:sascha.husa on Jun 27, 2024, 06:28

PhenomPv2 and all other precessing models do mix the modes in an inertial frame, e.g. the one of the observer. PhenomP/X/T essentially use the orbital plane rotation to rotate the waveform from a non-inertial frame to an inertial one. This procedure mixes the modes.

bilby-bot commented 3 months ago

In GitLab by @git.ligo:patricia-schmidt on Jun 27, 2024, 09:22

Hi @git.ligo:jacob.golomb @git.ligo:colm.talbot. As Pv2 only has the l = m = |2| modes in the co-precessing frame (from PhenomD in this case), it can be shown that the coalescence phase from PhenomD can be factored out in both the (inertial) plus and cross polarisation and hence both polarisations transform as

$h{+,\times} \rightarrow e^{-2i\delta \phi{c,D}} h_{+,\times}$.

This is the same relationship as for the aligned-spin case and applies to all precessing models that use te "twisting-up" methodology and only the l = m = |2| co-precessing frame modes (e.g. also XP). However, this relationship no longer holds when higher-order co-precessing frame modes are included to construct the inertial precessing waveform (e.g. as is the case for XPHM).

For the detailed maths, please refer to the Pv2 technical document, especially Sec. 3 and Eqn. (17) and (19).

Thanks!

bilby-bot commented 3 months ago

In GitLab by @git.ligo:jacob.golomb on Jun 28, 2024, 00:27

Thanks @git.ligo:patricia-schmidt. So I think based on the examples above and eq. 9 in the document Patricia linked above, the statement that phase marginalization is invalid for IMRPhenomPv2 can be removed?