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

Giant and controllable Goos–H¨anchen shift of a reflective beam off a hyperbolic metasurface of polar crystals

2024-01-25 07:13:02TianXue薛天YuBoLi李宇博HaoYuanSong宋浩元XiangGuangWang王相光QiangZhang張強ShuFangFu付淑芳ShengZhou周勝andXuanZhangWang王選章
Chinese Physics B 2024年1期

Tian Xue(薛天), Yu-Bo Li(李宇博), Hao-Yuan Song(宋浩元), Xiang-Guang Wang(王相光), Qiang Zhang(張強),Shu-Fang Fu(付淑芳),?, Sheng Zhou(周勝), and Xuan-Zhang Wang(王選章)

1Key Laboratory for Photonic and Electronic Bandgap Materials,Ministry of Education,and School of Physics and Electronic Engineering,Harbin Normal University,Harbin 150025,China

2Department of Basic Courses,Guangzhou Maritime University,Guangzhou 510725,China

Keywords: Goos–H¨anchen shift,black phosphorus,surface plasmon phonon polaritons,sensitivity,metasurfaces

1.Introduction

Real light beams with finite lateral width reflected at dielectric surfaces appear to be shifted in-plane and out-of-plane,which is not entirely consistent with geometrical optics predictions.The Goos–H¨anchen(GH)shift[1]refers to the beam shift parallel to the plane of incidence, while the Imbert–Fedorov(IF)shift[2]refers to the beam shift perpendicular to the incident plane.Artmann first theoretically described the GH shift in the regime of classical physics in 1948.[3]Optical losses due to optical absorption or scattering cannot be ignored in most materials, particularly metals.As a result, at an ordinary dielectric interface,the GH shift will be the same as the incident wavelength.The development of weak measurement technology has expanded the GH shift’s potential applications,[4]including the accurate measurement of optical properties.[5]Several studies have been conducted on various dielectric surfaces to enhance the GH shift, such as the metasurfaces,[6]topological insulators,[7]Weyl semimetals,[8]and metamaterials (MMs),[9,10]as well as two-dimensional(2D) atomic crystals,[11]and so on.A structure known as subwavelength gratings was recently used to improve the GH shift.[12]Wuet al., for example, demonstrated massive GH shifts in a compound grating waveguide structure that is aided by bound states in the continuum(BICs).[13,14]

Due to their special electronic and optical properties,2D atomic crystals have revealed a variety of physical phenomena since the discovery of graphene.[15]The hyperbolic properties of graphene-based MMs can be achieved by adjusting the chemical potential controlled by the gate voltage or the doping level.[16]As an alternative to graphene, black phosphorus (BP) has received more attention, showing extraordinary potential in many applications, including photodetectors,[17]phase shifters,[18]and absorbers.[19]Unlike graphene,BP processes the puckered honeycomb lattice structure, resulting in the significant in-plane anisotropy.Furthermore,BP has a direct and layer-sensitive bandgap ranging from 0.3 eV to 2 eV,making it a remarkable material capable of achieving a tunable optical response over a broad range of wavelengths.[20]It is worth noting that 2D BP monolayers have been used in the development of polarization-sensitive broadband photodetectors, as well as traditional terahertz (THz) absorbers.[21]Motivated by the emergence of stacked 2D materials with a twist angle, the spatial shifts of the reflective beam from the hBN covered by a rotated BP layer are investigated.[22]Reference[23]investigated the spin Hall effect(SHE)on the surface of a twisted few-layer BP film,and discovered that an in-plane shift occurs and is sensitive to the twist angle.Surface plasmon polaritons(SPPs)can also be generated on the surface of a 2D BP film, which makes it a viable option for designing optoelectronic devices.[24]

Ionic crystals(ICs),which constitute the majority of polar crystals,possess at least one reststrahlen band(RB)where the longitudinal and transverse components of permittivity have the opposite sign.[25]The optical phonon mode of lattice vibration in ICs strongly couples with electromagnetic waves to form surface phonon polaritons(SPhPs)in the THz range.Unlike photonic crystals and metamaterials,the excited SPhPs on the surface of ICs can travel longer distances than SPPs due to the much lower optical losses of ICs.[26]In the infrared to THz range, SPhPs can support sub-diffraction limited nearfield imaging,extraordinary transmission,and so on.[27]

