From 9cbba893e543a4b0389883af2c95fdf80663bea5 Mon Sep 17 00:00:00 2001
From: Kevin Kuns <kevin.kuns@ligo.org>
Date: Thu, 14 Jan 2021 16:16:41 -0500
Subject: [PATCH] add infinite volume residual gas damping

---
 gwinc/ifo/aLIGO/__init__.py |  9 +++--
 gwinc/ifo/aLIGO/ifo.yaml    |  2 +
 gwinc/ifo/noises.py         | 80 +++++++++++++++++++++++++++++++------
 gwinc/noise/residualgas.py  | 43 ++++++++++++++++++--
 4 files changed, 114 insertions(+), 20 deletions(-)

diff --git a/gwinc/ifo/aLIGO/__init__.py b/gwinc/ifo/aLIGO/__init__.py
index 821e4887..961f88fa 100644
--- a/gwinc/ifo/aLIGO/__init__.py
+++ b/gwinc/ifo/aLIGO/__init__.py
@@ -30,9 +30,12 @@ class ExcessGas(nb.Budget):
     )
 
     noises = [
-        ExcessGasH2,
-        ExcessGasN2,
-        ExcessGasH2O,
+        ExcessGasScatteringH2,
+        ExcessGasScatteringN2,
+        ExcessGasScatteringH2O,
+        ExcessGasDampingH2,
+        ExcessGasDampingN2,
+        ExcessGasDampingH2O,
     ]
 
 
diff --git a/gwinc/ifo/aLIGO/ifo.yaml b/gwinc/ifo/aLIGO/ifo.yaml
index eca5e990..28666445 100644
--- a/gwinc/ifo/aLIGO/ifo.yaml
+++ b/gwinc/ifo/aLIGO/ifo.yaml
@@ -37,11 +37,13 @@ Infrastructure:
       pressure: 4.0e-7              # Pa
       mass: 3.35e-27                # kg; Mass of H_2 (ref. 10)
       polarizability: 7.8e-31       # m^3
+      TestMassPressure: 4e-7
 
     N2:
       pressure: 1.33e-8
       mass: 4.65e-26
       polarizability: 1.79e-30
+      TestMassPressure: 4.65e-8
 
     H2O:
       pressure: 1.33e-8
diff --git a/gwinc/ifo/noises.py b/gwinc/ifo/noises.py
index 17f83d4a..f3467595 100644
--- a/gwinc/ifo/noises.py
+++ b/gwinc/ifo/noises.py
@@ -764,7 +764,7 @@ class ExcessGas(nb.Noise):
 
     def calc(self):
         cavity = arm_cavity(self.ifo)
