Self-compression of soliton-like laser pulses in the process of self-focusing

We study the possibility of efficient self-compression of femtosecond laser pulses in nonlinear media with anomalous dispersion of group velocity during the self-focusing of wave packets with a power several times greater than the critical self-focusing power. The results of qualitative analysis of the evolution of three-dimensional wave packets with the quasi-soliton field distribution are confirmed by the computer simulation. The simulation proves that the considered regime of compression of high-power laser pulses with initial durations of about ten optical cycles is stable relative to filamentation instability due to the influence of the nonlinear dispersion. We demonstrate the possibility of self-compression of laser pulses at a multi-millijoule energy level and up to one optical cycle with an energy efficiency of more then 50%.


I. INTRODUCTION
The generation of high-energy few-cycle optical pulses is a challenge in contemporary laser physics that has important implications for high-field science and as well for many other extreme light applications [1]. Conventional laser systems, based on broadband active media, e.g., such as Ti:sapphire crystals and/or parametric amplification technology, are able to provide pulses of ultrashort durations with high enough energies. However, to get high-energy pulses of few-cycle or even single cycle duration require additional spectral broadening with subsequent compression technique. For example, usual way is to broaden spectrum of the pulse in a high pressure gas and then to compress it by using gratings or chirp mirrors. Alternative but more refined way is to use selfcompression mode employing Kerr nonlinearity [2,3], for high-power pulses ionization [4][5][6] and relativistic [7] nonlinearities. It should be noted two important features on which this way relies. First, anomalous dispersion is required to realize self-compression mode. Second, waveguiding of pulse propagation is used for long interaction distances, which imposes some restrictions on the pulse energy. It is worth to emphasize that at least for powers about the critical power for self-focusing P cr when single filamentation occurs, self-guided pulses can be stable and their self-compressing favorably used for obtaining pulses of extremely short durations [8][9][10]. Of course, selfcompression of optical pulses in anomalous group velocity dispersion media with the help of transversely cummulating laser energy, i.e., self-focusing, looks very attractive idea for getting high-energy few-cycle pulses. Nevertheless, for high energy pulses with peak power much exceeding the threshold of self-focusing, the filamentation instability is the main obstacle that completely prevents using free space propagation in nonlinear bulk media.
It should be noted that most of the media in the visible frequency range are characterized by normal dispersion * sksa@ufp.appl.sci-nnov.ru of group velocities. However, in recent years, considerable success has been achieved in the field of generation of high-energy laser pulses in the mean IR range. In this frequency range, most of the media are characterized by anomalous group velocity dispersion. The interest in ultrashort high-energy laser pulses in the 2-8-µm range is directly related to their specific features connected with generation of higher harmonics [11,12], filamentation [13,14], and particle acceleration.
Numerical simulation and experiments demonstrated that filamentation of optical radiation in the mean IR range follows a different scenario [8,15,16]. A peculiar feature of filamentation dynamics is possible formation of a set of soliton structures in the optical field, which have rather short durations (two or three field periods) [10]. This regime of laser pulse shortening is limited by the powers about of the critical self-focusing power. This take place due to development of the filamentation instability, which leads to fast decomposition of the beam to many filaments, and the final loss of coherence makes the sources of such pulses inapplicable in practice.
In this work, basing on the results described in [17], we present detailed analytical and numerical studies, which demonstrate that nonlinear dispersion of the medium (dependence of the group velocity of the wave packet on its intensity) leads to stabilization of filamentation instability for ultrashort pulses. It is shown that at anomalous medium dispersion, the filamentation-free regime of laser pulse self-focusing allows one to realize self-compression up to one optical cycle for soliton-like pulses having powers much larger than critical one for the self-focusing.
The structure of this work is as follows. In Section II, the basic equation is formulated for describing self-action of ultrashort laser pulses comprising several optical cycles in a medium with a Kerr-type inertia-free nonlinearity. Sections III and IV present a qualitative study of the method for adiabatic compression of the laser pulse duration. In Section V, stability of the radiation self-focusing is studied in relation to the filamentation instability. In the last part, SectionVI, the results of numerical simulation are presented, and the initial parameters of the laser arXiv:1702.01965v1 [physics.optics] 7 Feb 2017 pulse for optimal self-compression of the wave packet are determined.

