Reconstruction of higher-order differential operators by their spectral data

This paper is concerned with inverse spectral problems for higher-order ($n>2$) ordinary differential operators. We develop an approach to the reconstruction from the spectral data for a wide range of differential operators with either regular or distribution coefficients. Our approach is based on the reduction of an inverse problem to a linear equation in the Banach space of bounded infinite sequences. This equation is derived in a general form that can be applied to various classes of differential operators. The unique solvability of the linear main equation is also proved. By using the solution of the main equation, we derive reconstruction formulas for the differential expression coefficients in the form of series and prove the convergence of these series for several classes of operators. The results of this paper can be used for constructive solution of inverse spectral problems and for investigation of their solvability and stability.


Historical background
Inverse problems of spectral analysis consist in the recovery of differential operators from their spectral information. Such problems arise in practice when one needs to determine certain physical parameters of a system from some measured data or to construct a model with desired properties. The majority of physical applications are concerned with linear differential operators of form (1.1) with n = 2, 3, 4. For n = 2, expression (1.1) turns into the Sturm-Liouville (Schrödinger) operator − ℓ 2 (y) = −y ′′ + q(x)y, (1.2) which models string vibrations in classical mechanics, electron motion in quantum mechanics, and is widely used in other branches of science and engineering. The third-order linear differential operators arise in the inverse problem method for integration of the nonlinear Boussinesq equation (see [10,11]), in mechanical problems of modeling thin membrane flow of viscous liquid and elastic beam vibrations (see [12] and references therein). Inverse spectral problems for the fourth-order linear differential operators attract much attention of scholars because of applications in mechanics and geophysics (see [13][14][15][16][17][18][19][20] and references therein).
The classical results of the inverse problem theory have been obtained for the Sturm-Liouville operator (1.2) with integrable potential q(x) in 1950th by Marchenko, Levitan, and their followers (see [21,22]). They have developed the transformation operator method, which reduces the nonlinear inverse Sturm-Liouville spectral problem to the linear Fredholm integral equation of the second kind. However, the transformation operator method appeared to be ineffective for the higher-order differential operators y (n) + n−2 k=0 p k (x)y (k) , n > 2. (1.3) Note that the differential expression (1.1) can be transformed into (1.3) in the case of sufficiently smooth coefficients {τ ν } n−2 ν=0 . Thus, the development of inverse spectral theory for the higher-order operators (1.3) required new approaches. Relying on the ideas of Leibenson [23,24], Yurko has created the method of spectral mappings. This method allowed him to construct inverse problem solutions for the higher-order differential operators (1.3) with regular (integrable) coefficients on the half-line x > 0 and on a finite interval x ∈ (0, T ) (see [25,26]). The case of Bessel-type singularities also was considered [27,28]. Later on, the ideas of the method of spectral mappings were applied to a wide range of inverse spectral problems, e.g., to inverse problems for the first-order differential systems [29], for differential operators on graphs [30], for quadratic differential pencils [31]. This method is based on the theory of analytic functions and mainly on the contour integration in the complex plane of the spectral parameter. The method of spectral mappings reduces a nonlinear inverse problem to a linear equation in a suitable Banach space. This space is constructed in different ways for different operator classes. In particular, for differential operators on a finite interval, the main equation is usually derived in the space m of infinite bounded sequences. It is also worth mentioning that an approach to inverse scattering problems for higher-order differential operators (1.3) on the full line was developed by Beals et al [32,33].
During the last 20 years, the inverse problems are actively investigated for the secondorder differential operators with distributional potentials (see, e.g., [34][35][36][37][38][39][40][41][42][43]). In particular, Hryniv and Mykytyuk [34][35][36] transferred the transformation operator method to the Sturm-Liouville operators (1.2) with potential q(x) of class W −1 2 (0, 1) and so generalized the basic results of inverse problem theory to this class of operators. Note that the space W −1 2 contains the Dirac δ-function and the Coulumb potential 1 x , which are used for modeling particle interactions in quantum mechanics [44]. The method of spectral mappings has been extended to the Sturm-Liouville operators with potentials of W −1 2 in [37,43,45]. This opens the possibility of constructing the inverse spectral theory for higher-order differential operators with distribution coefficients. However, till now, only the first steps have been taken in this direction. In [9,46], the uniqueness of recovering the higher-order differential operators with distribution coefficients on a finite interval and on the half-line has been studied. The goals of this paper are to derive the linear main equation of the inverse problem, to prove its unique solvability, and to obtain reconstruction formulas for the coefficients {τ ν } n−2 ν=0 of various classes.

