Semiclassical Spectral Series Localized on a Curve for the Gross–Pitaevskii Equation with a Nonlocal Interaction

: We propose the approach to constructing semiclassical spectral series for the generalized multidimensional stationary Gross–Pitaevskii equation with a nonlocal interaction term. The eigenvalues and eigenfunctions semiclassically concentrated on a curve are obtained. The curve is described by the dynamic system of moments of solutions to the nonlocal Gross–Pitaevskii equation. We solve the eigenvalue problem for the nonlocal stationary Gross–Pitaevskii equation basing on the semiclassical asymptotics found for the Cauchy problem of the parametric family of linear equations associated with the time-dependent Gross–Pitaevskii equation in the space of extended dimension. The approach proposed uses symmetries of equations in the space of extended dimension.


Introduction
The collective modes of coherent quantum atom ensembles in the Bose-Einstein condensate (BEC) confined by a trap field in the mean-field approximation are described by the nonlinear Schrödinger equation, the Gross-Pitaevskii equation (GPE) (see, e.g., [1][2][3][4]): Here, ∂ t = ∂ ∂t , ∆ is the Laplace operator, m is the particle mass, κ is the nonlinearity parameter, V tr ( x) is the external field (trap) potential, the function Ψ( x, t) is the macroscopic wavefunction of the condensate (the order parameter) such that |Ψ( x, t)| 2 is proportional to the BEC density and its phase gradient is proportional to the BEC velocity. For the sake of simplicity, we omit the argumenth of the function Ψ( x, t,h) in formulas where it does not cause the confusion.The nonlinear interaction operator W[Ψ]( x, t) in (1) is usually chosen as and it describes the local interaction with g = 4πh 2 a/m where a is the s-wave scattering length. However, to take into account long-range interaction, the nonlocal operators are considered here W( x, y)|Ψ( y, t)| 2 d y.
The local operator W local of (2) can be treated as the limit of the nonlocal version (3) with a delta-like function W( x, y) that corresponds to the short-range interaction. Some applications in physics, specifically the BEC with the dipolar-to-dipolar interaction [5][6][7][8], requiring accounting the long-range interaction do not allow one to get rid of nonlocality for the sake of simplicity. In these cases, the combination of the local and nonlocal terms is often used. We call Equation (1) with the term W local the local GPE, and with W nonlocal the nonlocal GPE.
Substitution of Ψ( which describes the stationary states of the BEC. If solutions ψ( x) of (4) are restricted to the unit norm ||ψ( then eigenvalues µ of the nonlinear eigenvalue problem (4) and (5) correspond to the chemical potential and eigenfunctions ψ( x) minimize the energy functional where E is the energy per the BEC particle. The ground state solution corresponds to the global minimum of the functional (6) while the exited states meet the local minima [9,10]. For the weak interaction (a small parameter κ), the stationary solutions can be found using the perturbation theory as a correction to solutions of the so-called associated linear Schödinger equation that inherits the GPE linear part and has not the interaction term, i.e., it is the Equation (4) with κ = 0. Solutions of the GPE obtained in this way are termed solutions with linear counterpart in [11]. The opposite limit of strong nonlinearity was also considered using the Thomas-Fermi approximation when the kinetic term can be ignored and the local GPE (4) takes the form of an algebraic equation. For the nonlocal case, the GPE (4) becomes an integral equation. Moreover, it was shown that there are stationary solutions of the local GPE that do not correspond to any eigenfunctions of the associated linear Schödinger equation for multiwell external potential in [12] (the GPE solutions without linear counterpart).
The minima of the functional (6) can be found numerically by iterative procedure (see, e.g., the review in [13] and references therein). The issue is that the iterative procedure needs a proper initial approximation for the respective stationary states (4). For the weak nonlinearity, the stationary states of the associated linear Schrödinger equation can be used as the initial data, and, for the strong nonlinearity, the initial data can be obtained from the Thomas-Fermi approximation [9,14,15]. Numerical methods are usually limited to calculating only the ground state [13] and do not allow one to construct spectral series in the general case.
Note that the solutions with linear counterpart inherit the symmetries of the external field potential while these symmetries are broken for solutions without linear counterpart.
In order to avoid the limitations of the weak or strong nonlinearity, the semiclassical approximation can be used. For example, in [21], the semiclassical approximation was used to study the non-stationary local GPE with the special type of external field and different orders of the nonlinear term.
Here, we propose an approach to the spectral problem (4) that deals with the arbitrary strength of the nonlinearity and is applicable for a wide range of the external fields in n-dimensional nonlocal GPE. Moreover, we study the nonlocal GPE solutions localized on curves reflecting the BEC topological features. We apply the well-known WKB-Maslov theory of semiclassical approximation [22][23][24] to the problem under consideration. The approach to the eigenvalue problem (4) and (5) is based on [25], where the leading term of semiclassical asymptotics for the Cauchy problem has been constructed explicitly in a class of functions localized in a neighborhood of a curve in the phase space of a dynamic system of moments of the nonlocal GPE solution termed the Hamilton-Ehrenfest system. The last one for the n-dimensional manifolds in the phase space was studied in [26][27][28]. As opposed to the linear case (κ = 0 in (4)), the variance matrix of the solutions constructed in [25] may be bounded due to the nonlinearity. Such solutions of the Cauchy problems are termed the semiclassical soliton-type solutions [29], and we use exactly such particular solutions for the construction of the spectral series of the nonlocal GPE.
Differing from the associated linear Schödinger equation in the sense of [11], our method uses a parametric family of the associated linear Gross-Pitaevskii equations that include terms corresponding to both the external field and the interaction. Among the solutions of these equations, the solution of the initial nonlinear GPE (4) can be found.
Our consideration is not limited to the convolution nonlinearity with the kernel W( x, y) = W( x − y), which naturally arise in the BEC application. Although the convolution nonlinearity is usually simpler for calculations, it is unsubstantial for our approach. Moreover, the consideration of the wider class of kernels W( x, y) allows one to apply this approach to the GPE with the modulated nonlinearity that was studied in [18,19].
The WKB-Maslov method sometimes leads to quite complex classical equations (the Hamilton-Ehrenfest system) that need semi-analytical considerations. However, the development of numerical-analytical math packages, such as Wolfram Mathematica, breathed new life into the semiclassical method [30]. In [31], it was applied to the nonlocal generalization of the Fokker-Plank equation with solutions localized on a zero-dimensional manifold (a point). In [32,33], the WKB-Maslov method allowed to construct asymptotic solutions concentrated on curves for the nonlocal Fisher-Kolmogorov-Petrovskii-Piskunov equation. Semiclassical spectral series were constructed in [34,35] for the Hartree type equation with the special type of the external potential V tr .
Our method focuses on formal mathematical structures, so it is applicable for constructing solutions of a wider class of equations. In particular, the approach can be generalized for nonlinear Schrödinger-like equations such as Lugiato-Lefever equation, which is an amplitude equation for an optical pulse in the pumped ring resonator [36]. The interest in this equation was encouraged by discovery of the frequency comb associated with the Kerr temporal solitons [37][38][39]. The results of our paper may be of value for the theoretical study of the Bose-Einstein condensation phenomena that occurs in a gas of atoms and of various quasiparticles [3,[40][41][42][43][44]. Note that the conception of condensation is even involved in the study of dark matter [45,46].
The paper is structured as follows. In Section 2, we describe the asymptotic particular solutions of the non-stationary nonlocal GPE concentrated on curves. In Section 3, we construct the asymptotic solutions of the stationary nonlocal GPE based on the Section 2. In Section 3, the semiclassical spectral series are obtained. In Section 4, our formalism is illustrated by an example for the specific 2D stationary nonlocal GPE. The semiclassical solutions of the spectral problem for this equation are obtained in the analytical form. In Section 5, the concluding remarks are given. In Appendix, the auxiliary equations involved are given.

Nonlocal Gross-Pitaevskii Equation
The semiclassical approach to the spectral problem (4) can be developed for the nonlocal stationary GPE of a more general form Here,ẑ = (ˆ p, x), x ∈ R n ,ˆ p = −ih∂ x , and V(ẑ) is a linear partial differential operator depending onẑ with a symbol V(z) = V( p, x) smoothly depending on p and x ( p ∈ R n ). Our approach to the problem (7) relies on the method proposed in [25] for solution of the Cauchy problem in the semiclassical approximation for the non-stationary nonlocal GPE: This method allows one to construct the asymptotical solutions to the problem (8) in the class of functions semiclassically concentrated on curves of the phase space of the dynamical system of moments of the desired solution.
The core idea of the approach proposed is as follows. The method of [25] provides a practical possibility to select from the set of asymptotical solutions Ψ( x, t) to the Cauchy problem (8), found for various initial functions ψ( x), those Ψ( x, t) that have the form In the process of finding the functions (9), we come to ψ( x) and µ which are the eigenfunctions and eigenvalues of the problem (7). Consider solutions of (8) concentrated in a vicinity of the one-dimensional manifold Λ 1 t (the curve) in the 2n-dimensional phase space M 2n = R n × R n , p ∈ R n , x ∈ R n : Here, s is the curve Λ 1 t parameter and: For the functions Ψ semiclassically concentrated on the curve Λ 1 t , the dynamics of the curve in the limitsh → 0 is determined by the (1, 1)-type Hamilton-Ehrenfest system (the term (1, 1) implies that the parameter s is a one-dimensional variable and the vector Z(s, t) corresponds to the first order moments of the function Ψ( x, t), see [25] for details): Here, ||Ψ|| is the L 2 -norm, which is conserved for exact solutions of the non-stationary nonlocal GPE (8). Hereinafter, for the sake of simplicity, we will use the unit normalization ||Ψ|| = 1. The class J τ h of functions semiclassically concentrated on curves Λ 1 is given by where P t h is a class of trajectory-concentrated functions [47]: In this definition, the function ϕ( ξ, t, s,h) belongs to a Schwarz space S in the variable ξ ∈ R n , smoothly depends on t and s, and regularly depends on √h ash → 0; the function S(t, s,h) together with functions Z(t, s) = ( P(t, s), X(t, s)) define the class and are to be determined; the family of hypersurfaces s = ωτ( x, t) are given by relation [23,24] Hereinafter, we assume that rank X s (t, s) = 1.
This condition ensures the nondegeneracy of the localization curve (10) and the nonsingularity of the BEC density on the curve.
The construction of solutions of the form to the non-stationary GPE (8) in the class of functions J τ h in the semiclassical approximation with the accuracy of O(h 3/2 ) can be reduced to the solution of the associated linear GPE: (18) with the additional algebraic constraint Here, we denoted , are the solutions of the second order Hamilton-Ehrenfest system with integration constants determined from the initial condition for the function χ( x, t, s). In general case, these functions are given by: Here, Â is the mean value of the operator with a Weyl symbol A(z, t) in the class J τ h : where ξ is (n − 1)-dimensional vector of variables that complement the variable s = ωτ( x, t) to form the coordinate system in R n , D is the domain of the variables ξ, and the operatorẑ acts subject to the condition s = const . The particular solution of (18) with auxiliary algebraic constraint (19) can be written in terms of solutions of the variational system, which is the linear Hamiltonian system: where a(t, s) ∈ C 2n , J = 0 −I n×n I n×n 0 is unit symplectic matrix, (I n×n is the unit n × nmatrix). The function a 0 (t, s) =Ż(t, s) is the particular solutions of (23). Let (n − 1) linear independent solutions of (23) be known and satisfy the conditions Here, δ lk is the Kronecker symbol, {·, ·} is the skew-orthogonal product. Denote a k (t, s) = W k (t, s), Z k (t, s) . Then, the particular solutions of (18) reads [47] where N is the normalization coefficient, and The conditions (16) and (24) ensures the non-degeneracy of matrix C(t, s). In order to make sure that function (40) satisfies the auxiliary constraint (19) it is enough to check the validity of the identity which results from the property following from the definition of the matrix C(t, s). Let us introduce first order operators of the form Due to (24), the following commutation relations hold: These operators are the symmetry operators of (18) that preserve the property (19) of the solutions. Thus, the operatorsâ + k , k = 1, n − 1 are creation operators of (18), and the operatorŝ a k , k = 0, n − 1 are annihilation operators. The function |0, t (40) is the "vacuum" solution.
Let us construct the "occupation number representation" for (18) by determination of the following countable set of solutions: Here, In the next Section, we will construct stationary solutions of the GPE based on the solutions (31) and (25).

Stationary Solutions of the GPE
In order to obtain solutions that do not depend on time, we construct time-invariant curves. Let us consider curves determined by the following form of functions Z(t, s): Here, a parameter ω is the ratio between the periods by s and by t. The curve Λ 1 t determined by such the functions is mapped into itself with change of t, i.e., Λ 1 t = Λ 1 does not depend on t. By the change of variables (t + s/ω) → t and (t + r/ω) → r, the system (12) can be written in the following form: In this case, the equation of hypersurfaces s = ωτ( x, t) has the form s = ωτ( x) + t or Let the solutions of the second order Hamilton-Ehrenfest satisfy the relations g(t, s) = g(t + s/ω), g(t + T) = g(t), i.e.
The second order Hamilton-Ehrenfest system and the way of obtaining the functions satisfying the relations (36) are considered in the Appendix A.
For the solutions of the (1, 1)-type Hamilton-Ehrenfest system satisfying (33), the variational system (23) where we put a(t, s) = a(t + s/ω), and made the change of the variable (t + s/ω) → t. Since H zz (t) is the periodic function, we can pose the Floquet problem for the variational system (37). We know at least one solution of the Floquet problem for (37), namely a 0 (t) =Ż(t) (the T-periodic solution). Let there be (n − 1) more linear independent solutions of the Floquet problem for (37) that satisfy the following Floquet conditions: and the skew-orthogonality condition (24). Since only linearly stable Λ 1 are of interest for us, we consider the solutions of (37) with the real Floquet exponent Ω k . If conditions (33), (36) and (38) are met, the particular solutions (25), (31) of (18) read Then, considering (9) and (35), the relation (17) is as follows: The transition (42) from the extended (n + 1)-dimensional space to the n-dimensional configuration space can give the multi-valued function ψ( x). However, the function ψ( x) is the single-valued one if the following periodicity condition holds: In the next Section, we will study what constraints are caused by this condition.

Semiclassical Spectral Series
In view of the theorems given in Appendix B and relation (39), the condition (43) for solutions (40) and (41) reads [24] T 0 and assume µ(h) = µ (0) +hµ (1) where µ (1) (h) regularly depends onh. We consider the periodic trajectories of the (1, 1)-type Hamilton-Ehrenfest system with the period T = T(µ (0) ) determined by the following integral of motion of (34): Then, the quantization condition (44) is satisfied if The first equation in (48) implicitly determines µ (0) . The energy per particle for constructed semiclassical spectral series is given by Note that the quantization condition (48) describes the spectrum well in the linear theory when l ∼ 1/h.

Example
In this section, we will illustrate the formalism of our approach by a simple but nontrivial example that admits the analytical integration of the Hamilton-Ehrenfest system. We consider the following two-dimensional nonlocal Gross-Pitaevskii equation: Here, the units of length and energy are given by such the scales that the term V(ẑ) corresponds to −¯h 2 2m ∆ + 1 2 mω 2 harm x 2 in (1) in those notations. This equation can be treated as the nonlocal generalization of the model that describes the BEC in a radially symmetric trap. When κ > 0, the nonlinear term in (50) corresponds to the repulsive interaction. For the unit normalization of the wave function ψ( x), the coefficient κ is proportional to the number of condensate particles. Since the equation has the axial symmetry, the (1, 1)-type Hamilton-Ehrenfest system (34) admits solutions corresponding to the circumference in the configuration space. We look for a solution to the Equation (12) of the form X(t + s/ω) = R cos(ωt + s), sin(ωt + s) , where the parameter s is the angle in polar coordinates. The (1, 1)-type Hamilton-Ehrenfest system (34) for the GPE (50) reads From the first equation of the system (52) and (51), we have For the function (51), the integrals in the Equations (52) are as follows: where I 0 , I 1 are modified Bessel functions of the first kind, Γ = R 2 γ 2 . Here, we have taken into consideration that ω = 2π T and have used the property I n (z) = 1 2π π −π e z cos ϕ cos nϕdϕ, n = 0, ∞.
As a result, the (1, 1)-type Hamilton-Ehrenfest system yields the following equation for R: Consider the approximation R γ. It meets the strongly local interaction. The following asymptotic estimate holds for the modified Bessel functions: Then, we have This equation has the unique solution. Thus, the right-hand side of (53) satisfies the solvability condition with respect to R for at least sufficiently small γ (the strongly local interaction). Note that the Equation (53) degenerates into ω = 1 for κ = 0, i.e., the period T depends on the chemical potential µ (0) only in the presence of the nonlinearity in (50). The first quantization condition in (48) for the solutions (53) reads The energy per particle the zeroth order approximation with respect toh is as follows: The Equation (15) for the function t + s/ω = τ( x) for the curve (51) reads −x sin(ωt + s) + y cos(ωt + s) = 0, and it has the solution The semiclassical solutions of (50) for ν = 0 and ν = 1 are obtained in Appendices C and D, respectively.
The solution for ν = 0 (the vacuum solution) reads where θ is the positive constant defined in Appendix C. The solution for ν = 1 is given by In Figures 2 and 3, the plots of |ψ ν ( x)| 2 and Arg [ψ ν ( x)] are shown for ν = 0 (56) and ν = 1 (57) when γ = 1, κ = 1,h = 0.1, and l = 10 (it corresponds to ω ≈ 0.883, R ≈ 1.06, θ ≈ 0.528). The phase incursion of 20π along the localization curve (Figures 2 and 3) indicates that we have the vortex solution that corresponds to the solutions of the GPE with 10 quantum vortices [48]. However, since we construct asymptotics only in the neighborhood of the circumference, it do not allows us to obtain the complete information about the phase features in the center.

Conclusions
We applied the method proposed in [25] to constructing the solutions of the spectral problem for the stationary nonlocal GPE (7). The considered equation has rather general form with the limitation of the smoothness of the kernel of the nonlocal nonlinearity and of the symbol of the linear part. The solutions of the spectral problem were obtained by choosing the special solutions of the non-stationary nonlocal GPE semiclassically concentrated on time-invariant curves. Such solutions are found in terms of periodic solutions of the second order Hamilton-Ehrenfest system (A2) and (A5) and the set of linear independent skew-orthogonal Floquet solutions of the variational system (37).
The generalization of the method [25] for the stationary GPE is based on the following ideas. The time independence of the Equation (8) operator yields the symmetry (33) with respect to the tand s-shift for the Equation (12) if the system (34) admits periodic solutions. The property (33) provides to the same symmetry for the second order Hamilton-Ehrenfest system, which leads to the Equations (A2) and (A5), and to the periodicity of coefficients in the variational system (23), that allows us to pose the Floquet problem for (37). Then, the periodic solutions of (A2) yield the symmetry for the associated linear GPE (18) and (20) that effectively reduces the number of independent variables in (18) from (n + 2) (n variables x, t, and s) to (n + 1) (n variables x, and (t + s/ω)). This symmetry together with the Floquet solutions of (37) allows us to obtain the particular solutions (40) and (41) of the associated linear GPE (18) and (20) that give the solutions (42) of the stationary GPE (7) subject to the quantization conditions (48).
The feature of our method is that we do not reduce the initial n-dimensional problem to the one-dimensional. It means that, in general, our approach does not require the GPE to have some special symmetries. Nevertheless, various concepts of symmetries may be required for analytical integration of the auxiliary Equations (34), (37) and (A5) (e.g., [49]).
The solution constructed for the simplest form of the external field reminds the giant hole solution that was studied in [50]. It is known that vortices in the BEC are quantized so that each vortex corresponds to the phase incursion of 2π (see, e.g., [48]). The giant hole solutions are characterized by the proximity of vortex centres so that they form one density pit of the BEC [50] like our example. The indistinguishability of the vortex centres is the consequence of the semiclassical approximation. The centres of vortices likely can be identified within our approach by the consideration of a less trivial localization curve that is an open problem.
The solutions obtained here within the semiclassical approach proposed have the limit κ = 0 (the linear case) and decrease rapidly (exponentially) with distance from the localization curve. In the sense of these properties, the solutions of the nonlocal GPE correspond to ones studied in [11] for the one-dimensional local GPE. However, although the solutions localized on a curve cannot be directly reduced to the one-dimensional solutions found in [11], they explicitly reflect topological features of the BEC.
Note that the numerical calculations for the local and nonlocal GPEs show that, for the delta-like interaction kernel, the solutions of these type weakly change with an increase in the nonlocality coefficient γ, and qualitatively do not differ from the solutions of the local GPE when γ is small. The solution is sensitive to the variation of γ when the localization domain of the interaction potential becomes larger than the localization area of the BEC density. It caused by the sharp decrease in the interaction energy. Therefore, the nonlocal GPE can be treated not only as a separate problem but as the deformation of the local GPE for using our approach.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. The Second Order Hamilton-Ehrenfest System
We can reduce the nonlocal GPE (8) to the parametric family of the associated linear GPEs (18) and (20) due to the possibility of obtaining the dispersion matrix σ xx and the first order correction with respect toh for the first order moment X (1) from the auxiliary system of equations in the semiclassical approximation, which is called the second order Hamilton-Ehrenfest system. This system was derived in the general case in [25]. In this Appendix, we consider the second order Hamilton-Ehrenfest system for the special type of solutions (9) to (8).
The matrix of second order moments ∆ 2 (t) accurate to O(h 3/2 ) is determined by the leading term of asymptotics since it corresponds to an operator with the estimate O(h). The formula for ∆ 2 (t) reads where l, k = 1, 2n, ξ is (n − 1)-dimensional vector of variables that complement the variable t = τ( x) to form the coordinate system in R n , and D is the domain of the variables ξ. Here, the operator ∆ẑ acts subject to the condition t = const . Note that the function H(t,h)[g ν (t)] included in the relation (40) affects only the phase of χ ν ( x, t) and commutes with ∆ẑ. Hence, the matrix σ xx (t) does not affect the result of the formula (A1), i.e., this formula is not recursive. The function given by (A1) satisfy the equatioṅ In general case, the solution of (A2) reads [25] where A(t) is the solution of the matrix variational system and the matrix D 0 = ∆ 2 (t) t=0 is given by (A1). In view of (33) and (36), the system of equations for Z (1) (t) = P (1) (t), X (1) (t) can be written as dr h W y X(t), X(r) , X (1) (r) + where the operator ∂ z = ∂ ∂z = ∂ ∂ p ; ∂ ∂ x acts on the arguments of the function V(z) and W( x, ·), and the function π 0 (t) is given by It was shown in [25] thatπ 0 (t,h) = O(h 3/2 ). Hence, in the system (A5), we can put π 0 (t) = const = π 0 (t 0 ) = π 0 .
In view of (A7), the system (A5) is the system of linear integro-differential equations. Note that, for every ν, the function ∆ 2 (t)[χ ν ] and π 0 [χ ν ], which are included in (A5), are different. As opposed to ∆ 2 (t), the function Z (1) (t) is not determined by the leading term of asymptotics. Thus, our formalism relies on the assumption that system (A5) for every ν admits T-periodic solution under some initial conditions for Z (1) (t).
Proof. In view of (26), (24) can be written as This relation is also valid for B(t) Then, Theorem A2. For the matrix C(t) given by (26), (24) and (37), the following relation holds: Proof. The statement of this theorem directly follows from (A9).
Appendix C. The Solution of (50) for ν = 0 The variational system for the function a(t) = W(t), Z(t) in view of (53) takes the form   ˙ Z(t) = W(t), For (A13), we seek the Floquet solutions satisfying (24). The system (A13) can be written as where g = κ − ω 2 , κ = −κ Γ 2γ 2 e −Γ I 0 (Γ) − 2I 1 (Γ) + I 2 (Γ) are scalar coefficients. Let us make the following change of variables corresponding to the rotation of the vector Z(t) by the angle ωt: Then, the variational system reads We seek U(t) in the form In view of g = κ − ω 2 , this relation can be written as The characteristic equation has the roots Substitution of g = κ − ω 2 into this relation yields It is clear that we have only aperiodic solutions for κ ≥ 2ω 2 , which do not satisfy the Floquet condition (24). Hence, we consider the case when κ < 2ω 2 . Denote θ = ω 2 − κ 2 .
The coefficient N is obtained from the unit normalization condition: Note that V zzz (z, t) = 0 in our example. Hence, it is not necessary to obtain the whole matrix ∆ 2 (t) for the system (A5) but it is enough to obtain only σ xx (t), which is given by Substituting (40) into (A6), the formula for π 0 is as follows: In view of (27), we have Substitution of (45) into this formula yields The identity (A14) is the integral of motion of the system (A5) for ν = 0. The system (A5) for the coordinate part of Z (1) (t) yieldṡ For the pulse part Z (1) (t), we havė The integral of the system (A14) can be written as T 0 dr exp[Γ cos(ω(r − t))] X(r) − X(t) , X (1) (r) + +ω 2 X(t), X (1) We seek the solution of (A15) in the form Then, for (A16), we get The constant vector U (1) (t) = U (1) satisfies (A17). Let us show that it also satisfies the Equation (A15). The integral (A17) with U (1) (t) = U (1) is as follows:  . Finally, we obtain X (1) (t) = U (1) R · X(t), P (1) (t) = U (1) R · P(t).