II. BASIC EQUATIONS
To describe adequately the spatio-temporal evolution of ultrashort, circularly polarized laser pulses (E = E(x 0 + iy 0 ), where E x and E y are the corresponding components of the electric-field intensity) in a medium with cubic nonlinearity at the self-focusing process. Let us turn now directly to the wave equation Here, P nl is the nonlinear medium response, and ε is the linear dielectric permittivity, which satisfies the fundamental Kramers Kronig relation where ε r and ε i are the real and imaginary parts of the dielectric permittivity ε. Let us apply Eq. (2) to weakly absorbing media, i.e., assume that one can neglect the imaginary part of the dielectric permittivity within the frequency range being of interest for us. Let's assume that the weak-absorption region spreads in a wide frequency range from ω 1 to ω 2 , and consider the frequencies ω, such that ω 1 ω ω 2 .
Then, we will limit our consideration to the case, where ε 0 ω 2 D /ω 2 bω 2 . This condition is fulfilled well for the majority of media in the mean IR range. As a result, the dielectric permittivity will take on the following form: In the case of a nonresonance medium with the instantaneous nonlinear response P nl = χ (3) |E| 2 E (here, χ (3) is the cubic nonlinear susceptibility), wave equation (1) will acquire the following form [18,19]: This equation describes the spatio-temporal evolution of the field intensity E of the laser pulse allowing for all considered important physical factors, such as medium dispersion, beam diffraction, Kerr nonlinearity, selfsteepening of the wave packet, and self-focusing of the transverse distribution. For detailed analysis of a dynamic problem, it is convenient to use the evolution equation for the field in the simplest form of a reduced wave equation. Let's suppose that the spatio-temporal structure of the wave field varies smoothly in the process of one-way pulse propagation, i.e., let's neglect the reflection effects. In the quasimonochromatic case, this approach corresponds to passing over to an envelope equation. Using the approximation of one-way propagation of the wave field along the z axis, i.e., supposing smallness of |c∂ z E + ∂ t E| |c∂ z E − ∂ t E|, we will represent the equation in dimensionless variables in the following form [18][19][20][21]: Here In the case of the monochromatic wave packet Eq. (6) yields easily an equation, which generalizes the nonlinear Schrödinger equation (NSE) for the envelope [22] i ∂Ψ ∂z Here, k z = −1/ω,τ = τ − z/ω 2 0 is the time in the accompanying system of coordinates, ω 0 = 1 is carrier frequency,T = 1 − i ∂ ∂τ is the operator correcting the approach of slowly changing amplitudes [23], and For the field distributions u(z, r ⊥ , τ ), which are localized in time and space, the following values are constant: Formula (8a) is responsible for the absence of the zeroth harmonic in the field distribution. The integral (8b) determines the conservation of the "number of quanta". Formula (8c) is a Hamiltonian. Using the "Hamiltonian" nature of Eq. (6), one can obtain the relationship which describes the variation of the efficient transverse width of the wave field in the process of self-focusing, In what follows, we will consider the initial distributions of laser pulses with the negative Hamiltonian H full < 0, that are definitely collapsing in the transverse direction along a finite propagation path (see Eq. (9)). In the case, where spatial effects are insignificant (∆ ⊥ ≡ 0), there exists a class of stable soliton solutions [24]. Wave solitons of Eq. (6) can be represented by the two-parameter family of solutions having the form where ω s is the characteristic carrier frequency, γ is the parameter determining the group velocity of the soliton, and ξ = ω s (τ − γz). The soliton amplitude G(ξ) and the nonlinear phase φ(ξ) obey the following equations: where F (G 2 ) = G 2 3/2(1 + δ 2 ) − (4 − 5G 2 )/4(1 − G 2 ) 2 , G m is the maximum amplitude of the soliton, and ξ 0 is the integration constant corresponding to the position of the maximum of the soliton amplitude. As seen from Eq. (12b), solutions for the soliton amplitude G(ξ) depend only on the parameter δ 2 = 1/(ω 2 s γ) − 1 and exist at 0 ≤ δ ≤ δ cr ≡ 1/8. The critical value gives the single-cycle soliton [24]. An important feature of the wave solitons is a semi-bounded spectrum of their admissible solutions, i.e., the presence of a boundary solution corresponding to the limiting soliton with the minimum possible pulse duration and, hence, the maximum possible amplitude. It should be noted that the existence of the limiting soliton is defined by the constraint ∞ −∞ udτ = 0, which is one of the integrals of Eq. (6). At δ = δ cr the shortest duration is equal to τ * s = 2.31ω −1 s . For the sake of comparison, Figs. 1(a,c) present exact soliton solutions for two different values of the parameter δ: (a) for δ = 0.06 and (b) for δ = 0.3. One can see that as the value of δ increases, the duration of the soliton decreases. Figures 1(b,d) show the spectral intensities for various δ. It follows from the figure that as δ increases, the soliton spectrum becomes wide and asymmetric due to a strong frequency modulation in the pulse (12a). Now let us demonsrate the connection between the solution of Eq. (12b) with the well-known envelope solitons, which exist within the framework of NSE and its generalizations. Passing over to long quasi-monochromatic wave packets corresponds to the case of small values of the amplitudes G 1. Keeping the terms of the order of G 4 in Eq. (12b), we obtain the solitary solution for the envelope Typical distributions of the field and spectral intensity of the soliton at small δ = 0.06 are shown in Fig. 1(a,b). Keeping the terms of the next order of smallness, G 6 , in Eq. (12b), we obtain the NSE solution allowing for nonlinear dispersion (amplitude dependence of the group velocity) A distinctive feature of this solution is the presence of a sufficiently strong frequency modulation in a laser pulse. Typical distributions of the field and spectral intensity of the soliton at δ = 0.3 are shown in Fig. 1(c,d).

