Lomb-Scargle Periodogram方法研究耀變體3C 454.3長周期光變特性

陸 林，張皓晶，張 歡，馬凱旋

(云南師范大學(xué)物理與電子信息學(xué)院，云南 昆明 650092)

耀變體的多波段光變時標(biāo)是一個重要的物理參數(shù)，耀變體的長周期光變時標(biāo)可以幫助我們研究天體的中心黑洞質(zhì)量、內(nèi)部結(jié)構(gòu)和輻射區(qū)域等問題。目前研究此類問題的方法很多，有結(jié)構(gòu)函數(shù)法、離散相關(guān)函數(shù)法(Discrete Correlation Function, DCF)、period4方法、功率譜(Power Spectral Density, PSD)、Jurkevich方法、加權(quán)小波Z變換法(Weighted Wavelet Z-transform, WWZ)[1]和LSP方法[2]等。由于天體光學(xué)波段和紅外波段的觀測受天氣、設(shè)備等多種因素的影響，觀測數(shù)據(jù)在時間序列上不連續(xù),在光變周期的研究中引入較大的誤差,結(jié)果的精確度一直是該研究領(lǐng)域的熱點(diǎn)問題。文[3]對這個天體光學(xué)和紅外波段的數(shù)據(jù)采用功率譜分析方法，得到了與本文相同的結(jié)果，但是在功率譜方法中要求數(shù)據(jù)為等間隔，所以文[3]對相鄰的數(shù)據(jù)進(jìn)行了平均插值處理，這樣可能損害了原始數(shù)據(jù)的真實(shí)性，引入人為因素，還可能帶來其他疊加對周期的干擾。本文從SMARTS和文[3]中收集了耀變體3C 454.3[3]近10年光學(xué)波段的準(zhǔn)同時性數(shù)據(jù)，利用LSP方法對耀變體3C 454.3的非均勻數(shù)據(jù)進(jìn)行分析。LSP方法不需要對數(shù)據(jù)進(jìn)行插值，但是在功率譜結(jié)果中可能出現(xiàn)虛假的峰值。為了進(jìn)一步評估功率譜峰值的可信度，通過模擬大量光變曲線，并采用LSP方法計算模擬光變曲線的功率譜，對模擬光變曲線的功率譜結(jié)果進(jìn)行抽樣，估計相應(yīng)波段的置信度。通過置信度評估所出現(xiàn)峰值的可靠性，我們可以得出更為準(zhǔn)確的準(zhǔn)周期，并且在分析結(jié)果中得到的周期只有一個，而沒有其他的干擾周期，這種模擬大量光變曲線的方法稱為蒙特卡洛模擬方法。為了驗(yàn)證蒙特卡洛模擬方法的可行性和結(jié)果的準(zhǔn)確性，我們根據(jù)5個波段的數(shù)據(jù)，利用加權(quán)小波Z變換方法進(jìn)行長周期光變分析，也得到相同的結(jié)果，并通過文[4]中的標(biāo)準(zhǔn)差公式計算相應(yīng)波段的置信度，得出了可靠的準(zhǔn)周期,從而驗(yàn)證了蒙特卡洛模擬方法的可行性和結(jié)果的準(zhǔn)確性。本文首次將耀變體3C 454.3光學(xué)B，V，R波段和紅外J，K波段的離散數(shù)據(jù)應(yīng)用加權(quán)小波Z變換法和蒙特卡洛模擬方法來研究相應(yīng)波段的光變周期，為耀變體的長周期光變研究找到了一種新方法，通過加權(quán)小波Z變換法和蒙特卡洛模擬方法可以得到非平穩(wěn)信號的準(zhǔn)周期，大大提高了周期的計算精度，從而得到更加準(zhǔn)確的準(zhǔn)周期。我們使用這兩種方法得到5個波段都存在約454天的準(zhǔn)光變周期。從長周期的研究中我們可以得到3C 454.3中心黑洞質(zhì)量和輻射區(qū)半徑，從而為活動星系核(Active Galactic Nucleus, AGN)物理模型的研究提供重要參數(shù)。

1 觀測數(shù)據(jù)和研究

