Modifying Feynman-Wheeler absorber theory to work with arbitrary potentials?

I’m trying to consider relativistic multi-body dynamics in special relativity. In classical mechanics, it’s easy to write a simple $n$-body system with arbitrary potential $V$:

begin{equation}
m ddot{x}_ i=sum_ j – nabla_ {x_ i} V(|x_ i-x_ j|). tag{1}
label{1}
end{equation}
In special relativity, it’s tempting to just replace this with the retarded potential, where $x_ j$ is evaluated at the time where $c |Delta t|=|x_ i-x_ j|$. However this ends up in solutions blowing up over time. I want to find an action for a 2-body system which reduces to equation ref{1} in the limit $vll c$, but which also has correct and physically meaningful conservation laws.

Since this is all within the realm of radiation reaction, I figure a surefire starting point is to consider things from a Lagrangian Feynman-Wheeler type system (Classical Electrodynamics in Terms of Direct Interparticle Action), since its symmetries will pretty directly give conservation laws (albeit with some speed of light delays). I label the two particles $a$ and $b$, and I’m working with $c=1$, unit charges and masses, signature $(- + + +)$, and $t$ an arbitrary parameter labeling the world lines. Then the action is:

$$A=-sum_{i=a,b}int dt sqrt{-dot x_i^mu dot x_{imu}} – iint delta((x_a-x_b)^2) dot x_a^mu dot x_{bmu}dt_1 dt_2 label{2}tag{2}$$

Note that $dt sqrt{-dot x_i^mu dot x_{imu}}$ should really be considered as $sqrt{-dx_i^mu dx_{imu}}$, and that the double integral should really be considered as $dx_a^mu dx_{bmu}$. So we really are reparameterization invariant, and we really are integrating with respect to the world lines. (Also note: "$x^2$" in the delta function means $x^mu x_mu$.)

It’s easy to see that this gives the Coulomb force: Fix particle $b$ to the origin so that $x_b^mu(t)=(t,vec{0})$. Then for $x_a^mu(t)=(t,vec{x}_a(t))$, we find $dot x_a^mu dot x_{bmu}=1$. Apply the delta function identity $delta(g(x))=sum_{g(x_0)=0} delta(x-x_0)/|g'(x_0)|$ and integrate with respect to $t_2$ to get

$$iint delta((x_a-x_b)^2) dot x_a^mu dot x_{bmu}dt_1 dt_2 =int dt_1 sum_{t_2=t_a,t_r}frac{1}{|2(x_a^mu-x_b^mu) dot x_{bmu}|}=int dt_1 sum_{t_2=t_a,t_r}frac{1}{|2Delta t|}.label{3}tag{3}$$

$t_a$ and $t_r$ are the advanced and retarded times with $|Delta t|=|Delta x|$, so summing over the two we get the action of a single particle in a Coulomb potential $$int dt_1 frac{1}{|Delta x|}$$

So the term $|(x_a^mu-x_b^mu) dot x_{bmu}|$ turned into a vector difference $|Delta vec{x}|$. This leads to the idea: just multiply the interaction term by terms like that. The corrected action term might look something like this:

$$iint F(|(x_a^mu-x_b^mu) dot x_{bmu} /sqrt{- dot x_b^nudot x_{bnu}}|) delta((x_a-x_b)^2) dot x_a^mu dot x_{bmu}dt_1 dt_2. label{4}tag{4}$$

If $F(x)=xV(x)$ and particle $b$ is fixed at the origin, this gives the correct limit, and is Lorentz covariant and reparameterization invariant (that’s what the $sqrt{-ldots}$ term is for), but it also favors $x_a$ over $x_b$! Symmetrizing with respect to $a$ and $b$ also seems OK, because for $|frac{d}{dt} vec{x}_a| ll 1$ we should have $dot x_{amu} /sqrt{- dot x_a^nudot x_{anu}}approx (1,vec{0})$, but it feels like there should be a more simple route to go down.

Does anyone know of a way to do this, or have any better ideas on how to modify the interaction term?

Lorentz covariance and reparameterization invariance put some heavy restrictions on the action, so maybe it’s not possible to get a very elegant action with the desired properties.