Family of asymptotic solutions to the two-dimensional kinetic equation with a nonlocal cubic nonlinearity

We apply the original semiclassical approach to the kinetic ionization equation with the nonlocal cubic nonlinearity in order to construct the family of its asymptotic solutions. The approach proposed relies on an auxiliary dynamical system of moments of the desired solution to the kinetic equation and the associated linear partial differential equation. The family of asymptotic solutions to the kinetic equation is constructed using the symmetry operators acting on functions concentrated in a neighborhood of a point determined by the dynamical system. Based on these solutions, we introduce the nonlinear superposition principle for the nonlinear kinetic equation. Our formalism based on the Maslov germ method is applied to the Cauchy problem for the specific two-dimensional kinetic equation. The evolution of the ion distribution in the kinetically enhanced metal vapor active medium is obtained as the nonlinear superposition using the numerical-analytical calculations.


I. INTRODUCTION
Kinetic equations are the theoretical footing for the dynamic phenomena of various nature that occur in physical systems of many interacting elements (particles). Examples of such system are diluted gases, gas discharges, a plasma, processes of coagulation [1], and biological systems such as, e.g., population systems [2][3][4]. In some systems, nonlocal collective (averaged) interactions of elements substantially contribute to dynamics. Such interactions are modeled by integral terms in kinetic equations that become integro-differential. In spatially heterogeneous kinetic phenomena, the interelement interactions occur along with the diffusion. Then, the model kinetic equation belongs to the class of reaction-diffusion (RD) equations. The study of RD equations with both local and nonlocal terms have formed an independent branch of mathematical physics. * shpv@phys.tsu.ru † aek8@tpu.ru ‡ ssaykmh@yandex.ru Due to the mathematical complexity of the study of RD equation with nonlocal interactions, methods of computer modeling prevail here. However, the demand for the analytical methods stimulates the development of approximate and asymptotically exact solutions. For a number of RD kinetic equations with nonlocal interactions, one can succeed using the WKB-Maslov theory of semiclassical approximation or the Maslov complex germ method [5][6][7]. Based on the WKB-Maslov theory, the method of semiclassical asymptotics was developed for a generalized Fisher-Kolmogorov-Petrovskii-Piskunov equation (Fisher-KPP) with a quadratic nonlocal term in [8,9] and for the nonlocal Gross-Pitaevskii equation in [10,11].
In this work, using the results of [9,12], we construct semiclassical asymptotics for the model kinetic equation with the nonlocal cubic nonlinearity of the form y, z, t)u( y, t)u( z, t). (1.1) Here, t is a time, u( x, t) is a distribution function (e.g., the particle density in a system), In a general case, the method under consideration is applicable for n-dimensional space, x = (x 1 , x 2 , . . . , x n ) = (x i ) ∈ R n . The nonlinearity parameter κ and the small diffusion parameter D are introduced explicitly for the sake of convenience. The n-dimensional Laplace operator in the Cartesian space x ∈ R n is denoted by ∆ x . The coefficients a( x, t) and b( x, y, z, t) are smooth functions of their spatial arguments that grow not faster than polynomially at each point t.
In the physical two-dimensional or three-dimensional space, Equation (1.1) is considered as a model of the optical metal vapor active medium (MVAM) excited by an electrical discharge [12]. The MVAM is a mixture of a buffer inert gas and metal vapors in a gas discharge tube (GDT) (see [13,14] and references therein). In the active medium excited by an electrical discharge, the ionization and recombination processes are mainly caused by the inelastic electron impact. For typical pressures of a buffer gas and metal vapors, preferentially metal atoms are ionized in the mixture. The process of triple recombination of an ion with two electrons is responsible for the deionization (see, e.g., [15]). Such dense plasma formed by metal ions and electrons can be considered as quasineutral. The contracted electrical discharge generates ions and electrons localized in the neighborhood of the GDT center. It means that the concentration of the charges rapidly decreases with the distance from the GDT center. In [12,16], the description of the plasma kinetics under assumptions made was based on the following equation: ∂ t n i = D a (t)∆ x n i + q i n e n neut − q tr n i (n e ) 2 , (1.2) where x are Cartesian coordinates of a point in R 2 or R 3 depending on the problem statement.
The quantity q i = q i ( x, t) is the kinetic coefficient for the electron impact ionization of neutral atoms with a concentration n neut = n neut ( x, t) (n e = n e ( x, t) is an electron concentration).
In the same sense, the coefficient q tr = q tr ( x, t) meets the process of triple recombination of ions with a concentration n i = n i ( x, t). We assume the plasma to be dense so that the triple recombination dominates over the dielectronic recombination. The coefficient D a (t) is an ambipolar diffusion coefficient. The dependence of q i , n neut , q tr , and D a on x and t is due to their dependence on electron temperature that substantially depends on x and t.
Assuming the quasineutrality of plasma, the concentrations of ions and electrons are the same, i.e., Then, for given a( x, t) = q i ( x, t)n neut ( x, t) and q tr ( x, t), Equation (1.2) becomes closed and determines the concentration n i ( x, t) for the given initial and boundary conditions.
To apply the method of semiclassical asymptotics borrowed from papers [9,12], we write Equation (1.2) in the nonlocal form (1.1). For the space R 2 or R 3 , the function b( x, y, z, t) is the probability density of a triple recombination due to the collision of an ion with two where D is the asymptotic small parameter.
In this work, following the method of semiclassical asymptotics [8,9], we have constructed where Φ( x, t, D) is a generic element of the class P t D ; ∆ x = x − X(t, D); the real function ϕ( η, t, D) belongs to the Schwartz space S in variables η, smoothly depends on t, and regularly depends on √ D as D → 0. The real smooth functions S(t, D) and X(t, D), characterizing the class P t D , regularly depend on √ D as D → 0 and are to be determined.
Note that the approach proposed can be useful for other models based on nonlinear equations similar to (1.2). Nonlinear kinetic equations arise in various areas such as cosmology models (see review [17]), superfluidity models [18], etc.
In the next section, we expound the main ideas of the our approach and basic notations, netic equation based on the nonlinear superposition principle is proposed. In Section VI, the specific physically motivated example of the two-dimensional kinetic equation is considered. We illustrate the general formalism of our semiclassical approach by constructing the evolution of the initial ion distribution in the relaxing kinetically enhanced active medium.
In Section VII, we conclude with some remarks.

