MODELING ULTRASHORT FILAMENTS OF LIGHT

Laser sources nowadays deliver optical pulses reaching 
few cycles in duration and 
peak powers exceeding several terawatt (TW). When such 
pulses propagate in transparent 
media, they first self-focus in space, until 
they generate a tenuous plasma by photo-ionization. These pulses 
evolve as self-guided objects, resulting from successive equilibria between the Kerr focusing process, the defocusing action of 
the electron plasma and the chromatic dispersion of the medium. Discovered ten years ago, this self-channeling mechanism reveals 
a new physics, having direct applications in the long-distance propagation of TW beams in air, supercontinuum emission as well as pulse 
self-compression. This review presents the major progress in this field. Particular emphasis is laid to the derivation of the propagation 
equations, for single as well as coupled wave components. Physics is discussed from 
numerical simulations and explained by analytical arguments. 
Attention is also paid to the multifilamentation instability, which breaks up broad beams into 
small-scale cells. Several experimental data validate theoretical descriptions.

1. Introduction. Progress in laser technology have permitted the generation of light wave packets with smaller and smaller durations, down to the femtosecond scale (1 fs = 10 −15 sec.), and optical intensities of several hundreds of terawatt (TW) per cm 2 . These ultrashort pulses enable researchers to probe ultrafast processes on never-before-accessed time scales and study light-matter interactions at high intensity levels, from which the pulse strength overcomes the Coulomb barrier and triggers optical-field ionization.
Ten years ago, infrared pulses with femtosecond durations and gigawatt (GW) powers were surprisingly observed to become confined to a long-living, self-confined tube of light capable of covering several tens of meters along the propagation axis [27]. The mechanism supporting this "light bullet" results from the balance between Kerr focusing, which increases the local optical index with the laser intensity, and self-induced ionization. Through this balance, the beam spatial profile exhibits a narrow extent along the optical path, which gives rise to a nonlinear object, called "femtosecond filament". In the diffraction plane, this filament is characterized by a white-light spot, surrounded by concentric rainbows of colors. The high Figure 1. Scheme for producing a femtosecond filament in air at the laser wavelength of 800 nm. Photograph at the right-hand side shows a transverse cut of the filament profile.
nonlinearities competing through the filamentation process produce an impressive supercontinuum leading to white-light emission, which is illustrated in Fig. 1.
Femtosecond filaments have been the topic of intense research activities, first because of their capability to convey high intensities over long distances, second because of their large spectrum that can be exploited in remote sensing experiments [37,74]. For appropriate beam configurations, filaments can be created not only in the atmosphere, but also in noble gases, liquids and dielectrics, as long as the electron plasma density remains at subdense levels. These optical structures are subject to strong modifications of their temporal profile caused, e. g., by the ionization front. Their duration can then be drastically shortened, potentially down to the optical cycle limit. Besides, femtosecond pulses with broad spatial extents create several filaments, whose mutual interactions support the self-guiding of the beam upon several kilometers. This recently led to develop ultrashort LIght Detection And Ranging (LIDAR) facilities, that exploit the white light emitted by filaments, in order to identify several pollutants within a single laser shot.
The present paper reviews current models describing the dynamics of femtosecond filaments. An accurate understanding of the filamentation requires a rigorous derivation of the propagation equations. For this purpose, Section 2 addresses the model governing the long range propagation of ultrashort laser pulses in any optically-transparent medium for isolated and two coupled pulses as well. Section 3 reviews the basic phenomenon of wave self-focusing and its limitation by potential players such as plasma generation or chromatic dispersion. Section 4 lists the major processes driving femtosecond filaments while Section 5 presents experimental results related to this field.
2. Propagation equations. The derivation of the equations modeling the propagation of ultrashort optical pulses follows the conventional description of nonlinear optics. Straightforward combination of the Maxwell's equations yields [1,66,136] where ǫ 0 , µ 0 and c denote the electric permittivity, magnetic permeability and the speed of light in vacuum, respectively. The optical electric field E, the polarization vector P , the carrier density ρ and the current density J are real valued. The current density J describes the motions of free electrons created by ionization of the ambient atoms. The polarization vector P originates from the bounded electron response driven by the laser radiation. It contains a linear part P L ≡ P (1) related to the firstorder susceptibility tensor ↔ χ (1) and a nonlinear one P NL satisfying | P (1) | ≫ | P NL |. Because we assume isotropic, homogeneous media and spectral ranges far from any material resonance, all susceptibility tensors ↔ χ (j) with even index j vanish due to inversion symmetry. P thus decomposes into an odd series of E such as P = P (1) ( r, ω) + P (3) ( r, ω) + P (5) ( r, ω) + . . .

Forward wave equations. The Fourier transform of Eqs. (1) expresses as
where ∇ 2 ⊥ ≡ ∂ 2 x + ∂ 2 y stands for transverse diffraction and F NL ≡ P NL + i J/ω (6) gathers all nonlinear contributions. Because Eqs. (5) are usually difficult to integrate in the full space-time domain, assumptions are requested to simplify them into more tractable form. The most fundamental of those consists in supposing that the wavefield keeps a transverse extension always fulfilling i. e., for k(ω) located around ω 0 , the transverse wave number of the beam (the inverse of the beam waist) is much smaller than the longitudinal one. The second assumption is that we deal with small nonlinearities, i. e., Conditions (7) and (8) have important impacts. First, they make vectorial effects negligible. Indeed, Eqs. (5) can be combined, with the help of the continuity equation ∂ t ρ + ∇ · J = 0, into the form whose last term scrambles nonlinear vectorial components. When we project all vectors onto the transverse and longitudinal axis, the nonlinear coupling of vectorial components described by the scrambling term behaves as O(k 2 ⊥ /k 2 ) [51,52,107]. So, vectorial couplings become important in the limit k ⊥ → k(ω) only. Otherwise, the influence of ∇( ∇ · F NL )/k 2 is negligible. This property justifies the use of a scalar description for linearly-polarized beams having, e. g., E x = 0, E y = 0. Second, the limit (7) implies that the wave components forming the angle θ = arcsin (k ⊥ /k) between the transverse and longitudinal directions mostly propagate forwards since θ ≪ π/2 [49]. In other words, the amount of backscattered photons is weak, compared with those propagating along the forward direction. The reason can be seen from the scalar version ( E = E e x ) of Eq. (9) expressed as where D ± (ω) ≡ ∂ z ∓ik(ω). By substituting the solution E = U + e ik(ω)z + U − e −ik(ω)z into Eq. (10), where U + and U − represent the Fourier components of the forward and backward running fields, and Taylorizing U ± over one optical cycle, Fibich et al. demonstrated that the backscattered component has a weak influence on the beam dynamics if Eq. (7) is verified [54]. Third, current models require that the second-order derivative in z of the envelope function U + for the forward component must be small compared with |k(ω)∂ z U + | [60,70]. This approximation, expressed as |∂ z U + | ≪ |k(ω) U + |, refers to the "paraxiality" assumption. It holds if the field envelope U + does not significantly change over propagation distances of the order of λ. Paraxiality is again linked to the weakness of both the ratio k ⊥ /k and the nonlinearities: If we assume the nonlinear function F NL clamped at a maximal constant level, F NL = F max E, the forward component then goes like This solution fits that of the paraxial model discarding second derivatives in z, i. e., U + ∼ e −i(k 2 ⊥ −µ0ω 2 Fmax)z/2k(ω) , as long as the two constraints (7) and (8) apply. Under the previous approximations, E is mainly described by U + e ik(ω)z and the so-called Forward Maxwell Equation (FME) emerges from Eq. (10). Alternatively, we may decompose the same equation (10) upon the linear modes, E lin = e ±i √ k 2 (ω)+∇ 2 ⊥ z , and project the nonlinearities along the forward direction only [82,83]. Rewriting Eq. (10) as with D ⊥ ± (ω) ≡ (∂ z ∓ i k 2 (ω) + ∇ 2 ⊥ ) and retaining the forward running component ⊥ z immediately yields under the unidirectional limit Models (12) and (14) are analogous when the condition k 2 ⊥ /k 2 ≪ 1 applies. Their major advantage is to elude the formal use of a central optical frequency and correctly describe the complete spectrum of nonlinear pulses, even when they develop very large bandwidths.
Let us now introduce the complex version of the electric field where c 1 ≡ ω 0 µ 0 /2k 0 defines an optical intensity I ≡ |E| 2 expressed in W/cm 2 and Θ(x) denotes the Heaviside function. Because E satisfies E * (ω) = E(−ω) * ( * means complex conjugate), it is then sufficient to treat the FME model (12) in the frequency domain ω > 0 only. When a central frequency ω 0 is imposed, Eq. (12) restitutes the Nonlinear Envelope Equation (NEE), earlier derived by Brabec and Krausz [26]. By virtue of the Taylor expansion whereω = ω − ω 0 is the envelope frequency shift, k ′ = ∂k/∂ω| ω=ω0 is the inverse of the group velocity and k (n) ≡ ∂ n k/∂ω n | ω=ω0 , Eq. (12) can be developed as where terms with k(ω) in their denominator are expanded up to first order inω only. Since E(ω) is the Fourier transform of E(t)e iω0t inω,ω corresponds to i∂ t acting on the field envelope by inverse Fourier transform. This naturally leads us to introduce the complex-field representation involving the novel envelope function U . Next, the new time variable t → t − z/v g shifts the pulse into the frame moving with the group velocity v g = k ′−1 . Eq. (17) then restores the NEE model where D ≡ +∞ n≥2 (k (n) /n!)(i∂ t ) n contains group-velocity dispersion (GVD) with coefficient k (2) = k ′′ and the operator T = 1 + i ω0 ∂ t accounts for deviations from the slowly-varying envelope approximation. This description holds whenever the difference between group and phase velocity is small, i. e., |k 0 − ω 0 k ′ |/k 0 ≪ 1. The operator T −1 introduces space-time focusing in front of the diffraction term The nonlinearities with the envelope function F env NL (U ) are also affected by the operator T , which refers to self-steepening.

