Integrals of motion for one-dimensional Anderson localized systems

All single particle states of Type-1 Hamiltonians (6), except possibly the ground state for $y\gt 0$ or the highest energy state for $y\lt 0$ are localized, see e.g. figure 1.

Zoom In Zoom Out Reset image size

This can be understood in more detail from the exact solution for the spectrum of these models [27]. Exact un-normalized single particle eigenstates of the Hamiltonian (6) are

$$\begin{eqnarray}&&| E\rangle =\displaystyle \sum _{i=0}^{N-1}\displaystyle \frac{{\gamma }_{i}{c}_{i}^{\dagger }}{E-{\epsilon }_{i}}| 0\rangle \end{eqnarray} \tag{ 7 }$$

and the corresponding eigenvalues E (energies) are solutions of the equation

$$\begin{eqnarray}&&\displaystyle \sum _{i=0}^{N-1}\displaystyle \frac{{\gamma }_{i}^{2}}{E-{\epsilon }_{i}}=-\displaystyle \frac{1}{y}.\end{eqnarray} \tag{ 8 }$$

Suppose ${\epsilon }_{i}$ are ordered in the ascending order. By plotting the left hand side of equation (8) as a function of E, one can verify that it has $N-1$ real roots ${E}_{1},{E}_{2},\ldots {E}_{N-1}$ located between consecutive ${\epsilon }_{i}$, i.e. ${\epsilon }_{i-1}\lt {E}_{i}\lt {\epsilon }_{i}$. The remaining root E₀ is also real and is below ${\epsilon }_{0}$ (ground state) for $y\gt 0$ and above ${\epsilon }_{N-1}$ for $y\lt 0$ (highest excited state).

Equations (7) and (8) also provide an exact solution for one fermion (Cooper) pair and one spin flip sectors of the BCS and Gaudin models, respectively

$$\begin{eqnarray}{H}_{\mathrm{BCS}} & = & \displaystyle \sum _{i,\sigma =\uparrow ,\downarrow }{\epsilon }_{i}{c}_{i\sigma }^{\dagger }{c}_{i\sigma }-y\displaystyle \sum _{{ij}}{c}_{i\downarrow }^{\dagger }{c}_{i\uparrow }^{\dagger }{c}_{j\uparrow }{c}_{j\downarrow },\\ & & {H}_{i}(y)={s}_{i}^{z}-y\displaystyle \sum _{j\ne i}\displaystyle \frac{{\vec{s}}_{i}\cdot {\vec{s}}_{j}}{{\epsilon }_{i}-{\epsilon }_{j}},\end{eqnarray} \tag{ 9 }$$

where ${c}_{i\sigma }$ are spin-full fermions and ${\vec{s}}_{i}$ are quantum spins of arbitrary magnitudes s_i, see [27] for details. For the BCS (Gaudin) model one needs to replace ${c}_{i}^{\dagger }\to {c}_{i\downarrow }^{\dagger }{c}_{i\uparrow }^{\dagger }$ ${\rm{[}}{s}_{i}^{+}{\rm{]}}$ in equation (7), set ${\gamma }_{i}=1$ ${\rm{[}}\sqrt{2{s}_{i}}{\rm{]}}$ and the corresponding eigenvalue is equal to $2E$ ${\rm{[}}2{s}_{i}{(E-{\epsilon }_{i})}^{-1}{\rm{]}}$ rather than E. Our results for the PR of Type-1 Hamiltonians therefore also apply to these sectors of these models.

The PR defined through equation (14) reads

$$\begin{eqnarray}&&{\mathrm{PR}}_{E}=\displaystyle \frac{{\left[{\displaystyle \sum }_{i}\displaystyle \frac{{\gamma }_{i}^{2}}{{(E-{\epsilon }_{i})}^{2}}\right]}^{2}}{{\displaystyle \sum }_{i}\displaystyle \frac{{\gamma }_{i}^{4}}{{(E-{\epsilon }_{i})}^{4}}}.\end{eqnarray} \tag{ 10 }$$

For concreteness we take $y\gt 0$. Then, the ground state is ${E}_{0}\lt {\epsilon }_{0}$. We assume that most ${\gamma }_{i}$ are of the same order of magnitude and consequently the vector with components ${\gamma }_{i}$ is delocalized. Further, we take ${\epsilon }_{i}$ to lie in a fixed interval that does not scale with N, e.g. from $-w$ to w.

For excited states E_k is between ${\epsilon }_{k-1}$ and ${\epsilon }_{k}$. The summations in the numerator and denominator of equation (10) both come from ${\epsilon }_{i}$ in a small vicinity of ${\epsilon }_{k}$ for large N and converge as ${\sum }_{n}{n}^{-2}$ and ${\sum }_{n}{n}^{-4}$, respectively, where $n=| i-k|$. The numerator and the denominator scale as ${[{\gamma }_{k}^{2}/{\delta }^{2}]}^{2}$ and ${\gamma }_{i}^{4}/{\delta }^{4}$, where $\delta \propto 1/N$ is the mean level spacing between ${\epsilon }_{i}$ in the vicinity of ${\epsilon }_{k}$. Therefore, ${\mathrm{PR}}_{{E}_{k}}$ is of order 1 (much smaller than N) meaning excited states are always localized. Figure 1 shows PR for $N={10}^{3}$ uncorrelated random ${\epsilon }_{i}$ uniformly drawn from an interval $(-1,1)$ and the same distribution of ${\gamma }_{i}$.

