[proofplan] The covariance sequence is first shown to be positive definite by writing its finite quadratic forms as variances of finite linear combinations of the centered process. The Herglotz representation theorem for positive-definite sequences on $\mathbb Z$ then gives the unique spectral measure on the torus. Under absolute summability, the Fourier series defining the proposed density converges uniformly; positivity follows by approximating it with nonnegative Fejer means coming from the same positive-definite quadratic forms. Finally, termwise integration recovers the covariance sequence, so uniqueness of the Herglotz measure identifies the spectral measure with this absolutely continuous measure. [/proofplan]

[step:Show that the autocovariance sequence is positive definite]Define the centered process $Y_t:\Omega\to\mathbb C$ by \begin{align*} Y_t(\omega):=X_t(\omega)-m, \qquad t\in\mathbb Z. \end{align*} For arbitrary $N\in\mathbb N$, arbitrary indices $t_1,\dots,t_N\in\mathbb Z$, and arbitrary coefficients $c_1,\dots,c_N\in\mathbb C$, define the square-integrable random variable $Z:\Omega\to\mathbb C$ by \begin{align*} Z(\omega):=\sum_{j=1}^{N} c_jY_{t_j}(\omega). \end{align*} Since $\mathbb E[|Y_t|^2]<\infty$ for every $t\in\mathbb Z$, the finite linear combination $Z$ belongs to $L^2(\Omega,\mathcal F,\mathbb P)$. Expanding the nonnegative number $\mathbb E[|Z|^2]$ gives \begin{align*} 0 \leq \mathbb E[|Z|^2] &=\mathbb E\left[\left(\sum_{j=1}^{N} c_jY_{t_j}\right) \overline{\left(\sum_{k=1}^{N} c_kY_{t_k}\right)}\right] \\ &=\sum_{j=1}^{N}\sum_{k=1}^{N} c_j\overline{c_k}\, \mathbb E[Y_{t_j}\overline{Y_{t_k}}] \\ &=\sum_{j=1}^{N}\sum_{k=1}^{N} c_j\overline{c_k}\,\gamma(t_j-t_k). \end{align*} Thus $\gamma:\mathbb Z\to\mathbb C$ is positive definite.[/step]

[guided]We must connect the probabilistic object, the covariance function, to the harmonic-analysis input used later. The required input is positive definiteness of the sequence $\gamma$. Define the centered process $Y_t:\Omega\to\mathbb C$ by \begin{align*} Y_t(\omega):=X_t(\omega)-m, \qquad t\in\mathbb Z. \end{align*} Since $(X_t)_{t\in\mathbb Z}$ is weakly stationary and $\mathbb E[|X_0|^2]<\infty$, every $Y_t$ is square-integrable. Now fix $N\in\mathbb N$, indices $t_1,\dots,t_N\in\mathbb Z$, and coefficients $c_1,\dots,c_N\in\mathbb C$. Define \begin{align*} Z:\Omega&\to\mathbb C\\ \omega&\mapsto \sum_{j=1}^{N} c_jY_{t_j}(\omega). \end{align*} This random variable belongs to $L^2(\Omega,\mathcal F,\mathbb P)$ because it is a finite linear combination of square-integrable random variables. Therefore $\mathbb E[|Z|^2]\geq 0$. Expanding the square and using linearity of expectation gives \begin{align*} 0 \leq \mathbb E[|Z|^2] &=\mathbb E\left[\left(\sum_{j=1}^{N} c_jY_{t_j}\right) \overline{\left(\sum_{k=1}^{N} c_kY_{t_k}\right)}\right] \\ &=\sum_{j=1}^{N}\sum_{k=1}^{N} c_j\overline{c_k}\, \mathbb E[Y_{t_j}\overline{Y_{t_k}}]. \end{align*} By the definition of the autocovariance function, \begin{align*} \mathbb E[Y_{t_j}\overline{Y_{t_k}}]=\gamma(t_j-t_k), \end{align*} because $\gamma(h)=\mathbb E[Y_{t+h}\overline{Y_t}]$. Hence \begin{align*} \sum_{j=1}^{N}\sum_{k=1}^{N} c_j\overline{c_k}\,\gamma(t_j-t_k)\geq 0. \end{align*} This is exactly positive definiteness of the sequence $\gamma$ on $\mathbb Z$.[/guided]

