On universality of critical behaviour in the focusing nonlinear Schr\"odinger equation, elliptic umbilic catastrophe and the {\it tritronqu\'ee} solution to the Painlev\'e-I equation

TL;DR

This paper demonstrates that the critical behavior near gradient catastrophe in the focusing nonlinear Schrödinger equation is approximately described by a special solution to the Painlevé-I equation, revealing universality in such nonlinear wave phenomena.

Contribution

It establishes a connection between the critical behavior of the focusing nonlinear Schrödinger equation and the tritronquée solution of the Painlevé-I equation, highlighting a universal pattern.

Findings

Critical behavior is described by Painlevé-I tritronquée solution.

Universality of critical phenomena in nonlinear Schrödinger equations.

Analytic initial data lead to Painlevé-I based asymptotics.

Abstract

We argue that the critical behaviour near the point of ``gradient catastrophe" of the solution to the Cauchy problem for the focusing nonlinear Schr\"odinger equation $ i\epsilon \psi_t +\frac{\epsilon^2}2\psi_{xx}+ |\psi|^2 \psi =0$ with analytic initial data of the form $\psi(x,0;\epsilon) =A(x) e^{\frac{i}{\epsilon} S(x)}$ is approximately described by a particular solution to the Painlev\'e-I equation.

Figures

Equations

$$i{\Psi }_t+\frac{1}{2}{\Psi }_{xx}+|\Psi {|}^2\Psi =0$$

$${\Psi }_t+\{\Psi (x),\mathcal{H}\}=0$$

$$\{\Psi (x),{\Psi }^{\ast }(y)\}=i\delta (x-y)$$

$$\mathcal{H}=\frac{1}{2}\int \left(|{\Psi }_x{|}^2-|\Psi {|}^4\right)dx$$

$$i{\Psi }_t+\frac{1}{2}{\Psi }_{xx}-|\Psi {|}^2\Psi =0$$

$$\Psi =A{e}^{i(kx-\omega t)},\omega =\frac{1}{2}{k}^2-A^2$$

$$\Psi (x,0;ϵ)=A(x)e^{\frac{i}{ϵ}S(x)}$$

$$iϵ{\Psi }_t+\frac{{ϵ}^2}{2}{\Psi }_{xx}+|\Psi {|}^2\Psi =0.$$

$$u=|\Psi {|}^2,v=\frac{ϵ}{2i}\left(\frac{{\Psi }_x}{\Psi }-\frac{{\Psi }_x^{\ast }}{{\Psi }^{\ast }}\right)$$

$$u_t+(uv{)}_x=0$$

$$v_t+v{v}_x-u_x+\frac{{ϵ}^2}{4}{\left(\frac{1}{2}\frac{u_x^2}{u^2}-\frac{u_{xx}}{u}\right)}_x=0.$$

$$u(x,0)=A^2(x),v(x,0)=S^{\prime }(x).$$

$$\left.\begin{array}[]{lr}u_{t}&+v\,u_{x}+u\,v_{x}=0\\ &\\ v_{t}&-~{}~{}u_{x}+v\,v_{x}=0\end{array}\right\}.$$

$$\left(\begin{array}[]{cc}v&u\\ -1&v\end{array}\right)$$

$${\Omega }_{\zeta \zeta }=6{\Omega }^2-\zeta .$$

$$\left.\begin{array}[]{r}u_{t}+\{u(x),H\}=0\\ \\ v_{t}+\{v(x),H\}=0\end{array}\right\}$$

$$\{u(x),v(y)\}={\delta }^{\prime }(x-y),$$

$$H=\int \left[\frac{1}{2}\left(u{v}^2-u^2\right)+\frac{{ϵ}^2}{8u}{u}_x^2\right]dx.$$

$${\left(u_{I,x}\right)}^2+{\left(v_{I,x}\right)}^2\ne 0.$$

$$u(x,0)=u_I(x),v(x,0)=v_I(x)$$

$$\left.\begin{array}[]{l}x=v\,t+f_{u}\\ \\ 0=u\,t+f_{v}\end{array}\right\}$$

$$f_{vv}+u{f}_{uu}=0.$$

$$F=\int f(u,v)dx,\frac{d}{dt}F=0.$$

$$\left.\begin{array}[]{r}u_{s}+\{u(x),F\}\equiv u_{t}+\left(f_{v}\right)_{x}=0\\ \\ v_{s}+\{v(x),F\}\equiv v_{t}+\left(f_{u}\right)_{x}=0\end{array}\right\}$$

$${\left(u_t\right)}_s={\left(u_s\right)}_t,{\left(v_t\right)}_s={\left(v_s\right)}_t.$$

$$g=-\frac{1}{2}{v}^2+u(\log u-1).$$

$$\left.\begin{array}[]{r}u_{s}+v_{x}=0\\ \\ v_{s}=\frac{u_{x}}{u}\end{array}\right\}$$

$$u_{ss}+{\left(\log u\right)}_{xx}=0.$$

$$\left.\begin{array}[]{l}x=v\,t+f_{u}\\ \\ s=u\,t+f_{v}\end{array}\right\}$$

$$\begin{array}[]{l}x_{0}=v_{0}t_{0}+f_{u}(u_{0},v_{0})\\ \\ s_{0}=u_{0}t_{0}+f_{v}(u_{0},v_{0})\\ \\ f_{uu}(u_{0},v_{0})=f_{vv}(u_{0},v_{0})=0,\quad f_{uv}(u_{0},v_{0})=-t_{0}\end{array}$$

$$f_{uu}^0=f_{vv}^0=0,f_{uv}^0=-t_0.$$

$$f_{uuv}^0\ne 0.$$

Full text

1 Introduction

The focusing nonlinear Schrödinger (NLS) equation for the complex valued function $\Psi=\Psi(x,t)$

[TABLE]

has numerous physical applications in the description of nonlinear waves (see, e.g., the books [47, 35, 36]). It can be considered as an infinite dimensional analogue of a completely integrable Hamiltonian system [48], where the Hamiltonian and the Poisson bracket is given by

[TABLE]

(here $\Psi^{*}$ stands for the complex conjugate function). Properties of various classes of solutions to this equation have been extensively studied both analytically and numerically [5, 6, 7, 8, 16, 21, 25, 26, 29, 32, 34, 43, 44]. One of the striking features that distinguishes this equation from, say, the defocusing case

