Dimensionality in Bistritzer MacDonald Hamiltonian I was trying to understand the origin of flat bands for twisted bilayer graphene and had a basic misunderstanding.
The starting Hamiltonian is $$H = \begin{bmatrix} -iv_0 \sigma_{\theta/2}\nabla && T(r) \\ T^{\dagger}(r) && -iv_0\sigma_{-\theta/2}\nabla \end{bmatrix} $$
where $\sigma_{\theta/2} = e^{-\frac{i\theta}{4}\sigma_z}(\sigma_x,\sigma_y)e^{\frac{i\theta}{4}\sigma_z}$ and $\nabla = (\partial_x, \partial_y)$, while $$ T(r) = \sum_{j=1}^3 T_je^{-i\mathbf{q_j}\mathbf{r}}$$
with $T_{j+1} = w_0\sigma_0+w_1[\cos(\phi j )\sigma_x+\sin(\phi j)\sigma_y]$, for $\phi = 2\pi/3$ and $ q_1 = k_{\theta}(0, −1), q_{2,3} = k_\theta(\pm
\sqrt3/2, 1/2)$, $k_{\theta} = 2k_D\sin(\theta/2)$.
My question is, this is a $4\times4$ Hamiltonian so diagonalizing this would give me 4 bands. How does this Hamiltonian generate the bands in the Moire BZ which would be the bands folded in from the unit cell BZ? How is the Moire Hamiltonian being constructed from this Hamiltonian?
 A: You can have more bands than the size of the matrix. Think of Harper's equation:
$$
\psi_{n+1}+\psi_{n-1} + (\lambda \cos q n)\psi_n = E\psi_n
$$
This is a nearest-neighbour hopping problem with an position dependent on-site potential.  The unit cell is not one site, but instead depends on the periodicity determined by $q$. If you have to go $N$ sites to get back the same $\lambda \cos qn $ potential there will be $N$ bands.
In your case the number of bands  will be four time the number of sites in the Moire unit cell.
A: After I posted my original answer, I realized that there's a better explanation, and it explains the name "continuum model": (for simplicity, I'll omit sublattice degree of freedom)
Suppose that the eigenstate is:
$$
\widetilde{\psi}(k_0)=
\begin{pmatrix}
\widetilde{\Psi}_1(k_0)\\
\widetilde{\Psi}_2(k_0)
\end{pmatrix}
=
\begin{pmatrix}
\sum_{k_1} \widetilde{\psi}_1(k_1)|k_1\rangle\\
\sum_{k_2} \widetilde{\psi}_2(k_2)|k_2\rangle
\end{pmatrix}
$$
1, 2 stand for the first and second layer. $\{k_1\}$ is a set of crystal momentum that live in the first Brillouin zone of the first layer of graphene, and they fold to the the same point $k_0$ in Moire Brillouin zone. The same applies to $\{k_2\}$ as well.
The Hamiltonian is:
$$\hat{H}(k_0)=
\begin{pmatrix}
\hat{H}_{11}(k_0) & \hat{H}_{12}(k_0)\\
\hat{H}_{12}^{\dagger}(k_0) & \hat{H}_{22}(k_0)\\
\end{pmatrix}
\\
\hat{H}_{11}(k_0)=\sum_{k_1}^{BZ_1}\hat{h}_{11}(k_1)|k_1\rangle\langle k_1|
\\
\hat{H}_{22}(k_0)=\sum_{k_2}^{BZ_2}\hat{h}_{22}(k_2)|k_2\rangle\langle k_2|
\\
\hat{H}_{12}(k_0)=\sum_{k_1}^{BZ_1}\sum_{k_2}^{BZ_2}\hat{h}_{12}(k_1, k_2)|k_1\rangle\langle k_2|
$$
Then, by defining $\psi_i(x)=\sum_{k_i}^{BZ_i}\widetilde{\psi}_i(k_i)e^{ik_i x} / \sqrt{N}, (i=1, 2)$, we have:
$$
\begin{align}
\hat{H}_{11}(k_0)\widetilde{\Psi}_1(k_0) & =
\sum_{k_1}^{BZ_1}\hat{h}_{11}(k_1)\widetilde{\psi}_1(k_1)|k_1\rangle
\\
& =\sum_m\sum_{k_1}^{BZ_1}\hat{h}_{11}(k_1)\widetilde{\psi}_1(k_1)\frac{e^{ik_1 m}}{\sqrt{N}}|m\rangle
\\
& =\sum_m \hat{h}_{11}(-i\nabla) \left[\sum_{k_1}^{BZ_1}\widetilde{\psi}_1(k_1)\frac{e^{ik_1 x}}{\sqrt{N}}\right]_{x=m} |m\rangle
\\
& =\sum_m \left[ \hat{h}_{11}(-i\nabla) \psi_1(x) \right]_{x=m} |m\rangle
\end{align}
$$

