Here's one way to do it without getting your hands too dirty. It's not constructive.
Let ${\cal A}$ be the union of a countable dense set in $C_0({\mathbb R}^n)$, and indicators of all balls of integer radius centered at the origin (not necessary, but helpful later).
Let $\{X_1,X_2,\dots \}$ be an IID $\mu$-distributed sequence.
Let $\mu_n = \frac{1}{n} \sum_{j\le n} \delta_{X_j}.$
Then $\mu_n$ is a random probability measure. Observe that for any bounded and measurable function,
$$\int f d\mu_n= \frac {1}{n} \sum_{j\le n } f(X_j).$$
Then by LLN
$$ \lim_{n\to\infty} \int f d\mu_n = \int f d\mu ,~\forall f \in {\cal A}\quad (*),$$
a.s. (here is where we use the countability of ${\cal A}$).
Let $\{\bar \mu_n\}$ be a realization of $\{\mu_n\}$ (or equivalently, $\{X_1,\dots\}$) for which $(*)$ holds. It remains to show that $\{\bar \mu_n\}$ converges weakly.
Let $f\in C_b({\mathbb R}^d)$, by a standard procedure, for every $M$, there exists $f_M\in C_0({\mathbb R}^d)$ such that $f_M$ coincides with $f$ on $\{x:|x|<M\}$ and $|f_M |\le |f|$. Fix $\epsilon>0$ and choose an integer $M$ such that $\mu(\{x:|x|>M\})<\epsilon$
There exists a continuous $g_M\in {\cal A}$ such that $\|f_M-g_M\|<\epsilon$.
We then have
$$ \int f d \bar \mu_n = \int g_M d \bar \mu_n + \int (f_M -g_M) d \bar \mu_n +
\int (f- f_M) d \bar \mu_n.$$
The first integral on the RHS converges to $\int g_M d\mu$ by construction. The second is bounded in absolute value by $\epsilon$. The third is bounded by $2\|f\|\mu_n (\{x:|x|>M\})$, which by assumption (recall the indicator functions in ${\cal A}$) converges to $2\|f\|\mu(\{x:|x|>M\})=2\|f\|\epsilon$.
Therefore
$$ \limsup_{n\to\infty} |\int f d\bar \mu_n - \int g_M d\mu| \le \epsilon (1+2\|f\|).$$
But
$$\begin{align*} |\int g_M d\mu - \int f d \mu| &\le \int |g_M -f_M| d \mu + \int |f_M -f| d\mu\\
& \le \epsilon + 2\|f\|\mu (\{x:|x|>M\})\le \epsilon (1+2\|f\|)\end{align*}.$$
The result follows.
Best Answer
The measure $\rho$ admits a density $f\colon X^2\to\mathbb R$ with respect to the measure $\mu\times \delta_0$. Letting $$ \nu(A):= \int_X f\left(x,x_0\right)d\mu(x), $$ we can show that the measures $\nu\times\delta_{x_0}$ and $\rho$ coincide on finite disjoint unions of rectangles and since all the involved measures are finite, we derive that $\nu\times \delta_{x_0}=\rho$.