How to separate electric and magnetic effects
The QED lagrangian has the form
$$
\newcommand{\opsi}{\overline\psi}
L = -\frac{1}{4e^2}F^{ab}F_{ab}
+ \opsi\gamma^a A_a\psi
+ L_\psi
\tag{1}
$$
where $L_\psi$ is the kinetic term for the matter field $\psi$. The coupling constant $e$ has been moved from the interaction term to the gauge-field kinetic term by rescaling the gauge field. The gauge-field kinetic term has the form
$$
-\frac{1}{4e^2}F^{ab}F_{ab}
=\frac{E^2-B^2}{2e^2}.
\tag{2}
$$
To separate electric and magnetic effects from each other, we can temporarily generalize (2) to
$$
\frac{E^2}{2e_E^2}-\frac{B^2}{2e_B^2}
\tag{3}
$$
where $e_E^2$ and $e_B^2$ are independent parameters. This ruins Lorentz symmetry, but only temporarily. After finishing the calculation and seeing how the electric/magnetic distinction manifests itself in the result, we can restore Lorentz symmetry by setting $e_E^2=e_B^2=e^2$.
The scattering amplitude
Since we've temporarily compromised Lorentz symmetry anyway, we might as well use the temporal gauge $A_0=0$. Then the action is
\begin{align}
L &= \frac{1}{2}\sum_{jk}\int_{\omega,p} A_j
\left(\frac{\omega^2\delta_{jk}}{e_E^2}
-\frac{p^2\delta_{jk}-p_j p_k}{e_B^2}
\right)A_k
\\
&+ \int_x \big(\opsi\gamma^k A_k\psi
+ L_\psi\big).
\tag{4}
\end{align}
with $j,k\in\{1,2,3\}$. The four-momentum is $(\omega,p_1,p_2,p_3)$. From (4), we can read off the form of the propagator for the gauge field:
\begin{align}
D_{jk}(\omega,p)
&=
\left(\frac{\omega^2}{e_E^2}-\frac{p^2}{e_B^2}\right)^{-1}
\left(\delta_{ij} - \frac{p_j p_k/e_B^2}{\omega^2/e_E^2}\right)
\\
&=
e_E^2\frac{e_B^2\delta_{jk}-e_E^2 p_j p_k/\omega^2}{
e_B^2\omega^2-e_E^2p^2}.
\tag{5}
\end{align}
The scattering amplitude is given by the usual expression, but with this propagator in place of the usual one.
Interpretation in the pure electric case ($e_B=0$)
When $e_B=0$ in the lagrangian, a nonzero magnetic field is infinitely costly, so the magnetic field is effectively constrained to be zero. Therefore, when $e_B=0$, we would expect (5) to describe the pure Coulomb interaction. This is not obvious the way (5) is written, because setting $e_B=0$ gives
$$
D_{jk}(\omega,p)=
e_E^2\frac{ p_j p_k}{p^2\omega^2}.
\tag{6}
$$
Despite appearances, this is consistent with a pure Coulomb interaction. To see why, recall that the interaction term in the original lagrangian is $J^a A_a$, where $J^a=\opsi\gamma^a\psi$ is the current. Gauge invariance implies current conservation $\partial_a J^a=0$, whose momentum-space version is $\omega J^0+p_k J^k=0$. Use this with (6) to get
$$
J^j D_{jk} J^k = e_E^2 \frac{J^0 J^0}{p^2}
\tag{7}
$$
in momentum space. This shows that the right-hand side of (7) is the expected pure Coulomb interaction. Throughout this answer, the current $J$ is an operator. The scattering amplitude (with $e_B=0$) is obtained by sandwiching (7) between the initial and final two-particle states.
Interpretation in the pure magnetic case ($e_E=0$)
When $e_E=0$ in the lagrangian, a nonzero electric field is infinitely costly, so the electric field is effectively constrained to be zero. Therefore, when $e_E=0$, we might have expected (5) to describe a purely magnetic interaction... but setting $e_E=0$ in (5) makes the propagator zero, which means no interaction at all! What's going on here?
In hindsight, setting $e_E=0$ in (5) gives exactly what we should have expected. In the temporal gauge, constraining the electric field to be zero is the same as constraining the gauge field to be independent of time. With this constraint, the current is unable to influence the gauge field, so the gauge field is unable to mediate any interaction between currents.
Another perspective: The electric and magnetic fields are not entirely independent of each other, because they're both written in terms of the gauge field. Constraining $B=0$ only constrains the transverse components of the gauge field, but constraining $E=0$ constrains the time-dependence of all components of the gauge field (in the temporal gauge).
How the amplitude reproduces magnetic forces
From now on, let's focus on the part of (5) that is not included in the purely electric part (6). I'll avoid calling this the purely magnetic part, for the reason given in the preceding paragraphs, but we can at least say that it includes the magnetic effects that the question is asking about. The difference between (5) and (6) is
\begin{align}
\delta D_{jk}(\omega,p)
&=
e_E^2\frac{e_B^2\delta_{jk}-e_E^2 p_j p_k/\omega^2}{
e_B^2\omega^2-e_E^2p^2}
-
e_E^2\frac{ p_j p_k}{p^2\omega^2}.
\tag{8}
\end{align}
Now that we've isolated the part of interest, we can set $e_E^2=e_B^2=e^2$ to get
\begin{align}
\delta D_{jk}(\omega,p)
&=
e^2\frac{\delta_{jk}-p_j p_k/\omega^2}{
\omega^2-p^2}
-
e^2\frac{ p_j p_k}{p^2\omega^2}
\\
&=
e^2\frac{\delta_{jk}-p_j p_k/p^2}{\omega^2-p^2}.
\tag{9}
\end{align}
The corresponding contribution to the electron-electron interaction is
$$
J^j \delta D_{jk} J^k.
\tag{10}
$$
The $p$ in (9) is the momentum transferred from one current to the other, and the currents $J$ here are operators. The scattering amplitude is obtained by sandwiching (10) between the initial and final two-particle states.
Now consider the Fourier transform of (10) in the center-of-mass frame, where the energy-transfer is $\omega = 0$. The result (shown in this Wikipedia page) is proportional to
$$
e^2 J^j(x)\frac{\delta_{jk}+r_j r_k/r^2}{r}J^k(y)
\tag{11}
$$
with $\vec r=\vec x-\vec y$. Insofar as the (quantum) particle has any well-defined location and velocity, the current $J(x)$ is proportional to the velocity and $x$ is the location, so (11) reproduces the Darwin lagrangian that was mentioned in the question.
The last paragraph was a bit loose with its interpretation of the currents $J$. Again, the currents in (10) are operators. By treating the currents as though they were classical currents, the last paragraph glossed over some important features of the scattering amplitude, including information about the particles' spins.
