999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

Al+離子3s21S0→3s3p3,1P躍遷同位素偏移的理論研究?

2018-03-27 06:12:10張婷賢李冀光劉建鵬
物理學報 2018年5期
關鍵詞:關聯質量模型

張婷賢 李冀光 劉建鵬

1)(北京應用物理與計算數學研究所,北京 100088)

2)(西北師范大學物理與電子工程學院,蘭州 730070)

3)(國防科技大學文理學院,長沙 410073)

(2017年10月19日收到;2017年12月14日收到修改稿)

1 引 言

在核物理研究中,原子核的均方根半徑〈r2〉是研究原子核性質的重要物理參數,對同一元素該參數沿同位素鏈的變化反映原子核內相互作用的變化規律.因此,對不同同位素均方半徑的確定不僅有助于發現奇異核,而且有助于理解原子核內的多體相互作用從而發展核結構理論[1].從原子物理的角度出發,結合對原子(或離子)同位素偏移的精密測量和原子結構高精度計算是提取原子核均方根半徑的有效途徑之一[2].該方法的優勢在于其提取得到的核參數與核模型無關,因此可以用于檢驗核模型.特別是近年來共線激光光譜技術的發展已經可以測量產量較低的不穩定核素的同位素偏移.但是在理論計算方面,由于同位素偏移因子對電子關聯效應的敏感性較強,在有限的計算資源下難以對多電子體系的計算精度進行準確地評估.

近年來,對某些體系的大規模計算中發現,在同位素質量偏移因子(包擴正規質量偏移因子和特殊質量偏移因子)和躍遷能隨電子關聯的收斂過程中,這兩個物理量的收斂性之間存在線性關系.以Carette和Godefroid對C原子和C?離子的研究為例,當參考組態基組固定時,隨著活動空間中關聯軌道的增加,C原子的2p23P態的能量本征值和特殊質量偏移因子均逐漸減小,并近似成線性關系.在活動空間大小一樣、參考組態基組逐漸擴大的情況下,特殊質量偏移因子逐漸增大,能量本征值逐漸減小,二者的收斂過程仍然近似成線性關系.這反映出C原子的2p23P態的能量本征值和特殊質量偏移因子隨電子關聯效應的收斂過程近似成線性關系[3].類似的情況還在S?離子、Cl原子,Cl?離子中被發現[4,5].后來,他們又在對中性Cu的4s躍遷的研究工作中發現,該躍遷的質量偏移因子與躍遷能隨電子關聯的收斂性之間存在線性關系,并且指出在他們的計算模型下,計算所得的質量偏移因子的誤差與躍遷能和實驗結果之間的誤差有一定的關聯[6].最近,Liu等[7]在B+體系中同樣發現了相似的線性關系,而且用該線性關系對計算結果的不確定度進行了評估,其評估結果與通過評估電子關聯效應的收斂性及Breit相互作用所得不確定度基本一致.但是,到目前為止,還沒有對這些線性關系普適證明的報道.另外,對于大部分體系,不論是實驗還是理論,給出躍遷能的結果要相對容易一些,如果同位素偏移因子與躍遷能的收斂性之間確實存在線性關系,這必然有助于對今后理論計算結果不確定度的評估.

本文在之前對Al+離子兩個躍遷的同位素偏移因子精細計算的基礎上,進一步細致地討論了同位素偏移因子隨著電子關聯的收斂過程.發現在大規模計算中,同位素偏移因子的收斂過程和躍遷能的收斂過程存在線性關系.進一步,我們嘗試利用此線性關系對之前的計算結果給出不確定度.

2 理論方法

2.1 MCDHF方法和RCI方法

多組態Dirac-Hartree-Fock方法,簡稱MCDHF方法[8],用具有相同宇稱P、總角動量J以及J的z方向的分量MJ的組態波函數Φ的線性組合表示原子態的波函數Ψ,即

