Your (unbounded) operator is in the limit circle case at $+\infty$. This means that self-adjoint realizations are obtained only when a boundary condition at $\infty$ is imposed, and the eigenvalues will depend on this boundary condition. (In fact, for any $\lambda\in\mathbb R$, there is exactly one boundary condition that makes this $\lambda$ an eigenvalue.)

To see this, let's make the (obvious) try $u_n=q^{-2n}g_n$. Then $Tg=\lambda g$ (viewed as a difference equation) is equivalent to
$$
q^{1/2}u_{n+1}+q^{-3/2}u_{n-1}+\frac{1+q}{q}u_n = \lambda q^{2n} u_n .
$$
The coefficient $\lambda q^{2n}$ is summable on $n\ge 1$, so can be ignored as far as the asymptotic behavior is concerned; to prove this carefully, we would have to refer to a discrete version of Levinson's theorem.

If we set the RHS equal to zero, then two basis solutions are $u_n=(-q)^{-n/2}$ and $u_n=(-q)^{-3n/2}$. Going back to $g$, we thus have two solutions $g_n\simeq q^{3n/2}$ and $g_n\simeq q^{n/2}$, both of which are exponentially decaying and in $\ell^2$. All solutions are square integrable, as claimed.

We can also observe that the coefficients of the left half line operator $T_-$ go to zero, so $\sigma_{ess}(T_-)=\{ 0\}$. This shows that (for any boundary condition at $\infty$) $\sigma_{ess}(T)=\{ 0\}$.