引力波天文II：探测响应与灵敏度

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

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

本文介绍通过激光干涉探测引力信号的基本原理，具体关注探测器对信号的响应函数以及灵敏度曲线。

基本探测原理

长波近似结论

TT transverse-traceless.

$$h(t, \vec{r}) = h_{\vec{r}=0}(t-\hat{k}\cdot \vec{r}/c) = \int \tilde{h}(f) e^{i2 \pi f(t-\hat{k}\cdot\frac{\vec{r}}{c})}$$

当$L \ll \lambda/2\pi, f \ll \frac{c}{2 \pi L}$，假设探测器位于坐标原点，有$2 \pi f \hat{k}\cdot\frac{\vec{r}}{c} \ll 1$，从而$h(t, \vec{r})$可近似简单替换为$h(t)$，忽略同一时刻探测器不同位置的信号差异，即引力波的传播效应，被称为长波近似。但光子传播过程，不同时刻的引力波信号差异仍需要考虑。

下面考虑激光沿干涉臂传播过程中受到的影响，为便于推导这里先关注单个频率成分：

$$\Delta t = \frac{1}{2c}\int_{s=0}^{L} h(s) d s= \frac{1}{2c}\int_{s=0}^{L} \tilde{h}(f) e^{i\omega t} ds$$

其中，$\vec{r}$ 和 $t$对参量$s$ 的依赖关系由激光传播决定：

$$t(s) = t|_{\vec{r}_1} + s/c = t|_{\vec{r}_2} - L/c + s/c$$

这里$\vec{r_1}, \vec{r_2}$分别为激光发射端与接收端，$\hat{x}$由激光发射端指向接收端，$t|_{\vec{r}_1}, t|_{\vec{r}_2}$ 分别为激光发射和接收时刻。注意，上式中只出现了$\hat{x}$而没有$\hat{k}$，因为上述关系完全由激光传播决定，而与引力波无关，激光是沿$\hat{x}$传播。

$$\begin{aligned} \Delta t &= \frac{1}{2c} ~\tilde{h}(f) \int_{s=0}^{L} e^{i\omega (t|_{\vec{r}_1} + \frac{s}{c})} ds\\ &= \frac{1}{2c} ~\tilde{h}(f)~ e^{i\omega t|_{\vec{r}_1}} \int_{s=0}^{L} e^{i\omega \frac{s}{c}} ds\\ &= \frac{1}{2c} ~\tilde{h}(f)~ e^{i\omega t|_{\vec{r}_1}} {\small \frac{c}{i\omega} } e^{i\omega \frac{s}{c}}\Big|_0^L\\ &= \frac{1}{2} ~\tilde{h}(f)~ e^{i\omega t|_{\vec{r}_1}} \small \frac{1}{i\omega} \left[e^{i\omega \frac{L}{c}}-1\right]\\ &= \frac{1}{2} ~\tilde{h}(f)~ e^{i\omega t|_{\vec{r}_1} } e^{i\omega \frac{2L}{c}} ~ \small \frac{2\sin\left[\omega \frac{L}{2c}\right] }{\omega}\\ &= \tilde{h}(f)~ e^{i\omega (t|_{\vec{r}_1} + \frac{L}{2c}) } ~{\small \frac{L}{2c}} ~ {\rm sinc}\left[\omega {\textstyle \frac{L}{2c}}\right]\\ \end{aligned}$$

这里$t|_{\vec{r}_1} + \frac{L}{2c}$ 对应激光传播到干涉臂中间的时间。因此对于单频引力波，单向传播的激光受到的总影响正比于干涉臂中间时的应变强度乘以${\rm sinc}$项，高频信号受到压制。${\rm sinc}$项傅里叶逆变换为$\rm rect$函数，乘积变卷积。

$$\begin{aligned} \Delta t (t|_{\vec{r}_1}) &= \frac{1}{2} h(t|_{\vec{r}_1} + {\textstyle \frac{L}{2c}}) * {\rm rect}\small \left(\frac{c}{L} t|_{\vec{r}_1} \right) \end{aligned}$$

引力波信号与矩形窗卷积，相当于信号的滑动平均，窗口宽度为$\frac{L}{c}$，测量效应为引力波信号的平均。

代入$t|_{\vec{r}_1} = t|_{\vec{r}_2} - \frac{L}{c}$，用接收端的时间和位矢表示有：

$$\begin{aligned} \Delta t &= \tilde{h}(f) e^{i\omega (t|_{\vec{r}_2} - \frac{L}{2c})} ~ {\small \frac{L}{2c}} ~ {\rm sinc}\left[\omega {\textstyle \frac{L}{2c}}\right] \end{aligned}$$

除了时间(相位)变化，对引力波效应的另一个常用描述是分数多普勒频移$\frac{\Delta \nu}{\nu} = \frac{\nu_{\rm recv} - \nu_{\rm emit}}{\nu_{\rm emit}}$。PTA相关文献中通常记为符号$z$，而在空间引力波探测中更常用$y$表示。根据相位关系：

$$\phi_{\rm recv} = \phi_{\rm emit} + 2 \pi \nu_{\rm emit}\frac{L}{c} + 2 \pi \nu_{\rm emit}\Delta t$$

求导可得$\frac{\Delta \nu}{\nu} = \dot{\Delta t} = i\omega \Delta t$，从而：

