Transverse Evolution Operator for the Gross-Pitaevskii Equation in Semiclassical Approximation

The Gross-Pitaevskii equation with a local cubic nonlinearity that describes a many-dimensional system in an external field is considered in the framework of the complex WKB-Maslov method. Analytic asymptotic solutions are constructed in semiclassical approximation in a small parameter $\hbar$, $\hbar\to 0$, in the class of functions concentrated in the neighborhood of an unclosed surface associated with the phase curve that describes the evolution of surface vertex. The functions of this class are of the one-soliton form along the direction of the surface normal. The general constructions are illustrated by examples.


Introduction
The Gross-Pitaevskii equation (GPE), derived independently by Gross [1] and Pitaevskii [2], arises in various models of nonlinear physical phenomena. This is a Schrödinger-type equation with an external field potential U ( x, t) and a local cubic nonlinearity: Here x = (x j ) ∈ R n ; i, j, k, l = 1, . . . , n; t ∈ R 1 ; ∂ t = ∂/∂t; p = −i ∇; ∇ is a gradient operator in x; κ and are real parameters. The complex function Ψ( x, t, ) determines the system state, |Ψ| 2 = |Ψ( x, t, )| 2 = Ψ( x, t, )Ψ * ( x, t, ); the function Ψ * ( x, t, ) is complex conjugate to Ψ( x, t, ). For U = 0, κ < 0, and n = 1, equation (1) is known to be integrable by the Inverse Scattering Transform (IST) method [3,4] and to have exact soliton solutions. The solutions are written in terms of space localized functions which conserve their form during the evolution. Solitons have numerous physical applications; for example, they serve to describe the propagation of optical pulse in nonlinear media [5,6].
The solutions of equation (1) that are not completely localized are also of interest, in particular, in the theory of nonlinear deep-water waves where the two-dimensional (n = 2) solutions of plane wave solitonlike type are localized only along the wave propagation direction [7]. The Gross-Pitaevskii equation (1) in physical dimensions (n = 2, 3) is used in the meanfield quantum theory of Bose-Einstein condensate (BEC) formed by ultracold bosonic coherent atomic ensembles [8]. As a rule, the function U ( x, t) is the potential of an external field of a magnetic trap and laser radiation. The wave function Ψ( x, t, ) corresponds to a condensate state. The local nonlinearity term κ|Ψ( x, t, )| 2 arises from an assumption about the deltashape interatomic potential. For multidimensional cases, no method of exact integration of the GPE is known, to our knowledge, and equation (1) is usually studied for some special cases by computer simulation. However, this approach has natural restrictions and also has nontrivial specific features as long as localized solutions of equation (1) are known to be unstable and to collapse when n > 1, U = 0, and κ < 0 [9]. Therefore, to develop methods for constructing analytical solutions of the GPE (1) for many-dimensional cases is critical.
In this paper we construct a class of asymptotic solutions for the (1 + n)-dimensional GPE with a focusing local cubic nonlinearity L(Ψ)( x) = −i ∂ t + H(ˆ p, x, t) − g 2 |Ψ( x, t, )| 2 Ψ( x, t, ) = 0. ( Here g is a real nonlinearity parameter and is an asymptotic parameter, → 0. The pseudodifferential operator H(ˆ p, x, t) is Weyl-ordered [10] and its symbol H( p, x, t) is quadratic in where a, b = n l=1 a l b l ; the real functions -the (n × n)-matrix H pp (t), the vector H( x, t), and the scalar H 0 ( x, t) -smoothly depend on their arguments. These functions model the external fields imposed on the system described by equation (2). A formalism of semiclassical asymptotics was developed for the GPE with a nonlocal nonlinearity in [11,12]. This equation was named the Hartree type equation. The semiclassical asymptotics of the nonstationary Hartree type equation were also studied [13,14,15,16,17,18,19,20]. This formalism is based on the WKB-Maslov complex germ method, the idea of which was first put forward in [21] and then comprehensively treated in [22,23,24].
Further developments of the Maslov method [25,26] (see also [27] and references therein) allow one to construct the localized solutions of the Cauchy problem in the class of trajectoryconcentrated functions.
Here we apply the ideas of the complex germ method and the approach used in [11,12] to construct the analytic solutions asymptotic in a small parameter , → 0, for equation (2). The asymptotic solutions are constructed in the class of functions localized in the vicinity of an unclosed surface associated with the phase curve that describes the evolution of the surface vertex.
The functions of this class have the one-soliton form of the one-dimensional nonlinear Schrödinger equation along the surface normal direction. A semiclassical linearization of equation (2) is performed accurate to O( 3/2 ), → 0. The asymptotic solutions constructed are illustrated by examples.

The class of paraboloid-concentrated solitonlike functions
Let us describe a class of functions Ψ( x, t, ) in which the asymptotics for the GPE (2) are constructed. Following [11,12], consider equation (2) for the (1 + 1)-dimensional case with no external field (U = 0) where x ∈ R 1 and ∂ x = ∂/∂x. Equation (4) has a one-soliton solution [3,4] Ψ(x, t, where ξ, η, x 0 , and ϕ 0 are the real parameters. Earlier [28] we obtained a semiclassical solution of equation (2) with the leading term It was assumed that the complex function Ω( x, t, ) = S( x, t, ) + iσ( x, t, ) regularly depends on and its high-order part, i.e., the function Ω( , is a complex solution of the Hamilton-Jacobi equation The functions (6) are localized in the neighborhood of a t-parameter family of surfaces Assume that the rank of Hesse matrix of the function σ (0) ( x, t) is n at any point x ∈ Γ t . The construction of global solutions of equation (7) is beyond the scope of this work. Here we consider the evolution of the functions localized in the neighborhood U x 0 of the surface point x 0 (8) at the initial time t = 0.
Denote by a smooth curve passing through the point x 0 (= X(0)) in the coordinate space x ∈ R n such that X(t) ∈ Γ t . This curve plays the role of the "classical trajectory" corresponding to the solution of the quantum equation (2). In the neighborhoodŨ X(t) ⊂ Γ t of the point X(t) ∈ γ, the hypersurface Γ t can be approximated by a simpler surface, for example, by the tangent plane. It is well known that a tangent plane is determined uniquely by its normal and the point of contact. Let π(t) be a normal vector to the surface Γ t at a point x = X(t).
To construct asymptotic solutions of equation (2) one needs the complex germ [22], i.e. the n-dimensional complex space associated with the equation under consideration and functions defined in the neighborhood of the hypersurface (8).
It is more precise to use a second order surface instead of the tangent plane. The surface (9) with a proper choice of the matrix Q 2 (t) is a quadratic approximation of (8) in the neighborhoodŨ X(t) . Introduce the class N t of complex functions with the generic element Ψ( x, t, ) = Φ(θ, x, t, ) = g −1 π(t), H pp (t) π(t) exp i S( x, t, ) The functions (10) are localized in the neighborhood U X(t) of the point x = X(t). This point is the "vertex" of the surface (9). In equation (10), θ = σ( x, t, )/ is a "fast" variable; the real functions u = u(θ, x, t, ) and v = v(θ, x, t, ) are regular in their arguments and are bounded in θ, and In equations (11) and (12), the real function S 0 (t), the real n-dimensional vector functions π(t), P (t), and X(t), the real (n × n)-matrix functions Q 1 (t) and Q 2 (t), the real functions Note that in constructing the class N t , a specific basis, "orthogonal" in some way to the vector π(t), is used 1 . As a result, the classical trajectory depends on the complex germ as well.
The functions (10) are not normalizable in the space L 2 (R n , d x) since they are localized on an unclosed surface, equation (9). Therefore, the surface Γ t (8) (compact in some physical applications) is substituted by its quadratic approximation (9) in the neighborhood U X(t) . As we are interested in finding functions localized in U X(t) , we come to the problem of small perturbations on the "background" of the hypersurface (9). To eliminate the "background", consider the functioñ The functionsΨ( x, t, ), equation (13), are normalizable in L 2 (R n , d x). Estimating the solutions of equation (2), we use the norm of the functionsΨ( x, t, ). Taking into account this normalization condition, let us assume that the functionsσ 1 ( ξ, t, ) andS 1 ( ξ, t, ) belong to the Schwartz space S(R n−1 ξ , dµ), where dµ is the measure on the hypersurface, determined by equation (9).
The class N t of the form (10) with n = 1 includes the exact one-soliton solution (5) of onedimensional nonlinear Schrödinger equation (4) with a special choice of the functional parameters of the class. If the rank of matrix Q 2 (t) is n, then in one-dimensional case equation (10) can not be transformed to (5) since the argument θ of the hyperbolic cosine contains the terms quadratic in x. The sufficient conditions for such a transformation are where (Q 2 (t), π(t)) is the augmented matrix of order n×(n+1). Note that under conditions (14), the surface (9) has the paraboloidal shape. Let us expand the operator L in equation (2) in the neighborhood U X(t) under the conditions The leading term of the asymptotic in the class of functions (10) under conditions (15) has the form (6). It is determined by the phase curve z = Z(t, ) = ( P (t, ), X(t, )), vector π = π(t), and functions σ( x, t, ) and S( x, t, ) of the form (11) and (12), respectively. We call these functions the paraboloid-concentrated solitonlike functions.

The linear associated Schrödinger equation
The asymptotic leading term is obtained when the asymptotic solution of the equation (2) is constructed accurate to O( 2 ) [28].
Let us substitute (10) in (2), take into account the terms of order O( 2 ), separate the equations with respect to the "fast" variable, and find the solutions of these equations that decrease as θ → ±∞. We than obtain for the function Ω = S + iσ With the substitution Ω = −i ln ψ( x, t, ) π, H pp (t) π −1/2 , equation (16) yields where the operator H(ˆ p, x, t) has the form of (3), We call the equation the linear associated Schrödinger equation corresponding to the GPE (2). Therefore, the leading term of the asymptotic solution of the nonlinear equation (2) is constructed by using a solution of the linear equation (18) with conditions (6) and (17).
Note that for the function (6) to satisfy (14) the solution (17) must be chosen in a specific way. Denote the class of these functions by S t .
Let us seek the solutions to equation (18) that possess this property in the class of functions Substitution of (19) in (18) yields Let us expand the operator of equation (20) in a power series in ∆ x accurate to O 3/2 in terms of the accuracy of (15). We then have Therefore, for the solution (19) in the approximation under consideration, the linear associated Schrödinger equation (18) takes the form Here the operatorĤ 0 (t) (quadratic with respect to the operators ∆ x and ∆ˆ p ) is obtained from (21). Thus, the construction of asymptotic leading term (6) in the class of functions (19) for the nonlinear equation (2) by solving equation (22) is complete.

Solutions of the Gross-Pitaevskii equation
The solutions of equation (22) are well-known (see, e.g., [29,30,31]). In particular, the evolution operatorÛ H 0 (t, s) is found in explicit form, whose action on the function ψ( x, ) referred to a time t = s is given by the relation Here G H 0 ( x, y, t, s) is the Green function of the linear associated Schrödinger equation (22), which is determined by the following conditions: To construct a subclass of the class N t (10) of solutions of equation (2), let us find a particular solution in the form where the functions S 0 (t), S 1 (t), and σ 1 (t), the n-dimensional vectors P (t), X(t), and the complex (n × n)-matrix Q(t)=Q 1 (t) + iQ 2 (t) are to be determined; N φ is a constant.
Substituting (25) in (21) and setting˙ X = H p , Q(t) = Q 1 (t) + iQ 2 (t), we obtain the determining system of equations for P (t), π(t), X(t), Q(t) Consider the following system of equationṡ Equations (28) are called a system in variations in vector form [27].
In general, it has n complex linear independent solutions, which can be written in the form of 2n-dimensional vector columns For equation (22) the vectors (29) set the symmetry operators which are linear in coordinates and momentâ since the conditions are satisfied because of the validity of equation (28); [Ĥ 0 (t),â j (t)] =Ĥ 0 (t)â j (t) −â j (t)Ĥ 0 (t) is the commutator of the linear operators. Let B(t) and C(t) denote the (n × n)-matrices whose columns are constructed of the vectors of the solutions of system (28): Then Q(t) = B(t)C −1 (t) and W (t) = Q(t) Z(t), and the system in variations (28) is equivalent to (27) (see [27]). The normalization factors N a j are chosen as follows. The matrices B(t) and C(t) satisfy the condition Here C + is Hermitian adjoint to C. In addition, if the matrix Q(0) = B(0)C −1 (0) is symmetric at the initial time t = 0, Q(0) = Q(0) ⊺ , then Q(t) = Q(t) ⊺ , t > 0. Let Then the operators (30) satisfy the ordinary commutation relations [â j (t),â + k (t)] = δ jk , j, k = 1, n − 1, and all the other commutators are equal to zero.
Let us denote the function ψ( x, t, ) in (19) as where the function φ( x, t, ) of the form (25) is denoted as |0, t . The normalization factor N 0 is determined from the normalization condition for the function Ψ( x, t, ) that is a solution of the nonlinear equation (2). Let us substitute (32) in (22); then Here ( P (t), X(t)) is a solution of the first pair of equations (26), the (n × n)-matrix C(t) is composed of the columns Z j (t). Note in addition that ∆ x, ImQ(t)∆ x ≥ 0 and a j (t)|0, t = 0, j, k = 1, n − 1.
Using expressions similar to (35)-(37), we obtain the function for the function (6) corresponding to ψ  ν ( x, t, ). The generic element of class K t is given by The conditioñ a n (t)ψ K ( x, t, ) =â n (t)φ K ( x, t, ) = 0 (42) is valid for the functions of the class K t . This follows immediately from (33), (34), and (38). Let φ K ( x, s, ) be a function of the class K s referred to a start time s, andÛ H 0 (t, s) is the evolution operator given by (23) and (24).
Then the function also belongs to the class K t , since the function is represented as The above relation follows from the uniqueness of solution of the Cauchy problem for equation (22). On the other hand, we have for the symmetry operatorâ n (t) of equation (22) the relation a n (t) =Û H 0 (t, s)â n (s)Û −1 H 0 (t, s) following from (23) and (31). If the function φ K ( x, s, ) belongs to the class K s at the initial time s and satisfies the condition a n (s)φ K ( x, s, ) = 0, then condition (42) is also satisfied for the function (43) at a time t > s, i.e., we havê a n (t)φ K ( x, t, ) = 0.
Therefore, the function φ K ( x, t, ) belongs to the class K t . The evolution operatorÛ H 0 (t, s) of equation (22) induces the evolution operatorÛ tr (t, s, ·) for the nonlinear equation (2) in the class S t .
Let a function ψ( x, ) be referred to an initial time s. According to (17), the function Ψ s ( x, ) corresponds to ψ( x, ). Then the function Ψ( x, t, ), determined by (23), corresponds to ψ( x, t, ). This correspondence can be written as a result of the action of the evolution operatorÛ tr (t, s, ·) on the initial function Ψ s ( x, ), Using the evolution operatorsÛ H 0 (t, s, ·) andÛL(t, s, ·) of the forms (23) and (45) respectively, we can define the symmetry operators for the nonlinear equation (2) in the class of functions S t .
To that end, let us take the function ϕ( x, ) from the class K s which satisfies condition (44) at an initial time t = s. Let ϕ( x, t, ) be a function obtained from Eq. (17) and Ψ( x, t, ) be the solution of the nonlinear equation (2) related to ϕ( x, t, ) according to (17). Consider the operatorÂ,Â : K s → K s such that [Â,â n (s)] = 0, and a function ϕ A ( x, ) =Âϕ( x, ) ∈ K s . Then the function Ψ A ( x, t, ) related to ϕ A ( x, ) by (17) can be treated as a result of the action of the symmetry operatorÂ nl of the nonlinear equation (2) Ψ A ( x, t, ) =Â nl Ψ( x, t, ).
In conclusion, we note that for a nonlinear Schrödinger equation with a focusing nonlinearity, the many-dimensional solutions localized at the initial time are unstable. This leads to the phenomenon of collapse in the course of evolution. The semiclassical asymptotics (39) behave in a similar manner. They can be constructed for special external fields within finite time intervals where singularities typical of collapse appear.

5
The three-dimensional anisotropic oscillator In case of a harmonic oscillator field, the linear operatorĤ, equation (3), reads where K = diag{k 1 , k 2 , k 3 } with k j ∈ R 1 , j = 1, 3. The GPE (2) then takes the form To construct the asymptotic solutions (6) for equation (47) we solve the dynamic system (26) and (28), which is reduced tȯ in the case under consideration. Here W j (t) and Z j (t), j = 1, 3, are the linear independent solutions of the system in variations (49), which determine the matrices To write down the solution of systems (48) and (49) and the function Ψ ν ( x, t, ), we introduce the notation Ω + j = k j /m for k j > 0, Ω − j = −k j /m for k j < 0, a j (t) = W j (t), Z j (t) ⊺ , and {a j (t), a k (t)} denoting the skew-scalar product of the 6-vectors a j (t), a k (t), j, k = 1, 3.
Define the operatorÛ (0) tr by its action on the function ϕ( Here S 0 ± is given by (50) with P (t) = P (t) ± and X(t) = X(t) ± , respectively. We call the operatorÛ as the initial function, i.e., at t = 0. Below we put m = 1, g = 1, K = E for Ψ + = Ψ + 00 , and K = −E for Ψ − = Ψ − 00 , where E is the unit matrix and Ψ ± 00 is the "vacuum" function (52). Then, for the (1+1)-dimensional case, the module of function (52) under the initial conditions P (0) = 0, x(0) = 0, π 0 = 1, σ 1 (0) = 0 take the form In Fig. 2a, the function |Ψ + | (56) is shown for the (1 + 1) space-time and the potential well (kx 2 /2, k > 0) of the oscillator form in (46) with = 1. It can be seen that the solution, being localized at the initial time, is focused within the time interval (0, π/2) (in the process of evolution) and collapses. The GPE with an external field was used to describe the dispersion and diffraction of nonlinear waves [32]. The collapse phenomenon was studied in [32] by numerical simulations. The numerical results qualitatively correlate with those obtained by analytic asymptotic methods. The asymptotic solution (53) describes the system for any T within a finite time interval [0, T ]. If the time interval is short enough, the evolution can end before the collapse occurs.
The function (56) is π-periodic and its graph is given in Fig. 2b. The term (− /2) ln | cos(t)| in the argument of hyperbolic cosine in equation (56) shifts the function centroid (56) within  the time interval 0 < t < π/2. Throughout the evolution time, this gives rise to oscillations of the function maximum about the plane x = 0.
For the potential hill of the oscillator form (kx 2 /2, k < 0) in equation (46), the state of the system in the (1 + 1) space-time is described by Ψ − . Let us choose the initial condition for the function Ψ − the same as for Ψ + . Then the dynamics of the system is characterized by defocusing and exponential damping (see Fig. 3). The term (− /2) ln(ch(t)) in the argument of hyperbolic cosine in (57) shifts the function centroid in the negative direction of x.
Consider the (2 + 1)-dimensional case. For the functions (52) put P (0) = 0, x(0) = 0, π 0 = 1, σ 1 (0) = 0, then we have The functions (58) and (59) coincide at t = 0. The initial graph of the functions (58) and (59) is presented in Fig. 4a. Figs. 4b,c show the graphs of function (58) for t = 1.2 and of the function (59) for t = 2, respectively. The solution considered is localized about the parabola with the vertex at the point P (t) = 0, X (t) = 0 of phase space. In the course of evolution of |Ψ + |, the parabola branches disperse and at t = π/2 it transforms into a straight line. The inverse behavior of the parabola takes place for the evolution of |Ψ − |: the parabola branches converge for an infinite time. The amplitude of the functions |Ψ ± | behaves, in the direction along the coordinate x, as in the (1 + 1)-dimensional case.

Conclusion remarks
The asymptotics obtained, equation (6), can be regarded as a necessary step in the construction of a global semiclassical asymptotic solution to the GPE (2). It may be supposed that the functions (6) describe the behavior of an element of the global solution in the neighborhood of a ray along the normal π ( Fig. 1) to the (closed) surface about which the global solution is concentrated. Substantiation of these global asymptotics for finite times t ∈ [0, T ], T = const, is a special nontrivial mathematical problem. This problem is concerned with obtaining a priori estimates for the solution of a nonlinear equation, which are uniform in parameter ∈]0, 1] and is beyond the scope of the present work. Note that, in view of the heuristic considerations given in [15], it seems that the difference between the exact and the constructed formal asymptotic solution can be found with the use of method developed in [15,23].
The technique of construction of semiclassical asymptotics developed in Section 4 provides a way for solving the problem of correspondence between the classical and quantum results for quantum systems described by nonlinear equations, namely, via finding solutions to the dynamic system (26) and (28). For nonlinear systems this problem differs from the relevant problem in the linear quantum mechanics. For linear quantum systems, a correct formal transition from the quantum theory to the classical one requires imposing special limitations on the quantum states (their semiclassical concentration). The states not satisfying these limitations are regarded as "essentially quantum", and those satisfying them are considered "near to classical". Thus, the dynamics of the classical objects obtained is described by the classical Hamiltonian equations no matter the domain where the wave function is concentrated (whether it be a point, a curve, or a surface). The study of the dynamic system (26) and (28) is a separate mathematical subject for research. For example, the Hamiltonian or Poissonian formalisms as applied to this system are of interest.
The construction of "transverse" evolution operator given in Section 4 allows one not only to obtain semiclassical asymptotics but also to construct approximate symmetry operators of special type acting in the class of S t functions under consideration. Such operators can be naturally referred to as semiclassical symmetry operators [33] (see also [12]).