I don't see how you get mixed terms. Your $U$ acts on every other site. For instance (as operators acting on different sites commute) $$\ldots+ U \sigma^{2i}_x \sigma^{2i+1}_x U^{-1} +\ldots= \ldots \sigma^{2i}_z \sigma^{2i}_x \sigma^{2i}_z \sigma^{2i+1}_x \ldots=\ldots i \sigma^{2i}_y \sigma^{2i}_z \sigma^{2i+1}_x \ldots=\ldots-\sigma^{2i}_x \sigma^{2i+1}_x\ldots,$$ and similarly for the Pauli $y$-matrices.
EDIT:
Maybe you made a mistake multiplying the Pauli matrices?
Here is the full multiplication table for the Pauli matrices:
http://www.its.caltech.edu/~hmabuchi/Ph125b/HW7-solutions.pdf
$J$ can be omitted because it is not zero compared to $\Delta$ and you can divide (rescale) the Hamiltonian by $J$ and absorb it into $\Delta/J \rightarrow \Delta'$. If $\Delta \rightarrow \infty$ then you just recover the Ising model.
Remember, this transformation is unitary, which means it's just a change of basis in Hilbert space. It can't change anything "physically" and that tells us that it's only the relative $\pm$ - sign between $J$ and $\Delta$ that is physical.
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.
Best Answer
1.) The easiest way to count the energy is as that of domain walls. Both of the configurations you have drawn have two domain walls. Each domain wall costs energy $E\Delta/2$, so both configurations have energy $E\Delta$ above the ground state. I don't understand either how one would count the spin of a domain as 1. Perhaps that is the spin/site?
2.) To see that a domain wall between the two antiferromagnetic ground states has spin $ 1/2$, look at a spin configuration with a domain wall and divide the spin at each site evenly between the adjacent links. For example, if there is a spin $S_z=1/2$ at site $i$ and a spin $S_z=-1/2$ at site $i+1$, you'd say there is $0$ spin on the link $<i,i+1>$. Likewise, if there is a spin $S_z=1/2$ at site $i$ and a spin $S_z=1/2$ at site $i+1$, you'd say there is spin $S_z=1/2$ on the link $<i,i+1>$. In this way you end up with 0 everywhere except where there are two spins of the same type in a row (i.e., where there is a domain wall). At the domain wall you get $S_z=\pm 1/2$.