2.2.
Optical nonlinearities. We henceforth assume a linearly polarized field supporting a scalar description. For simplicity, the Kerr response induced by the nonlinear polarization vector is limited to the cubic order in Eq. (2), since higher-order saturation concerns small corrections with, e. g., χ (5) /χ (3) ∼ 10 −12 in gases [25,136]. For centro-symmetric materials, only one relevant component remains in the cubic polarization P (3) [1]. Considering the susceptibility tensor χ (3) as keeping a constant value around ω 0 , the cubic polarization simplifies with a single component, such as This expression is valid whenever we suppose an instantaneous response of the medium. Strictly speaking, however, the phenomenon of Raman scattering comes into play when the laser field interacts with anisotropic molecules. Molecular scatterers have rotational eigenstates, which are excited through three atomic levels.
2.3. Plasma nonlinearities. Plasma generation implies the creation of free electrons by ionization processes. Free electrons induce a current density J = q e ρ v e , which depends on the electron charge q e = −1.6 × 10 −19 C, the electron density ρ and the electron velocity v e . J is computed from the fluid equations [47,143] where S represents external plasma sources and ν e is the electron collision frequency. Equations (24) can be combined to yield up to ponderomotive forces acting on slowly-varying time scales. These forces are negligible, as long as the peak intensities remain below 10 15 W/cm 2 , which will always be satisfied in the present scope [119,145]. Assuming electrons born at rest, the electron density is only governed by the source term S, i. e., that involves photo-ionization processes with rate W (I), collisional ionization with cross-section σ [99,100], and a function describing electron recombination with neighboring ions denoted by f (ρ). Here, ρ nt and U i are the density of neutral species and the ionization potential, respectively. Below, we shortly discuss the meaning of these three contributions.
(i) W (I) is the rate for photo-ionization. It is evaluated from perturbative theories valid as long as the electric field E is weaker than the atom field strength. This rate usually follows from theories applying to atoms or molecules and derived either by Keldysh, by Perelomov, Popov and Terent'ev (PPT), or later by Ammosov, Delone and Krainov (ADK) [9,77,121,122,123]. These stress two major limits bounded by the so-called Keldysh parameter, namely, the limit for multiphoton ionization (MPI, γ ≫ 1) and the tunnel limit (γ ≪ 1) concerned with high intensities, from which electrons can tunnel out the Coulomb barrier. Here, E p denotes the peak optical amplitude. For laser intensities I = |E| 2 < 10 13 W/cm 2 , the MPI limit applies and the ionization rate has the simple dependency where K = mod(U i / ω 0 ) + 1 is the number of photons necessary to liberate one electron. Despite the complexity of the available ionization formulas, all of them exhibit common dependencies on the laser field. Figure 2 illustrates the PPT ionization rate applied to a gas of xenon atoms and the ionization rate for crystals derived by Keldysh and applied to SiO 2 (silica) glass. Note that as long as the beam saturates below 10 13 W/cm 2 , the MPI limit alone can be retained. To obtain qualitative behaviors, this approximation can formally be used at higher intensities, while keeping the interaction physics valid.
Because free carriers are generated by photo-ionization, corresponding losses must be included into the propagation equations. These are determined by a local version of the Poynting theorem, i. e., d dt where ∂ t ρ PI ≡ W (I)(ρ nt − ρ) yields the amount of energy lost per time and volume units by single ionization events [71,126]. (ν e + iω) (ρ E).  Figure 2. Ionization rates. The solid curve corresponds to the PPT theory applied to a Xenon gas with U i = 12.1 eV at the wavelength λ 0 = 800 nm. The dashed curve follows from Keldysh's theory for crystals applied to a silica sample with potential gap U i = 7.8 eV at the same laser wavelength.
The current density term in Eq. (1a) then transforms as after introducing the critical plasma density ρ c ≡ ω 2 0 m e ǫ 0 /q 2 e , at which the laser wave number vanishes. The cross-section then provides a frequency-dependent collisional rate, which characterizes avalanche (cascade or impact) ionization in Eq. (26) [48,79,165]. For comparison, electron collision times ∼ 1/ν e are about 20 fs in silica and ten times larger in gases at atmospheric pressure.
(iii) Electron recombination leads to a decrease of the plasma density when freed electrons re-attach with ions. This varies with the medium density. In Eq. (26), the function f (ρ) in gases has a quadratic dependency on ρ [109,144,154] and recombination times belong to the nanosecond scale. In dielectrics, much shorter recombination times are involved (τ recomb = 50 − 150 fs) and the density linearly decreases like f (ρ) = ρ/τ recomb [13,118,155].
By combining Eqs. (12) and (23a), the propagation equation within the FME description reads in the Fourier domain as where For practical use, the collision cross-section σ is stated for a central frequency ω 0 . We furthermore assume ν 2 e /ω 2 0 ≪ 1 and ǫ(ω 0 )/ǫ(ω) ≈ 1. The link to the NEE model (19) is then straightforward, whenever the dispersion relation supports a Taylor expansion around a central frequency, ω 0 . By retaining only waveforms beating at this frequency, the nonlinear envelope equation for the forward component U is directly inferred from (19) as Both models (33) and (35) are closed with the equation (26) governing the evolution of the plasma density ρ. They both describe wave diffraction, Kerr self-focusing, plasma generation and related losses, chromatic dispersion, space-time focusing and self-steepening. Their last term includes multi-photon absorption (MPA). They are usually integrated numerically by considering, e. g., a super-Gaussian beam with the envelope function [Eq. (18)] as initial condition, where r = x 2 + y 2 . Pulses can be focused through a lens of focal length f and be temporally chirped if C = 0. For Gaussian beams (N = 1), U 0 = 2P in /πw 2 0 involves the input power P in , w 0 is the beam waist and t p the 1/e 2 pulse half-width, such that its full-width-at-half-maximum (FWHM) is ∆t = √ 2 ln 2t p . The initial level (t = −∞) of plasma density is zero. The pulse dynamics depends on the material parameters. Some of them, used throughout the present review, are given in Table 1 for different laser wavelengths.
In linear, parallel (f = +∞) propagation, such Gaussian pulses diffract over the well-known Rayleigh distance In nonlinear propagation regime, the interplay between all competitors in Eq. (35) prevents the immediate diffraction of the beam central core and can maintain the most intense components of the pulse in self-guided state over longer distances.
The major differences between Eqs. (33) and (35) are illustrated in Fig. 3. We performed three different series of simulations at 800 nm. The electron plasma is created from O 2 molecules (see Table 1), whose ionization potential is lower than that of nitrogen. The first simulation concerns a direct integration of Eq. (33). The second refers to the same pulse described by Eq. (35), which eludes the production of high-order harmonics. The third approach relies on Eq. (35), in which temporal steepening is omitted, i. e., we set T = T −1 = 1. Insets show the on-axis spectra at maximal extent. Spectra are here computed from the Fourier transform in time of the pulse intensity at the center r = 0. Following the FME description, conversion of the pump into third harmonic (TH) is limited, whereas the fundamental is widely extending towards the blue/UV wavelengths [ Fig. 3(a)]. In the NEE description, there is no harmonic generation. However, spectral broadening (or supercontinuum) is so amplified in the blue region by temporal steepening (T, T −1 operators), that it overlaps the third-harmonic zone and simply hides it A first observation can be drawn from Fig. 3: The maximal intensity reached in the FME model is lower and the frequency variations are diminished compared with the NEE spectra for the pump wave alone. These differences are caused by TH generation, which limits the increase of the pump wave. Apart from this difference, no significant other change is visible between both models, so that NEE seems equivalent to the FME description without the self-generated harmonics. In contrast, omitting the temporal derivatives of the operators T, T −1 implies more serious discrepancies, as can be seen from Fig. 3(c).