[TABLE]

is the phenomenon of modulation instability [1, 11, 16]. Namely, slow modulations of the plane wave solutions

[TABLE]

develop fast oscillations in finite time.

The appropriate mathematical framework for studying these phenomena is the theory of the initial value problem

[TABLE]

for the $\epsilon$-dependent NLS

[TABLE]

Here $\epsilon>0$ is a small parameter, $A(x)$ and $S(x)$ are real-valued smooth functions. Introducing the slow variables

[TABLE]

the equation can be recast into the following system:

[TABLE]

The initial data for the system (1) coming from (1.3) do not depend on $\epsilon$:

[TABLE]

The simplest explanation of the modulation instability then comes from considering the so-called dispersionless limit $\epsilon\to 0$. In this limit one obtains the following first order quasilinear system

[TABLE]

This is a system of elliptic type because of the condition $u>0$. Indeed, the eigenvalues of the coefficient matrix

[TABLE]

are complex conjugate, $\lambda=v\pm i\sqrt{u}$. So, the Cauchy problem for the system (1.8) is ill-posed in the Hadamard sense (cf. [33, 7]). Even for analytic initial data the life span of a typical solution is finite, $t<t_{0}$. The $x$- and $t$-derivatives explode at some point $x=x_{0}$ when the time approaches $t_{0}$. This phenomenon is similar to the gradient catastrophe of solutions to nonlinear hyperbolic PDEs [2].

For the full system (1) the Cauchy problem is well-posed for a suitable class of $\epsilon$-independent initial data (see details in [17, 46]). However, the well-posedness is not uniform in $\epsilon$. In practical terms that means that the solution to (1) behaves in a very irregular way in some regions of the $(x,t)$-plane when $\epsilon\to 0$. Such an irregular behaviour begins near the points $(x=x_{0},t=t_{0})$ of the “gradient catastrophe” of the solution to the dispersionless limit (1.8). The solutions to (1.8) and (1) are essentially indistinguishable for $t<t_{0}$; the situation changes dramatically near $x_{0}$ when approaching the critical point. Namely, when approaching $t=t_{0}$ the peak near a local maximum111Regarding initial data with local minima we did not observe cusps related to minima in numerical simulations. We believe they do not exist because of the focusing effect in the NLS that pushes maxima to cusps but seems to smoothen minima. of $u$ becomes more and more narrow due to self-focusing; the solution develops a zone of rapid oscillations for $t>t_{0}$. They have been studied both analytically and numerically in [8, 16, 21, 23, 25, 26, 34, 43, 44]. However, no results are available so far about the behaviour of the solutions to the focusing NLS at the critical point $(x_{0},t_{0})$.

The main subject of this work is the study of the behaviour of solutions to the Cauchy problem (1), (1.7) near the point of gradient catastrophe of the dispersionless system (1.8). In order to deal with the Cauchy problem for (1.8) we will assume analyticity222We believe that the main conclusions of this paper must hold true also for non analytic initial data; the numerical experiments of [8] do not show much difference in the properties of solution between analytic and non analytic cases. However, the precise formulation of our Main Conjecture has to be refined in the non analytic case. of the initial data $u(x,0)$, $v(x,0)$. Then the Cauchy problem for (1.8) can be solved for $t<t_{0}$ via a suitable version of the hodograph transform (see Section 2 below). An important feature of the gradient catastrophe for this system is that it happens at an isolated point of the $(x,t)$-plane, unlikely the case of KdV or defocusing NLS where the singularity of the hodograph solution takes place on a curve. We identify this singularity for a generic solution to (1.8) as the elliptic umbilic singularity (see Section 4 below) in the terminology of R.Thom [42]. This codimension 2 singularity is one of the real forms labeled by the root system of the $D_{4}$ type in the terminology of V.Arnold et al. [3].

Our main goal is to find a replacement for the elliptic umbilic singularity when the dispersive terms are added, i.e., we want to describe the leading term of the asymptotic behaviour for $\epsilon\to 0$ of the solution to (1) near the critical point $(x_{0},t_{0})$ of a generic solution to (1.8).

Thus, our study can be considered as a continuation of the programme initiated in [13] to study critical behaviour of Hamiltonian perturbations of nonlinear hyperbolic PDEs; the fundamental difference is that the non perturbed system (1.8) is not hyperbolic! However, many ideas and methods of [13] (see also [12]) play an important role in our considerations.

The most important of these is the idea of universality of the critical behaviour. The general formulation of the universality suggested in [13] for the case of Hamiltonian perturbations of the scalar nonlinear transport 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 in [13] via a particular solution to the fourth order analogue of the Painlevé-I equation (the so-called $P_{1}^{2}$ equation). The existence of the needed solution to $P_{1}^{2}$ has been rigorously established in [9]. Moreover, it was argued in [13] that this behaviour is essentially independent on the choice of the Hamiltonian perturbation. Some of these universality conjectures have been partially confirmed by numerical analysis carried out in [21].

The main message of this paper is the formulation of the Universality Conjecture for the critical behaviour of the solutions to the focusing NLS. Our considerations suggest the description of the leading term in the asymptotic expansion of the solution to (1) near the critical point via a particular solution to the classical Painlevé-I equation (P-I)

[TABLE]

The so-called tritronquée solution to P-I was discovered by P.Boutroux [4] as the unique solution having no poles in the sector $|\arg\zeta|<4\pi/5$ for sufficiently large $|\zeta|$. Remarkably, the very same solution333It is interesting that the same tritronquée solution (for real $\zeta$ only) appeares also in the study of certain critical phenomena in plasma [41]. In the theory of random matrices and orthogonal polynomials a different solution to P-I arises; see, e.g., [14]. arises in the critical behaviour of solutions to focusing NLS!

