Fix conversion from z to f and add noise figures

CFT-HY · Nov 25, 2024 · e3354e4 · e3354e4
1 parent 41cc55f
commit e3354e4
Show file tree

Hide file tree

Showing 13 changed files with 228 additions and 32 deletions.
diff --git a/docs/conf.py b/docs/conf.py
@@ -144,6 +144,7 @@
     # Other articles
     "ai_2023": ("https://arxiv.org/abs/2303.10171%s", "Ai et al., 2023%s"),
     "borsanyi_2016": ("https://arxiv.org/abs/1606.07494%s", "Borsanyi et al., 2016%s"),
+    "smith_2019": ("https://arxiv.org/abs/1908.00546%s", "Smith & Caldwell, 2019%s"),
     "caprini_2020": ("https://arxiv.org/abs/1910.13125%s", "Caprini et al., 2020%s"),
     "giese_2020": ("https://arxiv.org/abs/2004.06995%s", "Giese et al., 2020%s"),
     "giese_2021": ("https://arxiv.org/abs/2010.09744%s", "Giese et al., 2021%s"),

diff --git a/examples/const_cs_gw.py b/examples/const_cs_gw.py
@@ -54,6 +54,9 @@ def main():
             spectrum_kwargs={
                 "r_star": r_star
                 # "z": z
+                # "Tn": 100,
+                # "g_star": 100,
+                # "gs_star": 100
             }
             # bubble_kwargs={"allow_invalid": False}, allow_bubble_failure=True
         )
@@ -150,9 +153,9 @@ def main():
     fig.tight_layout()
     fig2.tight_layout()
     fig3.tight_layout()
-    utils.save_and_show(fig, "const_cs_gw.png")
-    utils.save_and_show(fig2, "const_cs_gw_v.png")
-    utils.save_and_show(fig3, "const_cs_gw_omgw0.png")
+    utils.save_and_show(fig, "const_cs_gw")
+    utils.save_and_show(fig2, "const_cs_gw_v")
+    utils.save_and_show(fig3, "const_cs_gw_omgw0")
 
 
 if __name__ == "__main__":

diff --git a/examples/noise.py b/examples/noise.py
@@ -0,0 +1,44 @@
+from matplotlib import pyplot as plt
+import numpy as np
+
+from examples.utils import save_and_show
+from pttools.omgw0 import noise
+
+
+def main():
+    fig: plt.Figure = plt.figure()
+    axs = fig.subplots(2, 2)
+
+    f = np.logspace(-4, -1, 50)
+    ax1 = axs[0, 0]
+    P_oms = noise.P_oms()
+    ax1.plot([f.min(), f.max()], [P_oms, P_oms], label=r"$P_\text{oms}$")
+    ax1.plot(f, noise.P_acc(f), label=r"$P_\text{acc}$")
+    ax1.set_ylabel("$P(f)$")
+
+    ax2 = axs[0, 1]
+    ax2.plot(f, noise.S_AE(f), label="$S_{A,E}$")
+    ax2.plot(f, noise.S_AE_approx(f), label=r"$S_{A,E,\text{approx}}$")
+    ax2.plot(f, noise.S_gb(f), label=r"$S_\text{gb}$")
+    ax2.set_ylabel("$S(f)$")
+
+    ax3 = axs[1, 0]
+    ax3.plot(f, noise.omega_ins(f), label=r"$\Omega_\text{ins}$")
+    ax3.plot(f, noise.omega_eb(f), label=r"$\Omega_\text{eb}$")
+    ax3.plot(f, noise.omega_gb(f), label=r"$\Omega_\text{gb}$")
+    ax3.plot(f, noise.omega_noise(f), label=r"$\Omega_\text{noise}$")
+    ax3.set_ylabel(r"$\Omega(f)$")
+
+    for ax in axs.flat:
+        ax.set_xlabel(r"$f(\text{Hz})$")
+        ax.set_xscale("log")
+        ax.set_yscale("log")
+        ax.legend()
+    fig.tight_layout()
+
+    return fig
+
+
+if __name__ == '__main__':
+    fig = main()
+    save_and_show(fig, "noise")
diff --git a/pttools/models/bag.py b/pttools/models/bag.py
@@ -35,6 +35,7 @@ class BagModel(AnalyticModel):
     DEFAULT_LABEL_LATEX = "Bag model"
     DEFAULT_LABEL_UNICODE = DEFAULT_LABEL_LATEX
     DEFAULT_NAME = "bag"