[step:Apply Herglotz representation to obtain the spectral measure] The sequence $\gamma:\mathbb Z\to\mathbb C$ is positive definite by the previous step, and its value at zero is finite because \begin{align*} \gamma(0)=\mathbb E[|X_0-m|^2]<\infty. \end{align*} By the Herglotz representation theorem for positive-definite sequences on $\mathbb Z$, there exists a unique finite positive Borel measure $\mu$ on $\mathbb T=\mathbb R/(2\pi\mathbb Z)$ such that \begin{align*} \gamma(h)=\int_{[-\pi,\pi]} e^{ih\lambda}\,d\mu(\lambda), \qquad h\in\mathbb Z. \end{align*} This proves the existence and uniqueness part of the theorem. [/step]

[step:Construct the Fourier series density under absolute summability]Assume \begin{align*} \sum_{h\in\mathbb Z}|\gamma(h)|<\infty. \end{align*} For each $h\in\mathbb Z$, define $\phi_h:[-\pi,\pi]\to\mathbb C$ by \begin{align*} \phi_h(\lambda):=\gamma(h)e^{-ih\lambda}. \end{align*} Since $|\phi_h(\lambda)|=|\gamma(h)|$ for every $\lambda\in[-\pi,\pi]$, the Weierstrass uniform convergence criterion applies to the series $\sum_{h\in\mathbb Z}\phi_h$. Hence the function \begin{align*} f:[-\pi,\pi]&\to\mathbb C\\ \lambda&\mapsto \frac{1}{2\pi}\sum_{h=-\infty}^{\infty}\gamma(h)e^{-ih\lambda} \end{align*} is well-defined and continuous. For $N\in\mathbb N$, define the Fejer mean $f_N:[-\pi,\pi]\to\mathbb C$ by \begin{align*} f_N(\lambda):=\frac{1}{2\pi}\sum_{|h|<N}\left(1-\frac{|h|}{N}\right)\gamma(h)e^{-ih\lambda}. \end{align*} For every $\lambda\in[-\pi,\pi]$, \begin{align*} f_N(\lambda) &=\frac{1}{2\pi N}\sum_{j=0}^{N-1}\sum_{k=0}^{N-1}\gamma(j-k)e^{-i(j-k)\lambda}. \end{align*} Positive definiteness of $\gamma$, applied with $t_j=j$ and $c_j=e^{-ij\lambda}$ for $j=0,\dots,N-1$, gives $f_N(\lambda)\geq 0$. Since $\sum_{h\in\mathbb Z}|\gamma(h)|<\infty$, dominated convergence for series gives $f_N\to f$ uniformly on $[-\pi,\pi]$. Therefore $f(\lambda)\geq 0$ for every $\lambda\in[-\pi,\pi]$, and $f$ is real-valued.[/step]