The paper is organized as follows. In Section 2 we develop a version of the hodograph transform for integrating the dispersionless system (1.8) before the catastrophe $t<t_{0}$. We also establish the shape of the singularity of the solution near the critical point; the latter is identified in Section 4 with the elliptic umbilic catastrophe of Thom. In Section 3 we develop a method of constructing formal perturbative solutions to the full system (1) before the critical time. In Section 5 we collect the necessary information about the tritronquée solution of P-I and formulate the Main Conjecture of this paper. Such a formulation relies on a much stronger property of the tritronquée solution: namely, we need this solution to be pole-free in the whole sector $|\arg\zeta|<4\pi/5$. Numerical evidence for the absence of poles in this sector is given in Section 6. In Section 7 we analyze numerically the agreement between the critical behaviour of solutions to focusing NLS and its conjectural description in terms of the tritronquée solution restricted on certain lines in the complex $\zeta$-plane. In the final Section 8 we give some additional remark and outline the programme of future research.

Acknowledgments. The authors thank K.McLaughlin for a very instructive discussion. One of the authors (B.D.) thanks R.Conte for bringing his attention to the tritronquées solutions of P-I. The results of this paper have been presented by one of the authors (T.G.) at the Conference “The Theory of Highly Oscillatory Problems”, Newton Institute, Cambridge, March 26 - 30, 2007. T.G. thanks A.Fokas and S.Venakides for the stimulating discussion after the talk. The present work is partially supported by the European Science Foundation Programme “Methods of Integrable Systems, Geometry, Applied Mathematics” (MISGAM), Marie Curie RTN “European Network in Geometry, Mathematical Physics and Applications” (ENIGMA). The work of B.D. and T.G. is also partially supported by Italian Ministry of Universities and Researches (MUR) research grant PRIN 2006 “Geometric methods in the theory of nonlinear waves and their applications”.

2 Dispersionless NLS, its solutions and critical behaviour

The equations (1) are a Hamiltonian system

[TABLE]

with respect to the Poisson bracket originated in (1)

[TABLE]

other brackets vanish, with the Hamiltonian

[TABLE]

Let us first describe the general analytic solution to the dispersionless system (1.8).

Lemma 2.1

Let $u_{I}(x)$, $v_{I}(x)$ be two real valued analytic functions of the real variable $x$ satisfying

[TABLE]

Then the solution $u=u(x,t)$, $v=v(x,t)$ to the Cauchy problem

[TABLE]

for the system (1.8) for sufficiently small $t$ can be determined from the following system

[TABLE]

where $f=f(u,v)$ is an analytic solution to the following linear PDE:

[TABLE]

Conversely, given any solution to (2.5) satisfying $u^{2}f_{uu}^{2}+f_{uv}^{2}\neq 0$ at some point $(u=\tilde{u},v=\tilde{v})$ such that $f_{v}(\tilde{u},\tilde{v})=0$, the system (2.4) determines a solution to (1.8) defined locally near the point $x=\tilde{x}:=f_{u}(\tilde{u},\tilde{v})$ for sufficiently small $t$.

Remark 2.2

The solutions to the linear PDE (2.5) correspond to the first integrals of dispersionless NLS:

[TABLE]

Taking them as the Hamiltonians

[TABLE]

yields infinitesimal symmetries of the dispersionless NLS:

[TABLE]

One of the first integrals will be extensively used in this paper. It corresponds to the Hamiltonian density

[TABLE]

The associated Hamiltonian flow reads

[TABLE]

Eliminating the dependent variable $v$ one arrives at the elliptic version of the long wave limit of Toda lattice:

[TABLE]

Due to commutativity (2.8) the systems (1.8) and (2.10) admit a simultaneous solution $u=u(x,t,s)$, $v=v(x,t,s)$. Any such solution can be locally determined from a system similar to (2.4)

[TABLE]

where $f=f(u,v)$, as above, solves the linear PDE (2.5).

The system (2.11) determines a solution $u=u(x,t,s)$, $v=v(x,t,s)$ provided applicability of the implicit function theorem. The conditions of the latter fail to hold at the critical point $(x_{0},s_{0},t_{0},u_{0},v_{0})$ such that

[TABLE]

In sequel we adopt the following system of notations: the values of the function $f$ and of its derivatives at the critical point will be denoted by $f^{0}$ etc. E.g., the last line of the conditions (2.12) will read

[TABLE]

Definition 1. We say that the critical point is generic if at this point:

[TABLE]

Let us the introduce real parameters $r$, $\psi$ determined by the third derivatives of the function $f$ evaluated at the critical point,

[TABLE]

Due to the genericity assumption

[TABLE]

In order to describe the local behaviour of a solution to the dispersionless NLS/Toda equations we define a function $R(X;S,\psi)$ of real variables $X$, $S$ depending on the real parameter $\psi$ satisfying

[TABLE]

by the following formula

[TABLE]

Put

[TABLE]

Observe that $P_{0}(X,S,\psi)$ and $Q_{0}(X,S,\psi)$ are smooth functions of the real variable $X$ provided validity of the inequality (2.14).

Lemma 2.3

Given an analytic solution $u(x,s,t)$, $v(x,s,t)$ to the dispersionless NLS/Toda equations with a generic critical point $(x_{0},s_{0},t_{0},u_{0},v_{0})$, and arbitrary real numbers $X$, $S$ satisfying (2.14), $T<0$, then there exist the following limits

[TABLE]

where the parameters $r$, $\psi$ are defined by (2.13).

*Proof  *From the linear PDE (2.5) it follows that

[TABLE]

Using these formulae we expand the implicit function equations (2.11) near the critical point in the form

[TABLE]

where we introduce the shifted variables

[TABLE]

The rescaling

[TABLE]

in the limit $\lambda\to 0$ yields the quadratic equation

[TABLE]

where the complex independent and dependent variables $z$ and $w$ read

[TABLE]

and the complex constant $a$ is defined by

[TABLE]

therefore

[TABLE]

The substitution

[TABLE]

reduces the quadratic equation to

[TABLE]

For $\bar{t}<0$ we choose the following root

[TABLE]

where the branch of the square root is obtained by the analytic continuation of the one taking positive values on the positive real axis. Equivalently,

[TABLE]

where

[TABLE]

This gives the formulae (2.15).

The result of the lemma describes the local structure of generic solutions to the dispersionless NLS/Toda equations near the critical point. It can also be represented in the following form

[TABLE]

where

[TABLE]

We want to emphasize that the approximation (2) works only near the critical point. Indeed, for large $x\to\pm\infty$ the function $u(x,s,t)$ and $v(x,s,t)$ have the following behaviour

[TABLE]

So, for sufficiently large $|x|$ the function $u(x,s,t)$ defined by (2) becomes negative.