Problem statement and methods
Our treatment of the differential expression (1.1) is based on the regularization approach. Namely, we will assume that the differential equation ℓ n (y) = λy, x ∈ (0, 1), (1.4) where λ is the spectral parameter, can be equivalently transformed into the first-order system Y ′ (x) = (F (x) + Λ)Y (x), x ∈ (0, 1), (1.5) where Y (x) is a column vector function of size n, Λ is the (n × n)-matrix whose entry at the position (n, 1) equals λ and all the other entries are zero, and F (x) = [f k,j (x)] n k,j=1 is a matrix function with the following properties: We denote the class of (n × n) matrix functions satisfying (1.6) by F n . By using any matrix F ∈ F n , one can define the quasi-derivatives and the domain Definition 1.1. A matrix function F (x) ∈ F n is called an associated matrix of the differential expression ℓ n (y) if ℓ n (y) = y [n] for any y ∈ D F . We call a function y a solution of equation (1.4) if y ∈ D F and y [n] = λy, x ∈ (0, 1).
In this paper, we assume that ℓ n (y) is any differential expression that has an associated matrix in terms of Definition 1.1. We do not impose any additional restrictions on {τ ν } n−2 ν=0 , since we are interested to formulate the abstract results which can be applied to various classes of differential operators. Certain restrictions on {τ ν } n−2 ν=0 will be imposed below when necessary. Let us proceed to the inverse problem formulation. Suppose that we have a differential expression of form (1.1) and an associated matrix F (x) = [f k,j ] n k,j=1 . By using the corresponding quasi-derivatives (1.7), define the linear forms where p s,a ∈ {0, . . . , n − 1}, p s,a = p k,a for s = k, and u s,j,a are some complex numbers. In addition, introduce the matrices U a = [u s,j,a ] n s,j=1 , u s,j,a := δ j,ps,a+1 for j > p s,a , a = 0, 1. Here and below, δ j,k is the Kronecker delta. We call by the problem L the triple (F (x), U 0 , U 1 ). Below we introduce various characteristics related to the problem L.
It has been proved in [9,Section 4] that, for all λ ∈ C except for a countable set, equation (1.4) has the so-called Weyl solutions {Φ k (x, λ)} n k=1 satisfying the boundary conditions Define the matrix function Φ(x, λ) = [ Φ k (x, λ)] n k=1 . The columns of the matrices C(x, λ) and Φ(x, λ) form fundamental solution systems of (1.5). Consequently, the following relation holds: where the matrix function M(λ) is called the Weyl matrix of the problem L (see [9]). The notion of Weyl matrix generalizes the notion of Weyl function for the second-order operators (see [21,26]). Weyl functions and their generalizations play an important role in the inverse spectral theory for various classes of differential operators. In particular, Yurko [25][26][27][28] has used the Weyl matrix as the main spectral characteristics for the reconstruction of the higher-order differential operators (1.3) with regular coefficients. The analogous inverse problem for the differential expression of form (1.1) can be formulated as follows. The uniqueness of Problem 1.2 solution has been proved in [9] for the Mirzoev-Shkalikov case: n = 2m, τ 2k+j ∈ W −(m−k−j) 2 (0, 1) and n = 2m + 1, τ 2k+j ∈ W −(m−k−j) 1 (0, 1), j = 0, 1. In [46], the uniqueness of recovering the boundary condition coefficients from the Weyl matrix has been studied.
It has been shown in [9,Section 4] that the Weyl matrix M(λ) = [M j,k (λ)] n j,k=1 is unit lower-triangular, and its non-trivial entries have the form where ∆ k,k (λ) := det[U s,1 (C r )] n s,r=k+1 and ∆ j,k (λ) is obtained from ∆ k,k (λ) by the replacement of C j by C k . The functions C r (1, λ), r = 1, n, s = 0, n − 1, are entire analytic in λ, so do the functions ∆ j,k (λ), 1 ≤ k ≤ j ≤ n. Hence, M(λ) is meromorphic in λ, and the poles of the k-th column of M(λ) coincide with the zeros of ∆ k,k (λ). At the same time, the zeros of the entire functions ∆ j,k (λ) , 1 ≤ k ≤ j ≤ n, coincide with the eigenvalues of some boundary value problems for equation (1.4), and the inverse problem by the Weyl matrix (Problem 1.2) is related to the inverse problem by n(n+1) 2 spectra (see [9] for details). We will say that the problem L belongs to the class W if all the zeros of ∆ k,k (λ) are simple for k = 1, n − 1. Then, in view of (1.12), the poles of M(λ) are simple. In general, the function ∆ k,k (λ) can have at most finite number of multiple zeros. The latter case can be treated by developing the methods of Buterin et al [52,53], who considered the non-self-adjoint Sturm-Liouville operators (n = 2) with regular potentials. However, the case of multiple zeros is much more technically complicated, so, in this paper, we always assume that L ∈ W .
Denote by Λ the set of the Weyl matrix poles. Consider the Laurent series We call the collection {λ 0 , N (λ 0 )} λ 0 ∈Λ the spectral data of the problem L. Obviously, the spectral data are uniquely specified by the Weyl matrix M(λ), so Problem 1.2 can be reduced to the following problem. Problem 1.3. Given the spectral data {λ 0 , N (λ 0 )} λ 0 ∈Λ , find the coefficients {τ ν } n−2 ν=0 . It is more convenient to study the reconstruction question for Problem 1.3. It is worth mentioning that, in fact, the Weyl matrix and the spectral data can be constructed according to the above definitions for any matrix function F (x) of class F n , not necessarily associated with any differential expression of form (1.1). But, in general, the matrix F (x) is not uniquely specified by the Weyl matrix (see Example 4.5 in [46]). Therefore, in this paper, the solution of Problem 1.3 is divided into the two steps: ν=0 . The recovery of the Weyl solutions {Φ k (x, λ)} n k=1 from the spectral data is studied for a matrix F (x) of general form, and then reconstruction formulas are derived for {τ ν } n−2 ν=0 of certain classes.
For a fixed F ∈ F n , we define the quasi-derivatives (1.7), the expression ℓ n (y) := y [n] , the problem L = (F (x), U 0 , U 1 ), its spectral data {λ 0 , N (λ 0 )} λ 0 ∈Λ as above, and focus on the following auxiliary problem. Problem 1.4. Given the spectral data {λ 0 , N (λ 0 )} λ 0 ∈Λ , find the Weyl solutions {Φ k (x, λ)} n k=1 . Let us briefly describe the method of solution. Along with L, we consider another problem L = (F (x),Ũ 0 ,Ũ 1 ) of the same form but with different coefficients. Similarly to Φ(x, λ), definẽ Φ(x, λ) forL. An important role in our analysis is played by the matrix of spectral mappings: For each fixed x ∈ [0, 1], the matrix function P(x, λ) is meromorpic in λ with poles at the eigenvalues Λ ∪Λ. The method is based on the integration of some functions by a special family of contours enclosing these eigenvalues. Applying the Residue theorem, we derive an infinite system of linear equations. Further, that system is transformed into a linear equation in the Banach space m of infinite bounded sequences. The main equation of the inverse problem has the form where, for each fixed x ∈ [0, 1], ψ(x) andψ(x) are elements of m,R(x) is a linear compact operator in m, and I is the unit operator. The elementψ(x) and the operatorR(x) are constructed by the model problemL and by the spectral data of the two problems L,L, respectively, while the unknown element ψ(x) is related to the desired functions {Φ k (x, λ)} n k=1 . We prove that the operator (I −R(x)) has the bounded inverse, and so the main equation is uniquely solvable (see Theorem 3.6). This implies the uniqueness of solution for Problem 1.4. Using the main equation, we obtain a constructive procedure for solving Problem 1.4 (see Algorithm 3.7). These results can be applied to a wide range of differential operators (1.1) with associated matrices of class F n .
The reconstruction formulas have the form of series, and the main difficulties in our analysis are related to studying the convergence of those series. These difficulties increase for the case of non-smooth and/or distribution coefficients. In order to prove the series convergence, we use the Birkhoff-type solutions constructed by Savchuk and Shkalikov [2] and the precise asymptotic formulas for the spectral data obtained in [54]. For the cases (ii) and (iii), we reconstruct the functions τ ν step-by-step for ν = n − 2, n − 3, . . . , 1, 0. The similar approach can be used in the case of odd n, which requires technical modifications.
By using the reconstruction formulas, one can develop numerical methods for solving inverse spectral problems (see [55] for the second-order case). However, this issue requires an additional work. In this paper, we obtain theoretical algorithms, which in the future can be used for investigation of existence and stability of the inverse problem solution.
It is worth mentioning that our method of inverse problem solution is the first one for higher-order differential operators with distribution coefficents. The obtained main equation and reconstruction formulas generalize the results of [45] for the Sturm-Liouville operators with distribution potential. The other methods applied to the second-order operators (see, e.g., [34,39]), to the best of the author's knowledge, appear to be ineffective for higher orders.
The paper is organized as follows. In Section 2, we provide preliminaries and study the properties of the spectral data. Section 3 is devoted to the contour integration and to the derivation of the main equation of the inverse problem in a Banach space. The unique solvability of the main equation is also proved. As a result, an algorithm for solving the auxiliary Problem 1.4 is obtained for arbitrary F ∈ F n . In Section 4, we derive the reconstruction formulas for the coefficients {τ ν } n−2 ν=0 and study the convergence of the obtained series. Section 5 contains a brief summary of the main results.