由SMARTS數(shù)據(jù)庫和文[3]中耀變體3C 454.3光學(xué)B，R，V波段和紅外J，K波段的數(shù)據(jù)得到的光變曲線如圖1，橫軸為儒略日(MJD)，縱軸為星等值(Mag)。

圖1 光學(xué)B，R，V波段和紅外J，K波段的光變曲線

根據(jù)流量與星等的轉(zhuǎn)化關(guān)系,5個波段的星等轉(zhuǎn)化為流量的公式為[5]

(1)

其中，F(xiàn)為每個波段對應(yīng)的流量值，單位為Jy；m為對應(yīng)波段的星等值。

1.1 光變曲線的長周期分析

從SMARTS 數(shù)據(jù)庫和文[3]獲得星等的數(shù)據(jù)之后，通過流量與星等的轉(zhuǎn)化公式得到各個波段的流量數(shù)據(jù)。接下來本文采用蒙特卡洛模擬[6]和加權(quán)小波Z變換法計算各個波段的周期和各自的置信度。即用圖1的光變曲線獲得流量數(shù)據(jù)后，利用LSP方法得到相應(yīng)的光變周期和蒙特卡洛模擬方法得到相應(yīng)的置信度估計。為了檢驗(yàn)蒙特卡洛模擬方法置信度評估的可行性和結(jié)果的準(zhǔn)確性，我們通過加權(quán)小波Z變換法計算相應(yīng)波段的光變周期，通過標(biāo)準(zhǔn)差公式，我們計算了相應(yīng)的加權(quán)小波Z變換周期的置信度。

1.2 LSP方法分析長周期光變

計算隱藏在噪聲中的周期信號是天文數(shù)據(jù)時間序列分析中的一個重要目標(biāo)。LSP方法可以對非均勻采樣的時間序列周期圖進(jìn)行相位修正處理，能夠在一定范圍內(nèi)對非均勻采樣的時間間隔引起的誤報周期進(jìn)行修正，并且LSP方法不要求數(shù)據(jù)是均勻的時間序列，不用對其進(jìn)行插值，從而減少了人為因素對數(shù)據(jù)的干擾。因此，LSP方法能尋找隱藏在噪聲中的準(zhǔn)周期振蕩光變。現(xiàn)在假定一組時間序列x(tj)，(j=1, 2, 3, …,N)，則時間序列的LSP功率譜為

(2)

(3)

圖1的光變曲線是一個非均勻的時間序列，要得到隱藏在其中的真實(shí)周期非常困難。我們利用LSP方法分析真實(shí)數(shù)據(jù)的功率譜結(jié)果，但是由于這些結(jié)果中出現(xiàn)虛假峰值，很難判別真實(shí)的周期頻率，所以為了解決這個問題，我們需要對這些功率譜的峰值進(jìn)行置信度評估。

1.3 蒙特卡洛模擬方法和置信度評估

活動星系核的功率譜一般呈紅噪聲冪律分布[6]，即P∝f-α，為了獲得功率譜指數(shù)α，對所得到的功率譜結(jié)果取對數(shù),并進(jìn)行一元線性回歸擬合，可以獲得功率譜指數(shù)α。光學(xué)B，R，V波段和紅外J，K波段的功率譜對數(shù)坐標(biāo)擬合結(jié)果如圖2。

圖2 光學(xué)B，R，V波段和紅外J，K波段的功率譜對數(shù)坐標(biāo)擬合

圖2中的橫坐標(biāo)是頻率f的對數(shù)，縱坐標(biāo)為功率譜的對數(shù)，紅色直線為相應(yīng)的一元線性擬合結(jié)果，圖中直線的斜率為功率譜指數(shù)α。通過一元線性回歸擬合，我們得到光學(xué)和紅外波段的各個冪律指數(shù)α?；诠β首V指數(shù)α和真實(shí)數(shù)據(jù)，我們模擬了5 000條光變曲線，對模擬的光變曲線利用LSP方法逐一進(jìn)行計算，獲得模擬光變曲線的功率譜結(jié)果，對這些模擬光變曲線的功率譜進(jìn)行置信度抽樣，得到各個波段不同的置信度曲線，最終將真實(shí)數(shù)據(jù)的功率譜和置信度曲線整合，即可得到準(zhǔn)確的周期頻率。

1.4 LSP方法的置信度評估結(jié)果

