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

旋轉二維激子極化激元凝聚渦旋疊加態的動力學特性*

2020-12-14 04:57:46吳昊任元劉通王元欽刑朝洋
物理學報 2020年23期
關鍵詞:體系

吳昊 任元? 劉通 王元欽 刑朝洋

1) (航天工程大學宇航科學與技術系, 北京 101400)

2) (航天工程大學, 激光推進及其應用國家重點實驗室, 北京 101400)

3) (北京航天控制儀器研究所, 北京 100094)

研究了二維激子極化激元凝聚正反渦旋疊加態在半導體微腔極化激元波色愛因斯坦凝聚(Bose–Einstein condensate, BEC)體系旋轉情形下的穩定性和動力學特性. 在體系旋轉情形下對單分量Gross-Pitaevskii 方程進行重構, 利用四階龍格庫塔方法和時域有限差分方法構建數值模型. 利用實時演化方法研究在體系旋轉的情況下, 不同拓撲荷數的正反渦旋疊加態的實時演化過程及穩態局域粒子數和體系旋轉角速率之間的關系. 研究了渦旋疊加態激發區域的旋轉速率與體系旋轉速率的關系, 并闡明了體系的旋轉速率對渦旋疊加態相位穩定性的影響機理.研究表明, 半導體微腔極化激元BEC 體系的旋轉速率對激子極化激元凝聚疊加態的演化過程及其動力學特性有重要影響.

1 引 言

對于激子極化激元, 尤其是對于微腔激子極化激元的研究, 近來發展十分迅速[1?7]. 微腔激子極化激元是一種光激發, 是部分腔光子和部分量子阱激子的量子態. 半導體平板微腔引入, 實現了實驗上對激子極化激元的反交叉色散行為的觀測[8,9].觀測結果表明, 色散曲線在零波失附近區域有一個明顯的下陷, 這表明微腔激子極化激元的有效質量非常輕, 很容易實現玻色-愛因斯坦凝聚(Bose–Einstein condensate, BEC). 而微腔激子極化激元發生凝聚的臨界溫度可以達到幾開爾文, 而且有可能實現更高溫度下的凝聚, 尤其是近年來人們發現半導體微腔中的激子極化激元凝聚體系(exciton-polariton condensates)有望在室溫下實現BEC[10?13],更給BEC 的應用從實驗室走向工程實踐提供了可能. 而激子極化激元凝聚體系中自發及可操控的量子化渦旋可用于量子導航、量子計算、量子傳輸等諸多領域, 具有十分廣闊的研究價值[14?21].

其中, 基于激子極化激元BEC 量子化渦旋構建波-粒渦旋陀螺儀在量子導航領域有極大的潛在應用價值. 基于此, 本文提出了實現以BEC 中量子化渦旋疊加態的Sagnac 干涉效應為基礎的波-粒渦旋陀螺儀的科學設想[22]. 本質上波-粒渦旋陀螺是一種量子渦旋陀螺, 它利用光與物質相互作用形成的渦旋疊加態作為陀螺效應產生的載體. 而實現這一科學設想的首要前提是對BEC 中渦旋疊加態的陀螺效應加以驗證, 因此本文對旋轉狀態下激子極化激元BEC 中渦旋疊加態的動力學特性進行了研究.