$$\begin{aligned} \Delta t &= \frac{1}{2}~\tilde{h}(f)~ e^{i\omega t|_{\vec{r}_1}} \small \frac{1}{i\omega } \left[e^{i\omega\frac{L}{c}}-1\right]\\ \frac{\Delta \nu}{\nu} & = \frac{1}{2}~\tilde{h}(f)~ e^{i\omega t|_{\vec{r}_1}} \left[e^{i\omega\frac{L}{c}}-1\right]\\ & = \frac{1}{2} \left[\tilde{h}(f)~e^{i\omega t|_{\vec{r}_2} } - \tilde{h}(f)~e^{i\omega t|_{\vec{r}_1}} \right] \end{aligned}$$

傅里叶逆变换

$$\begin{aligned} y(t) &= \frac{1}{2} \left[ h(t|_{\vec{r}_2}) - h(t|_{\vec{r}_1}) \right] \end{aligned}$$

这里$\vec{r_1}, \vec{r_2}$分别对应激光发射端与接收端，$t|_{\vec{r}_1}, t|_{\vec{r}_2}$ 则对应激光的发射和接受时间。即，引力波导致的多普勒频移仅由接收时刻(接收端)与发射时刻(发射端)的应变强度之差决定！

除了这里通过对光子路径积分获得距离变化，对应相位改变，求导可得到频率变化。还可以直接推导频率影响，之后积分得到相位改变。
高阶速度项的影响可忽略(Cornish & Rubbo 2003)

最后，对于单臂往返，发射与接收端重合，将最终接收时间$t_{\rm recv}$简记为$t$，考虑往返的总效应，有：

$$\begin{aligned} \Delta t &= \frac{1}{2}~\tilde{h}(f)~ e^{i\omega (t - \frac{L}{c}) } \small \frac{1}{i\omega} \left[e^{i\omega \frac{L}{c}} - 1\right] + \frac{1}{2}~\tilde{h}(f)~ e^{i\omega (t - \frac{2L}{c})} \small \frac{1}{i\omega } \left[e^{i\omega \frac{L}{c}}-1\right]\\ &= \frac{1}{2}~\tilde{h}(f)~ e^{i\omega t } \small \frac{1}{i\omega} \left[ 1 - e^{-i\omega \frac{2L}{c}} \right]\\ &= \frac{1}{2}~\tilde{h}(f)~ e^{i\omega (t -\frac{L}{c}) } {\textstyle \frac{2L}{c}} {\rm sinc}\left[\omega {\textstyle \frac{L}{c}}\right] \end{aligned}$$

$$\begin{aligned} \tilde{y}(f) &= \frac{1}{2}~\tilde{h}(f)~ \left[1- e^{-i\omega\frac{2L}{c}} \right] e^{i\omega t}\\ &= \frac{1}{2} ~\tilde{h}(f)~ e^{i\omega t} ~ 2i~ \sin\left({\textstyle \omega\frac{L}{c} }\right)~e^{-i\omega\frac{L}{c}} \\ \end{aligned}$$

$$\begin{aligned} y(t) = \frac{1}{2} \left[ h(t) - h(t-{\textstyle \frac{2L}{c} }) \right]\\ \end{aligned}$$

$$\boxed{\begin{aligned} \Delta t &= \frac{1}{2}~\tilde{h}(f)~ e^{i\omega (t - \frac{L}{c})} {\textstyle \frac{2L}{c} } {\rm sinc}\left[\omega {\textstyle \frac{L}{c}}\right] \end{aligned}}$$

如果将上述时间变化量转换为对应的无量纲应变响应，有：

$$\begin{aligned} s(f, t, \hat{x}, \hat{k}, \vec{r}_1) &\equiv {\small \frac{\Delta L}{L}} = \frac{c \Delta t/2}{L} = \frac{1}{2}~\tilde{h}(f)~ e^{i 2\pi f (t - \frac{L}{c})} {\rm sinc}\left[\omega {\textstyle \frac{L}{c}}\right] \end{aligned}$$

注意，这里$\Delta L = c \Delta t/2$，因为$\Delta t$为往返时长差，因此光程差对应于2倍的臂长差。
$s(f,t)$就是激光臂对单个频率成分 $\tilde{h}(f)e^{i\omega t}$ 的响应。可以看到响应是线性的，即$s \propto \tilde{h}(f)e^{i\omega t}$，因此总的引力波响应就是所有频率分量贡献的线性叠加，即对频率积分。用$\mathscr{D}$表示探测器响应：

$$\begin{aligned} s(t) &= \mathscr{D}\left(h(t)\right) = \mathscr{D}\left(\int \tilde{h}(f)e^{i 2\pi f t} df\right)\\ &= \int \mathscr{D}\left(\tilde{h}(f)e^{i 2\pi f t}\right) df = \int s(f, t) df\\ &= \int df ~ \frac{1}{2}~\tilde{h}(f)~e^{i 2\pi f t} ~\mathcal{T}(f, \hat{k} \cdot \hat{x})\\ &= \int \left[ \frac{1}{2}~\tilde{h}(f)~~\mathcal{T}(f, \hat{k} \cdot \hat{x}) \cdot \hat{x})\right] ~e^{i 2\pi f t} df \end{aligned}$$

$$\begin{aligned} \tilde{s}(f) &= \mathcal{F}\left\{s(t)\right\} = \frac{1}{2}~\tilde{h}(f)~ e^{-i \omega\frac{L}{c}} {\rm sinc}\left[\omega {\textstyle \frac{L}{c}}\right] \end{aligned}$$

至此还只是单臂往返，其中$h(t)$是引力波信号在探测臂方向的投影，依赖其具体指向，张量形式为$\hat{x}^i\hat{x}^j h_{ij}(t)$。从而对于双臂干涉