式中ci是混合系數,γi是描述組態波函數Φ所需其余量子數的統一標記.組態波函數Φ又是單電子Dirac軌道波函數的乘積的線性組合.理論上,(1)式中組態波函數Φ的個數應該是無數個時才能得到真實的原子態波函數.但是,在實際計算中只能選取有限多個組態波函數Φ.因此,組態波函數Φ的個數以及選取條件就決定著對電子關聯效應的描述程度.在MCDHF方法中,基于變分原理,通過自洽場計算同時優化混合系數和軌道波函數使原子態的能量達到最低.在得到單電子軌道基后,可以通過進一步擴大組態空間,利用相對論組態相互作用(RCI)的方法包含更多的電子關聯,但此時單電子軌道波函數是固定的,只有混合系數可以被繼續優化.另外,在低頻近似下的Breit相互作用也可以在RCI中被包括[9]:

2.2 同位素偏移

真實的原子核是由質子和中子構成的具有一定質量M的有限量子體系,其核電荷分布由原子核電荷均方半徑描述.因為不同同位素的中子數不同,所以不同同位素具有不同的質量數和核電荷分布,從而導致原子譜線的偏移,這一偏移量被稱之為同位素偏移(IS).根據起因,將同位素偏移分為由于原子核的有限質量引起的質量偏移(MS)和由于原子核的有限體積引起的場偏移(FS)[10,11].對于質量偏移部分,近似到α2(m/M)階,包含相對論的原子核核反沖效應的質量偏移的哈密頓量如(3)式,其單體算符和兩體算符分別對應著正規質量偏移(NMS)和特殊質量偏移(SMS)[12?15],

在此基礎上,將正規質量偏移因子KRNMS和特殊質量偏移因子KRSMS定義為[16]:

對于一個躍遷而言,同位素A和A′之間的質量偏移因子ΔKRMS為躍遷中上(u)下(l)能級之間質量偏移因子差值[17],即

3 計算模型

物理量的計算精度主要取決于對電子關聯的描述程度.在多組態Dirac-Hartree-Fock理論框架下,我們采用活動空間方法捕獲電子關聯效應,因此,組態空間的構建是計算的關鍵所在.在本文的計算中,組態空間是由從參考組態的單雙激發所產生的組態構成.在表1和表2中給出了各個模型的參考組態的基組.為了系統地描述電子關聯,根據微擾理論將電子關聯劃分為一階關聯和高階關聯,其中一階電子關聯效應用從單參考組態單雙激形成的組態空間描述.將1s,2s,2p軌道看作原子芯軌道(C),3s和3p軌道看作價軌道(V),相應地,一階電子關聯可以劃分為價電子之間的關聯(VV關聯),價電子與原子芯電子之間的關聯(CV關聯),以及原子芯電子之間的關聯(CC關聯).待一階關聯描述充分后,將一階關聯組態空間中權重較大的組態加入參考組態基組,用從多參考組態基組的單雙激發形成的組態空間描述高階電子關聯效應.因為從多參考組態的單、雙激發就相當于從單參考組態的限制性的三、四激發,而且我們選擇的是權重較大的組態加入參考組態基組,實質是優先考慮了較重要的高階電子關聯效應,所以在有限的計算資源下大大提高了計算效率.另外,為了更真實地反映各個能級的波函數,計算中對奇偶宇稱的波函數分別進行優化.詳細的計算模型在之前的工作[18]中已經介紹過,下面僅做簡要說明.

首先,在自洽場計算中同時包含了VV和CV關聯.這里的組態空間是由從單參考組態(即3s2和3s3p分別作為奇偶宇稱的參考組態)限制性的單雙激發到由關聯軌道(虛軌道)構成的活動空間產生組態形成,其中對單雙激發的限制條件為:原子芯軌道上至多有一個電子能被激發到關聯軌道.自洽場計算的起點是用Dirac-Hartree-Fock方法優化光譜軌道(占據軌道),并在隨后的計算中保持固定.為了同位素偏移因子能夠穩定收斂,如表1和表2中所示,關聯軌道是逐層加入的,而且每次只優化最新加入的軌道.在前期的測試計算中,我們發現部分高角動量的關聯軌道,如7i,7h,8g等軌道,對同位素偏移因子的影響很小,因此為了提高計算效率,這些軌道不包含在活動空間中.我們把上述計算模型被標記為C1sV-nl,其中n和l分別代表最新加入的關聯軌道的最大主量子數和軌道角量子數.

