低频准周期振荡是存在于黑洞双星暂现源的普遍现象，表现为X射线光变功率谱密度曲线上的宽峰，峰的中心频率从几毫赫兹到30 Hz^[1-3]。根据峰的品质因子、噪声和相位延迟等，低频准周期振荡通常分为A，B和C 3种。其中，C型低频准周期振荡是黑洞暂现源中最常见的类型，常在低硬态和硬中间态观测到，频率范围在几毫赫兹到10 Hz，但在高软态也能观测到，频率可达到约30 Hz。C型低频准周期振荡伴随着较强的平台噪声，连续谱呈现近似截断幂律谱轮廓，且通常由一个主峰和几个附加峰组成复杂的结构。主峰的中心频率与振荡的基频有关，附加峰是指分别在基频的半频和二次谐频处出现的副谐频峰和谐频峰^[2-5]。在少数情况下，谐频峰(或副谐频峰)的值与基频峰相当，形成一种特殊的双峰低频准周期振荡。如在XTE J1550-564和XTE J1859+226两个黑洞暂现源功率谱的相关频率处出现了两个强而窄的峰^[6-7]。

虽然对低频准周期振荡的观测和研究取得了很大的进展，但是其物理起源一直是一个有争议的问题，其中，一个可能的起源是吸积盘的振荡。在吸积盘中，流体元受扰动的恢复力作用时，在水平方向作惯性振荡，在赤道面附近垂直方向作简谐振荡，且不同方向的振荡是耦合的。这些振荡按照恢复力的性质可以分为g-模(gravity modes)、p-模(pressure modes)和c-模(corrugation modes)3种类型。研究发现，这些振荡被禁闭在吸积盘中一定区域内传播，即所谓陷波振荡(trapped oscillation)^[8-9]。禁闭振荡的频率与黑洞质量成反比，可能是黑洞暂现源中观测到的准周期振荡。其中，g-模和p-模振荡频率在高频准周期振荡的频率范围内，而c-模是低频振荡，被禁闭在内半径和内部垂向共振半径之间，其振荡频率与低频准周期振荡频率相当^[2]。此外，文[10]提出了吸积盘全局的垂向振荡模型，该模型假设整个盘的垂向振荡可以产生0.1 Hz左右的持续低频准周期振荡。

以上吸积盘的振荡模型并不能从观测的光谱上拟合低频准周期振荡，且从C型低频准周期振荡的性质及宽带噪声的形状和强度可以看出，准周期振荡应该与含有附加噪声的随机系统有关。因此，文[11-12]提出用广义朗之万方程描述吸积盘的随机振荡，并通过数值计算得到振荡光度的功率谱密度，解释了一些观测结果。文[13]利用此模型计算了3种情况的光变曲线和功率谱，并通过幂律拟合功率谱来获得不同情况下的光谱指数，发现观测天体光谱指数变化是由盘振荡在不同情况之间的转换引起的。对于这样的随机系统，在一定条件下能产生随机共振现象。而随机共振是乘性在随机动力系统中的非线性现象，已广泛应用于天文领域的研究^[14]。文[15-17]讨论了吸积盘随机振荡光度功率谱密度曲线中的随机共振现象，认为曲线中的共振峰为黑洞双星系统中准周期振荡现象的物理解释。在这些工作中，模型中的随机项仅存在于粘滞系数和外部作用力，吸积盘的振荡频率固定不变，在功率谱密度曲线中没有发现谐频峰。基于此，本文提出用包含频率涨落和随机力作用的广义朗之万方程模型描述黑洞吸积盘的垂向振荡，计算吸积盘随机振荡光度的功率谱密度，讨论出现的双峰随机共振现象，为解释C型低频准周期振荡提供理论依据。

1 吸积盘随机振荡的功率谱密度 1.1 吸积盘振荡模型方程

一个包围致密天体的标准吸积薄盘在受到外界随机力和内部粘滞力作用时，盘中单位质量粒子的垂向振荡可以用广义朗之万方程描述^[12]：