The function $u$ has a maximum at the point $X=S\,\tan{\psi}$, so locally

[TABLE]

At the critical point $(x_{0},s_{0},t_{0},u_{0},v_{0})$ the function $u$ develops a cusp. Let us consider only the particular case $S=0$ in order to avoid complicated expressions. In this case the local behaviour of the function $u$ near the critical point is given by

[TABLE]

(here $\hat{x}=\sqrt{u}_{0}(\bar{x}-v_{0}\bar{t})$). Thus the parameters $r$, $\psi$ describe the shape of the cusp at the critical point.

3 First integrals and solutions of the NLS/Toda equations

Let us first show that any first integral (2.6) of the dispersionless equations can be uniquely extended to a first integral of the full equations.

Lemma 3.1

Given a solution $f=f(u,v)$ to the linear PDE (2.5), there exists a unique, up to a total derivative, formal power series in $\epsilon$

[TABLE]

such that the integral

[TABLE]

commutes with the Hamiltonian of NLS equation:

[TABLE]

at every order in $\epsilon$. Explicitly,

[TABLE]

Here we use short notations

[TABLE]

Example 1. Taking $f=\frac{1}{2}(u\,v^{2}-u^{2})$ one obtains the Hamiltonian of the NLS equation

[TABLE]

In this case the infinite series truncates. It is easy to see that the series in $\epsilon$ truncates if and only if $f(u,v)$ is a polynomial in $u$. Polynomial in $u$ solutions to the linear PDE (2.5) correspond to the standard first integrals of the NLS hierarchy.

Example 2. Taking $g=-\frac{1}{2}v^{2}+u(\log u-1)$ (cf. (2.9)) one obtains the Hamiltonian of Toda equation

[TABLE]

written in terms of the function $\phi=\log u$ in the form

[TABLE]

Lemma 3.2

Any solution to the NLS/Toda equations in the class of formal power series in $\epsilon$ can be obtained from the equations

[TABLE]

where $f=f(u,v;\epsilon)$ is an arbitrary admissible solution to the linear PDE (2.5) in the class of formal power series in $\epsilon$,

[TABLE]

Now, we can apply to the system (3.2) the rescaling (2.19) accompanied by the transformation

[TABLE]

At the limit $\lambda\to 0$ we arrive at the following system of equations

[TABLE]

Using the complex variables $z$, $w$ defined in (2.21) we can rewrite the system in the following form:

[TABLE]

The last observation is that the Toda equations generated by the Hamiltonian $H_{g}=\int h_{g}\,dx$ (see Example 2 above) after the scaling limit (2.19), (3.4) yield the Cauchy - Riemann equations for the function $w=w(z)$,

[TABLE]

Therefore the system (3) can be recast into the form equivalent to the Painlevé-I (P-I) equation (see (5.1) below)

[TABLE]

Choosing

[TABLE]

we eliminate $\epsilon$ from the equation.

In the Section 5 below we will write explicitly the reduction of (3.7) to the Painlevé-I equation and give a conjectural characterization of the particular solution of the latter.

4 Critical behaviour and elliptic umbilic catastrophe

Separating again the real and complex parts of (3.6) one obtains a system of ODEs

[TABLE]

that can be identified with the Euler - Lagrange equations

[TABLE]

with the Lagrangian

[TABLE]

In the “dispersionless limit” $\epsilon\to 0$ the Euler - Lagrange equations reduce to the search of stationary points of a function (let us also set $\bar{t}=0$)

[TABLE]

where we redenote

[TABLE]

At $a_{+}=a_{-}=0$ the function $F$ has an isolated singularity at the origin $U=V=0$ of the type $D_{4,\,-}$ also called elliptic umbilic singularity, according to R.Thom [42]. This singularity appears in various physical problems; we mention here the caustics in the collisionless dark matter [40] to give just an example. The parameters $a_{+}$ and $a_{-}$ define two particular directions on the base of the miniversal unfolding of the elliptic umbilic; the full unfolding depending on 4 parameters reads

[TABLE]

It would be interesting to study the properties of the modified Euler - Lagrange equations for the Lagrangian

[TABLE]

This deformation does not seem to arrive from considering solutions to the NLS hierarchy.

5 The tritronquée solution to the Painlevé-I equation and the Main Conjecture

In this section we will select a particular solution to the Painlevé-I (P-I) equation

[TABLE]

Recall [22] that an arbitrary solution to this equation is a meromorphic function on the complex $\zeta$-plane. According to P. Boutroux [4] the poles of the solutions accumulate along the rays

[TABLE]

Boutroux proved that, for each ray there is a one-parameter family of particular solutions called intégrales tronquées whose lines of poles truncate for large $\zeta$. He proved that the intégrale tronquée has no poles for large $|\zeta|$ within two consecutive sectors of the angle $2\pi/5$ around the ray, and, moreover it has the asymptotic behaviour of the form

[TABLE]

for a suitable choice of the branch of the square root (see below) and a sufficiently small $\varepsilon>0$.

Furthemore, if a solution truncates along any two of the rays (5.2) then it truncates along three of them. These particular solutions to P-I are called tritronquées. They have no poles for large $|\zeta|$ in four consecutive sectors; their asymptotics for large $\zeta$ is given by (5.3). It suffices to know the tritronquée solution $\Omega_{0}(\zeta)$ for the sector

[TABLE]

In this case the branch of the square root in (5.3) is obtained by the analytic continuation of the principal branch taking positive values on the positive half axis $\zeta>0$. Other four tritronquées solutions are obtained by applying the symmetry

[TABLE]

The properties of the tritronquées solutions in the finite part of the complex plane were studied in the important paper of N.Joshi and A.Kitaev [24].

A. Kapaev [27] obtained a complete characterization of the tritronquées solutions in terms of the Riemann - Hilbert problem associated with P-I. We will briefly sketch here the main steps of his construction.

The equation (5.1) can be represented as the compatibility condition of the following system of linear differential equations for a two-component vector valued function $\Phi=\Phi(\lambda,\zeta)$

[TABLE]

The canonical matrix solutions $\Phi_{k}(\lambda,\zeta)$ to the system (5.9) - (5.13) are uniquely determined by their asymptotic behaviour

[TABLE]

in the sectors

[TABLE]

Here

[TABLE]