$$\boxed{\begin{aligned} \tilde{y}(f) &= \frac{1}{2} \left[1- e^{-i\omega\frac{2L}{c}} \right] \left[\hat{x}^i \hat{x}^j - \hat{y}^i \hat{y}^j \right] \tilde{h}_{ij}(f)\\ &= \frac{1}{2} ~2i~ e^{-i\omega\frac{L}{c}}~\sin\left({\textstyle \omega\frac{L}{c} }\right) \left[\hat{x}^i \hat{x}^j - \hat{y}^i \hat{y}^j \right] \tilde{h}_{ij}(f) \end{aligned}}$$

$$\boxed{\begin{aligned} \tilde{s}(f) &= \frac{1}{2} ~{\small \frac{1}{i\omega}~ \frac{1}{2L/c}}~ \left[1- e^{-i\omega\frac{2L}{c}} \right] \left[\hat{x}^i \hat{x}^j - \hat{y}^i \hat{y}^j \right] \tilde{h}_{ij}(f)\\ &= \frac{1}{2} ~ e^{-i \omega\frac{L}{c}}~ {\rm sinc}\left[\omega {\textstyle \frac{L}{c}}\right] \left[\hat{x}^i \hat{x}^j - \hat{y}^i \hat{y}^j \right] \tilde{h}_{ij}(f) \end{aligned}}$$

对比两者，区别在于$\sin$变为$\rm sinc$，没有虚数单位$\rm i$，同时因为$\Delta L/2L$消去了原本的2倍。
如果进一步做长波近似，取$f \rightarrow 0$，则可得到更常见的长波近似结果：

$$\boxed{\begin{aligned} \tilde{y}(f) &= {\small i \frac{\omega L}{c} } \left[\hat{x}^i \hat{x}^j - \hat{y}^i \hat{y}^j \right] \tilde{h}_{ij}(f)\\ \tilde{s}(f) &= \frac{1}{2} ~\left[\hat{x}^i \hat{x}^j - \hat{y}^i \hat{y}^j \right] \tilde{h}_{ij}(f) \end{aligned}}$$

---

引力波对时空的扰动会影响光束经过测底线上两点的时间间隔，进而反映在相位测量上，这是现代引力波测量的基本原理。在平直时空中任取两点，连线方向设为$x$轴，对应空间坐标分别为$(0,0,0)$和$(L,0,0)$。只考虑$x$方向($dy=dz=0$)，假设引力波引入的时空微扰在$x$方向的投影为$h(t)$，所对应的线元可表示为：

$$ds^2 = g_{\mu\nu} dx^\mu dx^\nu = -c^2dt^2 + [1+h(t)]dx^2$$

$$ds^2 = g_{\mu\nu} dx^\mu dx^\nu = -c^2dt^2 + [1+h_+(t)]dx^2 + [1+h_\times(t)]dx^2 + dz^2$$

这里暂时只考虑$x$方向的引力波观测效应，并暂时将$h_+(t)$简记为$h(t)$。光线沿零测地线(null geodesic)移动$ds^2 = 0$，从而对于沿$x$方向传播的光$c dt = \sqrt{1 + h_+(t)}dx$。这就是用光进行(干涉)测量的基本原理，即对光而言时空距离的改变会反映在坐标时间中，并最终通过相位观测记录。由发射点到接收点积分：

$$c(t_{\rm recv} -t_{\rm emit}) = \int_{t_{\rm emit}}^{t_{\rm recv}} c dt = \int_0^L \sqrt{1 + h[t(x)]}dx \approx \int_0^L \left\{1 +\frac{1}{2}h[t(x)]\right\}dx$$

这里利用$h(t) \ll 1$做了一阶近似展开。最终光信号到达接收点时间为：

$$t_{\rm recv} \approx t_{\rm emit} + \frac{L}{c} + \frac{1}{2c} \int_0^L h[t(x)] dx$$

最后一项为引力波引入的修正，但积分反过来又依赖自身$t(x)$，通常需迭代。一阶近似取$t(x) = t_{\rm emit} + \frac{x}{c}$，接收端对应时间为：

$$t_{\rm recv} \approx t_{\rm emit} + \frac{L}{c} + \frac{1}{2c} \int_0^L h\left(t_{\rm emit} + \frac{x}{c}\right) dx = t_{\rm emit} + \frac{L}{c} + \frac{1}{2} \int_{t_{\rm emit}}^{t_{\rm emit} + \frac{L}{c}} h(t) dt$$

最后一项为引力波引入的一阶修正，进一步迭代所引入的高阶修正量为$o(h^2)$量级，这里暂时忽略。最终，引力波效应为：

$$\boxed{\begin{aligned} \Delta t &= t - t_{\rm emit} - \frac{L}{c} = \frac{1}{2} \int_{t-\frac{L}{c}}^{t} h(\tau) d\tau\\ \Delta L &= c\Delta t =\frac{c}{2} \int_{t-\frac{L}{c}}^{t} h(\tau) d\tau \end{aligned}}$$

考虑到相位测量是在接收端进行，这里调整了积分限，并将$t_{\rm recv}$简记为$t$，并不影响结果。最终相位测量：

$$\varphi_{\rm recv} \approx \varphi_{\rm emit} + 2\pi\nu_0\frac{L}{c} + 2\pi\nu_0\Delta t$$

其中$2\pi \nu_0 \Delta t$为引力波引入的相位偏移，$\nu_0$为激光发射频率。通过激光干涉技术，可以对该相位偏移做出精确测量，这是现代引力波测量的核心。地基引力波探测处于长波极限，$L \ll \lambda$，波长覆盖整个探测器，引力波振幅在两点间变化不大，积分内的$h(t)$可近似为常数，从而$\Delta L = \frac{1}{2} L h(t) \propto h(t)$。

