Let
$$f(z)=\frac{e^{iz}-1}{iz}.$$
This function is in the Hardy class for any upper half-plane,
and has these properties: $f(0)=1,$ $f(2\pi n)=0$,
$$|f(z)|\leq C\frac{e^y+1}{|z|+1},$$
(this evidently holds for large and small $|z|$, therefore there is a constant $C$ so that this holds everywhere).
Since the $L^2$ norm has the property $\| f(kx)\|^2=\| f(x)\|^2/k,$ and the norm is invariant under real translations, the following function $g$ belongs to the Hardy spaces for upper half-planes:
$$g(z)=\sum_{n=0}^\infty f(n^2(z-a_n)),$$
where $a_n\in\{2\pi k: k\in {\mathbf{Z}}\}$.
Moreover, it is entire, if the sequence $a_n$ grows fast enough. (It follows from the estimate for $f$ that the series is uniformly convergent on any compact subset of the plane).
Now $g(a_n)=1$, so $h(z)=(g(z-i)-1)/(z-i)$, where has all properties that you required.
Some extra labor is needed to show that $h$ can be made of finite order.
What is clear here is that $H_2$ has spectrum above $4$, and that $\lambda\in\sigma_{\textrm{ess}}(H_2)$ also, with $\lambda:=\max\sigma(H_2)=\|H_2\|>4$. What $\lambda$ is actually equal to seems a rather delicate question. Certainly $\lambda<20$, and I'm also not at all sure that $\lambda\ge 12$; that may already be too ambitious.
To discuss these claims, let me switch to the continuous counterpart of your problem, the operator $H=-\Delta-V(x,y)$ on $L^2(\mathbb R^2)$. Everything I'm going to say has a precise analog in your (discrete) setting, but the details become considerably less tedious in the continuous case. Even so, I will only sketch most of them.
We now assume that $V(x,y)=16$ on $|x|\le 1$, $y\ge 0$, and $V=0$ otherwise. We are interested in the part of the spectrum below zero. (We clearly have $[0,\infty)\subseteq\sigma (H)$.)
First of all, there actually is negative spectrum. This would not be clear in dimension $3$ or higher, but here it still works. Compare my answer here. We can find $\lambda=\min\sigma (H)$ as the infimum of the quadratic form $\lambda=\inf Q(f)$,
$$
Q(f)=\int_{\mathbb R^2} \left(|\nabla f|^2-V|f|^2\right) , \tag{1}
$$
taken over $f\in C_0^{\infty}$ (say) with $\|f\|_{L^2}=1$. Clearly, we can only hope to make $Q(f)$ small if we concentrate $f$ where $V=16$, but then $Q(f)>-16$ because we can not avoid paying a (potentially steep) price in the first term. A more chunky potential that is equal to $16$ on a large disk rather than a thin strip would be more effective at pulling down the bottom of the spectrum.
In any event, we can try concrete test functions in (1) and get upper bounds on $\lambda$, but finding $\lambda$ exactly seems a formidable task.
Finally, there are $f\in L^2$ that make $\|(H-\lambda)f\|$ arbitrarily small. A carefully cut off version (compare the discussion in the comments to the linked answer) of $f$ will still keep $(H-\lambda)f$ small, and then we obtain a Weyl sequence by also shifting and considering $g(x,y)=f(x,y-L)$, $L\gg 1$.
Best Answer
Here is a suggestion (but I have not worked through the details). What follows is given for the upper-half plane (UHP for short) but of course a trivial rotation will convert the example to one on the right-half plane as you requested.
There is a weighted composition operator giving an isometry of Hilbert spaces from $H^2$ of the UHP to $H^2$ of the open unit disc. One can find an explicit formula on the Wikipedia page for Hardy spaces. If $M: H^2({\rm UHP})\to H^2({\bf D})$ is the particular isometry described at that link, then inspecting the formula there shows that for each $f\in (H^2\cap H^\infty)({\rm UHP})$, the function $M(f)$ belongs to $H^2({\bf D})$ and is essentially bounded as $z\to -1$.
It therefore suffices to pick any $g\in H^2({\bf D})$ which is not essentially bounded as $z\to -1$, and then $M^{-1}(g)$ will be a function in $H^2({\rm UHP})$ that is not essentially bounded. One such $g$ is given by $$ \displaystyle g(z)=\sum_{n=1}^\infty \frac{(-1)^{n-1}}{n}z^n=\log(1+z). $$