BEC 量子化渦旋的存在形式有單個渦旋, 渦旋陣列和正反渦旋疊加態. 渦旋疊加態是由具有相同拓撲荷數量子數的正反渦旋疊加形成渦旋疊加態. 2012 年Thanvanthri 等[23]研究了基于原子BEC的正反渦旋疊加態在超穩物質波陀螺中的應用, 并首次提出了渦旋干涉陀螺的概念. 2015 年, Padhi等[24]研究了激子極化激元凝聚體系中單渦旋在旋轉參考系中形成渦旋陣列中渦旋個數隨體系參數的變化規律. 2016 年, Dai 等[19]進一步研究了扁平勢和無序勢中渦旋疊加態的穩態結構及其Sagnac效應, 理論闡明了渦旋干涉陀螺中Sagnac 干涉相位, 拓撲荷數量子數和半導體微腔BEC 體系旋轉角速率三者之間的關系. 在此基礎上, 2018 年, 任元等[25]通過數值計算的方法初步驗證了在量子渦旋陀螺應用場景下BEC 單支渦旋的動力學特性. 由于波-粒渦旋陀螺是一個嶄新的前沿性概念, 因此對于其陀螺效應的實現載體—激子極化激元渦旋疊加態的旋轉動力學特性的研究很少, 且人們并未給出渦旋疊加態的演化規律與體系旋轉速率間的關系, 也未深入探討體系旋轉角速率、渦旋拓撲荷數與渦旋疊加態的瞬時角速率之間的關系, 而上述問題正是實現波粒渦旋陀螺所必須解決的基礎性理論問題.

針對上述問題, 本文利用Runge-Kutta 有限差分方法求解耗散體系的單分量Gross-Pitaevskii(GP)方程, 研究了處于無序勢中的二維激子極化激元凝聚正反渦旋疊加態在半導體微腔BEC 體系旋轉情形下的穩定性和動力學特性. 通過引入具有時間分量的旋轉坐標系, 對GP 方程進行重構, 得到具有科里奧利項的刻畫旋轉狀態的GP 方程. 利用四階龍格庫塔方法(Runge-Kutta method)和時域有限差分方法(finite difference time domain method, FDTD)構建數值模型. 利用實時演化方法研究在半導體微腔BEC 體系旋轉的情況下, 不同拓撲荷數的正反渦旋疊加態的實時演化過程及穩態局域激子極化激元粒子數和半導體微腔BEC 體系旋轉角速率之間的關系. 研究了渦旋疊加態激發區域的旋轉速率與半導體微腔BEC 體系旋轉速率的關系, 并給出了半導體微腔BEC 體系的旋轉速率對渦旋疊加態相位穩定性的影響. 研究表明, 半導體微腔BEC體系的旋轉速率對體系的演化過程和動力學特性有重要影響. 這對于構建極化激元凝聚系統以及后續開發相關陀螺儀設備有重要的理論指導意義.

本文的結構安排如下: 第2 節給出了描述旋轉情形下耗散體系的無量綱化形式的GP 方程, 局域激子極化激元粒子數, 渦旋疊加態場分布, 相位分的數值計算方法; 第3 節討論了正反渦旋疊加態的演化過程與體系會旋轉角速率的一系列關系; 第4 節討論了體系轉速在超過107rad/s 量級時體系渦旋疊加態無法保持穩定這一現象在機理層面的探究; 第5 節給出了半導體微腔BEC 體系的旋轉速率對極化基元BEC 體系的演化過程和動力學特性的重要影響的相關結論.

2 理論與數值模型建立

單分量模型是直接以下支激子極化激元為研究對象的, 適合于非共振激發, 并探討激子極化激元的Bose-Einstein 凝聚[26]. 考慮單分量模型:

其中,ψ(r) 是所研究的下支激子極化激元場,P(r)是恒定的外加激勵項(泵浦項),g為激子極化激元間的非線性相互作用,V(r) 是物理場所感受到的結構勢阱,γ是體系的耗散參數,η是衡量體系飽和的參數, 該方程也稱為GP 方程, 是描述Bose物理場的主要手段. 這是由于Bose 體系一般是多粒子體系, 在含有相互作用的情況下, 嚴格的全量子方法難以求解. 上面所述GP 方程需要有非零的初始物理場條件. 其中, 位置r可以表示為r=本文研究了將整個體系放置在繞體系幾何中心勻速旋轉的空間中, 渦旋疊加態的演化特性. 為得到可以描述當整個體系處于非慣性空間、圍繞激發區域幾何中心勻速旋轉的情況下, 體系演化的GP方程, 需要對(1)式進行變換. 此時, 將x,y用旋轉坐標系表示[27], 即