$ \frac{\mathrm{d}^{2} z}{\mathrm{~d} t^{2}}+\xi \frac{\mathrm{d} z}{\mathrm{~d} t}+\omega_{\perp}^{2} z=f(t), $

(1)

其中，z和ω_⊥分别为粒子垂向振荡的位移和频率；ξ为吸积盘内部的粘滞耗散作用系数，是一个常量；外部作用力f(t)是指数形式关联的湍流作用随机力，满足统计性质

$ \langle f(t)\rangle=0,\left\langle f(t) f\left(t^{\prime}\right)\right\rangle=D{\mathrm{e}}^{-\lambda\left|t-t^{\prime}\right|}, $

(2)

其中，〈〉表示取平均值；D为常数，表示随机力作用的强度；λ为随机力关联时间的倒数，即时间关联率，其大小和盘振荡频率的关系影响吸积盘振荡与湍流作用的耦合^[18]。在此，我们设定λ小于盘的振荡频率。

假设吸积盘的垂向振荡频率受外界因素的影响存在涨落现象，其表达式为^[19]

$ \omega_{\perp}^{2}=\omega_{0}^{2}+\eta(t). $

(3)

频率涨落η(t)是高斯白噪声随机量，统计特性为

$ \langle\eta(t)\rangle=0,\left\langle\eta(t) \eta\left(t^{\prime}\right)\right\rangle=\varepsilon \omega_{0}^{2} \delta\left(t-t^{\prime}\right), $

(4)

其中，ε为噪声强度，设0≤ε≤1，因为若ε值较大，频率噪声会迫使功率谱密度曲线上的峰较宽，而降低峰的品质因子Q值^[18]；ω₀为吸积盘垂向振荡的特征角频率，表达式为^[10]

$ \omega_{0}=2 {\rm{ \mathsf{ π} }} \nu_{0}=\frac{4.4 {\rm{ \mathsf{ π} }} \times 10^{3}}{m}\left\{\frac{2-\gamma}{r_{\mathrm{in}}^{3}} \frac{1-\gamma /\left[(\gamma+1) r_{\mathrm{adj}}\right]}{r_{\mathrm{adj}}^{\gamma} r_{\mathrm{out}}^{2-\gamma}}\right\}^{1 / 2} \mathrm{~Hz}, $

(5)

其中，m=M/M_⊙，M和M_⊙分别为中心天体和太阳的质量；r_in为吸积盘内边缘的半径，r_in=R_in/3R_S；R_S为史瓦西半径，R_S=2GM/c²；r_out=R_out/R_in和r_adj=R_adj/R_in分别为盘的外边缘半径和调整半径；面密度指数γ可以取3/5或3/4^[10]。

为了计算简便，假定β=ξ/2作为吸积盘内部的粘滞阻尼系数，则(1)式可以写成

$ \frac{\mathrm{d}^{2} z}{\mathrm{~d} t^{2}}+2 \beta \frac{\mathrm{d} z}{\mathrm{~d} t}+\left[\omega_{0}^{2}+\eta(t)\right] z=f(t). $

(6)

假设两个随机量不相关联，满足〈f(t)η(t′)〉=0。通过对(6)式作拉普拉斯变换，可以求得任意时刻粒子振荡的位移和速度：

$ z(t)=\langle z(t)\rangle+\int_{0}^{t} H(t-s)[f(s)-\eta(s) z(s)] \mathrm{d} s, $

(7)

$ v(t)=\langle v(t)\rangle+\int_{0}^{t} \dot{H}(t-s)[f(s)-\eta(s) z(s)] \mathrm{d} s, $

(8)

其中，$ \dot H $表示$\frac{{{\rm{d}}H}}{{{\rm{d}}t}} $。假设初始位移和速度为z(0)=0，v(0)=v₀，则位移和速度的平均值〈z(t)〉和〈v(t)〉为

$ \langle z(t)\rangle=v_{0} H(t),\langle v(t)\rangle=v_{0} \dot{H}(t). $

(9)

在ω₀ > β的情况下, 响应函数H(t)和$ \dot H $(t)为

