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

    Third-Order Optical Nonlinearity in Two-Dimensional Transition Metal Dichalcogenides?

    2018-09-10 06:39:50SinaKhorasani
    Communications in Theoretical Physics 2018年9期

    Sina Khorasani

    School of Electrical Engineering,Sharif University of Technology,Tehran,Iran

    Ecole Polytechnique F′ed′erale de Lausanne,CH-1015,Lausanne,Switzerland

    AbstractWe present a detailed calculation of the linear and nonlinear optical response of four types of monolayer twodimensional(2D)transition-metal dichalcogenides(TMDCs),having the formula MX2with M=Mo,W and X=S,Se.The calculations are based on 6-band tight-binding model of TMDCs,and then performing a semi-classical perturbation analysis of response functions.We numerically calculate the linearand nonlinear surface susceptibility tensorswith ωΣ = ωr+ωs+ωt.Both non-degenerate and degenerate cases are studied for thirdharmonic generation and nonlinear refractive index,respectively.Computational results obtained with no external fitting parameters are discussed regarding two recent reported experiments on MoS2,and thus we can confirm the extraordinarily strong optical nonlinearity of TMDCs.As a possible application,we demonstrate generation of a π/4?rotated squeezed state by means of nonlinear response of TMDCs,in a silica micro-disk resonator covered with the 2D material.Our proposed method will enable accurate calculations of nonlinear optical response,such as four-wave mixing and highharmonic generation in 2D materials and their heterostructures,thus enabling study of novel functionalities of 2D photonic integrated circuits.

    Key words:2D materials,nonlinear optics,transition metal dichalcogenides,quantum optics

    1 Introduction

    Following the discovery of the celebrated material,graphene in 2004,the field of two-dimensional(2D)materials has been rapidly expanding over the recent years.[1?9]Despite being only a monolayer thick,2D materials exhibit extraordinary electronic,mechanical,thermal,spintronic,and in particular,optical properties.These o ff er novel and unprecedented applications,which were not foreseen earlier such as new nonlinear capacitors[10]for using in cryogenic parametric amplifiers,circulators,and mixers,as well as a new type of capacitive qubit referred to as cubit[11]in quantum nonlinear superconducting circuits.Of primary importance,is the nonlinear optics of the 2D materials and structures made out of them,which is a matter of current deep investigations.

    A number of researches discuss the second-order nonlinear susceptibility as well as second-harmonic generation in various 2D materials.[12?16]Similarly,the third-order nonlinear optical properties of graphene has been studied in many works[17?25]and the two-dimensional nonlinear sheet susceptibility χ(3)2Dor the conductivity σ(3)2Dtensor of graphene have been calculated.However,recent experiments[26]still disagree with the expected numerical estimates from theory within an order of magnitude.The nonlinearity of graphene is remarkably large,however,it needs to be gated to adjust its Fermi level to the resonance conditions.That implies for ungated structures,there is still a tendency to examine other 2D materials as well.Similar discrepancies between measured sheet optical nonlinearity of alternative 2D materials and theoretical estimates are noticed by other researchers as well.

    In this work,we present a detailed study of the nonlinear optical properties of 2D transition metal dichalcogenides(TMDCs),[27]with the general formula MX2with M being a transition metal,here being either Molybdenum Mo or Tungsten W,and X being a chalcogen such as Sulphur S or Selenium Se.Based on a six-band Tight-Binding Hamiltonian,we obtain the nonlinear sheet susceptibility χ(3)2Dof the TMDCs through semi-classical perturbation approach.

    The semi-classical method retains its validity for low illumination intensities,where Zener tunneling and semimetal transitions could be well ignored and dismissed from the carrier dynamics,[28?30]leaving only multiphoton processes as important.This approximation also is limited by the finite span of optical wavelength and nonzero extent of atomic bonds,that is,the medium has to be effectively treated as continuous and its microscopic granular structure shall be ignored.For ultrashort wavelength optical excitations beyond the deep ultraviolet spectrum,this treatment is no longer valid.

    Such nonlinear interactions[31]are of primary importance in study of valley-spin dynamics[32]and Kerr spectroscopy[33]of TMDCs.In all these last works,the TMDC monolayer has been sandwiched in protecting 2D Boron Nitride shields,which has greatly contributed to the visibility and sharpness of emission spectra.[31?34]

    The reported data on χ(3)2Dof MX2are very scattered up to four orders of magnitude,as it is shown to be strongly dependent on the growth method and the posttreatment process.This makes a conclusive evaluation out of reach at present.However,we notice agreement roughly within an order of magnitude between the computed numerical figure and the stronger one of the experimental values.

    We study the nonlinear optical properties of an ultralow-loss amorphous silica microdisk covered by WSe2,and define an effective nonlinear susceptibilityfor the structure.We show through detailed COMSOL calculations that the effective nonlinearity marked byis improved up to two orders of magnitude by placing the TMDC on top of the silica microdisk.This platform has been shown to be useful for experimental study of optical emission and excitonic properties of TMDCs.[35]

    As an application example,we consider the squeezing property of light[36]through four-wave mixing under unpumped and pumped configurations and discuss these scenarios.While calculations reveal the dominance of twophoton absorption over nonlinear Kerr self-phase modulation in 2D TMDCs,we show through detailed theory that an unorthodox π/4-rotated squeezed state of light is produced with an elongated elliptical Wigner distribution and theoretically unlimited squeezing.Without consideration of this π/4-rotation,both quadratures appear to be desqueezed.However,unlimited squeezing is possible only along one π/4-rotated quadrature and the orthogonal quadrature will always be strongly desqueezed.

    2 Theory&Results

    2.1 Band Structure

    In order to compute the linear and nonlinear susceptibility from first principles,we would need to have accurate knowledge of the electronic transitions,valleys,and spinorbit interactions.The correct way to tackle this problem is to have an efficient code to derive the electronic band structure of the material,and since this information is going to be called upon quite frequently inside integration and summation loops,the computation has to be both accurate and very efficient.For this purpose,tight-binding(TB)scheme is very appealing since with correct implementation it could meet both criteria.

    The method TB is quite popular for low-dimensional carbon structures such as graphene and carbon nanotubes,[37?38]and six-band TB has been used for studying the band structure of hydrogenated graphane.[39?40]In recent years,two-band[41?43]expansion based on k ·p perturbation and L?wdin partitioning of Hamiltonian,[44]three-band,[45]four-band,[46]six-band,[47?50]sevenband,[42]and eventually eleven-band[51]TB models have been developed to calculate the band structure of TMDCs.However,six-band TB based[50]on Slater-Koster expressions[52]is ultimately chosen since it is a fairly accurate low-energy approximation to the extensive density functional theory(DFT)calculations,and thus expected to be both sufficiently accurate and efficient for our purpose.Moreover,the 6-band TB method,considers the effects of valley-dependent spin-orbit interaction of 2D TMDCs,[53]or better known as trigonal warping.[41?42]This method employs an orthonormal basis,and therefore does need inversion of the overlap matrix,[37,39]which eases out coding and efficiency at the same time.Since the details of this scheme is rather comprehensive,the reader is referred to the existing literature for further information.[47?50]However,in Appendix A we present the necessary ingredients brie fly.The orthonormal bases here are

    In general for honeycomb lattice with C3vsymmetry such as graphene[38]and in absence of chirality,it is normally sufficient to compute the band structure over the irreducible Brillouin zone,which is only 1/12 of the first Brillouin zone.However,monolayer TMDCs are noncentrosymmetric and satisfy a different spatial group denoted as,[54]which causes spin-valley coupling.As a result,the irreducible zone is only 1/6 of the first Brillouin zone.

    The calculated band structures are shown in Fig.1 for the four basic types of TMDCs MX2with M=Mo,W and X=S,Se.Here,red and blue curves correspond to respectively down and up spin polarizations.Interestingly,the calculations preserve valley-dependent spin-orbit interaction or the trigonal warping of TMDCs,which is strongly present in the valence bands at K and K′.To illustrate this,we calculate the entire first Brillouin zone of WSe2,in which spin-orbit interaction is highly different at K and K′.As a result of this symmetry breaking,the band structure shows a C3vpoint-group symmetry,instead of the expected C6v,and by changing the direction of spin the band structure rotates by π/3,thereby interchanging the roles of K and K′.The conduction and valence bands of WSe2are illustrated in Fig.2.

    Fig.1 (Color online)Band structures of four basic TMDCs,based on the 6-band TB model.[50]Red and blue lines correspond to the spin down and up bands.

    Fig.2 Illustration of valley-spin coupling in WSe2as an example,based on the 6-band TB model.[50]The effect of trigonal warping is seen to be strongly pronounced for the valence band,where K and K′behave very differently.

    2.2 Linear Susceptibility

    The concept of an interface occupied by a dielectric having finite surface electric and magnetic susceptibilities has been originally investigated in a pair of papers published in 1994.[55?56]Independently and still a few years before celebrated discovery of the first 2D material,graphene,the optical properties,possible modulation and switching applications of such media,referred to as the conducting interfaces were explored by a series of papers by the author.[57?62]The transfer matrix method has been reformulated so far independently by many authors to tackle the problem of optical wave refraction from layered structures containing 2D materials[56,58,63?65]and even more recently by Morano[66?67]for measurement of the surface conductivity of graphene.

    Based on the theory of conducting interfaces,an ultrathin 2D material could be treated by a discontinuity in tangential electromagnetic fields across the interface.However,the corresponding surface susceptibility or χ2Dwith the dimension of meter,or equivalently,surface conductivity σs=j?0ωχ2Dwith the dimension of ??1should be known.

    Here,we are able to compute the surface susceptibility tensor elements χ2D(ω)= χ(1)2D(?ω;ω)of TMDCs using the expression

    where

    is the Fermi-Dirac distribution with EFbeing the Fermi energy(here equal to zero for the intrinsicsemiconductor),and kBT being the thermal energy.The matrix elements(k)are related to the current operator(k)as

    in which kμare elements of the vector k.Moreover,q is the electronic charge and H(k;σ)is the spin-dependent Hamiltonian in k-space,which in our TB formalism appears as a 6×6 matrix.More on the details of this operators could be found in Ref.[68].In Eq.(3),?0is the permittivity of vacuum,AUCis the area of unit cell,and ABZ=(2π)2/AUCis the area of first Brillouin zone.[69?71]That implies the expression NAAUCis the area occupied by one mole of the 2D material,where NAis Avogadro’s constant.[69?75]One can easily verify thatwhere a is the 2D crystal’s lattice constant,being equal to the distance between the two nearest metal atoms.

    In the third-order nonlinear processes,the location of Fermi-energy has very little observable effect on the Kerr coefficient,and much less on the third-harmonic generation coefficient.This is a fact observed through extensive numerical experiments.Only if the Fermi level is deeply into the conduction or valence bands,that is,the semiconductor is made degenerate,the result would become different.Normally,since chemical doping is not in practice for TMDCs,any shift in the Fermi level would be triggered by electrostatic gating.This method of charge depletion or accumulation,however,is typically unpractical for making a 2D TMDC degenerate.The author believes that there is no practical reason to be concerned about the effect of extrinsic Fermi level as opposed to the intrinsic case with EF=0.

    When the TMDC is not under strain,then linear susceptibility χ(1)2Dμν(?ω;ω)=χ(1)2D(ω)δμνis a simple scalar and not a tensor quantity.This fact together with the spin-valley coupling of carriers can be used to simplify(3)a bit as

    Here,the integration on the reciprocal plane is taken only over the irreducible Brillouin zone,and summation on valley index τ= ± is for selection between K and K′valleys.Since valley contributions from opposite spins are simply equal,one may take advantage of symmetry considerations and drop the summation on τ,just multiplying the whole expression on the right by a factor of 2,thus arriving at

    Results of computations for the linear surface susceptibilities of the four basic TMDCs in the wavelength range(760–790)nm is illustrated in Fig.3.Although the absolute values are here shown,imaginary parts of χ(1)2D(?ω;ω)in the desired range is effectively zero.

    Fig.3 Typical values of first-order linear susceptibility χ(1)(?ω;ω)for different TMDCs.Contribution to the absolute value is entirely from the real part.

    2.3 Nonlinear Susceptibility

    Folllowing the standard perturbation scheme,[76?88]we have developed a code in Mathematica to compute the non-degenerate four-wave mixing(FWM)third-order nonlinear susceptibilitywhere ωΣ=ωr+ ωs+ ωtfor any of the widely used TMDCs.When studying 3rd harmonic generation,must be calculated,while for degenerate nonlinear propagation,should be found.Expressions are given in the Appendix B.It is worth mentioning that the method of equations of motions[68,89?90]could also be used to tackle the dynamics of nonlinear optics and fourwave mixing in semiconducting materials,however,the formalism does not directly give in an expression for the nonlinear susceptibility.

    Figure 4 represents typical calculated values from semi-classical perturbation theory of nonlinear interac-tions.Similarly,Fig.5 presents the typical nonlinear coefficient in the telecommunication window around the wavelength 1550 nm.The uniform convergence of the code has been demonstrated in Fig.6 versus computational grid.These set of figures present absolute value of the tensor componentcorresponding to the third-harmonic generation.Calculation of other non-zero tensor components has also been done,but nor presented for the sake of brevity.Further discussions on the relationship among non-zero tensor components is given in the following.

    Fig.4 Typical values of third-order nonlinear susceptibility χ(3)(?3ω;ω,ω,ω)for different TMDCs.Contribution to the absolute value is entirely from the real part.

    Fig.5 The calculated third-order nonlinear susceptibility χ(3)(?3ω;ω,ω,ω)for MoS2in the infrared telecommunication window.Contribution to the absolute value is entirely from the real part.

    Fig.6 Typical convergence of the code at the wavelength of 1560 nm versus resolution of the computational grid,corresponding to Fig.5.

    Computations on the Kerr nonlinear susceptibilityreveals that it is entirely imaginary in the wavelength range of interest,that is,it is the twophoton absorption,which is the dominant Kerr nonlinear effect.Furthermore,this coefficient turns out to be a scalar like χ(3)(?ω;ω,?ω,ω),as is discussed later below.Typical values for different types of TMDCs are illustrated in Fig.7.Comparing to Fig.4,much less oscillations are seen,and that is because of the much less number of available resonances with transitions among band edges.For the case of MoS2,in particular,many resonances occur as shown in Fig.7,which can be attributed to the bandgap of this material being the smallest of the four studied TMDCs.As a result,for a given photon energy above the energy gap,much more possible triple transitions required for a third-order nonlinear process become available.Evidently,selection of a photon energy sufficiently small,that is the case in Fig.5,no such behavior is seen,simply because of the absence in availability of corresponding possible transitions.

    Fig.7 Typical values of Kerr nonlinear susceptibility χ(3)(?ω;ω,?ω,ω)for different TMDCs.Contribution to the absolute value is entirely from the imaginary part.

    Furthermore,the magnitude of Kerr nonlinear index is significantly larger than the third-harmonic generation.It is thus to be noticed well that these two parameters of Kerr nonlinearity and third-harmonic generation,although related,but are essentially much different in both their physics and behavior.Actually,these two values could be quite far apart in phase and many orders in magnitude.

    We here point out again that the Kerr nonlinear susceptibility χ(3)(?ω;ω,?ω,ω)and Third-harmonic generation susceptibility χ(3)(?3ω;ω,ω,ω)are not the same.While both are having the third-order,these correspond to entirely different processes.This difference can be understood as follows.

    The Kerr nonlinear susceptibility only describes a process in which the local permittivity of the medium is linearly dependent on the local intensity of light,such aswhereis the phasor of the electromagnetic field.Hence,this effect can for instance cause local increase or decrease of the effective refractive index depending on the sign of ?[χ(3)],which correspondingly lead to self-focusing or self-defocusing phenomena.While higher harmonics inevitably may weakly arise under such conditions,because of the inherent nonlinearity in the propagation,it is only the first harmonic,which undergoes self- field propagation issues.More accurately,

    However,the third-order nonlinear susceptibility directly is connected to the local amplitude of the harmonic,as a result of locally present electromagnetic field.Hence,one could imagine that the third-harmonic be related to the first harmonic as

    Since the physical phenomena behind these two mechanism are not identical,one would naturally expect thatSometimes,this fact is not sufficiently made clear in the literature.

    (i)Critical Field

    De fining a critical electric field asas a measure of required lighs electrical field to reach the onset of third-order nonlinearity,we observe that for the TMDCs this value falls within the typical range of(7–8)kV/m,regardless of the choice of the particular material. For all the studied cases in the wavelength range of(750-790)nm,the approximations,as well asapply well to the linear and nonlinear surface polarizabilities.We use the conducting interface formulation[57?58]to carry out any electromagnetic analysis of monolayer 2D materials.

    (ii)Analogy with Bulk Values

    In this wavelength range of interest and for the linear response,the typical strength of linear susceptibility χ(1)is within the range of(3.5–4.4) ×10?8m,so that having an effective monolayer thickness of t2Dwe may assign an effective refractive index ofto the ultrathin TMDC layer.Since,t2Dis only of the order of typically(7–10) ?,then effective refractive index is expected to be about only n ≈ 6–8 in the wavelength range presented here.However,we notice that at longer infrared wavelengths,this value dramatically increases to higher values.

    Referring to Fig.7,TMDCs have an extremely large surface nonlinearity χ(3)2Din the range of 5 ×10?24m3·V?2to 2×10?23m3·V?2,[76,91]which when divided by the typical thickness of a monolayer of the order of 0.5 nm,give rise to an effective nonlinearity of the order of 10?14m2·V?2to 10?13m2·V?2.Given that this value is resulting from pure electronic polarization,this is quite remarkably large when put into perspective of expected values of the order of 10?22m2·V?2.[77]

    This places the expected nonlinearity of 2D TMDCs,to many orders of magnitude stronger the range of IIIV semiconductors such as GaAs 1.4 × 10?22m2·V?2,and than that of Diamond 2.5 × 10?21m2·V?2,fused silica 2.5 × 10?22m2·V?2,and even GaP at 577 nm 2.93 × 10?18m2·V?2,[80]which is known to have an extremely large nonlinear index of refraction.

    We first define a bulk-equivalent susceptibility aswhere t2Dis the physical thickness of the monolayer. For instance then MoS2at 577 nm and 1560 nm respectively exhibits an expected bulk-equivalent,absolute nonlinear Kerr susceptibilityof 2.11 × 10?13m2·V?2and 1.17 ×10?12m2·V?2.In a similar manner,the corresponding expected bulk-equivalent,absolute third-harmonic susceptibilitiesare 2.96 × 10?17m2·V?2and 2.11×10?15m2·V?2at 577 nm and 1560 nm respectively.

    These numbers could reasonably well explain the recently observed ultrastrong high harmonic generation in 2D TMDCs.[76,91?92]There are certain classes of materials or media,which could o ff er significantly stronger nonlinearity,however,either their slow response times(such as polymers)or complexity of formation(such as cold atomic gases),render them of limited use at optical frequencies.

    The unusually large nonlinearity of 2D materials compared to bulk 3D structures is hard to explain.However,we believe that it is a matter of geometrical confinement dimensions that sets this strength.At least it has been rigorously established that the third order nonlinear optical response of spin density wave insulators is much stronger in 2D than 3D and 1D structures.[93]

    (iii)Analogy with Experimental Values

    The validity of numerical results could be furthermore verified against three recent experimental data[76,91,94]on 3rd harmonic generation from MoS2at the wavelength of 1560 nm.While our developed code estimates a value of 1.6×10?25m3·V?2for,as shown in Fig.7,the reported experimental data are 3.9×10?15m2·V?2×0.75 nm=2.93×10?24m3·V?2,[76]as well as the very different values of 1.7 × 10?28m3·V?2,[91]as well as 1.2×10?19m2·V?2×0.65 nm=7.8×10?29m3·V?2as the pre-treatment and 1.95×10?28m3·V?2as the posttreatment values,[94]and also 6.5×10?29m3·V?2.[95]

    This shows that while experimental results[76,91,94]vary within four orders of magnitude,our numerical estimate is much closer to the first measurement[76]reporting the larger value.The difference between theoretical results with experimental ones,therefore,should not be surprising considering the remarkably scattered numerical values among experimental observations.While the nature of such differences is not exactly known yet,it could be due to material growth and transfer issues,defect concentrations,substrate effects,as well as the optical method of measurement.Given the large sensitivity of TMDCs to the fabrication process and even resilience under ambient conditions,there remains a question of how to unify experiments alike.It could be speculated that for a variety of reasons,the experimental reports actually either underestimate or overestimate the nonlinear susceptibility over the true theoretical value.

    (iv)Behavior of Tensor Components

    With regard to the individual tensor components recognized by the set of indicesμνζη,there should be 24=16 elements.However,half of them are zero and the only nonzero tensor elements are yyyy,xxxx,xyxy,yxyx,xyyx,yxxy,xxyy,and yyxx.Even though,all tensor elements are not quite independent and the following identities hold for all four TMDCs due to the crystal symmetries

    This limits the maximum number of independent tensor elements to three.

    For the case of Kerr nonlinearity where the parameteris concerned,one may verify that the elements yyxx and yxyx also identically vanish.Hence,we haveand get a fairly simple scalar form for the nonlinear Kerr susceptibility as

    Similarly,we could write for the third-harmonic generation the following

    3 Non-classical State of Light

    In this section,we discuss production of a non-classical state of light with elliptical Wigner distribution,due to two-photon absorption.Since the method of analysis is based on the theory of squeezing,we need to initially present an overview of light squeezing schemes.

    3.1 Methods of Squeezing

    Squeezing of light is usually done through either of the following general methods:[96?102]

    (i)Squeezed light by parametric down conversion.

    (ii)Squeezed light in optical fibers.

    (iii)Squeezed light in atomic ensembles.

    (iv)Squeezed light in semiconductor lasers.

    Out of the above four methods,the first three are mostly implemented in a non-monolithic experiments and normally require large optical setups.The second one is based on the χ(3)effect of fused silica in optical fibers,and thus requires long propagation paths over fibers to achieve squeezing,but squeezing as large as 9dB has been achieved using this scheme.The third one requires sophisticated techniques of atom vapor trapping and condensation at ultralow temperatures,and the largest observed squeezing of 13dB has been so far achieved this way.The fourth method[103]is relatively easy to achieve and based on monolithic integrated photonics,but is limited for various practical reasons to only 3-4dB of reduction in shot noise and thus squeezing.

    Recent advances in optomechanics,has brought the possibility of optomechanical squeezing of light into perspective as well.[104?109]Optomechanical squeezing of light in homodyne detection within a small amount also occurs,and can be observed using a novel quantum feedback control scheme,which has been recently reported.[110]Monolithic approaches to squeezing by optical parametric oscillators[111]have been shown to be feasible as well.

    As an alternative route,the possibility of using Si3N4micro-ring resonators for squeezing has been demonstrated in a multi-mode optical parametric oscillator[112]with 1.7 dB squeezing.This has also apparently been verified experimentally at room temperature by pumping a continuous laser,[113]and 0.5 dB squeezing below shot noise was observed in a self-homodyne setup.So,it could be safely claimed that using silica microdisks covered with 2D TMDCs,certainly measurable squeezing could be obtained.It has been furthermore recently shown that high-quality silica micro-ring resonators could provide a versatile platform for study of emission properties of 2D materials.[35]

    3.2 χ(3)Squeezing on Micro-Resonators

    It is known that the 3rd order nonlinearity could be employed to generate squeeze light. Depending on the optical setup,that whether the squeezed state is produced through unitary transformation in vacuum,or an interferometric setup,it leads to either of the Hamiltonians[114?115]as

    The first squeezed state(10)being referred to as the quadrature squeezing,has an elliptical Wigner distribution,while the second one(11)being referred to as the photon-number squeezing results in a kidney-or crescentshaped squeezing and is thus quite different and non-Hermitian when ?[ξ] ≠0.Here, ω is the frequency of light and ξ is a dimensionless parameter defined as[115]

    where u(r)is the vector mode pro file of a single photon in the cavity with components um(r),and a factor 3/2 is already included for standing wave,as it is going to be the case under study.Furthermore,since the optical wave is strongly confined in the micro-disk and undergoes a much larger effective propagation path,an extra dimensionless finesse factor F is also included.The finesse F is a direct measure of how many times the light pulse circulates the ring.[116?117]This is while the single photon mode u(r)having the physical dimension m?3/2must be normalized as

    Here,χ(1)(r)= ?r(r)?1 is the relative susceptibility of the micro-disk resonator,including the cladding and substrate.Hence,Eq.(13)is effectively taken on the mode volume of the cavity.Obviously,Eq.(12)could be written for an unnormalized mode as

    In presence of a 2D material,and using a conducting interface approximation,[57]these two latter expressions should be changed slightly as follows

    where it is supposed that the 2D TMDC is placed on the z=0 plane.Following Eq.(8)for Kerr nonlinearity and two-photon absorption in unstrained TMDCs,both χ(3)and χ(1)are scalar quantities.This greatly simplifies Eqs.(15)and(16)as

    (i)Effective Nonlinearity and Mode Volume

    Alternatively,Eq.(17)can be written as

    where V is mode volume[118]andis the effective nonlinear index defined as

    It should be mentioned here that ξ and thereforeby definition cannot be independent of the micro-disk radius r as well as wavelength λ.

    Another way to look into this is to view the circulating optical power in disk resonator as an optical field going through a long straight path,[119?120]thus experiencing an overall phase retardation.This is the basis of light squeezing in optical fibers,[121?124]and this point of view gives a more clear and straightforward measure for the squeezing parameter s as

    where r and n are disk’s radius and index of refraction,respectively,λ is the optical wavelength,and E(r)is the electric field inside the disk.Moreover,N is the number of photons inside the disk,and thus Nω is the total optical energy confined in the cavity.

    Now,plugging in Eq.(17)and further simplification gives the final expression as

    (ii)Noise Squeezing/Desqueezing

    As it was discussed in the above,χ(3)is purely imaginary in the wavelength range of interest between 760 nm to 790 nm,and thus it is the two-photon absorption,which wins over the Kerr nonlinear index.Despite this fact,this is unimportant for our particular application of producing a non-classical state of light with non-circular Wigner distribution.

    Some studies have only considered the absolute value of|χ(3)|to the squeeze parameter s=|ζ|,[121]where ζ is the complex squeeze parameter given by the squeeze operator as[97,121,125?127]

    Ignoring the real part of ζ and noting that its imaginary part is positive,we may rewrite the above as

    This is actually directly corresponding to the interaction Hamiltonian(10)in the above,noting that s ∝ |ξ|.As it will be discussed later below,the fact that χ(3)is imaginary does not disallow production of non-classical states.

    It should be mentioned that actual ratio of noise squeezing for the case of Eq.(10),here being denoted respectively by ? and φ for the two orthogonal quadratures is not the same as s,since ζ is not real valued.As we can actually show in Appendix C,when ?[ζ]=0 the correct relationship is

    As it has been demonstrated in the Appendix C,the measure of non-classicality of the resulting Wigner distribution is limited to 6.02 dB at high intensities or high

    3.3 Experimental Issues

    One therefore may produce a non-classical light using a micro-disk resonator covered with a 2D TMDC,such as WSe2,or MoS2.We would like to investigate this possibility by taking advantage of the very large χ(3)coefficient of 2D TMDCs.Combined with the very strong confinement of light in low loss silica disk resonators,we would expect a significant anharmonicity ξ as defined in Eq.(15).This has yet to be calculated.If ξ is hopefully found to be reasonably large,and leading to an observable value of desqueezing asymmetry in light,we would move ahead with experiments to characterize the properties of output emission and shot noise.The enhanced effective nonlinear susceptibilityshould in principle allow production of non-classical elongated state at much lower intensities.

    (i)Ultrashort Solitons

    A possible advantage of this scheme in case of successful design and experiments,would be relatively ease of fabrication and operation,as well as compatibility with monolithic integrated photonics. Usage of ultrashort solitons[128?129]together with strong confinement in micro-disks increases the overall nonlinear interactions and thereby squeezing as well.This has been already demonstrated in squeezing via optical fibers.[96]In that case,no extra pumping is needed,and the pulse undergoes squeezing by itself through self-pumping.Evidently,the squeezed output is pulsed.Since the quality factors of cavities is not in finite,a small dissipation at high optical power densities is unavoidable,[118,130?133]which can be treated anyhow as either a Lugiato-Lefever equation[134]or using a perturbative expansion.[135]

    (ii)Separate Pumping

    It is in practice bene ficial if the optical power required for excitation of nonlinearity is maintained by a pump held at a slightly different frequency such as ωp= ω ± FSR,where FSR is the free spectral range of the micro-disk resonator.The susceptibility then should be calculated from the non-resonant near-degenerate pumped value χ(3)(?ω;ωp,?ωp,ω). Since FSR ? ω,then χ(3)(?ω;ωp,?ωp,ω)stays within the same order of magnitude of χ(3)(?ω;ω,?ω,ω),thus leaving the presented analysis and discussions basically intact.A closely related scheme with two pumps with frequencies=ω±?? where?? is a non-zero integer multiple of FSR has been proposed in fiber[136?137]and recently used in integrated[130]ring resonators for squeezing and random number generation,respectively.It is worthwhile to mention that the reverse of this scheme where two identical pump photons decompose into two different signal and idler photons in higher and lower frequency neighboring sidebands has been shown[112]to be useful for on-chip optical squeezing,too.

    A numerical evaluation of near-degenerate χ(3)(?ω;ωp,?ωp,ω),as shown in Fig.8 for an FSR of 3 nm at the wavelength of 775 nm,for instance,in case of WSe2reveals that χ(3)is actually reduced from the degenerate value of|χ(3)|=1.2 × 10?23m3/V2roughly by 49.3%,to the new non-degenerate pumped value of|χ(3)|=6.05 × 10?24m3/V2at a pump wavelength of 778 nm,which is again purely imaginary.This drop could be easily accounted for by a proportional increase in pump power.Similarly,if pumping is done at 772 nm,then|χ(3)|=5.95 × 10?24m3/V2,which is a 50.2%reduction again.Therefore,when ?λ = λp? λ with λpbeing the pump wavelength is sufficiently small,then the approximation

    holds.This type of behavior is more or less the same for all other sorts of TMDCs,as shown in Fig.8.

    This scheme,as opposed to the ultrashort solitons,enables continuous pumping and signal feeding,and therefore a continuous squeezed output may be expected as well.

    Fig.8 Computation of near-dengenerate susceptibility χ(3)(?ω;ωp,?ωp,ω)versus pump wavelength difference?λ=λp?λ at fixed λ=775 nm for various TMDCs.

    4 Excitonic Effects

    Without doubt,the prominent role of excitons in light emission from TMDCs could be considered as a subject of deep study.This is due to the fact that all of the four basic types of TMDCs discussed in this paper support both of the dark and bright excitons,which is further complicated by the trignoal warping property of these materials and presence of charged excitons(also known as trions)and biexcitons.[35,138?139]

    In general,exciton binding energies in 2D materials cannot be explained by a simple hydrogenic model,because of the very different radial distribution of the wavefunction,which could easily be extend to a few nm in radius.It only can be investigated by DFT GW-BSE calculations,which exhibits numerous lines between 1s and 2p states.While the exciton binding energies in 2D materials can be large,these are directly dependent on the substrate screening effects.When there are insulating encapsulation or separation,for instance,using bilayer graphene(or BN in other works),the exciton binding energy may drastically reduce by a factor of 3 to 6.

    At room temperature,exciton emission peaks are significantly broadened,because of large electron-phonon coupling rates in typical TMDCs.Biexcitons in 2D materials are unstable at room temperature and dissociate,since their binding energies are on the order of only a few meV.The exciton-exciton annihilation(EEA),which is a four-particle process can occur under high exciton population(demanding high illumination intensities)in principle even at the room temperature,but it can be expected that at such high rates,non-radiative recombinations through defect states and impurities would be dominant over the entire EEA process.Dark excitons in most 2D materials can actually emit light.Since they are optically forbidden because of vanishing dipole.However,the true dipole should be complex valued instead of a real one,so that they can emit non-linearly polarized light(not necessarily exactly circular).The fact that the energy different between dark and bright excitons is on the same order of exciton binding energy complicates the correct interpretation of light emission.

    While the existence of biexcitons at room temperature is unlikely because of the small binding energy,there are strong reasons to believe that emission from charged excitons could be easily mistaken with defect emissions because of interface effects.[32,140?141]Quite possibly,the existence of higher-wavelength emission peaks in photoluminescence spectra of TMDCs could be due completely due to interface defects,and in that case trions/biexciton emissions can have no in fluence at room-temperature as they normally appear only at cryogenic conditions.[142]If ture,then the observed trion emission peaks could be actually a defect-related emission,and for this reason is entirely absent in boron-nitride(BN)encapsulated monolayer TMDCs.[32,142]Similarly,there exists experimental evidence for brightened dark excitons/trions in WSe2,where measurements are all done at low temperatures.[143]

    Another recent study[144]is citing the fact that the two emission peaks for MoS2on thermally grown SiO2,correspond to A-exciton and interface Defect.The defect state disappears at room temperature,which is somehow connected to observations reported in Ref.[32].Furthermore,they propose an elegant way to identify a defect state.

    The emission of a defect is almost always unpolarized regardless of the pump incidence angle,intensity,and polarization,and that is why it cannot form polaritons in cavity quantum electrodynamics(CQED)experiments.A careful polarization measurement on the emission spectra could unanimously reveal whether the lower energy peak is a defect or excitonic emission.

    In summary,the origin and nature of peaks in the emission spectra of TMDCs has been a matter of unresolved debate[32,140?157]and still remains with no conclusive agreement so far.

    5 Conclusions

    We have presented a detailed theoretical analysis of linear and third-order nonlinear optical response of twodimensional(2D)monolayer transition metal dichalgonides(TMDCs). Based on a rigorous six-band tightbinding model,we have calculated first order and third order susceptibility tensors,and observed reasonable consistency with experimental data.We predict an elongated non-classical light using high quality silica micro-disk resonators covered with the monolayer TMDCs.

    Appendix A Six-Band Tight-Binding

    The 6× 6 Hamiltonian can be decomposed as[47?50]

    The vectors ajform the lattice basis vectors and extend from every metal to the nearest neighbors.They are all equal in length,given as the lattice constant a=|aj|.The vectors δjhave equal length,given bybut make right angles with the basis vectors.A desirable choice for these vectors areandThen we haveand δ2= δ(0,1).With these conditions,the coordinates of high-symmetry reciprocal lattice points are

    In Table 1,the lengths of these vectors for various TMDCs are introduced.

    Table 1 Lattice constants a= of TMDCs.[50]

    Table 1 Lattice constants a= of TMDCs.[50]

    MoS2 MoSe2 WS2 WSe2 3.160? 3.288? 3.153? 3.260?

    The 3×3 submatrices in the above are given as

    The rest of the matrices are introduced as follows.Starting with TMMjwe have

    Table 2 presents the Slater-Koster parameters needed for this analysis.Data are compiled from Ref.[50]with minor corrections,and numerical results from the model presented in this Appendix are basically identical to the ones reported therein.Comparison to experimental values are already done in many of works,however,where the interested reader is referred to.

    Table 2 Slater-Koster parameters for TMDCs.All parameters are in units of electron-volts(adapted from Ref.[50]with minor corrections).

    Appendix B Expressions forχ(3)

    Following the standard perturbation method to calculate the nonlinear susceptibility tensor,[76?88]the χ(3)is composed of two paramagnetic Π and diamagnetic? parts,as

    where ωΣ= ωr+ωs+ωt.For the case of third-harmonic generation,we have ω = ωr= ωs= ωtand ωΣ=3ω.For the Kerr nonlinearity we have ω = ωr= ?ωs= ωt= ωΣ.

    For the paramagnetic contribution,we have

    The operator P represents all possible intrinsic permutations among frequencies ωr,ωs,ωt.That implies no permutation when all three frequencies are equal,three permutations when only one frequency is different as(ω1,ω2,ω2),(ω2,ω1,ω2),and(ω2,ω2,ω1),and six permutations otherwise as(ω1,ω2,ω3),(ω1,ω3,ω2),(ω2,ω1,ω3),(ω2,ω3,ω1),(ω1,ω2,ω3),and(ω1,ω3,ω2).

    For the diamagnetic contribution,we have

    In all the above relationships,a small positive imaginary part is normally added to all three frequencies,in order to preserve causality as well as to avoid numerical over flow at resonances.Needless to say,coding these relations are not straightforward and need extra care.The double integration over the reciprocal lattice could be done by evaluation of the integrand over a discrete triangular grid and multiplying by element sizes.Employing adaptive numerical integration methods is not practical because of the extremely large numerical burden.

    Appendix C Imaginary Squeezing

    If the third-order nonlinear susceptibility χ(3)is imaginary,then the two-photon absorption is dominant over the Kerr effect.While the Kerr effect could in principle produce unlimited squeezing,two-photon absorption may cause a limited squeezing of shot noise.This effect has been noticed by a number of authors in the past.[158?164]However,in all these works the temporal evolution of squeeze parameter is considered,in which squeezing generally increases with propagation,reaching an ultimate value in the limit of in finite propagation in a two-photon absorbing medium.That type of analysis is useful for nonlinearly lossy long fibers.

    Here,for the case of a ring resonator with constant propagation length,it is the intensity,or equivalently,the cavity photon number N,of the input beam,which could be varied.In what follows,we show that squeezing generally increases with the input power, first within the approximation of negligible loss over propagation,that is constant ζ.

    (i)Non-Rotated Squeezing

    The wavefunction corresponding to the squeezed coherent statewith α being the complex coherent state number,generated by an interaction of the type(10)or equivalently produced from a coherent state such asby application of the operator(ζ)(24),is given by[165?167]

    in which x is in the units of zero-point fluctuations xzp,x0= ?[α],p0= ?[α],and

    It is straightforward to verify via Fourier transformation that the momentum representation would be given by

    Since we have assumed that ζ=is,we have

    Therefore,the probablity distribution in position is given by

    where we have noticed that h is real-valued.Similarly,for the probablity distribution in the momentum representation we obtain

    Comparing to the Gaussian distribution of a simple coherent state,[97]we may deduce the noise squeezing ratio ? in position and φ in momentum,also known as the Fano factor,[164]respectively by solving the equations

    Interestingly,for a purely real ζ with ?[ζ]=0,we haveh=0,and this expression takes the simple solution ? =exp(ζ)and φ =exp(?ζ).This corresponds to a minimum uncertainty squeezed packet

    It is easy to observe that for large values of s,we getwhich represents unlimited desqueezing and shot noise increase.A plot of ? and φ versus s is illustrated in Fig.9.The momentum quadrature exhibits a minimum of φmin= ?2.73 dB at smin=0.797.

    Fig.9 Variation of squeezing amplitudes versus absolute value of purely imaginary squeeze parameter.Both quadratures desqueeze for large s.

    (ii)π/4-Rotated Squeezing

    While non-rotated quadratures both appear to be strongly desqueezed,the fact is that the true squeezing actually happens along the π/4-rotated quadratures.To illustrate this,we need to rigorously calculate the Wigner function of the squeezed vacuum first.This is given by

    where θ=x+ip, β =a+ib, μ =cosh|ζ|,and ν =exp(i∠ζ)sinh|ζ|,where ∠ζ stads for the phase of ζ.Using the fact that ζ=is,and some straightforward but significant algebra,we get the Wigner distribution of purely-imaginary squeezed vacuum as

    The probablity distributions obtained as[97]

    from the above Wigner function both appear to be desqueezed almost equally,as discussed in the above and shown in Fig.9 for large s.Here,the displacement operator(α),which produces a coherent state simply transforms x → x? ?[α]and p→ p??[α].In the the timedependent case,we only need to replace α =Aexp(iωt)where A is the amplitude of the coherent source and ω is the angular frequency of light.The typical resulting Wigner function for the squeezed coherent state is illustrated in Fig.10.

    Fig.10 Wigner distribution of the squeezed coherent state with α =1 and ζ=i.Both non-rotated orthogonal quadratures(x,p)appear to be desqueezed almost to the same amount.

    Now,we can define the π/4-rotated coordinates as

    In this new system of coordinates,we obtain the Wigner distribution of the squeezed stateas

    with the respective squeeze ratios

    plotted in Fig.11,and giving the uncertainty product

    which refers to the minimum-uncertainty packet.

    Fig.11 Variation of squeezing amplitudes in the π/4-rotated system versus absolute value of purely imaginary squeeze parameter.One rotated quadrature is squeezed and the other is strongly desqueezed.

    Appendix D Units of Surface Quantities

    While the linear bulk susceptibility χ(1)should be dimensionless,the surface linear susceptibility χ(1)2Dis not,actually having the dimension of length.But it is preferable to write the dimension as m ·?or ??where the redundant dimensionless square notation?emphasizes a 2D surface quantity.Similarly,the dimension of nonlinear susceptibility χ(3)2Dwould be preferably m3·V?2?instead of m3·V?2.It follows then,that the preferable dimension of linear and nonlinear sheet conductivity σ(3)2Dmust be respectively ??1?and ??1m2·V?2?.

    Acknowledgments

    The author would like to thank Dr.Habib Rostami at Istituto Italiano di Tecnologia(IIT)as well as Dr.Rafael Rold′an and Dr.Pablo San-Jose at Instituto de Ciencia de Materiales de Madrid(CSIC)for discussions and assistance on the Tight-Binding method.Enlightening discussions with Prof.Steven Louie from University of California at Berkeley,Dr.Bernhard Urbaszek at Laboratoire de Physique et Chimie des Nano-objets,Dr.Igor Aharonovich at University of Technology Sydney,as well as Dr.Anshuman Kumar and Cl′ement Javerzac-Galy at′Ecole Polytechnique F′ed′erale de Lausanne(EPFL)on excitonic emission effects is highly appreciated.

    The author is highly indebted to Hiwa Mahmoudi at Institute of Electrodynamics,Microwave and Circuit Engineering in Technische Universit?t Wien,and in particular,the Laboratory for Quantum Foundations and Quantum Information on the Nano-and Microscale in Vienna Center for Quantum Science and Technology(VCQ)at Universit?t Wien for their warm and receptive hospitality during the preparation of this article.

    This paper is dedicated to the celebrated artist,Anastasia Huppmann.

    av国产免费在线观看| 最近在线观看免费完整版| 欧美丝袜亚洲另类 | 99热精品在线国产| 亚洲专区国产一区二区| 在线观看一区二区三区| 欧美一区二区亚洲| 中文字幕精品亚洲无线码一区| 日本在线视频免费播放| 老熟妇乱子伦视频在线观看| 日韩欧美精品v在线| 免费高清视频大片| 久久久色成人| 亚洲欧美日韩高清在线视频| 麻豆成人av在线观看| 美女免费视频网站| 最近视频中文字幕2019在线8| 国产精品嫩草影院av在线观看 | 国内精品久久久久精免费| 国内精品一区二区在线观看| 午夜影院日韩av| 精品久久国产蜜桃| 亚洲国产色片| 啦啦啦观看免费观看视频高清| 老熟妇乱子伦视频在线观看| www日本黄色视频网| or卡值多少钱| 一本一本综合久久| 国产成人a区在线观看| 免费av不卡在线播放| 婷婷色综合大香蕉| 天美传媒精品一区二区| 美女 人体艺术 gogo| 国产成人a区在线观看| 久久国产精品影院| 欧美午夜高清在线| 99久久无色码亚洲精品果冻| 亚洲avbb在线观看| 亚洲av五月六月丁香网| 男女下面进入的视频免费午夜| 免费观看人在逋| 日韩 亚洲 欧美在线| 精品久久久久久久久久免费视频| 亚洲一区二区三区色噜噜| 色综合亚洲欧美另类图片| 久久6这里有精品| 性欧美人与动物交配| 赤兔流量卡办理| 精品一区二区三区视频在线观看免费| 搡女人真爽免费视频火全软件 | 人人妻人人澡欧美一区二区| 舔av片在线| 免费在线观看亚洲国产| 久久久久免费精品人妻一区二区| 2021天堂中文幕一二区在线观| 色综合站精品国产| 久久久国产成人免费| 国产成人a区在线观看| 欧美性猛交黑人性爽| 可以在线观看的亚洲视频| 看免费av毛片| 日韩精品青青久久久久久| 国产一区二区激情短视频| 麻豆成人午夜福利视频| 中文字幕av在线有码专区| 好男人电影高清在线观看| 亚洲内射少妇av| 久久午夜亚洲精品久久| 日韩大尺度精品在线看网址| 97碰自拍视频| 亚洲性夜色夜夜综合| 99久久精品国产亚洲精品| netflix在线观看网站| 男人舔奶头视频| 国产亚洲欧美在线一区二区| 成年免费大片在线观看| 老熟妇乱子伦视频在线观看| 亚洲18禁久久av| 麻豆国产av国片精品| 日韩大尺度精品在线看网址| 国产乱人视频| 搡老妇女老女人老熟妇| 别揉我奶头 嗯啊视频| 免费观看的影片在线观看| 可以在线观看的亚洲视频| 日韩欧美国产在线观看| av福利片在线观看| 久久精品国产自在天天线| 欧美日韩乱码在线| 欧美色欧美亚洲另类二区| 亚洲久久久久久中文字幕| 老司机福利观看| 久久久久久久久中文| 永久网站在线| 亚洲av美国av| 久久午夜亚洲精品久久| 美女被艹到高潮喷水动态| 午夜福利在线观看免费完整高清在 | 美女免费视频网站| 有码 亚洲区| 少妇人妻一区二区三区视频| 18禁在线播放成人免费| 一边摸一边抽搐一进一小说| 亚洲avbb在线观看| 久久久久国产精品人妻aⅴ院| 国产av在哪里看| 九色成人免费人妻av| 日韩高清综合在线| 中文在线观看免费www的网站| 日韩精品中文字幕看吧| ponron亚洲| 国产精品,欧美在线| 日本a在线网址| 国产成人aa在线观看| 人妻丰满熟妇av一区二区三区| 国产视频一区二区在线看| 男女下面进入的视频免费午夜| 最近中文字幕高清免费大全6 | 十八禁人妻一区二区| 日本三级黄在线观看| 在线a可以看的网站| 中文字幕人成人乱码亚洲影| 美女cb高潮喷水在线观看| 乱码一卡2卡4卡精品| 亚洲专区国产一区二区| 国产精品一区二区三区四区久久| 日韩有码中文字幕| 看黄色毛片网站| 欧美又色又爽又黄视频| 色综合站精品国产| 欧美3d第一页| 国产精品国产高清国产av| 日韩中文字幕欧美一区二区| 噜噜噜噜噜久久久久久91| 国产麻豆成人av免费视频| 精品人妻一区二区三区麻豆 | 简卡轻食公司| 18+在线观看网站| 亚洲 国产 在线| 男女之事视频高清在线观看| 日本在线视频免费播放| 中文字幕av在线有码专区| 亚洲国产精品999在线| 观看美女的网站| 欧美黄色淫秽网站| 波多野结衣巨乳人妻| 亚洲成人中文字幕在线播放| 最近最新免费中文字幕在线| 又黄又爽又刺激的免费视频.| 身体一侧抽搐| 亚洲内射少妇av| 日韩有码中文字幕| 日韩精品中文字幕看吧| 97超视频在线观看视频| 热99国产精品久久久久久7| 老女人水多毛片| 亚洲精品乱码久久久久久按摩| 久久久久久久午夜电影| 美女高潮的动态| 亚洲国产精品999| 91aial.com中文字幕在线观看| 国产精品精品国产色婷婷| 韩国av在线不卡| 赤兔流量卡办理| 亚洲av欧美aⅴ国产| 99久久精品热视频| 日本欧美国产在线视频| 青春草视频在线免费观看| 国产永久视频网站| 深爱激情五月婷婷| 可以在线观看毛片的网站| 国内精品美女久久久久久| 成人欧美大片| 亚洲一区二区三区欧美精品 | 中文资源天堂在线| 亚洲国产日韩一区二区| 毛片女人毛片| 亚洲精品国产av成人精品| 国产亚洲精品久久久com| 国产一区亚洲一区在线观看| 午夜福利在线在线| 日韩制服骚丝袜av| 免费电影在线观看免费观看| 免费观看在线日韩| 国内少妇人妻偷人精品xxx网站| 日本色播在线视频| av免费在线看不卡| 午夜免费鲁丝| 成年人午夜在线观看视频| 亚洲精品乱码久久久v下载方式| 国产亚洲av片在线观看秒播厂| 男女那种视频在线观看| 97超碰精品成人国产| 亚洲,欧美,日韩| 成人免费观看视频高清| 国产探花极品一区二区| 久久久久久九九精品二区国产| 日韩制服骚丝袜av| 少妇的逼好多水| 久久6这里有精品| 久久久久网色| 99久久人妻综合| 赤兔流量卡办理| 人人妻人人澡人人爽人人夜夜| 真实男女啪啪啪动态图| 国产成人a区在线观看| 狠狠精品人妻久久久久久综合| 国产中年淑女户外野战色| 国语对白做爰xxxⅹ性视频网站| 国国产精品蜜臀av免费| 色哟哟·www| av播播在线观看一区| 亚洲av欧美aⅴ国产| 人人妻人人澡人人爽人人夜夜| 女人久久www免费人成看片| 国产一区有黄有色的免费视频| 亚洲图色成人| 女人被狂操c到高潮| 成人无遮挡网站| 性插视频无遮挡在线免费观看| 男女边吃奶边做爰视频| 七月丁香在线播放| 午夜亚洲福利在线播放| 久久久久性生活片| 麻豆成人av视频| 免费av观看视频| 亚洲电影在线观看av| 亚洲天堂av无毛| 在线观看美女被高潮喷水网站| 免费av不卡在线播放| 99久久精品一区二区三区| 国产一区有黄有色的免费视频| 别揉我奶头 嗯啊视频| 国产亚洲av片在线观看秒播厂| 白带黄色成豆腐渣| 成人黄色视频免费在线看| av女优亚洲男人天堂| 黄色视频在线播放观看不卡| 超碰97精品在线观看| 久久精品熟女亚洲av麻豆精品| 久久久亚洲精品成人影院| 夫妻午夜视频| 国产老妇伦熟女老妇高清| 18+在线观看网站| 你懂的网址亚洲精品在线观看| 男女下面进入的视频免费午夜| 偷拍熟女少妇极品色| 精品一区二区三区视频在线| 日韩 亚洲 欧美在线| 色视频www国产| 国产 一区 欧美 日韩| 国产一区二区三区综合在线观看 | 亚洲国产欧美人成| 一区二区三区免费毛片| 韩国高清视频一区二区三区| 精品视频人人做人人爽| 在线精品无人区一区二区三 | 日韩视频在线欧美| 老司机影院毛片| 亚洲精品成人久久久久久| 天美传媒精品一区二区| 国产一级毛片在线| 99久久精品一区二区三区| 高清午夜精品一区二区三区| 欧美bdsm另类| 免费大片18禁| 热99国产精品久久久久久7| 亚洲久久久久久中文字幕| 国产一区二区三区综合在线观看 | 一区二区av电影网| 国产极品天堂在线| 亚洲一级一片aⅴ在线观看| 精品久久久久久久末码| 直男gayav资源| 亚洲国产精品999| av福利片在线观看| eeuss影院久久| 国产一区二区亚洲精品在线观看| 综合色丁香网| 最后的刺客免费高清国语| 亚洲精品第二区| 九九爱精品视频在线观看| 人妻 亚洲 视频| 大香蕉久久网| av在线蜜桃| 看非洲黑人一级黄片| 爱豆传媒免费全集在线观看| 在线观看国产h片| 性色avwww在线观看| 麻豆乱淫一区二区| 国产色爽女视频免费观看| 精品少妇黑人巨大在线播放| 人妻夜夜爽99麻豆av| 国产精品秋霞免费鲁丝片| 日韩av免费高清视频| 国产 精品1| 91精品一卡2卡3卡4卡| 国产高潮美女av| 一级a做视频免费观看| 99热这里只有是精品在线观看| 少妇的逼水好多| 大陆偷拍与自拍| 亚州av有码| 97热精品久久久久久| 精品人妻熟女av久视频| 国产日韩欧美在线精品| 国产有黄有色有爽视频| 日韩不卡一区二区三区视频在线| 一级片'在线观看视频| 久久99热这里只频精品6学生| 国产中年淑女户外野战色| 真实男女啪啪啪动态图| 欧美区成人在线视频| 久久久久久伊人网av| 亚洲精品影视一区二区三区av| 少妇 在线观看| 欧美亚洲 丝袜 人妻 在线| 国产极品天堂在线| 国产成人免费观看mmmm| 亚洲欧洲日产国产| 精品人妻熟女av久视频| 国产久久久一区二区三区| 国产欧美日韩精品一区二区| 国产一区二区三区综合在线观看 | 亚洲电影在线观看av| 91aial.com中文字幕在线观看| 亚洲自偷自拍三级| 日本爱情动作片www.在线观看| 一本一本综合久久| 欧美日韩亚洲高清精品| 中文字幕久久专区| 国产精品久久久久久久电影| 国产 一区精品| 男女边吃奶边做爰视频| 国产成人a∨麻豆精品| 免费黄色在线免费观看| 一级爰片在线观看| 亚洲欧美日韩东京热| 九草在线视频观看| 午夜福利在线在线| 自拍欧美九色日韩亚洲蝌蚪91 | 超碰av人人做人人爽久久| 嫩草影院入口| 亚洲精品第二区| 男女无遮挡免费网站观看| 黄色日韩在线| 麻豆成人午夜福利视频| 在线免费观看不下载黄p国产| 黄色欧美视频在线观看| 久久鲁丝午夜福利片| 22中文网久久字幕| 亚洲欧洲国产日韩| 男女国产视频网站| 精品一区二区三区视频在线| 欧美极品一区二区三区四区| 国产午夜精品一二区理论片| av福利片在线观看| 视频区图区小说| 狂野欧美激情性xxxx在线观看| 国产精品久久久久久精品电影小说 | 日韩欧美一区视频在线观看 | 国产成人精品久久久久久| 日日啪夜夜撸| 午夜福利高清视频| 亚洲综合色惰| a级毛片免费高清观看在线播放| 美女内射精品一级片tv| 又爽又黄无遮挡网站| 久久久欧美国产精品| 肉色欧美久久久久久久蜜桃 | 午夜福利视频精品| 麻豆精品久久久久久蜜桃| 久久久精品免费免费高清| 免费观看在线日韩| 九九久久精品国产亚洲av麻豆| 青春草亚洲视频在线观看| 好男人在线观看高清免费视频| 丝袜脚勾引网站| 联通29元200g的流量卡| 亚洲欧美精品自产自拍| 日韩av在线免费看完整版不卡| 各种免费的搞黄视频| 精品久久久精品久久久| 99热这里只有精品一区| 在线播放无遮挡| 久久久久久九九精品二区国产| 日本色播在线视频| 国国产精品蜜臀av免费| 国产美女午夜福利| 成人一区二区视频在线观看| 全区人妻精品视频| 亚洲欧美成人综合另类久久久| 深爱激情五月婷婷| 黄色视频在线播放观看不卡| 午夜免费男女啪啪视频观看| 国内精品美女久久久久久| 欧美日韩视频精品一区| 欧美一区二区亚洲| 久久精品国产亚洲网站| 久久这里有精品视频免费| 国产老妇女一区| 大话2 男鬼变身卡| 精品一区二区免费观看| av卡一久久| 99热这里只有精品一区| 亚洲欧美日韩东京热| 亚洲精品,欧美精品| 97在线视频观看| 亚洲欧美日韩另类电影网站 | 日本色播在线视频| 久久人人爽人人爽人人片va| 丰满少妇做爰视频| 亚洲国产精品国产精品| 亚洲欧美日韩东京热| 久久99蜜桃精品久久| 欧美性猛交╳xxx乱大交人| 少妇熟女欧美另类| 免费看不卡的av| 日韩av免费高清视频| 搡老乐熟女国产| 欧美高清性xxxxhd video| 亚洲色图av天堂| 欧美少妇被猛烈插入视频| 国产综合懂色| 欧美区成人在线视频| 91精品伊人久久大香线蕉| 欧美高清成人免费视频www| 2018国产大陆天天弄谢| av在线老鸭窝| 我的女老师完整版在线观看| 精品99又大又爽又粗少妇毛片| 97精品久久久久久久久久精品| 亚洲人成网站在线播| 永久网站在线| 嫩草影院入口| 国产成人一区二区在线| 少妇的逼好多水| 国产探花在线观看一区二区| 午夜激情久久久久久久| 久久久久久久久久久丰满| 99热国产这里只有精品6| 国产精品99久久99久久久不卡 | 国产伦在线观看视频一区| 深爱激情五月婷婷| 91狼人影院| 午夜亚洲福利在线播放| 久久精品国产亚洲av涩爱| 久久人人爽人人爽人人片va| 国产色婷婷99| 成人国产av品久久久| 校园人妻丝袜中文字幕| 新久久久久国产一级毛片| 久久99精品国语久久久| 午夜老司机福利剧场| 国产精品不卡视频一区二区| 最新中文字幕久久久久| 麻豆成人午夜福利视频| 人妻制服诱惑在线中文字幕| 黑人高潮一二区| 女人被狂操c到高潮| 欧美xxⅹ黑人| 日韩 亚洲 欧美在线| 人妻制服诱惑在线中文字幕| 久久久成人免费电影| 国产人妻一区二区三区在| 久久鲁丝午夜福利片| 美女主播在线视频| 久久午夜福利片| 亚洲精品乱码久久久v下载方式| 九九在线视频观看精品| 大码成人一级视频| 午夜爱爱视频在线播放| 中文资源天堂在线| 久久精品久久久久久噜噜老黄| 亚洲精品影视一区二区三区av| 在线观看人妻少妇| 国产精品一区二区性色av| 一级a做视频免费观看| 国产一区二区在线观看日韩| 舔av片在线| 国产老妇伦熟女老妇高清| 尾随美女入室| 婷婷色综合大香蕉| 日韩欧美精品免费久久| 男人添女人高潮全过程视频| 亚洲精品自拍成人| 观看美女的网站| 日日啪夜夜爽| 网址你懂的国产日韩在线| 免费不卡的大黄色大毛片视频在线观看| 2021少妇久久久久久久久久久| 国产色婷婷99| 一区二区av电影网| 又爽又黄无遮挡网站| 亚洲av欧美aⅴ国产| 日日摸夜夜添夜夜爱| 99久久精品一区二区三区| 少妇丰满av| 国产美女午夜福利| 最近中文字幕2019免费版| 一边亲一边摸免费视频| 国产v大片淫在线免费观看| 边亲边吃奶的免费视频| 国产精品久久久久久av不卡| 青春草视频在线免费观看| 交换朋友夫妻互换小说| 色播亚洲综合网| 国产v大片淫在线免费观看| 国内精品美女久久久久久| 在线观看国产h片| 亚洲av成人精品一区久久| 日韩av在线免费看完整版不卡| 国精品久久久久久国模美| 另类亚洲欧美激情| 一本色道久久久久久精品综合| 亚洲成人久久爱视频| 波野结衣二区三区在线| 国产一区亚洲一区在线观看| 亚洲精品456在线播放app| av在线蜜桃| 午夜免费鲁丝| 涩涩av久久男人的天堂| 成年女人看的毛片在线观看| 少妇猛男粗大的猛烈进出视频 | 国产一区二区三区av在线| 欧美日韩国产mv在线观看视频 | 伊人久久国产一区二区| 26uuu在线亚洲综合色| 国产v大片淫在线免费观看| 亚洲成人精品中文字幕电影| 午夜爱爱视频在线播放| 身体一侧抽搐| 在现免费观看毛片| 国产午夜福利久久久久久| 国内精品宾馆在线| 国内少妇人妻偷人精品xxx网站| 高清欧美精品videossex| 欧美日韩一区二区视频在线观看视频在线 | 国产一区二区亚洲精品在线观看| 女的被弄到高潮叫床怎么办| 亚洲精品视频女| 免费观看无遮挡的男女| 亚洲,欧美,日韩| 亚洲精品乱久久久久久| 丝瓜视频免费看黄片| 亚洲成人中文字幕在线播放| 女的被弄到高潮叫床怎么办| 中文天堂在线官网| 亚洲va在线va天堂va国产| 成年版毛片免费区| 成人漫画全彩无遮挡| 国产毛片在线视频| 日韩av免费高清视频| 亚洲在线观看片| 久久人人爽人人爽人人片va| 国产精品久久久久久久久免| 男的添女的下面高潮视频| 99久国产av精品国产电影| 波野结衣二区三区在线| 啦啦啦中文免费视频观看日本| 亚洲无线观看免费| 搡老乐熟女国产| a级一级毛片免费在线观看| 精品午夜福利在线看| av黄色大香蕉| 91久久精品电影网| 男插女下体视频免费在线播放| 网址你懂的国产日韩在线| 舔av片在线| 久久亚洲国产成人精品v| 国产精品爽爽va在线观看网站| 小蜜桃在线观看免费完整版高清| 夫妻午夜视频| 最近最新中文字幕大全电影3| 91久久精品电影网| 老司机影院成人| 亚洲av日韩在线播放| 2021少妇久久久久久久久久久| 高清午夜精品一区二区三区| 久久午夜福利片| .国产精品久久| 久久韩国三级中文字幕| 人人妻人人爽人人添夜夜欢视频 | 3wmmmm亚洲av在线观看| 嘟嘟电影网在线观看| 美女主播在线视频| 极品教师在线视频| 插逼视频在线观看| 亚洲高清免费不卡视频| 国产毛片a区久久久久| 亚洲av一区综合| 男人添女人高潮全过程视频| 国产精品成人在线| 国产老妇伦熟女老妇高清| 日韩中字成人| 乱系列少妇在线播放| 精品国产一区二区三区久久久樱花 | 看非洲黑人一级黄片| 欧美3d第一页| 免费高清在线观看视频在线观看| 狂野欧美白嫩少妇大欣赏| 你懂的网址亚洲精品在线观看| 毛片一级片免费看久久久久| 国产一区二区三区综合在线观看 | 最近2019中文字幕mv第一页| 久久99蜜桃精品久久| 成人鲁丝片一二三区免费| 草草在线视频免费看| 高清在线视频一区二区三区| 丰满人妻一区二区三区视频av| 国产一区有黄有色的免费视频| 成人黄色视频免费在线看| 蜜桃亚洲精品一区二区三区| 七月丁香在线播放| tube8黄色片|