+    TEMPERATURE_IS_PHYSICAL = False
 
     def __init__(
             self,

diff --git a/pttools/models/base.py b/pttools/models/base.py
@@ -13,14 +13,20 @@
 
 
 class BaseModel(abc.ABC):
-    """The base for both Model and ThermoModel"""
+    """The base for both Model and ThermoModel
+
+    All temperatures must be in units of GeV for the frequency conversion in Spectrum to work.
+    """
     DEFAULT_LABEL_LATEX: str = None
     DEFAULT_LABEL_UNICODE: str = None
     DEFAULT_NAME: str = None
     # Zero temperature would break many of the equations
     DEFAULT_T_MIN: float = 1e-3
     DEFAULT_T_MAX: float = np.inf
 
+    #: Whether the temperature is in proper physics units
+    TEMPERATURE_IS_PHYSICAL: bool = None
+
     def __init__(
             self,
             name: str = None,
@@ -29,13 +35,15 @@ def __init__(
             label_latex: str = None,
             label_unicode: str = None,
             gen_cs2: bool = True,
-            gen_cs2_neg: bool = True):
+            gen_cs2_neg: bool = True,
+            temperature_is_physical: bool = None):
         self.name = self.DEFAULT_NAME if name is None else name
         self.label_latex = self.DEFAULT_LABEL_LATEX if label_latex is None else label_latex
         self.label_unicode = self.DEFAULT_LABEL_UNICODE if label_unicode is None else label_unicode
         self.T_min = self.DEFAULT_T_MIN if T_min is None else T_min
         self.T_max = self.DEFAULT_T_MAX if T_max is None else T_max
         self.restrict_to_valid = restrict_to_valid
+        self.temperature_is_physical = self.TEMPERATURE_IS_PHYSICAL if temperature_is_physical is None else temperature_is_physical
 
         if self.name is None:
             raise ValueError("The model must have a name.")

diff --git a/pttools/models/const_cs.py b/pttools/models/const_cs.py
@@ -29,6 +29,7 @@ class ConstCSModel(AnalyticModel):
     DEFAULT_LABEL_LATEX = "Constant $c_s$ model"
     DEFAULT_LABEL_UNICODE = "Constant cₛ model"
     DEFAULT_NAME = "const_cs"
+    TEMPERATURE_IS_PHYSICAL = False
 
     def __init__(
             self,

diff --git a/pttools/models/const_cs_thermo.py b/pttools/models/const_cs_thermo.py
@@ -13,6 +13,7 @@ class ConstCSThermoModel(ThermoModel):
     DEFAULT_LABEL_LATEX = "Constant $c_s$ thermo-model"
     DEFAULT_LABEL_UNICODE = "Constant cₛ thermo-model"
     DEFAULT_NAME = "const_cs_thermo"
+    TEMPERATURE_IS_PHYSICAL = False
 
     GEFF_DATA_LOG_TEMP = np.linspace(-1, 3, 1000)
     GEFF_DATA_TEMP = 10**GEFF_DATA_LOG_TEMP

diff --git a/pttools/models/full.py b/pttools/models/full.py
@@ -49,7 +49,9 @@ def __init__(
             V_s=V_s, V_b=V_b,
             T_min=thermo.t_min, T_max=thermo.t_max,
             name=name, label_latex=label_latex, label_unicode=label_unicode,
-            gen_critical=False, gen_cs2=False, gen_cs2_neg=False, implicit_V=True)
+            gen_critical=False, gen_cs2=False, gen_cs2_neg=False, implicit_V=True,
+            temperature_is_physical=thermo.TEMPERATURE_IS_PHYSICAL
+        )
 
         self.temp_spline_s = splrep(
             np.log10(self.w(self.thermo.GEFF_DATA_TEMP, Phase.SYMMETRIC)), self.thermo.GEFF_DATA_LOG_TEMP

diff --git a/pttools/models/model.py b/pttools/models/model.py
@@ -43,6 +43,7 @@ def __init__(
             gen_cs2: bool = True,
             gen_cs2_neg: bool = True,
             implicit_V: bool = False,
+            temperature_is_physical: bool = None,
             allow_invalid: bool = False):
 
         if implicit_V:
@@ -62,6 +63,13 @@ def __init__(
             # if V_s == V_b:
             #     logger.warning("The bubble will not expand, when V_s <= V_b. Got: V_b = V_s = %s.", V_s)
 
+        self.temperature_is_physical = self.TEMPERATURE_IS_PHYSICAL if temperature_is_physical is None else temperature_is_physical
+        if self.temperature_is_physical is None:
+            raise ValueError(
+                "It has not been specified whether the temperature scale for the model is physical. "
+                "Please specify it in the model definition."
+            )
+
         self.T_ref: float = T_ref
         self.V_s: float = V_s
         self.V_b: float = V_b

diff --git a/pttools/models/sm.py b/pttools/models/sm.py
@@ -21,6 +21,8 @@ class StandardModel(ThermoModel):
     DEFAULT_LABEL_LATEX = "Standard Model"
     DEFAULT_LABEL_UNICODE = DEFAULT_LABEL_LATEX
     DEFAULT_NAME = "standard_model"
+    TEMPERATURE_IS_PHYSICAL = True
+
     # Copied from the ArXiv file som_eos.tex
     GEFF_DATA = np.array([
         [0.00, 10.71, 1.00228],

diff --git a/pttools/omgw0/noise.py b/pttools/omgw0/noise.py
@@ -13,7 +13,8 @@ def signal_to_noise_ratio(
     $$\rho = \sqrt{T_{\text{obs}} \int_{f_\text{min}}^{f_\text{max}} df \frac{
     h^2 \Omega_{\text{gw},0}^2}{
     h^2 \Omega_{\text{n}}^2}}
-    :gowling_2021:`\ ` eq. 3.12
+    :gowling_2021:`\ ` eq. 3.12,
+    :smith_2019:`\ ` p. 1
 
     :param f: frequencies
     :param signal: signal array
@@ -27,22 +28,33 @@ def signal_to_noise_ratio(
 def ft(L: th.FloatOrArr = const.LISA_ARM_LENGTH) -> th.FloatOrArr:
     r"""Transfer frequency
     $$f_t = \frac{c}{2\pi L}$$
+    :gowling_2021:`\ ` p. 12
     """
     return const.c / (2*np.pi*L)
 
+FT_LISA = ft()
+_ft_func = ft
 
-# def N_AE(f: th.FloatOrArr, L: th.FloatOrArr = const.LISA_ARM_LENGTH) -> th.FloatOrArrNumba:
-#     r"""A and E channels of LISA instrument noise
-#     $$N_A = N_E = ...$$
-#     :gowling_2021:`\ ` eq. 3.4
-#     """
-#     f_frac = f/ft
-#     cos_f_frac = np.cos(f_frac)
-#     W = 1 - np.exp(-2j*f_frac)
-#     return ((4 + 2*cos_f_frac)*P_oms(L) + 8*(1 + cos_f_frac + cos_f_frac**2 * P_acc(f, L))) * np.abs(W)**2
+
+def N_AE(f: th.FloatOrArr, ft: th.FloatOrArr = FT_LISA, L: th.FloatOrArr = const.LISA_ARM_LENGTH, W_abs2: th.FloatOrArr = None) -> th.FloatOrArrNumba:
+    r"""A and E channels of LISA instrument noise
+    $$N_A = N_E = ...$$
+    :gowling_2021:`\ ` eq. 3.4
+    """
+    cos_f_frac = np.cos(f/ft)
+    if W_abs2 is None:
+        W_abs2 = np.abs(W(f, ft))**2
+    return ((4 + 2*cos_f_frac)*P_oms(L) + 8*(1 + cos_f_frac + cos_f_frac**2) * P_acc(f, L)) * W_abs2
 
 
 def omega(f: th.FloatOrArr, S: th.FloatOrArr) -> th.FloatOrArr:
+    r"""Convert an effective noise power spectral density (aka. sensitivity) $S$
+    to a fractional GW energy density power spectrum $\Omega$
+    $$\Omega = \frac{4 \pi^2}{3 H_0^2} f^3 S(f)$$
+    :gowling_2021:`\ ` eq. 3.8,
+    :gowling_2023:`\ ` eq. 3.8,
+    :smith_2019:`\ ` p. 1
+    """
     return 4*np.pi**2 / (3*const.H0_HZ**2) * f**3 * S
 
 
@@ -58,10 +70,11 @@ def omega_ins(f: th.FloatOrArr) -> th.FloatOrArr:
     r"""LISA instrument noise
     $$\Omega_\text{ins} = \left( \frac{4 \pi^2}{3 H_0^2} f^3 S_A(f)$$
     """
-    return omega(f, S_AE(f))
+    return omega(f=f, S=S_AE(f))
 
 
 def omega_noise(f: th.FloatOrArr) -> th.FloatOrArr:
+    r"""$$\Omega_\text{noise} = \Omega_\text{ins} + \Omega_\text{eb} + \Omega_\text{gb}$$"""
     return omega_ins(f) + omega_eb(f) + omega_gb(f)
 
 
@@ -79,13 +92,46 @@ def P_oms(L: th.FloatOrArr = const.LISA_ARM_LENGTH) -> th.FloatOrArr:
     $$P_\text{oms}(f) = \left( \frac{1.5 \cdot 10^{-11} \text{m}}{L} \right)^2 \text{Hz}^{-1}$$
     :gowling_2021:`\ ` eq. 3.2
     This is white noise and therefore independent of the frequency.
+    Note that there is a typo on :gowling_2021:`\ ` p. 12:
+    the correct $L = 2.5 \cdot 10^9 \text{m}$.
     """
     return (1.5e-11 / L)**2
 
 
-def S_AE(f: th.FloatOrArr, L: th.FloatOrArr = const.LISA_ARM_LENGTH, both_channels: bool = True) -> th.FloatOrArr:
+def R_AE(f: th.FloatOrArr, ft: th.FloatOrArr = FT_LISA, W_abs2: th.FloatOrArr = None) -> th.FloatOrArr:
+    r"""Gravitational wave response function for the A and E channels
+    :gowling_2021:`\ ` eq. 3.6
+    """
+    if W_abs2 is None:
+        W_abs2 = np.abs(W(f, ft))**2
+    return 9/20 * W_abs2 / (1 + (3*f/(4*ft))**2)
+
+
+def S(N: th.FloatOrArr, R: th.FloatOrArr) -> th.FloatOrArr:
+    r"""Noise power spectral density
+    $$S = \frac{N}{\mathcal{R}}$$
+    :gowling_2021:`\ ` eq. 3.1
+    """
+    return N / R
+
+
+def S_AE(f: th.FloatOrArr, ft: th.FloatOrArr = FT_LISA, L: th.FloatOrArr = const.LISA_ARM_LENGTH, both_channels: bool = True) -> th.FloatOrArr:
+    r"""Noise power spectral density for the LISA A and E channels
+    $$S_A = S_E = \frac{N_A}{\mathcal{R}_A}$$
+    :gowling_2021:`\ ` eq. 3.7
+    """
+    # The W_abs2 cancels and can therefore be set to unity
+    ret = S(N=N_AE(f=f, ft=ft, L=L, W_abs2=1), R=R_AE(f=f, ft=ft, W_abs2=1))
+    if both_channels:
+        return 1/np.sqrt(2) * ret
+    return ret
+
+
+def S_AE_approx(f: th.FloatOrArr, L: th.FloatOrArr = const.LISA_ARM_LENGTH, both_channels: bool = True) -> th.FloatOrArr:
     r"""Approximate noise power spectral density for the LISA A and E channels
-    $$S_A = S_E$
+    $$S_A = S_E = \frac{N_A}{\mathcal{R}_A}
+    \approx \frac{40}{3} (P_\text{oms} + 4P_\text{acc}) \left( 1 + \frac{3f}{4f_t} \right^2$$
+    :gowling_2021:`\ ` eq. 3.7
     """
     ret = 40/3 * (P_oms(L) + 4*P_acc(f, L)) * (1 + (3*f/(4*ft(L)))**2)
     if both_channels:
@@ -104,3 +150,11 @@ def S_gb(
         d: float = 1680) -> th.FloatOrArr:
     r"""Noise power spectral density for galactic binaries"""
     return A * (1e-3 / f)**(-7/3) * np.exp(-(f/f_ref_gb)**a - b*f*np.sin(c*f)) * (1 + np.tanh(d*(fk - f)))
+
+
+def W(f: th.FloatOrArr, ft: th.FloatOrArr) -> th.FloatOrArr:
+    r"""Round trip modulation
+    $$W(f,f_t) = 1 - e^{-2i \frac{f}{f_t}}$$
+    :gowling_2021:`\ ` p. 12
+    """
+    return 1 - np.exp(-2j * f / ft)