Solution of an ODE upon singular perturbation

Question

The following question originates from a Physics problem, so I apologize if I am not using a suitable mathematical jargon.

The original system involves $N$ massless electric charges at position $\boldsymbol{r}_1$, $\boldsymbol{r}_2$, ..., $\boldsymbol{r}_N$ which can move on a plane pierced by a uniform and constant transverse magnetic field $\boldsymbol{B}$. The motion equation governing the dynamics of the $j$-th charge is: \begin{equation} \label{eq:Motion_eq} \boldsymbol{0}=q_j\,\dot{\boldsymbol{r}}_j\times\boldsymbol{B}\, +\, q_j\boldsymbol{E}(\boldsymbol{r}_j) \end{equation} where $q_j$ is the charge of the $j$-th electric charge and $\boldsymbol{E}(\boldsymbol{r}_j)$ is the electric field generated at position $\boldsymbol{r}_j$ by all the remaining charges. Moreover, at any time, the electric field at position $\boldsymbol{E}(\boldsymbol{r}_j)$ can be written as $$ \boldsymbol{E}(\boldsymbol{r}_j) = \sum_{i=1\\i\neq j}^N q_i\frac{\boldsymbol{r}_j-\boldsymbol{r}_i}{|\boldsymbol{r}_j-\boldsymbol{r}_i|^3}. $$ The motion equation above is manifestly of first order in time and the only initial conditions one needs are the initial positions of the massless charges, i.e. $\boldsymbol{r}_1(t=0)$, $\boldsymbol{r}_2(t=0)$, ..., $\boldsymbol{r}_N(t=0)$.

Now let us switch from massless to massive charges. This means that the aforementioned motion equation transforms as: \begin{equation} M_j \ddot{\boldsymbol{r}}_j=q_j\,\dot{\boldsymbol{r}}_j\times\boldsymbol{B}\, +\, q_j\boldsymbol{E}(\boldsymbol{r}_j) \end{equation} We assume that all the masses $M_j$ are very small. As far as I understand, the introduction of a non-zero mass constitutes a singular perturbation as it manifestly alters the order of the differential equations. Accordingly, the number of initial conditions which one should set doubles.

My question is: if one knows the solution $\{\boldsymbol{r}_1(t),\,\boldsymbol{r}_2(t),\,\dots,\,\boldsymbol{r}_N(t)\}$ of the massless problem (i.e. of the system of first-order differential equations), what can be said about the solution of the massive problem (i.e. of the system of second-order differential equations)?

For simplicity we can assume $M_j=M\,\forall j$ and we can also make additional assumptions, if needed. E.g. that the energy $$H_\mathrm{massive}=T+V=\sum_{j=1}^N\frac{1}{2}M_j\dot{\boldsymbol{r}}_j^2 + V$$ of the massive system is limited and does not significantly depart from that of the massless system $$H_\mathrm{massless}=V=\frac{1}{2}\sum_{j=1}^N\sum_{i\neq j}^N \frac{q_iq_j}{|\boldsymbol{r}_i-\boldsymbol{r}_j|}\,.$$ Another possible and reasonable assumption would be that the initial velocities of the 2nd-order problem deviate only little from the fixed initial velocities of the 1st order problem.

I guess that a term $\mathcal{O}(M^{-1})$ will show up in the solutions of the perturbed system. Is this correct? Apart from this term, is it possible to write a sort of Taylor expansion involving powers of $M$? Is there a general framework to approach this singular-perturbation problem?

Answer 1

Let me change the notations to fit the mathematical literature. I will denote by $x(t) \in \mathbb{R}^{3N}$ the positions of the particles, by $y(t) \in \mathbb{R}^{3N}$ their velocities and by $0 < \varepsilon \ll 1$ their common small mass. Your massive problem can be written as a first order ODE as: $$ \begin{cases} \dot{x} = f(x,y), \\ \varepsilon \dot{y} = g(x,y), \end{cases} $$ with $f(x,y) = y$ and $g(x,y) = q y \times B + q E(x)$. The question now becomes whether solutions to the degenerate system $\dot{x} = f(x,y)$ and $g(x,y) = 0$ correctly approximate solutions with small $\varepsilon > 0$.

This question has been studied by many authors, probably starting with Tikhonov's paper:

Andrei Nikolaevich Tikhonov. "Systems of differential equations containing small parameters in the derivatives." Mat. Sb.(NS) 31.73 (1952): 575-586.

As a starting point, I would suggest the book:

Wolfgang Wasow, Asymptotic expansions for ordinary differential equations. Reprint of the 1976 edition. Dover Publications, Inc., New York, 1987. x+374 pp. ISBN: 0-486-65456-7

In particular, Section 39 contains an introduction to the topic (and surveys Tikhonov's results), and Section 40 discusses series expansions (with respect to $\varepsilon$), which correspond to the framework you are looking for.

Eventually, let me mention two difficulties:

1. As already mentioned in the comments, a key difficulty in such singularly perturbed problems lies in the discrepancy between the boundary conditions, which can lead to the presence of boundary layers in the solutions with $\varepsilon > 0$ (in your case, boundary layers in time, located near $t = 0$). The study of boundary layers is a whole field by itself. If you are not familiar with it, Wikipedia's page on singular perturbations already features a nice example with the ODE $\varepsilon \ddot{x}(t) + \dot{x}(t) = - e^{-t}$.
2. In your physical problem, there is another source of singularity, since the electric field becomes singular as two particles become very close to each other. This difficulty is typically not included in the above mentioned literature, which always requires some kind of regularity on $f$ and $g$ (at least continuous for instance). Note that this difficulty is somehow independent, since it is already present with $\varepsilon = 0$ or $\varepsilon = 1$. In order to extend the known results for the general setting explained above to your case, you would first typically have to derive a priori estimates proving that the particles won't be collapsing together on the considered time interval.

---

Following your request, here are two (hopefully more "practical") examples:

a) Concerning boundary layers, let us look at $\dot{x} = 1$ on $\mathbb{R}$, with initial condition $x(0) = 0$. Then the solution is $x(t) = t$. Note that, in this massless case, $\dot{x}(0) = 1$.

Now the massive problem is $\varepsilon \ddot{x}_\varepsilon + \dot{x}_\varepsilon = 1$, with boundary conditions $x_\varepsilon(0) = 0$ and $\dot{x}_\varepsilon(0) = y^0$ given. The explicit solution is $$ x_\varepsilon(t) = t + \varepsilon (y^0 - 1) (1-e^{-t/\varepsilon}). $$ In particular, when $y^0 = 1$ (the initial speed according to the massless problem), one has $x_\varepsilon(t) = t$, so the massive solution is (in this particular case) exactly the same.

On the contrary, when $y^0 \neq 1$, there is a short period of time, say $t \in [0,5\varepsilon]$ where the solution $x_\varepsilon(t)$ does not resemble the massless solution $t$. The particle quickly transitions from its initial speed to the one predicted by the massless model. After this initial period called "boundary layer", one has $x_\varepsilon(t) \approx t + \varepsilon (y^0-1)$ so resembles the massless solution $t$, with a slight correction due to the initial phase.

b) Concerning expansions, let us consider the massless model $\dot{x} + x = 0$, with initial condition $x(0) = 1$, so that $x(t) = e^{-t}$ is the solution, with $\dot{x}(0) = -1$ the initial speed in this model. Now consider the massive version $\varepsilon \ddot{x}_\varepsilon + \dot{x}_\varepsilon + x_\varepsilon = 0$ with $x_\varepsilon(0) = 1$ and let us use $\dot{x}_\varepsilon(0) = -1$ as the initial condition matching the one of the massless problem, to avoid the boundary layers mentioned above.

You can then indeed look (at least formally) for a solution under the form $$ x_\varepsilon(t) = \sum_{k=0}^{+\infty} \varepsilon^k x_k(t) $$ with $\dot{x}_0 + x_0 = 0$, $x_0(0) = 1$ and, for $k \geq 0$, $$ \dot{x}_{k+1} + x_{k+1} = - \ddot{x}_k. $$ Here $x_0$ is the solution to the massless problem so $x_0(t) = e^{-t}$, and you can solve iteratively using the variation of constant formula. For example, $x_1(t) = - t e^{-t}$, and so on.

In this case, you can prove that if you truncate the series, you have a solution which approximates $x_\varepsilon$ at any precision.