Chasing down minus sign in the expression $p_\mu - e A_\mu$ for the momentum In his Quantum Field Theory book, Folland says that the prescription for incorporating a 4-potential $A$ into the Dirac equation $(i \gamma^{\mu} \partial_\mu - m)\psi = 0$ is to replace $i\partial_\mu \to i\partial_\mu - eA_\mu$. I'm trying to understand where this comes from, in a principled way.
Let's use the $(1, -1, -1, -1)$ metric convention. The action for a particle of mass $m$ and charge $e$ in Minkowski space $M^4$ with the 4-potential $A$ used by Folland is  $$\int(-m d\tau - eA)$$
Writing in coordinates the path $\mathbf{x} = \mathbf{x}(t)$, this is the integral with respect to coordinate time of the Lagrangian
$$\mathcal{L} = -m \sqrt{1 - \mathbf{v}^2} - e A_j \dot{x}^j - eA_0$$
This gives us the canonical momenta
$$\frac{\partial \mathcal{L}}{\partial \dot{x}^j} = p_j - e A_j$$
where $$p_j = m \frac{ \dot{x}^j }{ \sqrt{1 - \mathbf{v}^2} }$$ is the relativistic momentum.
We also have the energy
\begin{split} H & = \frac{\partial \mathcal{L}}{\partial \dot{x}^j} \dot{x}^j - \mathcal{L}\\
& = \left( m \frac{ \dot{x}^j }{ \sqrt{1 - \mathbf{v}^2} } - eA_j \right) \dot{x}^j - (-m \sqrt{1 - \mathbf{v}^2} - e A_j \dot{x}^j - eA_0)  \\
&= m \left( \frac{\mathbf{v}^2}{\sqrt{1 - \mathbf{v}^2}} + \sqrt{1 - \mathbf{v}^2} \right) + eA_0 \\
&= \frac{m}{\sqrt{1 - \mathbf{v}^2}} + eA_0 \\
&= p_0 + eA_0 \end{split}
where $p_0$ is the 0-momentum i.e. relativistic energy.
I understand that to quantize energy-momentum we replace $p^\mu \to i\partial^\mu$, so $p_0 \to i \partial_0$ and $p_j \to -i \partial_j$. So shouldn't $p_j - eA_j$ be replaced by $-i\partial_j - eA_j$? And similarly shouldn't $p_0 + eA_0$ be replaced by $i\partial_0 + eA_0$?
I reckon that people might prefer to think of the relativistic momentum as a vector and not a covector. So maybe what I've been calling $p_j$ is really $-p_j$. This would correctly get us to $i\partial_j - eA_j$ for the spatial indices, despite seeming kind of arbitrary to me. But we still have the problem with the $0$ index, for which I can't see a way around. Moreover, I'm interested in some principled way to understand what's the correct sign, rather than fiddling around with (what seem to me like) ad hoc rules until I get a matching formula.
 A: The second-to-last paragraph in the question starts with
$$
 p^\mu\to i\partial^\mu. 
$$
This implies
$$
 p_\mu\to i\partial_\mu,
$$
because the index on both sides is lowered using the same metric (regardless of the sign convention for the metric). This equation can be separated into
\begin{align*}
 p_0 &\to i\partial_0
\\
 p_j &\to i\partial_j.
