引力波天文III:空间探测技术

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

本文介绍引力波的空间探测技术。

空间观测概述

相关项目介绍

低频引力波
主流毫赫兹(Milli-hertz, mHz),有欧洲和美国合作的LISA以及国内的太极、天琴项目。
毫赫兹频段两边的分赫兹(Deci-hertz, dHz)及微赫兹(Micro-hertz, µHz)频段同样受到关注。
为什么要做空间探测

主流mHZ频段
激光干涉空间天线(Laser Interferometer Space Antenna)是

https://en.wikipedia.org/wiki/Laser_Interferometer_Space_Antenna#History
引力波的空间探测实验探索始于上世纪80年代,并于90年代正式提出LISA项目。1997年NASA加入,LISA项目已基本成型,被设计为呈三角分布的三颗卫星,臂长500万公里。2011年NASA退出项目,LISA项目降级为eLISA,三臂设计变为两臂设计,臂长缩水为100万公里。2015年引力波地面探测的成功,以及LISA Pathfinder的超预期结果,NASA再次加入项目,eLISA重新升级为LISA,恢复三臂设计,但臂长选择为250万公里。

三臂、500万公里
两臂、100万公里
三臂、250万公里

Pathfinder 2015年12月3日发射
无拖曳控制、重力感应、激光外差测量

地月距离36-40万公里,LISA/Taiji臂长约地月距离的10倍,天琴臂长约地月距离的一半。
Taiji 300万公里
天琴 17万公里

其它频段项目

卫星轨道简介

无拖曳卫星编队(satellite formation)
无拖曳(zero-drag or drag-free)卫星是指目标载荷(测试质量)不受引力之外的其它力作用的卫星。为实现该目标,卫星会将检验质量或测试质量(proof mass)置于真空壳层内,测试质量在壳层内自由悬浮,并实时监测其与外壁的距离。任何距离的改变都代表卫星受到了非重力因素影响,相对测试质量出现了运动位移,卫星会启动微推进器做出补偿调整,恢复两者相对位置。在外壳的保护下,太阳风、光压等非重力因素的影响被屏蔽,测试质量仅在重力作用下沿测地线运动。而通过追踪测试质量,卫星本身运动也始终追随测地线。

惯性传感器+微推进系统,其中惯性传感系统被称为重力参考传感器GRS(Gravitational Reference Sensor),包括测试质量+电极室+真空室及相关传感器。

无拖曳卫星需对测试质量在三个方向的平移及旋转,共计六个自由度上,对测试质量实现追踪。但对三臂的空间干涉阵而言,每颗卫星都有两个测试质量,分别用于与其余两颗卫星进行连接。显然卫星无法在所有自由度上同时对两者保持追踪,因此实际中,只会在与其它卫星连接的视线方向上保持测试质量完全自由。即对每个测试质量只在一个方向保持追踪,其余自由度全部由电容制动装置加以控制。若两个追踪方向相互正交相对易理解,追踪60°夹角的两个方向虽然不直观,但仍是可行的。

对于LISA而言,测试质量的精度可达到3×1015m/s2/Hz1/23\times 10^{-15}\rm m/s^2/{Hz}^{1/2}1mHz1{\rm mHz}

最简单的想法是类似地面,由同一光源分光后经过不同路径反射回本地干涉,但对于卫星干涉测量:

  • 一方面,数千至数百万公里以上的臂长,意味激光强度衰减严重,通过反射往返测量不现实。以LISA为例,取激光出射功率为1W,激光接受功率仅约250pW,有近10个量级的衰减。因此,实际中会采用单程测量,即各卫星分别用本地激光指向其它卫星,三臂往返一共六路激光测量数据。
  • 另一方面,三个卫星各自独立控制,考虑到月球及其它行星的引力作用,卫星的臂长在不断波动(“呼吸”),测试质量间存在相对运动。这会导致:
    • 不同卫星间臂长不相等。以LISA为例,臂长波动在1%量级,即几万公里的相对位移。激光频率的不稳定性在1Hz、kHz、MHz量级?
    • 卫星间往返路径不相等,除了不同卫星间
    • 测试质量间存在相对运动,引入激光频率的多普勒偏移。
      其“呼吸”频率在接近年的量级,低于目标引力波的频率

由同一光源分光经过相同长度的路径后自动干涉,无需考虑频率波动、无需测量初始相位。

考虑到

  • 臂长不相等:时间延迟干涉技术
  • 存在多普勒效应:外差干涉测量

空间引力波每条路径的光束与标准的参考激光对比,精确记录各自的相位,一共记录三臂往返的六路激光数据。不同路径的激光干涉是在数据(预)处理阶段离线进行。因为探测器响应记录到的是每路激光的相位延迟,并未直接提供不同路径光的干涉信息。

事实上,每路激光数据也是由三个干涉器共同完成测量。

通过轨道的合理设置,星间距离保持稳定。忽略干扰因素,三颗卫星各自沿测底线做开普勒运动,轨道平面与黄道面夹角一致,但升交点依次相错2π3\frac{2\pi}{3},在黄道上均匀排布。最终,理想情况下:

  • 三颗卫星的中心位于黄道面,沿地球轨道以年为周期公转,与地球保持约20°;
    20°是地球潮汐扰动与卫星通信需求两相权衡的折中选择
  • 三颗卫星构成的平面公转过程中始终面向太阳,与黄道面保持约60°倾斜;
    保持与黄道倾角的稳定,有助于保持卫星系统的温度稳定
  • 三颗卫星保持60°夹角/等边三角形,整体在卫星平面中绕中心以年为周期(反向)“自转”;
    转动引入的额外信号调制有助于提升引力波源系统的定位

但实际中,会受到各种扰动

  • 臂长存在约2%的波动(5万公里),非等臂干涉激光自身频率噪声无法在干涉测量时直接抵消
  • 臂长波动对应卫星视向相对速度约10m/s,激光频率偏移不再适用零差干涉,需用外差干涉
  • 臂长波动意味不再保持等边,卫星夹角将偏离60°,偏差在1°左右,激光指向需相应进行调整

数据链路简介

  • 六路激光多普勒数据
  • 卫星轨道

最直观的想法是激光在测试质量上反射后干涉,这也确实是上世纪LISA最初提出时的设计方案,但之后在2005年左右切换为现在的分立干涉(split interferometry)的方案,即将测试质量到测试质量的距离测量拆分为测试质量到光学平台、光学平台到光学平台、光学平台到测试质量三部分,如下图所示,注意与干涉测量中的分光(split)区分。

每颗卫星上有两个科学单元,如下图所示,每个单元的主要组件有:

  • 1个测试质量,位于惯性传感器(重力参考传感器)内部
  • 1个激光系统,红外1064nm Nd:YAG 稳频稳功率激光
  • 1个望远镜,用于发射和接受激光光束,口径30cm
  • 1个光学平台,位于望远镜与GRS之间,承载干涉光路

科学单元正式名称为可动光学子组件(Movable Optical Sub-Assembly, MOSA),每个组件包含一个望远镜、一个光学平台(Optical Bench, OB)以及一个重力参考传感器(Gravitational Reference Sensor, GRS),GRS内部则是自由悬浮的测试质量(Test Mass, TM)。

测试质量悬浮于电容电极室(capacitive electrode housing),通过电容实时监测测试其位置,并通过DFACS()系统实现现卫星的无拖曳控制

重力参考传感器GRS通过激光?/真空室电容?对测试质量位置进行监测,实现卫星的无拖曳控制。每颗卫星上的两个科学单元固定在同一支架上,各自有±1°的旋转自由,用于卫星漂移时保持激光指向对齐。两个科学单元的光学平台之间由光纤连接,用于两个本地激光之间的干涉。光学平台上共有3个光电探测器,分别监测3路干涉拍频信号,经低通滤波和ADC转换后由相位计读取相位。

载荷概念设计,两个科学单元挂载在同一支架,呈60°夹角(Amaro-Seoane et al. 2017)

GRS与DFACS连接

每个科学单元(光学平台)的3个光电探测器分别用于三个干涉测量:

分立干涉光路简化示意图,主要包含科学干涉、参考干涉及测试质量干涉
  • 科学干涉:本地激光与远程激光干涉,含GW信息
  • 参考干涉:本地激光与相邻光学平台的激光干涉
  • 本地干涉:本地激光经TM反射后与相邻激光干涉

其中科学干涉用于光学平台到光学平台测距,光学平台到测试质量的距离可通过对比参考干涉和本地干涉得到,两者结合最终实现测试质量到测试质量的测距。此外,为实现干涉测量需要在相邻的本地激光信号间建立关联,以构建干涉环路,这同样依赖参考干涉。作为对比,早期的两路干涉光路设计为:

主要区别在于本地干涉与科学干涉合并,远程激光会直接在测试质量上反射。这里科学干涉测量的是光学平台到测试质量的距离,仍需要参考干涉辅助实现测试质量到测试质量的距离测量。

  • 参考干涉测量的作用:
    • 与本地干涉测量对比抵消噪声,获取光学平台相对于测试质量的位移,即光学平台的位移噪声
    • 通过双向参考干涉与科学干涉测量对比,实现科学数据中本地与相邻激光频率噪声的转换
    • 前面一点也就是在全部六个激光之间的实现锁相

科学干涉,也称星间干涉(Interspacecraft Interferometer, ISI),包含本地及远程的激光频率噪声,两个光学平台位移噪声,引力波信息以及读出噪声。参考干涉(Reference Interferometer, RFI)包含两个本地激光频率噪声、连接连个光学平台的光纤所引入的相位延迟以及读出噪声。本地干涉,也称测试质量干涉(Test-mass Interferometer, TMI),包含两个本地激光频率噪声,光纤相位延迟、光学平台的位移噪声、测试质量的加速噪声及读出噪声。这些就是空间引力观测所涉及的主要噪声,除此之外,星间信息传递(频率同步)时信息失真对应于调制噪声,相位计读出采样时还会引入时钟噪声,这里暂不考虑。

所有其它噪声如时钟噪声或指向噪声的效应都要在TDI之后才能显现,相应的噪声抑制或修正也都建立在TDI的输出之上,而TDI本身又可能对这些噪声产生影响。
经过TDI处理后,数据是由三颗卫星各自的固有时?坐标时?表示,需要进一步将其转换到统一的时间参考系,如质心坐标时TCB,以方便后续分析。

其它噪声

  • 指向噪声
    指向噪声即TTL(tilt-to-length)噪声,是指卫星指向角抖动引入的光程变化。在探测器运行过程中会对望远镜指向进行实时的校准修正,在硬件级别直接抑制TTL噪声。但残余噪声仍需要在数据预处理阶段进一步修正。具体的,波前倾角可由差分波前感知技术(Differential Wavefront Sensing, DWS)进行监测,噪声抑制的关键是确定角度与光程的耦合系数。考虑到TDI对多路观测数据进行延时组合,指向噪声也会相应的通过线性组合混杂在一起。
  • 时钟噪声
    TDI可以同时消除时钟噪声和激光噪声

为使用数学符号表示,首先对标号做出约定:

这里采用LISA项目的最新记号规则(Hartwig et al. 2021, Bayle et al. 2021)。如上图所示,卫星沿顺时针依次记为1、2、3;卫星上的科学单元用两个下标指示,前者对应所在卫星,后者对应所指向的卫星;连接臂标号与其所指向的科学单元一致,沿其传播的激光最终会在对应科学单元上进行干涉。

s(t),ε(t),τ(t)s(t), \varepsilon(t), \tau(t)分别代表科学、本地及参考三个干涉的相位信号,有:

s12(t)=p21(tL12/c)p12(t)2πν21Δ21(tL12/c)Δ12(t)c+H12(t)+N12s(t)ε12(t)=p13(t)p12(t)2 × 2πν13δ12(t)Δ12(t)c+μ12(t)+N12ε(t)τ12(t)=p13(t)p12(t)+μ12(t)+N12τ(t)\begin{aligned} s_{12}(t) &= p_{21}(t-L_{12}/c) - p_{12}(t) - 2\pi\nu_{21}{\small \frac{\Delta_{21}(t-L_{12}/c) - \Delta_{12}(t)}{c}} + H_{12}(t) + N_{12}^{s}(t) \\ \varepsilon_{12}(t) &= p_{13}(t) - p_{12}(t) - 2{\small~\times~} 2\pi\nu_{13}{\small \frac{\delta_{12}(t) - \Delta_{12}(t)}{c}} + \mu_{12}(t) + N_{12}^{\varepsilon}(t)\\ \tau_{12}(t) &= p_{13}(t) - p_{12}(t) + \mu_{12}(t) + N_{12}^{\tau}(t)\\ \end{aligned}

其中p(t)p(t)为激光相位噪声;ν\nu为激光基准频率,不同激光会有所区别;Δ(t)\Delta(t)为光学平台相对测试质量的位移在连接臂方向的投影;δ(t)\delta(t)为测试质量对测地线的偏离在连接臂方向的投影;H(t)H(t)为引力波所引入的相位变化;N(t)N(t)为光电探测器读出噪声(激光散粒噪声),三个干涉测量读出噪声相互独立;μ(t)\mu(t)为连接两个光学平台的光纤所引入的相位变化,下标与激光的目标光学平台一致。

注意,参考干涉中的位移项δ12(t)Δ12(t)c\frac{\delta_{12}(t) - \Delta_{12}(t)}{c},其正负号与光路设计有关。具体的跟在TM上反射的激光以及激光打向TM的方向有关。前者有两种选择,常见的是相邻光学平台的激光在TM上反射,而最新的Bayle & Hartwig 2022中的光路图是本地激光在TM上反射,这将直接决定相位做差时,位移项属于减数(本地激光p12p_{12})还是被减数(相邻激光p13p_{13})。同时还要考虑激光指向,当激光打向TM的方向与远程激光传播方向一致时,位移项对相应激光相位的贡献为正,方向相反时,位移项对相应激光贡献为负。这两点的具体选择似乎没有定论,需要小心,比如这里是本地激光在TM上反射,位移项属于减数,同时激光指向与远程激光传播方向一致,相位贡献为正,最终整体上取负号。

引入沿激光臂方向的延时算符及单位矢量,上述公式可简记为:

s12=D12p21p122πν21n^12c(D12Δ21Δ12)+H12+N12sε12=p13p122 × 2πν13n^12c(δ12Δ12)+μ12+N12ετ12=p13p12+μ12+N12τ\begin{aligned} s_{12} &= \mathcal{D}_{12} p_{21} - p_{12} - 2\pi\nu_{21} \frac{\hat{\bm{n}}_{12}}{c} \cdot (\mathcal{D}_{12} \bm{\Delta}_{21} - \bm{\Delta}_{12}) + H_{12} + N_{12}^{s}\\ \varepsilon_{12} &= p_{13} - p_{12} - 2{\small~\times~} 2\pi\nu_{13}\frac{\hat{\bm{n}}_{12}}{c} \cdot (\bm{\delta}_{12} - \bm{\Delta}_{12}) + \mu_{12} + N_{12}^{\varepsilon}\\ \tau_{12} &= p_{13} - p_{12} + \mu_{12} + N_{12}^{\tau} \end{aligned}

ϕsci12=D12p21p122πν21n^12c(D12Δ21Δ12)+H12+N12sϕloc12=p13p122 × 2πν13n^12c(δ12Δ12)+μ12+N12εϕref12=p13p12+μ12+N12τ\begin{aligned} \phi_{ {\rm sci}_{12} } &= \mathcal{D}_{12} p_{21} - p_{12} - 2\pi\nu_{21} \frac{\hat{\bm{n}}_{12}}{c} \cdot (\mathcal{D}_{12} \bm{\Delta}_{21} - \bm{\Delta}_{12}) + H_{12} + N_{12}^{s}\\ \phi_{ {\rm loc}_{12} } &= p_{13} - p_{12} - 2{\small~\times~} 2\pi\nu_{13}\frac{\hat{\bm{n}}_{12}}{c} \cdot (\bm{\delta}_{12} - \bm{\Delta}_{12}) + \mu_{12} + N_{12}^{\varepsilon}\\ \phi_{ {\rm ref}_{12} } &= p_{13} - p_{12} + \mu_{12} + N_{12}^{\tau} \end{aligned}

s12=D12p21p122πν21n^12c(D12Δ21Δ12)+H12+νbeat+N12sε12=p13p122 × 2πν13n^12c(δ12Δ12)+μ12+N12ετ12=p13p12+μ12+N12τ\begin{aligned} s_{12} &= \mathcal{D}_{12} p_{21} - p_{12} - 2\pi\nu_{21} \frac{\hat{\bm{n}}_{12}}{c} \cdot (\mathcal{D}_{12} \bm{\Delta}_{21} - \bm{\Delta}_{12}) + H_{12} + \frac{\nu_{\rm beat}}{} + N_{12}^{s}\\ \varepsilon_{12} &= p_{13} - p_{12} - 2{\small~\times~} 2\pi\nu_{13}\frac{\hat{\bm{n}}_{12}}{c} \cdot (\bm{\delta}_{12} - \bm{\Delta}_{12}) + \mu_{12} + N_{12}^{\varepsilon}\\ \tau_{12} &= p_{13} - p_{12} + \mu_{12} + N_{12}^{\tau}\\ \end{aligned}

在TDI分析中会通过延时操作构建虚拟等臂干涉,在这一过程通常默认每颗卫星上只有一个激光噪声源。而上述干涉测量涉及六路激光,需要消除其中三个,得到等效的分光干涉测量,每颗卫星只有一个激光。

中间变量,这是后续TDI分析的基础。

η12=\begin{aligned} \eta_{12} &= \\ \end{aligned}

注意,科学干涉中本地激光在W量级,而远程激光仅0.1nW量级,约10个量级的差别;参考及本地干涉中的相邻激光也只有mW量级。此外参考干涉是双向的:当前科学单元的本地激光会与经光纤过来的相邻科学单元激光干涉,相邻科学单元的激光与会经光纤过去的本地激光干涉作为参考干涉。假设光纤的相位延迟是对易的,即从左向右传播与从右向左传播的相位延迟相等,则双向参考干涉相减,光纤噪声抵消,可得到两个科学单元激光相位(频率噪声)间的关联。

TDI 技术需要高精度的绝对距离数据。结合激光相位伪码调制,LISA 的 TDI 系统需在 5x106 km 臂长下实现绝对距离测量精度 30cm
相对误差6x10-11量级?引力波在pm量级???
TDI需要时间以及臂长信息用于确定延迟时间,时间精度要达到纳秒量级,而测距需要达到米多的量级
需要在卫星间同步时钟,Kalman滤波
需要控制拍频,频率同步locking、频率规划
次级噪声:TTL噪声、时钟噪声

在空间引力波探测中,两卫星之间最大速度约为 10m/s,对应相位计带宽范围约是 2MHz~25MHz,pm量级测量要求鉴相精度达到一个圆周的一百万分之一

涉及到的物理量之间的转换
用于测量的激光信号可用电场强度表征:

E(t)=E0(t)ejφ(t)=E0(t)ej[2πν0t+ϕ(t)]E(t) = E_0(t) e^{j\varphi(t)} = E_0(t) e^{j[2\pi\nu_0 t + \phi(t)]}

这里E0(t)E_0(t)为信号的(瞬时)振幅,受激光功率不稳定影响,振幅存在波动;φ(t)\varphi(t)为信号的瞬时相位,其导数为瞬时频率。ν0\nu_0为激光的中心频率,ϕ(t)\phi(t)则是瞬时相位φ(t)\varphi(t)相对线性增长的基准相位2πν0t2\pi\nu_0 t的波动残差。主要的噪声以及引力波信号都位于ϕ(t)\phi(t)中。

ω(t)=φ˙(t),   ν(t)=φ˙(t)2π=ν0+ϕ˙(t)2π;   φ(t)=2πt0tν(τ)dτ\omega(t) = \dot{\varphi}(t), ~~~ \nu(t) = \frac{\dot{\varphi}(t)}{2\pi} = \nu_0 + \frac{\dot{\phi}(t)}{2\pi}; ~~~ \varphi(t) = 2\pi \int_{t_0}^t\nu(\tau)d\tau

除了相位与频率外,另两个常用的描述相位演化的物理量是计时抖动(timing jitter)x(t)x(t)及其导数——相对频偏(relative frequency deviation)或分数多普勒频移(fractional doppler shift)y(t)y(t)。其中计时抖动是以中心频率ν0\nu_0作为基准,将相位残差ϕ(t)\phi(t)转换为时间,反映了瞬时相位φ(t)\varphi(t)所对应的时间与实际(线性)时间之间的偏离,即所谓抖动,其导数对应于激光瞬时频率ν(t)\nu(t)与中心频率ν0\nu_0的相对偏离。

x(t)=ϕ(t)2πν0=φ(t)2πν0t;     y(t)=x˙(t)=ν(t)ν0ν0x(t) = \frac{\phi(t)}{2\pi\nu_0} = \frac{\varphi(t)}{2\pi\nu_0} - t; ~~~~~ y(t) = \dot{x}(t) = \frac{\nu(t)-\nu_0}{\nu_0}

在激光的接收端,考虑到信号的衰减及时延:

Erecv(t)=E0recv(t)ej[2πν0t+ϕrecv(t)]=ρEemit(tL/c)E^{\rm recv}(t) = E^{\rm recv}_0(t) e^{j[2\pi\nu_0 t + \phi^{\rm recv}(t)]} = \rho E^{\rm emit}(t - L/c)

ϕrecv(t)=ϕ(tL/c)2πν0L/c;     E0recv(t)=ρE0(t)\phi^{\rm recv}(t) = \phi(t - L/c) - 2\pi\nu_0 L/c; ~~~~~ E^{\rm recv}_0(t) = \rho E_0(t)

yrecv(t)=  ϕ˙recv(t)2πν0=ϕ˙(tLc)2πν0(1L˙(t)/c)L˙(t)c=  ϕ˙(tLc)2πν0(1+ϕ˙(tLc)2πν0)L˙(t)c=  y(tL/c)[1+y(tL/c)]L˙(t)c\begin{aligned} y^{\rm recv}(t) = ~&~ \frac{\dot{\phi}^{\rm recv}(t)}{2\pi\nu_0} = \frac{\dot{\phi}(t-\frac{L}{c})}{2\pi\nu_0}\left(1-\dot{L}(t)/c\right) - \frac{\dot{L}(t)}{c}\\ = ~&~ \frac{\dot{\phi}(t-\frac{L}{c})}{2\pi\nu_0} - \left(1 + \frac{\dot{\phi}(t-\frac{L}{c})}{2\pi\nu_0}\right)\frac{\dot{L}(t)}{c}\\ = ~&~ y(t-L/c) - \left[1 + y(t-L/c)\right]\frac{\dot{L}(t)}{c} \end{aligned}

这里考虑了臂长的变化,第一项为发射端激光本身的相对频率波动,后一项为臂长变化所引入的相对频移,对应实际的多普勒效应。注意多普勒效应源于L˙\dot{L},因此在ϕrecv(t)\phi^{\rm recv}(t)是无法直接体现的,只有在νrecv(t)\nu^{\rm recv}(t)y(t)y(t)中才有对应的额外项。考虑到卫星相对速度10m/s\sim 10 \rm m/sL˙/c3×1081, y(t)1\dot{L}/c \sim 3\times 10^{-8}\ll 1, ~y(t) \ll 1,忽略两者相乘的高阶项(Bayle中似乎是1m/s??):

yrecv(t)=y(tL/c)L˙(t)cy^{\rm recv}(t) = y(t-L/c) - \frac{\dot{L}(t)}{c}