$ H(t)=\frac{\mathrm{e}^{-\beta t}}{\omega} \sin \omega t, \dot{H}(t)=\frac{\mathrm{e}^{-\beta t}}{\omega}(\omega \cos \omega t-\beta \sin \omega t), $

(10)

其中，$ \omega = \sqrt {\omega _0^2 - {\beta ^2}} $。

1.2 振荡位移和速度的关联函数

根据文[20]给出的广义朗之万方程中位移和速度关联函数的计算方法，定义位移的方差σ²(t)、及位移和速度的时间关联函数κ_z(t, τ)和κ_v(t, τ)为

$ \sigma^{2}(t)=\left\langle[z(t)-\langle z(t)\rangle]^{2}\right\rangle, $

(11)

$ \kappa_{z}(t, \quad \tau)=\langle[z(t+\tau)-\langle z(t+\tau)\rangle][z(t)-\langle z(t)\rangle]\rangle, $

(12)

$ \kappa_{v}(t, \quad \tau)=\langle[v(t+\tau)-\langle v(t+\tau)\rangle][v(t)-\langle v(t)\rangle]\rangle . $

(13)

根据随机量η(t)和f(t)统计过程的独立性，它们满足以下统计关系

$ \begin{gathered} \left\langle f\left(t_{1}\right) \eta\left(t_{2}\right) z\left(t_{2}\right)\right\rangle=\left\langle f\left(t_{1}\right) \eta\left(t_{1}\right) z\left(t_{2}\right)\right\rangle=0,\langle f(t) z(t)\rangle=0, \\ \left\langle z\left(t_{1}\right) z\left(t_{2}\right) \eta\left(t_{1}\right) \eta\left(t_{2}\right)\right\rangle=\varepsilon \omega_{0}^{2} \delta\left(t_{1}-t_{2}\right)\left\langle z^{2}\left(t_{2}\right)\right\rangle, \end{gathered} $

(14)

则位移和速度的关联函数可以表示为

$ \begin{aligned} \kappa_{z}(t, \tau)=& \int_{0}^{t} \int_{0}^{t+\tau} H\left(t_{1}\right) H\left(t_{2}+\tau\right)\left\langle f\left(t_{1}\right) f\left(t_{2}\right)\right\rangle \mathrm{d} t_{1} \mathrm{~d} t_{2}+\\ & \int_{0}^{t} H\left(t_{1}\right) H\left(t_{1}+\tau\right)\left\langle\eta\left(t_{1}\right) \eta(t)\right\rangle \sigma^{2}\left(t-t_{1}\right) \mathrm{d} t_{1}, \end{aligned} $

(15)

$ \begin{aligned} \kappa_{v}(t, \tau)=& \int_{0}^{t+\tau} \int_{0}^{t} \dot{H}\left(t_{1}\right) \dot{H}\left(t_{2}+\tau\right)\left\langle f\left(t_{1}\right) f\left(t_{2}\right)\right\rangle \mathrm{d} t_{1} \mathrm{~d} t_{2}+\\ & \int_{0}^{t} \dot{H}\left(t_{1}\right) \dot{H}\left(t_{1}+\tau\right)\left\langle\eta\left(t_{1}\right) \eta(t)\right\rangle \sigma^{2}\left(t-t_{1}\right) \mathrm{d} t_{1} . \end{aligned} $

(16)

当t→∞时，以上公式变为

$ \kappa_{z a}(\tau)=D \int_{0}^{\infty} \int_{0}^{\infty} H\left(t_{1}\right) H\left(t_{2}\right) \mathrm{e}^{-\lambda\left|t_{2}+\tau-t_{1}\right|} \mathrm{d} t_{1} \mathrm{~d} t_{2}+\bar{\sigma} \int_{0}^{\infty} H\left(t_{1}\right) H\left(t_{1}+\tau\right) \mathrm{d} t_{1}, $

(17)

