I would like to show that the Wasserstein distance between a probability measure $\mu$ (in $\mathbb{R}^d$) and a mollified version of it $\mu_{\epsilon}:=\mu\ast\rho_\epsilon$ where $\mu\ast\rho_\epsilon(x):=\int \rho_{\epsilon}(x-y)d\mu(y)$ goes to 0, i.e. $$\lim\limits_{\epsilon\to 0}W_p(\mu,\mu_{\epsilon}).$$ I was told that to do this I can show that $$W_p(\mu,\mu_{\epsilon})\le\epsilon m_p(\rho)$$ where $m_p(\rho):=\left(\int d(x,x)^p\rho(x)dx\right)^{1/p}$ denotes the $p$-th moment. I think I have to come up with a particular transport plan between $\mu$ and $\nu$. I tried with the most obvious one to me, i.e. the product measure $\mu\times\mu_{\epsilon}$ but I didn't arrive to the desired conclusion.
edit
Using the plan $\mu\times\mu_{\epsilon}$ I get \begin{align} W_p^p(\mu,\mu_{\epsilon})&\le \int d(x,y) d\mu(x)d\mu_\epsilon(y)\\ &= \int d(x,y)\mu\ast\rho_{\epsilon}(y)d\mu(x) dy\\ &= \int d(x,y)\left(\int\rho_{\epsilon}(y-z)d\mu(z)\right)d\mu(x) dy \end{align} But now I a am stuck
The problem with using the product measure $\mu \otimes \mu_{\epsilon}$ is it does not use the form of $\mu_{\epsilon}$ at all. Convolutions have a very specific probabilistic interpretation, as do couplings, and the $p$-th moment you wrote also frequently appears in probabilistic computations. Hence a second guess is to exploit this connection.
Recall that if $Y$ is a random vector with distribution $\mu$ and $X_{\epsilon}$ is a random variable with distribution $\rho_{\epsilon}(x) \, dx$, and if $Y$ and $X_{\epsilon}$ are independent, then $X_{\epsilon} + Y$ has distribution $\rho_{\epsilon} * \mu = \mu_{\epsilon}$. Accordingly, if $\nu$ is the distribution of the pair $(X_{\epsilon} + Y,Y)$, then the first marginal $\pi_{1*} \nu$ equals $\mu_{\epsilon}$ and the second marginal $\pi_{2*}\nu$ equals $\mu$. In particular, $\nu$ is a coupling of $\mu_{\epsilon}$ and $\mu$.
Also, I'm assuming that $d(x,y) = \|x - y\|$, where $\|\cdot\|$ is the Euclidean norm. However, all that is really necessary is that $d$ is translationally invariant (i.e. $d(x + z, y + z) = d(x,y)$ for all $x,y, z \in \mathbb{R}^{d}$). In particular, other norms would also work.
With the coupling $\nu$ in hand, let's try to bound $W_{p}(\mu_{\epsilon},\mu)$. I will use the random variables to do the computation since it is much clearer. If that's not your style, it would be a good exercise to write down a purely analytic proof. We find: \begin{align*} W_{p}(\mu_{\epsilon},\mu)^{p} &\leq \int_{\mathbb{R}^{d} \times \mathbb{R}^{d}} \|x - y\|^{p} \, \nu(dx \otimes dy) \\ &= \mathbb{E}(\|(X_{\epsilon} +Y) - Y\|^{p}) \\ &= \mathbb{E}(\|X_{\epsilon}\|^{p}) \\ &= \int_{\mathbb{R}^{d}} \|x\|^{p} \rho_{\epsilon}(x) \, dx \\ &= \epsilon^{p} \int_{\mathbb{R}^{d}} \|\epsilon^{-1} x\|^{p} \rho(\epsilon^{-1}x) \, \frac{dx}{\epsilon^{d}} \\ &= \epsilon^{p} M_{p}(\rho), \end{align*} where the $p$th moment of $\rho$ is $M_{p}(\rho) = \int_{\mathbb{R}^{d}} \|y\|^{p} \rho(y) \, dy$.
It is worth noting (and can be used to simplify the computation above somewhat) that $X_{\epsilon}$ has the same distribution as $\epsilon X$, where $X$ has distribution $\rho(x) \, dx$. In this sense, what we are doing in the last three lines above is simply computing $\mathbb{E}(\|\epsilon X\|^{p}) = \epsilon^{p} \mathbb{E}(\|X\|^{p})$ --- another indication that the probabilistic interpretation is very convenient here.
Hint to "un-probability" the computation above: Let $\nu = \nu_{\epsilon}$ be the probability measure on $\mathbb{R}^{d} \times \mathbb{R}^{d}$ given by its action on $f \in C_{c}(\mathbb{R}^{d} \times \mathbb{R}^{d})$ by \begin{equation*} \int_{\mathbb{R}^{d} \times \mathbb{R}^{d}} f(x,y) \nu_{\epsilon}(dx \otimes dy) = \int_{\mathbb{R}^{d}} \left( \int_{\mathbb{R}^{d}} f(y' + x',y') \rho_{\epsilon}(x') \, dx' \right) \mu(dy'). \end{equation*} Check that $\nu_{\epsilon}$ as defined is really a coupling of $\mu$ and $\mu_{\epsilon}$. It is also worth verifying that $\nu_{\epsilon}$ can be rewritten as \begin{equation*} \int_{\mathbb{R}^{d} \times \mathbb{R}^{d}} f(x,y) \nu_{\epsilon}(dx \otimes dy) = \int_{\mathbb{R}^{d}} \left( \int_{\mathbb{R}^{d}} f(y' + \epsilon x', y') \rho(x') \, dx' \right) \mu(dy'). \end{equation*} Hence $\nu_{\epsilon}$ is looking very close to the diagonal. Finally, check that \begin{equation*} \int_{\mathbb{R}^{d} \times \mathbb{R}^{d}} \|x-y\|^{p} \nu_{\epsilon}(dx \otimes dy) = \int_{\mathbb{R}^{d}} \|x'\|^{p} \rho_{\epsilon}(x') \, dx' \end{equation*} and complete the argument as in the answer above.
Edit: The product measure will not give a useful estimate. To see this, consider the case when $\mu$ is Gaussian, that is, $\mu(dx) = (2 \pi)^{-d/2}\exp(-\|x\|^{2}/2) \, dx$ and let $\rho$ also be Gaussian: \begin{equation*} \rho(x) = (2 \pi)^{-d/2} \exp(-\|x\|^{2}/2). \end{equation*} One can then check that $\mu \otimes \mu_{\epsilon}$ is nothing other than two independent Gaussians, the second with variance $1 + \epsilon^{2}$: \begin{equation*} (\mu \otimes \mu_{\epsilon})(dx \otimes dy) = (2 \pi)^{-d} \epsilon^{-d} \exp \left(- \|x\|^{2}/2 - \|y\|^{2}/2(1+\epsilon^{2})\right) \, dx dy. \end{equation*} Now in the easy case when $p = 2$, writing $\|x - y\|^{2} = \|x\|^{2} + \|y\|^{2} - 2 \langle x,y \rangle$ and eliminating the cross term leads to \begin{equation*} \int_{\mathbb{R}^{d} \times \mathbb{R}^{d}} \|x - y\|^{2} (\mu \otimes \mu_{\epsilon})(dx \otimes dy) = \int_{\mathbb{R}^{d}} \|x\|^{2} \mu(dx) + \int_{\mathbb{R}^{d}} \|y\|^{2} \mu_{\epsilon}(dy) = 2 + \epsilon^{2}. \end{equation*}