Theoretical and Numerical Study of Self-Organizing Processes in a Closed System Classical Oscillator and Random Environment

: A self-organizing joint system classical oscillator–random environment is considered within the framework of a complex probabilistic process that satisﬁes a Langevin-type stochastic differential equation. Various types of randomness generated by the environment are considered. In the limit of statistical equilibrium (SEq), second-order partial differential equations (PDE) are derived that describe the distribution of classical environmental ﬁelds. The mathematical expectation of the oscillator trajectory is constructed in the form of a functional-integral representation, which, in the SEq limit, is compactiﬁed into a two-dimensional integral representation with an integrand: the solution of the second-order complex PDE. It is proved that the complex PDE in the general case is reduced to two independent PDEs of the second order with spatially deviating arguments. The geometric and topological features of the two-dimensional subspace on which these equations arise are studied in detail. An algorithm for parallel modeling of the problem has been developed.


Introduction
From ancient sources, we know that Pythagoras (l. c. 571-497 BC) and the students of his school considered complex numerical and geometric constructions, as well as complicated natural phenomena. Nonetheless, the metaphorization of the word "system" seems to have been first proposed by Democritus (born c. 460 BCE-died c. 370), which meant the formation of complex bodies from atoms, similar to the formation of words from syllables and syllables from letters. In addition, in ancient Greek philosophy, the "system" characterized the orderliness and integrity of natural objects. Sometime later, Plato (427-347 BC) formulated the thesis that the whole is greater than the sum of its parts. Aristotle (384-322 BC), being in a polemic with Plato, formulated the opposite thesis, saying that the whole can be decomposed and studied separately, and then put back together again without losing anything (see for example [1,2]). In those distant times, the Aristotelian concept became more popular due to its simplicity, and for the next almost 2500 years, all research on the problems of natural science was carried out within the framework of this particular concept. Despite its idealized nature, this concept remained the main method of cognition in science for a long time, stimulated its development and contributed to the creation of a huge number of new technologies. Nevertheless, the triumph of Aristotle's concept at the beginning of the 20th century was the creation of a logically perfect theory-the classical mechanics of closed systems, which after a short time experienced a deep crisis, which led to the development of a more general physical representation of quantum mechanics.
In the 20th century, Einstein and Smoluchowski developed the theory of Brownian particle motion [3,4], which stimulated the emergence of a new approach in the study of nature-the theory of classical and quantum open systems [5].
Recall that an open physical system, interacting with the environment, unlike an isolated system, exchanges mass, energy, information, etc. with it. At present, the science of open systems is intensively developing both theoretically and experimentally, especially in the fields of modern quantum physics, chemical physics, etc. [6][7][8]. Note that the open systems approach is not equivalent to Plato's concept, but it is closer to it in spirit. Nevertheless, as shown by numerous studies, this approach also has serious difficulties that do not allow one to describe a number of important phenomena in nonequilibrium thermodynamics, manyparticle quantum systems, etc. Recall that the main drawback of all representations of open systems is that when describing various physical processes, a certain part of information is inevitably lost, especially when it comes to physical systems under extreme conditions. This is due to ignoring the influence of the system on the medium, which excludes the possibility of formation of a small environment (SE) self-consistent with it, which, most likely, can be considered as its integral part or continuation [9].
The purpose of this study is to combine two opposite concepts of nature cognition, namely, to supplement an open system with its environment and to describe the joint system JS as a closed system. Note that just such a statement of the problem would be equivalent to Plato's concept! In this paper, we demonstrate the implementation of this idea on the example of the well-known problem of a classical oscillator immersed in a random environment and under the action of an external force. Recall that this problem has been studied in sufficient detail within the framework of various models of Brownian motion [10], and its results are widely used in solving a number of important applied problems (see [11]). In particular, a Brownian particle moving in a viscous medium exchanges energy and momentum with the environment, which affects the motion of the particle itself. In other words, such mutual influence leads to the appearance of memory in the Brownian particle, i.e., its behavior becomes dependent on the entire previous history of the process. However, we know that a Markov model, by definition, cannot describe random processes with memory, and therefore, taking into account the entrainment of particles of the medium imparts a non-Markov character to the Brownian motion. To overcome the difficulties of describing non-Markovian Brownian motion, it is necessary to reformulate the standard statement of the problem, making significant changes to the mathematical apparatus. In connection with this, a number of authors have proposed the so-called generalized Langevin equation, in which instead of a resistance force proportional to velocity, an integral operator of the convolution type is used [12]. Note that despite the adaptation of the theory of Brownian motion to emerging new scientific and technical problems, this is the theory of open systems and has all the above limitations and disadvantages. In other words, the standard theories of Brownian motion are unsuitable for describing the properties of systems far from equilibrium or in critical states.
Thus, the solution to the problem lies in the development of a fundamentally new mathematical representation, which makes it possible to study the process of self-organization of the entire system, consisting of finite and infinite subsystems, as a closed system. The article is organized as follows: In Section 2, we present a statement of the problem and derive complex stochastic differential equations describing the movement of fields of a random environment for three different cases: • The oscillator frequency is random and there is no external field; • The oscillator frequency is random and the external field is a regular function; • The frequency of the oscillator is a regular function and the external force is random.
Note that for all three cases, we have obtained the distribution equations for the fields of the environment using the systems of the corresponding stochastic differential equations.
In Section 3, we represent in detail a method for constructing a function space measure. Formulas are derived for the mathematical expectation of the oscillator trajectory in the form of a double integral representation for the three cases indicated above, where the integrand is a solution of the corresponding complex second-order partial differential equation PDE.
In Section 4, we obtain a 4th order algebraic equation with coefficients depending on two variables and time, which generates a two-dimensional compactified subspace on which second-order complex PDEs are defined. The geometric and topological features of a two-dimensional manifold are studied in detail. We prove that the formation of topological manifolds with the first Betti [13] number is less than or equal to 4, depending on the interaction constants of the oscillator with the environment.
In Section 5, we present a complex second-order PDE as a system of two coupled real PDEs and formulate the Neumann initial-boundary value problem for this system. We analyze in detail the PDEs system for the case of symmetry or asymmetry as well as the absence of any symmetry of the desired solutions. Using the symmetry properties of the problem, we proved that in the case of symmetry or antisymmetry of solutions, the system of equations reduces to two independent PDEs. Finally, we have shown that when the solutions do not have a certain symmetry, the PDEs system again reduces to two independent PDEs but with deviant arguments.
In Section 6, we construct the time-dependent Shannon entropy for a classical oscillator as an open system affected by the environment. In the same section, we constructed the generalized Shannon entropy for a closed self-organizing system oscillator and random environment. In Section 7, we develop algorithms for numerical simulation of the PDE system as well as analyze, interpret, and visualize the results of various numerical experiments. In Section 8, we discuss the obtained theoretical and numerical results and outline directions for future research.