For conventional types of sensors that typically use noble metal like gold or silver to improve sensor sensitivity,[28]it is still a challenging to improve sensor sensitivity due to the low binding ability of gold with biomolecules and the poor chemical stability of silver in the air.Monitoring GH shift signals by controlling the structural parameters and material characteristics is an efficient way to realize the sensing performance.Based on this concept, a sensor that works in the THz range by measuring the GH shift effect is promised.A new dissolved oxygen sensor based on surface plasmon resonance (SPR) and GH shift has been proposed, where the hemoglobin is immobilized on the sensor chip through self-assembled monolayers.[29]A highly sensitive short-range mode resonance sensor made of multilayer structured hyperbolic MMs with a maximum sensitivity of 330 μm/RIU in the near-infrared band has recently been designed.[30]Reference [31] describes a graphene–MoS2heterostructure-coated Au GH shift sensor with a sensitivity of 5.545×105λ0/RIU.The sensibility of a bimetallic sensor based on graphene–hBN heterostructure can be improved by enhancing the GH shift in the infrared band, and a sensibility of 202×105λ0/RIU can be achieved.[32]By precisely controlling the GH shift in the graphene–substrate system at the optical communication band, Zhouet al.created a refractive index sensor with adjustable sensitivity coefficient, which can be increased to±1×108λ0/RIU by adjusting the Fermi energy.[33]In this work,we propose a metasurface structure made of polar crystals with BP-patches.The properties of the metasurface can be altered by adjusting parameters of BP-patches, such as inplane anisotropy,size,carrier density,and layer number.Another important feature is that the surface plasmon phonon polaritons (SPPPs) will be excited at the interface between the BP-patches and the polar crystal, which is coupled by SPPs in BP and SPhPs in the polar crystal.It has been proven that the generation of SPPPs is beneficial to enhance the spatial shift.[34,35]Based on these advantages,the proposed novel sensing method based on the enhanced GH shift exhibits superior sensibility, namely, the maximum sensitivity is up to 6.43×108λ0/RIU,which is almost two orders of magnitude higher than the conventional sensor.We believe that this work is conducive to the design of GH shift sensors in THz range,and may also play an important role in highly sensitive detection.

This paper is organized as follows.Section 2 describes the metasurface structure created by the BP-patches,ZnS crystal,and sensor medium.Using the transfer matrix and equivalent circuit models,we provide a theoretical analysis of the GH shift and the proposed sensor device.Section 3 discusses how numerical simulations are used to verify the accuracy of theoretical models and provide additional physical insights.Finally,we summarize our findings in Section 4.

2.Theoretical model and method

In Fig.1(a), a hyperbolic metasurface based on BPpatches is constructed.The thickness of IC isdand the period of BP-patches is set toDalong thex- andy-axes.A sensing layer(SL)under the IC is utilized to detect the slight changes in its refractive index based on the GH shift.LxandLyindicate the lengths of the BP-patch along thex- andydirections, respectively.The case where the armchair (AC)and the zigzag (ZZ) along thex- andy-directions is referred to as the model-I.The model-II is defined when the directions of the armchair and zigzag are switched, as shown in Fig.1(b).The corresponding equivalent circuit model is plotted in Fig.1(c).We assume that a radiation beam impinges on the hyperbolic metasurface with incident angleθ.With the help of Kubo’s formula, the Drude model of BP conductivity can be described by[36]

wherenBis the carrier density,η,mj,andsjrepresent the relaxation time,the effective mass of electrons,and the different strengths of the interband component along thejdirection,respectively.ωjis the frequency of the interband transitions for thejcomponent.Θis the step function.It is worth noting that the optical conductivity of an ultrathin BP film containing several monolayers(N <6)can be approximated asNσj,[37]whereNrepresents the layer number of BP film.The dielectric constant of the IC as the function of frequency is normally expressed as[38]

whereεhandεlare the high and low frequency permittivity,fTis the transverse lattice vibration frequency, andτdrepresents the damping of optical loss.The RB of IC is shown byif the damping responsible for losses is ignored.The dielectric constant of the SL can be

where Δnrefers to the change of refractive index of sensing medium.[31]Since the dimension of the unit-cell is assumed to be subwavelength, the surface with the patch-array can be characterized by a homogeneous surface impedance.For TMwaves,the surface impedance of BP-patches can be expressed by[39]

A serial relationship exists between the IC layer and the semiinfinite SL,and the serial impedance is

and the symbol “‖” represents the parallel relationship betweenZBPandZzs.The reflective coefficient is expressed by

here the free space impedance for different incident angles is given asZ0=η0cosθ, whereη0=120πis the plane wave impedance in free space.The GH shift is then obtained by the well-known result[42]

