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.
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
Best Answer
It is important to note that there are two different boundary conditions, the first is the boundary condition for real spin model, the second is for fermion model. In fact, for real spin model with a certain boundary condition, there may be different corresponding boundary conditions for fermion operator. When we make a simulation in reality, we need to deal with this "boundary conditions problem" carefully.
Jordan-Wigher transformation
For a Ising model with spin-$\frac{1}{2}$, we can maop it into a spinless fermion model via Jordan-Wigher transformation: $$\sigma_i^z=1-c_i^\dagger c_i$$ which means we map the spin-up state to empty state, and map spin-down state to occupied state: $$|\uparrow\rangle \longmapsto |0\rangle\\|\downarrow\rangle \longmapsto |1\rangle$$ Naively speaking, we can unify the creation and annihilate operator for fermion with the ladder operator for spin: $c_i^\dagger\sim \sigma_i^-$, $c_i\sim \sigma_i^+$. However, the fermion operator follows the anti-communication relation and spin operator behaves like hard-core boson, which follows communication relation. To unify such two different communication relation, we need to introduce an additional string: $$\sigma_i^+=P_{i-1}c_i\\\sigma_i^-=P_{i-1}c_i^\dagger$$ where $P_{i}$ is the so-called string operator, which is defined as $$P_{i}\prod_{j=1}^i(1-c_j^\dagger c_j)$$. The effect of string operator is to measure the parity, even or odd, of the number of fermion which on the left-side of i-site.
For simplicity, we will rotate the axis: $$\sigma_z \longmapsto \sigma_x\\\sigma_x \longmapsto -\sigma_z$$
Open boundary condition(OBC)
For the real spin system, i.e. Ising model, with open boundary condition: $$H_I=-J\sum_{i=1}^N g\sigma_i^x-J\sum_i^{N-1}\sigma_i^z\sigma_{i+1}^z$$ we can re-express it via fermion operator: $$H_I=-Jg\sum_{i=1}^N (1-2c_i^\dagger c_i)-J\sum_i^{N-1}P_{i-1}(c_i^\dagger+c_i)P_i(c_{i+1}^\dagger+c_{i+1}) \\=-Jg\sum_{i=1}^N (1-2c_i^\dagger c_i)-J\sum_i^{N-1}P_{i-1}(c_i^\dagger+c_i)P_{i-1}(1-2n_i)(c_{i+1}^\dagger+c_{i+1}) \\=-Jg\sum_{i=1}^N (1-2c_i^\dagger c_i)-J\sum_i^{N-1}(c_i^\dagger+c_i)(1-2n_i)(c_{i+1}^\dagger+c_{i+1}) \\=-Jg\sum_{i=1}^N (1-2c_i^\dagger c_i)-J\sum_i^{N-1}(c_i^\dagger-c_i)(c_{i+1}^\dagger+c_{i+1})$$
We can find that the OBC of real spin model will be consistent with the OBC of spinless fermion model after Jordan-Wigher transformation.
Periodic boundary condition(PBC)
For the real spin system, i.e. Ising model, with periodic boundary condition: $$\sigma_{N+1}=\sigma_1$$ there is an additional term: $$-J\sigma_N^z\sigma_{N+1}^z \\=-J\sigma_N^z\sigma_{1}^z \\=-JP_{N-1}(c_N^{\dagger}+c_N)(c_1^{\dagger}+c_1) \\=-JP_{N}(c_N^{\dagger}-c_N)(c_1^{\dagger}+c_1)$$ the $P_N$ operator measure the parity of the number of fermion in the whole system. It is important to note that for the whole system with odd fermion number, $P_N=1$, we can unify this additional term to the normal expression of fermion model via $c^\dagger_{N+1}=c_i^\dagger$, which means the PBC for fermion model. On the other hands, for the whole system with even fermion number, $P_N=-1$, we can unify this additional term to the normal expression of fermion model via $c^\dagger_{N+1}=-c_i^\dagger$, which means the APBC for fermion model. As the result, for the real spin system with periodic boundary condition, the corresponding periodic boundary condition of fermion model has two situation: - PBC: if total number of fermion is odd, the spin-less fermion follows PBC. The effect of PBC will constrain the value of momentum. Namely, after Fourier transformation, $c_{k}=\frac{1}{\sqrt{N}} \sum_{j=1}^{N} c_{j} \exp (-i k j)$, where $k=2 \pi n / N$, the PBC restricts $n$ can only be integer. - APBC: if total number of fermion is even, the spin-less fermion follows APBC. which restricts $n$ can only be half-integer.
As the result, the spin-less fermion Hamiltonian can be written in the momentum space: $$H=J \sum_{k}\left[2(g-\cos k) c_{k}^{\dagger} c_{k}+i \sin k\left(c_{-k}^{\dagger} c_{k}^{\dagger}+c_{-k} c_{k}\right)-g\right]$$ where the value of $k$ depends on the boundary condition.
After Bogoliubov transformation: $$\gamma_{k}=u_{k} c_{k}-i v_{k} c_{-k}^{\dagger}$$ we can obtain the final diagonal Hamiltonian: $$H_{I}=\sum_{k} \varepsilon_{k}\left(\gamma_{k}^{\dagger} \gamma_{k}-1 / 2\right)$$ where the dispersion is: $$\varepsilon_{k}=2 J\left(1+g^{2}-2 g \cos k\right)^{1 / 2}$$
Unifying boundary condition
Now, there will emerge a contradiction. For a real spin-$\frac{1}{2}$ model with PBC, assuming there are $N$ sites, the number of state is $2^N$. However, for the spin-less fermion with certain boundary condition, $k$ have $N$ values, thus, there are totally $2^N$ states for a certain boundary condition. In the other words, if we consider two boundary conditions, i.e. PBC and APBC for fermion model, there are $2\times2^N$ states, which means half of them are redundant.
In fact, not all the states which calculated from fermion model for a certain boundary condition are physical, we only need to take the states with correct fermion number. For example, when we choose $k=2 \pi n / N, n\in Z$, which satisfies the PBC, the number of total fermion should be odd, thus, we need to leave out all the states with even fermion number. As the result, for each boundary condition, we just need to take half of them.
Till now, we have discussed how to deal with boundary condition problem for real spin model with PBC. For real spin model with APBC, the methods is similar.