The Cauchy problem for the generalized hyperbolic Novikov-Veselov equation via the Moutard symmetries

We begin by introducing a new procedure for construction of the exact solutions to Cauchy problem of the real-valued (hyperbolic) Novikov-Veselov equation which is based on the Moutard symmetry. The procedure shown therein utilizes the well-known Airy function $\Ai(\xi)$ which in turn serves as a solution to the ordinary differential equation $\frac{d^2 z}{d \xi^2} = \xi z$. In the second part of the article we show that the aforementioned procedure can also work for the $n$-th order generalizations of the Novikov-Veselov equation, provided that one replaces the Airy function with the appropriate solution of the ordinary differential equation $\frac{d^{n-1} z}{d \xi^{n-1}} = \xi z$.


I. INTRODUCTION
The Moutard transformation [1] is a very interesting form of discrete symmetry of second-order linear equations with variable coefficients. Consider the differential equation where ∆ is the two-dimensional Euclidean Laplacian. Let ψ = ψ(ξ 1 , ξ 2 ) and φ = φ(ξ 1 , ξ 2 ) be two partial solutions of this equation, i.e. its solutions for two different Cauchy problems (with different initial conditions). The function φ we will call the prop solution because it will play the role of a foundation on which we will build the Moutard's mathematical apparatus. This apparatus, known as the Moutard symmetry is to be defined as the following transformation: where and we used the following (standard) tensor notations: µ ∈ {1, 2}; ∂ µ = ∂/∂ξ µ ; ε µν is a fully antisymmetric tensor with ε 12 = 1; and, as usual, we have a summation over the repeated indices. Note that the one-form being integrated in (3) is a closed one when ψ and φ are both solutions of (2). Hence, the shape of the contour of integration Γ in (3) is irrelevant. If we switch from (ξ 1 , ξ 2 ) to the cone variables x and y, then the equation (1) will take a new form: To define the Moutard symmetry in the new cone variables one has to override the closed one-form θ. Namely, let Ψ = Ψ(x, y) and Φ = Φ(x, y) be two solutions of the (4): Define a differential form dθ[Ψ; Φ] such that and Note that since by definition both Ψ and Φ are solutions of (5), the one-form is closed too, i.e.
and thus the shape of the contour of integration Γ is irrelevant. The Moutard symmetry has the form It means that It is worth pointing out at this step that rather then zero. The Moutard symmetry in usual variables can be formally reduced in a one-dimensional limit to the Darboux transformation (originally introduced in 1882 in article [2] and further developed by Crum in [3]). In other words, the Moutard symmetry group includes the Darboux transformation group as a subgroup. This is potentially a very significant observation since the group structure of the Darboux transformation is associated with a subalgebra of the Kats-Moody algebra of the group SU(2) by means of the following commutation relation: where c µν ρ are the structure constants of the group SU(2) (i.e. the group of unitary 2 × 2 matrices with determinant 1, whose generators are the famous Pauli matrices.) [4] . It is important to note, however, that while the Darboux transformation is indeed a limit case of the Moutard transformation in the one-dimensional limit, there still exist a crucial difference between the two. While both require two eigenfunctions to jump start the process of transformation, the Darboux transformation allows them to belong to different eigenvalues, as long as the prop function is positively defined everywhere (to ensure the regularity of a new dressed potential). In contrast to that, the Moutard symmetries demands a pair of functions to belong to the same eigenvalue, which in this paper for simplicity is set equal to zero.
Another interesting property of the Moutard symmetries is that they ((8), (11) , (9) or (2), (3)) can be iterated several times, and the result of their application can be expressed via the corresponding Pfaffian forms [5]. The Moutard symmetries have a number of interesting physical applications, ranging from the Maxwell equations in a 2D inhomogeneous medium [6] to the problems of quantum cosmology and the Wheeler-DeWitt equation [7], as well as the quantum (phantom) modified gravity models [8]. But the Moutard transformations are most effective in the theory of integrable evolutionary equations in (1 + 2) dimensions.
In 1984 two Soviet mathematicians, Sergey Novikov and Aleksander Veselov, made a very important discovery [9]. At that time they were studying the class of two-dimensional Schrödinger operators L with a property of a "finitezoneness w.r.t. a single energy level", originally introduced in 1976 by Dubrovin, Krichever and Novikov in [10], and is a property of a Schrödinger operator L whose Bloch functions (i.e. the eigenfunctions that L shares with the periodic operators of spatial translations) on a single energy level is meromorphic (i.e. holomorphic everywhere except for a few isolated poles) on a Riemann surface of a finite genus ( [10], see also [11]). More specifically, Veselov and Novikov were interested in those "finite-zone" 2D Schrödinger operators L that are simultaneously real and purely potential. For that end, a theorem has been proven that for eigenfunction ψ of such an operator to be meromorphic everywhere except at two points and possess the required asymptotes (see [9] for details), those eigenfunctions necessarily have to satisfy the following system of equations: where the overhead bar denotes complex conjugate, the operators and the coefficients a j and V are uniquely defined by the aforementioned asymptotes. The pair of operators L and A coupled by the system (12) deeply resembles the famous Lax pair L, A [12]: two time-dependent operators that satisfy the condition which in turn produces a new (1 + 1)-dimensional differential equation -with such notable examples as, Korteweg-de Vries (KdV), sine-Gordon and the nonlinear Schrödinger equations (for more details cf. [12] and [5], see also [13]). However, in the case of (12) the proper "entanglement" between the L and A was a tad more complex, and required an additional operator B n (similar in form to A n−1 albeit with its own coefficients): In other words, the system (12) was essentially a letter of acquaintance from an infinitely diverse family of novel (1 + 2) dimensional equations. And the very first (and the simplest) member of that very family was the n = 1 equation (with A 1 = ∂ 3 + 3ω(x, y, t)∂) that has been henceforth known as the Novikov-Veselov equation (NV) [9]: Subsequent studies of NV equation has produced a lot of very interesting observation: for example, in a onedimensional case (i.e. when both u and ω are independent of y variable) the (14) reduces to the KdV equation; whereas if we add to (14) one additional term, λ∂ω, and then take a limit λ → ±∞, we will instead wound up with either one of two Kadomtsev-Petviashvili equations [14] -another famous (1 + 2) generalization of KdV! In fact, a huge and ever-growing body of works related to the study of NV equations has been established (see, for example, [15], [16], [17] and also [18] for a rather extensive review of a recent literature on the subject). In particular, a lot of spotlight has been focused on the solutions of (14) and (15). For example, the article [19] shows how the method of superposition originally proposed in [20] can be used to obtain a 2-solitary wave solution of the Novikov-Veselov equation. The more general N-solitons solutions were subsequently constructed in [21]. A conspicuous absence of exponentially localized solitons for NV equation with a positive energy was explored and explained in [22], whereas the impossibility of such solutions for the negative energy NV was proven in [23]. Many articles were dedicated to unusual and fascinating properties of the multi-dimensional solutions, including those for seemingly ordinary flat waves. In particular, in [24] it has been shown that plane wave soliton solutions of NV equation are not stable for transverse perturbations; the paper [18] demonstrates that NV equation permits such interesting solutions as multi-solitons, ring solitons, and the breathers; while the authors of [25] construct a Mach-type soliton of the NV equation. One of the most effective mathematical tools for studying the NV equation is the inverse scattering method. It was developed and applied in many articles, such as, for example, [26] and [16]. We must also mention an important paper [17], which looked at a zero-energy Novikov-Veselov equation for the initial data of conductivity type. Taking into account that the (1 + 2) nonlinear equations to this day remains mostly "terra incognita", it generates a lot of attention when someone manages to establish a relationship between the various members of a small group of currently known integrable models. As one such example we can refer to the article [27] which has uncovered a curious relationship between the hyperbolic NV equation and the stationary Davey-Stewartson II equation -here an adjective "hyperbolic" simply means that the L equation (4) is hyperbolic, i.e. that both x and y variables are real (accordingly, since the "original" NV equation is associated with the system (12), it can be called "elliptic"). Finally, a lot of literature has been written on the subject of various generalizations of NV equations, of their properties and of their solutions [28], [29], [30], and see also [31], where the analogue of NV equations is shown to arise in the nonlinear optics in a dispersion-free limit.
The observant reader will of course notice that one of the most prominent aspects of the majority of the articles we have mentioned is an almost universal adoption of an inverse scattering method as a primary tool for conducting the research and finding the exact solutions of NV. However, in this article we wish to discuss an alternative method of solving the Cauchy problem for NV (the hyperbolic version). This method, albeit simple in principle, appears to be deep enough to be applicable to a very broad class of equations, NV being just the first one -just as it is the but a first member of the Novikov-Veselov hierarchy.
The article is organized as follows. In Sec. II we introduce all the necessary ingredients of our proposed method, namely: the Lax pair for the hyperbolic (real-valued) NV equation, the Moutard transformation and the Airy functions -and describe how to use them to produce the exact solutions to the NV equation. In the next section, Sec. III, we up the ante by adding the initial conditions into the mix -and show how to make sure the new solution satisfies those conditions. Finally, in Sec. IV we discuss the generalization of the proposed method to the higher-order equations from the Novikov-Veselov hierarchy.

