[proofplan]
We test the autocovariance kernel against an arbitrary finite vector of real coefficients. The key construction is the centered linear combination $Y=\sum_{i=1}^n a_i(X_i-\mu)$, whose square has non-negative expectation. Expanding $\mathbb{E}[Y^2]$ and using weak stationarity identifies that expectation with the required quadratic form.
[/proofplan]
[step:Form the centered linear combination associated to the coefficient vector]
Fix $n\in\mathbb{N}$ and fix real numbers $a_1,\dots,a_n\in\mathbb{R}$. Let $(\Omega,\mathcal{F},\mathbb{P})$ denote the probability space on which the process $(X_t)_{t\in\mathbb{Z}}$ is defined, let $\mathcal{B}(\mathbb{R})$ denote the Borel $\sigma$-algebra on $\mathbb{R}$, and let $L^2(\Omega,\mathcal{F},\mathbb{P})$ denote the space of real-valued square-integrable random variables on this probability space, modulo equality $\mathbb{P}$-almost surely. Let $\mu\in\mathbb{R}$ denote the common mean $\mu=\mathbb{E}[X_t]$, which is independent of $t$ by weak stationarity. For each $i\in\{1,\dots,n\}$, define the centered random variable
\begin{align*}
Z_i:(\Omega,\mathcal{F})&\to(\mathbb{R},\mathcal{B}(\mathbb{R}))\\
\omega&\mapsto X_i(\omega)-\mu.
\end{align*}
Since $X_i$ has finite second moment, $Z_i\in L^2(\Omega,\mathcal{F},\mathbb{P})$. Define the finite linear combination
\begin{align*}
Y:(\Omega,\mathcal{F})&\to(\mathbb{R},\mathcal{B}(\mathbb{R}))\\
\omega&\mapsto \sum_{i=1}^{n}a_i Z_i(\omega).
\end{align*}
Because $L^2(\Omega,\mathcal{F},\mathbb{P})$ is closed under finite linear combinations, $Y\in L^2(\Omega,\mathcal{F},\mathbb{P})$.
[/step]
[step:Use the non-negativity of the second moment]
Since $Y^2\ge 0$ pointwise on $\Omega$ and $Y\in L^2(\Omega,\mathcal{F},\mathbb{P})$, the expectation $\mathbb{E}[Y^2]$ is finite and satisfies
\begin{align*}
\mathbb{E}[Y^2]\ge 0.
\end{align*}
[/step]
[step:Expand the square into the covariance quadratic form]
Using finite distributivity and linearity of expectation, we compute
\begin{align*}
\mathbb{E}[Y^2]
&=\mathbb{E}\left[\left(\sum_{i=1}^{n}a_i Z_i\right)\left(\sum_{j=1}^{n}a_j Z_j\right)\right]\\
&=\mathbb{E}\left[\sum_{i=1}^{n}\sum_{j=1}^{n}a_i a_j Z_i Z_j\right]\\
&=\sum_{i=1}^{n}\sum_{j=1}^{n}a_i a_j\,\mathbb{E}[Z_i Z_j].
\end{align*}
For each pair $i,j\in\{1,\dots,n\}$, the definition of $Z_i$ and $Z_j$ gives
\begin{align*}
\mathbb{E}[Z_i Z_j]
=\mathbb{E}[(X_i-\mu)(X_j-\mu)].
\end{align*}
[guided]
The point of introducing $Y$ is that its square automatically produces all pairwise covariance terms. We first expand the square as a finite algebraic identity:
\begin{align*}
Y^2
&=\left(\sum_{i=1}^{n}a_i Z_i\right)\left(\sum_{j=1}^{n}a_j Z_j\right)\\
&=\sum_{i=1}^{n}\sum_{j=1}^{n}a_i a_j Z_i Z_j.
\end{align*}
There is no convergence issue here because both sums are finite. Since each $Z_i\in L^2(\Omega,\mathcal{F},\mathbb{P})$, the products $Z_iZ_j$ are integrable by the Cauchy-Schwarz inequality for random variables:
\begin{align*}
\mathbb{E}[|Z_iZ_j|]\le \mathbb{E}[Z_i^2]^{1/2}\mathbb{E}[Z_j^2]^{1/2}<\infty.
\end{align*}
Therefore linearity of expectation applies to the finite sum, giving
\begin{align*}
\mathbb{E}[Y^2]
&=\mathbb{E}\left[\sum_{i=1}^{n}\sum_{j=1}^{n}a_i a_j Z_i Z_j\right]\\
&=\sum_{i=1}^{n}\sum_{j=1}^{n}a_i a_j\,\mathbb{E}[Z_i Z_j].
\end{align*}
Finally, because $Z_i=X_i-\mu$ and $Z_j=X_j-\mu$, each expectation in the last expression is exactly
\begin{align*}
\mathbb{E}[Z_i Z_j]
=\mathbb{E}[(X_i-\mu)(X_j-\mu)].
\end{align*}
[/guided]
[/step]
[step:Identify each covariance term with the autocovariance function]
By weak stationarity, the covariance between $X_i$ and $X_j$ depends only on the lag $i-j$. Since $\gamma(h)=\mathbb{E}[(X_{t+h}-\mu)(X_t-\mu)]$ for every $h,t\in\mathbb{Z}$, taking $h=i-j$ and $t=j$ gives
\begin{align*}
\mathbb{E}[(X_i-\mu)(X_j-\mu)]
=\mathbb{E}[(X_{j+(i-j)}-\mu)(X_j-\mu)]
=\gamma(i-j).
\end{align*}
Substituting this identity into the expansion of $\mathbb{E}[Y^2]$ yields
\begin{align*}
\mathbb{E}[Y^2]
=\sum_{i=1}^{n}\sum_{j=1}^{n}a_i a_j\gamma(i-j).
\end{align*}
Combining this equality with $\mathbb{E}[Y^2]\ge0$ proves
\begin{align*}
\sum_{i=1}^{n}\sum_{j=1}^{n}a_i a_j\gamma(i-j)\ge0.
\end{align*}
This is the desired positive definiteness of the autocovariance function.
[/step]