-        n = noise.residualgas.residual_gas_cavity(
+        n = noise.residualgas.residual_gas_scattering(
             self.freq, self.ifo, cavity, self.ifo.Infrastructure.ResidualGas)
         # FIXME HACK: it's unclear if a phase noise in the arms like
         # the excess gas noise should get the same dhdL strain
@@ -779,55 +779,109 @@ class ExcessGas(nb.Noise):
         return n * 2 / sinc_sqr
 
 
-class ExcessGasH2(nb.Noise):
-    """Excess gas for H2
+class ExcessGasScatteringH2(nb.Noise):
+    """Excess gas scattering for H2
 
     """
     style = dict(
-        label='H2',
+        label='H2 scattering',
         color='xkcd:red orange'
     )
 
     def calc(self):
         cavity = arm_cavity(self.ifo)
         species = self.ifo.Infrastructure.ResidualGas.H2
-        n = noise.residualgas.residual_gas_cavity(
+        n = noise.residualgas.residual_gas_scattering(
             self.freq, self.ifo, cavity, species)
         dhdl_sqr, sinc_sqr = dhdl(self.freq, self.ifo.Infrastructure.Length)
         return n * 2 / sinc_sqr
 
 
-class ExcessGasN2(nb.Noise):
-    """Excess gas for N2
+class ExcessGasScatteringN2(nb.Noise):
+    """Excess gas scattering for N2
 
     """
     style = dict(
-        label='N2',
+        label='N2 scattering',
         color='xkcd:emerald'
     )
 
     def calc(self):
         cavity = arm_cavity(self.ifo)
         species = self.ifo.Infrastructure.ResidualGas.N2
-        n = noise.residualgas.residual_gas_cavity(
+        n = noise.residualgas.residual_gas_scattering(
             self.freq, self.ifo, cavity, species)
         dhdl_sqr, sinc_sqr = dhdl(self.freq, self.ifo.Infrastructure.Length)
         return n * 2 / sinc_sqr
 
 
-class ExcessGasH2O(nb.Noise):
-    """Excess gas for H2O
+class ExcessGasScatteringH2O(nb.Noise):
+    """Excess gas scattering for H2O
 
     """
     style = dict(
-        label='H2O',
+        label='H2O scattering',
         color='xkcd:water blue'
     )
 
     def calc(self):
         cavity = arm_cavity(self.ifo)
         species = self.ifo.Infrastructure.ResidualGas.H2O
-        n = noise.residualgas.residual_gas_cavity(
+        n = noise.residualgas.residual_gas_scattering(
             self.freq, self.ifo, cavity, species)
         dhdl_sqr, sinc_sqr = dhdl(self.freq, self.ifo.Infrastructure.Length)
         return n * 2 / sinc_sqr
+
+
+class ExcessGasDampingH2(nb.Noise):
+    """Residual gas damping for H2
+
+    """
+    style = dict(
+        label='H2 damping',
+        color='xkcd:red orange',
+        linestyle='--',
+    )
+
+    @nb.precomp(sustf=precomp_suspension)
+    def calc(self, sustf):
+        species = self.ifo.Infrastructure.ResidualGas.H2
+        n = noise.residualgas.residual_gas_damping(
+            self.freq, self.ifo, species, sustf)
+        return 4 * n
+
+
+class ExcessGasDampingN2(nb.Noise):
+    """Excess gas damping for N2
+
+    """
+    style = dict(
+        label='N2 damping',
+        color='xkcd:emerald',
+        linestyle='--',
+    )
+
+    @nb.precomp(sustf=precomp_suspension)
+    def calc(self, sustf):
+        species = self.ifo.Infrastructure.ResidualGas.N2
+        n = noise.residualgas.residual_gas_damping(
+            self.freq, self.ifo, species, sustf)
+        return 4 * n
+
+
+class ExcessGasDampingH2O(nb.Noise):
+    """Excess gas damping for H2O
+
+    """
+    style = dict(
+        label='H2O damping',
+        color='xkcd:water blue',
+        linestyle='--',
+    )
+
+    @nb.precomp(sustf=precomp_suspension)
+    def calc(self, sustf):
+        species = self.ifo.Infrastructure.ResidualGas.H2O
+        n = noise.residualgas.residual_gas_damping(
+            self.freq, self.ifo, species, sustf)
+        return 4 * n
diff --git a/gwinc/noise/residualgas.py b/gwinc/noise/residualgas.py
index 61cb962f..bdf31f5e 100644
--- a/gwinc/noise/residualgas.py
+++ b/gwinc/noise/residualgas.py
@@ -7,18 +7,18 @@ from numpy import sqrt, log, pi
 from .. import const
 
 
-def residual_gas_cavity(f, ifo, cavity, species):
-    """Residual gas noise strain spectrum in an evacuated cavity.
+def residual_gas_scattering(f, ifo, cavity, species):
+    """Residual gas noise strain spectrum due to scattering.
 
     Noise caused by the passage of residual gas molecules through the
-    laser beams in a cavity.
+    laser beams in a cavity due to scattering.
 
     :f: frequency array in Hz
     :ifo: gwinc IFO structure
     :cavity: arm cavity structure
     :species: molecular species structure
 
-    :returns: arm displacement noise power spectrum at :f:
+    :returns: arm strain noise power spectrum at :f:
 
     The method used here is presented by Rainer Weiss, Micheal
     E. Zucker, and Stanley E. Whitcomb in their paper Optical
@@ -55,3 +55,38 @@ def residual_gas_cavity(f, ifo, cavity, species):
     zint[zint < 0] = 0
 
     return zint
+
+
+def residual_gas_damping(f, ifo, species, sustf):
+    """Noise due to residual gas damping for one test mass
+
+    :f: frequency array in Hz
+    :ifo: gwinc IFO structure
+    :species: molecular species structure
+    :sustf: suspension transfer function structure
+
+    :returns: displacement noise
+    """
+    kT = ifo.Infrastructure.Temp * const.kB
+    mass = species.mass
+    radius = ifo.Materials.MassRadius
+    thickness = ifo.Materials.MassThickness
+    thermal_vel = sqrt(kT/mass)  # thermal velocity
+
+    # pressure near the test mass
+    # possibly different from the pressure in the arms due to outgassing near
+    # the test mass
+    if 'TestMassPressure' in species:
+        pressure = species.TestMassPressure
+    else:
+        pressure = species.pressure
+
+    # infinite volume viscous damping coefficient for a cylinder
+    # table 1 of https://doi.org/10.1016/j.physleta.2010.06.041
+    beta_inf = pi * radius**2 * pressure/thermal_vel
+    beta_inf *= (1 + thickness/(2*radius) + pi/4)
+
+    # convert force to displacement noise using the suspension susceptibility
+    noise = 4 * kT * beta_inf * abs(sustf.tst_suscept)**2
+
+    return noise
-- 
GitLab