II. THE MOUTARD TRANSFORMATION
Let us start by introducing the hyperbolic NV equation: where from now on the indices will denote the partial derivatives w.r.t. the corresponding variables. This system allows for a Lax pair of the following type: If one knows two linearly independent solutions Ψ 1 (x, y, t) and Ψ 2 (x, y, t) for (16), then one can utilize the famous Moutard transformation to construct a new function Ψ[1](x, y, t) that will serve as a solution to the same equation (16) albeit with a new potential u[1](x, y, t). The new potential will then satisfy the relation Let us assume that u = a = b = 0. Then the entire system (16) simplifies to The equation (18) can be resolved by separating the variables. The resulting solution will be of a form: where A, B are two arbitrary functions that are continuously differentiable by x and y, correspondingly. Substituting (20) into (17) yields a following post-Moutard form of function u[1](x, y, t): As follows from (21), our next goal should lie in ascertaining the exact forms of the functions A(x, t) and B(y, t). This task can be accomplished by looking at the equation (19) which we have ignored so far. We will rewrite it as a standard Cauchy problem by introducing the initial conditions for A(t, x), B(t, y) and rewriting the (19) as a system where T = T (t) is an arbitrary time-dependent function. The apparently symmetric nature of (23) allows us to restrict our attention on just one of the equations therein, namely -the first one. We begin by introducing the Fourier transformÃ(p, t) of the function A(x, t): This transformation is handy because of the identity which, after being substituted into (23), yields the equation The equation (25) must be satisfied for all x and p, and therefore leads to: where C(p) is a function, determinable from the initial conditions (22). Using the inverse Fourier transform (24) we come to the following conclusion: According to (22), so the unknown C(p) is an inverse Fourier transform of the initial condition φ(x), i.e.: and we subsequently end up with the following general formula for the function A(x, t): The (29) can be further simplified by pointing out the similarity between the integrals with respect to variable p and the Airy function Ai(ξ). The Airy function is a particular solution of the eponymous Airy equation: that has a following integral representation: Using this fact together with the apparent identity: together with the equation (29) and the similar one written for B(y, t) finally yields: So, we end up with both the solution Ψ 1 = A + B of the Lax pair (18), (19), and, as a courtesy of Moutard transform (17), with a solution u[1] of the NV equation (15) as well. In other words, to find a non-zero solution of the NV equation, it will suffice to start with u ≡ 0, impose the boundary conditions (22) on the Lax pair (18), (19), use (32) to find its solution and conclude the calculations by finding a function u[1] via the Moutard transformation (17). As straightforward as it is, there is one question we should ask: what would happen should we try to invert the process and instead start out with he boundary conditions for the NV equation itself?

