You can't prove it, because it isn't true. Given any 4 by 4 matrices that satisfy the Dirac algebra, you can make 8 by 8 matrices $\gamma_i$ which are equal to the four by four matrices in their upper diagonal, and also equal to the same thing in the lower diagonal. The goal is to look for the lowest dimensional representation, the smallest possible matrices. In this stupid trick, the upper 4 components of the spinor transforms exactly the same as the lower 4 components, meaning that you can reduce the representation by setting the top and bottom components to be equal. You are looking for an irreducible representation, which just means you can't do that.
To understand the size of the lowest dimensional representation Dirac matrices, there is a nice trick described in an article of Scherk's which does this for arbitrary dimensions. In Euclidean signature, make complex dimensions out of pairs of real dimensions
$$ z_1 = x_1 + i x_2 $$
$$ z_2 = x_3 + i x_4 $$
and then linearly combine the Dirac matrices as indicated by this coordinate change:
$$ \gamma'_1 = \gamma_1 + i \gamma_2 $$
$$ \gamma'_2 = \gamma_3 + i \gamma_4 $$
Then the $\gamma$ algebra in terms of the new $\gamma'$ matrices and their conjugates turns into the commuting Fermionic raising and lowering operators. These have a minimal representation which begins by defining the state $|0>$ which is annihilated by all lowering operators, and then acting the raising operators at most once, to produce states. This produces $2^n$ different states, where n is the dimension divided by 2, and this is the basic starting point, ignoring three annoying complications.
These states produced by raising and lowering give you the dimension of the spinor (ignoring the two complications) This starting point tells you that the Dirac algebra should be represented by 4 by 4 matrices in 4d, by 8 by 8 matrices in 6d, by 16 by 16 in 8d, 32 by 32 in 10d. Two to the power of half the dimensions.
The three annoying complications are: odd dimensions, Weyl Fermions, and Majorana Fermions, which each are a separate not so long discussion, but you can see how big the Dirac matrices are supposed to be from the argument above, at least approximately, and you can figure out the useful forms in 2d,3d, and 4d from just piddling around, as Pauli, Dirac, Weyl and Majorana did. The generalization to higher dimensions for string theory is the only time you need to get systematic about it, so it isn't discussed well outside of string theory literature.
That's essentially what we mean by a quasiparticle: a degree of freedom that has been otherwise decoupled from the rest of the system and which behaves exactly like a particle as far as quantum mechanics can tell. Since the (anti)commutation relations are the crucial part of the quantum mechanics of the relevant particles, those need to be carried as well.
Moreover, the canonical anticommutation relations are the best you could hope for anyways. This is because the transformation,
\begin{align}
b_j &= \sum_{k=1}^n \left(U_{jk}c_k + V_{jk}c_k^\dagger\right)\\
b_j^\dagger &= \sum_{k=1}^n \left( U^*_{jk}c_k^\dagger + V_{jk}^*c_k\right)
,
\end{align}
requires the anticommutator to read
\begin{align}
\{b_i^\dagger,b_j\}
&=
\left\{
\sum_{k=1}^n \left( U_{ik}^*c_k^\dagger + V_{ik}^*c_k\right),
\sum_{l=1}^n \left(U_{jl}c_l + V_{jl}c_l^\dagger\right)
\right\}
\\&=
\sum_{k=1}^n \sum_{l=1}^n \left[
U_{ik}^*U_{jl}
\left\{c_k^\dagger,c_l \right\}
+
V_{ik}^*V_{jl}
\left\{c_k,c_l^\dagger\right\}
\right]
\\&=
\sum_{k=1}^n \sum_{l=1}^n \left[
-U_{ik}^*U_{jl}
+
V_{ik}^*V_{jl}
\right]
\delta_{kl}
\\&=
\sum_{k=1}^n\left[ V_{ik}^*V_{jk} - U_{ik}^*U_{jk} \right]
\end{align}
in the most general case. Thus, in the worst-case scenario, you've bungled up the (anti)commutation relations, and you're going to need a nonzero (anti)commutator any time you need to slide a $b_j$ past a $b_i^\dagger$, even if $i\neq j$. Thus, you require at the very least that the anticommutator be diagonal.
Moreover, since $\{b,b^\dagger\}=bb^\dagger + b^\dagger b$ is the sum of two positive semidefinite hermitian operators, the diagonal anticommutator needs to be real and non-negative, so you've got two possibilities:
- it can be positive, in which case it can be rescaled to one, or
- it can be zero, in which case you're essentially losing information, probably via a degenerate matrix, and the algebra spanned by the $b_i,b_i^\dagger$ will be unable to fully span all of the $c_i,c_i^\dagger$.
More broadly speaking, you're obliged to "conserve the quantumness" of the problem (unless you explicitly want to make mean-field approximations or similar) and this requires you to maintain a full set of generators for the operator algebra and therefore to have nontrivial anticommutators between your new quasiparticle fermionic operators.
Best Answer
The matrix is just a discretized version of $-\partial_x^2$. Define a vector ${\bf x}=x_n$ where $n$ in the site, then in ${\bf y}=H{\bf x}$ the vector $y$ has components $$ y_n=(x_{n+1}- x_n) - (x_n-x_{n-1})= x_{n+1}-2x_n + x_{n-1}, $$ which being the difference of two differences approximates the second derivative.