Statement of the Problem
The classical action of a one-dimensional oscillator immersed in a random environment can be represented as (see [14]): where t i and t f are the moments of time when the interaction of the oscillator with the environment turns on and off, respectively; in addition, L ẋ(t), x(t), t is the Lagrangian describing the oscillator with a random environment: Recall that {f} and {g} are complex probabilistic processes (stochastic sources or generators), whose properties will be refined below. Obviously, the presence of stochastic generators in the Lagrangian makes the action also stochastic. Despite this, one can still require the fulfillment of the minimization condition: Performing the standard procedure for varying the expression (3), taking into account the conditions δx(t i ) = 0 and δx(t f ) = 0 (see [15]), we obtain the following second-order differential equation describing the motion of a test particle, i.e., classical oscillator in a random environment (thermostat): It is important to note that the randomness of the action S[x, t i , t f ) and, accordingly, the Lagrangian L ẋ(t), x(t), t does not affect the variation procedure, as a result of which a stochastic Equation (4) is found. However, as we well know, the second-order Equation (4) is still undefined, since the stochastic equation must be first order.
For definiteness, we will assume that random generators satisfy the white noise correlation relations: where E ... denotes the expectation of the random variable, the set of random generators { f (r) , g (r) } and { f (i) , g (i) }; then, they characterize, respectively, elastic and inelastic collisions of the oscillator with a random environment, while { g } are sets of constants that describe the powers of these processes.
We consider two different cases: 1. When randomness in a JS generates a complex process {f} = 0, and the second source of the random process is absent {g} ≡ 0, and, accordingly, 2.
when {f} ≡ 0 and randomness in a JS generates the generator {g} = 0, which has a complex character.
In the case when the external force is an arbitrary regular function of time, i.e., F 0 (t) = F(t; {g}) {g}≡ 0 , the solution of Equation (4) can be formally represented as (see [16]): where the symbol " * " denotes the complex conjugation of a function, ξ(t) is the solution of the homogeneous Equation (4), i.e., when F t; {g} ≡ 0, in addition, the following notations are made; Ω − 0 = lim t→−∞ Ω t, {f} = const − and Ω + 0 = lim t→+∞ Ω t, {f} = const + . Note that in the general case, the asymptotic states (in) at t → −∞ and (out) at t → +∞ can be different and, accordingly, const − = const + . Below, for definiteness, we will use the model of the regular frequency Ω 0 (t), which has the following form: where γ, ν > 0 are some constants.

Derivation of Environmental Fields Distribution Equations
Theorem 1. If we assume that Equation (4) for the case F(t; {g}) ≡ 0 reduces to a complex Langevin SDE, and the complex force {f} is a Gauss-Markovian random process (5), then the conditional probability distribution of the environmental fields in the limit of statistical equilibrium will obey the Fokker-Planck-type equation.
Proof. We can represent the solution of the classical oscillator Equation (4) as: where ξ 0 (t) is the solution of the classical oscillatory Equation (4) in the case when the frequency is a regular function of time Ω 0 (t) = Ω(t; {f}) {f}≡ 0 and the external regular force is identically equal to zero F 0 (t) ≡ 0; in addition, t 0 = 0 denotes the time of switching on a random environment. As for the function φ(t), it denotes a complex probabilistic process. Substituting (8) into (4), taking into account the regular equation: we obtain the following non-linear stochastic differential equation (SDE) of Langevin type: where . For further study of the problem, it is convenient to represent the complex probabilistic process φ(t) as a sum of fields: Using Equation (10) and representation (11), we can write the following system of non-linear SDEs [9]: where Note that the environment fields satisfy the following initial conditions: Let us consider the following functional describing the distribution of the conditional probability of fields: Differentiating expression (13) with respect to the time "t", taking into account Equation (10), we obtain: where u t = ∂ t u and u ≡ u(t ), in addition: Taking into account that the vector probabilistic process {f} satisfies the correlation relations (5), we can calculate the second term in expression (14). Using Wick's theorem for an arbitrary functional N u, t; { f })|u , t of the argument { f } (see [17]), we can obtain: Since u 1 (t) and u 2 (t) are stochastic functions, the variational derivatives with respect to the independent random forces f (r) (t) and f (i) (t) can be defined as follows: (17) After carrying out the regularization procedure in the sense of the Fourier expansion, we find its value at the time t = t : f . Considering the equalities in (17) for the conditional probability, we obtain the following Fokker-Planck equation: In Equation (18), the evolution operatorL(u, t) has the form: where k 0 = 4u 1 . In addition, in Equations (18) and (19), the variables u 1 and u 2 denote the coordinates of the distribution of the environmental fields in the state of quasi-equilibrium.
In the case when t = t 0 , the conditional probability P(u, t) ≡ P(u, t|u 0 , t 0 ) describes the distribution of classical fields without taking into account the influence of the oscillator on the environment. An important condition for the exact formulation of the problem is the definition of the type of two-dimensional space, on which, let me remind you, the equation for the distribution of environmental fields is set. The latter, in particular, implies writing the conditional probability Equations (18) and (19) in tensor form.
For simplicity, below, we assume that the Fokker-Planck Equations (18) and (19) are defined on a two-dimensional Euclidean space and solve this equation as Neumann's initial-boundary value problem [9]. The numerical study of the free fields of the environment P(u 1 , u 2 , t) is carried out using the mathematical algorithm-difference Equation (87) developed in Listing 1 (Section 7). To illustrate the calculations, graphs of the distribution of fields for various media depending on time are plotted (see Figures 1-3 of Section 7.1). Now, consider the case when the frequency of the oscillator is regular Ω 0 (t), while the external force, on the contrary, is random and can be represented as the sum: where F (r) t; g (υ) (t) = F 0 (t) + (r) gḡ (t) and F (i) t; g (υ) (t) = (i) gḡ (t); in addition,ḡ(t) is a real Gauss-Markov random process, which will be clearly defined below. In particular, using the definition (20) and given that the frequency is regular, we can write the Equation (4) as follows:ẍ Theorem 2. If the oscillator trajectory obeys Equation (21), and the random functionḡ(t) satisfies the Gauss-Markov random process (5), then the distribution of fields of the environment in the limit of statistical equilibrium will be described by PDE of the second order, which, in the (out) asymptotic state or in the limit t → +∞, transforms into the PDE of the Fokker-Planck type.
Proof. Let us represent the solution of the Equation (21) in the form: where x 0 (t) is the solution of a classical oscillator with a regular non-stationary frequency under the action of an external non-stationary force and satisfies the following second-order differential equation:ẍ Substituting (22) into Equation (21) and assuming that θ(t) = w 1 (t) + iw 2 (t), we obtain the following system of stochastic integro-differential equations: where As for the functions σ 1 (t) and σ 2 (t), they are singly differentiable, i.e., belong to class L 1 and are represented as: Assuming that the random functionḡ(t) satisfies the white noise correlation relations: we can use a system of stochastic differential Equations (66) and obtain the following equation for the conditional probability of the environmental fields: where the evolution operatorL has the form: in addition: Obviously, if the medium is turned on in the time interval [t i , t f ], then after t ≥ t f , Equation (25) becomes the Fokker-Planck-type equation.