$ \kappa_{v a}(\tau)=D \int_{0}^{\infty} \int_{0}^{\infty} \dot{H}\left(t_{1}\right) \dot{H}\left(t_{2}\right) \mathrm{e}^{-\lambda\left|t_{2}+\tau-t_{1}\right|} \mathrm{d} t_{1} \mathrm{~d} t_{2}+\bar{\sigma} \int_{0}^{\infty} \dot{H}\left(t_{1}\right) \dot{H}\left(t_{1}+\tau\right) \mathrm{d} t_{1}, $

(18)

其中，σ为渐近状态下位移方差的时间齐次部分，t→∞时，κ_za(0)=σ，则可以得到

$ \bar{\sigma}=\frac{2 D(2 \beta+\lambda)}{\omega_{0}^{2}(4 \beta-\varepsilon)\left(\omega_{0}^{2}+2 \beta \lambda+\lambda^{2}\right)}. $

(19)

1.3 随机振荡光度的功率谱密度

吸积盘中单位质量粒子振荡的总能量为

$ E=\frac{1}{2} v^{2}+\frac{1}{2} \omega_{0}^{2} z^{2}. $

(20)

由于外部随机力、内部粘滞耗散和频率涨落噪声的作用，粒子的总能量部分损失，损失的能量即为粒子振荡的光度L(t)，由(6)式和(20)式可以得到

$ L(t)=-\frac{\mathrm{d} E}{\mathrm{~d} t}=2 \beta v^{2}+\eta(t) z v-f(t) v. $

(21)

当t→∞时，光度的稳态平均自相关函数定义为

$ C_{L L}(\tau)=\lim \limits_{t \rightarrow \infty}\langle[L(t+\tau)-\langle L(t+\tau)\rangle][L(t)-\langle L(t)\rangle]\rangle. $

(22)

将(17), (18)和(21)式代入(22)式，得到

$ \begin{gathered} C_{L L}(\tau)= \lim \limits_{t \rightarrow \infty}\left[4 \beta^{2}\left\langle\varTheta^{2}(t+\tau)\right\rangle\left\langle\varTheta^{2}(t)\right\rangle+8 \beta^{2}\langle\varTheta(t+\tau) \varTheta(t)\rangle^{2}-2 \beta\left\langle\varTheta^{2}(t+\tau)\right\rangle\langle f(t) \varTheta(t)\rangle\right.\\ -4 \beta\langle\varTheta(t+\tau) \varTheta(t)\rangle\langle f(t) \varTheta(t+\tau)\rangle-2 \beta\langle f(t+\tau) \varTheta(t+\tau)\rangle\left\langle\varTheta^{2}(t)\right\rangle+\\ \langle f(t+\tau) \varTheta(t+\tau)\rangle\langle f(t) \varTheta(t)\rangle+\langle\eta(t+\tau) \eta(t)\rangle\langle\varTheta(t+\tau) \varTheta(t) \varPhi(t+\tau) \varPhi(t)\rangle]。\end{gathered} $

(23)

其中，Θ(t)=$ \int_0^t $$ \dot H $(t-t′)[f(t′)-η(t′)z(t′)]dt′；Φ(t)=$ \int_0^t $H(t-t′)[f(t′)-η(t′)z(t′)]dt′。

通过傅里叶变换，(23)式转换成振荡光度的功率谱密度

$ P S D(\nu)=\int_{-\infty}^{+\infty} C_{L L}(\tau) \mathrm{e}^{-\mathrm{i} 2 {\rm{ \mathsf{ π} }} \nu \tau} \mathrm{d} \tau. $

(24)

通过复杂的计算，并忽略不含频率ν的常数项，可以得到PSD(ν)的表达式为

