But I'm very new to this so I could be wrong, please feel free to correct me. I have added the fact that the MPS is PBC and other things for clarity. Thank you!

The boundary here is PBC, which is why I took the trace. Because there's only one tensor for all sites, I think whether I contract it from the left or right leg it would still be normalized. I understood an MPS being left-canonical (satisfying the above condition) mean if I contract the tensors at each site, with new tensors being added to the right of the old ones, then the state would stay normalized. But when I take the inner product, the norm is 1.25.

Thank you for helping me! I think my main confusion is whether the degeneracy is of classical or quantum nature. So I understand that if we have a system evolving to some Hamiltonian, and we let it interact with a reservoir (canonical ensemble), it will thermalize following the Boltzmann distribution. This degeneracy is to me of classical nature. But let's say we prepare an initial state, and let it evolve solely with that Hamiltonian without a reservoir (microcanonical), and then measure to find it in a degenerate level, would the state end up in a quantum superposition? Is this the answer?

But wouldn't that mean the physics of a Hamiltonian is highly sensitive to the initial state, which seems like something that would break the assumption of statistical mechanics?

I'm sorry, but could you clarify it? From what I understand, it's similar to the answer by Lukas below, stating that the wave function collapses to the eigenspace, but it states nothing about the specific distribution of those eigenstates and why is it so (even though this distribution certainly exists, as we can just do another measurement). My question is what is the state of the system in-between the two measurements (let's say 1st energy and 2nd in spin basis), when the system is in that eigenspace; what is its state/distribution in the basis of the second measurement?

I think the problem is that this ratio is totally arbitrary. For the same Hamiltonian, it could be 50/50, 70/30, etc. and according to Wikipedia, it's 50/50. If true, why is that the case and if not, what's the actual state of the system here?

Yes I see, but this doesn't seem to answer the question. Obviously I don't know the state prior to measurement, so I don't know those $\alpha_{\lambda}$. Measuring to $\tilde{E}$ would collapse to the subspace of the degenerate states, so I now know $\alpha_{\lambda} = 0 \ \forall \lambda \neq 1,2$, but still I don't know the distribution for $1,2$ or whether it's 50/50. Am I correct in saying that?

I see, thank you. So if I understand correctly, the free particle dispersion is actually $\epsilon^{\text{free}}_k=\frac{\hbar^2 k^2}{2m}$, and the tight-binding $\epsilon^{\text{tb}}_k=-2t\cos k$ is an approximation of the 1D free fermi gas. In fact, with the correct $\epsilon^{\text{free}}_k$, far-hopping is possible, but we ignore this to get $\epsilon^{\text{tb}}_k$. Am I correct? But why do we do this?

Ah I see, the hopping terms show up again with negative coefficients, and the number operator terms have positive coefficients, meaning the Hamiltonian favors hopping and disfavor sitting. But I tried for 3 sites, and there are non-zero coefficients for far-hopping like $c^\dagger_1 c_3$. It still doesn't reproduce the Hubbard model exactly. Is this physical?