Start from a stochastic differential equation of the form: $$dr_t=a(b-r_t)dt+\sigma\sqrt{r_t}dB_t\tag{1}$$ with $\left(B_t\right)_{t\geq0}$ denoting a Brownian motion on the filtered probability space $\left(\Omega,\mathcal{F},\mathcal{F}_n,\mathbb{P}\right)$. By property of Brownian motion, $(1)$ can be rewritten as: $$dr_t=a(b-r_t)dt+\sigma\sqrt{r_tdt}Z\tag{2}$$ with $Z=\mathsf{N}(0,1)$ denoting a standard normal distribution, mean $0$ and standard deviation $1$.
I would like to numerically simulate the solution to $(2)$. To this, I could use for example Euler-Maruyama method, with simulation step-size equal to $\Delta t$, discretizing $(2)$ as follows: $$r_{t+\Delta t}=r_{t}+a(b-r_t)\Delta t+\sigma\sqrt{r_t\Delta t}\mathsf{N}(0,1)\tag{3}$$
I would like to know the probability associated to the event (for a given $t$):
$$\mathbb{P}\left(r_{t+\Delta t}<0|r_{t}=\varepsilon\right)\hspace{1.5cm}\text{for }\varepsilon>0\text{ arbitrarily small}$$
In other terms, if I suppose that the value of the solution process at previous time is quite small, which is the probability that at next simulation step the solution process will lie below $0$?
After several passages, I get that:
\begin{align}
\mathbb{P}\left(r_{t+\Delta t}<0\right)&=\mathbb{P}\left(\varepsilon+a(b-\varepsilon)\Delta t+\sigma\sqrt{\varepsilon \Delta t}Z<0\right)\\
&=1-\mathbb{P}\left(Z>-\dfrac{\sqrt{\varepsilon \Delta t}}{\sigma}\left(\dfrac{1}{\Delta t}+a\left(\dfrac{b}{\varepsilon}-1\right)\right)\right)\nonumber\\
&\underbrace{\leq}_{\text{by Chebyshev inequality}}1-\dfrac{\sigma^2 \Delta t \varepsilon}{\varepsilon^2+\left(\Delta t a b\right)^2+\left(a\varepsilon \Delta t\right)^2+2ab\Delta t\varepsilon-2a\varepsilon^2 \Delta t-2a^2\Delta t^2b\varepsilon}
\end{align}
At this point I would like to know if I can investigate the behaviour of $\mathbb{P}\left(r_{t+\Delta t}<0|r_{t}=\varepsilon\right)$ as $\varepsilon$ and/or $\Delta t$ vary (varies). Hence, in other terms, I would like to find an alternative (simplified) form for:
$$1-\dfrac{\sigma^2 \Delta t \varepsilon}{\varepsilon^2+\left(\Delta t a b\right)^2+\left(a\varepsilon \Delta t\right)^2+2ab\Delta t\varepsilon-2a\varepsilon^2 \Delta t-2a^2\Delta t^2b\varepsilon}\tag{4}$$
just to understand how $\mathbb{P}\left(r_{t+\Delta t}<0|r_{t}=\varepsilon\right)$ is affected by $\varepsilon$ and $\Delta$ (e.g. probability rises as $\varepsilon$ increases or probability lowers as $\Delta t$ increases, etc.)