III. THE CAUCHY PROBLEM FOR THE NOVIKOV-VESELOV EQUATION
In the previous chapter we have shown that there shall exist a solution u[1](x, y, t) to the NV equation, whose exact form can be derived via the Moutard transformation (21) from the solutions of the system (18,19), provided we are given the initial conditions (22). But what would happen if the exact forms of the functions φ(x) and Φ(y) are unknown and we are instead given the initial conditions for the NV equation itself, and would it still be possible to find the required u [1]? In other words, is it possible to find an analytic solution to the Cauchy problem for the NV equation provided we only know that the solution has a general structure (21)? In short, the answer is "yes".
Let us start by introducing the set of initial and boundary conditions for the NV equation: where C ∈ R is some constant that is given to us together with the boundary conditions A 1 and B 1 . Since we know that u[1] satisfies the Moutard transformation, we also know that: where φ and Φ are defined as in Sec II, and and · denote the partial derivatives with respect to x and y variables correspondingly. From (34) and (33) it immediately follows that The first two differential equations in (35) can be easily integrated; for example, the first one after the integration with respect to the variable x yieldsΦ which leads us to the following conclusion: where we have introduced the notation: In a similar fashion, the boundary condition Φ(y) and its derivative will satisfy: (37) The system (36, 37) depends on four constants: φ 0 , Φ 0 , φ 0 andΦ 0 . Three of them can be chosen arbitrarily, whereas the fourth one would have to satisfy the equation (35), namely: Curiously, this choice does not affect the Cauchy problem of the NV equation in the least, for it can be shown by direct substitution into (34) that: i.e. the initial condition u 0 (x, y) depends only on the known initial boundary conditions A 1 (x), B 1 (y) and C.
We are now ready to answer the question posed in the beginning of this section: provided we know the initial conditions (38), how do we solve the corresponding Cauchy problem of the hyperbolic real-valued Novikov-Veselov equation? The answer lies in repeating the Moutard transformation process we described in Sec. II! Indeed, since the unknown functions φ(x) and Φ(y) satisfy the relations (36) and (37), all we really have to do is substitute them into the system (32), derive A(x, t) and B(y, t), and substitute them in equation (21) to find out the sought after u[1](x, y, t), which will conclude the problem.
Lets summarize everything we have said so far. In order to find an exact solution u(x, y, t) to the hyperbolic real-valued Novikov-Veselov equation with the given initial boundary conditions u(x, 0, 0) = A 1 (x), u(0, y, 0) = B 1 (y), u(0, 0, 0) = C, that correspond to the initial condition one shall: 1. Choose a differentiable function T (t) and four numbers α, β, γ and δ that satisfy the condition, 2. Find two support function φ(x) and Φ(y) via the formulas

