How to find $n$ of the $n$-fold rotation symmetry explicitly for a $4times 4$ Hamiltonian?

Given a 2D crystal Hamiltonian $H(k_x,k_y)$ in terms of momenta $k_x$ and $k_y$, how do I use the rotation operator $R_n(theta)=e^{-i theta vec{J}cdotvec{n}}$ (or otherwise) to find the Hamiltonian’s possible $n$-fold rotational symmetry.
I know that if $H(k_x,k_y)$ has $n$-fold rotational symmetry, then $R_n(theta) H(k_x,k_y) R_n^{-1}(theta)=H(k_x,k_y)$.
Since this is a 2D crystal, I expect this $vec{n}$ to be in the $k_z$ direction, perpendicular to the $(k_x,k_y)$ plane.
I expect to plug in various $n$ to $theta=2pi/n$ in $R_n(theta)$ to find which $n$ would work out.

This post discusses this for a $2times2 $ system, but I am lost on which matrices to use for $vec{J}$ in $R_n(theta)$. This is because I have a bunch of $k$-space Hamiltonians written in the basis of valence band and conduction band states, like the one below from this paper:

I doubt I can use a general rotation matrix such as the one below that rotates about the $xy$ plane, because it is independent of $k_x$ and $k_y$:
$$
R_n(theta)=
begin{bmatrix}
cos{theta} & -sin{theta} & 0 & 0
sin{theta} & cos{theta} & 0 & 0
0 & 0& 1 &0
0 & 0& 0 & 1
end{bmatrix}.
$$

Or do I use a general expression for $L_z$?
$$
L_z
=
– i h (x frac{partial}{partial y}-y frac{partial}{partial x})
$$?

Or instead of matrices, do I simply use $R=e^{i theta (k_x + k_y)}$?

Any suggestions or resources to help me tackle this?