$ \begin{aligned} {PSD}(\nu) &=\frac{\sqrt{2}}{\sqrt{{\rm{ \mathsf{ π} }}} a_{3}}\left[\frac{4 D^{2} \lambda^{2}\left(2 a_{2} a_{4}^{2}+a_{4} a_{7} a_{5}-a_{5}^{2} a_{6}\right)}{\omega^{2}\left(a_{1}^{2}-4 \beta^{2} \lambda^{2}\right)^{2}}+\frac{\varepsilon^{2} \omega_{0}^{4} \bar{\sigma}^{2}\left(2 a_{6} \omega^{2}-2 a_{6} \beta^{2}+3 a_{7} \beta \omega-4 a_{2} \beta^{2}\right)}{2 \omega^{2}}\right] \\ &+\frac{(\beta+\lambda)\left(a_{1}+2 \beta \lambda+4 {\rm{ \mathsf{ π} }}^{2} \nu^{2}\right)}{\sqrt{2 {\rm{ \mathsf{ π} }}} a_{8}}\left[\frac{D^{2} \lambda a_{5}}{\beta \omega\left(a_{1}-2 \beta \lambda\right)^{2}}+\frac{D \varepsilon^{2} \omega_{0}^{4} \bar{\sigma}}{2 \beta}\right]+\frac{\varepsilon^{2} \omega_{0}^{4} \bar{\sigma}^{2} \beta\left(\beta^{2}+\omega_{0}^{2}\right)}{\sqrt{2 {\rm{ \mathsf{ π} }}} \omega^{2}\left(4 \beta^{2}+4 {\rm{ \mathsf{ π} }}^{2} \nu^{2}\right)} \\ &+\frac{\omega\left(a_{1}+2 \beta \lambda-4 {\rm{ \mathsf{ π} }}^{2} \nu^{2}\right)}{\sqrt{2 {\rm{ \mathsf{ π} }}} a_{8}}\left[\frac{D^{2} \lambda a_{4}}{2 \beta \omega\left(a_{1}-2 \beta \lambda\right)^{2}}-\frac{D \varepsilon \omega_{0}^{2} \bar{\sigma}\left(a_{1}^{2}+12 \beta^{2} \lambda^{2}\right)}{2 \omega\left(a_{1}^{2}-4 \beta^{2} \lambda^{2}\right)}\right] \\ &+\frac{8 \sqrt{2} D^{2} \lambda^{3} a_{1}^{2}}{\sqrt{{\rm{ \mathsf{ π} }}}\left(4 \lambda^{2}+4 {\rm{ \mathsf{ π} }}^{2} \nu^{2}\right)\left(a_{1}^{2}-4 \beta^{2} \lambda^{2}\right)}, \end{aligned} $

(25)

其中，a₁=ω₀²+λ²；a₂=βω²(12π²ν²-4ω₀²)；a₃=[16π⁴ν⁴+32(β²-ω²)π²ν²+16ω₀⁴](4β²+4π²ν²)；a₄=β(ω₀²-λ²)；a₅=ω(ω₀²+λ²)；a₆=β[4(β²+π²ν²)²+24β²ω²-8ω²π²ν²+8ω⁴]；a₇=16ω(β²+π²ν²)(π²ν²-ω₀²)；a₈=(a₁+4π²ν²)²-16ω²π²ν²。

2 模型参数对功率谱密度曲线中低频准周期振荡的影响

根据解析表达式(25)，我们可以计算不同参数作用下盘振荡光度的功率谱密度，并发现曲线中出现了单峰和双峰的随机共振现象。这些随机共振峰即为低频准周期振荡。下面主要讨论系统参数(ε, λ, β, ν₀)对低频准周期振荡的影响。

2.1 低频准周期振荡随特征振荡频率ν₀的变化

在我们的模型中，通过计算发现随机力的强度D不改变谱线的轮廓，因此，在以下的分析中取D=1，且吸积盘垂向振荡的特征频率ν₀=ω₀/2π，近似为C型低频准周期振荡的基频。由(5)式可知，ω₀依赖于吸积盘的结构和中心天体的质量。选择不同的盘参数，使振荡频率在低频准周期振荡的范围内。在其他参数保持不变时，不同特征频率下的功率谱密度曲线均呈现与低频准周期振荡相吻合的低频截断幂律谱轮廓(图 1)。由图 1可以看出，当ν₀取不同值时，小于截断处频率的低频部分功率变化很小，并且高频部分的谱斜率基本保持不变。然而，曲线中的随机共振峰对ν₀十分敏感，只有当ν₀达到一定的值时，才出现一个基频峰和一个附加的二次谐频峰，双峰的中心频率分别对应ν₀和2ν₀。随着ν₀减小，两个峰变低变宽，且基频峰变化更快。当ν₀减小到0.07 Hz时，基频峰几乎消失，在曲线上只留下一个较宽的二次谐频峰。随着ν₀进一步减小，谐频峰也越来越弱，直至完全消失。在ν₀=0.02 Hz时，功率谱密度随频率的增加而单调减小，不存在低频准周期振荡。这说明吸积盘振荡的频率越高，由随机共振导致的准周期振荡现象越明显。