光學(xué)B，R，V波段和紅外J，K波段的置信度評估結(jié)果如圖3。圖3分別為各波段的蒙特卡洛模擬LSP方法結(jié)果，橫坐標(biāo)為頻率f的對數(shù)，縱坐標(biāo)為相應(yīng)波段功率譜的對數(shù)，紅色曲線代表真實(shí)數(shù)據(jù)的功率譜結(jié)果，黑色曲線代表置信度為99.7%，紫色曲線代表置信度為99%，橙色曲線代表置信度為95%。我們?nèi)〕^置信度99.7%的峰值為所對應(yīng)的周期頻率，如圖中黑色箭頭所指的位置。

圖3 光學(xué)B，R，V波段和紅外J，K波段的置信度評估結(jié)果

由圖3我們看到，在光學(xué)B，R，V波段和紅外J，K波段只存在一個很明顯的峰值超過99.7%，如圖3中箭頭所指的位置，相應(yīng)的峰值頻率保留4位小數(shù)都為f= 0.002 2(d-1)，所以光學(xué)和紅外波段都存在約454天的周期。通過這種方法我們得到了與文[3]相同的結(jié)果，并且這種方法不會引入其他干擾因素。為了檢驗(yàn)以上方法的可行性和結(jié)果的準(zhǔn)確性，我們通過加權(quán)小波Z變換法[1]計算相應(yīng)波段的周期。由于標(biāo)準(zhǔn)差和置信度的關(guān)系，我們計算了相應(yīng)的加權(quán)小波Z變換功率譜的置信度，本文采用文[4]給出的正態(tài)分布標(biāo)準(zhǔn)差公式。

2 加權(quán)小波Z變換法分析長周期光變

經(jīng)典的時間頻率分析一般采用傅里葉變換和小波分析方法，但是由于天文數(shù)據(jù)往往是非等間距的離散數(shù)據(jù)，在使用傅里葉變換和小波變換分析時，我們需要對數(shù)據(jù)進(jìn)行插值變成均勻的數(shù)據(jù)，這樣損害了數(shù)據(jù)的真實(shí)性，在使用傅里葉變換處理這些非等間隔數(shù)據(jù)時可能出現(xiàn)偽周期，所以在小波變換的基礎(chǔ)上文[1]提出了加權(quán)小波Z變換法，不僅可以更加有效地處理非等間隔的時間序列的周期，還可以一定程度上反映周期的穩(wěn)定性。加權(quán)小波Z變換是將時間序列投影到3個正交歸一的基函數(shù)上，即φ1(ti)=1,φ2(ti)=cos[ω(t-τ0)]和φ3(ti)=sin[ω(t-τ0)]，在投影上做統(tǒng)計加權(quán)ωi=exp[-cω2(ti-τ0)]，將不均勻的數(shù)據(jù)通過權(quán)重調(diào)節(jié)，避免周期分析時受到數(shù)據(jù)過密的影響，其中使用的母函數(shù)為Morlet小波[1]。加權(quán)小波Z變換的定義為

(4)

這個等式的分子分母分別滿足自由度為Neff-3和2的F分布。其中，Neff為有效數(shù)據(jù)點(diǎn)個數(shù)；Vx和Vy分別為觀測數(shù)據(jù)和模擬函數(shù)的加權(quán)變量，即

(5)

(6)

(7)

其中，c為衰減因子；τ0為時移；ω為尺度因子；ωλ(λ=1, 2, 3, …,n)為對應(yīng)的測試頻率。

圖4 光學(xué)B，R，V波段和紅外J，K波段的加權(quán)小波Z變換法分析結(jié)果和相應(yīng)的周期頻率

3 3C 454.3中心黑洞質(zhì)量和輻射區(qū)域半徑的估計

3.1 中心黑洞質(zhì)量

耀變體中心有一個超大質(zhì)量黑洞，很多輻射現(xiàn)象是由中心黑洞引起。通過長時標(biāo)光變能夠獲得中心黑洞的質(zhì)量，假設(shè)由薄吸積盤理論引起的長周期光變的公式為[7]

(8)

