The following is taken/paraphrased from page 18 of Geometry of Sets and Measures in Euclidean Spaces by Mattila.
Theorem. Let $\{\psi_\epsilon\}_{\epsilon>0}$ be an approximate identity and $\mu$ a Radon measure on $\mathbb{R}^n$. Then the functions $\psi_\epsilon * \mu$ converge weakly to $\mu$ as $\epsilon \to 0$, that is, $$\lim_{\epsilon \to 0} \int\varphi(\psi_\epsilon * \mu) \text{d}\mathcal{L}^n=\int \varphi \text{d}\mu$$ for all $\varphi \in C_0(\mathbb{R}^n)$.
Proof. We use Fubini's theorem, change of variable and the facts that spt $\psi_\epsilon \subset B(\epsilon)$ and $\int \psi_\epsilon \text{d}\mathcal{L}^n=1$ to compute $$\begin{equation} \begin{aligned} \int\varphi(\psi_\epsilon * \mu) \text{d}\mathcal{L}^n - \int \varphi \text{d}\mu &= \int \varphi(x) \int \psi_\epsilon(x-y) \text{d}\mu y\text{d}x - \int \varphi(y) \int \psi_\epsilon(x) \text{d}x\text{d}\mu y \textbf{ (1)} \\ &= \cdots \text{ (Applying Fubini's theorem and arithmetic)} \\ &= \iint_{B(\epsilon)}[\varphi(x+y)-\varphi(y)]\psi_\epsilon(x)\text{d}x\text{d}\mu y. \end{aligned} \end{equation}$$ Then, as $\varphi$ is uniformly continuous, and $\int \psi_\epsilon \text{d}\mathcal{L}^n=1$, this goes to zero, so we have our result.
I am concerned about the equality marked $\mathbf{(1)}$ above. Why is it that the composition of the convolution and $\varphi$ ends up just looking like the product of the two integrals instead?
[Edit: My question here is why don't we have something like
$$\int\varphi(\psi_\epsilon * \mu) \text{d}\mathcal{L}^n - \int \varphi \text{d}\mu = \int \varphi(\int \psi_\epsilon(x-y)) - \int \varphi(y) \int \psi_\epsilon(x)$$
instead in step $\textbf{(1)}$? This might just be a basic arithmetic mistake on my part.]
Also, I don't really understand the deeper meaning of this theorem. What is its meaning?