从0开始的EMRI

目前，EMRI 波形的理论建模主要有两种路径：

• 基于黑洞微扰理论，求解由小质量天体系统的近测地运动构成源项的 Teukolsky 方程。
• 基于小质量天体运动轨迹而直接利用引力波辐射近似公式得到 Kludge 数值波形。

前者基于微扰的理念正好适配极端质量比的情况，能够比较全面的对系统进行半解析建模；后者采用数值方法，能够快速的生成波形。而在未来的引力波探测中，往往需要快速的生成波形模板进行数据分析，因此如何快速的生成EMRI波形成为当下的主要问题，而 Kludge 波形正是应运而生。

接下来先展示如何通过黑洞微扰理论来得到EMRI的波形snapshot。

OverLook

黑洞微扰理论描述EMRI系统主要分为三大步：

• 构建小质量天体系统的运动轨迹，并给出对应的能动张量。
• 求解该源项对应的Teukolsky方程，得到静态的EMRI快照。
• 加入引力波辐射等因素造成的轨道衰减，重新迭代计算引力波形得到演化的EMRI波形。

本文内容描述了前两部分的流程，具体包含：

• 概述：
  • 从 Kerr 度规利用 NP 形式得到 Teukolsky 方程
  • 齐次径向方程的解法，得到了格林函数法的生成函数；齐次解的得到有如下几类方法：
    • MST法：利用超几何函数的级数展开构造（本文理论上概述了）
    • 基于 Heun 函数的数值方法（本文数值计算使用，因为速度比MST快近十倍）
    • SN法：考虑SN变换将径向方程变为具有短程势的形式
    • JH法：在无穷远和视界处级数展开，最后用常点处的级数展开进行拼接
• 详细推导：
  • 有源径向方程的格林函数解，最终应变波形的构造
  • 测地运动的 频率与轨迹解析解、数值计算理论
  • 利用BlackHoleToolKit 程序包数值计算实现 EMRI波形
• 下一步：
  • 学习方向：代码构建/齐次解推导/能流修正->绝热演化->EOB理论/MPD方程->自力理论
  • 研究想法：
    • 一般性的非真空轴对称Type-D度规（暗物质）/球对称薄盘度规 （Meng Kun 2026）
    • b-EMRI 有效单体描述 or 其他类似方法 （韩文标老师组大概在干这个）
    • 电磁对应：拉氏量加个电磁项/考虑 黑洞-中子双星的 EMRI/b-EMRI（报告听来的）
    • 记忆效应修正（和辐射反应的高阶修正一样弱）

其中对于Teukolsky方程的求解已经颇为成熟，主流方法是求解齐次方程后利用格林函数法求得非齐次方程的解。对于齐次Teukolsky方程，黑洞微扰理论中发展了多种方法，其中能够满足EMRI精度要求的方法有三种：MST法，SN法，Jiang-Han法。其中MST法为完全解析求解，精度最高但计算效率较慢；SN法在于将径向方程变换为类薛定谔方程形式来求解，轻巧便捷但精度较为一般；最新提出的Jiang-Han法则建立在连续分数法的基础上，在保持计算精度的前提下计算速度大幅提升。

对于小质量天体系统的动力学构建上，最早是考虑测试无自旋中性点粒子的类时测地运动，随后的MPD方程描述了带自旋的测试质量运动，以其为基础发展的Dixon形式描述了带自旋且具有内部多级矩结构的质量运动。需要注意的是，从点粒子到带自旋多极矩的测试质量，运动轨迹是从测地线到测地偏离运动发展，对应的运动方程也从具有解析解到目前只能数值求解，实际上对于EMRI系统，最难的部分就是运动轨迹的建模。

NOTE：本文表达式均遵循 BlackHoleToolKit 程序。

Teukolsky 方程与齐次情况求解

采取NP形式可以直接推导出Teukolsky方程，对NP方程进行形式理论推导后的主要计算量在于对各个非零旋量繁杂的显式计算，不过可以用已有的 Mathematica 程序完成，此处暂时掠过（至少对于Kerr黑洞以及非真空 TYPE-D 度规而言计算过程都没有很大的困难）：

$$\begin{array}{rl}& [\frac{(r^2+a^2{)}^2}{\mathrm{\Delta }}-a^2{\sin }^2\theta ]\frac{{\partial }^2{}_s\mathrm{\Psi }}{\partial t^2}+\frac{4Mar}{\mathrm{\Delta }}\frac{{\partial }^2{}_s\mathrm{\Psi }}{\partial t\partial \phi }+[\frac{a^2}{\mathrm{\Delta }}-\frac{1}{{\sin }^2\theta }]\frac{{\partial }^2{}_s\mathrm{\Psi }}{\partial {\phi }^2}\\ & -{\mathrm{\Delta }}^{-s}\frac{\partial }{\partial r}\left({\mathrm{\Delta }}^{s+1}\frac{\partial {}_s\mathrm{\Psi }}{\partial r}\right)-\frac{1}{\sin \theta }\frac{\partial }{\partial \theta }(\sin \theta \frac{\partial {}_s\mathrm{\Psi }}{\partial \theta })\\ & -2s[\frac{a(r-M)}{\mathrm{\Delta }}+\frac{i\cos \theta }{{\sin }^2\theta }]\frac{\partial {}_s\mathrm{\Psi }}{\partial \phi }-2s[\frac{M(r^2-a^2)}{\mathrm{\Delta }}-r-ia\cos \theta ]\frac{\partial {}_s\mathrm{\Psi }}{\partial t}\\ & +(s^2{\cot }^2\theta -s){}_s\mathrm{\Psi }=4\pi \mathrm{\Sigma }\hat{T}\end{array}$$ $$\begin{array}{rl}& \mathrm{\Psi }(t,r,\theta ,\phi )=\int \mathrm{d}\omega \sum _{l,m}{R}_{lm}(r){}_s{S}_{lm}(\theta )e^{-i\omega t+im\phi }\\ & 4\pi \mathrm{\Sigma }\hat{T}=\int \mathrm{d}\omega \sum _{l,m}{\mathcal{T}}_{lm}(r){}_s{S}_{lm}(\theta )e^{-i\omega t+im\phi }\end{array}$$

利用自旋加权球谐函数的正交归一性质，可以得到有源情况下的 Teukolsky 分离方程组：

$$\begin{array}{rl}& {\mathrm{\Delta }}^{-s}\frac{\mathrm{d}}{\mathrm{d}r}\left({\mathrm{\Delta }}^{s+1}\frac{\mathrm{d}{R}_{lm}}{\mathrm{d}r}\right)+(\frac{K^2-2is(r-M)K}{\mathrm{\Delta }}-\lambda )R_{lm}={\mathcal{T}}_{lm}\\ & \frac{1}{\sin \theta }\frac{\mathrm{d}}{\mathrm{d}\theta }[\sin \theta \frac{\mathrm{d}{}_s{S}_{lm}}{\mathrm{d}\theta }]+[-\frac{(m+s\cos \theta {)}^2}{{\sin }^2\theta }-a^2{\omega }^2{\sin }^2\theta -2a\omega s\cos \theta \\ & +s+2ma\omega +{\lambda }_{lm}]{}_s{S}_{lm}=0\\ & \mathrm{\Delta }=r^2-2Mr+a^2,K=(r^2+a^2)\omega -am,\mathrm{\Sigma }=r^2+a^2{\cos }^2\theta \end{array}$$

注意到角向方程就是自选加权椭谐函数，对于该函数的研究较为成熟，一般直接调用现有程序包。因此绝大部分注意力都在径向方程的求解。最主要的是齐次情况的径向方程求解，对于黑洞要求满足最基本的边界条件：

$$\begin{array}{rl}& R_{lm}\sim e^{-\mathrm{i}k{r}_{\ast }},r\to r_{+}\\ & R_{lm}\sim e^{\mathrm{i}\omega r_{\ast }},r\to \mathrm{\infty }\end{array}$$

在目前主流的方法是半解析的 MST，主要使用超几何函数展开得到解：

$$\begin{array}{rl}& R^{in}(x)=e^{\mathrm{i}ϵ\kappa x}(-x{)}^{-s-\mathrm{i}\frac{(ϵ+\tau )}{2}}(1-x{)}^{-\mathrm{i}\frac{(ϵ+\tau )}{2}}\sum _{n=-\mathrm{\infty }}^{\mathrm{\infty }}{f}_n{}_2{F}_1(n+\nu +1-\mathrm{i}\tau ,-n-\nu -\mathrm{i}\tau ;1-s-\mathrm{i}(ϵ+\tau );x)\\ & R^{up}(x)=2^{\nu }{e}^{-\pi ϵ}{e}^{-\mathrm{i}\pi (1+s+\nu )}{e}^{\mathrm{i}\hat{z}}{\hat{z}}^{\nu +\mathrm{i}{ϵ}_{+}}(\hat{z}-ϵ\kappa {)}^{-s-{ϵ}_{+}}\\ & \times \sum _{n=-\mathrm{\infty }}^{\mathrm{\infty }}{\mathrm{i}}^n\frac{\nu +1+s-\mathrm{i}ϵ}{\nu +1-s+\mathrm{i}ϵ}{f}_n(2\hat{z}{)}^{n}U(n+\nu +1+s-\mathrm{i}ϵ,2n+2\nu +2;-2\mathrm{i}z),\hat{z}=\omega (r-r_{-})\end{array}$$

也有更轻便的替代方法：SN形式，其将Teukolsky方程进行变换为具有短程势垒的形式而求解。
近期也有具有更快速度和较好精度的JH 方法，其分别在无穷远与视界处进行级数展开，最后在常点处进行级数展开来进行匹配拼接。而数值上一般利用 Heun 函数去数值求解。总之，不管是哪种方法，最后都是产生一对满足相同边界条件的线性无关的解函数 $R^{in}$ , $R^{up}$ ：

