Does particle-hole symmetry always imply half-filling and real correlations $\langle c^\dagger_n c_{n+1} \rangle$? Suppose we had a lattice Hamiltonian $H$ which was symmetric under the particle-hole transformation
$$ c_n \mapsto U^\dagger c_nU=(-1)^nc^\dagger _n$$
such that $[H,U] = 0$, where $c_n$ are Fermionic operators obeying $\{ c_n , c_m^\dagger \} = \delta_{nm}$ and $\{ c_n,c_m \} = \{ c^\dagger_n, c^\dagger_m \} = 0$ where the indices label the lattice sites of the system.
If the ground state $|\Omega \rangle$ of $H$ is unique, then it would be an eigenstate of the particle-hole transformation as $U |\Omega\rangle = e^{i \theta} |\Omega\rangle$, where the eigenvalue is some phase as $U$ is unitary. From this, if I calculated the nearest-neighbour correlation, I would find
\begin{align*}
\langle \Omega | c_n^\dagger c_{n+1} |\Omega \rangle & = \langle \Omega | U^\dagger c_n^\dagger c_{n+1} U |\Omega \rangle \\
& = \langle \Omega |( U^\dagger c_n^\dagger U)(U^\dagger c_{n+1} U) |\Omega \rangle \\
&= (-1)^{2n+1}\langle \Omega| c_n c_{n+1}^\dagger |\Omega \rangle \\
&=- \langle \Omega |(-c_{n+1}^\dagger c_n) |\Omega \rangle \\
&=  \langle \Omega| c^\dagger_{n+1} c_n |\Omega \rangle \\
&= \langle \Omega | c_n^\dagger c_{n+1} |\Omega \rangle^*
\end{align*}
therefore I find it is real. A similar calculation would also show that this implies half-filling $\langle \Omega| c^\dagger_n c_n |\Omega \rangle = \frac{1}{2}$ as shown in this answer.
This calculation relied upon the fact that $|\Omega\rangle $ is unique, however many systems will have ground state degeneracy in which case $U |\Omega \rangle \neq e^{i\theta} |\Omega\rangle$ in general. One can always find a state in the ground state subspace that is an eigenstate however, and the above results will apply.
(Edit) Example Hamiltonian
An example of an interacting Hamiltonian that has particle-hole symmetry is given by
$$ H = -t \sum_n c_n^\dagger c_{n+1} + \lambda \sum_n c^\dagger_n c_{n+1} c_{n+2}^\dagger c_{n+3} + \mathrm{h.c.} = H_0 + H_\mathrm{int} $$
where $t,\lambda \in \mathbb{R}$. We have
\begin{align*}
U^\dagger H_0 U & = -t \sum_n (-1)^{2n+1} c_n c^\dagger_{n+1} + \mathrm{h.c.} \\
& = t \sum_n c_n c^\dagger_{n+1} + \mathrm{h.c.} \\
& = -t\sum_n c_{n+1}^\dagger c_n + \mathrm{h.c} \\
& = H_0 
\end{align*}
And for the interaction part
\begin{align*}
U^\dagger H_\mathrm{int} U & = \lambda \sum_n (-1)^{4n + 6} c_n c^\dagger_{n+1} c_{n+2} c^\dagger_{n+3} + \mathrm{h.c.}\\
& = \lambda \sum_n c_{n+3} c^\dagger_{n+2} c_{n+1} c_n^\dagger + \mathrm{h.c.} \\
& = H_\mathrm{int}
\end{align*}
where I have interchanged parts with the hermitian conjugate and anti-commuted fermions past each other.
As this Hamiltonian is interacting, I am not sure how I would derive the correlation matrix of this, but all I would like to know is if this system has half-filling and real nearest-neighbour correlations in the ground state.
My question
In the case of no degneracy we can say that particle-hole symmetry implies half filling $\langle c^\dagger_n c_n \rangle = \frac{1}{2}$ and the correlation $\langle c_n^\dagger c_{n+1} \rangle \in \mathbb{R}$. However, in the case of degeneracy the above derivation would not apply unless we chose a state in the ground state subspace that is an eigenstate of $U$. My question is, am I allowed to choose the particle-hole symmetric ground state over another as "special" and can I always assume that particle-hole symmetry implies half-filling and real nearest-neighbour correlations?
 A: I think this is quite a standard question in condensed matter. Particle-hole symmetry always implies half-filling, but whether it implies $\langle c^{\dagger}_{n+1}c_n \rangle$ to be real depends on what your particle-hole symmetry really is. If the transformation is what you defined in the question, the answer is yes.