III. QUALITATIVE ANALYSIS OF WAVE PACKET SELF-COMPRESSION
Now, let us turn to the problem of self-compression of wave packets in the self-focusing process. Equation (6) is much more complicated than the corresponding equation for the quasi-monochromatic radiation. Therefore, in order to obtain analytical relationships, first we turn to simplify Eq. (7), where we neglect such important effects as high-order dispersion, nonlinear dispersion (amplitude dependence of the group velocity) and spatio-temporal self-focusing (T = 1) [25] i ∂Ψ ∂z Here, the second term describes anomalous media dispersion, the third term describes pulse diffraction, and the fourth term allows for the cubic nonlinearity of the medium. In this equation, the "number of quanta" and the Hamiltonian are also saved for localized distributions: Note that the issue of pulse self-focusing in media with anomalous and normal dispersion is addressed in many papers [26][27][28][29][30]. When a wave packet propagates in a nonlinear medium, it is affected simultaneously by dispersion and diffraction. At the same time, however, these two effects become interconnected due to the nonlinearity of the medium. This connection leads to the possibility of a spatio-temporal collapse. For equation (15), one can also write down a relationship for variation in the efficient scales (width and duration) of the wave packet (similar to Eq. (9)), which characterizes the global behavior of the system: where I q ρ 2 ⊥ = r 2 ⊥ |Ψ| 2 dτ dr ⊥ is the efficient transverse width of the wave packet, and I q ρ 2 ⊥ + τ 2 = (r 2 ⊥ + τ 2 )|Ψ| 2 dτ dr ⊥ is the total efficient scale of the wave packet. From here, as in the case of relationship (9) at the negative Hamiltonian H q < 0, one can make a certain conclusion about the collapse of not only efficient transverse width (17a), but also the total effective scale of the laser pulse (17b) in the case of its propagation in a nonlinear medium.
For a more detailed study of the evolution of the wave packet in a nonlinear medium, we turn to the variation approach [31], which allows one to find an approximate solution of Eq. (15). The essence of the method consists in finding solutions for this class of the functions Ψ = f (r; σ(z)), where the set of parameters σ(z) depends on the evolution variable and is determined basing on the solutions of the corresponding system of differential equations. This system of equations is found on the basis of requirements of minimization of the action functional, in which the integration with respect to spatial variables (r ⊥ , τ ) is performed: As a result, the "abbreviated" Lagrangiañ yields a system of equations in a usual way: Note that the main difficulty of this method consists in the choice of the class of functions. If it is not well chosen, then the resulting approximate solution will be a little like a true solution of the original equation. Let us use the variational approach in order to find the approximate solution of Eq. (15). We will seek the solution in the class of Gaussian functions, i.e., in the so-called aberration-free approximation where a ⊥ (z), a (z), α(z) and β(z) are the parameters of the function, which depend on the evolution variable z.
The advantages of the Gaussian profile of solution (21) consist in its well-localization and in the absence of singularity at center. Substituting expression Eq. (21) to the action functional (18) with Lagrangian we obtain an "abbreviated" LagrangianL Variation with respect to the variables σ = {a ⊥ , a , α, β} yields the following equations • from variation with respect to α: • from variation with respect to β: • from variation with respect to a ⊥ : • from variation with respect to a : It is seen from Eqs. (24c)-(24d) that the decrease in width and duration of the laser pulse in the process of radiation self-focusing will take place in the case, where the initial packet duration a (0) satisfies the following relationship: In the case of a great difference in the scales (a ⊥ a ), as follows from Eqs. (24c)-(24d), the compression will occur only along the transverse coordinate. In this case, when the initial width and duration are not strongly different, the scales along two coordinates will line up in the process of nonlinear dynamics [26]. The regime of spherically symmetric collapse, when the longitudinal and transverse scales are equal (a ⊥ = a ), was studied in [32]. Now, we will consider the most interesting special case of the self-action regime, where the pulse envelope have longitudinal scale being much less than the transverse ones (a a ⊥ ). In this case, equation (24d) is an equation with small coefficient at highest derivative. Its "slow" motion is At "slow" motion the soliton-like law for pulse duration is fulfilled: d 2 a /dz 2 ≈ 0. This corresponds to an adiabatic decrease in the duration of the soliton-like wave packet. Substituting Eq. (26) into Eq. (24c) we obtain: It follows from Eq. (26) that the pulse duration decreases in proportion to square of pulse width, i.e., it is inversely proportional with the field intensity at the pulse axis. Note that for oblate pulses (a a ⊥ ), the self-focusing condition is fulfilled automatically (the right-hand part of Eq. (27) should be negative), since the second term in Eq. (27), which is responsible for the medium nonlinearity, exceeds the first term describing the pulse diffraction, In what follows, we will neglect the first term in Eq. (27) for such distributions. This allows us to find the variation laws for the longitudinal and transverse pulse scales depending on z. The solution of Eq. (27) allowing for Eq. (28) can be represented in quadratures under the following initial conditions at z = 0: To solve this equation, we will use the fact that as the pulse width decreases (a ⊥ < a ⊥0 ), the second term in Eq. (29) becomes less than the first one and, correspondingly, it can be neglected. As a result, we obtain approximate solutions for the pulse width a ⊥ and duration a where a ⊥0 is the initial pulse width, and a 0 is the initial pulse duration. It should be noted that, as follows from Eqs. (30), the pulse duration a always stays smaller than the pulse width a ⊥ (a ⊥ a ), i.e., the anisotropy of the wave packet distribution is retained.
Using formulas (26) and (30b), one can evaluate the compression length z comp , at which the pulse duration will turn to zero, [a (z comp ) = 0]: Let us rewrite this formula using dimensional units, It is seen from Eq. (32) that the length of the medium, at which the duration of the soliton turns to zero, is proportional to initial duration and width of the wave packet and decreases with decreasing of media linear dispersion.

