How to understand the time-reversal symmetry in graphene?

Two-band model for graphene. To simplify the discussion of Dirac points, it is sufficient to consider a nearest-neighbor tight-binding 2-band model for graphene. This is a spinless model because spin-orbit is negligible in graphene. The Bloch Hamiltonian is:

$$ hat{H}(mathbf{k})= begin{pmatrix} 0 & h(mathbf{k}) \ h^{dagger}(mathbf{k}) & 0 end{pmatrix}, $$

where $h(mathbf{k})=-tsum_j e^{imathbf{k}cdotmathbf{u}_j}$, with $t$ the hopping parameter, and $j=1,2,3$ running over the three $mathbf{u}_j$ nearest neighbors. In the pseudospin representations, this becomes:

$$ hat{H}(mathbf{k})=-tsum_jleft[sigma_1cos(mathbf{k}cdotmathbf{u}_j)-sigma_2sin(mathbf{k}cdotmathbf{u}_j)right]. $$

The energy spectrum is $E_{pm}(mathbf{k})=pm|h(mathbf{k})|$, and to obtain Diract points we must find a $mathbf{k}_0$ for which $h(mathbf{k}_0)=0$. This can generally happen in graphene because we have two terms in the Hamiltonian (one proportional to $sigma_1$ and one proportional to $sigma_2$) and two variables ($k_x$ and $k_y$) that we can tune to enforce the vanishing of $h(mathbf{k})$ at $mathbf{k}=mathbf{k}_0$. In graphene, the Dirac points occur at the K and $-$K points in the Brillouin zone.

General 2-level system in 2 dimensions. The Hamiltonian of a general $2$-level system can always be written in the basis of the three Pauli matrices. In two dimensions, this means that we have three terms ($sigma_1$, $sigma_2$, and $sigma_3$) but only two parameters ($k_x$ and $k_y$) to tune, so that a general 2-dimensional system has no crossing points. The reason graphene has crossing points is because the $sigma_3$ term is not present in the Hamiltonian. Why is is not present? Because there are some symmetries in the system that force it to be zero. As you correctly identified, they are time reversal symmetry and inversion symmetry.

Time reversal symmetry. Time reversal symmetry means that the Bloch Hamiltonian must obey:

$$ hat{H}(mathbf{k})=hat{H}^{ast}(-mathbf{k}). $$

You should convince yourself that the tight-binding Hamiltonian above obeys this. Imagine we add a term $f_3(mathbf{k})sigma_3$ to the Hamiltonian, to get:

$$ hat{H}_1(mathbf{k})= begin{pmatrix} f_3(mathbf{k}) & h(mathbf{k}) \ h^{dagger}(mathbf{k}) & -f_3(mathbf{k}) end{pmatrix}. $$

In this case,

$$ hat{H}^{ast}_1(-mathbf{k})= begin{pmatrix} f_3^{ast}(-mathbf{k}) & [h(-mathbf{k})]^{ast} \ [h^{dagger}(-mathbf{k})]^{ast} & -f_3^{ast}(-mathbf{k}) end{pmatrix} = begin{pmatrix} f_3^{ast}(-mathbf{k}) & h(mathbf{k})\ h^{dagger}(mathbf{k}) & -f_3^{ast}(-mathbf{k}) end{pmatrix}, $$

where we have used the fact that for the original Hamiltonian $hat{H}(mathbf{k})=hat{H}^{ast}(-mathbf{k})$ in the last step. If the system has time reversal symmetry, then the expression for $hat{H}^{ast}_1(-mathbf{k})$ must equal $hat{H}_1(mathbf{k})$, which implies $f_3^{ast}(-mathbf{k})=f_3(mathbf{k})$.

Inversion symmetry. Inversion symmetry implies $sigma_1hat{H}(-mathbf{k})sigma_1=hat{H}(mathbf{k})$, and you should again check that the original Hamiltonian obeys this. You can then also add a term $f_3(mathbf{k})sigma_3$ to the Hamiltonian, and proceeding in the same way I did for time reversal symmetry (I encourage you to try), you will find that this imposes the condition $f_3(-mathbf{k})=-f_3(mathbf{k})$.

Combined time reversal and inversion symmetries. If we have the term $f_3(mathbf{k})sigma_3$ in the Hamiltonian, and insist that both symmetries are obeyed then we have the following simultaneous conditions on $f_3$:

$$ begin{cases} f_3^{ast}(-mathbf{k})=f_3(mathbf{k}),\ f_3(-mathbf{k})=-f_3(mathbf{k}). end{cases} $$

Putting them together, they imply that $f_3^{ast}(-mathbf{k})=-f_3(-mathbf{k})$, which means that $f_3(mathbf{k})$ is purely imaginary. But for a Hermitian operator, $f_3(mathbf{k})$ must be purely real. These two conditions can only be obeyed simultaneously if $f_3(mathbf{k})=0$. Therefore, graphene has no term proportional to $sigma_3$ in the Hamiltonian because of a combination of time reversal and inversion symmetries. This in turn implies that graphene has Dirac points because of a combination of time reversal and inversion symmetries.

Opening a gap in graphene. If you break either of these two symmetries, then you no longer have $f_3=0$, and a gap develops. As you correctly identified, a term that breaks inversion symmetry is $f_3(mathbf{k})sigma_3=msigma_3$, which effectively places inequivalent on-site potentials on the two atoms in the primitive cell. This is what happens for example in monolayer hexagonal BN, which indeed has a gap. A term that you can add to break time reversal symmetry is to add a complex next-nearest-neighbor hopping $t^{prime}e^{ivarphi}$. This is the famous term that Haldane added to kick-start the field of topological materials, and you can read the original paper here. If you go through the maths, you will find that the Bloch Hamiltonian acquires the term

$$ f_3(mathbf{k})sigma_3=-2t^{prime}sum_jsin(mathbf{k}cdotmathbf{v}_j)sigma_3, $$

where I have chosen the complex phase to be $varphi=pi/2$ and $mathbf{v}_j$ are the three next-nearest-neighbors.

Spin-orbit coupling. The discussion is different if spin-orbit coupling is included. Although it is completely negligible in graphene ($mu$eV scale), it played a very important role historically. With spin-orbit coupling we have a 4-band model, and you can add a spin-orbit term that conserves time reversal and inversion symmetries, but still opens a gap. This is precisely what Kane and Mele did in their 2005 paper that led to the discovery of topological insulators.