當自洽場計算的結果收斂之后,將單電子軌道波函數固定,通過進一步擴大組態空間在RCI計算中包括剩余的電子關聯效應.首先在RCI計算中包括的是2s和2p軌道的電子關聯.此時組態空間中被加入的是從單參考組態中2s,2p軌道上雙激發到最大活動空間所形成的組態.相應地,此時的計算模型用CC2s標記.至此,Al+離子中除CC1s關聯外的一階電子關聯效應均已被包括.隨后,逐步把一階關聯中權重大于0.05,0.02,0.01的組態加入參考組態基組再次逐步擴大組態空間,從而包擴了2s,2p軌道和3s,3p軌道之間的高階關聯效應.將這部分計算模型依次標記為MR1,MR2,MR3.在MR1模型中,{2s22p63s2;2s22p63p2}和{2s22p63s3p;2s22p63p3d}分別作為奇偶宇稱的參考組態.這里的關聯軌道仍然是逐層添加,直至該模型下的計算結果達到收斂.待MR1模型下的計算結果收斂后,將{2s22p43s24p2;2s22p43s23d2;2s22p53s3p3d}和{2s22p63s5p;2s22p63s4p;2s22p63d5p}分別加入奇偶宇稱的參考組態基組中,構成MR2模型中的參考組態基組. 同理,MR2模型下的計算結果收斂后,將{2s22p43s23p4p;2s22p43s23d4d;2s2p53s24s4p;2s22p43s24p5p}和{2s22p43s3p4p2;2s22p43s3p4d2}加入參考組態基組,即MR3.最后,用RCI計算評估了Breit作用對Al+離子的同位素偏移因子的貢獻,并將其作為計算結果的不確定度的一部分.本文計算用Grasp2K程序包[19?21]完成.

圖1 Al+離子3s3p3,1P→3s21S0躍遷的躍遷能(單位:cm?1)隨電子關聯的變化Fig.1. Transition energies(in cm?1)of the 3s3p 3,1P→3s21S0transitions in Al+ion as a function of the computational models.

為了反映計算得到的原子態波函數的品質,在表1和表2中給出了各個模型下組態波函數的個數(NCSF)以及能量本征值(單位:cm?1).同時在圖1中展示了本文所涉及的兩個躍遷的躍遷能(單位:cm?1)隨電子關聯的變化.從表1和表2中三個能級的能量本征值隨組態空間的收斂可以看出,最后一步得到的能量本征值最

低,因此,此時得到的原子態波函數與真實情況最接近.但在圖1中最后計算的的躍遷能與虛線標注的NIST[22]推薦值相比并不是最好的,其誤差為1.35%.這是由于初末態能級的關聯效應存在細微的不平衡,這些誤差需要用更大規模的計算來改進,但這已經超出了我們目前的計算資源.通過系統地考慮電子關聯效應,可以準確地給出物理量的誤差范圍,并呈現在表3中.

表1 各個計算模型下3s21S0態的參考組態基組,活動軌道(AO),組態波函數的個數(NCSF)以及能量本征值(單位:cm?1)Table 1. The reference configurations,the active orbitals(AO),the number of CSFs(NCSF)and the energy(in cm?1)for the 3s21S0state in various correlation models.

表2 各個計算模型下3s3p3,1P態的參考組態基組,活動軌道(AO),組態波函數的個數(NCSF)以及能量本征值(單位:cm?1)Table 2.The reference configurations,the active orbitals(AO),the number of CSFs(NCSF),and the energy(in cm?1)for the 3s3p3,1Pstates in various correlation models.

