引力波天文IV：数据分析方法

引力波天文系列文章将介绍引力波能观测到什么，如何观测，以及基本的数据分析方法。侧重于引力波的空间探测和分析，但核心的测量原理、仪器响应以及数据分析方法基本都是通用的。

• 研究目标：波源概述(基本波源及频谱划分)
• 探测技术：响应与灵敏度、空间探测技术
• 数据分析：数据分析方法、数据分析实践

本文介绍引力波数据的分析方法。

数据分析基础

上面推导的频域响应是严格的，但只针对特定的时刻：

$$s(t)|_{t_{\rm orbit}} = \mathcal{F}^{-1}\{ \tilde{s}(f)|_{t_{\rm orbit}}\}$$

只有$t=t_{\rm orbit}$的值才对应实际观测。随着$t_{\rm orbit}$改变，探测器位置改变，$\tilde{s}(f)$改变。
实际观测数据$\bar{s}(t_{\rm orbit})$是一系列不同位置处$s(t)$在对应时刻取值所组成的序列，对它做傅里叶变换得到的是？怎么得到前面推导的$\tilde{s}(f)$？

注意：对于双星绕转或并和等系统，$s(t)$是有傅里叶变换的。但探测器噪声$n(t)$以及随机背景信号$s(t)$都是随机过程，属于能量信号，严格讲是没有傅里叶变换的(积分发散)，只有功率谱，这也是为什么会引入特征应变$h_c(f)$。这里只是个记号，不求严谨。

$$\begin{aligned} d(t) = s(t) + n(t); ~~~ s(t) \equiv F^P(t) * h_P(t)\\ \tilde{d}(f) = \tilde{s}(f) + \tilde{n}(f); ~~~ \tilde{s}(f) = F^P(f) \tilde{h}_P(f) \end{aligned}$$

信噪比

$$\begin{aligned} d(t) = F^P(t) * h_P(t) + n(t)\\ \tilde{d}(f) = F^P(f) \tilde{h}_P(f) + \tilde{n}(f) \end{aligned}$$

在描述引力波信号时，不同偏振模式不能简单合并，有各自不同的转移函数，而且一般地$F^P$是复值函数。但不同引力波信号的功率(能量)可以直接相加：

$$\begin{aligned} S_h(f) = \sum_P |\tilde{h}_P(f)|^2 \leftrightarrow S_h(f) = \frac{1}{2}\sum_P |\tilde{h}_P(f)|^2 \end{aligned}$$

对各态历经的平稳随机过程，有两种等价的自(互)相关：

• 一种是系综平均：$E\{ d_I(t)d_J(t+\tau) \} ~~~ \big\langle d_I(t)d_J(t+\tau)\big\rangle$
• 一种是时间平均：$\lim_{T\rightarrow\infty} \frac{1}{T}\int_{-T/2}^{T/2} d_I(t)d_J(t+\tau) dt ~~~~ d_I(t) \star d_J(t+\tau)$
  相比函数的互相关$\int_{-\infty}^{\infty} f(t) g(t+\tau) dt$，多了对时间的平均？？为什么
  不同于一般互相关定义，这里互相关是包含时间平均的，具体参考Creighton & Anderson. 2011.

$$\langle x(t)x(t+\tau)\rangle = r^x(\tau), ~~~~ \langle \tilde{x}(f)\tilde{x}^*(f')\rangle = S(f)\delta(f-f')$$

$$\langle x^2(t)\rangle = r^x(0) = \int S(f)df, ~~~~ \langle |\tilde{x}(f)|^2\rangle = S(f)\delta(0)$$

通过探测信号自相关或互相关计算(互)功率谱，与
自相关$r^d(t_1, t_2) = E\{ d(t_1) d(t_2) \}; ~~ r^s(t_1, t_2) = E\{ s(t_1) s(t_2) \}$
互相关$r^d_{IJ}(t_1, t_2) = E\{ d_I(t_1) d_J(t_2) \}; ~~ r^s_{IJ}(t_1, t_2) = E\{ s_I(t_1) s_J(t_2) \}$
其中期望对应系综平均

对平稳随机过程自相关函数只与时间差有关，不随具体计算的时间点变化。
对各态历经过程，系综平均等于时间平均

自相关$r^s(t) = r^s(t_1-t_2) = \left\langle s(t) s(t) \right\rangle,~~~ S(f) = \tilde{s}(f)\tilde{s}^*(f) = |\tilde{s}(f)|^2$
互相关$r^s_{IJ}(t) = r^s_{IJ}(t_1 - t_2) = \left\langle s_I(t) \star s_J(t) \right\rangle,~~~ S_{IJ}(f) = \tilde{s}_I(f)\tilde{s}_J^*(f)$

$$r^s(t) = \int S(f) e^{i 2 \pi f t} df,~~~ \left\langle s^2(t) \right\rangle = r^s(0) = \int N(f) df$$

$$\left\langle s(t) \star s(t) \right\rangle \leftrightarrows \left\langle \tilde{s}(f)\tilde{s}^*(f) \right\rangle \leftrightarrows S(f)$$