這樣, 將x′,y′代入(1)式, 通過計算, 可以得到體系以Ω為角速率旋轉情形下的GP 方程:

對(3)式進行無量綱處理, 可以得到

由于使用的數值模型是基于柱坐標系的, 因此使用了如下變換關系:

對于柱坐標系下的拉普拉斯算符?2, 數值模型中將其用差分形式改寫, 其中角向與徑向分量分別用ρ,φ表示, 具體的差分形式是:

由此可以構建基于四階龍格庫塔方法的GP 方程的差分形式, 即可以對以轉速為Ω旋轉的體系的渦旋疊加態的場分布、相位分布和局域粒子總數的含時演化特性進行分析.

3 渦旋疊加態在旋轉狀態下的動力學演化特性

研究激子極化激元凝聚體系中渦旋疊加態的陀螺效應的前提在于分析非慣性系下體系的動力學特性. 簡化模型起見, 考慮渦旋疊加態體系圍繞其圓心做勻速旋轉的情形, 如圖1 所示. 此時, 以拓撲荷數為±l的雙支渦旋的渦旋疊加態作為GP 方程中激子場分布ψX(r,t) 的初始解. 因此在含時演化的t = 0?/meV 時刻(下文中時刻t的單位均為?/meV)微腔中存在相干疊加圖樣(圖樣呈現花瓣狀, 后文簡稱這種干涉圖樣為“干涉花瓣”,用以描述激子場的分布)為2l的“渦旋疊加態解”.通過參數的改變, 研究“渦旋疊加態解”隨時間演化的動力學特性.

圖1 旋轉狀態下的激子極化激元渦旋疊加態體系Fig. 1. System of exciton polariton condensates on the rotational state.

3.1 體系旋轉對渦旋疊加態激子場密度分布的影響

首先考慮拓撲荷數分別為±2 的雙支渦旋的疊加情形. 雙支渦旋在半導體微腔中發生干涉, 形成關于圓心對稱的4 個花瓣狀自發輻射場分布極大值區域.現令Ω′=3.81×1010Ωrad/s,其中Ω′為體系的旋轉角速率, 而Ω為其無量綱形式, 本文中體系的旋轉角速率均由無量綱的Ω給出. 圖2(a)為Ω= 0 即體系靜止時, 當渦旋疊加態演化至穩態時的激子場分布, 可以發現, 在雙渦旋的干涉效應下, 激子場分布在φ= ±45°和φ= ± 135°位置出現4 個極大區域. 如果使激子極化激元體系處于角 速 率Ω約為 4.0×10?7, 8.0×10?7和1.2×10?6的旋轉狀態下, 則體系的演化狀態則會發生很大改變. 從圖2(b)—圖2(d)的激子場分布可以發現, 體系的旋轉影響了渦旋疊加態的干涉效應, 隨著旋轉角速率的增大, 雙渦旋的干涉作用變弱, 而整個勢阱內被激發的粒子數卻不斷增大, 導致干涉相干與相消位置區域連通, 干涉花瓣逐漸模糊.

圖2 l=±2 時不同旋轉角速率對激發場分布的影響 (a) Ω=0 , t=80?/meV ; (b) Ω=4.0×10?7 , t=80?/meV ; (c) Ω =8.0×10?7 , t=80?/meV ; (d) Ω=1.2×10?6 , t=80?/meVFig. 2. Effects of different angular velocities of the rotation on the exciton field when l=±2 : (a) Ω=0 , t=80?/meV ; (b) Ω =4.0×10?7 , t=80?/meV ; (c) Ω=8.0×10?7 , t=80?/meV ; (d) Ω=1.2×10?6 , t=80?/meV .