IV. ULTIMATE POSSIBILITIES OF LASER PULSE SELF-COMPRESSION IN THE PROCESS OF RADIATION SELF-FOCUSING
The qualitative analysis performed in the previous section on the basis of the NSE equation (15) shows that the pulse duration will decrease adiabatically down to the zero in the self-focusing process. The question about a minimal pulse duration at self-compression arises. Evidently, as the duration of the wave packet decreases, additional effects start to show themselves, which can limit the shortening of the pulse duration, specifically, the dependence of the group velocity of the packet on the field amplitude i|Ψ| 2 ∂ τ Ψ, dispersion of the group velocity of a higher order n>2 (1/i n n!)(∂ n k z /∂ω n )∂ n Ψ/∂τ n , where k z is the wave number, and the spatio-temporal self-focusing (see Eq. (15)). In this case, we should return to initial equation (6) for analysis of the ultimate self-compression of laser pulses in the process of radiation self-focusing. As it has been already noted, the initial distributions of the wave field at H full < 0 will undergo self-focusing in the transverse direction, since collapse condition (9) is fulfilled for such distributions.
For qualitative studies of the ultimate self-compression of laser pulses in the framework of initial equation (6), it is convenient to pass over from the laboratory system of coordinates to the system of coordinates, which collapses towards a certain point (r ⊥ = 0, z = z 0 ). Let us represent the field u(z, r ⊥ , τ ) as where the new variables are Here, ζ is the new reference scale for the evolution variable, for which the moment of singularity formation is shifted to infinity. The function ρ(z) describes the variation in the transverse width of the field. This representation of the solution of Eq. (6) allows for two processes simultaneously, namely, radiation self-focusing and formation of a characteristic "horsehoe" structure of the field distribution, which is determined by the variable θ.
As a result of applying transformations (33) and (34) to equation (6), we arrive at the following equation: (35) Transformation to the "collapsing" system of coordinates allows us, as in the case of quasimonochromatic radiation (Eq. (15)), to segregate the self-focusing process in the system and reduce the problem to studying the quasione-dimensional longitudinal evolution of the pulse with respect to the variable θ. The characteristic transverse scale of the quasi-waveguide structure in new variables is of the order of unity. Now, we represent the field in the near-axis region (η ≈ 0) as = A(ζ, θ) · 1 − η 2 /4 . Substituting this formula to Eq. (35) and setting the coefficients in front of η 0 and η 2 equal to zero, we find equations for A and ρ: When obtaining Eqs. (36b), we averaged |A| 4 with respect to the pulse shape |A| 4 = 1 the characteristic scale of the field ρ(z) is a function of z only, by assumption. Formula (36a) is valid, when the threshold for self-focusing is exceeded significantly, which is valid for oblate wave packets, as we have mentioned in the previous section.
One can see that Eq. (36a), which describes the field dynamics in the near-axis region (η ≈ 0), coincides with Eq. (6) at ∆ ⊥ ≡ 0 (in the absence of spatial effects). The second term in Eq. (36a) describes weakening of the medium dispersion role with the decreasing of the pulse width ρ(z). As it was mentioned at the end of Section II, within the framework of Eq. (36a) at a constant value of ρ, there is a family of soliton solutions (11)(12) having the shape, which is similar to that of Schrödinger soliton (13). An important difference of these solutions is the presence of a strong frequency modulation (14) in solitons at short durations.
Let's suppose now that the function ρ(ζ) varies smoothly in new variables where ∆ω is the spectral width of the laser pulse. So, the soliton parameters will be adjusted smoothly in the selffocusing process and an adiabatic increase in the soliton amplitude u ∝ A/ρ will take place as the pulse width ρ(z) will decrease Eq. (36b). The solution of Eq. (36a) can be written for A(ζ, θ) = B exp(iφ) in the following form depending on the current soliton duration τ ∝ 1/δ: Allowing for preservation of the total "number of quanta", it is convenient to rewrite the solution of Eq. (37) using this quantity. Thus, we manage to get rid of the parameter δ, which is now related via the energy I full in the laser pulse and the current pulse width ρ(z). As a result, we obtain the final system of equations for the dynamics of the wave packet: where α is a number of the order of unity. In this case, the frequency modulation in the soliton φ(θ) depending on the current duration of the soliton obeys formula (37).
One can see from formula (39a) that as the transverse pulse width ρ decreases, the soliton duration τ p decreases in the near-axis region (η ≈ 0). Thus, due to the process of pulse self-focusing, the soliton duration will decrease adiabatically. Additionally, the soliton velocity 1/γ ω 2 /ρ 2 will tend to the velocity of light c/ √ ε 0 . One can see that for the pulse width, Eq. (39b) coincides up to numerical coefficient with Eq. (27).
Here, if a Schrödinger-like soliton with δ 1 and, correspondingly, without a frequency modulation φ(θ) 0 is sent to the entrance of the nonlinear medium, then, as the pulse duration decreases, a strong frequency modulation φ(θ) 3 2 θ −∞ |A(θ )| 2 dθ will arise in the pulse, which will manifest itself in the wave packet spectrum (see Fig. . 1(d)).
As it has been noted in section III, a distinctive feature of considered wave solitons (within the framework of Eq. (6) at ∆ ⊥ ≡ 0) is the semi-bounded spectrum of their admissible solution, i.e., the presence of the bounding solution 0 ≤ δ ≤ δ cr = 1/8. Therefore, the maximum degree of self-compression will be determined by this bounding soliton, whose duration is comparable with the optical cycle. This is the difference between the final regime of laser pulse self-compression and the regime, which we considered on the basis of nonlinear Schrödinger equation (15).