Coupled pulses.
The previous strategy for deriving the NEE model can be applied to the coherent coupling of two pulse components characterized by different frequencies and wavenumbers, (ω 0 , k 0 ) and (ω q , k q ). For simplicity we here consider configurations where ω q is an integer multiple of the pump frequency ω 0 (ω q = qω 0 ). Both frequencies are far from each other. We start from the scalar equation of Eq. (5a) with ∇ · E = 0 and decompose the total real electric field E into two components This yields the set of coupled equations where | kj ,ωj means that we select the nonlinearity with oscillations in e ikj z−iωj t . We shall use the Taylor expansion k(ω) ≃ k j +k ′ jω + n≥2 (k (n) j /n!)ω n withω = ω −ω j , k ′ j = ∂k/∂ω| ωj being the inverse of the jth group velocity and k The inverse Fourier transform of Eq. (39a) can then be re-expressed in the frame moving with the pump wave, by applying the change of variables z → ξ = z, 0 are omitted. The standard approximation |1 − ω 0 k ′ 0 /k 0 | ≪ 1 is also applied, such that Concerning the harmonic field, Eq. (39b) can be similarly manipulated. Reframed in the group-velocity frame of the pump field, it reduces to where and we suppose paraxiality on the harmonic field envelope, We next perform the following steps Contributions in ∂ τ ∂ ξ U q coming from the last term of Eq. (44a) can be discarded, compared with the third term of Eq. (42), by virtue of (43) [112]. Also, we assume that phase and group velocities of the harmonic field are close to each other ( The propagation equation for the harmonics then reads The Kerr nonlinearities are computed from the cubic response (20) of the total field (38), for which dispersion in the susceptibility tensor is ignored and only an instantaneous contribution is considered. In practical uses, the resulting formulation applies whenever the pulse U q remains of small amplitude compared with the pump wave, i. e., |U q | ≪ |U 0 |. In that case, plasma nonlinearities can be derived within an MPI-like formulation, such as where W j (I j ) = σ Kj I Kj j and I j ≡ |U j | 2 , discarding electron recombination. Applied to third harmonic (TH) generation, i. e., ω q = 3ω 0 , the polarization vector P has two distinct contributions, one beating at ω 0 , the other one at 3ω 0 . In the propagation equations, after normalizing the intensity with c 1 = ω 0 µ 0 /2k 0 , we can perform the change U 3ω0 →Ũ 3ω0 e i∆kz in the equation for the harmonic field, where ∆k ≡ 3k 0 − k 3ω0 yields the phase mismatch in wavenumber. Omitting tilde symbols for notational convenience and introducing the group-velocity mismatch parameter ∆v = (k ′ 3ω0 − k ′ 0 ) −1 , equations for TH generation finally read In the above model, ionization losses apply to each component, separately. Intensitydependent couplings (I j E k with j = k) refer to cross-phase modulation (XPM). Phase-dependent nonlinearities correspond to four-wave mixing (FWM). To make the reading easier, the meaning of the principal abbreviations used in this paper is recalled in appendix.
3. Self-focusing and arrest of collapse.

3.1.
Wave blow-up. For beams with no temporal dispersion, Eq. (35) with only an instantaneous cubic Kerr term (x K = 0) describes the self-focusing of optical wave packets with n 2 > 0 [11,105]. Self-focusing makes the field intensity increase and causes a compression of the beam in the diffraction plane. This dynamics leads to "wave collapse" when the nonlinearity is not saturated. A necessary condition for the collapse is that the input power P in = |E| 2 d r exceeds some critical value, computed on the Townes mode (see below). The beam waist decreases more and more, as the field amplitude |U | diverges ( Figure 4). For purely Gaussian beams, the blow-up distance, also termed as "nonlinear focus", is given by the well-known semi-empirical Marburger formula [105] Below we briefly review properties of wave collapse and how self-focusing is affected by temporal dispersion and plasma generation. In order to argue conveniently on these mechanisms, we rescale Eqs. (35) and (26) in dimensionless form by using the substitutions r → w 0 r, t → t p t, z → 4z 0 z, U → √ c 2 ψ and ρ → (n 2 0 ρ c /2z 0 k 0 )ρ, where c 2 ≡ λ 2 0 /8π 2 n 0 n 2 w 2 0 . They result into the extended nonlinear Schrödinger equation (NLSE)  where ∂ t ρ = Γ|ψ| 2K when we neglect collisional ionization and recombination. F (ψ) includes "perturbations" of a critical collapse, such as GVD, MPI, MPA and pulse steepening. GVD and MPA have the normalized coefficients δ With a purely cubic nonlinearity, solutions to the Cauchy problem (51a) can blow up at finite distance [78,148] .., formally accounts for the dispersion of a wave-packet along D orthogonal axes [ r = (x, y, t, ...)], where D = 2 or 3. The wavefunction ψ evolves from the spatially-localized initial datum ψ( r, 0) ≡ ψ 0 ( r), assumed to belong to the Hilbert space H 1 with finite norm ψ H 1 = ( ψ 2 2 + ∇ψ 2 2 ) 1/2 , where f p ≡ ( |f | p d r) 1/p . Two invariants are associated with ψ, namely, the L 2 norm (power) P and Hamiltonian H: The following "virial" equality can be established [62,157] P d 2 where r 2 = r 2 |ψ| 2 d r/P denotes the mean-squared radius of the solution ψ. By a double integration in z, Eq. (53) shows that, whenever D ≥ 2, there exist initial conditions for which r 2 vanishes at finite distance. For finite norms P , the inequality P ≤ (2/D) 2 r 2 × ∇ψ 2 2 thus implies that the gradient norm diverges in collapse regimes. As H is finite, the collapse dynamics makes the L 4 norm ψ 4 4 blow up in turn and max r |ψ| diverges accordingly, by virtue of the mean-value theorem |ψ| 4 d r ≤ max r |ψ| 2 × P [87]. This leads to a finite-distance blow-up, at which the solution ψ stops to exist in H 1 [129]. This mathematical singularity reflects the ultimate issue of the nonlinear self-focusing in the absence of saturation of the Kerr response and under the paraxiality assumption.
While H < 0 arises from Eq. (53) as a sufficient condition for collapse, sharper requirements can be derived by means of the Cagliardo-Nirenberg inequality In the critical case D = 2, this inequality is used to bound H from below, so that the gradient norm blows up only if P fulfills the constraint P > P c . The best constant in Eq. (54) is exactly C best = 2/P c . It involves the quantity P c = R 2 d r = 11.68, where R is the radially-symmetric soliton solution, called the "Townes mode", of [36,161]. P c justifies the existence of a critical power for the 2D self-focusing of optical beams in nonlinear Kerr media. In the supercritical case D = 3, a criterion for collapse, sharper than H < 0, follows from a combination of Eqs. (53) and (54) as H < P 2 c /P for gradient norms initially above 3P 2 c /P [88]. Here, P c again corresponds to the mass of the 3D soliton satisfying −R + r −2 ∂ r r 2 ∂ r R + R 3 = 0.
In collapse regime, the solution focuses self-similarly near the singularity point z c as [135] ψ( r, z) where ξ = r/L(z), ζ(z) ≡ z 0 du/L 2 (u) and the parameter λ is positive for making the new wavefunction φ localized. The function L(z) represents the scale length that vanishes as collapse develops, and φ converges to an exactly self-similar form φ( ξ) fulfilling ∂ ζ φ → 0. For radial solutions Eq. (51a) transforms into where ] is viewed as a complex turning point, with β ≡ − 1 4 L 3 L zz . As L(z) → 0, φ can be treated by means of quasi-self-similar techniques [15]. The solution φ is split into a nonlinear core, φ c , extending in the range ξ ≪ ξ T , and a linear tail, φ T ∼ e −λπ/β /ξ 1+iλ/β , defined in the complementary spatial domain ξ ≫ ξ T where the nonlinearity vanishes. The length L(z) is then identified from the continuity equation describing the mass exchanges between the core and tail parts of ψ. The dynamics of self-similar collapses vary with the space dimension number as follows.
• For D = 2, the compression scale L(z) has a twice-logarithmic correction: [58,92,104]. As z → z c , the exponential contribution of the tail decreases to zero, while the core converges to the Townes mode R like The power P thus relaxes to the critical value P c and stays mostly located around the center. • For D = 3, β attains a fixed point β 0 = 0, leading to the scaling law L(z) The power is no longer preserved self-similarly in space, since P = L(z) |φ| 2 d ξ. This integral behaves as P ≃ P core (z) + P tail (z), where P core (z) ∼ L(z) vanishes, while P tail (z) contains almost all the initial mass as L(z) → 0. A 3D collapse thus promotes an expulsion of mass towards the large distances where it keeps a stationary density r 2 |ψ| 2 → const [85,95,158,168].
3.2. Saturation by temporal dispersion. Let us now consider Eqs. (51) in which F only contains group-velocity dispersion (F = −δ∂ 2 t ψ). A critical collapse in the transverse plane is then accompanied by a defocusing in time for k ′′ > 0 (δ > 0, normal dispersion) or a temporal compression for k ′′ < 0 (δ < 0, anomalous dispersion). The interplay of these processes results in the symmetric splitting of In all cases, plasma response is omitted.
the pulse along the time axis with normal GVD [34,35] and to a 3D spatiotemporal collapse with anomalous GVD [88].
With δ > 0, normal GVD transfers power towards non-zero instants, symmetrically located with respect to t = 0 [55,103]. It "splits" a focusing pulse into two regular, symmetric spikes at powers < 2P c . For higher powers, the peak edges develop shock profiles and disintegrate into ripplelike cells [61,56]. Because of the hyperbolicity of the operator ∇ 2 ⊥ − δ∂ 2 t , one splitting event transforms the optical field distribution in the (r, t) plane into an X-shaped waveform [18,43]. In addition, the operator T in front of the Kerr term (self-steepening) induces a shock dynamics: The field develops a singular profile with |ψ t | → +∞ in the trail of the pulse [10]. This dynamics is reinforced by space-time focusing [133]. Self-steepening and spacetime focusing produce a transfer of power from the leading pulse (corresponding to the region t < 0 in the time domain) to the trailing portion of the pulse (located at t > 0). This asymmetrizes the temporal distribution, which was retrieved by direct experiments [128]. Figure 5 depicts temporal profiles of pulses undergoing steepening effects in normally dispersive regime (silica glass at 790 nm), together with a transverse collapse. Note that when the pulse is only subject to self-focusing, it naturally undergoes on-axis compression in time (dotted curve) [105].
GVD is able to halt the wave collapse at powers moderately above critical. The stronger the GVD coefficient, the larger the power interval in which the collapse is arrested by dispersion. By solving the cubic NLSE, a boundary δ crit (p), function of the ratio of input power over critical,p = P in /P cr , can be calculated in such a way that initial conditions fulfilling δ > δ crit (p) may inhibit the divergence of the solution, through, e. g., GVD splitting for normal dispersion [102]. For anomalous GVD (k ′′ < 0), a mapping |δ| > δ crit (p) can again be constructed on the basis of virial-type arguments [21]. It predicts pulse spreading when the dispersion length ∼ t 2 p /|k ′′ | is short enough to prevail over diffraction and Kerr nonlinearity. Higherorder dispersion together with steepening effects (T, T −1 ) modify the theoretical curves in the sense of shifting the self-focusing threshold to higher powers [53,137,140]. Figure 6 summarizes the zone of collapse and no-collapse in the plane (δ,p) for both normally and anomalously dispersive media computed from Eq. (35). Open circles represent initial conditions that do not collapse; full circles indicate a strong divergence. Three zones clearly occur: (I) A dispersion-dominated domain leading to pulse spreading; (II) A transition zone in which GVD and steepening operators inhibit the self-focusing; (III) A Kerr-dominated region, in which chromatic dispersion is unable to stop the wave blow-up. Here, the field intensity can increase by several decades before reaching the ionization threshold.
3.3. Saturation by plasma defocusing. Plasma generation depletes the pulse temporal profile through the emergence of an ionization front near the focus point z c . In the MPI limit, the plasma equation ∂ t ρ = Γ|ψ| 2K can be integrated with the ansatz (58) completed by the Gaussian shape e −t 2 as where Erf(x) = (2/ √ π) x 0 e −u 2 du denotes the error function. As the beam width vanishes (L → 0), the plasma coupling term F = −ρψ can then sharply compete with the Kerr contribution. The plasma response arises like a defocusing plateau with the step function Erf( √ 2Kt) + 1, as the field intensity reaches its maximum. All time slices belonging to the interval t > 0 are defocused. At negative times, the pulse continues to self-focus and feeds plasma defocusing until it forms a thin time slice located near t = t * ≃ −(ln P in /P cr ) 1/2 at powers close to critical [17,68]. Thus, plasma generation is accompanied by a sharp duration shortening of the pulse. When GVD and MPA come into play, the leading peak becomes unstable and refocusing of the trail can occur for high enough input powers. With MPA, the leading component is partly damped in intensity. Consequently, the electron density attains much lower levels and permits the emergence of a trail. To sum up, plasma defocusing has two characteristic signatures: (i) MPI shortens the pulse duration near the nonlinear focus; (ii) Because ρ scales as R 2K (i. e., an almost Gaussian profile up to the power K > 1), the spatial zone of plasma defocusing takes place inside a narrow region of the intensity distribution, which creates spatial rings. This "self-guiding" mechanism is not static and pulse components are defocused to the benefit of others. As an example, Figure 7 depicts focusing/defocusing cycles affecting the peak intensity, electron density, beam radius and the distribution in space and time of a 150-fs pulse for atmospheric propagation. We can observe the rings formed just at the stage of plasma defocusing [ Fig. 7(e)] and the occurrence of a sharp trail when steepening effects (T, T −1 ) are included [ Fig. 7(d)].
Note that high-order effects may smooth the balance between Kerr and plasma nonlinearities. For instance, the Raman response in Eqs. (23) becomes all the weaker as the pulse is short. This affects the critical power for self-focusing like [151]. The function G(t) follows from computing the integral h(t − t ′ )|ψ(t ′ )| 2 dt ′ ∼ G(t)|ψ| 2 over the pulse profile. Unlike long pulses, for which G(t) → 1, the noninstantaneous nonlinearity reduces the effective value of the Kerr index and results in delaying the nonlinear focus from which a filament emerges [5,38,97]. Besides, the cubic polarization ∼ E 3 allows to generate third harmonic [Eq. (38)], which can soften the growth of the pump component. If we impose the limits T 0 → 1 and T q → 1, Eqs. (48a) and (48b) describe TH generation for slowlyvarying envelopes having narrow-spectral bandwidths ∆ω j /ω j ≪ 1 (j = ω 0 , 3ω 0 ). The coupling is insured by cross-phase modulation (XPM) and four-wave mixing (FWM) induced by the cubic nonlinearity. In dimensionless form, the propagation equations then read like Eqs. (51) where the function F for the two components U ω and U 3ω must be adapted as with ∆v ′ = t p ∆v/4z 0 and ∆k ′ = 4z 0 ∆k. In self-focusing regimes, the third harmonic is mainly issued from the source (FWM) term ψ 3 ω and it usually contributes by a little percentage to the overall beam fluence [4]. Despite the smallness of the TH field, this component can act as a saturable nonlinearity for the pump component when ∆k < 0. The term containing this mismatch parameter has an order of magnitude comparable with the FWM and XPM terms, ψ 3ω ∼ −ψ 3 ω /3(|∆k ′ | + 2|ψ ω | 2 ), while all other contributions driving the TH component are neglected. Once inserted into the pump wave equation, this so-called "cascading" limit [29] introduces a saturable quintic nonlinearity that lowers the peak intensity of the fundamental pulse and enhances the propagation over longer distances [23]. Figure 8 illustrates the growth of TH intensity for a meter-range propagation of one filament (P in /P cr = 4, w 0 = 0.5 mm and t p = 127 fs) issued from an infrared Gaussian pulse propagated in parallel geometry in air. A "two-color" filament coupling the infrared (800 nm) and ultraviolet (266 nm) components propagates with a clamped intensity whose maxima are less than those of the carrier wave alone. This explains the lower maximum intensity reached through the FME model (12) in Fig. 3(a).

