Bochner's Theorem: If $\varphi : \mathbb{R}^d \to \mathbb{C}$ is positive definite, continuos and $\varphi(0)=1$ then it is the characteristic function of a probability measure, i.e. the Bochner's theorem.
Proof:
We can prove that if $f \in L^1(\mathbb{R}^d,\mathbb{C})$ then $\int_{\mathbb{R}^d \times \mathbb{R}^d}\varphi(t-s)f(t)\overline{f(s)}dtds \geq 0$ using the fact that $\varphi$ is positive definite.
Thus if $f(t)=e^{i \langle t,x \rangle-\varepsilon \|t\|^2}$ we can get that $h_{\varepsilon}(x)=\int_{\mathbb{R}^d}\varphi(u)e^{i \langle u,x \rangle -\frac {\varepsilon}{2}\|u\|^2}du \geq 0$.
If I could compute $\int_{\mathbb{R}^d}h_{\varepsilon}(x)dx$ I would modify $h_{\varepsilon}(x)$ in such a way it is a probability density.
Then I would compute the characteristic function of a random variable with this density and I will find that it is $\varphi_{\varepsilon}(u)=\varphi(u)e^{\frac{-\|t\|^2\varepsilon}{2}}$ and I will conclude thanks to Lévy's theorem with $\varepsilon \to 0$ that $\varepsilon$ is a characteristic fucntion of some probability.
My question is how can we evaluate $\int_{\mathbb{R}^d}h_{\varepsilon}(x)dx$?
It is true that $0\leq \int_{\mathbb{R}^d \times \mathbb{R}^d}\varphi(t-s)f(t)\overline{f(s)}dtds =\int_{\mathbb{R}^d}\varphi(u)e^{i \langle u,x \rangle - \frac{\varepsilon \|u\|^2}{2}}(\frac{\pi}{2\varepsilon})^{\frac d 2 }du$?
Below I will use the (inverse) Fourier transform - there is the never-ending question about the constants in the definition of the (inverse) Fourier transform, so let me fix the definitions which I am using: \begin{align*} (\mathcal{F})f(\xi) &:= \frac{1}{(2\pi)^d} \int_{\mathbb{R}^d} f(x) e^{-i \langle x,\xi \rangle} \, dx, \\ (\mathcal{F}^{-1} f)(\xi) &:= \int_{\mathbb{R}^d} f(x) e^{i \langle x,\xi \rangle} \, dx. \end{align*}
Set $\varphi_{\epsilon}(u) := \varphi(u) \exp \left(- \frac{\epsilon}{2} \|u\|^2 \right)$, then $$h_{\epsilon}(x) = \int \varphi_{\epsilon}(u) e^{i \langle u,x \rangle} \, du = (\mathcal{F}^{-1} \varphi_{\epsilon})(x). \tag{1}$$ Denote by
$$p_{t}(x) := \frac{1}{(2\pi t)^{d/2}} \exp \left(- \frac{|x|^2}{2t} \right)$$ the density of the Gaussian distribution with mean zero and covariance matrix $t \cdot \text{Id}_{d \times d}$. Note that the Fourier transform $\mathcal{F}p_t$ is given by $$(\mathcal{F}p_t)(x) = \frac{1}{(2\pi)^d} \exp \left(- t \frac{|x|^2}{2} \right). \tag{2}$$ Thus,
\begin{align*} \int h_{\epsilon}(x) e^{-\frac{t}{2} |x|^2} \, dx &\stackrel{(1)}{=} \int (\mathcal{F}^{-1} \phi_{\epsilon})(x) e^{-\frac{t}{2} |x|^2} \, dx \\ &= \int \phi_{\epsilon}(\xi) (\mathcal{F}^{-1} e^{-\frac{t}{2} |\bullet|^2})(\xi) \, d\xi \\ &\stackrel{(2)}{=} (2\pi)^d \int \phi_{\epsilon}(\xi) p_t(\xi) \, d\xi. \end{align*}
Since $\int_{\mathbb{R}^d} p_t(\xi) \, d\xi=1$ and $|\phi(\xi)| \leq \phi(0)$ (because $\phi$ is positive definite), we get
$$\int h_{\epsilon}(x) e^{-\frac{t}{2} |x|^2} \, dx \leq (2\pi)^d \phi(0).$$
Thus, by Fatou's lemma,
$$\int h_{\epsilon}(x) \, dx \leq (2\pi)^d \phi(0).$$
Since $h_{\epsilon} \geq 0$, this shows that $h_{\epsilon} \in L^1$ and $\|h_{\epsilon}\|_{L^1} \leq (2\pi)^d \phi(0)$.
If we define $\mu_{\epsilon}(dx) := h_{\epsilon}(x) \, dx$, then the so-defined (non-negative) measure is finite and, by (1), $\mathcal{F}\mu_{\epsilon} = \phi_{\epsilon} \to \phi$ as $\epsilon \to 0$. Applying Lévy's continuity theorem, it follows that there exists a finite measure $\mu$ such that $\mu_n \to \mu$ weakly and $\mathcal{F}\mu=\phi$. Finally, $$\frac{1}{(2\pi)^d} \stackrel{**}{=} \varphi(0) = (\mathcal{F}\mu)(0)=\frac{1}{(2\pi)^d} \mu(\mathbb{R}^d)$$ shows that $\mu$ is a probability measure. (Remark on $(**)$: I'm using a different normalization for the Fourier transform than you do; for the normalizing I am using the Fourier transform $\mathcal{F}\nu$ of a probability measure $\nu$ on $\mathbb{R}^d$ satisfies $(\mathcal{F}\nu)(0)=\frac{1}{(2\pi)^d}$; therefore the assumption $\phi(0)=1$ needs to be replaced by $\phi(0)=\frac{1}{(2\pi)^d}$. Sorry about that.)