I am writing programs to construct the spectra of models with known exact solutions, and soon noticed some subtleties that are not often mentioned in most references. These subtleties are not important in the thermodynamic limit, but since I am writing numerical recipes to finite size systems, they do matter. Now there is still one puzzle concerning Bogoliubov transformation (or so I believe) that I cannot resolve on my own.
Take for example, the 1D anisotropic XY model in a transverse field (including the case of Ising model), given by the following Hamiltonian
$$ H = – J \sum_{j=1}^N \left( g \sigma_j^z + \frac{1+\eta}{2} \sigma_j^x \sigma_{j+1}^x + \frac{1-\eta}{2} \sigma_j^y \sigma_{j+1}^y\right), $$
assuming periodic boundary condition $\vec{\sigma}_{N+1} = \vec{\sigma}_1$.
This model can be solved by applying Jordan-Wigner transformation, Fourier transform, and Bogoliubov transformation successively. The final diagonal form of the Hamiltonian, as given in many references, is
$$ H = \sum_k \varepsilon(k) (b^\dagger_k b_k – \frac{1}{2} ) $$
where $ \varepsilon(k) = 2 J \sqrt{(g – \cos(ka))^2 + \eta^2 \sin^2(ka)}$.
Two known subtleties:
- The original Hilbert space of periodic spin chain is mapped to TWO sectors of Fock space: odd/even parity sector ($n = \sum c^\dagger_j c_j$ odd/even) with periodic/anti-periodic boundary condition $c_1 = \pm c_{N+1}$ respectively.
- The momentum $k$ takes values of integer/half-integer multiples of $2\pi/Na$ in the odd/even parity sector respectively.
Another subtitles that (almost) immediately follows (but unresolved) is that the dispersion relation given above actually only gives the absolute value of a single fermion's energy. The true single fermion energy should be given by
$$ \varepsilon(k) = \pm 2 J \sqrt{(g – \cos(ka))^2 + \eta^2 \sin^2(ka)}. $$
This undetermined sign matters in a finite size system because there are two parity sectors with different choices of $k$, hence the "presence" and "absence" of fermions with "opposite energies" are not equivalent (in fact, they can't have exactly opposite energies in general). Even if this sign ambiguity does not affect total energy $E_n$ itself, it most certainly still disrupts the dispersion relation of the total spectrum $E_{n,k}$ (this is the case when $g=\eta=0$).
I believe this sign is determined in the Bogoliubov transformation, but I do not know how. The Hamiltonian right before the Bogoliubov transformation is
$$ H = J \sum_k \left\{ 2(g- \cos(ka)) c^\dagger_k c_k + i \eta \sin(ka) (c^\dagger_{-k} c^\dagger_k + c_{-k} c_k ) \right\} + \text{const.} $$
For example, in the isotropic case $\eta = 0$, by comparing test cases with results from exact diagonalization, I know the correct single fermion energy is given by $ \varepsilon(k) = 2J (g – \cos(ka)) $, not any other variants. However, when making the Bogoliubov transformation:
$$ c_k = \cos \frac{\theta_k}{2} b_k + i \sin \frac{\theta_k}{2} b^\dagger_{-k}, $$
and assuming $\cos \frac{\theta_k}{2} = \cos \frac{\theta_{-k}}{2}, \;
\sin \frac{\theta_k}{2} = -\sin \frac{\theta_{-k}}{2} $, we get
$$ H_{\eta=0} = J \sum_k \left\{ (g-\cos(ka) )
\left[ 2 \cos \theta_k b^\dagger_k b_k + i \sin \theta_k ( b^\dagger_k b^\dagger_{-k} – b_{-k} b_k ) \right] \right\} + \text{const.} $$
The only restriction that I can see is $\sin\theta_k = 0$. If that's the case, then $\cos\theta_k$ can take $\pm 1$ however it likes. This is NOT what happens! Only $\cos\theta_k = 1$ gives the correct solution!
If someone could point out what I am or might be missing here, or point me to some relevant reference (most references don't discuss finite size case in such details), that would be great. Thanks!
Best Answer
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:
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:
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