Consistent with our numerical results, we estimate the largest PR for excited states to scale as $\mathrm{ln}N$, i.e.

$$\begin{eqnarray}&&{\mathrm{PR}}_{{E}_{k}}^{\mathrm{max}}\approx \alpha \mathrm{ln}N,\end{eqnarray} \tag{ 11 }$$

for large N, where α depends on N much weaker than $\mathrm{ln}N$. Such values of PR come from clustering in ${\epsilon }_{i}$. Indeed, suppose spacings ${\delta }_{i}={\epsilon }_{i+1}-{\epsilon }_{i}$ between m of ${\epsilon }_{i}$ for i from k to k + m are all much smaller than ${\delta }_{k-1}$ and, moreover, ${\epsilon }_{k+m}-{\epsilon }_{k}\leqslant {\delta }_{k-1}$. It follows from equation (10) that ${\mathrm{PR}}_{{E}_{k}}\approx {\mathrm{PR}}_{{E}_{k+m+1}}\approx m$ because the above ${\epsilon }_{i}$ contribute most to these PRs. Normalized spacings ${s}_{i}={\delta }_{i}/\delta$ are distributed according to the Poisson distribution $P(s){\rm{d}}s={{\rm{e}}}^{-s}{\rm{d}}s$. The probability of having m spacings between 0 and ${s}_{0}\ll 1$ is then roughly s₀^m. We need ${{ms}}_{0}\leqslant 1$ and also ${{Ns}}_{0}^{m}=1$ so that at least one such clustering occurs⁸ . This implies $m\approx \mathrm{ln}N/\mathrm{ln}(\mathrm{ln}N)$ and equation (11) follows. Numerically we find that typical values of $\alpha \approx 1-3$ and averaged over disorder $\bar{\alpha }\approx 1.7$, at least for N = 2⁴ – 2¹². Note that according to this argument such large values of PR typically come in pairs spaced by $m+1$, roughly equal to the value of the PR itself. We also stress that, in contrast to the largest PR, a typical (and average) PR is something between one and three for any N (does not scale) as can be seen from figure 1.

It is interesting to compare this $\mathrm{ln}N$ behavior to the flat band localization studied earlier [44, 45]. The latter leads to a (weakly) divergent PR in the localized regime, a phenomenon that is viewed as corresponding to critical (power law type) localization. The Type-1 Hamiltonian kinetic energy may also be viewed as a 'flat band' model, with a flat dispersion for all except one state. Indeed, for ${t}_{{ij}}={\gamma }_{i}{\gamma }_{j}$ all but one eigenvalues of the second term in equation (2) are zero. The non-zero eigenvalue (ground state for $y\gt 0$) corresponds to the eigenstate ${\gamma }_{i}{c}_{i}^{\dagger }| 0\rangle$.

Let us consider limits $y\to 0$ and $y\to \infty$ separately. When $y\to 0$ all states are localized as expected. Indeed, equation (8) implies ${E}_{k}\to {\epsilon }_{k}$, summations in equation (10) are dominated by the i = k term and we obtain ${\mathrm{PR}}_{{E}_{k}}=1$ for all k. When $y\to \infty$ excited states are localized as before because E_k for $k\geqslant 1$ remains trapped in the interval $({\epsilon }_{k-1},{\epsilon }_{k})$. The ground state energy on the other hand diverges— equation (8) implies ${E}_{0}\to -y{\sum }_{i}{\gamma }_{i}^{2}$. Then, ${\epsilon }_{i}$ are negligible as compared to E₀ in equation (10) and

$$\begin{eqnarray}&&{\mathrm{PR}}_{{E}_{0}}=\displaystyle \frac{{\left[{\displaystyle \sum }_{i}{\gamma }_{i}^{2}\right]}^{2}}{{\displaystyle \sum }_{i}{\gamma }_{i}^{4}},\end{eqnarray} \tag{ 12 }$$

which is of order N according to our choice of ${\gamma }_{i}$. The ground state is therefore delocalized for $y\to \infty$. It undergoes a localization–delocalization crossover at a certain y_c, which we estimate below in this section.

It is possible to evaluate the PR analytically to leading order in $1/N$ for distributions of ${\epsilon }_{i}$ and ${\gamma }_{i}$ with negligible short range fluctuations (such that the spacing ${\delta }_{i}={\epsilon }_{i+1}-{\epsilon }_{i}$ changes slowly with i—$| {\delta }_{i+1}-{\delta }_{i}| /{\delta }_{i}$ is of order $1/N$ for all i—and similarly for ${\gamma }_{i}$). For simplicity, let us take constant ${\gamma }_{i}$, which we can set to one with no loss of generality, and equally spaced ${\epsilon }_{i}$, i.e. ${\delta }_{i}=\delta =2w/N$.