The Measure of the Functional Space
For further analytical constructions of the theory, it is necessary to determine the distance between functions in the functional space R {ξ} or, more precisely, the measure of the function space (see [18,19]). Let the conditional probability P(u, t|u , t ) satisfy the following limiting condition: The latter means that for small intervals, i.e., for ∆t = t − t 1, we can represent the solution of Equations (18) and (19) as: where is the second-rank matrix with elements; 11 = (r) , 22 = (i) and 12 = 21 = 0, while [· · ·] T denotes a vector transposition. Using (28), we can write an explicit expression for the distribution: where the coefficients k 1 and k 2 are defined from expression (15).
In the case when there is no dissipation in the environment, i.e., (i) = 0, the distribution (29) takes the following form: As can be seen from the expression (29), the evolution of the system in the functional space R {ξ} is characterized by a regular shift with a speed K(u, t) against the background of Gaussian fluctuations with the diffusion matrix ij . Concerning to the trajectory u(t) in the space R {ξ} , it is determined by the following equations (see for example [20]): As can be seen from (31), the trajectory is continuous everywhere, i.e., u(t + ∆t) ∆t→ 0 = u(t), but it is nonetheless everywhere non-differentiable due to the presence of the term ∼ ∆t 1/2 . If the time interval is represented as ∆t = t/N, where N → ∞, then expression (29) can be interpreted as the probability of transition from the vector u l (t) to the vector u l+1 (t) during ∆t in the Brownian motion model. Now, we can define the Fokker-Planck measure of the functional space: where dµ(u 0 ) = δ(u 1 − u 1(0) )δ(u 2 − u 2(0) )du 1 du 2 denotes the measure of the initial distribution; in addition, the following notations are made: Note that the measure Dµ(w), which describes the probability of a given trajectory in the functional space R {x} , can be constructed in a similar way using the equation distributions for classical fields of the environment (25).

The Stages of the Expected Trajectory Calculation
We can now rigorously calculate the expected trajectory of the oscillator for the three cases described above.
Definition 1. The functional integral along the random trajectory u 1 (t), u 2 (t), t will be called the mathematical expectation of the trajectory: where Let us consider the case when the oscillator is not subjected by a random force, i.e., F(t; {g}) ≡ 0. Then, the representation for the mathematical expectation of the trajectory, taking into account (8) and (33), will have the following form: Computing the functional integral (34) using the generalized Feynman-Kac theorem [19], one can find the following two-dimensional integral representation for the mathematical expectation of the trajectory: where the function Q(u 1 , u 2 , t) is the solution of the following second-order complex PDE: Since ξ 0 (t 0 ) is a constant, the main role in determining the expectation of the trajectory is played by the function Λ Q (t).
A numerical study of a complex PDE (36) is carried out using the developed algorithm in the form of a system of difference equations (see Listing 2 Section 7). The results of numerical simulation of the function Q(u 1 , u 2 , t) depending on the state of the environment and time are shown in Figures 4-6 of Section 7.2. The mathematical expectation of the trajectory, depending on various parameters and time, was calculated and represented on the graphs (for details see Section 7.3, Figure 7). Now, let us calculate the mathematical expectation of the trajectoryx(t) when the oscillator is acted upon by a regular external force F 0 (t).
Given Equation (6), the trajectory can be formally written as follows: where the symbol · · · denotes functional integration with respect to the Fokker-Planck measure (see expression (32)): The functional integral in (38) can be calculated and brought to a two-dimensional integral representation if we use the following auxiliary relation: The value under the derivative sign can be calculated using the generalized Feynman-Kac theorem: where the function Q λ (u 1 , u 2 , t) is the solution of the following second-order PDE: Differentiating Equation (40) with respect to the parameter λ, we find the following equation: . Now, introducing the notation D(u 1 , u 2 , t) = D λ (u 1 , u 2 , t) λ= 0 , we obtain the following two-dimensional integral representation: where the function D(u 1 , u 2 , t) satisfies the following complex second-order PDE: To solve Equation (40), it is necessary that the solution satisfies the following initial condition: As we can see, Equation (41) includes two external parameters; one of them is the regular function F 0 (t) (external force), and the other is the random oscillator trajectory ξ(t) ∈ L 1 , which is a solution to stochastic Equation (12).
Integrating Equation (41) with respect to the Fokker-Planck measure (32), one can obtain a new equation: whereξ(t) denotes a regular function representing the mathematical expectation of the oscillator trajectory without external influence (see (34)). The ∂ t D term in Equation (43) can be represented as: where M(t) denotes the exponent of the Fokker-Planck measure (32) in the limit N → ∞ when the sum goes into an integral. It is easy to see that the term ∂ t M(t) is a random function and, accordingly, the new averaging does not allow one to obtain a regular equation for the distribution functionD. However, since at large times ∂ t D → 0, then from Equation (43) in the asymptotic t → ∞, we can obtain the following regular stationary equation: Taking into account (45), we can write the expectation of the oscillator trajectory in the external field when the process enters the asymptotic channel (out): Considering that at t → ∞, the ergodic properties of JS will prevail, we expect that the oscillator trajectoryx(t) becomes a regular function of time.
In conclusion, we calculate the mathematical expectation of the oscillator trajectory under the action of an external random force F(t; {g}). Carrying out similar reasoning, we can write the following functional integral for it: Doing a similar calculation in the functional integral (47), we obtain: where Σ (2) w (t) is a two-dimensional manifold, the geometric and topological features of which must be studied specially. In addition, the function x 0 (t 0 ) is a solution to the regular Equation (23), and the function Q 1 (w 1 , w 2 , t) is a solution to the following complex PDE: As we can see, Equation (48) differs significantly from the parabolic complex PDE (36). It can go into the usual complex PDE in the t → ∞ limit when JS comes to statistical equilibrium.

Geometric and Topological Features of a Compactified Space
As we saw in the previous section, in the limit of statistical equilibrium, the functional space R {ξ} compactifies into the two-dimensional manifold. In particular, for a random frequency and no external force, the functional space R {ξ} compactifies into the twodimensional manifold Σ (2) u (t), and for a random external force and a regular frequency, the functional space, respectively, is compactified into another two-dimensional manifold; w (t). Thus, it is obvious that in this case, the JS in the limit of statistical equilibrium is described in three-dimensional space; R 3 w (t) are two-dimensional manifolds. Below, as an example, we will study in detail the topological and geometric features of the manifold Σ    u (t) with a doubly covariant tensor g µν defined on it, which we will call the generalized metric tensor. Theorem 3. If the motion of a dynamical system is described by stochastic differential equations of the Langevin type (12), then in the limit of statistical equilibrium, these equations generate a two-dimensional space with an antisymmetric metric (Riemann-Cartan manifold). (18) and (19) in tensor form (see for example [21]):

