Discrete and Continuous Models and Applied Computational ScienceDiscrete and Continuous Models and Applied Computational Science2658-46702658-7149Peoples' Friendship University of Russia named after Patrice Lumumba (RUDN University)3032610.22363/2658-4670-2022-30-1-52-61Research ArticleOn the many-body problem with short-range interactionGambaryanMark M.<p>PhD student of Department of Applied Probability and Informatics</p>gamb.mg@gmail.comhttps://orcid.org/0000-0002-4650-4648MalykhMikhail D.<p>Doctor of Physical and Mathematical Sciences, Assistant Professor of Department of Applied Probability and Informatics of Peoples’ Friendship University of Russia (RUDN University); Researcher in Meshcheryakov Laboratory of Information Technologies, Joint Institute for Nuclear Research</p>malykh-md@rudn.ruhttps://orcid.org/0000-0001-6541-6603Peoples’ Friendship University of Russia (RUDN University)Meshcheryakov Laboratory of Information Technologies Joint Institute for Nuclear Research01042022301526125022022Copyright © 2022, Gambaryan M.M., Malykh M.D.2022<p style="text-align: justify;">The classical problem of the interaction of charged particles is considered in the framework of the concept of short-range interaction. Difficulties in the mathematical description of short-range interaction are discussed, for which it is necessary to combine two models, a nonlinear dynamic system describing the motion of particles in a field, and a boundary value problem for a hyperbolic equation or Maxwell’s equations describing the field. Attention is paid to the averaging procedure, that is, the transition from the positions of particles and their velocities to the charge and current densities. The problem is shown to contain several parameters; when they tend to zero in a strictly defined order, the model turns into the classical many-body problem. According to the Galerkin method, the problem is reduced to a dynamic system in which the equations describing the dynamics of particles, are added to the equations describing the oscillations of a field in a box. This problem is a simplification, different from that leading to classical mechanics. It is proposed to be considered as the simplest mathematical model describing the many-body problem with short-range interaction. This model consists of the equations of motion for particles, supplemented with equations that describe the natural oscillations of the field in the box. The results of the first computer experiments with this short-range interaction model are presented. It is shown that this model is rich in conservation laws.</p>many-body problemGalerkin methodshort-range interactionзадача многих телметод Галёркинаблизкодействие<h1 id="interaction">Interaction</h1>
<p>Studying the motion of a beam of charged particles in an external electromagnetic field with the interaction of particles taken into account is one of the most important and popular problems in plasma electronics. In the framework of the generally accepted approach to its study , the time interval is divided into discrete steps of length <span class="math">\(\Delta t\)</span>. At each step, based on the current positions of the charges and their velocities that define the currents, the induced field is calculated as a solution of Maxwell’s equations. Then this ‘induced’ field is added to the external field and the new positions and velocities of the particles are calculated, which they acquire in this field under the action of the Lorentz force.</p>
<p>The described scheme allows for many variations . However, these details do not at all remove the division into processes: for one time step, first, the charges and their velocities generate a field, and then the field acts on the bodies through the Lorentz force.</p>
<p>It is quite obvious that what has been said gives a description of a numerical method for studying a certain mathematical model of the many-body problem with short-range interaction. The latter is explicitly taken into account in the model: at each step, the field is calculated and the interaction between particles is carried out through this field, which is described using the Maxwell’s equations, that is, hyperbolic equations that describe the propagation of signals with the speed of light <span class="math">\(c\)</span>. However, the model itself remains undescribed; moreover, the issue of the convergence of the described numerical method, i.e., the study of the limit at <span class="math">\(\Delta t \to 0\)</span>, is usually avoided.</p>
<p>From a mathematical point of view, it is necessary to combine two models into one system: a nonlinear dynamic system that describes the motion of charges based on the Lorentz law, and a linear system of Maxwell’s equations that describes the dynamics of the electromagnetic field. Separately, these models are well studied. Dynamical systems with analytic right-hand sides are solved in analytic functions, and the convergence of the finite-difference method is proved in the <span class="math">\(C\)</span> norm . Maxwell’s equations, as well as linear partial differential equations in general, are naturally solved in Sobolev spaces, and an approximate solution is also sought in one or another integral norm, for example, in <span class="math">\(L^2\)</span> over the space . However, when combining these models, we must consider dynamical systems, the right-hand sides of which are elements of Sobolev spaces, and Maxwell’s equations, in which currents and charges are combinations of <span class="math">\(\delta\)</span>-functions. We do not have a theorem on the existence of a solution for such problems.</p>
<p>A detailed description of the model, separated from the numerical method of its study, is very useful, firstly, in order to be able to assess the quality of the study in terms of closeness to the exact solution, and not in terms of closeness to the expectations of the experimenters. Secondly, good mathematical models always have a large number of symmetries, which correspond to conservation laws. Checking their performance provides another important criterion for assessing the quality of the numerical method. Finally, it cannot be ruled out that less obvious, but more effective numerical methods for studying this model can be found.</p>
<p>Thus, for example, by means of computer experiments it was found that the Boris difference scheme for solving the equations of motion corresponds to the expectations of experimenters more than others . To explain this effect, Hong Qin et al. showed that this scheme is the phase volume when integrating the equations of motion of one particle in an external electromagnetic field. The question of whether the Boris scheme inherits the properties of the original system in the many-many problem, which, we note, is not Hamiltonian, was not raised.</p>
<p>In this paper, we consider the simplest formulation of the many-body problem with short-range interaction described by the wave equation. To add boundary conditions to the wave equation, we consider the problem in a finite domain. The question of setting the radiation conditions in such a problem does not seem trivial to us, although to simplify the problem it is usually assumed that the field in the far zone should be equal to zero.</p>
<h1 id="short-range-interaction-mathematical-model">Short-range interaction mathematical model</h1>
<p>Let there be <span class="math">\(N\)</span> identical bodies of mass <span class="math">\(m\)</span>, under the assumption of short-range interaction between them, they produce a field with potential <span class="math">\(u\)</span> and move in it in accordance with the second Newton’s law. The simplest formulation can be written as follows: the dynamics of particles is described by the set of equations</p>
<p><span class="math">\[\label{eq:ode}
m\ddot{\vec r}_n = - \nabla u\big|_{\vec r = \vec r_n}, \quad n=1,2,\dots, N,\]</span></p>
<p>and the field dynamics is described by a wave equation</p>
<p><span class="math">\[\label{eq:pde}
\frac{1}{c^2}\frac{\partial^2u}{\partial t^2} = \Delta u + \rho,\]</span></p>
<p>where <span class="math">\(\rho\)</span> is the density of mass distribution:</p>
<p><span class="math">\[\label{eq:rho}
\rho = \gamma \sum \limits_{n=1}^N \delta_s (\vec r - \vec r_n).\]</span></p>
<p>Here it is reasonable to consider <span class="math">\(\delta\)</span> as a smoothed prototype of Dirac delta function that tends to the delta function in the limit <span class="math">\(s \to 0\)</span>.</p>
<p>By virtue of the Poisson formula and regardless of the boundary conditions imposed on the field, this problem becomes classical if we first proceed to the limit <span class="math">\(c \to \infty\)</span>, and then to the limit <span class="math">\(s \to 0\)</span>.</p>
<p>[th:1] Let there be a family of solutions to the system –, parameterized by two parameters <span class="math">\(c\geqslant 0\)</span> and <span class="math">\(s0\)</span>, and let it satisfy the condition <span class="math">\[u, u_t \in L^2(\mathbb R^3)\]</span> at <span class="math">\(t=0\)</span>. If we first proceed to the limit <span class="math">\(c \to \infty\)</span>, and then to the limit <span class="math">\(s \to 0\)</span>, then this solution becomes the solution of the classical many-body problem.</p>
<p>According to the Poisson formula</p>
<p><span class="math">\[\begin{gathered}
u(\vec r, t) = \frac{1}{4\pi} \iiint \limits_{|\vec r|ct}
\frac{\rho(\vec r', t-|\vec r-\vec r'|/c)}{|\vec r-\vec r'|}dv' + \\
+\frac{1}{4\pi c} \frac{\partial}{\partial t} \iint \limits_{|\vec
r|=ct} \frac{u(\vec r',0)}{|\vec r-\vec r'|}ds' +\frac{1}{4\pi c}
\iint \limits_{|\vec r|=ct} \frac{u_t(\vec r',0)}{|\vec r-\vec
r'|}ds'. \end{gathered}\]</span></p>
<p>Under the assumptions made about the initial conditions, the last two terms tend to zero as <span class="math">\(c \to \infty\)</span> and we get <span class="math">\[u(\vec r, t) = \frac{1}{4\pi} \iiint \limits_{\mathbb{R}^3}
\frac{\rho(\vec r', t)}{|\vec r-\vec r'|}dv'.\]</span></p>
<p>Substituting Eq. here yields an expression that, at <span class="math">\(s\to 0\)</span>, becomes <span class="math">\[u = \gamma \sum \limits_{n=1}^N \frac{1}{|\vec r - \vec r_n| }.\]</span></p>
<p>However, we cannot substitute it directly into because this would lead to dividing by zero. However, for <span class="math">\(s \not =0\)</span>, the expression for <span class="math">\(u\)</span> is a sum of terms of the form <span class="math">\(\phi_s(\vec r -\vec r_n)\)</span>, having an extreme at <span class="math">\(\vec r=\vec r_n\)</span>. So <span class="math">\[\left.\nabla \phi_s(\vec r -\vec r_n)\right|_{\vec r=\vec r_n}=0\]</span> and there is no division by zero: <span class="math">\[\nabla u \big|_{\vec r=\vec r_m} = \gamma \sum \limits_{n\not=m} \nabla
\phi_s(\vec r_m -\vec r_n).\]</span></p>
<p>Now, proceeding to the limit, we get in the right-hand side of equation exactly an expression that should be in the many-body problem <span class="math">\[m\ddot{\vec r}_n = - \gamma \nabla_{\vec r_n} \sum \limits_{m \not =n}
\frac{1}{|\vec r_n - \vec r_m| }.\]</span></p>
<p>The proved theorem allows us to hope that for large <span class="math">\(c\)</span> the solutions of the system under consideration resemble the classical many-body problem. However, it is important to emphasize that the order of proceeding to the limit is important.</p>
<p>We are interested in constructing a model of many-body motion, in which short-range interaction is explicitly taken into account, rather than in the classical limit itself. For our purpose, it is necessary to supplement the differential equations with initial and boundary conditions.</p>
<p>Let the bodies occupy fixed positions up to <span class="math">\(t0\)</span>, then for <span class="math">\(t0\)</span> we know <span class="math">\(u\)</span> as a solution to the Poisson equation <span class="math">\[\Delta u =- \rho.\]</span></p>
<p>At the moment <span class="math">\(t=0\)</span> the bodies are given initial velocities. Adding the initial condition <span class="math">\(u_t=0\)</span> to the wave equation, we get the classical initial value problem for finding the potential <span class="math">\(u\)</span>, if we assume that the density <span class="math">\(\rho\)</span> is known.</p>
<p>Let us turn to the boundary conditions. We assume that the bodies do not radiate waves that are noticeable in the far zone. To treat this problem numerically, we place the systems in a Dirichlet box <span class="math">\(G\)</span> and set the conditions <span class="math">\(u\big|_{\partial G}=0\)</span> on its boundary. This box will replace the boundary conditions at infinity.</p>
<p>The complete problem is formulated as follows. Given initial positions <span class="math">\(\vec r_n^{(0)}\)</span> and initial velocities <span class="math">\(\vec v_n^{(0)}\)</span> of the bodies, the solution is calculated to the boundary value problem</p>
<p><span class="math">\[\label{eq:u0}
\begin{cases}
\Delta u_0 = -\gamma \sum \limits_{n=1}^N \delta_s (\vec r - \vec r_n^{(0)}), \\
u\big|_{\partial G}=0.
\end{cases}\]</span></p>
<p>It is required to find the functions <span class="math">\(\vec r_n(t)\)</span> and <span class="math">\(u(x,y,z,t)\)</span> satisfying the initial and boundary value problem:</p>
<p><span class="math">\[\label{eq:sys}
\left\{
\begin{aligned}
m\ddot{\vec r}_n = - \nabla u\big|_{\vec r = \vec r_n}, \quad n=1,2,\dots, N,\\
\frac{1}{c^2}\frac{\partial^2u}{\partial t^2} = \Delta u + \gamma
\sum \limits_{n=1}^N \delta_s (\vec r - \vec r_n),
\end{aligned}
\right.\]</span></p>
<p>with the initial conditions <span class="math">\[\vec r_n=\vec r_n^{(0)}, \quad \dot{\vec r}_n=\vec v_n^{(0)}, \quad
u=u_0, \quad u_t=0 \quad (t=0)\]</span> and the boundary conditions <span class="math">\(
u\big|_{\partial G} =0.
\)</span></p>
<p>We believe that this problem has a unique solution for small <span class="math">\(t\)</span>. However, the proof of this assertion requires a more careful description of the class of functions in which the solution is sought. We confine ourselves to a few computer experiments with this model.</p>
<h1 id="galerkin-method">Galerkin method</h1>
<p>A natural method for solving the oscillation equation in a finite domain is the Galerkin method . Let <span class="math">\(\phi_n\)</span> be the normalized eigenfunctions of the Laplace operator in <span class="math">\(G\)</span>, and let <span class="math">\(\alpha_n^2\)</span> be the corresponding eigenvalues. We seek the solution of the wave equation in the form</p>
<p><span class="math">\[\label{eq:series}
u(x,y,z,t)= \sum \limits_{j=1}^\infty u_j(t) \phi_j(x,y,z),\]</span></p>
<p>where <span class="math">\(u_j\)</span> are coefficients yet unknown. Then <span class="math">\[\frac{1}{c^2} \frac{d^2 u_j}{dt^2} + \alpha_j^2 u_j = \gamma \sum
\limits_{n=1}^N \iiint \limits_{G} \delta_s (\vec r - \vec r_n(t))
\phi_j dxdydz , \quad j=1,2, \dots\]</span> and <span class="math">\[m\frac{d^2\vec r_n}{dt^2} = - \sum \limits_{j=1}^\infty u_j(t) \nabla
\phi_j \big|_{\vec r = \vec r_n}, \quad n=1,2,\dots, N.\]</span></p>
<p>If we truncate the sum over <span class="math">\(j\)</span> to any finite number of terms <span class="math">\(J\)</span>, then the system has a unique solution, taking into account the initial conditions. Inthelimit <span class="math">\(s \to 0\)</span> we get</p>
<p><span class="math">\[\label{eq:1}
\frac{1}{c^2} \frac{d^2 u_j}{dt^2} + \alpha_j^2 u_j = \gamma \sum \limits_{n=1}^N \phi_j \big|_{\vec r = \vec r_n} , \quad j=1,2, \dots, J\]</span></p>
<p>and</p>
<p><span class="math">\[\label{eq:2}
m\frac{d^2\vec r_n}{dt^2} = - \sum \limits_{j=1}^J u_j(t) \nabla \phi_j \big|_{\vec r = \vec r_n}, \quad n=1,2,\dots, N.\]</span></p>
<p>The initial conditions for <span class="math">\(\vec r_n\)</span> are given and for <span class="math">\(u_j\)</span> they are found from equation using the explicit formulae <span class="math">\[\alpha_j^2 u_j(0) = \gamma \sum \limits_{n=1}^N \iiint \limits_{G} \delta_s(\vec r - \vec r_n) \phi_j dxdydz\]</span> or</p>
<p><span class="math">\[\label{eq:3}
u_j(0)=\frac{\gamma}{\alpha_j^2} \sum \limits_{n=1}^N \phi_j\big|_{\vec r = \vec r_n}\]</span></p>
<p>and</p>
<p><span class="math">\[\label{eq:4}
\dot u_j(0)=0.\]</span></p>
<p>By virtue of the Weyl lemma , the eigenfunctions of the Laplace operator are twice continuously differentiable in the domain considered. Therefore, the system of ordinary differential equations , falls under the conditions of the classical Cauchy theorem. This means that the initial value problem for equations , with initial conditions , has a unique solution, at least in the vicinity of the initial data. Moreover, standard numerical methods can be applied to this problem, for example, the Runge-Kutta method .</p>
<p><img src="plot-0" alt="image"><img title="fig:" src="plot-1" alt="The first body trajectory at c=1 and c=10"> [fig:1]</p>
<p>For example, let us take the box in the form of a cube <span class="math">\([0,L]^3\)</span>. Then the eigenfunctions are expressed as <span class="math">\[\sin \frac{\pi m x}{L} \sin \frac{\pi n y}{L} \sin \frac{\pi k z}{L}, \quad n,m,k \in \mathbb{N},\]</span> with the corresponding eigenvalues <span class="math">\[\alpha_{mnk}^2= \frac{\pi^2}{L^2} (m^2+n^2+k^2).\]</span></p>
<p>Taking the first <span class="math">\(J\)</span> functions from this set, the initial positions and velocities of the bodies, we uniquely determine the initial problem –, which we will solve by the classical Runge-Kutta method of the 4th order.</p>
<p>Let us take, for example, <span class="math">\( c=1\)</span>, <span class="math">\(m=1\)</span>, <span class="math">\(\gamma=1\)</span>, <span class="math">\(L=10\)</span> and consider the problem of two bodies. We place the first body at the point <span class="math">\((6,5,5)\)</span> and the second one at the point <span class="math">\((4,5,5)\)</span>. Let the first body be at rest, and the second one have an initial velocity <span class="math">\(\vec v_2 =(0,1,0)\)</span> In the classical case, this leads to the rotation of bodies along ellipses around their center of gravity <span class="math">\((5,5,5)\)</span>, and the motion occurs in the <span class="math">\(xy\)</span> plane. Our computer experiment shows that in the case of short-range interaction, the motion also turns out to be planar, but instead of ellipses, more complex non-closed curves are obtained. If we set the velocity in the direction of the <span class="math">\(Oz\)</span> axis, the motion still remains flat, only the plane itself changes. Therefore, our system is rich in integrals of motion.</p>
<h1 id="conclusion">Conclusion</h1>
<p>The initial value problem – can and should be considered as a mathematical model describing the many-body problem with short-range interaction. Equation has a very simple physical meaning of a mechanical equation of motion (the second Newton law), and equation describes the natural oscillations of the field in the resonator <span class="math">\(G\)</span>. The transformation of the box <span class="math">\(G \)</span> into a resonator seems quite natural in the framework of the theory of short-range interaction.</p>
<p>A few computer experiments that we have managed to perform demonstrate that this system is rich in conservation laws. However, it is not yet clear to us how to study them analytically. We hope that further experiments with this new problem will clarify the issue.</p>
<p>With respect to the system , , this problem is a simplification, however, a simplification different from that leading to classical mechanics. By virtue of theorem [th:1], we will pass to classical mechanics if we first proceed to the limit <span class="math">\(c \to \infty\)</span> (long-range interaction), and then to the limit <span class="math">\(s \to 0\)</span> (narrowing the charge density to <span class="math">\(\delta \)</span>-functions). When deriving the system –, we restrict the number of oscillations in the box to a finite number of modes (Galerkin method) and immediately proceed to the limit <span class="math">\(s \to 0\)</span>. In this case, the limit <span class="math">\(c \to \infty\)</span> makes the singularity problem perturbed, and, from the point of view of the Tikhonov and Vasilieva theory , slow variables correspond to the bodies, and fast variables correspond to the field.</p>[C. K. Birdsall and L. A. Bruce, Plasma physics via computer simulation. Bristol, Philadelphia and New York: Adam Hilger, 1991.][R. W. Hockney and J. W. Eastwood, Computer simulation using particles. Bristol, Philadelphia and New York: Adam Hilger, 1988.][V. P. Tarakanov, User’s manual for code KARAT. Springfield, VA: Berkley Research, 1999.][E. Hairer, G. Wanner, and S. P. Nørsett, Solving Ordinary Differential Equations, 3rd ed. New York: Springer, 2008, vol. 1.][G. Duvaut and J. L. Lions, Inequalities in Mechanics and Physics. Berlin: Springer-Verlag, 1976.][H. Qin, S. Zhang, J. Xiao, J. Liu, and Y. Sun, “Why is Boris algorithm so good?” Physics of Plasmas, vol. 20, p. 084503, 2013. DOI: 10.1063/1.4818428.][A. G. Sveshnikov, A. N. Bogolyubov, and K. V. V., Lectures on Mathematical Physics [Lektsii po matematicheskoy fizike]. Moscow: MGU, 1993, in Russian.][I. I. Vorovich, “On some direct methods in the nonlinear theory of oscillations of shallow shells [O nekotorykh pryamykh metodakh v nelineynoy teorii kolebaniy pologikh obolochek],” Izvestiya Akademii Nauk USSR, Seriya Matematicheskaya, vol. 21, no. 6, pp. 747-784, 1957, in Russian.][P. G. Ciarlet, The finite element method for elliptic problems. NorthHolland, 1978.][N. G. Afendikova, “The history of Galerkin’s method and its role in M.V. Keldysh’s work [Istoriya metoda Galerkina i yego rol’ v tvorchestve M.V.Keldysha],” Keldysh Institute preprints, no. 77, 2014, in Russian.][G. Hellwig, Partial differential equations. An introduction. Leipzig: Teubner, 1960.][G. Hellwig, Differential operators of Mathematical Physics. Reading, MA: Addison-Wesley, 1967.][A. N. Tikhonov, “Systems of differential equations containing small parameters at derivatives [Sistemy differentsial’nykh uravneniy, soderzhashchiye malyye parametry pri proizvodnykh],” Mat. Sb., vol. 31, no. 3, pp. 575-586, 1952, in Russian.][A. B. Vassilieva and V. F. Butuzov, Asymptotic methods in singular perturbation theory [Asimptoticheskiye metody v teorii singulyarnykh vozmushcheniy]. Moscow: Vysshaya shkola, 1990, in Russian.]