耀变体(Blazar)是活动星系核(Active Galactic Nuclei, AGN)中比较特殊的一个子类，是最活跃的天体之一，其辐射主要来自喷流中的高能粒子，由于相对论聚束效应，在伽玛波段可以观测到辐射活动。耀变体可分为平谱射电耀变体(FSRQ)和蝎虎座BL型天体(BL Lac)两个子类^[1]，观测特征包括高光度、高偏振、快速光变以及非热辐射等^[2]；能谱分布(Spectra Energy Distribution, SED)呈现双峰结构，低能峰是由于同步辐射产生，高能峰是由于电子的逆康普顿辐射过程或者强子辐射产生^[3-4]。多普勒因子的估算方法有很多种，其中，通过射电光变估算多普勒因子是比较简单的一种，因而被广泛应用^[5-6]。光变研究能探索耀变体内部的辐射过程和物理机制，因此，研究天体的光变及其可能的周期性，对了解活动星系核具有重要意义。

耀变体CGRaBS J0835 + 6835(BWE 0830 + 6845)是一颗红移z=1.414的天体(α=08 35 47.608 086 678 5，δ=+68 35 11.493 386 746)^[7]。文[8]测得这个源在R波段的流量为0.000 13 mJy，文[9]在研究GHz峰频的样本源时，给出了此源在325 MHz，609 MHz和5 GHz等不同频率下的流量密度。本文利用欧文斯谷射电天文台的长时间监测数据(15 GHz)研究CGRaBS J0835 + 6835的射电光变，检验是否存在准周期振荡，还使用双指数函数拟合爆发过程，检验是否存在明显的聚束效应。目前，分析耀变体光变曲线准周期振荡的方法主要有两种：LSP方法^[10]和加权小波Z变换^[11]，这两种方法非常适合处理非等间隔数据。本文利用这两种方法对耀变体CGRaBS J0835 + 6835的15 GHz射电波段的光变周期进行研究，并对结果的可靠性和置信度进行评估。

1 光变曲线

如图 1是来自欧文斯谷射电天文台40 m望远镜观测的耀变体CGRaBS J0835 + 6835在15 GHz射电波段的光变曲线，去除少量的坏点后，采用460个数据点，约化儒略日MJD范围约为54 908~58 871。采用光变指数A^[12-13]判断天体的活跃程度：

$ A = \frac{{{f_{\max }} - {f_{\min }}}}{{{f_{\max }} + {f_{\min }}}}. $

(1)

其中，f_max和f_min分别为流量的最大值和最小值，光变指数A值越大表示天体越活跃。通过计算可以得到CGRaBS J0835 + 6835在射电波段的光变指数为0.98，这表明CGRaBS J0835 + 6835是个非常活跃的天体。

2 光变曲线的周期分析 2.1 LSP方法

检测隐藏在噪声中的周期信号是天文数据时间序列分析的一个重要目标。LSP方法可以对非均匀采样的时间序列周期图进行相位修正，能够在一定范围内对非均匀采样的时间间隔引起的误报周期进行修正^[10]。因此，LSP方法能很好地寻找隐藏在噪声中的准周期振荡光变。LSP方法的原理是基于傅里叶变换，使用正弦函数模型，基本公式为^[14-15]

$ Px(\omega ) = \frac{1}{2}\left\{ {\frac{{\sum\limits_{j = 1}^{{N_0}} {\left[ {x\left( {{t_j}} \right) - \bar x} \right]} \cos \omega \left( {{t_j} - \tau } \right)}}{{\sum\limits_{j = 1}^{{N_0}} {{{\cos }^2}} \omega \left( {{t_j} - \tau } \right)}} + \frac{{\sum\limits_{j = 1}^{{N_0}} {\left[ {x\left( {{t_j}} \right) - \bar x} \right]} \sin \omega \left( {{t_j} - \tau } \right)}}{{\sum\limits_{j = 1}^{{N_0}} {{{\sin }^2}} \omega \left( {{t_j} - \tau } \right)}}} \right\}, $

(2)

其中，τ为对应时间t的相位修正，计算公式为

$ \tan (2\omega \tau ) = \frac{{\sum\limits_{j = 1}^{{N_0}} {{{\sin }^2}} 2\omega {t_j}}}{{\sum\limits_{j = 1}^{{N_0}} {{{\cos }^2}} 2\omega {t_j}}}. $