Supercontinuum generation.
The simplest equation i∂ z ψ = −|ψ| 2 ψ admits the exact solution ψ = ψ 0 e i|ψ0| 2 z where ψ 0 ≡ ψ(z = 0). This describes a self-induced phase modulation experienced by the optical field along z, which is called "selfphase modulation" (SPM). This intensity-dependent phase shift is responsible for spectral broadening by virtue of the relation ∆ω = −∂ t arg(ψ) [1,136]. Because the frequency spectrum is expanded by the nonlinearity, SPM leads to supercontinuum and white-light emissions, as the wave intensity strongly increases through the selffocusing process. With plasma generation, this phenomenon can be described by employing the following approximations. First, we assume that the MPI response behaves like a static density plateau [see, e. g., Fig. 7(f)], so that the ratio ρ/|ψ| 2 is either zero or close to unity with ∂ z (ρ/|ψ| 2 ) ≃ 0. Second, omitting GVD and MPA we retain self-steepening to the detriment of space-time focusing, i. e., we neglect diffraction in Kerr-dominated regimes. When we model the pulse initially located at t = t 0 by A 2 = A 2 0 /cosh(τ ), τ = (t − t 0 )/τ 0 , the solving for the amplitude and phase of solutions ψ = Ae iϕ to Eq. (51a) yields the spectral broadening [166] ∆λ where Q ≡ zA 2 0 /t p ω 0 τ 0 is linked to self-steepening. As long as Q ≪ 1 and in the absence of MPI, variations in wavelength ∆λ/λ 0 ≈ −Qsinh(τ )/cosh 2 (τ ) starting with t 0 = 0 represent the early broadening through SPM (symmetric in frequency). When MPI forms a defocusing plateau, ∆λ/λ 0 then becomes larger in the region where ρ = 0 (non-defocused leading pulse), than when ρ → A 2 (defocused trail). Consequently, as the beam reaches the first focus point z c , the dominant part of the pulse is the front edge (ρ = 0) and MPI creates a primary redshift [94]. In the trail, full plasma coupling (ρ/A 2 → 1) limits the spectral enlargement to the opposite side. However, self-steepening induced at increasing Q inhibits the MPI redshifting and displaces more the spectrum to the blue wavelengths. This creates an asymmetric spectral broadening with a prominent blueshift ∆λ < 0 for τ > 0 [59,160,166]. These behaviors justify the differences in the insets of Figs. 3(b) and (c) for Gaussian pulses filamenting in air at 800 nm.
3.5. Multiple filamentation. In nonlinear optics, an intense wave propagating in a focusing Kerr medium can break up into several spots from perturbations of its initial distribution [24,30]. The resulting small-scale structures are usually named "filaments". Because it mostly concerns spatial distortions, this process of "multifilamentation" can be understood from Eqs. (51), where temporal variations are discarded. (2 + 1)-dimensional models freezing temporal dependencies are indeed able to yield rather good qualitative agreements with experimental measurements [22,44,46,50,73,124]. For instance, by assuming that MPI mainly counterbalances Kerr focusing at a dominant time slice t = t c (z), the field envelope can be decomposed as ψ =ψ(x, y, z) × e −[t−tc(z)] 2 /τ 2 with constant extent τ . Plugging the above expression into Eqs. (51), computing the expression of ρ and averaging the result over the entire time domain supplies the extended NLSE for the spatial profileψ: where the coefficient α averages the Raman delayed response. This model does not formally depend on the longitudinal location of the time slice t = t c (z). It faithfully reproduces experimental fluences of multifilamented pulses by choosing τ = 0.1 [139]. Multiple filamentation is currently described in its early stage by a perturbation theory involving a steady-state solution, expressed asψ s ( r, z) = φ( r)e iλz . Defined in the limit of no dissipation (ν → 0), φ satisfies the differential equation and λ = const. Stability of φ is investigated from perturbations v + iw with small real-valued components (v, w) acting against this stationary mode. Linearizing Eqs. (63) yields the eigenvalue problem [89] where L 0 and L 1 are self-adjoint operators with f ′ (φ 2 ) = ∂f (u)/∂u| u=φ 2 . Combining Eqs. (65), we then obtain ∂ 2 z v = −L 0 L 1 v, from which different kinds of instabilities may be investigated.
For instance, the modulational instability (MI) involves oscillatory perturbations with an exponential growth rate, that split the beam envelope approximated by a background uniform solution [24]. Perturbative modes are chosen as v, w ∼ cos (k x x) cos (k y y)e γz and they apply to a plane wave φ which satisfies ∇ 2 ⊥ φ = 0 and λ = f (φ 2 ). The growth rate γ is then given by Plane waves are unstable with γ 2 > 0 in the range 0 < k ⊥ < √ 2A and the maximum growth rate γ max = A is attained for k ⊥ = k max = √ A. This instability promotes the beam breakup into small-scale filaments periodically distributed in the diffraction plane with the transversal spacing λ mod ≃ 2π/k max and longitudinal length ∼ γ −1 max . For more accuracy, plane waves can be replaced by the soliton modes of the NLSE (64). The resulting instability appears when a soliton φ is perturbed by oscillatory modulations developing along one transverse axis. Perturbations are, e. g., local in x and they promote the formation of bunches periodically distributed over the y axis. Numerical computations are then required for solving the spectral problem (65) [2,134,169].
Besides, broad beams usually break up upon ring diffraction patterns after a stage of azimuthal instability [49]. To model this behavior, the Laplacian in Eq. (63a) must be rewritten as ∇ 2 ⊥ = r −1 ∂ r r∂ r + r −2 ∂ 2 θ , where θ denotes the azimuthal angle. Within a first approximation, unstable modes v, w ∼ cos(M θ)e γM z with index number M break up a spatial ring, which is modeled by a uniform background solution φ lying on a circular path with length s =rθ and mean radiusr. Eqs. (65) then yield a growth rate formally identical to Eq. (66), where k ⊥ must be replaced by the ratio M/r [12,142]. The maximum number of modulations on the ring is provided by the integer part of M max =r √ A. In this regards, solitary waveforms having a ring-shaped distribution may also exhibit a phase singularity with an integer number of windings, m (topological charge). Such structures are termed "optical vortices". They convey a constant orbital angular momentum [86] and are experimentally produced by means of phase masks and Figure 9. Multifilamentation patterns in Kerr focusing regime at powers above critical (P in ≤ 30P cr ) for (a) a Gaussian beam (x 0 = 0.35, y 0 = 0, ǫ = 0.1) and (b) an annular ring with flat phase (x 0 = 1.7, y 0 = 0, ǫ = 0.01). These two beams undergo MI from a local defect with relative amplitude ǫ located at the coordinates (x 0 , y 0 ). (c) Vortex-shaped beams with charge |m| = 2 initially perturbed by a 10% random noise.
holographic techniques [42,150]. For cubic media, they undergo azimuthal instability, which make vortex solitons decay into M max ∼ 2|m| + 1 filaments [156,159]. For cubic-quintic nonlinearities, optical vortices may be linearly stable for |m| = 1 [14,125,152], but they decay into M max ∼ 2|m| filaments otherwise [57]. Figure  9 shows the multifilamentation of several beam distributions including an m = 2 vortex that breaks up into 2|m| + 1 filaments in a cubic medium. MI is triggered either from an initial defect or from random noise [19]. Note the robustness of the Gaussian profile, whose shape differs from a plane wave.
Modulational instabilities rather affect the early Kerr stage in the beam propagation. At later distances, the resulting filaments become fully nonlinear and they can mutually interact [20,84]. For instance, in air, many collapsing spots can be nucleated from TW beams, in which nonlinear dissipation consumes small energy per filament [108]. The filaments are then able to merge and relaunch recurrent collapse events at further distances (see for instance Fig. 10) [69,153]. For broad pulses with enough power (e. g., > 100 P cr ), the self-focusing distance is close to the MI filamentation distance z fil ∼ 1/A ∼∼ 1/n 2 I 0 , which varies like 1/P in [31,50]. 4.1. Equilibrium and robustness. Summarizing the previous properties, the appropriate scenario for filament propagation is known as the "dynamical spatial replenishment". During the focusing stage, the beam generates a narrow plasma that strongly defocuses the trailing part of the pulse and creates one leading peak. Once the intensity decreases enough (via, e. g., MPA), plasma generation turns off. The back pulse can then refocus, which produces a two-spiked temporal profile [109,110]. In space, plasma defocuses only the inner part of the beam intensity, so that the trailing edge decays into spatial rings [72]. With increasing distances, the front pulse decreases, the spatial rings coalesce under Kerr compression and allow refocusing of the trail (Figs. 7). At high enough powers, self-channeling is supported by several focusing/defocusing cycles triggered from a dynamical interplay between stringent peaks alternating in the temporal pulse profile [16,38,39,116,155]. "Pushed" by steepening effects, these peaks are responsible for the formation of optical shocks in the medium [5,32].
In agreement with this fundamental scenario, estimations for peak intensities (I max ), electron densities (ρ max ) and filament radius (L min ) can be deduced from equating diffraction, Kerr and ionization responses. This yields the simple relations represents the effective Kerr index when the Raman-delayed response is maximum over the initial pulse profile and W (I max ) reduces to W (I max ) = σ K I K max in MPI regime. From the parameter values reported in Table 1 (see Sec. 2), the above relations supply filament diameters of the order of 150 µm in air at 800 nm and 10−20 µm in dense transparent media (dielectrics, water), which agrees with current data. Note that for laser wavelengths about 1µm, the requirement (7) is always fulfilled.
Nonetheless, there exist different scenarios having attracted attention in the past few years. For instance, Di Trapani and co-workers [44,45,124] emphasized the apparent absence of free electron emission in the self-guiding of 527 nm, 200 fs pulses with ∼ 100 µm waist and powers up to 10 P cr , propagating in water within a loosely focused geometry. Numerical simulations outlined the filament reshaping into a nonlinear X (biconical) wave, driven by the combination between linear diffraction, normal GVD, self-focusing and nonlinear losses [83]. Whereas GVD cannot halt the beam collapse even in these configurations (see Fig. 6), normal dispersion plays an important role by spreading out the most intense time slices of the pulse. Instead of an extensive production of free electrons, the nonlinear losses (MPA) dominate at laser wavelengths as low as 527 nm. The percentage of MPA over the emission of free electrons augments all the more as the laser wavelength is low, so that MPA saturates the collapse and becomes a key player in dense media whenever λ 0 < 600 nm. This mainly follows from the dependences over λ 0 of both the critical plasma density ρ c ≃ 1.11 × 10 21 /λ 2 0 [µm] cm −3 and the photo-ionization rate W (I). As a result, self-focusing is stopped at lower intensities and the free electron density is maintained at weaker levels [140].
Femtosecond filaments can propagate over several meters and simply reform after hitting an obstacle. In the atmosphere, water droplets, even as large as 100 µm in diameter, do not block the filamentation and limit the energy losses to only 15% of the filament energy [40]. The filament rapidly self-heals and seems unaffected by the droplet, so that an energy balance between the filament core and the unfocused part of the beam was conjectured to explain the rebuilding of the pulse. Unfocused parts of the beam refer to the "photon bath", i. e., the low-amplitude background surrounding the filament core. Three-dimensional simulations enclosing the filament inside a circular disk of 300 µm in diameter clearly showed, however, that the beam components outside this zone play no significant role in the filament reconstruction [138]. Some of these results are summarized in Fig. 11. Recent experiments using pinholes to block the photon bath confirmed these theoretical predictions. The energy reservoir actively feeding the filament core was measured within a diameter of 220−440 µm and was found to contain up to 50% of the pulse energy [98]. Similar properties were earlier reported from experiments in water [44], where selectively blocked portions of the beam allowed to clearly evidence this reservoir.