卫星的相对速度是10m/s,还是1m/s?
y(t)在什么量级??L˙(t)c\frac{\dot{L}(t)}{c}可以通过轨道监测消除,为什么能直接忽略??
δν1Hz,y(t)1015L˙c;δν10MHz,y(t)3×108L˙c\delta \nu \sim 1{\rm Hz}, y(t) \sim 10^{-15} \ll \frac{\dot{L}}{c}; \delta \nu \sim 10{\rm MHz}, y(t) \sim 3\times10^{-8} \sim \frac{\dot{L}}{c}

外差测量
E1(t)ejπ/2+E2(t)E_1(t)e^{j\pi/2} + E_2(t)pi/2pi/2是分光镜反射所引入的??

E1(t)+E2(t)E_1(t) + E_2(t)

I1+2E1,0E2,0E1,02+E2,02cos[2π(ν1ν2)t+ϕ1(t)ϕ2(t)]I \propto 1 + \frac{2 E_{1, 0}E_{2, 0}}{E^2_{1, 0} + E^2_{2, 0}} \cos[2\pi(\nu_1-\nu_2)t + \phi_1(t)-\phi_2(t)]

ν12=ν1ν2,   ϕ12=ϕ1(t)ϕ2(t)\nu_{12} = \nu_1 - \nu_2, ~~~ \phi_{12} = \phi_1(t) - \phi_2(t)

y12=ϕ˙122πν0=y1(t)y2(t)y_{12} = \frac{\dot{\phi}_{12}}{2\pi\nu_0} = y_1(t) - y_2(t)

引力波频率用ff,电磁波频率用ν\nu,TM是测试质量,OMS是光学测量系统,BBS日心系,GRS、ADC、TDI
optical metrology system.

引力波信号对于天线臂长(卫星间距离)的扰动可表示为:

h(t)=h+SSBξ+(u,v,n)+h×SSBξ×(u,v,n)h(t) = h^{\rm \tiny SSB}_+ \xi_+(\bm{u}, \bm{v}, \bm{n}) + h^{\rm \tiny SSB}_× \xi_×(\bm{u}, \bm{v}, \bm{n})

其中h+SSB,h×SSBh^{\rm \tiny SSB}_+, h^{\rm \tiny SSB}_×为引力波振幅在太阳质心(SSB)坐标系的表示,ξ+,ξ×\xi_+, \xi_×为天线响应函数。

t1t0+Lc12c0LH(t0kx0(t0)c+1kn^cx)dxt_1 \approx t_0 + \frac{L}{c} - \frac{1}{2c} \int_0^L H\left(t_0 - \frac{\bm{k}\cdot\bm{x}_0(t_0)}{c} + \frac{1-\bm{k}\cdot\bm{\hat{n}}}{c} x \right) d x

yGW(t0)12kn^[H(t0kx0(t0)c)H(t0kx1(t1)c+Lc)]y^{\rm GW}(t_0) \approx \frac{1}{2\bm{k}\cdot\hat{\bm{n}}}\left[ H\left(t_0 - \frac{\bm{k}\cdot\bm{x}_0(t_0)}{c}\right) - H\left(t_0 - \frac{\bm{k}\cdot\bm{x}_1(t_1)}{c} + \frac{L}{c}\right)\right]

x1(t1)=x0(t0)+Ln^\bm{x}_1(t_1) = \bm{x}_0(t_0) + L\hat{\bm{n}}

yGW(t1)12kn^[H(t1Lckx0(t0)c)H(t1kx1(t1)c)]y^{\rm GW}(t_1) \approx \frac{1}{2\bm{k}\cdot\hat{\bm{n}}}\left[ H\left(t_1 - \frac{L}{c} - \frac{\bm{k}\cdot\bm{x}_0(t_0)}{c}\right) - H\left(t_1 - \frac{\bm{k}\cdot\bm{x}_1(t_1)}{c}\right)\right]

只考虑引力波的效应,接收端与发射端的频率相对偏移量为:

ν(t)νemit(tL/c)νemit(tL/c)=12[h(t)h(tL/c)]\frac{\nu(t)-\nu_{\rm emit}(t-L/c)}{\nu_{\rm emit}(t-L/c)} = \frac{1}{2} \Big[h(t)-h(t-L/c)\Big]

激光不稳定性10KHz量级、预稳频后在30Hz量级

数据噪声
只考虑引力波的效应,接收端与发射端的频率相对偏移量为:

ν(t)ν(tL/c)ν012[h(t)h(tL/c)]\boxed{\frac{\nu(t)-\nu(t-L/c)}{\nu_0} \approx \frac{1}{2} \Big[h(t) - h(t - L/c)\Big]}

由于激光波本身的频率不稳定性,这里发射端频率不再稳定为ν0\nu_0,而是tL/ct-L/c(发射时刻)的ν\nu
νn\nu_n表示频率本身不稳定性所贡献的频率噪声,有

ν(t)ν(tL/c)ν0νn(t)νn(tL/c)ν012[h(t)h(tL/c)]\frac{\nu(t)-\nu(t-L/c)}{\nu_0} - \frac{\nu_n(t)-\nu_n(t-L/c)}{\nu_0} \approx \frac{1}{2} \Big[h(t) - h(t - L/c)\Big]

加速度噪声a=FTMmωs2XTMa = \frac{F_{\rm TM}}{m} - \omega_s^2X_{\rm TM}
加速噪声对应(分数)频率导数,对前面的公式求导加入加速噪声有:

ν˙(t)ν˙(tL/c)ν0 12[h˙(t)h˙(tL/c)]+1ν0[ν˙n(t)ν˙n(tL/c)]  +1c[a(t)a(tL/c)]\begin{aligned} \frac{\dot{\nu}(t)-\dot{\nu}(t-L/c)}{\nu_0} \approx ~ & \frac{1}{2} \Big[\dot{h}(t) - \dot{h}(t - L/c)\Big] + \frac{1}{\nu_0}\Big[\dot{\nu}_n(t)-\dot{\nu}_n(t-L/c)\Big] \\ ~&~ + \frac{1}{c}\Big[a(t) - a(t-L/c)\Big] \end{aligned}

散粒噪声是光电感应测量时引入的?5pm/Hz5\rm pm/\sqrt{Hz} 功率的不稳定性? 104/Hz10^{-4}/\rm \sqrt{Hz}
激光频率抖动 2pm/Hz2\rm pm/\sqrt{Hz} 激光相对稳定性 等效1021/Hz10^{-21}/\rm \sqrt{Hz} 预稳频1013/Hz10^{-13}/\rm \sqrt{Hz}
激光指向抖动 1pm/Hz1\rm pm/\sqrt{Hz}
时钟稳定性 1013s/Hz10^{-13}\rm s/\sqrt{Hz}
光学平台光程噪声 1pm/Hz1\rm pm/\sqrt{Hz}

1pm的要求是什么意思??
干涉测距指标8pm/Hz8\rm pm/\sqrt{Hz}
加速度噪声指标3×1015ms2/Hz3\times10^{-15}\rm ms^{-2}/\sqrt{Hz}

卫星轨道分析

卫星编队的轨道分析通常使用CW参考系,即卫星编队的共动参考系:取参考点绕中心天体做圆轨道运动,以参考点为原点,坐标轴分别沿轨道径向、切向和法向,随参考点一起运动(非惯性系)。CW参考系由Clohessy & Wiltshire提出,由于卫星偏离参考点通常远小于其轨道半径,因此可在参考点附近做小量展开分析卫星运动。

在CW参考系下,满足轨道稳定性,同时与原点(参考点)距离恒定的粒子运动轨迹为与xyxy平面±60\pm 60^{\circ}夹角的平面内的圆周运动,具体可参考Dhurandhar et al. 2005。此时卫星轨道相当于大圆(参考点)叠加倾斜平面内的小圆,卫星编队位于该平面内即可实现构型稳定。对于LISA,三颗卫星依次相错2π3\frac{2\pi}{3},参考点对应卫星三角中心,运动过程中三角构型保持不变,但整体在卫星平面内旋转。

当然,上述CW参考系下的讨论相当于到偏心率一阶的近似,实际上卫星构型并不像上面所说的完全稳定,三颗卫星间距离存在一定波动,而根据Nayak et al. 2006的分析,当卫星平面倾斜角取π3+58ι\frac{\pi}{3} + \frac{5}{8}\iota时臂长更为稳定,其中ι\iota为卫星编队臂长与轨道长轴之比。

轨道运动描述

取参考点在黄道面内做圆轨道运动,半径为日地平均距离R1AU1.5×108kmR \sim 1AU \sim 1.5\times 10^8\rm km。卫星编队三角边长L2.5×106kmL\sim 2.5\times 10^6 \rm km,卫星平面与黄道面夹角ϵ\epsilon。卫星轨道半长轴同样为R=1AUR=1\rm AU,但有一定偏心率。
注意,地球轨道半长轴等于日地平均距离(相对偏差106\sim 10^{-6}),根据开普勒第三定律,参考点轨道半径及卫星轨道半长轴取RR,保证了它们周期为1年。

如上图所示,取轨道远心点为对应轨道最高点,利用简单几何关系,卫星相对黄道面高度为L3sinϵ\frac{L}{\sqrt{3}}\sin\epsilon,底边R+L3sinϵR+\frac{L}{\sqrt{3}}\sin\epsilon,从而:

tanε=L3sinϵ/(R+L3cosϵ), R(1+e)=(L3sinϵ)2+(R+L3cosϵ)2\textstyle \tan\varepsilon = \frac{L}{\sqrt{3}}\sin\epsilon \Big /(R+\frac{L}{\sqrt{3}}\cos\epsilon), ~ R(1+e) = \sqrt{\left(\frac{L}{\sqrt{3}}\sin\epsilon\right)^2 + \left(R+\frac{L}{\sqrt{3}}\cos\epsilon\right)^2}

由此卫星轨道倾角和偏心率为:

ε=arctan(23ιsinϵ1+23ιcosϵ)e=1+43ιcosϵ+43ι21\begin{aligned} \varepsilon &= \arctan\left(\frac{\frac{2}{\sqrt{3}}\iota\sin\epsilon}{1 + \frac{2}{\sqrt{3}}\iota\cos\epsilon}\right)\\ e &= \sqrt{1 + \frac{4}{\sqrt{3}}\iota \cos\epsilon + \frac{4}{3}\iota^2} - 1 \end{aligned}

其中ι=L2R0.0083\iota=\frac{L}{2R} \sim 0.0083。如果不考虑修正,取ϵ=π3\epsilon=\frac{\pi}{3}

ε=arctan(ι1+ι/3)e=1+23ι+43ι21\boxed{\begin{aligned} \varepsilon &= \arctan\left(\frac{\iota}{1 + \iota/\sqrt{3}}\right)\\ e &= \sqrt{1 + \frac{2}{\sqrt{3}}\iota + \frac{4}{3}\iota^2} - 1 \end{aligned}}

对小量ι\iota展开到一阶近似有:

ει3e\varepsilon \sim \iota \sim \sqrt{3}e

显然臂长LL越短,倾角越小,偏心率也越小。

LISA轨道偏心率e0.0096e \sim 0.0096,轨道倾角ε0.96\varepsilon \sim 0.96^\circ (Colpi et al. 2024)。
意味偏心率的高阶效应要比前一阶小100个量级,作为对比地球轨道偏心率0.0167\sim 0.0167。最后,取近心点为轨道最高点也有类似结论,此时对应卫星平面与黄道面60-60^\circ,偏心率和倾角略有区别,倾角更大,偏心率更小。

前面的分析限制了轨道形状R,eR, e及轨道倾角ε\varepsilon,同时远心点对应轨道最高点暗含了要求轨道近心点幅角(升交点与近心点夹角)为π2\frac{\pi}{2}(?)。在这些限制下,轨道具体指向仍可以绕黄道面法向旋转,对应于改变升交点经度,基准情况下暂取π2\frac{\pi}{2},相当于轨道面绕yy轴倾斜。前面计算倾角时取远心点对应最高点,以近心点为轨道参考方向,对应倾角ε-\varepsilon。下面具体分析卫星的轨道方程。在轨道平面内的卫星坐标可表示为:

x=R(cosψe)y=R1e2sinψ\begin{aligned} x &= R(\cos\psi - e)\\ y &= R\sqrt{1-e^2}\sin\psi \end{aligned}

其中ψ\psi是偏近点角。绕yy轴倾斜ε-\varepsilon,变换到黄道坐标系:

x=R(cosψe)cosεy=R1e2sinψz=R(cosψe)sinε\begin{aligned} x &= R(\cos\psi - e)\cos\varepsilon\\ y &= R\sqrt{1-e^2}\sin\psi\\ z &= -R(\cos\psi - e)\sin\varepsilon\\ \end{aligned}

为获得坐标随时间的演化,要建立偏近点角ψ\psi与时间的关系,这点可由开普勒方程得到:

ψ(t)esinψ(t)=ωˉ(tt0)+m0\psi(t) - e\sin\psi(t) = \bar{\omega} (t-t_0) + m_0

方程右侧为平近点角m(t)m(t),随时间均匀增加,ωˉ\bar{\omega}为轨道的平均角速度2π1yr\frac{2\pi}{1\rm yr}m0m_0t0t_0时刻的平近点角。注意这是个超越方程,没有解析解,需要数值求解。考虑到参考点也是匀圆周运动,因此平近点角与参考点相位直接对应,而初始平近点角与参考点初始位置对应。

最后,卫星实际轨道可由上述基准轨道绕zz旋转得到,三颗卫星轨道的升交点需相依次相错2π3\frac{2\pi}{3},同时相位(平近点角)也要依次差2π3\frac{2\pi}{3}(符号相反)。

(xkykzk)=(cosλksinλk0sinλkcosλk0001)(R(cosψke)cosεR1e2sinψkR(cosψke)sinε)λk=λ0+2π3(k1),    k=1,2,3ψkesinψk=ωˉ(tt0)+m02π3(k1)\boxed{\begin{aligned} \begin{pmatrix} x_k\\ y_k\\ z_k \end{pmatrix} &= \begin{pmatrix} \cos\lambda_k & -\sin\lambda_k & 0\\ \sin\lambda_k & \cos\lambda_k & 0\\ 0 & 0 & 1 \end{pmatrix} \begin{pmatrix} R(\cos\psi_k - e)\cos\varepsilon\\ R\sqrt{1-e^2}\sin\psi_k\\ -R(\cos\psi_k - e)\sin\varepsilon\\ \end{pmatrix}\\ \lambda_k &= \lambda_0 + \frac{2\pi}{3}(k-1), ~~~~ k = 1,2,3 \\ \psi_k &- e\sin\psi_k = \bar{\omega}(t-t_0) + m_0 - \frac{2\pi}{3}(k-1) \end{aligned}}

其中λ0\lambda_0为卫星1相对于基准轨道升交点经度的变化量,m0m_0为初始t0t_0时刻卫星1的轨道相位(平近点角),其余卫星的轨道方向及位置则完全由卫星1决定。

这里轨道形式上与Cornish & Rubbo 2003/Rubbo et al. 2004附录的公式不同,但两者等价。这里是从偏近点角出发推导,而相关文章是从真近点角出发推导。具体的,在轨道平面内由开普勒半径rr和真近点角γ\gamma可得到位置坐标,之后依次绕黄道坐标系yy轴和zz轴旋转,同时轨道倾角近似取3e\sqrt{3}e,即可得到相关文章中形式。

x=rcosγy=rsinγ   Ry(ε)   x=rcosγcosε  y=rsinγ          z=rcosγsinε   Rz(β)\begin{matrix} x = r\cos\gamma\\ y = r\sin\gamma \end{matrix} ~ ~ ~ \underrightarrow{\mathcal{R}_y(\varepsilon)} ~ ~ ~ \begin{matrix} x = r\cos\gamma \cos\varepsilon~~\\ y = r\sin\gamma~~~~~~~~~~\\ z = -r\cos\gamma \sin\varepsilon \end{matrix} ~ ~ ~ \underrightarrow{\mathcal{R}_z(\beta)}

这里开普勒半径r=R(1e2)1+ecosγr=\frac{R(1-e^2)}{1+e\cos\gamma}γ\gamma是真近点角,相关文中称之为方位角(azimuthal angle)或经度并不合适。其它符号上,β\beta对应上面λk\lambda_kλ\lambda对应λ0\lambda_0αβ\alpha-\beta对应轨道相位(平近点角),κ\kappa对应m0+λ0m_0 + \lambda_0

关于参考系
天球坐标系通常以北半球春分点(Spring/Vernal Equinox)为参考方向,对应黄道与天赤道交点。考虑到地轴的进动,赤道面是不断变动的,春分点也随之移动,为了建立固定的参考系,需要取特定时间的赤道面确定春分点。目前通常采用J2000.0纪元(epoch),即以2000年1月1日地球时12:00的赤道面为基准平面确定春分点。注意季节源于地轴倾斜,与地球轨道的近远日点无关!事实上,春分为每年的3月末,而近日点则在1月初(巧合,公历1月1日的选择源于历史而非天文)。
以太阳系质心(Solar System Barycentre, SSB)为原点,黄道面为参考平面,J2000.0纪元春分点为参考方向可建立黄道坐标系。注意与太阳系质心天球参考系(BCRS)/国际天球坐标系(ICRS)区分,虽然同样以SSB为原点,但后者属于赤道坐标系(但不随地球运动变化)。这里J2000.0纪元用于确定空间参考,而非时间起点。J2000.0本身属于儒略年计时,与之相对的有儒略日JDxxx、格里历(公历)、世界协调时UTC等计时标准。

计时通常采用太阳质心坐标时TCB,英文Barycentric Coordinate Time,缩写TCB来自法语Temps-coordonnée barycentrique。TCB属于质心天球参考系BCRS定义的一部分,但黄道坐标系也可以用。

轨道近似及臂长变化

x=R(cosψe)cosεy=R1e2sinψz=R(cosψe)sinε\begin{aligned} x &= R(\cos\psi - e)\cos\varepsilon\\ y &= R\sqrt{1-e^2}\sin\psi\\ z &= -R(\cos\psi - e)\sin\varepsilon\\ \end{aligned}

这里对上述基准轨道进行近似,轨道倾角和偏心率都展开到ι\iota一阶,有ε=3e\varepsilon=\sqrt{3}e。对于对于偏近点角ψ\psi,由开普勒方程:

ψ(t)esinψ(t)=m(t)\psi(t) - e\sin\psi(t) = m(t)

考虑到偏心率ee很小,迭代近似:

ψ0=m,   ψ1=m+esinψ0,   ψ2=m+esinψ1,...\psi_0 = m, ~~~ \psi_1 = m + e \sin \psi_0, ~~~ \psi_2 = m + e \sin \psi_1, ...

到偏心率的二阶有:

ψψ2=m+esin[m+esin(m)]=m+esin(m)cos[esin(m)]+ecos(m)sin[esin(m)]m+esin(m)+e2cos(m)sin(m)\begin{aligned} \psi &\approx \psi_2 = m + e \sin [m + e \sin (m)]\\ &= m + e \sin (m) \cos[e \sin (m)] + e \cos (m) \sin[e \sin (m)]\\ &\approx m + e \sin (m) + e^2 \cos (m)\sin (m)\\ \end{aligned}

由此,轨道方程中所有项都可展开到偏心率二阶:

sinψsin(m)+esin(m)cos(m)+e2sin(m)[cos2(m)12sin2(m)]cosψcos(m)esin2(m)32e2cos(m)sin2(m)cosε112ε2132e2sinεε3e1e2112e2\begin{aligned} \sin\psi &\approx \sin (m) + e \sin(m)\cos(m) + e^2 \sin(m)[\cos^2 (m) - \frac{1}{2}\sin^2 (m)]\\ \cos\psi &\approx \cos (m) - e \sin^2 (m) - \frac{3}{2}e^2 \cos(m)\sin^2 (m)\\ \cos\varepsilon &\approx 1 - \frac{1}{2}\varepsilon^2 \approx 1 - \frac{3}{2}e^2\\ \sin\varepsilon &\approx \varepsilon \approx \sqrt{3}e\\ \sqrt{1-e^2} &\approx 1-\frac{1}{2}e^2 \end{aligned}

带入轨道方程,保留到偏心率二阶项有:

x=Rcos(m)12Re[3cos(2m)]+18Re2[3cos(3m)10cos(m)5cos(m)]y=Rsin(m)+12Resin(2m)+18Re2[3sin(3m)10sin(m)+5sin(m)]z=3Recos(m)+3Re2[sin2(m)+1]\small\begin{aligned} x &= R\cos (m) - \frac{1}{2}Re [3- \cos(2m)] + \frac{1}{8}Re^2 [3\cos(3m) -10\cos(m) - 5\cos(m) ]\\ y &= R\sin (m) + \frac{1}{2}Re \sin(2m) + \frac{1}{8}Re^2 [3\sin(3m) - 10\sin(m) + 5\sin(m)]\\ z &= -\sqrt{3} R e \cos(m) + \sqrt{3} R e^2 [\sin^2(m) + 1]\\ \end{aligned}

考虑到mm对应αβ\alpha-\beta及绕zz轴旋转β\beta角,上述结果与Cornish & Rubbo 2003/Rubbo et al. 2004文中的一致。

近似到偏心1一阶,参考点轨道为黄道面内圆轨道,卫星臂长恒定为23Re2\sqrt{3}Re;近似到偏心率2阶,参考点仍为圆轨道,轨道面与黄道面平行,但纵向相对黄道面高332Re238L2R\frac{3\sqrt{3}}{2}Re^2\sim\frac{\sqrt{3}}{8}\frac{L^2}{R},。

探测响应调制

运动导致干涉臂指向(及臂长)变化会影响探测响应。对于地基观测,恒星级双黑洞观测时标t10st\sim 10s,地球自转约1°,该假设问题不大。但对于空间观测,超大质量双黑洞最终并和,观测时标在月的量级,卫星编队已有显著公转和自转,需要考虑探测器的运动。

投影部分,

臂长变化

L12(t)=L{1+e32[15sin(αλ+π6)cos(3α3λ)]}L13(t)=L{1+e32[15sin(αλπ6)cos(3α3λ)]}L23(t)=L{1e32[15cos(αλ)cos(3α3λ)]}\begin{aligned} L_{12}(t) &= L\left\{ 1 + \frac{e}{32}\left[15\sin(\alpha-\lambda+\frac{\pi}{6}) - \cos(3\alpha-3\lambda)\right] \right\}\\ L_{13}(t) &= L\left\{ 1 + \frac{e}{32}\left[-15\sin(\alpha-\lambda-\frac{\pi}{6}) - \cos(3\alpha-3\lambda)\right] \right\}\\ L_{23}(t) &= L\left\{ 1 - \frac{e}{32}\Big[15\cos(\alpha-\lambda) - \cos(3\alpha-3\lambda) \Big] \right\} \end{aligned}

L12(t)=L{1+e32[15sin(m+π6)cos(3m)]}L13(t)=L{1+e32[15sin(mπ6)cos(3m)]}L23(t)=L{1e32[15cos(m)cos(3m)]}m(t)=ωˉ(tt0)+m0\begin{aligned} L_{12}(t) &= L\left\{ 1 + \frac{e}{32}\left[15\sin(m +\frac{\pi}{6}) - \cos(3m)\right] \right\}\\ L_{13}(t) &= L\left\{ 1 + \frac{e}{32}\left[-15\sin(m -\frac{\pi}{6}) - \cos(3m)\right] \right\}\\ L_{23}(t) &= L\left\{ 1 - \frac{e}{32}\Big[15\cos(m) - \cos(3m) \Big] \right\}\\ m(t) &= \bar{\omega}(t-t_0) + m_0 \end{aligned}

探测响应

DP(t)=(x^12ax^12bx^13ax^13b)eabP\begin{aligned} D^P(t) = \left(\hat{x}^a_{12}\hat{x}^b_{12} - \hat{x}^a_{13}\hat{x}^b_{13}\right)e^P_{ab} \end{aligned}

