{'fraction_non_alphanumeric': 0.05652674552496598, 'fraction_numerical': 0.022266752403882102, 'mean_word_length': 3.9337464954994834, 'pattern_counts': {'":': 0, '<': 17, '<?xml version=': 0, '>': 39, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 26, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Variants of the focusing NLS equation. Derivation, justification and open problems related to filamentation

The focusing cubic NLS is a canonical model for the propagation of laser beams. In dimensions 2 and 3, it is known that a large class of initial data leads to finite time blow-up. Now, physical experiments suggest that this blow-up does not always occur. This might be explained by the fact that some physical phenomena neglected by the standard NLS model become relevant at large intensities of the beam. Many ad hoc variants of the focusing NLS equation have been proposed to capture such effects. In this paper, we derive some of these variants from Maxwell's equations and propose some new ones. We also provide rigorous error estimates for all the models considered. Finally, we discuss some open problems related to these modified NLS equations.


Introduction
The cubic, focusing, nonlinear Schrödinger equation in space dimension d is given by (1) i∂ It is a canonical model for the propagation of laser beams. From a result of Ginibre and Velo [15], equation (1) is locally well-posed in H 1 = H 1 (R d ) for d = 1, 2, 3, and thus, for v 0 ∈ H 1 , there exists 0 < T ≤ +∞ and a unique solution v(τ ) ∈ C([0, T ), H 1 ) to (1) and either T = +∞, we say the solution is global, or T < +∞ and then lim τ ↑T ∇v(t) L 2 = +∞, we say the solution blows up in finite time.
Referring to conservation of the L 2 norm by the flow, (1) is said to be L 2 −subcritical if s c < 0, L 2 −critical if s c = 0 and L 2 −supercritical if s c > 0. Thus, (1) is L 2 −subcritical if d = 1, L 2 −critical if d = 2, and L 2 −supercritical if d ≥ 3. In the subcritical case, global existence (in C([0, ∞), L 2 )) holds for arbitrarily large data in L 2 . It turns out that in this case, global existence (in C([0, ∞), H 1 )) also holds for arbitrarily large data in H 1 , due to the conservation of mass and energy. In the critical and supercritical cases however, there exist stable finite time blow-up dynamics. This has been known since the 60ies using global obstructive arguments based on the virial identity (see e.g. [42]). 1 There is however a discrepancy between the blow-up results predicted by (1) and physical observations. Indeed, while the blow-up signifies a break-down of the solution v, physical observations show in many cases that lasers begin to focus according to the scenarios associated to (1) but depart from this behavior slightly before the focusing time. The reason advanced by physicists is that some physical phenomena that have been neglected to derive (1) become relevant at high intensities, and therefore near focusing. This phenomenon is called filamentation: defocusing physical phenomena are triggered at high intensities and halt the collapse of the beam. This interplay between diffraction, self-focusing, and defocusing mechanisms allow for the beam to propagate along several times the focusing distance (called Rayleigh length in optics) and the resulting structure is called filament .
Many variants of (1) have been derived in optics to take into account these additional physical phenomena and reproduce the filamentation mechanism. In many cases, it is a mathematical open problem to prove whether these additional terms prevent focusing or not, and a fortiori to understand the modification of the dynamics induced by them.
Rather than adding as usual ad hoc modifications to (1) in order to take new physical effects into account, we choose here to rigorously derive such modifications from Maxwell's equations. We then comment on some of the most physically relevant open mathematical problems that these modified equations raise and that are natural milestones towards the understanding of filamentation. These variants can roughly be classified into two groups, depending on whether they take ionization processes into account or not. Notation 1.1. In the brief presentation below, we denote by z the direction of propagation of the laser, by X ⊥ = (x, y) the transverse variables, and by ∆ ⊥ = ∂ 2 x + ∂ 2 y the transverse Laplace operator. In dimension d = 2, the variable y is omitted (and hence ∆ ⊥ = ∂ 2 x ), while in dimension d = 1, functions depend only on z (so that ∆ ⊥ = 0).

Models without ionization processes
We give below a family of variants to (1) that incorporate many physical phenomena neglected by (1). It is of course possible to look at one or several of these additional effects simultaneously. We state the equations in their most general form, starting with a family of scalar NLS equations, and then give the corresponding vectorialand more general -form of these equations. Let us therefore consider (3) iP 2 (ε∇)∂ τ v + (∆ ⊥ + α 1 ∂ 2 z )v + iα 2 v + (1 + iεα 3 · ∇) 1 + f (ε r |v| 2 ) |v| 2 v = 0, 1 Much more is known about the finite time blow up dynamics for the focusing NLS and we refer the interested reader to [38] [26] [27] [28] [29] [30] [31] [32] and references therein.
where v is a complex-valued function. Here, ε > 0 is a (small) parameter; P 2 (ε∇) is a second (at most) order, self-adjoint, positive operator; α 1 = 0, ±1; α 2 ≥ 0; α 3 ∈ R d ; f : R + → R is a smooth mapping vanishing at the origin, and r > 0. The physical meaning of these terms is commented below: (1) Nonlinearity. The cubic nonlinearity in (NLS) corresponds to a first order approximation of the nonlinear optical phenomena. At high intensities, it is often worth including some next order terms captured here by the additional term f (ε r |v| 2 ). We consider here three situations: (a) Cubic nonlinearity: f = 0.
(c) Saturated nonlinearity: f is a smooth function on R + vanishing at the origin and such that (1 + f (r))r is bounded on R + (e.g. f (r) = − r 1+r ).  (4) Off-axis variations of the group velocity. The operator P 2 (ε∇) is here to account for the fact that self-focusing pulses become asymmetric due to the variation of the group velocity of off-axis rays 2 . The operator P 2 (ε∇) is a second order, self-adjoint, and positive operator in the sense that (P 2 (ε∇)u, u) ≥ C|u| 2 * , with |u| 2 * ≥ |u| 2 2 . The norm | · | * may also control derivatives of u; we consider three cases: (a) No off-axis dependence: P 2 (ε∇) = 1, and therefore | · | * = | · | 2 . (b) Full off-axis dependence: the norm | · | * controls all first order derivatives, |u| 2 * ∼ |u| 2 2 + ε 2 |∇u| 2 2 . (c) Partial off-axis dependence: the norm | · | * controls some but not all first order derivatives. More precisely, there exists j (j < d) linearly independent vectors v j ∈ R d such that |u| 2 * ∼ |u| 2 2 +ε 2 j k=1 |v k ·∇u| 2 2 . (5) Self-steepening of the pulse. The operator (1 + iεα 3 · ∇) in front of the nonlinearity accounts for off-axis dependence of the nonlinearity, responsible for the possible formation of optical shocks. Various cases are considered here: (a) No self-steepening. This corresponds to α 3 = 0 and to the usual situation where the nonlinearity does not contain any derivative. (b) Longitudinal self-steepening. When α 3 is colinear to e z , there is a derivative in the nonlinearity along the direction z of propagation of the laser. (c) Transverse self-steepening. When α 3 = 0 and α 3 · e z = 0, there is a derivative in the nonlinearity along a direction orthogonal to the direction of propagation.
(d) Oblique self-steepening. When α 3 is neither colinear nor orthogonal to e z .
As previously said, (3) stems from a more general vectorial equation. For the sake of simplicity, we give here the equation corresponding to the cubic case (or f = 0 in (3)): where v is now a C 2 -valued function.
Remark 1.2. Equation (3) is in fact a particular case of (3) vect corresponding to initial data living on a one-dimensional subspace of C 2 . Indeed, if the initial condition to (3) vect has the form v |τ=0 = v 0 (x)v 0 with v 0 ∈ R 2 and v 0 a scalar-valued function, then the solution to (3) vect takes the form v(τ,

Models with ionization processes
In addition to the physical phenomena taken into account in (3), it is necessary at high intensities to include ionization processes for a correct description of the laser pulse. The reason why this phenomenon is singled out here is because a system of two equations must be considered instead of the single equation (3). In the most simple case (i.e. P 2 = 1, f = 0, α 2 = 0, α 3 = 0; for a more general model, see (50)), this system is given by with α 4 , α 5 ≥ 0, c > 0, and where ρ is the density of electrons created by ionization, while c g = c g e z is the group velocity associated to the laser pulse.
The system (4) does not directly compare to (3) and (3) vect ; indeed, (4) is written in the fixed frame of the laboratory, while (3) and (3) vect are written in a frame moving at the group velocity c g = c g e z and with respect to a rescaled time τ = εt. Rather than (4), the NLS equation with ionization used in the physics community is its version written in the same variables as (3) and (3) vect . More precisely, if we set and τ = εt, the equations (4) are approximated 3 by the following ones, The ionization processes taken into account by the systems (4) and (5) are: 3 The approximation lies in the equation on ρ. In the new variables, the second equation of (4) is given by ε∂τρ − cg∂zρ = εα 4 |v| 2K + εα 5ρ |v| 2 . In the physics literature, the term ε∂τρ is neglected, and this corresponds to (5).
(1) Photo-ionization. This corresponds to α 4 > 0 and K > 0 (K is the number of photons necessary to liberate one electron). (2) Collisional ionization. When α 5 > 0, a term corresponding to collisional ionization is added to the evolution equation on ρ.
Remark 1.3. (1) The coupling with ρ can of course be added to any equation of the family (3).
(2) A vectorial variant of (4) and (5) can also be derived in the same lines as (3) vect .
The rest of the paper is as follows. In section 2, we recall the Maxwell equations, we give an abstract formulation and we discuss the spaces of initial data used for the Cauchy problem. In section 3, we prove our main result about the rigorous derivation of general abstract versions of the models (3), (3) vect , and (4), (5). In section 4, we analyze the role of the various parameters in (3) and (4) or (5) In particular, we consider whether they indeed prevent the breakdown in finite time or not. We also formulate a number of interesting open problems for these modified NLS equations. Finally, an appendix contains explicit computations for a physically relevant system of Maxwell equations which show that the abstract models derived in Section 3 take indeed the form of (3), (3) vect , and (4).
where D is given in terms of the electric field E and a polarization P -modeling the way the dipole moment per unit volume depends on the strength of the electric field -by the relation and where we used standard notations ǫ 0 and µ 0 for the electric permittivity and magnetic permeability in vacuum. The evolution equations (6) go along with two constitutive laws, where ρ is the electric charge density. As a consequence of the relation ∇ · D = ρ and the second equation of (6) we get the continuity equation coupling ρ to the current density J, Introducing the speed of light in vacuum the equations (6) can also be rewritten as a set of two evolution equations on the magnetic field B and the electric field E, In order to get a closed system of equations, we still need two physical informations: (1) A description of the polarization response to the electric field E.
(2) A description of the current density J. We first address the description of the polarization response in absence of current density and then proceed to describe the modification to be made when current density is included.
2.1.1. The polarization response to the electric field. Throughout this section, we assume that there is no charge nor current density (ρ = 0, J = 0). The general case will be handled in §2.1.2 below.
There exist various ways to describe the polarization P; we use here a simple and natural model called "nonlinear anharmonic oscillator", according to which the polarization is found by solving the second order ODE (11) ∂ 2 t P + Ω 1 ∂ t P + Ω 2 0 P − ∇V N L (P) = ǫ 0 bE, where b ∈ R is a coupling constant and Ω 0 , Ω 1 > 0 are frequencies, and where V N L accounts for nonlinear effects. When such effects are neglected, the description (11) goes back to Lorentz [23] and expresses the fact that electrons are bound to the nucleus by a spring force. Nonlinearities have been added to this description by Bloembergen [5] and Owyyoung [34] and the mathematical investigation of their influence was initiated by Donnat, Joly, Métivier and Rauch [11,16] (see also [20]).
Remark 2.1. In physics books, the polarization P is often sought as an expansion where the operator χ 1 is called the linear susceptibility of the material, while for j > 1, the operators χ j are the j-th order nonlinear susceptibilities. It is easy to check that the linear susceptibility corresponding to (11) is given by the nonlocal (in time) operator where we used the Fourier multiplier notation, (ii) Cubic/quintic nonlinearity: |P| 6 and therefore ∇V N L (P) = a 3 |P| 2 P − a 5 |P| 4 P.
(iii) Saturated nonlinearity: there exists a function v sat : and v ′′ sat bounded on R + and such that for instance, one can take v sat (r) = a 3 2 and is therefore the same at the origin as in (ii) above, up to higher order terms (seventh order terms here).
We show in Appendix A.1 that Maxwell's equations can be put under the following dimensionless form 4 for all the nonlinearities considered in Example 2.1, where γ, ω 0 , ω 1 , r and p are constants, 0 < ε ≪ 1 is a small parameter (the ratio of the duration of an optical cycle over the duration of the pulse, see Appendix A.1), while f is a smooth function vanishing at the origin.
2.1.2. The case with charge and current density. The main mechanism at stake in laser filamentation is certainly the local ionization of the medium: once a powerful self-focusing laser beam reaches high enough intensities, it ionizes the medium around itself. It leaves behind a narrow channel of plasma, hereby causing local defocusing that prevents blowup. Taking current density into account, we come back to the set of equations (10)- (11), and a physical description of the current density J is needed. This current density has the form where J e and J i are respectively the free electron and ionization current densities.
-Free electron current density. Partial ionization of the material medium by the laser generates free electrons, with charge q e (= −1.6 × 10 −19 C). This induces a free electron current density J e = q e ρ e v e , where ρ e is the electron density, and v e is the electron velocity. A rough 5 , but standard model in nonlinear optics is to take (see [4] and references therein), where k l and ω l are the laser wave number and pulsation respectively, with |∂ t E 01 | ≪ |ω l E 01 | (slowly varying envelope approximation); the polarization current is then taken under the form where m e is the electron mass. The drawback of this model is that it assumes that the electric field and the current density field can be written at leading order as wave packets (i.e. are given under the form (14)- (15)). In particular, it does not provide any relation between the current density J and the electric field E that could be used in Maxwell's equations (10). We therefore propose here such a relation, namely, where k l = |k l |, and H is the regularization of the Hilbert transform given by the Fourier multiplier Quite obviously, this is consistent with the usual model (15) since this latter is recovered at leading order when the electric field is a wave packet under the form (14).
Finally, the evolution of the electron density ρ e is given by a source term S representing external plasma sources. Taking into account photo-ionization and collisional ionization, but neglecting electron recombination (see for instance [3] for richer models), we have where the intensity is I = |E| 2 and ρ nt is the constant density of neutral species.
In the regime considered here 6 , ρ e is negligible compared to ρ nt and the photoionization rate W (I) takes the form for some constant coefficient σ K > 0 and with K > 1 the number of photons needed to liberate one electron. The collisional ionization cross-section σ depends 5 This approximation can be deduced formally by assuming that ions are at rest and that electron motion is described by the compressible Euler system (see for instance [3]). Neglecting electron collisions, such a model yields ∂tJe = q 2 e me ρeE which formally yields (15) assuming that E is as in (14) and that ρe is not oscillating at leading order. It would of course be interesting to provide a rigorous justification to these approximations. 6 For higher intensities, electrons can tunnel out the Coulomb barrier of atoms, and W (I) is modified.
on the laser frequency, and U i is the ionization potential. Summing up, we get the following expression for the free electron current J e and ρ = ρ e , -Ionization current density. It is also necessary to take into account losses due to photo-ionization. We therefore introduce a ionization current density J i such that J i · E represents the energy lost by the laser to extract electrons (per time and volume unit). More precisely, J i · E is equal to the energy necessary to extract one electron (given by the ionization potential U i ) multiplied by the number of electrons per time and volume unit (given by ∂ t ρ). Using the second equation of (17), this gives We show in Appendix A.2 that after nondimensionalization, the set of equations (10)-(13)-(11)-(17)-(18) (for the nonlinearities considered in Example 2.1) becomes, with the same notations as in (12) for the constants γ, ω 0 , ω 1 , r and p, the small parameter ε, and the function f . In addition, we have here constants c 0 , c 1 , c 2 ≥ 0, and we also recall that the definition of the regularized Hilbert transform H is given in (16 where U is a R n (n ≥ 1) valued function depending on the time variable t and the space variable x ∈ R d (d ≥ 1), The operator A(∂) is defined as

10ÉRIC DUMAS, DAVID LANNES AND JÉRÉMIE SZEFTEL
where ∂ j is the differentiation operator with respect to the j-th space coordinate. The matrix A 0 has size n × n, p is a positive number, and the following assumption is made on the matrices A j and E, and on the nonlinearity F .
The matrix E is a constant coefficient n×n, real valued, skew symmetric matrix. (iii) There exists a smooth mapping f : R + → R vanishing at the origin, a real number r > 0, a quadratic form Q : C n → R + and a trilinear symmetric mapping where, as in § 2.2.1, U is a R n -valued function, whereas W is an R-valued function of (t, x) ∈ R × R d . The matrices A j and E, as well as the nonlinearity F , satisfy Assumption 2.1. Concerning the other coefficients of the system, we assume the following.
The constant c is positive. (iii) There exists two real, positive constants c 1 and c 2 and an integer K ≥ 1 such that Remark 2.3. As in Remark 2.2 for the case without ionization, we can put Maxwell's equation with ionization terms (19) under the abstract form (21). Using the same notations as in Remark 2.2, the matrix C 1 is the projection matrix such that C 1 U = E, and w = ρ.
2.3. The Cauchy problem. We are considering initial conditions that correspond to laser pulses. In the case without charge nor current density (equation (20)), they are fast oscillating wave packets slowly modulated by an envelope, where k ∈ R d is the (spatial) wave-number of the oscillations. Taking charge and current density into account (equation (21)), we need to provide initial conditions for W; since we are interested here in the situation where this quantity is created by the laser when it reaches high intensities near self-focusing, we take these initial conditions to be initially zero 7 for the sake of clarity.
The evolution equation (20) (as well as (21)) being of semilinear nature, it is natural to work with Banach algebra in view of a resolution by Picard iterations. Throughout this article, we assume that u 0 ∈ B, with (the so called Wiener algebra, which is better adapted than H t0 (R d ) n to handle short pulses, see [9,21]). In both cases, B is stable by translations in Fourier space (this ensures that if u 0 ∈ B in (22) then U |t=0 ∈ B) and is a Banach algebra in the sense that ∀f, g ∈ B, f · g ∈ B and |f · g| B |f | B |g| B .
For all k ∈ N, we also define endowed with its canonical norm. We are interested in deriving asymptotics to the solution formed by (20)- (22), with initial envelope u 0 ∈ B, and more generally (21)-(22)-(23), if we want to be able to handle ionization processes in nonlinear optics. This requires a further assumption on the nonlinearity F , namely that F acts on B and is locally Lipschitz. Assumption 2.3. In addition to (iii) of Assumption 2.1, the mapping F satisfies, uniformly with respect to ε ∈ [0, 1), 3 is always satisfied (by Moser's inequality); when B = W (R d ) n , the assumption holds for analytic nonlinearities. 7 One could more generally and without supplementary difficulty consider non-oscillating initial conditions for the charge density W

Derivation of NLS-type equations
The Schrödinger approximation takes into account the diffractive effects that modify over large times the propagation along rays of standard geometrical optics. These diffractive effects are of linear nature and are known [11,17,16,36,20] to appear for time scales of order O(1/ε) for the initial value problem formed by the linear part of (20) and (22). This is the reason why we are interested in proving the existence and describing the solutions to the (nonlinear) initial value problem (20)-(22) over such time scales.
For the sake of simplicity, the initial value problem (20)- (22) (no ionization) is first considered. Up to minor modifications, the results of § §3.1, 3.2 and 3.3 are known [9]; we reproduce them here because they are necessary steps to derive the family of NLS equations (3) and, for the sake of clarity, their proof is sketched in a few words. The general idea to derive the Schrödinger equations of §3.5 with improved dispersion was introduced in [9] but the computations are carried further here. We then derive in §3.6 a new class of models with derivative nonlinearity for which a local well-posedness result is proved. When applied to Maxwell's equations, these derivative nonlinearities yield the so-called self-steepening operators; to our knowledge, this is the first rigorous explanation of these terms.
3.1. The profile equation. We show here that under reasonable assumptions on F , solutions to the initial value problem (20)-(22) exist for times of order O(1/ε) and that there can be written under a very convenient form using a profile U , with U (t, x, θ) periodic with respect to θ and for any ω ∈ R, provided that U solves the profile equation Here, we used the notation where U solves the profile equation (25).
Proof. The proof is a slight adaptation of the one given in [9] in the trilinear case; consequently, we just give the main steps of the proof. Quite obviously, a solution U to (20)-(22) is given by (24) ) is a Banach algebra; moreover the evolution operator S(t) associated to the linear part of (25), is unitary on H k (T; B) (thanks to point (i) and (ii) of Assumption 2.1). One can therefore construct a (unique) solution to (25) by a standard iterative scheme Indeed, one has thanks to Assumption 2.3, thanks to the ε in front of the integral. An estimate of the same kind is valid for a difference of iterates, by point (ii) of Assumption 2.3. By a fixed point argument, this ensures that the sequence converges to a solution on [0, T /ε] for some T > 0 independent of ε. Uniqueness then follows classically from an energy estimate on the difference of two solutions.

3.2.
The slowly varying envelope approximation. The slowly varying envelope approximation (SVEA) consists in writing the profile U under the form plugging this approximation into the profile equation (25) and keeping only the first harmonic in the Fourier expansion yields easily (writing u = u env ), Denoting D = −i∇, we observe that where the last notation is of course consistent with (26).

14ÉRIC DUMAS, DAVID LANNES AND JÉRÉMIE SZEFTEL
As a consequence of these computations, we see that in order for (28) to hold, it is necessary that u = u env satisfies the envelope equation As implicitly assumed by omitting the fast oscillating scale in the argument of the envelope function u env (t, x), the envelope must not contain any fast oscillation. However, • The singular part of the linear term in (30) creates fast oscillations with frequencies ω − ω j (k), where the ω j (k) stand for the eigenvalues of L(0, k). • The nonlinearity creates other oscillations that may resonate with the linear propagator. There is one way to avoid the catastrophic effects of these two scenarios. Choosing ω = ω j (k) for some j and assuming that, up to O(ε) terms, the initial envelope u 0 is contained in the corresponding eigenspace prevents the creation of oscillations by the linear propagator. This is the polarization condition. The nonlinearity will however create harmonics of the main oscillation k · x − ω j (k)t and it is necessary to make a non resonance assumption. What is called characteristic variety in the assumption below is the set C L ⊂ R d+1 defined as Let us also recall that we assumed that the nonlinearity is under the form with f (0) = 0, Q a quadratic form and T a trilinear symmetric mapping. If U is a monochromatic oscillation, the nonlinearity εF (ε, U ) creates third harmonic with size O(ε), a fifth harmonic (if f ′ (0) = 0) with size O(ε 1+r ), etc. The non-resonance condition stated below holds for the (2p + 3)-th harmonics, for all p ≥ 0 such that pr < 1 (the contribution of higher harmonics is small enough to be controlled even if it is resonant).
Notation 3.1. We denote by π j (k) (j = 1, . . . , m) the eigenprojectors of the eigenvalues ω j (k) of A(k) + 1 i E; in particular, we have Example 3.2. For Maxwell's equations, it is shown in Appendix B that Assumption 3.1 is satisfied with m = 7 in dimension d = 3, for ω = 0. Explicit expressions for the eigenprojectors π j (k) are also provided in Appendix B.
Let Assumptions 2.1, 2.3 and 3.1 be satisfied and assume moreover that u 0 = π 1 (k)u 0 + εr. Then (i) There exist T > 0 and, for all ε ∈ (0, 1], a unique solution u ∈ C([0, T /ε]; B (1) ) to (30) with initial condition u 0 . (ii) There exists ε 0 > 0 such that for all 0 < ε < ε 0 , the solution U to (20) provided by Theorem 3.1 exists on [0, T /ε] and Proof. Here again, the proof is a small variation of the one given in [9] for the trilinear case and B = W (R d ) n . We just indicate the main steps of the proof: Step 1. Existence and bounds of the solution u to (30) is established by a fixed point argument as in Theorem 3.1.
Step 3. Thanks to the assumption that ω = ω 1 (k) one gets from the equation obtained by applying Step 4. Using a non-stationary phase argument (on the semigroup formulation) relying on point (ii) of Assumption 3.1 and the bound on ∂ t u 1 established in Step 3, and taking advantage of Assumption 2.3, we get that 1 ε |u II (t)| B remains uniformly bounded on [0, T /ε].
Step 5. Using the non-resonance condition (iii) of Assumption 3.1, one can show that the third and higher harmonics created by the nonlinearity remain of order O(ε). More precisely, the solution U ∈ H 1 (T; B) to (25) provided by Theorem 3.1 can be written as Step 6. Since U (t) − U app (t) = εV (t), it follows from the above that sup and the theorem follows therefore from the observation that

16ÉRIC DUMAS, DAVID LANNES AND JÉRÉMIE SZEFTEL
3.3. The Full Dispersion model. The idea is to diagonalize (30) in order to work with a scalar equation. We project therefore (30) onto the eigenspace corresponding to the oscillating term. These naturally leads to introduce which naturally leads to the the full dispersion scalar 8 equation (writing u = u F D ) and with ω 1 (·) as in Assumption 3.1.
The following corollary shows that the full dispersion scalar equation yields an approximation of the same precision as the envelope equation for times t ∈ [0, T /ε]. to (31) with initial condition u 0 . (ii) There exists ε 0 > 0 such that for all 0 < ε < ε 0 , the solution U to (20) Remark 3.1. Equation (31) does not correspond exactly to the "full dispersion" model of [9,21], where the r.h.s is επ 1 (k)F env (ε, u) rather than επ 1 (k+εD)F env (ε, u) (but it can be found as an "intermediate model" in [8]). This change does not affect the estimate given in the Corollary, but it is important to keep track of the frequency dependence of the polarization of the nonlinear term to introduce the "self-steepening" operators in §3.6. Note also that the "full dispersion" model is related to the so-called "unidirectional pulse propagation equation" used in nonlinear optics [19,4].
Proof. This is actually a byproduct of the proof of Theorem 3.2 since u F D coincides with u 1 in Step 2 of the proof of Theorem 3.2.
3.4. The nonlinear Schrödinger (NLS) equation. The NLS equation is easily derived from (31) using Taylor expansions where c g = ∇ω 1 (k) and H k stands for the Hessian of ω 1 at k. Neglecting the O(ε 2 ) terms in (31) we define the N LS approximation u = u N LS as the solution to We then get easily (see [9,21]) the following justification of the NLS approximation.
Corollary 3.2 (Schrödinger approximation). Under the assumptions of Theorem 3.2, one has for all u 0 ∈ B (3) such that u 0 = π 1 (k)u 0 , 8 The linear propagator is a scalar operator but the equation remains vectorial because of the nonlinear term. Indeed, π 1 (k) is in general not of rank 1 (i.e. ω 1 (k) is in general not of multiplicity one) (i) There exists T > 0 and, for all ε ∈ (0, 1], a unique solution u ∈ C([0, T /ε]; B (3) ) to (32) with initial condition u 0 . (2) There exists ε 0 > 0 and c N LS > 0 such that for all 0 < ε < ε 0 , the solution U to (20) provided by Theorem 3.1 exists on [0, T /ε] and in the error estimate of the Corollary is due to the bad frequency behavior of the Schrödinger equation when the envelope of the oscillations ceases to be well localized in frequency. This is for instance the case for short pulses, chirped pulses ( [8,9,21]), and near the focusing point. To describe such extreme situations, the standard NLS approximation does a poor job, and this is why various variants have been derived in physics (e.g. [3]).
Remark 3.3. We assumed here that the polarization of the initial condition is exact (i.e. r = 0 in Theorem 3.2) for the sake of simplicity; indeed, the solution remains polarized along π 1 (k) for all times and computations on real physical models are much easier.
Example 3.3. In the frequent case where ω 1 (·) has a radial symmetry, and writing with a slight abuse of notation ω 1 (k) = ω 1 (k), with k = |k|, we can write, where we assumed without loss of generality that (0z) is the direction of the wave number k, k = ke z . In particular, (32) reads 3.5. The nonlinear Schrödinger equation with improved dispersion relation. We propose here to investigate further the Schrödinger equation with improved dispersion relation derived in [9]. As said in Remark 3.2, the NLS approximation has too bad dispersive properties to capture correctly the envelope of oscillating solutions to Maxwell's equations in extreme situations, where high frequencies are released. Indeed, the dispersion relation ω 1 (·) of (20) is approximated by the second order polynomial For Maxwell's equations and in dimension d = 1, Figure 1 shows that this dispersion relation is very poor when k ′ is not close to k. The idea introduced in [9] is to replace the linear part of the Schrödinger approximation by a linear operator that differs from the linear part of the Schrödinger approximation by O(ε 2 ) terms only, but whose dispersion relation is far better. More precisely, we consider an approximation under the form where u = u imp solves the nonlinear Schrödinger equation with improved dispersion relation where b ∈ C d , B ∈ M d×d (R) and C 3 (∇) is a third order homogeneous differential operator. We assume moreover that (36) B is symmetric positive, b ∈ Range(B), and 4 − b · (B −1 b) > 0 (note that even though B −1 b is not unique when B is not definite, the scalar b · (B −1 b) is uniquely defined). These assumptions ensure that the operator (1 − iεb · ∇ − ε 2 ∇ · B∇) is invertible. This new model allows one to replace (33) by .
A good choice of b, B and C 3 allows a much better approximation of ω 1 (·), as shown in Figure 1 for Maxwell's equation in dimension d = 1. Exactly as for Corollary 3.2 we get the following result, where the only difference in the error estimate with respect to Corollary 3.2 is the constant c imp (with is much smaller than c N LS for good choices of b, B and C). We refer to [9] for the proof and numerical validations of this model for the approximation of short pulses and chirped pulses. to (35) with initial condition u 0 .
(2) There exists ε 0 > 0 and c imp > 0 such that for all 0 < ε < ε 0 , the solution U to (20) provided by Theorem 3.1 exists on [0, T /ε] and A similar equation has been proposed in [22] §8.5.3 in the framework of water waves equations.
3.6. The NLS equation with frequency dependent polarization. The idea here is to improve the rough approximation π 1 (k + εD) ∼ π 1 (k) + O(ε) used to derive the NLS approximation (see §3.4). Indeed, when the description of the envelope of the laser pulse requires a broad band of frequencies as in the situations mentioned in §3.5, the variations of the polarization term π 1 (k + εD) in front of the nonlinearity in (31) become important and should be taken into account. We therefore make the following approximation, , where b and B are the same as in the NLS approximation with improved dispersion (35). In particular, if b = 0 and B = 0 (standard NLS), then the above approximation coincides with the first order Taylor expansion The general formula has the same accuracy as this Taylor expansion as ε → 0. The corresponding approximation is given by where u = u pol solves the nonlinear Schrödinger equation with frequency dependent polarization where b, B and C 3 (∇) are the same as in (31).
Contrary to all the previous models, the nonlinearity in (39) seems to be of quasilinear rather than nonlinear nature. It turns out however that the presence of the operator 1−iεb·∇−ε 2 ∇·B∇ in front of the time derivative plays a smoothing role allowing the control of one or several first order derivatives (see the discussion in point (4) in p. 3). If the first order derivatives involved in the nonlinearity are all controlled by this smoothing operator, then the nonlinearity remains semilinear in nature. As shown in the proof below, the component −iε(b·∇)π 1 (k)F env (ε, π 1 (k)u) of the nonlinearity is always semilinear in this sense. This is not the case for the component επ 1 (k)π ′ 1 (k) · DF env (ε, π 1 (k)u) that may be of quasilinear nature, in which case a symmetry assumption is needed on the nonlinearity to ensure local well-posedness. In order to state this assumption, it is convenient to introduce the norm | · | * defined as where e j is the unit vector in the j-th direction of R d and d v F env is the derivative at v of the mapping u → F env (ε, u). The approximation furnished by (39) is justified by the following corollary (the difference in the estimate with respect to Corollary 3.3 is just a better nonlinear constant, denoted C pol to insist on this point). For the sake of simplicity, we take B = H t0 (R d ) n (t 0 > d/2) here, but adaptation to Wiener spaces are possible. to (39) with initial condition u 0 . (2) There exists ε 0 > 0 and c imp > 0 such that for all 0 < ε < ε 0 , the solution U to (20) provided by Theorem 3.1 exists on [0, T /ε] and  (39)).
Remark 3.5. Note that in (39), we have applied π 1 (k) to the full nonlinearity (hence the term π 1 (k)π ′ 1 (k) · D instead of π ′ 1 (k) · D). This means that we keep the effect of the frequency dependent polarization on the main polarized space Range(π 1 (k)) only. Similarly, we have replaced F env (ε, u) by F env (ε, π 1 (k)u). This latter substitution would not change anything to the previous NLS models since we have seen that polarized initial conditions remain polarized. Its purpose in (39) is to make Assumption 3.2 much easier to check.
Proof. The justification of (39) is performed as for the other models. The only difference here is that local well-posedness for a time scale of order O(1/ε) must be established. We therefore show here that such a local well-posedness result holds if u 0 ∈ H s+1 (R d ) n with s > t 0 + 1. We just prove a priori energy estimates for (39); existence, uniqueness and stability with respect to perturbations can be deduced classically. The natural energy associated to (39) is given for all s ≥ 0 by Under the assumptions (36) on B and b, E s (u) 1/2 defines a norm that controls uniformly the H s -norm. It may also control a certain number of first order derivatives. The important point for the local well-posedness of (39) is that it always controls first order derivatives in the direction b · ∇; more precisely, we claim that there exists c > 0 independent of ε such that for all u,

A quick look on the Fourier side shows that it is equivalent to prove that
which is a consequence of (36). Multiplying Λ s (39) by Λ sū and integrating by parts, we get where we denoted F (u) = F env (ε, π 1 (k)u). It is straightforward to control I 1 , and Moser's estimate gives a control of I 2 , I 1 ε 1+p |u| 2 H s , I 2 ≤ εC(|u| ∞ )|u| 2 H s . In order to control I 3 , we must split it into two parts, the first component is controlled using Assumption 3.2 while the Kato-Ponce and Moser estimate give a control of the second one, Remarking that Λ s (b · ∇)F (u) = Λ s d u F (b · ∇)u and using the tame product estimate |f g| H s |f | L ∞ |g| H s + |f | H s |g| L ∞ and Moser's inequality, we get Gathering all the above estimates and using (40), we obtain Since moreover |u| 2 W 1,∞ E s (u) for s > d/2 + 1, we deduce from Gronwall's estimate that for such s, the energy E s (u) remains bounded for times of order O(1/ε).
Again, U (t, x, θ) and W (t, x, θ) are periodic with respect to θ, and we use any ω ∈ R. The action of the Fourier multiplier H ε k k 2 · D from (19) (and (16)) is transferred at the profile level into the operator H D θ + ε k k 2 · D , with k = |k|, where H D θ + ε k k 2 · D is the Fourier multiplier which acts continuously on H k (T, B), for any k ∈ N. In order to get a solution to the original problem, it is sufficient that (U, W ) solves where Z = (U, W ) solves the profile equation (43), with the initial condition (44).
Proof. Similar to the one of Theorem 3.1, by an iterative scheme in H k (T; B), with k ≥ 1.

3.7.2.
The slowly varying envelope approximation. In this section, as in in section 3.2, we shall assume that ω = ω 1 (k) is some characteristic frequency for L(·, k).
Postulating the Ansatz we obtain formally the following system for (u env , w env ) (denoted (u, w)), Here, F env is given by (29) and G env is defined in the same way, filtering oscillations, Remark 3.6. The equation one obtains for w from direct computations is actually in view of Assumption 2.2.
The approximation (45) is justified in the following theorem.
Step 1. As in the proof of Theorem 3.1 and Theorem 3.3, we have local in time T /ε (with T independent of ε) existence of (u, w), solution to (46), together with (uniform w.r.t. ε) bounds.
Step 2. Decomposing u = u 1 + u II as in Step 2 of the proof of Theorem 3.2, one obtains in the same way that |∂ t u 1 | L ∞ ([0,T /ε],A n ) is bounded w.r.t. ε. This is also the case for |∂ t w| L ∞ ([0,T /ε],A) (even of order O(ε)), as shows directly the third equation in (46).
Step 3. As in Step 4 of the proof of Theorem 3.2, we deduce that This is obtained by integration by parts in the integral formulation giving the u ′ j s of u II , using Step 2 and Assumption 3.1 to have non-stationary phase; we conclude by a Gronwall argument.
Step 4. Approximation of Z by Z SV EA . Compared to Theorem 3.2, the new point is the component w; for the sake of simplicity, we therefore set F = 0 throughout this proof. Denote where Z = (U, W ) is the solution to the profile equation (43) provided by Theorem 3.3. We estimate Z = ( U , W ) in H k (T; B) (defined in (27)), k ≥ 1.
• For W . We have W |t=0 = 0 and The terms in c 1 and c 2 can be treated similarly. For the sake of clarity, we therefore set c 1 = 0 and c 2 = 1 in this proof, so that the right-hand side is given by we easily deduce that Splitting u = u 1 + u II as in Theorem 3.2, we have a uniform (in ε) bound for Step 1 and Step 3. The only component left to control is therefore the one involving the product (C 1 u 1 ) · (C 1 u 1 ), for which we write Using the equation satisfied by w to control ∂ t w and Step 2 to control ∂ t u 1 , one readily deduces that • For U . We have U |t=0 = 0 and where, for the sake of simplicity, we also have taken c = 0 because the corresponding terms do not raise any difficulty. Since the Fourier multiplier H D θ + ε k k 2 · D − i has a norm of order ε when acting on A and because the difference of the last two terms in easily bounded in terms of εZ, we get that | U (t)| H k (T;A) satisfies the same upper bound as | W (t)| H k (T;A) in (47).
Gathering the upper bounds for W and U we therefore get By a Gronwall estimate, we finally obtain a bound on Z, Using Theorem 3.4, this model can be justified with a straightforward adaptation of the case without ionization.
3.7.4. The most general model. In (48), ionization effects have been added to the standard cubic NLS equation. It is of course possible to add them to a most sophisticated model that takes into account more general nonlinearities, improved frequency dispersion and frequency dependent polarization. We then obtain the following generalization of (48) (for which the same justification as in Corollary 3.5 holds), As shown in Appendix B.2, this system of equations takes the following form corresponding to (4) in the case of Maxwell's equations,

Analysis of the equations (3) and (4), and open problems
In this section, we analyze and formulate open problems for the NLS-type equations (3), as well as the NLS-type equations (4) taking the ionization processes into account. We first consider in section 4.1 the equations (3) in the case of respectively no and anomalous GVD (resp. α 1 = 0 and α 1 = 1), and also briefly discuss (4). Then, we discuss (3)  As recalled in the introduction, (1) is critical in dimension 2, and supercritical in dimension 3, and there exist finite time blow-up solutions in both cases. Let us also recall that some of these blow-up dynamics are stable (see for example [26] [27] [28] [29] [30] in the critical case, and [32] in a slightly supercritical case). Now, let us recall that the breakdown of solutions is not always observed in physical experiments on the propagation of laser beams and that instead a filamentation phenomenon occurs. In this section, we would like to analyze the possibility that the modified model (3) in the case of no or anomalous GVD prevents the formation of singularities. Below, we analyze the role of each parameter of (3) separately starting with the nonlinearity f . 4.1.1. The nonlinearity. In the case where P 2 = I, α 2 = 0, and α 3 = 0, (3) takes the following form, where the dimension is d = 2 in the case of no GVD, and d = 3 in the case of anomalous GVD.
As recalled in the introduction, standard modifications of the cubic nonlinearity consist either of the cubic/quintic nonlinearity, i.e. f (s) = −s, or a saturated nonlinearity, i.e. f is a smooth function on R + vanishing at the origin and such that (1 + f (s))s is bounded on R + (e.g. f (s) = − s 1+s ). Let us first consider the case of a saturated nonlinearity. In that case, the fact that (1 + f (s))s is bounded on R + implies the following control of the nonlinear term in (51): Thus, running a fixed point argument yields the fact that this equation is locally well posed in L 2 (R d ) for any integer d ≥ 1, with a time of existence only depending on the size of u 0 L 2 . Since equation (51) still satisfies the conservation of mass, this immediately implies global existence for any initial data in L 2 (R d ) and for any d ≥ 1. Therefore, modifying the cubic nonlinearity in the standard NLS equation by a saturated nonlinearity does indeed prevent finite-time breakdown of the solutions.
Next, we consider the modification by a cubic/quintic nonlinearity. In that case, (51) becomes, Note that (52) still satisfies the conservation of mass, and the conserved energy is now given by, Note also that, which together with the conservation of mass and energy yields, Let us first consider the case of dimension 2. In this case, both the cubic focusing NLS and the quintic defocusing NLS are H 1 subcritical in the sense that the equation is locally well posed in H 1 (R 2 ) with a time of existence only depending on the size of u 0 H 1 (see [15]). The proof extends to (52) which is thus locally well posed in H 1 (R 2 ) with a time of existence only depending on the size of u 0 H 1 . Together with the bound (53), this immediately implies global existence for any initial data in H 1 (R 2 ). The case of dimension 3 is more involved, since the quintic defocusing NLS is H 1 critical in the sense that the equation is locally well posed in H 1 (R 3 ) with a time of existence depending on the shape of u 0 (see [6]). Thus, one cannot rely solely on the bound (53) to prove global existence. However, the case of combined nonlinearities is addressed in [41], where global existence of solutions with H 1 initial data is proved, when the nonlinearity is a sum of two H 1 -subcritical powers, or the sum of an H 1 -subcritical power and of a defocusing H 1 -critical power. This includes equation (52). Therefore, modifying the cubic nonlinearity in the standard NLS equation by a quintic defocusing nonlinearity does indeed prevent finite-time breakdown of the solutions both in dimension 2 and 3. Let us mention an interesting phenomenon regarding the qualitative behavior of the solutions. Physical experiments suggest an oscillatory behavior of the solution which focuses, then defocuses, refocuses,... in an almost periodic fashion (see [3] and references therein). Such a behavior is also observed in numerical simulations and suggested by heuristic arguments (see for example [38] and [13]).
Another open problem concerns the behavior of the solution as ε → 0. Pick an initial data v 0 ∈ H 1 (R d ), d = 2, 3, leading to a finite time blow-up solution v at time T > 0 to the cubic NLS (1). Consider the solution v ε of (51) with the same initial data v 0 . In the case of dimension 2, it is shown in [24] that An interesting question is the understanding of the limit of v ε (t) as ε → 0 for t > T . Partial results have been obtained in this direction in [25]. There, it is proved that only few scenarios are possible, but one would like to establish whether all scenarios do occur or only some of them, and which scenarios are generic.
Open Problem 2. Let v 0 ∈ H 1 (R d ), d = 2, 3 be leading to a finite time blow-up solution v at time T > 0 to the cubic NLS (1). Consider the solution v ε of (51) with the same initial data v 0 . Describe the behavior of v ε (t) as ε → 0 for t > T .

4.1.2.
Taking the ionization process into account. Let us discuss equations (4) and (5) in the case of no or anomalous GVD which take the ionization process into account and also correspond to a modification of the nonlinearity in the focusing cubic NLS (1). Ionization is certainly the most important phenomenon leading to the formation of laser filaments. Recall that (4) is given by where α 4 , α 5 ≥ 0, c > 0, where ρ is the density of electron created by ionization, while c g = c g e z is the group velocity associated to the laser pulse, and where d = 2 in the case of no GVD, and d = 3 in the case of anomalous GVD. Also, recall that (5) is given by whereρ corresponds to ρ written in a frame moving at the group velocity c g = c g e z and with respect to a rescaled time τ = εt, i.e. ρ(t, X ⊥ , z) =ρ(εt, X ⊥ , z − c g t).
In view of the second equation of (55), a boundary condition forρ has to be prescribed at z = z 0 for some z 0 in order to obtain a well-posed problem. A natural choice, which ensures thatρ ≥ 0, is to prescribeρ at +∞: Let us first discuss the local well-posedness theory for equations (54) and (55), starting with the first one. One may obtain the existence of solutions to (54) over an interval of time with size O(ε −1 ). Indeed, for an integer N > d/2, differentiating N times with respect to space variables the equations both for v and ρ, multiplying (54) respectively by ∇ N v and ∇ N ρ, integrating on R d , taking the imaginary part and integrating by parts for the first equation, and summing both equations yields: where we used the fact that N > d/2 together with simple product rules and the Sobolev embedding. Integrating this differential inequality, we obtain a time of existence with size O(ε −1 ).
Next, we discuss the local well-posedness theory for equation (55) supplemented with the boundary condition (56). For any integer N , we introduce the space H N , For an integer N > d, differentiating N times with respect to space variables the equations for v, multiplying the first equation of (55) by ∇ N v, integrating on R d , taking the imaginary part and integrating by parts yields: where we used the fact that N > d together with the Sobolev embedding in the last inequality. Next, we estimate ρ H N . In view of the second equation of (55) and the boundary condition (56), we havẽ ρ(t, x, y, z) = +∞ z α 4 ε|v| 2K (t, x, y, σ) exp α 5 ε c σ z |v| 2 (t, x, y, s)ds dσ.
In view of the boundary condition (56), this yields x, y, s)ds dσ.
Hence, we deduce . By induction, we obtain Integrating this differential inequality, we obtain local existence for (v,ρ) in H N × H N for N > d.
Remark 4.1. Note that the factor ε in the second equation of (55), while present in the model, is not needed for the well-posedness theory.
Let us now come back to the issue of global well-posedness/finite time singularity formation. For equations (54) and (55), note that the L 2 norm of u (resp. v) is dissipated. Indeed, in the case of (55), multiplying by v, integrating on R d , taking the imaginary part and integrating by parts yields: from which we deduce that the L 2 norm is dissipated sinceρ ≥ 0. On the other hand, the energy is not conserved, nor dissipated. We thus can not carry out the analysis of section 4.1.1 for equation (52), even in the case of dimension 2 9 . An interesting problem would then be to prove that taking the ionization process into account (i.e. replacing the cubic focusing NLS (1) by either equation (54) or equation (55)) does indeed prevent finite-time breakdown of the solutions in dimensions 2 and 3. This is formulated in the following open problem.
Open Problem 3. Prove that the solutions to equations (54) and (55) are global in dimensions 2 and 3.

The damping.
In the case where P 2 = I, α 3 = 0 and f = 0, (3) takes the following form, where d = 2 in the case of no GVD, and d = 3 in the case of anomalous GVD. The mass and energy are not conserved quantities anymore. The mass decreases, Note that a simpler model of damped NLS where the ρ orρ term are not present has been investigated in [2]. The authors obtain global existence for a certain range of parameters by controlling a modified energy even if it is not conserved.
while for the energy we have, so that the energy is neither increasing nor decreasing. Equation (57) has been analyzed in several works (see e.g. [40], [38], [12] and [33]). In particular, from standard arguments using Strichartz estimates, there is in dimensions 2 and 3 a continuous and increasing function θ such that global existence holds for α 2 > θ( u 0 H 1 ) (see [38] and [33]). Also, in the 3-dimensional case, there is a function θ such that stable blow-up solutions exist for 0 ≤ α 2 ≤ θ(u 0 ) (see [40] and [33]). Finally, the same behavior has been observed in numerical simulations and suggested by heuristic arguments in the 2-dimensional case (see [12]). Thus, a small α 2 does not prevent finite time blow up. Since the constant α 2 obtained in (4) is usually small, it appears that the damping cannot by itself explain the physical observation according to which the breakdown of solutions does not occur.

4.1.4.
Off-axis variation of the group velocity. In the case where α 2 = 0, α 3 = 0 and f = 0, (3) takes the following form, where d = 2 in the case of no GVD, and d = 3 in the case of anomalous GVD. The energy remains unchanged and is still conserved, while the mass is replaced by the following conserved quantity, Let us first consider the case of full off-axis dependence, i.e. the case where, The operator P 2 (ε∇) is a second order self-adjoint operator which is invertible in view of (62). We denote by P 2 (ε∇) −1 its inverse. We rewrite (60) in the following form, Using Duhamel's formula, and the semi-group e itP2(ε∇) −1 ∆ , we obtain, As a simple consequence of (62) and the Sobolev embedding in dimensions 2 and 3, we have, 3. Now, a fixed point argument based on the formulation (64) together with the estimate (65) and the fact that the semi-group e itP2(ε∇) −1 ∆ is unitary on H 1 , implies that (60) is locally well posed in H 1 for dimensions d = 2, 3, with a time of existence only depending on the size of u 0 H 1 . Together with the lower bound (62) and the conserved quantity (61) which yield a uniform in time bound on the H 1 norm of the solution, this immediately implies global existence for any initial data in H 1 when d = 2, 3. Therefore, modifying the focusing cubic NLS equation by adding a full off-axis dependence does indeed prevent the finite-time breakdown of the solutions.
Next, let us consider the case of partial off-axis dependence, i.e. the case where, for vectors v k in R d , k = 1, . . . , j, where j = 1 if d = 2, and j = 1 or j = 2 if d = 3 so that the right-hand side of (66) does not control the full H 1 norm. Notice that the well-posedness theory has yet to be investigated in this case. Indeed, consider the nonlinear term in the formulation (64), In the directions v k of (66), the operator P 2 (ε∇) −1 gains two derivatives, while the semi-group e itP2 ( [1] for more details on optical shocks). The reason we included this modification in the discussion is because it may account for the physical observation of the self-steepening of the pulse: an initial pulse which is symmetric becomes asymmetric after propagating over a large distance and its profile seems to develop a shock. This phenomenon has been widely observed and we refer the reader to [3] and references therein. Thus, an interesting open problem would be to exhibit solutions of (68) for which the corresponding profile develops a shock. We formulate a slightly more general open problem below.
Open Problem 5. Describe the blow-up scenarios for the finite time blow-up solutions to equation (68).
Equation (69) admits the following conservation laws in the energy space H 1 , Note that (69) is not invariant under the usual Galilean transform. However, it is invariant under a twisted Galilean transform. For β = (β 1 , β 2 , β 3 ) ∈ R 3 , we define, Then (69) is invariant under the following twisted Galilean transform, . Note also that the scaling symmetry u(t, x) = λ 0 v(λ 2 0 t, λ 0 x) leaves the space L 2 (R 2 ) invariant so that (69) is critical with respect to the conservation of mass in dimension 2, while it leaves the homogeneous Sobolev spaceḢ 1 2 (R 3 ) invariant so that (69) is supercritical with respect to the conservation of mass in dimension 3.
In contrast to the cubic focusing NLS, the existence or absence of finite time blow up solutions for (69) is widely open. While there exists a counterpart to the virial for (69), it is too weak to provide the existence of finite-time blow up dynamics (see the discussion in [38]). Also, there are no standing waves in the form v(t, x) = Q(x)e iωt , ω ∈ R in the energy space for the equation (69) (see [14]), while for the focusing cubic NLS, such an object exists and is fundamental in the analysis of the blow up dynamics (see e.g. [26] [27] [28] [29] [30] [31] in dimension 2, and [32] in dimension 3). Now, numerical simulations suggest that the solutions to (69) do not break down in finite time (see e.g. [38], [13] and references therein). In particular, the simulations exhibit a phenomenon called pulse splitting where the pulse focuses until it reaches a certain threshold. Once this threshold is attained, the pulse splits in two pulses of less amplitude moving away from each other. This phenomenon might repeat itself (multiple splitting) and ultimately prevent finite time blow up, but this has not been clearly backed up by numerics so far. Thus an interesting problem would be to prove that the cubic NLS equation (69) (i.e. in the case of normal GVD) does not have finite time blow up solutions, whether this is due to pulse splitting, or to some other phenomenon. This is formulated in the two open problems below.
Open Problem 6. Prove that the solutions to equation (69) are global in dimension 2 and 3.
Open Problem 7. Describe rigorously the phenomenon of pulse splitting for equation (69).
Remark 4.2. We only considered in this section the case where, P 2 = I, α 2 = 0, α 3 = 0 and f = 0 which is (69). While it is of interest to investigate the qualitative role of the parameters P 2 , α 2 , α 3 and f of equation (3) in the case of normal GVD (i.e. the case α 1 = −1), we chose not to investigate their effects since we are interested in this paper primary on the non existence of focusing dynamics, and equation (69) is expected to have only global solutions. Remark 4.3. In a recent preprint [18], standing waves and selfsimilar solutions have been exhibited for equation (69) in dimension 2. However, these solutions do not belong to the energy space. In particular, they are not in C 2 and do not decay along the two diagonals in R 2 .

4.3.
Mixing several phenomena. In section 4.1 and section 4.2, we have investigated the effect of each of the parameters P 2 , α 1 , α 2 , α 3 and f on the solutions to equation (3), and in particular whether these modifications of the cubic NLS (1) prevent the existence of finite time blow-up solutions. We have also formulated several open problems. Now, instead of studying all these effects separately, one may investigate the case when all these phenomenon are present at the same time, and ask which are the dominant ones? Consider for instance the case where there is partial off-axis variation of the group velocity in the direction z, and at the same time a potential self-steepening of the pulse in the same direction, i.e. P 2 = 1 − ∂ 2 z and α 3 = e z , Then, one may wonder whether the partial off-axis variation of the group velocity prevents the self-steepening from taking place. This is formulated in the open problem below.
where v is now a C 2 -valued function, and where we consider for simplicity the equation corresponding to the cubic case (i.e. f = 0 in (3)). In fact, (3) is a particular case of (3) vect corresponding to initial data living on a one-dimensional subspace of C 2 (see Remark 1.2). Now, one may of course consider the previous questions formulated for equation (3), and investigate the same problems for the vectorial case. This is formulated in the open problem below.
Open Problem 9. Investigate the vectorial counterparts for equation (3) Such a justification is far from being sharp. In particular, for the standard cubic focusing NLS equation, one has T < T foc , where T foc is the time when focusing occurs. Indeed, the above justification process requires that u app remains bounded on [0, T /ε] (after rescaling, this is equivalent to require that the solution v to (3) is bounded on [0, T ]). In the error estimate above, the constant C(T, |u app | L ∞ ([0,T /ε];B) ) therefore blows up as T → T foc . It follows that all the NLS-type equations derived here are justified far enough from the focusing time. Now, relevant differences between the different models can only be observed close enough to focusing (they all differ by formally O(ε 2 ) terms that become relevant only near the focusing point). So, roughly speaking, we have only proved that all the models are justified on a time interval for which they are basically identical.
As seen above, existence beyond the focusing time of the standard cubic NLS equation is proved or expected for many of the variants considered here. For such models, the above argument does not work, i.e., the constant C(T, |u app | L ∞ ([0,T /ε];B) ) of the error estimate does not blow up as T → T foc . This is for instance the case of the cubic/quintic NLS equation (52) that admits a global solution v ε . However (see Open Problem 2), we have so that there is no reason to expect the error term εC(T, |v ε | L ∞ ([0,T /ε];B) ) to be small near the focusing point and for small ε.
For the moment, the merits of the NLS-type models derived here to describe correctly the mechanisms at stake during focusing can only by assessed numerically. Hence the following interesting open problem, Open Problem 10. Rigorously justify one of the NLS-type models derived here on a time scale [0, T /ε] with T ≥ T foc . 10 In the error estimates of §3, the constant on the r.h.s. is C(T, |u 0 | B ) rather than C(T, |uapp| L ∞ ([0,T /ε];B) ). Since |uapp| L ∞ ([0,T /ε];B) = C(T, |u 0 | B ), this is of course equivalent, but the first form is more convenient for the present discussion.
Acknowledgements. The autors warmly thank Christof Sparber for pointing reference [41] out.
Appendix A. Nondimensionalization of the equations A.1. The case without charge nor current density. There are two characteristic times for the situation considered here. The first one, denotedt, is the inverse of the frequency of the laser pulse; the second one, denotedT is the duration of the pulse. In the regimes considered here 11 , one hasT ≫t, and the small parameter ε is defined as The time variable is naturally nondimensionalized withT , while the space variables are nondimensionalized by L = cT , namely, where dimensionless quantities are denoted with a tilde. We also nondimensionalize the unknowns E, B and P with typical orders E 0 , B 0 and P 0 . The typical scale for the electric and magnetic fields are directly deduced from P 0 , In order to choose the remaining P 0 , some consideration on the polarization equation is helpful. We recall that equation (11) is given by (iii) Saturated nonlinearity: All these examples can therefore be put under the form ∇PṼ N L (P) = 1 +f (ãε r |P| 2 ) |P| 2P , wheref : R + → R is a smooth function such thatf (0) = 0.
In dimensionless form, the polarization equation becomes therefore We therefore choose P 0 to be the typical nonlinear scale for which the nonlinear terms in the above equations are of the same order as the second order time derivative, namely, Finally, we need some information on the size of Ω 0 , Ω 1 and b to give our final dimensionless form of the equations. The resonance frequency of the harmonic oscillator is typically of the same order as the frequency of the laser pulse, the damping frequency is written under the form where the dimensionless damping frequency ω 1 has size O(1) while the coefficient α depends on the medium in which the laser propagates. Typically, α = 1 or α = 2. Writing the coupling constant b as the dimensionless polarization equations can therefore be written as Omitting the tildes for the sake of clarity, Maxwell's equations (10)-(11) read therefore (without charge nor current density), Introducing the auxiliary unknowns, and working with nonlinearities as those considered in Example A.1, this system can be written as a first order system, where we wrote f (x) =f (ã γ This system has the form (20) with n = 12, U = (B, E, Q ♯ , P ♯ ) T and while A 0 is the block diagonal matrix A 0 = diag(0, 0, ω 1 I, 0), so that the first two points of Assumption 2.1 are satisfied. For the third one, remark that the nonlinearity is given by A.2. The case with charge and current density. Now, we deal with the nondimensionalization of equations (10)-(11)- (17). We proceed in the same way as in the previous paragraph, simply adding the characteristic sizes J 0 and ρ 0 of the free electron current density J e and charge density ρ. Considerations on the polarization equation are the same. We choose ρ 0 from the equation for the charge density, Many configurations are possible, according to the numerical value of the various coefficients involved (see for instance [3] for experimental data). We choose here a configuration that leads to the richest model, where the coupling constants of the two nonlinear terms are of the same order; more precisely ∂tρ = εc 1 |Ẽ| 2K + εc 2ρ |Ẽ| 2 , for non-negative constants c 1 , c 2 . This corresponds to the following choice for ρ 0 , Knowing ρ 0 , we can then determine J 0 from the equation for the free electron current density, J e = q 2 e ω l m e ρ 0 E 0ρ H(εDz)Ẽ, which naturally leads to J 0 = q 2 e ω l m e ρ 0 E 0 = ε q 2 eT ρ 0 E 0 m e .
Going back to Maxwell's equations, we get    ∂tB + curlẼ = 0, Here again, many configurations can be found, and we choose one that leads to the richest models, namely, where c 0 , c 3 are dimensionless positive constants. ReplacingJ e by c 3Je ,ρ by c 3ρ and c 1 by c 1 c 3 , we can assume that c 3 = 1. Finally, dropping the tildes, introducing the unknowns Q ♯ and P ♯ , and using the same notations as in (73), we obtain Appendix B. Explicit computations for Maxwell equations B.1. The case without charge nor current density. We derive here the variants of the NLS equations derived in this paper in the particular case of the Maxwell equations (73). We derive all the versions of the NLS equations that do not take into account the frequency dependence of the polarization (this corresponds to (37)). This latter effect, which leads to the generalization (41) of the previous model, is examined separately. For the sake of simplicity, we only consider cubic nonlinearities here; we show that in this context, (41)  and ω ∈ R, u = (b, e, q ♯ , p ♯ ) ∈ C 12 . From the last two equations of (73), we can always write p ♯ and q ♯ in terms of e, Let us now consider several cases: • If e k then by the first equation, either ω = 0 and b k or b = 0 and ω 2 = ω 2 0 + γ.