V. STUDY OF THE STABILITY OF THE 3D WAVE PACKET RELATIVE TO THE FILAMENTATION INSTABILITY
The considered key idea of the laser pulse compression is an adiabatic decrease of the soliton duration in the self-focusing process. The initial amplitude of the soliton is related to its duration by the following formula: where γ is the soliton velocity, τ in p = 1/(2δ 0 ) is the initial duration, and δ 0 is the initial parameter of the soliton. The total energy in the spatio-temporal bounded quasi-soliton distribution of the field is determined by formula (38). For practical realization of the proposed quasi-soliton method of pulse self-compression, of interest are wide-aperture wave packets, in which the energy flow I s = |u| 2 dτ is completely determined by the nonlinearity and dispersion of the medium. In this connection, one faces an important problem of studying the stability of the considered regime of wave packet compression in relation to various spatio-temporal perturbations of the initial structure.
The conclusion about the character of filamentation instability of continuous radiation is usually made on the basis of analysis of the instability of a plane wave [33]. The transverse perturbations of the wave front in the interval 0 < κ ⊥ ≤ κ cr = √ 2B 0 increase exponentially ∝ exp(Γz), as the wave packet propagates along z. Here the growth rate Γ 2 = κ 2 ⊥ (2B 2 0 − κ 2 ⊥ ), where B 0 is the amplitude of the plane wave. The growth rate Γ has the maximum value Γ max = B 2 0 at κ ⊥m = B 0 . In the twodimensional case, this indicates that the wave field splits into a set of beams with a power being of the order of magnitude of the critical value for self-focusing. This conclusion is also valid for inhomogeneous wave structures with the soliton distribution along one of the coordinates [34]. A more general consideration on the basis of the NSE allowing for time dispersion also leads to the existence of an instability with the growth rate [25] In the case of the modulation frequency Ω being equal to zero (Ω = 0), formula (42) describes instability of the plane wave discussed above, where we talked about the consequences of wave evolution. For κ ⊥ = 0, the corresponding instability is known as the modulation instability. Generally, the instability has the spatio-temporal character.
However, in the case of a shorter soliton pulse and, therefore, a more intense peak power, additional nonlinear effects, such as the dependence of the group velocity on intensity, should be taken into account, which in the first approximation modifies the NSE equation into the so-called derivative nonlinear Schrödinger equation (7).
To analyze the stability of the plane wave relative to the perturbations, we will seek for the solution in the form where u 0 is the amplitude of the incident wave, Substituting formula (43) to Eq. (7), in the zeroth approximation with respect to v we obtain that φ = u 2 0 z +2u 2 0 τ . As a result, the linearized equation of the first order of smallness with respect to v will acquire the form: This equation differs from the usual equation for analysis of stability of the plane wave [33] and from more general equation with high-order dispersion terms [22] by the presence of the last term and is connected with the allowance for the nonlinear medium dispersion. We will seek for the solution of the equation in the form v = a+ib, where a, b ∝ e Γz+iΩτ −iκ ⊥ r ⊥ . As a result, we obtain a system of two homogeneous equations, which has a nontrivial solution in the case, where Ω and κ ⊥ satisfy the following dispersion relation: where K(Ω) = k z (ω 0 + Ω) − k z (ω 0 ) − Ω/v gr . As a result, we find that a drift of perturbations along the pulse takes place at a rate determined by the imaginary part of Γ in the pulse frame of reference, τ = ω 0 (t+z/v gr ). Indeed, in the structure of perturbations a, b ∝ e Γ(Ω,κ ⊥ )+iΩτ −iκ ⊥ r ⊥ along with the exponential increase in the perturbations at 0 ≤ K(Ω) + κ 2 ⊥ /(1 + Ω) ≤ √ 2u 0 (which is described by the second term in formula (45)). This is easily seen, when one considers the time dependence of v: where S(Ω) is the perturbation spectrum. It is seen from Eq. (46) that the perturbations move with the relative velocity 2u 2 0 . Consequently, the homogeneous solution persists to be unstable. However, the instability changes its type and becomes convective with group velocity ∂ Im Γ/∂Ω = −2u 2 0 . Hence, for laser pulses with a duration less than a certain value τ p < τ cr , filamentation instability has no time to develop.
Next, we estimate the critical pulse duration τ cr at which the perturbation growth is stabilized due to the drift to the rear pulse front where its amplification becomes negligible. Dangerous perturbations slip down to the rear side of the pulse, because their group velocity is less than the soliton one. The "slipping" length z = z * at which the perturbation v shifts by half the pulse duration is determined as 2u 2 0 z * τ p /2. This length should be smaller than ln(δΨ/u 0 ) ≈ 15 . . . 20 of maximal perturbation growth length 1/Γ max = 1/u 2 0 . As the result, we obtain the following inequality on the critical length (in dimensional units): In other words, if the pulse duration is less than about ten optical cycles then the transverse modulation instability is suppressed. Moreover, this suppression doesn't depend on the pulse amplitude or media properties. Let us consider the results of numerical simulations using Eq. (7) (i.e., keeping lowest nonzero term inD) to confirm the qualitative analysis. Since we study the problem of structural stability, we restrict numerical simulation by the case (2D+1) to simplify the analysis, because the result depends weakly on dimension of the problem. Figure 2 shows the results of numerical simulation for the wave packet with input distribution at two different cases demonstrating different evolution of the filamentation instability. The evolution of pulse intensity |Ψ(z, x, τ )| 2 in the absence of the fourth term in Eq. (7), i.e., by neglecting pulse steepening, is shown in Fig. 2(a). The pulse evolution with nonlinear dispersion taken into account is shown for comparison in Fig. 2(b). As seen in Fig. 2(a), filamentation instability splits the wave packet into separate pulses in the transverse direction [33,34]. Here, the contour level lines are shown by the green color. However, if nonlinear dispersion is taken into account in Eq. (7), the instability is stabilized. It is clear from Fig. 2(b) that the inhomogeneities pass to the rear part of the pulse and cease to grow. This is the key difference between the considered pulse evolution mode and the laser pulse evolution within the framework of an ordinary nonlinear Schrödinger equation. A similar process of the filament instability stabilizing for few-cycle pulses occurs also for the laser pulse dynamics in the framework of the original wave equation (6) for (3D+1) case. Figure 3 shows the results of evolution of two different Gaussian laser pulses comprising 30 (a) and 10 (b) optical cycles with initial noise level of about 10 −4 of pulse amplitude. The process of pulse selffocusing demonstrated in the (x, y) plane occurs in both cases together with the strong pulse compression shown in the (y, τ ) plane. However, the longer pulse compression is accompanied by simultaneous development of the filamentation instability (see the (x, y) cross section in Fig. 3(a)), i.e., the pulse splits into separate filaments. For the shorter pulse ( Fig. 3(b)) the spatial structure of the pulse remains smooth in the process of adiabatic shortening of the pulse duration (see the (x, y) cross section). So, the transverse instability can be suppressed, and no violation of the pulse symmetry is observed. Correspondingly, we perform a thorough numerical study of the pulse dynamics for the axisymmetric case.
Thus, as shown by the results of the analytical and numerical study, laser pulses with durations of less than ten optical cycles are not a subject to filamentation instability due to medium nonlinear dispersion.