旋轉速率與激子極化激元渦旋疊加態的演化速率以及穩態時的局域激子極化激元粒子數有密切關聯. 雙渦旋的拓撲荷數取l=±2 , 旋轉角速率分別為(Ω1,···Ωk,···Ω9)=(1.0×10?5,···k×10?5,···9.0×10?5), 可以得到一組激子極化激元體系的含時演化局域激子極化激元粒子數曲線. 首先, 圖3(a),(b)反映了Ω1=1.0×10?5和Ω9=9.0×10?5時,局域激子極化激元粒子數處于t=80?/meV 時刻的穩態時, 激子極化激元渦旋疊加態的相位分布. 相位突變處為渦旋產生處, 當轉動角速率較小(Ω1=1.0×10?5)時, 相位突變發生于φ= 0°, 90°, 180°和270°處, 而當轉動角速率較大(Ω9=9.0×10?5)時, 體系中的渦旋發生分裂, 穩定的渦旋疊加態無法存在, 即雙渦旋干涉現象逐漸消失, 激子場聯通成一個環狀.

圖3 轉動角速率與渦旋疊加態演化的關系 (a) 當 Ω 1 =1.0×10?5 , t=80?/meV 時, 激子極化激元渦旋疊加態的相位分布;(b) 當 Ω 9 =9.0×10?5 , t=80?/meV 時, 激子極化激元渦旋疊加態的相位分布; (c) 當體系處于不同轉速時激發區域內局域粒子數隨時間的變化Fig. 3. The relationship between angular velocities of the rotation and the evolution of superposition state of vortexes: (a) The phase distribution of the superposition state of exciton polariton vortexes when Ω 1 =1.0×10?5 and t=80?/meV ; (b) the phase distribution of superposition state of exciton polariton vortexes when Ω 9 =9.0×10?5 and t=80?/meV ; (c) the description curve of the relationship between time and quasi-particle number on various speeds of rotation of the system.

又如圖3(c)所示, 進一步通過對局域激子極化激元粒子數隨時間變化的關系可以發現, 當旋轉角速率增大, 激子極化基元體系的含時演化速率也會增大, 旋轉速率與激子極化激元體系達到穩態所需時間負相關. 而在泵浦光形式、化學勢和能量一定的情況下, 體系旋轉速率越高, 激子場密度越大,這說明旋轉會促進激子極化激元的激發, 使穩態時體系的局域極化激元粒子數大大增加. 轉動角速率越大, 渦旋疊加態越容易在更短的時間內達到穩態. 這可能是因為轉動破壞了原有的渦旋疊加態結構, 使大渦旋裂化為更小的渦旋, 這一過程促使光-激子場能量耦合增強. 而在其他參數一定的情況下, 體系旋轉速率越高, 激子場密度越大, 這說明體系旋轉會促進激子極化激元的激發, 使穩態時局域極化激元粒子數增加, 且當旋轉角速率小于Ω5=5.0×10?5時, 激子極化激元體系能夠在相對久的時間內維持弱平衡狀態, 局域激子極化激元粒子數隨時間的漲落是非常微幅的. 而當旋轉速率超過Ω6=6.0×10?5時, 渦旋疊加態體系局域激子極化激元粒子數會迅速飽和, 此后粒子數會隨時間變化發生明顯的漲落, 這表明激子極化激元渦旋疊加態體系在較高的旋轉速率情形下難以維持長久的平衡狀態.

在較小的轉動角速率情況下, 局域粒子數不隨時間產生明顯的漲落, 但渦旋疊加態的干涉位置依舊會隨時間發生變化. 以旋轉角速率Ω=9.0×10?7、拓撲荷數l=±2 和泵浦光為半徑 6 00 μm 的高斯光束為計算條件, 研究處于ρ∈(5,10) μm 的環狀無序勢阱中渦旋疊加態干涉花瓣位置隨時間變化的情況. 如圖4, 依次是t為0, 15, 30, 45, 60,75?/meV 時刻渦旋疊加態的激子場分布情況. 可以發現如下規律:

1)干涉極大區域在ρ方向的分布不隨時間推移而變化, 此處干涉極大區域(花瓣中心位置)在演化過程中始終位于ρ=7.5 μm 處;

