For large $\rho_0$ you can use an asymptotic expansion of the Bessel function, which gives
$$p_{\rm appr}=\frac{1}{2\sqrt{\pi}}e^{-\rho^2-\rho_0^2}e^{2\rho\rho_0}.$$
Note that your distribution is normalised to $1/2$.

The approximation is already quite accurate for $\rho\gtrsim 3$, see the plot (blue = exact, orange = approximate).

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