Full Text

Interaction

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 \(\Delta t\). 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.

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.

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 \(c\). 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 \(\Delta t \to 0\), is usually avoided.

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 \(C\) 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 \(L^2\) 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 \(\delta\)-functions. We do not have a theorem on the existence of a solution for such problems.

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.

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.

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.

Short-range interaction mathematical model

Let there be \(N\) identical bodies of mass \(m\), under the assumption of short-range interaction between them, they produce a field with potential \(u\) 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

\[\label{eq:ode}
m\ddot{\vec r}_n = - \nabla u\big|_{\vec r = \vec r_n}, \quad n=1,2,\dots, N,\]

and the field dynamics is described by a wave equation

\[\label{eq:pde}
\frac{1}{c^2}\frac{\partial^2u}{\partial t^2} = \Delta u + \rho,\]

where \(\rho\) is the density of mass distribution:

\[\label{eq:rho}
\rho = \gamma \sum \limits_{n=1}^N \delta_s (\vec r - \vec r_n).\]

Here it is reasonable to consider \(\delta\) as a smoothed prototype of Dirac delta function that tends to the delta function in the limit \(s \to 0\).

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 \(c \to \infty\), and then to the limit \(s \to 0\).

[th:1] Let there be a family of solutions to the system –, parameterized by two parameters \(c\geqslant 0\) and \(s>0\), and let it satisfy the condition \[u, u_t \in L^2(\mathbb R^3)\] at \(t=0\). If we first proceed to the limit \(c \to \infty\), and then to the limit \(s \to 0\), then this solution becomes the solution of the classical many-body problem.

According to the Poisson formula

\[\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}\]

Under the assumptions made about the initial conditions, the last two terms tend to zero as \(c \to \infty\) and we get \[u(\vec r, t) = \frac{1}{4\pi} \iiint \limits_{\mathbb{R}^3}
\frac{\rho(\vec r', t)}{|\vec r-\vec r'|}dv'.\]

Substituting Eq. here yields an expression that, at \(s\to 0\), becomes \[u = \gamma \sum \limits_{n=1}^N \frac{1}{|\vec r - \vec r_n| }.\]

However, we cannot substitute it directly into because this would lead to dividing by zero. However, for \(s \not =0\), the expression for \(u\) is a sum of terms of the form \(\phi_s(\vec r -\vec r_n)\), having an extreme at \(\vec r=\vec r_n\). So \[\left.\nabla \phi_s(\vec r -\vec r_n)\right|_{\vec r=\vec r_n}=0\] and there is no division by zero: \[\nabla u \big|_{\vec r=\vec r_m} = \gamma \sum \limits_{n\not=m} \nabla
\phi_s(\vec r_m -\vec r_n).\]

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 \[m\ddot{\vec r}_n = - \gamma \nabla_{\vec r_n} \sum \limits_{m \not =n}
\frac{1}{|\vec r_n - \vec r_m| }.\]

The proved theorem allows us to hope that for large \(c\) 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.

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.

Let the bodies occupy fixed positions up to \(t<0\), then for \(t<0\) we know \(u\) as a solution to the Poisson equation \[\Delta u =- \rho.\]

At the moment \(t=0\) the bodies are given initial velocities. Adding the initial condition \(u_t=0\) to the wave equation, we get the classical initial value problem for finding the potential \(u\), if we assume that the density \(\rho\) is known.

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 \(G\) and set the conditions \(u\big|_{\partial G}=0\) on its boundary. This box will replace the boundary conditions at infinity.

The complete problem is formulated as follows. Given initial positions \(\vec r_n^{(0)}\) and initial velocities \(\vec v_n^{(0)}\) of the bodies, the solution is calculated to the boundary value problem

\[\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}\]

It is required to find the functions \(\vec r_n(t)\) and \(u(x,y,z,t)\) satisfying the initial and boundary value problem:

\[\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.\]

with the initial conditions \[\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)\] and the boundary conditions \(
u\big|_{\partial G} =0.
\)

We believe that this problem has a unique solution for small \(t\). 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.

Galerkin method

A natural method for solving the oscillation equation in a finite domain is the Galerkin method . Let \(\phi_n\) be the normalized eigenfunctions of the Laplace operator in \(G\), and let \(\alpha_n^2\) be the corresponding eigenvalues. We seek the solution of the wave equation in the form

\[\label{eq:series}
u(x,y,z,t)= \sum \limits_{j=1}^\infty u_j(t) \phi_j(x,y,z),\]

where \(u_j\) are coefficients yet unknown. Then \[\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\] and \[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.\]

If we truncate the sum over \(j\) to any finite number of terms \(J\), then the system has a unique solution, taking into account the initial conditions. In the limit \(s \to 0\) we get

\[\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\]

and

\[\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.\]

The initial conditions for \(\vec r_n\) are given and for \(u_j\) they are found from equation using the explicit formulae \[\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\] or

\[\label{eq:3}
u_j(0)=\frac{\gamma}{\alpha_j^2} \sum \limits_{n=1}^N \phi_j\big|_{\vec r = \vec r_n}\]

and

\[\label{eq:4}
\dot u_j(0)=0.\]

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 .

[fig:1]

For example, let us take the box in the form of a cube \([0,L]^3\). Then the eigenfunctions are expressed as \[\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},\] with the corresponding eigenvalues \[\alpha_{mnk}^2= \frac{\pi^2}{L^2} (m^2+n^2+k^2).\]

Taking the first \(J\) 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.

Let us take, for example, \( c=1\), \(m=1\), \(\gamma=1\), \(L=10\) and consider the problem of two bodies. We place the first body at the point \((6,5,5)\) and the second one at the point \((4,5,5)\). Let the first body be at rest, and the second one have an initial velocity \(\vec v_2 =(0,1,0)\) In the classical case, this leads to the rotation of bodies along ellipses around their center of gravity \((5,5,5)\), and the motion occurs in the \(xy\) 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 \(Oz\) axis, the motion still remains flat, only the plane itself changes. Therefore, our system is rich in integrals of motion.

Conclusion

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 \(G\). The transformation of the box \(G \) into a resonator seems quite natural in the framework of the theory of short-range interaction.

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.

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 \(c \to \infty\) (long-range interaction), and then to the limit \(s \to 0\) (narrowing the charge density to \(\delta \)-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 \(s \to 0\). In this case, the limit \(c \to \infty\) 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.

About the authors

Mark M. Gambaryan

Peoples’ Friendship University of Russia (RUDN University)

Author for correspondence.
Email: gamb.mg@gmail.com
ORCID iD: 0000-0002-4650-4648

PhD student of Department of Applied Probability and Informatics

6, Miklukho-Maklaya St., Moscow, 117198, Russian Federation

Mikhail D. Malykh

Peoples’ Friendship University of Russia (RUDN University); Meshcheryakov Laboratory of Information Technologies Joint Institute for Nuclear Research

Email: malykh-md@rudn.ru
ORCID iD: 0000-0001-6541-6603

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

6, Miklukho-Maklaya St., Moscow, 117198, Russian Federation; 6, Joliot-Curie St., Dubna, Moscow Region, 141980, Russian Federation

References