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}(\rho)=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)=\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.

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

<sub>

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