[proofplan] We use the **method of cumulants** (Brillinger's approach). Writing each ordinate in terms of its real and imaginary parts, we assemble the real random vector of all these parts and show that its joint cumulants converge to those of a centred Gaussian: the first cumulants vanish (zero mean), the second cumulants converge to a diagonal covariance matrix whose blocks are $\tfrac{1}{2}f(\omega_r)I_2$, and every joint cumulant of order $\ge 3$ tends to zero. The second-order analysis rests on summing geometric series of unimodular terms — bounded uniformly because the limiting frequencies and their pairwise sums and differences avoid $0\pmod{2\pi}$ on $(0,\pi)$ — together with the [Dominated Convergence Theorem](/theorems/4) on $\mathbb Z$; the higher-order analysis rests on the cumulant summability hypothesis, which forces an $n^{1-k/2}$ decay. Cumulant convergence yields convergence of all moments to the Gaussian moments, and since the multivariate normal is determined by its moments, the method of moments gives convergence in distribution; the [Continuous Mapping Theorem](/theorems/1847) then transfers the conclusion from the real coordinate vector to the complex ordinates, the diagonal covariance encoding both the independence of distinct ordinates and the independence of real and imaginary parts within each ordinate. [/proofplan]

[step:Set up the cumulant method and reduce to a real coordinate vector]Throughout, expectation is $\mathbb E[Y]=\int_\Omega Y\,d\mathbb P$ on the underlying probability space $(\Omega,\mathcal F,\mathbb P)$. We recall the notion of joint cumulant. For random variables $Y_1,\dots,Y_k$ with finite moments of all orders, the **$k$-th order joint cumulant** is \begin{align*} \operatorname{cum}(Y_1,\dots,Y_k)=\sum_{\pi\in\mathcal P(k)}(|\pi|-1)!\,(-1)^{|\pi|-1}\prod_{B\in\pi}\mathbb E\!\Big[\prod_{j\in B}Y_j\Big], \end{align*} where $\mathcal P(k)$ is the set of partitions of $\{1,\dots,k\}$ and $|\pi|$ is the number of blocks of $\pi$ (the Leonov–Shiryaev formula). In particular $\operatorname{cum}(Y_1)=\mathbb E[Y_1]$ and $\operatorname{cum}(Y_1,Y_2)=\mathbb E[Y_1Y_2]-\mathbb E[Y_1]\mathbb E[Y_2]$. Two structural facts are used repeatedly, both immediate from this formula: - **Multilinearity.** Each summand is a product of expectations of products of the $Y_j$, hence multilinear in $(Y_1,\dots,Y_k)$; the formula extends to complex linear combinations by $\mathbb C$-linearity, which is how we read cumulants of the complex-valued $d_n(\omega)$. - **Reality.** Since $X_t$ is real, $\overline{d_n(\omega)}=d_n(-\omega)$. For $r\in\{1,\dots,p\}$ define the real and imaginary parts of the ordinate at the perturbed frequency $\omega_{r,n}$, \begin{align*} A_{r,n}&=\operatorname{Re}d_n(\omega_{r,n})=\frac{1}{\sqrt{2\pi n}}\sum_{t=1}^{n}X_t\cos(\omega_{r,n} t),\\ B_{r,n}&=\operatorname{Im}d_n(\omega_{r,n})=-\frac{1}{\sqrt{2\pi n}}\sum_{t=1}^{n}X_t\sin(\omega_{r,n} t), \end{align*} so that $d_n(\omega_{r,n})=A_{r,n}-iB_{r,n}$, and assemble the real random vector \begin{align*} V_n=\big(A_{1,n},B_{1,n},\dots,A_{p,n},B_{p,n}\big)\in\mathbb R^{2p}. \end{align*} By the reality relation, $A_{r,n}=\tfrac12\big(d_n(\omega_{r,n})+d_n(-\omega_{r,n})\big)$ and $B_{r,n}=\tfrac{i}{2}\big(d_n(\omega_{r,n})-d_n(-\omega_{r,n})\big)$, so every coordinate of $V_n$ is a fixed complex-linear combination, with coefficients in $\{\pm\tfrac12,\pm\tfrac{i}{2}\}$, of ordinates $d_n(\lambda)$ with $\lambda\in\{\pm\omega_{r,n}\}$. Our goal is to compute the limiting joint cumulants of $V_n$: Step 2 handles orders $1$ and $2$, Step 3 handles orders $\ge 3$, and Step 4 concludes.[/step]

[guided]The strategy is the method of cumulants. Convergence in distribution of a sequence of random vectors can be established by showing that all of their cumulants (equivalently, all of their moments) converge to those of a target law, provided that target law is determined by its moments. The Gaussian is the prototypical moment-determinate law, and — crucially — it is the *unique* law whose cumulants of order $\ge 3$ all vanish. So the plan is: prove that $V_n$ has vanishing cumulants of order $\ge 3$ in the limit, and that its first two cumulants converge to those of a specific Gaussian. Why work with the real vector $V_n$ rather than directly with the complex ordinates? Convergence in distribution and the notion of a real Gaussian limit are cleanest for real random vectors; we will recover the complex statement at the very end through the continuous map $(a,b)\mapsto a-ib$. The two structural facts deserve emphasis. **Multilinearity** is what makes the entire method tractable: because each $d_n(\omega)$ is a *linear* functional of $(X_1,\dots,X_n)$, every joint cumulant of the $d_n$'s expands into a sum of joint cumulants of the $X_t$'s weighted by products of the unimodular factors $e^{-i\omega t}$. This converts a probabilistic question into the asymptotic analysis of exponential sums. **Reality**, $\overline{d_n(\omega)}=d_n(-\omega)$, holds because $X_t\in\mathbb R$ gives $\overline{X_t e^{-i\omega t}}=X_t e^{i\omega t}$; it lets us express the real and imaginary parts $A_{r,n},B_{r,n}$ as combinations of ordinates at $\pm\omega_{r,n}$, so that cumulants of the coordinates of $V_n$ reduce to cumulants of ordinates, which is the only object we will need to estimate.[/guided]

[step:Compute the second-order moments and identify the limiting covariance]Since $\mathbb E[X_t]=0$, multilinearity gives $\mathbb E[d_n(\omega)]=0$, hence every coordinate of $V_n$ has mean $0$ (all first cumulants vanish). For the second order, introduce the second-moment kernel \begin{align*} R_n(\alpha,\beta):=\mathbb E\!\big[d_n(\alpha)\,\overline{d_n(\beta)}\big]=\frac{1}{2\pi n}\sum_{s=1}^{n}\sum_{t=1}^{n}\gamma(s-t)\,e^{-i\alpha s+i\beta t}, \end{align*} where we used $\mathbb E[X_sX_t]=\gamma(s-t)$. Substituting $s=t+h$ and writing $e^{-i\alpha s+i\beta t}=e^{-i\alpha h}e^{-i(\alpha-\beta)t}$, \begin{align*} R_n(\alpha,\beta)=\frac{1}{2\pi n}\sum_{|h|<n}\gamma(h)\,e^{-i\alpha h}\sum_{t\in I_h}e^{-i(\alpha-\beta)t},\qquad I_h:=\{t: 1\le t\le n,\ 1\le t+h\le n\}, \end{align*} where $I_h$ is an integer interval with $|I_h|=n-|h|$. Note $\mathbb E[d_n(\alpha)d_n(\beta)]=R_n(\alpha,-\beta)$. **Diagonal case $\alpha=\beta=\omega_n\to\omega$.** Then $\alpha-\beta=0$, the inner sum equals $|I_h|=n-|h|$, and \begin{align*} R_n(\omega_n,\omega_n)=\frac{1}{2\pi}\sum_{|h|<n}\Big(1-\tfrac{|h|}{n}\Big)\gamma(h)\,e^{-i\omega_n h}. \end{align*} The summand, viewed as a function on $\mathbb Z$ with counting measure, is dominated by $|\gamma(h)|\in\ell^1(\mathbb Z)$ and converges pointwise to $\gamma(h)e^{-i\omega h}$ as $n\to\infty$ (since $1-|h|/n\to1$ and $e^{-i\omega_n h}\to e^{-i\omega h}$). By the [Dominated Convergence Theorem](/theorems/4) applied to the $\sigma$-finite counting measure on $\mathbb Z$, \begin{align*} R_n(\omega_n,\omega_n)\ \longrightarrow\ \frac{1}{2\pi}\sum_{h\in\mathbb Z}\gamma(h)e^{-i\omega h}=f(\omega). \end{align*} **Off-diagonal case $\alpha_n-\beta_n\to\rho\not\equiv 0\pmod{2\pi}$.** Then $e^{-i(\alpha_n-\beta_n)}\to e^{-i\rho}\ne 1$, so for all large $n$ we have $|1-e^{-i(\alpha_n-\beta_n)}|\ge\tfrac12|1-e^{-i\rho}|=:\delta>0$. Summing the geometric progression over the interval $I_h$, \begin{align*} \Big|\sum_{t\in I_h}e^{-i(\alpha_n-\beta_n)t}\Big|\le\frac{2}{|1-e^{-i(\alpha_n-\beta_n)}|}\le\frac{2}{\delta}\qquad(\text{all large }n), \end{align*} a bound independent of $h$ and $n$. Therefore \begin{align*} |R_n(\alpha_n,\beta_n)|\le\frac{1}{2\pi n}\cdot\frac{2}{\delta}\sum_{h\in\mathbb Z}|\gamma(h)|=O\!\big(n^{-1}\big)\ \longrightarrow\ 0. \end{align*} We now apply these two cases. All the relevant frequency combinations are nonzero modulo $2\pi$ because $\omega_1,\dots,\omega_p$ are distinct points of $(0,\pi)$: - **Variance.** $R_n(\omega_{r,n},\omega_{r,n})=\mathbb E\big[(A_{r,n})^2+(B_{r,n})^2\big]\to f(\omega_r)$. - **Complementary covariance.** $\mathbb E[d_n(\omega_{r,n})^2]=R_n(\omega_{r,n},-\omega_{r,n})$ has frequency difference $\alpha_n-\beta_n=2\omega_{r,n}\to 2\omega_r\in(0,2\pi)$, hence $\not\equiv 0\pmod{2\pi}$; so $\mathbb E[d_n(\omega_{r,n})^2]\to 0$. Expanding $d_n^2=(A-iB)^2=A^2-B^2-2iAB$, the real and imaginary parts give \begin{align*} \mathbb E\big[(A_{r,n})^2\big]-\mathbb E\big[(B_{r,n})^2\big]\to 0,\qquad \mathbb E\big[A_{r,n}B_{r,n}\big]\to 0. \end{align*} Together with the variance limit, this yields \begin{align*} \mathbb E\big[(A_{r,n})^2\big]\to\frac{f(\omega_r)}{2},\quad \mathbb E\big[(B_{r,n})^2\big]\to\frac{f(\omega_r)}{2},\quad \mathbb E\big[A_{r,n}B_{r,n}\big]\to 0. \end{align*} - **Cross-frequency covariances ($r\ne s$).** Both $R_n(\omega_{r,n},\omega_{s,n})$ (difference $\to\omega_r-\omega_s\in(-\pi,\pi)\setminus\{0\}$) and $R_n(\omega_{r,n},-\omega_{s,n})$ (difference $\to\omega_r+\omega_s\in(0,2\pi)$) have frequency differences $\not\equiv 0\pmod{2\pi}$, hence both tend to $0$. Expanding $\mathbb E[d_n(\omega_{r,n})\overline{d_n(\omega_{s,n})}]$ and $\mathbb E[d_n(\omega_{r,n})d_n(\omega_{s,n})]$ into real and imaginary parts as above shows that all four cross-covariances $\mathbb E[A_{r,n}A_{s,n}],\ \mathbb E[A_{r,n}B_{s,n}],\ \mathbb E[B_{r,n}A_{s,n}],\ \mathbb E[B_{r,n}B_{s,n}]$ tend to $0$. Hence the covariance matrix of $V_n$ converges to the block-diagonal matrix \begin{align*} \Sigma=\operatorname{diag}\!\Big(\tfrac{f(\omega_1)}{2},\tfrac{f(\omega_1)}{2},\ \dots,\ \tfrac{f(\omega_p)}{2},\tfrac{f(\omega_p)}{2}\Big)\in\mathbb R^{2p\times 2p}. \end{align*}[/step]

[guided]The whole step is one computation — the second moment $R_n(\alpha,\beta)=\mathbb E[d_n(\alpha)\overline{d_n(\beta)}]$ — evaluated in two regimes, "diagonal" (equal frequencies) and "off-diagonal" (separated frequencies). The substitution $s=t+h$ is the standard device for a stationary covariance: it replaces the double sum over $(s,t)$ by a sum over the lag $h=s-t$ times an inner sum over the base index $t$, because $\gamma(s-t)$ depends only on the lag. *Diagonal case.* When $\alpha=\beta$ the inner exponential $e^{-i(\alpha-\beta)t}$ is identically $1$, so the inner sum just counts the admissible $t$, namely $|I_h|=n-|h|$. This produces the Cesàro-type weights $1-|h|/n$, which are the only difference between $R_n(\omega_n,\omega_n)$ and the full spectral sum $f(\omega)$. To pass to the limit we want to send $n\to\infty$ inside the sum over $h$. The justification is the [Dominated Convergence Theorem](/theorems/4): we view the sum as an integral against counting measure on $\mathbb Z$, which is $\sigma$-finite; the dominating function is $|\gamma(h)|$, summable by Hypothesis (1); and the integrand converges pointwise because $1-|h|/n\to1$ for each fixed $h$ and $e^{-i\omega_n h}\to e^{-i\omega h}$ by continuity. The limit is exactly $f(\omega)$. This is where Hypothesis (1) is consumed. *Off-diagonal case.* When the two frequencies stay apart, the inner sum is a genuine geometric series in $e^{-i(\alpha_n-\beta_n)}$. A geometric series of unimodular terms is bounded by $2/|1-\text{ratio}|$ regardless of how many terms it has — the partial sums never grow, they orbit on a circle. The ratio stays bounded away from $1$ precisely because $\rho\not\equiv 0\pmod{2\pi}$, which is guaranteed by the geometry of $(0,\pi)$. With the inner sum $O(1)$ uniformly in $h$, the prefactor $1/(2\pi n)$ kills everything: $R_n=O(1/n)\to 0$. Why do all the needed frequency combinations avoid $0\pmod{2\pi}$? Three checks, all using $\omega_r\in(0,\pi)$. First, $2\omega_r\in(0,2\pi)$, so it is never a multiple of $2\pi$ — this kills the *complementary* covariance $\mathbb E[d_n^2]$ and is exactly what forces the real and imaginary parts to have equal variance and zero correlation, i.e. produces a *proper* (circularly symmetric) complex Gaussian. Second, for distinct $r,s$ the difference $\omega_r-\omega_s\in(-\pi,\pi)$ is nonzero. Third, the sum $\omega_r+\omega_s\in(0,2\pi)$ is again never a multiple of $2\pi$. These last two kill the cross-frequency covariances and are the source of asymptotic *independence* across frequencies. Had we allowed $\omega=0$ or $\omega=\pi$, the complementary covariance would not vanish and the limit would be a real (degenerate) Gaussian — which is exactly why those endpoints are excluded from the hypothesis. Finally, the bookkeeping converting the complex identities $\mathbb E[d_n^2]\to 0$ and the cross-moments into statements about $A_n,B_n$ is just expanding $d=A-iB$ and matching real and imaginary parts. The upshot is that, in the limit, every coordinate of $V_n$ becomes orthogonal to every other, and within a frequency the two parts share the variance $f(\omega_r)$ equally: the limiting covariance is the diagonal matrix $\Sigma$.[/guided]

[step:Show that all joint cumulants of order at least three vanish asymptotically]Fix $k\ge 3$ and frequencies $\lambda_1,\dots,\lambda_k$, each of which we allow to be of the form $\pm\omega_{r,n}$. By multilinearity of cumulants and the definition of $d_n$, \begin{align*} \operatorname{cum}\!\big(d_n(\lambda_1),\dots,d_n(\lambda_k)\big)=\frac{1}{(2\pi n)^{k/2}}\sum_{t_1=1}^{n}\cdots\sum_{t_k=1}^{n}\Big(\prod_{j=1}^{k}e^{-i\lambda_j t_j}\Big)\operatorname{cum}\!\big(X_{t_1},\dots,X_{t_k}\big). \end{align*} By translation invariance of the cumulants, $\operatorname{cum}(X_{t_1},\dots,X_{t_k})=c_k(t_1-t_k,\dots,t_{k-1}-t_k)$. Change variables to $u_j=t_j-t_k$ for $j=1,\dots,k-1$ and $v=t_k$; then $t_j=u_j+v$ and \begin{align*} \sum_{j=1}^{k}\lambda_j t_j=\sum_{j=1}^{k-1}\lambda_j u_j+v\sum_{j=1}^{k}\lambda_j . \end{align*} Writing $\Lambda:=\sum_{j=1}^k\lambda_j$ and letting $J_u\subseteq\{1,\dots,n\}$ be the set of $v$ for which all $t_j=u_j+v$ (and $t_k=v$) lie in $\{1,\dots,n\}$, \begin{align*} \operatorname{cum}\!\big(d_n(\lambda_1),\dots,d_n(\lambda_k)\big)=\frac{1}{(2\pi n)^{k/2}}\sum_{u\in\mathbb Z^{k-1}}c_k(u)\,e^{-i\sum_{j<k}\lambda_j u_j}\sum_{v\in J_u}e^{-i\Lambda v}. \end{align*} We use only the crude bound $\big|\sum_{v\in J_u}e^{-i\Lambda v}\big|\le|J_u|\le n$ together with $|e^{-i\sum_{j<k}\lambda_j u_j}|=1$, giving \begin{align*} \Big|\operatorname{cum}\!\big(d_n(\lambda_1),\dots,d_n(\lambda_k)\big)\Big|\le\frac{n}{(2\pi n)^{k/2}}\sum_{u\in\mathbb Z^{k-1}}|c_k(u)|=(2\pi)^{-k/2}\,n^{\,1-k/2}\sum_{u\in\mathbb Z^{k-1}}|c_k(u)|. \end{align*} By Hypothesis (2) the sum is finite, and for $k\ge 3$ the exponent $1-k/2\le-\tfrac12<0$, so the right-hand side tends to $0$. Hence \begin{align*} \operatorname{cum}\!\big(d_n(\lambda_1),\dots,d_n(\lambda_k)\big)\ \longrightarrow\ 0\qquad(k\ge 3), \end{align*} uniformly over the finite set of frequency choices $\lambda_j\in\{\pm\omega_{1,n},\dots,\pm\omega_{p,n}\}$. Now any joint cumulant of order $k\ge 3$ of the coordinates of $V_n$ is, by multilinearity and the representations $A_{r,n}=\tfrac12(d_n(\omega_{r,n})+d_n(-\omega_{r,n}))$, $B_{r,n}=\tfrac{i}{2}(d_n(\omega_{r,n})-d_n(-\omega_{r,n}))$ from Step 1, a finite linear combination — with fixed coefficients of modulus $2^{-k}$ — of cumulants of the form $\operatorname{cum}(d_n(\lambda_1),\dots,d_n(\lambda_k))$ with $\lambda_j\in\{\pm\omega_{r,n}\}$. Each such term tends to $0$, so every joint cumulant of order $\ge 3$ of $V_n$ tends to $0$.[/step]

[guided]This is the engine of the proof, and it is short once the cumulant formula is written down. Multilinearity lets us pull the linear coefficients $e^{-i\lambda_j t_j}/\sqrt{2\pi n}$ out of the cumulant, leaving the joint cumulant of the raw variables $X_{t_1},\dots,X_{t_k}$. Stationarity (translation invariance of the cumulants) then says this depends only on the $k-1$ lags $u_j=t_j-t_k$, not on the absolute position $v=t_k$. The change of variables is the same idea as in Step 2, now in $k$ dimensions: separate the $k-1$ "shape" variables $u$ from the single "translation" variable $v$. The exponential factorises as $e^{-i\sum_{j<k}\lambda_j u_j}\cdot e^{-iv\Lambda}$, with $\Lambda=\sum_j\lambda_j$. The first factor pairs with $c_k(u)$ and is summable; the second is summed over $v$. Here is the key point about the counting of powers of $n$. There are $k$ summation variables $t_1,\dots,t_k$, but after the change of variables only the *translation* variable $v$ ranges freely over $O(n)$ values for each fixed shape $u$, while the shape variables $u$ range over a set on which $|c_k(u)|$ is summable — a set of effectively *bounded* total mass, independent of $n$. So the sum contributes a factor $\le n$ (from $v$) times a finite constant (from $u$), against a normalisation $(2\pi n)^{-k/2}$. The net size is $n^{1-k/2}$. Why is this decisive only for $k\ge 3$? For $k=2$ we would get $n^{0}=O(1)$ — consistent with Step 2, where second cumulants converge to nonzero limits. For $k\ge 3$ the half-power $n^{1-k/2}$ is a negative power of $n$, so the cumulant is crushed to zero. This dichotomy is the analytic heart of the Gaussian limit: the only cumulants that survive the $(2\pi n)^{-k/2}$ normalisation are the second-order ones, and a law with vanishing cumulants of order $\ge 3$ is exactly a Gaussian. Notice we did not even need the sum over $v$ to be small — the crude bound $|J_u|\le n$ suffices — so, unlike in Step 2, we required no non-degeneracy of the frequencies here. Hypothesis (2), summability of the $k$-th order cumulants, is precisely what guarantees $\sum_u|c_k(u)|<\infty$ and is consumed once for each $k$. The final paragraph just propagates the estimate from ordinates to the coordinates of $V_n$: since each $A_{r,n},B_{r,n}$ is a fixed two-term combination of ordinates at $\pm\omega_{r,n}$, multilinearity expands any order-$k$ cumulant of coordinates into at most $2^k$ ordinate-cumulants, each vanishing in the limit.[/guided]

[step:Conclude via the method of moments and the continuous mapping theorem]By Hypotheses (1)–(2) the process has finite moments of all orders, so $V_n$ has finite moments of every order. From Steps 2 and 3, the joint cumulants of $V_n$ converge: \begin{align*} \text{order }1:\ \to 0,\qquad \text{order }2:\ \to\Sigma,\qquad \text{order }\ge 3:\ \to 0, \end{align*} with $\Sigma$ the diagonal matrix of Step 2. These are exactly the cumulants of the centred Gaussian law $\mathcal N(0,\Sigma)$ on $\mathbb R^{2p}$ (a Gaussian has vanishing cumulants of all orders $\ge 3$, first cumulant equal to its mean, and second cumulant equal to its covariance). Because moments are fixed finite polynomials in the cumulants of equal or lower order (the inverse Leonov–Shiryaev relations), convergence of all joint cumulants of $V_n$ implies convergence of all joint moments of $V_n$ to the corresponding moments of $\mathcal N(0,\Sigma)$. The multivariate normal distribution is determined by its moments: for a vector $W\sim\mathcal N(0,\Sigma)$ and any $\theta\in\mathbb R^{2p}$, the scalar $\theta\cdot W$ is $\mathcal N(0,\theta^\top\Sigma\theta)$, whose even moments $(2m-1)!!\,(\theta^\top\Sigma\theta)^m$ satisfy Carleman's condition $\sum_m \mathbb E[(\theta\cdot W)^{2m}]^{-1/(2m)}=\infty$, so each one-dimensional projection — and hence, by the Cramér–Wold principle, the joint law — is moment-determinate. Therefore, by the method of moments (the multivariate Fréchet–Shohat theorem: if the moments of $V_n$ converge to those of a moment-determinate law, then $V_n$ converges in distribution to that law), \begin{align*} V_n=\big(A_{1,n},B_{1,n},\dots,A_{p,n},B_{p,n}\big)\ \xrightarrow{\ d\ }\ \mathcal N(0,\Sigma). \end{align*} *(Citing results not yet in the wiki: the Fréchet–Shohat method-of-moments theorem and Carleman's moment-determinacy criterion.)* Because $\Sigma$ is diagonal, the limiting Gaussian vector $(A_1,B_1,\dots,A_p,B_p)\sim\mathcal N(0,\Sigma)$ has uncorrelated, jointly Gaussian, hence **independent** coordinates; in particular the pairs $(A_r,B_r)$ are independent across $r$, and within each pair $A_r$ and $B_r$ are independent with $A_r,B_r\sim\mathcal N\big(0,\tfrac12 f(\omega_r)\big)$. Finally apply the continuous map \begin{align*} \Psi:\mathbb R^{2p}&\to\mathbb C^{p}, & \Psi(a_1,b_1,\dots,a_p,b_p)&=(a_1-ib_1,\dots,a_p-ib_p), \end{align*} which is continuous, with $\Psi(V_n)=(d_n(\omega_{1,n}),\dots,d_n(\omega_{p,n}))$. By the [Continuous Mapping Theorem](/theorems/1847), \begin{align*} \big(d_n(\omega_{1,n}),\dots,d_n(\omega_{p,n})\big)\ \xrightarrow{\ d\ }\ \big(\zeta_1,\dots,\zeta_p\big),\qquad \zeta_r=A_r-iB_r, \end{align*} where the $\zeta_r$ are independent (images of independent blocks under the coordinatewise map $\Psi$) and each $\zeta_r=A_r-iB_r$ has independent real part $A_r\sim\mathcal N(0,\tfrac12 f(\omega_r))$ and imaginary part $-B_r\sim\mathcal N(0,\tfrac12 f(\omega_r))$, i.e. $\zeta_r\sim\mathcal N_{\mathbb C}(0,f(\omega_r))$. This is precisely the claimed limit, and taking $p=1$ recovers the marginal statement $d_n(\omega_n)\xrightarrow{d}\mathcal N_{\mathbb C}(0,f(\omega))$. $\qquad\blacksquare$[/step]