the branch cut on the complex $\lambda$-plane for the fractional powers of $\lambda$ is chosen along the negative real half-line.

The Stokes matrices $S_{k}$ are defined by

[TABLE]

They have the triangular form

[TABLE]

and satisfy the constraints

[TABLE]

where

[TABLE]

Due to (5.27) two of the Stokes multipliers $s_{k}$ determine all others; they depend neither on $\lambda$ nor on $\zeta$ provided $\Omega(\zeta)$ satisfies (5.1).

In order to obtain a parametrization of solutions to the P-I equation (5.1) by Stokes multipliers of the linear differential equation (5.9) one has to reformulate the above definitions as a certain Riemann - Hilbert problem. The solution of the Riemann - Hilbert problem depends on $\zeta$ through the asymptotics (5.18). If the Riemann - Hilbert problem has a unique solution for the given $\zeta_{0}\in\mathbb{C}$ then the canonical matrices $\Phi_{k}(\lambda,\zeta)$ depend analytically on $\zeta$ for sufficiently small $|\zeta-\zeta_{0}|$; the coefficient $\Omega=\Omega(\zeta)$ will then satisfy (5.1). The poles of the meromorphic function $\Omega(\zeta)$ correspond to the forbidden values of the parameter $\zeta$ for which the Riemann - Hilbert problem admits no solution.

We will now consider a particular solution to the P-I equation specified by the following Riemann - Hilbert problem. Denote four oriented rays $\gamma_{0}$, $\gamma_{\pm 1}$, $\rho$ in the complex $\lambda$-plane defined by

[TABLE]

directed towards infinity. The rays divide the complex plane in four sectors. We are looking for a piecewise analytic function $\Phi(\lambda,\zeta)$ on

[TABLE]

depending on the parameter $\zeta$ continuous up to the boundary with the asymptotic behaviour at $|\lambda|\to\infty$ of the form (5.18) satisfying the following jump conditions on the rays

[TABLE]

Here the plus/minus subscripts refer to the boundary values of $\Phi$ respectively on the left/right sides of the corresponding oriented ray, the jump matrices are given by

[TABLE]

The following result is due to A.Kapaev444Our solution $\Omega_{0}(\zeta)$ coincides with $y_{3}(x)\equiv y_{-2}(x)$, $x=-\zeta$, of [27] (see eq. (2.73) of [27]; Kapaev uses a different normalization $y^{\prime\prime}=6y^{2}+x$ of the P-I equation). .

Theorem 5.1

The solution to the above Riemann - Hilbert problem exists and it is unique for

[TABLE]

for a sufficiently large positive number $R$. The associated function

[TABLE]

is analytic in the domain (5.31), it satisfies P-I and enjoys the asymptotic behaviour

[TABLE]

Moreover, any solution of P-I having no poles in the sector (5.31) for some large $R>0$ coincides with $\Omega_{0}(\zeta)$.

Joshi and Kitaev proved that the tritronquée solution has no poles on the positive real axis. They found a numerical estimate for the position of the first pole $\zeta_{0}$ of the tritronquée solution $\Omega_{0}(\zeta)$ on the negative real axis:

[TABLE]

(cf. also [10]). Besides this estimate very little is known about the location of poles of $\Omega_{0}(\zeta)$. Our numerical experiments (see below) suggest the following

Main Conjecture. Part 1. The tritronquée solution $\Omega_{0}(\zeta)$ has no poles in the sector

[TABLE]

We are now ready to describe the conjectural universal structure behind the critical behaviour of generic solutions to the focusing NLS. For simplicity of the formulation let us assume $\cos\psi>0$.

Main Conjecture. Part 2. Any generic solution to the NLS/Toda equations near the critical point behaves as follows

[TABLE]

*where $\Omega_{0}(\zeta)$ is the tritronquée solution to the Painlevé-I equation (5.1). *

The above considerations can actually be applied replacing the NLS time flow by any other flow of the NLS/Toda hierarchy. The local description of the critical behaviour remains unchanged.

Remark 5.2

Note that the angle of the line $\zeta(\bar{x})$ in (5) for $\bar{t}$ fixed is equal to $\psi/5+\pi/2$, $\psi\in[-\pi,\pi]$. Thus the maximal value of $\mbox{arg}\zeta$ is equal to $7\pi/10<4\pi/5$. The lines in (5) consequently do not get close to the critical lines of the tritronquée solution.

6 Numerical analysis of the tritronquée solution of P-I

In this section we will numerically construct the tritronquée solution $\Omega_{0}$, i.e. the tritronquée solution with asymptotic behavior (5.3). We will drop the index [math] in the following. The solution will be first constructed on a straight line in the complex plane. In a second step we will then explore global properties of these solutions within the limitations imposed by a numerical approach555Cf. [15] where a similar technique was applied to solve numerically the Painlevé-II equation in the complex domain..

Let the straight line in the complex plane be given by $\zeta=ay+b$ with $a,b\in\mathbb{C}$ constant (we choose $a$ to have a non-negative imaginary part) and $y\in\mathbb{R}$. The asymptotic conditions are

[TABLE]

for $y\to\pm\infty$. The root is defined to have its cut along the negative real axis and to assume positive values on the positive real axis. This choice of the root implies the following symmetry for the solution:

[TABLE]

Thus $\Omega$ is real on the real axis, see [24].

Numerically it is not convenient to impose boundary conditions at infinity. We thus assume that the wanted solution can be expanded in a Laurent series in $\sqrt{\zeta}$ around infinity. Such an asymptotic expansion is possible for the considered tritronquée solution in the sector $|\arg\zeta|<4\pi/5$. The formal series can be written there (see [24]) in the form

[TABLE]

where $a_{0}=1$, and where the remaining coefficients follow from the recurrence relation for $k\geq 0$

[TABLE]

This formal series is divergent, the coefficients $a_{k}$ behave asymptotically as $((k-1)!)^{2}$, see [24] for a detailed discussion.

It is known that divergent series can be used to get numerically acceptable approximations to the values of the sum by properly truncating the series. Generally the best approximations for the sum result from truncating the series where the terms take the smallest absolute values (see e.g. [18]). Since we work in Matlab with a precision of 16 digits and with values of $|\zeta|\geq 10$, we typically consider up to 10 terms in the series. In this case the terms corresponding to $a_{10}$ are of the order of machine precision ($10^{-14}$ and below).