4.2.
White light generation. A major feature of femtosecond filaments is the broad supercontinuum developed by these nonlinear objects. When discarding steepening operators, frequency variations are dictated by where ϕ( r, t) is the phase of the field envelope, which varies with the superimposed actions of the Kerr and plasma responses. Near the focus point, only the front leading edge survives from this interplay and a redshift occurs. At later distances, second focusing/defocusing sequences relax the spectra to the blue side while the redshifted components decrease in intensity. A salient blueshift around the central frequency follows from self-steepening that creates a shock edge at the back of the pulse and amplifies a "blue shoulder" in the spectrum [5]. This spectral dynamics readily follows from Eq. (62) and was earlier reported in Sec. 2. Besides Kerr and plasma effects, the material dispersion can play a relevant role in condensed media, for which the GVD coefficient is usually high. Kolesik et al. [80,81] compared the extension of supercontinuum wings generated in air (k ′′ = 0.2 fs 2 /cm at 800 nm) and in water (k ′′ = 500 fs 2 /cm at 400 nm). Numerical simulations accounting for the complete linear dispersion relation displayed evidence of that chromatic dispersion becomes a key player (for dense materials above all) in determining the spectral extent of supercontinuum generation.
On the other hand, spectral broadening becomes also enhanced by harmonic generation. At 800 nm, the coupling of TH with the infrared pump produces a "twocolored" femtosecond filament from the threshold intensity I(ω) ≥ 2 × 10 13 W/cm 2 [7,167]. Along this process, the pump wave injects part of its energy into the third harmonic. The amount of pump energy depends on the linear wave vector mismatch parameter ∆k = [3k(ω) − k(3ω)] −1 fixing the coherence length L c = π/|∆k|. The smaller the coherence length, the weaker TH fields, since most of the pump energy is periodically injected and depleted by the pump. In self-channeling regime, the balance between TH, pump wave nonlinearities and the linear mismatch parameter makes the phase difference ∆ϕ = 3ϕ(ω) − ϕ(3ω) be clamped at a constant value equal to π [4]. Along meter-range distances, the two-colored filament is capable of covering longer scales than an infrared pulse alone (see, e. g., Fig. 8), as the TH component limits the growth of the pump wave [23].
Manipulating intense, powerful filaments centered at 800 nm results in supercontinuum generation extending down to 230 nm in air and yielding a continuous spectral band of UV-visible wavelengths [149]. Numerical simulations evidenced that UV-visible spectral broadening is created by the overlap of the supercontinuum of the pump and TH broadening. In this process, the enlargement of the pump to the blue wavelengths caused by the T, T −1 operators prevails [see Fig. 3(a)].
4.3. Filamentation vs. laser wavelength. Let us here comment on the role of the laser wavelength. Besides the Rayleigh length z 0 , several physical parameters become modified when varying λ 0 . Among those, the Kerr index n 2 decreases when λ 0 is increased [1,113]. The photo-ionization rate and the number of photons consumed by MPI transitions ∼ U i / ω 0 also vary with λ 0 . In addition, the value and sign of the GVD parameter, k ′′ = ∂ 2 k/∂ω 2 | ω=ω0 , are modified by changing λ 0 .
Changes in supercontinuum generation have numerically been addressed for laser wavelengths as large as 1550 nm [3]. The latter corresponds to eye-safety and is particularly important for remote-sensing applications. For these wavelengths, the GVD parameters for both pump and TH frequencies are close to each other. However, their wavevector mismatch strongly decreases from −5 cm −1 (800 nm) to −0.65 cm −1 (1550 nm), whereas the temporal walk-off parameter ∆v increases from 0.4 to 2.0 cm/s, respectively. Hence, the coherence length augments at large wavelengths, which allows higher TH intensities and better conversion efficiency.
More surprising features occur when the laser wavelength is selected in such a way that k ′′ becomes negative and leads to anomalous dispersion. In that case, the pulse undergoes temporal compression, in addition to the spatial Kerr focusing. This happens, e. g., in fused silica at λ 0 = 1550 nm, for which k ′′ ≃ −280 fs 2 /cm (see Table 1). Recent experiments at this wavelength have shown that collapse events looked "extended" along the z axis, unlike the localized events promoted by normal GVD. Plasma halts the collapse at powers above critical, but anomalous GVD continues to transfer energy into the collapse region, resulting in the formation of longer filaments before the beam eventually defocuses [111]. The pulse can thus remain confined along several Rayleigh lengths. Strong temporal compression produces optical spots whose duration is shrunk to the few-cycle limit [21,96]. For comparison, Fig. 12 displays the temporal distributions of a femtosecond filament created in fused silica at 790 nm (k ′′ = 370 fs 2 /cm > 0) and at 1550 nm (k ′′ < 0). At 790 nm, normal GVD stretches the pulse along the time direction. Reversely, at 1550 nm, anomalous GVD compresses it temporally while the pulse is shifted to positive instants through third-order dispersion and self-steepening. This latter property follows from the evolution of the centroid integral which is computed from Eq. (35) under the limits x K = 0, T −1 ≃ 1 − (i/ω 0 )∂ t and up to fourth-order dispersion. Here, the first three terms, which correspond to space-time focusing, self-steepening and third-order dispersion, respectively, are strictly positive.