对于单频信号，取$h(t) = h_0 e^{-i\omega t}$，带入前面的公式有：

$$\begin{aligned} \Delta t &= \frac{1}{2} \int_{t-\frac{L}{c}}^{t} h_0 e^{-i\omega \tau} d\tau = \frac{h_0}{2} \frac{1}{-i\omega} e^{-i\omega \tau}\Big|_{t-\frac{L}{c}}^{t} \\ &= \frac{h_0}{2} \frac{1}{-i\omega} e^{-i\omega t} \left[1 - e^{i\omega \frac{L}{c}}\right]\\ & = h_0 e^{-i\omega t} e^{i\omega \frac{L}{2c}} \frac{\sin(\omega \frac{L}{2c})}{\omega} \end{aligned}$$

简单整理有：

$$\boxed{\Delta t = {\small \frac{L}{2c}} ~ h(t-{\textstyle \frac{L}{2c}}) ~ {\rm sinc}(\omega {\textstyle \frac{L}{2c}})}$$

这里${\rm sinc}(x) = \frac{\sin x}{x}$，在长波极限下$L \ll \lambda, \omega \frac{L}{2c} \ll \pi$，${\rm sinc}(\omega \frac{L}{2c}) \sim 1$。而随着引力波频率增加，$\omega \frac{L}{2c} \gtrsim \pi$，引力波信号将比显著抑制。注意，这里测量信号依赖$h(t- {\textstyle \frac{L}{2c}})$，如果是沿$x$方向往返，$L \rightarrow 2L$，则对应于在远端反射时刻的信号强度。

另一个方面，对前面公式两边求导有：

$$2\pi\nu(t) \approx 2\pi\nu_0 + 2\pi\nu_0\dot{\Delta t} = 2\pi\nu_0 + 2\pi\nu_0 \frac{1}{2} \Big[h(t) - h(t - {\textstyle \frac{L}{c}})\Big]$$

这里利用了变限积分的导数，假设激光发射频率稳定为$\nu_0$，并定义瞬时频率$\nu(t) = \frac{\dot{\varphi}(t)}{2\pi}$。简单整理有：

$$\boxed{\frac{\nu(t)-\nu_0}{\nu_0} \approx \frac{1}{2} \Big[h(t) - h(t - {\textstyle \frac{L}{c}})\Big]}$$

即接收到的激光瞬时频率相对发射端基准频率的相对偏移正比于引力波在两点的振幅差。直观理解，引力波的波峰或波谷的到达时间受到引力波调制，最终的接收频率也将随之发生变化。长波极限下，$h(t)$在光传播过程中近似常数，$\frac{\nu(t)-\nu_0}{\nu_0} \approx 0$，因此地基引力波探测无需考虑频率偏移，但对于空间引力波探测，激光频率的相对偏移是经常讨论的物理量，后面会具体展开。值得注意的，这里激光的频率偏移仅依赖引力波在两端点时刻的引力波振幅，而与激光传播过程中的引力波演化无关(不限制$h(t)$的形式)！

激光干涉测量

https://en.wikipedia.org/wiki/Interferometry
https://en.wikipedia.org/wiki/Sagnac_effect
https://en.wikipedia.org/wiki/Common-path_interferometer#Zero-area_Sagnac
后面会提到，激光的频率噪声远高于引力波测量所需精度，因此上述单臂单向响应在实际的激光测量中并不现实，需要等臂干涉抑制激光频率噪声。而反过来说，如果计时精度足够，则无需干涉，这正是脉冲星计时用于引力波测量的原理。注意，脉冲星观测中实现计时所依靠的并非射电辐射本身的频率精度，而是射电脉冲间隔(脉冲到达时间)精度。因此，对脉冲星计时测量而言，多普勒频移是脉冲间隔的频移，而非电磁波本身频移（虽然两者受引力影响相同？？），后者对应多普勒追踪(Doppler tracking)测量，属于另一种引力波测量方法。最后，脉冲计时测量随机背景时，需要在不同脉冲星测量间进行互相关，注意区分相关分析与干涉测量，干涉仪测量背景引力波时同样需要不同探测器间进行相关分析。

对于引力波测量而言，两种常见的激光干涉主要有迈克尔逊干涉(Michelson)和萨格纳克干涉(Sagnac)。如下图所示，两者的主要区别在于光线的传播路径：前者经过分光的两束光线走不同路径进行干涉，而后者则是经过相同路径进行干涉。

迈克尔逊干涉更为知名，当前运行的aLIGO, Adv-Virgo等都属于迈克尔逊干涉。
早在狭义相对论提出(1905)之前，Michelson为了检验以太理论，于1881年首次构想并实施了激光干涉测量，并在之后通过改进的迈克尔逊-莫雷实验(1887)否定了“以太”理论，而被载入史册，直接推动了洛伦兹变换的提出及狭义相对论的发展。迈克尔逊1886年就提出过激光干涉的萨格纳克构型，随后又被萨格纳克重新提出(Sagnac 1913)，而广义相对论的提出则要到更晚的1915年。

萨格纳克效应用于激光或光纤陀螺仪

迈克尔逊干涉实验是非常著名的相对论检验实验，但萨格纳克效应

而将萨克纳克干涉用于引力波测量的想法很早就有人提出(R. W. P. Drever 1983, R.Weiss 1987, Sun et al. 1996)，
为此人们提出了零面积的萨格纳克干涉与LIGO的构型非常接近，只是在分光镜附近增加了一个反射镜，使得两束光经相同路径后再干涉。
在讨论下一代引力波探测器时，除了延续迈克尔逊干涉思路的CE、ET，基于萨格纳克干涉的设想同样被不断提及。

