After days of thinking, searching, discussions, and testing, I can finally answer my own question now. The answer is much more involved than I expected from such a "simple" XY model (even just for the Ising model)!
All "correct solutions/spectrum" stated below are checked against results from exact diagonalization.
Simply put, there's nothing wrong with choosing either sign for the single fermion energy $\varepsilon(k)$, and it doesn't matter (nor help) whether we stick to a single branch of $\theta_k = \arctan \xi$. The crucial thing to take into account (and make things right) is that, after making the Bogoliubov transformation, the $b$-fermion vacuum state $|0\rangle_b$ is different from the $c$-fermion vacuum state $|0\rangle_c$, and may have different parity $P_0$ and momentum $k_0$!
To further illustrate my statement above, let me first clarify and fix my terminologies and notations:
- $c_j$ are the fermion operators after the Jordan-Wigner transformation, $c_k$ are the fermion operators after the Fourier transform, $b_k$ are the fermion operators after the Bogoliubov transformation.
- $a$-fermion vacuum state refers to the zero particle state $|0\rangle_a$ annihilated by fermion operators $a$'s. This is NOT necessarily the ground state of the system.
- Parity always refers to the operator $P \equiv \prod_j (-\sigma_j^z) = (-1)^{\sum_j c_j^\dagger c_j}$ or its corresponding eigenvalues. Note that it is defined in the spin basis, and has a simple representation in the $c$-fermion basis, but not necessarily in the $b$-fermion basis.
- $P_0$ and $k_0$ refer to the parity and the momentum of a vacuum state.
As stated in the question, after the Jordan-Wigner transformation, the original Hilbert space of the periodic spin chain is mapped into two sectors of Fock space with $P=\pm 1$ (even/odd) respectively. The $c$-fermion vacuum states $|0\rangle_c$ in both sectors have $P_0=1$ and $k_0=0$. The Fourier transform respects parity and particle number conservation, and does not change the vacuum state. If the Hamiltonian is already diagonal in this Fourier basis (i.e. the isotropic case $\eta=0$), we can directly build up the entire spectrum (including dispersion relation) using $n$-particle state $c_{k_1}^\dagger c_{k_2}^\dagger \dots c_{k_n}^\dagger |0\rangle_c$ whose parity and momentum are given by $P = (-1)^n$ and $k = k_1 + k_2 + \dots + k_n$. With appropriate choice of $n$ (odd/even) and $k$ (integers/half-integers) in the two sectors, there's no more subtleties.
Demonstration with the $\eta=0, \theta_k = \pi$ case
Things complicate with the Bogoliubov transformation, because $|0\rangle_b \neq |0\rangle_c$, and the latter state has different $P_0$ and $k_0$, and may not even be the same in the two sectors! For demonstration purpose, let's still take the $\eta=0$ isotropic case, and make a nontrivial Bogoliubov transformation with $\theta_k = \pi$:
$$ c_k = i b_{-k}^\dagger, \quad c_k^\dagger = -i b_{-k}. $$
Nothing is wrong with this transformation. We just need to be careful about the new vacuum state $|0\rangle_b$. In this case, there's a simple relation between the two vacuum states:
$$ |0\rangle_b = \prod_k c_k^\dagger |0\rangle_c. $$
This means that, in the $b$-fermion basis:
- If total number of sites $N$ is even, $P_0=1$ in both sectors, $k_0 = 0$ or $\pi/2$ in the parity even/odd sector respectively.
- If total number of sites $N$ is odd, $P_0=-1$ in both sectors, $k_0 = \pi/2$ or $0$ in the parity even/odd sector respectively.
The parity and momentum of a state with $n$ $b$-fermions is given by $P = (-1)^n P_0$ and $k = k_0 + k_1 + \dots + k_n$. Using these formulas and choosing the appropriate $P$ and $k$'s in the two sectors, I can still build up the correct spectrum and dispersion relation with $\theta_k = \pi$.
Discussion about the general case
The principle is the same in the general case ($\eta \neq 0$). All we need to do is to determine $P_0$ and $k_0$ of $|0\rangle_b$ in the two sectors, and use $P = (-1)^n P_0$ and $k = k_0 + k_1 + \dots + k_n$ to determine the parity and momentum of a $n$ $b$-fermion state. However, the determination of $P_0$ and $k_0$ becomes very nontrivial, because there is no simple relation between $|0\rangle_b$ and $|0\rangle_c$ when $\eta \neq 0$ (however, see reference at the end).
As correctly calculated in @cesaruliana's answer, we have
$$ \tan\theta_k = \frac{ \eta \sin(ka)}{ g - \cos(ka)}, $$
all branches of $\theta_k = \arctan\xi$ give valid transformations, but it is not obvious which one(s) is the most helpful. Our goal is to find $P_0$ and $k_0$ of the ground state. One strategy we can use is choosing
$$ \cos \theta_k = \frac{ g - \cos(ka) }{ \sqrt{ (g-\cos(ka))^2 + \eta^2 \sin^2(ka)} },
\quad
\sin \theta_k = \frac{ \eta \sin(ka) }{ \sqrt{ (g-\cos(ka))^2 + \eta^2 \sin^2(ka)} }, $$
hence making all the single $b$-fermion energies positive:
$$ \varepsilon(k) = +2J \sqrt{ (g-\cos(ka))^2 + \eta^2 \sin^2(ka)}. $$
Now the vacuum state $|0\rangle_b$ is the ground state of the $b$-fermion system, making the ground state (of either the fermion or the spin system) much easier to track. Next, we analyze the ground state of the original spin model in different limits.
For example, in the Ising weak field limit
$ H = -J \sum \sigma_j^z \sigma_{j+1}^z, $
there are two degenerate ground states with definite parities:$ |\pm\rangle = \prod_j \left| \leftarrow \right>_j \pm \prod_j \left| \rightarrow \right>_j$, $ P |\pm \rangle = \pm | \pm \rangle$ and both have momentum $k_0=0$. We note that this degeneracy is the characteristic feature of a symmetry breaking phase. So, until a phase transition occurs (at $g=1$ or $\eta=1$), the two-fold degeneracy should always be present and exact, and the quantum number $P$ and $k$ of the two ground states should be unchanged. To have such ("robust") two-fold degeneracy for all values of $0\leq g\leq 1$ and $0\leq \eta \leq 1$, the two ground states must be the two vacuum states $|0\rangle_b$ of the two sectors. This means that, for $0 \leq g\leq 1$, $0 \leq \eta \leq 1$, $P_0 = + 1$ in the parity even sector and $P_0 = -1$ in the parity odd sector, and $k_0=0$ for both sectors.
Therefore, in both sectors, only states with even number of $b$-fermions should be included in the final spectrum.
The correct spectrum is indeed built up from this construction.
$P_0$ and $k_0$ may change after a phase transition. Similar analysis must be done again in other phases.
The mistake I was making was that I calculated parity and momentum naively using $b$-fermions, always taking states with even/odd number of $b$-fermions in the parity even/odd sector, essentially always assuming $P_0=1, k_0=0$. This is NOT true in general.
Acknowledgements and references
I would like to thank my fellow student Wen Wei for very helpful discussion. He first found out that the correct solution can be constructed if we take even number of particles in both sectors in the symmetry breaking phase. A very helpful note written by Prof. Andreas Schadschneider and Prof. Götz S. Uhrig that ultimately resolved the confusion can be found here:
Strongly Correlated Systems in Solid State Physics
Update: After I posted this answer, I found another reference (master thesis) by Erik Eriksson where he rigorously calculated in detail the relation between $|0\rangle_b$ ("Bogoliubov vacuum") and $|0\rangle_c$ in the case of Ising model $\eta=1$. I believe this can be easily generalized to the case of $\eta < 1$. His thesis can be found here:
Quantum Phase Transitions in an Integrable One-Dimensional Spin Model