2)隨時間推移, 二維空間中的激子極化激元向環狀無序勢阱中心集中, 干涉區域隨時間演化局域激子極化激元粒子數不斷增加, 而勢阱外區域的場密度分布不斷下降, 至穩態時降為零;

3)干涉花瓣在一段時間內維持初始位置, 在一定時刻后開始旋轉, 且旋轉方向在每一時刻都保持一致.

3.2 體系旋轉速率與渦旋疊加態旋轉速率的關系

由3.1 節所述現象之三(干涉花瓣在一段時間內維持初始位置, 在一定時刻后開始旋轉, 且旋轉方向在每一時刻都保持一致)可以引申出一組研究:渦旋疊加態干涉花瓣的轉動角速率是否與體系轉速相關. 為降低計算誤差, 使用高階激發, 取拓撲荷 數l=±6 ,計 算 了Ω=2.0×10?7, 4.0×10?7,8.0×10?7. 記渦旋疊加態干涉花瓣的瞬時角速率為Ωvortex. 圖5(a)所描述的為體系角速率Ω=2.0×10?7時,t=0?/meV 時 刻 和t=100?/meV 時 刻 渦旋疊加態干涉花瓣的位置對比, 而圖5(b)所描述的為體系角速率Ω=8.0×10?7時,t=0?/meV 時刻和t=100?/meV 時刻渦旋疊加態干涉花瓣的位置對比. 在Ω=8.0×10?7時,t=100?/meV 時刻的干涉花瓣相對于初態旋轉了更多的角度, 這一點可以在“花瓣暈影”的空間分布上得到體現.

圖4 旋轉狀態下不同時刻的激子場分布 (a) t=0?/meV ; (b) t=15?/meV ; (c) t=30?/meV ; (d) t=45?/meV ; (e) t=60?/meV ; (f) t=75?/meVFig. 4. The exciton field distribution at different moments in the rotational state: (a) t=0?/meV ; (b) t=15?/meV ; (c) t=30?/meV ; (d) t=45?/meV ; (e) t t= 6 0?/meV ; (f) t=75?/meV .

由于當t≥200?/meV 時體系完全失穩而失去研究意義, 因此計算的演化時間限定在200?/meV以內. 圖5(c)反映了三種不同旋轉角速率下, 渦旋疊加態在演化過程的各時刻的瞬時旋轉角速率. 計算表明, 渦旋疊加態從初態到一個中間態的過程中, 干涉花瓣的轉動角速率為零. 從該中間態開始,干涉花瓣開始隨體系轉動而轉動, 且Ωvortex的總體趨勢是逐漸增大的. 體系旋轉角速率Ω越大, 則渦旋疊加態干涉花瓣越早開始旋轉, 且Ωvortex的增速越快. 而這一過程也與圖4 中的計算結果相吻合.

3.3 體系旋轉對不同拓撲荷數渦旋疊加態的影響

最后, 本文研究了體系旋轉角速率對不同拓撲荷數的渦旋疊加態的影響. 利用尋找干涉極大區域的中心位置的方法, 可以找出始末時刻“花瓣解”極大值的位置, 從而測算出體系在經歷了一段時間的旋轉后, 激子極化激元渦旋疊加態旋轉的角度. 圖6(a)反映了當體系旋轉角速率Ω ∈(2.0×10?7,1.0×10?6) 時, 在t=80?/meV 時刻, 拓撲荷數分別為l=±4 和l=±12 的渦旋疊加態相對于初態轉過的角度. 如前文所得出的結論, 體系轉速越高, 到達穩態時渦旋疊加態相較于初態轉過的角度就 越大.當l=±4 時,渦旋疊加態最終轉動了14.1°, 而l=±12 時渦旋疊加態最終轉動了8.3°.顯然, 拓撲荷數越大, 其渦旋疊加態的位置受體系轉動影響越小. 圖6(b)給出了不同拓撲荷數情況下, 體系旋轉角速率對演化過程產生的影響, 其Y軸表示體系到達穩態所需的時間. 可見, 一定體系轉速情況下, 渦旋疊加態拓撲荷數與其容易受激發的程度成反比, 因此可以推測, 當拓撲荷數較大時渦旋疊加態更易因體系的旋轉而過度激發, 從而失穩. 相反地, 拓撲荷數越小, 渦旋疊加態在旋轉狀態下的穩定性越好. 這種穩定性來源于相位分布的穩定性, 也即是渦旋的穩定性. 當拓撲荷數足夠大時, 渦旋相對更容易被破壞并分裂, 這就表現為體系容易被激發達到飽和狀態.