Thus we have constructed approximations to the numerical values of the tritronquée solution for large values of $|\zeta|$. These can be used as in [24] to set up an initial value problem for the P-I equation and to solve this with a standard ODE solver. In fact the approach works well on the real axis starting from positive values until one reaches the first singularity on the negative real axis. It is straightforward to check the results of [24] with e.g. ode45, the Runge-Kutta solver in Matlab corresponding to the Maple solver used in [24]. If one solves P-I on a line that avoids the sector $|\arg\zeta|>4\pi/5$, one could integrate until one reaches once more large values of $|\zeta|$ for which the asymptotic conditions are known. This would provide a control on the numerical accuracy of this so-called shooting approach. Shooting methods are problematic if the second solution to the initial value problem has poles as is the case for P-I. In this case the numerical errors in the initial data (here due to the asymptotic conditions) and in the time integration will lead to a large contribution of the unwanted solution close to its poles which will make the numerical solution useless. It is obvious that P-I has such poles from the numerical results in [24] and the property (5.5). In [24] the task was to locate poles in the tritronquée solution, and in this case the shooting approach seems to be the only available. Here we are studying, however, the solution on a line in the complex plane where we know the asymptotic conditions for the affine parameter tending to $\pm\infty$.

Thus we use as in [15] the asymptotic conditions on lines avoiding the sector $|\arg\zeta|>4\pi/5$ to set up a boundary value problem for $y=\pm y_{0}$, $y_{0}\geq 10$. The solution in the interval $[-y_{0},y_{0}]$ is numerically obtained with a finite difference code based on a collocation method. The code bvp4c distributed with Matlab, see [39] for details, uses cubic polynomials in between the collocation points. The P-I equation is rewritten in the form of a first order system. With some initial guess (we use $\Omega=-\sqrt{\zeta/6}$ as the initial guess), the differential equation is solved iteratively by linearization. The collocation points (we use up to 10000) are dynamically adjusted during the iteration. The iteration is stopped when the equation is satisfied at the collocation points with a prescribed relative accuracy, typically $10^{-10}$. The solution for $a=i$ and $b=0$ is shown in Fig. 1.

The values of $\Omega$ in between the collocation points are obtained by interpolation via the cubic polynomials in terms of which the solution has been constructed. This interpolation leads to a loss in accuracy of roughly one order of magnitude with respect to the precision at the collocation points. To test this we determine the numerical solution via bvp4c for P-I on Chebychev collocation points and check the accuracy with which the equation is satisfied via Chebychev differentiation, see e.g. [45]. It is found that the numerical solution with a relative tolerance of $10^{-10}$ on the collocation points satisfies the ODE to roughly the same order except at the boundary points where it is of the order $10^{-8}$, see Fig. 2 where we show the residual $\Delta$ by plugging the numerical solution into the differential equation for the above example. It is straightforward to achieve a prescribed accuracy by requiring a certain value for the relative tolerance. Notice that we are not interested in a high precision solution of P-I here, but in a comparison of solutions to the NLS equation close to the point of gradient catastrophe of the semiclassical system with an asymptotic solution in terms of P-I transcendents. For this purpose an accuracy of the solution of the order of $10^{-4}$ will be sufficient in all studied cases.

The quality of the used boundary conditions via the asymptotic behavior can be checked by computing the solution for different values of $y_{0}$. One finds that the difference between the asymptotic square root and the tritronquée solution is only visible near the origin, see Fig. 3. For large $x$ it can be seen that the difference between the square root asymptotics and the tritronquée solution reaches quickly values below the aimed at threshold of $10^{-4}$. It is interesting to note that this difference is actually smaller than the difference between the tritronquée solution and the truncated formal asymptotic series except at the boundary, where the latter condition is implemented (see Fig. 3).

The dominant behavior of the square root changes if one approaches the critical lines $a=\exp(4\pi i/5)$, $b=0$. As can be seen from Fig. 4, the solution shows oscillations on top of the square root. The closer one comes to the critical lines, the slower is the fall off of the amplitude of the oscillations. We conjecture that these oscillations will have on the critical lines only a slow algebraic fall off towards infinity.

The above approach thus allows the computation of the tritronquée solution for a line avoiding the sector $|\arg\zeta|>4\pi/5$ with high accuracy. The picture one obtains by computing $\Omega$ along several such lines is that there are indeed no singularities in the sector $|\arg\zeta|<4\pi/5$, and that the square root behavior is followed for large $|\zeta|$. To obtain a more complete picture, we compute the tritronquée solution for $|\arg\zeta|<4\pi/5-0.05$ and $|\zeta|<R$ ( we choose $R=20$). The boundary data for $|\zeta|=R$ follow as before from the truncated asymptotic series, the data for $\arg\,\zeta=\pm 4\pi/5-0.05$ are obtained by computing the tritronquée solution on the respective lines as above.

To solve the resulting boundary value problem for the P-I equation is, however, computationally expensive since we have to solve an equation in 2 real dimensions iteratively. Since the solution we want to construct is holomorphic there, we can instead solve the harmonicity condition (the two dimensional Laplace equation) for the given boundary conditions. To this end we introduce polar coordinates $r$, $\phi$ and use a spectral grid as described in [45]: the main point is a doubling of the interval $r\in[0,R]$ to $[-R,R]$ to allow for a better distribution of the Chebychev collocation points. Since we work with values of $\phi<\phi_{0}$, we cannot use the usual Fourier series approach for the azimuthal coordinate. Instead we use again a Chebychev collocation method. The found solution in the considered domain is shown in Fig. 5 and Fig. 6.

The quality of the solution can be tested by plugging the found solution to the Laplace equation into the P-I equation. Due to the low resolution and problems at the boundary, the accuracy is considerably lower in the two dimensional case than on the lines. This is, however, not a problem since we need only the one dimensional solutions for quantitative comparisons with NLS solutions. The two dimensional solutions give nonetheless strong numerical evidence for the conjecture that the tritronquée solution has globally no poles in the sector $|\arg\zeta|<4\pi/5$.

7 Critical behavior in NLS and the tritronquée solution of P-I: numerical results

