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

    The propagation dynamics and stability of an intense laser beam in a radial power-law plasma channel

    2022-01-10 14:51:12XuerenHONG洪學(xué)仁DeshengZHANG張德生JimingGAO高吉明RonganTANG唐榮安PengGUO郭鵬andJukuiXUE薛具奎
    Plasma Science and Technology 2021年12期

    Xueren HONG (洪學(xué)仁),Desheng ZHANG (張德生), Jiming GAO (高吉明),Rongan TANG (唐榮安),?, Peng GUO (郭鵬) and Jukui XUE (薛具奎)

    1 Key Laboratory of Atomic Molecular Physics&Functional Materials of Gansu Province,and College of Physics and Electronic Engineering,Northwest Normal University,Lanzhou 730070,People?s Republic of China

    2 School of Mathematics and Physics,Lanzhou Jiaotong University,Lanzhou 730070,People?s Republic of China

    Abstract By containing ponderomotive self-channeling,the propagation behavior of an intense laser beam and the physical conditions are obtained theoretically in a radial power-law plasma channel.It is found that ponderomotive self-channeling results in the emergence of a solitary wave and catastrophic focusing,which apparently decreases the region for stable propagation in a parameter space of laser power and the ratio of the initial laser spot radius to the channel radius (RLC).Direct numerical simulation confirms the theory of constant propagation,periodic defocusing and focusing oscillations in the parameter space, and reveals a radial instability which prevents the formation of bright and dark solitary waves.The corresponding unstable critical curve is added in the parameter space numerically and the induced unstable region above the unstable critical curve covers that of catastrophic focusing, which shrinks the stable region for laser beams.For the expected constant propagation, the results reveal the need for a low RLC.Further study illustrates that the channel power-law exponent has an obvious effect on the final stable region and laser propagation, for example increasing this exponent can enlarge the stable region significantly,which is beneficial for guiding of the laser and increases the lowest RLC for constant propagation.Our results also show that the initial laser amplitude has an apparent influence on the propagation behavior.

    Keywords: laser beam, power-law channel, propagation dynamics, stability

    1.Introduction

    The propagation of laser beams [1, 2] has received considerable attention due to their important applications, such as plasma-based acceleration [3–10], harmonic generation[11, 12], x-ray lasers [13, 14] and advanced laser fusion schemes [15].A laser beam can propagate in a vacuum by about a Rayleigh length due to natural diffraction.For the above applications, however, more than several Rayleigh lengths are required.Propagation in a plasma may effectively increase the traveling distance.In a uniform plasma, the laser suffers relativistic self-focusing, which can overcome natural diffraction and ensure a long propagation distance if the laser power is larger than a critical value.When the laser power is smaller than this critical value, however, relativistic selffocusing fails to prevent diffraction and other focusing methods are needed.A preformed plasma channel is an effective technique for providing a powerful focusing effect and can play the role of a waveguide tube [7]; preformed plasma channels are widely used for laser guiding in related works.

    As we know,preformed channels are usually designed to be parabolic in the radial direction[7,16,17].Other channels have previously received less attention than the usual parabolic channel.Nevertheless,some special advantages have been reported for non-parabolic channels.For example, a class of Raman-type instabilities may be stabilized in a leaky density channel[18,19]and good stable laser transmission with high quality in inhomogeneous plasmas can be achieved in a hollow channel[20–25].In addition, recent work on the development ofq-Gaussian, super-Gaussian and quadruple Gaussian lasers [26–32] may inspire investigation of laser guiding in non-parabolic channels.In experiments, various channel structures can be achieved by controlling some of the experimental parameters [4, 7, 33–35].Therefore,the investigation of laser propagation in non-parabolic channels is necessary and crucial.Quasi-matched propagation of laser pulses in a plasma channel with a general radial profile has recently been studied [36].Attention to the corrugated channel has revealed rich laser beam behaviors, such as aperiodic oscillation, resonance, beat-like waves and periodic oscillation with multipeaks[37].Our recent work on laser propagation in a radial power-law channel [38], which is a reasonable and general extension of a parabolic density channel, has revealed that constant propagation is sensitive to the power-law exponent.

    In preliminary work [38] on radial power-law plasma channels, however, ponderomotive self-channeling was not considered because of the difficulty of theoretical analysis.Ponderomotive self-channeling is related to the ponderomotive force, which can be viewed as the radiation pressure, i.e.the gradient of the electromagnetic energy density [7].Under the action of this force, electrons are expelled transversely away from regions of high laser intensity in the long pulse length limit and ponderomotive self-channeling emerges.Previous studies on ponderomotive self-channeling[39–41]in parabolic channels revealed a more complex parameter space than in[38].After comparisons, it can be found that ponderomotive self-channeling can result in a parameter region for catastrophic focusing, and apparently decreases the traditional propagation domain in parameter space for a parabolic channel.Therefore,for an investigation into a radial power-law channel ponderomotive self-channeling must be contained.

    After some effort we overcame the theoretical difficulty of treating ponderomotive self-channeling in a radial power-law plasma channel and developed the theory presented in this work.Meanwhile, we also investigated the propagation dynamics by direct simulation of the wave equation with the finite difference method.The laser behaviors and their physical conditions are analyzed.As expected, ponderomotive self-channeling induces the emergence of solitary waves and catastrophic focusing,which obviously decreases the area for stable laser propagation.Some interesting findings were obtained.(1) Theoretical results indicate that the increase in the channel power-law exponent may dramatically decrease the domain for catastrophic focusing and be conducive to laser propagation.(2) Direct numerical simulation confirms the theory of constant propagation,periodic defocusing and focusing oscillations in parameter space, and reveals a radial instability which prevents the formation of bright and dark solitary waves.We added the corresponding unstable critical curve to the parameter space numerically and the induced unstable region above the unstable curve covers that of catastrophic focusing, which further shrinks the stable region for laser beams.(3) For the usually expected constant propagation,the results reveal the need for a lowest RLC(ratio of the initial laser spot radius to the channel radius).(4) Increasing the channel power-law exponent can enlarge the final stable region significantly and also increase the lowest RLC for constant propagation.It should be noted that findings(2)and(3)recover the previous knowledge even for the case of a parabolic channel,such as in a previous study of solitary waves [39] and other studies associated with constant propagation.Our studies also demonstrate the obvious effects of the initial amplitude of the laser on the propagation behavior.These findings are important for the guiding of lasers and other related applications.

    This paper is organized as follows.The basic model and evolution equations are given in section 2.The propagation types and characteristics are obtained theoretically in section 3.Numerical simulation of the wave equation and further analysis are done in section 4.A summary is provided in section 5.

    2.Basic model and evolution equations

    Taking into account ponderomotive self-channeling and the weakly relativistic limit (|a|2?1, where a is the normalized vector potential and its normalization is presented after equation (2)), an intense laser beam propagating in an underdense plasma (n/nc<1, wherenis the density of the underdense plasma andncis the critical density) is described by the following two equations [40–44]:

    where a(r,z,t) is the vector potential and is normalized bym0c2/e,m0andeare the electron rest mass and charge,respectively,cis the velocity of light in a vacuum,ωp=(4πn0e2m0)12is the plasma frequency and δnis the perturbed electron density.n(r) denotes the power-law density channel

    wheren0is the initial axial density of the electron,rchis the channel radius andd>0 is the power-law exponent.Whend=2,equation(3)is reduced to the familiar parabolic case.The Coulomb gauge ?·a=0 is adopted in deriving equations (1)and (2).Obviously, equation (1) is the wave equation for the laser field and equation (2) describes the perturbed electron density associated with the excitation of the wakefield.For simplicity, we consider the dynamics in the long pulse length limit ωpτl?1,with τlthe laser pulse width.Thus,equation(2)is simplified as [40]

    wherekp=ωp/crepresents the plasma wave number.Equation (4) shows the formation of ponderomotive self-channeling [7], which was ignored in our earlier work [38] for simplicity.

    Here, we study a circularly polarized laser beam propagating along thezdirection in a plasma.So, the form of the vector potential is set as

    wherea(r,z,t)represents a complex slowly varying envelope,ω0represents the laser center frequency,k0=vgω0/c2represents the wave number of the laser center,vgrepresents the laser group velocity and ?xand ?yare unit coordinate vectors;c.c.is the complex conjugate.

    By using equations (3)–(5) and the coordinate transformationz=z, τ=t?z/vgwith paraxial treatment and the slowly varying envelope assumption |?a/?τ|?|ωpa|,equation (1) is changed to

    In equation (6), the fifth term is caused by ponderomotive self-channeling.Without this term, equation (6) reduces to equation (4) in [38].Settingd=2 corresponding to the parabolic case and adopting the assumption of the linear dispersion relationequation (6) is simplified to become equation (3) in [39] and equation (5)in [40].

    To get the solution to equation (6), we adopt the following trial function:

    wherear,rs,band φ represent the amplitude, spot radius,spatial chirp parameter and phase shift of the laser beam in plasmas at any distancez, respectively.By applying the variational method [38, 40, 41, 44–46] to equation (6), one can obtain the following equations forar,rs,band φ:

    3.Propagation types and characteristics

    Equation (8) shows thatas2rs2=a02r02, witha0being the initial (z=0) laser amplitude andr0the initial laser spot radius.Usingarrs=a0r0and combining equations (9) and(10), the normalized evolution equation of the laser spot radius becomes

    Integrating equation (12) with the initial conditionrs|z=0=1 and drs/dz|z=0=0, one obtains

    where the expression for the Sagdeev potential is

    with

    The form of equation (13) is similar to the energy conservation equation for a particle with ‘timez’ and ‘positionrs’.Thus, we can investigate equation(13)by analyzing the Sagdeev potential.

    To get the features of the Sagdeev potential,we begin with the investigation of the two equationsV(rs)=0 andV′ (rs)=0.For the parabolic case, it is easy to list the roots ofV(rs)=0[40, 41, 44], so that one can get the propagation characteristics directly.However, it seems that there is no algebraic theory to give their roots directly for a generald.Nevertheless,we can get some important properties of the roots under the appropriate conditions 00 andd>1.A detailed analysis of the properties is included in the Appendix and the main results are

    (1) forV(rs)=0, the value 1 is one root and there are at most two other positive roots

    (2)V′(rs)=0may have no positive root,one positive root or two positive roots.

    According to the two properties, considering the asymptotic behaviors ofV(rs), i.e.V(rs)→∞whenrs→∞ andV(rs)→?∞whenrs→0, the qualitative figures of the Sagdeev potential aroundV(rs)=0 can be drawn as in figure 1,which includes nine kinds of profile.Figures 1(a)–(c)correspond to the case of three roots ofV(rs)=0.Figures 1(d)–(g)represent the case when there are two roots.Figures 1(h)and(i)are for the case of only one root.For the one-root case,it should be mentioned that there may exist other profiles besides figures 1(h) and (i).However, it is not necessary to investigate these other mathematical profiles because the basic physical behaviors in the one-root case are adequately included in figures 1(h) and (i).

    Figure 1.Diagram of V(rs) around V(rs)=0.Parts (a), (b) and (c)correspond to the case when there are three roots of V(rs)=0.Parts(d), (e), (f) and (g) correspond to the case of two roots of V(rs)=0.Parts(h)and(i)correspond to the case of only one root.The marked dots represent the initial position of a test particle with rs|z=0=1 and drs/dz|z=0=0.

    Figure 2.Parameter regions of laser propagation types.The black,green and red curves represent the conditions for constant propagation, bright solitary wave and dark solitary wave respectively.The dashed curves are for d=2 and the types of curves for other d values are labeled in the figure.For a fixed d, the areas beneath the green and black curves, between the black and red curves, and above the green and red curves correspond to periodic defocusing oscillation,periodic focusing oscillation and catastrophic focusing, respectively.The small black, red, green and blue rectangles correspond to d=1.5, 2, 3 and 4, respectively, and represent the phase transition points for the division of the parameter regions.The empty circles are parameter points used for simulation in figure 6 of section 4.Here, a0=0.3.

    Figure 3.Display of the laser propagation types shown in figure 2.Parts (a), (b) and (c) represent constant propagation, dark solitary waves and bright solitary waves respectively.Parts (d) shows the variation of laser behavior with increase in laser power p for d=4,and it can be seen that periodic defocusing oscillation, constant propagation, periodic focusing oscillation, a dark solitary wave and catastrophic focusing may occur with the increase in p.Here,a0=0.3.

    Figure 4.The effect of a0 on the parameter regions.The black,green and red curves represent the conditions for constant propagation,bright solitary wave and dark solitary wave, respectively.The dashed curves are for a0=0.3 and the types of curves for other a0 are labeled in the figure.For a fixed a0, the areas beneath the green and black curves, between the black and red curves, and above the green and red curves correspond to periodic defocusing oscillation,periodic focusing oscillation and catastrophic focusing,respectively.The small black, red, green and blue rectangles correspond to a0=0.2, 0.3, 0.4 and 0.5, respectively, and represent the phase transition points for the division of the parameter regions.The empty circle is a parameter point used for simulation in figure 9 of section 4.Here, d=4.

    From figure 1, one can directly get the following propagation types and characteristics of the laser beam.

    (1) Catastrophic focusing.In the cases shown in figures 1(a), (f), (h) and (i), the particle will move to the positionrs=0 considering the asymptotic behavior ofV(rs),i.e.V(rs)→?∞whenrs→0.In this case,the laser beam will catastrophically focus at the position Obviously, the pointrs=0 will never be reached because the intensity of laser will large enough to exceed the weakly relativistic limit when the laser spot radius continuously decreases.The condition of this case isV′ (rs)∣rs=1>0and ?ra?(0, 1),V(ra)<0 orV′ (rs)∣rs=1=0andV″(rs)∣rs=1=0, i.e.

    (2) Periodic defocusing oscillation.In the case shown in figure 1(b), the motion of the particle is bounded between the positions 1 andr1, wherer1>1.It can be seen thatrs≥1 in this kind of evolution ofrs.The condition for this case isi.e.

    The corresponding spatial period of the oscillation is

    (3) Periodic focusing oscillation.In the case shown in figure 1(c), the motion of the particle is bounded between the positionsr2and 1, wherer2<1.It can be seen thatrs≤1 in this kind of evolution ofrs.The condition of this case isV′(rs)∣rs=1>0and ?rb?(0,1)satisfiesi.e.

    The corresponding spatial period is

    (4) Constant propagation.In the case shown in figure 1(g),the particle locates at the minimum point ofV(rs), initially with with zero velocity.So the particle will stay at rest.This means that the laser spot radius does not change and the type of propagation is constant propagation.The corresponding power can be obtained by solvingV′ (rs)∣rs=1=0underV″ (rs)∣rs=1>0, which yields

    with

    (5) Solitary wave.In the figures 1(d) and (e), it can be concluded that the maximum pointsrs=1 in figure 1(d)andrs=r3in figure 1(e) are all saddle points by linear stability analysis.According to the corresponding nonlinear theory[44],equation(13)has solitary wave solutions in these cases.Case figure 1(d) corresponds to a bright solitary wave and the condition isV′ (rs)∣rs=1=0underV″ (rs)∣rs=1<0, i.e.

    with

    Case figure 1(e)corresponds to a dark solitary wave and the condition ?r3?(0,1)satisfeisV(rs)∣rs=r3=0,V′ (rs)∣rs=r3=0 andV″ (rs)∣rs=r3<0, i.e.

    As we known,solitary waves are typical solutions for a variety of nonlinear systems, such as in the deep ocean[47], plasmas [48] and others.Here, a bright solitary wave means a single‘hump’ofrs,while a dark solitary wave represents a single‘pit’ofrs.One should note that these are quite different.

    Obviously,the propagation types in this power-law density channel are very rich.Meanwhile, the condition of each propagation type strongly depends ond.To make a further study of the propagation characteristics and the channel effects, we plot the propagation characteristics map in parameter space (p,r0/rch), i.e.figure 2, in which the parameter regions corresponding to different types are classified using the above conditions.The correctness of figure 2,i.e.the above theory based on equation (12), is well verified by solving this ordinary differential equation (12) numerically and we display some of the solutions in figure 3.With the parameters presented in figure 2,the behaviors of the constant propagation, dark solitary waves and bright solitary waves are respectively shown in figures 3(a),(b)and(c)for differentd.All the above five types are displayed in figure 3(d)by settingd=4 as an example.With the increase in laser powerp, periodic defocusing oscillation, constant propagation,periodic focusing oscillation,dark solitary waves and catastrophic focusing are illustrated in figure 3(d).For a fixedd,it can seen that figure 3(b) displays a string of dark solitary waves instead of a single one,which is related to the instability of the saddle point in the nonlinear system.

    Let us first observe and analyze the parameter regions for a givendin figure 2.It can be seen that there exists a phase transition pointmarked with a rectangle for a givend.Whenr0/rchis larger than the horizontal ordinate of this point, there are five regions which correspond to periodic defocusing oscillation, constant propagation, periodic focusing oscillation, a dark solitary wave and catastrophic focusing, respectively.Whenhowever,there are only three regions for periodic defocusing oscillation,a bright solitary wave and catastrophic focusing,respectively.That is to say, the motions of constant propagation and periodic focusing oscillation are impossible when the injected spot radius is less than a small value.We should emphasize that this finding was omitted and has not been considered in previous work,even in the parabolic channel.Comparing figure 2 with figure 1 in[38],it is found that the parameter regions are much more complex than in[38]and this is induced by containing the ponderomotive self-channeling.Besides solitary waves, we think that the main contribution of ponderomotive self-channeling is catastrophic focusing because it occupies a large area which dramatically decreases the region for stable propagation.

    Now, we attend to the influence ofdon parameter areas in figure 2.Obviously, the phase transition point moves towards largerr0/rchwhendincreases.Increasingdflattens the central channel part, which decreases preformed channel focusing whenr0/rchis small.To recover from the loss of the channel focusing effect, a relatively largerr0/rchis naturally needed for the new phase transition point with largerd.For constant propagation,there exists a critical valuer0/rch|1(labeled in figure 2).Whenr0/rchr0/rch|1;this was stated in[38]so we pay no attention to it here.For the new discovery in this work, there is also another critical valuer0/rch|2(≈1), below whichpfor the dark solitary wave apparently increases withd.Obviously,r0/rch|2>r0/rch|1.We think that this is because the laser spot radius of the dark solitary wave is always less than or equal to that for the constant propagation.Thus, for the same real critical radial length of the laser with regard to the variation of channel focusing withd, a relatively larger initial spot radius is needed for a dark solitary wave, andr0/rch|2>r0/rch|1is reasonable.In addition, it should be mentioned that the intersection points atr0/rch|1andr0/rch|2are both very small intersection regions.So the two critical valuesr0/rch|1andr0/rch|2should strictly be two very small regions aroundr0/rch|1[38]andr0/rch|2.Nevertheless,it is natural to think that our adoption of the two critical valuesr0/rch|1andr0/rch|2is a very good approximation and helps to illustrate the physics.We should possibly pay more attention to the region of catastrophic focusing, which is sensitive to the channel power-law exponentd.With increasingd, the area of catastrophic focusing shrinks significantly and thus the region for stable propagation enlarges.

    Besidesd, it can be found from equations (16)–(25) that the parameter regions also depend ona0.In figure 4, we plot the parameter regions for different values ofa0.It can be seen that with increase ina0the phase transition point moves to the bottom right, causing an obvious change in the phase diagram.Correspondingly, the curve of the solitary wave (containing bright and dark solitary waves) moves down,increasing the area of catastrophic focusing; the area of periodic oscillation (containing periodic defocusing oscillation and periodic focusing oscillation) gradually decreases.Meanwhile, the constant propagation curve moves up, inducing an increase in the area of periodic defocusing oscillation and a decrease in the area of periodic focusing oscillation.Compared with the dark solitary wave curve,a0has a more obvious effect on the curves for constant propagation and the bright solitary wave.Now we analyze the reason for the main variations of the parameter regions.According to equation (12), whena0increases and other parameters are fixed, the ponderomotive self-channeling effect is enhanced while the preformed channel focusing effect is weakened.The variation of parameter regions witha0is the result of competition between these two effects.For the parameters around the constant propagation curve,rsvaries about 1.In this case,the reduction of the preformed channel focusing effect induced bya0is more obvious than the enhancement of the ponderomotive self-channeling effect.This causes the constant propagation curve to move up.For the values around the minimumrsof a dark solitary wave, however, the enhancement of the ponderomotive self-channeling effect is relatively obvious this time, and the dark solitary curve moves down.

    In figures 2 and 4, the laser propagation behaviors are analyzed in the phase diagram(p,r0/rch),wherepandr0/rchare dimensionless parameters and their variation may mean a complete change of the laser beam and plasma background.Now we study the propagation behaviors from another perspective, i.e.paying attention to the laser behaviors in the parameter space(a0,r0)of the initial laser under the condition of a given plasma backgroundn0andrch.In order to display this simply,we actually use(a0,r0/rch)in the following part.Therefore, by usingn0=0.02ncandrch=30 μm, a corresponding phase diagram is drawn in figure 5(a).As expected,the parameter space is divided into five regions through a phase transition point,and the distribution of these parameter regions is also similar to figure 2 from a topological point of view.However, figure 5(a) reveals some important information that cannot be found from figure 2.Firstly, one cannot simply think that constant propagation can be achieved at any the initial spot radiusr0for a given plasma background.Actually, figure 5(a) only shows a small allowed range, i.e.r0/rch?(0.24, 0.38).It is worth noting that this conclusion does not mean that the range(0,1)for constant propagation in figure 2 is irrational, because it is always possible to move and adjust the range ofr0/rchfor constant propagation by changing the plasma background.Secondly, on the curve for constant propagation,a0increases rapidly with the decrease ofr0, In other words,a0is sensitive to changes inr0for constant propagation, especially for smalla0.Finally, the above two limitations of the constant propagation condition reveal the importance of focusing and defocusing oscillations.Figure 5(a)shows that there is a considerable area supporting the two behaviors.Figures 5(b)and(c)show the variations of parameter regions with the changes ofn0andrch, respectively.Asn0andrchincrease separately, parameter regions move to the bottom left.Meanwhile,we note that the changes do not induce variation of the main characteristics revealed by figure 5(a).

    Figure 5.(a) Parameter regions of laser propagation types in the space(a0,r0/rch)for rch=30 μm and n0=0.02nc.The black,green and red curves represent the conditions for constant propagation,bright solitary wave and dark solitary wave, respectively.The areas beneath the green and black curves, between the black and red curves and above the green and red curves correspond to periodic defocusing oscillation,periodic focusing oscillation and catastrophic focusing,respectively.(b)The variation of parameter regions with n0 for rch=20 μm.From right to left, the curves correspond to n0=0.02nc, 0.04nc and 0.06nc, respectively.(c) The variation of parameter regions with r0/rch for n0=0.02nc.From right to left,the curves correspond to rch=20 μm, 30 μm and 40 μm, respectively.Here, d=4 and nc ≈9.84×1026 m?3 for a laser with wavelength λ=1064 nm.

    Figure 6.The modulus of the slowly varying laser envelope |a(r, z)| obtained by numerical simulations with parameters marked by empty circles in figure 2 for d=4, a0=0.3 and r0/rch=0.4.From (a) to (e), p=0.10, 0.23, 0.29, 0.40, 0.58 and 0.59, respectively.

    Figure 7.Image profile of |a(r, z)| obtained by the numerical simulation with laser power p=0.585, d=4 and r0/rch=0.4 around the unstable critical value for laser propagation failure.Part(a)presents the image profile and part(b)presents the corresponding modulus of the laser amplitude along the radial direction at z=7.388.

    Figure 8.Parameter regions of laser propagation types(shown in figure 2)with unstable critical curves obtained by numerical simulation for d=1.5 in (a), 2 in (b), 3 in (c) and 4 in (d).The curves with black rectangles, black circles, black triangles and white circles represent the unstable critical curves corresponding to different values of d.In the area above an unstable critical curve,the laser beam fails to propagate.Other curves and labels are identical to figure 2.Here, a0=0.3.

    Figure 9.The modulus of the slowly varying laser envelope|a(r,z)|obtained by numerical simulation with parameters marked by the empty circle in figure 4 for d=4, r0/rch=0.4 and p=0.29.From (a) to (d) a0=0.2, 0.3, 0.4 and 0.5, respectively.

    Figure 10.Sketch map for the properties of the Sagdeev potential.

    Considering the complexity of the real propagation and the obvious effects of the channel power-law exponentdand the initial laser amplitudea0, it is necessary to make a numerical simulation of equation (6) to verify the theoretical results in the next section.We will see that the numerical results reveal a new unstable critical curve corresponding to radial instability [49, 50] in parameter space (p,r0/rch).The theoretical predictions below this curve are verified,while the formation of solitary waves is prevented because of the radial instability.Many other interesting findings will be found in the next simulation.

    4.Numerical simulation

    To illustrate the propagation characteristics of the laser beam and verify the above theoretical results,we will make a direct numerical simulation of the wave equation, i.e.partial differential equation (6), in this section.Therefore, with the assumption of an initial linear dispersion relationwe can simplify equation (6) under axisymmetric condition and normalize it as

    In this equation,ris normalized byr0andpandzare normalized quantities identical to those in equation (12).The injected initial laser beam is set to bea(r, 0) =a0exp( ?r2).Then equation (26) can be numerically solved using the finite difference method as successfully performed in[28, 37, 38, 51].

    To begin with we make the simulation with parameters marked by empty circles in figure 2 ford=4.The results are presented in figure 6 and the corresponding physical parameters are also labeled.It can be seen that the the simulation results shown in figures 6(a)–(e), corresponding top=0.1–0.58, are identical to the theoretical results shown in figure 2, i.e.figures 6(a) and (b) show periodic defocusing oscillations, figure 6(c) shows constant propagation and figures 6(d) and (e) show periodic focusing oscillations.The measured periods of 4.5 in figure 6(a), 3.3 in figure 6(b), 3.2 in figure 6(d) and 3.4 in figure 6(e) are also approximately equal to the calculated values of the theoretical period formulas given by equations (18) and (20) as 4.2, 3.4, 3.0 and 2.8, where the difference may be relatively large for figure 6(e); we think this is because the corresponding laser power of 0.58 is much nearer the 0.585 of the following radial instability.This good consistency confirms the theoretical results in section 3 and the numerical simulation.However,when the laser powerpexceeds 0.58 and increases to 0.59,the laser beam fails to propagate.Here, we think the appearance of the irregular behaviors (except constant propagation, periodic focusing and defocusing oscillations)means a failure of propagation.With a clear observation of the simulation atp=0.585 shown in figure 7, we find the failure of propagation is related to instability in the radial direction [49, 50] (radial shape change and some microstructures related to the instability can be seen in figure 7(b)).

    Numerical simulations with other values ofr0/rchconfirm this phenomenon.That is,for a fixedd,there should exist an unstable critical curve in figure 2, above which the laser beam fails to propagate stably.Ford=1.5, 2, 3 and 4, we plot these unstable critical curves plus the corresponding parameter regions of figure 2 in figures 8(a)–(d),respectively.The critical curves of instability are obtained through a large number of numerical simulations.The discrete values ofr0/rchare 0.01,0.05,0.10,0.15,...,and 1.0 respectively.For a fixedr0/rch,we find the criticalpof failure propagation by changingpin steps of 0.005.The critical curves are obtained by connecting the different critical points.Three important findings can be observed from figure 8.(1)The region above the unstable critical curve covers that of catastrophic focusing and shrinks the region for stable laser propagation further;and because of this, the behaviors of bright and dark solitary waves will never occur and the laser beam can only show periodic defocusing oscillation, propagation with a constant spot radius or periodic focusing oscillation.(2)Whenr0/rchis less than a certain value(the intersection point of the unstable critical curve and the curve of constant propagation) for a fixedd, even the usual constant propagation and periodic focusing oscillation cannot exist and the laser beam can only show periodic defocusing oscillation;in other words,and it is only above this lowest value forr0/rchthat the laser beam can propagate with constant spot radius.(3) With the increase ofd, the lowestr0/rchincreases apparently, and the unstable region above the unstable critical curve shrinks significantly,which enlarges the area for stable propagation.

    In the previous theoretical part it was revealed thata0has an obvious influence on the laser propagation parameter region.In order to further verify and visually show the influence ofa0on laser propagation, we make numerical simulations withr0/rch=0.4 andp=0.29, marked with the empty circle in figure 4 for differenta0, and show the results in figure 9.We can see that whena0=0.2, the propagation behavior of the laser is periodic focusing oscillation.Whena0=0.3, the laser propagates with a constant spot radius.Whena0=0.4 and 0.5, propagation takes the form of periodic defocusing oscillations.This is exactly consistent with the theoretical prediction in figure 4.In section 3, figure 5(a)revealed that constant propagation occurs only in a narrow ranger0/rch?(0.24, 0.38) for the given plasma backgroundrch=30 μm andn0=0.02nc.To verify this result, numerical simulation was performed for some parameter points along two linesr0/rch=0.23 andr0/rch=0.4 on the two sides of the range (0.24,0.38) and also for some parameters of the constant propagation shown in figure 5(a).Atr0/rch=0.23,a0is increased from 0.05 to 0.50 in steps of 0.05.The simulation shows periodic defocusing oscillations.Atr0/rch=0.4,a0is increased from 0.03 to 0.30 in steps of 0.03.The numerical simulation gives periodic focusing oscillations.Constant propagation is also observed for the selected parameter points (r0/rch=0.36,a0=0.23), (0.37,0.16)and(0.38,0.06)on the constant propagation curve.The numerical simulation results confirm the theoretical results in figure 5(a).Due space limitations,the results of the numerical simulation are not presented in this paper.

    According to the above analysis,the numerical simulation successfully verifies the theoretical results and finds an unstable critical curve above which the laser beam fails to propagate stably.The unstable critical curve amends the parameter space for laser propagation and gives rise to some interesting results even for the parabolic channel.The modified parameter space varies significantly with the channel power-law exponentd.

    5.Summary

    Containing the effect of ponderomotive self-channeling, the propagation theory of a laser beam traveling in a power-law channel is developed.Ponderomotive self-channeling induces a large variation of the parameter space(p,r0/rch),which means a much more complex laser behavior.In the variation, more attention should be paid to catastrophic focusing because it occupies a relatively large region and shrinks the domain of laser stable propagation in the parameter space.Increasing the power-law exponent can decrease the catastrophic region,which is beneficial for laser propagation.Increasing the initial laser amplitude, however, may enlarge the catastrophic region moderately and induce obvious variations in the parameter regions for other stable propagation behaviors.The effects of the initial amplitude on the propagation behavior are also presented and studied in another parameter space(a0,r0/rch)for a given plasma background.By directly solving the wave equation with the finite difference method, the propagation dynamics and stability are also studied numerically.The numerical simulation confirms the main theoretical results containing the constant propagation, periodic focusing and defocusing oscillation.Meanwhile,the numerical simulation reveals a radial unstable critical curve in the parameter space(p,r0/rch), above which stable propagation fails.The unstable critical curve brings some interesting findings, for example,it prevents the formation of bright and dark solitary waves and the unstable region above the unstable critical curve covers that of catastrophic focusing and further shrinks the stable domain for the laser beam; this indicates that the laser beam can propagate with constant spot radius only above a certain low value ofr0/rch.The unstable critical curve is sensitive to the channel power-law exponent and this results in the apparent decrease of the unstable region above the unstable critical curve and the increase in the lowestr0/rchvalue for constant propagation with this exponent.It can be seen that the radial instability corresponding to the unstable critical curve is one of the important discoveries in our work, and a detailed theoretical and numerical investigation is necessary.However, this is a systematic work with many calculations and much analysis,and we will address it further in the near future.To summarize,the obtained results recover our previous knowledge and reveal some valuable information about conditions for stable laser propagation, it may therefore help bring insight into the future studies related to laser propagation.

    Acknowledgments

    This work is supported by National Natural Science Foundation of China (Nos.11 765 017, 11 865 014, 12 047 574,11 847 304,11 764 039 and 12 165 018), and by the Scientific Research Project of Gansu Higher Education(No.2019B-034).

    Appendix: properties of the Sagdeev potential

    Under the appropriate condition 00 andd>1,we investigateV(rs)=0 andV′ (rs)=0:

    (1) The properties of the roots ofV(rs)=0.Firstly,rs=1 is one of the roots.Secondly, we rewriteV(rs)=0 in the formthen set the expression on the left-hand side asf(rs) and that on the righthand side asg(rs).Obviously, the curvedecreases first and reaches a minimum value then increases finally whenrsvaries from 0 to infinity (this is directly concluded by analyzing the expression).Meanwhile, the curvedecreases quickly, then slowly and finally varies linearly.These are the geometric properties of the two curves.For illustration, an example is shown in figure 10(a).Considering the geometric properties and the asymptotic behaviors∣′f(rs)∣ < ∣g′(rs)∣forrs→0 andf′(rs)>g′(rs)forrs→∞; making an elementary geometric analysis,one can find that there are at most three intersection points between the two curvesf(rs) andg(rs) in the range ofrs>0.This means thatV(rs)=0 has three positive roots at most.In summary, one root ofV(rs)=0 is equal to 1, and there are at most two other positive roots.

    (2) The properties of the roots ofV′ (rs)=0.Differentiating equation (14) once, one obtains(1?p) rs2.Considering the conditions 00 andd>1, and making a similar geometric analysis (an example is shown in figure 10(b)),one may conclude that the number of intersection points of the curvesand(1?p) rs2may be zero, one or two.That is,V′(rs)=0may have no positive root,one positive root or two positive roots.We have done some numerical solutions ofV(rs)=0 andV′ (rs)=0.The results are all consistent with the above properties.

    亚洲成av人片在线播放无| 久久精品91蜜桃| 国产单亲对白刺激| 卡戴珊不雅视频在线播放| 麻豆国产97在线/欧美| 最新中文字幕久久久久| 神马国产精品三级电影在线观看| 人人妻人人澡人人爽人人夜夜 | 精品人妻视频免费看| 亚洲成人av在线免费| 不卡视频在线观看欧美| 日韩欧美国产在线观看| 韩国av在线不卡| 精品免费久久久久久久清纯| 99在线人妻在线中文字幕| 91aial.com中文字幕在线观看| 中文亚洲av片在线观看爽| 少妇高潮的动态图| 亚洲国产日韩欧美精品在线观看| 国产色婷婷99| 国产色爽女视频免费观看| 欧美3d第一页| 国产精品美女特级片免费视频播放器| 一区二区三区乱码不卡18| 青春草亚洲视频在线观看| 国产男人的电影天堂91| 精品一区二区免费观看| 在线观看一区二区三区| 国产真实伦视频高清在线观看| 免费看a级黄色片| 边亲边吃奶的免费视频| 99久久九九国产精品国产免费| 亚洲欧美日韩东京热| 日本熟妇午夜| 中文字幕久久专区| 高清毛片免费看| 亚洲精品456在线播放app| 国产69精品久久久久777片| 国内精品宾馆在线| 69av精品久久久久久| 搡老妇女老女人老熟妇| 在线观看美女被高潮喷水网站| 深爱激情五月婷婷| 18禁在线播放成人免费| 18禁在线播放成人免费| 校园人妻丝袜中文字幕| 亚洲av熟女| 久久韩国三级中文字幕| 岛国在线免费视频观看| 女人十人毛片免费观看3o分钟| 天美传媒精品一区二区| АⅤ资源中文在线天堂| 久久久亚洲精品成人影院| 免费观看在线日韩| 亚洲av日韩在线播放| 国产极品精品免费视频能看的| 国产av码专区亚洲av| 免费av不卡在线播放| 天天一区二区日本电影三级| 欧美高清成人免费视频www| 国产日韩欧美在线精品| 又黄又爽又刺激的免费视频.| 国产大屁股一区二区在线视频| 又爽又黄a免费视频| 亚洲久久久久久中文字幕| 女的被弄到高潮叫床怎么办| 亚洲不卡免费看| 国产成人精品一,二区| 日韩一区二区视频免费看| 久久久久久久久中文| 啦啦啦韩国在线观看视频| 亚洲丝袜综合中文字幕| 久久久精品大字幕| 中文在线观看免费www的网站| 国产成人a区在线观看| 国产淫语在线视频| 少妇的逼水好多| 大又大粗又爽又黄少妇毛片口| 一级毛片aaaaaa免费看小| 日日摸夜夜添夜夜添av毛片| 亚洲国产日韩欧美精品在线观看| 久久久欧美国产精品| 国产精品99久久久久久久久| av福利片在线观看| 亚洲电影在线观看av| 成人毛片a级毛片在线播放| 午夜日本视频在线| 欧美日韩在线观看h| 26uuu在线亚洲综合色| 成人国产麻豆网| 七月丁香在线播放| 久久精品影院6| 色综合站精品国产| 一级毛片久久久久久久久女| 嫩草影院新地址| 亚洲精品国产av成人精品| 少妇人妻精品综合一区二区| 亚洲av不卡在线观看| av黄色大香蕉| 久久人妻av系列| 级片在线观看| 秋霞伦理黄片| 成人午夜精彩视频在线观看| 欧美日韩国产亚洲二区| 中国美白少妇内射xxxbb| 18+在线观看网站| 内射极品少妇av片p| 少妇熟女aⅴ在线视频| 国产精品1区2区在线观看.| 美女脱内裤让男人舔精品视频| 免费观看在线日韩| 日本-黄色视频高清免费观看| 午夜福利网站1000一区二区三区| 在线观看av片永久免费下载| 亚洲自拍偷在线| 看非洲黑人一级黄片| 国产 一区精品| 搡女人真爽免费视频火全软件| 色吧在线观看| 亚洲国产欧美在线一区| 丰满乱子伦码专区| 深爱激情五月婷婷| 插阴视频在线观看视频| 日本一二三区视频观看| 国产人妻一区二区三区在| 麻豆av噜噜一区二区三区| 久热久热在线精品观看| 乱码一卡2卡4卡精品| 成年av动漫网址| 国产亚洲最大av| 青春草视频在线免费观看| 你懂的网址亚洲精品在线观看 | 亚洲av不卡在线观看| 亚洲av电影在线观看一区二区三区 | 女人久久www免费人成看片 | 免费电影在线观看免费观看| 国产一区亚洲一区在线观看| 成人特级av手机在线观看| 国产黄色小视频在线观看| 午夜老司机福利剧场| 免费看日本二区| 亚洲国产精品国产精品| 久久这里有精品视频免费| 精品酒店卫生间| 国内精品美女久久久久久| 汤姆久久久久久久影院中文字幕 | 国产一级毛片在线| 国产探花在线观看一区二区| 免费看光身美女| 精品酒店卫生间| 免费观看在线日韩| 在线观看av片永久免费下载| 亚洲av成人精品一二三区| 成人特级av手机在线观看| 一级毛片我不卡| 在线免费观看不下载黄p国产| 欧美另类亚洲清纯唯美| 亚洲在线自拍视频| 看片在线看免费视频| 在线a可以看的网站| ponron亚洲| 日韩欧美精品v在线| 亚洲伊人久久精品综合 | 少妇猛男粗大的猛烈进出视频 | .国产精品久久| 国产高潮美女av| 深夜a级毛片| 长腿黑丝高跟| 校园人妻丝袜中文字幕| 亚洲最大成人手机在线| 国产一级毛片在线| 中文字幕久久专区| 搞女人的毛片| 一级黄色大片毛片| 人人妻人人澡欧美一区二区| 爱豆传媒免费全集在线观看| 亚洲av一区综合| 最近手机中文字幕大全| 亚洲成人久久爱视频| 国模一区二区三区四区视频| 国产成人aa在线观看| 中文字幕免费在线视频6| 国产激情偷乱视频一区二区| 网址你懂的国产日韩在线| 永久网站在线| 久久精品人妻少妇| 蜜臀久久99精品久久宅男| 精品国产三级普通话版| 在线免费观看的www视频| 18禁在线播放成人免费| 一区二区三区乱码不卡18| 国产乱来视频区| 日本一本二区三区精品| 高清视频免费观看一区二区 | 亚洲精品乱码久久久v下载方式| 亚洲在线观看片| 亚洲性久久影院| 白带黄色成豆腐渣| 欧美一区二区精品小视频在线| 欧美日韩在线观看h| 日日撸夜夜添| 深爱激情五月婷婷| 男人的好看免费观看在线视频| 水蜜桃什么品种好| 亚洲美女搞黄在线观看| 观看免费一级毛片| 日本熟妇午夜| 久久热精品热| 国产精品一及| 深爱激情五月婷婷| 在现免费观看毛片| 亚洲久久久久久中文字幕| 亚洲欧美成人精品一区二区| 亚洲第一区二区三区不卡| 日韩精品有码人妻一区| 色噜噜av男人的天堂激情| 九九热线精品视视频播放| 国产亚洲av嫩草精品影院| 亚洲图色成人| 国产伦一二天堂av在线观看| 国产爱豆传媒在线观看| 国产 一区精品| 色综合站精品国产| 亚洲婷婷狠狠爱综合网| 桃色一区二区三区在线观看| 丝袜喷水一区| 小说图片视频综合网站| 日韩中字成人| or卡值多少钱| 中文乱码字字幕精品一区二区三区 | 日韩中字成人| 国产精品三级大全| 少妇裸体淫交视频免费看高清| 春色校园在线视频观看| 晚上一个人看的免费电影| 国产淫语在线视频| 高清午夜精品一区二区三区| 亚洲国产精品sss在线观看| 色吧在线观看| 日本爱情动作片www.在线观看| 99在线视频只有这里精品首页| 亚洲人成网站在线播| a级一级毛片免费在线观看| 午夜福利在线观看免费完整高清在| 成人国产麻豆网| 舔av片在线| 国产人妻一区二区三区在| 国产成人a区在线观看| 成年女人看的毛片在线观看| 国产精品av视频在线免费观看| 少妇猛男粗大的猛烈进出视频 | 国产又黄又爽又无遮挡在线| 日本三级黄在线观看| 久久久久久久久中文| 人体艺术视频欧美日本| 2021天堂中文幕一二区在线观| 国产欧美另类精品又又久久亚洲欧美| 欧美xxxx黑人xx丫x性爽| 成年av动漫网址| 麻豆成人av视频| 国产亚洲5aaaaa淫片| 九色成人免费人妻av| 欧美成人一区二区免费高清观看| 欧美丝袜亚洲另类| 成人午夜精彩视频在线观看| 国产爱豆传媒在线观看| 九九热线精品视视频播放| 精品人妻熟女av久视频| 国产精品麻豆人妻色哟哟久久 | 国产精品.久久久| 国产亚洲5aaaaa淫片| 一级毛片电影观看 | av线在线观看网站| .国产精品久久| 日韩制服骚丝袜av| 成人午夜高清在线视频| 男人狂女人下面高潮的视频| 亚洲,欧美,日韩| 日本与韩国留学比较| 久久久精品欧美日韩精品| 亚洲欧美清纯卡通| 国产av一区在线观看免费| 国产精品国产高清国产av| 亚洲国产精品sss在线观看| 一夜夜www| 91精品国产九色| 人妻制服诱惑在线中文字幕| eeuss影院久久| 91久久精品国产一区二区成人| 国产v大片淫在线免费观看| 国产精品1区2区在线观看.| 日韩三级伦理在线观看| 99久久无色码亚洲精品果冻| 亚洲久久久久久中文字幕| 久久久久久久久久成人| 蜜臀久久99精品久久宅男| 91久久精品国产一区二区三区| 简卡轻食公司| 可以在线观看毛片的网站| 我的老师免费观看完整版| 亚洲国产欧美人成| 国产人妻一区二区三区在| 免费av观看视频| 国产精品无大码| 深爱激情五月婷婷| 久久精品国产亚洲网站| 成人综合一区亚洲| 亚洲电影在线观看av| 国产一区亚洲一区在线观看| 中文乱码字字幕精品一区二区三区 | 亚洲自偷自拍三级| 老司机影院成人| 欧美潮喷喷水| 亚洲精品日韩在线中文字幕| 久久久久久久午夜电影| 国产黄色视频一区二区在线观看 | 国产精品国产三级国产av玫瑰| 卡戴珊不雅视频在线播放| 久热久热在线精品观看| 亚洲,欧美,日韩| 日韩欧美精品v在线| 六月丁香七月| 国产久久久一区二区三区| 人人妻人人澡人人爽人人夜夜 | 国语自产精品视频在线第100页| 26uuu在线亚洲综合色| 丝袜美腿在线中文| 麻豆乱淫一区二区| 91精品一卡2卡3卡4卡| 亚洲激情五月婷婷啪啪| 欧美一区二区国产精品久久精品| 国产黄色视频一区二区在线观看 | 99在线人妻在线中文字幕| 最近2019中文字幕mv第一页| 亚洲av中文字字幕乱码综合| 亚洲真实伦在线观看| 中文字幕熟女人妻在线| 三级毛片av免费| 国产伦在线观看视频一区| 黄片wwwwww| 嫩草影院入口| 国产综合懂色| 久久久亚洲精品成人影院| 两个人视频免费观看高清| 欧美成人免费av一区二区三区| 成年版毛片免费区| 美女黄网站色视频| 简卡轻食公司| 日日摸夜夜添夜夜添av毛片| 精品人妻视频免费看| 色播亚洲综合网| 如何舔出高潮| 精品酒店卫生间| 日韩制服骚丝袜av| 老司机影院成人| 亚洲av中文字字幕乱码综合| 国产一区二区三区av在线| 国产成人免费观看mmmm| 色播亚洲综合网| 国产免费又黄又爽又色| 亚洲av中文av极速乱| 国产黄片视频在线免费观看| 午夜久久久久精精品| 我的老师免费观看完整版| 亚洲最大成人av| 99久国产av精品| 免费av不卡在线播放| 麻豆国产97在线/欧美| 国产视频内射| 国产毛片a区久久久久| 七月丁香在线播放| 成人性生交大片免费视频hd| 天天一区二区日本电影三级| 久久这里只有精品中国| 一本一本综合久久| 国产成人午夜福利电影在线观看| 国产一区二区在线观看日韩| 国产视频首页在线观看| 国产精品久久久久久精品电影小说 | 青春草国产在线视频| 久久久久国产网址| av免费在线看不卡| 亚洲国产色片| av在线亚洲专区| 成人国产麻豆网| 国产成人aa在线观看| 精品少妇黑人巨大在线播放 | 国产精品伦人一区二区| 亚洲经典国产精华液单| 爱豆传媒免费全集在线观看| 久久精品国产99精品国产亚洲性色| 国产老妇女一区| 国产成人精品婷婷| 免费看av在线观看网站| 女人久久www免费人成看片 | 欧美高清性xxxxhd video| 欧美三级亚洲精品| 一级毛片电影观看 | 免费看光身美女| 高清日韩中文字幕在线| 精品人妻熟女av久视频| 国产午夜精品论理片| 久久婷婷人人爽人人干人人爱| 欧美xxxx黑人xx丫x性爽| 国产精品电影一区二区三区| 国产乱人视频| 少妇猛男粗大的猛烈进出视频 | 有码 亚洲区| 免费看av在线观看网站| av福利片在线观看| 国产一区二区亚洲精品在线观看| 一级av片app| 亚洲成人av在线免费| 一级毛片电影观看 | 一个人看视频在线观看www免费| 日韩大片免费观看网站 | 久久久久久久亚洲中文字幕| 国产三级在线视频| 美女cb高潮喷水在线观看| 国产精品人妻久久久久久| 3wmmmm亚洲av在线观看| 一本久久精品| 国产精品乱码一区二三区的特点| 最近最新中文字幕免费大全7| 欧美丝袜亚洲另类| 99视频精品全部免费 在线| 久久久久久国产a免费观看| 成年av动漫网址| 中文字幕av成人在线电影| 卡戴珊不雅视频在线播放| 免费观看精品视频网站| 男女下面进入的视频免费午夜| 国产精品久久久久久精品电影小说 | 成人亚洲精品av一区二区| 在线天堂最新版资源| 中国国产av一级| 真实男女啪啪啪动态图| 久久国产乱子免费精品| 国产片特级美女逼逼视频| 精品酒店卫生间| 国产精品爽爽va在线观看网站| 18禁动态无遮挡网站| 午夜福利高清视频| 精品久久久久久久末码| 国产精品不卡视频一区二区| 最近最新中文字幕免费大全7| 国产精品伦人一区二区| 男女啪啪激烈高潮av片| 乱人视频在线观看| videos熟女内射| 国产真实伦视频高清在线观看| 好男人在线观看高清免费视频| 精品不卡国产一区二区三区| 国产黄片视频在线免费观看| 天堂av国产一区二区熟女人妻| 九九久久精品国产亚洲av麻豆| 久久99热6这里只有精品| 欧美三级亚洲精品| 欧美zozozo另类| 亚洲欧美成人精品一区二区| 免费电影在线观看免费观看| 午夜激情福利司机影院| 综合色av麻豆| 桃色一区二区三区在线观看| 永久网站在线| 99国产精品一区二区蜜桃av| 国产成人午夜福利电影在线观看| 国产高清视频在线观看网站| 久久精品国产亚洲av天美| 亚洲精品自拍成人| 美女cb高潮喷水在线观看| 亚洲国产精品合色在线| 欧美性猛交黑人性爽| 国产精品永久免费网站| 自拍偷自拍亚洲精品老妇| 日本欧美国产在线视频| 国语自产精品视频在线第100页| 欧美一区二区国产精品久久精品| 丝袜喷水一区| 免费av毛片视频| 晚上一个人看的免费电影| 久久久午夜欧美精品| 男人和女人高潮做爰伦理| 国产黄片美女视频| 日韩成人av中文字幕在线观看| 91精品一卡2卡3卡4卡| 国产精品三级大全| 日本一本二区三区精品| 成人国产麻豆网| 建设人人有责人人尽责人人享有的 | 精品久久久久久久人妻蜜臀av| 亚洲国产精品专区欧美| 美女高潮的动态| 中文字幕av成人在线电影| 最近手机中文字幕大全| 国产精品国产高清国产av| 偷拍熟女少妇极品色| 少妇的逼水好多| 国产一区二区在线av高清观看| 国产亚洲av片在线观看秒播厂 | 亚洲色图av天堂| 综合色av麻豆| 久久人妻av系列| 亚洲av成人精品一二三区| 久久精品国产99精品国产亚洲性色| 亚洲精品亚洲一区二区| 国产精品嫩草影院av在线观看| 一个人看的www免费观看视频| 69人妻影院| 久久久久久久亚洲中文字幕| 男人舔女人下体高潮全视频| 国产单亲对白刺激| 日韩精品青青久久久久久| 男人的好看免费观看在线视频| 天天一区二区日本电影三级| 一级二级三级毛片免费看| 午夜视频国产福利| 天堂√8在线中文| 久久人妻av系列| 久久久久久伊人网av| 晚上一个人看的免费电影| 国产亚洲5aaaaa淫片| 亚洲成人久久爱视频| av.在线天堂| 亚洲av熟女| 日本熟妇午夜| 日韩成人av中文字幕在线观看| 欧美+日韩+精品| 激情 狠狠 欧美| 午夜免费男女啪啪视频观看| 少妇的逼好多水| 精品久久久噜噜| 少妇的逼好多水| 国产探花在线观看一区二区| 日韩三级伦理在线观看| 国产综合懂色| 七月丁香在线播放| 91精品国产九色| 成人三级黄色视频| 日本黄色片子视频| 欧美日本亚洲视频在线播放| 99热6这里只有精品| 边亲边吃奶的免费视频| 国产在视频线精品| 国产极品精品免费视频能看的| 97热精品久久久久久| 高清午夜精品一区二区三区| 毛片一级片免费看久久久久| 99久久九九国产精品国产免费| 日本免费a在线| 亚洲,欧美,日韩| 我要搜黄色片| 色5月婷婷丁香| 99九九线精品视频在线观看视频| 免费av观看视频| 亚洲最大成人av| 欧美高清成人免费视频www| 中国国产av一级| 亚洲无线观看免费| 我要看日韩黄色一级片| 91精品一卡2卡3卡4卡| 狂野欧美激情性xxxx在线观看| 18+在线观看网站| 日韩强制内射视频| 亚洲,欧美,日韩| 女人十人毛片免费观看3o分钟| 亚洲av二区三区四区| 一本久久精品| 亚洲第一区二区三区不卡| 99在线视频只有这里精品首页| 三级国产精品欧美在线观看| av女优亚洲男人天堂| 午夜福利网站1000一区二区三区| 国产淫片久久久久久久久| 久久精品人妻少妇| 中文字幕精品亚洲无线码一区| 国产亚洲5aaaaa淫片| 久热久热在线精品观看| av国产免费在线观看| 别揉我奶头 嗯啊视频| 免费av毛片视频| 国产一区有黄有色的免费视频 | 亚洲电影在线观看av| 国产精品嫩草影院av在线观看| 欧美一区二区精品小视频在线| 国产美女午夜福利| 午夜福利在线在线| 99热全是精品| 国产午夜精品久久久久久一区二区三区| 国产一级毛片在线| 男女下面进入的视频免费午夜| 精品久久久久久久久久久久久| 日日撸夜夜添| 中文字幕久久专区| 久久久久久大精品| 欧美人与善性xxx| 99久久无色码亚洲精品果冻| 人人妻人人澡人人爽人人夜夜 | 搡女人真爽免费视频火全软件| 国产av一区在线观看免费| 国产精品国产三级国产av玫瑰| 国产精品电影一区二区三区| or卡值多少钱| 免费看光身美女| 国产精品久久视频播放| 九九在线视频观看精品| av女优亚洲男人天堂| 欧美最新免费一区二区三区| 最近手机中文字幕大全| 纵有疾风起免费观看全集完整版 | av在线播放精品| 欧美最新免费一区二区三区| 国内精品宾馆在线| 国产精品综合久久久久久久免费| 久久精品夜色国产|