whereλ0is the vacuum wavelength of the incident light.Obviously,the GH shift can reach a peak when the phase difference experiences a sharp variation with the incident angle.Finally,the phase (φ) and reflectivity (R) of reflection coefficient are obtained from the formulae

We define ΔGHas the change of GH shift with the changing of the refractive index and Δn=0.0002.Therefore,the sensitivity can be defined as[32]

Fig.1.(a) Schematic diagram of the spatial shift on the metasurface of BP–ZnS–SL.(b) Model-I is defined as the armchair (AM) direction along the x-axis and model-II means the zigzag(ZZ)along x-axis.(c)The equivalent circuit diagram of this structure.

3.Results and discussion

For subsequent theoretical and numerical calculations,both the BP-patch period (D) and the thickness of ZnS (d)are 500 nm, andLxandLyare fixed at 100 nm.Of course,the structure parameters of BP can be adjusted to flexibly control the GH shift.The ZnS crystal is selected as an example wherefT=8.22 THz,εh=5.0 andεl=8.3.[43]The damping of ZnS is supposed to beτd=3 cm?1.Figure 2 presents the permittivity of ZnS versus the frequency.In the region offT<f <1.2884fTsituated in the mid-infrared range, the ZnS crystal possesses a negative permittivity which leads to the hyperbolic properties.For the conductivity of BP,a particular set of parameters for model-I is used,namely,mx=0.2m0andmy=0.7m0(m0is the static electron mass),η=0.01 eVωx=1.0 eV,sx=1.7,ωy=0.35 eV, andsy=3.7.[44]For model-II,we only need to exchange the electron mass alongxandy-directions.The real and imaginary parts of BP conductivity as a function of the frequency under the differentnBare illustrated in Fig.3.Atf=38.48 THz,a transition takes place for the real part ofσACand a dip is observed for the imaginary part.It is well-known that the BP exhibits the strong structure anisotropy, which means that the puckered lattice results in two inequivalent directions in-plane.As a result,the transition point for the real part ofσZZoccurs at about 13.74 THz,which is close to the RB of ZnS.In addition,the values of conductivity increase with the increasingnB.The electromagnetic response of the anisotropic 2D material is normally classified into two district regions.In the lower frequency region, the conductivity is considered as a pure Drude type, where the in-plane anisotropy results from the effective mass of the electrons traveling in various directions.Yet,the contribution from interband electron transitions may become dominant at higher frequencies.In the following discussions,we mainly focus on the two typical models.

Fig.2.(a)Real and(b)imaginary parts of permittivity of ZnS slab versus frequency.The colored region indicates the RB.

Fig.3.Conductivity of BP as a function of the frequency under the different nB and the transition occurs at(a) f =38.48 THz for σAC and(b)f =13.48 THz for σZZ.

The relative impedance and the GH shift versus the frequency atN=1 andnB=1.0×1013cm?2are presented in Fig.4.For model-I shown in Fig.4(a),the relative impedance amplitude undergoes two transitions located in the lower and higher frequency regions,respectively.The first transition occurs at 11.2 THz,which is just close to the right boundary of the ZnS’s RB(see Fig.2).The other transition near 38.48 THz should be caused by the BP layer if we check its conductivity.An obvious GH shift is observed near the first transition atf=11.12 THz where the impedance approaches to zero,as shown in Fig.4(b).As a result of the zero reflection (see Eq.(8)),the incident angle will be close to the Brewster angle.When the BP’s zigzag direction is rotated to the direction of the armchair, only a transition in the lower frequency region is preserved and shifts from the higher frequency to the lower frequency, as illustrated in Fig.4(c).The inset in Fig.4(c)also indicates that the significant GH shift normally is excited near Brewster angles,as shown in Fig.4(d).It is reasonable to predict that the impedance transition is the result of the joint action of BP and ZnS.The in-plane anisotropy of BP will play a crucial role in the generation of large GH shifts.

Fig.4.Impedance of the BP–ZnS–SL structure and the GH shift versus the frequency at N=1,nB=1.0×1013 cm?2,Lx=Ly=100 nm and(a)θ =60.89° for model-I,(b)θ =60.91° for model-II.

Fig.5.Variation of (a) GH shift, (b) reflectivity and (c) phase with respect to angle of incidence under different nB at f =13.74 THz, N =1 and Lx=Ly=100 nm for model-I.