For excited states, we write ${E}_{k}={\epsilon }_{k}-{\alpha }_{k}\delta$, where $0\lt {\alpha }_{k}\lt 1$, and solve equation (8) for ${\alpha }_{k}$ to the leading order in $1/N$ as described in appendix B of [46]. This yields

$$\begin{eqnarray}&&\mathrm{cot}\pi {\alpha }_{k}=\displaystyle \frac{\delta }{\pi y}+\displaystyle \frac{1}{\pi }\mathrm{ln}\displaystyle \frac{{\epsilon }_{k}+w}{w-{\epsilon }_{k}}\equiv f({\epsilon }_{k}).\end{eqnarray} \tag{ 13 }$$

We note that $\lambda =y/\delta$ is the proper dimensionless coupling constant in the sense that it must stay finite in the $N\to \infty$ limit. This is because the second summation in equation (2) scales as N² for ${t}_{{ij}}={\gamma }_{i}{\gamma }_{j}$ and our choice of ${\gamma }_{i}$. Therefore, we need $y\propto \delta \propto 1/N$ so that both terms in equation (2) are extensive in the thermodynamic limit. For the BCS Hamiltonian in equation (9), so defined λ is the dimensionless superconducting coupling [47].

Equations (10) becomes to leading order in $1/N$

$$\begin{eqnarray}&&{\mathrm{PR}}_{{E}_{k}}=\displaystyle \frac{{\left[{\displaystyle \sum }_{n=0}^{\infty }\left(\displaystyle \frac{1}{{(n+{\alpha }_{k})}^{2}}+\displaystyle \frac{1}{{(n+1-{\alpha }_{k})}^{2}}\right)\right]}^{2}}{{\displaystyle \sum }_{n\quad =\quad 0}^{\infty }\left(\displaystyle \frac{1}{{(n+{\alpha }_{k})}^{4}}+\displaystyle \frac{1}{{(n+1-{\alpha }_{k})}^{4}}\right)},\end{eqnarray} \tag{ 14 }$$

which evaluates to

$$\begin{eqnarray}&&{\mathrm{PR}}_{{E}_{k}}=\displaystyle \frac{3}{1\quad +\quad 2{\mathrm{cos}}^{2}\pi {\alpha }_{k}}=\displaystyle \frac{3\quad +\quad 3{f}^{2}({\epsilon }_{k})}{1\quad +\quad 3{f}^{2}({\epsilon }_{k})}.\end{eqnarray} \tag{ 15 }$$

This answer is in good agreement with numerics already for N = 20, see also figure 1. Note that $1\leqslant {\mathrm{PR}}_{{E}_{k}}\leqslant 3$.

We saw above that the ground state energy ${E}_{0}\to -\infty$ as $y\to \infty$, while ${E}_{0}\to {\epsilon }_{0}$ for $y\to 0$. Let y be large enough that E₀ is well separated from ${\epsilon }_{0}$. Then, we can replace summation in equation (8) with integration and obtain

$$\begin{eqnarray}&&\mathrm{ln}\displaystyle \frac{{E}_{0}-w}{{E}_{0}+w}=\displaystyle \frac{\delta }{y}=\displaystyle \frac{2w}{{Ny}}.\end{eqnarray} \tag{ 16 }$$

Performing the same replacement in equation (10) and using equation (16), we derive

$$\begin{eqnarray}&&{\mathrm{PR}}_{{E}_{0}}=\displaystyle \frac{3N}{1+2\mathrm{cosh}(\delta /y)}.\end{eqnarray} \tag{ 17 }$$

Note that in the limit $y\to \infty$, ${\mathrm{PR}}_{{E}_{0}}=N$ in agreement with equation (12). This expression also allows us to estimate the value y_c beyond which the ground state becomes extended. We obtain ${\lambda }_{{\rm{c}}}={y}_{{\rm{c}}}/\delta \approx 1/\mathrm{ln}N$. This also corresponds to the coupling for which the gap in the spectrum ${\rm{\Delta }}={E}_{1}-{E}_{0}\approx -w-{E}_{0}$ becomes comparable to the spacing δ. For a superconductor described by the BCS model (9) this localized-extended crossover translates into a normal-superconducting one[48, 49]. As $N\to \infty$ this crossover becomes a quantum phase transition at $\lambda =0$, i.e. any infinitesimal coupling is sufficient to make the ground state extended (superconducting). The localized character of the excited states for the specific case of ${\gamma }_{i}=1$ has been demonstrated in a previous work as well [43].

Proceeding as for the case of Type-1 Hamiltonians, we focus on the conserved charge Q₀, corresponding to the site i = 0, which to lowest order is equal to n₀. However, in this case Q₀ is not simply linear in y. In fact, it can be argued that the an expansion of Q₀ in the hopping does not truncate at any finite order in the thermodynamic limit. Indeed, as explained in the Introduction, conserved charges are generally infinite power series in y. We thus assume Q_i of the form