Preliminaries
Throughout the paper, we use the following notations.
1. I is the (n × n) unit matrix, e k is the k-th column of I, k = 1, n.
2. The sign T denotes the matrix transpose.
6. The notations ⌊x⌋ and ⌈x⌉ are used for rounding a real number x down and up, respectively.
7. The binomial coefficients are denoted by C k n = n! k!(n − k)! .
8. Along with L, we will consider the problemsL, L ⋆ ,L ⋆ of the same form but with different coefficients. We agree that, if a symbol γ denotes an object related to L, then the symbols γ, γ ⋆ ,γ ⋆ will denote the analogous objects related toL, L ⋆ ,L ⋆ , respectively. Note that the quasi-derivatives for the problemsL, L ⋆ ,L ⋆ are defined by using the matricesF (x), , respectively, which may be different from F (x).

The notation y [k]
is used for quasi-derivatives defined by (1.7) (or analogously by using the entries ofF (x), F ⋆ (x), orF ⋆ (x)). The notation y(x) is used for the column vector of the quasi-derivatives 10. In estimates, the symbol C is used for various positive constants independent of x, l, k, etc.

Problems
as follows: Let F (x) be a fixed matrix function of class F n . Suppose that y ∈ D F and z ∈ D F ⋆ , the quasi-derivatives for y are defined via (1.7) by using the elements of F (x), the quasi-derivatives for z are defined as and Lemma 2.1. The following relation holds: From (2.2) and (1.7), we obtain Substituting the latter relations into (2.4), we get Note that Taking (2.1) into account, we arrive at (2.3).
If y and z satisfy the relations ℓ n (y) = λy and ℓ ⋆ n (z) = µz, respectively, then ) by using the corresponding quasi-derivatives (1.7) and (2.2), and the matrix J : For a = 0, 1, let U a = [u s,j,a ] n s,j=1 be an (n × n) matrix such that u s,j,a = δ j,ps,a+1 for j > p s,a , where p s,a ∈ {0, . . . , n − 1}, p s,a = p k,a for s = k. The matrices U a define the linear forms U s,a via (1.8).
Along with U a , consider the matrices The matrices U ⋆ a are chosen is such a way that the following relation holds: for any y ∈ D F , z ∈ D F ⋆ . Indeed, the right-hand side of (2.8) can be represented in the matrix form [U ⋆ a z(a)] T J a U a y(a), Taking (2.6) and (2.7) into account, we arrive at (2.8).
By virtue of Theorem 1.1 in [54], the spectrum of L k is a countable set of eigenvalues Λ k := {λ l,k } l≥1 having the following asymptotics (counting with multiplicities): (2.14) where {κ l,k } ∈ l 2 and χ k are constants which depend only on n, k, and {p s,a }. Hence, for a fixed k ∈ {1, . . . , n − 1} and sufficiently large l, the eigenvalues λ l,k are simple. Assume that L ∈ W , that is, all the zeros of ∆ k,k (λ) are simple for k = 1, n − 1. Then, in view of (1.12) and (2.11), the poles of M(λ) and M ⋆ (λ) are simple. It follows from (1.11) and (2.10) that the matrix functions Φ(x, λ) and Φ ⋆ (x, λ) for each fixed x ∈ [0, 1] also have only simple poles.
Proof. This lemma is proved similarly to Lemma 2.3.1 in [26], so we outline the proof briefly. If λ 0 ∈ Λ k , then Φ k, −1 (x, λ 0 ) = 0. On the other hand, it follows from (2.17) that Applying the linear forms U s,0 to this relation for s = k + 1, n, we conclude that N s,k (λ 0 ) = 0, s = k + 1, n. Thus, the assertion (i) is proved for j = k. The proof for j = k − 1, . . . , 2, 1 can be obtained by induction.
In order to prove (ii), we suppose that In view of the asymptotics (2.14), we have λ l,k = λ r,k+1 for sufficiently large l and r. Therefore, Lemma 2.4 implies the following corollary.
Define the weight numbers β l,k := N k+1,k (λ l,k ). It is worth considering β l,k only for sufficiently large l. It follows from (1.13) and (1.12) that .
Consequently, Theorem 6.2 from [54] yields the asymptotics where the constants β 0 k depend only on n, k, and {p s,a }.