Substitute φ(x) and Φ(y) into the equations
4. Substitute the new functions A(x, t) and B(y, t) into the equation The resulting function u(x, y, t) will be a proper solution of the Cauchy problem since by construction it will satisfy both the NV equation (39), and the initial conditions u(x, y, 0) = u 0 (x, y). We would like to emphasize here that this procedure does not involve anything more complicated than partial differentiation and integration and can therefore be used for both the analytic study of the properties of the solutions of NV equation and the corresponding numerical calculations.
Before we conclude this section, we would like to offer two interesting examples of Cauchy problems that might be used for the algorithm we have described above and that serve as the proof that even the seemingly simple case of (18) with u = 0 can lead to some rather interesting problems.
The first is based on the solution of (20) of the form: Ψ 1 (x, y, 0) = cosh(x − x 0 ) + cosh(y − y 0 ), where x 0 , y 0 = 0. According to (21) it corresponds to the following dromion solution: depicted on Fig. 1. It is important to note that x 0 and y 0 must be non-zero constants, otherwise the functions φ(x) = u 0 (x, 0) and Φ(y) = u 0 (0, y) in (40) will be identically zero. Any other choice for x 0 , y 0 ∈ R, however, would be fine and will result in nontrivial solutions of the NV equation. Also, a little note is in order. The solution we have just produced is the exponentially exponentially localized soliton localized on the 2D plane. The solutions of this have been previously constructed for the Davey-Stewartson-I (DS) equation in [32], [33] and [34], while the term "dromion" itself stems from the 1989 paper by Fokas ans Santini. [35]. The DS equations describe two interacting fields, and in the case of the dromion on of them describes a certain exponentially localized (on a plane) structure, while the other has orthogonal equipotential lines. It is thanks to this very property that we are at liberty to call the solution (42) a "dromion". The second solution is based on the function: Ψ 1 (x, y, 0) = e (x−x0) 2 + e −(y−y0) 2 , where x 0 , y 0 = 0. The solution u 0 (x, y) then has the form: depicted on Fig. 2, and it is easy to see that on the plane 0xy it describes an exponentially localized structure, centered around the point (x 0 , y 0 ).