Next,we investigate the regulation of BP on the GH shift by varying the carrier density, layer number, and size of BP.At first,we consider the effect ofnBon the GH shift since it is easily controlled by changing the gate voltage.For model-I, the GH shift, reflectivity, and phase with respect toθat 13.74 THz are presented in Fig.5.With the increasing ofnB, the GH shift turns from the positive to the negative one,and a large GH shift about Δxmax≈?378.98λ0is obtained atnB=3×1013cm?2.The lowest points of reflectivity are closer to zero and the phase becomes sharper, indicating that the incident angle is located near the Brewster angle, i.e.,θ=60.89°.AsnBincreases further, the GH shift decreases rapidly and maintains a relatively stable value about?73.7λ0whennB≥7×1013cm?2.On the other hand, the GH shifts atf=11.2 THz and 38.48 THz(see Fig.4(b))are very small,while it has been proven that the GH shifts are independent ofnB.For model-II,figure 6 illustrates the variation of GH shift,reflectivity, and phase withnBat the two special frequencies(see Fig.4(c)),respectively.For both cases,the maximum values of GH shift are obtained atnB=1×1013cm?2and rapidly decrease with the increasing ofnB.In Fig.6(a), the corresponding Brewster angle is fixed atθ=61.21°.However, in Fig.6(b),the GH shift exhibits the red shift phenomenon asnBincreases, resulting in the different Brewster angles for each curve.Moreover,the phase curve is steeper,implying that the conditions for obtaining a larger GH shift are relatively more stringent.These findings confirm the validity of regulating GH shift via gate voltage,which may provide an effective way to flexibly increase the GH shift, leading to the optimum sensitivity of the structure.

Fig.6.Variation of GH shift,reflectivity,and phase with respect to angle of incidence under different nB at(a) f =11.12 THz,(b) f =13.609 THz,N=1 and Lx=Ly=100 nm for model-II.

The effect of layer numberNof BP film on the GH shift is investigated further, and the maximum values of GH shift are screened out by precisely adjusting the size of BP while remainingnB=3×1013cm?2.We limitLxandLyto vary from 10 nm to 500 nm.Figure 7(a)depicts the maximum GH shift caused by a blue shift atf=11.2 THz ifN <5,resulting in different Brewster angles for each GH shift.The largest GH shift is found atN=3,which is about?7565.58λ0.It can also be seen that the size of BP required to obtain the maximum GH shift should satisfy the condition ofLx >Ly.In the inset of Fig.7(a),we present the distribution of the electric field at the maximum GH shift by using Comsol Multiphysics software.The localization effect is clearly visible on the edges of the BP-patches.Figure 7(b)demonstrates that large GH shift is found at the small size of BP-patch whenf=13.74 THz,such as,Δxmax=3150.76λ0atLx=11 nm,Ly=357 nm andN=2.Unlike Fig.7(a),the Brewster angle is normally fixed at 60.889°and the length of BP satisfiesLx <Ly.SinceLx ?Dthe interaction between the BP-patches obviously disappears,as shown in the inset of Fig.7(b).The stronger light localization only exists around each BP-patch.Finally,figure 7(c)plots the GH shift versus the incident angle atf=38.48 THz.It is apparent that whenN >2,the GH shift becomes negative and increases with the increase ofN.The largest GH shift can be obtained atN=5 and the BP-patch is close to the square.Although the interaction between the adjacent BP-patches is clearly enhanced compared to Fig.7(b), the GH shift is reduced due to the overall reduction of the optical localization.To summarize,if a large GH shift is pursued,then we should consider the optimum combination of BP-patch size and layer number given in Fig.7(a).However,if a small BP-patch size is expected,then we have to sacrifice the GH shift and consider the recommended dimensions in Fig.7(b).

Fig.7.Maximum of GH shift under the different size of BP and layer numbers at nB = 3.0×1013 cm?2, and (a) f = 11.12 THz, (b) f =13.609 THz,and(c)nB=1.0×1013 cm?2, f =38.48 THz for model-I.The insets show the distributions of electric fields for the maximum GH shifts.