Main equation
This section is devoted to the constructive solution of the auxiliary Problem 1.4, that is, to the recovery of the Weyl solutions {Φ k (x, λ)} n k=1 from the spectral data {λ 0 , N (λ 0 )} λ 0 ∈Λ . We consider this problem for L = (F (x), U 0 , U 1 ) ∈ W with an arbitrary F ∈ F n . Thus, the results of this section can be applied to a wide class of differential expressions (1.1) with associated matrix of F n .
Along with L, we consider another problemL = (F (x),Ũ 0 ,Ũ 1 ) of the same form but with different coefficients. Assume thatF ∈ F n , p s,a =p s,a , s = 1, n, a = 0, 1. The quasi-derivatives forL are defined by the matrixF (x), so they are different from the quasi-derivatives of the problem L. The problemL ⋆ is defined similarly to L ⋆ . For simplicity, we assume thatL ∈ W . The caseL ∈ W requires technical modifications (see Remark 3.4). Denote I := Λ ∪Λ.
In Subsection 3.1, we reduce the studied problem to the infinite system (3.31) of linear equations with respect to some entries of φ 0 (x, λ 0 ), λ 0 ∈ I. Our technique is based on the contour integration in the λ-plane and on the Residue theorem. In Subsection 3.2, the system (3.31) is transformed into the main equation (3.43) in the Banach space m of infinite bounded sequences. The unique solvability of the main equation is proved. Finally, we arrive at the constructive Algorithm 3.7 for finding {Φ k (x, λ)} n k=1 by the spectral data. This algorithm will be used in the next section for solving the inverse spectral problem.