Consider the partition function $e^{-\beta \hat{H}}$ where $\beta>0$. The correlation function you care is
$$\langle c^{\dagger}_{n+1}c_n \rangle = {\text{Tr}\big(c^{\dagger}_{n+1}c_n e^{-\beta \hat{H}}\big) \over \text{Tr}\big(e^{-\beta \hat{H}}\big)}$$
Note that when $\beta \rightarrow +\infty$, we have the ground state limit given whether or not the ground state is non-degenerate. Since the Hamiltonian $\hat{H}$ is invariant under particle-hole transformation, we have
\begin{align}
\text{Tr}\big(c^{\dagger}_{n+1}c_n e^{-\beta \hat{H}}\big)^{\dagger} & = \text{Tr}\big(e^{-\beta \hat{H}}c^{\dagger}_nc_{n+1}\big) \\
& = \text{Tr}\big(c^{\dagger}_nc_{n+1} e^{-\beta \hat{H}}\big) \\
& = \text{Tr}\big(\hat{U}^{\dagger}c^{\dagger}_nc_{n+1} e^{-\beta \hat{H}}\hat{U}\big) \\
& = \text{Tr}\big(\big(\hat{U}^{\dagger}c^{\dagger}_n\hat{U}\big)\big(\hat{U}^{\dagger}c_{n+1} \hat{U}\big)\big(e^{-\beta \hat{U}^{\dagger}\hat{H}\hat{U}}\big)\big) \\
& = -\text{Tr}\big(c_nc^{\dagger}_{n+1} e^{-\beta \hat{H}}\big) \\
& = \text{Tr}\big(c^{\dagger}_{n+1}c_n e^{-\beta \hat{H}}\big)
\end{align}
This implies $\langle c^{\dagger}_{n+1}c_n \rangle$ is real. You can also use the same approach to show $\langle c^{\dagger}_nc_n \rangle={1 \over 2}$ and $\text{Re}\big(\langle c^{\dagger}_{n+2}c_n \rangle\big)=0$.
So it is a good time to answer your question. Whether we have $\text{Tr}\big(c^{\dagger}_{n+1}c_n e^{-\beta \hat{H}}\big)^{\dagger}=\text{Tr}\big(c^{\dagger}_{n+1}c_n e^{-\beta \hat{H}}\big)$ depends on the form of your particle-hole transformation. If you have $c_n$ transformed by $\hat{U}^{\dagger}c_n\hat{U}=c^{\dagger}_n$, we can no longer show $\langle c^{\dagger}_{n+1}c_n \rangle$ is real ($\text{Re}\big(\langle c^{\dagger}_{n+1}c_n \rangle\big)$ is in fact $0$). However, we still can prove $\langle c^{\dagger}_nc_n \rangle={1 \over 2}$ and $\text{Re}\big(\langle c^{\dagger}_{n+2}c_n \rangle\big)=0$. The very particular case is $\langle c^{\dagger}_nc_n \rangle$. Since the phase given by the particle-hole transformation in $c^{\dagger}_nc_n$ will be cancelled out, no matter which transformation is given, we keep the same result $\langle c^{\dagger}_nc_n \rangle={1 \over 2}$. However, the particle-hole transformation is defined in your question is useful when the lattice of the Hamiltonian is bipartite, including the hexagonal lattice and the square lattice.
A: While I realise your question is in the second quantized picture and this may not be the answer you're looking for, I just wanted to offer a simple argument in the first quantized picture that would demonstrate the symmetry of PHS Hamiltonians about $E = 0$ (half filling).
Define $\hat{\mathcal{P}}$ as the particle-hole operator in the second quantized picture. The action of $\hat{\mathcal{P}}$ is defined as
\begin{equation}
    \hat{\mathcal{P}}\hat{c}^{}_{A}\hat{\mathcal{P}}^{-1} = \sum_{B} \hat{c}_{B}^{\dagger}(U_{P})_{B, A}