$$r^s_{IJ}(t) = \int S_{IJ}(f) e^{i 2 \pi f t} df,~~~ \left\langle s_I(t)s_J(t) \right\rangle = r^s_{IJ}(0) = \int S_{IJ}(f) df$$

对于随机背景如果不做时间延迟，相关就是信号相乘并对时间积分，由于是随机过程，需要求系综平均，即多次测量求平均。噪声不相关求平均为0。

$$\left\langle d_I(t) d_J(t) \right\rangle = \left\langle s_I(t) s_J(t) \right\rangle$$

$$\begin{aligned} \left\langle d_I(t) d_J(t) d_I(t') d_J(t') \right\rangle =& \left\langle s_I(t)s_I(t') s_J(t)s_J(t') \right\rangle + \left\langle n_I(t) n_I(t') n_J(t)n_J(t') \right\rangle\\ &+ \left\langle s_I(t)s_I(t') n_J(t)n_J(t') \right\rangle + \left\langle n_I(t)n_I(t') s_J(t)s_J(t') \right\rangle \end{aligned}$$

$$\left\langle s_I(t)s_I(t') s_J(t)s_J(t') \right\rangle \leftrightarrow \left\langle s_I(t) s_J(t) \right\rangle\left\langle s_I(t') s_J(t') \right\rangle \leftrightarrow \left\langle s_I(t) s_I(t') \right\rangle\left\langle s_J(t) s_J(t') \right\rangle$$

$$\begin{aligned} &\left\langle d_I(t) d_J(t) d_I(t') d_J(t') \right\rangle - \left\langle s_I(t) s_J(t) \right\rangle\left\langle s_I(t') s_J(t') \right\rangle \\ &=\left\langle n_I(t) n_I(t') n_J(t)n_J(t') \right\rangle+ \left\langle s_I(t)s_I(t') n_J(t)n_J(t') \right\rangle + \left\langle n_I(t)n_I(t') s_J(t)s_J(t') \right\rangle\\ &=\left\langle n_I(t) n_I(t')\right\rangle \left\langle n_J(t)n_J(t') \right\rangle + \left\langle s_I(t)s_I(t')\right\rangle \left\langle n_J(t)n_J(t') \right\rangle + \left\langle n_I(t)n_I(t')\right\rangle \left\langle s_J(t)s_J(t') \right\rangle\\ \end{aligned}$$

在弱信号假设下$\left\langle s(t)s(t')\right\rangle \ll \left\langle n(t)n(t') \right\rangle$

$$\left\langle d_I(t) d_J(t) d_I(t') d_J(t') \right\rangle - \left\langle d_I(t) d_J(t) \right\rangle^2 \approx \left\langle n_I(t) n_I(t') \right\rangle \left\langle n_J(t)n_J(t') \right\rangle$$

考虑到引力波信号与噪声不相关：

$$r(t) = d \star d = s \star s + n \star n; ~~~ r_{IJ}(t) = d_I \star d_J = s_I \star s_J + n_I \star n_J$$

如果进一步假设不同探测器噪声不相关：

$$r(t) = s \star s + n \star n; ~~~ r_{IJ}(t) = s_I \star s_J$$

$$\left\langle d_1 d_2 \right\rangle = \left\langle s_1 s_2\right\rangle + \left\langle n_1 n_2\right\rangle + \left\langle s_1 n_2\right\rangle + \left\langle n_1 s_2\right\rangle$$

轨道平均 v.s. 方位平均
地基探测器：地球自转+公转
空间探测器：探测器自转+公转，卫星三角在探测器平面内自转，

信号功率谱描述

$$\begin{aligned} S(f) = |\tilde{s}(f)|^2 = \tilde{s}(f)\tilde{s}^*(f), ~~~S_h(f) = \textcolor{red}{\frac{1}{2}??} \sum_P |\tilde{h}_P(f)|^2 \end{aligned}$$

其中，$S_{\rm \tiny GW}(f)$为引力波信号原始功率为不同偏振模式信号功率之和。

由此可定义引力波的功率响应函数：

$$\boxed{\mathcal{R}(f) = \frac{S(f)}{S_h(f)} = \frac{\tilde{s}(f)\tilde{s}^*(f)}{\textcolor{red}{\frac{1}{2}??} \sum_P |\tilde{h}_P(f)|^2}}$$

通常说的引力波响应函数指的就是功率转移函数(power transfer function)

这里要有一个R(f)的曲线图，覆盖从PTA到地基观测的所有波段，参考Romano & Cornish 2017 P74

显然这里$\tilde{s}(f), S(f), \mathcal{R}(f)$都是依赖方位角及偏振角的函数。一般情况下，难以化简，但若对偏振角求平均

多探测器联合

$$\begin{aligned} \begin{pmatrix} F^{+}\\ F^{\times} \end{pmatrix} = \mathcal{R}_{2\psi} \begin{pmatrix} D^{+} \\ D^{\times} \end{pmatrix} = \begin{pmatrix} \cos 2\psi & -\sin 2\psi \\ \sin 2\psi & \cos 2\psi \end{pmatrix} \begin{pmatrix} D^{+} \\ D^{\times} \end{pmatrix} \end{aligned}$$