其中，β為黏度參數(shù)；M6=M/(106M⊙)為中心黑洞質(zhì)量約化單位。吸積率M·和廣義應(yīng)力張量參數(shù)μ，當(dāng)μ=0.5時，M·=ME/ε，ME為愛丁頓吸積率，ε為吸積效率。爆發(fā)時間間隔主要取決于μ。 文[7]建議μ=0.5時，磁場的逃逸速率更低。此時,熱極限周期時間為tcyc=2tburts，

(9)

其中，β0.1=β/0.1，所以當(dāng)β=0.1，μ=0.5，M·=0.2M·C時，我們采用5個波段1.244年的周期得到中心黑洞的質(zhì)量約為M≈0.24×106M⊙。由于利用薄吸積盤理論的分析方法得到這個結(jié)果，沒有考慮黑洞自旋的影響，所以比一般的平譜射電類星體的中心黑洞質(zhì)量偏小。

3.2 輻射區(qū)半徑

對于耀變體的長周期光變現(xiàn)象，目前還沒有最佳的物理模型解釋，現(xiàn)有的理論為雙黑洞模型[8]、薄盤的熱不穩(wěn)定性模型[8]和螺旋噴流等模型[9-10]。如果是薄盤不穩(wěn)定性造成的，從搜集的數(shù)據(jù)中獲得各個波段的星等變化情況如表1。

由表1可以看出，耀變體3C 454.3在光學(xué)和紅外波段都具有較為劇烈的變化，并且紅外波段的變化幅度比光學(xué)波段更加劇烈。采用薄盤不穩(wěn)定性分析長周期T[8]，通過長周期T可以確定熱不穩(wěn)定性產(chǎn)生的區(qū)域，公式為

表1 每個波段的星等變化

(10)

其中，x=R/Rg，Rg為施瓦西半徑；β為黏滯系數(shù)；M8=M/(108M⊙)。這里我們?nèi)≈ 1.244 year，中心天體黑洞質(zhì)量為M≈0.24×106M⊙，所以我們得到輻射區(qū)域半徑為R≈3.68×108km。

4 結(jié) 果

對耀變體3C 454.3光學(xué)波段和紅外波段的研究中，本文采用LSP方法得到了各個波段的周期，評估了峰值頻率的置信度，得到更為可靠的結(jié)果，并采用加權(quán)小波Z變換法驗(yàn)證蒙特卡洛模擬評估置信度方法的結(jié)果，我們發(fā)現(xiàn)：

(1)兩種方法都在光學(xué)B，R，V波段和紅外J，K波段發(fā)現(xiàn)454天，約1.244年的光變周期，通過蒙特卡洛模擬評估置信度方法，我們清楚地看到對于5個波段都只有一個峰值頻率超過99.7%，并且沒有引入其他干擾結(jié)果。為了進(jìn)一步驗(yàn)證這種方法的可行性，我們對5個波段采用加權(quán)小波Z變換的分析方法，也得到了與蒙特卡洛模擬方法一致的結(jié)果，都為454天，約1.244年。

(2)對于中心的超大質(zhì)量黑洞，通過文獻(xiàn)給出的長周期與黑洞質(zhì)量的公式，結(jié)合本文中5個波段共同的長周期數(shù)據(jù)，我們得到中心黑洞的質(zhì)量為M≈0.24×106M⊙。由于這個結(jié)果是利用薄吸積盤理論的分析方法，沒有考慮黑洞自旋的影響，所以比一般平譜射電類星體的中心黑洞質(zhì)量偏小，輻射區(qū)半徑為R≈3.68×108km，計算得到的輻射區(qū)半徑比文[11]給出的輻射區(qū)半徑在相應(yīng)的施瓦西半徑下要小，由于文[11]采用的黑洞質(zhì)量比本文中所用的黑洞質(zhì)量大，所以此處我們得到的輻射區(qū)半徑是合理的。

(3)通過以上的計算，我們預(yù)測耀變體3C 454.3在2021年6月左右將再次爆發(fā)。我們將用云南天文臺麗江觀測站2.4 m光學(xué)望遠(yuǎn)鏡做進(jìn)一步實(shí)測驗(yàn)證，為蒙特卡洛模擬方法研究耀變體的長周期光變尋求觀測證據(jù)。

致謝：真誠致謝SMARTS團(tuán)隊(duì)和吳月承等人在文[3]中提供的光學(xué)和紅外波段的星等數(shù)據(jù)。