Consider a nonrelativistic free fermion Hamiltonian given by $$H = \sum_{i} (\epsilon_i - \mu) \psi^\dagger_i \psi_i \ ,$$ where we have explicitly included the chemical potential $\mu$. The ground state is given by the "Fermi surface" $$|\text{GS}\rangle= \prod_{i|\epsilon_i < \mu} \psi^\dagger_i|0\rangle \ ,$$ and it is with respect to this nonempty ground state that correlation functions, e.g. the propagator $$\left<\text{GS} | T \psi(r,t) \psi^\dagger(r',t')|\text{GS}\right>$$ are computed.
Why do we never consider the same thing in relativistic QFT?
For example, a real scalar field Hamiltonian is given by $$H_0 = \int d^d r \left(\frac{m^2}{2}\phi^2 + \frac{1}{2}\pi^2 + \frac{1}{2}(\nabla\phi)^2\right)$$ with $\pi = \partial_t\phi$.
This is diagonalized by setting $$\phi \propto \sum_k a_k e^{i(kx-\omega t)} + a^\dagger_k e^{-i(kx-\omega t)}$$ with $\omega=\sqrt{k^2+m^2}$, and the ground state is given by the vacuum $|\text{GS}\rangle = |0\rangle$.
What if we explicitly add a chemical potential, i.e. $$H = H_0 - \sum_k \mu a^\dagger_k a_k \ .$$ I suppose this breaks Lorentz invariance? If so, how about instead considering a complex scalar field, and setting $$H=H_0 - \sum_k \mu (a_k^\dagger a_k - b_k^\dagger b_k) \ .$$ Is this Lorentz invariant?