Peierls Substitution with Time-Dependent Vector Potential My question is whether Peierls substitution really holds true for time-dependent electromagnetic (EM) potentials and, if yes, why.
To implement an electromagnetic field in a condensed matter system described by a Bloch Hamiltonian, I have often seen people make the following substitution of the hopping integrals
$t_{12}\rightarrow t_{12}' = e^{i\frac{q}{\hbar}\int_{\vec{R}_1}^{\vec{R}_1}\vec{A}(\vec{r})d\vec{r}} t_{12} \tag{1}\label{subst}$
with the EM vector potential $\vec{A}$ and connected lattice sites $\vec{R}_i$. This is correct as a result of the Peierls substitution (as long as closed integrals over $\vec{A}$ can be assumed $0$). The proof on the Wikipedia page is very nice, actually.
My problem is that I have also seen substitution $\ref{subst}$ for time-dependent EM potentials $\vec{A}(\vec{r},\tau)$ with time $\tau$. If one checks the proof on Wikipedia (every other proof I have seen is conceptually equivalent), one can easily see that for a time-independent EM potential, the Schrödinger equation is still satisfied: Assume
$
H(\vec{r})=\frac{(\vec{p})^2}{2m}+U(\vec{r})
\qquad\text{and}\qquad
t_{12}=-\int\phi_{\vec{R}_1}(\vec{r})^\dagger H\phi_{\vec{R}_2}(\vec{r})
\tag{2}\label{tDef}
$
with Wannier functions $\phi_{\vec{R}_i}$. If we now have eigenfunctions
$
\Psi_k(\vec{r})=\frac{1}{\sqrt{N}}\sum_{\vec{R}_i} e^{i\vec{k}\vec{R}_i}\phi_{\vec{R}_i}(\vec{r})
\qquad\text{with}\qquad
H\Psi_k(\vec{r})=E(\vec{k})\Psi_k(\vec{r}),
\tag{3}\label{energy}
$
we can use them to obtain
$
i\hbar \frac{d}{d\tau}
\left(
e^{-i\frac{E(\vec{k})}{\hbar}\tau}
\Psi_k(\vec{r})
\right)
=
E(\vec{k})
\left(
e^{-i\frac{E(\vec{k})}{\hbar}\tau}
\Psi_k(\vec{r})
\right)
=
H
\left(
e^{-i\frac{E(\vec{k})}{\hbar}\tau}
\Psi_k(\vec{r})
\right)
\tag{4}\label{satisfies}
$
and to thus satisfy the time-dependent Schrödinger equation.
Now we introduce an EM potential, such that
$
\tilde{H}=\frac{(\vec{p}-q\vec{A}(\vec{r}))^2}{2m}+U(\vec{r}).
\tag{5}\label{newHami}
$
Upon Peierls substitution
$
\phi_{\vec{R}_i}(\vec{r})\rightarrow
\tilde{\phi}_{\vec{R}_i}(\vec{r})=
e^{i\frac{q}{\hbar}\int\vec{A}(\vec{r})d\vec{r}}
\phi_{\vec{R}_i}(\vec{r}),
\tag{6}\label{subst2}
$
equation $\ref{satisfies}$ is still satisfied as
$
i\hbar \frac{d}{d\tau}
\left(
e^{-i\frac{E(\vec{k})}{\hbar}\tau}
\tilde{\Psi}_k(\vec{r})
\right)
=
\tilde{H}
\left(
e^{-i\frac{E(\vec{k})}{\hbar}\tau}
\tilde{\Psi}_k(\vec{r})
\right)
\tag{7}\label{satisfies2}
$
with $\tilde{\Psi}(\vec{r})$ like in definition $\ref{energy}$, but with all $\phi_{\vec{R}_i(\vec{r})}$ replaced with $\tilde{\phi}_{\vec{R}_i}(\vec{r})$ and with the original $E(\vec{k})$ (easy to check).
If one now wants to express this in the basis of the $\tilde{\phi}_{\vec{R}_i}(\vec{r})$ (like e.g. for a numerical evaluation where the hopping integrals are known), one indeed makes the substitution $\ref{subst}$. However, equation $\ref{satisfies2}$ remains true only if
$
\frac{d}{d\tau}\int\vec{A}(\vec{r},\tau)d\vec{r}=0
\tag{8}\label{question}
$
for all the integrals inside the $\tilde{\phi}_{\vec{R}_i}(\vec{r})$. Otherwise, on the left hand side of equation $\ref{satisfies2}$, countless obstructive prefactors arise, but they do not appear on the right hand side.
I personally want to use Floquet formalism (not immediately relevant for this discussion) to desribe a time-dependent system. However, I wonder whether I am truly allowed to use Peierls substitution. So far, I am using another gauge instead, such that $\vec{A}(\vec{r},\tau)$ vanishes and I am left with a scalar potential $\Phi(\vec{r},\tau)$. This, however, leads to terms which are not diagonal in $\vec{k}$ which is not nice either.
 A: Ok since I have to give a full, valid answer to the question, that can stand on its own:
The Peierls substitution does NOT work for time-dependent Vector potentials. The proof at Wikipedia is wrong, as it stands.
Here is my proof:
We start from the field operator representation of the Hamiltonian, for simplicity we assume that the Vector potential is spatial homogeneous (long wave length approximation)
\begin{align}
\label{Eq.ONE}
H(t) = \int \mathrm{d} r \Psi^\dagger(r) \left[ \frac{(-i \hbar \nabla - e A (t) )^2}{2m} + V(r) \right] \Psi(r)
\end{align}
We express the Field operators in time-independent atomic basis which are site diagonal (i.e. Wannier orbitals)
\begin{align}
\Psi(r) = \sum\limits_l \phi(r-R_l) c_l
\end{align}
Substituting yields
\begin{align}
H(t) = \sum\limits_{l l'} J_{l l'} (t) c_l^\dagger c_{l'}
\end{align}
with
\begin{align}
\label{Eq.Hop}
J_{l l'} (t) = \int \mathrm{d} r \phi^*(r-R_l) \left[ \frac{(-i \hbar \nabla - e A (t) )^2}{2m} + V(r) \right] \phi(r-R_{l'})
\end{align}
Note that the operators $c_l$ are time independent, the time dependence is only in the hopping integral.
We introduce a new set of Bloch basis states,
\begin{align}
\tilde{\psi}_k = \frac{1}{\sqrt{N}} \sum\limits_{R_l} e^{i\int_{R_l}^{r} e A \mathrm{d}r} e^{ikR_l} \phi_{R_l}(r)
\end{align}
and the corresponding inverse transformation
\begin{align}
\phi_{R_l}(r) = \frac{1}{\sqrt{N}} \sum\limits_k e^{-ik{R_l}} e^{-i\int_{R_l}^r e A \mathrm{d}r}  \tilde{\psi}_k(r) 
\end{align}
where $\phi_{R_l}(r) = \phi(r-{R_l})$.
Note that $\tilde{\psi}_k$ are instantaneous Eigenstates of the operator
\begin{align}
H(t) = \left[ \frac{(-i \hbar \nabla - e A (t) )^2}{2m} + V(r) \right],
\end{align}
they are time-dependent, however they only enter in the calculation of the hopping matrix elements. They are not used to solve the physical time-dependent Schroedinger equation.
The eigenequation reads
\begin{align}
H(t) \tilde{\psi}_k = \tilde{E}_k \tilde{\psi}_k
\end{align}
We use the representation of the Wanier orbitals in the diagonal Bloch states
to rewrite
\begin{align}
J_{l l'} &= \int \mathrm{d}r \sum\limits_{k,k'} \frac{1}{N} e^{ikR_l} e^{i\int_{R_l}^r eA \mathrm{d}r'} \tilde{\psi}^*_k(r) H(t) \tilde{\psi}_k(r) e^{-ik'R_{l'}} e^{-i \int_{R_l'}^r eA \mathrm{d}r'} \\
&= e^{ i \int_{R_l}^{R_{l'}} eA \mathrm{d}r'} \sum\limits_{k,k'}  e^{-ik'R_{l'}}  e^{ikR_{l}}  \int \mathrm{d}r\tilde{\psi}_k^*(r) H(t) \tilde{\psi}_{k'}(r) \\
&= e^{ i \int_{R_l}^{R_{l'}} eA \mathrm{d}r'} \sum\limits_{k} e^{ik(R_l-R_{l'})} \tilde{E}_k \\
\end{align}
Here comes the main problem. Only if $\tilde{E}_k = E_k$ holds we can deduce that the Peierls substitution works (i.e. its just a phase factor). But that is in general not true.
If $A$ is not time dependent, we can introduce a new set of Wannier Orbitals and indeed proof that the peierls substition works, as you do in your derivation. But if $A$ is time dependent, we should not introduce time depentent Wannier Orbitals, as that changes the form of the time dependent Schrödinger equation.
A: It holds but the Wikipedia, wannier orbitals and hamiltonian formalism unnecessarily blurs this very obvious fact.
First assume that we have an action of Dirac fermions(this can also be non relativistic fermions as in your case it won't effect the argument) in $2+1$ flat space-time.
$$S_F=\int  d^3x\bar{\psi}\gamma^{\mu}i(\partial_\mu)\psi+m_0\bar{\psi}\psi$$
now this is invariant under global $U(1)$ symmetry that is if I let
$$\psi(x)\to\psi(x)e^{i\alpha}$$ for some alpha and substitute it to above action, the action would not change. However if $\alpha$ depends on $x$ than the above action would change. So we write a new action which would not change under such transformations with $x$ dependent $\alpha$ that is
$$S_F=\int  d^3x\bar{\psi}\gamma^{\mu}i(\partial_\mu+iA_\mu)\psi+m_0\bar{\psi}\psi$$
such that we defined a new quantity $A$ which transforms as
$$A_\mu(x)\to A_\mu(x)+\partial_\mu\alpha(x)$$
so this is how we actually add an electromagnetic coupling to any system, just make the system invariant under $x$ dependent gauge transformation, it automatically adds the gauge coupling.
Ok, now the next step is just doing this in the lattice. First we wick rotate, than we enforce the local gauge invariance and then, we discretize, only the  spatial coordinates, and leave the temporal coordinate continuous then  the action is now,
\begin{equation}
S_{WF}[A]=\int d\tau \bigg[\sum_{\vec{r}}\bar{\psi}(r)(\gamma^{0}_E(\partial_0+iA^{E}_0)\psi(r)+(m_0+2R)\sum_{\vec{r}}\bar{\psi}(r)\psi(r)-\frac{1}{2}\sum_{\vec{r}i}[\bar{\psi}(r)(R-\gamma^i_E)\psi(r+\hat{i})U_{ir}+\bar{\psi}(r+\hat{i})(R+\gamma^i_E)\psi(r)U^{\dagger}_{ir} ]\bigg] 
\end{equation}
where $r=(\tau,\vec{r})$, $R$ is the wilson parameter which is dropped in continuum limit and as long as it is finite, the particular value is not important and $U_{ir}=\exp(iA_{ir})$
Here the $U$ is actually what you called  peierls substitution, but the actual name is Schwinger integral. Somehow they usually teach this very fundamental stuff in a very obscure weird and unecesarrly useless peirls point of view.
Edit To be more clear:
Let us assume we have a three dimensional Hamiltonian  in lattice that is,
$$H =a^2\bigg[(m_0a+3R)\sum_n\bar{\psi}_n\psi_n -\frac{1}{2}\sum_{n\mu}[\bar{\psi}_n(R-\gamma_\mu)\psi_{n+\hat{\mu}} +\bar{\psi}_{n+\hat{\mu}}(R+\gamma_\mu)\psi_n]\bigg], $$
where $a$ is the lattice spacing,  $\gamma^\mu$ are the usual pauli matrices $\mu$ runs from $1$ to $3$ which represents coordinates of 3 dimensional spatial space $x,y,z$. Now we want to make this locally gauge invariant under transformations,
$$\psi_n\to G_n\psi_n$$
where $G\in U(1)$ thus notice that if $\psi$ defined to transform as above we would also have the following transformation (just substitute above to LHS of below and you get the RHS of below),
$$\bar{\psi}_n\psi_{n+\hat{\mu}} +\bar{\psi}_{n+\hat{\mu}}\psi_n\to \bar{\psi}_nG^{-1}_{n}G_{n+\hat{\mu}}\psi_{n+\hat{\mu}} +\bar{\psi}_{n+\hat{\mu}}G^{-1}_{n+\hat{\mu}}G_{n}\psi_n$$ so if we define the locally gauge invariant Hamiltonian as,
$$H=a^2\bigg[(m_0a+3R)\sum_n\bar{\psi}_n\psi_n -\frac{1}{2}\sum_{n\mu}[\bar{\psi}_n(R-\gamma_\mu)U_{n,n+\hat{\mu}}\psi_{n+\hat{\mu}} +\bar{\psi}_{n+\hat{\mu}}(R+\gamma_\mu)U_{n,n+\hat{\mu}}^\dagger\psi_n]\bigg]$$
and define the transformation rule on $U$ as,
$$U_{n,n+\hat{\mu}}\to G(n)U_{n,n+\hat{\mu}}G^{-1}(n+\hat{\mu})$$
we see that the above Hamiltonian is indeed locally gauge invariant (Just apply the transformation rules we have just shown on $\psi$ and $U$ in that hamiltonian, and you see that it is indeed invariant and stay same). So we have derived the locally invariant hamiltonian. Finally, we can define another field $A$ to express $U$ more conveniently,
$$U=e^{iaA_\mu(n)}$$
A: I have been informed in the comments that I made a mistake here. However, since I consider it interesting, I want to keep my answer and only add an additional comment: In the following, one needs to be able to cancel the vector potential in the Hamiltonian with that from the phase. For that to be possible, one needs an indefinite integral of the vector potential in the phase. The existence of an indefinite integral, however, means that the vector potential is rotation-free and thus does not give rise to a magnetic field.
I interpret Peierls substitution as a means to handle an electric vector potential $\vec{A}$ in a Schrödinger equation (I explicitly do not write Hamiltonian, since further down, the time differential will be changed in the process, too, and one should not overlook that part). Assume the simplest SDE for this case, i.e.
$i\frac{\partial}{\partial t}\psi(\vec{r},t)
 =
 \left(
 \frac{\left(\vec{p}-\vec{A}(\vec{r},t)\right)^2}{2m}+U(\vec{r},t)
 \right)\psi(\vec{r},t).
$
To simplify it, one modifies the wave function via
$\psi(\vec{r},t)
 \rightarrow
 \psi(\vec{r},t)'=
 \psi(\vec{r},t)\cdot
 \exp\left(
  i\frac{q}{\hbar}
  \int_{\vec{r}_0}^{\vec{r}}
   \vec{A}(\vec{r}',t)
  d\vec{r}'
 \right)
$
with some $\vec{r}_0$ that does not matter for now. Upon this substitution, one obtains
$H(\vec{r},t)\psi(\vec{r},t)'\\
 =
 \left(
 \frac{\left(\vec{p}-\vec{A}(\vec{r},t)\right)^2}{2m}+U(\vec{r},t)
 \right)
 \exp\left(
 i\frac{q}{\hbar}
 \int_{\vec{r}_0}^{\vec{r}}
 \vec{A}(\vec{r}',t)
 d\vec{r}'
 \right)
 \psi(\vec{r},t)\\
 =
 \exp\left(
 i\frac{q}{\hbar}
 \int_{\vec{r}_0}^{\vec{r}}
 \vec{A}(\vec{r}',t)
 d\vec{r}'
 \right)
 \left(
 \frac{\left(\vec{p}+\vec{A}(\vec{r},t)-\vec{A}(\vec{r},t)\right)^2}{2m}+U(\vec{r},t)
 \right)
 \psi(\vec{r},t)\\
 =
 \exp\left(
 i\frac{q}{\hbar}
 \int_{\vec{r}_0}^{\vec{r}}
 \vec{A}(\vec{r}',t)
 d\vec{r}'
 \right)
 \left(
 \frac{\vec{p}^2}{2m}+U(\vec{r},t)
 \right)
 \psi(\vec{r},t)
$
and
$i\frac{\partial}{\partial t}\psi(\vec{r},t)'=
 \exp\left(
 i\frac{q}{\hbar}
 \int_{\vec{r}_0}^{\vec{r}}
 \vec{A}(\vec{r}',t)
 d\vec{r}'
 \right)
 \left(-q\Phi(\vec{r},t)\right)~
 i\hbar\frac{\partial}{\partial t}\psi(\vec{r},t).
$
with the new scalar potential
$
\Phi(\vec{r},t)=
 \int_{\vec{r}_0}^{\vec{r}}
 \frac{\partial\vec{A}(\vec{r}',t)}{\partial t}
 d\vec{r}'.
$
Plugging together the equations above, one obtains the new SDE
$
\quad\exp\left(
 i\frac{q}{\hbar}
 \int_{\vec{r}_0}^{\vec{r}}
 \vec{A}(\vec{r}',t)
 d\vec{r}'
 \right)\cdot
 i\hbar\frac{\partial}{\partial t}\psi(\vec{r},t)\\
 =
 \exp\left(
 i\frac{q}{\hbar}
 \int_{\vec{r}_0}^{\vec{r}}
 \vec{A}(\vec{r}',t)
 d\vec{r}'
 \right)\cdot
 \left(
  \frac{\vec{p}^2}{2m}
  +\Phi(\vec{r},t)
  +U(\vec{r},t)
 \right)
 \psi(\vec{r},t),
$
where the new scalar potential is ultimately a gauge transformation originating from the vector potential. For a time-independent potential $\vec{A}$, the new $\Phi$ is zero. For a time-dependent one, however, this transformation just changes an SDE involving a vector potential into one involving a scalar potential. And yes, this scalar potential needs to be taken into account when working with this substitution.
To apply this to a condensed matter system, one should take a look at the Wikipedia proof mentioned in the initial statement. Basically, it adds a sum over lattice sites and gives the lower integration limit a meaning. This will then lead to the modified hopping terms.
As a final note, J. M. Luttinger, Phys. Rev. 84, 814 (1951) shows that Peierls Substitution is a good approximation for slowly varying fields. For fast variation, however, further precautions are necessary. So people have been aware of this issue back then already.