PROBLEM SOLUTION
In this section, we recapitulate the general scheme for the method of constructing the leading term of semiclassical asymptotics for Equation (1.1). The detailed description of this method can be found in [12].
According to [12], the following asymptotic estimates hold for functions u( x, t, D) from the class P t D : ·, · is the scalar product of vectors, and k is the non-negative integer. In particular, (2.1) For simplicity, we will omit the parameter D in expressions where it does not cause confusion.
In view of the estimates (2.1), the asymptotic expansion of the coefficient in Equation ( in the class P t D , where M ≥ 2 is the highest power of ∆x i accounted in the expansion. The leading term of asymptotics of the solution to (1.1) is determined by the expansion of up to O(D 3/2 ). Following [12], the respective expansions in matrix notations can be written as Here, X = X(t) = X(t, D), ∆ x = x − X(t), ∆ y = y − X(t), and z = z − X(t) are column vectors; (·) is a transposed matrix; a x , b x , b y , and b z are row vectors of , row vectors b y , b z have the analogous form; a xx and b xx are symmetric matrices of the form a xx = ∂ 2 a ∂x i ∂x j x= X(t) , , matrices b yy , b zz , b xy , b yx , b zy ,b yz , and b xz have the analogous form.
The key point of the considered approach [12] is that the nonlocal nonlinearity enters into the approximate kinetic equation obtained with the help of an asymptotic expansion of (1.1) in the form of the moments of the solution u( x, t, D), and dynamical equations that determine the evolution of these moments can be solved separately.
In order to construct the leading term of asymptotics, we need corresponding moments of up to the second order that are defined as follows Here, α where I n is the identity matrix of size n, and the expression a xx α (2) u implies a product of matrices a xx and a (2) u . Let us consider Equation (2.5) as the system where the aggregate of moments is substituted for the set of independent variables σ(t), x(t), α (2) (t) , and α (2) = α (2) ij (t) that are not related to the function u( x, t) in the general case. Thus, the resulting system can be treated as an independent dynamical system. According to [12], this system is termed the Einstein-Ehrenfest (EE) system of the second order for the nonlocal kinetic Equation Let the general solution of this system be where C is a set of arbitrary integration constants. Then, the substitution of the expansion (2.2) into Equation (1.1) with the replacement of moments (2.6) by the general solution (2.7) yields the following linear equation: where the linear operatorL( x, t, C) is given bŷ Equations (2.8) and (2.9) in [12] are termed the associated linear equation for Equation (2.10) Next, we impose a restriction on the integration constants C involved in the general solutions (2.7). The restriction is given by the following algebraic condition: which yields C = C ϕ . Here, g ϕ is the aggregate of the moments (2.6) that is determined by the initial condition ϕ( x, D) as (2.12) Let us consider the Cauchy problem for the associated linear Equations (2.8) and (2.9) with According to [9,12], the solution of the Cauchy problem for Equation (1.1) with the initial condition (2.10) and the solutions of the Cauchy problems (2.8), (2.13) in the class of TCF are related as follows: where C ϕ is given by the algebraic condition (2.11).

III. SEMICLASSICAL ASYMPTOTICS IN A TWO-DIMENSIONAL PLANE-PARALLEL CASE
In this section, based on the method proposed in the previous section, we will obtain an explicit expression for a family of asymptotic solutions for the special case of Equation (1.1) analogous to the one considered in [12].
Let us consider the problem in the plane orthogonal to the GDT axis. Let x = (x 1 , x 2 ) ∈ R 2 be Cartesian coordinates in this plane and the coefficients in (1.1) be given by where the functions a(t) and b(t) are assumed to be monotone decreasing and increasing, respectively, and the parameter µ characterizes the nonlocality of the nonlinearity kernel For functions (3.1), Equations (2.4) and (2.5) yield˙ X(t) = 0. The identity˙ X = 0 also holds in a more general case for the problem with the symmetric configuration of a GDT.
We choose the origin of coordinates so that X(t) t=0 = 0. It leads to X(t) = 0. Ions are usually localized on the GDT axis, which is taken as the origin of coordinates in our case.
In the case under consideration, Equation (2.5) readṡ Here, we denoted where d is a diagonal 2 × 2-matrix, diagonal elements of the matrix 2Dd characterize the degree of localization (dispersion) of the initial axial ion distribution with respect to x 1 , and x 2 , I 2 is the identity 2 × 2-matrix.
The general solution of Equation (3.2) is given by Thus, the relation (3.2) for α (2) for σ(t) = σ(t, C) yield the general solution of the EE system with arbitrary integration constant C: where the set of integration constants reads Let us proceed to the construction of the family of particular solutions for the kinetic equation. For the two-dimensional case under consideration, x ∈ R 2 , with the coefficients (3.1), Equation (1.1) can be written as where ∆ x is a two-dimensional Laplace operator. In view of (2.1) and (2.2), Equation (3.7) u (t) of the required solution u( x, t) by the general solutions of the dynamical EE system (3.5).
In view of (3.2), the relation (2.9) reads which allows us to write the associated linear Equation (2.8) aŝ (3.10) We are looking for the particular solutions of Equation (3.10) in the form of the Gaus- where N 0 is a normalization constant that is related to the initial number of ions. The mul- The substitution of (3.11) into (3.10) yieldṡ (3.12) Equation (3.12) determines the functions S(t, D) and φ(t, D) through quadratures: (3.13) The matrix Riccati equation in (3.12) can be represented as the linear system by the substitution: where B(t) and C(t) are nondegenerate matrices that satisfy the following matrix linear system of differential equations: Since the coefficients in the system (3.15) are scalar, the matrices B(t) and C(t) can be sought in the diagonal form: where W i (t) and Z i (t), i, j = 1, 2, are scalar functions.
Note that the matrix Q(t) is also diagonal is this case: For the functions in (4.1), the system (3.15) readṡ The equations for functions W 1 (t), Z 1 (t) and W 2 (t), Z 2 (t) are identical. Hence, these functions can differ only due to their different initial condition. Let us consider the system of equations for the functions W (t) and Z(t), whose solutions for different initial conditions determine the matrices B(t) and C(t) and, correspondingly, the matrix Q(t). The system (4.4) was coined "the variational system" in [9].
The system (4.4) has two linearly independent solutions. Let us denote them by the following formulae: These solutions are determined by the following initial conditions: For the functions W ± i (t), Z ± i (t) , i, j = 1, 2, the numbers β i that determine the initial conditions by can be different in a general case. Now, we can construct the family of solutions to Equation (3.10) based on the particular solution (3.11). For this purpose, we use the well-known quantum mechanics method widely used in various problems [19][20][21].
Denote the two-dimensional symplectic identity matrix as J = and the skew scalar product which is a conserved quantity in view of (4.4): (4.10) Therefore, for t = 0, we have i.e., Define the operatorsâ (4.13) These operators satisfy the following commutation relations: (4.14) We define the normalization multiplier as N a = 1 √ 2βD . Then, we havê For the diagonal matrix Q(t) (4.2), the solution (3.11) can be written as , and the functions W For the two-dimensional case, the operators (4.15) read a (−) 18) and the commutators are as follows: It can be shown that the operatorsâ  (4.20) and that operatorsâ (±) (t) commute with the operatorL( x, t, C) of Equation (3.10): L ( x, t, C),â Here, n = (n 1 , n 2 ) is the two-dimensional multi-index, n i ∈ {0, 1, 2, . . .}, |n| = n 1 + n 2 , n! = n 1 !n 2 !, N n are normalization coefficients.
The solutions to (4.22) can be written in the explicit form. Define In view of the formula for the Hermitian polynomials [22] the relation (4.22) reads v n ( x, t, C) = N n N 0 (4.25) The set (4.25) is the parametric family of solutions to the associated linear Equation Next, we obtain the constants (5.2) in an explicit form. Let us write the solutions to In order to obtain the moments for the whole set of solutions v n ( x, t, C), we use the representation of the Hermitian polynomial through the generating function [22]: Then, the generating function V ( x, t, t 1 , t 2 , C) for solutions (5.3) is given as follows, t n 1 1 t n 2 2 n! .  u (t) · σ u (t) are linear functionals with respect to u( x, t). This property allows us to obtain moments σ n (t) and α Introduce the following functions: ii,n (t)σ n (t) t n 1 1 t n 2 2 n! .
Here, the subscripts ii indicate the number of a matrix element.
Straightforward calculations of integral in the formulae (5.6) yield (5.7) The formula for A 22 (t, t 1 , t 2 ) can be obtained from one for A 11 (t, t 1 , t 2 ) (5.7) by the formal interchanging 1 ↔ 2 in all subscripts. In view of (5.6), the expansion of functions (5.7) in powers of t 1 , t 2 yields the following expression for moments: 11,(2n 1 ,2n 2 ) (t) = −D Note that functions (5.8) are particular solutions to the EE system (3.2). Hence, the solutions to the EE system (3.2) are determined by the solutions to the variational system (4.3). It is a corollary of the fact shown in [12] that the leading term of asymptotics for the function u( x, t) uniquely defines the functions σ(t), α (2) (t) within accuracy of O(D 3/2 ).
Substituting t = 0 into the formulae (5.8) and taking into account (4.7) and (3.13), we obtain the following initial conditions for moments included in the integration constants C n : while the integration constants (5.2) read The functions v n ( x, t) (5.3) can be written as v n ( x, t) = Υ n (t)Ψ n ( x, t), n = (n 1 , n 2 ), x i .
The substitution of t = 0 into (5.11) yields It can be seen that the following orthogonality condition holds for the functions ψ n ( x): Let us expand the initial condition (2.13) in functions ψ n ( x): k n ψ n ( x). (5.14) Then, the initial condition ϕ( x) corresponds to the following solution to the associated linear