Proof. Let us represent Equations
where the following notations are made: u 1 = u 1 and u 2 = u 2 .
To find the elements of the metric tensor, we write the two-dimensional Laplace-Beltrami operator ∇ 2 in explicit form: Comparing (50) with (19) and requiring the corresponding terms in the equations to be equal, we find: As can be seen from (51), the metric tensor of the subspace Σ (2) u (t) is antisymmetric and, therefore, the corresponding geometry is non-commutative. Note that spaces with such properties arise both in mathematics and in quantum physics, and all of them are associated with non-commutative algebras [22]. For a deep understanding of non-commutative spaces-non-commutative versions of vector bundles, connections, curvature, etc.-as a rule, the operator algebra is used [23]. Recall that a non-commutative algebra is an associative algebra in which the multiplication is not commutative; i.e., xy does not always equal yx. More generally, this algebraic structure, in which one of the basic binary operations is not commutative, furthermore allows for additional structures, such as topology or a norm, that can be carried over to a non-commutative functional algebra.
The first feature of the generalized metric is that the non-symmetric part does not contribute to the definition of the length, sinceȳ ij v i v j = 0, and therefore: Recall that the tensor λ ij defines the Euclidean geometry of the plane tangent to the manifold v ∈ Σ Based on this definition and taking into account the antisymmetry property of the off-diagonal element of the metric tensorḡ 12 (u 1 , u 2 , t) = −ḡ 21 (u 1 , u 2 , t), it is easy to obtain two expressions for the cosine of an angle [24]: where (v 1 ) i denotes the projection of the v 1 vector onto the u i axis (see (13)).
After doing some simple calculations, we find: where ∆ϑ = ϑ 2 − ϑ 1 is the Euclidean angle between the vectors v 1 and v 2 . As for the angle δ, it is determined from the following relations: From (55) also follows the important conditions for the Euclidean angles ϑ + = ∆ϑ + δ and ϑ − = ∆ϑ − δ. In particular, it follows from the definition of Euclidean geometry that the angles must satisfy the following constraint conditions: Recall that two different values of the angle ψ + and ψ − between the vectors v 1 = v 1 (u 1 , u 2 ) and v 2 = v 2 (u 1 , u 2 ) (see (55)) depending on the direction of rotationto the right or to the left-is a characteristic peculiarity of Kozyrev's theory [25].
Taking into account the antisymmetry of the metric of the two-dimensional space Σ (2) u (t), it is easy to prove that its Gaussian curvature is equal to zero. However, following Cartan [26,27], one can introduce a generalized linear connection: where Γ i jk = 1 2 λ il λ lj;k + λ lk;j − λ jk;l , (λ lj;k = ∂λ lj /∂ū k ) is the Christoffel symbol and K i jk denotes the curvature tensor formed by the interaction of an oscillator with a random medium. This tensor can be defined as follows: where δ ij denotes the Kronecker symbol. Note that since the tensor K ijk (ū 1 ,ū 2 , t;ȳ) is antisymmetric with respect to the first pair of indices, then the connection G ijk (ū 1 ,ū 2 , t;ȳ) is consistent with the metric. Now, we can write the equation of motion of a quasi-particle or excitation of an environment, which, taking into account the representation (53) will have the following form: where s = λ ij du i du j = dū i dū i andu i = ∂ū i /∂s.
Thus, Equation (58) can be considered as the equation of a geodesic line in space with connection K i jk (ū 1 ,ū 2 , t;ȳ). To solve this equation, it is necessary to know the form of the contortion tensor K i jk (ū 1 ,ū 2 , t;ȳ) as a function of coordinates and time.

Topology of Two-Dimensional Subspace Σ
(2) u (t) Theorem 4. The generalized metric g ij (u 1 , u 2 , t) satisfying (49)-(51) is defined by the antisymmetric metric element g 12 (u 1 , u 2 , t) = −g 21 (u 1 , u 2 , t), which satisfies a fourth-degree algebraic equation generating a topological manifold with possibly non-trivial singularities and a first Betti number less than or equal to 4.