\end{align*}
The second-to-last paragraph writes $p_j\to -i\partial_j$ instead, which is incorrect (I mean, inconsistent with $p^\mu\to i\partial^\mu$). That's one sign error.
By the way, the replacement $p_\mu\to -i\partial_\mu$ is a more common convention. Both conventions are legitimate (they're just conventions), but when combined with the result derived below, the more common convention seems to be the one that agrees with Folland's convention $i\partial_\mu\to i\partial_\mu-eA_\mu$.

Independently of any relationship between the canonical momenta and the partial derivatives $\partial_\mu$, the question also asks about the inconsistency between the expressions $p_j-e A_j$ and $p_0+eA_0$. To resolve that, let's use a Lorentz-covariant approach.
First, the $A$ in the action is an abbreviation for the one-form
$$
 A = A_\mu dx^\mu = A_\mu \dot x^\mu\,d\tau,
$$
so the action is
$$
 S = -\int (m\,+e\dot x^\mu A_\mu)d\tau,
$$
where a dot means a derivative with respect to the proper time $\tau$. To compute variational derivatives in a covariant way, we can switch to an arbitrary worldline parameter $\lambda$:
$$
 S = -\int \left(m\,(\dot x_\mu\dot x^\mu)^{1/2} 
  +e\dot x^\mu A_\mu\right)d\lambda,
$$
where now a dot means a derivative with respect to $\lambda$. The lagrangian is
$$
 L = -\left(m\,(\dot x_\mu\dot x^\mu)^{1/2} 
 + e\dot x^\mu A_\mu\right),
$$
so the conjugate momenta are
$$
 p_\mu\equiv \frac{\delta L}{\delta \dot x^\mu}
 = -\frac{m}{(\dot x_\mu\dot x^\mu)^{1/2}}\dot x_\mu
  - eA_\mu.
$$
Now that we're done computing the variational derivatives, we can impose the constraint $\dot x_\mu\dot x^\mu=1$ so that the parameter $\lambda$ becomes the proper time $\tau$. This gives
$$
 p_\mu = -m\dot x_\mu - e A_\mu,
$$
which implies
\begin{align*}
 p_0 &= -m\dot x_0 - e A_0 \\
 p_j &= -m\dot x_j - e A_j.
\end{align*}
The sign inconsistency is gone, of course, because it never had any opportunity to arise: the whole calculation was Lorentz covariant.

Now that we have this result, we can look back at the approach used in the question and try to understand what might have gone wrong. Two steps in that approach are suspect:

*

*In going from the action to the lagrangian, the question tacitly set the worldline parameter equal to the time coordinate $t\equiv x^0$. The action is indeed reparameterization invariant, but setting the integration parameter equal to one of the dynamic variables $x^0$ is not an obviously-legitimate thing to do. If we want to do it, then we need to think carefully about what it really means and how to do it right.


*Even if we gloss over issue #1, how is the expression for the Hamitonian $H$ shown in the question related to the canonical momentum $p_0$ that we defined above? We can't just assume that the $0$th component of the canonical momentum defined as above is the same thing as the Hamiltonian (or that $H=p_0+eA_0$). The Hamiltonian may turn out to be related to the canonical momentum in some simple way in this case, but it's clearly not true in general: the canonical momenta are specific to the particle/field $F$ in the denominator of $\delta L/\delta \dot F$, but the Hamiltonian is the generator of time-translations for the whole system, which in general may include many particles/fields.
I'm betting that both of these interesting conceptual issues are already addressed elsewhere on Physics Stack Exchange. If I can find them, then I'll come back later and add the relevant links.
A: I remember that I have also struggled with the sign of the minimal coupling term several times in the past. Therefore, I feel compassionate with your attempts to derive it.
Would you accept the non-relativistic Schrödinger equation for the bare electric field of the nucleus (i.e. without magnetic potential) as

"some principled way to understand what's the correct sign"?

If the Schrödinger equation got the sign wrong, a bound state would turn into a free state, so there cannot be any ambiguity.
$$i\hbar\frac{\partial}{\partial t}\psi=\hat H\psi = (\hat T+V)\psi=\frac{1}{2m}\left(-i\hbar\frac{\partial}{\partial \vec x}\right)^2\psi +q\Phi\psi$$
Note that if $\Phi>0$ represents the nucleus field and $q<0$ the electron, we have $q\Phi\sim -1/|\vec r|$.
Writing the Schrödinger equation by 4-components yields
$$i\hbar\frac{\partial}{\partial x^0}\psi=\hat T\psi+qA_0\psi$$
Note, that the electric potential is the zero component of the electromagnetic 4-potential $A^\mu$, which in your metric convention is not different between covariant and contravariant, i.e. $A_0=A^0$, so no danger to make a sign mistake here.
Moving the potential term to the left-hand side we get
$$(i\hbar\partial_0-qA_0)\psi=\hat T\psi$$
No matter how fancy it is to justify the relativistic generalization of the Schrödinger equation (Dirac, Klein-Gordon), especially the kinetic term on the right hand side, the only way to claim equal rights for time and space is to generalize the "minimal coupling substitution" $i\hbar\partial_0\to i\hbar\partial_0-qA_0$ for space to
$$i\hbar\partial_k\to i\hbar\partial_k-qA_k$$
Be aware that the particular charge $q$ might include a sign, while $e$ might also be considered strictly positive (elementary charge constant), depending on context/author.
A: Let $C$ be a field such that one had a covariant derivative $\partial_\mu+C_\mu$ satisfying:
\begin{align*}
&(D_\mu \Psi)'=e^{i\alpha(x)}D_\mu \Psi
\\
&\Longrightarrow (\partial_\mu+C'_\mu)[e^{i\alpha(x)}\Psi]=e^{i\alpha(x)}(\partial_\mu+C_\mu)\Psi
\\
&\Longrightarrow C'_\mu\Psi=C_\mu\Psi-i\partial_\mu \alpha(x)\Psi
\end{align*}
So one can see that our field $C$ is gauge transformed exactly as $ieA$. Then our covariant derivative is $D_\mu=\partial_\mu+ieA_\mu$. So the Dirac equation becomes:
\begin{equation*}
(i\gamma^\mu D_\mu-m)\Psi=0 \Longrightarrow (\gamma^\mu (i\partial_\mu-e A_\mu)-m)\Psi=0
\end{equation*}
