Fermionic occupation for an inhomogeneous tight-binding model

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.

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!