$$\begin{eqnarray}&&{Q}_{i}={P}_{i0}+{{yP}}_{i1}+{y}^{2}{P}_{i2}+\cdots ,\end{eqnarray} \tag{ 19 }$$

where ${P}_{i0}={n}_{0}$ and ${P}_{i1},{P}_{i2}\ldots$ are operators to be determined in terms of the microscopic parameters subject to the condition $[{Q}_{i},H]=0$. For concreteness, we first take our one-dimensional system to be a finite-sized ring of $N+1$ sites going from 0 to N.

Since the Hamiltonian H and all the zero order charges n_i are quadratic in the creation and annihilation operators, we take all the operators ${P}_{i1},{P}_{i2}...$ to be similarly quadratic, i.e.

$$\begin{eqnarray}&&{P}_{{im}}=\displaystyle \sum _{{jk}}{\eta }_{{jk}}^{(m)}(i){c}_{k}^{\dagger }{c}_{j},\end{eqnarray} \tag{ 20 }$$

where the symmetric coefficients ${\eta }_{{jk}}^{(m)}(i)={\eta }_{{kj}}^{(m)}(i)$ are to be determined. We have

$$\begin{eqnarray*}&&[{Q}_{i},H]=[{P}_{i0},{H}_{0}]+\displaystyle \sum _{m}{y}^{m+1}([{P}_{{im}},{H}_{1}]+[{P}_{{im}+1},{H}_{0}]).\end{eqnarray*}$$

The requirement that the commutator vanishes to all orders in y requires

$$\begin{eqnarray}&&[{P}_{{im}},{H}_{1}]+[{P}_{{im}+1},{H}_{0}]=0\end{eqnarray} \tag{ 21 }$$

and yields a recursion relation among $\eta ^{\prime} s$