表2 各個計算模型下3s3p3,1P態的參考組態基組,活動軌道(AO),組態波函數的個數(NCSF)以及能量本征值(單位:cm?1)Table 2.The reference configurations,the active orbitals(AO),the number of CSFs(NCSF),and the energy(in cm?1)for the 3s3p3,1Pstates in various correlation models.

Energy/cm?12 ?53111869 ?53079345 CV1s-4f 4s,4p,3d,4f 600 ?53114653 ?53089871 CV1s-5g 5s,5p,4d,5f,5g 2326 ?53117018 ?53094133 CV1s-6h 6s,6p,5d,6f,6g,6h 5573 ?53117958 ?53095104 CV1s-7g 7s,7p,6d,7f,7g,6h 9860 ?53118256 ?53095425 CV1s-8f 8s,8p,7d,8f,7g,6h 14292 ?53118331 ?53095553 CV1s-9f 9s,9p,8d,9f,7g,6h 19594 ?53118377 ?53095619 CV1s-10f 10s,10p,9d,10f,7g,6h 25766 ?53118399 ?53095643 CV1s-11d 11s,11p,10d,10f,7g,6h 30696 ?53118405 ?53095660 CV1s-12d 12s,12p,11d,10f,7g,6h 36146 ?53118411 ?53095667 CV1s-13d 13s,13p,12d,10f,7g,6h 42116 ?53118413 ?53095670 CC2s 13s,13p,12d,10f,7g,6h 223468 ?53188733 ?53163561 MR1 +{2s22p63p3d} 8s,8p,8d,8f,7g,6h 352702 ?53189063 ?53165717 MR2 +{2s22p63s5p; 5s,7p,5d,5f,5g 400690 ?53189141 ?53165917 2s22p63s3p;2s22p63d5p}MR3 +{2s22p43s3p4p2; 4s,6p,5d,4f 1327012 ?53189741 ?53166574 2s22p43s3p4d2}DF {2s22p63s3p}

4 結果與討論

4.1 單個能級的收斂

圖2 3s2(1S0)和3s3p(1,3Po1)的質量偏移因子(KRNMS和KRSMS(單位:GHz·u))與能量本征值E(單位:cm?1)隨電子關聯的收斂Fig.2.The convergency of the NMS and SMS factor(in GHz·u)of the states 3s2(1S0)and 3s3p(1,3P)versus the eigenvalue energy(in cm?1).

首先研究了3s21S0和3s3p1,3Po1能級的同位素偏移因子和能量本征值隨電子關聯的收斂過程.圖2給出了上述能級的正規質量偏移因子KRNMS(左)、特殊質量偏移因子KRSMS(右)和能量本征值的收斂過程.圖中的δK=|K(a)?K(b)|,δE=|E(a)?E(b)|,其中a代表每一步計算模型,b代表最后一步計算模型,即MR3.圖中實心方點對應的是在自洽場計算中物理量隨VV和CV關聯的收斂過程;空心方點對應在RCI計算中隨CC2s關聯的收斂過程;空心三角則代表在RCI中隨高階關聯的收斂過程.對于能量本征值,我們的計算結果已收斂到0.001%,對于質量偏移因子KRNMS和KRSMS則大約收斂到0.1%,這反映出同位素偏移因子對電子關聯的敏感性較高,對于能量本征值影響很小的電子關聯,對同位素偏移的影響不一定可以忽略.另外,從圖中可以看出,在自洽場過程中與MR過程同位素偏移因子隨能量本征值的收斂斜率不同,這反映出在不同的計算模型下同位素偏移因子和能量本征值的收斂速度不同,但是隨著越來越多的關聯被包括后,能量本征值和同位素偏移因子的收斂速度趨于一致.