Contour integration
In order to formulate and prove the main lemma of this subsection (Lemma 3.2), we first need some preliminaries. Introduce the notations and similarly defineD α (x, λ 0 , λ).
Lemma 3.2. The following relations hold: where the series converge in the sense uniformly by λ, µ on compact sets of (C \ I).
Proof. In this proof, a crucial role is played by the matrix of spectral mappings It follows from (2.12) and (3.13) that 14) The proof consists of three steps.
Step 2. Contour integration. In view of (3.14), the matrix function P(x, λ) is meromorpic in λ with the poles I. Hence, P(x, λ) is analytic in Ξ ± R . Let P 1 (x, λ) be the first row of P(x, λ). The Cauchy formula implies Consequently, dξ. where It follows from (3.16) that Step 3. Residues. Using the first row of (3.14): and the Residue theorem, we obtain 1 2πi  Using analytic continuation, we conclude that these relations hold for λ, µ ∈ (C \ I).
The relations (3.31) can be considered as an infinite linear system with respect to ϕ l,k,ε (x), (l, k, ε) ∈ V . However, it is inconvenient to use (3.31) as the main equation system for the inverse problem, because the series in (3.31) converges only "with brackets": Therefore, in the next section, we transform the system (3.31) to a linear equation in a suitable Banach space. The relation (3.32) will be used to prove the unique solvability of the main equation.
Remark 3.4. IfL ∈ W , that is, the poles ofM(λ) are not necessarily simple, then this influences on the calculation of the residues in (3.22). Consequently, we obtain the following relation instead of (3.11): where m λ 0 is the multiplicity of λ 0 ∈Λ. Using (3.33), one can derive an infinite system analogous to (3.31), containing not only entries of the vectors φ 0 (x, λ 0 ) but also of φ k (x, λ 0 ) for k = 1, m λ 0 − 1.
Thus, we arrive at the following algorithm for solving Problem 1.4.