Pulse shortening in gases.
In the last decade, powerful techniques of ultrabroadband dispersion control have been developed, in order to compress pulses to durations of a few optical cycles only (τ o.c. = λ 0 /c). Nisoli and co-workers [115] manipulated the spectral broadening of ∼ 1 mJ, 20-fs pulses along hollow fiber filled with atomic gases. Through an appropriate combination of SPM and gas dispersion managed by external chirped compensation systems, pulses as short as 4.5 fs could be produced at 800 nm. The right combination could at that time be satisfied by a delicate tuning of the pressure parameter. For a gaseous medium with pressure p, the coefficients affecting Eq. (35) indeed vary as follows so that an optimal compression of the pulse was reached by varying the cell length and the local pressure p. Laser parameters were kept with an input power below the collapse threshold, in order to avoid plasma formation. However, plasma generation may not constitute a drawback for pulse compression. For instance, Chen et al. [33] achieved self-compression from 50 to 20 fs, by making 800 nm pulses with energy < 1 mJ pass through a BK7 glass plate. Measurements revealed that compression takes place as long as the front pulse focuses while plasma cuts off the trail part.
Self-compression experiments are nowadays directed to the single-cycle regime. This requires an accurate control of the pulse temporal distribution, in such a way that one peak compressed down to a few fs can be isolated [32]. To realize this challenge, Hauri et al. used a configuration of two cascaded gas-filled cells with intermittent dispersion compensation and different pressures for producing sub-mJ light pulses with durations down to 5.7 fs at ∼ 800 nm [65]. Alternatively, Stibenz et al. recently demonstrated an efficient pure self-compression obliterating the need for any kind of dispersion compensation, pressure gradients, or capillaries for beam guiding [147]. The experimental setup involves beams in convergent geometry with energy up to 5 mJ, w 0 = 11 mm, f = 1.5 m, durations of 45 fs and input powers equal to 5 times critical. The beam propagates inside an 1-m long cell whose center is positioned at the location of the geometrical focus. Output pulses reach durations down to about 10 fs and develop a strong blueshift, as reported in Figs. 13(a,b). Figure 13(c) details the underlying mechanism, reproduced numerically from Eq. (35). First, near z c ≃ 1.45 m, the front pulse focuses and plasma depletes its rear part. Second, at the linear focus (z ≃ 1.5 m), the plasma density decreases, and the back pulse rises again. Finally, from z ≃ 1.6 m, only the components close to center, keeping intensity levels above 10 13 W/cm 2 , are effectively trapped in the filament, whereas the temporal wings diffract rapidly. A robust, temporallycompressed structure of 11 fs forms in the filament core [ Fig. 13(d)]. The pulse center relaxes to a narrow waveguide, preserving its energy at low plasma density levels (< 10 14 cm −3 ) after the linear focus. Fig. 13(e) shows the on-axis intensity spectra, corresponding to this temporal compression [141]. This experimental scheme revealed for the first time the existence of ultrashort "light bullets", which should further provide novel few-cycle optical sources.

