Function: simulate_sb()

0

void simulate_sb(int n, float delta_t, float omega_0, float f_low, float f_high, double *gamma12, double *whiten1, double *whiten2, float *out1, float *out2, int *pseed)
This function simulates the generation of an isotropic and unpolarized stochastic background of gravitational radiation having a constant frequency spectrum: $\Omega_{\rm gw}(f)=\Omega_0$ for $f_{\rm low}\le f \le f_{\rm high}$. The outputs are (not necessarily continuous-in-time) whitened data stream $o_1(t)$ and $o_2(t)$ representing the detector outputs when only a stochastic background signal is present.

The arguments of simulate_sb() are:

n: Input. The number $N$ of data points corresponding to an observation time $T:=N\ \Delta t$, where $\Delta t$ is the sampling period of the detectors, defined below. $N$ should equal an integer power of 2.

delta_t: Input. The sampling period $\Delta t$ (in sec) of the detectors.

omega_0: Input. The constant value $\Omega_0$ (dimensionless) of the frequency spectrum $\Omega_{\rm gw}(f)$ for the stochastic background:

$$\Omega_{\rm gw}(f)=\left\{ \begin{array}{cl} \Omega_0 & \qua... ... f_{\rm high}\0 & \quad {\rm otherwise.} \end{array}\right.$$

$\Omega_0$ should be greater than or equal to zero.

f_low: Input. The frequency $f_{\rm low}$ (in Hz) below which the spectrum $\Omega_{\rm gw}(f)$ of the stochastic background is zero. $f_{\rm low}$ should lie in the range $0\le f_{\rm low}\le f_{\rm Nyquist}$, where $f_{\rm Nyquist}$ is the Nyquist critical frequency. (The Nyquist critical frequency is defined by $f_{\rm Nyquist}:=1/(2\Delta t)$, where $\Delta t$ is the sampling period of the detectors.) $f_{\rm low}$ should also be less than or equal to $f_{\rm high}$.

f_high: Input. The frequency $f_{\rm high}$ (in Hz) above which the spectrum $\Omega_{\rm gw}(f)$ of the stochastic background is zero. $f_{\rm high}$ should lie in the range $0\le f_{\rm high}\le f_{\rm Nyquist}$. It should also be greater than or equal to $f_{\rm low}$.

gamma12: Input. gamma12[0..n/2-1] is an array of double precision variables containing the values of the overlap reduction function $\gamma(f)$ for the two detector sites. These variables are dimensionless. gamma12[i] contains the value of $\gamma(f)$ evaluated at the discrete frequency $f_i=i/(N\Delta t)$, where $i=0,1,\cdots,N/2-1$.

whiten1: Input. whiten1[0..n-1] is an array of double precision variables containing the values of the real and imaginary parts of the spectrum $\tilde W_1(f)$ of the whitening filter of the first detector. These variables have units ${\rm rHz}/{\rm strain}$ (or ${\rm sec}^{-1/2}$), which are inverse to the units of the square root of the noise power spectrum $P_1(f)$. whiten1[2*i] and whiten1[2*i+1] contain, respectively, the values of the real and imaginary parts of $\tilde W_1(f)$ evaluated at the discrete frequency $f_i=i/(N\Delta t)$, where $i=0,1,\cdots,N/2-1$.

whiten2: Input. whiten2[0..n-1] is an array of double precision variables containing the values of the real and imaginary parts of the spectrum $\tilde W_2(f)$ of the whitening filter of the second detector, in exactly the same format as the previous argument.

out1: Output. out1[0..n-1] is an array of floating point variables containing the values of the whitened data stream $o_1(t)$ representing the output of the first detector when only a stochastic background signal is present. $o_1(t)$ is the convolution of detector whitening filter $W_1(t)$ with the gravitational strain $h_1(t)$. The variables out1[] have units of rHz (or ${\rm sec}^{-1/2}$), which follows from the definition of $h_1(t)$ as a strain and $\tilde W_1(f)$ as the ``inverse'' of the square root of the noise power spectrum $P_1(f)$. out1[i] contains the value of $o_1(t)$ evaluated at the discrete time $t_i=i\Delta t$, where $i=0,1,\cdots,N-1$.

out2: Output. out2[0..n-1] is an array of floating point variables containing the values of the whitened data stream $o_2(t)$ representing the output of the second detector when only a stochastic background signal is present, in exactly the same format as the previous argument.

pseed: Input. A pointer to a seed value, which is used by the random number generator routine.

simulate_sb() simulates the generation of an isotropic and unpolarized stochastic background of gravitational radiation having a constant frequency spectrum $\Omega_{\rm gw}(f)=\Omega_0$ for $f_{\rm low}\le f \le f_{\rm high}$ in the following series of steps:

(i)

It first constructs random variables $\tilde h_1(f_i)$ and $\tilde h_2(f_i)$ in the frequency domain that have zero mean and satisfy:

$$\langle \tilde h_1^*(f_i)\tilde h_1(f_j)\rangle = {1\over 2}\ T\ \delta_{ij}\ {3H_0^2\over 10\pi^2}\ f_i^{-3}\ \Omega_0$$ (11.10.231)

$$\langle \tilde h_2^*(f_i)\tilde h_2(f_j)\rangle = {1\over 2}\ T\ \delta_{ij}\ {3H_0^2\over 10\pi^2}\ f_i^{-3}\ \Omega_0$$ (11.10.232)

$$\langle \tilde h_1^*(f_i)\tilde h_2(f_j)\rangle = {1\over 2}\ T\ \delta_{ij}\ {3H_0^2\over 10\pi^2}\ f_i^{-3}\ \Omega_0 \ \gamma(f_i)\ ,$$ (11.10.233)

where $\langle\ \rangle$ denotes ensemble average. Here $\tilde h_1(f_i)$ and $\tilde h_2(f_i)$ are the Fourier components of the gravitational strains $h_1(t)$ and $h_2(t)$ at the two detectors. The above equations are the discrete frequency versions of equation (3.17) of Ref. [36], with $\Omega_{\rm gw}(f)=\Omega_0$ for $f_{\rm low}\le f \le f_{\rm high}$. They can be realized by setting

$$\tilde h_1(f_i) = {1\over 2}\sqrt{T}\ \left({3H_0^2\over 10\pi^2}\right)^{1/2} f_i^{-3/2}\ \Omega_0^{1/2}\ (x_{1i} + i y_{1i})$$ (11.10.234)

$$\tilde h_2(f_i) = \tilde h_1(f_i)\ \gamma(f_i)+$$ (11.10.235)

$${1\over 2}\sqrt{T}\ \left({3H_0^2\over 10\pi^2}\right)^{1/2} f_i^{-3/2}\ \ \Omega_0^{1/2}\sqrt{1-\gamma^2(f_i)}\ (x_{2i}+i y_{2i})\ ,$$ (11.10.236)

where $x_{1i}$, $y_{1i}$, $x_{2i}$, and $y_{2i}$ are statistically independent (real) Gaussian random variables, each having zero mean and unit variance. (Note: The $x_{1i}$, $y_{1i}$, $x_{2i}$, and $y_{2i}$ random variables are statistically independent of the $u_i$ and $v_i$ random variables defined in Sec. .) These Gaussian random variables are produced by calls to the Numerical Recipes in C random number generator routine gasdev(). Note also that the second term of $\tilde h_2(f_i)$ (which is proportional to $\sqrt{1-\gamma^2(f_i)}$) is needed to obtain equation (). Without this term, $\langle\tilde h_2^*(f_i)\tilde h_2(f_j)\rangle$ would include an additional (unwanted) factor of $\gamma^2(f_i)$.

(iii)

simulate_sb() then whitens the data in the frequency domain by multiplying $\tilde h_1(f_i)$ and $\tilde h_2(f_i)$ by the frequency components $\tilde W_1(f_i)$ and $\tilde W_2(f_i)$ of the whitening filters of the two detectors:

$$\tilde o_1(f_i) := \tilde h_1(f_i)\ \tilde W_1(f_i)$$ (11.10.237)

$$\tilde o_2(f_i) := \tilde h_2(f_i)\ \tilde W_2(f_i)\ .$$ (11.10.238)

This (complex) multiplication in the frequency domain corresponds to the convolution of $h_1(t)$ and $W_1(t)$, and $h_2(t)$ and $W_2(t)$ in the time domain. By convention, the DC (i.e., zero frequency) and Nyquist critical frequency components of $\tilde o_1(f_i)$ and $\tilde o_2(f_i)$ are set to zero.

(iii)

The final step consists of Fourier transforming the frequency components $\tilde o_1(f_i)$ and $\tilde o_2(f_i)$ into the time domain to obtain the whitened data streams $o_1(t_i)$ and $o_2(t_i)$. Here $t_i=i\Delta t$ with $i=0,1,\cdots,N-1$. Since $\tilde o_1^*(f_i)$ and $\tilde o_2^*(f_i)$ are the Fourier transforms of real data sets, the two Fourier transforms can be performed simultaneously.

Authors: Bruce Allen, ballen@dirac.phys.uwm.edu, and Joseph Romano, romano@csd.uwm.edu

Comments: Although it is possible and more efficient to write a single function to simulate the generation of a stochastic background and intrinsic detector noise simultaneously, we have chosen--for the sake of modularity--to write separate functions to perform these two tasks separately. (See also the comment at the end of Sec. .)