VI. RESULTS OF NUMERICAL SIMULATION
It was shown in the previous section that the spatial modulation instability can be suppressed, and, correspondingly, the pulse symmetry will not be violated. To perform detailed numerical analysis of soliton selfcompression on the basis of Eq. (6), in what follows we will turn to studying the dynamics of axisymmetric circular polarized laser pulses.
As the initial distribution of the laser pulse, we will specify the soliton-like distribution (Eq. (12)) along τ and the transverse Gaussian distribution with the characteristic scale a: Here, γ is the soliton velocity, and G and φ are the amplitude and phase distribution of the soliton, which are found by solving system of equations (12), and the parameter N characterizes the number of solitons, into which the initial laser pulse will be decomposed at the asymptotic stage within the framework of the onedimensional problem (∆ ⊥ ≡ 0).
A. Single-soliton dynamics (N 1) Figure 4 presents the results of numerical simulation with the initial distribution (49) with N = 1, δ 0 = 0.03, ω s = 1, and a = 400. In this case, the duration of the wave packet corresponds to ten optical cycles (τ in p = 10 T 0 , where T 0 is the field period). It should be emphasized, that the distribution of the laser pulse with the longitudinal scale being much smaller than the transverse one is specified at the input to the nonlinear medium, i.e., the dispersion length of the wave packet is much shorter than the diffraction length in the dimensionless equation (6). One can see from Fig. 4(a-f ) that self-focusing of the radiation in the transverse direction is accompanied by the adiabatic decrease in the soliton duration. The evolution of the field in the pulse on the beam axis u x (r = 0) = Re(u) is shown as the blue line. It is seen in the figure that the pulse scales decrease significantly. For further adjustment of the numerical simulation results and the above-presented qualitative analysis of the problem, it is convenient to determine the integral pulse width ρ ⊥ (10). Analysis shows that the average pulse width ρ 2 ⊥ decreases by 3 times, from ρ 2 ⊥ = 400 to ρ 2 ⊥ 133, and the intensity of the field in the compressed pulse increases by 230 times.
The dashed line in Fig. 4(g) represents the initial distribution of the pulse envelope |u| = u 2 x + u 2 y at the beam axis and the field distribution in the compressed pulse for u x = Re(u) at the output of the nonlinear medium, z = 2700. One can see that at the half-intensity level, the laser pulse is compressed by 14 times, from τ in p = 10 T 0 to τ out p = 0.71 T 0 , which corresponds to the duration being slightly less than the optical cycle, while the r.m.s. duration of the pulse τ pulse , which is calculated on the basis of the second-order momentum, amounts to τ pulse = 1.1 T 0 . Here, τ is the center of mass of the wave packet.
It is seen that the pulse duration τ pulse averaged with respect to the transverse distribution of the field intensity, is 1.5 times larger than the duration of the wave packet pulse at the beam axis. Then, we use formula (26) to estimate the duration of the compressed laser pulse. The average wave packet width have decreased by 3 times. Therefore, the pulse duration determined by using Eq. (26) should decrease to a = 1.05 T 0 , which is a little different from the r.m.s. duration τ pulse of the laser pulse. Therefore, the results of numerical simulation agree with the above-presented qualitative analysis.
It should be emphasized that in the process of nonlinear dynamics of the laser pulse, the longitudinal scale is always smaller than the transverse one, i.e., the distribution of the wave packet is not symmetrized and, therefore, self-focusing of the pulse does not pass over to the regime of spherically symmetric collapse.
Evidently, such a great decrease in the duration of the laser pulse should be accompanied by a great widening of the wave packet spectrum. The dashed line in Fig. 4(h) shows the spectral intensity of the input pulse at the beam axis, and the solid line is the spectral intensity of the compressed pulse. One can see in this figure that the spectrum of the wave packet at the output from the nonlinear medium is asymmetric and alike the spectrum presented in Fig. 1(d) for the precise soliton solution is found within the framework of the one-dimensional problem (∆ ⊥ ≡ 0). This asymmetry of the spectrum intensity takes place due to the fact that, as the pulse duration decreases, the term being responsible for steepening of the wave packet profile starts manifesting itself. As it has been already mentioned, in the case of short durations, the soliton solutions have a sufficiently strong frequency modulation (see Eq.(12a)), which is reflected as significant widening of the short-wave part of the spectrum. Note, that the pulse compression in the beam axis is stronger than the average one, since the field intensity is larger in the near-axis region of the pulse.
The solid red line in Fig. 4(i) shows the dependence of the pulse duration, which is determined at the halfmaximum of the intensity normalized with respect to the initial value, on the evolution variable z. This plot indicates that the dependence of the wave packet duration has two scales. It follows, in particular, from the system of equations for the wave packet duration (26) within the framework of NSE. In the case, where the wave packet width is much smaller than the initial value (a ⊥ a ⊥0 ), which corresponds to a significant decrease in the wave packet duration, the behavior of the pulse duration is described by Eq. (30b). Let's consider now the case, where the pulse width is slightly smaller, a ⊥ − a ⊥0 = δa ⊥ 1, which corresponds to an insignificant decrease in the wave packet duration.
The dashed line in Fig. 4(i) shows the dependence of the part of the initial energy in the compressed pulse on the evolution variable z. It is seen from the figure presented that the compressed wave packet contains more than 55% of the initial energy and, since the pulse duration is ten times shorter, the peak power of the compressed pulse increases by 5 times. The remaining energy is located in the halo around the central peak.