5.2.
Pulse propagation in dense media. Investigations on ultrashort pulse propagation in dielectrics (e. g., fused silica) started from the early 2000's. Tzortzakis et al. [155] observed the 1-cm long self-channeling of a 160-fs focused pulse conveying 3 critical powers at 800 nm. The pulse traveled across a fused silica sample like a narrow waveguide with ∼ 20 µm waist. Local heating and damage by local intensities > 10 TW/cm 2 were avoided by making the sample move in the (x, y) plane. Earlier, Brodeur and Chin reported a similar behavior accompanied by a wide supercontinuum in glasses [28]. As an example, Fig. 14(a) shows a transverse photograph of the self-guided filament measured at input energy of 2 µJ in fused silica. Fig. 14 Besides solid materials, femtosecond filaments can also be self-guided in water, as examined by Dubietis and collaborators a few years ago [45]. By launching a Gaussian beam centered at 527 nm into a water-filled cuvette in focused geometry, a single filament formed and was capable of covering up to 4 cm along the propagation axis, while keeping a diameter of a few tens of µm. The experiment revealed that the filament dynamics was not sustained by a balance between Kerr-induced self-focusing and plasma-induced defocusing. As already signaled in the previous section, numerical simulations outlined, instead, the spontaneous reshaping of the beam into a Bessel-type X-wave fulfilling the requirements of minimum nonlinear losses, maximum stationarity and localization. The X shape of the beam lines appears when looking at the angular frequency spectrum.