圖3 3s2(1S0)和3s3p(1,3P)的質量偏移因子(KRNMS和KRSMS(單位:GHz·u))與能量本征值E(單位:cm?1)相鄰兩步的差值隨電子關聯的變化Fig.3.Convergency of the NMS and SMS factor(in GHz·u)of the states 3s2(1S0)and 3s3p(1,3P)versus the eigenvalue energy(in cm?1)with the adjacent steps.

圖3給出了3s21S0和各能級的正規質量偏移因子KRNMS和特殊質量偏移因子KRSMS在計算過程中每相鄰兩步計算結果的差值隨能量本征值每相鄰兩步的差值的變化.圖中的δK=|K(α)?K(β)|,δE=|E(α)?E(β)|,此時α代表每一步計算模型,β代表α之前一步的計算模型.與圖2的標示方法相同,即實心方點部分對應在自洽場計算中物理量隨VV和CV關聯的收斂過程,空心方點對應在RCI計算中隨CC2s關聯的收斂過程,空心三角代表隨高階關聯的收斂過程.從圖3可以看出,對于單個能級,其同位素偏移因子和能量本征值隨電子關聯的收斂性之間近似存在線性關系,這進一步說明這兩種物理量在每一步收斂過程中都是相關的.當電子關聯對能量本征值的影響減小時對同位素偏移因子的影響也相應地減小.相比于自洽場計算部分,δK和δE之間線性關系在RCI計算中隨CC2s關聯的收斂過程中更明顯.類似的情況同樣出現在Liu等[7]對B+的研究中.在圖3中,物理量隨高階關聯的收斂過程中δK和δE之間線性關系不明顯,主要是由于我們的計算模型中所包含的高階關聯有限.在Liu等對B+的計算模型中直接采用了從單參考組態的三、四激發形成組態空間的方法描述高階關聯,因此對高階關聯的描述很充分.但是受限于計算資源,我們的計算中僅包含了重要的高階關聯,這導致得出的線性關系不及B+的明顯.

圖4 3s21S0→3s3p3,1P躍遷的質量偏移因子(ΔKRNMS和ΔKRSMS(單位:GHz·u))與躍遷能ΔE(單位:cm?1)相鄰兩步的差值隨電子關聯的變化Fig.4.Convergency of the NMS and SMS factor(in GHz·u)of the transitions 3s21S0→ 3s3p3,1Pversus the transition energy(in cm?1)with the adjacent steps.

4.2 躍遷的收斂

圖4給出了3s21S0→3s3p3,1P兩個躍遷的同位素質量偏移因子ΔKRNMS和ΔKRSMS每相鄰兩步的差值隨躍遷能ΔE的每相鄰兩步差值的變化規律.對于躍遷而言,其同位素偏移因子和躍遷能都是該躍遷初末態的相應物理量的差值,由于初末態之間存在相互抵消作用,因此相比于單個能級,躍遷所對應的收斂過程中的線性關系較差.另外,由于在我們的計算中對奇偶宇稱的波函數是分別優化的,不可避免地造成對初末態收斂的不平衡,這種不平衡會破壞躍遷同位素偏移因子和躍遷能之間的線性關系.從圖4中可以看出,在大規模的計算中,躍遷能與同位素偏移因子之間確實存在線性關系,相比于自洽場計算和CC2s關聯部分的收斂過程,高階關聯部分的收斂過程中的線性關系更明顯.對于正規質量偏移因子KNMS在非相對論框架下有標度關系,ΔKNMS=?ν/1823,其中ν是躍遷頻率[23].對于Al+離子我們的計算結果表明相對論的原子核反沖效應對質量偏移因子的影響不超過4%,因此在正規質量偏移因子與躍遷能隨電子關聯的收斂速度基本一致的情況下,也就是電子關聯描述較充分的情況下,正規質量偏移因子與躍遷能的收斂過程呈現線性關系并不意外.對于特殊質量偏移因子,單從其算符形式很難看出線性關系的存在,但它本質上是與電子的動能相關的,因此推測可能存在線性關系,并通過大規模數值計算予以驗證.

