Source code for tgr.massive_graviton

import numpy as np
import scipy.constants
from scipy.integrate import quad
import pycbc.cosmology



[docs]
def mg_to_lambda_g(mg):
    '''
    Convert graviton mass to Compton wavelength.

    Parameters
    ----------
    mg : float
        Graviton mass in eV/c^2

    Returns
    -------
    lambda_g : float
        Compton wavelength of graviton in km
    '''
    h_eV_s, _, _ = scipy.constants.physical_constants['Planck constant in eV/Hz']
    return h_eV_s / (mg / scipy.constants.c) / 1000 # convert from m to km





[docs]
def lambda_g_to_mg(lambda_g):
    '''
    Convert graviton Compton wavelength to mass.
    
    Parameters
    ----------
    lambda_g : float
        Compton wavelength of graviton in km
    
    Returns
    -------
    mg : float
        Graviton mass in eV/c^2
    '''
    h_eV_s, _, _ = scipy.constants.physical_constants['Planck constant in eV/Hz']
    lambda_g_m = lambda_g * 1000 # convert from km to m
    return h_eV_s / (lambda_g_m / scipy.constants.c)



def _get_cosmology(cosmology=None, **kwargs):
    return pycbc.cosmology.get_cosmology(cosmology=cosmology, **kwargs)


def _effective_distance_from_redshift(z, cosmo):
    def integrand(zp):
        return 1.0 / ((1.0 + zp)**2 * cosmo.efunc(zp))

    integral, _ = quad(integrand, 0, z)
    c_km_s = scipy.constants.c / 1000.0
    return c_km_s * (1.0 + z) / cosmo.H0.value * integral




[docs]
def effective_distance(luminosity_distance, cosmology=None, **kwargs):
    '''
    Calculate effective distance D_eff for massive graviton phase correction.
    D_eff = (1+z) * c / H0 * integral_0^z dz' / ((1+z')^2 * E(z'))

    Parameters
    ----------
    luminosity_distance : float
        Luminosity distance to the source in Mpc
    
    Returns
    -------
    deff : float
        Effective distance in Mpc
    '''
    cosmo = _get_cosmology(cosmology=cosmology, **kwargs)
    z = pycbc.cosmology.redshift(luminosity_distance, cosmology=cosmo)
    return _effective_distance_from_redshift(z, cosmo)





[docs]
def mg_phase_correction(mg, distance, frequencies, cosmology=None, **kwargs):
    '''
    Calculate the massive graviton phase correction for given frequencies.
    delta_phi = -pi * D_eff * c^3 * m_g^2 / h^2 / (1 + z) / f

    Parameters
    ----------
    mg : float
        Graviton mass in eV/c^2
    distance : float
        Luminosity distance to the source in Mpc
    frequencies : array-like
        Array of frequencies in Hz

    Returns
    -------
    delta_phi : array-like
        Phase correction for each frequency in radians
    '''
    cosmo = _get_cosmology(cosmology=cosmology, **kwargs)
    z = pycbc.cosmology.redshift(distance, cosmology=cosmo)
    deff = _effective_distance_from_redshift(z, cosmo) * 1e6 * scipy.constants.parsec # in unit of m
    c = scipy.constants.speed_of_light
    h_eV_s, _, _ = scipy.constants.physical_constants['Planck constant in eV/Hz']
    delta_phi = - np.pi * deff * mg**2 / c / h_eV_s**2 / (1 + z) / frequencies
    return delta_phi





[docs]
def gen_mg_waveform(**kwds):
    '''
    Generate frequency-domain waveform with massive graviton phase correction.
    Based on Will (1997) PRD 57, 2061. Assumes no modification to binary dynamics.

    Parameters
    ----------
    kwds : dict
        Must contain ``base_gr_approximant`` for the GR approximant name,
        ``mg`` for the graviton mass in eV/c^2, and ``distance`` for the
        luminosity distance in Mpc.

        Other parameters (masses, spins, etc.) passed to get_fd_waveform.

    Returns
    -------
    hp: pycbc.types.FrequencySeries
        Plus polarization time series
    hc: pycbc.types.FrequencySeries
        Cross polarization time series
    '''
    from pycbc.waveform import get_fd_waveform

    # sanity checks
    for key in ['base_gr_approximant', 'mg', 'distance']:
        if key not in kwds or kwds[key] is None:
            raise ValueError(f"Missing required argument: {key}")

    # Generate GR waveforms
    if 'approximant' in kwds:
        kwds.pop("approximant")
    hp, hc = get_fd_waveform(approximant=kwds['base_gr_approximant'], **kwds)

    # Apply massive graviton phase correction
    delta_phi = mg_phase_correction(kwds['mg'], kwds['distance'], hp.sample_frequencies[1:])

    # slicing with index 1 to avoid dividing zero frequency
    hp[1:] *= np.exp(1j * delta_phi)
    hc[1:] *= np.exp(1j * delta_phi)
    
    return hp, hc