$$\begin{eqnarray}&&{\eta }_{{ab}}^{(m+1)}(i)={\delta }_{{ab}}{R}_{a}^{(m+1)}(i)+\displaystyle \frac{1-{\delta }_{{ab}}}{{\epsilon }_{a}-{\epsilon }_{b}}\displaystyle \sum _{j}[({t}_{{aj}}{\eta }_{{jb}}^{(m)}(i)-{\eta }_{{aj}}^{(m)}(i){t}_{{jb}}],\end{eqnarray} \tag{ 22 }$$

with initial conditions ${\eta }_{{ab}}^{(0)}(i)={\delta }_{{ia}}{\delta }_{{ib}}$. The diagonal term ${R}_{a}^{(m+1)}$ represents a 'gauge' freedom, since the corresponding term in ${P}_{{im}}$ commutes trivially with H₀. We further discuss this freedom below. Specializing to the case of nearest neighbor hopping equation (18) and with i = 0, it can be verified that terms present in ${P}_{0m}$ are of the form

• ${\eta }_{0\;m}^{(m)}({c}_{0}^{\dagger }{c}_{m}+{c}_{m}^{\dagger }{c}_{0}),$
• ${\eta }_{0,N-(m-1)}^{(m)}[{c}_{0}^{\dagger }{c}_{N-(m-1)}+{c}_{N-(m-1)}^{\dagger }{c}_{0}],$
• ${\sum }_{{ij}\ne m,| i-j| =\mathrm{even}\leqslant m}{\eta }_{{ij}}^{(m)}({c}_{i}^{\dagger }{c}_{j}+{c}_{j}^{\dagger }{c}_{i})$ (if m is even),
• ${\sum }_{{ij}\ne m,| i-j| =\mathrm{odd}\leqslant m}{\eta }_{{ij}}^{(m)}({c}_{i}^{\dagger }{c}_{j}+{c}_{j}^{\dagger }{c}_{i})$ (if m is odd).

This is shown schematically in figure 2 for the first few P_0m.

Zoom In Zoom Out Reset image size

The ${Q}_{i}^{\prime }s$ are related to each other by translating all site indices in the above relations by an appropriate number. By construction, they all commute with H. Since H is generally non-degenerate, this implies Q_i also commute among themselves, $[{Q}_{i},{Q}_{j}]=0$ ∀ $i,j$. To see this, first recall that for Hermitian matrices $[A,B]=[A,C]=0$ implies $[B,C]=0$ as long as eigenvalues of A are non-degenerate. All operators involved in the above construction of Q_i are of the form $\hat{A}={\sum }_{{ij}}{A}_{{ij}}{c}_{i}^{\dagger }{c}_{j}$, where A_ij is a Hermitian N × N matrix, which represents operator $\hat{A}$ in the sector with total particle number n = 1. Moreover, the commutativity of any two such operators is equivalent to that of the underlying matrices. Eigenvalues of the Hamiltonian in the n = 1 sector at y = 0 are ${\epsilon }_{i}$, which are assumed to be distinct, i.e. the corresponding matrix is non-degenerate at y = 0. By continuity of the eigenvalues in y, it remains non-degenerate in some finite interval (until the first level crossing) of the real axis containing y = 0. Thus, $[{Q}_{i},{Q}_{j}]=0$ ∀ $i,j$ in this interval of y. But, as can be seen e.g. from the above construction of Q₀, commutativity of Q_i on any finite interval of values of y implies that they commute for all y.

We noted above that commutation relations (21) and consequently recursion relations (22) do not constrain the diagonal part of the coefficients ${\eta }^{(m)}$, i.e. ${R}_{r}^{(m)}$, for $m\geqslant 1$. The choice of ${R}_{a}^{(m)}(i)$ however does affect the off-diagonal part of ${\eta }^{(k)}$ for $k\gt m$. In our construction of Q_i we set ${R}_{a}^{(m)}(i)=0$ for all $m\geqslant 1$, since this leads to the most compact description of these objects. we will refer to this as the standard gauge. Conserved charges ${\tilde{Q}}_{i}$ resulting from any other choice${R}_{a}^{(m)}(i)$ uniquely relate to our standard gauge Q_i's, a brief calculation shows their relationship is

$$\begin{eqnarray}&&{\tilde{Q}}_{i}={Q}_{i}+\displaystyle \sum _{m}{y}^{m}\displaystyle \sum _{r}{R}_{r}^{(m)}(i){Q}_{r}.\end{eqnarray} \tag{ 23 }$$

Another advantage of our choice of a gauge is in a simple relationship between the Hamiltonian (18) and the conserved charges, namely

$$\begin{eqnarray}&&H=\displaystyle \sum _{i}{\epsilon }_{i}{Q}_{i}.\end{eqnarray} \tag{ 24 }$$

To see this, consider the difference

$$\begin{eqnarray}&&H-\displaystyle \sum _{i}{\epsilon }_{i}{Q}_{i}={{yW}}_{1}\quad +\quad {y}^{2}{W}_{2}\quad +\quad {y}^{3}{W}_{3}+\cdots ,\end{eqnarray} \tag{ 25 }$$

where W_i are y-independent operators. Note that the zeroth order term cancels in the difference. Since H commutes with all Q_i, the right -hand side (rhs) of equation (25) must also commute. This implies in particular $[{W}_{1},{n}_{i}]=0$ for all i (from the coefficient at the lowest power of y in the commutator of the rhs with Q_i), which in turn means that ${W}_{1}={\sum }_{i}{r}_{i}^{1}{n}_{i}$. Now note that the left-hand side (lhs) has zero diagonal matrix elements, i.e. no terms of the form ${c}_{r}^{\dagger }{c}_{r}$. This is because the zeroth order term is absent, while higher order terms have no diagonal matrix elements since ${\eta }_{{rr}}^{(m)}=0$ for all $m\geqslant 1$ in our gauge (and similarly the diagonal is absent in other Q_i). Then, the diagonal matrix elements must vanish on the rhs as well, to all orders in y. In particular, ${r}_{i}^{1}=0$, i.e. ${W}_{1}=0$ and

$$\begin{eqnarray}&&H-\displaystyle \sum _{i}{\epsilon }_{i}{Q}_{i}={y}^{2}{W}_{2}\quad +\quad {y}^{3}{W}_{3}+\cdots .\end{eqnarray} \tag{ 26 }$$

Applying the same argument to the rhs of this equation we similarly obtain ${W}_{2}=0$ etc., until we finally arrive at equation (24).

We have seen above that the conserved charges are power series in the hopping. This is unlike the case of Type-1 Hamiltonians, where the power series truncates after the first term. The gauge where the series truncates corresponds to having distinct terms for m = 1, one can see in equation (3) (the gauge terms are indicated in the lower braces).

It is an amusing exercise to determine the correct gauge terms that lead to truncation, starting from the recursion relations equation (22). To obtain Type-1 Hamiltonians we set ${t}_{{ij}}={\gamma }_{i}{\gamma }_{j}$, so that the recursions simplify to

$$\begin{eqnarray}{\eta }_{{ab}}^{(m+1)}(i) & = & {\delta }_{{ab}}{R}_{a}^{(m+1)}(i)-\displaystyle \frac{1-{\delta }_{{ab}}}{{\epsilon }_{a}-{\epsilon }_{b}}({Y}_{{ab}}^{(m)}(i)-{Y}_{{ba}}^{(m)}(i))\\ {Y}_{{ab}}^{(m)}(i) & = & \displaystyle \sum _{j}{\eta }_{{aj}}^{(m)}(i){\gamma }_{j}\;{\gamma }_{b}.\end{eqnarray} \tag{ 27 }$$

With the initial condition ${\eta }_{{ab}}^{(0)}(i)={\delta }_{{ia}}{\delta }_{{ib}}$, we obtain at the first level

$$\begin{eqnarray}&&{\eta }_{{ab}}^{(1)}(i)={\delta }_{{ab}}{R}_{a}^{(1)}(i)+\displaystyle \frac{{\gamma }_{a}{\gamma }_{b}(1-{\delta }_{{ab}})}{{\epsilon }_{a}-{\epsilon }_{b}}({\delta }_{{ib}}-{\delta }_{{ia}}).\end{eqnarray} \tag{ 28 }$$

At this point we pause and ask if we can choose the gauge term ${R}_{a}^{(1)}(i)$ such that ${\eta }_{{ab}}^{(2)}(i)$ can be made to vanish identically, so that the iterations stop at the first level. From equation (27) we see that the relevant condition is the vanishing of $({Y}_{{ab}}^{(1)}(i)-{Y}_{{ba}}^{(1)}(i))$. Using equation (27) compute

$$\begin{eqnarray}&&{Y}_{{ab}}^{(1)}(i)={\gamma }_{a}{\gamma }_{b}\left\{{R}_{a}^{(1)}(i)+\displaystyle \frac{{\gamma }_{i}^{2}}{{\epsilon }_{a}-{\epsilon }_{i}}(1-{\delta }_{{ia}})-{\delta }_{{ia}}\displaystyle \sum _{j}^{\prime }\displaystyle \frac{{\gamma }_{j}^{2}}{{\epsilon }_{j}-{\epsilon }_{i}}\right\}.\end{eqnarray} \tag{ 29 }$$

We may choose ${R}^{(1)}$ so that the term in braces vanishes, thus leading to the truncation of the iterations. From equation (28) we have the complete first order term, and we can proceed to construct the charge (denoting the currents by the symbol $\tilde{Q}$)

$$\begin{eqnarray}&&{\tilde{Q}}_{i}={n}_{i}+y\displaystyle \sum _{{ab}}{\eta }_{{ab}}^{(1)}(i){c}_{a}^{\dagger }{c}_{b},\end{eqnarray} \tag{ 30 }$$

which is identical to that in equation (3).

The use of the gauge term here is very special, and guided by our understanding of this model. On the other hand, we could by default set all the gauge terms ${R}^{(m)}$ to zero, giving us the irreducible (i.e. standard gauge) currents. These no longer truncate even for Type-1 Hamiltonians. For completeness we note the second order term for the current in this (standard) gauge

$$\begin{eqnarray}{Q}_{i} & = & {n}_{i}+(y+{y}^{2}\displaystyle \sum _{j}^{\prime }\displaystyle \frac{{\gamma }_{j}^{2}}{{\epsilon }_{j}-{\epsilon }_{i}})\times \displaystyle \sum _{i\ne j}\displaystyle \frac{{\gamma }_{i}{\gamma }_{j}}{{\epsilon }_{j}-{\epsilon }_{i}}({c}_{i}^{\dagger }{c}_{j}+{c}_{j}^{\dagger }{c}_{i})\\ & & +{y}^{2}{\gamma }_{i}^{2}\times \displaystyle \sum _{a,b}^{\prime }\displaystyle \frac{{\gamma }_{a}{\gamma }_{b}}{({\epsilon }_{a}-{\epsilon }_{i})({\epsilon }_{a}-{\epsilon }_{j})}{c}_{i}^{\dagger }{c}_{j}+O({y}^{3}).\end{eqnarray} \tag{ 31 }$$

Thus the Type-1 Hamiltonians allow for variety of expressions of the constants of motion. To establish their equivalence in general is a subtle problem, where some surprising results have been found quite recently in [30].

This type of gauge choice, made explicit in our construction could be exploited further to test the possibility that the series can take simpler forms, as compared to a brute force expansions to infinite order. We leave this interesting question for future investigation.

A natural question that arises is the relationship between the currents found above and those found from a brute force expansion of the projection operators of the Anderson model in powers of the coupling constant y. The model has a formal single particle eigenfunction expansion in the form

$$\begin{eqnarray}&&| {\rm{\Psi }}(y)\rangle =\displaystyle \sum _{k}{u}_{0k}(y){c}_{k}^{\dagger }| 0\rangle ,\end{eqnarray} \tag{ 32 }$$

with an initial condition localized say at the site 0 as ${u}_{0k}(0)={\delta }_{k0}$. The projector $Q=| {\rm{\Psi }}(y)\rangle \langle {\rm{\Psi }}(y)|$ can be expanded in a series in y

$$\begin{eqnarray}&&\hat{Q}=\displaystyle \sum _{j,k}{u}_{0j}{u}_{0k}^{*}{c}_{j}^{\dagger }{c}_{k}=\hat{P}(0)+y\hat{P}(1)+{y}^{2}\hat{P}(2)+\cdots \end{eqnarray} \tag{ 33 }$$

so that the basic expansion of the wave functions in a Rayleigh–Schrödinger (RS) series in y generates the conserved currents. We can use the standard result in text books⁹ to write a perturbative expansion for the state at site 0 with standard normalization to ${u}_{00}=1$ as

$$\begin{eqnarray}&&| {\rm{\Psi }}(y)\rangle ={c}_{0}^{\dagger }| 0\rangle +\displaystyle \sum _{k\ne 0}{u}_{0k}{c}_{k}^{\dagger }| 0\rangle ,\end{eqnarray} \tag{ 34 }$$

with a power series expansion for u_0k

$$\begin{eqnarray}&&{u}_{0k}=-y\displaystyle \frac{{t}_{0k}}{{\epsilon }_{0}-{\epsilon }_{k}}+{y}^{2}\displaystyle \sum _{l\ne 0}\displaystyle \frac{{t}_{{kl}}{t}_{l0}}{({\epsilon }_{0}-{\epsilon }_{k})({\epsilon }_{0}-{\epsilon }_{l})}-{y}^{2}\displaystyle \frac{{t}_{00}{t}_{k0}}{{({\epsilon }_{0}-{\epsilon }_{k})}^{2}}+O({y}^{3}).\end{eqnarray} \tag{ 35 }$$

Using this expansion, we may generate the series equation (33), the result is explicitly stated below in equation (38). From this series we can verify to second order, that this series differs from that in the standard gauge equation (19) by specific gauge terms. The advantage of equation (19) is that this gauge invariance is manifest in the construction by the nested commutators. On the other hand, equations (33)–(35), corresponds to a particular gauge picked out by the R–S method, and the currents found here are some linear combinations of the ones in equation (19) as in equation (23).

It seems to us that the series in equation (19) possesses an essential simplicity relative to the Rayleigh–Schrödinger series equations (33)–(35). The R–S perturbation expansion simultaneously determines the energy eigenvalue, and for this purpose very specific gauge terms are needed. On the other hand all terms in equation (19) are generated by completely off diagonal terms, those terms that avoid multiple visits to any site. This leads to simpler recursion relations, as in equation (22), relative to the RS series. For this reason our numerical work in this paper uses the series in equation (19).

The Rayleigh–Schrödinger series can also be constructed for Type-1 Hamiltonians using the exact eigenstates $| E\rangle$ with eigenvalues E as given in equations (7) and (8). The projector $| E\rangle \langle E| ={\sum }_{{ij}}\frac{{\gamma }_{i}{\gamma }_{j}{c}_{i}^{\dagger }{c}_{j}}{(E-{\epsilon }_{i})(E-{\epsilon }_{j})}$ can be expanded in y as shown in equation (33). In the limit $y\to 0$, the roots of equation (8) tend to ${\epsilon }_{i}$. We take the root $E\to {\epsilon }_{0}$ to obtain ${\tilde{Q}}_{0}$ the conserved charge corresponding to site 0 calculated using the Rayleigh–Schrödinger gauge. Other roots yield other ${\tilde{Q}}_{i}$. Expanding equation (8) for E in y near $E={\epsilon }_{0}$, we get

$$\begin{eqnarray}&&E={\epsilon }_{0}-y{\gamma }_{0}^{2}\quad +\quad {y}^{2}{\gamma }_{0}^{2}\displaystyle \sum _{i\ne 0}\displaystyle \frac{{\gamma }_{i}^{2}}{{\epsilon }_{0}-{\epsilon }_{i}}+O({y}^{3}).\end{eqnarray} \tag{ 36 }$$

Since, the projector diverges in $y\to 0$ limit, we define our conserved charge as ${\tilde{Q}}_{0}=\frac{{(E-{\epsilon }_{0})}^{2}}{{\gamma }_{0}^{2}}| E\rangle \langle E|$ to make it well behaved. ${\tilde{Q}}_{0}$ is given by

$$\begin{eqnarray}&&{\tilde{Q}}_{0}={n}_{0}+\displaystyle \frac{E-{\epsilon }_{0}}{{\gamma }_{0}}\displaystyle \sum _{j\ne 0}\displaystyle \frac{{\gamma }_{j}({c}_{0}^{\dagger }{c}_{j}+{c}_{j}^{\dagger }{c}_{0})}{E-{\epsilon }_{j}}+\displaystyle \frac{{(E-{\epsilon }_{0})}^{2}}{{\gamma }_{0}^{2}}\displaystyle \sum _{i,j\ne 0}\displaystyle \frac{{\gamma }_{i}{\gamma }_{j}{c}_{i}^{\dagger }{c}_{j}}{(E-{\epsilon }_{i})(E-{\epsilon }_{j})}.\end{eqnarray} \tag{ 37 }$$

Then, replacing $\frac{{(E-{\epsilon }_{0})}^{2}}{{\gamma }_{0}^{2}}\to {y}^{2}{\gamma }_{0}^{2}$ and then $E\to {\epsilon }_{0}$, we have obtained ${\tilde{Q}}_{0}$ as a combination of Q_i (see equation (3)) as follows

$$\begin{eqnarray}&&{\tilde{Q}}_{0}={Q}_{0}-y\displaystyle \sum _{i\ne 0}\displaystyle \frac{{\gamma }_{0}^{2}{Q}_{i}+{\gamma }_{i}^{2}{Q}_{0}}{{\epsilon }_{0}-{\epsilon }_{i}}+O({y}^{3}).\end{eqnarray} \tag{ 38 }$$

Other ${\tilde{Q}}_{k}$ can be obtained with the replacement $0\to k$. Unlike Q_k, there is no indication of the series truncating at any finite order for ${\tilde{Q}}_{k}$.

The conserved charges constructed above depend on the microscopic parameters of the Hamiltonian, i.e. the hopping and on-site energies. As we shall show later, the same Hamiltonian can have a localized and delocalized phase depending on the values of these parameters. It is thus important to understand if and how the conserved charges themselves differ in the two phases. More precisely, how do the conservation laws 'know' whether a particular choice of microscopic parameters produces a localized or delocalized phase?

The answer has to do with their convergence since they are expressed as power series in the microscopic parameters and particle operators. We thus need to state in what sense the power series are convergent. A reasonable condition for convergence is a sufficiently rapid decay of the coefficients ${\eta }_{{ij}}^{m}$ with increasing m. However, this is complicated by the fact that there are energy difference denominators in the coefficients ${\eta }_{{ij}}^{m}$ that can cause them to blow up when the on-site energies at two different sites are equal. To avoid this, we restrict ourselves to a particular type of disorder that may be termed 'non-resonant'. By this we mean any ensemble of ${\epsilon }_{i}$, which shows 'level repulsion', i.e. the probability of finding ${\epsilon }_{i}$ very close to each other is very small.

From the random matrix theory, we know that the eigenvalues of a generic matrix display level repulsion in their eigenvalues of various degree, the Gaussian orthogonal ensemble (GOE) [50] of real symmetric matrices has the least level repulsion. This condition ensures that perturbative resonances from small denominators, that would otherwise cause individual terms in the expansions of the conserved charges to diverge, are prohibited. This choice is similar to the one involving limited level attraction recently adopted in the context of many-body localization [51].

We have verified that this distribution of onsite energies gives us localization (as indicated from a calculation of the PR) immediately upon switching on the hopping term. Thus, this particular choice of onsite energies, which is of great convenience from the point of view of calculations, is also not unphysical. The on-site energies ${\epsilon }_{i}$ are drawn from the eigenvalues of a real symmetric matrices whose elements are taken from a Gaussian random distribution with fixed variance. The eigenvalues of these matrices are assigned randomly to different sites. Different random assignments then constitute different realizations of disorder, which can then be averaged over to check for convergence. The result of this procedure is shown in figure 3, where ${\epsilon }_{i}$ are drawn from the eigenvalues of real symmetric matrices whose elements are taken from a Gaussian distribution of variance σ = 0.1, 0.25 and 0.4. It can be seen that the ${\eta }^{m}$ decrease rapidly with increasing order of power series m indicating convergence. We have also checked the convergence of the power series for ${\epsilon }_{i}$ drawn from the eigenvalues of non-integrable $t-{t}^{\prime }-V$ model, which also follow a GOE distribution [33, 52].

Zoom In Zoom Out Reset image size

Since ${\eta }_{{ij}}^{m}$ contain more than one term for each m, we checked the convergence of a typical term, which is of the form $\frac{{t}^{m}}{({\epsilon }_{{a}_{1}}-{\epsilon }_{{b}_{1}})({\epsilon }_{{a}_{2}}-{\epsilon }_{{b}_{2}})\ldots ({\epsilon }_{{a}_{m}}-{\epsilon }_{{b}_{m}})}$. Recall that the ${m}^{{\rm{th}}}$ order term in the calculation of ${\tilde{Q}}_{0}$ involves sites with labels between $N-(m-1)$ and m as can be seen from figure 2. Thus, the only values of ${\epsilon }_{i}$ involved are are those chosen from $[{\epsilon }_{N-(m-1)},{\epsilon }_{m}]$ (${\epsilon }_{0}$ is at the center) such that ${\epsilon }_{{a}_{i}}\ne {\epsilon }_{{b}_{i}}$, $\forall i$ and $\mathrm{max}| {a}_{i}-{b}_{i}| =m$.

As the aim of this work is to construct conserved charges in localized systems, it is legitimate to ask whether this slightly non-standard choice of disorder distribution produces localization. We have verified this through numerical exact diagonalization by calculating the PR. We find that the PR for different eigenstates is indeed close to zero for systems of size N = 500 as shown in figure 4, consistent with localization. We thus conclude that our model with on-site energies taken from a GOE distribution does indeed produce a localized phase. A similar exercise to construct the conservation laws for the above model has been carried out in [53]. In that work too, the conserved charges have been constructed as infinite operator series but whose coefficients correspond to the amplitudes of a particle to be on the sites of a square lattice whose sides are the physical one-dimensional lattice. The recursion relation obtained is between conserved charges on different sites and the convergence of the series is assumed to follow from the exponential decay of the eigenfunctions of the Hamiltonian. In our calculations, we construct the conserved charges directly in terms of the microscopic parameters of the Hamiltonian and our convergence criterion is not based on any assumption about the nature of the eigenstates of the Hamiltonian. In fact, as we show in the next section, the convergence of the series for the conserved charges can be used to identify the delocalized and localized phases instead of the eigenfunctions.

Zoom In Zoom Out Reset image size