The standard way to treat a tight binding method in a magnetic is to replace the hopping matrix element:
$t_{i,j}\rightarrow e^{i\int_i^j\mathbf{A(x)}.d\mathbf{x}}$
the so called "Peierls substitution", How is this justified?
The usual way people try to justify this is to perform a consistency check (i.e. to show it's gauge invariant). Is there a more rigorous way to drive this? What are it's limits of validity?
The Hamiltonian is given by
\begin{equation}
H=\frac{\mathbf{p}^2}{2m}+U\left(\mathbf{r}\right),
\end{equation}
where
$U\left(\mathbf{r}\right)$
is the potential landscape due to the crystal lattice. The Bloch
theorem asserts that the solution to the problem
\begin{equation}
H\Psi_{\mathbf{k}}=E\left(\mathbf{k}\right)\Psi_{\mathbf{k}},
\end{equation}
is to be sought in the Bloch sum form
\begin{equation}
\Psi_{\mathbf{k}}=\frac{1}{\sqrt{N}}\sum_{\mathbf{R}}e^{i\mathbf{k}\cdot\mathbf{R}}\phi\left(\mathbf{r}-\mathbf{R}\right),
\end{equation}
where $N$ is the number of unit cells, and $\phi$ are atomic orbitals, also known as Wannier states. The corresponding eigenvalues $E\left(\mathbf{k}\right)$, which form bands
depending on the crystal momentum $\mathbf{k}$, are obtained by calculating
the matrix element $\langle\Psi_{\mathbf{k}}|H|\Psi_{\mathbf{k}}\rangle$
\begin{equation} \frac{1}{N}\sum_{\mathbf{R}\mathbf{R}^{\prime}}e^{i\mathbf{k}\left(\mathbf{R}^{\prime}-\mathbf{R}\right)} \int d\mathbf{r}\phi^*\left(\mathbf{r}-\mathbf{R}\right)H\phi\left(\mathbf{r}-\mathbf{R}^{\prime}\right) \end{equation} and ultimately depend on material-related hopping integrals $t_{12}=-\int d\mathbf{r}\phi^*\left(\mathbf{r}-\mathbf{R_1}\right)H\phi\left(\mathbf{r}-\mathbf{R_2}\right)$.
In the presence of the magnetic field the Hamiltonian changes to
\begin{equation} H=\frac{\left(\mathbf{p}-q\mathbf{A}\right)^2}{2m}+U\left(\mathbf{r}\right), \end{equation} where $q$ is the charge of the particle. The additional term introduces complications, and the original Bloch sum becomes inadequate. As it turns out, simply adding a phase term
\begin{equation} \Psi_{\mathbf{k}}=\frac{1}{\sqrt{N}}\sum_{\mathbf{R}}e^{i\left(\mathbf{k}\cdot\mathbf{R}+\frac{q}{\hbar}G_{\mathbf{R}}\right)}\phi\left(\mathbf{r}-\mathbf{R}\right), \end{equation} where
\begin{equation} G_{\mathbf{R}}=\int_{\mathbf{R}}^{\mathbf{r}}\mathbf{A}\cdot d\mathbf{l}, \end{equation} resolves the issue. The hopping matrix elements now read
\begin{align}
\langle\Psi_{\mathbf{k}}|H|\Psi_{\mathbf{k}}\rangle=&\frac{1}{N}\sum_{\mathbf{R}\mathbf{R}^{\prime}}e^{i\mathbf{k}\left(\mathbf{R}^{\prime}-\mathbf{R}\right)}\int d\mathbf{r}e^{-i\frac{q}{\hbar}G_{\mathbf{R}}}\phi^*\left(\mathbf{r}-\mathbf{R}\right)\left[\frac{\left(\mathbf{p}-q\mathbf{A}\right)^2}{2m}+U\left(\mathbf{r}\right)\right]e^{i\frac{q}{\hbar}G_{\mathbf{R}^{\prime}}}\phi\left(\mathbf{r}-\mathbf{R}^{\prime}\right)\nonumber\\
=&\frac{1}{N}\sum_{\mathbf{R}\mathbf{R}^{\prime}}e^{i\mathbf{k}\left(\mathbf{R}^{\prime}-\mathbf{R}\right)}e^{i\frac{q}{\hbar}\int_{\mathbf{R}^{\prime}}^{\mathbf{R}}\mathbf{A}\cdot d\mathbf{l}}\nonumber\\&\times\int d\mathbf{r}e^{i\frac{q}{\hbar}\Phi\left(\mathbf{r}\right)}\phi^*\left(\mathbf{r}-\mathbf{R}\right)\left[\frac{\left(\mathbf{p}-q\mathbf{A}+q\nabla G_{\mathbf{R}^{\prime}}\right)^2}{2m}+U\left(\mathbf{r}\right)\right]\phi\left(\mathbf{r}-\mathbf{R}^{\prime}\right)\nonumber\\
=&\frac{1}{N}\sum_{\mathbf{R}\mathbf{R}^{\prime}}e^{i\mathbf{k}\left(\mathbf{R}^{\prime}-\mathbf{R}\right)}e^{i\frac{q}{\hbar}\int_{\mathbf{R}^{\prime}}^{\mathbf{R}}\mathbf{A}\cdot d\mathbf{l}}\int d\mathbf{r}\phi^*\left(\mathbf{r}-\mathbf{R}\right)\left[\frac{\mathbf{p}^2}{2m}+U\left(\mathbf{r}\right)\right]\phi\left(\mathbf{r}-\mathbf{R}^{\prime}\right).
\end{align}
The relation $\nabla G_{\mathbf{R}^{\prime}}=\mathbf{A}$ holds for the tight
binding condition and in the case when the magnetic field is invariant at the
scale of the crystal lattice. On the other hand, the flux
$\Phi\left(\mathbf{r}\right)=\oint_{\mathbf{R}^{\prime}\rightarrow\mathbf{r}\rightarrow\mathbf{R}}\mathbf{A}\cdot
d\mathbf{l}$ is larger when the integrand $\mathbf{r}$ is further from the two
vectors $\mathbf{R}$ and $\mathbf{R}^{\prime}$, where the atomic orbitals Wannier states are
effectively zero, while the flux is vanishing where the hopping integral is
nonzero. Having these two things in mind helps explain the transition from the
second to the third line.
Now it becomes obvious that the matrix elements are the same as in the case without magnetic field, apart from the phase factor picked up, which is denoted the Peierls phase. This is tremendously convenient, since then we get to use the same material parameters regardless of the magnetic field value, and the corresponding phase is computationally trivial to take into account. For electrons it amounts to replacing the hopping term $t_{ij}$ with $t_{ij}e^{i\frac{e}{\hbar}\int_i^j\mathbf{A}\cdot d\mathbf{l}}$. Finally, note that a beautiful and elucidating explanation for this phase can also be found in Feynman's Lectures (Vol. III, Chapter 21).
