Does the Darwin Lagrangian neglect deviations from the Coulomb field? The Darwin Lagrangian is said to describe the interaction between two charges to order $(v/c)^2$, and consists of a free part
$$L_0 = \sum_{i = 1, 2} \frac12 m_i v_i^2 + \frac{1}{8c^2} m_i v_i^4$$
and an interaction
$$L_{\text{int}} = - \frac{q_1 q_2}{r} + \frac{q_1 q_2}{r} \frac{1}{2 c^2} (\mathbf{v}_1 \cdot \mathbf{v}_2 + (\mathbf{v}_1 \cdot \hat{\mathbf{r}})(\mathbf{v}_2 \cdot \hat{\mathbf{r}})).$$
However, it appears to me that this interaction doesn't capture all of the $O(v^2)$ corrections. The first term is the $O(v^0)$ Coulomb potential, while the second term presumably captures the Lorentz force on one charge due to the magnetic field produced by the other; each of these are $O(v)$, giving an $O(v^2)$ effect.
However, another effect is that the electric field produced by a moving charge with uniform velocity differs from the Coulomb field at $O(v^2)$. This seems to imply that $L_{\text{int}}$ should include terms proportional to $v_1^2$ and $v_2^2$, but it doesn't.
I can think of a couple things that could be going on here.


*

*The Darwin Lagrangian just explicitly excludes this interaction. But then what do people mean when they say it's accurate to order $(v/c)^2$?

*Even though it looks like the required term isn't present in $L_{\text{int}}$, the desired contribution appears after one computes the complicated Euler-Lagrange equations. But this didn't seem to occur when I did it. The accelerations of the charges are quite complicated, but they don't contain any terms proportional to $q_1 q_2$ besides the magnetic term.

*This interaction cannot be included without accounting for retardation or radiation effects, which the Darwin Lagrangian explicitly excludes, since otherwise one would need to keep track of the field configuration. But the interaction is present even for charges which have had uniform velocity forever.

*The magnetic interaction somehow also accounts for the effect. I don't see how this could happen, because the magnetic force on charge $i$ vanishes when $v_i = 0$, and this effect doesn't. 


What's going on here?
 A: 
… what do people mean when they say it's accurate to order $(v/c)^2$?

Lagrangian of relativistic particle carrying a charge $q_a$ in external EM field is given by:
$$
L_a = -m_a c^2  \sqrt{1-\frac{v_a^2}{c^2}} - q_a \varphi + \frac{q_a}{c} \mathbf{A}\cdot \mathbf{v}_a.
$$
As usual, we assume that the fields $\varphi$ and $\mathbf{A}$ are created by the rest of the charges in the system. We are formally expanding this Lagrangian in powers of $1/c$ (the actual small dimensionless quantity of the perturbation theory would be $v/c$). Radiative effects appear at order $c^{-3}$, so we could expect that if we keep only terms up to $1/c^2$ the resulting approximate Lagrangian would include relativistic corrections but would still correspond to a conservative system. 
This means that in the interaction Lagrangian it is enough to calculate the scalar potential $\varphi$ up to $1/c^2$ terms, while vector potential  $\mathbf{A}$ only up to $1/c$ (since corresponding term in the Lagrangian already has $1/c$).

… another effect is that the electric field produced by a moving charge with uniform velocity differs from the Coulomb field at $O(v^2)$. This seems to imply that $L_{\text{int}}$ should include terms proportional to $v_1^2$ and $v_2^2$, but it doesn't.

Not necessarily. Note that the Lagrangian does not has “electric field” but only scalar and vector potentials. By imposing the Coulomb gauge condition ($\nabla \cdot \mathbf{A}=0$) the equation for scalar potential reads:
$$
\Delta \varphi = - 4 \pi \rho.
$$
This means that the scalar potential in this gauge for the system of moving charges is independent of velocity and do not contribute any $1/c^2$ terms to $L_\text{int}$:
$$
\varphi (\mathbf{r}) = \int \frac{\rho(\mathbf{r}') }{|\mathbf{r}'-\mathbf{r}|} \,d^3\mathbf{r}'
$$
Vector potential in Coulomb gauge satisfies the equation:
$$
\Delta \mathbf{A} - \frac{1}{c^2} \frac{\partial^2 \mathbf{A} }{\partial t^2} = -\frac{4 \pi}{c}\left(\mathbf{j}-\frac{1}{4 \pi}\nabla \frac{\partial \varphi}{\partial t}\right).
$$
Second time derivative on the left is proportional to $1/c^2$ so in our approximation it could be dropped: the vector potential could be taken as quasi-static. We also make use of charge conservation law:$$\frac{\partial \rho }{\partial t} + \nabla \cdot \mathbf{j} =0.$$
This allows us to write solution for the vector potential:
$$
\mathbf{A}(\mathbf{r})=\frac{1}{2c}\int \frac{\left[\mathbf{j}(\mathbf{r}')+\hat{\mathbf{r}}(\hat{\mathbf{r}}\cdot \mathbf{j}(\mathbf{r}'))\right]}{|\mathbf{r}'-\mathbf{r}|}d^3\mathbf{r}'.
$$
As expected, this expression is proportional to $1/c$.
If we assume that $\rho$ and $\mathbf{j}$ correspond to charge and current density  of all charges except $q_a$ the potentials would now read:
$$
 \varphi = \sum_b' \frac{q_b}{r_{ab}}, \qquad \mathbf{A} =\frac{1}{2c}\sum_b' 
\frac{q_b[\mathbf{v}_b+\hat{\mathbf{r}}(\hat{\mathbf{r}}\cdot \mathbf{v}_b)]}{r_{ab}}.
$$
Inserting these expressions back into interaction Lagrangians would indeed produce the Darwin Lagrangian.
The takeaway here, is that this is wrong to call terms of $L_\text{int}$ magnetic and electric. Instead, they correspond to scalar and vector potential. Of course, for moving charge the vector potential does contribute to the perturbation of electric field due to the motion of the charge.