[guided]The displayed formula for $f$ is an infinite Fourier series, so the first task is to justify that it defines an honest function. For each integer $h$, define \begin{align*} \phi_h:[-\pi,\pi]&\to\mathbb C\\ \lambda&\mapsto \gamma(h)e^{-ih\lambda}. \end{align*} The bound \begin{align*} |\phi_h(\lambda)|=|\gamma(h)| \end{align*} holds for every $\lambda\in[-\pi,\pi]$. Because $\sum_{h\in\mathbb Z}|\gamma(h)|<\infty$, the Weierstrass uniform convergence criterion implies that $\sum_{h\in\mathbb Z}\phi_h$ converges uniformly. Thus \begin{align*} f:[-\pi,\pi]&\to\mathbb C\\ \lambda&\mapsto \frac{1}{2\pi}\sum_{h=-\infty}^{\infty}\gamma(h)e^{-ih\lambda} \end{align*} is well-defined and continuous, since it is the uniform limit of continuous partial sums. It remains to prove that this continuous function is nonnegative. The raw symmetric partial sums need not be nonnegative, so we use Fejer means. For $N\in\mathbb N$, define \begin{align*} f_N:[-\pi,\pi]&\to\mathbb C\\ \lambda&\mapsto \frac{1}{2\pi}\sum_{|h|<N}\left(1-\frac{|h|}{N}\right)\gamma(h)e^{-ih\lambda}. \end{align*} The weights count pairs $(j,k)$ with $j-k=h$, so \begin{align*} \sum_{j=0}^{N-1}\sum_{k=0}^{N-1}\gamma(j-k)e^{-i(j-k)\lambda} = \sum_{|h|<N}(N-|h|)\gamma(h)e^{-ih\lambda}. \end{align*} Therefore \begin{align*} f_N(\lambda) &=\frac{1}{2\pi N}\sum_{j=0}^{N-1}\sum_{k=0}^{N-1}\gamma(j-k)e^{-i(j-k)\lambda}. \end{align*} Now apply positive definiteness of $\gamma$ with indices $t_j=j$ and coefficients $c_j=e^{-ij\lambda}$ for $j=0,\dots,N-1$. It gives \begin{align*} \sum_{j=0}^{N-1}\sum_{k=0}^{N-1}c_j\overline{c_k}\gamma(j-k) = \sum_{j=0}^{N-1}\sum_{k=0}^{N-1}\gamma(j-k)e^{-i(j-k)\lambda} \geq 0. \end{align*} Hence $f_N(\lambda)\geq 0$ for every $\lambda\in[-\pi,\pi]$. Finally, the absolute summability of $\gamma$ implies \begin{align*} \sup_{\lambda\in[-\pi,\pi]}|f_N(\lambda)-f(\lambda)| \leq \frac{1}{2\pi}\sum_{h\in\mathbb Z} \left|a_N(h)-1\right|\,|\gamma(h)|, \end{align*} where $a_N:\mathbb Z\to[0,1]$ is the coefficient map defined by \begin{align*} a_N(h):= \begin{cases} 1-\frac{|h|}{N}, & |h|<N,\\ 0, & |h|\geq N. \end{cases} \end{align*} For each fixed $h\in\mathbb Z$, $a_N(h)\to 1$, and $|a_N(h)-1||\gamma(h)|\leq |\gamma(h)|$. Since $\sum_{h\in\mathbb Z}|\gamma(h)|<\infty$, dominated convergence for series gives uniform convergence $f_N\to f$. A uniform limit of nonnegative real-valued functions is nonnegative and real-valued, so $f(\lambda)\geq 0$ for every $\lambda\in[-\pi,\pi]$.[/guided]

[step:Recover the covariance sequence from the density] Define a finite positive Borel measure $\nu$ on $[-\pi,\pi]$ by \begin{align*} \nu(A):=\int_A f(\lambda)\,d\mathcal L^1(\lambda), \qquad A\in\mathcal B([-\pi,\pi]). \end{align*} For every $h\in\mathbb Z$, uniform convergence of the defining Fourier series for $f$ permits termwise integration: \begin{align*} \int_{[-\pi,\pi]} e^{ih\lambda}\,d\nu(\lambda) &=\int_{[-\pi,\pi]} e^{ih\lambda}f(\lambda)\,d\mathcal L^1(\lambda)\\ &=\frac{1}{2\pi}\sum_{k\in\mathbb Z}\gamma(k) \int_{[-\pi,\pi]} e^{i(h-k)\lambda}\,d\mathcal L^1(\lambda). \end{align*} The elementary orthogonality relation for exponentials on $[-\pi,\pi]$ gives \begin{align*} \int_{[-\pi,\pi]} e^{i(h-k)\lambda}\,d\mathcal L^1(\lambda) = \begin{cases} 2\pi, & k=h,\\ 0, & k\neq h. \end{cases} \end{align*} Consequently, \begin{align*} \int_{[-\pi,\pi]} e^{ih\lambda}\,d\nu(\lambda)=\gamma(h), \qquad h\in\mathbb Z. \end{align*} The measure $\nu$ and the Herglotz spectral measure $\mu$ have the same Fourier coefficients. By the uniqueness part of the Herglotz representation theorem on $\mathbb T$, $\mu=\nu$. Therefore $\mu$ is absolutely continuous with respect to $\mathcal L^1$, with density $f$. [/step]