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).

<IMG SRC="https://i.sstatic.net/4rWiW.png" WIDTH="400"/>