sinΥ32,Υς\sin\Upsilon \leftarrow \frac{\sqrt{3}}{2}, \Upsilon\varsigma

D+(t)=sinΥ32{36sin2θ sin(2α2λ)+(3+cos2θ)×[cos2φ(9sin2λsin(4α2λ))+sin2φ(cos(4α2λ)9cos2λ)]43sin2θ[sin(3α2λφ)3sin(α2λ+φ)]}D×(t)=sinΥ8{cosθ [9cos(2λ2φ)cos(4α2λ2φ)]23sinθ [cos(3α2λφ)+3cos(α2λ+φ)]}\begin{aligned} \small D^+(t) = \frac{\sin\Upsilon}{32} \bigg\{& \small -36\sin^2\theta ~ \sin(2\alpha -2\lambda) + (3+\cos2\theta) \times \\ &\small \Big[ \cos2\varphi \Big(9\sin2\lambda - \sin(4\alpha-2\lambda) \Big) + \sin2\varphi \Big(\cos(4\alpha -2\lambda) - 9 \cos2\lambda\Big)\Big]\\ &\small - 4\sqrt{3} \sin2\theta\Big[\sin(3\alpha -2\lambda-\varphi) - 3\sin(\alpha - 2\lambda + \varphi)\Big] \bigg\} \\ \small D^\times(t) = \frac{\sin\Upsilon}{8} \bigg\{ & \small \cos\theta ~ \Big[9\cos(2\lambda-2\varphi) - \cos(4\alpha-2\lambda-2\varphi)\Big]\\ &\small - 2\sqrt{3}\sin\theta ~ \Big[\cos(3\alpha - 2\lambda-\varphi) + 3\cos(\alpha -2\lambda + \varphi)\Big]\bigg\} \end{aligned}

D+(t)=364{36sin2θ sin(2α2λ)+(3+cos2θ)×[cos2φ(9sin2λsin(4α2λ))+sin2φ(cos(4α2λ)9cos2λ)]43sin2θ[sin(3α2λφ)3sin(α2λ+φ)]}D×(t)=116{3cosθ [9cos(2λ2φ)cos(4α2λ2φ)]6sinθ [cos(3α2λφ)+3cos(α2λ+φ)]}\begin{aligned} \small D^+(t) = \frac{\sqrt{3}}{64} \bigg\{& \small -36\sin^2\theta ~ \sin(2\alpha -2\lambda) + (3+\cos2\theta) \times \\ &\small \Big[ \cos2\varphi \Big(9\sin2\lambda - \sin(4\alpha-2\lambda) \Big) + \sin2\varphi \Big(\cos(4\alpha -2\lambda) - 9 \cos2\lambda\Big)\Big]\\ &\small - 4\sqrt{3} \sin2\theta\Big[\sin(3\alpha -2\lambda-\varphi) - 3\sin(\alpha - 2\lambda + \varphi)\Big] \bigg\} \\ \small D^\times(t) = \frac{1}{16} \bigg\{ & \small \sqrt{3}\cos\theta ~ \Big[9\cos(2\lambda-2\varphi) - \cos(4\alpha-2\lambda-2\varphi)\Big]\\ &\small - 6\sin\theta ~ \Big[\cos(3\alpha - 2\lambda-\varphi) + 3\cos(\alpha -2\lambda + \varphi)\Big]\bigg\} \end{aligned}

其中 α=2πfmt+κ\alpha=2 \pi f_m t + \kappa 是随时间变化的,其余角度都是固定的,θ,ϕ\theta, \phi是波源方位角,κ,λ\kappa, \lambda定义了卫星编队的初始位置和指向

\begin{aligned} \end{aligned}

坐标变换

考虑到探测器轨道运动,在黄道参考系对系统进行描述更为方便。前面推导的系统转移函数为矢量(张量)形式,不依赖具体参考系。因此形式上,转移函数保持不变,只是此时x^,y^\hat{x}, \hat{y}都是随时间变化的。

虽然转移函数形式上保持不变,张量缩并、矢量点乘结果也都与坐标系选择无关,但实际计算时要将上述表达式展开为具体的方位角表示,是依赖坐标选择的,需要在不同坐标系的角度之间进行转换。

θS,ψS,θL,ψL\theta_S, \psi_S, \theta_L, \psi_L

波源的方位、波源轨道角动量

利用黄道坐标系(太阳系质心SSB参考系)建立引力波源与探测器的关系。
从波源到黄道坐标系:
从黄道坐标到探测器:

λ,β\lambda, \beta

经度、维度

信号调制
除了空间角度的不同,从探测器参考系到黄道坐标系,信号的到达时间也不同:

hdet(t)=hSSB(tk^rdetc)=hSSB(t)eiωk^rdetc\bm{h}_{\rm det}(t) = \bm{h}_{\rm \tiny SSB}\left(t - \hat{k}\cdot\frac{\vec{r}_{\rm det}}{c}\right) = \bm{h}_{\rm \tiny SSB}(t) e^{- i\omega\hat{k}\cdot\frac{\vec{r}_{\rm det}}{c}}

注意,虽然考虑到探测器运动,x^\hat{x}rdet\vec{r}_{\rm det}是含时的,但这并不影响积分结果。因为rdet\vec{r}_{\rm det}位于积分之外,而激光传播方向x^\hat{x}在激光发出后的传播过程中并不随探测器运动而改变,最终对于单向传播,只需要考虑时间延迟(调制)项即可。对于往返干涉情况,激光的反射及接收位置与发射位置的关系将受轨道运动影响:

12x^i(tLc)x^j(tLc)eijPeiω[1+k^x^(tLc)]L2c sinc{ω[1+k^x^(tLc)]L2c}eiωk^r2(tLc)Lx^(tLc)c+12x^i(t2Lc)x^j(t2Lc)eijPeiω[3+k^x^(t2Lc)]L2c sinc{ω[1k^x^(tL2c)]Lc}eiωk^r1(t2Lc)c\begin{aligned} &\small \frac{1}{2}\hat{x}^i(t-{\textstyle \frac{L}{c}}) \hat{x}^j(t-{\textstyle \frac{L}{c}}) e_{ij}^P e^{-i\omega [1 + \hat{k}\cdot\hat{x}(t-\frac{L}{c})]\frac{L}{2c}} ~ {\rm sinc}\left\{\omega [1 + \hat{k}\cdot\hat{x}(t-{\textstyle \frac{L}{c}})] {\textstyle \frac{L}{2c}}\right\} e^{-i\omega \hat{k} \cdot \frac{\vec{r}_2(t-\frac{L}{c}) - L\hat{x}(t-\frac{L}{c})}{c}}\\ &\small + \frac{1}{2}\hat{x}^i(t-{\textstyle \frac{2L}{c}}) \hat{x}^j(t-{\textstyle \frac{2L}{c}}) e_{ij}^P e^{-i\omega [3+\hat{k} \cdot \hat{x}(t-\frac{2L}{c})] \frac{L}{2c}} ~ {\rm sinc}\left\{\omega [1 - \hat{k}\cdot\hat{x}(t-{\textstyle \frac{L}{2c}})] {\textstyle \frac{L}{c}}\right\} e^{-i\omega \hat{k} \cdot \frac{\vec{r}_1(t-\frac{2L}{c})}{c}} \end{aligned}

由于卫星运动,激光指向其实要有一定提前量

x^ij(ti)rj(tj)ri(ti)ij(ti)ij(ti)rj(tj)ri(ti)=rj[ti+ij(ti)c]ri(ti)Lij(ti)rj(ti)ri(ti)\begin{aligned} \hat{x}_{ij}(t_i) &\equiv \frac{\vec{r}_j(t_j) - \vec{r}_i(t_i)}{\ell_{ij}(t_i)}\\ \ell_{ij}(t_i) &\equiv \|\vec{r}_j(t_j) - \vec{r}_i(t_i)\| = \left\|\vec{r}_j\left[t_i + {\textstyle \frac{\ell_{ij}(t_i)}{c}}\right] - \vec{r}_i(t_i)\right\|\\ L_{ij}(t_i) &\equiv \|\vec{r}_j(t_i) - \vec{r}_i(t_i)\| \end{aligned}

这里 ij\ell_{ij} 是从发射到接收的实际光程距离,注意与通常所说的臂长 LijL_{ij} 区分。将 ij,rj\ell_{ij}, \vec{r}_jvc\frac{v}{c}展开,有:

rj[ti+ij(ti)c]=rj(ti)+titi+ij(ti)cvj(τ)dτrj(ti)+ij(ti)cvj(ti)\vec{r}_j\left[t_i + {\textstyle \frac{\ell_{ij}(t_i)}{c}}\right] = \vec{r}_j(t_i) + \int_{t_i}^{t_i + \frac{\ell_{ij}(t_i)}{c}} \vec{v}_j(\tau) d\tau \approx \vec{r}_j(t_i) + {\textstyle \frac{\ell_{ij}(t_i)}{c}} \vec{v}_j(t_i)

ij(ti)=rj(ti)ri(ti)+ij(ti)cvj(ti)=\ell_{ij}(t_i) = \left\|\vec{r}_j(t_i) - \vec{r}_i(t_i) + {\textstyle \frac{\ell_{ij}(t_i)}{c}} \vec{v}_j(t_i)\right\| =

ij(ti)Lij1+[vj(ti)c]2Lij{1+12[vj(ti)c]2}\ell_{ij}(t_i) \approx L_{ij} \sqrt{1 + \left[\frac{v_j(t_i)}{c}\right]^2} \approx L_{ij} \left\{1 + \frac{1}{2}\left[\frac{v_j(t_i)}{c}\right]^2\right\}

ij(ti)=Lij(ti)[1+x^ij(ti)vj(ti)c+o(v2c)]\ell_{ij}(t_i) = L_{ij}(t_i)\left[1 + \hat{x}_{ij}(t_i)\cdot\frac{\vec{v}_j(t_i)}{c} + o\left(\frac{v^2}{c}\right)\right]

考虑卫星运动的提前量导致的光程变化在vc1c2πfmR104\frac{v}{c} \approx \frac{1}{c}2\pi f_m R \approx 10^-4量级,作为对比卫星轨道运动本身所导致的臂长变化在1/1001/100量级。对于5×109m5\times10^9 \rm m的平均臂长,前者为100km100\rm km,后者为104km10^4\rm km。如果从角度上考虑,这个提前量对应多大角度?

对于目前的空间探测,激光沿干涉臂传播时间仅约1010秒,暂时忽略该过程中探测器位置及干涉臂方向的改变,最终唯一的区别是多了个时间延迟(调制)项。

FSSBP(f,t)=FP(f,t)ei2πfk^rdet(t)c\boxed{ F_{\rm \tiny SSB}^P(f, t) = F^P(f, t) e^{- i 2\pi f\hat{k}\cdot\frac{\vec{r}_{\rm det}(t)}{c}} }

注意,这里的响应调制项是频率依赖的。

刚性绝热近似 Rigid Adiabatic Approximation cf/f˙Lc f/\dot{f} \ll L
对比 长波近似 f<c2πLf < \frac{c}{2\pi L}

时间延迟干涉

通常讨论TDI是针对外差相位(距离变化)测量的,实际中到底采用相位(距离),还是频率(速度)或者啁啾度(加速度)还是一个开放问题,这涉及相位计的内部设计和处理流程,包括滤波、量化等过程的影响。在臂长(相等?)不变的情况,不同物理量的TDI可简单实现转换,当考虑到臂长变化不同物理量的TDI需要具体讨论。Bayle et al. 2021中就给除了频率的TDI,Spadaro et al. 2023提出了加速度TDI,但是在不考虑臂长变化的前提下的,可以扩展。

相位、相对频移、

暂时不考虑多普勒效应,对于分光干涉系统

似乎有两种思路:

一种是双臂干涉直接对应X、Y、Z通道TDI的响应
另一种是利用单臂响应从TDI推导出X、Y、Z通道
如果考虑到不等臂肯定是后者,如果是等臂的两者是不是就是等价的。

在臂长相等情况下,原则上并不需要TDI,但也可以按TDI进行组合,此时时间延迟算子是同一个,可以等效视为0代TDI
在臂长不同情况下,此时时间延迟算子不同的,为通常所谓的1代TDI
在臂长变化情况下,此时时间延迟算是怎么处理,含时吗?,为通常所谓的2代TDI
通常在科学分析时,并不真正考虑激光频率噪声,只是用TDI组合
根据Hartwig & Muratore 2021的讨论,常见的0代TDI与1代TDI差别?!

相应TDI组合的噪声曲线?

以激光接收端为参考,x^\hat{x} 反向(指向发射端)

sij(f,t)=12x^ax^bh~ab(f)eiω(trik^ric)1iω(1k^x^ij)[1eiω(1k^x^ij)Lc]=12x^ax^bh~ab(f)eiω(trik^ric)eiω(1k^x^ij)L2c Lc sinc[ω(1k^x^ij)L2c]Δνν=x^ax^b2(1+k^x^)h~ab(f)eiω(trik^ric)[1eiω(1k^x^)Lc]=x^ax^b2(1k^x^)[h~ab(f)eiω(trik^ric)h~ab(f)eiω(trjk^rjc)]\begin{aligned} s_{ij} (f, t) &= \frac{1}{2}\hat{x}^a\hat{x}^b \tilde{h}_{ab}(f) e^{i\omega (t|_{\vec{r}_i} - \hat{k}\cdot\frac{\vec{r}_i}{c})} \small \frac{1}{i\omega (1 - \hat{k}\cdot\hat{x}_{ij})} \left[1 - e^{-i\omega (1 - \hat{k}\cdot\hat{x}_{ij})\frac{L}{c}}\right]\\ &= \frac{1}{2}\hat{x}^a\hat{x}^b \tilde{h}_{ab}(f) e^{i\omega (t|_{\vec{r}_i} - \hat{k}\cdot\frac{\vec{r}_i}{c})} e^{-i\omega (1- \hat{k} \cdot \hat{x}_{ij}) \frac{L}{2c} } ~ \small \frac{L}{c} ~ {\rm sinc}\left[\omega(1-\hat{k} \cdot \hat{x}_{ij}) {\textstyle \frac{L}{2c}}\right]\\ \frac{\Delta \nu}{\nu} &= \frac{\hat{x}^a\hat{x}^b}{2 (1 + \hat{k}\cdot\hat{x})} \tilde{h}_{ab}(f) e^{i\omega (t|_{\vec{r}_i} - \hat{k}\cdot\frac{\vec{r}_i}{c})} \left[1 - e^{-i\omega (1 - \hat{k}\cdot\hat{x})\frac{L}{c}}\right]\\ & = \frac{\hat{x}^a\hat{x}^b}{2(1 - \hat{k}\cdot\hat{x})} \left[\tilde{h}_{ab}(f) e^{i\omega (t|_{\vec{r}_i} - \hat{k}\cdot\frac{\vec{r}_i}{c})} - \tilde{h}_{ab}(f) e^{i\omega (t|_{\vec{r}_j} - \hat{k}\cdot\frac{\vec{r}_j}{c}) } \right] \end{aligned}

s~ij(f)=12x^ax^bh~ab(f)eiωk^ric1iω(1k^x^ij)[1eiω(1k^x^ij)Lc]=12x^ax^bh~ab(f)eiωk^riceiω(1k^x^ij)L2c Lc sinc[ω(1k^x^ij)L2c]y~ij(f)=x^ax^b2(1+k^x^)h~ab(f)eiωk^ric[1eiω(1k^x^)Lc]=x^ax^b2(1k^x^)h~ab(f)[eiωk^riceiω(Lc+k^rjc)]\begin{aligned} \tilde{s}_{ij} (f) &= \frac{1}{2}\hat{x}^a\hat{x}^b \tilde{h}_{ab}(f) e^{-i\omega \hat{k}\cdot\frac{\vec{r}_i}{c}} \small \frac{1}{i\omega (1 - \hat{k}\cdot\hat{x}_{ij})} \left[1 - e^{-i\omega (1 - \hat{k}\cdot\hat{x}_{ij})\frac{L}{c}}\right]\\ &= \frac{1}{2}\hat{x}^a\hat{x}^b \tilde{h}_{ab}(f) e^{-i\omega \hat{k}\cdot\frac{\vec{r}_i}{c}} e^{-i\omega (1- \hat{k} \cdot \hat{x}_{ij}) \frac{L}{2c} } ~ \small \frac{L}{c} ~ {\rm sinc}\left[\omega(1-\hat{k} \cdot \hat{x}_{ij}) {\textstyle \frac{L}{2c}}\right]\\ \tilde{y}_{ij}(f) &= \frac{\hat{x}^a\hat{x}^b}{2 (1 + \hat{k}\cdot\hat{x})} \tilde{h}_{ab}(f) e^{-i\omega \hat{k}\cdot\frac{\vec{r}_i}{c}} \left[1 - e^{-i\omega (1 - \hat{k}\cdot\hat{x})\frac{L}{c}}\right]\\ & = \frac{\hat{x}^a\hat{x}^b}{2(1 - \hat{k}\cdot\hat{x})} \tilde{h}_{ab}(f) \left[e^{-i\omega \hat{k}\cdot\frac{\vec{r}_i}{c}} - e^{-i\omega ( \frac{L}{c} + \hat{k}\cdot\frac{\vec{r}_j}{c}) } \right] \end{aligned}

X=[D31D13D21η12+D31D13η21+D31η13+η31][η21+D21η12+D21D12η31+D21D12D31η13]=(D31D13D21D21)η12+(D31D131)η21+(D31D21D12D31)η13+(1D21D12)η31\begin{aligned} X = & [\mathcal{D}_{31} \mathcal{D}_{13} \mathcal{D}_{21} \eta_{12} +\mathcal{D}_{31} \mathcal{D}_{13} \eta_{21} +\mathcal{D}_{31} \eta_{13} + \eta_{31}]\\ & - [\eta_{21} +\mathcal{D}_{21} \eta_{12} +\mathcal{D}_{21} \mathcal{D}_{12} \eta_{31} +\mathcal{D}_{21} \mathcal{D}_{12} \mathcal{D}_{31} \eta_{13}]\\ = & (\mathcal{D}_{31} \mathcal{D}_{13} \mathcal{D}_{21} -\mathcal{D}_{21}) \eta_{12} + (\mathcal{D}_{31} \mathcal{D}_{13} - 1) \eta_{21} \\ & + (\mathcal{D}_{31} -\mathcal{D}_{21} \mathcal{D}_{12} \mathcal{D}_{31}) \eta_{13} + (1 -\mathcal{D}_{21} \mathcal{D}_{12}) \eta_{31} \end{aligned}

Dij\mathcal{D}_{ij} 是时间延迟算子(时域)
Δij\Delta_{ij} 是相位延迟因子(频域),Δij=eiωLijc\Delta_{ij} = e^{-i\omega\frac{L_{ij}}{c}} v.s. Δij=eiωLijc\Delta_{ij} = e^{i\omega\frac{L_{ij}}{c}}

s~X(f)=(Δ21Δ13Δ31Δ21)s~12+(Δ13Δ311)s~21+(Δ31Δ31Δ12Δ21)s~13+(1Δ12Δ21)s~31y~X(f)=(Δ21Δ13Δ31Δ21)y~12+(Δ13Δ311)y~21+(Δ31Δ31Δ12Δ21)y~13+(1Δ12Δ21)y~31\begin{aligned} \tilde{s}_X(f) =& (\Delta_{21} \Delta_{13} \Delta_{31} - \Delta_{21}) \tilde{s}_{12} + (\Delta_{13} \Delta_{31} - 1) \tilde{s}_{21}\\ & + (\Delta_{31} - \Delta_{31} \Delta_{12} \Delta_{21})\tilde{s}_{13} + (1 - \Delta_{12} \Delta_{21}) \tilde{s}_{31}\\ \tilde{y}_X(f) =& (\Delta_{21} \Delta_{13} \Delta_{31} - \Delta_{21}) \tilde{y}_{12} + (\Delta_{13} \Delta_{31} - 1) \tilde{y}_{21}\\ & + (\Delta_{31} - \Delta_{31} \Delta_{12} \Delta_{21})\tilde{y}_{13} + (1 - \Delta_{12} \Delta_{21}) \tilde{y}_{31} \end{aligned}

FijP(f)=12x^ax^beabPeiωk^ric1iω(1k^x^ij)[1eiω(1k^x^ij)Lc]=12x^ax^bh~ab(f)eiωk^riceiω(1k^x^ij)L2c Lc sinc[ω(1k^x^ij)L2c]FijP(f)=x^ax^b2(1+k^x^)h~ab(f)eiωk^ric[1eiω(1k^x^)Lc]=x^ax^beabP2(1k^x^)[eiωk^riceiω(Lc+k^rjc)]\begin{aligned} F^P_{ij} (f) &= \frac{1}{2}\hat{x}^a\hat{x}^b e^P_{ab} e^{-i\omega \hat{k}\cdot\frac{\vec{r}_i}{c}} \small \frac{1}{i\omega (1 - \hat{k}\cdot\hat{x}_{ij})} \left[1 - e^{-i\omega (1 - \hat{k}\cdot\hat{x}_{ij})\frac{L}{c}}\right]\\ &= \frac{1}{2}\hat{x}^a\hat{x}^b \tilde{h}_{ab}(f) e^{-i\omega \hat{k}\cdot\frac{\vec{r}_i}{c}} e^{-i\omega (1- \hat{k} \cdot \hat{x}_{ij}) \frac{L}{2c} } ~ \small \frac{L}{c} ~ {\rm sinc}\left[\omega(1-\hat{k} \cdot \hat{x}_{ij}) {\textstyle \frac{L}{2c}}\right]\\ F^P_{ij}(f) &= \frac{\hat{x}^a\hat{x}^b}{2 (1 + \hat{k}\cdot\hat{x})} \tilde{h}_{ab}(f) e^{-i\omega \hat{k}\cdot\frac{\vec{r}_i}{c}} \left[1 - e^{-i\omega (1 - \hat{k}\cdot\hat{x})\frac{L}{c}}\right]\\ & = \frac{\hat{x}^a\hat{x}^b e^P_{ab}}{2(1 - \hat{k}\cdot\hat{x})} \left[ e^{-i\omega \hat{k}\cdot\frac{\vec{r}_i}{c}} - e^{-i\omega ( \frac{L}{c} + \hat{k}\cdot\frac{\vec{r}_j}{c}) } \right] \end{aligned}

FXP(f)=(Δ21Δ13Δ31Δ21)F12P+(Δ13Δ311)F21P+(Δ31Δ31Δ12Δ21)F13P+(1Δ12Δ21)F31P\begin{aligned} F^P_X(f) =& (\Delta_{21} \Delta_{13} \Delta_{31} - \Delta_{21}) F^P_{12} + (\Delta_{13} \Delta_{31} - 1) F^P_{21}\\ & + (\Delta_{31} - \Delta_{31} \Delta_{12} \Delta_{21})F^P_{13} + (1 - \Delta_{12} \Delta_{21}) F^P_{31} \end{aligned}

Seto 2002
η12,η13\eta_{12}, \eta_{13} 是单臂单向的响应
η31+Δ31η13,η21+Δ21η12\eta_{31} + \Delta_{31}\eta_{13}, \eta_{21} + \Delta_{21}\eta_{12} 是单臂双向的响应
η31+Δ31η13(η21+Δ21η12)\eta_{31} + \Delta_{31}\eta_{13} - (\eta_{21} + \Delta_{21}\eta_{12}) 是等臂迈克尔逊的响应
这里假设臂长相等,所有的时间延迟算符作用效果相等,相位延迟因子均为Δ=eiωLc=eix\Delta = e^{-i\omega\frac{L}{c}}= e^{-i x}。此时TDI 1.0及1.5的X通道数据为