IV. GENERALIZATION OF THE METHOD: THE HIGHER-ORDER EQUATIONS
Let us now say a few words about the more general problem. First of all, let us look at the following operator-type Lax pair: where n ∈ N + in some non-zero natural number. This system will correspond to a family of Lax equations, with the special case n = 3 corresponding to the hyperbolic NV equation. It will still allow for the Moutard transformation, and therefore the crux of our discussion would still be applicable for arbitrary n. However, one thing that must change is the exact form of the equations for A(x, t) and B(y, t). This is due to a very simple reason: the Fourier transform A(p, t) for the function A(x, t) had to satisfy the equation (26), where one of the terms had a specific factor ip 3 . It is this very factor that had allowed us to invoke the Airy function and with its aid derive the exact formula for A(x, t).
Unfortunately, this very factor is endemic to the Veselov-Novikov equation, and in the more general case of (44) it must be replaced with the factor −(ip) n -destroying any relationship with the Airy function proper. As a result, the entire formula (40) becomes no longer applicable for the general case and thus should be properly replaced. In order to find out the suitable replacement, we shall separately consider two alternative cases: when n is odd and when n is even.
Case 1: Odd n. Let n = 2m + 1, where m ≥ 0. Following our previous discussion, let us consider the special case u ≡ 0. Then the system (44) turns into Since the first equation requires that Ψ = A(x, t)+B(y, t), the system subsequently splits into the following equations: where for simplicity we have omitted the arbitrary function T (t). Using the Fourier transformatioñ we end up with the differential equation Solving it and returning back to A(x, t) as described in Sec.II yields As we know, in the special case m = 1 (i.e. n = 3) the inner integral in (48) can be rewritten in terms of the Airy function which serves as a solution to the Airy equation and is easily derived using either Fourier or Laplace transformation; in case of the Laplace transformation the contour of integration must be chosen lying inside of a sector where the polar angle θ satisfies the condition cos(3θ) > 0. Similarly, it is easy to show that one of a solutions to a more general equation will be a higher-order generalization of the Airy function: which means that the required functions A 2m+1 and B 2m+1 can be derived from the initial conditions φ(x) and Φ(y) by the following formulas: SIDE NOTE. We would like to remind the reader that in literature the term generalized Airy function is commonly assigned to the solutions of the second order O.D.E. w (x) = x n w(x); hence the addition of the term -order in our case is necessary to avoid a possible confusion.
Case 2: Even n. Let n = 2m, where m ≥ 0. This time let us utilize not a Fourier but a Laplace transform: where we have introduced the equation forÃ(p, t) is so the required function A(x, t) will satisfy the equation It is not difficult to show that the Laplace transformation method applied to the ordinary differential equation will yield a following solution and so the even case produces the formulas that are quite similar to the old ones, namely: Note the appearance of a negative sign under the root in (50), which serves as a indication of an ill-posedness of our problem for t > 0. SIDE NOTE The choice of the lower limit in the Laplace transform being equal to −∞ is neither random nor capricious -it is necessary so that during the integration by parts (which is required to get the solution (49)) all the boundary terms will vanish.
In order to conclude this section, we have to produce a sample of a differential equation that allows for a Moutard transformation and whose Lax pair reduces to (44) when u = 0. It is possible to demonstrate that no such equation exists for n = 4, which gives a very strong indication that the method is intricately tied up to the Novikov-Veselov hierarchy. And indeed, if we consider the following Lax pair: Ψ xy = −uΨ Ψ t = Ψ 5x + Ψ 5y + 5 (ω yy Ψ yy ) y + 5 (ω xx Ψ xx ) x + aΨ y + bΨ x , with the functions a, b and ω satisfying the conditions: ω xy = u a x = 5u yyy + 5 (uω yy ) y + 10u y ω yy b y = 5u xxx + 5 (uω xx ) x + 10u x ω xx , then we end up with the second equation from the Novikov-Veselov hierarchy, which has the following form: u t = u 5x + u 5y + 5 (ω yy u y ) yy + (ω xx u x ) xx + (au) y + (bu) x , which exactly satisfies both of our underlying assumptions whenever we set ω = 0, or, more generally, choose ω = c 1 x + c 2 y + c 3 , where c 1 , c 2 , c 3 ∈ R. Subsequent application of the proposed method to this equation (plus the initial conditions) for the case u = 0 we leave to the reader.

V. CONCLUSION
At last, let us reflect on the results we have gained. The main aim of the article was to demonstrate that the Moutard symmetry is not only well-suited to construct explicit partial solutions of nonlinear partial differential equations such as Novikov-Veselov equation, but it is also versatile enough to solve the Cauchy problems like (33). This is a very important and useful property of the Moutard symmetry, because the common way to solve the Cauchy problem invokes the inverse scattering method whose applications for the (1 + 2)-dimensional equations remains far from being completely understood. Another interesting application of the Moutard symmetry may be connected with the construction of so-called dressing chains of discrete symmetries. Such chains connect different integrable hierarchies and in the long run suggest that potentially all the nonlinear integrable equations might be different manifestations of one unique integrable differential equation, but simply written in different calibrations (see, for example, [36]).
Another interesting open question is the evolution of initially exponentially localized structures. There is a set of theorems that explicitly forbids the existence of exponentially localized solitons of NV equation (see the discussion in Sec. I), so we shall expect that the evolution of such structures as (42) and (43) will destroy the localization and transform the solution into some other type of structure. Our approach opens the window of opportunity to study these complex processes analytically.