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
This is a very good question. The same operator algebra does not imply that $H(J,h)$ and $H(h,J)$ have the same spectrum. As has been mentioned in Dominic's answer, even the ground state degeneracy is different under the interchange of $J$ and $h$ ($J\gg h$: symmetry-broken two-fold degeneracy, and $J\ll h$ unique ground state), therefore it is impossible to establish a one-to-one mapping between the eigenstates of $H(J,h)$ and $H(h,J)$. One must keep in mind that the duality transformation only preserves local dynamics but the global (topological) properties will be lost. This statement becomes sharper in higher dimensions. Like in 2D, the $\mathbb{Z}_2$ lattice gauge theory is dual to the quantum Ising model, however the topological order of the gauge theory (most prominently the topology-dependent ground state degeneracy) is completely lost in the dual Ising model, even though there is a beautiful correspondence between their local excitations (e.g. charge and vison).
Nevertheless, duality is still very useful if we only focus on the local excitations. Many important problems like low-energy dynamics, phase transitions and criticality are only related to local excitations, then the duality transformation can help us a lot in understanding these things.
To appreciate the duality in the 1D transverse field Ising model, it is better to just look locally and try to figure out the local correspondence of the bulk excitations, without caring too much about the global properties such as boundary conditions, infinite strings, ground state degeneracy etc. The idea of duality is actually simple: one can describe a 1D Ising chain either by the spin variables living on each site, or by the kink variables living on each link. A kink on the spin chain is a link across which the Ising spins are opposite. Then every link $l$ can only have two possible states:
$$\tau_l^z=\left\{\begin{array}{cc}+1 & \text{unkinked,}\\ -1 & \text{kinked,}\end{array}\right.$$
If we have specified all the kink configuration $\tau_l^z$ on each link $l$, we can actually determine the spin configuration $\sigma_i^z$ on each site $i$, with only one additional piece of knowledge about the left-most spin $\sigma_0^z$. The trick is just to accumulate the kink configurations from left all the way to the right,
$$\sigma_i^z=\sigma_0^z\prod_{0<l<i}\tau_l^z.$$
So the spin configuration is uniquely determined by the kink configuration (up to the left-most spin). If we count the number of states in the Hilbert space, the Hilbert space dimension will be $2^{N_\text{site}}$ in the spin language and $2^{N_\text{link}}$ in the kink language, where $N_\text{site}$ and $N_\text{link}$ are the number of sites and links respectively, which are equal (apart from the left-most site) on the 1D lattice, so the Hilbert space dimension is actually the same in both languages. In this sense, we can say that the correspondence between the spin and the kink descriptions is almost one-to-one (especially in the thermodynamic limit), even though there might be some complication arising from the left-most boundary (which however is to be ignored in the duality transform, as duality only cares local properties).
Now we can rephrase the original transverse field Ising model
$$H=-J\sum_{i}\sigma_i^z\sigma_{i+1}^z-h\sum_{i}\sigma_i^x,$$
in the kink language. It is not hard to see that the coupling between the Ising spins is just the chemical potential for the kink
$$-J\sum_{i}\sigma_i^z\sigma_{i+1}^z=-J\sum_l\tau_l^z,$$
which basically follows from the physical meaning of the kink variable $\tau_l^z$. Admittedly this equality may run into a little trouble on the boundary, where some sites or links might be missing, but they do not affect the local properties deep in the bulk, so we just ignore. The translation of transverse field term is more involved. In the spin language $\sigma_i^x$ operator just flips the spin on site $i$, which would correspond to simultaneous creation or annihilation two kinks on the links adjacent to that site, or moving an existed kink across the site.
In either case, flipping a spin would correspond to simultaneously changing the kink variables $\tau^z$ on adjacent links, which can be carried out by $\tau^x_l\tau^x_{l+1}$, s.t.
$$-h\sum_i\sigma_i^x=-h\sum_{l}\tau^x_l\tau^x_{l+1}.$$
Obviously the relation between $\tau^x$ and $\tau^z$ is exactly the same as our familiar $2\times 2$ Pauli matrices $\sigma^x$ and $\sigma^z$, s.t. $\tau^x|\tau^z=+1\rangle=|\tau^z=-1\rangle$ and $\tau^x|\tau^z=-1\rangle=|\tau^z=+1\rangle$, and the algebraic relations like $\{\tau^x,\tau^z\}=0$ are just consequences that follow. As the OP have pointed out, the algebraic relations can not guarantee the representation to be fundamental, it is actually the above physical picture that guarantees the representation of the kink operators.
Putting together the above results, we arrive at the Hamiltonian in terms of the kink operators $\tau_l^x$ and $\tau_l^z$
$$H=-h\sum_{l}\tau^x_l\tau^x_{l+1}-J\sum_l\tau_l^z.$$
One might have been wondering for a while that why I kept using the symbol $\tau$ instead of $\mu$ in the original post. Now it is clear that $\tau$ is still one step away from $\mu$ by a basis transformation on each link that redefines $\tau^x=\mu^z$ and $\tau^z=\mu^x$ (relabeling $x\leftrightarrow z$). Such a unitary transform will not change any physics, but just to bring back the standard form of the transverse field Ising model to accomplish the duality transform,
$$H=-h\sum_{l}\mu^z_l\mu^z_{l+1}-J\sum_l\mu_l^x.$$
So the duality between $J$ and $h$ is now manifest, but what is the physical meaning of $\mu_l^z$ indeed? To answer this question, one should first understand that the relation between $\tau^z$ and $\tau^x$ is just like that between the coordinate and the momentum. They are related by a $\mathbb{Z}_2$ version of the Fourier transformation. If we treat the two states of $|\tau^z=\pm 1\rangle$ as two position eigenstates in a two-site system, then the $\tau^x$ eigenstates $|\tau^z=+1\rangle\pm|\tau^z=-1\rangle$ are nothing but the momentum eigenstates with the momentum = $0$ and $\pi$ respectively. In this sense, we can say $\mu_l^z\equiv\tau_l^x$ is the conjugate momentum of the kink variable $\tau_l^z$ on each link. In fact, this concept is so important that people invent a name for $\mu^z$, i.e. the vison variable, s.t.
$$\mu_l^z=\left\{\begin{array}{cc}+1 & \text{vison off,}\\ -1 & \text{vison on.}\end{array}\right.$$
By saying that the vison is the conjugate momentum of the kink, we mean that if there is a vison sitting on a link then the kinked and the unkinked configurations across that link will be differed by a minus sign in the wave function. Unlike the kink which is just another way to encode the spin configuration, the vison does not have a correspondent spin configuration. In fact, the vison configuration is encoded in the relative sign between different spin configurations in the wave function. It represents the inter-relation among the spin configurations other than any particular spin configuration itself, or in other words, the quantum entanglement in the spin wave function.
Mathematically this can be seen from the fact that the vison operator $\mu_l^z$ is non-local in terms of the spin operator
$$\mu_l^z=\prod_{i<l}\sigma_i^x=\tau_l^x,$$
which flips all the spins to its left to create (or annihilate) a kink. Now let us discuss more about this infinite string of $\sigma^x$ stretching all the way to the left. A first question is that can we "gauge" this string to right? This can be done by applying the operator $S\equiv\prod_i\sigma_i^x$, as
$$\mu_l^z\to S\mu_l^z= \prod_{i}\sigma_i^x\prod_{i<l}\sigma_i^x = \prod_{i>l}\sigma_i^x.$$
One can see the operator $S$ simply flips all the spins in the system, meaning that it actually implement the global Ising symmetry transformation to the spins $\sigma_i^z\to-\sigma_i^z$. Because $S$ is a symmetry of the Hamiltonian (as $[S,H]=0$), the eigenstates of $H$ are spitted into the even ($S=+1$) and the odd ($S=-1$) sectors. In the even sector, we can gauge the string to the right; while in the odd sector, gauging the string to the right will induce a $\mathbb{Z}_2$ gauge transformation of the vison $\mu_l^z\to-\mu_l^z$. So combined with the gauge transformation of the vison, the vison string can be made invisible indeed. Then what is the significance of the this vison string? Recall that $\mu_l^z=\tau_l^x$ is also the creation/annihilation operator of the the kink. So applying a string of $\sigma^x$ operators will actually create two kinks at both ends of the string
$$\tau_{l_1}^x\tau_{l_2}^x=\prod_{l_1<i<l_2}\sigma_i^x.$$
The kink is a local excitation of the system (in the $J>h$ phase), so it can be viewed as a particle. From this perspective, we can see emergent particles at the ends of the string, which is exactly one of the central theme of Levin-Wen model and string-net condensation.
In fact, there is a very interesting relation between the 1D transverse field Ising model and the 2D Levin-Wen model (with $\mathbb{Z}_2$ topological order), that the former can be considered as a synthetic dislocation of the latter by anyon condensation, which was described in a paper (arXiv:1208.4109) that I wrote with my friend Chao-Ming and Prof. Wen. We basically showed that the 1D transverse field Ising model can emerge in the 2D $\mathbb{Z}_2$ topological ordered system as some kind of line defect, where some particular type of the strings between the anyons in 2D will naturally degrade to the vison strings along the emergent 1D Ising chain. So in this sense, the invisible string in the Ising model is really the same invisible string in the string-net condensate (but just confined to the 1D system).
Best Answer
After thinking about it I must say it is not as simple as I thought it would be. The JW transformation on the transverse Ising model contains quite a few subtleties.
So to proceed,
1) Take your ground state for ANY $h$ expressed in the spinless fermion language. I stress ANY because this condition is true always - it's not just for $h<1$. Now this is the vacuum, specified by $| 0 \rangle$ s.t. $a_k |0\rangle = 0$. This is a non-trivial condition written in the spin-language, i.e. we take the operators $a_k$, and do the following:
2) Apply the reverse Bogoliubov transform on $a_k$: $\{a_k\}\to \{b_k\}$.
3) Apply an inverse Fourier transform: $\{b_k\} \to \{b_i\}$
4) Apply the inverse Jordan Wigner transformation: $b_i = f(\sigma^x,\sigma^y,\sigma^z)$ .
All these transformations are invertible (see this pdf for JW and inverse JW transformation), which is why you can do that. So composing all the maps one can express $a_k = g(\sigma^x, \sigma^y, \sigma^z)$, where $g$ is the highly non-trivial function.
Then one has to find the kernel of $g(\sigma^x, \sigma^y, \sigma^z)$, i.e. $|\psi\rangle $ s.t. $g(\sigma^x, \sigma^y, \sigma^z) |\psi\rangle = 0$. $|\psi \rangle $ is the ground state written in the spin-basis.
You can write a program to do it for you symbolically, but for all your effort what you'll end up with is a highly non-local ground state in the spin-basis because of all the Jordan-Wigner strings.
Remark:
There are many subtleties associated with this transformation. What is very often not mentioned when one derives the spectrum $\epsilon(k,h)$ is that the JW transformation MUST be performed separately on states with different parity in the Hilbert space.
This is because the imposition of periodic boundary conditions in spin space implies the imposition of periodic boundary conditions for ODD number of fermions but anti-periodic boundary conditions for EVEN number of fermions. This affects the Fourier transform. In the calculation of any macroscopic quantity in the thermodynamic limit, there is no difference, and many books/resources just discard talking about the two cases. But this distinction must be made if one wants to be careful about it.
One question that arose to me when I was thinking about this problem was: hm, in one limit, $h\to\infty$, the ground states is unique, while in the other limit $h \to 0$ the ground states is 2-fold degenerate. Can I see that in the fermion language easily?
There are two resolutions I can think to the problem.
1) It could be that the ground state $|0\rangle_k$ for each $k$ is not unique. That is, instead of the irreducible 2-dimension representation of the fermionic CAR that we usually assume $a_k$ acts in, $a_k$ could be operators in a $2 \times d$ (reducible) dimensional representation, with $d$ 'ground states'.
2) The even and odd sectors of the full Hilbert space give rise to two conditions on the ground state: $a_k$ in the even sector gives a condition $g(\sigma^x, \sigma^y, \sigma^z) |\psi \rangle = 0$ while $a_k$ in the odd sector gives another condition $g'(\sigma^x, \sigma^y, \sigma^z)|\psi'\rangle$ = 0.
It could perhaps be the case that when $h\to \infty$, $|\psi\rangle = |\psi'\rangle$, while when $h \to 0$, $|\psi\rangle \neq |\psi'\rangle$.
It would seem more likely to me that 2) is the correct analysis, though it's going to be one tough assertion to prove.