萨格纳（Sagnac）原理：也称萨氏效应（相位差正比于输入角速率）。该原理适用于光纤角速率陀螺；激光角速率陀螺等。Sagnac法国物理学家（1869年～1926年），居里夫妇的朋友。1913年发明萨氏效应。

As noted above, Sagnac interferometers are, to first order, insensitive to any static or low-frequency displacement of their optical components, nor are the fringes affected by minor frequency variation in the lasers or birefringence.
This variant of the Sagnac interferometer is hence insensitive to rotation or low frequency drift of its optical components,

insensitive to rotations while remaining sensitive to time-dependent displacement in the two arms.
Sun et al. have investigated the possibility of detecting a broad-band GW detector with a zeroarea Sagnac interferometer.

外差干涉测量

两个频率接近的电磁波信号相干叠加，其电场强度可表示为(这里假设两个信号偏振方向相同)：

$$E_s\cos(\omega_s t + \phi) + E_{\rm LO}\cos(\omega_{\rm LO} t)$$

在探测器进行探测时，光电感应器的响应信号与光强而非场强成正比(非线性响应)，而光强

$$\begin{aligned} I \propto ~&~ [E_s\cos(\omega_s t + \phi) + E_{\rm LO}\cos(\omega_{\rm LO} t)]^2\\ =~&~ E_s^2\cos^2(\omega_s t + \phi) + E_{\rm LO}^2\cos^2(\omega_{\rm LO} t) + 2E_sE_{\rm LO}\cos(\omega_s t + \phi)\cos(\omega_{\rm LO} t)\\ =~&~ \frac{1}{2}E_s^2[1+\cos 2(\omega_s t + \phi)] + \frac{1}{2}E_{\rm LO}^2[1+\cos(2\omega_{\rm LO} t)]\\ ~&~ + E_sE_{\rm LO}\big\{cos[(\omega_s-\omega_{\rm OL}) t + \phi] + \cos[(\omega_s+\omega_{\rm OL}) t + \phi]\big\} \end{aligned}$$

探测器只能探测到信号在响应时间内的积分(平均强度)，考虑到光电感应器响应时间≫信号本身周期，上式中的高频项积分响应为0，从而最终的响应

$$\propto \frac{1}{2}(E_s^2+E_{\rm LO}^2) + E_sE_{\rm LO}\cos[(\omega_s-\omega_{\rm OL}) t + \phi]$$

也可以解释为，通过低通滤波移除了所有高频信号。在目标信号与参考输入频率差值足够小时，差频(拍频)项$\omega_s-\omega_{\rm OL}$将得以保留。而当$\omega_s$与$\omega_{\rm LO}$相等，系统响应

$$\propto \frac{1}{2}(E_s^2+E_{\rm LO}^2) + E_sE_{\rm LO}\cos \phi$$

外差干涉与零差干涉

• 探测器实际响应的是$E_sE_{\rm LO}$项，相比信号本身($E_s^2$)，通过增强$E_{\rm LO}$可实现弱信号的增益放大。
• 相位差值信息得到了完美保留，可以实现相位的线性响应
• 拍频降低频率，可以实现相位差值或频率偏移的精确测量

拍频越低，测量精度越高，但为什么不能太低？
外差干涉需要细致的频段规划(frequency plan)，同时通过自适应技术保持六路拍频尽可能低，但又不为零(或跨越零)。而为了确保卫星漂移的多普勒效应是频率偏移的唯一来源，所有激光的频率需要相互锁定，并最终与稳定的参考基准保持锁定。

同时由于卫星运动的多普勒效应，在MHz量级，为了拍频不跨越零，

激光功率不稳定性在-5Hz~5Hz内占据主导，因此拍频信号需要避开该频段，即频率下限为5Hz。

下限为6MHz，
通过合适的锁频算法，拍频够可以保持在7~23MHz

考虑到为本身多普勒效应在MHz量级，为保证拍频不跨越零，也需要至少为MHz量级，在5~25MHz之间波动

仪器响应函数

$$\Delta t =$$

正交性

$$\int \frac{d\hat{k}}{4\pi} F_+F_\times = 0$$

平均响应

$$\int \frac{d\hat{k}}{4\pi} F_+ \neq \int \frac{d\hat{k}}{4\pi} F_\times; ~~~~ \int \frac{d\psi}{2\pi} F_+ = \int \frac{d\psi}{2\pi} F_\times$$

$$\int \frac{d\psi}{2\pi} \int \frac{d\hat{k}}{4\pi} F_+ = \int \frac{d\psi}{2\pi} \int \frac{d\hat{k}}{4\pi} F_\times$$

对给定的引力波信号$h_{ij}$，时空伸缩模式是确定的，有确定的特殊标架方向。
对于单个探测器而言，$\psi$是冗余，改变张量极化基矢，最终响应旋转不变，也就是无法测偏振？

天线模式函数

探测器坐标系、日心坐标系、潮汐效应、相对论效应
需要对轨道进行数值积分，以达到所需精度

探测器对于信号的响应依赖空间方位，被称为天线响应模式(antenna response pattern)或简单的天线模式(antenna pattern)，也可称为天线的指向性函数(antenna directivity)或天线波束(antenna beam pattern)，后者在射电领域更常见。对于引力波探测而言，天线模式具体依赖天线构型和引力波的方位及偏振。