[η31+Δ31η13(η21+Δ21η12)]ΔΔ[η31+Δ31η13(η21+Δ21η12)]=(1ei2x)[η31+Δ31η13(η21+Δ21η12)]\begin{aligned} &\Big[\eta_{31} + \Delta_{31}\eta_{13} - (\eta_{21} + \Delta_{21}\eta_{12})\Big] - \Delta\Delta\Big[\eta_{31} + \Delta_{31}\eta_{13} - (\eta_{21} + \Delta_{21}\eta_{12})\Big]\\ &= \left(1-e^{-i 2 x}\right)\Big[\eta_{31} + \Delta_{31}\eta_{13} - (\eta_{21} + \Delta_{21}\eta_{12})\Big] \end{aligned}

(1ei2x)2 12[r^12ar^12br^13ar^13b]h~ab(f)\left(1-e^{-i 2 x}\right)^2 ~ \frac{1}{2}\left[\hat{r}^a_{12}\hat{r}^b_{12} - \hat{r}^a_{13}\hat{r}^b_{13}\right] \tilde{h}_{ab}(f)

其中(1ei2x)2(1-e^{-i 2 x})^2为TDI所引入的额外因子,进一步展开4x2\sim 4 x^2

1ei2x=eix(eixeix)=2 i eixsinx1-e^{-i 2 x} = e^{-i x}\left(e^{i x} - e^{-i x}\right) = 2 ~i~ e^{-i x} \sin x

4ei2xsin2x r^12ar^12br^13ar^13b2h~ab(f)-4 e^{-i 2x} \sin^2 x ~ \frac{\hat{r}^a_{12}\hat{r}^b_{12} - \hat{r}^a_{13}\hat{r}^b_{13}}{2} \tilde{h}_{ab}(f)

ei2xsin2xeixxsinxsin2xxsinxx2e^{-i 2x} \sin^2 x \sim e^{-i x} x \sin x \sim \sin^2 x \sim x \sin x \sim x^2

Cornish & Shuman 2020
长波近似,忽略k^\hat{k}项目,记xxxx
(f=c/2πL,x=2πfL/cf_* = c/2 \pi L, x = 2\pi f L/c)

T(f)=sinc(x/2) ei3x/2+sinc(x/2) eix/2=sinc(x/2)eix(eix/2+eix/2)=sin(x/2)x/2eix 2cos(x/2)=2sinxxeix     2\begin{aligned} \mathcal{T}(f) &= {\rm sinc} (x/2)~ e^{-i3x/2} + {\rm sinc}(x/2)~e^{-ix/2}\\ & = {\rm sinc}(x/2) e^{-ix} \left(e^{-ix/2} + e^{ix/2} \right)\\ & = \frac{\sin(x/2)}{x/2} e^{-ix} ~ 2\cos(x/2)\\ & = 2 \frac{\sin x}{x} e^{-ix} ~~~~~ {\color{red} \rightarrow 2}\\ \end{aligned}

X(f)=x eixsinx[r^12ar^12br^13ar^13b]h~ab(f) T(f)=x eixsinx[r^12ar^12br^13ar^13b]h~ab(f) 2sinxxeix=2ei2xsin2x[r^12ar^12br^13ar^13b]h~ab(f)\begin{aligned} X(f) &= -x ~ e^{-ix} \sin x \left[\hat{r}^a_{12}\hat{r}^b_{12} - \hat{r}^a_{13}\hat{r}^b_{13}\right]\tilde{h}_{ab}(f) ~ \mathcal{T}(f) \\ &= -x ~ e^{-ix} \sin x \left[\hat{r}^a_{12}\hat{r}^b_{12} - \hat{r}^a_{13}\hat{r}^b_{13}\right]\tilde{h}_{ab}(f) ~ 2 \frac{\sin x}{x} e^{-ix} \\ &= - 2 e^{-i 2 x} \sin^2 x \left[\hat{r}^a_{12}\hat{r}^b_{12} - \hat{r}^a_{13}\hat{r}^b_{13}\right]\tilde{h}_{ab}(f) \\ \end{aligned}

X(f)2sin2x[r^12ar^12br^13ar^13b]h~ab(f)2xsinx[r^12ar^12br^13ar^13b]h~ab(f)4xsinx32x^ax^by^ay^b2h~ab(f)\begin{aligned} X(f) &\approx - 2 \sin^2 x \left[\hat{r}^a_{12}\hat{r}^b_{12} - \hat{r}^a_{13}\hat{r}^b_{13}\right]\tilde{h}_{ab}(f) \\ &\approx - 2 x \sin x \left[\hat{r}^a_{12}\hat{r}^b_{12} - \hat{r}^a_{13}\hat{r}^b_{13}\right]\tilde{h}_{ab}(f) \\ &\approx 4 x \sin x \frac{\sqrt{3}}{2} \frac{\hat{x}^a\hat{x}^b - \hat{y}^a\hat{y}^b}{2}\tilde{h}_{ab}(f) \\ \end{aligned}

Babak et al. 2021

X(f)4ωLsin(ωL)32(F+h~++F×h~×)\begin{aligned} X(f) &\approx 4\omega L \sin(\omega L) \frac{\sqrt{3}}{2}(F_+\tilde{h}_+ + F_\times\tilde{h}_\times)\\ \end{aligned}

与上面Cornish & Shuman 2020的结果对比,F+F_+对应直角迈克尔逊天线模式。 对应到E通道

E(f)2xsinx32x^ax^by^ay^b2h~ab(f)\begin{aligned} E(f) &\approx 2x \sin x \frac{3}{\sqrt{2}} \frac{\hat{x}^a\hat{x}^b - \hat{y}^a\hat{y}^b}{2}\tilde{h}_{ab}(f)\\ \end{aligned}

Seto 2006 (少了2倍?)

E32x2x^ax^by^ay^b2h~ab(f)\begin{aligned} E &\approx \frac{3}{\sqrt{2}} x^2 \frac{\hat{x}^a\hat{x}^b - \hat{y}^a\hat{y}^b}{2}\tilde{h}_{ab}(f) \\ \end{aligned}

Krolak et al. 2004

LDC(2a), Vallisneri & Galley 2012

SX=16sin2x[SOMS+2(1+cos2x)Sacc],   SXY=8sin2xcosx(SOMS+4Sacc)S_X = 16 \sin^2 x\left[ S_{\rm OMS} + 2(1+\cos^2 x)S_{\rm acc}\right], ~~~ S_{XY} = -8\sin^2 x \cos x \left( {S_{\rm OMS} + 4 S_{\rm acc}} \right)

A=ZX2,   E=X2Y+Z6,   T=X+Y+Z3A = \frac{Z - X}{\sqrt{2}}, ~~~ E = \frac{X - 2Y + Z}{\sqrt{6}}, ~~~ T = \frac{X + Y + Z}{\sqrt{3}}

最早出现于Prince et al. 2002

SA=12(SX+SZSXZSZX)=SXSXYS_A = \frac{1}{2}\left(S_X + S_Z - S_{XZ} - S_{ZX}\right) = S_X - S_{XY}

SE=SXSXY,      ST=SX+2SXYS_E = S_X - S_{XY}, ~~~~~~ S_T = S_X + 2S_{XY}

最终有:

SA=SE=16sin2x[(1+12cosx)SOMS+2(1+cosx+cos2x)Sacc]S_A = S_E = 16\sin^2 x\left[(1+\frac{1}{2}\cos x)S_{\rm OMS} + 2(1+\cos x + \cos^2 x)S_{\rm acc}\right]

ST=16sin2x[(1cosx)SOMS+2(1cosx)2Sacc]S_T = 16\sin^2 x\left[(1-\cos x)S_{\rm OMS} + 2(1-\cos x)^2 S_{\rm acc}\right]

对比Cornish 2002中的噪声谱

Sn=4[SOMS+2(1+cos2x)Sacc]S_n = 4\left[S_{\rm OMS} + 2(1+\cos^2 x)S_{\rm acc}\right]

少了因子(2sinx)2(2\sin x)^2,因为这里是等臂迈克尔逊噪声谱,相比X通道少了TDI因子1ei2x2sinx|1-e^{-i2 x}| \sim 2 \sin x
上述噪声谱是针对光程变化测量,Robson et al. 2019则将其对应到应力:

Sn=1L2[SOMS+2(1+cos2x)Sacc]S_n = \frac{1}{L^2}\left[S_{\rm OMS} + 2(1+\cos^2 x)S_{\rm acc}\right]

如果将AE通道的光程变化噪声谱对应到应力,有:

Sn=sin2xL2[(1+12cosx)SOMS+2(1+cosx+cos2x)Sacc]S_n = \frac{\sin^2 x}{L^2}\left[\left(1+\frac{1}{2}\cos x\right)S_{\rm OMS} + 2(1+\cos x + \cos^2 x)S_{\rm acc}\right]

Cornish & Shuman (2020)

hab=h+ϵab++h×ϵab×h_{ab} = h_+\epsilon_{ab}^+ + h_\times \epsilon_{ab}^\times

h+=h01+cos2ι2cosΦ,   h×=h0cosι sinΦh_+ = h_0 \frac{1+\cos^2\iota}{2} \cos \Phi, ~~~ h_\times = h_0 \cos\iota ~ \sin\Phi

h+=h01+cos2ι2eiΦ,   h×=h0cosι ei(Φπ2)h_+ = h_0 \frac{1+\cos^2\iota}{2} e^{i\Phi}, ~~~ h_\times = h_0 \cos\iota ~ e^{i(\Phi-\frac{\pi}{2})}

hab=h0eiΦ[1+cos2ι2ϵab+icosι ϵab×]=h(t)[A+ϵab++iA×ϵab×]h_{ab} = h_0 e^{i\Phi} \left[ \frac{1+\cos^2\iota}{2} \epsilon_{ab}^+ - i \cos\iota ~ \epsilon_{ab}^\times\right] = h(t) \left[ A_+ \epsilon_{ab}^+ + i A_\times \epsilon_{ab}^\times\right]

δνν0\frac{\delta \nu}{\nu_0}为观测量转换到以δLL\frac{\delta L}{L}为观测测量时只需要移除4f/f=4x4f/f_* = 4x因子。
长波近似

X(f)=2xsinx[r^12ar^12br^13ar^13b]h~ab(f)X(f) = - 2x \sin x \left[\hat{r}^a_{12}\hat{r}^b_{12} - \hat{r}^a_{13}\hat{r}^b_{13}\right] \tilde{h}_{ab}(f)

Y(f)=2xsinx[r^23ar^23br^21ar^21b]h~ab(f)Y(f) = - 2x \sin x \left[\hat{r}^a_{23}\hat{r}^b_{23} - \hat{r}^a_{21}\hat{r}^b_{21}\right] \tilde{h}_{ab}(f)

Z(f)=2xsinx[r^31ar^31br^32ar^32b]h~ab(f)Z(f) = - 2x\sin x \left[\hat{r}^a_{31}\hat{r}^b_{31} - \hat{r}^a_{32}\hat{r}^b_{32}\right] \tilde{h}_{ab}(f)

