# Copyright (C) 2022 Shichao Wu, Alex Nitz
# This program is free software; you can redistribute it and/or modify it
# under the terms of the GNU General Public License as published by the
# Free Software Foundation; either version 3 of the License, or (at your
# option) any later version.
#
# This program is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
# Public License for more details.
#
# You should have received a copy of the GNU General Public License along
# with this program; if not, write to the Free Software Foundation, Inc.,
# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
#
# =============================================================================
#
# Preamble
#
# =============================================================================
#
"""
This module provides (semi-)analytical PSDs and sensitivity curves for space
borne detectors, such as LISA. Based on LISA technical note
<LISA-LCST-SGS-TN-001>, LDC manual <LISA-LCST-SGS-MAN-001>,
and paper <10.1088/1361-6382/ab1101>.
"""
import numpy as np
from astropy.constants import c
from pycbc.psd.read import from_numpy_arrays
[docs]def psd_lisa_acc_noise(f, acc_noise_level=3e-15):
""" The PSD of LISA's acceleration noise.
Parameters
----------
f : float or numpy.array
The frequency or frequency range, in the unit of "Hz".
acc_noise_level : float
The level of acceleration noise.
Returns
-------
s_acc_nu : float or numpy.array
The PSD value or array for acceleration noise.
Notes
-----
Please see Eq.(11-13) in <LISA-LCST-SGS-TN-001> for more details.
"""
s_acc = acc_noise_level**2 * (1+(4e-4/f)**2)*(1+(f/8e-3)**4)
s_acc_d = s_acc * (2*np.pi*f)**(-4)
s_acc_nu = (2*np.pi*f/c.value)**2 * s_acc_d
return s_acc_nu
[docs]def psd_lisa_oms_noise(f, oms_noise_level=15e-12):
""" The PSD of LISA's OMS noise.
Parameters
----------
f : float or numpy.array
The frequency or frequency range, in the unit of "Hz".
oms_noise_level : float
The level of OMS noise.
Returns
-------
s_oms_nu : float or numpy.array
The PSD value or array for OMS noise.
Notes
-----
Please see Eq.(9-10) in <LISA-LCST-SGS-TN-001> for more details.
"""
s_oms_d = oms_noise_level**2 * (1+(2e-3/f)**4)
s_oms_nu = s_oms_d * (2*np.pi*f/c.value)**2
return s_oms_nu
[docs]def lisa_psd_components(f, acc_noise_level=3e-15, oms_noise_level=15e-12):
""" The PSD of LISA's acceleration and OMS noise.
Parameters
----------
f : float or numpy.array
The frequency or frequency range, in the unit of "Hz".
acc_noise_level : float
The level of acceleration noise.
oms_noise_level : float
The level of OMS noise.
Returns
-------
low_freq_component, high_freq_component :
The PSD value or array for acceleration and OMS noise.
"""
low_freq_component = psd_lisa_acc_noise(f, acc_noise_level)
high_freq_component = psd_lisa_oms_noise(f, oms_noise_level)
return low_freq_component, high_freq_component
[docs]def omega_length(f, len_arm=2.5e9):
""" The function to calculate 2*pi*f*LISA_arm_length.
Parameters
----------
f : float or numpy.array
The frequency or frequency range, in the unit of "Hz".
len_arm : float
The arm length of LISA.
Returns
-------
omega_len : float or numpy.array
The value of 2*pi*f*LISA_arm_length.
"""
omega_len = 2*np.pi*f * len_arm/c.value
return omega_len
[docs]def analytical_psd_lisa_tdi_1p5_XYZ(length, delta_f, low_freq_cutoff,
len_arm=2.5e9, acc_noise_level=3e-15,
oms_noise_level=15e-12):
""" The TDI-1.5 analytical PSD (X,Y,Z channel) for LISA.
Parameters
----------
length : int
Length of output Frequencyseries.
delta_f : float
Frequency step for output FrequencySeries.
low_freq_cutoff : float
Low-frequency cutoff for output FrequencySeries.
len_arm : float
The arm length of LISA, in the unit of "m".
acc_noise_level : float
The level of acceleration noise.
oms_noise_level : float
The level of OMS noise.
Returns
-------
fseries : FrequencySeries
The TDI-1.5 PSD (X,Y,Z channel) for LISA.
Notes
-----
Please see Eq.(19) in <LISA-LCST-SGS-TN-001> for more details.
"""
len_arm = np.float64(len_arm)
acc_noise_level = np.float64(acc_noise_level)
oms_noise_level = np.float64(oms_noise_level)
fr = np.linspace(low_freq_cutoff, (length-1)*2*delta_f, length)
s_acc_nu, s_oms_nu = lisa_psd_components(
fr, acc_noise_level, oms_noise_level)
omega_len = omega_length(fr, len_arm)
psd = 16*(np.sin(omega_len))**2 * (s_oms_nu+s_acc_nu*(3+np.cos(omega_len)))
return from_numpy_arrays(fr, psd, length, delta_f, low_freq_cutoff)
[docs]def analytical_psd_lisa_tdi_2p0_XYZ(length, delta_f, low_freq_cutoff,
len_arm=2.5e9, acc_noise_level=3e-15,
oms_noise_level=15e-12):
""" The TDI-2.0 analytical PSD (X,Y,Z channel) for LISA.
Parameters
----------
length : int
Length of output Frequencyseries.
delta_f : float
Frequency step for output FrequencySeries.
low_freq_cutoff : float
Low-frequency cutoff for output FrequencySeries.
len_arm : float
The arm length of LISA, in the unit of "m".
acc_noise_level : float
The level of acceleration noise.
oms_noise_level : float
The level of OMS noise.
Returns
-------
fseries : FrequencySeries
The TDI-2.0 PSD (X,Y,Z channel) for LISA.
Notes
-----
Please see Eq.(20) in <LISA-LCST-SGS-TN-001> for more details.
"""
len_arm = np.float64(len_arm)
acc_noise_level = np.float64(acc_noise_level)
oms_noise_level = np.float64(oms_noise_level)
fr = np.linspace(low_freq_cutoff, (length-1)*2*delta_f, length)
s_acc_nu, s_oms_nu = lisa_psd_components(
fr, acc_noise_level, oms_noise_level)
omega_len = omega_length(fr, len_arm)
psd = (64*(np.sin(omega_len))**2 * (np.sin(2*omega_len))**2 *
(s_oms_nu+s_acc_nu*(3+np.cos(2*omega_len))))
return from_numpy_arrays(fr, psd, length, delta_f, low_freq_cutoff)
[docs]def analytical_csd_lisa_tdi_1p5_XY(length, delta_f, low_freq_cutoff,
len_arm=2.5e9, acc_noise_level=3e-15,
oms_noise_level=15e-12):
""" The cross-spectrum density between LISA's TDI channel X and Y.
Parameters
----------
length : int
Length of output Frequencyseries.
delta_f : float
Frequency step for output FrequencySeries.
low_freq_cutoff : float
Low-frequency cutoff for output FrequencySeries.
len_arm : float
The arm length of LISA, in the unit of "m".
acc_noise_level : float
The level of acceleration noise.
oms_noise_level : float
The level of OMS noise.
Returns
-------
fseries : FrequencySeries
The CSD between LISA's TDI-1.5 channel X and Y.
Notes
-----
Please see Eq.(56) in <LISA-LCST-SGS-MAN-001(Radler)> for more details.
"""
len_arm = np.float64(len_arm)
acc_noise_level = np.float64(acc_noise_level)
oms_noise_level = np.float64(oms_noise_level)
fr = np.linspace(low_freq_cutoff, (length-1)*2*delta_f, length)
s_acc_nu, s_oms_nu = lisa_psd_components(
fr, acc_noise_level, oms_noise_level)
omega_len = omega_length(fr, len_arm)
csd = (-8*np.sin(omega_len)**2 * np.cos(omega_len) *
(s_oms_nu+4*s_acc_nu))
return from_numpy_arrays(fr, csd, length, delta_f, low_freq_cutoff)
[docs]def analytical_psd_lisa_tdi_1p5_AE(length, delta_f, low_freq_cutoff,
len_arm=2.5e9, acc_noise_level=3e-15,
oms_noise_level=15e-12):
""" The PSD of LISA's TDI-1.5 channel A and E.
Parameters
----------
length : int
Length of output Frequencyseries.
delta_f : float
Frequency step for output FrequencySeries.
low_freq_cutoff : float
Low-frequency cutoff for output FrequencySeries.
len_arm : float
The arm length of LISA, in the unit of "m".
acc_noise_level : float
The level of acceleration noise.
oms_noise_level : float
The level of OMS noise.
Returns
-------
fseries : FrequencySeries
The PSD of LISA's TDI-1.5 channel A and E.
Notes
-----
Please see Eq.(58) in <LISA-LCST-SGS-MAN-001(Radler)> for more details.
"""
len_arm = np.float64(len_arm)
acc_noise_level = np.float64(acc_noise_level)
oms_noise_level = np.float64(oms_noise_level)
fr = np.linspace(low_freq_cutoff, (length-1)*2*delta_f, length)
s_acc_nu, s_oms_nu = lisa_psd_components(
fr, acc_noise_level, oms_noise_level)
omega_len = omega_length(fr, len_arm)
psd = (8*(np.sin(omega_len))**2 *
(4*(1+np.cos(omega_len)+np.cos(omega_len)**2)*s_acc_nu +
(2+np.cos(omega_len))*s_oms_nu))
return from_numpy_arrays(fr, psd, length, delta_f, low_freq_cutoff)
[docs]def analytical_psd_lisa_tdi_1p5_T(length, delta_f, low_freq_cutoff,
len_arm=2.5e9, acc_noise_level=3e-15,
oms_noise_level=15e-12):
""" The PSD of LISA's TDI-1.5 channel T.
Parameters
----------
length : int
Length of output Frequencyseries.
delta_f : float
Frequency step for output FrequencySeries.
low_freq_cutoff : float
Low-frequency cutoff for output FrequencySeries.
len_arm : float
The arm length of LISA, in the unit of "m".
acc_noise_level : float
The level of acceleration noise.
oms_noise_level : float
The level of OMS noise.
Returns
-------
fseries : FrequencySeries
The PSD of LISA's TDI-1.5 channel T.
Notes
-----
Please see Eq.(59) in <LISA-LCST-SGS-MAN-001(Radler)> for more details.
"""
len_arm = np.float64(len_arm)
acc_noise_level = np.float64(acc_noise_level)
oms_noise_level = np.float64(oms_noise_level)
fr = np.linspace(low_freq_cutoff, (length-1)*2*delta_f, length)
s_acc_nu, s_oms_nu = lisa_psd_components(
fr, acc_noise_level, oms_noise_level)
omega_len = omega_length(fr, len_arm)
psd = (32*np.sin(omega_len)**2 * np.sin(omega_len/2)**2 *
(4*s_acc_nu*np.sin(omega_len/2)**2 + s_oms_nu))
return from_numpy_arrays(fr, psd, length, delta_f, low_freq_cutoff)
[docs]def averaged_lisa_fplus_sq_approx(f, len_arm=2.5e9):
""" An approximant for LISA's squared antenna response function,
averaged over sky and polarization angle.
Parameters
----------
f : float or numpy.array
The frequency or frequency range, in the unit of "Hz".
len_arm : float
The arm length of LISA, in the unit of "m".
Returns
-------
fp_sq_approx : float or numpy.array
The sky and polarization angle averaged squared antenna response.
Notes
-----
Please see Eq.(36) in <LISA-LCST-SGS-TN-001> for more details.
"""
from scipy.interpolate import interp1d
from astropy.utils.data import download_file
if len_arm != 2.5e9:
raise Exception("Currently only support 'len_arm=2.5e9'.")
# Download the numerical LISA averaged response.
url = "https://zenodo.org/record/7497853/files/AvFXp2_Raw.npy"
file_path = download_file(url, cache=True)
freqs, fp_sq = np.load(file_path)
# Padding the end.
freqs = np.append(freqs, 2)
fp_sq = np.append(fp_sq, 0.0012712348970728724)
fp_sq_interp = interp1d(freqs, fp_sq, kind='linear',
fill_value="extrapolate")
fp_sq_approx = fp_sq_interp(f)/16
return fp_sq_approx
[docs]def averaged_response_lisa_tdi_1p5(f, len_arm=2.5e9):
""" LISA's TDI-1.5 response function to GW,
averaged over sky and polarization angle.
Parameters
----------
f : float or numpy.array
The frequency or frequency range, in the unit of "Hz".
len_arm : float
The arm length of LISA, in the unit of "m".
Returns
-------
response_tdi_1p5 : float or numpy.array
The sky and polarization angle averaged TDI-1.5 response to GW.
Notes
-----
Please see Eq.(39) in <LISA-LCST-SGS-TN-001> for more details.
"""
omega_len = omega_length(f, len_arm)
ave_fp2 = averaged_lisa_fplus_sq_approx(f, len_arm)
response_tdi_1p5 = (4*omega_len)**2 * np.sin(omega_len)**2 * ave_fp2
return response_tdi_1p5
[docs]def averaged_response_lisa_tdi_2p0(f, len_arm=2.5e9):
""" LISA's TDI-2.0 response function to GW,
averaged over sky and polarization angle.
Parameters
----------
f : float or numpy.array
The frequency or frequency range, in the unit of "Hz".
len_arm : float
The arm length of LISA, in the unit of "m".
Returns
-------
response_tdi_2p0 : float or numpy.array
The sky and polarization angle averaged TDI-2.0 response to GW.
Notes
-----
Please see Eq.(40) in <LISA-LCST-SGS-TN-001> for more details.
"""
omega_len = omega_length(f, len_arm)
response_tdi_1p5 = averaged_response_lisa_tdi_1p5(f, len_arm)
response_tdi_2p0 = response_tdi_1p5 * (2*np.sin(2*omega_len))**2
return response_tdi_2p0
[docs]def sensitivity_curve_lisa_semi_analytical(length, delta_f, low_freq_cutoff,
len_arm=2.5e9,
acc_noise_level=3e-15,
oms_noise_level=15e-12):
""" The semi-analytical LISA's sensitivity curve (6-links),
averaged over sky and polarization angle.
Parameters
----------
length : int
Length of output Frequencyseries.
delta_f : float
Frequency step for output FrequencySeries.
low_freq_cutoff : float
Low-frequency cutoff for output FrequencySeries.
len_arm : float
The arm length of LISA, in the unit of "m".
acc_noise_level : float
The level of acceleration noise.
oms_noise_level : float
The level of OMS noise.
Returns
-------
fseries : FrequencySeries
The sky and polarization angle averaged semi-analytical
LISA's sensitivity curve (6-links).
Notes
-----
Please see Eq.(42-43) in <LISA-LCST-SGS-TN-001> for more details.
"""
len_arm = np.float64(len_arm)
acc_noise_level = np.float64(acc_noise_level)
oms_noise_level = np.float64(oms_noise_level)
fr = np.linspace(low_freq_cutoff, (length-1)*2*delta_f, length)
fp_sq = averaged_lisa_fplus_sq_approx(fr, len_arm)
s_acc_nu, s_oms_nu = lisa_psd_components(
fr, acc_noise_level, oms_noise_level)
omega_len = 2*np.pi*fr * len_arm/c.value
sense_curve = ((s_oms_nu + s_acc_nu*(3+np.cos(2*omega_len))) /
(omega_len**2*fp_sq))
fseries = from_numpy_arrays(fr, sense_curve/2,
length, delta_f, low_freq_cutoff)
return fseries
[docs]def sensitivity_curve_lisa_SciRD(length, delta_f, low_freq_cutoff):
""" The analytical LISA's sensitivity curve in SciRD,
averaged over sky and polarization angle.
Parameters
----------
length : int
Length of output Frequencyseries.
delta_f : float
Frequency step for output FrequencySeries.
low_freq_cutoff : float
Low-frequency cutoff for output FrequencySeries.
Returns
-------
fseries : FrequencySeries
The sky and polarization angle averaged analytical
LISA's sensitivity curve in SciRD.
Notes
-----
Please see Eq.(114) in <LISA-LCST-SGS-TN-001> for more details.
"""
fr = np.linspace(low_freq_cutoff, (length-1)*2*delta_f, length)
s_I = 5.76e-48 * (1+(4e-4/fr)**2)
s_II = 3.6e-41
R = 1 + (fr/2.5e-2)**2
sense_curve = 10/3 * (s_I/(2*np.pi*fr)**4+s_II) * R
return from_numpy_arrays(fr, sense_curve, length, delta_f, low_freq_cutoff)
[docs]def sensitivity_curve_lisa_confusion(length, delta_f, low_freq_cutoff,
len_arm=2.5e9, acc_noise_level=3e-15,
oms_noise_level=15e-12,
base_model="semi", duration=1.0):
""" The LISA's sensitivity curve with Galactic confusion noise,
averaged over sky and polarization angle.
Parameters
----------
length : int
Length of output Frequencyseries.
delta_f : float
Frequency step for output FrequencySeries.
low_freq_cutoff : float
Low-frequency cutoff for output FrequencySeries.
len_arm : float
The arm length of LISA, in the unit of "m".
acc_noise_level : float
The level of acceleration noise.
oms_noise_level : float
The level of OMS noise.
base_model : string
The base model of sensitivity curve, chosen from "semi" or "SciRD".
duration : float
The duration of observation, between 0 and 10, in the unit of years.
Returns
-------
fseries : FrequencySeries
The sky and polarization angle averaged
LISA's sensitivity curve with Galactic confusion noise.
Notes
-----
Please see Eq.(85-86) in <LISA-LCST-SGS-TN-001> for more details.
"""
if base_model == "semi":
base_curve = sensitivity_curve_lisa_semi_analytical(
length, delta_f, low_freq_cutoff,
len_arm, acc_noise_level, oms_noise_level)
elif base_model == "SciRD":
base_curve = sensitivity_curve_lisa_SciRD(
length, delta_f, low_freq_cutoff)
else:
raise Exception("Must choose from 'semi' or 'SciRD'.")
if duration < 0 or duration > 10:
raise Exception("Must between 0 and 10.")
fr = np.linspace(low_freq_cutoff, (length-1)*2*delta_f, length)
f1 = 10**(-0.25*np.log10(duration)-2.7)
fk = 10**(-0.27*np.log10(duration)-2.47)
sh_confusion = (0.5*1.14e-44*fr**(-7/3)*np.exp(-(fr/f1)**1.8) *
(1.0+np.tanh((fk-fr)/0.31e-3)))
fseries_confusion = from_numpy_arrays(fr, sh_confusion,
length, delta_f, low_freq_cutoff)
fseries = from_numpy_arrays(base_curve.sample_frequencies,
base_curve+fseries_confusion,
length, delta_f, low_freq_cutoff)
return fseries
[docs]def semi_analytical_psd_lisa_confusion_noise(length, delta_f, low_freq_cutoff,
len_arm=2.5e9, duration=1.0,
tdi="1.5"):
""" The TDI-1.5/2.0 PSD (X,Y,Z channel) for LISA Galactic confusion noise,
no instrumental noise.
Parameters
----------
length : int
Length of output Frequencyseries.
delta_f : float
Frequency step for output FrequencySeries.
low_freq_cutoff : float
Low-frequency cutoff for output FrequencySeries.
len_arm : float
The arm length of LISA, in the unit of "m".
duration : float
The duration of observation, between 0 and 10, in the unit of years.
tdi : string
The version of TDI, currently only for 1.5 or 2.0.
Returns
-------
fseries : FrequencySeries
The TDI-1.5/2.0 PSD (X,Y,Z channel) for LISA Galactic confusion
noise, no instrumental noise.
"""
fr = np.linspace(low_freq_cutoff, (length-1)*2*delta_f, length)
if tdi == "1.5":
response = averaged_response_lisa_tdi_1p5(fr, len_arm)
elif tdi == "2.0":
response = averaged_response_lisa_tdi_2p0(fr, len_arm)
else:
raise Exception("The version of TDI, currently only for 1.5 or 2.0.")
fseries_response = from_numpy_arrays(fr, np.array(response),
length, delta_f, low_freq_cutoff)
f1 = 10**(-0.25*np.log10(duration)-2.7)
fk = 10**(-0.27*np.log10(duration)-2.47)
sh_confusion = (0.5*1.14e-44*fr**(-7/3)*np.exp(-(fr/f1)**1.8) *
(1.0+np.tanh((fk-fr)/0.31e-3)))
fseries_confusion = from_numpy_arrays(fr, sh_confusion,
length, delta_f, low_freq_cutoff)
psd_confusion = 2*fseries_confusion.data * fseries_response.data
fseries = from_numpy_arrays(fseries_confusion.sample_frequencies,
psd_confusion, length, delta_f,
low_freq_cutoff)
return fseries
[docs]def analytical_psd_lisa_tdi_AE_confusion(length, delta_f, low_freq_cutoff,
len_arm=2.5e9, acc_noise_level=3e-15,
oms_noise_level=15e-12,
duration=1.0, tdi="1.5"):
""" The TDI-1.5 PSD (A,E channel) for LISA
with Galactic confusion noise.
Parameters
----------
length : int
Length of output Frequencyseries.
delta_f : float
Frequency step for output FrequencySeries.
low_freq_cutoff : float
Low-frequency cutoff for output FrequencySeries.
len_arm : float
The arm length of LISA, in the unit of "m".
acc_noise_level : float
The level of acceleration noise.
oms_noise_level : float
The level of OMS noise.
duration : float
The duration of observation, between 0 and 10, in the unit of years.
tdi : string
The version of TDI, currently only for 1.5.
Returns
-------
fseries : FrequencySeries
The TDI-1.5 PSD (A,E channel) for LISA with Galactic confusion
noise.
"""
if tdi != "1.5":
raise Exception("The version of TDI, currently only for 1.5.")
psd_AE = analytical_psd_lisa_tdi_1p5_AE(length, delta_f, low_freq_cutoff,
len_arm, acc_noise_level,
oms_noise_level)
psd_X_confusion = semi_analytical_psd_lisa_confusion_noise(
length, delta_f, low_freq_cutoff,
len_arm, duration, tdi)
# Here we assume the confusion noise's contribution to the CSD Sxy is
# negligible for low-frequency part. So Sxy doesn't change.
fseries = psd_AE + psd_X_confusion
return fseries