\documentclass[a4paper,11pt]{article}
\parskip=1.5ex
\usepackage{verbatim,a4margins,amsmath,amssymb,bm,hyperref}

\newcommand{\dd}[2]{\frac{\mathrm d #1}{\mathrm d #2}}
\newcommand{\dds}[2]{\frac{\mathrm d^2 #1}{{\mathrm d #2}^2}}
\newcommand{\pdd}[2]{\frac{\partial #1}{\partial #2}}
\newcommand{\pdds}[2]{\frac{\partial^2 #1}{\partial{#2}^2}}
\newcommand{\pddd}[3]{\frac{\partial^2 #1}{\partial #2\partial #3}}

\newcommand{\rowfour}[4]{\left[\begin{array}{cccc}#1&#2&#3&#4\end{array}\right]}
\newcommand{\vectwo}[2]{\left[\begin{array}{c}#1\\#2\end{array}\right]}
\newcommand{\mxtwo}[4]{\left[\begin{array}{cc}#1&#2\\#3&#4\end{array}\right]}

\newcommand{\sq}[1]{\left(#1\right)^2}
\newcommand{\floor}[1]{\left\lfloor#1\right\rfloor}
\newcommand{\abs}[1]{\left|#1\right|}

\newcommand{\rem}[1]{}

\newcommand{\A}{\mathbf A}
\newcommand{\E}{\mathbf E}
\newcommand{\F}{\mathbf F}
\newcommand{\B}{\mathbf B}
\newcommand{\J}{\mathbf J}
\newcommand{\0}{\mathbf 0}

\newcommand{\av}{\mathbf a}
\newcommand{\f}{\mathbf f}
\newcommand{\uv}{\mathbf u}
\newcommand{\vv}{\mathbf v}
\newcommand{\w}{\mathbf w}
\newcommand{\x}{\mathbf x}
\newcommand{\y}{\mathbf y}
\newcommand{\z}{\mathbf z}

\newcommand{\p}{\mathbf p}
\newcommand{\s}{\mathbf s}
\newcommand{\rv}{\mathbf r}
\newcommand{\cv}{\mathbf c}

\newcommand{\e}{\mathbf e}
\newcommand{\m}{\mathbf m}
\newcommand{\M}{\mathbf M}

\newcommand{\de}{\,\mathrm d}

\newcommand{\ii}{\mathrm i}

\DeclareMathOperator{\cosmh}{cos\{h\}}
\DeclareMathOperator{\sinmh}{sin\{h\}}

\begin{document}

\title{Exact Timestep for a Pairwise Coulomb Collision}
\author{Stephen Brooks}
\date{February 24, 2025}
\maketitle

\section{Introduction}
Standard numerical integrators work well for many-body Coulomb repulsion problems when the timestep is much shorter than the timescale of relative position changes.  However, for `hard' collisions in which two particles have a near miss and exchange a lot of momentum within one timestep, they understandably struggle.  This note proposes using the exact solution of Keplerian two-body orbits (usually hyperbolic) to calculate the momentum exchange with other particles: either a selection of the `closest' ones or all of them.

\section{Pairwise Collision in Central Frame}
Take the masses of two particles to be $m_{1,2}$ and their charges to be $q_{1,2}$.  Adopt the frame where the total momentum $\p=m_1\vv_1+m_2\vv_2=\0$ and also the centre of mass $m_1\x_1+m_2\x_2=\0$ is at the origin.  From this, it is always true that $\x_2=-\frac{m_1}{m_2}\x_1$ and $\x_1-\x_2=\frac{m_1+m_2}{m_2}\x_1$.  Similar relations hold for $\vv_1$, meaning the problem is restricted to a two dimensional subspace spanned by $\x_1$ and $\vv_1$.  The radial Coulomb force on the first particle is
\[ \F_1(|\x_1|) = \frac{q_1q_2}{4\pi\epsilon_0}\frac{\x_1-\x_2}{|\x_1-\x_2|^3} = \frac{q_1q_2}{4\pi\epsilon_0}\frac{\frac{m_1+m_2}{m_2}\x_1}{(\frac{m_1+m_2}{m_2})^3|\x_1|^3}
 = \frac{q_1q_2}{4\pi\epsilon_0}\left(\frac{m_2}{m_1+m_2}\right)^2\frac{\x_1}{|\x_1|^3}. \]
Defining $r=|\x_1|$, this is a central force of magnitude
\[ F_1(r) = \frac{q_1q_2}{4\pi\epsilon_0}\left(\frac{m_2}{m_1+m_2}\right)^2\frac1{r^2}. \]

\section{Binet Equation}
The trajectory shape can be conveniently found from the Binet equation.  Defining $u=1/r$ and using polar coordinates, it states, for a central force $F$,
\[ F = -mh^2u^2\left(\dds{u}{\theta}+u\right), \]
where $h=r^2\dot\theta$ is a conserved angular-momentum-like quantity.  For this problem, $F=ku^2$ where $k=\frac{q_1q_2}{4\pi\epsilon_0}(\frac{m_2}{m_1+m_2})^2$.  Thus,
\[ ku^2 = -m_1h^2u^2\left(\dds{u}{\theta}+u\right)
\qquad \Rightarrow \qquad
\dds{u}{\theta}+u = -\frac{k}{m_1h^2} \]
\[ \Rightarrow \qquad u(\theta) = A\cos\theta + B\sin\theta -\frac{k}{m_1h^2} \]
for some constants $A$, $B$.  The origin for $\theta$ is free to choose at this point, so it simplifies the calculation to choose it such that $B=0$ and $A\ge0$.  The solution is now
\[ u(\theta) = A\cos\theta +C  \]
where $C=-\frac{k}{m_1h^2}$ is a constant that is negative for repelling particles with $q_1q_2>0$ and positive for attracting particles with $q_1q_2<0$.

\section{Initial Conditions}
Call the value of $\theta$ corresponding to the initial condition $\theta_0$, so that 
\[ u(\theta_0) = \frac1{|\x_1|} = A\cos\theta_0 + C, \]
where in this section $\x_1$ and $\vv_1$ are evaluated at the initial condition.

The first derivative satisfies
\[ u' = \dd{u}{\theta} = \dd{u}{r}\dd{r}{t}\dd{t}{\theta} = \frac{-1}{r^2}\dot r \frac1{\dot\theta}. \]
Initially we have $\dot r = \vv_1\cdot\x_1/|\x_1|$ and $\dot\theta = |\x_1\times\vv_1|/|\x_1|^2$, so 
\[ u'(\theta_0) = \frac{-1}{|\x_1|^2}\frac{\vv_1\cdot\x_1}{|\x_1|}\frac{|\x_1|^2}{|\x_1\times\vv_1|} = \frac{-\x_1\cdot\vv_1}{|\x_1||\x_1\times\vv_1|}. \]
These initial conditions also determine
\[ h = r^2\dot\theta = |\x_1|^2 \frac{|\x_1\times\vv_1|}{|\x_1|^2} = |\x_1\times\vv_1|, \]
which is a conserved quantity proportional to angular momentum.  Using the solution for $u(\theta)$ gives
\[ u'(\theta_0) = \frac{-\x_1\cdot\vv_1}{|\x_1|h} = -A\sin\theta_0. \]

Taking the sum of squares gives 
\[ (A\cos\theta_0)^2+(-A\sin\theta_0)^2 = A^2 = \left(\frac1{|\x_1|}-C\right)^2 + \left( \frac{-\x_1\cdot\vv_1}{|\x_1|h} \right)^2 \]
\[ \Rightarrow \qquad A = \frac1{|\x_1|}\sqrt{\left(1-|\x_1|C\right)^2 + \left( \frac{\x_1\cdot\vv_1}{h} \right)^2}. \]
Taking the ratio gives
\[ \frac{A\sin\theta_0}{A\cos\theta_0} = \tan\theta_0 = \frac{ \frac{\x_1\cdot\vv_1}{|\x_1|h} }{ \frac1{|\x_1|}-C }
 = \frac{ \x_1\cdot\vv_1 }{ (1-|\x_1|C)h }. \]
But also \[ \cos\theta_0 = \frac{ \frac1{|\x_1|}-C}{A} = \frac{1-|\x_1|C}{|\x_1|A} \] is useful.

\section{Time Dependence and the Eccentric Anomaly}
As $h=r^2\dot\theta$ is a constant found in the previous section, $\dot\theta=\dd{\theta}{t}=hu^2$.  The equation to solve is
\[ \dd{\theta}{t}=hu^2=h(A\cos\theta+C)^2. \]
Kepler's eccentric anomaly $E$ is defined via
\[ \cosmh E = \frac{e+\cos\theta}{1+e\cos\theta} \]
for some eccentricity $e$ that will be found later.  Terms in \{brackets\} are used for the hyperbolic form when $|e|>1$.  The sign of $\sinmh E$ is defined to be the same as $\sin\theta$.  The variable of the differential equation can be changed to $E$ via the chain rule:
\begin{eqnarray*}
\dd{E}{t} &=& \dd{E}{\cosmh E}\dd{\cosmh E}{t} \\
 &=& \frac{\{-\}1}{-\sinmh E}\left(\frac1{1+e\cos\theta}\dd{\cos\theta}{t}+(e+\cos\theta)\frac{-1}{(1+e\cos\theta)^2}e\dd{\cos\theta}{t}\right) \\
 &=& \frac{\{-\}1}{-\sinmh E}\left(\frac1{1+e\cos\theta}+\frac{-e(e+\cos\theta)}{(1+e\cos\theta)^2}\right)\dd{\cos\theta}{t} \\
 &=& \frac{\{-\}1}{-\sinmh E}\frac{1+e\cos\theta-e^2-e\cos\theta}{(1+e\cos\theta)^2}(-\sin\theta)h(A\cos\theta+C)^2 \\
 &=& \frac{\{-\}\sin\theta}{\sinmh E}(1-e^2)h\frac{(A\cos\theta+C)^2}{(1+e\cos\theta)^2} \\
 &=& \frac{\sin\theta}{\sinmh E}|1-e^2|h\frac{(A\cos\theta+C)^2}{(1+e\cos\theta)^2}.
\end{eqnarray*}
Next, $\sinmh E$ needs to be evaluated:
\[ \sinmh E = \pm\sqrt{\{-\}(1-\cosmh^2 E)} = \pm\sqrt{\{-\}\left(1-\frac{(e+\cos\theta)^2}{(1+e\cos\theta)^2}\right)} \]
\[ = \pm\sqrt{\{-\}\left(\frac{1+2e\cos\theta+e^2\cos^2\theta-e^2-2e\cos\theta-\cos^2\theta}{(1+e\cos\theta)^2}\right)} \]
\[ = \pm\frac{\sqrt{\{-\}(1-e^2-\cos^2\theta+e^2\cos^2\theta)}}{1+e\cos\theta}
= \pm\frac{\sqrt{|1-e^2|(1-\cos^2\theta)}}{1+e\cos\theta} = \frac{\sqrt{|1-e^2|}\sin\theta}{|1+e\cos\theta|}. \]
Substituting this back in,
\[ \dd{E}{t} = \frac{\sin\theta}{\sinmh E}|1-e^2|h\frac{(A\cos\theta+C)^2}{(1+e\cos\theta)^2}
 = \frac{|1+e\cos\theta|}{\sqrt{|1-e^2|}} |1-e^2|h\frac{(A\cos\theta+C)^2}{(1+e\cos\theta)^2} \]
\[ = \sqrt{|1-e^2|}h\frac{(A\cos\theta+C)^2}{|1+e\cos\theta|}
= \sqrt{|1-e^2|}hC^2\frac{(1+\frac{A}{C}\cos\theta)^2}{|1+e\cos\theta|}. \]
Setting $e=\frac{A}{C}$ gives
\[ \dd{E}{t} = \sqrt{|1-e^2|}hC^2|1+e\cos\theta|. \]
Note that this $e$ has the same sign as $C$, so is the usual definition of eccentricity for attracting particles, but negative for repelling particles.

To express $\dd{E}{t}$ back in terms of $E$, note that
\[ \cosmh E - \frac1e = \frac{e-\frac1e}{1+e\cos\theta} \]
\[ \Rightarrow \qquad 1+e\cos\theta = \frac{e-\frac1e}{\cosmh E - \frac1e} = \frac{e^2-1}{e\cosmh E -1}. \]
Therefore 
\[ \dd{E}{t} = \sqrt{|1-e^2|}hC^2\frac{|1-e^2|}{|1-e\cosmh E|} = \frac{|1-e^2|^{3/2}hC^2}{|1-e\cosmh E|}. \]
Kepler's equation suggests an implicit solution of the form \[ qt= E - e\sinmh E \] for some constant $q$.  Differentiating both sides with respect to $t$ gives 
\[ q = \dd{E}{t} - e\cosmh E\dd{E}{t} = (1-e\cosmh E)\dd{E}{t} \]
and substituting the $\dd{E}{t}$ found above gives
\[ q = (1-e\cosmh E)\frac{|1-e^2|^{3/2}hC^2}{|1-e\cosmh E|} = \sigma |1-e^2|^{3/2}hC^2, \]
which is indeed a constant.  The sign $\sigma$ of $1-e\cosmh E$ is
\begin{itemize}
\item $1$ for elliptic orbits ($0\le e<1$), 
\item $-1$ for hyperbolic attractive orbits ($e>1$, $C>0$) and 
\item $1$ for hyperbolic repulsive orbits ($e<-1$, $C<0$).
\end{itemize}

\section{Time Step}
With $A$ and $\theta_0$ known, calculate $e=\frac{A}{C}$ and proceed via
\[  \cosmh E_0 = \frac{e+\cos\theta_0}{1+e\cos\theta_0}, 
\qquad
 qt_0= E_0 - e\sinmh E_0 \]
to find the initial time $t_0$ in Kepler's equation.

Suppose the position is required at $t=t_0+\delta t$ for some time step $\delta t$.  A solution of $qt= E - e\sinmh E$ is needed.  This does not have a closed-form solution but Newton--Raphson iterations exist (see Appendices).  Then
\[ \cos\theta = \frac{e-\cosmh E}{e\cosmh E - 1}, \qquad
\sin\theta = \frac{|1+e\cos\theta|}{\sqrt{|1-e^2|}} \sinmh E \]
 and \[ u(\theta)=A\cos\theta+C. \]

A convenient way of handling the frames of reference is to rotate $\x_1$ by angle $\theta-\theta_0$ about the $\x_1\times\vv_1$ axis and then scale by $r/r_0=u(\theta_0)/u(\theta)$.

To find the velocity, it is easier to use the known differential equations $\dot\theta=hu^2$ and $u'(\theta)=-A\sin\theta$ to get
\[ \dd{r}{t} = \dd{r}{u}\dd{u}{\theta}\dd{\theta}{t} = \frac{-1}{u^2}(-A\sin\theta)hu^2 = Ah\sin\theta. \]

\end{document}