$$
\begin{align}
\hat{H}_{12}(k_0) \widetilde{\Psi}_2(k_0) & =
\sum_{k_1}^{BZ_1}\sum_{k_2}^{BZ_2} \hat{h}_{12}(k_1, k_2) \widetilde{\psi}_2(k_2) |k_1\rangle
\\
& =\sum_m \sum_{k_1}^{BZ_1}\sum_{k_2}^{BZ_2} \hat{h}_{12}(k_1, k_2) \widetilde{\psi}_2(k_2) \frac{e^{ik_1 m}}{\sqrt{N}} |m\rangle
\\
& \approx\sum_m \sum_{q=k_1-k_2} \hat{h}_{12}(q) e^{iqm} \sum_{k_2}^{BZ_2} \widetilde{\psi}_2(k_2) \frac{e^{i k_2 m}}{\sqrt{N}} |m\rangle
\\
& =\sum_m \left[ \left(\sum_{q=k_1-k_2} \hat{h}_{12}(q) e^{iqx} \right) \psi_2(x) \right]_{x=m} |m\rangle
\\
& =\sum_m \left[ \hat{T}(x) \psi_2(x) \right]_{x=m} |m\rangle
\end{align}
$$
There're similar equations for $\hat{H}_{12}^{\dagger}(k_0)$ and $\hat{H}_{22}(k_0)$. In the derivation above, I've made an approximation $\hat{h}_{12}(k_1, k_2) \approx \hat{h}_{12}(k_1-k_2)=\hat{h}_{12}(q)$. I haven't rigorously analysed this approximation, but I feel that it's fine as long as $k_2$ is close to a specific valley of single graphene layer. Now, the problem of finding eigenstate and eigenvalue for tBLG now becomes solving differential equations:
$$
\hat{h}_{11}(-i\nabla) \psi_1(x) + \hat{T}(x) \psi_2(x) = E \psi_1(x)
\\
\hat{h}_{22}(-i\nabla) \psi_2(x) + \hat{T}^{\dagger}(x) \psi_1(x) = E \psi_2(x)
$$
One thing that might be asked of these differential equations is: where is the dependence of $k_0$, the momentum in Moire Brillouin zone? Note that the sum $\sum_{k_i}^{BZ_i}$ doesn't involve all $k_i$ terms in the first Brillouin zone $BZ_i$, only those that can fold to $k_0$ is considered. As a consequence, $\psi_i(x)$, a solution of the differential equations with an arbitrarily chosen $E$, might not be converted to the form $\sum_{k_i}^{BZ_i}\widetilde{\psi}_i(k_i)e^{ik_i x} / \sqrt{N}$, and thus, might not be our solution. So, given $k_0$, eigenvalue $E$ can only take a discrete set of values, and then we get our band structure.
I'm not sure how this set of differential equtions is actually solved and how eigenvalue $E$ is correctly chosen in practice. Or perhaps I have some misunderstanding in this continuum model? Please discuss with me if you have any idea.


I was having exactly the same problem. The following is my own understanding, please discuss with me if there's any mistake.
I think $T(r)$ is just a symbol for interlayer hopping, and $\mathit{H}$ is just a symbol showing you the necessary ingredients to model tBLG, not the actual Hamiltonian to diagonalize.
Note that $T(r)=\sum_{j=1}^{3}T_{j}e^{-i\pmb{q}_{j}\pmb{r}}$ has the form of Fourier transformation (the sum only has 3 terms due to approxiamtion). $T_{j}$ corresponds to an interlayer hopping process with momentum transfer $\pmb{q}_{j}$, and one can recover $T_{j}$ by doing integral $\frac{1}{\Omega}\int_{\Omega}d\pmb{r}\, T(\pmb{r})e^{i\pmb{q}_j\pmb{r}}$ , where $\Omega$ is the area of Moire superlattice. So, $T(r)$ can indeed be a symbol representing interlayer hopping.