(3)

LSP方法有几个优点是离散傅里叶变换所不具备的：首先，公式中包含时间项t，使得周期计算不随时间原点的变化而变化；其次，这种形式使得周期图分析完全等价于正弦曲线对数据的最小二乘拟合。为了验证LSP方法计算准周期的正确性，如图 2，先对周期图进行幂率拟合，得到蒙特卡罗方法计算置信度所需参数的幂率指数β^[16]，然后用Python编写程序计算CGRaBS J0835 + 6835的置信度，结果如图 3。图 3中绿色线表示准周期图，峰值就是准周期结果；蓝色、红色、黑色线分别表示蒙特卡罗模拟的95%、99%、99.7%置信度。图 3中的绿线有两个明显的峰值分别为562天和1 020天，并且两个峰的置信度都超过99.7%，说明两个峰值结果都是可靠的，但是1 020天约为562天的两倍，所以取最小正周期即562天作为CGRaBS J0835 + 6835的准周期结果。

2.2 加权小波Z变换方法

小波分析(Wavelet Aanlysis)是一种新的分析方法，是继傅里叶分析之后纯粹数学和应用数学完美结合的又一典范。利用Morlet小波作为小波母函数进行变换，复Morlet小波为^[17]

$ \psi (t) = {{\rm{e}}^{ - {t^{2/2}}}}\left( {{{\rm{e}}^{{\rm{i}}{w_0}t}} - {{\rm{e}}^{ - {w_0}/2}}} \right), $

(4)

文[18]定义了加权小波变换

$ WWT = \frac{{\left( {{N_{{\rm{eff }}}} - 1} \right){V_y}}}{{2{V_x}}}, $

(5)

其中，

$ {N_{{\rm{eff }}}} = \frac{{{{\left( {\sum {{\omega _\alpha }} } \right)}^2}}}{{\sum {\omega _\alpha ^2} }};{V_x} = \frac{{\sum\limits_\alpha {{\omega _\alpha }} {x^2}\left( {{t_\alpha }} \right)}}{{\sum\limits_\beta {{\omega _\beta }} }} - {\left[ {\frac{{\sum\limits_\alpha {{\omega _\alpha }} x\left( {{t_\alpha }} \right)}}{{\sum\limits_\beta {{\omega _\beta }} }}} \right]^2};{V_y} = \frac{{\sum\limits_\alpha {{\omega _\alpha }} {y^2}\left( {{t_\alpha }} \right)}}{{\sum\limits_\beta {{\omega _\beta }} }} - {\left[ {\frac{{\sum\limits_\alpha {{\omega _\alpha }} y\left( {{t_\alpha }} \right)}}{{\sum\limits_\beta {{\omega _\beta }} }}} \right]^2}. $

由于小波形状的变化，低频部分的有效数据N_eff 比高频部分的有效数据多，因此，加权小波变换的值向高频部分偏移，导致结果出现偏差，于是文[18]根据Z统计量定义了加权小波Z变换为^[10]

$ Z = \frac{{\left( {{N_{{\rm{eff}}}} - 3} \right){V_y}}}{{2\left( {{V_x} - {V_y}} \right)}}. $

(6)

利用(6)式可以得到如图 4的准周期图，(a)是时间周期平面中加权小波Z变换功率彩色缩放的分布，(b)是时间平均功率(黑色曲线)作为周期的函数图，其中蓝色虚线表示置信度在99.7%和95%的位置。根据加权小波Z变换图的极大值可以确定光变曲线的周期。图 4显示了一个明显的极大值560天且置信度超过99.7%，也就是说560天的准周期结果是可靠的。

3 多普勒因子分析

多普勒因子可以通过喷流中物质流速度V以及视角θ两个本征参数定义，然而这两个量都是难以直接观测的，因此，用间接方法估算多普勒因子。在各种估算多普勒因子的方法中，通过射电光变估算相对比较简单，因而被广泛应用。

文[19-20]详细介绍了亮温度和双指数函数公式及估算多普勒因子的过程。双指数函数公式为

