In physics, this is a resonance. Indeed, you can "follow" the energy eigenvalue when it has an analytic dependence in the parameter. However, it will have a complex energy, the imaginary part being proportional to the lifetime of the "state." This can only happen in infinite dimensional Hilbert spaces.
Mathematically, you can understand this using the resolvent as you suggested. Recall that since $H$ is hermitian, $R$ is analytic on the lower and upper complex half planes, with isolated poles and branch cuts on the real line (isolated eigenvalues and continuum bands respectively). The trick is to analytically continue the resolvent across the cuts, and view the usual definition of the resolvent as the principal branch of a multivalued function on the complex plane. Resonances are precisely poles that appear in the other branches. When the bounds state "disappears", it can therefore survive by crossing the cut and moving to another branch.
Spectrally, you can also understand this by considering the density of states. When a bound state disappears to become a resonance, it gives rise to an intensive contribution to the density of states. You can therefore view it as a continuous version of degeneracy lifting.
Physically, these resonances have many observable consequences. Typically, you can observe these resonances in the scattering properties (enhancement of cross sections etc.) which is why it is relevant experimentally. Since the resolvent is the time-frequency Fourier transform of the impulse response, the life time is observable in the evolution of a localised state.
1D
To illustrate this, consider for example the 1D finite well:
$$
H = -\partial_x^2+V \\
V = \begin{cases} v & |x|<1 \\ 0 & |x|>1 \end{cases}
$$
and focus on the even parity states. You can explicitly compute the resolvent in the real space basis:
$$
R(x,y) = \langle x|(H-\omega)^{-1}|y\rangle.
$$
For the the density of states $\rho(\omega)$, you just need the diagonal terms, using:
$$
\rho(\omega) = \frac1\pi\Im\operatorname{Tr}R(\omega+i0^+)
$$
In this case, the extensive contribution of the trace is:
$$
\operatorname{Tr}R = \frac12\frac{1+(\kappa+1)\frac{\tanh\lambda}\lambda}{\kappa+\lambda\tanh\lambda}+\frac1{(2\kappa)^2}\frac{\kappa-\lambda\tanh\lambda}{\kappa+\lambda\tanh\lambda}
$$
(modulo some mistakes...) with:
$$
\kappa = \sqrt{-\omega} \quad \lambda = \sqrt{v-\omega}
$$
It is important to note that for $R$, you need to use the usual principle branch of the square root for $\kappa$ (cut on $\mathbb R_-$), but for $\lambda$ it does not matter as $R$ is really a function of $\lambda^2$ (purely for convenience).
There are naturally two extreme cases: $v=0^-$ and $v=-\infty$. In the former, you only have resonances, while in the latter you only have bound states. As you decrease $v$, the resonances move from the second branch to the principal branch through the branch point $0$ and are converted to proper bound states. If you have issues picturing it, you can use a complex function plotter.
Explicitly, bound states are in the principal branch where:
$$
\kappa+\lambda\tanh\lambda = 0
$$
whereas the resonances are in the second branch where:
$$
-\kappa+\lambda\tanh\lambda = 0
$$
General case
The same applies in higher dimensions and higher angular momentum, you just need to use the appropriate Bessel functions. I don't think that there is any trick, you just need to roll up your sleeves and solve for the inhomogeneous Schrödinger equation. In $d$ dimensions and angular momentum $l$, the equation for the resolvent is:
$$
-r^{1-d}\partial_r(r^{d-1}\partial_r\psi)+\frac{l(l+d-2)}{r^2}\psi+V\psi = \omega\psi +s^{1-d}\delta(r-s)
$$
you can revert to $d=2$ by changing function $\psi = r^{1-d/2}\phi$, setting $\alpha = l+d/2-1$:
$$
-\phi''-\frac1r\phi'-\frac{\alpha^2}{r^2}\phi+(V-\omega)\phi = s\delta(r-s)
$$
with the new inner product:
$$
\langle\psi|\psi\rangle = \int_0^\infty |\psi|^2r^{d-1}dr = \int_0^\infty|\phi|^2rdr
$$
Using the same notation:
$$
\kappa = \sqrt{-\omega}\quad\lambda = \sqrt{v-\omega}
$$
you obtain:
- for $s>1$:
$$
R(r,s) = \begin{cases}
\frac{I_\alpha(\lambda r)}BK_\alpha(\kappa s) & r\leq 1\\
\left(I_\alpha(\kappa r)+\frac{A_>}BK_\alpha(\kappa r)\right)K_\alpha(\kappa s) & 1\leq r\leq s\\
\left(I_\alpha(\kappa s)+\frac{A_>}BK_\alpha(\kappa s)\right) K_\alpha(\kappa r) & s\leq r
\end{cases}
$$
where:
$$
A_> = I_\alpha(\lambda)\kappa I_\alpha'(\kappa)-\lambda I_\alpha'(\lambda)I_\alpha(\kappa) \\
B = \lambda I_\alpha'(\lambda)K_\alpha(\kappa)-I_\alpha(\lambda)\kappa K_\alpha'(\kappa)
$$
- for $s<1$:
$$
R(r,s) = \begin{cases}
\left(K_\alpha(\lambda s)+\frac{A_<}BI_\alpha(\lambda s)\right)I_\alpha(\lambda r) & r\leq s\\
\left(K_\alpha(\lambda r)+\frac{A_<}BI_\alpha(\lambda r)\right)I_\alpha(\lambda s) & s\leq r\leq 1\\
\frac{K_\alpha(\kappa r)}BI_\alpha(\lambda s) & 1\leq r
\end{cases}
$$
where:
$$
A_> = K_\alpha(\lambda)\kappa K_\alpha'(\kappa)-\lambda K_\alpha'(\lambda)K_\alpha(\kappa) \\
B = \lambda I_\alpha'(\lambda)K_\alpha(\kappa)-I_\alpha(\lambda)\kappa K_\alpha'(\kappa)
$$
You recover the previous case when $\alpha = -1/2$. It is a bit trickier to check this time, but you for fixed $r,s$, you only have a branch cut at $E=0$ due to $\kappa$. Even if the formulas suggest a branch cut at $E=v$ due to $\lambda$, $R$ is a function of $\lambda^2$ like in the detailed 1D case. The trace is obtained by integrating:
$$
\begin{align}
\text{Tr }R &= \int_0^\infty R(r,r)rdr\\
&= \int_0^\infty \psi_<\psi_>rdr \\
\psi_> &= \begin{cases}
I_\alpha(\lambda r) & r\leq 1 \\
I_\alpha(\kappa r)+\frac{A_>}BK_\alpha(\kappa r) & 1\leq r
\end{cases} \\
\psi_< &= \begin{cases}
K_\alpha(\lambda r)+\frac{A_<}BI_\alpha(\lambda r) & r\leq 1 \\
K_\alpha(\kappa r) & 1\leq r
\end{cases}
\end{align}
$$
Formally, the integral diverges, since you have an extensive contribution (proportional to the infinite volume of space), due to the unbounded states. The finite contribution corresponds to the resonances.
The previous discussion of resonances still apply. $R_l$ when viewed as a function of $E$ is a double branched analytic function with a branch point at $E=0,\infty$. The bound states are simple poles on $E\in\mathbb R_-^*$ of the principal branch, and the complex poles in the second branch correspond to resonances. As you vary $v$, they will cross the branch point and be converted into proper bound states.
Up to now, only the Hamiltonian restricted to angular momentum $l$ was considered. The full trace is obtained by summing all the contributions:
$$
\text{Tr }R = \sum_{l=0}^\infty (l+1)R_l
$$
On a side note, there is a link between resolvent formalism and the usual scattering approach. The trick is to compute the discontinuity of the the the resolvent at the cut:
$$
\text{Disc }R = \lim_{\epsilon\to0^+}R(E+i\epsilon)-R(E-i\epsilon)
$$
This will give you the scattered state and allow you to recover the $S$ matrix. You can check this explicitly by using the relation between the modified Bessel function $K_\alpha$ with the Hankel functions $H_\alpha$.
Hope this helps.