圖5 激子極化激元渦旋疊加態瞬時轉動角速率與體系轉速的關系Fig. 5. The relationship between the instantaneous angular rate of the superposition state of exciton polariton vortexes and the rotate speed of the system.

圖6 轉動角速率對不同拓撲荷數激子極化激元渦旋疊加態的影響 (a)體系旋轉角速率Ω ∈(2.0×10?7,1.0×10?6) , t= 80?/meV 時 刻, 拓 撲 荷 數分別為l =±4和 l=±12 的渦旋疊加態相對于初態轉過的角度; (b)不同拓撲荷數情況下, 體系旋轉角速率對演化過程產生的影響Fig. 6. Effects of the angular velocities on the superposition state of exciton polariton vortexes with different topological charge number: (a) Angle of rotation of superposition state vortexes to the initial state at the moment of t= 80?/meV with different topological charge of l=±4 and l=±12 in the angular rate range of Ω ∈(2.0×10?7, 1 .0×10?6) ;(b) effect of the system angular rate on the evolution process with different topological charges.

4 討 論

當Ω小于 1 0?7數量級時, 若泵浦強度、體系耗散γ、環形無序勢阱Vext和衡量飽和參數μ選取得當, 則渦旋疊加態會保持長時間的穩態, 即粒子總數的漲落處于微幅區間中, 渦旋疊加態干涉位置不隨時間改變.而激子極化激元凝聚的超流特性也從理論上支撐在慣性系中“干涉花瓣”不隨體系轉動而發生偏轉的設想.然而種種計算表明當體系旋轉角速率Ω大于10?7數量級時,渦旋疊加態的演化終究會導致激子極化激元干涉圖樣發生偏轉.可能的解釋是,弱平衡下的激子極化激元體系經過足夠長時間的演化,其初始的渦旋結構會被破壞掉,而渦旋結構的分裂會導致干涉作用漸漸不明顯,最后激子場均勻分布于環形無序勢阱中,為化學勢與泵浦光所束縛.而體系轉動加速了渦旋結構裂化衰減的過程,半導體微腔的材料特性,如無序性和一些晶格特性,使得渦旋結構的裂化沿著旋轉的切線方向進行,從而導致在有限的計算時間內就發現了干涉圖樣的旋轉.

5 結 論

研究表明,半導體微腔的旋轉速率對激子極化激元渦旋疊加態的演化過程及其動力學特性有重要影響.微腔的旋轉會加快激子極化激元渦旋疊加態的演化速度,顯著提升結構勢阱中激子場密度分布.過快的旋轉角速率會破壞激子極化激元凝聚體系的激發-耗散平衡,使體系中的局域激子極化激元粒子數發生大幅漲落,進而破壞原有的渦旋結構,導致渦旋疊加態失穩.研究表明,處于旋轉狀態下的激子極化激元凝聚體系,其渦旋疊加態并非從初始態起就跟隨半導體微腔一同旋轉,而是在演化至結構勢阱中場密度分布趨于飽和后才開始發生旋轉,且渦旋疊加態的轉動速率會不斷增加,而這種增加與體系轉速是正相關的.研究表明,不同拓撲荷數的渦旋疊加態,演化特性受體系旋轉的影響不同,高拓撲荷數渦旋疊加態的旋轉效應明顯弱于低拓撲荷數的情形,拓撲荷數越高渦旋疊加態隨體系旋轉的現象越不明顯.然而,相同體系轉速情況下高拓撲荷數渦旋疊加態的演化速度明顯大于低拓撲荷數的情形,具有高不穩定性.而當體系的轉速ω≤103rad/s量級時,激子極化激元渦旋疊加態會在旋轉的體系中保持長時間的穩定且幾乎不隨體系發生轉動,這就意味著在絕大多數的量子渦旋陀螺儀的使用場景中,疊加態渦旋都具備陀螺效應.