$$\begin{aligned} &\small F^{+} {F^{+}}^* = \cos^2 2\psi D^{+}{D^{+}}^* - \frac{1}{2}\sin 4\psi (D^{+} {D^{\times}}^* + {D^{+}}^*D^{\times} ) + \sin^2 2\psi D^{\times} {D^{\times}}^*\\ &\small F^{\times} {F^{\times}}^* = \sin^2 2\psi D^{+}{D^{+}}^* + \frac{1}{2}\sin 4\psi (D^{+} {D^{\times}}^* + {D^{+}}^*D^{\times} ) + \cos^2 2\psi D^{\times} {D^{\times}}^*\\ &\small F^{+} {F^{\times}}^* = \frac{1}{2}\sin 4\psi \left( D^{+}{D^{+}}^* - D^{\times}{D^{\times}}^* \right) + \cos^2 2\psi D^{+} {D^{\times}}^* - \sin^2 2\psi {D^{+}}^* D^{\times}\\ &\small {F^{+}}^* F^{\times} = \frac{1}{2}\sin 4\psi \left( D^{+}{D^{+}}^* - D^{\times}{D^{\times}}^* \right) + \cos^2 2\psi {D^{+}}^* D^{\times} - \sin^2 2\psi D^{+} {D^{\times}}^*\\ \end{aligned}$$

对于长周期信号，随着轨道运动
对轨道/探测器方位求平均$\left\langle D^{+} {D^{\times}}^* \right\rangle = \left\langle {D^{+}}^* D^{\times} \right\rangle$，从而有：

$$\begin{aligned} &\small \left\langle F^{+} {F^{+}}^*\right\rangle = \cos^2 2\psi \left\langle D^{+}{D^{+}}^*\right\rangle - \sin 4\psi \left\langle D^{+} {D^{\times}}^*\right\rangle + \sin^2 2\psi \left\langle D^{\times} {D^{\times}}^*\right\rangle\\ &\small \left\langle F^{\times} {F^{\times}}^*\right\rangle = \sin^2 2\psi \left\langle D^{+}{D^{+}}^*\right\rangle + \sin 4\psi \left\langle D^{+} {D^{\times}}^*\right\rangle + \cos^2 2\psi \left\langle D^{\times} {D^{\times}}^*\right\rangle\\ &\small \left\langle F^{+} {F^{\times}}^*\right\rangle = \left\langle {F^{+}}^* F^{\times}\right\rangle = \frac{1}{2}\sin 4\psi \left( \left\langle D^{+}{D^{+}}^*\right\rangle - \left\langle D^{\times}{D^{\times}}^*\right\rangle \right) + \cos 4\psi \left\langle D^{+} {D^{\times}}^*\right\rangle \end{aligned}$$

对长周期的单频信号(双白矮星)

进一步的，如果对偏振角求平均有：

