On Critical Behaviour in Systems of Hamiltonian Partial Differential Equations

We study the critical behaviour of solutions to weakly dispersive Hamiltonian systems considered as perturbations of elliptic and hyperbolic systems of hydrodynamic type with two components. We argue that near the critical point of gradient catastrophe of the dispersionless system, the solutions to a suitable initial value problem for the perturbed equations are approximately described by particular solutions to the Painlevé-I (PI\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_I$$\end{document}) equation or its fourth-order analogue PI2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_I^2$$\end{document}. As concrete examples, we discuss nonlinear Schrödinger equations in the semiclassical limit. A numerical study of these cases provides strong evidence in support of the conjecture.

we discuss nonlinear Schrödinger equations in the semiclassical limit. A numerical study of these cases provides strong evidence in support of the conjecture.

Introduction
Critical phenomena in the solutions of partial differential equations (PDEs) are important from various theoretical and applied points of view since such phenomena generally indicate the appearance of new behaviours as the onset of rapid oscillations, the appearance of multiple scales, or a loss of regularity in the solutions. Some of the most powerful techniques in the asymptotic description of such phenomena are due to the theory of completely integrable systems which were so far restricted to integrable PDEs. In Dubrovin (2006), this restriction was overcome by introducing the concept of approximate integrability up to a finite order of some small parameter . This has allowed to apply techniques from the theory of integrable systems to a large class of non-integrable equations and to obtain asymptotic descriptions of solutions to such equations in the vicinity of critical points of these PDEs. The scalar case was studied along these lines in Dubrovin (2006). Basically it was shown that solutions to dispersive regularizations of a nonlinear transport equation near a point of gradient catastrophe (for the transport equation itself) behave like solutions of the celebrated Korteweg-de Vries equations, which at such point can be asymptotically expressed in terms of a particular solution to a fourth-order ordinary differential equation from the Painlevé-I family. In Dubrovin et al. (2009), this concept was generalized to the study of the semiclassical limit of the integrable focusing cubic nonlinear Schrödinger equation (NLS) which can be seen as a perturbation of 2 × 2 elliptic system and in Dubrovin (2008) to a certain class of integrable Hamiltonian perturbation of 2 × 2 elliptic and hyperbolic systems. The idea that integrable behaviour persists in certain non-integrable cases has been already developed in the study of long-time behaviour of solutions to several non-integrable equations, like the perturbed NLS equation (Deift and Zhou 2002), see also Tao (2009) for a general overview about the soliton resolution conjecture.
The persistence of integrability for a rather general class of infinite-dimensional Hamiltonian systems and in particular for perturbed integrable equations has been also considered in the framework of the KAM theory. For example, as it was shown in Kuksin (1988), under suitable regularity conditions on the perturbation to the KdV equation, there exists a Cantor set of invariant tori supporting linearly stable solutions periodic in space and quasi-periodic in time. A similar result has been proven for NLS-type equations in Kuksin and Poeschel (1996). However, the phenomena studied in the present paper seem to be of a different nature as they describe local asymptotics near the point where the trajectory of the infinite-dimensional Hamiltonian system switches from a family of zero-dimensional to one-dimensional invariant tori.
In this paper, we consider general two-component Hamiltonian systems which contain a small dispersion parameter . When = 0, the Hamiltonian system reduces to a 2 × 2 quasilinear system of elliptic or hyperbolic type so that the Hamiltonian system can be considered as a perturbation of the elliptic or hyperbolic systems. We study the behaviour of solutions to such Hamiltonian systems when the parameter tends to zero. The fundamental question we address is how does a solution to Hamiltonian equations behave near the point where the solution of the unperturbed elliptic or hyperbolic system breaks up.
We consider Hamiltonian PDEs obtained as perturbations of systems of hydrodynamic type of the form (1.1) is a smooth function of u and v. Such perturbations can be written in the form .
( 1.2) where H is the perturbed Hamiltonian, H = H 0 + H 1 + 2 H 2 + · · · . By definition, the kth-order term of the perturbative expansion must have the form . . . , u (k) , v (k) . . . , u (k) , v (k) is a graded homogeneous polynomial of degree k in the variables u x , v x , …, u (k) , v (k) ; i.e. it satisfies the identity . . . , u (k) , v (k) for an arbitrary λ. The Hamiltonian system (1.2) can be considered as a weakly dispersive perturbation of the first-order quasilinear system (1.1). After certain simplification of the system (1.2) by a suitable class of -dependent canonical transformations generated by a Hamiltonian F (see Sect. 2) the system can be spelled out as follows up to terms of order 3 . Here a = a (u, v), b = b(u, v), and c = c (u, v) (1.4) The family of equations of the form (1.3) contains important examples such as the generalized nonlinear Schrödinger (NLS) equations (also in a non-local version), the long-wave limit of lattice equations like the Fermi-Pasta-Ulam or Toda lattice equation, Boussinesq equation, two-component Camassa-Holm equation (Falqui 2006), and many others. For certain choices of the functions h (u, v), a(u, v), b(u, v), and c(u, v), the system of Eq. (1.3) is integrable up to the order 3 (Dubrovin 2008). However, the complete classification of integrable cases in the class of equations of the form (1.3) remains open; see Degasperis (2009), Dubrovin et al. (2006), Kodama and Mikhailov (1997) for the current state of the art in this context.
The study of scalar weakly dispersive equations of the form similar to (1.2), (1.4) in the limit → 0 in the strongly nonlinear regime was initiated by the seminal paper by Gurevich and Pitaevskii (1973) about "collisionless shock waves" described by KdV equation (see also the book Novikov et al. 1984 and references therein). Rigorous mathematical results in this direction were obtained by Lax and Levermore (1983), Lax et al. (1993), Venakides (1990), and Deift et al. (1997) (see also Klein 2007, 2012 for numerical comparison). For twocomponent systems (1.3), an analogous line of research was started with the works Carles (2008), Gérard (1993), Grenier (1998) on the semiclassical limit of generalized defocusing nonlinear Schrödinger equation in several space dimensions for times less than the critical time t 0 of the cusp catastrophe. It was studied in more details for arbitrary times for the integrable case (Zakharov and Shabat 1972), namely for the spatially one-dimensional cubic defocusing NLS in Jin et al. (1994Jin et al. ( , 1999, DiFranco and Miller (2008). Another system that is included in the class (1.2) is the long-wave limit of the Toda lattice equation that has been studied in detail for arbitrary times in Deift and McLaughlin (1998), and in the context of Hermitian random matrix models with exponential weights by many authors, see the book (Deift 1999) and references therein. Interesting results, in the spirit of the original Gurevich and Pitaevsky setting, have been obtained for certain non-integrable cases in El (2005), Hoefer and Ilan (2012). Possible relations between integrable and non-integrable behaviour have been also analysed in the framework of the long-wave limit of the Fermi-Pasta-Ulam system by Zabusky and Kruskal (1965) and, more recently, in Bambusi and Ponno (2008), Lorenzoni and Paleari (2006), Benettin and Ponno (2011). The study of solutions to Hamiltonian systems of the form (1.3) in the limit → 0 with the leading term (1.1) of elliptic type was initiated by the analysis of the semiclassical limit of the focusing cubic nonlinear Schrödinger equation (Kamvissis et al. 2003;Tovbis et al. 2004); see also Bronski and Kutz (2002), Ceniceros and Tian (2002), Lyng and Miller (2007), Tovbis et al. (2006). Other interesting Hamiltonian systems not included in the class (1.3) have been considered in the limit → 0 in Miller and Xu (2012), Buckingham and Miller (2012).
Our study can be considered as a continuation of the programme initiated in Dubrovin (2006Dubrovin ( , 2008 and Dubrovin et al. (2009) aimed at studying critical behaviour of Hamiltonian perturbations of quasilinear hyperbolic and elliptic PDEs. The most important of the concepts developed in these papers is the idea of universality of the critical behaviour. We borrow this notion from the theory of random matrices where various universality types of critical behaviour appear in the study of phase transitions in random matrix ensembles; see, for example, Bleher and Its (1999), Bertola and Tovbis (2011), Deift et al. (1999a, b), Duits and Kuijlaars (2006),  for mathematically oriented references. The description of the critical behaviour for generalized Burgers equation with small viscosity was found by Il'in (1992); for more general weakly dissipative equations, see Dubrovin and Elaeva (2012), Arsie et al. (2013).
In the present paper solutions, u(x, t; ), v(x, t; ) to the Cauchy problem u(x, 0; ) = u 0 (x), v(x, 0; ) = v 0 (x) , (1.5) for the system (1.3) with -independent smooth initial data in a suitable functional class will be under consideration. The contribution of higher-order terms is believed to be negligible as long as the solution (u(x, t; ), v(x, t; )) remains a slowly varying function of x and t; that is, it changes by O(1) on the space-and timescale of order O −1 . A rigorous proof of such a statement would justify existence, for sufficiently small values of the parameter , of the solution to the Cauchy problem (1.2), (1.5) on a finite time interval 0 ≤ t ≤ t 0 depending on the initial condition but not on . This was proven by Lax and Levermore (1983), Lax et al. (1993) for the particular case of the Korteweg-de Vries (KdV) equation with rapidly decreasing initial data. In a more general setting of a certain class of generalized KdV equations with no integrability assumption, the statement was proven more recently in Masoero and Raimondo (2013). Actually, we expect validity of a more bold statement that, in particular, gives an efficient upper bound for the lifespan of a solution to (1.2) with given initial data (1.5). Namely, we start with considering the solution (u(x, t), v(x, t)) to the Cauchy problem for the unperturbed system (1.1) with the same -independent smooth initial Such a solution exists for times below the time t 0 of gradient catastrophe. 2 We expect that the lifespan of the perturbed solution (u(x, t; ), v(x, t; )) for sufficiently small is at least the interval [0, t 0 ]. More precisely, we have the following Main Conjecture consisting of three parts.
In the Main Conjecture, we do not specify the class of boundary conditions for the smooth (or even analytic, in the elliptic case) initial data u 0 (x), v 0 (x). We believe that the statement is applicable to a wide class of boundary conditions like rapidly decreasing, step-like, periodic. Moreover, the shape of the universal critical behaviour at the point of catastrophe should be independent of the choice of boundary conditions. The Cauchy problem for the elliptic system is ill-posed for non-analytic initial data, while the Cauchy problem for the corresponding dispersive regularization is generically well posed, at least locally. This is the case, for example, for the semiclassical limit of the focusing NLS equation with initial data with compact support and with discontinuities. The behaviour of the solution in the semiclassical limit has been stud-1 Analyticity of the initial data will be assumed in case the quasilinear system (1.1) is of elliptic type. The precise formulation of our Main Conjecture has to be refined in the non-analytic case. 2 As was shown in Kodama and Mikhailov (1997), existence of a point of gradient catastrophe for a secondorder system of quasilinear PDEs is a very general phenomenon, provided that the system is not linearly degenerate (see also the book Majda 1984). Here we choose the first catastrophe point; due to our genericity assumptions, it is assumed to be isolated and satisfy certain non-degeneracy conditions. Our local analysis is applicable also to subsequent generic catastrophe points. ied in this case in Jenkins et al. (2014), where it is shown that the solution develops oscillations for t > 0, without developing a point of elliptic umbilic catastrophe.
The last statement of the Main Conjecture refers to the behaviour of a generic perturbed solution near the point of gradient catastrophe of the unperturbed one. Our main goal is to find an asymptotic description for the dispersive regularization of the elliptic umbilic singularity or the cusp catastrophe when the dispersive terms are added; i.e. we want to describe the leading term of the asymptotic behaviour for → 0 of the solution to (1.3) near the critical point, say (x 0 , t 0 ), of a generic solution to (1.1).
At the point of catastrophe, the solutions u (x, t), v(x, t) to the Cauchy problem (1.1), (1.6) remain continuous, but their derivatives blow up. The generic singularities of solutions to the quasilinear systems (1.1) are classified as follows (Dubrovin 2006(Dubrovin , 2008. • If the system (1.1) is hyperbolic, h uu h vv > 0, then the generic singularity is a point of cusp catastrophe or more precisely the Whitney W3 (Whitney 1955) singularity. • If the system (1.1) is elliptic, h uu h vv < 0, then the generic singularity is a point of elliptic umbilic catastrophe. This codimension 2 singularity is one of the real forms labelled by the root system of the D 4 type in the terminology of Arnold et al. (1993).
Elliptic umbilic singularities appear in experimental and theoretical studies of diffraction in more than one spatial dimension (Berry et al. 1979), in plasma physics (Slemrod 2002;Sikivie 1999), in the Hele-Shaw problem (Martínez-Alonso and Medina 2009), and also in random matrices (Fokas et al. 1991;Bertola and Tovbis 2011). Formation of singularities for general quasilinear hyperbolic systems in many spatial dimensions has been considered in Alinhac (1995), Majda (1984) (see Manakov and Santini 2011 for an explicit example). For the particular case of 2 × 2 systems we are mainly dealing with, the derivation of the cusp catastrophe was obtained for C 4 initial data in Kong (2002), see also Dubrovin (2008) for an alternative derivation.
Let us return to the Cauchy problem for the perturbed system (1.3) with the same initial data (1.5). The fundamental idea of universality first formulated in Dubrovin (2006) for scalar Hamiltonian PDEs suggests that, at the leading order of asymptotic approximation, such behaviour does depend neither on the choice of generic initial data nor on the choice of generic Hamiltonian perturbation. One of the goals of the present paper is to give a precise formulation of the universality conjecture for a quite general class of systems of Hamiltonian PDEs of order two (for certain particular subclasses of such PDEs the universality conjecture has already been formulated in Dubrovin 2008).
The general formulation of universality introduced in Dubrovin (2006) for the case of Hamiltonian perturbations of the scalar nonlinear transport equation and in Dubrovin (2008) for Hamiltonian perturbation of the nonlinear wave equation says that the leading term of the multiscale asymptotics of the generic solution near the critical point does not depend on the choice of the solution, modulo Galilean transformations, and rescalings. This leading term was identified via a particular solution to the fourth-order analogue of the Painlevé-I (P I ) equation (the so-called P 2 I equation). Earlier the particular solution to the P 2 I equation proved to be important in the theory of random matrices (Moore 1990;Brézin et al. 1990); in the context of the so-called Gurevich-Pitaevsky solution to the KdV equation, it was derived in Kudashev and Suleimanov (1996). The existence of the needed smooth solution to P 2 I has been rigorously established in . Moreover, it was argued in Dubrovin (2006Dubrovin ( , 2008 that the shape of the leading term describing the critical behaviour is essentially independent of the particular form of the Hamiltonian perturbation. Some of these universality conjectures have been supported by numerical experiments carried out in Grenier (1998), Dubrovin et al. (2011). The rigorous analytical proof of this conjecture has been obtained for the KdV equation in Claeys and Grava (2009) for analytic initial data decreasing at infinity sufficiently fast so that inverse scattering is applicable.
In Dubrovin et al. (2009), the universality conjecture for the critical behaviour of solutions to the focusing cubic NLS has been formulated, and in Dubrovin (2008), the universality conjecture has been extended to other integrable Hamiltonian perturbations of elliptic systems. The universality conjecture in this case suggests that the description of the leading term in the asymptotic expansion of the solution to the focusing NLS equation in the semiclassical limit, near the point of elliptic umbilic catastrophe, is given via a particular solution to the classical Painlevé-I equation (P I ), namely the tritronquée solution first introduced by Boutroux (1913) one hundred years ago; see Joshi and Kitaev (2001), Kapaev (1995Kapaev ( , 2004 regarding some important properties of the tritronquée solution and its characterization in the framework of the theory of isomonodromy deformations. The smoothness of the tritronquée solution in a sector of the complex z-plane of angle | arg z| < 4π/5 conjectured in Dubrovin et al. (2009) has only recently been proved in Costin et al. (2014). Other arguments supporting the universality conjecture for the cubic focusing NLS case were found in Bertola and Tovbis (2013). Namely, the validity of a modified version of the conjecture has been established in the important paper (Bertola and Tovbis 2013) where the authors have considered -dependent (in the slow variables obtained after the Madelung transform) initial data built from the ad hoc semiclassical asymptotics of the spectral data. For particular initial data, namely u(x, 0, ) = sech 2 x and v(x, 0, ) = −μ log cosh x where the slow variables u (x, t, ) and v(x, t, ) are defined in (5.2) the conjecture has been proved in its original form (Bertola and Tovbis 2013). A proof of the original conjecture of Dubrovin et al. (2009) with -independent generic initial data remains an open problem, to the best of our knowledge. In this paper, we extend these ideas to the more general class of systems of the form (1.3). Our main goal is a precise formulation of the following conjectural statement.
• The solution of the generic system (1.3) with generic -independent smooth initial data near a point of cusp catastrophe of the unperturbed hyperbolic system (1.1) is described in the limit → 0 by a particular solution to the P 2 I equation.
• The solution of the generic system (1.3) with generic -independent analytic initial data near a point of elliptic umbilic catastrophe of the unperturbed elliptic system (1.1) in the limit → 0 is described by the tritronquée solution to the P I equation.
An important aspect of the above conjectures is the existence of the solution of the perturbed Hamiltonian systems (1.3) for times t up to and slightly beyond the critical time t 0 for the solution of the unperturbed system (1.1). The study of the local or global well-posedness of the Cauchy problem for the full class of Eq. (1. 3) remains open even though a large class of equations has been studied; see, for example, Ginibre and Velo (1979) or Tao (2006), Linares andPonce (2009), Bourgain (1999) for a survey of the state-of-the-art. For finite , it is known that the solution of the Cauchy problem of certain classes of equations of the form (1.3) develops blow-up in finite time; see, for example, Sulem and Sulem (1999), Kenig and Merle (2006). For the class of equations of the form (1.2) and initial data such that the solution develops a blow-up in finite time t B , we consequentially conjecture that, for sufficiently small , the blow-up time t B is always larger than the critical time t 0 of the dispersionless system. The blow-up behaviour of solutions to certain class of equations, like the focusing NLS equation, has been studied in detail in Merle and Raphael (2004); however, the issue of the determination of the blow-up time remains open. For the particular case of the quintic focusing NLS equation, we claim that the blow-up time t B which depends on is close in the limit → 0 to the time of elliptic umbilic catastrophe, more precisely the ratio (t B − t 0 )/ 4 5 is asymptotically equal to a constant that depends on the location of the first pole of the P I tritronquée solution on the negative real axis.
This paper is organized as follows. In Sect. 2, we single out the class of Hamiltonian systems (1.2) and we recall the procedure of obtaining solutions of the system (1.1) by a suitable form of the method of characteristics. In Sect. 3, we study the generic singularity of the solutions to (1.1) and describe the conjectural behaviour for the generic solution of a Hamiltonian perturbation (1.3) of the hyperbolic system (1.1) in the neighbourhood of such singularity. The same programme is realized in Sect. 4 for Hamiltonian perturbations of an elliptic system of the form (1.1). In Sect. 5, we consider in more details the above results for the generalized nonlinear Schrödinger equation (NLS) and the non-local NLS equation, and in Sect. 6, we study analytically some particular solutions of the system (1.1) up to the critical time t 0 for the generalized NLS equation. In Sects. 7-9, we present numerical evidences supporting the validity of the above conjectures.

Hamiltonian Systems
In this section, we identify the class of Hamiltonian equations we are interested in. Let us consider the class of systems of Hamiltonian PDEs of the form where we are taking the sum over repeated indices. The system of coordinates on the space of dependent variables can be chosen in such a way that the Poisson bracket takes the standard form (Dubrovin and Novikov 1989) where η i j is a constant symmetric non-degenerate matrix. Choosing a Hamiltonian in the form one obtains the following representation of the system (2.1) This yields, in particular, that Let us observe that a nonlinear change of dependent variables brings the Poisson bracket (2.2) to the form (Dubrovin and Novikov 1989) where the symmetric tensorg is a (contravariant) metric of zero curvature (not necessarily positive definite) and is expressed via the Christoffel coefficients of the Levi-Civita connection for the metric i j Any Hamiltonian system with n dependent variables can be locally reduced to the standard form (2.4), (2.3) by the action of the group of generalized Miura transformations (Degiovanni et al. 2005;Getzler 2002) changing the dependent variables as follows We will now concentrate on the case of a second-order Hamiltonian system, n = 2. It will be assumed that the metric η i j in the coordinates (u, v) has the canonical antidiagonal form Thus, the Hamiltonian system with a Hamiltonian H = H [u, v] reads . (2.10) A general perturbation of degree 2 of the Hamiltonian H 0 takes the form u, v) are some smooth functions. A simple calculation yields the following explicit form of the Hamiltonian flow where The linear terms in can be eliminated from Eq. (2.12) by a canonical transformation, as it follows from defined by the time-shift generated by the Hamiltonian (2.14) transforms the Hamiltonian system (2.10) to a system of the same form (2.15) The new HamiltonianH defined bỹ Proof By definition, one has So the HamiltonianH has the form (2.16) with The statement of Lemma readily follows from these expressions.
Thus, any Hamiltonian of the form (2.11) can be reduced to the form where the terms of order have been eliminated by a canonical transformation. The system of the form (2.12) can then be reduced to the form (1.3) (see above). Let us compute the general solution to the leading term of (1.3) obtained by setting (2.18) We will consider systems for which the eigenvalues h uv ± √ h uu h vv of the above matrix are distinct, namely h uu h vv = 0.
We will deal with smooth initial data only.
The following version of the classical hodograph transform will be used for the local description of non-degenerate solutions. such that on a sufficiently small neighbourhood of this point the following two equations hold identically true, (u(x, t), v(x, t)) . (2.21)

Conversely, given any solution f
locally defined by the system of the implicit function theorem holds true at the point (u 0 , v 0 ) such that With the help of (2.18), one arrives at (2.25) This system can be recast into the form Hence, there locally exists a pair of functions φ = φ (u, v), ψ = ψ(u, v) such that This implies Therefore, a function f = f (u, v) locally exists such that Thus, The linear PDE (2.20) as well as the implicit function Eq. (2.21) readily follows. The proof of the converse statement can be obtained by a straightforward computation using the expressions derived with the help of the implicit function theorem (here D is the determinant (2.23)).

Remark 2.3
Observe invariance of the implicit function Eq. (2.21) with respect to transformations of the dependent variables (u, v) preserving the antidiagonal form (2.9) of the metric η.

Hyperbolic Case
In this section, we study solutions to the system (1.3) when the unperturbed systems (2.18) is hyperbolic. We will restrict our analysis to smooth initial data. We first derive the generic singularity of the solution to the hyperbolic systems of the form (2.18), and then, we study the local behaviour of the solution of the system (1.3) with > 0 near such a singularity. Our first observation is that, in a suitable system of dependent and independent coordinates, the system of equations (1.3) decouples in a double scaling limit near the singularity into two equations: one ODE and one PDE equivalent to the Korteweg de Vries equation. We then argue that the local behaviour of the solution of (1.3) near the singularity of the solution to (2.18), in such a double scaling limit is described by a particular solution to the P 2 I equation. The system (2.18) is hyperbolic if the eigenvalues of the coefficient matrix The proof of the following statement is straightforward.

Lemma 3.1 The hodograph Eq. (2.22) can be rewritten in the form
Denoting by r ± the Riemann invariants of the system, we get for their differentials where κ ± = κ ± (u, v) are integrating factors. The leading order system (2.18) becomes diagonal in the coordinates r + , r − ; i.e.
It is convenient to write the hodograph Eq. (3.3) in terms of the Riemann invariants r = (r + , r − ) where the functions μ ± = μ ± (r ) must satisfy the linear system equivalent to (2.20). The functions μ + (r ), μ − (r ) have to be determined from the system (3.8) along with the conditions at t = 0 for given C ∞ initial data u(x, 0), v(x, 0) for the system (2.18). It is easy to see that such a solution is determined uniquely, and it is smooth on any interval of monotonicity of both initial Riemann invariants r + (x), r − (x) provided the values of the characteristic velocities λ ± (r (x)) := λ ± (r + (x), r − (x)) on the initial curve are distinct then the hyperbolic system is called linearly degenerate. In this case, there exists (Majda 1984) a global solution r ± (x, t) for all t ≥ 0. In the present paper, it is always assumed that the system is not linearly degenerate. Furthermore, as in the scalar case, in order to have a point of gradient catastrophe, we need to assume that the initial data r ± (x) is not monotone (increasing or decreasing depending on the sign of the characteristics speeds λ ± ) (Kong 2002). In Kong (2002), solutions to strictly hyperbolic system of the form (3.6) with C 1 initial data are considered such that the first point of gradient catastrophe (x 0 , t 0 ) occurs at x 0 < ∞. The class of initial data satisfying this requirement is quite big including, for example, C ∞ periodic initial data, compactly supported initial data, or initial data such that d dx r ± (x) → 0 as |x| → ∞.
Below we will assume that the smooth initial data r ± (x) are such that the solution of the Cauchy problem (3.6) has its first point of gradient catastrophe for the Riemann invariant r − (x, t).
Our first goal is to derive a normal form of the system (3.6) near a point of gradient catastrophe of the leading term (3.6). The limiting values of the solutions r ± (x, t) to (3.6) at the point of gradient catastrophe (x 0 , t 0 ) will be denoted Let us also introduce the shifted dependent variables denoted as (3.10) and the notation etc., for the values of the coefficients and their derivatives at the point of catastrophe.
In the generic situation, the x-derivative of only one of the Riemann invariants becomes infinite at the point of catastrophe. To be more specific, let us assume We say that the point of catastrophe (3.11) is generic if and, moreover, the graph of the function r − (x, t 0 ) has a non-degenerate inflection point at x = x 0 . Introduce characteristic variables at the point of catastrophe. One can represent the functions r ± = r ± (x, t) as functions of (x + , x − ). Let us redenoter ± = r ± (x + , x − )−r 0 ± the resulting transformed functions. It will be convenient to normalize 3 the Riemann invariants in such a way that (3.14) Here and below, we will use notations for partial derivatives with respect to Riemann invariants r ± similar to those in (3.12) Choose a pair of sufficiently small real numbers X + , X − satisfying Lemma 3.2 For a generic solution to the system (3.6) and for arbitrary sufficiently small real numbers X + , X − satisfying (3.15), there exist the limits (3.16) The limiting functions satisfy the system Proof A generic solution to (3.6) for t < t 0 is determined from the implicit function Eq. (3.7). At the point of catastrophe of the Riemann invariant r − , one has The point is generic if, along with the condition (3.12), one also has Expanding Eq. (3.7) in Taylor series near the point (r 0 + , r 0 − ) and using (3.8) one obtains, after the rescaling the relations Applying a similar procedure directly to the system (3.6), one obtains the following Lemma 3.3 The limiting functions (3.16) satisfy the following system of PDEs where the constant β is defined in (3.18).

Remark 3.4
The study of solutions of integrable partial differential equation in the limit → 0 can be tackled via a Riemann-Hilbert (RH) formulation of the Cauchy problem and then (Deift and Zhou 1993) steepest descent analysis. Technically such an analysis can be quite involved, and so far it has been rigorously performed just for few integrable equations which have a hyperbolic limit as → 0, like the defocusing nonlinear Schrödinger equation (Jin et al. 1994;DiFranco and Miller 2008), the Korteweg-de Vries equation (Deift et al. 1997), the Toda lattice (Deift et al. 1999b;Deift and McLaughlin 1998), and few others. Let us remind, for experts in the field, that the point of gradient catastrophe (3.19) and (3.20) corresponds, in the RH analysis, to a type III singularity for a complex function of a complex variable called g-function (Deift et al. 1999b). In this case, the g-function has a zero of order 7 2 at one of the end points of its support.
Let us proceed to the study of solutions to the perturbed system (1.3). Choosing the Riemann invariants r ± = r ± (u, v) of the leading term as a system of coordinates on the space of dependent variables, we obtain the system (1.3) in the form (3.25) We are now ready to prove the first result of this section.
Theorem 3.5 Let r ± = r ± (x, t; ) be a solution to the system (3.24) defined for • there exists a time t 0 satisfying 0 < t 0 < τ such that for any 0 ≤ t < t 0 and sufficiently small |x − x 0 |, the limits exist and satisfy the system (3.6). Let us consider the solution r ± (x, t) represented in the hodograph form (3.7), and assume that • it has a gradient catastrophe at the point (x 0 , t 0 ) of the form described in Lemma 3.3; • there exist the limits (3.26) • the constants α, β, γ in (3.18) do not vanish and β γ > 0; The limiting function R + = R + (X + , X − ; ) is given by the formula A solution r ± (x, t; ) to the system (3.24) with a hyperbolic leading term satisfying the assumption (3.12) along with will be called generic.

Conjecture 3.6 A generic solution to the -independent
with α, β, γ , and ρ defined in (3.18) and (3.27), respectively, and where U = U (X, T ) is the smooth solution to the P 2 I equation for each fixed T ∈ R.
The existence of such solution to the P 2 I equation has been conjectured in Dubrovin (2006) (for T = 0, such a conjecture has already been formulated in Brézin et al. 1990) and proved in . See Fig. 1 below for a plot of such solution in the (X, T ) plane.
Lemma 3.7 The function U (X, T ) satisfies also the KdV equation Proof First, it is easy to check that Eqs. (3.33) and (3.35) are compatible (cf Moore 1990; Kudashev and Suleimanov 1996;Dubrovin 2006). That means that, given a solution U (X, T ) to the KdV equation (3.35) such that, for some T 0 the function U (X, T 0 ) satisfies Eq. (3.33) with T = T 0 , then U (X, T ) satisfies Eq. (3.33) for all values of the parameter T . Choose an arbitrary real value T 0 ; denote U 0 (X ) the (unique) smooth solution to the Eq. (3.33) with T = T 0 such that Due to results of Menikoff (1972), there exists a unique smooth solution U (X, T ) to the KdV Eq. (3.35) satisfying the initial condition U (X, it will satisfy the asymptotic (3.34). Due to the above compatibility statement, the solution to KdV will satisfy Eq. (3.33) for all real T .
Thus, the asymptotic formulae (3.32) meet the following two conditions: • for t < t 0 , the solution (3.32) tends to the hodograph solution (3.17) as → 0; • near the point of break-up, the rescaled Riemann invariant r − approximately satisfies the KdV Eq. (3.28) while the rescaled Riemann invariant r + admits an approximate representation (3.29). Indeed, choosing So, after rescaling of the characteristic variables Moreover, for largeX − and negativeX + it behaves like the root of the cubic equation is a solution to KdV satisfying these properties. Returning to the original variablesr − , x ± , one arrives at the formula (3.32).

Elliptic Case
In this section, we study solutions to the system (1.3) when the unperturbed systems (2.18) is elliptic. We will restrict our analysis to analytic initial data. We first derive the generic singularity of the solution to the elliptic systems of the form (2.18), and then, we study the local behaviour of the solution of the system (1.3) with > 0 near such a singularity. We argue that such behaviour in a double scaling limit is described by the tritronquée solution to the P I equation.
Let us now proceed to considering the elliptic case for the system (2.18), namely The initial data u(x, 0) and v(x, 0) are analytic functions. The Riemann invariants and the characteristic speeds are complex conjugate (the asterisk will be used for the complex conjugation), At the point of elliptic break-up of a solution, written in the form (3.7), the following two complex conjugated equations hold (4.4) In Sect. 6, we provide several examples of initial data for which the Eq. (4.4) have a solution. However, to the best of our knowledge, the problem of characterizing a class of initial data such that the solution of the elliptic system (2.18) has a point of elliptic umbilic catastrophe, is still open. The characteristic variables at the point of catastrophe are defined as and are also complex conjugate. One can represent the functions r ± = r ± (x, t) as functions of (x + , x − ). Let us redenoter ± = r ± (x + , x − ) − r 0 ± the resulting shifted and transformed Riemann invariants.
exist and satisfy the quadratic equation In the sequel, it will be assumed that a ± = 0 (4.9) (this condition will be added to the genericity assumptions).
Proof Differentiating the hodograph relations (3.7), one obtains Moreover, differentiating (3.8) one finds that Hence, due to (4.4), all these combinations of the second derivatives vanish at the break-up point. Expanding the hodograph Eq. (3.7) in Taylor series near the point of catastrophe, one easily arrives at (4.7).
Remark 4.2 Also in this case as in remark 3.4, the study of solutions of integrable partial differential equation in the limit → 0 can be tackled via a Riemann-Hilbert formulation of the Cauchy problem and then with (Deift and Zhou 1993) steepest descent analysis. Such an analysis can be quite involved, and it has been rigorously performed to the best of our knowledge only for the cubic focusing NLS equation (Kamvissis et al. 2003;Tovbis et al. 2004). In this case, the point of elliptic umbilic catastrophe (4.7) corresponds to a singularity of a certain type of a complex phase function related to the so-called g-function. This complex phase function has a zero of order 5 2 at the point of elliptic umbilici catastrophe (Bertola and Tovbis 2013).
Choosing Riemann invariants r ± = r ± (u, v) of the leading term as a system of coordinates on the space of dependent variables, and x ± as independent variables, the system (1.3) takes the form (4.11) As above we will denoter ± =r ± (x + , x − ; ) a shifted generic solution to the system (4.10) with -independent initial data written as functions of the complex conjugated linearized characteristic variables (4.5). Like above, we will be interested in the multiscale expansion of these complex conjugated functions (4.12) We will now show that the existence of such expansions implies that the leading term is a holomorphic/antiholomorphic function satisfying an ODE.
Theorem 4.3 Let r ± (x, t, ; ) be a solution of the system (4.10) such that there exist the limits Then, the function R + = R + (X + , X − ; ε) satisfies the Cauchy-Riemann equation and also the equation where c + is a holomorphic function of X + such that , as X + → ∞ and δ > 0. (4.16) Here C + + has been defined in (4.11). The function R − = R − (X − ; ε) is antiholomorphic and satisfies the complex conjugate of (4.15). The function R + (X + , X − ; ε) satisfies the equation where C + − has been defined in (4.11). The function R − (X + , X − ; ε) satisfies the complex conjugate of the above equation.
Proof In order to prove the theorem, it is sufficient to plug the expansion (4.12) and (4.13) into Eq. (4.10) giving the following expansions Since λ 0 + = λ 0 − , from the leading term, it readily follows that Separating holomorphic and antiholomorphic parts in the terms of order O(1), one arrives at Eqs. (4.15), (4.17), and their complex conjugates. Equation (4.15) must have a solution with asymptotic behaviour determined by (4.7), namely This immediately gives that c + is an analytic function of X + with asymptotic behaviour at infinity (4.21) Assuming c + = const, we arrive at an ODE for the function R + = R + (X + ) equivalent to the P I equation, with asymptotic behaviour (4.20). The complex conjugate of the above equation gives the corresponding P I equation for R − = R − (X − ). If we linearize the increments of the Riemann invariants, we obtain (4.23) For simplicity, we normalize the constant κ 0 ± to (4.24) From (4.22) and (4.23), we arrive at the following.

Conjecture 4.4 The functions u(x, t, )
and v(x, t, ) that solves the system (1.3) admit the following asymptotic representation in the double scaling limit x → x 0 , t → t 0 and → 0 in such a way that Here a + and C + + (u 0 , v 0 ) have been defined in (4.8) and (4.11), respectively, and = (ξ ) is the tritronquée solution to the P I equation (4.29) The smoothness of the solution of (4.28) with asymptotic condition (4.29) in a sector of the complex z-plane of angle | arg z| < 4π/5 conjectured in Dubrovin et al. (2009) has only recently been proved in Costin et al. (2014). For a plot of such solution in the complex plane, see Dubrovin et al. (2009).
is given by the complex conjugate of (4.26). Remark 4.6 We write the constant a + in the form with q > 0 and ψ ∈ [−π, π]. One can check that when ψ = 0 and t = t 0 , the quantity ξ defined in (4.27) has to be purely imaginary, and this gives a rule for the selection of the fifth root, namely Note that the angle of the line ξ = ξ(x − x 0 ) for fixed t is equal to π 2 + ψ 5 with ψ ∈ [−π, π]; thus, the maximal value of arg ξ is equal to 7 10 π < 4 5 π .
In the next subsection, we consider an alternative derivation of the P I equation for a subclass of Hamiltonian PDEs having the structure of a generalized nonlinear Schrödinger equation.

P I Equation and Approximately Integrable PDEs
In this subsection, we give an alternative derivation of the Conjecture 4.4 for the nonlinear wave equation u tt − ∂ 2 x P (u) = 0. It can be represented in the form (2.18) as a second-order system with the Hamiltonian after eliminating the dependent variable v. Here P(u) is a smooth function satisfying P (u) < 0. More generally the arguments given below will work for any Hamiltonian system with the Hamiltonian commuting with H nlin . Its density must satisfy We will see below that, in particular, a very general family of nonlinear Schrödinger equations belongs to this subclass. The condition P (u) < 0 (4.34) guarantees that the unperturbed quasilinear system is elliptic. A local solution of the system (1.1) with h(u, v) satisfying (4.33) for given analytic initial data u 0 (x), v 0 (x) takes the form where the function f = f (u, v) satisfies equation and the condition (4.37) The equation for determining the point of elliptic umbilic catastrophe characterized by Eq. (4.4) in the variables u and v takes the form and the constants a ± take the form To study the critical behaviour of solutions of (2.12), we first restrict ourselves to approximately integrable cases in the sense of Dubrovin (2008). Recall that a perturbation of the Hamiltonian system (2.12) with a Hamiltonian H of the form (2.17) is called integrable up to corrections of order 3 if, for any first integral F 0 = f (u, v) dx of the unperturbed system (2.12) there exists a perturbed functional with Hamiltonian H h = Dh dx and Hamiltonian density Dh given by Furthermore, a class of solutions of the system (4.41) characterized by an analogue of the string equation is given by the following theorem.
Theorem 4.8 (Dubrovin 2008) The solutions to the string equation also solve the Hamiltonian equations is another solution to f uu = P (u) f vv , and h := ∂h ∂v .
We remark that (4.44) is a system of coupled ODEs for u and v having t has a parameter. We can apply to the system (4.44) the rescaling (4.12). Let us first introduce the Riemann invariants for the Hamiltonians H 0 satisfying (4.33) Choosing the Riemann invariants r ± = r ± (u, v) as a systems of coordinates on the space of dependent variables, one can write the string Eq. (4.44) in the form with λ ± as in (4.3) and where the coefficientsC ± ± are as in (4.11) with a = a (u, v), b = b(u, v) and c = c(u, v) obtained by comparing the Hamiltonian H f − t H h to the general form (2.17).

Proposition 4.9
The string Eq. (4.44) in the scaling (4.12) reduces to the P I equation Proof Using the Riemann invariants as a system of dependent coordinates, the string Eq. (4.44) takes the form (4.46). Changing the independent coordinates (x, t) to (x + , x − ) defined in (4.5) and performing the scalings (4.12), one obtains for k → 0 where P 0 = P(u 0 , v 0 ). Requiring the compatibility of the leading order expansion of the string equation with the leading order expansion of the system (4.10), we get that (4.19) has to be compatible with (4.48), namely which is equivalent to the P I equation. We observe that the quantity ρ 0 u + i |P 0 |ρ 0 v can be rewritten in the form with C + + as in (4.11). In a similar way, one can write the complex conjugate. Therefore, Eq. (4.49) can be written in the form (4.47). Proof It is sufficient to show that for given a = a (u, v), b = b(u, v), and c = c(u, v), one can find ρ u (u, v), ρ v (u, v), and s 3 (u, v) such that at the critical point (u 0 , v 0 ) the following identities hold: (4.52) The constants ρ 0 u and ρ 0 v can be chosen in an arbitrary way since they solve the secondorder Eq. (4.43), and s 3 (u, v) is an arbitrary function. The system (4.52) is solvable for ρ 0 u , ρ 0 v , and s 0 3 as a function of a 0 , b 0 , c 0 . For a given initial datum, the solutions of two different Hamiltonian perturbations of the form (2.12) with the same unperturbed Hamiltonian density h (u, v) satisfying h uu = h vv P (u) have the same approximate solution for t < t 0 . From our Conjecture 4.4, it follows that the solutions near the critical point have the same leading asymptotic expansion if the coefficients of the two systems satisfy at the critical point the relation (4.52).

An Example: Generalized Nonlinear Schrödinger Equations
Let us now consider the example of generalized nonlinear Schrödinger (NLS) equa- is a complex variable and V is a smooth function monotone increasing on the positive real axis. The case V (u) = u is called cubic NLS, and the case V (u) = u 2 /2 is called quintic NLS, and so on. The case with positive sign in front of the potential V is the so-called focusing NLS, while the negative sign corresponds to the defocusing NLS. For sufficiently regular V and for finite > 0, the initial value problem of the defocusing NLS equation is globally well posed in some suitable functional space, see Ginibre and Velo (1979), Bourgain (1999), and references therein, while the solution of the initial value problem of the focusing case is globally well posed when the nonlinearity V |ψ| 2 = |ψ| 2 (Ginibre and Velo 1979). Equation (5.1) can be rewritten in the standard Hamiltonian form (2.10) with two real-valued-dependent functions, the so-called Madelung transform (the star stands for the complex conjugation). Then, Eq. (5.1) reduces to the system of equations 3) The above system can be written in the Hamiltonian form 5 (5.5) The semiclassical limit of this system is of elliptic or hyperbolic type, respectively, provided V (u) is a monotonically increasing smooth function on the positive semiaxis.
Another interesting NLS-type model is given by the non-local NLS equation (Conti et al. 2009;Rasmussen et al. 2005;Ghofraniha et al. 2007), where η is a positive constant. In the slow variables u and v, this non-local NLS model can be equivalently written as Writing θ from the last equation as the formal series θ = u + 2 η u x x + 4 η 2 u x x x x + · · · and keeping terms up to order 2 only, one arrives at a system of the above class

The non-local NLS can be written in the Hamiltonian form (5.4) with the Hamiltonian H = h dx and
(5.11) The above Hamiltonian coincides with the one of the cubic NLS when η = 0.
We are going to study the critical points of the solutions of the system (5.6) for some initial data and then the solutions of Eqs. (5.1) or (5.7) for the same data near the critical points of the solution of (5.6). We first treat the hyperbolic case.

Defocusing Generalized NLS
The Riemann invariants and the characteristic velocities of Eq. (5.6), in the hyperbolic case, are The general solution to (5.6) can be represented in the implicit form where the function f = f (u, v) solves the linear PDE of the form (2.20) (5.14) The coordinates (u 0 , v 0 ) of the point of a generic break-up of the second Riemann invariant r − can be determined from the system In the Riemann invariants, the system (5.3) reads The asymptotic representation of the shifted Riemann invariants is given as a function of the shifted characteristic variables (5.17) In particular for the non-local defocusing NLS equation, the shifted Riemann invariants as functions of the shifted characteristic variables behave in the vicinity of the point of gradient catastrophe as in (3.32) with α, β, and γ as in (5.17) with V 0 = 1 and ρ and σ given by

Focusing Generalized NLS
The Riemann invariants and the characteristic velocities of system (5.6) in the elliptic case are

The general solution of (5.6) is obtained via the hodograph equations
where the function f (u, v) solves the linear equation The point of elliptic umbilic catastrophe is determined by the equations (5.19) and the conditions f uu = 0, t + f uv = 0. (5.21) The asymptotic formula (4.26) near the point of elliptic umbilic catastrophe takes the form Remark 5.1 In the formula (5.22), the convention for choosing the fifth root is defined by the following condition: For symmetric initial data and t = t 0 , the argument of the tritronquée solution has to be purely imaginary. So, defining (5.25) Remark 5.2 In the focusing non-local NLS model (5.7), the behaviour of the solution near the point of elliptic umbilic catastrophe is given by the expression (5.22) with a + and λ 0 +,+ as in (5.23) and For η = 0, such a formula was derived in an equivalent form in Dubrovin et al. (2009).

Studying Particular Solutions
The present section is devoted to the comparison of solutions to the defocusing and focusing NLS equations with their unperturbed counterparts near the critical points of solutions of the unperturbed system with, respectively, the asymptotic formula (3.32) and (5.22). We consider various examples of nonlinear potentials V and initial data. Let us consider the Cauchy problem If the initial data ϕ ± (x) are bounded analytic functions of x, then in virtue of the Cauchy-Kowalevskaya theorem (see Bressan 2000) r ± (x, t) are analytic functions in the variable x up to the time t < t 0 where t 0 is the time of gradient catastrophe. The implicit solution of (6.1) is given by the hodograph equations as where μ ± solves the system of linear PDEs equivalent to (5.14) with the constraint (6.4)

Defocusing Cubic NLS
The cubic NLS equation written as corresponds to the case V (u) = u, and the Riemann invariants and the characteristics velocities (5.12) take the form Let us consider an initial datum rapidly going to a constant value at infinity The solution of the corresponding quasilinear system (3.6) is obtained as described below. Let us suppose that the initial datum ϕ + (x) has a single-positive hump at x M and that ϕ − (x) has a single-negative hump at x m ≤ x M , and denote by h + L/R (r + ), the inverse of the increasing and decreasing part of ϕ + (x) and by h − L/R (r − ), the inverse of the decreasing and increasing part of ϕ − (x), respectively. Since λ + > λ − , it follows that x M (t) ≥ x m (t) for all t ≥ 0. In order to obtain the quantities μ ± (r + , r − ), we use the formula by Tian and Ye (1999): dτ dx (6.7) For different choices of initial data, more complicated relations can be obtained. Within the interval of monotonicity of the function ϕ ± , the solution (6.5) can be written also in the equivalent form (Tian and Ye 1999) where x(τ ) is the inverse function of ϕ ± (x) in the interval of monotonicity. For the particular case v(x, 0) = 0, u(x, 0) = a 2 sech 2 x, one has The critical point is obtained by the two Eqs. (6.2) together with which give the equations Solving the above two equations together with (6.2) yields (6.12) The constants b = 12 ρ β , β, γ , and α defined in (5.17) at the critical time are given

Defocusing Quintic NLS
Let us now proceed to the case V (u) = u 2 /2. The Riemann invariants of the quintic defocusing NLS are given by The Eq. (5.6) reduce to the two decoupled Riemann wave equations ∂ t r ± + r ± ∂ x r ± = 0, which can be solved by the method of characteristics. For the initial data r ± (x, 0) = ρ ± (x), one has the solution in implicit form The point of gradient catastrophe is determined by the conditions where F ± is the inverse of the decreasing part of the initial data ρ ± (x). The constants b = 12 ρ β , β and σ defined in (5.17) at the critical time are given by (6.14) The constants α and γ in (5.17) depend on the initial data and are evaluated for several initial data below. Symmetric Initial Data This name will be applied to the class of NLS initial data ψ 0 (x) := ψ(x, 0) satisfying the condition (the asterisque stands for complex conjugation) or, equivalently, The initial values of the Riemann invariants then satisfy If none of the conditions (6.15) holds true, then the solution will be called asymmetric. We begin with considering the following symmetric initial data with A a positive constant. For such initial data, both r ± have a point of gradient catastrophe. The evolution in time of the decreasing part of r + (x, t) gives The point of gradient catastrophe is given by The second Riemann invariant r − (x 0 The constants γ and α in (3.18) take the form The evolution in time of the decreasing part of r − (x, t) gives with F − (r − ) as in (6.17). The point of gradient catastrophe is given by The constants γ and α in (3.18) take the form as in (6.16). "Dark" Initial Data. We consider the initial data In the evolution of this initial data, two points of gradient catastrophe occur, one at x 0 + < 0 for the Riemann invariant r + and one at x 0 − > 0 for the Riemann invariant r − . For these initial data, the Riemann invariant r + (x, t) for x < x m , where x m is the point of the minimum of u, is determined by with critical point The constants γ and α in (3.18) take the form The evolution of r − (x, t) for x > x m , where x m is the point of minimum, is determined by the equation The point of gradient catastrophe occurs at The point r 0 The constants γ and α in (3.18) take the form

Focusing Cubic NLS
The case of the focusing cubic NLS equation was considered extensively in Dubrovin et al. (2009).
For the initial data ψ(x, t = 0) = A 0 sech x, or equivalently u = A 2 0 sech 2 x, v = 0, (6.18) the solution of the Eq. (5.6) in the elliptic case is given by (6.20) and the function f (u, v) takes the form The point of elliptic umbilic catastrophe is given by so that a + in (5.23) becomes a + = −i

Focusing Quintic NLS
The Riemann invariants of the equation are given by The Eq. (5.6) reduce to two uncoupled Riemann wave equations ∂ t r ± + r ± ∂ x r ± = 0.
For the symmetric initial datum the solution is given by where F is the inverse of the increasing part of the initial data (6.18), namely An equivalent result can be obtained considering the decreasing part of the initial data. Comparing (6.22) with (5.19), one has (6.23) and it easily follows that f uu + f vv = 0. The point of elliptic umbilic catastrophe is determined by the equations (6.22) and the condition The solution is given by The constants r and ψ in (5.24) are given by Asymmetric initial data Let us first consider symmetric initial data u(x) = sech x and v(x) = − tanh x. The solution defined by the hodograph transform takes the form where F is the inverse of the increasing part of the initial data r + (x, The breaking condition implies that the critical point is given by The constants r and ψ in (5.24) are given by The quantities in (5.23) take the form To obtain an initial datum which is manifestly not symmetric, we use the fact that, if F(r + ) is an analytic function, also d dr + F(r + ) is an analytic function, and therefore, solves the Laplace equation (5.20). We choose asymmetric initial data of the form (6.27) where F is given in (6.26) and F (r + ) = F(r + ), namely The time evolution of r + (x, t) is given by the hodograph equation In order to determine the point of elliptic umbilic catastrophe, it is sufficient to consider the solution of (6.28) together with the condition The real and imaginary parts of the Eqs. (6.28) and (6.29) give The solution of the above system determines the critical point (x 0 , t 0 ) and the values . The constants r and ψ in (5.24) are given by "Dark" Initial Data. We consider the initial data u(x, t = 0) = tanh 4 x and v = 0. For such initial data, the hodograph equations are where r + = v + iu. The break-up point is determined by the above complex equation together with the condition As in this case, it is not possible to obtain a simple analytic expression for the point of elliptic umbilic catastrophe (x 0 , t 0 ), and for r 0 + , r 0 − , they are determined numerically. The constants r and ψ that appear in (5.24) are given by

Numerical Methods
The numerical task in treating the semiclassical limit of the NLS equations consists in solving the NLS equations, the numerical evaluation of implicit solutions to certain ODEs, and the direct solution of ODEs of Painlevé type for a given asymptotic behaviour. The present section provides a summary of how these different tasks are solved numerically, and how the numerical accuracy is controlled.

NLS Equations
Critical phenomena are generally believed to be independent of the chosen boundary conditions. Thus, we study a periodic setting in the following. This also includes rapidly decreasing functions which can be periodically continued as smooth functions within the finite numerical precision. This allows to approximate the spatial dependence via truncated Fourier series which leads for the studied equations to large systems of ordinary differential equations (ODEs). Fourier methods are convenient because of their excellent approximation properties for smooth functions (the numerical error in approximating smooth functions decreases faster than any power of the number N of Fourier modes) and for minimizing the introduction of numerical dissipation which is important in the study of the purely dispersive effects considered here. In Fourier space, Eq. (5.1) have the formψ whereψ denotes the (discrete) Fourier transform of ψ, and L and N denote linear and nonlinear operators, respectively. The resulting system of ODEs consists in this case of stiff equations. A stiff system is essentially a system for which explicit numerical schemes as explicit Runge-Kutta methods are inefficient, since prohibitively small time steps have to be chosen to control exponentially growing terms. The standard remedy for this is to use stable implicit schemes, which require, however, the iterative solution of a system of nonlinear equations at each time step which is computationally expensive. In addition, the iteration often introduces numerical errors in the Fourier coefficients. The stiffness appears here in the linear part L (it is a consequence of the distribution of the eigenvalues of L), whereas the nonlinear part is free of derivatives. In the semiclassical limit, this stiffness is still present despite the small term in L. This is due to the fact that the smaller is, the higher wave numbers are needed to resolve the strong gradients. A possible way to deal with stiff systems are so-called implicit-explicit (IMEX) methods. The idea of IMEX is the use of a stable implicit method for the linear part of Eq. (7.1) and an explicit scheme for the nonlinear part which is assumed to be non-stiff. In Kassam and Trefethen (2005), such schemes did not perform satisfactorily for dispersive PDEs which is why we consider a more sophisticated variant here. Driscoll's idea (see Driscoll (2002)) was to split the linear part of the equation in Fourier space into regimes of high and low wave numbers. He used the fourth-order Runge-Kutta (RK) integrator for the low wave numbers and the lineary implicit RK method of order three for the high wave numbers. He showed that this method is in practice of fourth order over a wide range of step sizes. In Klein (2008), we showed that this method performs best for the focusing case. We use it here also for the defocusing case where it was very efficient, but slightly outperformed by so-called time-splitting schemes as in Bao et al. (2002Bao et al. ( , 2003. For a discussion of exponential integrators in this context, see Berland and Skaflestad (2005), Berland et al. (2007), Klein (2008). Numerical approaches to the semiclassical limit of NLS can be also found in Ceniceros (2002), Ceniceros and Tian (2002).
The accuracy of the numerical solution is controlled via the numerically computed conserved energy of the solution which is an exactly conserved quantity for NLS equations. Numerically, the energy E will be a function of time due to unavoidable numerical errors. We define E := |(E(t) − E(0))/E(0)|. It was shown in Klein (2008) that this quantity can be used as an indicator of the numerical accuracy if sufficient resolution in space is provided. The quantity E typically overestimates the precision by two to three orders of magnitude.
Since we are interested in an accuracy at least of order , we will always ensure that the Fourier coefficients of the final state decrease well below 10 −5 , and that the quantity E is smaller than 10 −6 (in general it is of the order of machine precision; i.e. 10 −14 ). Focusing NLS equations have a modulational instability due to the fact that they can be seen as a hyperbolic regularization of an elliptic semiclassical system for which initial value problems are ill-posed. In our context, this instability shows up in the form of spurious growing modes for high wave numbers. To address this problem, we use a Krasny filter (Krasny 1986), which means we put the Fourier coefficients with modulus below some threshold (typically 10 −12 ) equal to zero. Thus, the effect of rounding errors is reduced. In Klein (2008), it was pointed out that sufficient spatial resolution has to be provided to resolve the maximum of the solution close to the critical time to avoid instabilities. Thus, we use 2 14 to 2 16 Fourier modes, and 10 4 to 10 5 time steps for the computations.

Numerical Solution of the Semiclassical Equations
The solutions to the semiclassical equations are obtained in implicit form via hodograph techniques. These equations are of the form where the S i denote some given real function of the y i and x, t. The task is to determine the y i in dependence of x and t. To this end, we determine the y i for given x and t as the zeros of the function S := M i=1 S 2 i . This is done numerically via a Newton iteration which is very efficient for a sufficiently good initial iterate. This iteration has the advantage that it can be done for all values of x at the same time, i.e. in a vectorized way. Alternatively, we use the algorithm (Lagarias et al. 1988) pointwise to solve (7.3). We calculate the zeros to the order of machine precision. The residual of the equations provides a check of the numerical accuracy.

Painlevé Transcendents
The asymptotic solutions near the break-up point are given by pole-free solutions with a given asymptotic behaviour for x → ±∞ to the P I and the P 2 I equation. The standard way to solve these equations for large |x| is to give a series solution to the respective equation with the imposed asymptotics that is generally divergent. These divergent series are truncated at finite values of x, x l < x r at the first term that is of the order of machine precision. 6 The sum of this truncated series at these points is then used as boundary data, and similarly for derivatives at these points. Thus, the problem is translated to a boundary value problem on the finite interval [x l , x r ].
In Grava and Klein (2008), we used for the P 2 I solution a collocation method with cubic splines distributed as bvp4 with MATLAB, and the same approach in Dubrovin et al. (2009) for the tritronquée solution of P I . Note that the tritronquée solutions are constructed on lines in the complex plane in the sector where the solution is conjectured (see Dubrovin et al. 2009) to have no poles. As in Grava and Klein (2012), we use here a Chebyshev collocation method for both equations. The solution of the ODEs is sampled on Chebyshev collocation points x j , j = 0, . . . , N c which can be related to an expansion of the solution in terms of Chebyshev polynomials. The action of the derivative operator is in this setting equivalent to the action of a Chebyshev differentiation matrix on this space, see for instance (Trefethen 2000). The ODE is thus replaced by N c + 1 algebraic equations. The boundary data are included via a so-called τ method: The equations for j = 0 and for j = N c (for the fourth-order equation j = 0, 1, N c − 1, N c ) are replaced by the boundary conditions. The resulting system of algebraic equations is solved with a standard Newton method with relaxation which is necessary for the oscillatory P 2 I solution (there is no good initial iterate for the oscillatory solutions). The convergence of the solutions is in general very fast. We always stop the Newton iteration when machine precision is reached. Again the highest Chebyshev coefficients are taken as an indication of sufficient resolution of the solutions (they have to reach machine precision). An efficient solution of the ODE is especially important in the P 2 I case where the asymptotic solution to (3.33) has to be computed for many values of the parameter t. It can be seen in Fig. 1. For a more detailed discussion of this special P 2 I solution, also in the complex plane, see Kapaev et al. (2013). are compared to the corresponding semiclassical ones and for t ∼ t 0 to an asymptotic description in terms of a special solution to the second equation in the Painlevé-I hierarchy. We will consider the cubic and the quintic version of these equations. The cubic NLS is the only completely integrable equation studied in this paper. Since the results for both cubic and quintic are very similar in this case, we present a more detailed investigation for the non-integrable quintic NLS. We also study a non-local variant of the cubic NLS equation. Unless otherwise noted, the considered critical point is always at the centre of the figures.

Sech x Initial Data for the Cubic Defocusing NLS Equation
We will study the initial data ψ 0 (x) = sech x for several values of . In this case, there are two break-up points at ±x c with x c =∼ 2.2093 at the same time t 0 =∼ 1.5244. We will consider in the following always the break-up for negative values of x where the Riemann invariant r − = v − 2 √ u has a gradient catastrophe. In Fig. 2, the NLS solution, the semiclassical solution, and the P 2 I solution (3.32) can be seen at the critical time close to the critical point of the semiclassical solution.
The corresponding Riemann invariants can be seen in Fig. 3. For smaller , the agreement of NLS and semiclassical solution becomes better. We show the Riemann invariants r ± for = 10 −3 in Fig. 3. Notice that there are also oscillations in the invariant r + which stays smooth at this point in the semiclassical limit.

Sech x Initial Data for the Defocusing Quintic NLS Equation
We will first study the initial data ψ 0 (x) = sech x for values of of 0.1, 0.09,…, 0.01, 0.009,…, 0.001. In this case, there are two break-up points at ±x c with x c = ln(( √ 3 + 1)/ √ 2) + √ 3/2) ∼ 1.5245 at the same time t 0 = 3 √ 3/4 ∼ 1.2990. The solution up to the critical time can be seen in Fig. 4, where the defocusing effect of the equation can be recognized. The critical value of the Riemann invariants at the respective break-up point is ±2/3. We will consider in the following always the break- At the critical time, the difference of the Riemann invariants r − between the semiclassical solution and the solution to the focusing quintic NLS scales roughly as 2/7 . More precisely we find via a linear regression analysis for the logarithm of the difference − between NLS and semiclassical solution a scaling of the form ∝ a with a = 0.2952 (2/7 ∼ 0.2857) with standard deviation σ a = 0.0017 and correlation coefficient r = 0.9999. At the same point, the difference + between the Riemann invariants r + = v + u between the semiclassical and the NLS solution scales roughly as 4/7 as predicted by the theory. A linear regression analysis for the logarithm of the difference + gives a scaling of the form ∝ a with a = 0.5988 (4/7 ∼ 0.5714) with standard deviation σ a = 0.0053 and correlation coefficient r = 0.9998.
In Fig. 5, the NLS solution, the semiclassical solution, and the P 2 I solution (3.32) can be seen at the critical time close to the critical point of the semiclassical solution.
The corresponding Riemann invariants can be seen in Fig. 6. For smaller , the agreement of NLS and semiclassical solution becomes better. We show the Riemann invariants r ± for = 10 −3 in Fig. 6. Note that there are also oscillations in the invariant r − which stays smooth at this point in the semiclassical limit.
The P 2 I solution (3.32) gives a much better agreement with the NLS solution close to the critical point as can be seen in Figs. 5 and 6. The agreement is in fact so good that the difference of the solutions has to be studied. The P 2 I solution only gives locally an asymptotic description, at larger distances from the critical point the semiclassical solution provides a better description as can be also seen from Fig. 7.
We can identify the regions where each of the asymptotic solutions gives a better description of NLS than the other by identifying the values x l , x r such that for all x l < x < x r the P 2 I solution provides a better asymptotic description than the semiclassical solution. Due to the oscillatory character of the NLS and the P 2 I solution (3.32), such a definition leads to ambiguities and oscillations also in the boundaries of these zones for r ± . No clear scaling could thus be identified for these limits. The oscillatory character of the solution also implies there is no obvious scaling of the maximal error in the asymptotic description for the values of we could treat.
The matching procedure nonetheless clearly improves the asymptotic description near the critical point. In Fig. 8, we see the difference between this matched asymptotic solution and the NLS solution for two values of . Visibly the zone, where the solutions are matched, decreases with (note the rescaling of the x axes with a factor 6/7 ). The same procedure can be carried out for the invariant r + which stays smooth at this point. Obviously, the P 2 I solution (3.32) provides a description of higher order at this point as can be seen in Fig. 9. Thus, the P 2 I solution (3.32) provides as expected an Modulus of the difference af the Riemann invariants for the defocusing quintic NLS equation for the initial data ψ 0 (x) = sech x for = 0.01 at the critical time t 0 and the semiclassical solution in blue, and the difference between the corresponding P 2 I solution (3.32) and the NLS solution in green; on the left the invariant that has a break-up in the semiclassical limit, and on the right the invariant that stays smooth (Color figure online) asymptotic description of the oscillations for the Riemann invariant which remains smooth in the semiclassical limit.
The P 2 I solution (3.32) holds for small |x − x c | and |t − t 0 |. To illustrate the latter effect, we compare it with the NLS solution for the times t ± = t 0 ± 0.0027. Note that t − t 0 appears in the formula (3.32) for the P 2 I solution at several places with different In the upper part of the left figure, one can see the modulus of the difference + of the Riemann invariant for the defocusing quintic NLS equation for the initial data ψ 0 (x) = sech x at the critical time t 0 and the semiclassical solution for = 0.01. The lower part shows the same difference, which is replaced close to the critical point by the difference between NLS solution and the P 2 I solution (3.32) (in red where the error is smaller than the one shown above). The right figure shows the same situation as the lower figure on the left for = 0.01 above and = 0.001 below. The x axes are rescaled by a factor 6/7 (Color figure online) powers of . Thus in contrast to the elliptic case (5.24), there is no simple dependence on t in the hyperbolic case. In Fig. 10, we show the quantities r ± at the time t ± . It can be seen that the P 2 I solution gives again a clearly better asymptotic description near the break-up point than the semiclassical solution.

"Dark" Initial Data for the Defocusing Quintic NLS
It is well known that the defocusing cubic NLS equation has exact solutions called dark solitons, i.e. solutions that do not tend to zero for |x| → ∞. Such solutions are physically problematic since they have infinite energy and are mathematically difficult to handle, but they are nonetheless of importance in applications. Therefore, we will here also study initial data which do not decay to zero at spatial infinity. We will consider the example ψ 0 (x) = tanh 2 x in the following. The time evolution of the solution up to the critical time t 0 ∼ 1.3448 can be seen in Fig. 11. The steepening of the two fronts of the pulse can be seen as well as the formation of a small oscillation on each side. For times t t 0 , each of the initial oscillations develops into an oscillatory zone which will eventually overlap.
Clearly, there will be two regions with strong gradients symmetric in x. We will concentrate on positive values of x where the Rieman invariant r − breaks in the semiclassical solution. In Fig. 12, the Riemann invariants for the NLS solution, the corresponding semiclassical solution, and the P 2 I asymptotics (3.32) can be seen close to x c ∼ 0.5476 for = 0.001.

Defocusing Non-local NLS
We will study the small dispersion limit of the non-local NLS (5.7) close to the breakup of the corresponding semiclassical solutions. We will concentrate on values of 1 for all studied values of . For both cases, we will consider the initial data ψ 0 = sech x. In the defocusing variant of the non-local NLS equation (5.7), the non-locality has the effect to reduce the defocusing effect of the equation. The dispersion and the steepening of the gradient close to the break-up of the corresponding semiclassical solution are reduced as can be seen in Fig. 13. This also suppresses the formation of dispersive shocks, i.e. the oscillations close to the gradient catastrophe of the semiclassical solution (see Ghofraniha et al. 2007). Due to the possible sign change of the quantity ρ in (3.27), an other effect can be observed in Fig. 13: for large enough η, the oscillations appear on the other side of the critical point. We again consider the initial data ψ 0 = sech x at the critical time t 0 ∼ 1.5244 near the break-up of the Riemann invariant r − at x c ∼ −2.2094 in the semiclassical limit.
For larger times, this implies for ρ < 0 that there is just one oscillation to the right of −x c as described asymptotically by the P 2 I solution, and many small oscillations on the other side of the critical point as can be seen in Fig. 14. The situation is similar to the one of certain Kawahara solutions in the small dispersion limit as discussed in Dubrovin et al. (2011). In the case ρ = 0 in (3.27), the P 2 I asymptotics cannot be used. In the present example, this is the case for η ∼ 1.3060. The solution at the critical time for this value of η can be seen in Fig. 15.
For smaller η, the non-local NLS behaves qualitatively like the defocusing cubic NLS close to the critical time as can be seen in Fig. 16 for the Riemann invariant breaking in the semiclassical limit. For smaller values of , the same behaviour can be seen, but on smaller scales. Again there are two different scales in the P 2 I asymptotics (3.32) which means there is no clear scaling in the coordinates x and t. For the representation, we nonetheless rescale x by a factor of 6/7 to be able to compare the case = 0.001 with = 0.01. The y axes are rescaled to optimally use the space of the figure. The approximation visibly gets better with smaller . The Riemann invariant staying smooth in the semiclassical limit can be seen for the same situation in the right part of Fig. 16. The asymptotic description again improves clearly with smaller .
For larger η, the smoothing out of the gradients near the shock of the semiclassical equations implies that the semiclassical solution only provides a valid asymptotic description for larger |x − x c | than is the case for smaller η. The P 2 I asymptotics (3.32) catches this behaviour as can be seen for η = 100 in Fig. 17 on the left for the invariant breaking in the semiclassical limit. There are essentially no oscillations in this case.
The invariant r + can be seen on the right part of Fig. 17. There is essentially only one oscillation to the right of the critical point in this case. The P 2 I asymptotics has an oscillation close to the oscillation of the non-local NLS and thus catches this behaviour in an asymptotic sense.

Numerical Study of Focusing Generalized and Non-local NLS Equations
In this section, we will study numerically solutions to the focusing NLS before and close to the break-up of the corresponding semiclassical solutions. Since the case of the focusing cubic NLS was studied in detail in Dubrovin et al. (2009), we concentrate here on the not integrable quintic NLS. We compare solutions to NLS and semiclassical equations and for t ∼ t 0 to an asymptotic solution in terms of the tritronquée solution of the P I equation. The same is done for a non-local variant of the cubic NLS equation.  Fig. 18. The focusing effect can be clearly recognized.
For times much smaller than the critical time, one finds that the difference between semiclassical and NLS solution scales as 2 . For instance for t = t 0 /2 t 0 , we obtain for = |u N L S − u sc | via a linear regression analysis for the logarithm of a scaling of the form ∝ a with a = 1.985 with standard deviation σ a = 0.0018 and correlation coefficient r = 0.999998.
At the critical time, the difference between the semiclassical solution and the solution to the focusing quintic NLS scales roughly as 2/5 . More precisely we find via a linear regression analysis for the logarithm of the difference between NLS and semiclassical solution a scaling of the form ∝ a with a = 0.403 with standard deviation σ a = 0.001 and correlation coefficient r = 0.99998. As can be seen in Fig. 19, the semiclassical solution has a cusp. Thus, the maximal difference between semiclassical and NLS solution is always observed for the critical point.
For smaller , the agreement of NLS and semiclassical solution becomes better, but the biggest difference is always at the critical point as can be seen in the bottom of Fig. 19.
The P I solution (5.24) gives a much better agreement with the NLS solution close to the critical point as can be seen in Fig. 19. The agreement is in fact so good that the difference of the solutions has to be studied. The P I solution only gives locally Modulus of the difference of the solution to the focusing quintic NLS equation for the initial data ψ 0 (x) = sech x for = 0.1 at the critical time t 0 and the difference between the corresponding P I solution (5.24) for several values of ; on the left the difference for u, and on the right the difference v for v. The x axes are rescaled with a factor 4/5 an asymptotic description, at larger distances from the critical point the semiclassical solution provides a better description as can be also seen from Fig. 20.
We can identify the regions where each of the asymptotic solutions gives a better description of NLS than the other by identifying the value of x r such that for all x > x r the semiclassical solutions give a better asymptotic description than the multiscales solution (since the solution is symmetric with respect to x, we only consider positive values of x here). We find that the width of this zone scales roughly as 3/5 . A linear regression analysis for the dependence of log 10 x r on log 10 yields a = 0.634 with standard deviation σ a = 0.0036 and correlation coefficient r = 0.99993.
This matching procedure clearly improves the NLS description near the critical point. In Fig. 21, we see the difference between this matched asymptotic solution and the NLS solution for two values of . Visibly the zone, where the solutions are matched, decreases with (note the rescaling of the x axes by a factor 4/5 ).
A linear regression analysis for the logarithm of the difference between NLS and multiscales solution in the matching zone gives a scaling of the form ∝ a with a = 0.6659 with standard deviation σ a = 0.032 and correlation coefficient r = 0.995. The found scaling is thus in the whole interval clearly better than the 2/5 of the semiclassical solution, but does not reach the expected 4/5 scaling in the whole interval. This indicates that transition formulae between the multiscales and the semiclassical solution have to be established as in Grava and Klein (2012) for KdV, which is, however, beyond the scope of the present paper.
The P I solution (5.24) holds for small |x − x c | and |t − t 0 |. To illustrate the latter effect, we compare it with the NLS solution for the times t ± ( ) = t 0 ± 0.01 4/5 where we take care of the scaling of t in (4.26). In Fig. 22, we show the quantity for 2 values of at the times t − ( ). The x axes are rescaled by a factor 4/5 . It can be seen that the quality of the asymptotic description is slightly lower than at the critical time, but that the error is of a similar order. The situation is similar at the time t + = t 0 + 0.01 4/5 as can be seen also in Fig. 22.

Non-Symmetric Initial Data for the Focusing Quintic NLS
To study solutions to the focusing quintic NLS for the asymmetric initial data (6.27), we first have to solve equations (6.27) numerically for α = 0.2. This is done for values of |x| < 15 in a standard way by solving (6.30) on some Chebyshev collocation points with a Newton iteration. The choice of this interval is determined by the fact that the residual of the Newton iterate is smaller than 10 −10 on the whole interval. We choose N c = 512 collocation points to ensure that the coefficients of an expansion of the solution decrease to machine precision and that the solution is thus numerically fully resolved. For values of |x| > 15, we solve (6.27) asymptotically, for x → +∞ and r = 1+i exp(x)2 1+2α +2 2+4α exp(2x) −0.5 + 2α 2 ln(2) + αx − α +O(exp(3x)) (9.2) for x → −∞. Machine precision is reached for |x| > 15 for this asymptotic solution. Initial data for α = 0.2 can be seen in Fig. 23.
To obtain initial data for the NLS equation from r = v + iu in the form ψ = √ u exp(i x x 0 v(x )dx / ), we have to integrate the real part of r with respect to x. This is done by using an expansion of the solution for |x| < 15 in terms of Chebyshev polynomials via a discrete cosine transform (this is the reason why the solution was computed on Chebyshev collocation points) and applying the well-known formula for the integral of Chebyshev polynomials. For values of |x| > 15, the asymptotic formulae (9.1) and (9.2) are integrated analytically by choosing the integration constants to obtain a continuous matching with the numerically integrated v. This way we obtain initial data with an accuracy of better than 10 −10 . We put the Krasny filter to the order of this threshold and thus obtain initial data resolved up to the level of the Krasny filter.
For = 0.1, the solution to the focusing quintic NLS equation for the asymmetric initial data as well as the semiclassical and the P I asymptotics (5.24) can be seen in Fig. 24. As expected, the P I asymptotics gives a much better description of the NLS solution close to the critical point of the semiclassical solution. The error in the approximation is, however, also not symmetric here.
The agreement gets even better for smaller . We can reach values as low as = 0.02. For smaller , the blow-up singularity of quintic NLS solutions seems to be too close to the critical time of the semiclassical solution which breaks the code. The case = 0.02 is, however, numerically fully resolved. As can be seen in the lower row of Fig. 24, the agreement is as expected. Note that also in this case the x axes of the bottom figures have been rescaled by a factor 4/5 .

"Dark" Initial Data
Focusing NLS equations do not have dark solitons as exact solutions, i.e. solutions which tend asymptotically to a nonzero constant and which vanish for finite values of x. But it is mathematically interesting to study how initial data of this form lead to a breakup of the semiclassical equations, and how the corresponding NLS solution behaves in the vicinity of the critical point. We consider here initial data of the form ψ 0 = tanh 2 x. The solution breaks here in the form of two cusps symmetric with respect to x = 0. The critical time is at t 0 = 0.9041 . . ., the cusps form at x c = ±1.8723 . . .. The corresponding solution can be seen in Fig. 25. For = 0.1, the solution to the focusing quintic NLS equation for the dark initial data as well as the semiclassical and the P I asymptotics (5.24) can be seen in Fig. 26. As expected, the P I asymptotics gives a much better description of the NLS solution close to the critical point of the semiclassical solution. The agreement gets better for smaller . We can reach values as low as = 0.04, where the modulation instability leads to problems for smaller values of because of the asymptotically non-vanishing solution. The case = 0.04 is, however, numerically accessible. As can be seen in the bottom figures of Fig. 26, the agreement is as expected.

Blow-Up
For the cubic focusing NLS, solutions in the semiclassical limit for times t t 0 develop a zone of rapid modulated oscillations as can be seen for instance in Fig. 27. The central hump close to the critical time splits into several humps of smaller amplitude. For the quintic NLS on the other hand, it is known, see e.g., Merle and Raphael (2004), that initial data with negative energy have a blow-up in finite time. For the NLS with the semiclassical parameter we consider in this paper, this will be always the case for sufficiently small . Thus, the solution of the quintic NLS looks for small very differently from the solution to the cubic NLS for the same initial data and the same value of as can be seen in Fig. 27. The central hump develops in this case into a blow-up. For obvious reasons, it is impossible to treat a blow-up exactly numerically, but the numerical solution can get sufficiently close to this case. Driscoll's composite Runge-Kutta method produces an overflow error close to the L ∞ blow-up encountered here because of the term |ψ| 4 ψ. We stop the code when this happens and note the last time with finite value of ψ as a lower bound t B for the blow-up time. The error in the determination of the blow-up time with this method is largest for larger . Using linear regression, we find for ln(t B − t 0 ) = a ln + b for values of = 0.01, 0.02, . . . , 0.1 the value a = 0.83 close to 4/5 with standard deviation σ a = 0.0439, b = −0.1267 with standard deviation σ b = 0.0138 correlation coefficient r = 0.999, see Fig. 28.
As expected from the P I solution (5.24), the time scales with 4/5 . Since we expect the error in the determination of the blow-up time to decrease with , a slightly stronger decrease with of the time t B than predicted is no surprise.
It is an interesting question whether the blow-up time in the limit → 0 is related to the first pole of the tritronquée solution on the negative real axis. In Joshi and Kitaev (2001), it was shown that the first pole is located at ξ pole = −2.3841687 . . . . x one can see that for quintic NLS and sechx initial data, the point of elliptic umbilic catastrophe is at x 0 = 0, and for symmetry reasons, the blow-up is at x B = 0. Using the above formula, with V (u) = u 2 2 so that V 0 = u 0 , V 0 = 1 and a + = − i re iψ = −i with u 0 = 1.5858 determined in (6.24) for this specific example, the blowup time t B is then conjectured to satisfy the equation which gives a value of |b| = ln(2.3841/2.0324) = 0.1596, in reasonable agreement with the numerically found value |b| ∼ 0.1267.

Focusing Non-local NLS
We will study the small dispersion limit of the non-local NLS (5.7) close to the breakup of the corresponding semiclassical solutions. We will concentrate on values of η such that η 2 1 for all studied values of . For both cases, we will consider the initial data ψ 0 = sech x. The effect of the non-locality in (5.7) is to reduce the focusing effect of the focusing NLS. This means the larger η, the smaller the value for the maximum at the critical time of the corresponding semiclassical solution, and the less pronounced the focusing of the maximum, i.e. smaller gradients in the solution. This effect can be clearly seen in Fig. 29. At the critical time, the tritronquée solution to P I gives as expected a much better description of the non-local NLS solution than the semiclassical solution as can be seen for η = 0.1 for u in Fig. 31. The quality of the approximation increases visibly for smaller . Note that the x axes are rescaled with a factor 4/5 .
The corresponding plots for v can be seen in Fig. 31. The same behaviour as for u is visible (Fig. 32).
For larger values of η, the agreement is less good for both the semiclassical and the P I asymptotics. This is clear for the former since the semiclassical solution is independent of η, and since the focusing effect of the non-local NLS is less pronounced for larger values of η. The P I asymptotics takes this into account, the value of its maximum is also reduced, but more so than for the non-local NLS which implies that the agreement between the two solutions is best for η = 0, i.e. the cubic NLS. The approximation gets, however, better for smaller as can be seen for η = 1 in Fig. 33.
The corresponding plots for v can be seen in Fig. 34.