## 3 Direct Numerical Methods for Strongly Disordered Systems
### Transfer Matrix Method  
The transfer matrix method (TMM) is used extensively on various disorder problems [39; 40; 41]. Unlike brute force diagonalization methods, the TMM can handle rather large system sizes. When combined with finite-size scaling, this method is very robust for detecting the localization transition and its corresponding exponents. Most of the accurate estimates of critical disorder and correlation length exponents for disorder models in the literature are based on this method [40; 41].  
The simplifying assumption of the TMM is that the system can be decomposed into many slices (Figure 10), and each slice only connects to its adjacent slice. Precisely for this reason, the TMM is not ideal for models with long-range hopping, or long-range disorder potentials or interactions.  
We can understand the computational scaling of the TMM by a simple 3D example without an explicit interaction. We assume the system has a width and height equal to \(M\) for each slice of a \(N\)-slice cuboid, forming a \"bar\" of length \(N\). The Hamiltonian can be decomposed into the form  
\[H=\sum_{i}H_{i}+\sum_{i}(H_{i,i+1}+H.c.), \tag{8}\]  
where \(H_{i}\) describes the Hamiltonian for slice \(i\) and \(H_{i,i+1}\) contains the coupling terms between the \(i\) and \(i+1\) slices. The Schrodinger equation can be written as  
\[H_{n,n+1}\psi_{n+1}=(E-H_{n})\psi_{n}-H_{n,n-1}\psi_{n-1}\,, \tag{9}\]  
Figure 10: Schematic of a transfer matrix method (TMM) calculation. Assuming the system has a width and height equal to \(M\) for each slice of a \(N\)-slice cuboid, forming a “bar” of length \(N\), the amplitude of the wavefunction in the 0-th slice can be related to that in the N-th slice via the transfer matrix, Equation (10).  
where \(\psi_{i}\) is a vector with \(M^{2}\) components which represent the wavefunction of the slice \(i\). This may be reinterpreted as an iterative equation  
\[\left[\begin{array}{c}\psi_{i+1}\\ \psi_{i}\end{array}\right]=T_{i}\times\left[\begin{array}{c}\psi_{i}\\ \psi_{i-1}\end{array}\right]. \tag{10}\]  
where the transfer matrix  
\[T_{i}=\begin{bmatrix}H_{i,i+1}^{-1}(E-H_{i})&-H_{i,i+1}^{-1}H_{i,i-1}\\ 1&0\end{bmatrix}. \tag{11}\]  
The goal of the TMM is to calculate the localization length, \(\lambda_{M}(E)\) for a system with linear size \(M\) at energy \(E\), from the product of \(N\) transfer matrices  
\[\tau_{N}\equiv\prod_{i=1}^{N}T_{i}. \tag{12}\]  
The Lyapunov exponents, \(\alpha\), of the matrix \(\tau_{N}\) is given by the logarithm of its eigenvalues, \(Y\), at the limit of \(N\rightarrow\infty\), \(\alpha=\lim_{N\rightarrow\infty}\frac{\ln(Y)}{N}\). The smallest exponent corresponds to the slowest exponential decay of the wavefunction and thus can be identified as corresponding to the localization length, \(\lambda_{M}(E)=1/\alpha_{\min}\)[128; 129; 130; 131; 132; 133; 134].  
Since the repeated multiplication of \(T_{i}\) is numerically unstable, periodic reorthogonalization is needed in the numerical implementation [39; 40; 41]. For the 3D Anderson model, the reorthogonalization is done for about every 10 multiplications. This is the major bottleneck for the TMM method, as reorthogonalization scales as the third power of the matrix size. Therefore, the method in general scales as \(M^{3}\).