4.3 不確定度評估

表3中給出了用兩種評估方法得到的同位素偏移因子的不確定度.一種是通過評估電子關聯效應和Breit作用所得結果,一種是用圖4中高階關聯所對應的數值點線性擬合得到的直線結合我們計算的躍遷能與實驗結果誤差所得的不確定度.在本文的計算模型下,兩種方法所得不確定度差距較大,其原因主要由于在我們的計算中對高階關聯的捕獲有限,因此高階關聯所對應的數值個數有限導致圖4中線性關系的擬合誤差較大.但是,在計算資源更優越的條件下,我們可以對高階關聯的描述程度更充分,從而得到的線性關系也會更準確,相應地,用該線性關系評估得到的不確定度也更為可靠.

表3 兩種不確定度評估方式下所得結果對比,line代表用本文線性關系所得結果,correlation代表之前工作中通過評估各個電子關聯效應貢獻所得結果Table 3.Comparing of the uncertainties which are evaluated by the different method,the “line”means the results evaluated using the linearity in this work,the“correlation”means the uncertainties in our previous work.

5 結 論

通過研究Al+中3s21S0→3s3p3,1P躍遷的同位素偏移因子和躍遷能隨電子關聯的收斂性,發現在Al+體系中同位素偏移因子和躍遷能隨電子關聯的收斂性之間確實存在線性關系.而且通過這一線性關系可以根據理論計算的躍遷能和實驗結果的符合程度實現對同位素偏移因子的理論計算結果的不確定度予以評估.由于計算資源有限,在本文的計算模型下用兩種方法評估所得的不確定度相差較大,而該偏差可以通過在計算模型中包含更充分的高階關聯的方法來減小.對于暫時沒有同位素偏移實驗結果的多電子體系,上述線性關系提供了一種有效評估不確定度的方法.鑒于同位素偏移因子與能量本征值或躍遷能的收斂性之間的線性關系已經在多個體系中被發現,我們猜測這種線性關系并非偶然,但是需要在更多體系中予以驗證.

[1]Campbell P,Moore I D,Pearson M R 2016Prog.Part.Nucl.Phys.86 127

[2]Cheal B,Cocolios T E,Fritzsche S 2012Phys.Rev.A86 042501

[3]Carette T,Godefroid M R 2011Phys.Rev.A83 062505

[4]Carette T,Godefroid M R 2010Phys.Rev.A81 042522

[5]Carette T,Godefroid M R 2013J.Phys.B:At.Mol.Opt.Phys.46 095003

[6]Carette T,Godefroid M R 2016arXiv:1602.06574.

[7]Liu J P,Li J G,Zou H X 2017Chin.Phys.B26 23104

[8]Grant I P 2007Relstivistic Quantum Theory of Atom and Molecules(New York:Spinger)

[9]Li J G,J?nsson P,Godefroid M R,Dong C Z,Gaigalas G 2012Phys.Rev.A86 052523

[10]Tupitsyn I I,Shabaev V M,Crespo López-Urrutia J R,Dragani? I,Soria Orts R,Ullrich J 2003Phys.Rev.A68 022511

[11]Filippin L,Beerwerth R,Ekman J,Fritzsche S,Godefroid M R,J?nsson P 2016Phys.Rev.A94 062508

[12]Palmer C W P 1987J.Phys.B:At.Mol.Opt.Phys.20 5987

[13]Torbohm G,Fricke B,Rosén A 1985Phys.Rev.A31 2038

[14]Filippin L,Godefroid M R,Ekman J,J?nsson P 2016Phys.Rev.A93 062512

[15]Filippin L,Bieroń J,Gaigalas G,Godefroid M R,J?nsson P 2017Phys.Rev.A96 042502

[16]Gaidamauskas E,Nazé C,Rynkun P,Gaigalas G,J?nsson P,Gaigalas G 2011J.Phys.B:At.Mol.Opt.Phys.44 175003

