I don't know the article you refer to, but I believe the Hamiltonian you discuss should get a $\pi$-phase shift after one turn around a (2D) lattice cell. So I guess it should read $H=F^{\dagger}\cdot H_{\pi}\cdot F$ with
$$H_{\pi}=t\left(\begin{array}{cccc}
0 & e^{\mathbf{i}\pi/4} & 0 & e^{-\mathbf{i}\pi/4}\\
e^{-\mathbf{i}\pi/4} & 0 & e^{\mathbf{i}\pi/4} & 0\\
0 & e^{-\mathbf{i}\pi/4} & 0 & e^{\mathbf{i}\pi/4}\\
e^{\mathbf{i}\pi/4} & 0 & e^{-\mathbf{i}\pi/4} & 0
\end{array}\right)$$
and $F^{\dagger}=\left(\begin{array}{cccc}
f_{1}^{\dagger} & f_{2}^{\dagger} & f_{3}^{\dagger} & f_{4}^{\dagger}\end{array}\right)$. Then, one has
$$H_{\pi}=\dfrac{t}{\sqrt{2}}\left[\left(1+\tau_{x}\right)\otimes\eta_{x}-\left(1-\tau_{x}\right)\otimes\eta_{y}\right]$$
where the $\eta$ and $\tau$ are the usual Pauli matrices.
Time reversal symmetry operator -- when it exists -- is defined as an anti-unitary operator which commutes with the Hamiltonian. Such an operator can be defined as $T=\mathscr{K}\tau_{z}\otimes\mathbf{i}\eta_{y}$ and thus $H$ is time reversal symmetric. $\mathscr{K}$ is the anti-unitary operator $\mathscr{K}\left[\mathbf{i}\right]=-\mathbf{i}$ and thus $\mathscr{K}\left[\eta_{y}\right]=-\eta_{y}$. One verifies that $\left[H_{\pi},T\right]=0$ as it must.
Please tell me if I started with the wrong Hamiltonian.
A few words about the definition (as follow from the comment below): The time-reversal operator is defined as I did, i.e. one applies it to the Hamiltonian $H_{\pi}$, (call it the Hamiltonian density if you wish, since in my way of writing $H=F^{\dagger}\cdot H_{\pi}\cdot F$, the dots should include summation(s) over phase-space-time [delete as appropriate]). You could prefer to define the action of an operator as transforming the operators (or the wave-function). But you should not use both definitions at the same time. It is clear that you can not do both, since otherwise you transform $H=F^{\dagger}\cdot H_{\pi}\cdot F \rightarrow F^{\dagger}\cdot U^{\dagger}\cdot \left(U \cdot H_{\pi} \cdot U^{\dagger}\right) \cdot U\cdot F = H$ trivially, whatever (anti-)unitary transformation $U$ you choose. It is clear that what your are looking for is something like $H=F^{\dagger}\cdot H_{\pi}\cdot F \rightarrow F^{\dagger}\cdot U^{\dagger}\cdot H_{\pi} \cdot U\cdot F \sim H$ and you see what I just said: apply the transformation to the Hamiltonian (density) or to the fields, but not both. In condensed matter we usually choose the convention I gave to you: we transform the Hamiltonian. One of the reasons is that the operators (especially the fermionic creation/annihilation ones) are seen as encoding the statistics of the fields, whereas the Hamiltonian encodes the dynamics, and it is simple imagination to change the dynamics.
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
Basically, the answer is yes: $H$ is TRI because it is real. Reality condition really means that the Hamiltonian obeys a certain anti-unitary symmetry. In this case, the time-reversal operation is simply $T=K$ where $K$ is the complex conjugation. It is not the usual one($T=K\prod_i i\sigma^y_i$), and in particular $T^2=1$, so there is no Kramers' theorem and the spectrum is not doubly degenerate. The fact that level statistics follows GOE of course is a consequence of the reality condition. In fact, I think if there was a $T^2=-1$ time-reversal symmetry, the statistics would follow a different ensemble (the symplectic one, I believe).
You asked what if one changes the transverse field $g\sum_i \sigma^x_i$ to the $y$ direction. In that case, the two Hamiltonians are unitarily related (i.e. a $\pi/2$ spin rotation $U_z$ around $z$ would bring it back). Let me call the Hamiltonian with $x$ transverse field $H_x(g)$ where $g$ is the transverse field, and with $y$ transverse field $H_y(g)$. Define $U_z=e^{i\pi \sum_i\sigma^z_i/4}$, then it is easy to check that $U_z H_y(g) U_z^{-1}= H_x(g)$. Since we know $H_x^*(g)=H_x(g)$, we can easily find $H_y^*=U_z^2 H_y U_z^{-2}$. Therefore one just has to redefine the time-reversal symmetry to be $T=K U_z^2$. If you really want to break the reality condition, in a way that can not be fixed by additional unitary transformations, then one needs to turn on transverse fields along all three dimensions.
Last comment on your "Arguments for yes": the first argument you gave, namely one also flips the external parameters, does not work. In this way, there would be no time-reversal symmetry breaking, except the CP violation in the fundamental processes! When we talk about the symmetry of a Hamiltonian, we should just treat the system on its own, not with all the external devices that generate the various terms -- unless you want to consider the dynamics of these devices, but then it is a different Hamiltonian.