For model-II,according to the discussion in Fig.6 we fixnB=1×1013cm?2and further adjust the size and layer number of BP to obtain the maximum of GH shift,as illustrated in Fig.8.Atf=11.12 THz in Fig.8(a),the GH shift presents a red shift phenomenon and changes from positive to negative.The maximum GH shifts achieve atN=1 whenLx=191 nm andLy=10 nm.The Brewster angle for each GH shift also drifts to the right asNincreases.In addition,Lx ?Lyif the maximum GH shift is desired.The strength and distribution of the electric field are presented in the insets.For the case atf=13.609 THz shown in Fig.8(b),the Brewster angle is fixed atθ ≈60.91°independent ofN.In addition,whenLx ≈LyandN=1 the GH shift is largest.However,onceN >1,the conditionLx >Lyhas to be selected to obtain the maximum value of GH shift.Compared with Fig.8(a),the maximum value of GH shift is relatively smaller,which may be understood if we carefully examine the strength of the electric field shown in the inset of Fig.8(b).Based on the above discussion,the large GH shift can be realized by artificially designing the size and the layer number of the BP-patch in the actual device fabrication.BP is a 2D anisotropic material,which can be deposited on the surface of ZnS by the chemical vapor deposition(CVD)technique[45]and molecular beam epitaxy(MBE).[46]

Fig.8.Maximum of GH shift under the different size of BP and layer numbers at nB =1×1013 cm?2, and (a) f =11.12 THz and (b)f =13.609 THz for model-II.The insets show the distributions of electric fields for the maximum GH shifts.

Fig.9.Highest sensibility with the largest GH shift for the different layer numbers of BP.The refractive index varies from 1.75 RIU to 1.85 RIU.The other parameters are the same as in Figs.7 and 8 for the different models, respectively.The inset in panel(c)shows the S is very small at N=1.

Finally, based on the controllable GH shift, the possibility of its application in sensors is discussed.We use a sensitivity factorSto characterize the performance of the sensor.Both the layer numberNand carrier densitynBof the BP film greatly determine the GH shift.Therefore, the control variable method is used to analyze the performance of the sensor based on GH shift under the differentNandnBat Brewster angles.The highest sensibility with largest GH shift (see Figs.7 and 8) is sought when the refractive index varies from 1.75 RIU to 1.85 RIU with a scanning interval of 0.0002.Thus, it will make the detection of GH shift more accurate and convenient, and the sensitivity of the sensor can be improved correspondingly.The condition to produce the highest sensitivities under the differentNis shown in Fig.9,and the other parameters,including the size of the BP film,the carrier density, and the frequency, are the same as Fig.7 for the model-I and Fig.8 for model-II, respectively.The variation of the refraction index for eachSis indicated at the top of the bars.The sensibilityScorresponding to the GH shift in Figs.7(a)and 8(a)is presented in Fig.9(a).The maximum value ofScan reach 3.95×107λ0/RIU atN=2 for model-I and 1.62×107λ0/RIU atN=1 for model-II.Moreover, the order of magnitude ofSremains 107λ0/RIU for both models,with only a little difference whenN >2.Figure 9(b)presents the largest sensibility related to the GH shift in Figs.7(b)and 8(b).It can be seen that the sensibility atN=2 can reach about 1.5×107λ0/RIU for model-I and 1.72×107λ0/RIU atN=4 for model-II.Compared with Fig.9(a), the sensibility is obviously smaller.In particular, the smallestSis about 7.78×106λ0/RIU atN=5 for model-II.Finally, we check the sensibility atf=38.48 THz only for model-I related to Fig.7(c).As shown in Fig.9(c), the largest sensibility will reach 8.36×106λ0/RIU atN=5, which is smaller than the largest ones obtaining for the other two frequencies by an order of magnitude.The inset gives the sensibilitySforN=1 although the GH shift is not shown in Fig.7(c).It is clearly seen that the order of magnitude of the correspondingSis 10?7λ0/RIU,which may imply that a small GH shift has the potential to accompany a small sensitivityS.To summarize,it can be deduced that not only large GH shifts but also the flexible adjustment of sensibility are achieved in the BP–ZnS–SL sensor.It is noteworthy that the carrier densitynBhas a great influence on GH shift, which can be easily controlled by the external field.Thus,we will consider the sensibility regulation ofnBwhen the proposed structure is applied to the sensor.

Based on the simulation results withN, we investigated the effect ofnBon the sensibility, as shown in Fig.10.The physical parameters are same as those in Fig.9,including the selection of the refraction index for each curve.In general,theSof model-I is higher than that of model-II,which can be the order of 108λ0/RIU;for example,S=6.43×108λ0/RIU atN=5,nB=2.98×1013cm?2andf=11.2 THz as well asS=5×108λ0/RIU atN=4,nB=0.997×1013cm?2,andf=38.48 THz.It can be seen from the above discussion that the precise control of the carrier density helps to obtain the higher sensibility in BP–ZnS–SL sensor.Compared with that in the symmetrical graphene-cladding waveguide and in a graphene–hBN heterostructure,[32]the sensibility obtained in the BP–ZnS–SL metasurface structure is enhanced at least two orders of magnitude,which suggests the potential application of sensors based on the GH shift in THz range.