In this section we will compare the numerical solution of the focusing NLS equation for two examples of initial data for values of $\epsilon$ between $0.1$ and $0.025$ with the asymptotic solutions discussed in the previous sections, the semiclassical solution up to the breakup and the tritronquée solution to the Painlevé I equation. The numerical approach to solve the NLS equation is discussed in detail in [29]. For values of $\epsilon$ below $0.04$ we have to use Krasny filtering [30] (Fourier coefficients with an absolute value below $10^{-13}$ are put equal to zero to avoid the excitation of unstable modes). With double precision arithmetic we could thus reach $\epsilon=0.025$, but could not go below.

7.1 Initial data

We consider initial data where $u(x,0)$ has a single positive hump, and where $v(x,0)$ is monotonously decreasing. For initial data of the form $u(x,0)=A^{2}(x)$ and $v(x,0)=0$ where the function $A(x)$ is analytic with a single positive hump with maximum value $A_{0}$, the semiclassical solution of NLS follows from (2.11) with $f(u,v)$ given by

[TABLE]

where

[TABLE]

and where $x^{\pm}(\eta)$ are defined by $A(x^{\pm}(\eta))=i\eta$. The formula (7.1) follows from results by S.Kamvissis, K.McLaughlin and P.Miller in [26].

From $f(u,v)$ it is straightforward to recover the initial data from the equations

[TABLE]

Numerically we study the critical behavior of two classes of initial data, one symmetric with respect to $x$ which were used in [34], and initial data without symmetry with respect to $x$ which are built from the initial data studied in [43]. For the former class the corresponding exact solution of focusing NLS is known in terms of a determinant. Nonetheless we integrate the NLS equation for these initial data numerically since this approach is not limited to special cases, but can be used for general smooth Schwartzian initial data as in the latter case.

7.1.1 Symmetric initial data

We consider the particular class of initial data data

[TABLE]

Introducing the quantity

[TABLE]

we find that the semiclassical solution for these initial data follows from (2.11) with

[TABLE]

where

[TABLE]

For $\mu=0$ we recover the Satsuma-Yajima [37] initial data that were studied numerically in [34]. The function $f(u,v)$ takes the form

[TABLE]

which can also be recovered from (7.1) by setting $\rho=i$. The critical point is given by

[TABLE]

Furthermore we have

[TABLE]

where $r$, $\psi$ are defined in (2.13). For $A_{0}=1$ the initial data (7.3) coincides with the one studied in [43] by A.Tovbis, S.Venakides and X.Zhou. In the particular case $\mu=2$, $A_{0}=1$ the function $f(u,v)$ in (7.1) simplifies to

[TABLE]

In this case the critical point is given by

[TABLE]

Furthermore,

[TABLE]

7.1.2 Non-symmetric initial data

Recall that we are interested here in Cauchy data in the Schwartz class of rapidly decreasing functions. The above initial data are symmetric with respect to $x$, $u$ is an even and $v$ an odd function in $x$. To obtain a situation which is manifestly not symmetric, we use the fact that if $f$ is a solution to (2.5), the same holds for derivatives and anti-derivatives of $f$ with respect to $v$ and for any linear combination of those. If $f_{v}$ is an even function in $v$, this will obviously not be the case for a linear combination of $f$ and $f_{v}$.

As a specific example, we consider the linear combination

[TABLE]

where $f_{1}$ coincides with (7.8) and

[TABLE]

The function $f_{2}$ is obtained from the integration $f_{2,v}=f_{1}-v$. The critical point is given in this case by

[TABLE]

thus we have $|\alpha|<1/4$. Furthermore,

[TABLE]

such that

[TABLE]

We determine the initial data corresponding to $f$ for a given value of $|\alpha|<1/4$ by solving (7.2) for $u$, $v$ in dependence of $x$. This is done numerically by using the algorithm of [31] which is implemented as the function fminsearch in Matlab. The algorithm provides an iterative approach which converges in our case rapidly if the starting values are close enough to the solution, which is achieved by choosing $u$ and $v$ corresponding to $f_{0}$ as an initial guess. For $\alpha$ close to $1/4$ we observe numerically a steepening of the initial pulse which will lead to a shock front in the limit $\alpha\to 1/4$. For the computations presented here, we consider the case $\alpha=0.1$ which leads to the initial data shown in Fig. 7.

The initial data are computed in the way described above to the order of the Krasny filter on the interval $x\in[-15,11]$ on Chebychev collocation points. Standard interpolation via Chebychev polynomials is then used to interpolate the resulting data to a Fourier grid. To avoid a Gibbs phenomenon at the interval ends due to the non-periodicity of the data, we use a Fourier grid on the interval $[-10\pi,10\pi]$ to ensure that the function $u$ takes values of the order of the Krasny filter. For $x<-15$ and $x>11$, the function $u$ is exponentially small which implies the zero-finding algorithm will no longer provide the needed precision. Thus we determine the exponential tails of the solution to leading order analytically. We find for $x\to-\infty$

[TABLE]

and for $x\to+\infty$

[TABLE]

where $v_{\pm}=(\sqrt{1\pm\alpha}-1)/\alpha$. The initial data for the NLS equation in the form $\Psi=\sqrt{u}\exp(iS/\epsilon)$ are then found by integrating $v$ on the Chebychev grid by standard integration of Chebychev polynomials. The exponential tails for $S$ follow from (7.10) and (7.11). The matching of the tails to the Chebychev interpolant is not smooth and leads to a small Gibbs phenomenon. The Fourier coefficients decrease, however, to the order of the Krasny filter which is sufficient for our purposes. Thus we obtain the non-symmetric initial data with roughly the same precision as the analytic symmetric data.

7.2 Semiclassical solution

For times $t\ll t_{c}$, the semiclassical solution gives a very accurate asymptotic description for the NLS solution. The situation is similar to the Hopf and the KdV equation [19]. We find for the symmetric initial data for $t=t_{c}/2$ that the $L_{\infty}$ norm of the difference between the solutions decreases as $\epsilon^{2}$. More precisely a linear regression analysis in the case of symmetric initial data (for the values $\epsilon=0.03,0.04,\ldots,0.1$) for the logarithm of this norm leads to an error proportional to $\epsilon^{a}$ with $a=1.94$, a correlation coefficient $r=0.9995$ and standard error $\sigma_{a}=0.03$. In the non-symmetric case, we find $a=1.98$, $r=0.999996$ and $\sigma_{a}=0.003$.