广义相对论框架下，采用横向无迹(TT)规范，引力波为横波且只有两种偏振模式，在垂直引力波传播平面内

$$\begin{aligned} \bm{h} = \begin{pmatrix} h_{+} & h_{\times} \\ h_{\times} & -h_{+} \end{pmatrix} = h_{+} \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix} + h_{\times}\begin{pmatrix} 0 & 1\\ 1 & 0 \end{pmatrix} = h_{+} \bm{e}^{+} + h_{\times}\bm{e}^{\times} \end{aligned}$$

其中$\bm{e}^{+}, \bm{e}^{\times}$代表引力波的两种张量偏振模式，$h_{+}, h_{\times}$为对应的强度。在广义相对论之外，引力波还可以有横向的呼吸模式、沿传播方向的纵向模式(纵波)以及横纵交叉的矢量模式。直角坐标系下，引力波的张量模式可通过矢量外积(张量积)描述：

$$\begin{aligned} \bm{e}^{+} &= \hat{u} \otimes \hat{u} - \hat{v}\otimes \hat{v}\\ \bm{e}^{\times} &= \hat{u}\otimes \hat{v} + \hat{v}\otimes \hat{u}\\ \bm{e}^{\circ} &= \hat{u} \otimes \hat{u} + \hat{v}\otimes \hat{v}\\ \bm{e}^{l} &= \hat{k}\otimes \hat{k} \\ \bm{e}^{x} &= \hat{u} \otimes \hat{k} - \hat{k}\otimes \hat{u}\\ \bm{e}^{y} &= \hat{v}\otimes \hat{k} + \hat{k}\otimes \hat{v} \end{aligned}$$

这里$\hat{u}, \hat{v}, \hat{k}$为直角坐标系基矢，其中$\hat{k}$为引力波传播方向，$\hat{u}, \hat{v}$沿$+$模式伸缩方向，$\times$模式相对旋转45°，此时引力波偏振度为0。下面暂不考虑张量模式之外的偏振。
在垂直传播方向的平面内取$\hat{p}, \hat{q}$，相对$\hat{u}, \hat{v}$顺时针旋转$\psi$角，或者说引力波$+$模式主轴相对参考系$\hat{p}, \hat{q}$逆时针旋转$\psi$角则：

$$\begin{aligned} \bm{h} &= \bar{h}_{+} \bm{\epsilon}^{+} + \bar{h}_{\times}\bm{\epsilon}^{\times}\\ \bm{\epsilon}^{+} &= \hat{p}\otimes \hat{p} - \hat{q}\otimes \hat{q}\\ \bm{\epsilon}^{\times} &= \hat{p}\otimes \hat{q} + \hat{q}\otimes \hat{p} \end{aligned}$$

$$\begin{aligned} \begin{pmatrix} \hat{p} \\ \hat{q} \end{pmatrix} = \mathcal{R}_{\psi} \begin{pmatrix} \hat{u} \\ \hat{v} \end{pmatrix} ; ~~~~~ \mathcal{R}_{\psi} = \begin{pmatrix} \cos \psi & -\sin \psi \\ \sin \psi & \cos \psi \end{pmatrix} \end{aligned}$$

$$\begin{aligned} \begin{pmatrix} \hat{p} \\ \hat{q} \end{pmatrix}\otimes \begin{pmatrix} \hat{p} \\ \hat{q} \end{pmatrix} = \left[\mathcal{R}_{\psi} \begin{pmatrix} \hat{u} \\ \hat{v} \end{pmatrix}\right] \otimes \left[\mathcal{R}_{\psi} \begin{pmatrix} \hat{u} \\ \hat{v} \end{pmatrix}\right] = \mathcal{R}_{\psi} \left[\begin{pmatrix} \hat{u} \\ \hat{v} \end{pmatrix} \otimes \begin{pmatrix} \hat{u} \\ \hat{v} \end{pmatrix}\right] \mathcal{R}^T_{\psi} \end{aligned}$$

$$\begin{aligned} \mathcal{R}_{\psi} \begin{pmatrix} h_+ & h_\times\\ h_\times & - h_+ \end{pmatrix} \mathcal{R}^T_{\psi} = \mathcal{R}_{\psi}\mathcal{R}_{\psi} \begin{pmatrix} h_+ & h_\times\\ h_\times & - h_+ \end{pmatrix} = \mathcal{R}_{2\psi} \begin{pmatrix} h_+ & h_\times\\ h_\times & - h_+ \end{pmatrix} \end{aligned}$$

这个对普通矩阵不成立，是这种特殊情况下的结果？

$$??\begin{aligned} \begin{pmatrix} \bar{h}_{+} \\ \bar{h}_{\times} \end{pmatrix} = \mathcal{R}_{2\psi} \begin{pmatrix} h_{+} \\ h_{\times} \end{pmatrix}, ~ \begin{pmatrix} \bm{\epsilon}^{+} \\ \bm{\epsilon}^{\times} \end{pmatrix} = \mathcal{R}_{2\psi} \begin{pmatrix} \bm{e}^{+} \\ \bm{e}^{\times} \end{pmatrix} \end{aligned}??$$

这里$\psi$就是引力波在 $\hat{p}, \hat{q}$ 参考系下的偏振角。注意，引力波是四极辐射，因此实际旋转矩阵角度为$2\psi$，旋转$\frac{\pi}{4}$两种偏振模式互换(伸缩反相)，旋转$\frac{\pi}{2}$回到初始偏振(伸缩反相)。

