The Gross–Pitaevskii Equation with a Nonlocal Interaction in a Semiclassical Approximation on a Curve

: We propose an approach to constructing semiclassical solutions for the generalized multidimensional Gross–Pitaevskii equation with a nonlocal interaction term. The key property of the solutions is that they are concentrated on a one-dimensional manifold (curve) that evolves over time. The approach reduces the Cauchy problem for the nonlocal Gross–Pitaevskii equation to a similar problem for the associated linear equation. The geometric properties of the resulting solutions are related to Maslov’s complex germ, and the symmetry operators of the associated linear equation lead to the approximation of the symmetry operators for the nonlocal Gross–Pitaevskii equation.


Introduction
The development of new mathematical methods for studying Bose-Einstein condensates (BECs) with complex geometric and topological properties is motivated by the progress in experimental physics. New types of traps with toroidal, ring-shaped, and other topologies allow one to produce condensates with the density localized in the vicinity of a closed or unclosed curve like a circumference, a cigar shape, an elongated line, etc (see, e.g., the review paper [1], the recent paper [2], and the references therein). The exploration of BECs with various topological features contributes to a deeper understanding of the BEC properties and the mechanisms of its creation, and opens new prospects for developing precise measurement technologies using coherent properties of matter waves [2].
Theoretical approaches to modeling diluted coherent quantum ensembles of Bose particles at low temperatures in the BEC phase use mean-field approximations. Within the framework of a mean-field theory, many of the effects characteristic of a BEC and its evolution can be described on the basis of the Gross-Pitaevskii equation (GPE) [11,12] (see also [4]) and its modifications that take into where ∂ t = ∂/∂t, ∇ is the gradient operator andh is the Planck constant. The mean-field theory involves a single effective particle characterized by a function Ψ( x, t), which is the order parameter (classical field or wave function) of the condensate depending on the space and time variables ( x and t, respectively). The external potential U( x, t) (trap potential) in Equation (1) can be almost arbitrarily designed, and κ|Ψ( x, t)| 2 is the inter-atomic short-range and isotropic contact (local) interaction potential. The parameter κ is proportional to the s-wave scattering length [11,12]. The BEC phenomena emerging in various trapping geometries were studied using analytic and numerical methods for the 3D GPE with a local nonlinearity term (1) by a number of researchers.
Brand and Reinhardt [13] numerically simulated the three-dimensional BEC in a toroidal trap geometry using a stochastic GPE and showed that the BEC formation can spontaneously generate quantized circulation of the newborn condensate. Zurek studied the formation of grey soliton states in a BEC in elongated traps [14].
As the BEC properties substantially depend on the interactions, it is of interest to elucidate which changes in the BEC properties can be related to the nonlocal inter-atomic interaction. The robust way of modeling a BEC with non-local interactions can be based on the mean-field theory of many-particle quantum mechanics and the Hartree-Fock approximation (for details, see, e.g., [15,16], and references therein).
Arguments in favor of non-local generalizations of the Gross-Pitaevsky equation are discussed in [17], where, in particular, the orbital stability of the solutions in a potential of a special type has been revealed.
A detailed consideration of the interparticle interactions and long-range correlations arising in the ensemble's interacting atoms leads to a nonlocal generalization of the GPE (1), which can be written as where the nonlocal term, W[Ψ(t)]( x) = Ψ( y, t)|W( x, y)|Ψ( y, t) = R n W( x, y) |Ψ( y, t)| 2 d y, can be treated as an expectation of a real weight function W( x, y) depending on spatial points x and y. The function W( x, y) is maximal at x = y and decreases with increasing distance between x and y.
The size of the domain where the function W( x, y) substantially differs from zero characterizes the degree of nonlocality of the interaction. The concept of the nonlocal interaction arises naturally in theoretical mean-field studies of dipolar quantum gases, since the dipolar interaction between dipoles described by an appropriate function W( x, y) is anisotropic and has a long-range character (see reviews [18,19]). The creation of the dipolar BEC was reported in [20,21].
For analysis of models built on the nonlocal GPE, exact and approximate analytical methods are required. Some important results for non-local versions of the GPE were obtained in the last decade. For instance, the conditions of the Cauchy problem for the nonlocal GPE (2) were analyzed for U = 0 and different forms of W [22]. Analytical approaches for the nonlocal GPE were studied in [23], particular solutions were found in [24], and the collapse problem of the localized waves described by (2) was discussed in [25].
However, analytical methods for constructing solutions of the nonlocal GPE in a multidimensional space-time when the solution has rather complex geometric properties were not developed due to the mathematical complexity of the problem. This motivated us to seek an approach for constructing solutions of the Cauchy problem for the nonlocal GPE describing a BEC with complex topological features to contribute to non-linear methods of mathematical physics.
In the search for a possible approach to solving this problem, we note that the expectations of operators, calculated by analogy with Equation (3) from the solutions Ψ( y, t) of Equation (2), should satisfy special conditions resulting from the properties of Equation (2). Consider the problem of obtaining a closed extended system of equations describing the expectations of W( x, y) and including the original equation. As a result, the original Cauchy problem will be reduced to solving the resulting extended system. Should this procedure be implemented, it may be possible to find the exact solution of the original nonlinear Equation (2) (see [26][27][28]).
For the general case, asymptotic techniques can yield an approximate solution of the equation with a prescribed accuracy. Note that the deduction of this "extended" system and its solution should be mutually consistent within the accepted approximations. In view of this, the use of the WKB-Maslov method is justified, as it captures all of the core items and is well-proven in solving a wide class of quantum-mechanical problems.
Our aim is to apply the well-known semiclassical WKB-Maslov method to construct solutions concentrated in the neighborhood of a certain manifold of a nonzero dimension. To do this, we extend the dimension of the original nonlocal GPE in order that the basic theorems of the Maslov theory, which hold for solutions concentrated around a point (zero-dimensional manifold), can be used. The finding of solutions of the nonlocal GPE in extended dimensions has some specific features, the analysis of which constitutes the essence of the approach proposed.
Semiclassical asymptotics concentrated on linearly stable closed geodesics (Λ 1 ) were first obtained in [29] for linear equations. These asymptotics are local wave packets which are gaussian along the direction transverse to Λ 1 and oscillatory along Λ 1 . Afterwards, this case was studied in [30][31][32].
The rigorous theory of the semiclassical quantization of nonintegrable Hamilton systems for linearh −1 -(pseudo)differential operators is based on the Maslov complex germ method [33][34][35]. In this theory, the problem of the construction of semiclassical asymptotics is reduced to the construction of geometric objects in a 2n-dimensional phase space (family of Lagrangian manifolds Λ k with a complex germ r n ). Each manifold Λ k is generated by a solution of the classical Hamilton system and has a dimension k, 0 ≤ k < n. The complex germ r n is generated by a special set of n linearly independent complex solutions of the variational system (Hamilton system linearized in the neighborhood of Λ k ).
Generalization of the above constructions for nonlinear quantum systems is nontrivial because it is not clear which classical equations correspond to the semiclassical limit (h → 0) of a given nonlinear quantum system. The form of these equations depends on the class of functions in which the asymptotic expansion is constructed and on the properties of nonlinear operators in this class.
Asymptotics of the Hartree-type equation with a Coulomb self-action were studied in [36][37][38], where the authors obtained solutions concentrated on low-dimensional manifolds using a "singular" version of the WKB method. For a Hartree-type operator with smooth symbols in its linear part and a self-action term, semiclassical solutions concentrated on n-dimensional manifolds Λ n were obtained in [39][40][41], where other classical equations were used.
In this paper, we partially generalize these results to a one-dimensional manifold Λ 1 and propose a new method for constructing the Cauchy problem's approximate solutions concentrated on a curve for a multidimensional nonlocal GPE.
The paper has the following structure. In Section 2, we introduce a concentration manifold in a phase space and describe its evolution in a semiclassical approximation. In addition, a class of functions semiclassically concentrated on the manifold is defined and its main properties are described. In Section 3, we show that in the given class, the semiclassical approximation allows us to reduce the nonlocal GPE to a partial differential equation with an additional algebraic condition. In Section 4, we deduce auxiliary differential equations, which must be solved to linearize the nonlocal GPE. In Section 5, the linearized GPE is considered. We obtain the evolution operator of the linearized GPE and clarify its connection with the asymptotic solutions of the nonlocal GPE. In Section 6, the symmetry properties of the associated linear GPE are discussed. In Section 7, concluding remarks are given.