[17]Li J G,Godefroid M R,Wang J G 2016J.Phys.B:At.Mol.Opt.Phys.49 115002

[18]Zhang T X,Xie L Y,Li J G,Lu Z H 2017Phys.Rev.A96 012514

[19]J?nsson P,He X,Fischer C F,Grant I 2007Comput.Phys.Commun.177 597

[20]J?nsson P,Gaigalas G,Bierón J,Fischer C F,Grant I 2013Comput.Phys.Commun.184 2197

[21]Nazé C,Gaidamauskas E,Gaigalas G,Godefroid M R,J?nsson P 2013Comput.Phys.Commun.184 2187

[22]Kramida A,Ralchenko Y,Reader J 2016NIST Atomic SpectraDatabase(Version 5)https://www.nist.gov/pml/atomic-spectra-database

[23]Godefroid M R,Fischer C F,J?nsson P 2001J.Phys.B:At.Mol.Opt.Phys.34 1079

猜你喜歡
關聯質量模型
一半模型
“質量”知識鞏固
“苦”的關聯
當代陜西(2021年17期)2021-11-06 03:21:36
質量守恒定律考什么
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
做夢導致睡眠質量差嗎
奇趣搭配
智趣
讀者(2017年5期)2017-02-15 18:04:18
3D打印中的模型分割與打包
主站蜘蛛池模板: 无码高清专区| 久久综合干| 亚洲女同一区二区| 日本高清有码人妻| 亚洲手机在线| 天天躁夜夜躁狠狠躁躁88| 在线视频亚洲欧美| 中国黄色一级视频| 狠狠久久综合伊人不卡| 久久毛片免费基地| 日本精品视频一区二区| 欧美翘臀一区二区三区| 国产成人免费高清AⅤ| 欧美a在线看| 四虎影视国产精品| 99精品国产高清一区二区| 亚洲人成网站在线播放2019| 午夜毛片免费观看视频 | 国产91色在线| 亚洲男女天堂| 中文字幕人妻无码系列第三区| 在线观看国产精美视频| 色国产视频| 亚洲欧美成人网| 亚洲综合中文字幕国产精品欧美| 亚洲系列中文字幕一区二区| 日韩123欧美字幕| 亚洲欧洲日产国码无码av喷潮| 美女无遮挡免费网站| 国产精品欧美在线观看| www.99在线观看| 日韩精品成人网页视频在线| 日本三级黄在线观看| 一级片一区| 69综合网| 全部毛片免费看| 欧洲免费精品视频在线| 国产成人综合日韩精品无码不卡| 狠狠亚洲五月天| 国产成人亚洲综合a∨婷婷| 五月婷婷导航| av在线5g无码天天| 亚洲69视频| 国产精品999在线| 色婷婷视频在线| 亚洲天堂伊人| 成年人久久黄色网站| 十八禁美女裸体网站| 国产精品黄色片| 亚洲欧美日本国产综合在线| 大学生久久香蕉国产线观看| 奇米影视狠狠精品7777| 刘亦菲一区二区在线观看| 少妇高潮惨叫久久久久久| 四虎永久免费地址在线网站| 久久国产黑丝袜视频| 成年午夜精品久久精品| 亚洲男人的天堂在线观看| 日本高清有码人妻| 在线亚洲天堂| 亚洲激情区| 亚洲水蜜桃久久综合网站 | 欧美www在线观看| 国产大全韩国亚洲一区二区三区| 伊人成人在线视频| 就去吻亚洲精品国产欧美| 欧美亚洲香蕉| 国产高清在线精品一区二区三区| AV无码一区二区三区四区| 六月婷婷精品视频在线观看| 免费播放毛片| 亚洲区视频在线观看| 中文字幕乱妇无码AV在线| 国产av一码二码三码无码| 国产美女在线观看| 欧美成在线视频| 精品91在线| 99精品一区二区免费视频| 亚洲大尺码专区影院| 精品1区2区3区| 三区在线视频| 亚洲免费毛片|