下面讨论探测器对引力波的响应，如下图所示，考虑直角的迈克尔逊干涉仪，沿引力波探测臂建立坐标系，引力波沿$\hat{k}$方向到达探测器，对应方位角$(\theta, \varphi)$。这里引力波由$x,y$平面下向上传播，由上至下的传播对应$\theta >\pi/2$。

We denote by $\theta, \varphi$ the angles that define the propagation direction $\hat{k}$ of the GW from the source to us, so the polar angles of the source, as seen from the Earth, are $\theta_s = \pi - \theta$ and $\varphi_s = \pi + \varphi$.

首先考虑单臂情况，探测信号对应应变张量在探测臂方向的投影，对$x$轴方向的探测臂：

$$\boxed{h_{xx} = \frac{1}{2}~\hat{x}\otimes\hat{x} : \bm{h} = \frac{1}{2}~\hat{x}\cdot \bm{h} \cdot \hat{x}}$$

其中$\hat{x}\otimes \hat{x}$为$x$方向上的投影算符。代入引力波应变张量：

$$\begin{aligned} h_{xx} = \frac{1}{2}~\hat{x}\otimes\hat{x} : \begin{pmatrix} \bar{h}_{+} & \bar{h}_{\times} \end{pmatrix} \begin{pmatrix} \bm{\epsilon}^{+}\\ \bm{\epsilon}^{\times} \end{pmatrix} = \frac{1}{2}\begin{pmatrix} \bar{h}_{+} & \bar{h}_{\times} \end{pmatrix} \begin{pmatrix} \hat{x}\cdot \bm{\epsilon}^{+} \cdot \hat{x}\\ \hat{x} \cdot \bm{\epsilon}^{\times} \cdot \hat{x} \end{pmatrix} \end{aligned}$$

同时，探测信号用响应函数表示有：

$$h = F^P h_P = \begin{pmatrix} h_{+} & h_{\times} \end{pmatrix} \begin{pmatrix} F^{+} \\ F^{\times} \end{pmatrix} = \begin{pmatrix} \bar{h}_{+} & \bar{h}_{\times} \end{pmatrix} \begin{pmatrix} D^{+} \\ D^{\times} \end{pmatrix}$$

其中$A$对应偏振模式，这里仅考虑广义相对论框架下的两种偏振，对比两式有：

$$\boxed{ F_{x}^P = \frac{ \hat{x}\cdot \bm{e}^P \cdot \hat{x} }{2};~~~ D_{x}^P = \frac{ \hat{x}\cdot \bm{\epsilon}^P \cdot \hat{x} }{2} }$$

$$\begin{aligned} D_{x}^{+} &= \frac{\hat{x}\cdot \bm{\epsilon}^{+} \cdot \hat{x}}{2} = \frac{(\hat{p}\cdot \hat{x})^2 - (\hat{q}\cdot \hat{x})^2}{2} \\ D_{x}^{\times} &= \frac{\hat{x} \cdot \bm{\epsilon}^{\times} \cdot \hat{x}}{2} = (\hat{p}\cdot \hat{x})(\hat{q}\cdot \hat{x}) \end{aligned}$$

根据球坐标系与直角坐标基矢关系 $\hat{\theta}\cdot\hat{x} = \cos\theta\cos\varphi, ~\hat{\varphi}\cdot\hat{x} = {\footnotesize -}\sin\varphi$ 有：

$$\hat{p}\cdot\hat{x} = \cos\theta\cos\varphi, ~~~ \hat{q}\cdot\hat{x} = \sin\varphi$$

从而：

$$\begin{aligned} D_x^{+}(\theta, \varphi) &= \frac{\cos^2\theta \cos^2\varphi - \sin^2\varphi}{2}\\ D_x^{\times}(\theta, \varphi) &= \cos\theta \cos\varphi \sin\varphi\\ \end{aligned}$$

指向函数体现了天线对不同方向上信号的相对灵敏度：引力波沿探测臂方向($\theta=\frac{1}{2}\pi, \varphi=0, \pi$)，探测器完全无响应；引力波垂直探测臂方向($\varphi=\frac{1}{2}\pi,\frac{3}{2}\pi$)时，对$+$模式响应取最大值$\frac{1}{2}$，但对$\times$模式同样无响应。疑问：当$\theta=0, \pi$，$\varphi$取值不改变引力波传播方向，但响应函数依赖$\varphi$，这点如何理解？

类似的，对$y$方向的探测器 $h_{yy}(t) = \frac{1}{2}~\hat{y}\otimes\hat{y} : \bm{h}(t)$

$$\hat{p}\cdot\hat{x} = \cos\theta\sin\varphi ; ~~~ \hat{q}\cdot\hat{y} = -\cos\varphi$$

$$\begin{aligned} D_{y}^{+} &= \frac{\hat{y}\cdot \bm{\epsilon}^{+} \cdot \hat{y}}{2} = \frac{(\hat{p}\cdot \hat{y})^2 - (\hat{q}\cdot \hat{y})^2}{2} \\ D_{y}^{\times} &= \frac{\hat{y} \cdot \bm{\epsilon}^{\times} \cdot \hat{y}}{2} = (\hat{p}\cdot \hat{y})(\hat{q}\cdot \hat{y}) \end{aligned}$$

从而：

$$\begin{aligned} D_y^{+}(\theta, \varphi) &= \frac{\cos^2\theta \sin^2\varphi {\footnotesize -} \cos^2\varphi}{2}\\ D_y^{\times}(\theta, \varphi) &= -\cos\theta \sin\varphi \cos\varphi \end{aligned}$$

对于迈克尔逊双臂干涉，总体响应为两个方向响应之差，即：