Close to the critical time the semiclassical solution only provides a satisfactory description of the NLS solution for large values of $|x-x_{c}|$. In the breakup region it fails to be accurate since it develops a cusp at $x_{c}$ whereas the NLS solution stays smooth. This behavior can be well seen in Fig. 8 for the symmetric initial data.

The largest difference between the semiclassical and the NLS solution is always at the critical point. We find that the $L_{\infty}$ norm of the difference scales roughly as $\epsilon^{2/5}$ as suggested by the Main Conjecture. More precisely we find a scaling proportional to $\epsilon^{a}$ with $a=0.38$ and $r=0.999997$ and $\sigma_{a}=4.2*10^{-4}$. For the non-symmetric initial data, we find $a=0.36$, $r=0.9999$ and $\sigma_{a}=0.002$. The corresponding plot for $u$ can be seen in Fig. 9.

The function $v$ for the same situation as in Fig. 8 is shown in Fig. 10. It can be seen that the semiclassical solution is again a satisfactory description for $|x-x_{c}|$ large, but fails to be accurate close to the breakup point.

The phase for the non-symmetric initial data can be seen in Fig. 11.

In the following we will always study the scaling for the function $u$ without further notice.

7.3 Multiscales solution

It can be seen in Fig. 8 and Fig. 10 that the multiscales solution (5) in terms of the tritronquée solution to the Painlevé I equation gives a much better asymptotic description to the NLS solution at breakup close to the breakup point than the semiclassical solution for the symmetric initial data. For larger values of $|x-x_{c}|$, the semiclassical solution provides, however, the better approximation. The rescaling of the coordinates in (5) suggests to consider the difference between the NLS and the multiscales solution in an interval $[-\gamma\epsilon^{4/5},\gamma\epsilon^{4/5}]$ (we choose here $\gamma=1$, but within numerical accuracy the result does not depend on varying $\gamma$ around this value). These intervals can be seen in Fig. 12. We find that the $L_{\infty}$ norm of the difference between these solutions in this interval scales roughly like $\epsilon^{4/5}$. More precisely we have a scaling $\epsilon^{a}$ with $a=0.76$ ($r=0.998$ and $\sigma_{a}=0.019$).

For the non-symmetric initial data, the situation at the critical point can be seen in Fig. 9 and Fig. 11. Again the multiscales solution (5) gives a much better description close to the critical point than the semiclassical solution. However, the approximation is here much better on the side with weak slope for $u$ than on the side with strong slope. We consider again the $L_{\infty}$-norm of the difference between the multiscales and the NLS solution in the interval $[-\gamma\epsilon^{4/5},\gamma\epsilon^{4/5}]$. The scaling behavior of the solution can be seen in Fig. 13. For $\gamma=1$ we find $a=0.71$, $r=0.998$ and $\sigma_{a}=0.02$. These values do not change much for larger $\gamma$. For smaller $\gamma$ there are not enough points to provide a valid statistics. The value of $a$ smaller than the predicted $4/5$ is seemingly due to the strong asymmetry in the quality of the approximation of NLS by the multiscales solution as can be seen from Fig. 9. In the considered interval, the deviation is already so big that the scaling no longer holds as in the symmetric case. To study the scaling with a reliable statistics would, however, require the use of a considerably higher resolution which would be computationally too expensive.

Going beyond the critical time, one finds that the real part of the NLS solution continues to grow before the central hump breaks up into several humps. Notice that the multiscales solution always leads to a function $u$ that is smaller than the corresponding function of the NLS solution at breakup and before. This changes for times after the breakup as can be inferred from Fig. 14 which shows the time dependence of the NLS and the corresponding multiscales solution for the non-symmetric initial data. The approximation is always best at the critical time.

To study the quality of the approximation (5), we use rescaled times. The scaling of the coordinates in (5) suggests to consider the NLS solution close to breakup at the times $t_{\pm}(\epsilon)$ with

[TABLE]

where $\beta$ is a constant (we consider $\beta=0.1$). We will only study the symmetric initial data in this context. Before breakup we obtain the situation shown in Fig. 15. It can be seen that the multiscales solution always provides a better description close to $x_{c}$ than the semiclassical solution, and that the quality improves in this respect with decreasing $\epsilon$. We find that the $L_{\infty}$ norm of the difference scales in this case as $\epsilon^{a}$ with $a=0.55$ ($r=0.994$ and $\sigma_{a}=0.03$).

The situation for times after breakup can be inferred from Fig. 16. Close to the central region the multiscales solution shows a clear difference to the NLS solution. But it is interesting to note that the ripples next to the central hump are well approximated by the Painlevé I solution. The $L_{\infty}$ norm of the difference between the two solutions scales roughly like $\epsilon$. More precisely we find a scaling $\epsilon^{a}$ with $a=1.02$ ($r=0.9999$ and $\sigma_{a}=7.7*10^{-3}$).

8 Concluding remarks

In this paper we have started the study of the critical behavior of generic solutions of the focusing nonlinear Schrödinger equation. We have formulated the conjectural analytic description of this behavior in terms of the tritronquée solution to the Painlevé-I equation restricted to certain lines in the complex plane. We provided analytical as well as numerical evidence supporting our conjecture. In subsequent publications we plan to further study the Main Conjecture of the present paper by applying techniques based, first of all, on the Riemann - Hilbert problem method [26, 43, 44] and the theory of Whitham equations (see [19] for the numerical implementation of the Whitham procedure in the analysis of oscillatory behavior of solutions to the KdV equations). The latter will also be applied to the asymptotic description of solutions inside the oscillatory zone. Furthermore we plan to study the possibility of extending the Main Conjecture to the critical behavior of solutions to the Hamiltonian perturbations of more general first order quasilinear systems of elliptic type. Last but not least, it would be of interest to study the distribution of poles of the tritronquée solution in the sector $|\arg\zeta|>\frac{4\pi}{5}$ and to compare these poles with the peaks of solutions to NLS inside the oscillatory zone. The elliptic asymptotics obtained by Kitaev [28] might be useful for studying these poles for large $|\zeta|$.

In this paper we did not study the behaviour of solutions to NLS near the boundary $u=0$. Such a study is postponed for a subsequent publication.