Proof.
Comparing the operators (19) and (50) and taking into account the coordinate transformations (52), we can obtain the following first-order partial differential equations for the antisymmetric element of the metric tensor y = g 12 (u 1 , u 2 , t): where χ = y/(a + y 2 ) and a = (r) (i) . However, as follows from (59), the system of equations is redefined with respect to the desired function y. In this regard, a reasonable question arises: under what conditions are these equations compatible?
From Equation (59), it is easy to find expressions for two different derivatives of the off-diagonal component of the metric tensor: Using Equation (60), we can find the following two expressions for the mixed second derivatives of the metric tensor element y = g 12 (u 1 , u 2 , t): where χ ;j = ∂χ/∂u j and k i;j = ∂k i /∂u j , (i, j = 1, 2). It is important to note that the antisymmetry of the off-diagonal elements of the metric tensor arises at the stage of choosing the coordinate system and, accordingly, the orientation of the considered sub-manifold.
As for the question of the symmetry of mixed second derivatives on any oriented manifolds, then, based on the basic requirement of mathematical analysis, the following identity must hold at any time at each point of the manifold (Schwarz's theorem, see [28]): which is a necessary condition for a twice continuously differentiable function. In the context of partial differential equations, it is called the Schwarz integrability condition. If we write this equality (62) explicitly, taking into account expressions (61), then it will look like this: Finally, given (60) from Equation (63), we obtain the following 4th degree algebraic equation: where the coefficients of the algebraic equation A n (u 1 , u 2 , t) are defined by the expressions: (64) is solved exactly by the Ferrari method and has four solutions, some of which may be complex [29]. Since the coefficients of the Equation (64) are functions of two coordinates (u 1 , u 2 ) and time, the solutions must form a continuum of sets in two-dimensional Euclidean space. We will be interested in those sets of solutions that are complex. Obviously, if we cut out and remove from the Euclidean space all the domains on which the solutions of the algebraic Equation (64) are complex, then the remaining space will have topological singularities. As the numerical solution of the algebraic Equation (64) shows, depending on the parameters ε (r) and ε (i) , the manifold Σ Thus, we have proved the necessary condition for the compatibility of the system of Equation (59). To prove whether this condition is sufficient for the compatibility of the system of Equation (64), it is necessary to substitute the solution of the algebraic Equation (64) into this system. Unfortunately, this assertion verification procedure turns out to be a very difficult task. However, due to the fact that all sequential mathematical constructions are performed analytically by algebraic methods, there is every reason to believe that inverse transformations also take place; i.e., up to linear functions u 1 and u 2 , the function y is a solution to the Equation (59).
To illustrate the above, below are the results of visualization of a series of calculations (see Figures 8-11), which allow obtaining a detailed idea of the topological features of the two-dimensional manifold Σ u (t), which arises after the compactification of the function space R {x} (see Section 7.4 of Section 7 for details). It is also important to note that a similar analysis for the complex Equation (36) describing the solution Q(u 1 , u 2 , t) proves that the exact geometry for solving this problem is the manifold Σ

Statement of the Initial-Boundary Value Problem for the Complex PDE
For definiteness, we will study the initial-boundary value problem in the case when an external force does not act on the oscillator in the thermostat. In this case, the mathematical expectation of the trajectory is described by a two-dimensional integral representation (35), where the function Q(u 1 , u 2 , t) is the solution of a complex second-order PDE (36). We will consider the problem when Σ (2) u (t) is a two-dimensional Euclidean space, that is: Theorem 5. If a two-dimensional complex PDE has the form (36), then using the internal symmetry of the equation, it can be reduced to two independent PDFs belonging to the class of PDFs with a deviated argument given by affine transformations such as reflection.
Proof. Representing the solution of Equation (36) as a sum of real and imaginary parts: one can obtain the following system of PDEs: The functions Q (i) (u 1 , u 2 , t) and Q (r) (u 1 , u 2 , t) can be normalized and given the meaning of the probability density: and, obviously, the functions must satisfy the normalization condition: Q (i) (u 1 , u 2 , t) +Q (r) (u 1 , u 2 , t) du 1 du 2 = 1.
It is easy to see that when changing the coordinates (u 1 , u 2 ) → (u 1 , −u 2 ), the system of Equation (66) becomes: If we assume that Q (i) (u 1 , −u 2 , t) = Q (r) (u 1 , u 2 , t) and, accordingly, Q (r) (u 1 , −u 2 , t) = Q (i) (u 1 , u 2 , t), then the system of Equation (68) takes the original form (66), i.e., the first equation goes over into the second, and the second goes into the first. Using these obvious symmetry properties, we can write the system of coupled PDEs (68) as two independent PDEs: As we can see in the system (69), the equations are independent, and each of them belongs to the PDE class with a spatially deviated argument given by affine transformations such as reflection. By solving one of the equations, we can obtain the solution of the second one using a 180-degree rotation in the two-dimensional Euclidean space R 2 . Before proceeding to the numerical solution of these PDEs, we consider three possible scenarios: (a) When the solutions of Equation (69) are even functions with respect to the u 2 coordinate, i.e., Q (i) (u 1 , −u 2 , t) = Q (i) (u 1 , u 2 , t) and Q (r) (u 1 , −u 2 , t) = Q (r) (u 1 , u 2 , t); (b) When the solutions of Equation (69) are odd functions with respect to the coordinate u 2 , i.e., Q (i) (u 1 , −u 2 , t) = −Q (i) (u 1 , u 2 , t) and Q (r) (u 1 , −u 2 , t) = −Q (r) (u 1 , u 2 , t), and, accordingly, the case; (c) When the indicated functions do not have definite parity. In the first (a) case from (69), we obtain two unrelated PDEs: In the second (b) case, we again obtain a system of uncoupled differential equations, only in this case, it is necessary to replace Q (i) (u 1 , u 2 , t) → Q (r) (u 1 , u 2 , t) in the first equation and Q (r) (u 1 , u 2 , t) → Q (i) (u 1 , u 2 , t) in the second one, respectively.
In the third (c) case, the functions Q (i) (u 1 , u 2 , t) and Q (r) (u 1 , u 2 , t) are described by the system of Equation (69). Note that this is the most general and difficult case for numerical simulation, which will be considered in detail below.
We will consider a more complicated case where the PDE is a deviant argument. Our task will be to formulate an initial-boundary value problem for solving one of the PDEs of the equations system (69).
For definiteness, let us consider the second equation in (69). We can represent it as the following system: Recall that the second equation in (71) is obtained from the first one by replacing u 2 → −u 2 .
As an initial condition, we assume that the probability distribution of the environmental fields Q (r) (u 1 , u 2 , t) at time t 0 is described by the Dirac delta function: As for the boundary conditions, we define the Neumann boundary conditions on the perpendicular axes u 1 and u 2 , respectively: ;2 (u 1 , 0, t), where Q (r) ; i = ∂Q (r) /∂u i , (i = 1, 2). We will determine their values based on a number of physical considerations. As we will see below, the conditions (73) lead to two different equations.
First, we consider the behavior of the second equation in (69) near the u 1 ∈ (−∞, +∞) axis. The solution of the equation near this axis can be represented as: where a 0 , a 1 and a 2 are some unknown constants that will be determined based on physical considerations. Substituting (74) into the second equation in (69), in the limit of u 2 → 0, one obtains the following second-order PDE: In particular, for the first Neumann condition, we obtain: Since the equation for the function Q (r) (u 1 , u 2 , t) is not symmetric with respect to the change u 2 → −u 2 , the first Neumann boundary condition cannot be equal to zero and, accordingly, a 1 = 0. In addition, the function Q (r) ;2 (u 1 , t) ∼ Q (r) 1 (u 1 , t) has the meaning of the probability density on the u 1 axis; we can normalize it to one. The latter means that a 1 = 1 and the normalization of the solution Q (r) 1 (u 1 , t) is equal to one, that is: where c 0 (t) is the normalization constant.
Proceeding from the fact that the boundary conditions on the perpendicular axes u 1 and u 2 describe the probability distributions of elastic and inelastic collisions independent of each other, it is natural to assume that the term of inelastic collision in Equation (75) should be identically equal to zero. In other words, we can require that the equality 2a 2 /a 0 − 1 = 0 becomes satisfied, and thus, Equation (75) can be written as: To solve this equation, it is necessary to formulate the following initial-boundary value problem: and, correspondingly, where s 1 and s 2 denote points far enough from the origin 0, which are located to the left and right, respectively. Now, to establish the second Neumann boundary condition, we consider the solution Q (r) (u 1 , u 2 , t) near the axis u 2 ∈ (−∞, +∞). To do this, we represent the solution in the form: where b 0 , b 1 and b 2 are some constants that we will define below. Substituting (80) into (69), in the limit of u 1 → 0, we obtain the following equation: Since only inelastic processes are taken into account on the u 2 axis, we can require that the following condition becomes fulfilled: 2b 2 /b 0 − 1 = 0. In addition, for definiteness, we can set the constant b 0 = 1 and b 1 = 1/2. Taking into account the clarifications made, Equation (82) can be simplified and presented as: Because (82) is a second-order PDE with a deviant argument, it can be solved together with the same equation but after replacing the argument u 2 → −u 2 . For Equation (82), we can formulate the following initial-boundary conditions: and, correspondingly, Considering the above, for the second boundary condition, we can write the following expression: As in the case of the solution Q (r) 1 (u 1 , t), we can interpret the solution Q (r) 2 (u 2 , to) as the density probability and normalize it to one.

Entropy of a Self-Organizing System
As known, for a classical dynamical system, an important characteristic is the nonstationary Shannon entropy [30]. In particular, the entropy production rate is a quantitative measure of a non-equilibrium processes, and knowledge of its quantity indicates information about the dissipated heat [31,32], the difference in free energy between two equilibrium states [33,34], and also about the efficiency, if the considered non-equilibrium system is an engine [35][36][37]. It should be noted that the rate of entropy production provides important information for systems with hidden degrees of freedom [38,39], as well as for interacting subsystems, where the amount of information plays a key role [40][41][42][43].
In philosophy, physics and mathematics, the term negentropy is often used, which has a negative value and, therefore, has the opposite to entropy sense. Note that if entropy characterizes the measure of orderliness and organization of the system, then negentropy is the possibility of reducing entropy or an effort toward order. Recall that this concept was first proposed by Schrödinger [44] when explaining the behavior of living systems: in order not to die, a living system struggles with the surrounding chaos and with the entropy it produces, organizing and ordering the latter by introducing negentropy. This, in particular, explains the behavior of self-organizing systems.
For definiteness, let us consider the question of the change in the entropy of a classical oscillator in the case when its frequency has a random component. We first calculate the dynamics of the oscillator entropy without taking into account its influence on the random environment, when, as in the second case, we will consider this influence consistently and strictly. In particular, in the first case, we can define the non-stationary entropy in the standard way: It is often important to know the change in entropy or entropy production over a certain period of time: Since processes of different nature occur in the problem under consideration, we can also introduce the concept of partial entropy, which characterizes the production of entropy of a particular process: Note that the partial entropies S    Finally, we can determine the entropy of a closed self-organizing system classical oscillator + thermostat. In this case, the expression for the generalized Shannon entropy S gen (t) can be represented as: As can be seen from the graphs (see Figure 11), at intermediate times, the behavior of the entropies calculated by the formulas (84) and (86) differ greatly in value and in character. Moreover, in some time intervals, the generalized entropy (86), regardless of the intensity of the processes in the environment, takes a negative value, which is quite natural for a self-organizing system. Finally, as follows from the simulation and the corresponding visualizations, the behavior of both ordinary and generalized entropy tends to a constant value in the limit of large times.
Details of entropy calculations are discussed in Sub-section E of Section 7.

Numerical Methods for Solving the Problem
Numerical simulation of the self-organization processes in the joint system (classical oscillator and thermostat), even for a simple case, i.e., in the absence of an external field F t; {g} ≡ 0, requires large computational resources. This is due to the fact that the complex probabilistic processes described by the three key PDEs, Equations (18) and (19), as well as the system of Equation (66), are numerically difficult solvable tasks. Recall that in [9], we have already considered a PDE of this type. In the same work, we considered various finite-difference methods of solution. Taking into account the analysis performed and test calculations for the numerical solution, we chose an explicit finite-difference scheme of the second order of accuracy in coordinates and the first order in time. Despite the external simplicity of this method, we believe that this scheme satisfies the goals of our work in terms of efficiency and accuracy. Note that Equations (18), (19) and (66) are second-order PDEs of the parabolic type with a source term. The main difficulty that we face in the calculations arises in connection with the convective transfer, which must be taken into account when increasing the size of the computational domain.
Below, we consider a numerical algorithm for solving the initial-boundary value problems for these PDEs in the two-dimensional Euclidean space R 2 , as shown in Listings 1 and 2. For a PDE system (88) with appropriate conditions on the coordinate axes, the finite difference scheme is similar to Listing 1.
Note that in this work, the proposed difference scheme for a linear two-dimensional equation with a source term, which can be attributed to the parabolization schemes for equations of the convection-diffusion type, has no independent meaning. The stability conditions for such schemes have been considered in sufficient detail in various monographs; see, for example, [45]. Let us present the necessary stability conditions for scheme (87). The conditions for scheme (88) look similar. Let us rewrite expression (87) in the following form, ignoring the source term 4∆tP n j,k : P n+1 j,k = P n j,k + r 1 P n j+1,k − 2P n j,k + P n j−1,k + r 2 P n j,k+1 − 2P n j,k + P n j,k−1 + C 1 P n j+1,k − P n j−1,k + C 2 P n 1,k+1 − P n j,k−1 , Listing 1. Numerical algorithm for solving PDEs (18) and (19) after coordinate transformations (52).

1.
The continuous area R 2 for Equations (18) and (19) is replaced by the calculated area: In the computational domain, a uniform difference grid is set in time t and in spatial coordinatesū 1 andū 2 : where ∆ū 1 and ∆ū 2 are steps in spatial coordinates, and ∆t is a step in time.

2.
In these notations, for the constructed grid, we have the following difference equation for the PDEs (18) and (19): P n+1 j,k = P n j, k + r 1 P n j+1,k − 2P n j,k + P n j−1,k + r 2 P n j,k+1 − 2P n j,k + P n j,k−1 + 4∆tP n j,k + ∆t where the following notations are introduced:

3.
To calculate Equations (18) and (19), it is also necessary to set two boundary conditions in the form of difference equations on the coordinate axesū 1 andū 2 , respectively, which can be easily found by approximating Equation (87). Note that these difference equations must be solved taking into account the Dirichlet boundary conditions; P l (x)| x= ± |s| = 0, (|s| 1), where the index (l = 1, 2) denotes the first and second boundary conditions, respectively. 4.
The condition is set at the center of the coordinate axes: In addition, the Dirichlet condition P(ū 1 ,ū 2 , t) ∂G = 0 is specified on the boundaries of the computational domain, where ∂G denotes the boundary.

5.
As an initial condition, instead of the Dirac delta function, we use the Gaussian distribution: Note that the parameters σ = 500, ω = πσa 2 and a = b = 1/2 included in the function were chosen in such a way to normalize the initial distribution to unit.
In this case, the following inequalities must hold: Such necessary conditions are quite sufficient in determining the upper limit of the time step ∆t crit : ∆t ≤ ∆t crit for all options presented in the table below. To estimate the sufficient stability condition, we can use the approximation taken from the one-dimensional case: In real test calculations, the original differential equation behaves more complicated than the model linear equation. However, the real ∆t crit obtained from the calculations is close to the value determined from the given difference scheme stability inequalities.
Formally, we can represent the system of Equations (87) and (88) in matrix form, which in the general case looks like: where A l denotes the transition matrix for the variable f , in which connection A 1 is the transition matrix for the variable f n+1 j,k 1 = P (see Equation (87)), and the matrices A 2 and A 3 are for the variables f n+1 j,k 2 = Q (i) and f n+1 j,k 3 = Q (r) (see Equation (88)). To save space, we do not write out these matrices in detail. It is important to note that the calculation process itself is clear from Equations (87) and (88).

Listing 2.
Numerical algorithm for solving the system of PDEs (52) after coordinate transformation.

1.
The continuous region R 2 for the PDEs system (66) is replaced by a discrete grid, as described in Listing 1.

3.
Similarly, as in Listing 1, for the solutions Q (i)ū 1 ,ū 2 , t) and Q (r) (ū 1 ,ū 2 , t), boundary conditions are set on theū 1 andū 2 axes in the form of difference equations, which can be obtained by approximating Equation (88) on the same axes. Note that we solve each equation obtained for the boundary conditions as an internal Dirichlet problem with a zero value of the solution at the boundary.