2.2 低频准周期振荡随ε, λ和β等参数的变化

由(21)式可以发现，吸积盘的振荡光度主要由粘滞耗散、随机力和频率涨落噪声共同作用产生，这3种作用在模型中分别用粘滞阻尼系数β、时间关联率λ和噪声强度ε表示。为了讨论各参数对低频准周期振荡的影响，我们取同文[10]一致的吸积盘几何参数，即r_in=1，r_out=10⁴，γ=3/5，r_adj=3，选择频率为0.1~10 Hz的黑洞候选体XTE J1550-564的质量作为天体质量^[7]。根据文[1]的评估，该候选体质量范围8.4M_⊙~10.8M_⊙，取m=9.6。将这些参数代入(5)式，计算得到吸积盘振荡的特征频率为ν₀=0.289 Hz。当ν₀固定不变，λ, ε和β取不同值时，振荡光度的功率谱密度曲线如图 2。由图 2可以看出，曲线上都出现了期望的单峰或双峰低频准周期振荡，且参数对谱线的影响各不相同。

图 2清晰地表明，当吸积盘的振荡频率固定时，低频准周期振荡基频峰的中心频率不随其他参数发生变化。这些参数只影响低频准周期振荡峰的高度和宽度。如随机力时间关联率λ对基频的峰影响较大，而对二次谐频峰几乎不产生影响(如图 2(a))，且随着λ的增大，曲线上移，说明随机力的关联时间越短，谱密度越大，但是基频对应的峰值却在减小，峰变得更宽更平，直至消失，此时随机共振峰只包含一个中心频率为2ν₀的单峰。相反，频率噪声强度ε和粘滞阻尼系数β只对二次谐频峰有较大的影响，对基频峰的影响却很小(如图 2(b)和(c))。当增大ε，或减小β，功率谱密度曲线上移，二次谐频峰变得更加尖锐。而ε取较小值，或β取较大值时，二次谐频峰都消失。以上说明随机力关联时间越短，频率噪声强度越大，阻尼越小，吸积盘振荡与湍流作用、内部噪声和粘滞作用更容易产生共振，低频准周期振荡现象越明显。

3 总结

本文考虑频率涨落噪声和指数形式关联随机力的作用，改进了吸积盘垂向随机振荡的朗之万方程模型，发现在振荡光度功率谱密度曲线上产生了广义的单峰和双峰随机共振，用来解释包含一个基频峰和二次谐频峰的低频准周期振荡。低频准周期振荡峰是由粘滞耗散、频率噪声和外界湍流随机力共同作用产生共振的结果，而且基频峰对应的中心频率为吸积盘振荡的特征频率ν₀，随机力关联时间决定了基频峰的高度和宽度，频率噪声强度和粘滞阻尼对二次谐频峰产生影响。

研究结果证实了吸积盘的随机振荡模型对于描述X射线光变是近似有效的，能较好地解释低频准周期振荡现象。在模型中，低频准周期振荡的中心频率仅由吸积盘的振荡频率决定。可以想象，如果通过调整吸积盘的几何结构，使振荡频率在瞬时获得足够大，将产生高频准周期振荡，在合适的动力学条件下，可以通过共振解释观测到的3∶2高频准周期振荡对，这需要我们下一步建立更加复杂的吸积盘随机振荡模型拟合观测到的高频准周期振荡。