[guided]We now harvest the cumulant computations. The logical chain is: cumulant convergence $\Rightarrow$ moment convergence $\Rightarrow$ (because the limit is moment-determinate) convergence in distribution of the real vector $V_n$ $\Rightarrow$ (by a continuous map) convergence of the complex ordinates. *From cumulants to moments.* Cumulants and moments are linked by fixed combinatorial identities: each joint moment of order $k$ is a polynomial in the joint cumulants of order $\le k$, and vice versa. These polynomials do not depend on $n$. Hence "all cumulants of $V_n$ converge" is equivalent to "all moments of $V_n$ converge," and the limiting moments are computed from the limiting cumulants (order $1\to 0$, order $2\to\Sigma$, order $\ge 3\to 0$) — which are exactly the cumulants, hence the moments, of $\mathcal N(0,\Sigma)$. *Why moment convergence yields convergence in distribution.* In general, matching all moments is not enough — one needs the limiting law to be *determined* by its moments (the [moment problem](/theorems/1216) must be determinate). For the Gaussian this holds. The cleanest verification is via projections: by the Cramér–Wold principle it suffices that each one-dimensional projection $\theta\cdot W$, $W\sim\mathcal N(0,\Sigma)$, be moment-determinate. But $\theta\cdot W\sim\mathcal N(0,\sigma^2)$ with $\sigma^2=\theta^\top\Sigma\theta$, whose $2m$-th moment is $(2m-1)!!\,\sigma^{2m}$; since $((2m-1)!!)^{1/(2m)}$ grows like $\sqrt{2m/e}$, the series $\sum_m \big((2m-1)!!\,\sigma^{2m}\big)^{-1/(2m)}$ behaves like $\sum_m c/\sqrt{m}=\infty$, so Carleman's condition holds and the projection is moment-determinate. The method-of-moments theorem (Fréchet–Shohat) then upgrades moment convergence to convergence in distribution: $V_n\xrightarrow{d}\mathcal N(0,\Sigma)$. We flag that Fréchet–Shohat and Carleman's criterion are not yet in the wiki; everything else is self-contained. *Reading off independence.* The diagonality of $\Sigma$ is doing double duty. For a *jointly Gaussian* vector, zero covariance is equivalent to independence — a property special to the Gaussian. So the off-diagonal zeros from the cross-frequency analysis (Step 2) become independence of the ordinates at different frequencies, and the within-block diagonality $\tfrac12 f(\omega_r)I_2$ becomes independence of the real and imaginary parts of a single ordinate, each carrying half the spectral mass. *Back to the complex world.* We have proved the real statement; the map $\Psi(a,b)=a-ib$, applied coordinatewise, is continuous and sends $V_n$ to the tuple of ordinates. The [Continuous Mapping Theorem](/theorems/1847) preserves convergence in distribution under continuous maps, so the complex ordinates converge to $\Psi$ of the Gaussian limit. Since $\Psi$ acts blockwise, it preserves independence across $r$, and within each block it produces a circularly symmetric complex Gaussian $\zeta_r=A_r-iB_r$ with $\mathbb E|\zeta_r|^2=\tfrac12 f(\omega_r)+\tfrac12 f(\omega_r)=f(\omega_r)$ and $\mathbb E[\zeta_r^2]=\mathbb E[(A_r)^2-(B_r)^2]-2i\mathbb E[A_rB_r]=0$ — exactly the law $\mathcal N_{\mathbb C}(0,f(\omega_r))$. Setting $p=1$ gives the single-ordinate convergence, completing the proof.[/guided]