4.
As in the case of the probability density P(ū 1 ,ū 2 , t) (see Listing 1), the solutions Q (i) (ū 1 ,ū 2 , t) and Q (r) (ū 1 ,ū 2 , t) are subject to similar conditions at the center of the coordinate axes: In addition, we will assume that at the boundary of the computational domain, the solutions are subject to the following conditions:
Note that each line of this table characterizes different states of the environment: a. When the processes going on in the environment, both elastic and inelastic, are weak; b. When elastic processes are strong and inelastic processes are weak; c. When both elastic and inelastic processes are strong in the environment. Figure 1. The different stages of the evolution of the free fields of the environment, respectively, at t 1 = 1.5, t 2 = 10 and t 3 = 20. Note that these distributions were calculated using the data from the first row of the table, which corresponds to weak elastic and inelastic processes in the environment.
Comparing the distributions at different times, it is easy to see that as t ∼ 10, the distribution P(u 1 , u 2 , t) normalized to unity is established or tends to its stationary limit.

Figure 2.
A series of graphs of the distribution of free fields of the environment at time points t 1 = 1.5, t 2 = 10, t 3 = 20. Recall that we used the data of the second line table, which characterizes strongly elastic and weakly inelastic processes occurring in the environment. An analysis of the graphs shows that the distribution tends to the stationary limit already at t ∼ 7.