Fig.10.Sensibility versus nB for the different layer numbers based on Fig.9.

4.Conclusion

In conclusion, we propose a strategy to achieve tunable and larger GH shifts with the high sensibility in THz range via the excitation of SPPPs based on the BP–ZnS–SL metasurface structure.Under differentN,the largest GH shifts are predicted by dynamically adjusting the BP-patch size and the electron doping level.A comparison was made for the GH shift between model-I and model-II near the RB of ZnS and the transition of BP’s conductivity.It has been demonstrated that the largest GH shift can reach about?7565.58λ0atf=11.2 THz andθ=61.095°whenN=3,nB=3×1013cm?2,Lx=308 nm andLy=100 nm for model-I.The localized electric field around BP-patches can enhance the GH shift significantly.As a SPPP sensor based on the changing GH shift in THz range, the sensibility can reach at 6.43×108λ0/RIU when the index refraction of the senor medium is in the range of 1.75 RIU–1.85 RIU and the scan interval is 0.0002.The sensibility is at least two orders of magnitude larger than that of the conventional SPR sensor, which is extremely sensitive to the parameters of BP,including the size,carrier density,and the layer number near Brewster angles.These findings can serve as a theoretical foundation for designing SPPP sensors with high sensibility accompanied by a large GH shift in THz frequency.

Acknowledgments

Project supported by the Natural Science Foundation of Heilongjiang Province of China(Grant No.LH2020A014)and the Graduate Students’Research Innovation Project of Harbin Normal University(Grant No.HSDSSCX2022-47).

主站蜘蛛池模板: 欧美日韩专区| 久久久亚洲色| 99热线精品大全在线观看| 丰满人妻中出白浆| 精品99在线观看| 玖玖免费视频在线观看| 亚洲成人一区二区三区| 亚亚洲乱码一二三四区| 中文字幕啪啪| 一区二区在线视频免费观看| 88av在线| 国产凹凸视频在线观看| 午夜无码一区二区三区在线app| 久久国产黑丝袜视频| 国产h视频在线观看视频| 91精品国产情侣高潮露脸| 欧美成人日韩| 国产人免费人成免费视频| 国产精品漂亮美女在线观看| 亚洲国产日韩视频观看| 91小视频在线| 91精品国产综合久久不国产大片| 夜色爽爽影院18禁妓女影院| 日本欧美视频在线观看| 久久精品中文字幕少妇| 毛片网站免费在线观看| 日本不卡在线视频| 91无码人妻精品一区| 国产精品视频999| 亚洲无码37.| 国产在线观看91精品| 色网站在线视频| 91视频免费观看网站| 亚洲黄色成人| 玖玖精品在线| 日本伊人色综合网| 欧美另类一区| 久久一本日韩精品中文字幕屁孩| 中文字幕va| 最新加勒比隔壁人妻| 久久综合丝袜日本网| 精品国产免费观看| 亚洲天堂网2014| 中文字幕无码电影| 亚洲女同一区二区| 亚洲成人动漫在线观看 | 国产精品一区不卡| 国产欧美在线观看一区| 亚洲视频免费在线看| 91小视频在线播放| 久久午夜夜伦鲁鲁片无码免费| 国产麻豆另类AV| 米奇精品一区二区三区| 日韩国产黄色网站| 狠狠色成人综合首页| 最新亚洲av女人的天堂| 91系列在线观看| 国产97视频在线| 一级毛片在线播放免费| 91青青草视频| 欧洲亚洲欧美国产日本高清| 激情视频综合网| a级毛片毛片免费观看久潮| 国产精品视频3p| 一级做a爰片久久毛片毛片| 九九热视频在线免费观看| 精品成人一区二区三区电影 | 青青国产视频| 亚洲va欧美va国产综合下载| 午夜a级毛片| 国产a v无码专区亚洲av| 亚洲有无码中文网| 国产精品一区不卡| 九色视频最新网址| 欧美三級片黃色三級片黃色1| 亚洲AV电影不卡在线观看| 国产av色站网站| 久草视频精品| av无码久久精品| 老司国产精品视频91| 亚洲男女天堂| 黄色污网站在线观看|