You wish to approximate $$p_{\rm exact}(\rho)=e^{-\rho^2-\rho_0^2} \rho I_0(2\rho \rho_0).$$ For $\rho_0\ll 1$ a Taylor series in powers of $\rho_0$ is accurate. Let me consider the opposite regime of large $\rho_0$. An asymptotic expansion of the Bessel function gives $$p_{\rm approximate}(\rho_0)=\frac{1}{2\sqrt{\pi}}e^{-\rho^2-\rho_0^2}e^{2\rho\rho_0}.$$ This large-$\rho_0$ approximation to $\int_0^x p(\rho)\,d\rho$ is already quite accurate for $\rho_0= 3$, see the plot (blue = exact, orange = approximate, almost indistuishable for $\rho_0=3$).