$ F(t) = {F_{\rm{c}}} + {F_0}{\left( {{{\rm{e}}^{\frac{{{t_0} - t}}{{{t_{\rm{r}}}}}}} + {{\rm{e}}^{\frac{{t - {t_0}}}{{{t_{\rm{d}}}}}}}} \right)^{ - 1}}. $

(7)

利用亮温度公式估算CGRaBS J0835 + 6835的射电亮温度分布在3.7 × 10¹³~2.3 × 10¹⁴ K，超过逆康普顿极限温度或平衡亮温度几个量级, 表明此源的射电辐射存在明显的聚束效应。双指数函数拟合的6个爆发过程如图 5。类似文[21]选取平衡亮温度作为内禀亮温度，本文进一步估算对应的光变多普勒因子。每个爆发过程的拟合参数如表 1，第1列为约化儒略日范围；第2列为拟合所得的约化最小残差平方和${\hat \beta }$，详细定义见文[22]；第3列为归一化后的卡方值，定义为χ²/dof，${\chi ^2} = \sum\limits_{i = 1}^N {{{\left( {\frac{{{F_i} - {f_i}}}{\sigma }} \right)}^2}} $，其中F_i为实测流量，f_i为拟合流量，σ为残差绝对值的标准偏差，自由度dof=N-5；第4列为拟合峰值对应的约化儒略日及误差；第5列为基底流量及误差；第6列为爆发的幅度及误差；第7列和第8列分别为指数上升时标和下降时标及误差；第9列为每个爆发过程对应的多普勒因子，δ_V的值分布在8.87 ± 0.68~14.45 ± 0.23，平均值约为10.76。

4 讨论与结论

本文通过收集欧文斯谷射电天文台的长时间监测数据，利用LSP方法和加权小波Z变换方法对CGRaBS J0835 + 6835的射电波段光变曲线的准周期振荡进行分析。LSP方法分析得到的准周期约为562天和1 020天，并且置信度都大于3σ(99.7%)，1 020天大约是562天的两倍，故可以认为CGRaBS J0835 + 6835的准周期约为562天。加权小波Z变换方法得到的准周期约为560天，置信度也超过了3σ(99.7%)。两种检测准周期的方法得到的结果都大致约为560天，检验的结果相互印证，所以，我们认为CGRaBS J0835 + 6835的准周期约为560天。同时，本文利用双指数函数拟合了6个爆发过程，得到δ_V的值分布在8.87 ± 0.68~14.45 ± 0.23，平均值约为10.76。结果表明，CGRaBS J0835 + 6835的射电辐射存在明显的聚束效应，支持相对论性喷流模型。目前，耀变体长周期光变的物理模型仍然是悬而未决的问题，为了解决这些问题，天文学家提出了很多物理模型，例如双黑洞模型^[23]、螺旋喷流模型^[24]和薄吸积盘理论^[25]等。

根据CGRaBS J0835 + 6835 560天的准周期，其物理过程可能是喷流的螺旋运动^[26]。文[27]在给出的观测准周期P与实际物理驱动周期P_d之间的关系公式为

$ {P_{\rm{d}}} \approx \frac{{\gamma _{\rm{b}}^2}}{{1 + z}}P, $

(8)

其中，红移z=1.414^[7]；γ_b为体洛伦兹因子，约等于7.5^[28]。基于(8)式可以得到实际物理驱动周期P_d≈35.75年，另外，如果γ_b=15作为参数^[29]，P_d≈143.00年。当主黑洞和次黑洞的质量比R ≤ 1/3，称为主并合；如果3 ≤ R ≤ 10⁴，称为次并合^[30]。无论质量比如何，主黑洞的质量都可以用公式^[31]

$ M \approx P_{\rm{d}}^{\frac{8}{5}}{R^{\frac{3}{5}}}{10^6}{M_ \odot } $

(9)

估计，其中，P_d以年为单位，对于超大质量双黑洞系统的主并合，质量比可以假设为R=3/2，代入(9)式，得到CGRaBS J0835 + 6835主黑洞质量大约为M=3.7 × 10⁸M_⊙。同理，如果使用γ_b=15作为参数，次并合的超大质量双黑洞系统的主黑洞质量为M=3.5 × 10⁹M_⊙，故CGRaBS J0835 + 6835的主黑洞质量约为3.7 × 10⁸M_⊙~3.5 × 10⁹M_⊙。