Figure 3.
Distributions of free fields of the environment, respectively, at the time points t 1 = 1.5, t 2 = 3, t 3 = 10 and t 4 = 20. The data of the third row table, corresponding to strong elastic and inelastic processes occurring in the environment, were used for the calculation. As can be seen from the figures, the greater the constants that determine the powers of elastic and inelastic processes, the faster the distribution of environmental fields is established.

Distributions of Environmental Fields Taking into Account the Influence of the Oscillator
We now present figures illustrating the evolution of the normalized functionsQ (r) (u 1 , u 2 , t) andQ (i) (u 1 , u 2 , t) depending on time. Using the algorithm developed in Listing 1, we can calculate and visualize all of these solutions. Below are the normalized distributions for three different cases.   In particular, as we see from these figures, with an increase in the constants of interactions with the environment, the time for establishing distributions is reduced.

Mathematical Expectation of the Oscillator Trajectory
Using all the above calculations, one can simulate the expected value of the oscillator trajectory in the absence of an external field based on Equation (35). We modeled the expected oscillator trajectory for both its real and imaginary parts for five different environmental states and visualized their behavior as a function of time (see Figure 7). In particular, as can be seen from the figures, in all the cases under consideration, the mathematical expectation of the trajectory not only has a non-trivial oscillatory character, but it decreases with time, taking both positive and negative values. Mathematical expectations of the real (blue) and imaginary (red) parts of the oscillator trajectory for five sets of parameters. Left to right, front row: (ε (r) = ε (i) = 0.01), (ε (r) = 1, ε (i) = 0.001) and (ε (r) = 1, ε (i) = 0.01); in the second row: (ε (r) = 1, ε (i) = 0.5) and (ε (r) = ε (i) = 1).