ZX=2xsinx[r^31ar^31br^32ar^32br^12ar^12b+r^13ar^13b]h~ab(f)=2xsinx 32 [y^ay^bx^ax^b]h~ab(f)X2Y+Z=2xsinx[r^12ar^12br^13ar^13b2r^23ar^23b+2r^21ar^21b+r^31ar^31br^32ar^32b]h~ab(f)=2xsinx 3 [r^21ar^21br^23ar^23b]h~ab(f)X+Y+Z=2xsinx[r^12ar^12br^13ar^13b+r^23ar^23br^21ar^21b+r^31ar^31br^32ar^32b]h~ab(f)=0\footnotesize \begin{aligned} Z - X &= - 2x\sin x \left[\hat{r}^a_{31}\hat{r}^b_{31} - \hat{r}^a_{32}\hat{r}^b_{32} - \hat{r}^a_{12}\hat{r}^b_{12} + \hat{r}^a_{13}\hat{r}^b_{13}\right] \tilde{h}_{ab}(f)\\ &= - 2x\sin x ~\frac{3}{2}~\left[\hat{y}^a\hat{y}^b - \hat{x}^a\hat{x}^b\right] \tilde{h}_{ab}(f)\\ X - 2Y + Z &= - 2x\sin x \left[ \hat{r}^a_{12}\hat{r}^b_{12} - \hat{r}^a_{13}\hat{r}^b_{13} - 2\hat{r}^a_{23}\hat{r}^b_{23} + 2\hat{r}^a_{21}\hat{r}^b_{21} + \hat{r}^a_{31}\hat{r}^b_{31} - \hat{r}^a_{32}\hat{r}^b_{32}\right] \tilde{h}_{ab}(f)\\ &= - 2x\sin x ~3~\left[\hat{r}^a_{21}\hat{r}^b_{21} - \hat{r}^a_{23}\hat{r}^b_{23}\right] \tilde{h}_{ab}(f)\\ X + Y + Z =& -2x\sin x \left[ \hat{r}^a_{12}\hat{r}^b_{12} - \hat{r}^a_{13}\hat{r}^b_{13} + \hat{r}^a_{23}\hat{r}^b_{23} - \hat{r}^a_{21}\hat{r}^b_{21} + \hat{r}^a_{31}\hat{r}^b_{31} - \hat{r}^a_{32}\hat{r}^b_{32}\right] \tilde{h}_{ab}(f) = 0 \end{aligned}

ZX2=32 xsinxx^ax^by^ay^b2h~ab(f)X2Y+Z6=32 xsinxx^ay^b+y^ax^b2h~ab(f)SA=SE=16sin2x[(1+12cosx)SOMS+2(1+cosx+cos2x)Sacc]\begin{aligned} \frac{Z - X}{\sqrt{2}} &= 3\sqrt{2} ~ x\sin x \frac{\hat{x}^a\hat{x}^b - \hat{y}^a\hat{y}^b}{2} \tilde{h}_{ab}(f)\\ \frac{X - 2Y + Z}{\sqrt{6}} &= 3\sqrt{2} ~ x\sin x \frac{\hat{x}^a\hat{y}^b + \hat{y}^a\hat{x}^b}{2} \tilde{h}_{ab}(f) \\ S_A = S_E &= 16\sin^2 x\left[(1+\frac{1}{2}\cos x)S_{\rm OMS} + 2(1+\cos x + \cos^2 x)S_{\rm acc}\right] \end{aligned}

这种归一化下,信号放大了约3/21.2\sim \sqrt{3/2} \sim 1.2倍。取噪声谱中x=0x=0可知,低频噪声相应放大为原先的3/23/2倍。调整归一化因子可使信号及噪声,在长波极限下,都保持与XYZ通道相同的水平。

E=ZY3,   A=2XYZ3,   T=X+Y+Z3E = \frac{Z - Y}{\sqrt{3}}, ~~~ A = \frac{2X - Y - Z}{3}, ~~~T = \frac{X + Y + Z}{3}

最早出现于Adams & Cornish 2010

X(f)=2xsinx[r^12ar^12br^13ar^13b]h~ab(f)=2xsinxsinα[x^ay^b+y^ax^b]h~ab(f)=23 xsinxx^ax^by^ay^b2h~ab(f)\begin{aligned} X(f) &= - 2x\sin x \left[\hat{r}^a_{12}\hat{r}^b_{12} - \hat{r}^a_{13}\hat{r}^b_{13}\right] \tilde{h}_{ab}(f)\\ % = 2x\sin x 2\sin\frac{\alpha}{2}\cos\frac{\alpha}{2} \left[\hat{x}^a\hat{y}^b + \hat{y}^a\hat{x}^b\right] \tilde{h}_{ab}(f)\\ &= 2x\sin x \sin \alpha \left[\hat{x}^a\hat{y}^b + \hat{y}^a\hat{x}^b\right] \tilde{h}_{ab}(f)\\ &= 2\sqrt{3}~ x\sin x \frac{\hat{x}'^a\hat{x}'^b - \hat{y}'^a\hat{y}'^b}{2} \tilde{h}_{ab}(f)\\ \end{aligned}

XY=2xsinx[r^31ar^31br^32ar^32br^23ar^23b+r^21ar^21b]h~ab(f)=2xsinx 32 [y^ay^bx^ax^b]h~ab(f)2XYZ=2xsinx[2r^12ar^12b2r^13ar^13br^23ar^23b+r^21ar^21br^31ar^31b+r^32ar^32b]h~ab(f)=2xsinx 3 [r^12ar^12br^13ar^13b]h~ab(f)=2xsinx332[x^ay^b+y^ax^b]h~ab(f)\footnotesize \begin{aligned} X - Y &= - 2x\sin x \left[\hat{r}^a_{31}\hat{r}^b_{31} - \hat{r}^a_{32}\hat{r}^b_{32} - \hat{r}^a_{23}\hat{r}^b_{23} + \hat{r}^a_{21}\hat{r}^b_{21}\right] \tilde{h}_{ab}(f)\\ &= - 2x\sin x ~\frac{3}{2}~\left[\hat{y}^a\hat{y}^b - \hat{x}^a\hat{x}^b\right] \tilde{h}_{ab}(f)\\ 2X - Y - Z &= - 2x\sin x \left[ 2\hat{r}^a_{12}\hat{r}^b_{12} - 2\hat{r}^a_{13}\hat{r}^b_{13} - \hat{r}^a_{23}\hat{r}^b_{23} + \hat{r}^a_{21}\hat{r}^b_{21}- \hat{r}^a_{31}\hat{r}^b_{31} + \hat{r}^a_{32}\hat{r}^b_{32}\right] \tilde{h}_{ab}(f)\\ &= - 2x\sin x ~3~\left[\hat{r}^a_{12}\hat{r}^b_{12} - \hat{r}^a_{13}\hat{r}^b_{13}\right] \tilde{h}_{ab}(f) = 2x\sin x \frac{3\sqrt{3}}{2} \left[\hat{x}^a\hat{y}^b + \hat{y}^a\hat{x}^b\right] \tilde{h}_{ab}(f)\\ \end{aligned}

ZY3=23 xsinxx^ax^by^ay^b2h~ab(f)2XYZ3=23 xsinxx^ay^b+y^ax^b2h~ab(f)SA=SE=23×16sin2x[(1+12cosx)SOMS+2(1+cosx+cos2x)Sacc]\begin{aligned} \frac{Z - Y}{\sqrt{3}} &= 2\sqrt{3}~ x\sin x \frac{\hat{x}^a\hat{x}^b - \hat{y}^a\hat{y}^b}{2} \tilde{h}_{ab}(f)\\ \frac{2X - Y - Z}{3} &= 2\sqrt{3}~ x\sin x \frac{\hat{x}^a\hat{y}^b + \hat{y}^a\hat{x}^b}{2} \tilde{h}_{ab}(f) \\ S_A = S_E &= \frac{2}{3} \times 16\sin^2 x\left[(1+\frac{1}{2}\cos x)S_{\rm OMS} + 2(1+\cos x + \cos^2 x)S_{\rm acc}\right] \end{aligned}

Babak et al. 2021
这里应该是分数频移响应

X~1.5=4ffsinx32[F+h~++F×h~×]\tilde{X}_{1.5} = 4 \frac{f}{f_*} \sin x \frac{\sqrt{3}}{2} [F_+\tilde{h}_+ + F_\times \tilde{h}_\times]

变为应力响应(除以xx还是4x4x??)

X~1.5=sinx32[F+h~++F×h~×]\tilde{X}_{1.5} = \sin x \frac{\sqrt{3}}{2} [F_+\tilde{h}_+ + F_\times \tilde{h}_\times]

Sh,X1.5LW(f)=103Sn,X1.5(f)16sin2x x2S^{LW}_{h,X_{1.5}} (f) = \frac{10}{3} \frac{S_{n,X_{1.5}}(f)}{16\sin^2x ~ x^2}

考虑到长波近似下X+Y+Z=0X + Y + Z=0,在Cornish报告中还见到另一种形式的A, E, T,但找不到文章出处。

A=32X,    E=12(X+2Y),    T=13(X+Y+Z)A = \frac{\sqrt{3}}{2} X, ~~~~ E = \frac{1}{2}(X + 2Y), ~~~~ T = \frac{1}{3}(X + Y + Z)

[1+exp(iω 2L/c)][1+exp(iω 2L/c)][1 + \exp (- i\omega ~ 2L/c)] [1+ \exp(i\omega ~ 2L/c)]

1+exp(iω 2L/c)+exp(iω 2L/c)+11 + \exp (i\omega ~ 2L/c) + \exp (-i\omega ~ 2L/c) + 1

2+2cos(ω 2L/c)2 + 2 \cos(\omega ~ 2L/c)

因为TDI引入的e指数因子,ORF最终为虚数,计算信噪比时取模。

相位读出(采样)与插值

以上所有分析都是基于连续时间函数,默认任意时刻的数据都可直接得到,但实际数据是离散采样的。进行时延操作时,若延迟时间与采样率不匹配,就无法直接得到所需数据,需进行插值操作。

cosϕ(t)\cos \phi(t)
利用光电感应器的强度变化可以记录瞬时相位,其中包含频率波动ω(t)\omega(t)以及引力波引入的相位差变动ϕ(t)\phi(t)。其中频率波动来源包括激光频率本身的不稳定性以及卫星运动引入的多普勒频移。

在5MHz到25MHz波动的拍频下准确读出相位
相位读出精度 10pm/Hz10 \rm pm/\sqrt{Hz}
相位噪声要求 1pm/Hz1 \rm pm/\sqrt{Hz}

自适应激光锁频、激光噪声抑制、星间激光通信

激光不稳定性除了涉及相位/频率还涉及功率,被称为相对强度噪声,会耦合进光电探测器的散粒噪声,影响观测。目前的激光SRIN104Hz\small \sqrt{S_{\rm RIN}} \sim 10^{-4} \sqrt{\rm Hz}

激光的相对强度噪声、电子噪声也会影响光学测量读出噪声,低频抬升?

位移噪声:

  • 光程噪声(optical pathlength noise):光学平台位移、望远镜位移,(TTL噪声)
  • 读出噪声(read out noise):载波读出噪声、边带读出噪声
    主要噪声来源为散粒噪声(shoot noise)、相对强度噪声(relative intensity noise, RIN)、电子噪声
  • 时钟噪声(clock noise):使用导频信号抑制时钟噪声,同时会引入新噪声
    • 读出:边带读出噪声
    • 传输:导频传输路径上存在电子元件、电缆、光电调制(相位)、光纤、光纤放大等一些列噪声源
  • 测量噪声(metrology system):激光频率残余噪声(TDI之后)、相位计噪声

位移噪声在六路激光测量间相对独立,加速度噪声对于共享测试质量的连接臂是相关的

时钟噪声抑制需要将时钟信号调制到载波信号,消耗10%功率
测距需要将伪随机码调制到载波信号,消耗1%功率

数据处理技术

除了TDI激光频率噪声外,空间数据分析还有一些的不同于地面观测的问题。
数据间断(gap)、非稳噪声(glitch)
地面的信号相对短,数据间断的影响更少,空间信号相对长,数据间断的影响更严重

白矮星前景扣除以及重叠信号的分离重构

数据异常处理

数据间断(gap)、非稳噪声(glitch)

混叠信号分离

白矮星前景扣除以及重叠信号的分离重构

Global Fit
LDASoft: LISA Data Analysis Software
https://github.com/tlittenberg/ldasoft

LISA data analysis using MCMC methods
LISA Gravitational Wave Sources in a Time-varying Galactic Stochastic Background
https://www.ipam.ucla.edu/abstract/?tid=17067

GPU加速计算

Katz et al. 2020