$$\overline{ \left\langle F^P {F^{P'}}^* \right\rangle } = {\small \frac{ \left\langle {\left|D^{+}\right|}^2\right\rangle + \left\langle {\left|D^{\times}\right|}^2\right\rangle}{2} } \delta_{PP'}$$

$$\overline{ \left\langle F^+ {F^+}^* \right\rangle } = \overline{ \left\langle F^\times {F^\times}^* \right\rangle } = {\small \frac{ \left\langle {\left|D^{+}\right|}^2\right\rangle + \left\langle {\left|D^{\times}\right|}^2\right\rangle}{2} }$$

$$\overline{ \left\langle F^+ {F^\times}^* \right\rangle } = \overline{ \left\langle F^\times {F^+}^* \right\rangle } = 0$$

对方位角和偏振角平均后所有偏振交叉项都为零，且两种偏振模式的(平均)功率响应相同。

$$\tilde{s}(f)\tilde{s}^*(f) = \sum_P \overline{ \left\langle F^P {F^{P'}}^* \right\rangle } |\tilde{h}_P(f)|^2$$

由于对偏振角以及方位角平均后有$\overline{ \left\langle F^+ {F^+}^* \right\rangle } = \overline{ \left\langle F^\times {F^\times}^* \right\rangle }$，从而响应函数：

$$\mathcal{R}(f, \hat{k}) = \frac{\tilde{s}(f)\tilde{s}^*(f)}{\sum_P |\tilde{h}_P(f)|^2} = \overline{\left\langle F^+ {F^+}^*\right\rangle} = \overline{\left\langle F^\times {F^\times}^*\right\rangle} = {\small \frac{ \left\langle {\left|D^{+}\right|}^2\right\rangle + \left\langle {\left|D^{\times}\right|}^2\right\rangle}{2} }$$

$$\mathcal{R}(f, \hat{k}) = \frac{\tilde{s}(f)\tilde{s}^*(f)}{\textcolor{red}{\frac{1}{2}??} \sum_P |\tilde{h}_P(f)|^2} = \sum_P \overline{\left\langle F^P {F^P}^*\right\rangle} = \left\langle {\left|D^{+}\right|}^2\right\rangle + \left\langle {\left|D^{\times}\right|}^2\right\rangle$$

引力波频谱描述：$S_h, Ω, h_c$

检测与参数推断

最佳滤波/匹配滤波
信号探测能力/灵敏度曲线
Fisher信息矩阵
MCMC

探测器联合观测

联合分析()：多台探测器数据合并，最佳信噪比及对应的参数限制
利用所有探测器的数据纳入同一统计量，计算最佳匹配模板。“综合孔径”(Aperture Synthesis)
对于模板确定的波源，如致密双星并和系统，联合观测的信噪比就是各探测器之和，Fisher信息矩阵也是各探测器之和？。

$$\rho^2 = \sum_I \rho_I^2, \mathcal{I}_{ij} = \sum_I \mathcal{I}^I_{ij}$$

N个(全同)探测器的网络相比单探测器，最佳匹配滤波的信噪比提升为$\sqrt{N}$倍

定义

$$\mathbb{S}^{\rm tot} = \frac{\sum \lambda_i \mathbb{S}^{(i)}}{\sum_i \lambda_i}$$

其中 $\lambda_i$ 的取值需要使得整体信噪比 $\frac{\mu_{\mathbb{S}_{\rm tot}}}{\sigma_{\mathbb{S}_{\rm tot}}}$ 最大化。

$$\mu_{\mathbb{S}^{(i)}} = T^{(i)}\int df S_h(f) \Gamma_{IJ}(f), ~~~ \sigma_{\mathbb{S}^{(i)}} = T^{(i)}\int df N_I^{(i)}(f)N_J^{(i)}(f) \tilde{Q}(f)\tilde{Q}^*(f)$$

$$\left\langle \left(\sum \lambda_i \mathbb{S}^{(i)} \right)^2\right\rangle \rightarrow \sum \lambda^2_i \left\langle { \mathbb{S}^{(i)}}^2\right\rangle$$

$$\mu_{\mathbb{S}^{\rm tot}} = \frac{\sum \lambda_i \mu_{\mathbb{S}^{(i)}} }{\sum_i \lambda_i}, ~~~ \sigma_{\mathbb{S}^{\rm tot}} = \frac{\sqrt{\sum \lambda^2_i \sigma^2_{\mathbb{S}^{(i)}} } }{\sum_i \lambda_i}$$

$$\rho = \frac{ \mu_{\mathbb{S}^{\rm tot}} }{ \sigma_{\mathbb{S}^{\rm tot}} } = \frac{ \sum_i \lambda_i \mu_{\mathbb{S}^{(i)}} }{ \sqrt{\sum_i \lambda^2_i \sigma^2_{\mathbb{S}^{(i)}} } } = \frac{\left\langle \lambda \big| \sigma^{-2} \right\rangle}{\sqrt{\left\langle \lambda \big| \lambda \right\rangle}}$$

其中$\left\langle a \big| b \right\rangle \equiv \sum_i a_i b_i \sigma_i^2 T^{(i)}$，当$\lambda_i = \sigma_i^{-2}$时总信噪比最大：

$$\rho = \sqrt{\sum_i \frac{ T^{(i)}\mu }{ \sigma^2_i}} = \sqrt{\sum_i {\rho_{\tiny IJ}^{(i)}}^2}$$

三角定位：
直观定性的，两个探测器的空间定位是天球面上的一个环带，三台探测器时两个环带交叉，对应天球面上两个不同方位的小区域。
具体定量的，从Fisher参数推断角度理解，Fisher信息矩阵是参数协方差矩阵的倒数，总参数协方差为

$$\Sigma_{ij} = \frac{\prod \Sigma^I_{ij}}{\sum \Sigma^I_{ij}}$$

在参数空间中参数等概率曲面为椭球，参数观测误差越大，对应方向的轴越长。
协方差是相乘关系，因此只要不同探测器的参数不存在兼并，即参数长轴方向有区别，联合分析就可以显著提升对参数误差的限制。

一致性分析(Coincidence Analysis)：独立分析各探测数据，对比确认探测器中一致的候选事件。

相关分析：相干coherence v.s. 相关correlation
上面讨论的都是数据与模板匹配，这里是数据与数据匹配

$$\int_{-T/2}^{T/2} dt_1 \int_{-T/2}^{T/2} dt_2 ~ d_I(t_1) d_J(t_2) Q(t_1 - t_2)$$

参考随机噪声，最佳滤波信噪比：

$$\rho_{\tiny IJ} = \langle \mathcal{R} \rangle \sqrt{T} \sqrt{\int df \frac{ \gamma_{IJ}^2(f) S^2_h(f)}{N_I(f) N_J(f)} }$$

$$S_h(f) = f |\tilde{h}_c(f)|^2 = \frac{1}{\langle \mathcal{R} \rangle} f |\tilde{s}_c(f)|^2$$

$$\rho_{\tiny IJ} = \sqrt{T} \sqrt{\int df ~ f^2 \gamma_{IJ}^2(f) \frac{ |\tilde{s}_c(f)|^2}{N_I(f)} \frac{ |\tilde{s}_c(f)|^2}{N_J(f)} }$$

$f^2$项开根号后与$\sqrt{T}$也抵消不了，量纲不对？？？
这里是对角向积分过后的结果，对于点源而言是不需要的！！
这里的相关分析和前面的联合分析的区别和联系，信噪比会比匹配滤波显著低吗？？？
同样需要知道$S_h(f)$的先验，以确定最佳滤波；重叠约化函数
对比匹配滤波的信噪比

$$\rho_{\tiny IJ} = \sqrt{\int df \frac{ |\tilde{s}(f)|^2}{N(f)} }, |\tilde{s}(f)|^2 = \mathcal{R}(f)S_h(f)$$

随机过程的相关函数
对各态历经的平稳随机过程，有两种等价的自(互)相关：

• 一种是系综平均：$E\{ n_I(t)n_J(t+\tau) \}, \big\langle n_I(t)n_J(t+\tau)\big\rangle$
• 一种是时间平均：$\lim_{T\rightarrow\infty} \frac{1}{T}\int_{-T/2}^{T/2} n_I(t)n_J(t+\tau) dt$
  注意，相比确定性函数的互相关$\int_{-\infty}^{\infty} f(t) g(t+\tau) dt$，这里随机过程互相关多了对时间平均。从各态历经角度理解，系综平均等于时间平均，这里的$\frac{1}{T}$来自对时间平均；而直观地，不同于通常函数，对于随机过程，当$T\rightarrow\infty$积分是发散的，除以$T$整体结果才收敛。随机过程是功率信号，能量积分发散，需要除以时间得到功率；相对的，对确定性的能量信号，能量有限，反而不能除以(无限长的)时间。

考虑到随机过程能量发散，并不存在傅里叶变换，但形式上可记为$\tilde{n}(f)$。

随机过程能量发散(功率信号)，$\tilde{n}(f)$是发散的，但功率谱密度有限

$$\left\langle \tilde{n}(f) \tilde{n}^*(f') \right\rangle = \delta(f-f')S_n(f)$$

$$\left\langle \tilde{n}(f) \tilde{n}^*(f) \right\rangle = \delta(0)S_n(f) \rightarrow \infty$$

对于确定信号？？？

信噪比与误报率
假设目标信号服从正态分布，对于信噪比阈值，在统计上有：

$$\rho_{\rm threshold} = \sqrt{2}\Big({\rm erfcinv}(2\alpha_{\rm FAR})-{\rm erfcinv}[2(1-\beta_{\rm FNR})]\Big)$$

这里$\alpha_{\rm FAR},\beta_{\rm FNR}$分别为预设的误报率(False alarm rate)和漏报率(False negative rate)，$1-\beta_{\rm FNR}$对应于信号探测率。信号探测的误报和漏报率都取千分之一，则对应信噪比阈值约6.2，引力波观测中信噪比阈值通常取6 ~ 8。

对于随机背景，多探测器互相关，通常取信噪比阈值2，对应95%以上的置信度(误报率5%以下)，单探测器自相关，信噪比阈值则约定为5。

特征应变强度

$\delta$函数与$\rm sinc$函数

$$\begin{aligned} \delta(x) &= \lim_{\varepsilon\rightarrow\infty} \frac{1}{|\varepsilon|\sqrt{\pi} } e^{-\frac{x^2}{\varepsilon^2}} = \lim_{\varepsilon\rightarrow\infty} \frac{1}{2 \varepsilon} e^{-\frac{|x|}{\varepsilon}} = \lim_{\varepsilon\rightarrow0} \frac{\sin\left(\frac{x}{\varepsilon}\right)}{\pi x} = \lim_{\varepsilon\rightarrow 0} \frac{\sin\left(\frac{\pi x}{\varepsilon}\right)}{\pi x}\\ &= \lim_{\varepsilon\rightarrow 0} \frac{1}{\varepsilon}{\rm sinc}\left(\frac{x}{\varepsilon}\right) = \lim_{\varepsilon\rightarrow \infty} \frac{1}{\varepsilon} ~ {\rm rect}\left(\frac{x}{\varepsilon}\right) = \lim_{\varepsilon\rightarrow 0} \varepsilon ~ {\rm rect}\left(\varepsilon x\right) \end{aligned}$$

$$\boxed{\begin{gathered} \int_{-T/2}^{T/2} e^{-i 2\pi f t} dt = \int {\rm rect}\left({\small \frac{t}{\footnotesize T}}\right) e^{-i 2\pi f t} dt = T {\rm sinc}(Tf)\\ \delta(x) = \lim_{T\rightarrow\infty} \int_{-T/2}^{T/2} e^{-i 2\pi f t} dt = \lim_{T\rightarrow\infty} T {\rm sinc}(Tf) \end{gathered}}$$

$$\begin{gathered} \delta_T(f) \equiv T {\rm sinc}(Tf) = \frac{\sin(\pi Tf)}{\pi f}\\ \delta_T(0) = T, ~~~ \lim_{T\rightarrow\infty}\delta_T(f) = \delta(f) \end{gathered}$$

$$\begin{aligned} \sigma_{\mathbb{S}} &= \int_{-T/2}^{T/2} dt \int_{-T/2}^{T/2} dt' \int df \int df' e^{i 2\pi f (t-t')} e^{-i 2\pi f' (t-t')} N_I(f) N_J(f')\\ &= \int df \int df' N_I(f) N_J(f') \int_{-T/2}^{T/2} dt \int_{-T/2}^{T/2} dt' e^{i2\pi f(t-t')} e^{i2\pi f'(t-t')} \\ &= \int df \int df' ~ N_I(f) N_J(f') ~ \delta^2_T(f-f')\\ &\approx \int df \int df' ~ N_I(f) N_J(f') ~ \delta(f-f') ~ \delta_T(f-f')\\ &= \int df ~ N_I(f) N_J(f) ~ \delta_T(0)\\ &= T \int df ~ N_I(f) N_J(f)\\ \end{aligned}$$

噪声功率谱作为频率的函数，在Hz的尺度上变动，而信号的观测时间$T$在月的量级，对应$\delta_T(f) = T {\rm sinc}(Tf)$中心主峰宽度$\sim \frac{1}{T}$，远小于噪声功率谱变化尺度。在上面的计算中其中的一个$\delta_T(f)$被近似替换为$\delta(f)$。

$$\int_{-T/2}^{T/2} dt \int_{-T/2}^{T/2} dt' e^{i 2\pi f (t-t')} e^{-i 2\pi f' (t-t')} \leftrightarrow \delta^2_T(f-f')???$$

准周期连续源

波源简介

Monochromatic
旋转中子星、双白矮星、EMRI系统(FastEMRIWaveform)

方法概述

对于单频信号，理想情况下，直接傅里叶变换，频谱峰值即对应信号频率，实际情况则复杂很多。假设信号很弱，积累足够信噪比需要1年时间，为保证DFT后信号落入单个频率区间，观测时间内信号频率演化需小于DFT的频率分辨率，即$\dot{f}T_{\rm obs}\lesssim 1/T_{\rm obs} \sim 30\rm nHz$，$\dot{f}\lesssim 1/T^2_{\rm obs} \sim 10^{-15}\rm Hz/s$。即便不考虑探测器运动引入的多普勒调制，波源本身由于引力波辐射导致的频率演化通常也不满足该要求，而且信号越强，演化越明显。显然，只有对多普勒调制及波源演化都进行建模和修正之后，简单的傅里叶分析才可行，而这需要对波源位置及参数有足够精确先验，如电磁观测精确建模的波源。

对周期性信号的搜索，根据观测目标可分为靶向搜索、定向搜索以及全天搜索，分别对应波源已知(方位及频率)、方位已知以及波源未知三种情况。根据计算方法又可分为全相干法和半相干法：全相干法一次性分析所有观测数据，考虑信号频率和相位跨越全部观测时间的演化以及探测器噪声的变化；半相干法则是将数据分段，每段内进行相干分析(同时考虑频率和相位演化)，各数据段之间则忽略相位的演化，非相干叠加。

显然，全相干法具有更好的灵敏度，但通常并不现实，下面从波源的模板和定位两方面估算。对于波源本身，将其相位演化一般性地展开：

$$\phi(t) \approx \phi_0 + 2\pi \left[f_s(t-t_0) + \frac{1}{2}\dot{f}_s(t-t_0)^2 + \frac{1}{6}\ddot{f}_s(t-t_0)^3\right]$$

假设要求在相干时间$T \sim t-t_0$内，模型对相位的预测偏差不超过$\Delta\phi$，则粗略地，对频率及其导数误差要求为：

$$\Delta f_s \sim \frac{1}{2\pi}\frac{\Delta\phi}{T}, ~~~~ \Delta \dot{f}_s \sim \frac{2}{2\pi}\frac{\Delta\phi}{T^2}, ~~~~ \Delta \ddot{f}_s \sim \frac{6}{2\pi}\frac{\Delta\phi}{T^3}$$

这大致对应于模板精度的需求，从而近似到二阶导，所需模板数量将$\sim T^6$，考虑到模板匹配的计算量同样正比于相干时间$T$，总计量算$\sim T^7$。这只是粗略估算，实际中并不会如此极端，低阶效应主导时无需考虑更高阶效应，比如对于二阶导数$\ddot{f}_s$，只有当相干时间足够长时影响才显著。除了波源本身演化，探测器运动引入的多普勒调制是另一个重要因素，且依赖于具体波源方位，不同方位需要使用不同修正，大致地：

$$\Delta\theta \sim \frac{\frac{1}{2}\Delta f_s}{f_s v_{\perp}/c} \sim \frac{1}{2 v_{\perp}/c} \frac{1}{f_s T}, ~~~~ \Delta\Omega \sim (\Delta\theta)^2$$

其中$v_{\perp}$是探测器速度在垂直波源方向上的投影。由此多普勒修正的模板数$\sim f_s^2 T^2$，结合前面波源本身模板的需求，随着相干时间增加，计算量将迅速上升，因此全相干分析只有在电磁观测或半相干法提供了足够精确的参数先验情况下才有可能实现。

实际中通常是基于半相干法的层级搜索(Hierarchical Search)，即先在完整参数空间中使用相对粗略的模板快速搜索，之后在潜在目标区域精细化搜索，以平衡搜索精度和计算量需求，具体地可大致分为搜索、否决和证认三个阶段：

• 搜索：基于半相干法全局参数搜索(Search Pipeline)
  显然半相干分析会降低探测灵敏度，且不同的搜索策略会有所区别，但大致地，对给定的观测时间，假设数据分为$N$段，对信号应力强度的灵敏度损失$\sim N^\frac{1}{4}$。
  • $\mathcal{F}$-统计量：GCT、Wave、
  • 频谱分析：弱相干法PowerFlux/Falcon、互相关法Cross-Correlation
  • 时频谱分析：霍夫变换SkyHough/FrequencyHough
    频谱分析基于频谱构建统计量，时频谱法则更关注频率演化，通过对二值化时频谱进行霍夫变换，定位频率演化曲线。
  • 维特比算法：Viterbi search
  • 机器学习
• 否决：通过后处理(Post-Processing)分析，消除误报
  后处理的核心是参数聚类(Clustering)以及一致性分析(Coincidence)，其中一致性分析可在同一探测器不同时间段或不同探测器之间进行。
  • 一致性否决
  • $\chi^2$卡方否决
  • 窄带噪声否决
  • 零假设否决
• 证认：增加相干时间，对候选波源进行后随(Follow-up)证认
  后随证认本质是用更高精度对候选信号进行分析，基本的策略是增加相干时间，提升参数网格精度，或基于候选区间对目标参数进行采样(非网格化方法)。具体实现可分为单阶段或多阶段(精度依次提升)，数据可复用旧数据或使用新观测数据。
  • 旧数据 v.s. 新数据：证认阶段可复用旧数据或使用新的观测数据。
  • 单阶段 v.s. 多阶段：
  注意，后随认证阶段增加了分析结果的信噪比和可靠性，但无法提升探测灵敏度，灵敏度完全取决于半相干分析阶段的搜索策略。增加相干时间在增加相位一致性要求的同时也提升了频谱的频率分辨率？[两者是一个东西吗？]。

搜索算法

相干分析：外差法、重采样、狄利克雷滤波、时域分析、五矢量法、$\mathcal{F}$-统计量
半相干法：stack-slide、PowerFlux、HoughTransform、stacked $\mathcal{F}$-statistic

$\mathcal{F}$-统计量
服从卡方分布，而非F分布，与通常的F检验不同？

滑动堆栈Stack-slide
Stack：对功率平方求和
霍夫变换Hough Transform：搜索频率演化曲线

时频谱：短时傅里叶变换
“白化”：噪声功率谱估计并非通常的Welch法，而是对周期图进行移动平均，近似于Daniell法。
二值化：阈值、峰值
变换：拉东变换(Radon Transform)、霍夫变换(Hough Transform)

随机背景信号

波源简介

能谱/分数能量密度谱

假设背景是稳态、各向同性、无偏振；高斯
稳态：随机信号的统计性质不随时间演化，具体的，互相关信号只与时间差有关而于具体时间选择无关。
各向同性：信号源自分布相对均匀的遥远天体或早期宇宙过程，与本地物质不存在关联。
无偏振：各方向偏振强度相等，没有明确的“主”偏振方向。
圆偏振涉及的是两种偏振模式的相位差，这里所谓的“无偏振”是指偏振角随机分布，或者说只是无线偏振，两者并不矛盾？

高斯性：多个探测器信号$s(t)$的联合概率密度函数为多元正态分布。对于高斯稳态随机过程，假设均值为零，统计性质完全由其功率谱(自相关函数/二阶距)描述。
对于稳定高斯随机过程，其性质完全由均值(一阶矩)和互相关/自相关(二阶矩)决定。也意味着N点关联可表示为两点关联乘积之和(均值通常设为零)。
其合理性源自大数定律，
根据大数定律，任意足够多的独立随机变量之和都渐近服从正态分布，而与随机变量本身分布无关。因此对于大量独立天体源叠加而成的随机背景，高斯假设是合理的，对于相互关联的早期宇宙源？或少数天体源的叠加则需要谨慎。

在长波近似下，地基观测
在短波近似下，就是Hellings-Downs(HD)曲线。
PTA中观测量通常为多普勒红移(对应频率变化)，而非应变响应(对应相位变化)。
此时$z(t) = $

对于单个探测器而言，随机背景与探测器噪声原则是没法区分的，因此对于随机背景的探测通常需要多个探测器互相关，或者对于仪器噪声有准确的测量或模型描述。

随机背景分析

• 多探测器互相关：多个探测器进行联合观测
• 探测器噪声先验：通过零通道获取噪声信息

$$h_{ab}(t, \vec{r}) = \int df \int \tilde{h}_{ab}(f, \hat{k})e^{-i 2\pi f (t-\hat{k}\cdot\frac{\vec{r}}{c})}$$

$\tilde{h}_{ab}(f, \hat{k}) = \sum_P \tilde{h}_P(f, \hat{k})e_{ab}^P(\hat{k})$
$\sum_P e_{ab}^P e_{ab}^P = 4$

$$\overline{\left\langle \tilde{h}_P(f,\hat{k})\tilde{h}^*_{P'}(f',\hat{k}') \right\rangle} = \frac{1}{4\pi} \delta(f-f')~\delta^2(\hat{k},\hat{k}')~\delta_{PP'}~S_h(f)$$

$\delta^2(\hat{k},\hat{k}') = \delta(\cos\theta-\cos\theta')\delta(\varphi-\varphi')$

$$\int d^2\hat{k}~d^2\hat{k}' \overline{\left\langle \tilde{h}_P(f,\hat{k})\tilde{h}^*_{P'}(f',\hat{k}') \right\rangle} = \delta(f-f')~\delta_{PP'}~S_h(f)$$

$d^2\hat{k}' = d\cos\theta d\varphi$

$$\frac{1}{2}\int d^2\hat{k}~d^2\hat{k}'\left[\overline{\left\langle \tilde{h}_+(f,\hat{k})\tilde{h}^*_+(f',\hat{k}') \right\rangle} + \overline{\left\langle \tilde{h}_\times(f,\hat{k})\tilde{h}^*_\times(f',\hat{k}') \right\rangle}\right] = \delta(f-f')~ S_h(f)$$

$$\left\langle \tilde{h}_P(f)\tilde{h}^*_{P'}(f') \right\rangle = \int d^2\hat{k}~d^2\hat{k}'\overline{\left\langle \tilde{h}_P(f,\hat{k})\tilde{h}^*_{P'}(f',\hat{k}') \right\rangle} = \delta(f-f')~\delta_{PP'} S_h(f)$$

