Wavefunction magnitudes being degenerate everywhere on parameter space even though energy degeneracies occur at isolated points?

Question

Cross-posted here.

Consider the usual simple 2-level gapless graphene Hamiltonian in momentum-space where the energy dispersion is degenerate/gapless at a Dirac point:

\begin{equation}\tag{1} {\small H(k,M)=-t \sum_{\delta} [\cos(k\cdot\delta)\sigma_x-\sin(k\cdot\delta)\sigma_y+M\sigma_z],} \end{equation}

where $t$ is the hopping, $\sigma_i$ are Pauli matrices, $\delta$ lattice vectors and $M$ an on-site potential ($M=0$ for graphene).

When I solve for its wavefunctions both numerically and analytically, the eigenvectors/wavefunctions themselves are complex numbers with an arbitrary complex phase (as expected, depending on the numerical software used). The complex values of the wavefunction at different points on parameter space are different, also as expected.

However, I noticed that the complex magnitude of the eigenvectors are $(\frac{1}{\sqrt2},\frac{1}{\sqrt2})^T$ for all values of $(k_x,k_y)$. By 'complex magnitude', I mean: Each eigenvector $|\psi_i\rangle\equiv z$ is a complex-valued vector with two row elements. Being a complex number, each element may be written as $|z|e^{i\theta}$ for the complex magnitude $|z|$ and complex phase $\theta$. So, in this case, the wavefunctions share the same normalization $(\frac{1}{\sqrt2},\frac{1}{\sqrt2})^T = \frac{1}{\sqrt{2}}(1,1)^T$.

This does not sit well with me because, if I add a small on-site potential/mass to the graphene Hamiltonian to gap the system ($M\neq0$ above; so that the dispersion is non-degenerate everywhere on k-space), these wavefunction magnitudes change depending on $(k_x,k_y)$. So, as expected, there is a higher magnitude at the Dirac points (as in the first picture below). But, the magnitude is constant everywhere when $M=0$ (graphene, second picture), even though I expected it to be a Dirac delta function-like peak at the Dirac points.

I know that Berry-quantities such as the Berry curvature cannot be defined because of the degeneracy, but I want to understand the root and significance of the degenerate wavefunction magnitudes everywhere in parameter space. Besides the common probability interpretation (integration squared = probability = $\frac{1}{2}$ in this case), is there anything else about quantum mechanics / Hilbert spaces / numerical work / physics that this magnitude-observation sheds light on, keeping the issue of gapped vs gapless in perspective? Or is the observed magnitude indicative of the wavefunctions not being well-defined away from the gapless point? What I don't understand is the reason all this happens, and any insight/significance of this.

First picture: Magnitude of one row element of one wavefunction in k-space when $M\neq0$:

Second picture: Magnitude of one row element of one wavefunction in k-space when $M=0$:

Answer 1

This is a somewhat handwaving answer but I think it gets the main point across:

1. We should start by pointing out that the two components of your eigenstates represent the wave function amplitudes on the A/B sublattices of graphene. If you would like to plot the eigenstates in real space, besides a dependency on $(k_x,k_y)$ there will also be a dependence on the lattice vector $\delta$ but let's forget about that for the moment as it is not required for the discussion.
2. In perfect graphene, the two sublattices are symmetric. In the absence of spontaneous symmetry breaking, the electronic density corresponding to your eigenstates (=observable) will therefore be as well. This holds for all energies.
3. As explained in equation (20) of your source, the onsite potential $M$ you're introducing is applied with positive sign on sublattice A and with negative sign on sublattice B, breaking sublattice symmetry. For example, for $M>0$, the lowest-energy states in the spectrum will now have more weight on sublattice B than on sublattice A, while the highest-energy states in the spectrum will have more weight on sublattice A.