You wish to approximate $$p_{\rm exact}(\rho)=e^{-\rho^2-\rho_0^2} \rho I_0(2\rho \rho_0).$$
For small $\rho_0$ a Taylor series in powers of $\rho_0$ is accurate, $$p_{\rm small}=e^{-\rho^2} \rho (1+\rho^2 \rho_0^2-\rho_0^2).$$ For large $\rho_0$ an asymptotic expansion of the Bessel function gives $$p_{\rm large}(\rho_0)=\frac{1}{2\sqrt{\pi}}e^{-\rho^2-\rho_0^2}e^{2\rho\rho_0}.$$
In the plots below I compare $\int_0^x p(\rho)\,d\rho$ for the exact expression (blue) and the approximation (orange). It turns out the approximations $p_{\rm small}$ and $p_{\rm large}$ are already quite accurate for $\rho_0\lesssim 0.5$ and $\rho_0\gtrsim 3$, respectively.
Left plot: comparison of $\int_0^x p_{\rm exact}(\rho)\,d\rho$ and $\int_0^x p_{\rm small}(\rho)\,d\rho$ for $\rho_0=0.5$. Right plot: comparison of $\int_0^x p_{\rm exact}(\rho)\,d\rho$ and $\int_0^x p_{\rm large}(\rho)\,d\rho$ for $\rho_0=3$. The approximations (orange) are almost indistinguishable from the exact answer (blue).