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