B. Multi-soliton dynamics (N ≥ 2)
Let's consider now the case of higher N . Our interest in this case is connected with the fact that we can operate with high energies in the laser pulse. This is important in the problem of optimization of the laser pulse selfcompression, since the length of the nonlinear medium, at which the laser pulse duration reaches the minimal value, increases in proportion to the increase in the initial transverse width a (Eq. (32)). It will be shown in what follows that for N > 1, the length of the laser pulse compression will become several times shorter as compared with the initial soliton-like distribution at N = 1.
Within the one-dimensional problem (∆ ⊥ ≡ 0), i.e., in the absence of spatial effects, the nonlinear dynamics of the laser pulse is determined entirely by the parameter N . It has been already noted that the initial distribution of the pulse coincides exactly with the soliton solution, Eq. (11) in the case of N = 1. At N ≥ 2, the initial wave packet at the asymptotic stage will split into a sequence of solitons with the parameters δ n = (2n − 1)δ 0 , where n = 1, . . . , [N ] is a sequence of integer numbers [24]. Here [N ] is the integer part of N . An important feature of the wave solitons under consideration (Eqs. (11) and (12)) is the semi-bounded spectrum of their permissible solutions, i.e., the presence of the bounding parameter δ cr corresponding to the limiting soliton with the minimum possible pulse duration and, correspondingly, the maximum permissible amplitude. Therefore, the number of solitons, into which the initial pulse splits actually, is an integer number [N ], obeying effectively the following inequality (δ ([N ]) < δ cr ): The pulse dynamics will be more complicated, but in the long run, the solitons with δ ([N ]) < δ cr will be formed.
We turn now to the initial problem and analyze the dynamics of the laser pulse (Eq. (49)) within the framework of Eq. (6) allowing for the spatial effects. Figure 5(a-f ) presents the spatio-temporal evolution of the laser pulse at N = 2.05 and a = 400, when the distribution was specified in the longitudinal direction at δ = 0.03 and ω s = 1. It is seen from this figure that at the initial stage (z 350), the laser pulse is strongly compressed according qualitative law (26). At this the pulse width decreases not so strongly.
Then, just as within the framework of the onedimensional problem, the wave packet starts blurring in the longitudinal direction (z ∼ 660), which leads to a decrease in the pulse self-focusing rate, since the amplitude of the field has decreased. As seen from the figure, at z = 660, a horseshoe-shaped structure starts forming. Then, as follows from the figure, at z ∼ 850, the pulse starts splitting in the longitudinal direction into two solitons. The soliton with a short duration and, hence, a high amplitude, is located at the rear front of the pulse (τ ∈ [−25, 25]), and the soliton with a greater duration and, hence, a lower amplitude is located at the leading front of the time distribution (τ ∈ [−150, 0]). Note that in the near-axis part of the pulse the field becomes again stronger, than in the edge region r ∼ a/2. Further, Figure 5. (a-f ) Dynamics of the circularly polarized field |u(z, τ, r)| for N = 2.05, δ0 = 0.03, and ωs = 1 with the Gaussian distribution in the transverse direction with width a = 400. The blue line in inplots shows the pulse evolution at the beam axis (r = 0) for one of the field components, ux = Re(u). Here, the field is normalized with respect to the maximum value. The (g) dashed blue line is the distribution of the field envelope of the input laser pulse at the beam axis |u(τ, r = 0)|, the red line is the distribution of the compressed-pulse field at the beam axisux = Re(u), the dashand-dotted magenta line is the distribution of the time mask, the dashed magenta line is the distribution of the spectrum intensity after application the time mask over the compressed pulse, the (h) dashed blue line is the initial spectrum, the red line is the spectrum of the compressed pulse, and the (i) red line is the dependence of the wave packet duration normalized with respect to the initial duration for z.
due to the process of pulse self-focusing, the duration of these solitons will decrease adiabatically. Evidently, the short soliton will be compressed faster, since its energy is greater, I 2 12πδ 0 a 2 0 /a 2 cur , as compared with the longer one I 1 4πδ 0 a 2 0 /a 2 cur , where a cur is the current beam width, and a 0 is the initial beam width. As a result of the laser beam self-focusing, a soliton with a shorter duration will be segregated which is shown in Fig. 5(a-f ) for z ∼ 1010.
The blue dashed line in Fig. 5(g) shows the initial distribution of the envelope of the pulse |u| = u 2 x + u 2 y at the beam axis and the red line shows the distribution of the field in the compressed pulse for u x = Re(u) at the output of the nonlinear medium, z = 1010. It is seen from the figure that the laser pulse at the half-maximum of the field intensity is compressed by 19 times, from τ in p = 10 T 0 to τ out p = 0.52 T 0 , which corresponds to a value being slightly less than the optical cycle. As is seen in the figure, a long soliton, which is almost unnoticeable and is represented by a small pedestal at the level 0.02, is located at the leading edge of the time distribution of the laser pulse τ ∈ [−125, −25]. Therefore, this pedestal will be totally unobservable in the intensity distribution against the background of the main signal.
Such a significant shortening of the output pulse duration, as it has been already mentioned, should be accompanied with a significant widening of the spectrum of the compressed pulse. The blue dashed line in Fig. 5(h) represents spectral intensity of the input pulse at the beam axis, and the solid red line shows the spectral intensity of the compressed pulse. Note that the ruggedness of spectral intensity of the laser pulse at the output of the nonlinear medium z 1010 appears due to the interference of two temporally separated wave structures.
To develop spectral intensity of just one soliton with a short duration, we applied the time-interval mask M(τ ) to the compressed pulse, in order to remove the soliton with a longer duration, which is located on the leading edge. The distribution of the time-interval value M(τ ), Eq. (53) is shown in Fig. 5(g) as a dash-and-dotted magenta line, whereas the dash-and-dotted magenta line in Fig. 5(h) rerpresents spectral intensity of the output laser pulse at the beam axis after application of the timeinterval value M(τ ). As seen in Fig. 5(i), the resulting spectrum has become smooth and more asymmetric, and looks like the spectra shown in Figs. 1(d) and 4(h). In this case, the spectral intensity of the short soliton is wider than the spectral intensity of the compressed laser pulse for the case, where only one soliton contains in the initial wave packet (N = 1). The solid red line in Fig. 5(i) shows the dependence of the pulse duration, which is determined at half-maximum of the field intensity normalized with respect to the initial value τ 0 , on the evolution variable z. This figure demonstrates two stages in the evolution of a laser pulse, which we discussed earlier. Specifically, as follows from Fig. 5(i), the duration of the wave packet reaches the intermediate minimum τ p 0.18τ 0 at z 380. Presumably, one can restrict consideration to this medium length for compression of the initial laser pulse. However, in this case, self-compression of the wave packet is rather sensitive to the length of the nonlinear medium, since, as seen in Fig. 5(i), at z ∼ 400 the pulse duration starts increasing again, just as within the framework of the one-dimensional problem, and reaches τ p = 0.52 τ 0 at z 650. It has been already mentioned that the pulse splits into two solitons, which start compressing adiabatically. So, from the practical point of view, it is preferable to compress a pulse at z 800 due to following reasons. First, the amplitude of the tailing soliton exceeds significantly the amplitude of the leading soliton. Second, the decrease in the duration of the wave packet occurs monotonically, as the pulse propagates in the medium.
Note that at z 1010, the duration of the laser pulse is 3.9 times shorter, than at z 380. As it follows from the comparison of Figs. 4(i) and 5(i), the initial laser pulse at N = 2.05 compresses along the path z 1010, which is significantly shorter, than the path z 2700 for the case N = 1. It should be noted that such a significant decrease in the length of the wave packet compression happens due to the pulse splitting into two structures, the duration of a shorter pulse decreases by three times compared with the initial duration, which in the long run will lead to a decrease in the compression length.
In the process of the further increase of the parameter N , self-compression of the laser pulse in the process of self-focusing of the transverse distribution is preserved (Fig. 6, 7). The figure presents the evolution of the laser pulse for two different values of the parameter N : (a) -N = 2.5, (b) -N = 3.02. It follows from the figure that at z ∼ 650 (a), the laser pulse splits into two wave structures, and the soliton with a lower amplitude becomes more pronounced in the background of a shorter-duration soliton, unlike the case shown in Fig. 5(a-f ). At N = 3.02, the laser pulse splits into three structures at z ∼ 450 (see Fig. 7). It is evident that for self-compression of the input laser pulse, it is desirable to restrict consideration to the length of the nonlinear medium being z ∼ 295 for N = 2.5 and z ∼ 200 for N = 3.02, since further the time structure of the laser pulse becomes rather complicated.
Thus, in the case of self-compression of laser pulse under the conditions of self-focusing of spatial distribution, when the linear dispersion length is much shorter than the diffraction length, the key role in the pulse dynamics is played by the solitons found within the framework of the one-dimensional problem [24]. In the case of N > 1, the prevalent role is played by the modulation instability of the wave field, rather than the filamentation one. It leads to splitting of the laser pulse into a set of soliton-like structures [24], which further, due to the process of selffocusing, will compress individually and monotonically in the longitudinal direction. As shown by the numerical analysis, in order to obtain extremely short laser pulses with high time contrast, it is preferable to specify wave field distribution (49) with N 2.5 to the input of the nonlinear medium.

