The model
Consider the simple one-dimensional fermionic tight-binding chain of $N$ sites with inhomogenous hopping couplings $t_n$:
$$ H = – \sum_n t_n c^\dagger_{n+1} c_n + \text{h.c.} \equiv \sum_{n,m} c^\dagger_n h_{nm} c_m$$
where $h_{nm} = -t_m \delta_{n,m+1} + \text{h.c.}$. Examples for $t_n$ could be something strange like $t_n = 1 + 0.1n \sin(2n \pi/N)$.
This system does not have translational symmetry due to the inhomogenous hoppings $t_n$, so cannot be diagonalised with a Fourier transform. The single-particle Hamiltonian $h$ can be diagonalised numerically using a unitary transformation $h = UDU^\dagger$ to give
$$ H = \sum_k E(k) c^\dagger_k c_k $$
where $c^\dagger_k = \sum_n U_{nk} c_n^\dagger$ etc. We define the ground state as the state with all negative energy states filled:
$$ |\Omega \rangle = \prod_{k: E(k) < 0} c^\dagger_k |0\rangle $$
where the notation "$k: E(k) < 0$" means all modes $k$ whose energy $E(k)$ is negative and $|0\rangle $ is the vacuum annihilated by all $c_k$.
Fermionic occupation
Suppose I want to evaluate the fermionic occupation $n_i = c^\dagger_i c_i$ on the ground state $|\Omega\rangle$. If we define the correlation matrix $C_{ij} = \langle \Omega | c^\dagger_i c_j |\Omega \rangle$, then the occupations are simply the diagonal elements of the correlation matrix. The correlation matrix is given by
$$ C_{ij} = \sum_{k:E(k) < 0} U^*_{ik} U_{jk} $$
(see Eq. (5) of this paper).
My question
What I find, quite surprisingly, is that even with inhomogeneity in the system, the quantity $\langle n_i \rangle \equiv \langle \Omega | n_i | \Omega \rangle$ is completely uniform throughout the system, taking the value of $\langle n_i\rangle = 1/2$. Assuming that my Python code is correct, this seems strange to me. As translational symmetry has been broken, I see no reason for $n_i$ to be translationally symmetric in the ground state. Why is this occuring?
I can update my question with my Python code if requested, however I appreciate the computational physics tag is not for answering programming or debugging questions.
Best Answer
You can perform a particle-hole transformation $c_i = (-1)^i d^\dagger_i$ where $(-1)^i$ alternates between $+1$ and $-1$. This transformation preserves fermionic statistics, and your hopping term becomes $$c_{i+1}^\dagger c_i + h.c. = (-1) d_{i+1} d_i^\dagger + h.c. = d_i^\dagger d_{i+1} + d_{i+1}^\dagger d_i.$$
This implies that your Hamiltonian is particle-hole symmetric and the only ground state configuration that satisfies particle-hole symmetry will have $\langle n_i \rangle = 1/2$ on each site.
Next, you could explore the wave function overlap $\langle c_{i+1}^\dagger c_i \rangle$ which should exhibit inhomogeneous behaviour … Or add a (inhomogeneous?) chemical potential to destroy PHS!