For large $\rho_0$ you can use an asymptotic expansion of the Bessel function, which gives $$p_{\rm exact}=e^{-\rho^2-\rho_0^2} \rho I_0(2\rho \rho_0),$$ $$p_{\rm approximate}=\frac{1}{2\sqrt{\pi}}e^{-\rho^2-\rho_0^2}e^{2\rho\rho_0}.$$
The approximation to $\int_0^x p(\rho)\,d\rho$ is already quite accurate for $\rho_0= 3$, see the plot (blue = exact, orange = approximate).