5.3.
Filaments in the atmosphere. Air was historically the "birthplace" for the science of femtosecond filaments, whose illuminating example has been given in diameter and ∼ 1 mJ in energy at infrared wavelengths (800 nm) [27]. Emitting an impressive spectrum, they couple with electron densities of a few times 10 16 cm −3 [91] for peak intensities of about 5 × 10 13 W/cm 2 [75]. We review below major achievements of their applications.
The propagation of high-power (TW) femtosecond laser pulses in air has attracted considerable attention from the pioneering observation of the white-light supercontinuum beyond 10 km [164]. A number of important practical applications followed, including remote sensing [127], directed energy delivery and artificial lightning [90]. To achieve high intensities at remote distances, a negative frequency chirp can be introduced in the laser pulse [6,162]. The pulse then undergoes a temporal compression as it compensates dispersion along the propagation axis. Chirping effects can be measured by evaluating differences in the supercontinuum emission with and without pulse chirping from altitudes beyond 18 km [63,117,131]. The onset distance and longitudinal extent of ultrashort filaments can also be conditioned by the initial focusing geometry (f < +∞). Most experiments emphasize the use of rather large ratios f /w 0 , in order to avoid an immediate plasma defocusing in tightly focused configuration.
Open-air experimental campaigns make it possible to observe optical focal spots formed by many filamentary cells for input powers containing several thousands of critical powers in air. The 2D reduced model (63) using the experimental fluence as initial condition actually reproduces the evolution of the multifilamentation pattern over long scales, saving computational resources when one neglects the temporal dimension. As an example, Fig. 15 shows the overall envelope of the Teramobile bundle launched in parallel geometry. The initial pulse contains 700 critical powers and its transverse profile is scanned at different propagation distances. Filaments rise from the initial beam defects, form a crown of dots growing from the diffraction ring and then excite clusters of cells. Some of these clusters are identified by the labels (1), (2) and (3). They are faithfully reproduced by the numerics. The filament number remains in the order of P in /P fil , where P fil ≃ 3 − 5P cr is the power in one filament. At large distances, the primary brightest spots decay into secondary filaments by exchanging power through the energy reservoir formed by the background field [22].
However, the filamentation pattern becomes severely modified when the focusing geometry is changed. Figure 16    produce more than one hundred filaments in parallel geometry. Numerical computations for this configuration indicate that many filaments are created before the linear focus, but they fuse into three strings of light acquiring a high directivity afterwards [139].
The white-light continuum generated in air by ultrashort laser pulses has first been characterized in the visible [8,114]. J. Kasparian et al. next investigated the infrared region [76]. Experimental spectra around λ 0 = 800 nm exhibit universal features, some of which can be refound in Fig. 3, mostly at low wavelengths. As shown in Figure 17, the continuum is very broad, extending at least to 4.5 µm. An almost exponential decay over 4 orders of magnitude up to 2.5 µm is observed, followed by a slower decay of one order of magnitude only. Variations in the input energy can make spectral intensity change by one decade in some spectral regions. Nevertheless, the spectral shape developed from various setups remains similar and their large width finds direct uses in LIDAR technology [106,163]. In this scope, a laser pulse is launched into the atmosphere. The backscattered light is collected on a telescope and detected as a function of time, with a typical resolution of 1-10 nanoseconds. This temporal window yields a high spatial resolution, since the flight time of the detected photons is directly proportional to the distance where they have been backscattered. Such a spatial resolution, combined with the possibility of sweeping the laser beam, provides two-and three-dimensional maps of measured atmospheric species.
LIDARs using femtosecond pulses (or so-called "FemtoLIDARs") combine several advantages. The spectral bandwidth covered by fs filaments spans at least from 230 nm in the ultraviolet to 4.5 µm in the mid-infrared (Fig. 17). This includes the absorption band of water between 1.8 and 2.5 µm together with that of volatile organic compounds (VOCs) between 3 and 3.5 µm. In addition, the flat continuum spanning from the visible down to 230 nm provides a promising light source for the measurement of trace gases that absorb in the blue or the UV, such as ozone, toluene, benzene, SO 2 or nitrogen oxides.
The ability of filaments to remotely deliver intensities as high as 50 TW/cm 2 also opens the way to innovative exploration techniques, such as the "remote filamentinduced breakdown spectroscopy" [146]. This consists of a combination of LIDAR and laser-induced breakdown spectroscopy [41], which is a versatile tool allowing the analysis of numerous materials. It relies on the local ionization of their surface by a strongly focused pulsed laser. Here, the low energy of subpicosecond laser pulses limits the heating of the sample [132], while self-guided filaments can deposit intensities higher than the ablation threshold at distances of hundreds of meters or even kilometers. In this novel technique, laser filaments are launched on a remote target, and the light emitted by an excited plasma plume is collected and analyzed through a detection system comparable to a LIDAR-based setup. The characteristic spectra provide direct information about the nature of the material. 6. Outlook. Femtosecond lasers enable robust filamentary structures to propagate high intensities over long distances in a rich variety of transparent media. These "femtosecond filaments" are characterized by micrometric sizes and they can reach ultrashort durations down to the optical cycle limit. They design new objects, that can be exploited for delivering few-cycle pulses at high power levels in the future. Ionization processes do not break this property, so that a precise control of plasma generation should push new frontiers in nonlinear optics. The basic scenario of a Figure 17. Spectrum of the white light continuum emitted by 800 nm pulses and assembled from five spectral regions. Main atmospheric absorption bands are marked above the curve. dynamical balance between Kerr focusing and plasma defocusing explains the existence and the long life of femtosecond filaments in most of propagation media. This scenario must, however, be revised at visible wavelengths in condensed materials, in which chromatic dispersion and nonlinear absorption can drive the self-guiding mechanism through an important conical emission managed by X-shaped waveforms in normally dispersive regimes. This aspect requires more investigations. Also, propagation regimes of anomalous dispersion, which promote sharp temporal compression, should be deepened from a mathematical point of view.
In the last years, impressive progresses have been made in the manipulation of high-power laser pulses launched over long distances in air. Filaments can be used to remotely deliver high intensities generating in-situ nonlinear effects. The spectrum of the white-light continuum is much broader than that expected a few years ago, especially in the ultraviolet where the mixing between the fundamental and third harmonic wavelengths together with spectral amplification caused by pulse steepening give rise to a spectral plateau extending down to 230 nm. These progresses open the way to novel applications such as remote sensing of gaseous pollutants and aerosols by ultrafast LIDAR techniques as well as probing solid samples by remote spectroscopy analysis.