$$\begin{array}{rl}& R_{lm}^{in}⟶\{\begin{array}{c}{B}_{lm}^{trans}{\mathrm{\Delta }}^{-s}{e}^{-\mathrm{i}k{r}_{\ast }},r\to r_{+}\\ \\ r^{-1-2s}{B}_{lm}^{ref}{e}^{\mathrm{i}\omega r_{\ast }}+r^{-1}{B}_{lm}^{inc}{e}^{-\mathrm{i}\omega r_{\ast }},r\to \mathrm{\infty }\end{array}\\ \\ & R_{lm}^{up}⟶\{\begin{array}{c}{C}_{lm}^{up}{e}^{\mathrm{i}k{r}_{\ast }}+C_{lm}^{ref}{\mathrm{\Delta }}^{-s}{e}^{-\mathrm{i}k{r}_{\ast }},r\to r_{+}\\ \\ r^{-1-2s}{C}_{lm}^{trans}{e}^{\mathrm{i}\omega r_{\ast }},r\to \mathrm{\infty }\end{array}\end{array}$$

其中 "in" 解满足在视界处纯入射、"up" 解满足在无穷远纯出射 的边界条件，并且能够高效精准的计算渐近振幅： $B^{trans},B^{ref},B^{inc};C^{trans},C^{ref},C^{up}$ 。
note：一般会令 $C^{trans}=1$ 来归一化。

对于需要考虑的EMRI系统，只要不对度规进行修改或者考虑高阶扰动，那么所需要的齐次解是不变的，因此以上仅先做简单介绍，具体推导细节暂时搁置。

有源Teukolsky 方程的求解

实际上是对非齐次径向方程求解，采用Green 函数法。因此首先概述所用到的Green函数法：

考虑格林函数满足的二阶常微分方程：

$$\begin{array}{rl}& \frac{\mathrm{d}}{\mathrm{d}r}(p(r;\xi )\frac{\mathrm{d}G}{\mathrm{d}r})+q(r;\xi )=\delta (r-\xi )\end{array}$$

由黑洞视界纯入射、无穷远纯出射波的边界条件，可以假设方程具有如下形式的解：

$$\begin{array}{rl}& G(r;\xi )=\{\begin{array}{c}{c}_1(\xi )y_1(r),r_{+}<r<\xi \\ c_2(\xi )y_2(r),r>\xi \end{array}\end{array}$$

其中 $y_1,y_2$ 是满足齐次情况的解，利用 $\delta$ 函数的连续性方程容易得到：

$$\begin{array}{rl}& c_1=\frac{y_2(\xi )}{p(\xi )w(\xi )}=\frac{y_2(\xi )}{W(\xi )},c_1=\frac{y_1(\xi )}{p(\xi )w(\xi )}=\frac{y_1(\xi )}{W(\xi )}\end{array}$$

其中 $W(\xi )$ 为朗斯基行列式，由边界条件可以给出：

$$\begin{array}{rl}& G(r;\xi )=\{\begin{array}{c}\frac{1}{W}{R}_{lm}^{up}(\xi )R_{lm}^{in}(r),r_{+}<r<\xi \\ \frac{1}{W}{R}_{lm}^{in}(\xi )R_{lm}^{up}(r),r>\xi \end{array}\end{array}$$

可以证明，朗斯基行列式是变量无关的常数，即可以证明等式 ${[p(r)(y_1{y}_2^{\prime }-y_1^{\prime }{y}_2)]}_a^b=0$ 成立，方法是考虑 $y_1,y_2$ 均是齐次情况方程解，联立整理后求分步积分即可。因为对于引力波辐射最终要的是无穷远处的波形，因此通常在无穷远处计算该行列式，稍作计算并整理后结果为：

$$\begin{array}{rl}W(\mathrm{\infty })& ={[p(r)(y_1\frac{\mathrm{d}{y}_2}{\mathrm{d}r}+y_2\frac{\mathrm{d}{y}_1}{\mathrm{d}r})]}_{r\to \mathrm{\infty }}\\ & =2\mathrm{i}\omega B_{lm}^{inc}{C}_{lm}^{trans}\end{array}$$

最后，可以证明对于一般性的边界条件，格林函数具有对称性 $(r;\xi )\leftrightarrow (\xi ,r)$ ，由此可以得到非齐次径向方程的格林函数解：

$$\begin{array}{rl}{R}_{lm}(r)& ={\int }_{r_{+}}^{\mathrm{\infty }}G(r;\xi ){\mathrm{\Delta }}^s{\mathcal{T}}_{lm}\mathrm{d}\xi ={\int }_{r_{+}}^{\mathrm{\infty }}G(\xi ;r){\mathrm{\Delta }}^s{\mathcal{T}}_{lm}\mathrm{d}\xi \\ & =\frac{1}{W_{lm}}\{R_{lm}^{up}(r){\int }_{{}_{+}}^r{R}_{lm}^{in}{\mathrm{\Delta }}^s(\xi ){\mathcal{T}}_{lm}(\xi )\mathrm{d}\xi +R_{lm}^{in}(r){\int }_{{}_r}^{\mathrm{\infty }}{R}_{lm}^{up}{\mathrm{\Delta }}^s(\xi ){\mathcal{T}}_{lm}(\xi )\mathrm{d}\xi \}\end{array}$$

注意到由此便有：

$$\begin{array}{rl}& R_{lm}(r\to \mathrm{\infty })=\frac{R_{lm}^{up}(r\to \mathrm{\infty })}{W_{lm}}{\int }_{{}_{+}}^r{R}_{lm}^{in}{\mathrm{\Delta }}^s(\xi ){\mathcal{T}}_{lm}(\xi )\mathrm{d}\xi ={\tilde{Z}}_{lm}^{in}{R}_{lm}^{up}(r\to \mathrm{\infty })\\ & R_{lm}(r\to r_{+})=\frac{R_{lm}^{in}(r\to r_{+})}{W_{lm}}{\int }_{{}_{+}}^r{R}_{lm}^{up}{\mathrm{\Delta }}^s(\xi ){\mathcal{T}}_{lm}(\xi )\mathrm{d}\xi ={\tilde{Z}}_{lm}^{up}{R}_{lm}^{in}(r\to \mathrm{\infty })\end{array}$$

接下来计算 源项 $T_{lm}$ ，由傅里叶逆变换与自旋椭谐函数的正交归一性有：

$$\begin{array}{rl}& {\mathcal{T}}_{lm}=\frac{1}{2\pi }{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}\mathrm{d}t\int \mathrm{d}\mathrm{\Omega }4\pi \mathrm{\Sigma }\hat{T}{e}^{\mathrm{i}\omega t-\mathrm{i}m\phi }{}_s{S}_{lm}(\theta )\end{array}$$

此处选取 $\int {|{}_s{S}_{lm}(\theta )|}^2\mathrm{d}\theta =\frac{1}{2\pi }$ (与 BlackHoleToolkit 保持一致)，而源项 $4\pi \mathrm{\Sigma }\hat{T}$ 为：

$$\begin{array}{rl}4\pi \mathrm{\Sigma }\hat{T}& =4\pi \mathrm{\Sigma }(2{\rho }^{-4}{T}_4)\\ T_4& =B_2^{\prime }+B_2^{{\ast }^{\prime }}\\ & =-\frac{{\rho }^8\overline{\rho }}{2}{L}_{-2}[{\rho }^{-4}{L}_0({\rho }^{-2}{\overline{\rho }}^{-1}{T}_{nn})]+\frac{{\mathrm{\Delta }}^2{\rho }^8\overline{\rho }}{2\sqrt{2}}{L}_{-1}[{\rho }^{-4}{\rho }^{-2}{J}_{+}({\rho }^{-2}{\overline{\rho }}^{-2}{\mathrm{\Delta }}^{-1}{T}_{n\overline{m}})]\\ & +\frac{{\mathrm{\Delta }}^2{\rho }^8\overline{\rho }}{2\sqrt{2}}{J}_{+}[{\rho }^{-4}{\rho }^{-2}{\mathrm{\Delta }}^{-1}{L}_{-1}({\rho }^{-2}{\overline{\rho }}^{-2}{T}_{n\overline{m}})]-\frac{{\mathrm{\Delta }}^2{\rho }^8\overline{\rho }}{4}{J}_{+}[{\rho }^{-4}{J}_{+}({\rho }^{-2}\overline{\rho }{T}_{\overline{m}\overline{m}})]\\ \rho & =-\frac{1}{r-\mathrm{i}a\cos \theta },\overline{\rho }=-\frac{1}{r+\mathrm{i}a\cos \theta },\mathrm{\Delta }=r^2-2rM+a^2,\mathrm{\Sigma }=r^2+a^2\cos \theta \\ J_{+}& =\frac{\partial }{\partial r}+\frac{\mathrm{i}K(r)}{\mathrm{\Delta }},L_s=\frac{\partial }{\partial \theta }+m\csc \theta -a\omega \sin \theta +s\cot \theta \end{array}$$

其中 $T_{ab}=T^{\alpha \beta }{n}_{\alpha }{m}_{\beta }$ 便是测试粒子的能动张量在NP标架两个方向的投影。对于点粒子，其能动张量可以构造为：

$$\begin{array}{rl}{T}^{\alpha \beta }& =\mu \int u^{\alpha }{u}^{\beta }{\delta }^{(4)}[x-z(\tau )]\mathrm{d}\tau \\ & =\mu \int u^{\alpha }{u}^{\beta }(-g{)}^{-\frac{1}{2}}\delta [t-\tilde{t}(\tau )]\delta [\theta -\tilde{\theta }(\tau )]\delta [r-\tilde{r}(\tau )]\delta [\phi -\tilde{\phi }(\tau )]\mathrm{d}\tau \end{array}$$

其中 $g$ 为度规对应的行列式，考虑到粒子的坐标时间与固有时间是对应的（局域上实际上就是一回事），因此能动张量化为：

$$\begin{array}{rl}& T^{\alpha \beta }=\mu \frac{u^{\alpha }{u}^{\beta }}{\dot{t}\mathrm{\Sigma }\sin \theta }\delta [\theta -\tilde{\theta }(t)]\delta [r-\tilde{r}(t)]\delta [\phi -\tilde{\phi }(t)]\end{array}$$

而 NP 标架为：

$$\begin{array}{rl}& n_{\alpha }=\frac{1}{2}[\frac{\mathrm{\Delta }}{\mathrm{\Sigma }},1,0,-\frac{a\mathrm{\Delta }{\sin }^2\theta }{\mathrm{\Sigma }}],l_{\alpha }=\frac{1}{\mathrm{\Delta }}[\mathrm{\Delta },-\mathrm{\Sigma },0,a\mathrm{\Delta }{\sin }^2\theta ]\\ & m_{\alpha }=-\frac{\overline{\rho }}{\sqrt{2}}[,-\mathrm{\Sigma },0,a\mathrm{\Delta }{\sin }^2\theta ],\overline{m}\end{array}$$

因此：

$$\begin{array}{rl}& T_{ab}=\mu \frac{C_{ab}}{\sin \theta }\delta [\theta -\tilde{\theta }(\tau )]\delta [r-\tilde{r}(\tau )]\delta [\phi -\tilde{\phi }(\tau )]\\ & C_{nn}=\frac{\mathrm{d}\lambda }{\mathrm{d}t}\frac{1}{4\mathrm{\Sigma }}{[\mathcal{E}(r^2+a^2)-a{L}_z+\frac{\mathrm{d}r}{\mathrm{d}\lambda }]}^2\\ & C_{n\overline{m}}=\frac{\mathrm{d}\lambda }{\mathrm{d}t}\frac{\rho }{2\sqrt{2}\mathrm{\Sigma }}[\mathcal{E}(r^2+a^2)-a{L}_z+\frac{\mathrm{d}r}{\mathrm{d}\lambda }][\frac{\mathrm{d}\theta }{\mathrm{d}\lambda }+\mathrm{i}\sin \theta (a\mathcal{E}-\frac{L_z}{{\sin }^2\theta })]\\ & C_{\overline{m}\overline{m}}=\frac{\mathrm{d}\lambda }{\mathrm{d}t}\frac{{\rho }^2}{2\mathrm{\Sigma }}{[\frac{\mathrm{d}\theta }{\mathrm{d}\lambda }+\mathrm{i}\sin \theta (a\mathcal{E}-\frac{L_z}{{\sin }^2\theta })]}^2\end{array}$$

又由于：

$$\begin{array}{rl}\int f(r,\theta )\delta [\phi -\tilde{\phi }(t)]e^{-\mathrm{i}m\phi }\mathrm{d}\mathrm{\Omega }& ={\int }_0^{\pi }f(r,\theta )\mathrm{d}\theta {\int }_0^{2\pi }\delta [\phi -\tilde{\phi }(t)]e^{-\mathrm{i}m\phi }\mathrm{d}\phi \\ & ={\int }_0^{\pi }f(r,\theta )e^{-\mathrm{i}m\tilde{\phi }(t)}\mathrm{d}\theta \end{array}$$

接下来就是处理 $\delta [r-\tilde{r}],\delta [\theta -\tilde{\theta }]$ ，因为 $\delta$ 函数在能动张量中，而源项表达式中能动张量被算符 $J_{+},L_s$ 包裹，因此需要将其运用分部积分提取出来；因为径向部分要在非齐次解中进行积分，因此首先处理角向部分，步骤为：首先对于角向部分只需要考虑将其从 $L_s$ 算符中提取出来，其次将被积函数全部展开化简，最终所有的角向算符 $L_s,L_s^{+}$ 均只作用在自选加权椭谐函数 ${}_{-2}{S}_{lm}(\theta )$ （后续均是 $s=-2$ 的结果）；主要运用到了三个分部积分等式：

$$\begin{array}{rl}& \int g(r,\theta )L_s[h(r,\theta )]\sin \theta \mathrm{d}\theta =-\int h(r,\theta )L_{1-s}^{+}[g(r,\theta )]\sin \theta \mathrm{d}\theta \\ & J_{+}(gh)=g{J}_{+}(h)+h\frac{\partial }{\partial r}g\\ & L_s(gh)=g{L}_s(h)+h\frac{\partial }{\partial \theta }g\end{array}$$

以及一个对易关系：

$$\begin{array}{rl}& [J_{+},L_s]=0\end{array}$$

以上四个关系中有：

$$\begin{array}{r}{L}_s^{+}=\frac{\partial }{\partial \theta }-m\csc \theta +a\omega \sin \theta +s\cot \theta \end{array}$$

其中 $L_s^{+},L_s$ 本质是一致的，满足一样的等式关系。最后经过近十页的整理并且根据径向导数的阶数整理，可以得到结果：

$$\begin{array}{rl}{\mathcal{T}}_{lm}& =\mu \int e^{i[\omega t-m\tilde{\phi }(t)]}\mathrm{d}t{\mathrm{\Delta }}^2\{\\ & [-\frac{2{\rho }^{-3}{\overline{\rho }}^{-1}}{{\mathrm{\Delta }}^2}{C}_{nn}(L_1^{+}{L}_2^{+}({}_{-2}{S}_{lm})+2ia\sin \theta L_2^{+}({}_{-2}{S}_{lm}))\\ & -\frac{2\sqrt{2}{\rho }^{-3}}{\mathrm{\Delta }}{C}_{n\overline{m}}[(\frac{iK}{\mathrm{\Delta }}-\rho -\overline{\rho })L_2^{+}({}_{-2}{S}_{lm})+\frac{iK}{\mathrm{\Delta }}(\rho -\overline{\rho })ia\sin \theta {}_{-2}{S}_{lm}]\\ & -{\rho }^{-3}\overline{\rho }{C}_{\overline{m}\overline{m}}{}_{-2}{S}_{lm}[{\left(\frac{iK}{\mathrm{\Delta }}\right)}^2+2\rho \frac{iK}{\mathrm{\Delta }}+\frac{\partial }{\partial r}\left(\frac{iK}{\mathrm{\Delta }}\right)]]\delta [r-\tilde{r}(t)]\\ & +\frac{\mathrm{d}}{\mathrm{d}r}[-\frac{2\sqrt{2}{\rho }^{-3}}{\mathrm{\Delta }}{C}_{n\overline{m}}(L_2^{+}({}_{-2}{S}_{lm})+ia\sin \theta (\rho -\overline{\rho }){}_{-2}{S}_{lm})\\ & -2{\rho }^{-3}{\overline{\rho }}^{-1}(\frac{iK}{\mathrm{\Delta }}-\rho ){}_{-2}{S}_{lm}{C}_{\overline{m}\overline{m}}]\delta [r-\tilde{r}(t)]\\ & +\frac{{\mathrm{d}}^2}{\mathrm{d}{r}^2}[-{\rho }^{-3}{\overline{\rho }}^{-1}{}_{-2}{S}_{lm}{C}_{\overline{m}\overline{m}}\delta [r-\tilde{r}(t)]]\}\\ & =\mu \int e^{i[\omega t-m\tilde{\phi }(t)]}\mathrm{d}t{\mathrm{\Delta }}^2\{A_0\delta [r-\tilde{r}(t)]+\frac{\mathrm{d}}{\mathrm{d}r}[A_1\delta [r-\tilde{r}(t)]]+\frac{{\mathrm{d}}^2}{\mathrm{d}{r}^2}[A_2\delta [r-\tilde{r}(t)]]\}\end{array}$$

最后便可以计算非齐次径向方程的格林函数解：

$$\begin{array}{rl}{R}_{lm}(r)& ={\int }_{r_{+}}^{\mathrm{\infty }}G(r;\xi ){\mathrm{\Delta }}^s{\mathcal{T}}_{lm}\mathrm{d}\xi ={\int }_{r_{+}}^{\mathrm{\infty }}G(\xi ;r){\mathrm{\Delta }}^s{\mathcal{T}}_{lm}\mathrm{d}\xi \\ & =\frac{1}{W_{lm}}\{R_{lm}^{up}(r){\int }_{{}_{+}}^r{R}_{lm}^{in}{\mathrm{\Delta }}^s(\xi ){\mathcal{T}}_{lm}(\xi )\mathrm{d}\xi +R_{lm}^{in}(r){\int }_{{}_r}^{\mathrm{\infty }}{R}_{lm}^{up}{\mathrm{\Delta }}^s(\xi ){\mathcal{T}}_{lm}(\xi )\mathrm{d}\xi \}\end{array}$$

最关心的是无穷远处的出射波，因此仅计算：

$$\begin{array}{rl}& R_{lm}(r\to \mathrm{\infty })=\frac{R_{lm}^{up}(r\to \mathrm{\infty })}{W_{lm}}{\int }_{{}_{+}}^r{R}_{lm}^{in}{\mathrm{\Delta }}^s(\xi ){\mathcal{T}}_{lm}(\xi )\mathrm{d}\xi ={\tilde{Z}}_{lm}^{in}{R}_{lm}^{up}(r\to \mathrm{\infty })=r^3{e}^{\mathrm{i}\omega r_{\ast }}{\tilde{Z}}_{lm}^{in}\end{array}$$ $$\begin{array}{rl}{\tilde{Z}}_{lm}^{in}& =\frac{\mu }{2\mathrm{i}\omega B_{lm}^{inc}}{\int }_{r_{+}}^{\mathrm{\infty }}{R}_{lm}^{in}(\xi ){\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{e}^{\mathrm{i}[\omega t-m\tilde{\phi }(t)]}\mathrm{d}t\{[A_0(\xi ,\tilde{\theta }(t))\delta [r-\tilde{r}(t)]]\\ & +\frac{\mathrm{d}}{\mathrm{d}r}[A_1(\xi ,\tilde{\theta }(t))\delta [r-\tilde{r}(t)]]+\frac{{\mathrm{d}}^2}{\mathrm{d}{r}^2}[A_2(\xi ,\tilde{\theta }(t))\delta [r-\tilde{r}(t)]]\}\mathrm{d}{r}^{\prime }\end{array}$$

利用 $\delta$ 函数及其导数的性质

$$\begin{array}{r}\frac{\mathrm{d}\delta (x)}{\mathrm{d}x}=0,x\ne 0\end{array}$$

最终得到振幅表达式：

$$\begin{array}{rl}{\tilde{Z}}_{lm}^{in}& =\frac{\mu }{2\mathrm{i}\omega B_{lm}^{inc}}{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{e}^{\mathrm{i}[\omega t-m\tilde{\phi }(t)]}\mathrm{d}t\{A_0{R}_{lm}^{in}-A_1{\frac{\mathrm{d}{R}_{lm}^{in}}{\mathrm{d}\xi }|}_{\xi =\tilde{r}(t)}+A_0{\frac{{\mathrm{d}}^2{R}_{lm}^{in}}{\mathrm{d}{\xi }^2}|}_{\xi =\tilde{r}(t)}\}\\ & =\frac{\mu }{2\mathrm{i}\omega B_{lm}^{inc}}{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{e}^{\mathrm{i}[\omega t-m\tilde{\phi }(t)]}\mathrm{d}t{I}_{lm}\end{array}$$

由测地线运动可知： $t=\mathrm{\Gamma }\lambda +\mathrm{\Delta }t[\tilde{r}(\lambda ),\tilde{\theta }(\lambda )],\tilde{\phi }={\mathrm{{\rm Y}}}_{\phi }\lambda +\mathrm{\Delta }\phi [\tilde{r}(\lambda ),\tilde{\theta }(\lambda )]$ ，由此振幅化为：

$$\begin{array}{rl}{\tilde{Z}}_{lm}^{in}& =\frac{\mu }{2\mathrm{i}\omega B_{lm}^{inc}}{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{e}^{\mathrm{i}[\omega t-m\tilde{\phi }(t)]}\mathrm{d}t{I}_{lm}\\ & =\frac{\mu }{2\mathrm{i}\omega B_{lm}^{inc}}{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{e}^{\mathrm{i}\lambda [\omega \mathrm{\Gamma }-m\mathrm{{\rm Y}}\phi ]}\mathrm{d}\lambda e^{\mathrm{i}(\omega \mathrm{\Delta }t-m\mathrm{\Delta }\phi )}{I}_{lm}\frac{\mathrm{d}t}{\lambda }\\ & =\frac{\mu }{2\mathrm{i}\omega B_{lm}^{inc}}{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{e}^{\mathrm{i}\lambda [\omega \mathrm{\Gamma }-m\mathrm{{\rm Y}}\phi ]}\mathrm{d}\lambda J_{lm}\end{array}$$

对于被积分函数 $J_{lm}$ 同样由于角向与径向的可分离性，能够进行傅里叶展开：

$$\begin{array}{rl}& J_{lm}=\sum _{k,n}{J}_{lmnk}{e}^{\mathrm{i}(n{\mathrm{{\rm Y}}}_r+k{\mathrm{{\rm Y}}}_{\theta })\lambda },{\omega }_r=n{\mathrm{{\rm Y}}}_r,{\omega }_{\theta }=n{\mathrm{{\rm Y}}}_{\theta }\\ & J_{lmnk}=\frac{1}{4{\pi }^2}{\int }_0^{2\pi }\mathrm{d}{\omega }_r{\int }_0^{2\pi }\mathrm{d}{\omega }_{\theta }{e}^{\mathrm{i}(n{\omega }_r+k{\omega }_{\theta })}{J}_{lm}[r({\omega }_r),\theta ({\omega }_{\theta }),r,\lambda ]\end{array}$$

最后借助 $\delta$ 函数的傅里叶变换方程得到振幅最终表达式：

$$\begin{array}{rl}{\tilde{Z}}_{lm}^{in}& =\frac{\mu }{2\mathrm{i}\omega B_{lm}^{inc}}{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{e}^{\mathrm{i}\lambda [\omega \mathrm{\Gamma }-m\mathrm{{\rm Y}}\phi ]}\mathrm{d}\lambda J_{lm}\\ & =\frac{\mu }{2\mathrm{i}\omega B_{lm}^{inc}}{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}\sum _{n,k}{e}^{\mathrm{i}\lambda [(\omega \mathrm{\Gamma }-m\mathrm{{\rm Y}}\phi )-(n{\mathrm{{\rm Y}}}_r+k{\mathrm{{\rm Y}}}_{\theta })]}{J}_{lmnk}\mathrm{d}\lambda \\ & =\frac{\mu }{2\mathrm{i}\omega B_{lm}^{inc}}{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}\sum _{n,k}{e}^{\mathrm{i}\lambda \mathrm{\Gamma }[\omega -(m\frac{\mathrm{{\rm Y}}\phi }{\mathrm{\Gamma }}+n\frac{{\mathrm{{\rm Y}}}_r}{\mathrm{\Gamma }}+k\frac{{\mathrm{{\rm Y}}}_{\theta }}{\mathrm{\Gamma }})]}{J}_{lmnk}\mathrm{d}\lambda \\ & =\frac{\mu }{2\mathrm{i}\omega B_{lm}^{inc}}\sum _{n,k}\frac{2\pi J_{lmkn}}{\mathrm{\Gamma }}\delta (\omega -{\omega }_{mnk}),{\omega }_{mnk}=m\frac{\mathrm{{\rm Y}}\phi }{\mathrm{\Gamma }}+n\frac{{\mathrm{{\rm Y}}}_r}{\mathrm{\Gamma }}+k\frac{{\mathrm{{\rm Y}}}_{\theta }}{\mathrm{\Gamma }}\end{array}$$

最终无穷远处的波动函数写为：

$$\begin{array}{rl}{\psi }_4^{r\to \mathrm{\infty }}& ={\rho }^4\sum _{l,m,n,k}\int \mathrm{d}\omega R_{lm}(r){}_s{S}_{lm}(\theta )e^{-\mathrm{i}\omega t+\mathrm{i}m\phi }\\ & =\frac{\mu }{r}\sum _{l,m,n,k}\frac{e^{-\mathrm{i}{\omega }_{mnk}(t-r_{\ast })}}{2\mathrm{i}{\omega }_{mnk}{B}_{lm}^{inc}}\frac{2\pi J_{lmkn}}{\mathrm{\Gamma }}{}_{-2}{S}_{lm}^{a\omega }(\theta ,\phi )\end{array}$$

对应的引力波应变为：

$$\begin{array}{rl}{h}_{+}-\mathrm{i}{h}_{\times }& =-\frac{2{\psi }_4^{r\to \mathrm{\infty }}}{{\omega }_{mnk}^2}\\ & =-\frac{2\mu }{r}\sum _{l,m,n,k}\frac{e^{-\mathrm{i}{\omega }_{mnk}(t-r_{\ast })}}{{\omega }_{mnk}^2}\frac{2\pi J_{lmkn}}{2\mathrm{i}{\omega }_{mnk}{B}_{lm}^{inc}\mathrm{\Gamma }}{}_{-2}{S}_{lm}^{a\omega }(\theta ,\phi )\\ & =\frac{2\mu }{r}\sum _{l,m,n,k}\frac{e^{-\mathrm{i}{\omega }_{mnk}(t-r_{\ast })}}{{\omega }_{mnk}^2}{\tilde{Z}}_{lmnk}^{in}{}_{-2}{S}_{lm}^{a\omega }(\theta ,\phi )\end{array}$$

至此便完成了对于无演化的波形快照建模，可以发现在数值计算上需要做的就是计算那个产生于傅里叶展开的二重积分，也是这个积分的计算需要调用径向方程的齐次解，因此该积分的计算速度几乎是等于波形的计算速度，该积分的收敛性基本决定了波形的收敛性，对于径向上当高偏心可能出现的积分发散也是体现在这个二重积分上。

最后关于初值的问题，当不在转折点开始时，需要加上一个初始相位，体现在时间上便是，如果从转折点起步，则是 $\lambda =0$ 开始运动，而如果具有初始相位，则不再从转折点起步，此时计时应该是从 $\lambda ={\lambda }_0$ 开始，令 $\lambda =\tilde{\lambda }+{\lambda }_0$ ，其中起始有 $\tilde{\lambda }=0$ ，对应到波形计算过程便是存在一个相位修正：

$$\begin{array}{r}{\tilde{Z}}_{lm}^{in}⟹{\tilde{Z}}_{lm}^{in}{e}^{\mathrm{i}[\omega (q_{t_0}-\mathrm{\Delta }{t}_0-n{\omega }_{r_0})-m(q_{{\phi }_0}-\mathrm{\Delta }{\phi }_0+k{\omega }_{r_0})]}\end{array}$$

不过一般考虑初始相位等于 0 。

测试粒子的轨迹 —— 能动张量的构造

Kerr 度规

在BL系下Kerr度规写为：

$$\begin{array}{rl}\mathrm{d}{s}^2& =-(1-\frac{2Mr}{{\rho }^2})\mathrm{d}{t}^2+\frac{{\rho }^2}{\mathrm{\Delta }}\mathrm{d}{r}^2+{\rho }^2\mathrm{d}{\theta }^2\\ & +[(r^2+a^2){\sin }^2\theta +\frac{2Mr{a}^2{\sin }^4\theta }{{\rho }^2}]\mathrm{d}{\phi }^2-\frac{4Mra{\sin }^2\theta }{{\rho }^2}\mathrm{d}t\mathrm{d}\phi \\ & {\rho }^2\equiv r^2+a^2{\cos }^2\theta ,\mathrm{\Delta }\equiv r^2-2Mr+a^2\end{array}$$

其具有两个 Killing 矢量： $k^{\mu },m^{\mu }$ ，对应 $t,\phi$ 两个时空对称性，即能量 $E$ 与角动量 $L_z$ 守恒。

此外还存在一个隐藏的 Killing 矢量，对应 Carter 常数，后续可以在哈密顿-雅可比理论框架下利用分离变量将其给出。

运动方程

考虑一个测试点粒子，我们的目的便是得到其类时测地线方程。哈密顿-雅可比框架是一个简单快速的方法，并且可以直接得到 Carter 常数；接下来我们便采用这一方法导出测地线。

点源的拉格朗日量为（对应为测试粒子的拉氏密度，对应的守恒量也应为密度）：

$$\begin{array}{r}\mathcal{L}=\frac{1}{2}{g}_{\mu \nu }{\dot{x}}^{\mu }{\dot{x}}^{\nu }\end{array}$$

容易发现由于两个 Killing 矢量的存在，对应的 $t,\phi$ 是循环坐标，正是对应了两个守恒量：能量密度 $\mathcal{E}$ 与角动量密度： $L_z$ 。利用 Killing 矢量的性质可以直接给出,其中 $\mu =0,1,2,3$ ：

$$\begin{array}{r}\mathcal{E}=-u^{\mu }{k}_{\mu }=-u_t=-g_{t\mu }{u}^{\mu }=A\dot{t}+B\dot{\phi }\\ L_z=u^{\mu }{m}_{\mu }=u_{\phi }=g_{\phi \mu }{u}^{\mu }=-B\dot{t}+C\dot{\phi }\end{array}$$

自然可以反解出 $\dot{t},\dot{\phi }$ :

$$\begin{array}{rl}\frac{\mathrm{d}t}{\mathrm{d}\tau }& =\frac{C\mathcal{E}-B{L}_z}{\mathrm{\Delta }{\sin }^2\theta }=\frac{1}{\mathrm{\Delta }}[(r^2+a^2+\frac{2Mr{a}^2}{\mathrm{\Sigma }}{\sin }^2\theta )\mathcal{E}-\frac{2Mra}{\mathrm{\Sigma }}{L}_z]\\ \frac{\mathrm{d}\phi }{\mathrm{d}\tau }& =\frac{B\mathcal{E}+A{L}_z}{\mathrm{\Delta }{\sin }^2\theta }=\frac{1}{\mathrm{\Delta }}[(1-\frac{2Mr}{\mathrm{\Sigma }})\frac{L_z}{{\sin }^2\theta }+\frac{2Mra}{\mathrm{\Sigma }}\mathcal{E}]\end{array}$$

接着考虑 哈密顿-雅可比方程：

$$\begin{array}{r}\frac{\partial S}{\partial \lambda }+H(x^{\mu },\frac{\partial S}{\partial x^{\mu }})=0\end{array}$$

这是个四维方程，具有四个待求变量，因此如果解 $S$ 依赖于四个独立的常数，那么便构成完全积分，由此便有：

$$\begin{array}{r}\frac{\partial S}{\partial x^{\mu }}=p_{\mu }\end{array}$$

由欧拉-拉格朗日方程，可以得到粒子的正则动量：

$$p_{\mu }=\frac{\partial \mathcal{L}}{\partial {\dot{x}}^{\mu }}=g_{\mu \nu }{\dot{x}}^{\mu }$$

因此对应的哈密顿量（密度）可以通过勒让德变换得到：

$$H(x^{\mu },p_{\nu })=p_{\mu }{\dot{x}}^{\mu }(p_{\nu })-\mathcal{L}(x^{\mu },{\dot{x}}^{\mu }(p_{\nu }))=\frac{1}{2}{g}^{\mu \nu }{p}_{\mu }{p}_{\nu }$$

则容易由四速度的归一性： $g_{mu\nu }{\dot{x}}^{\mu }{\dot{x}}^{\nu }=-1$ 得到三个常数：

$$\begin{array}{rl}H(x^{\mu },p_{\nu })& =\frac{1}{2}\kappa =const\\ p_t& =-\mathcal{E}\\ p_{\phi }& =L_z\end{array}$$

为使得系统可积，接下来寻找第四个常数，将作用量泛函设为：

$$S=-\frac{1}{2}\kappa -\mathcal{E}t+L_z\phi +S^{r\theta }(r,\theta )$$

容易发现是符合已有守恒量要求的，接着假设可以进行分离变量（以此分离常数便是需要寻找的第四个常数）：

$$\begin{array}{r}{S}^{r\theta }(r,\theta )=S^r(r)+S^{\theta }(\theta )\end{array}$$

因此可以由：

$$\begin{array}{r}\frac{\partial S}{\partial \lambda }+H(x^{\mu },\frac{\partial S}{\partial x^{\mu }})=-\frac{1}{2}\kappa +\frac{1}{2}{g}^{\mu \nu }\frac{\partial S}{\partial x^{\mu }}\frac{\partial S}{\partial x^{\nu }}=0\end{array}$$

得到方程：

$$\begin{array}{r}-\kappa +\frac{\mathrm{\Delta }}{\mathrm{\Sigma }}(\frac{\mathrm{d}{S}^r}{\mathrm{d}r}{)}^2+\frac{1}{\mathrm{\Sigma }}(\frac{\mathrm{d}{S}^{\theta }}{\mathrm{d}\theta }{)}^2-\frac{1}{\mathrm{\Delta }}[r^2+a^2+\frac{2Mr{a}^2}{\mathrm{\Sigma }}{\sin }^2\theta ]{\mathcal{E}}^2+\frac{4Mra}{\mathrm{\Sigma }}\mathcal{E}{L}_z+\frac{\mathrm{\Delta }-a^2{\sin }^2\theta }{\mathrm{\Sigma }\mathrm{\Delta }{\sin }^2\theta }{L}_z=0\end{array}$$

这是可分离的，化简得到：

$$\begin{array}{rl}& \mathrm{\Delta }(\frac{\mathrm{d}{S}^r}{\mathrm{d}r}{)}^2-\kappa r^2-\frac{1}{\mathrm{\Delta }}{[\mathcal{E}(r^2+a^2)-a{L}_z]}^2+(L_z-a\mathcal{E}{)}^2=-Q\\ & (\frac{\mathrm{d}{S}^{\theta }}{\mathrm{d}\theta }{)}^2-{\cos }^2\theta [{\mathcal{E}}^2{a}^2+\kappa a^2-\frac{L_z^2}{1-{\cos }^2\theta }]=Q\end{array}$$

其中 $Q$ 称之为 Carter 常数（也可以定义： $q=Q+(L_z-a\mathcal{E}{)}^2$ 为 Carter 常数）。由此可以进一步利用： $\frac{\partial S}{\partial x^{\mu }}=p_{\mu }=g_{\mu \nu }{\dot{x}}^{\nu }$ 得到径向与角向的运动方程：

$$\begin{array}{rl}(g_{rr}\frac{\mathrm{d}r}{\mathrm{d}\tau }{)}^2& ={\left(\frac{\mathrm{\Sigma }}{\mathrm{\Delta }}\frac{\mathrm{d}r}{\mathrm{d}\tau }\right)}^2=\frac{1}{{\mathrm{\Delta }}^2}{[\mathcal{E}(r^2+a^2)-a{L}_z]}^2+\frac{\kappa r^2+-(L_z-a\mathcal{E}{)}^2-Q}{\mathrm{\Delta }}\\ (g_{\theta \theta }\frac{\mathrm{d}\theta }{\mathrm{d}\tau }{)}^2& ={\left(\mathrm{\Sigma }\frac{\mathrm{d}\theta }{\mathrm{d}\tau }\right)}^2={\cos }^2\theta [{\mathcal{E}}^2{a}^2+\kappa a^2-\frac{L_z^2}{1-{\cos }^2\theta }]+Q\end{array}$$

我们还希望把运动方程中对角向与径向的依赖函数分离开来，而引入 Mino time: $\lambda$ 便可以做到：

$$\begin{array}{r}\mathrm{d}\tau =\mathrm{\Sigma }\mathrm{d}\lambda \end{array}$$

由此运动方程为：

$$\begin{array}{rl}& {\left(\frac{\mathrm{d}r}{\mathrm{d}\lambda }\right)}^2={[\mathcal{E}(r^2+a^2)-a{L}_z]}^2+\mathrm{\Delta }[\kappa r^2+-(L_z-a\mathcal{E}{)}^2-Q]\\ & {\left(\frac{\mathrm{d}\theta }{\mathrm{d}\lambda }\right)}^2={\cos }^2\theta [{\mathcal{E}}^2{a}^2+\kappa a^2-\frac{L_z^2}{1-{\cos }^2\theta }]+Q\end{array}$$

为了方便后续展开，对角向方程做如下变换： $\mathrm{d}\cos \theta =-\sin \theta \mathrm{d}\theta$ ，则：

$$\begin{array}{rl}{\left(\frac{\mathrm{d}\cos \theta }{\mathrm{d}\lambda }\right)}^2& =Q-[Q-a^2({\mathcal{E}}^2+\kappa )+L_z^2]{\cos }^2\theta -a^2({\mathcal{E}}^2+\kappa ){\cos }^4\theta \\ & =\mathrm{\Theta }(\cos \theta )\\ {\left(\frac{\mathrm{d}r}{\mathrm{d}\lambda }\right)}^2& =R(r)\end{array}$$

已知讨论的是束缚椭圆轨道，因此存在运动转折点，径向体现为近心点 $r_p$ 与远心点 $r_a$ ，角向体现为升降点： ${\theta }_{min},\pi -{\theta }_{min}$ ，由于径向与角向的分离性，可以独立对待两个运动，因此在各自的转折点均应该有：

$$\begin{array}{rl}& {\frac{\mathrm{d}\cos \theta }{\mathrm{d}\lambda }|}_{{\theta }_{min},\pi -{\theta }_{min}}=0\\ & {\frac{\mathrm{d}r}{\mathrm{d}\lambda }|}_{r_a,r_p}=0\end{array}$$

这是转折点，运动方程在此发生变号，这也正是体现在平方号上。转折点的存在使得运动方程被迫分段，通常只能在半个周期内进行计算，需要手动改变正负号才能够计算整个周期的运动。

对于 $t,\phi$ 方向的运动，由于径向与角向的可分离性因此也可以写为径向部分与角向部分，进行这般处理后便可以得到点源的测地线方程：

$$\begin{array}{rl}& {\left(\frac{\mathrm{d}\cos \theta }{\mathrm{d}\lambda }\right)}^2=\mathrm{\Theta }(\cos \theta )=Q-[Q-a^2({\mathcal{E}}^2+\kappa )+L_z^2]{\cos }^2\theta -a^2({\mathcal{E}}^2+\kappa ){\cos }^4\theta \\ & {\left(\frac{\mathrm{d}r}{\mathrm{d}\lambda }\right)}^2=R(r)={[\mathcal{E}(r^2+a^2)-a{L}_z]}^2+\mathrm{\Delta }[\kappa r^2+-(L_z-a\mathcal{E}{)}^2-Q]\\ & \frac{\mathrm{d}t}{\mathrm{d}\lambda }=T_r+T_{\theta }=\frac{\mathcal{E}(r^2+a^2{)}^2}{\mathrm{\Delta }}-\frac{(r^2+a^2)a{L}_z}{\mathrm{\Delta }}+a{L}_z-a^2\mathcal{E}(1-{\cos }^2\theta )\\ & \frac{\mathrm{d}\phi }{\mathrm{d}\lambda }={\mathrm{\Phi }}_r+{\mathrm{\Phi }}_{\theta }=\frac{a\mathcal{E}(r^2+a^2)-a^2{L}_z}{\mathrm{\Delta }}-a\mathcal{E}+\frac{L_z}{1-{\cos }^2\theta }\end{array}$$

束缚轨道的Kerr类时测地运动的数值计算的办法由考虑对 $\theta ,r$ 进行变量替换消除远点与近点的变号问题而解决，而解析的办法直至2012年才由 Fujita & hikita 给出，为不失一般性，接下来将从数值方法开始介绍。

转折点的处理

不管是数值计算还是解析计算，对于转折点的处理本质上是相同的，均是将束缚轨道的振荡特性等价转换为周期函数的增加。
对 $\theta$ 方向，由于 $\theta \in [{\theta }_{min},\pi -{\theta }_{min}]$ ，因此 $\cos \theta$ 仅能在 $\cos {\theta }_{min}-0-\cos (\pi -{\theta }_{min})-0-\cos {\theta }_{min}$ 来回振荡，尽管变量是周期的，但是由于 $\cos {\theta }_{min}<1$ 使得函数只能在半个周期的区间内来回振荡，变好正是对应了”折返“；而考虑变换： $\cos {\theta }_{min}\cos \chi =\cos \theta$ ，则可以解出限制，新变量 $\chi \in [0,2\pi ]$ 且对应有 $\cos \chi$ : $1\to 0\to -1\to 0\to 1\to \cdots$ ，自然发现对于 N 个周期的运动，对应着 $\chi \in [0,N\pi ]$ 。至于运动方向，有：

$$\begin{array}{rl}& \mathrm{d}\chi /\mathrm{d}\theta =\frac{\sin \theta }{\cos {\theta }_{min}\sin \chi }\\ \\ & \frac{\mathrm{d}\theta }{\mathrm{d}\lambda }=\frac{\mathrm{d}\theta }{\mathrm{d}\chi }\frac{\mathrm{d}\chi }{\mathrm{d}\lambda }=\pm \sqrt{\mathrm{\Theta }}\end{array}\$\$\mathrm{变}\mathrm{换}\mathrm{前}，\mathrm{由}\mathrm{于}\mathrm{运}\mathrm{动}\mathrm{方}\mathrm{向}\mathrm{会}\mathrm{发}\mathrm{生}\mathrm{改}\mathrm{变}，\mathrm{而}\mathrm{始}\mathrm{终}\mathrm{有}\$\sqrt{\mathrm{\Theta }}>0\$，\mathrm{因}\mathrm{此}\mathrm{需}\mathrm{要}\mathrm{进}\mathrm{行}\mathrm{变}\mathrm{号}\mathrm{来}\mathrm{保}\mathrm{证}\mathrm{运}\mathrm{动}\mathrm{方}\mathrm{程}，\mathrm{变}\mathrm{换}\mathrm{后}\mathrm{可}\mathrm{以}\mathrm{发}\mathrm{现}：+\mathrm{当}\$\theta :{\theta }_{min}\to \pi -{\theta }_{min}\$，\$\mathrm{d}\theta /\mathrm{d}\lambda >0,\mathrm{d}\theta /\mathrm{d}\chi >0\$+\mathrm{当}\$\theta :\pi -{\theta }_{min}\to {\theta }_{min}\$，\$\mathrm{d}\theta /\mathrm{d}\lambda <0,\mathrm{d}\theta /\mathrm{d}\chi <0\$\mathrm{可}\mathrm{以}\mathrm{看}\mathrm{到}\$\mathrm{d}\theta /\mathrm{d}\chi \$\mathrm{将}\mathrm{使}\mathrm{得}\mathrm{运}\mathrm{动}\mathrm{方}\mathrm{程}\mathrm{始}\mathrm{终}\mathrm{可}\mathrm{以}\mathrm{保}\mathrm{持}\mathrm{符}\mathrm{号}\mathrm{不}\mathrm{发}\mathrm{生}\mathrm{改}\mathrm{变}（\mathrm{或}\mathrm{说}\mathrm{已}\mathrm{经}\mathrm{自}\mathrm{动}\mathrm{编}\mathrm{码}\mathrm{了}\mathrm{变}\mathrm{号}）。\mathrm{对}\mathrm{于}\$r\$\mathrm{方}\mathrm{向}，\mathrm{可}\mathrm{以}\mathrm{用}\mathrm{椭}\mathrm{圆}\mathrm{轨}\mathrm{道}\mathrm{常}\mathrm{见}\mathrm{的}\mathrm{形}\mathrm{式}\mathrm{来}\mathrm{代}\mathrm{换}：\$\$r=\frac{pM}{1-e\cos \psi }$$

可见 $(r_p,r_a)\to (\psi =0,\psi =\pi )$ ，同样有：

$$\begin{array}{r}\mathrm{d}r/\mathrm{d}\psi =\frac{pM}{(1-e\cos \psi {)}^2}e\sin \psi \end{array}$$

可见变号过程已经被编码在了 $\sin \psi$ ，对于一个周期： $\psi \in [0,2\pi ]$ 。

这两个变换在数值计算使用，而对于解析表达，则使用的是具有周期性的椭圆函数及其反函数来编码此类 “往返” 周期性运动，其中的道理是相同的，在后面将揭示这一点。

运动频率

傅里叶展开

由于径向与角向的分离性，且束缚轨道的周期特性，只需要考虑半个周期的运动便可以得到轨道的频率信息，因此直接有：

$$\begin{array}{rl}{\mathrm{\Lambda }}_r& =2{\int }_{r_p}^{r_a}\frac{\mathrm{d}r}{\sqrt{R}}=\frac{2\pi }{{\mathrm{{\rm Y}}}_r}\\ {\mathrm{\Lambda }}_{\theta }& =4{\int }_{0=\cos \frac{\pi }{2}}^{\cos {\theta }_{min}}\frac{\mathrm{d}\cos \theta }{\sqrt{\mathrm{\Theta }}}=\frac{2\pi }{{\mathrm{{\rm Y}}}_{\theta }}\end{array}$$

对于 $t,\phi$ 运动，由于 $r,\theta$ 的独立性，可以进行如下傅里叶展开：

$$\begin{array}{rl}\frac{\mathrm{d}t}{\mathrm{d}\lambda }& =T_r+T_{\theta }+a{L}_z=T(r,\theta )=\sum _{n,k}{T}_{nk}{e}^{-\mathrm{i}(n{\mathrm{{\rm Y}}}_r+k{\mathrm{{\rm Y}}}_{\theta })\lambda }\\ T_{nk}& =\frac{1}{(2\pi {)}^2}{\int }_0^{2\pi }\mathrm{d}{\omega }_r{\int }_0^{2\pi }\mathrm{d}{\omega }_{\theta }T(r,\theta )e^{\mathrm{i}(n{\omega }_r+k{\omega }_{\theta })}\end{array}$$

其中相位变量 ${\omega }_{·}={\mathrm{{\rm Y}}}_{·}\lambda$ ，展开系数 $T_{nk}$ 是一个确定的实数（有意义），由傅里叶基的正交归一性：

$$\begin{array}{r}{\int }_0^{2\pi }{e}^{\mathrm{i}n\omega }\mathrm{d}\omega =\{\begin{array}{c}0,n\ne 0\\ \\ 2\pi ,n=0\end{array}\end{array}$$

展开系数写为：

$$\begin{array}{r}{T}_{nk}=\frac{1}{(2\pi {)}^2}[{\int }_0^{2\pi }{T}_r{e}^{\mathrm{i}n{\omega }_r}\mathrm{d}{\omega }_r{\int }_0^{2\pi }{e}^{\mathrm{i}k{\omega }_{\theta }}\mathrm{d}{\omega }_{\theta }+{\int }_0^{2\pi }{e}^{\mathrm{i}n{\omega }_r}\mathrm{d}{\omega }_r{\int }_0^{2\pi }{T}_{\theta }{e}^{\mathrm{i}k{\omega }_{\theta }}\mathrm{d}{\omega }_{\theta }+a{L}_z{\int }_0^{2\pi }{e}^{\mathrm{i}n{\omega }_r}\mathrm{d}{\omega }_r{\int }_0^{2\pi }{e}^{\mathrm{i}k{\omega }_{\theta }}\mathrm{d}{\omega }_{\theta }]\end{array}$$

因此系数具有对称性：

$$\begin{array}{r}{T}_{-n,-k}={\bar{T}}_{nk}=T_{nk},T_{00}\ne 0,T_{0k},T_{n0}\ne 0\end{array}$$

注意由于周期性，因此自然有 $T_{-n,-k}=T_{nk}$ ，这直接等价于取 $-\omega$ ，由此容易得到：

$$\begin{array}{r}\sum _{n,k}^{\pm \mathrm{\infty }}{T}_{nk}{e}^{-\mathrm{i}(n{\omega }_r+k{\omega }_{\theta })}=\sum _{n,k=0}^{\mathrm{\infty }}[T_{nk}{e}^{-\mathrm{i}(n{\omega }_r+k{\omega }_{\theta })}+{\bar{T}}_{nk}{e}^{\mathrm{i}(n{\omega }_r+k{\omega }_{\theta })}]\end{array}$$

以及非零项表达式：

$$\begin{array}{rl}{T}_{00}& =\frac{1}{2\pi }{\int }_0^{2\pi }{T}_r[r({\omega }_r)]\mathrm{d}{\omega }_r+\frac{1}{2\pi }{\int }_0^{2\pi }{T}_{\theta }[\theta ({\omega }_{\theta })]\mathrm{d}{\omega }_{\theta }+a{L}_z\\ T_{n0}& =\frac{1}{2\pi }{\int }_0^{2\pi }{T}_r[r({\omega }_r)]e^{\mathrm{i}{\omega }_{r}n}\mathrm{d}{\omega }_r,T_{0k}=\frac{1}{2\pi }{\int }_0^{2\pi }{T}_{\theta }[\theta ({\omega }_{\theta })]e^{\mathrm{i}{\omega }_{\theta }k}\mathrm{d}{\omega }_{\theta }\end{array}$$

因此运动方程化为：

$$\begin{array}{rl}\mathrm{d}t& =T_{00}+[\sum _{n=1}^{+\mathrm{\infty }}{T}_{n0}{e}^{-\mathrm{i}{\omega }_{r}n}+\sum _{k=1}^{+\mathrm{\infty }}{T}_{0k}{e}^{-\mathrm{i}{\omega }_{\theta }k}+c.c.]\mathrm{d}\lambda \\ t& =T_{00}\lambda +[\sum _{n=1}^{+\mathrm{\infty }}\frac{\mathrm{i}{T}_{n0}}{n{\mathrm{{\rm Y}}}_r}{e}^{-\mathrm{i}n{\mathrm{{\rm Y}}}_r}+c.c.]+[\sum _{k=1}^{+\mathrm{\infty }}\frac{\mathrm{i}{T}_{0k}}{k{\mathrm{{\rm Y}}}_{\theta }}{e}^{-\mathrm{i}k{\mathrm{{\rm Y}}}_{\theta }}+c.c.]\\ & =\mathrm{\Gamma }\lambda +\mathrm{\Delta }t=\mathrm{\Gamma }\lambda +\mathrm{\Delta }{t}^{(r)}+\mathrm{\Delta }{t}^{(\theta )}\end{array}$$

其中 $\mathrm{\Gamma }$ 代表一个周期增长的相位，是时间频率， $\mathrm{\Delta }t$ 是周期运动过程中的小振动产生的相位差；最后注意到：

$$\begin{array}{r}\frac{1}{2\pi }{\int }_0^{2\pi }{T}_r[r({\omega }_r)]\mathrm{d}{\omega }_r=\frac{{\mathrm{{\rm Y}}}_r}{2\pi }{\int }_0^{{\mathrm{\Lambda }}_r}\frac{T_r}{\sqrt{R}}\mathrm{d}\lambda ={⟨T_r⟩}_{\lambda }=\frac{2}{{\mathrm{\Lambda }}_r}{\int }_0^{{\mathrm{\Lambda }}_r/2}\frac{T_r}{\sqrt{R}}\mathrm{d}\lambda \end{array}$$

因此也有表达式：

$$\begin{array}{r}\mathrm{\Gamma }=T_{00}={⟨T_r⟩}_{\lambda }+{⟨T_{\theta }⟩}_{\lambda }+a{L}_z=a{L}_z+{\mathrm{{\rm Y}}}_{t^{(r)}}+{\mathrm{{\rm Y}}}_{t^{(\theta )}}\end{array}$$

可令：

$$\begin{array}{rl}\frac{\mathrm{d}t}{\mathrm{d}\lambda }& ={⟨T_r⟩}_{\lambda }+{⟨T_{\theta }⟩}_{\lambda }+a{L}_z+\frac{\mathrm{d}{t}^{(r)}}{\mathrm{d}\lambda }+\frac{\mathrm{d}{t}^{(\theta )}}{\mathrm{d}\lambda }\\ \frac{\mathrm{d}{t}^{(r)}}{\mathrm{d}\lambda }& =\frac{\mathrm{d}\mathrm{\Delta }{t}^{(r)}}{\mathrm{d}\lambda }=T_r-{⟨T_r⟩}_{\lambda }\\ t^{(r)}& ={\int }^r\frac{T_r-{⟨T_r⟩}_{\lambda }}{\sqrt{R(r^{\prime })}}\mathrm{d}{r}^{\prime }\end{array}$$

对于 $\phi$ 方向完全一样。至此就可以开始数值计算了，但距离得到解析解还需要用椭圆积分来表示前面出现的几个积分。

角向与径向频率解析表达

对于解析计算，首先解析求解角向频率，令 $z=\cos \theta$ ：

$$\begin{array}{rl}\frac{\mathrm{d}\cos \theta }{\mathrm{d}\lambda }& =\pm \sqrt{\mathrm{\Theta }}=(z^2-z_{-}^2)[a^2(1-{\mathcal{E}}^2)z^2-a^2(1-{\mathcal{E}}^2)z_{+}^2]\\ & =(z^2-z_m^2)[\beta z^2-z_p^2]=(z^2-z_m^2)(\beta z^2-z_p^2)\end{array}$$

其中 $z_m^2=\cos {\theta }_{min},z_p^2=\beta z_{+}^2=\frac{Q}{z_{-}^2}=\beta +\frac{L_z^2}{1-z_m^2}$ ，需要注意，转折点仅对应 $z=z_{-}={\cos }^2{\theta }_{min}$ 这两个角度，另外两个根 $z=z_{+}$ 没有物理意义，一般而言是虚数或者大于0； 再令 $z=z_m{y}_{\theta }$ ，则角向周期可表为：

$$\begin{array}{rl}{\mathrm{\Lambda }}_{\theta }& =4{\int }_0^{z_m}\frac{\mathrm{d}z}{z_m{z}_p\sqrt{(1-\frac{z^2}{z_m^2})(1-\beta \frac{z^2}{z_p^2})}}=\frac{4}{z_m{z}_p}{\int }_0^1\frac{z_m\mathrm{d}{y}_{\theta }}{z_m{z}_p\sqrt{(1-y^2)(1-\beta \frac{z_m^2}{z_p^2}{y}_{\theta }^2)}}\\ & =\frac{4}{z_p}{\int }_0^{\frac{\pi }{2}}\frac{\mathrm{d}\psi }{\sqrt{1-k_{\theta }^2{\sin }^2\varphi }}=\frac{4}{z_p}K(k_{\theta }^2)\end{array}$$

其中 $K(m)$ 为第一类完全椭圆积分（其中 $m_r=k_r^2$ ，后面均用 $m_{r/\theta }$ 表示），倒数第二个等号使用了： $y_{\theta }=\sin \varphi$ ， $m=k_{\theta }^2=\beta \frac{z_m^2}{z_p^2}$ 。
对于径向同理，考虑如下代换：

$$\begin{array}{rl}{y}_r=\sqrt{\frac{r-r_2}{r-r_3}\frac{r_1-r_3}{r_1-r_2}},r_1& =\frac{pM}{1-e},r_2=\frac{pM}{1-e},r_3=\frac{A+B+\sqrt{(A+B{)}^2-4AB}}{2},r_4=\frac{AB}{r_3}\\ A+B& =\frac{2M}{1-{\mathcal{E}}^2}-(r_1+r_2),AB=\frac{a^{2}Q}{(1-{\mathcal{E}}^2)r_1{r}_2}\end{array}$$

以及将运动方程改写为（ $r_1>r_2>r_3>r_4$ ）：

$$\begin{array}{r}\frac{\mathrm{d}r}{\mathrm{d}\lambda }=\sqrt{R}=(1-{\mathcal{E}}^2)(r_1-r)(r-r_2)(r-r_2)(r-r_4)\end{array}$$

同样需要注意，按这个写法实际上有 $(r_2,r_1),(r_4,r_3)$ 两个束缚轨道区间，但是一般而言 $r_4$ 在视界内部，因此并不是有意义的物理束缚轨道；与角向积分相同，径向频率可以写为：

$$\begin{array}{rl}{\mathrm{\Lambda }}_r=2{\int }_{r_2}^{r_1}\frac{\mathrm{d}r}{\sqrt{R}}& =2\frac{2}{\sqrt{(1-{\mathcal{E}}^2)(r_1-r_3)(r_2-r_4)}}{\int }_0^1\frac{\mathrm{d}{y}_r}{\sqrt{(1-y_r^2)(1-k_r^2{y}_r^2)}},k_r^2=\frac{r_3-r_4}{r_2-r_4}\frac{r_1-r_2}{r_1-r_3}\\ & =\frac{2}{\sqrt{(1-{\mathcal{E}}^2)(r_1-r_3)(r_2-r_4)}}{\int }_0^{\frac{\pi }{2}}\frac{\mathrm{d}\psi }{\sqrt{1-k_r^2{\sin }^2\varphi }},y_r=\sin \varphi \\ & =\frac{2}{\sqrt{(1-{\mathcal{E}}^2)(r_1-r_3)(r_2-r_4)}}K(m_r)\end{array}$$

时间与轴向方向频率的解析表达

而 $t,\phi$ 方向的频率为：

$$\begin{array}{rl}\mathrm{\Gamma }& =\frac{2}{{\mathrm{\Lambda }}_r}{\int }_{r_2}^{r_1}\frac{T_r}{\sqrt{R}}\mathrm{d}r+\frac{4}{\mathrm{\Lambda }}{\int }_0^{\cos {\theta }_{min}}\frac{T_{\theta }}{\sqrt{\mathrm{\Theta }}}\mathrm{d}\cos \theta +a{L}_z={\mathrm{{\rm Y}}}_{t^{(r)}}+{\mathrm{{\rm Y}}}_{t^{(\theta )}}+a{L}_z\\ {\mathrm{{\rm Y}}}_{\phi }& =\frac{2}{{\mathrm{\Lambda }}_r}{\int }_{r_2}^{r_1}\frac{{\mathrm{\Phi }}_r}{\sqrt{R}}\mathrm{d}r+\frac{4}{\mathrm{\Lambda }}{\int }_0^{\cos {\theta }_{min}}\frac{{\mathrm{\Phi }}_{\theta }}{\sqrt{\mathrm{\Theta }}}\mathrm{d}\cos \theta -a\mathcal{E}={\mathrm{{\rm Y}}}_{{\phi }^{(r)}}+{\mathrm{{\rm Y}}}_{{\phi }^{(\theta )}}-a\mathcal{E}\end{array}$$

角向部分只存在 $z$ 的一次项，容易积分，直接计算即可，结果为：

$$\begin{array}{rl}{\mathrm{{\rm Y}}}_{t^{(\theta )}}& =\frac{2a^2\mathcal{E}{\mathrm{{\rm Y}}}_{\theta }}{\pi \beta z_p}[(z_p^2-\beta )K(m_{\theta })-z_p^{2}E(m_{\theta })]\\ {\mathrm{{\rm Y}}}_{{\phi }^{(\theta )}}& =\frac{2L_z{\mathrm{{\rm Y}}}_{\theta }}{\pi z_p}\mathrm{\Pi }(\frac{\pi }{2},-z_m^2,m_{\theta })\end{array}$$

其中 $E(m),\mathrm{\Pi }(\frac{\pi }{2},c,m)$ 依次为第二、三类完全椭圆积分；对于径向部分则由于二次项的存在而复杂得多，因此需要利用多项式展开将 $r^2$ 降阶为 $r$ 或 $r^{-1}$ ，具体的可以展开为：

$$\begin{array}{rl}{T}_r& =\mathcal{E}{r}^2+2Mr\mathcal{E}+(a^2+4M^2)\mathcal{E}-a{L}_z+\frac{2M}{r_{+}-r_{-}}\{\frac{(4M^2\mathcal{E}-a{L}_z)r_{+}-2M^2{a}^2\mathcal{E}}{r-r_{+}}-(+\leftrightarrow -)\},a<M\\ & =4M^2\mathcal{E}+\frac{2M(4M^2\mathcal{E}-a{L}_z)}{r-M}+\frac{2M^2(2M^2\mathcal{E}-a{L}_z)}{(r-M{)}^2},a=M\\ \\ {\mathrm{\Phi }}_r& =a\mathcal{E}+\frac{a}{r_{+}-r_{-}}(\frac{2Ma\mathcal{E}{r}_{+}-a{L}_z}{r-r_{+}}-\frac{2Ma\mathcal{E}{r}_{-}-a{L}_z}{r-r_{-}}),a<M\\ & =a\mathcal{E}+\frac{2Ma\mathcal{E}}{r-M}+\frac{2M^{2}a\mathcal{E}-a{L}_z}{(r-M{)}^2},a=M\end{array}$$

由此径向积分由以下几个积分组成：

$$\begin{array}{r}{\int }_{r_2}^{r_1}\frac{1}{\sqrt{R}}\mathrm{d}r,{\int }_{r_2}^{r_1}\frac{r}{\sqrt{R}}\mathrm{d}r,{\int }_{r_2}^{r_1}\frac{r^2}{\sqrt{R}}\mathrm{d}r,{\int }_{r_2}^{r_1}\frac{1}{(r-r_{\pm })\sqrt{R}}\mathrm{d}r,{\int }_{r_2}^{r_1}\frac{1}{(r-r_{\pm }{)}^2\sqrt{R}}\mathrm{d}r\end{array}$$

第一项就是径向周期，第二项与第四项为一阶项，也可以直接计算得到：

$$\begin{array}{rl}& {\int }_{r_2}^{r_1}\frac{r}{\sqrt{R}}\mathrm{d}r=\frac{2(r_2-r_3)}{\sqrt{(1-{\mathcal{E}}^2)(r_1-r_3)(r_2-r_4)}}\mathrm{\Pi }(\frac{\pi }{2},-h_r,k_r^2)+\frac{2r_3}{\sqrt{(1-{\mathcal{E}}^2)(r_1-r_3)(r_2-r_4)}}K(k_{}^{2}r)\\ & {\int }_{r_2}^{r_1}\frac{1}{(r-r_{\pm })\sqrt{R}}\mathrm{d}r=\frac{2}{\sqrt{(1-{\mathcal{E}}^2)(r_1-r_3)(r_2-r_4)}}[K(k_r^2)-\frac{r_2-r_3}{r_2-r_{\pm }}\mathrm{\Pi }(\frac{\pi }{2},-h_{\pm },k_r^2)]\\ & h_r=\frac{r_1-r_2}{r_1-r_3},h_{\pm }=h_r\frac{r_3-r_{\pm }}{r_2-r_{\pm }}\end{array}$$

对于二阶项，已知通过变量替换有：

$$\begin{array}{r}\sqrt{R}\sim \sqrt{\psi }=\sqrt{(1-y_r^2)(1-k_r^2{y}_r^2)}\end{array}$$

我们希望用三椭圆积分来表示目标积分，而第三类椭圆积分在特殊情况下可以转化为另外两类，而第三类椭圆积分为：

$$\begin{array}{r}\mathrm{\Pi }(\varphi ,-h,m)={\int }_0^{\varphi }\frac{\mathrm{d}y}{(1-h{y}^2)\sqrt{\psi }}\end{array}$$

由于面对的被积函数是二阶的，可以预计会出现二次项的分母 $\sim (1-h{y}^2{)}^2$ ，事实上正是如此，我们可以展开 $r^2$ ：

$$\begin{array}{rl}& r^2=r_3^2+\frac{2r_3(r_2-r_3)}{1-h^r{y}_r^2}+\frac{(r_2-r_3{)}^2}{(1-h_r{y}_r^2{)}^2}\end{array}$$

注意到其中可以进一步展开：

$$\begin{array}{r}\frac{1}{(1-h_r{y}_r^2{)}^2}=\frac{A}{(y_r-h_r^{-\frac{1}{2}}{)}^2}+\frac{B}{(y_r+h_r^{-\frac{1}{2}}{)}^2}\end{array}$$

这样分母便被降阶为了： $y_r^2$ 项，已经接近第三类椭圆积分了；最终便有：

$$\begin{array}{r}{r}^2=r_3^2+\frac{(r_2-r_3)(r_2+3r_3)}{2}\frac{1}{1-h^r{y}_r^2}+\frac{(r_2-r_3{)}^2}{4h_r}[\frac{1}{(y_r-h_r^{-\frac{1}{2}}{)}^2}+\frac{1}{(y_r+h_r^{-\frac{1}{2}}{)}^2}]\end{array}$$

对于 $(r-M{)}^{-2}$ 同理：

$$\begin{array}{rl}& {\left(\frac{r_3-M}{r-M}\right)}^2=\\ & 1+\frac{1}{2}(4-\frac{r_2-M}{r_2-M})\left(\frac{r_2-r_3}{r_2-M}\right)\frac{1}{1-h^r{y}_r^2}+\frac{(\frac{r_2-r_3}{r_2-M}{)}^2}{4h_r}[\frac{1}{(y_r-h_r^{-\frac{1}{2}}{)}^2}+\frac{1}{(y_r+h_r^{-\frac{1}{2}}{)}^2}]\end{array}$$

接下来的问题就是对展开式的最后一项积分，恰恰发现在椭圆积分的递推关系有一个一般表达式：

$$\begin{array}{r}{J}_n(c)=\int \frac{\mathrm{d}y}{(y-c{)}^n\sqrt{\psi (y)}}\end{array}$$

因此我们需要求的是： $J_2[c]+J_2[-c]$ ，这实际上是求： $\int \frac{\mathrm{d}y}{(y^2-c^2{)}^2\sqrt{\psi }}$ ，而这可以通过构造：

$$\begin{array}{r}\frac{\mathrm{d}}{\mathrm{d}y}\left(\frac{y\sqrt{\psi }}{(y^2-c^2)}\right)=\frac{\mathrm{d}}{\mathrm{d}y}\left(\frac{y}{(y^2-c^2)}\sqrt{\psi }\right)=\frac{\mathrm{d}}{\mathrm{d}y}\left(fg\right)\end{array}$$

来解决，经过分部积分便可以导出等式：

$$\begin{array}{rl}{J}_2[c]+J_2[-c]& =\\ & \frac{2}{\sqrt{\psi (c)}}\{[k_r^2(2c^2-1)-1]\mathrm{\Pi }(\varphi ,-c^{-2},m)+(1-c^2{k}_r^2)F(\varphi ,m)-E(\varphi ,m)-{\frac{y\sqrt{\psi }}{y^2-c^2}|}_0^y\}\end{array}$$

最终借此可以给出 ${\mathrm{{\rm Y}}}_{t^{(r)}},{\mathrm{{\rm Y}}}_{{\phi }^{(r)}}$ ，由此最终给出（ $a<M$ ）频率的解析表达式：

$$\begin{array}{rl}\mathrm{\Gamma }& =4M^2\mathcal{E}+\frac{2a^2\mathcal{E}{\mathrm{{\rm Y}}}_{\theta }{z}_p}{\pi \beta }[K(m_{\theta })-E(m_{\theta })]\\ & +\frac{2{\mathrm{{\rm Y}}}_r}{\pi \sqrt{(1-{\mathcal{E}}^2)(r_1-r_3)(r_2-r_4)}}\{2M\mathcal{E}[r_{3}K(m_r)+(r_2-r_3)\mathrm{\Pi }(\frac{\pi }{2},-h_r,m_r)]\\ & +\frac{\mathcal{E}}{2}[(r_3(r_1+r_2+r_3)-r_1{r}_2)K(m_r)+(r_1-r_3)(r_2-r_4)E(m_r)\\ & +(r_2-r_3)(r_1+r_2+r_3+r_4)\mathrm{\Pi }(\frac{\pi }{2},-h_r,m_r)]\\ & +[\frac{2M}{(r_{+}-r_{-})}((4M^2\mathcal{E}-a{L}_z)r_{+}-2M^2{a}^2\mathcal{E})(\frac{K(m_r)}{r_3-r_{-}}-\frac{(r_2-r_3)\mathrm{\Pi }(\frac{\pi }{2},-h_{+},m_r)}{(r_2-r_{+})(r_3-r_{+})})-(+\leftrightarrow -)]\}\\ \\ {\mathrm{{\rm Y}}}_{\phi }& =\frac{2L_z{\mathrm{{\rm Y}}}_{\theta }}{\pi z_p}\mathrm{\Pi }(\frac{\pi }{2},-z_m^2,m_{\theta })\\ & +\frac{2{\mathrm{{\rm Y}}}_r}{\pi \sqrt{(1-{\mathcal{E}}^2)(r_1-r_3)(r_2-r_4)}}\frac{a}{r_{+}-r_{-}}[\frac{2Ma\mathcal{E}{r}_{+}-a{L}_z}{r_3-r_{+}}(K(m_r)-\frac{r_2-r_3}{r_2-r_{+}}\mathrm{\Pi }(\frac{\pi }{2},-h_{+},m_r))\\ & -(+\leftrightarrow -)]\\ \\ {\mathrm{{\rm Y}}}_r& =\frac{\pi \sqrt{(1-{\mathcal{E}}^2)(r_1-r_3)(r_2-r_4)}}{2K(m_r)}\\ \\ {\mathrm{{\rm Y}}}_{\theta }& =\frac{\pi z_p}{2K(m_{\theta })}\end{array}$$

由于篇幅太长暂时就不考虑 $a=M$ 情况，日后再做补充。

运动轨迹

接下来可以求解运动轨迹的解析表达式，首先需要的是选定一个演化起始点，存在两个选择：上行区与下行区，由两个转折点分开，对应为： $\frac{\mathrm{d}r}{\mathrm{d}\lambda }>0,\frac{\mathrm{d}r}{\mathrm{d}\lambda }<0$ 。我们暂时只展示起始点位于上行区的情况，因为下行区与上行区本质上只存在一个相位差。此外由于周期性，我们只需要讨论最后一个周期内的运动就可以完成整个系统全时间运动的描述。

径向运动

选取初始点： $r_0^{(1)}$ ，其中 $(1)$ 代表选择的上行区的初值点。对于上行区，有 $sgn(L_z)>0,{\frac{\mathrm{d}r}{\mathrm{d}\lambda }|}_{\lambda =0}$ ，从运动方程可以得到累计的Mino时间为：

$$\begin{array}{r}\lambda ={\lambda }^{(r)}={\int }_{r_0^{(1)}}^r\frac{\mathrm{d}{r}^{\prime }}{\sqrt{R}}\end{array}$$

其中 ${\lambda }^{(r)}$ 用于标记这是由径向运动标记的时间；考虑运动的顺序为： $r_0^{(1)}\to r_1\to r_2\to r_0^{(1)}$ ，因此分三段来表示上述积分：

$$\begin{array}{r}{\lambda }^{(r)}=\{\begin{array}{c}({\int }_{r_2}^r-{\int }_{r_2}^{r_0^{(1)}})\frac{\mathrm{d}{r}^{\prime }}{\sqrt{R}},r_0^{(1)}\to r_1\\ (-{\int }_{r_2}^r+2{\int }_{r_2}^{r_1}-{\int }_{r_2}^{r_0^{(1)}})\frac{\mathrm{d}{r}^{\prime }}{\sqrt{R}},r_1\to r_2\\ ({\int }_{r_2}^r+2{\int }_{r_2}^{r_1}-{\int }_{r_2}^{r_0^{(1)}})\frac{\mathrm{d}{r}^{\prime }}{\sqrt{R}},r_2\to r_0^{(1)}\end{array}\end{array}$$

现在令：

$$\begin{array}{rl}{\lambda }_0^{(r)}(r)& =\frac{2}{\sqrt{(1-{\mathcal{E}}^2)(r_1-r_3)(r_2-r_4)}}F(\varphi ,m_r)=\frac{2}{\sqrt{(1-{\mathcal{E}}^2)(r_1-r_3)(r_2-r_4)}}{u}_r(r)={\int }_{r_2}^r\frac{\mathrm{d}{r}^{\prime }}{\sqrt{R}}\\ {\mathrm{\Lambda }}_0^{(1)}& ={\int }_{r_2}^{r_0^{(1)}}\frac{\mathrm{d}{r}^{\prime }}{\sqrt{R}},\alpha =\frac{2}{\sqrt{(1-{\mathcal{E}}^2)(r_1-r_3)(r_2-r_4)}}\end{array}$$

则可以将Mino时间用一系列周期函数表示为：

$$\begin{array}{r}{\lambda }^{(r)}=\{\begin{array}{c}{\lambda }_0^{(r)}-{\mathrm{\Lambda }}_r^{(1)},r_0^{(1)}\to r_1\\ -{\lambda }_0^{(r)}+{\mathrm{\Lambda }}_r-{\mathrm{\Lambda }}_r^{(1)},r_1\to r_2\\ {\lambda }_0^{(r)}+{\mathrm{\Lambda }}_r-{\mathrm{\Lambda }}_r^{(1)},r_2\to r_0^{(1)}\end{array}\end{array}$$

又从径向频率可以有：

$$\begin{array}{r}{\mathrm{\Lambda }}_r=2\alpha K(m_r)\Rightarrow \alpha =\frac{{\mathrm{\Lambda }}_r}{2K}\end{array}$$

因此可以定义一个新函数：

$$\begin{array}{r}{u}_r=F(\varphi ,m_r)=\frac{2K}{{\mathrm{\Lambda }}_r}{\lambda }_0^{(r)}\end{array}$$

当然 $u_r$ 也是分段表出的，将时间的分段表述反解出来即可，在此就先不多赘述。引入 $u_r$ 最重要的一个作用是得到坐标 $r$ 的解析表达，首先可以发现， $u_r$ 表征了从初始点开始积累的相位：

$$\begin{array}{r}{u}_r=\frac{K}{\pi }{\mathrm{{\rm Y}}}_r{\lambda }_0^{(r)}=\frac{K}{\pi }{q}_r\end{array}$$

可以通过椭圆函数的反函数来提取坐标表达式：

$$\begin{array}{rl}& u_r=F(\arctan y_r,m_r)\Rightarrow \arctan y_r=\mathrm{am}[u_r,m_r]\\ \\ & \Rightarrow y_r=\sqrt{\frac{r-r_2}{r-r_3}\frac{r_1-r_3}{r_1-r_2}}=\mathrm{sn}(u_r,m_r)=\mathrm{sn}(\frac{K}{\pi }{q}_r,m_r)\\ \\ & r(\lambda )=\frac{r_3(r_1-r_2){\mathrm{sn}}^2-r_2(r_1-r_3)}{(r_1-r_2){\mathrm{sn}}^2-(r_1-r_3)}\end{array}$$

其中 $q_r(\lambda )$ 便是相位，我们可以直接以此为新变量来计算运动轨迹。这里也可以看到， $\mathrm{sn}$ 函数是周期性的，其完整编码了径向的周期性往返运动！

角向运动

角向运动与径向运动类似，同样考虑上行区的起始点： ${\theta }_0^{(1)}$ ， ${\frac{\mathrm{d}\cos \theta }{\mathrm{d}\lambda }|}_{\lambda =0}>0$ ；
运动轨迹为： ${\theta }_0^{(1)}\to 0\to {\theta }_{min}\to 0\to \pi -{\theta }_{min}\to {\theta }_0^{(1)}$ ，同样由角向运动计量的Mino时间为：

$$\begin{array}{r}{\lambda }^{(\theta )}=\{\begin{array}{c}({\int }_0^{\cos \theta }-{\int }_0^{\cos {\theta }_0^{(1)}})\frac{\mathrm{d}\cos {\theta }^{\prime }}{\sqrt{\mathrm{\Theta }}},{\theta }_0^{(1)}\to {\theta }_{min}\\ (-{\int }_0^{\cos \theta }+2{\int }_0^{\cos {\theta }_{min}}-{\int }_0^{\cos {\theta }_0^{(1)}})\frac{\mathrm{d}\cos {\theta }^{\prime }}{\sqrt{\mathrm{\Theta }}},{\theta }_{min}\to \pi -{\theta }_{min}\\ ({\int }_0^{\cos \theta }+4{\int }_0^{\cos {\theta }_{min}}-{\int }_0^{\cos {\theta }_0^{(1)}})\frac{\mathrm{d}\cos {\theta }^{\prime }}{\sqrt{\mathrm{\Theta }}},\pi -{\theta }_{min}\to {\theta }_0^{(1)}\end{array}\end{array}$$

同样令：

$$\begin{array}{rl}& {\lambda }_0^{(\theta )}={\int }_0^{\cos \theta }\frac{\mathrm{d}\cos {\theta }^{\prime }}{\sqrt{\mathrm{\Theta }}}=\frac{F(\arcsin y_{\theta },m_{\theta })}{z_p},{\mathrm{\Lambda }}_{\theta }^{(1)}={\int }_0^{\cos {\theta }_0^{(1)}}\frac{\mathrm{d}\cos {\theta }^{\prime }}{\sqrt{\mathrm{\Theta }}}\end{array}$$

因此也有：

$$\begin{array}{r}{\lambda }^{(\theta )}=\{\begin{array}{c}{\lambda }_0^{(\theta )}-{\mathrm{\Lambda }}_{\theta }^{(1)},{\theta }_0^{(1)}\to {\theta }_{min}\\ -{\lambda }_0^{(\theta )}+\frac{{\mathrm{\Lambda }}_{\theta }}{2}-{\mathrm{\Lambda }}_{\theta }^{(1)},{\theta }_{min}\to \pi -{\theta }_{min}\\ {\lambda }_0^{(\theta )}+{\mathrm{\Lambda }}_{\theta }-{\mathrm{\Lambda }}_{\theta }^{(1)},\pi -{\theta }_{min}\to {\theta }_0^{(1)}\end{array}\end{array}$$

同样的可以通过引入函数 $u_{\theta }$ 来推导 $\cos \theta$ 的表达式：

$$\begin{array}{rl}& u_{\theta }=F(\arcsin y_{\theta },m_{\theta })=z_p{\lambda }_0^{(\theta )}=\frac{2K}{\pi }{\mathrm{{\rm Y}}}_{\theta }{\lambda }_0^{(\theta )}=\frac{2K}{\pi }{q}_{\theta }\\ \\ & y_{\theta }=\frac{\cos \theta }{\cos {\theta }_{min}}=\mathrm{sn}(u_{\theta },m_{\theta })\Rightarrow \cos \theta =\cos {\theta }_{min}\mathrm{sn}(\frac{2K}{\pi }{q}_{\theta },m_{\theta })\end{array}$$

初值选取

最后是关于初值的选取，首先需要说明的是，椭圆积分是允许相位超过 $\frac{\pi }{2}$ 的，对于一个周期而言，需要椭圆积分的积分相位为： $r:(0,\pi ),\theta :(0,2\pi )$ ，又因为椭圆积分本身的被积函数是大于等于0且周期为 $\pi$ 的且关于极值点对称的函数，因此椭圆积分本身具有 $\frac{\pi }{2}$ 周期。因此恰好可以描述所需要的轨道相位演化。

其次，因为坐标时间是统一的，因此有：

$$\begin{array}{rl}& \lambda ={\lambda }^{(r)}={\lambda }^{(\theta )}\end{array}$$

为了简单起见，选取转折点 $r_2,{\theta }_{min}$ 为初始点，由此便不必刻意区分上行或者下行，进一步取径向时间作为基准，则有：

$$\begin{array}{rl}& \lambda ={\lambda }^{(r)}={\lambda }_0^{(r)}\Rightarrow q_r={\mathrm{{\rm Y}}}_r{\lambda }_0^{(r)}={\mathrm{{\rm Y}}}_r\lambda \\ \\ & \lambda ={\lambda }^{(\theta )}={\lambda }_0^{(\theta )}-\frac{{\mathrm{\Lambda }}_{\theta }}{4}\\ \\ & \Rightarrow q_{\theta }={\mathrm{{\rm Y}}}_{\theta }{\lambda }_0^{(\theta )}={\mathrm{{\rm Y}}}_{\theta }\lambda +\frac{\pi }{2}\end{array}$$

由此确定了数值计算时径向与角向的相位变量。如果初始点不在转折点，那么便是从 $(r_0,{\theta }_0)$ 开始计时，变换到上述使用转折点开始只需要减去一个沿着轨道方向的相对 $(r_1,{\theta }_{min})$ 的相位差即可。

时间与轴向运动

从前面的傅里叶展开出发：

$$\begin{array}{rl}t& =T_{00}\lambda +[\sum _{n=1}^{+\mathrm{\infty }}\frac{\mathrm{i}{T}_{n0}}{n{\mathrm{{\rm Y}}}_r}{e}^{-\mathrm{i}n{\mathrm{{\rm Y}}}_r}+c.c.]+[\sum _{k=1}^{+\mathrm{\infty }}\frac{\mathrm{i}{T}_{0k}}{k{\mathrm{{\rm Y}}}_{\theta }}{e}^{-\mathrm{i}k{\mathrm{{\rm Y}}}_{\theta }}+c.c.]\\ & =\mathrm{\Gamma }\lambda +\mathrm{\Delta }t=\mathrm{\Gamma }\lambda +\mathrm{\Delta }{t}^{(r)}+\mathrm{\Delta }{t}^{(\theta )}=\mathrm{\Gamma }\lambda +t^{(r)}+t^{(\theta )}\\ \\ \phi & ={\mathrm{\Phi }}_{00}\lambda +[\sum _{n=1}^{+\mathrm{\infty }}\frac{\mathrm{i}{\mathrm{\Phi }}_{n0}}{n{\mathrm{{\rm Y}}}_r}{e}^{-\mathrm{i}n{\mathrm{{\rm Y}}}_r}+c.c.]+[\sum _{k=1}^{+\mathrm{\infty }}\frac{\mathrm{i}{\mathrm{\Phi }}_{0k}}{k{\mathrm{{\rm Y}}}_{\theta }}{e}^{-\mathrm{i}k{\mathrm{{\rm Y}}}_{\theta }}+c.c.]\\ & ={\mathrm{{\rm Y}}}_{\phi }\lambda +\mathrm{\Delta }\phi =\mathrm{{\rm Y}}\lambda +\mathrm{\Delta }{\phi }^{(r)}+\mathrm{\Delta }{\phi }^{(\theta )}=\mathrm{{\rm Y}}\lambda +{\phi }^{(r)}+{\phi }^{(\theta )}\end{array}$$

结合上一节对 $r,\cos \theta$ 解析表达式的求解，显然发现

$$\begin{array}{rl}& {\omega }_r=q_r,{\omega }_{\theta }=q_{\theta }+\frac{\pi }{2}\end{array}$$

因此至此就可以对 $t,\phi$ 运动数值求解了，但是为了得到解析表达，还需要进一步处理

$$\begin{array}{rl}& t^{(r)}={\int }_{r_0}^r\frac{T_r-{\mathrm{{\rm Y}}}_{t^{(r)}}}{\sqrt{R(r^{\prime })}}\mathrm{d}{r}^{\prime }={\int }_{r_0}^r\frac{T_r}{\sqrt{R(r^{\prime })}}\mathrm{d}{r}^{\prime }-\frac{2}{{\mathrm{\Lambda }}_{r_1}}{\int }_{r_2}^r\frac{T_r}{\sqrt{R(r^{\prime })}}\mathrm{d}{r}^{\prime }{\int }_{r_0}^r\frac{\mathrm{d}{r}^{\prime }}{\sqrt{R(r^{\prime })}}\\ & t^{(r)}={\int }_{r_0}^r\frac{T_r-{\mathrm{{\rm Y}}}_{t^{(r)}}}{\sqrt{R(r^{\prime })}}\mathrm{d}{r}^{\prime }={\int }_{r_0}^r\frac{T_r}{\sqrt{R(r^{\prime })}}\mathrm{d}{r}^{\prime }-\frac{2}{{\mathrm{\Lambda }}_{r_1}}{\int }_{r_2}^r\frac{T_r}{\sqrt{R(r^{\prime })}}\mathrm{d}{r}^{\prime }{\int }_{r_0}^r\frac{\mathrm{d}{r}^{\prime }}{\sqrt{R(r^{\prime })}}\end{array}$$

对于 $\phi$ 同理，经过复杂运算可以给出：

$$\begin{array}{rl}{t}^{(r)}& =\frac{\alpha \mathcal{E}}{2}(r_1-r_3)(r_2-r_4){\tilde{E}}_r(\varphi ,-h_r,m_r)+\frac{\alpha \mathcal{E}}{2}(r_2-r_3)[4M+(r_1+r_2+r_3+r_4)]{\tilde{\mathrm{\Pi }}}_r(\varphi ,-h_r,m_r)\\ & -\alpha [\frac{2M(2M{r}_{+}^2\mathcal{E}-a{L}_z{r}_{+})(r_2-r_3)}{(r_{+}-r_{-})(r_3-r_{+})(r_2-r_{+})}{\tilde{\mathrm{\Pi }}}_r(\varphi ,-h_{+},m_r)-(+\leftrightarrow -)]\\ {\phi }^{(r)}& =-\frac{a\alpha }{r_{+}-r_{-}}[\frac{2Ma\mathcal{E}{r}_{+}-a{L}_z}{r_3-r_{+}}\frac{r_2-r_3}{r_2-r_{+}}{\tilde{\mathrm{\Pi }}}_r(\varphi ,-h_{+},m_r)-(+\leftrightarrow -)]\\ t^{(\theta )}& =-\frac{a^2\mathcal{E}{z}_p}{\beta }{\tilde{E}}_{\theta }(\varphi ,m_{\theta }),{\phi }^{(\theta )}=\frac{L_z}{z_p}\widetilde{(}\varphi ,-z_m^2,m_{\theta })\end{array}$$

其中（以上行区的分段积分为例子，对于数值计算时取转折点为起始点则可以直接套用）：

$$\begin{array}{rl}& {\tilde{E}}_r(\varphi ,-h_r,m_r)=E_r(\varphi ,-h_r,m_r)-\frac{{\mathrm{{\rm Y}}}_r{\lambda }^{(r)}}{\pi }E(m_r)\\ & E_r(\varphi ,-h_r,m_r)={\int }_{{\varphi }_0}^{\varphi }\sqrt{1-m_r{\sin }^2{\varphi }^{\prime }}\mathrm{d}{\varphi }^{\prime }+{\frac{y_r\sqrt{\psi }}{y_r^2-h_r^{-1}}|}_{y_0}^{y_r}\\ & =(E(\varphi ,m_r)+{\frac{y_r\sqrt{\psi }}{y_r^2-h_r^{-1}}|}_0^{y_r})-(E({\varphi }_0,m_r)+{\frac{y_r\sqrt{\psi }}{y_r^2-h_r^{-1}}|}_0^{y_0})\\ & =E_r^{(0)}(\varphi ,-h_r,m_r)-E_r^{(0)}({\varphi }_0,-h_r,m_r),r_0\to r_1\\ & =2E_r^{(0)}(\frac{\pi }{2},-h_r,m_r)-E_r^{(0)}(\varphi ,-h_r,m_r)-E_r^{(0)}({\varphi }_0,-h_r,m_r),r_1\to r_2\\ & =2E_r^{(0)}(\frac{\pi }{2},-h_r,m_r)-E_r^{(0)}(\varphi ,-h_r,m_r)+E_r^{(0)}({\varphi }_0,-h_r,m_r),r_2\to r_0\\ & \\ & {\tilde{\mathrm{\Pi }}}_r(\varphi ,-h_{r/\pm },m_r)={\mathrm{\Pi }}_r(\varphi ,-h_{r/\pm },m_r)-\frac{{\mathrm{{\rm Y}}}_r{\lambda }^{(r)}}{\pi }{\mathrm{\Pi }}_r(\varphi ,-h_{r/\pm },m_r)\\ & {\mathrm{\Pi }}_r(\varphi ,-h_{r/\pm },m_r)=\mathrm{\Pi }(\varphi ,-h_{r/\pm },m_r)-\mathrm{\Pi }({\varphi }_0,-h_{r/\pm },m_r),r_0\to r_1\\ & =2\mathrm{\Pi }(\frac{\pi }{2},-h_{r/\pm },m_r)-\mathrm{\Pi }(\varphi ,-h_{r/\pm },m_r)-\mathrm{\Pi }({\varphi }_0,-h_{r/\pm },m_r),r_1\to r_2\\ & =2\mathrm{\Pi }(\frac{\pi }{2},-h_{r/\pm },m_r)-\mathrm{\Pi }(\varphi ,-h_{r/\pm },m_r)+\mathrm{\Pi }({\varphi }_0,-h_{r/\pm },m_r),r_2\to r_0\\ & \\ & {\tilde{E}}_{\theta }(\varphi ,m_{\theta })=E_{\theta }(\varphi ,m_{\theta })-\frac{2{\mathrm{{\rm Y}}}_{\theta }}{\pi }{\lambda }^{(\theta )}E(\frac{\pi }{2},m_{\theta })\\ & E_{\theta }(\varphi ,m_{\theta })=E(\varphi ,m_{\theta })-E({\varphi }_0,m_{\theta }),{\theta }_0\to {\theta }_{min}\\ & =-E(\varphi ,m_{\theta })+2E(\frac{\pi }{2},m_{\theta })-E_{(}{\varphi }_0,m_{\theta }),{\theta }_{min}\to \pi -{\theta }_{min}\\ & =E(\varphi ,m_{\theta })+4E(\frac{\pi }{2},m_{\theta })-E_{(}{\varphi }_0,m_{\theta }),\pi -{\theta }_{min}\to {\theta }_0\\ & {\tilde{\mathrm{\Pi }}}_{\theta }(\varphi ,-z_m^2,m_{\theta })={\mathrm{\Pi }}_{\theta }(\varphi ,-z_m^2,m_{\theta })-\frac{2{\mathrm{{\rm Y}}}_{\theta }}{\pi }{\lambda }^{(\theta )}\mathrm{\Pi }(\frac{\pi }{2},-z_m^2,m_{\theta })\\ & {\mathrm{\Pi }}_{\theta }(\varphi ,-z_m^2,m_{\theta })=\mathrm{\Pi }(\varphi ,-z_m^2,m_{\theta })-\mathrm{\Pi }({\varphi }_0,-z_m^2,m_{\theta }),{\theta }_0\to {\theta }_{min}\\ & =-\mathrm{\Pi }(\varphi ,-z_m^2,m_{\theta })+2\mathrm{\Pi }(\frac{\pi }{2},-z_m^2,m_{\theta })-{\mathrm{\Pi }}_{(}{\varphi }_0,-z_m^2,m_{\theta }),{\theta }_{min}\to \pi -{\theta }_{min}\\ & =\mathrm{\Pi }(\varphi ,-z_m^2,m_{\theta })+4\mathrm{\Pi }(\frac{\pi }{2},-z_m^2,m_{\theta })-{\mathrm{\Pi }}_{(}{\varphi }_0,-z_m^2,m_{\theta }),\pi -{\theta }_{min}\to {\theta }_0\end{array}$$

至此已经完成对于类时束缚轨道测地线的运动解析求解（以及提供了数值求解的理论基础）。
对于数值计算法，理论上有：直接四阶龙格库塔求解四个线性方程或者更粗暴直接数值积分，以及采取本文所讲述的计算傅里叶级数的半解析积分。

从观测量到守恒量

直接从运动方程出发，以及转折点的构造给出观测量 $(p,e,\cos \iota )$ 与 守恒量 $(\hat{E},{\hat{L}}_Z,\hat{Q})$ 之间的映射。
其中

$$\begin{array}{rl}& a=\frac{\hat{a}}{M},\mathcal{E}=\frac{\hat{E}}{\mu },L_z=\frac{\widehat{L_z}}{\mu M},Q=\frac{\hat{Q}}{{\mu }^2{M}^2}\end{array}$$

对于倾斜角 $\iota$ 定义为在赤道为 0，沿着转动方向时为 九十度 ，大于 九十度 时标志着运动方向与黑洞自旋方向相反，即逆行；因此有：

$$\begin{array}{rl}& \iota =\frac{\pi }{2}-sgn(L_z){\theta }_{min}\end{array}$$

且有 :

$$\begin{array}{rl}& \mathrm{\Theta }=0\Rightarrow Q={\cos }^2{\theta }_{min}(\beta +\frac{L_z^2}{1-{\cos }^2{\theta }_{min}})\end{array}$$

再带入径向方程，考虑转折点处的方程以及圆轨道情况，可以解出：

$$\begin{array}{rl}& {\mathcal{E}}^2=\frac{\kappa \rho +2ϵ\sigma +2(sgn{L}_z)\sqrt{\sigma (\sigma {ϵ}^2+\rho \kappa ϵ-\eta {\kappa }^2)}}{{\rho }^2+4\eta \sigma }\\ & L_z=-\frac{g_1\mathcal{E}}{h_1}+(sgn{L}_z)\sqrt{\frac{g_1^2{\mathcal{E}}^2}{h_1^2}+\frac{f_1{\mathcal{E}}^2-d_1}{h_1}}\end{array}$$

具体系数暂时懒得写了，反正考虑两个等式： $\frac{\mathrm{d}r}{\mathrm{d}\lambda }=0,\frac{{\mathrm{d}}^{2}r}{{\mathrm{d}}^2\lambda }=0$ 联立求解即可，对于偏心轨道，采用第一个等式分别在近远点的情况；对于用圆轨道，两个等式联立求解，详情见参考文献。

snapshot 的绘制（采取BlackHoleToolKit）

数值计算是EMRI除运动轨迹建模外另一个难点，计算的步骤为：

• 计算测地线方程，首先给出轨道四个方向的频率，其次给出径向与角向运动 $r,\theta$ 表达式，最终给出 $\mathrm{\Delta }t$ ， $\mathrm{\Delta }\phi$ 。
• 构造Teukolsky方程齐次解，获取无穷远/视界的渐近振幅
• 构造源项表达式，计算二重卷积积分，获得振幅
• 多模态叠加得到EMRI波形
  测地线的难度主要在理论建模，数值计算上并没有太大阻碍，真正困难的在于二重积分的计算，该积分的被积函数是振荡的，由运动轨迹和 Teukolsy 方程齐次解构成，因此齐次解的计算也直接决定了二重积分的计算，例如，如果采取MST方法，那么二重积分的计算速度将远慢于采取数值直接求解Teukolsky方程（并行计算下至少慢十倍）；
  计算设置：

• BlackHoleToolKit程序包实现，调用了 Teukolsky 、 KerrGeodesic 、 SpinWeightedSpheroidalHarmonics 。
• 轨道信息调用

• orbit = KerrGeodesicOrbit[a,p,e,x,"Method"->"Analytical"]
    
  1
  2
  3
  4
  5
  6
  7
  8
  	即采取前面讲述过的解析法，尽管计算轨道的时候更慢，但是在给出给定轨道的各类参数时计算要更快。
  - 调用 Teukolsky 包其中内置了 对于二重积分的计算，采取的是 ``` Tracy ``` 数值积分法；且集成了径向齐次方程与自旋加权椭谐函数的调用，官方程序包默认调用的齐次解为采用 Heun 方程的数值求解法： ``` Numerical  ``` ，不过我们修改了程序，调用了 ```MST``` 进行了同样的计算；注意，径向齐次解均已经对振幅进行了归一化。由此可以获得单模态振幅：
  - ```mathematica
    \[Psi] = TeukolskyPointParticleMode[s, l, m, n, k, orbit]]
    Z = \[Psi]["Amplitudes"]["\[ScriptCapitalI]"]
    S = \[Psi]["AngularFunction"][\[Pi]/3, 0]
    \[Omega] = \[Psi]["\[Omega]"]
    h_{+}-ih_{\times} = 2/D Z/\[Omega]^2 S Exp[-I \[Omega] (t-r_{*})]
  
  

• 多模态叠加：对 l,n,k 进行截断，通常截断在当固定其他参数，目标参数连续几位的振幅都保持衰减且均低于设置阈值，其中 l 衰减的比较快，一般采取截断 l=10,11 ，对于 n 其是径向运动产生的指标，对于偏心率越大的轨道，需要将其截断到更高的数值，因为高偏心率导致 远点与近点相距更远，辐射强度量级差距大，二重积分容易发散，这可以追溯到径向方程的格林函数解。对于 k 是角向运动的指标，一般对其截断要求较低，因为角向频率较低且角向运动产生的引力波辐射并不占主导。
  本文对于低偏心率轨道（e=0.1,0.3）采取 l=5,n=10,k=5 截断，高偏心率轨道（e=0.5）采取 l=5,n=15,k=5 截断，需要注意，理论上需要计算的模态数有： $(2l+1)(l+1)(2n+1)(2k+1)-1$ ，减 1 是因为量子数均为 0 时频率也为 0 ，没有意义。因此我们的参数配置理论上有 7388(低偏心);10908(高偏心) 个模态参与叠加。
• 结果展示（基于 ai agent 协助构建代码，EMRI波形的计算成本很大，即便采取并行计算以及计算速度最快的数值齐次解方法，在个人电脑上全部的16个核并行，以下结果也一共需要六小时左右，如果采取 MST 法那单个低偏心率也起码算几个小时）：
  
  
  

Note：目前最新进展有JH方法（Rust）、Heun数值计算法以及广义SN方法（Julia），下一步势必需要自己从BlackHoleToolKit中重构出一个专属的整合了不同齐次解解法的计算程序包。此外二重积分计算的收敛性问题也需要解决，BlackHoleToolKit程序内置的二重积分反而在计算角向振荡积分时难以达到设定收敛精度。

下一步

辐射反应——轨道衰减

能流计算

测地线修正

SFG 理论

Effective-One-Body

Lense

考虑到的几个潜在小项目：

• 天体物理环境下的 EMRI ：
  超大质量黑洞最常见的出现位置是星系中心，而ANG是一个得到天文观测证明的天体，吸积盘的物质密度比较高，具有丰富的恒星形成与演化过程，是天然的EMRI产生场所。

  尽管越靠近克尔视界，物质密度会降低，直到越过最小稳定轨道；因此直接把吸积盘作为一个摩擦流体考虑的话其对EMRI的影响较小（尤其是对于旋近后期），但是这没有考虑引力；如果考虑吸积盘对背景时空的形变，情况或许会不一样。

  因此可以考虑一个有吸积盘形变的时空度规，以便全相对论的考虑吸积盘对EMRI的影响，更近一步的或许可以研究其对双星并和信号的影响，以此来估计目前事件发生在吸积盘上的比例。

  同理暗物质尖峰/晕。

  对于非真空轴对称 Type-D 假设下，采取 NP 形式具有度规解：

$$\begin{array}{rl}& \mathrm{d}{s}^2=A(r,\theta )[\frac{L(r)}{\mathrm{\Sigma }}{(\mathrm{d}t-a{\sin }^2\theta \mathrm{d}\phi )}^2-\frac{\mathrm{\Sigma }}{L(r)}\mathrm{d}{r}^2\\ & -\mathrm{\Sigma }\mathrm{d}{\theta }^2-\frac{(r^2+a^2{)}^2{\sin }^2\theta }{\mathrm{\Sigma }}{(\mathrm{d}\phi -\frac{a}{r^2+a^2}\mathrm{d}t)}^2]\end{array}$$

对应的引力扰动的 Teukolsky 方程为：

$$\begin{array}{rl}& [\frac{(r^2+a^2{)}^2}{L(r)}-a^2{\sin }^2\theta ]\frac{{\partial }^2{}_{-2}\mathrm{\Psi }}{\partial t^2}+[\frac{2(r^2+a^2)}{L(r)}\frac{\partial L(r)}{\partial r}-8r-4\mathrm{i}a\cos \theta ]\frac{\partial {}_{-2}\mathrm{\Psi }}{\partial t}\\ & +[\frac{a^2}{L(r)}-\frac{1}{{\sin }^2\theta }]\frac{{\partial }^2{}_{-2}\mathrm{\Psi }}{\partial {\phi }^2}-L(r{)}^2\frac{\partial }{\partial r}(L(r{)}^{-1}\frac{\partial {}_{-2}\mathrm{\Psi }}{\partial r})-\frac{1}{\sin \theta }\frac{\partial }{\partial \theta }(\sin \theta \frac{\partial {}_{-2}\mathrm{\Psi }}{\partial \theta })\\ & +[\frac{2a{L}^{\prime }(r)}{L(r)}+\frac{4\mathrm{i}\cos \theta }{{\sin }^2\theta }]\frac{\partial {}_{-2}\mathrm{\Psi }}{\partial \phi }+[\frac{2a(r^2+a^2)}{L(r)}-2a]\frac{{\partial }^2{}_{-2}\mathrm{\Psi }}{\partial t\partial \phi }\\ & +(4{\cot }^2\theta +2){}_{-2}\mathrm{\Psi }=4\pi \mathrm{\Sigma }\hat{T}\end{array}$$

Note：该方程的导出需要进行分离规范选择，而Kerr情况的度规并不需要。

近期 Meng Kun el.ta 2026 研究的DM EMRI 所用的形变度规正是此形式。

评估：主要难点在于需要从度规开始全流程计算，实施成本较大（预估三个月），且物质对时空的形变较小，公式复杂，数值上对精度的需求可能比较大。

对于吸积盘情况就更为复杂，目前主要还在球对称黑洞薄盘模型，例如 `Che-Yu Chen 2023 :

$$\begin{array}{rl}& \mathrm{d}{s}^2=-f(r)e^{2{\nu }_{disk}}\mathrm{d}{t}^2+e^{2{\lambda }_{ext}-2{\nu }_{disk}}\frac{\mathrm{d}{r}^2}{f(r)}+r^2{e}^{-2{\nu }_{disk}}(e^{2{\lambda }_{ext}}\mathrm{d}{\theta }^2+{\sin }^2\theta \mathrm{d}{\phi }^2)\end{array}$$

不过对于球对称度规的话比较好处理，但是考虑到吸积盘结构的复杂性，预估与考虑DM的难度会差不多。

• 有效单体描述 b-EMRI
  目前关于相对论建模的 b-EMRI 采取的 Dixon 形式来描述作为扰动源项的内双星，通过一系列坐标变换从内双星平坦时空的运动映射到背景的BL时空，以此获得可以计算的源项能动张量，但是Dixon形式描述的是具有内部结构的扩展体运动，需要假设时空渐近平坦，且需要展开到足够高阶（八极矩）以揭露内双星运动的特性。

  但是可以发现有一个共性：先描述内双星轨道运动，通过坐标变换到背景时空BL系，以此获取内双星单个天体的轨迹，问题就在于如何更好的描述内双星运动。

  近期有EOB理论可以描述具有极端质量比的双星系统的 IMR 过程，因此可以考虑使用Kerr背景下的EOB运动描述内双星，采取质心在BL系做测地运动，再通过坐标变换得到内双星天体的运动轨迹（四速度，时空坐标），以此代替点粒子能动张量中的运动量来给出单个天体的源项，而内双星的源项由两个黑洞的线性叠加即可。这样还有一个好处，可能可以用EOB来描述内双星的 IMR 过程。

  评估：四月份太极数据训练营听韩文标老师的报告，他们也在想办法去弄可以 IMR 的b-EMRI，EOB 国内有谁能做的过韩老师组啊

• 多信使

  1. 辐射反应会导致测地偏运动，假设是带有磁场（中子星？）或者电（气体云？）的天体旋近，在辐射反应加剧的时候可能产生可以观测到的电磁信号。
     不过看到的几篇文献都是解 爱因斯坦-麦克斯韦方程来考虑电荷的。

  2. 黑洞-中子星旋近会发出X射线，bnb-EMRI 组合或许会有一些结果?

• 记忆效应在EMRI中的长期叠加效应
  引力波辐射中会伴随一种非振荡的长期/永久的应变，在BBH时间中或许难以观测，但是如果考虑EMRI的长期轨道演化，是否会对其相位与波形产生足以观测到的永久性偏离？
  补充：侯绍齐老师说效应太弱，对测试粒子轨迹的修正是个很小的量，基本不考虑。

---

• Next Post →