$$\begin{aligned} D^P = \frac{1}{2}\hat{l}_1\otimes\hat{l}_1 : \bm{e}^P - \frac{1}{2}\hat{l}_2\otimes\hat{l}_2 : \bm{e}^P = \frac{\hat{l}_1\otimes\hat{l}_1 - \hat{l}_2\otimes\hat{l}_2}{2} : \bm{e}^P \end{aligned}$$

注意，在这种形式本身不要求干涉臂方向 $\hat{l}_1, \hat{l}_2$ 垂直，只需要将
这里干涉臂方向矢量不用$\hat{x},\hat{y}$，用$\hat{p},\hat{q}$、$\hat{u},\hat{v}$或 $\hat{l}_1, \hat{l}_2$

对于直角的迈克尔逊双臂干涉而言，总体响应为两者之差：

$$\begin{aligned} D^P = \frac{\hat{x}\otimes\hat{x} - \hat{y}\otimes\hat{y}}{2} : \bm{e}^P \end{aligned}$$

$$\begin{aligned} D^{+}(\theta, \varphi) &= D_x^{+} - D_y^{+} = \frac{1}{2}(1+\cos^2\theta)\cos2\varphi\\ D^{\times}(\theta, \varphi) &= D_x^{\times} - D_y^{\times} = \cos\theta\sin2\varphi \end{aligned}$$

$$D_{\rm rms}(\theta, \varphi) = \sqrt{ {D^{+}}^2 + {D^{\times}}^2 } = \sqrt{\frac{1}{4}\sin^4\theta\cos^2 2\varphi + \cos^2\theta}$$

对单独的极化模式：
两种极化模式均方根响应：
垂直探测器平面时响应最强
探测器角平分线或平行探测器平面且垂直角平分线时响应最弱
立体图和天球图两个视角下理解

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

def plot_beam_pattern(F, theta, phi, fig, viewangle=()):
    # Convert spherical coordinates to Cartesian coordinates
    r = np.abs(F)
    x = r * np.sin(theta) * np.cos(phi)
    y = r * np.sin(theta) * np.sin(phi)
    z = r * np.cos(theta)

    # Create a 3D plot
    ax = fig.add_subplot(111, projection='3d')
    ax.view_init(10, 35) #view angle: inclination, azimuth

    # Plot the surface with the 'Blues' colormap to represent F values
    surf = ax.plot_surface(x, y, z, facecolors=plt.cm.Blues(F))#, shade=False)

    # Set labels for the axes
    ax.set_xlabel(r'$x$')
    ax.set_ylabel(r'$y$')
    ax.set_zlabel(r'$z$')
    ax.set_title(r'$F^{+}$')

    # Create a ScalarMappable to associate colormap with values
    mappable = plt.cm.ScalarMappable(cmap='Blues')
    mappable.set_array(F)

    # Add a color bar, specifying the axis (ax) to place it
    fig.colorbar(mappable, label='F(theta, phi)', ax=ax)
    ax.set_box_aspect([1, 1, 1.8])

# Define the range of theta and phi
theta_range = np.linspace(0, np.pi, 100)
phi_range = np.linspace(0, 2*np.pi, 100)

# Create a grid of theta and phi values
theta, phi = np.meshgrid(theta_range, phi_range)

# Calculate the corresponding values of F(theta, phi)
# F = np.cos(theta) * np.cos(phi) * np.sin(phi)
# F = (np.cos(theta)**2*np.cos(phi)**2 - np.sin(phi)**2)/2
# F_plus = (1+np.cos(theta)**2)/2*np.cos(2*phi)
# F_cross = np.cos(theta)*np.sin(2*phi)
F = np.sqrt(np.sin(theta)**4*np.cos(2*phi)**2/4 + np.cos(theta)**2)


fig = plt.figure()
plot_beam_pattern(F, theta, phi, fig)

# Show the plot
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from astropy.coordinates import SkyCoord

# Define the range of longitude and latitude
lon_range = np.linspace(-np.pi, np.pi,360)
lat_range = np.linspace(-np.pi/2., np.pi/2.,180)

# Create a grid of lon and lat values
lon, lat = np.meshgrid(lon_range,lat_range)
theta, phi = np.meshgrid(theta_range, phi_range)



# Convert lon, lat to theta, phi and Calculate F(theta, phi)
theta = np.pi/2 - lat
phi = (lon + 2*np.pi) % (2*np.pi)
F_plus = (1+np.cos(theta)**2)/2*np.cos(2*phi)
F_cross = np.cos(theta)*np.sin(2*phi)
F = np.sqrt(F_plus**2 + F_cross**2)

# Plot the function F on the Skymap
# fig = plt.figure(figsize=(8, 4.2))
ax = plt.subplot(111, projection="mollweide")
c = ax.pcolormesh(lon, lat, F, cmap='Blues')#, shading='auto')
plt.colorbar(c, shrink=0.5)

plt.title(r'Mollweide Projection of $F$', y=1.1)
plt.grid(color='white', linestyle='--', linewidth=0.5)  # Add grid lines
plt.show()

最后，如果考虑到偏振角$\psi$：

$$\begin{aligned} \begin{pmatrix} \bar{h}_+ & \bar{h}_\times \end{pmatrix} \begin{pmatrix} D^{+}\\ D^{\times} \end{pmatrix} = \begin{pmatrix} h_+ & h_\times \end{pmatrix} \mathcal{R}_{2\psi} \begin{pmatrix} D^{+}\\ D^{\times} \end{pmatrix} = \begin{pmatrix} h_+ & h_\times \end{pmatrix} \begin{pmatrix} F^{+}\\ F^{\times} \end{pmatrix} \end{aligned}$$