Reconstruction formulas
In this section, we use the solution ψ(x) of the main equation (3.43) to obtain the solution of Problem 1.3 for some classes of differential operators. We derive the reconstruction formulas in the form of series for the coefficients {τ ν } n−2 ν=0 of the differential expression (1.1). In Subsection 4.1, the general approach to obtaining reconstruction formulas is described. However, for certain classes of the coefficients {τ ν } n−2 ν=0 , the convergence of the obtained series has to be studied in the corresponding spaces. Therefore, in Subsection 4.2, we prove an auxiliary lemma on the series convergence. In Subsections 4.3-4.5, we study the three classes of operators: (i) n = 3, τ 0 ∈ W −1 2 (0, 1), τ 1 ∈ L 2 (0, 1); (ii) n is even, τ ν ∈ L 2 (0, 1), ν = 0, n − 2; (iii) n is even, τ ν ∈ W −1 2 (0, 1), ν = 0, n − 2. For each case, we provide the uniqueness theorem of the inverse problem solution in an appropriate statement, obtain reconstruction formulas and prove the convergence of the series, and so get constructive algorithms for solving Problem 1.3. For the cases (ii) and (iii), we recover the coefficients τ n−2 , τ n−3 , . . . , τ 1 , τ 0 one-by-one in order to achieve the convergence estimates for the corresponding series. The even order in (ii) and (iii) is considered for definiteness. The similar ideas can be applied to the odd-order differential operators. For simplicity, in all the three cases, we choose such boundary conditions that their coefficients cannot be uniquely recovered from the spectral data and so do not consider their reconstruction. However, for other types of boundary conditions, the recovery of their coefficients also can be studied similarly to the regular case (see Lemma 2.3.7 in [26]).
In this section, we use the following notations for an index v = (l, k, ε) ∈ V ′ : Additionally, define the scalar functions

Series convergence
In this subsection, we prove the following auxiliary lemma.
Here and below, the quasi-derivatives for φ v (x) are generated by the matrix F (x) and for g v (x), byF ⋆ (x). In order to prove Lemma 4.1, we need to formulate preliminary propositions.
Using the arguments above, we obtain the estimate |Z l,k (x)| ≤ Cl j 1 +j 2 ξ l .
Hence, in the case {l j 1 +j 2 ξ l } ∈ l 1 , the series l≥l 0 Z l,k (x) converges absolutely and uniformly with respect to x ∈ [0, 1]. Taking (4.19) into account, we arrive at the assertion of the lemma.
Note that the algorithms of this section are valid forL ∈ W . However, the caseL ∈ W requires only technical modifications due to Remark 3.4, which do not influence on the convergence of the series.

Conclusion
Let us briefly summarize the results of this paper. We have studied the inverse spectral problem which consists in recovering the coefficients {τ ν } n−2 ν=0 of the differential expression (1.1) from the spectral data {λ 0 , N (λ 0 )} λ 0 ∈Λ . An approach to constructive solution of the inverse problem is developed. Our approach can be applied to a wide class of differential expressions ℓ n (y), which admit regularization in terms of associated matrix. The inverse problem solution consists of the two steps. First, we consider the auxiliary problem of finding the Weyl solutions {Φ k (x, λ)} n k=1 by using the spectral data. This problem is reduced to the linear equation (3.43) in the Banach space m of bounded infinite sequences. Theorem 3.6 on the unique solvability of the main equation (3.43) is proved. Second, by using the solution of the main equation, we derive reconstruction formulas for the coefficients {τ ν } n−2