Calculation of Topological and Geometric Features of the Manifold Σ
(2) u (t) As we saw above, the off-diagonal term of the metric tensor y = g 12 (u 1 , u 2 , t) of a manifold Σ  (see data of the table), we obtain sets of surfaces for the offdiagonal term of the metric tensor (see , which allows us to study and understand the geometric and topological features of the manifold Σ (2) u (t). An analysis of the surfaces (see Figure 8) shows that when the oscillator is immersed in an environment with weak elastic and inelastic interactions, these surfaces do not have interesting topological features. However, there is an obvious singularity at the point u 1 = 0 − , (u 1 ∈ (−∞, 0] \ 0, > 0), because as → 0, the term of y = g 12 (u 1 , u 2 , t) tends to plus infinity in the upper half-space when as in the lower half-space, the term y = −g 21 (u 1 , u 2 , t) tends to minus infinity. It should be noted that in the case under consideration, the evolution of the environment characteristically does not change the geometry of space but only shifts the minimum point of the surface.  u (t) from the asymptotic state (in) with the environment data ε (r) = ε (i) = 0.01 and frequency Ω − 0 = 1, to the asymptotic state (out) with frequency Ω + 0 = 3. Figure 9. The first row on the right shows a topological space with the first Betti number equal to 2 in the (in) asymptotic state, which in the process of evolution to the (out) asymptotic state transits into two simply connected half-spaces with displaced holes. The second row shows three two-dimensional pictures-projections of the manifold Σ (2) u (t) onto the plane u 1 , u 2 , characterizing the evolution of the topological singularities of the manifold when transitioning between these asymptotic states defined by the data (Ω − 0 = Ω 0 (t) = 1 Ω 0 (t) = 2, Ω + 0 = Ω 0 (t) = 3) and (ε (r) = 1, ε (i) = 0.01), respectively.
In the presence of strong elastic and weak inelastic processes in the environment, the calculations lead to the formation of the following surfaces (see Figure 9). An analysis of the pictures shows that the arising manifold is topological space with the first Betti number is equal to 2, which during evolution transits to the manifold consisting from two sub-manifolds each of which is a simply connected topological space with spatially shifted typological singularities. Moreover, the shift between the topologies of the upper and lower sub-manifolds is present already in the (in) asymptotic state. In the course of evolution, this shift only increases, while the gap intersection decreases and already in the (out) state becomes equal to zero. Moreover, the shift between the topologies of the upper Σ (2)+ u (t) and lower Σ (2)− u (t) sub-manifolds is present in the (in) asymptotic state (this can be verified by analyzing three-dimensional graphs in Figure 9). In the course of evolution, this shift only increases, while the gap intersection decreases and already in the (out) state becomes equal to zero. In other words, the manifold Σ (2) u (t) can be represented as the union of two sub-manifolds: where Σ (2)+ u (t) denotes a submanifold of the upper half-space, while Σ (2)− u (t) is a submanifold of the lower half-space. As we will see below, this manifold can evolve into a disjoint union of two sub-manifolds in the course of evolution.  u (t) from the asymptotic state (in) with data Ω − 0 = 1, ε (r) = 1, ε (i) = 0.5 to the (out) asymptotic state with data Ω + 0 = 3, ε (r) = 1, ε (i) = 0.5 . The second row represents two-dimensional projections of this manifold onto the (u 1 , u 2 ) plane, showing the evolution of its topological features. (2)+ u (t) from the asymptotic state (in) with data Ω − 0 = 1, ε (r) = 1, ε (i) = 0.5 , to (out) an asymptotic state where the frequency is Ω + 0 = 3. In the second row (from left to right), there are figures showing a similar evolution for the sub-manifold Σ Finally, let us consider the case when strong elastic and non-weak inelastic processes occur in the environment. As calculations show (see Figure 10), in this case, the manifold also has the first Betti nuber equal to 2 but with different sizes and slot configurations. In particular, the first row in Figure 10 shows three-dimensional visualizations of the elements of the metric tensor g 12 (u 1 , u 2 , t) and g 21 (u 1 , u 2 , t) as functions of coordinates (u 1 , u 2 ) in asymptotic states (in) and (out), respectively. The second row ( Figure 10) shows the graphs showing the evolution of the topological features of the manifold Σ   u (t) is the union of two sub-manifolds (see (89)), a natural question arises: do these sub-manifolds retain any of their topological singularities during evolution? As the calculations and their two-dimensional visualizations ( Figure 11) show, each of the sub-manifolds (see the first and second rows of the graphs) in the process of evolution passes from a topological space with the first Betti number equal to 2 to a simply connected topology. However, as we can see, these topologies are displaced and do not have a common hole (see the graphs in the second line of Figure 10, obtained by combining the graphs of the first two lines of Figure 11). It is important to note that depending on the parameters, a topological manifold with the first Betti number equal to 3 also arises, which is a union of two oriented topological submanifolds Σ (2)+ u (t) and Σ (2)− u (t) with the first Betti number equal to 3. In the end, we note that the manifold Σ (2) u (t), regardless of the state of the environment, at the end of its evolution in the state (out) is a disjoint union of two sub-manifolds.

Features of Entropy Calculation
The integrals in the equations for the ordinary (84) and generalized (86) entropies were calculated using the trapezoidal method.
However, due to the fact that in the expressions for the entropy, there is a term of type A ln(A), which has a singularity at A = 0, the condition A ln(A) = 0 is introduced, which makes it possible to eliminate this singularity. We performed entropy calculations using Formulas (84) and (86) for three different states of the environment and visualized them (see Figure 12) taking into account the data of Table 1. As follows from these graphs, the usual entropy S(t) (left column) in the first two cases continuously increases with time and reaches a constant value in the (out) state. In the third case, when inelastic precessions in the environment are strong enough, as we observe, from some moment, the increase in entropy passes into the stage of its decrease, and already at the long times, when the system goes into (out) state, it takes a constant value. The right column contains graphs of the generalized entropy S gen (t), which, in addition to being non-monotonic, also contain areas with negative values. This behavior of entropy is quite explainable for a closed self-organizing system, which at some point can generate negentropy to stabilize the state of the joint system. In particular, the negative value of entropy can also be explained by the short-term capture or "swallowing" of a small or oscillator subsystem by a large subsystem or thermostat. As can be seen from the graphs of the right column, in the first two cases, the generalized entropy S gen (t) relatively fast tends a constant value, while in the third case, the process of stabilization of the combined system proceeds non-monotonically and takes a long time.

Conclusions
The main achievement of this work is the development of a mathematically rigorous representation that allows one to study the statistical properties of a classical oscillator with its random environment as a problem of self-organization of a closed self-consistent system. Note that such a statement in the philosophical sense corresponds to the consideration of the problem within the framework of Plato's concept, which excludes the loss of information in the JS. We have carried out a mathematical implementation of this idea within the framework of a complex probabilistic process that satisfies the SDE of the Langevin type. The paper considers three typical scenarios of a random environment or thermostat, for which, in the limit of statistical equilibrium, the kinetic equations for the distribution of the fields of the environment are derived (see Equations (18), (19), (25) and (26)). With the help of these equations, the measures of functional spaces are determined, which make it possible in the future to build mathematical expectations of the corresponding physical parameters of the problem. We use the generalized Feynman-Kac theorem [18,19] to compactify the infinite-dimensional functional integral describing the expectation of the oscillator trajectory and reduce it to the two-dimensional integral representation, where the integrand is the solution of the complex second-order PDE given on a two-dimensional manifold Σ  u (t) is generally described by a non-commutative geometry, which may also have important geometric and topological features. In the case of non-intensive random processes in the environment, i.e., when (ε (r) , ε (i) ) 1, the manifold Σ (2)− u (t) are Euclidean subspaces with one singular boundary (see Figure 8). As the power of random processes in the environment (ε (r) , ε (i) ) ∼ 1 increases, both the geometric properties of the manifold Σ (2)− u (t) and their topological features change strongly. In particular, as shown in Figure 10, the sub-manifolds Σ (2)+ u (t) and Σ (2)− u (t) have a non-Euclidean curvilinear geometry and topologies with the first Betti number (see for example [13]) equal to 2, which, in the (out) state, go over to a simply connected topology. In other words, conformational transformations of an additional subspace Σ (2) u (t) of a self-organizing classical system, taking into account its geometric and typological features, lead to radical differences in the description of the dynamics of a classical system without an environment and when it is immersed in a random environment. However, it is obvious that the correct formulation is to solve the problem not on the two-dimensional Euclidean space R 2 but on the manifold Σ (2) u (t), which is the union of two topological sub-manifolds Σ (2)+ u (t) and Σ (2)− u (t), respectively. We emphasize once again that this is due to the fact that the generated manifold Σ (2) u (t) is generally described by a non-commutative geometry that has non-trivial topological singularities. In this paper, we develop an efficient mathematical algorithm for the sequential and parallel calculation of various characteristic parameters of the problem, taking into account that the manifold Σ (2) u (t) is the Euclidean space. Note that this approximation takes place when the interaction constants of the oscillator with the environment are small.
In the near future, we plan to generalize the computational algorithm for performing calculations on just such a manifold. Recall that in this case, the distribution of the environmental fields will be described by the tensor Equations (49) and (50) where the off-diagonal element g 12 (u 1 , u 2 , t) = −g 21 (u 1 , u 2 , t) will be determined by the algebraic equation of the 4th degree (64). The latter will allow numerical methods to study important features of the dynamics of a classical system within the framework of an ideologically more consistent and accurate Platonic concept.
In conclusion, it should be noted that the study of quantum analogues of the considered classical models, taking into account the non-commutativity of the emerging geometries, will be extremely interesting and rich in new and unexpected results.