$$\left\langle \tilde{h}_+(f)\tilde{h}^*_+(f') \right\rangle = \left\langle \tilde{h}_\times(f)\tilde{h}^*_\times(f') \right\rangle = \delta(f-f')~ S_h(f)$$

$$\frac{1}{2}\left(\left\langle \tilde{h}_+(f)\tilde{h}^*_+(f') \right\rangle + \left\langle \tilde{h}_\times(f)\tilde{h}^*_\times(f') \right\rangle\right) = \delta(f-f')~ S_h(f)$$

$S_h(f)$是偏振平均的功率谱，如果便是为单边谱，则$1/2$抵消了，成了求和？如果考虑更多偏振模式？$\Omega_GW$表示随机引力波能量占比，应该是对偏振求和？

$$\boxed{ \left\langle h_{ab}(t)h^{ab}(t)\right\rangle = 4 \int df S_h(f) }???$$

$$\left\langle h_+^2(t) + h_\times^2(t)\right\rangle = 2 \int df S_h(f) ~~~~~~~ h_{ab}h^{ab}=2(h_+^2 + h_\times^2)$$

$$S_h(f) = f^2 h^2_c(f), \frac{\left\langle h_+^2(t) + h_\times^2(t)\right\rangle}{2} = \int df f h_c(f)$$

$$S_h(f) = \frac{3H_0^2}{4\pi^2 f^3} \Omega_{\rm \tiny GW}(f)$$

$S_h(f)$似乎有两种定义方式，一种是各种偏振模式的平均，一种是各种偏振模式的求和。对$+, \times$模式就相差两倍，因此会看到有的文章使用$S_h(f) = \frac{3H_0^2}{2\pi^2 f^3} \Omega_{\rm \tiny GW}(f)$。跟单边功率谱或双边功率谱似乎是没有关系的，因为如果这里$S_h$是单边的，双边的应该是1/2而非2倍。如果考虑更多偏振模式，这里应该就不再是2倍的关系？

$$S_h(f) = \left\langle \tilde{h}_P(f, \hat{k})\tilde{h}^*_P(f, \hat{k}) \right\rangle, ~~~ S_h(f) = 4\pi \overline{\left\langle \tilde{h}_P(f)\tilde{h}^*_P(f) \right\rangle}$$

$$\left\langle \tilde{n}(f)\tilde{n}^*(f')\right\rangle_T = N(f)\delta(f-f')$$

注意$\left\langle \tilde{n}(f) \tilde{n}^*(f) \right\rangle = \delta(0) N(f) \rightarrow \infty$ 趋向无穷大，而非等于$N(f)$！

$$\boxed{\left\langle n(t) n(t') \right\rangle \leftrightarrow \left\langle \tilde{n}(f) \tilde{n}^*(f') \right\rangle}$$

$$\begin{aligned} \left\langle n(t) n(t') \right\rangle &= \left\langle \int n(f) e^{i2\pi f t} df \int n^*(f') e^{-i2\pi f' t'} df' \right\rangle\\ &= \int df \int df'\left\langle n(f) n^*(f') \right\rangle e^{i2\pi f t} e^{-i2\pi f' t'}\\ &= \int df \int df'N(f)\delta(f-f') e^{i2\pi f t} e^{-i2\pi f' t'}\\ &= \int df N(f) e^{i2\pi f (t-t')}\\ \end{aligned}$$