• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    Blazar天體3C 66A光學(xué)波段準(zhǔn)周期光變分析

    2022-04-02 08:33:30孫中諾董富通
    天文學(xué)報(bào) 2022年2期
    關(guān)鍵詞:光變時(shí)標(biāo)噴流

    程 巖 劉 風(fēng) 孫中諾 董富通

    (1 德州機(jī)電工程學(xué)校信息工程系 德州 251200)

    (2 德州科技職業(yè)學(xué)院信息工程系 德州 251200)

    1 引言

    Blazar是活動(dòng)星系核(AGN)的一個(gè)子類, 它是一種活動(dòng)非常劇烈的天體. 光變是Blazar天體的一個(gè)重要特征, Blazar天體在全波段都會(huì)表現(xiàn)出劇烈的光變, 光變時(shí)標(biāo)從幾分鐘到幾十年不等[1–2]. 尋找和研究這類天體的光變特性對(duì)于認(rèn)識(shí)它們的內(nèi)部結(jié)構(gòu)有重要的意義. 比如根據(jù)Blazar天體的天內(nèi)光變(光變時(shí)標(biāo)小于1 d), 可以推測(cè)其輻射區(qū)的尺度以及中心黑洞的質(zhì)量等[3–4].

    大量研究發(fā)現(xiàn), Blazar天體的光變?cè)诙鄠€(gè)波段存在周期/準(zhǔn)周期的變化. 在光學(xué)波段, Fan等[5]整理了10個(gè)射電選BL Lac天體的光變數(shù)據(jù), 并在其中9個(gè)天體中找到了周期為1.4–17.85 yr的周期光變; Graham等[6]研究了PG 1302–102的光學(xué)光變數(shù)據(jù), 并發(fā)現(xiàn)了(1884±88) d的光變周期; Rani等[7]在Blazar天體S5 0716+714的光學(xué)R波段發(fā)現(xiàn)了時(shí)標(biāo)為~15 min的準(zhǔn)周期光變. 在X射線/γ射線波段人們也發(fā)現(xiàn)了大量的準(zhǔn)周期光變的證據(jù). 比如:Espaillat等[8]在3C 273的X射線光變數(shù)據(jù)中發(fā)現(xiàn)了時(shí)標(biāo)為~55 min的準(zhǔn)周期光變; Rani等[9]分別在AO 0235+164和1ES 2321+419的X射線數(shù)據(jù)中發(fā)現(xiàn)了時(shí)標(biāo)為17 d和420 d的光變周期; Ackermann等[10]在PG 1553+113的γ射線波段發(fā)現(xiàn)了2.2 yr的光變周期; Zhou等[11]在PKS 2247–131的γ射線光變數(shù)據(jù)中找到了~34.5 d的準(zhǔn)周期光變. 而在射電波段也有大量的Blazar天體被發(fā)現(xiàn)具有準(zhǔn)周期光變行為, 比如: AO 0235+164[12–13]、PKS 1510–089[14]、NRAO 530[15]和PKS 1156+295[16]等; 另外, Ciaramella等[17]和Fan等[18]分別研究了兩個(gè)樣本的Blazar天體多個(gè)波段的射電數(shù)據(jù)并發(fā)現(xiàn)了周期性光變.

    Blazar天體3C 66A從射電波段到γ射線波段都表現(xiàn)出明顯的光變. 在X射線波段, Ghosh等[19]發(fā)現(xiàn)了時(shí)標(biāo)為幾小時(shí)的天內(nèi)光變. 這種光變也存在于紅外波段和射電波段[20].γ射線波段存在時(shí)標(biāo)為1.5 d的光變[21]. 在光學(xué)波段, 人們發(fā)現(xiàn)了很多周期性光變的證據(jù). Lainela等[22]發(fā)現(xiàn)了65 d的光變周期, 但是這一周期在后續(xù)的觀測(cè)和研究中并沒(méi)有得到證實(shí). 而最新的觀測(cè)數(shù)據(jù)表明存在其他的光變周期, 比如: Fan等[23]使用13 yr的光學(xué)V波段觀測(cè)數(shù)據(jù)找到了~156.0 d和~1.96 yr的光變周期. Belokon等[24]和Kaur等[25]分別分析了3C 66A不同時(shí)間段的光變數(shù)據(jù), 都發(fā)現(xiàn)了時(shí)標(biāo)為2.5 yr的準(zhǔn)周期光變.Otero-Santos等[26]分析了3C 66A光學(xué)R波段和V波段的數(shù)據(jù), 找到了時(shí)標(biāo)為3 yr的光變周期. 為了進(jìn)一步確定3C 66A的光變周期,更好地認(rèn)識(shí)其內(nèi)部物理性質(zhì), 本文收集了2003—2021年間約18 yr的光學(xué)觀測(cè)數(shù)據(jù), 并使用了不同的方法對(duì)其進(jìn)行了分析.

    2 觀測(cè)數(shù)據(jù)

    3C 66A在光學(xué)波段有很長(zhǎng)時(shí)間的觀測(cè)歷史, 我們收集了光學(xué)V波段的星等數(shù)據(jù)進(jìn)行光變分析. 本文數(shù)據(jù)主要有3個(gè)來(lái)源, 第1個(gè)是Fan等[23]使 用 上 海 天 文 臺(tái)(ShAO)1.5 m望 遠(yuǎn)鏡觀測(cè)的數(shù)據(jù)(時(shí)間范圍2003—2009年), 第2個(gè)是通過(guò)AAVSO (The American Association of Variable Star Observers)數(shù)據(jù)庫(kù)獲得(2003—2021年)1https://www.AAVSO.org的數(shù)據(jù), 最后一個(gè)是來(lái)自于Steward天文臺(tái)的觀測(cè)數(shù)據(jù)(2009—2018年)[27]. 圖1給出了3種數(shù)據(jù)綜合之后的光變曲線, 其中橫坐標(biāo)為儒略日(Julian Date, J.D.). 圖中共2304個(gè)數(shù)據(jù)點(diǎn),總的時(shí)間跨度約為18 yr, 平均時(shí)間間隔為2.7 d, 最大時(shí)間間隔為237 d. 之前的研究發(fā)現(xiàn)3C 66A的典型的光變周期在1.96–3 yr之間, 所以這些數(shù)據(jù)的時(shí)間間隔對(duì)于我們尋找年量級(jí)的光變周期不會(huì)有很大影響.

    圖1 3C 66A光學(xué)V波段的歷史光變曲線, 圖中不同來(lái)源的數(shù)據(jù)使用了不同的符號(hào)及顏色來(lái)表示.Fig.1 The historical light curve of 3C 66A in the optical V band. Different colours and symbols denote data from different observations.

    3 光變周期分析

    3.1 Jurkevich方法

    為了研究3C 66A的準(zhǔn)周期光變行為, 我們首先使用了Jurkevich方法對(duì)其光變曲線進(jìn)行了分析來(lái)尋找光變周期. Jurkevich方法屬于“相位折疊方法”, 非常適用于處理非等間隔數(shù)據(jù), 并已經(jīng)被廣泛用于Blazar天體的周期分析中[28–29]. Jurkevich方法的原理是根據(jù)不同的測(cè)試周期, 按照相位不同將數(shù)據(jù)劃分到m個(gè)分組中, 計(jì)算每一組的離均差平方和, 進(jìn)而得到所有組的總離均差平方和V2m[30].V2m的極小值對(duì)應(yīng)的試驗(yàn)周期就接近樣本的真實(shí)周期. 但是V2m的值還會(huì)受到很多因素的影響, 比如分組數(shù)、數(shù)據(jù)樣本本身的采樣方式等,會(huì)導(dǎo)致V2m出現(xiàn)很多虛假極小值. 因此, Kidger等[31]提出了周期真實(shí)性的F判據(jù):

    一般情況下, 如果F≥0.5, 表明數(shù)據(jù)中有較強(qiáng)的周期, 而如果F≤0.25, 表明數(shù)據(jù)的周期比較弱.

    圖2給出了我們對(duì)3C 66A光學(xué)V波段星等數(shù)據(jù)的Jurkevich分析結(jié)果. 圖中有兩個(gè)比較明顯的極小值, 其中第1個(gè)極小值為V2m= 0.776, 對(duì)應(yīng)的F=0.29. 第2個(gè)極小值為V2m= 0.72, 對(duì)應(yīng)的F= 0.39.根據(jù)F判據(jù),F值的兩個(gè)極小值都大于0.25, 說(shuō)明它們所對(duì)應(yīng)的周期很可能是數(shù)據(jù)的真實(shí)周期. 這兩個(gè)周期值分別為(850±90)d和(1150±140)d, 其中周期誤差是由V2m曲線中這兩個(gè)極小值對(duì)應(yīng)的半峰全寬所得到的. 我們所得到的(850±90)d周期與Fan等[23]在V波段發(fā)現(xiàn)的(1.96±0.63)yr周期在誤差范圍內(nèi)是一致的, 但是我們的結(jié)果中并沒(méi)有出現(xiàn)他們所發(fā)現(xiàn)的~156.0 d光變周期. Otero-Santos等[26]分析了光學(xué)R波段的數(shù)據(jù)得到了2.28 yr和3.14 yr的光變周期, 而我們的結(jié)果與此非常一致.

    圖2 3C 66A光學(xué)V波段數(shù)據(jù)的Jurkevich方法分析結(jié)果Fig.2 The Jurkevich test result for the period search in the optical V band of blazar 3C 66A

    3.2 Lomb-Scargle方法

    Lomb-Scargle (LS)方法是分析天體準(zhǔn)周期光變的一種非常有效的方法, 它使用正弦波來(lái)模擬光變信號(hào), 并利用最小二乘法對(duì)光變數(shù)據(jù)進(jìn)行擬合,因此該方法非常適用于處理非等間隔數(shù)據(jù)[32–33].圖3給出了我們用LS方法計(jì)算的3C 66A光學(xué)V波段的能譜, 其中有兩個(gè)明顯的峰值, 對(duì)應(yīng)的周期值分別為(869±70)d和(1111±90)d. 其中的誤差是由LS能譜減去蒙特卡羅模擬的平均能譜后峰值的半峰全寬(FWHM)獲得(關(guān)于這一過(guò)程, 我們?cè)谙旅鏁?huì)做詳細(xì)說(shuō)明). 這一結(jié)果與我們用Jurkevich方法所得的結(jié)果是相同的, 另外在LS能譜中也沒(méi)有發(fā)現(xiàn)明顯的~156 d光變周期.

    圖3 3C 66A光學(xué)V波段數(shù)據(jù)的LS能譜, 圖中還給出了1000次蒙特卡羅模擬得到的95%和99%的置信水平.Fig.3 The LS power spectra of the optical V band of 3C 66A, along with that of the 95% and 99% confidence level generated by 1000 Monte Carlo simulations.

    Blazar天體的光變數(shù)據(jù)中通常包含一種隨機(jī)過(guò)程, 該過(guò)程的傅立葉變換能譜一般表現(xiàn)出冪律譜的形式:p(f)∝f-α(α≥1), 其中p(f)表示隨機(jī)過(guò)程的傅立葉變換能譜,f表示頻率. 這種隨機(jī)過(guò)程稱為“紅噪聲”. “紅噪聲”對(duì)Blazar天體的光變會(huì)產(chǎn)生重要的影響, 主要表現(xiàn)為在較長(zhǎng)的光變時(shí)標(biāo)下會(huì)產(chǎn)生較大振幅的隨機(jī)光變[34]. 這種隨機(jī)光變通過(guò)一般的周期分析方法很容易被誤認(rèn)為周期信號(hào). 因此為了評(píng)估所得周期的真實(shí)性, 我們需要在計(jì)算中考慮這種“紅噪聲”過(guò)程.

    我們使用了能譜響應(yīng)方法(PSRESP)[35]來(lái)評(píng)估“紅噪聲”對(duì)3C 66A光變數(shù)據(jù)的影響. 該方法已經(jīng)被廣泛用于活動(dòng)星系核的能譜分析中了(比如:Chatterjee等[36]、Bhatta等[37]). 它使用一個(gè)簡(jiǎn)單的冪律譜來(lái)模擬“紅噪聲”的傅立葉變換能譜, 通過(guò)蒙特卡羅方法模擬出光變曲線, 與真實(shí)的光變曲線進(jìn)行統(tǒng)計(jì)量比較, 得到“紅噪聲”的能譜, 進(jìn)而可以評(píng)估“紅噪聲”對(duì)光變數(shù)據(jù)的影響.

    本文使用冪律譜p(f)∝f-α來(lái)模擬3C 66A光變曲線中的“紅噪聲”,α的取值范圍為1–2.5共20個(gè)點(diǎn). 首先將光變曲線進(jìn)行傅立葉變換, 得到真實(shí)的傅立葉變換能譜.之后使用Timmer等[38]的方法,將不同α值的冪律譜進(jìn)行1000次的蒙特卡羅模擬, 并計(jì)算它們能譜的平均值. 通過(guò)比較真實(shí)能譜與這些平均能譜, 得到了斜率α的連續(xù)分布函數(shù)(如圖4所示), 其中最大值對(duì)應(yīng)的就是最佳擬合的α值. 我們得到的最佳擬合斜率α= 1.56±0.25, 其誤差是由圖4中斜率α的連續(xù)分布函數(shù)的半峰全寬得到的. 由此我們得到了最佳擬合能譜模型如圖5所示.

    圖4 3C 66A能譜斜率α的連續(xù)分布函數(shù)Fig.4 The success function of the power spectra slopes α of 3C 66A

    圖5 3C 66A光學(xué)V波段數(shù)據(jù)的能量譜及其最佳擬合模型Fig.5 The power spectra of 3C 66A in the V band, along with the best fitted model

    最后我們使用上面得到的最佳擬合模型對(duì)光變曲線進(jìn)行了1000次模擬, 對(duì)每一次模擬的光變曲線進(jìn)行重新采樣, 使其與真實(shí)的光變曲線的采樣率相同. 使用這些模擬的光變數(shù)據(jù), 我們建立了關(guān)于圖3中LS能譜峰值的置信水平. 其中(869±70)d的周期對(duì)應(yīng)的峰值具有>99%的置信水平, 而周期值為(1111±90)d的峰值置信水平為>95%, 而圖中其他峰值的置信水平都低于95%, 說(shuō)明并沒(méi)有其他在統(tǒng)計(jì)上比較顯著的周期存在.

    4 討論

    我們收集了來(lái)自于ShAO、AAVSO和Steward天文臺(tái)3個(gè)來(lái)源的Blazar天體3C 66A的光學(xué)V波段18 yr的光變數(shù)據(jù), 并使用Jurkevich方法和Lomb-Scargle方法兩種方法研究了其準(zhǔn)周期光變行為, 找到了~850 d (2.3 yr)和1150 d (3.2 yr)的光變周期.這一結(jié)果證實(shí)了之前對(duì)3C 66A準(zhǔn)周期光變的研究.Belokon等[24]研究了3C 66A在1972—1992年間光學(xué)B波段的光變數(shù)據(jù), 并發(fā)現(xiàn)了~2.5 yr的光變周期; Kaur等[25]分析了2005—2016年間的R波段數(shù)據(jù)同樣找到了~2.5 yr的光變周期;而我們使用2003—2021年間的V波段數(shù)據(jù)得到了~2.3 yr的光變周期,另外Otero-Santos等[26]使用了2005—2020年間R波段的光變數(shù)據(jù)同樣得到了~2.3 yr的光變周期, 這些結(jié)果足以說(shuō)明~2.3 yr的光變周期在3C 66A的歷史光變數(shù)據(jù)中是持續(xù)存在的. 而對(duì)于3.2 yr的光變周期, 在較早的觀測(cè)以及研究中并沒(méi)有發(fā)現(xiàn)該周期, 而且Otero-Santos等[26]研究發(fā)現(xiàn)該周期只存在于3C 66A的高活躍狀態(tài)(輻射流量較高的狀態(tài)), 說(shuō)明該周期并非穩(wěn)定周期.

    為了驗(yàn)證這一結(jié)論, 我們分別使用周期為2.3 yr和3.2 yr的正弦函數(shù)擬合了3C 66A的光學(xué)V波段光變曲線, 其結(jié)果如圖6所示. 為了更清楚地表示其長(zhǎng)周期光變趨勢(shì), 我們將光變曲線做了50 d的平滑處理. 首先在圖中可以看出3C 66A的長(zhǎng)周期光變并不是嚴(yán)格的正弦形式, 而且在整個(gè)時(shí)間范圍內(nèi)經(jīng)歷了高活躍狀態(tài)(J.D. 2456000之前, 平均星等值為14.38)和低活躍狀態(tài)(J.D. 2456000之后, 平均星等值為14.9)兩種狀態(tài). 因此我們使用單一的正弦函數(shù)并不能很好地?cái)M合其光變曲線, 但是通過(guò)比較光變曲線和正弦函數(shù)的變化趨勢(shì), 我們可以判定前面所獲得周期的可靠程度. 其中周期為3.2 yr的正弦函數(shù)在J.D. 2456290之前與光變曲線的變化趨勢(shì)基本相同, 而在J.D. 2456290之后其與光變曲線出現(xiàn)了比較大的偏差. 這驗(yàn)證了我們之前的結(jié)論,即3.2 yr的周期并不是穩(wěn)定周期, 并且只出現(xiàn)于3C 66A 的高活躍狀態(tài). 而周期為2.3 yr的正弦函數(shù)雖然其振幅與真實(shí)數(shù)據(jù)相差較大, 但是其總體變化趨勢(shì)與光變曲線是比較一致的. 這同樣也驗(yàn)證了我們之前的結(jié)論, 即2.3 yr的光變周期是一個(gè)穩(wěn)定周期. 雖然3.2 yr的光變周期并不穩(wěn)定, 但是它對(duì)3C 66A的長(zhǎng)周期光變有著重要影響. 實(shí)際上, 如果我們使用3.2 yr和2.3 yr兩個(gè)周期同時(shí)擬合3C 66A的光學(xué)V波段光變曲線, 可以得到較好的擬合效果.其結(jié)果如圖6所示, 雖然在某些部分?jǐn)M合曲線與光變曲線有較大的偏差(可能受到微光變或其他光變機(jī)制的影響), 但是總體上擬合曲線與光變曲線的變化趨勢(shì)基本一致.

    圖6 3C 66A光學(xué)V波段歷史光變曲線的正弦擬合結(jié)果, 黑線表示3C 66A光學(xué)V波段的光變曲線, 紅色和綠色的虛線分別表示周期為2.3 yr和3.2 yr的正弦曲線, 而藍(lán)色虛線表示這兩個(gè)周期的擬合結(jié)果.Fig.6 Sinusoidal fit to the historical light curve of 3C 66A in optical V band. The black line represents the light curve, the red and green dashed lines are sinusoidal functions with periods 2.3 yr and 3.2 yr respectively, and the blue dashed line is the result of fitting with the two periodicities.

    人們提出了各種模型來(lái)解釋Blazar天體的周期性光變, 其中比較成功的模型是超大質(zhì)量雙黑洞模型. Sillanp¨a¨a等[39]提出了該模型, 并解釋了OJ 287光學(xué)波段~12 yr的周期光變. 雙黑洞軌道運(yùn)動(dòng)可能導(dǎo)致周期性的吸積, 吸積盤(pán)的進(jìn)動(dòng)以及噴流的進(jìn)動(dòng)和噴流的扭曲, 進(jìn)而產(chǎn)生周期性光變[40–42]. 如果3C 66A的2.3 yr光變周期是由雙黑洞軌道運(yùn)動(dòng)造成的, 那么我們可以根據(jù)開(kāi)普勒運(yùn)動(dòng)公式估計(jì)其軌道半徑為:

    其中Tobs為觀測(cè)周期,G為萬(wàn)有引力常數(shù).3C 66A的中心黑洞質(zhì)量近似為M ~8.57×108M⊙[25], 其紅移值為z= 0.444[43]. 根據(jù)不同的觀測(cè), 其多普勒因子δ有不同的取值, 其變化范圍在2.6–15之間[23]. 如果我們?nèi)∑淦骄? 則可以得到雙黑洞軌道半徑R ~8.3×1016cm (~0.03 pc), 說(shuō)明3C 66A中的雙黑洞已經(jīng)演化到了“最后1 pc”階段(0.01–1 pc)[44].

    除此之外, 螺旋結(jié)構(gòu)的噴流模型也會(huì)導(dǎo)致Blazar天體的周期性光變[45]. 由強(qiáng)磁場(chǎng)或軌道運(yùn)動(dòng)導(dǎo)致的扭曲噴流, 在視角改變時(shí), 由于多普勒效應(yīng)會(huì)使得觀測(cè)的流量產(chǎn)生變化. 對(duì)于非均勻旋轉(zhuǎn)螺旋噴流, 其中發(fā)射物質(zhì)速度的方向與我們的視線夾角會(huì)周期性的變化. 由于多普勒增亮效應(yīng), 就會(huì)導(dǎo)致觀測(cè)的流量周期性變化[46–47]. 3C 66A是Blazar天體, 其輻射流量主要來(lái)自于噴流中的非熱輻射, 而甚長(zhǎng)基線陣列(VLBA)觀測(cè)發(fā)現(xiàn)其噴流呈現(xiàn)扭曲的形狀[48], 所以其周期性光變也可能是螺旋噴流旋轉(zhuǎn)造成的. 另外, 3C 66A同時(shí)也是BL Lac天體, 其中可能存在一個(gè)輻射效率很低的吸積盤(pán)[49], 而吸積盤(pán)的不穩(wěn)定傳遞到噴流當(dāng)中也會(huì)造成觀測(cè)流量的周期變化[50].

    致謝 我們感謝審稿人為本文提出的寶貴意見(jiàn), 本文使用了AAVSO數(shù)據(jù)庫(kù)和Steward天文臺(tái)提供的開(kāi)源數(shù)據(jù).

    猜你喜歡
    光變時(shí)標(biāo)噴流
    BL Lacs 天體1823+568 射電波段光變周期分析
    二階非線性中立型時(shí)標(biāo)動(dòng)態(tài)方程趨向于零的非振動(dòng)解的存在性
    時(shí)標(biāo)上具非正中立項(xiàng)的二階動(dòng)力方程的動(dòng)力學(xué)性質(zhì)
    “慧眼”發(fā)現(xiàn)迄今距離黑洞最近的高速噴流
    耀變天體3C 454.3 高能光變行為的研究
    噴流干擾氣動(dòng)熱數(shù)值模擬的若干影響因素
    100元人民幣上的“100”為什么會(huì)變顏色?
    耀變體噴流高能電子譜的形成機(jī)制
    發(fā)生在活動(dòng)區(qū)11931附近的重復(fù)噴流?
    基于奇異擾動(dòng)的永磁風(fēng)力發(fā)電機(jī)組雙時(shí)標(biāo)控制
    19禁男女啪啪无遮挡网站| 久久久国产欧美日韩av| 欧美午夜高清在线| 男人的好看免费观看在线视频| 美女高潮的动态| 一级毛片精品| 午夜福利高清视频| 国产成人福利小说| 国产成人一区二区三区免费视频网站| 黑人欧美特级aaaaaa片| 色老头精品视频在线观看| 嫩草影视91久久| 亚洲熟妇中文字幕五十中出| 高清毛片免费观看视频网站| 国产人伦9x9x在线观看| 精品一区二区三区视频在线观看免费| 国产精品久久久人人做人人爽| 亚洲欧美一区二区三区黑人| xxxwww97欧美| 天天一区二区日本电影三级| 日韩成人在线观看一区二区三区| 狠狠狠狠99中文字幕| 欧美黄色淫秽网站| 国产毛片a区久久久久| 午夜福利18| 久久久久久大精品| 最新在线观看一区二区三区| 91在线观看av| 国产真人三级小视频在线观看| 在线观看美女被高潮喷水网站 | 一二三四社区在线视频社区8| 国模一区二区三区四区视频 | 亚洲精品美女久久av网站| 久久久久九九精品影院| 午夜精品一区二区三区免费看| 日韩中文字幕欧美一区二区| 亚洲精品美女久久av网站| 欧美色视频一区免费| 久久国产精品人妻蜜桃| 亚洲精品在线美女| 国产高清视频在线观看网站| 午夜激情欧美在线| 一本综合久久免费| 久久精品aⅴ一区二区三区四区| 亚洲国产精品成人综合色| 91麻豆av在线| 一级黄色大片毛片| 五月玫瑰六月丁香| 热99re8久久精品国产| 看片在线看免费视频| 免费无遮挡裸体视频| 窝窝影院91人妻| 免费搜索国产男女视频| 夜夜爽天天搞| 亚洲人与动物交配视频| 久久伊人香网站| 日韩大尺度精品在线看网址| 国产伦在线观看视频一区| 级片在线观看| 国产午夜精品论理片| 99国产精品99久久久久| 白带黄色成豆腐渣| 日本五十路高清| 五月玫瑰六月丁香| 一个人免费在线观看电影 | 国产成人精品久久二区二区免费| 欧美乱码精品一区二区三区| 又黄又粗又硬又大视频| 国产精品99久久99久久久不卡| 动漫黄色视频在线观看| 国产精品久久久av美女十八| 国产成人福利小说| 毛片女人毛片| 亚洲男人的天堂狠狠| 9191精品国产免费久久| 午夜福利在线观看免费完整高清在 | 99精品久久久久人妻精品| 成人鲁丝片一二三区免费| 国产午夜精品论理片| 女同久久另类99精品国产91| 国产av在哪里看| 岛国视频午夜一区免费看| 午夜亚洲福利在线播放| 成人性生交大片免费视频hd| 久久久国产成人精品二区| 桃色一区二区三区在线观看| 午夜日韩欧美国产| 免费看光身美女| 99久久久亚洲精品蜜臀av| 国内揄拍国产精品人妻在线| 变态另类丝袜制服| 欧美在线一区亚洲| 午夜福利视频1000在线观看| 1024手机看黄色片| 99热这里只有是精品50| 麻豆一二三区av精品| 麻豆成人av在线观看| av天堂在线播放| 国产91精品成人一区二区三区| 欧美性猛交╳xxx乱大交人| 麻豆国产av国片精品| 法律面前人人平等表现在哪些方面| 亚洲在线观看片| 99热精品在线国产| 黄色丝袜av网址大全| 国产精品野战在线观看| 黑人欧美特级aaaaaa片| 黄色丝袜av网址大全| 搞女人的毛片| 又粗又爽又猛毛片免费看| 欧美黄色片欧美黄色片| 国产高清有码在线观看视频| 欧美乱妇无乱码| 亚洲欧洲精品一区二区精品久久久| 老司机在亚洲福利影院| 一本综合久久免费| 日本黄色视频三级网站网址| 国产三级在线视频| 国产午夜精品久久久久久| 哪里可以看免费的av片| 老汉色av国产亚洲站长工具| 一夜夜www| 亚洲人成伊人成综合网2020| 精品久久久久久久久久免费视频| 亚洲色图av天堂| 亚洲成人免费电影在线观看| 婷婷六月久久综合丁香| 一个人免费在线观看的高清视频| 国产精品久久视频播放| 成年人黄色毛片网站| 又爽又黄无遮挡网站| 国产成人av激情在线播放| 不卡av一区二区三区| 日韩中文字幕欧美一区二区| 国产乱人视频| 在线永久观看黄色视频| 听说在线观看完整版免费高清| 中亚洲国语对白在线视频| 成人亚洲精品av一区二区| 真人做人爱边吃奶动态| 欧美日韩瑟瑟在线播放| 亚洲中文字幕一区二区三区有码在线看 | 99国产极品粉嫩在线观看| 国产av不卡久久| 日本 欧美在线| 亚洲成a人片在线一区二区| 99热这里只有精品一区 | 人人妻人人澡欧美一区二区| 欧美黄色片欧美黄色片| 欧美性猛交╳xxx乱大交人| 成人av在线播放网站| 99热这里只有精品一区 | 日韩 欧美 亚洲 中文字幕| 亚洲五月天丁香| 国产高清三级在线| 在线a可以看的网站| 此物有八面人人有两片| 美女大奶头视频| 亚洲精品色激情综合| 欧美黑人巨大hd| 99久久国产精品久久久| 看免费av毛片| 在线观看66精品国产| 在线观看一区二区三区| 91在线精品国自产拍蜜月 | 亚洲人成伊人成综合网2020| 看黄色毛片网站| 午夜福利在线在线| 18禁国产床啪视频网站| 俄罗斯特黄特色一大片| 男人的好看免费观看在线视频| 精品久久久久久久久久久久久| 熟妇人妻久久中文字幕3abv| 999久久久国产精品视频| 长腿黑丝高跟| 欧美激情在线99| 狂野欧美激情性xxxx| 国产亚洲精品综合一区在线观看| 亚洲 欧美 日韩 在线 免费| 999精品在线视频| 久久久国产欧美日韩av| www.999成人在线观看| 色精品久久人妻99蜜桃| 国产97色在线日韩免费| 亚洲专区字幕在线| 国产毛片a区久久久久| 国产精品综合久久久久久久免费| 久久中文看片网| 免费大片18禁| 男人的好看免费观看在线视频| 欧洲精品卡2卡3卡4卡5卡区| 香蕉久久夜色| 国产乱人视频| 一二三四社区在线视频社区8| 男人和女人高潮做爰伦理| 国产aⅴ精品一区二区三区波| 亚洲av五月六月丁香网| www日本黄色视频网| 床上黄色一级片| 国产真实乱freesex| 色尼玛亚洲综合影院| 老司机在亚洲福利影院| 欧美在线黄色| 色精品久久人妻99蜜桃| 久久伊人香网站| 国产探花在线观看一区二区| 无限看片的www在线观看| 激情在线观看视频在线高清| 国产午夜福利久久久久久| АⅤ资源中文在线天堂| 99热6这里只有精品| 欧美最黄视频在线播放免费| 黄色日韩在线| 亚洲成人中文字幕在线播放| 亚洲av第一区精品v没综合| 可以在线观看的亚洲视频| 精品一区二区三区av网在线观看| 国产精品久久久久久久电影 | 村上凉子中文字幕在线| 国产三级黄色录像| bbb黄色大片| 99久久99久久久精品蜜桃| 免费在线观看影片大全网站| 久久久水蜜桃国产精品网| 性色avwww在线观看| 真实男女啪啪啪动态图| 亚洲中文字幕日韩| 99久久久亚洲精品蜜臀av| 久久天堂一区二区三区四区| 91在线精品国自产拍蜜月 | 久久伊人香网站| 在线观看免费视频日本深夜| 国内揄拍国产精品人妻在线| 亚洲精品国产精品久久久不卡| 日本一二三区视频观看| av黄色大香蕉| 亚洲人成网站在线播放欧美日韩| 黄片小视频在线播放| 精品久久久久久久末码| 99久久综合精品五月天人人| 欧美日本视频| 岛国在线观看网站| 日韩 欧美 亚洲 中文字幕| 国产精品久久视频播放| 色尼玛亚洲综合影院| 亚洲熟女毛片儿| 一个人免费在线观看的高清视频| 日本撒尿小便嘘嘘汇集6| 免费观看精品视频网站| 亚洲精品一卡2卡三卡4卡5卡| 一区福利在线观看| 91在线精品国自产拍蜜月 | 两人在一起打扑克的视频| 国产成+人综合+亚洲专区| 无遮挡黄片免费观看| 国产精品九九99| 国产真实乱freesex| 中文字幕av在线有码专区| 亚洲成人久久性| 国产三级黄色录像| 久9热在线精品视频| 精品久久久久久久人妻蜜臀av| 色哟哟哟哟哟哟| 黄色女人牲交| 精品欧美国产一区二区三| 黑人欧美特级aaaaaa片| 欧美又色又爽又黄视频| 俄罗斯特黄特色一大片| 精品午夜福利视频在线观看一区| 亚洲av成人一区二区三| 欧美+亚洲+日韩+国产| 日本免费a在线| 中出人妻视频一区二区| 国产伦一二天堂av在线观看| 国产又色又爽无遮挡免费看| 久久九九热精品免费| 国产伦精品一区二区三区四那| 国产三级中文精品| a级毛片在线看网站| 少妇的逼水好多| 天天躁狠狠躁夜夜躁狠狠躁| 最近最新中文字幕大全免费视频| 国产av一区在线观看免费| 免费在线观看影片大全网站| 久久中文字幕人妻熟女| 女人被狂操c到高潮| 白带黄色成豆腐渣| 深夜精品福利| 日本一二三区视频观看| 中文字幕久久专区| 久久精品人妻少妇| 精品国产乱码久久久久久男人| 免费在线观看成人毛片| 色老头精品视频在线观看| 少妇丰满av| 毛片女人毛片| 亚洲成av人片在线播放无| 欧美av亚洲av综合av国产av| 黄片大片在线免费观看| 国产亚洲精品综合一区在线观看| 18禁观看日本| 黄色视频,在线免费观看| 成年女人毛片免费观看观看9| 婷婷亚洲欧美| 国产精品,欧美在线| 国产69精品久久久久777片 | 亚洲熟妇熟女久久| 在线免费观看不下载黄p国产 | 宅男免费午夜| 两人在一起打扑克的视频| 三级男女做爰猛烈吃奶摸视频| 琪琪午夜伦伦电影理论片6080| 18禁黄网站禁片午夜丰满| 在线国产一区二区在线| 日韩欧美一区二区三区在线观看| 亚洲七黄色美女视频| 欧美一级a爱片免费观看看| 悠悠久久av| 成人午夜高清在线视频| 怎么达到女性高潮| 午夜影院日韩av| 日本 欧美在线| 在线国产一区二区在线| 九九久久精品国产亚洲av麻豆 | 人妻丰满熟妇av一区二区三区| 18禁观看日本| 男女视频在线观看网站免费| 国产探花在线观看一区二区| 日本成人三级电影网站| 亚洲av成人av| 国产成人欧美在线观看| 一区二区三区高清视频在线| а√天堂www在线а√下载| 特大巨黑吊av在线直播| 99国产综合亚洲精品| 日本黄色片子视频| 日本成人三级电影网站| 久久久精品欧美日韩精品| av黄色大香蕉| 男女午夜视频在线观看| 母亲3免费完整高清在线观看| 国产成人av激情在线播放| 午夜影院日韩av| 欧美成狂野欧美在线观看| 国产成+人综合+亚洲专区| 国产高清视频在线观看网站| 成人亚洲精品av一区二区| 又粗又爽又猛毛片免费看| 国产黄色小视频在线观看| 精品久久久久久久久久久久久| ponron亚洲| 国产伦人伦偷精品视频| 高清在线国产一区| 成人av一区二区三区在线看| 99久久无色码亚洲精品果冻| 久久久精品大字幕| 国产精品久久久久久亚洲av鲁大| 99riav亚洲国产免费| 亚洲av成人一区二区三| 亚洲美女视频黄频| 在线观看免费午夜福利视频| 国产免费男女视频| 亚洲精品色激情综合| 国产单亲对白刺激| 午夜两性在线视频| 亚洲精品粉嫩美女一区| 欧美一级毛片孕妇| 国产一级毛片七仙女欲春2| 嫩草影院精品99| 午夜成年电影在线免费观看| 国产午夜福利久久久久久| 露出奶头的视频| 国产精品一区二区免费欧美| 亚洲乱码一区二区免费版| 亚洲成av人片免费观看| 亚洲成av人片在线播放无| 色综合亚洲欧美另类图片| 日本在线视频免费播放| 日韩欧美一区二区三区在线观看| 黄色丝袜av网址大全| 色综合亚洲欧美另类图片| 高清在线国产一区| 亚洲av免费在线观看| 久久久久久久久免费视频了| 亚洲精品国产精品久久久不卡| 看片在线看免费视频| 成人特级av手机在线观看| 淫妇啪啪啪对白视频| 男女床上黄色一级片免费看| 麻豆一二三区av精品| a级毛片a级免费在线| 国产亚洲精品综合一区在线观看| 黄频高清免费视频| 一个人看的www免费观看视频| 国模一区二区三区四区视频 | 麻豆国产av国片精品| 午夜福利高清视频| 色视频www国产| 亚洲国产精品999在线| 免费看日本二区| 男人舔女人下体高潮全视频| 中文在线观看免费www的网站| 亚洲av美国av| 国产激情偷乱视频一区二区| 国内毛片毛片毛片毛片毛片| 特级一级黄色大片| 最近最新免费中文字幕在线| 精品国产亚洲在线| 18禁黄网站禁片午夜丰满| 午夜免费成人在线视频| 人妻久久中文字幕网| 国产伦精品一区二区三区视频9 | 欧美日本视频| 国产成人精品无人区| 免费av不卡在线播放| 久久久久九九精品影院| 国内少妇人妻偷人精品xxx网站 | 国产蜜桃级精品一区二区三区| 丰满人妻熟妇乱又伦精品不卡| 黄片小视频在线播放| 看免费av毛片| 叶爱在线成人免费视频播放| 国产乱人伦免费视频| 高清毛片免费观看视频网站| 国内少妇人妻偷人精品xxx网站 | 亚洲,欧美精品.| 99精品在免费线老司机午夜| 成在线人永久免费视频| 欧美不卡视频在线免费观看| 男人的好看免费观看在线视频| 两个人的视频大全免费| 在线看三级毛片| 精华霜和精华液先用哪个| 精品久久久久久成人av| 亚洲精品色激情综合| 亚洲av中文字字幕乱码综合| 亚洲av第一区精品v没综合| 国语自产精品视频在线第100页| 午夜免费激情av| 99久久久亚洲精品蜜臀av| 老司机午夜十八禁免费视频| 在线观看一区二区三区| 精品一区二区三区av网在线观看| or卡值多少钱| 欧美不卡视频在线免费观看| 两个人的视频大全免费| 国内精品久久久久久久电影| 日韩 欧美 亚洲 中文字幕| svipshipincom国产片| 免费电影在线观看免费观看| 一区二区三区国产精品乱码| 国产精品av久久久久免费| www日本在线高清视频| 国产精品久久视频播放| 老熟妇仑乱视频hdxx| 欧美色欧美亚洲另类二区| 日韩 欧美 亚洲 中文字幕| 午夜久久久久精精品| 哪里可以看免费的av片| 又黄又粗又硬又大视频| 一边摸一边抽搐一进一小说| 国产高清视频在线播放一区| 亚洲第一欧美日韩一区二区三区| 他把我摸到了高潮在线观看| 看片在线看免费视频| 黄频高清免费视频| 中文亚洲av片在线观看爽| 熟女人妻精品中文字幕| 国产真人三级小视频在线观看| 欧美大码av| 琪琪午夜伦伦电影理论片6080| 老司机深夜福利视频在线观看| 免费看十八禁软件| 伊人久久大香线蕉亚洲五| 久久久久国产一级毛片高清牌| 欧美午夜高清在线| 成人精品一区二区免费| 亚洲国产精品成人综合色| 成年女人看的毛片在线观看| 两人在一起打扑克的视频| 男女做爰动态图高潮gif福利片| 免费观看精品视频网站| 国产黄色小视频在线观看| 国产亚洲av高清不卡| 国产美女午夜福利| 男人舔女人下体高潮全视频| 看免费av毛片| 亚洲精品中文字幕一二三四区| 亚洲中文字幕日韩| 欧美色欧美亚洲另类二区| 久久精品aⅴ一区二区三区四区| 啦啦啦观看免费观看视频高清| 国产精品自产拍在线观看55亚洲| 亚洲精品一卡2卡三卡4卡5卡| 一级作爱视频免费观看| 美女免费视频网站| 夜夜躁狠狠躁天天躁| 国产精品精品国产色婷婷| 黄色成人免费大全| 热99re8久久精品国产| 男人舔女人下体高潮全视频| 啦啦啦韩国在线观看视频| 日本免费a在线| 成人特级黄色片久久久久久久| 很黄的视频免费| 黄色女人牲交| 人人妻人人看人人澡| 夜夜爽天天搞| 日韩人妻高清精品专区| 蜜桃久久精品国产亚洲av| 床上黄色一级片| 国内精品一区二区在线观看| 国产淫片久久久久久久久 | 色老头精品视频在线观看| 男人和女人高潮做爰伦理| 村上凉子中文字幕在线| xxx96com| 一个人免费在线观看电影 | 国产一区在线观看成人免费| 久久久精品大字幕| 精品免费久久久久久久清纯| 三级国产精品欧美在线观看 | 国产伦精品一区二区三区视频9 | 免费av不卡在线播放| 熟女人妻精品中文字幕| 国产主播在线观看一区二区| 免费一级毛片在线播放高清视频| 男女视频在线观看网站免费| 伊人久久大香线蕉亚洲五| 欧洲精品卡2卡3卡4卡5卡区| 视频区欧美日本亚洲| 国产精品1区2区在线观看.| 狂野欧美白嫩少妇大欣赏| 久久亚洲精品不卡| 国产激情偷乱视频一区二区| 亚洲国产欧美一区二区综合| 男人舔奶头视频| 国产高清videossex| 嫩草影院精品99| 99视频精品全部免费 在线 | 久久精品国产清高在天天线| 亚洲无线观看免费| 欧美激情在线99| 国产精品1区2区在线观看.| 国产亚洲av高清不卡| 国内毛片毛片毛片毛片毛片| 国产伦在线观看视频一区| 舔av片在线| 国产精品乱码一区二三区的特点| 黄色视频,在线免费观看| 中文资源天堂在线| 国产精品av久久久久免费| 特级一级黄色大片| 嫩草影院入口| 此物有八面人人有两片| 午夜久久久久精精品| 香蕉丝袜av| 三级国产精品欧美在线观看 | 欧美又色又爽又黄视频| 欧美中文综合在线视频| 一个人免费在线观看的高清视频| 国产视频一区二区在线看| 桃色一区二区三区在线观看| 五月玫瑰六月丁香| 国产精品一区二区三区四区免费观看 | 美女 人体艺术 gogo| 国产伦精品一区二区三区视频9 | 国产视频一区二区在线看| 最新在线观看一区二区三区| 亚洲熟妇中文字幕五十中出| 手机成人av网站| 亚洲人与动物交配视频| 丰满的人妻完整版| 成人永久免费在线观看视频| 久久久色成人| 久久亚洲精品不卡| 欧美乱码精品一区二区三区| 免费观看的影片在线观看| 五月伊人婷婷丁香| 少妇的逼水好多| 午夜免费激情av| 精华霜和精华液先用哪个| 天堂av国产一区二区熟女人妻| 一边摸一边抽搐一进一小说| 天天躁狠狠躁夜夜躁狠狠躁| 丁香六月欧美| 18禁裸乳无遮挡免费网站照片| 亚洲国产精品sss在线观看| 国产欧美日韩精品一区二区| 成人无遮挡网站| 亚洲欧美日韩无卡精品| 窝窝影院91人妻| 这个男人来自地球电影免费观看| 成人三级黄色视频| 国产欧美日韩精品一区二区| 这个男人来自地球电影免费观看| 97人妻精品一区二区三区麻豆| 国产伦精品一区二区三区四那| 又大又爽又粗| av国产免费在线观看| 国产成人精品久久二区二区免费| 日本一二三区视频观看| 一a级毛片在线观看| netflix在线观看网站| 后天国语完整版免费观看| 国产伦精品一区二区三区视频9 | 成人特级黄色片久久久久久久| 亚洲欧美精品综合久久99| 90打野战视频偷拍视频| 动漫黄色视频在线观看| 国产精品久久电影中文字幕| 亚洲九九香蕉| 日韩欧美 国产精品| 亚洲天堂国产精品一区在线| 两个人视频免费观看高清|