\end{equation}
for some unitary matrix $U_{P}$ (i.e. we turn an electron into a hole and vice-versa). $\mathcal{P}$ is also linear in the sense $\mathcal{P}i\mathcal{P}^{-1} = i$ The action of $\mathcal{P}$ on the Hamiltonian, $H$, is therefore
\begin{equation}
    \begin{split}
        \hat{\mathcal{P}}\hat{H}\hat{\mathcal{P}}^{-1} &= \sum_{A, B} \hat{\mathcal{P}} \hat{c}_{A}^{\dagger} \hat{\mathcal{P}}^{-1} \hat{\mathcal{P}} H^{}_{A, B} \hat{\mathcal{P}}^{-1} \hat{\mathcal{P}} \hat{c}^{}_{B} \hat{\mathcal{P}}^{-1}\\
        & =\sum_{A, B} \sum_{C, D} (U_{P}^{})^{\dagger}_{A, C} \hat{c}_{C} H_{A, B} \hat{c}_{D}^{\dagger} (U_{P}^{})^{}_{D, B} \\
        &= \sum_{A, B}\sum_{C, D} \delta^{}_{C, D}(U_{P})^{\dagger}_{A, C}H^{}_{A, B}(U_{P})^{}_{D, B} - \hat{c}^{\dagger}_{D}(U_{P})^{}_{D, B}H^{}_{A, B}(U_{P})^{\dagger}_{A, C}\hat{c}^{}_{C} \\
        & = \sum_{A, B}\sum_{C, D} (U_{P})^{\dagger}_{A, C}H^{}_{A, B}(U_{P})^{}_{C, B} - \hat{c}^{\dagger}_{D}(U_{P})^{}_{D, B}H^{T}_{B, A}(U_{P})^{\dagger}_{A, C}\hat{c}^{}_{C} \\
        & = \sum_{A, B}\left( \delta^{}_{A, B}H_{A, B} - \sum_{C, D}\hat{c}^{\dagger}_{D}(U_{P})^{}_{D, B}H^{T}_{B, A}(U_{P})^{\dagger}_{A, C}\hat{c}^{}_{C} \right) \\
        & = \text{Tr}(H) - \sum_{C, D} \hat{c}^{\dagger}_{D}H^{}_{D, C}\hat{c}^{}_{C},
    \end{split}
\end{equation}
where $H$ is the first quantized Hamiltonian. For the Hamiltonian to respect PHS we require $\hat{\mathcal{P}}\hat{H}\hat{\mathcal{P}}^{-1} = \hat{H}$ and therefore
\begin{equation}
    \begin{split}
        &\text{Tr}(H) = 0 \\
        U_{P}H^{*}&U_{P}^{\dagger} = -H.
    \end{split}
\end{equation}
These conditions are complementary in the sense that the second condition implies the spectrum is symmetric about $E=0$ (try it!) which implies the first. Indeed, diagonalizing the Hamiltonian, $H = \Lambda D \Lambda^{-1}$, we find $\text{Tr}(H) = \text{Tr}(\Lambda D \Lambda^{-1}) = \text{Tr}(D) = 0$. Therefore, our system has $n$ filled and $n$ empty bands and is at half-filling.