The Class of Functions Concentrated on a Curve
The nonlocal Gross-Pitaevskii Equation (2) can be written in a generalized form as Here,ẑ = (ˆ p, x),ˆ p = −ih∂ x , x ∈ R n , and V(ẑ, t) is a Weyl-ordered linear operator with a symbol V(z, t) [46,47]. The linear operator V(ẑ, t) in Equation (4) stands for 1 2m (−ih∇) 2 + U( x, t) in Equation (2). A more general form of V(ẑ, t) extends the applicability of Equation (4). In particular, Equation (4) can be considered in a non-inertial frame of reference. The function W( x, y, t) can describe both short-range and long-range interactions including the dipolar-dipolar interaction if the integral in Equation (4) is regularized in an appropriate way; for instance, the explicit cutoff and the single-mode approximation were used in [48].
Note that for Ψ in Equation (4), the L 2 -norm squared, if it exists, is conserved during the time evolution. In addition, in what follows, by GPE, we often imply the nonlocal GPE if this does not cause confusion. For constructing concentrated solutions to Equation (4), we introduce a k-dimensional manifold Λ k t in a 2n-dimensional phase space M 2n = R n × R n , p ∈ R n , x ∈ R n : Here, s = (s 1 , ..., s k ) is the Λ k t manifold parameter, k < n, and Z(s, t) are the functions defining the manifold.
Let us explain what is meant by the localization of functions Ψ( x, t,h) on the manifold Λ k t of the phase space M 2n .
for an arbitrary operatorÂ with a Weyl symbol A(z, t). Here, σ(s) is a weight function.
Here, we use the notations the norm squared, ||Ψ|| 2 , is given by (5) and it is assumed that the integrals in these expressions exist.
We also assume that the Weyl symbol A(z, t) of the operatorÂ is an infinitely differentiable function with respect to all its arguments and it grows no faster than a polynomial with |z| → ∞, z ∈ R 2n for every instant t.
Localized solutions of the GPE (4) are possible on manifolds (6) satisfying some special dynamic conditions.
The time evolution of the manifold Λ k t in (6) is given by a 2n-dimensional vector Z(s, t). To derive the equations governing the dynamics of Z(s, t), we consider the expectation determined by (8) for an operatorÂ(t) = A(ẑ, t) Hermitian with respect to the L 2 -scalar product (5).
Hereinafter, for the simplicity of notation, we will omit the variableh in formulas where this does not cause confusion.
Differentiating (8) with respect to t and taking ∂Ψ( x, t)/∂t from (4), we get where Â ,B =ÂB −BÂ is the commutator of the linear operatorsÂ andB. We call (9) the Ehrenfest equation for the nonlinear equation (4) by analogy with quantum mechanics.
In the limith → 0, Equations (7) and (9) yield Here, is a Poisson bracket [46], Z(s, t) is defined in (6), z, w = ∑ 2n j=1 z k w k is the inner product in R 2n , and J is a symplectic identity 2n × 2n matrix where I n×n and 0 are an identity and a zero n × n matrix, respectively.
The term "(k; 1)-type" implies that the parameter s in system (12) is a k-dimensional vector, and the vector Z(s, t) corresponds to the first moments of the functions Ψ( x, t) in (9), where Ψ( x, t) are the solutions of Equation (4).
Definition 1 has nothing to do with the GPE (4) itself. However, if the function Z(s, t) satisfies the (k; 1)-type Hamilton-Ehrenfest system (12) with the initial condition Z 0 (s) given by where the weight function σ(s) is determined by the manifold Λ k t geometry, then the function Ψ( x, t,h) in (7) is the asymptotic solution of (4).
In what follows, we will consider the one-dimensional manifolds Λ 1 t for k = 1, i.e., the curves in the phase space M 2n , for which (6) reads The curves Λ 1 t can be considered as a family of trajectories of the Hamilton-Ehrenfest (HE) system (12) parametrized by s ∈ [s 1 , The definition of the class of functions by relation (7) is not constructive. Therefore, we consider a class of functions with the property (7) that can be easily constructed. A sufficient condition for a function ϕ to be concentrated on the curve Λ 1 t is that it belongs to the class (14): ϕ ∈ J τ h . It is important that the asymptotic estimates that hold on this class allow us to construct approximate solutions of Equation (4).
We now describe the class J τ h of functions localized on curves (13). The class J τ h is defined as where P t h is a class of trajectory-concentrated functions [49], Here, the function ϕ( ξ, s, t,h) belongs to a Schwarz space S in the variable ξ ∈ R n , smoothly depends on t and s, and regularly depends on , the 2n-dimensional vector-function Z(s, t) = ( P(s, t), X(s, t)) is the solution of the Hamilton-Ehrenfest system (12), and, by definition, the function S(s, t,h) = S (0) (s, t) +hS (1) Note that the relation (16) indicates that the function S (0) (s, t) for the Hamilton-Ehrenfest system (12) is an analog of an action in classical mechanics.
Define the family of hypersurfaces s = τ( x, t) by the condition [34,35] Note that the implicit function theorem holds for τ( x, t) when, for every instant t, the relation is satisfied in the neighborhood of each solution ( x, s) of Equation (17). The introduction of the class of functions J τ h allows us to estimate the basic operators −ih∂ t , z = (ˆ p, x) entering into Equation (4) and, hence, represent (4) aŝ where Ψ ∈ J τ h , and the operatorL 1 can be considered small, i.e.,L 1 ∼Ô(h α ), α > 1. Hereinafter, where O(h γ ) is the uniform-in-time estimate for t ∈ [0, T], T = const, and the norm ||Ψ|| is defined in (5). Thus, constructing approximate solutions to Equation (4) is reduced to constructing exact solutions to the "unperturbed equation determining the principal term of the semiclassical asymptotics and the equations determining the higher-order corrections. The next section is devoted to the obtaining of the operatorL 0 ; its explicit form is given there. We call Equation (21) the reduced Gross-Pitaevskii Equation (4) in the semiclassical approximation. The properties of the class P t h , which were studied in [49], hold if s is an independent variable (s ∈ [s 1 , s 2 ]). So, we need to increase the dimension of the GPE to use them. In the class J τ h , a momentum operatorˆ p = −ih∇ can be represented as Here, we used relation (14), valid for all Φ( x, t) ∈ J τ h . Let us introduce the designationŝ To construct asymptotic solutions to the Gross-Pitaevskii Equation (4), we use the following asymptotic estimates. Theorem 1. The following asymptotic estimates hold in the class of functions J τ h : Here, {∆ẑ} α is the operator with the Weyl symbol (∆z) α , As mentioned earlier, all estimates through the text are uniform in time for t ∈ [0, T].
Proof. We give a proof only for the case k = 1, s ∈ [0, T], as the solutions concentrated on curves are considered in this paper. However, this proof can be generalized for multidimensional manifolds.
1. The symbol of the operator {∆ẑ} α can be written as and the mean value σ α (h) of the Weyl-ordered pseudo-differential operator {∆ẑ} α is given by The explicit form of a function from the class J τ h is Then, we have Here, ∆ y = y − X(τ( y, t),h).

Let us define an n × n-matrix A(t) as
where According to (30), the inverse of the matrix A(s, t) is given by Define the variables ξ = (ξ 1 , ξ 2 , . . . , ξ n ) ∈ R n by the relation In view of the definition of a hypersurface s = τ( x, t) (17), the variable ξ 1 is identically zero (ξ 1 ≡ 0).
Then, for the function χ( x, t, s), we have and hence The change of variables in the integral of (28) yields where ζ, ω is given by (17), (32), and (34). Then, we have . . dω n and we used the following relation: In particular, we have The function ϕ( ξ, s, t,h) regularly depends on √h , so M α (t,h) and M 0 (t,h) also regularly depend on Proof. 1. The relation (36) is just the particular case of (24).
2. According to the definition (15), the s-derivative of χ( x, s, t,h) reads As the function ϕ ∆ x √h , s, t,h smoothly depends on s, we havē The expectation Â (t) of the operator A(t) = A ẑ, τ( x, t), t is unsuitable for constructing semiclassical asymptotics on a manifold, as it does not provide any information about the behavior of solutions in the neighborhood of the manifold. Similarly to Z(s, t), we will define the expectation of the operator with the Weyl symbol A( x, s, t) in the class of functions J τ h . Let x −→ (s, ξ 2 , ..., ξ n ) be a change of variables in a coordinate space R n , where the variables ξ ∈ W ⊂ R n−1 complement the variable s to form the coordinate system in R n with the Jacobian J(s, Let us define the expectation A(t) of the operator A(t) with the Weyl symbol A( x, s, t) in the class of functions J τ h as It is easily seen that the following relations hold: and where f is an arbitrary function of the variables t and s.
In the next section, we focus on solving Equation (21), and the operatorL 1 in (19) is not considered. With the aid of estimates (36), we will obtain the reduced Equation (21) in explicit form and construct its solutions.

The Reduced Gross-Pitaevskii Equation
In this section, we seek a solution concentrated on a curve Λ 1 t for the Gross-Pitaevskii equation in the sense of (7). To this end, we involve the results of [26,28,[42][43][44], where asymptotic solutions to the nonlocal GPE were constructed in the class P t h of the trajectory-concentrated functions (15) accurate to O(h 3/2 ) using an associated linear equation of the Schrödinger type.
In the next section, we obtain the evolution equation for these moments.

Equations for the First and the Second Moments
The functions Z (1) (s, t) and α 0,2 (s, t,h) in Equations (55) and (57) depend on χ( x, s, t), and the Ehrenfest system for the moments of the function χ( x, s, t) includes an infinite number of equations. However, it is necessary to obtain only the first approximation of the functions Z (1) (s, t) and α 0,2 (s, t,h) in order for estimate (55) to be correct. As well as for Z(s, t), Equation (58) for the functions Z (1) (s, t) and α 0,2 (s, t,h) can be linearized in a semiclassical approximation. To obtain the required equations for the moments, we need some additional results.
Let A(t) be a Hermitian operator with respect to the scalar product (5) for the functions of the class ,t) , the equation for the expectation follows from (9) and (40): In view of estimates (36), the integral part of Equation (4) can be approximately expressed in terms of the moments of the function χ( x, s, t), given by (39).
The moments of the trajectory-concentrated functions can be obtained from an auxillary system of differential equations which are semiclassical approximations to the Ehrenfest equation (see [26][27][28][42][43][44]). In the class J τ h , the expectations are given by (58).
Proof. In (58), let us expand the functions V(ẑ, t) and W( x, y, t) in power series in ∆ẑ. Consider the series up to the third-order terms to obtain an equation accurate to O(h 3/2 ). Using the commutators where the operatorẑ is subject to the condition s = const , we obtain the third powers of the ∆ẑ coefficients accurate to O( √h ), the second powers of the ∆ẑ coefficients accurate to O(h), and the first powers of the ∆ẑ coefficients accurate to O(h 3/2 ) in the same way as we did in (49). These coefficients are as follows: dr σ(r) W X(s, t), X(r, t), t + W y X(s, t), X(r, t), t , X(r, t,h) − X(r, t) + Thus, we arrive at (62).
The general solution of the second-order Hamilton-Ehrenfest system (59), (61), and (62) will be denoted by g(t, s, C), where C is the set of integration constants.

The Associated Linear Gross-Pitaevskii Equation
Making a formal replacement of g[χ](t, s) by g(t, s, C) in Equation (66) Below, we denote the solution of the Cauchy problem for Equation (68) with the initial condition (67) as φ = φ( x, t, s, C).
Consider the following equation: Choose the initial conditions C in accordance with the initial condition for the function χ( x, s, t, C), . Then, the function χ( x, s, t, C[ϕ]) satisfies (55) and (56).
Proof. The operatorâ 0 (s, t) = a 0 (s, t), J∆ẑ is the symmetry operator of Equation (69) if it commutes with the equation operator. It can be shown that this condition is satisfied if a 0 (s, t) is a solution of the variational system a 0 (s, t) = JH zz (s, t)a 0 (s, t).
Taking into consideration the above theorems, we can write Equations (69) and (70) as Equation (72) Here, ∆y = y − X(s, 0), and the matrix M(s, t) = The solution χ( x, s, t, C) of the ALGPE can be obtained as whereÛ(s, t, C[ϕ 1 ]) is the exact evolution operator of the ALGPE. As the operatorâ 0 (s, t) commutes with the evolution operators of Equation (72) in the sensê a 0 (s, t)Û(s, t, C) =Û(s, t, C)â 0 (s, 0), the function χ( x, s, t, C) satisfies both equations of the systems (72) and (73). Let ϕ 1 ( x, s) = ϕ( x, s). Then, χ( x, s, t, C[ϕ]) s=τ( x,t) is the asymptotic solution of Equation (4) that corresponds to the initial condition ϕ( x, τ( x, t)) t=0 . So, the evolution operatorÛ(s, t, C) of the ALGPE is the key tool for constructing the solution to the Cauchy problem for the nonlocal Gross-Pitaveskii equation. In the class J τ h , the asymptotic solution can be obtained for any initial condition which satisfies the constraint (73). The construction of the evolution operatorÛ(s, t, C) consists of solving the systems (12), (59), (62), (61), and (75).
The initial condition (67) generates corresponding initial conditions for the second-order Hamilton-Ehrenfest system: Denote the solutions of the second order Hamilton-Ehrenfest system (77) as g[χ 0 (s)](t, s). From the identity In view of relation (42), the function χ( x, t, s) (79) yields the solution Ψ( x, t, s) of the GPE (7). This solution is obtained from the function χ( x, t, s) (79) using relation (42).

Symmetry Operators of the Associated Linear Gross-Pitaevskii Equation
In this section, we discuss the symmetry operators of the ALGPE as an alternative way of generating asymptotic solutions. The construction of the symmetry operators of the ALGPE is based on the solutions of the variational systeṁ a(s, t) = JH zz (s, t)a(s, t), where a(s, t) ∈ C 2n . The first-order symmetry operator corresponding to the solutions of the system (80) is given byâ (s, t) = N a(s, t), J∆ẑ , where N j = const. Obviously, the Hermitian conjugate operatorâ + (s, t) is the symmetry operator as well. We are interested only in symmetry operators which commute with the operatorâ 0 (s, t) = Z s (s, t), J∆ẑ , as they convert the ALGPE solution that satisfies the additional condition (56) to a similar solution. Consider (n−1) linearly independent complex solutions a j (s, t) of (80), which are skew-orthogonal to Z s (s, t): We have shown that the asymptotic solutions of the GPE that are concentrated on the manifold Λ 1 (13) in the sense of (7) can be constructed on manifolds (on the curves k = 1) whose dynamics are determined by the system of integro-differential equations (the first-order Hamilton-Ehrenfest system (12)). This integro-differential Hamilton-Ehrenfest system plays the role of the "classical" equations of motion corresponding to the quantum system described by the nonlocal GPE (2). The explicit dependence of the first-order Hamilton-Ehrenfest system on the dimension k of the manifold Λ k is the direct consequence of the nonlinearity of the original Equation (2). Thus, the evolution of the manifold corresponding (in the "classical" limit) to the quantum system described by Equation (2) depends on the manifold dimension k. Note that the first-order Hamilton-Ehrenfest system for the case k = 0 was obtained in [42,43].
The main idea of the method proposed is extending the dimension of the original GPE (2) by introducing the parameter s of the curve Λ 1 t as a new variable of the GPE according to relation (42). In the extended dimension, the solution is constructed using the technique of the trajectory-concentrated functions developed for the case k = 0 [42,43]. The passage from solving the Cauchy problem found for the GPE in a higher dimension to solving the Cauchy problem for the original GPE is carried out by truncating the domain of the solution to the family of hypersurfaces (17) orthogonal to the curve Λ 1 t using relation (42). Our method takes advantage of the nonlocality of the initial problem and can be applied to a wide range of kernels W( x, y, t) in (3) and (4), including the strongly localized one (e.g., Gaussian type), the long-range one (e.g., dipolar-dipolar), or their combination. A sufficient condition for the applicability of the method under consideration with respect to the specific kernel is the convergence of integrals over the curve Λ 1 t in (12), (57), and (62). There is the so-called giant hole solution of the GPE (a giant vortex absorbing phase singularities, see [50] for details) describing the fast-rotating condensate. Such a solution holds for both the local and nonlocal GPE and has rather semiclassical behavior, so it is a good example of semiclassical states localized on a curve arising in physical problems.
Note that we consider formal asymptotic solutions of the nonlocal GPE with the arbitrarily small right-hand-side residual with respect to the parameterh,h → 0. The substantiation of these asymptotics at a given finite interval t ∈ [0, T], T = const is the separate nontrivial mathematical problem. This problem is related to the obtaining of priori estimates of the nonlinear Equation (4) that are uniform with respect to the parameterh, and is not considered in this paper. For heuristic reasons given in [41], it seems that the estimate of the difference between exact and asymptotic solutions can be obtained using the methods that are developed in [33,34].
The equations of the GP model have a similar structure to those of reaction-diffusion systems [51][52][53], which have a number of specific features such as, for example, the formation of patterns, special kinds of waves, etc. Therefore, the methods applied to this model can be of more general interest from a mathematical point of view. In addition, the considered model catches characteristic features of the collective behavior of the condensate at very low atomic concentrations. Such features may be of interest in studies of properties of systems with low concentrations (e.g., [54]).
In conclusion, we note that a similar problem was studied for the nonlocal Fisher-Kolmogorov-Petrovskii-Piskunov equation [55,56].
However, although the method proposed here and the method developed in [55,56] are ideologically close, they have significant differences and require independent implementations.