猜你喜歡
體系
TODGA-TBP-OK體系對Sr、Ba、Eu的萃取/反萃行為研究
“三個體系”助力交通安全百日攻堅戰
杭州(2020年23期)2021-01-11 00:54:42
構建體系,舉一反三
探索自由貿易賬戶體系創新應用
中國外匯(2019年17期)2019-11-16 09:31:14
常熟:構建新型分級診療體系
中國衛生(2015年12期)2015-11-10 05:13:40
如何建立長期有效的培訓體系
現代企業(2015年1期)2015-02-28 18:43:18
E-MA-GMA改善PC/PBT共混體系相容性的研究
汽車零部件(2014年5期)2014-11-11 12:24:28
“曲線運動”知識體系和方法指導
加強立法工作 完善治理體系
浙江人大(2014年1期)2014-03-20 16:19:53
日本終身學習體系構建的保障及其啟示
主站蜘蛛池模板: 无码福利日韩神码福利片| 91精品网站| 久久超级碰| 欧美日韩另类在线| 婷婷色狠狠干| 欧美亚洲一二三区| 国产福利一区二区在线观看| 久青草免费在线视频| 多人乱p欧美在线观看| 国产激爽大片在线播放| 亚洲福利一区二区三区| 婷婷丁香在线观看| 久久综合亚洲鲁鲁九月天| 成年女人a毛片免费视频| 亚洲人成人无码www| 国产亚洲精久久久久久久91| 亚洲国产中文综合专区在| 色135综合网| 大香网伊人久久综合网2020| 黄色网址免费在线| h网站在线播放| 中文字幕日韩久久综合影院| 午夜一区二区三区| 一区二区三区四区日韩| 亚洲国产综合精品一区| 精品国产成人av免费| 日韩 欧美 国产 精品 综合| 人妻无码一区二区视频| 国产午夜不卡| 久久亚洲AⅤ无码精品午夜麻豆| 欧美丝袜高跟鞋一区二区| 亚洲成人动漫在线观看| 精品国产乱码久久久久久一区二区| а∨天堂一区中文字幕| 少妇精品网站| 欧美一级99在线观看国产| 无码一区二区三区视频在线播放| 久久精品一品道久久精品| 亚洲精品中文字幕午夜| 国产免费黄| 亚洲欧美极品| 91麻豆国产在线| 日本一区高清| аv天堂最新中文在线| 中文无码影院| 国产在线无码一区二区三区| 中文无码影院| 2021亚洲精品不卡a| 国产三级视频网站| 亚洲精品无码人妻无码| 欧洲亚洲欧美国产日本高清| h网址在线观看| 精品91视频| 精品自拍视频在线观看| 国产精品免费电影| 亚洲人成色在线观看| 日韩毛片免费| 亚洲中文字幕av无码区| 无码一区中文字幕| 日本精品一在线观看视频| 亚洲欧美成aⅴ人在线观看| 午夜a级毛片| 99久久精品免费观看国产| 成人在线观看一区| 国产激情在线视频| 亚洲最新在线| 波多野结衣无码视频在线观看| 国产精品亚洲一区二区三区z| 精品国产网| 国产极品美女在线播放| 午夜久久影院| 亚洲国产黄色| 精品人妻系列无码专区久久| 手机在线国产精品| 亚洲成A人V欧美综合| 亚洲综合婷婷激情| 中美日韩在线网免费毛片视频| 99久久99这里只有免费的精品| 不卡国产视频第一页| 国产精品分类视频分类一区| 国产毛片高清一级国语| 九九视频免费在线观看|