The standard Lagrangian $L$ written in local coordinates for a free, relativistic particle of mass $m > 0$ is given by $$L(q, \dot{q}) = -m \sqrt{-g_{\mu\nu}(q) \dot{q}^\mu \dot{q}^\nu}\tag{1}$$ where I have chosen to use the mostly pluses metric signature. Parameterizing the particle's world line by its proper time $\tau$, the corresponding action $S$ is $$S[q] = \int_{\tau_1}^{\tau_2} -m \sqrt{-g_{\mu\nu}(q) \frac{dq^\mu}{d\tau} \frac{dq^\nu}{d\tau}}d\tau.\tag{2}$$ We might consider extending this situation by introducing a (possibly 4-velocity-dependent) external potential $U$ so that $$S[q] = \int_{\tau_1}^{\tau_2} \bigg[-m \sqrt{-g_{\mu\nu}(q) \frac{dq^\mu}{d\tau} \frac{dq^\nu}{d\tau}} - U\Big(q, \frac{dq}{d\tau}\Big)\bigg]d\tau.\tag{3}$$ This is all nice and Lorentz-covariant so far. But what happens when we try to introduce multiple particles with masses $m_{(i)} > 0$ (Parentheses are to distinguish labels from indices) and their interaction potentials $$V_{(i, j)}(q_{(i)}, q_{(j)}, \dot{q}_{(i)}, \dot{q}_{(j)})~?\tag{4}$$ Let us assume that the form of the interactions is compatible with special relativity. Perhaps they are mediated by some field whose dynamics are specified by a separate Lagrangian density. What I would like to know is: What is the Lorentz-covariant Lagrangian for this new situation with multiple particles?
The first issue I encounter when considering this problem is what parameter to use in the action. To keep things Lorentz-covariant, we might try something like the proper time. But then which particle's proper time do we use? To avoid singling out any one particle, we might choose to parameterize each world line by its coordinate time. But then this artificially singles out a preferred reference frame.