(5.15)
In order to obtain the asymptotic solution of the original kinetic Equation (3.7), we must impose the algebraic condition C = C ϕ . The integration constants C ϕ for the function ϕ( x) (5.14) are given by where i = 1, 2.
Note that our functions v n ( x, t, C) contain two free positive parameters, β 1 and β 2 .
The change of these parameters yields us a new family of the asymptotic solutions ( where γ 1 > γ 2 > 0, 0 ≤ ε ≤ 1. For such constraints for the parameters γ 1 , γ 2 , ε, and N > 0, the function ϕ( x) does not take on negative values. If ε > γ 2 γ 1 , then the initial distribution (6.1) has a minimum at the center x = 0. The minimum of the ion distribution described by Equation (3.7) is physically realizable by the addition of hydrogen into the metal vapor active medium (by creating the so-called kinetically enhanced active medium) [24][25][26].
Let us apply the nonlinear superposition principle to the initial condition (6.1). Note that Since the function φ( x) (6.1) has the symmetry x 1 ↔ x 2 , we put In view of (5.14) and (6.3), the coefficient k n in (5.14) is given by The integral in (6.4) yields 2n = (2n 1 , 2n 2 ). (6.5) Let the value of the parameter β be considered optimal when the coefficients k 2n converge to zero, as |n| → ∞ is the most rapid. It corresponds to the minimum value of the follow- ing function: The function f (β) reaches its minimum at to the the most rapid convergence can be more complicated in the general case. However, it is not crucial to obtain the exact value of β i corresponding to the least |k n | for large |n|.
Hence, it can be done using some approximations for k n (e.g., [27]).
The initial conditions for the EE system (3.2) can be obtained either by the substitution of the coefficients (6.5) into the formulae (5.16) or by the substitution of the initial condition (6.1) into the relations (2.12). Both ways yield (6.8) The EE system (3.2) is integrable for the coefficients in Equation (3.7) of the form The general solutions to the EE system for the coefficients in (6.9) were obtained in [12].
This case is treated as the model of the plasma relaxation. Let us illustrate the solutions corresponding to the initial condition (6.1) for this case. Note that the variational system (4.4) is not integrable for such coefficients as far as we know. It is a remarkable fact, since the analytical solutions of the EE system, which are quite cumbersome though, can be expressed in terms of the solutions of the variational system by analogy with (5.8). For our example, we use the solution to the EE system from [12] and construct the solutions to the variational system (4.4) numerically. Note that the variational system (4.4) is the system of linear differential equations with constant-sign coefficients. The subsequent calculations are presented for κ = 1, Figure 1 shows the evolution of σ(t) for two values of ε. The physical meaning of this function is the total number of ions in the active medium. Figure 2 shows the solutions of the variational system. We provide the plot of these solutions just for one value of ε, since the difference between these plots for ε = 0.85 and ε = 1 is barely perceptible due to the little difference in σ(t). The sign of W (−) (t) was reverted for compactness of the figure. The values of the expansion coefficients are presented in Table I. Note that they tend to zero rapidly as |n| increases due to the optimal choice of β and they equal zero for odd n 1 or n 2 . Figure 3 shows the asymptotic evolution u( x, t) of the initial ion distribution (6.1) constructed as the proposed nonlinear superposition in the weak diffusion approximation with an accuracy of O(D 3/2 ).  Note that Z (+) (t) does not tend to nonzero constant. In a general case, this constant can be negative, which leads to the zero value of Z (+) (t) at some point t = t 0 . At this point, the condition of nondegeneracy of the matrix C(t) (4.1) is violated. It is worthy of discussion how it affects the solutions v n ( x, t, C n ). It can be shown that the functions (5.8) have the removable discontinuity at this point. The same is true for the solutions in (4.25). Hence, the asymptotic solutions u( x, t) = v n ( x, t, C n ) regularly depend on t in a neighborhood of Table I. The values of the coefficient k (n 1 ,n 2 ) for n 1 , n 2 = 0, 1, 2, 3, 4 and two values of ε. It means that if we construct the germ r 2 t from the point t = 0 to a point left of the point t = t 0 and then we construct the germr 2 t after the point t = t 0 , the resulting asymptotic solutions u( x, t), t < t 0 , andũ( x, t), t > t 0 , generated by the germ r 2 t andr 2 t , respectively, determine the continuous asymptotic evolution of the initial state u( x, t) t=0 for both t < t 0 and t > t 0 . Therefore, the formal constraint of the nondegeneracy of the matrix C(t) (4.1) can be bypassed. Note that the error of the semiclassical solutions usually grows over time. It can result in an adverse effect when the semiclassical solution jumps from the asymptotics for one exact solution to another at a sufficiently large time (see the work [28]). Therefore, the estimate of the period of time where the asymptotics are valid is the subject to study separately and is to be obtained in every particular case. Figure 1 illustrates that the ion distribution with the less pronounced minimum at the GDT center relaxes faster. However, Figure 3 shows that both distributions tend to the Gauss-like profile over time. It means that their evolution will be similar on large times. Hence, the difference is significant only at the initial stage of the relaxation.
Note that our method is applicable for the linear case (κ = 0). For κ = 0 and a( x, t) = a(t), the method proposed yields the exact solution to the kinetic equation.

VII. CONCLUSIONS
Within the framework of the approach presented in [12], we have constructed a countable family of asymptotic solutions to the two-dimensional kinetic Equation nonlocal cubic nonlinearity. The approach proposed can be generalized for n-dimensional space. We have considered the two-dimensional problem for the sake of simplicity in order to demonstrate our method and to give its physical interpretation. The constructed asymptotic solutions correspond to the weak diffusion approximation. The approach is based on the solutions to the Cauchy problem to the nonlinear dynamical system (the Einstein-Ehrenfest system) (2.5) and (3.2). These solutions generate the linear parabolic Equations (2.9) and (3.10) associated with the nonlinear kinetic equation. Solutions of this associated linear equation subject to the algebraic condition (2.11) yield the asymptotic solutions (4.25) and (5.1) to the original nonlinear kinetic equation. The family of the asymptotic solutions is constructed based on the set of skew-orthogonal, linearly indepen-dent solutions to the system of linear ODEs (4.4) (the variational system). Such asymptotic solutions form the orthogonal basis for the nonlinear superposition principle that allows one to construct the asymptotic evolution of arbitrary initial distribution ϕ( x) (2.13) expanded into the series with respect to the basis formed by the functions (4.25). The nonlinearity of the superposition principle is caused by the necessity of solving the nonlinear dynamical system (2.5) and (3.2) with initial conditions (2.12) determined by the initial distribution.
Note that the functions v n ( x, t, C) with C determined from (5.16) form an orthogonal basis for the nonlinear superposition principle, while the functions v n ( x, t, C n ) with C n given by (5.9) and (5.10), which determine independent asymptotic solutions to the nonlinear kinetic equation, are not orthogonal in the general case.
This work extends the results of our work [12] where the asymptotic evolution operator It allows us to study the properties of the asymptotic solutions through the properties of the variational system (4.4) that generates the Maslov germ r 2 t on the zero-dimensional manifold Λ 0 t and the symmetry operators. Moreover, we have found the relation (5.8) between the solutions to the Einstein-Ehrenfest system and the variational system. It means that our approach can yield analytical solutions expressed in terms of the solutions to the system of the linear ODEs (4.4) only. However, in the considered specific case, the nonlinear dynamical system admits the analytical solutions, while the variational system (4.4) is solved numerically. Note that the asymptotic evolution operator in [12] and the nonlinear superposition principle yield the exact solutions to the associated linear Equation (3.10), i.e., they yield the same asymptotic solutions for the given initial condition. However, the approach proposed here broadens the possibilities of analysis for these solutions. In addition, the nonlinear superposition principle leads to the expressions with free parameters that can be chosen so that the functional series (5.15) converges rapidly. Hence, in some cases, the study of the solutions (5.15) can be practically reduced to the study of the first few terms of the series. It makes the analysis of such solutions even easier, especially when the solutions are numerical-analytical.
The future prospects for this work are related to the study of a more complex twocomponent kinetic model. We plan to generalize the approach [29] for the latter problem.