VII. CONCLUSIONS
In this paper, we justify theoretically the promising method of self-compression of multi-millijoule laser pulses up to one optical cycle. Self-focusing of a wave field Figure 6. Dynamics of the circularly polarized field |u(z, τ, r)| for N = 2.5, δ0 = 0.03, and ωs = 1 with the Gaussian distribution in the transverse direction with width a = 400. The blue line in inplots shows the pulse evolution at the beam axis (r = 0) for one of the field components, ux = Re(u). Here, the field is normalized with respect to the maximum value. Coordinates r, z, τ are dimensionless according to Eq. (6). in a medium with the Kerr-type inertia-free nonlinearity and anomalous dispersion of group velocity leads to an adiabatic decrease of the wave packet duration to a duration comparable with the optical cycle for wide wave packets with the soliton-like field distribution along the longitudinal coordinate. Analysis shows that the soliton duration decreases in proportion to the square of characteristic width of the wave beam (26). It should be noted that the strongly oblate ellipsoidal distribution of the wave beam is preserved in the process of evolution, and no symmetrization occurs. Self-compression of the laser pulse proceeds under the conditions of a noticeable excess over the threshold value of the self-focusing power.
Thorough numerical studies of the evolution of a 3D axisymmetric wave packet have been performed. The results obtained on the basis of qualitative analysis in the aberration-free approximation have been confirmed by the numerical simulation. In the case of the quasi-soliton field distribution in the longitudinal direction, the pulse self-focusing is accompanied by a monotonic decrease of its duration down to the one optical cycle, corresponding to the duration of the limiting soliton (11).
In the case of the initial high amplitude of the laser pulse (N 2), the initial longitudinal distribution of the wave packet splits into a sequence of solitons, which further self-compress monotonically and diverge in the longitudinal direction. Numerical simulations show that it is preferable to use wave packets containing not more than two quasi-soliton structures.
We should denote that self-action of intense laser pulses having durations shorter than ten optical cycles is stable relative to the transverse filamentation instability [17]. This takes place due to allowance for the nonlinear medium dispersion (dependence of the group velocity on the amplitude), under which the type of self-focusing instability changes from the absolute to the convective one. As a result, in the case of sufficiently short pulses, the noise amplitudes shift relative to the pulse and have no time to reach to arbitrary noticeable values.
In recent paper [35], the possibility of laser pulse selfcompression at a wavelength of 3.9 µm from 94 fs to 30 fs was demonstrated for an input radiation power exceeding the critical self-focusing power by four orders of magnitude. The duration of the wave packet decreased by three times after it had passed through a YAG plate having a thickness of 2 mm. Simple evaluations show that the plate thickness was chosen to be less than the length, at which the filamentation instability develops, which leads to decomposition of the transverse field distribution to separate filaments. However, our study shows that much stronger compression down to 7 fs pulse duration (about optical cycle) can be achieved under the similar conditions, but with 7. . . 8 times larger medium length.
Let us place some estimates for the realistic experimental realization of proposed laser pulse compression. The typical initial parameters of an appropriate laser pulse are: the wavelength of 4µm, the initial duration τ in = 100 fs, the energy W in = 15 mJ, the initial