Why is every Gaussian process a linear process? In Section 4.2.4 of [1], the authors write

In this section we consider a causal linear process
$$
X_t = \sum_{j = 0}^\infty a_j \varepsilon_{t - j}, \quad t \in \mathbb{N},
$$
where, without loss of generality, $\sum_{j = 0}^\infty a_j^2 = 1$ and $\varepsilon_t$ $(t \in \mathbb{Z})$ are i.i.d. zero mean random variables with Var$(\varepsilon_1) = \sigma^2_\varepsilon < \infty$. Thus, Var$(X_1) = \sigma^2_X = \sigma^2_\varepsilon$. Note that
Gaussian processes are included in this definition, but the class is much more general.

I read the last sentence as "all Gaussian processes are causal linear processes". And [2] seems to support that. It argues that "...the Mallows MA closure is exhausted by three types of processes. The first type is the set of stationary Gaussian processes with mean zero...".
However, a counter-example can be constructed by a Gaussian process which is defined as $X(t) = Y$ for all $t \in \mathbb{N}$ where $Y$ is a zero-mean Gaussian random variable. This is an (admittedly degenerate) Gaussian process and I fail to see how it can be expressed as a linear process.
So, is my counter-example wrong? Or do we need additional assumptions for a Gaussian process to be a (causal) linear process? [1] often works with regularly varying auto-covariance functions resp. spectral densities. Has this maybe been assumed implicitly? If so, how do these assumptions help to prove the claim?
References
[1] J. Beran, Y. Feng, S. Ghosh, and R. Kulik. Long-memory processes. Springer, 2016.
[2] https://www.pnas.org/doi/pdf/10.1073/pnas.93.22.12128
 A: $\newcommand\ep\varepsilon\newcommand\si\sigma\newcommand\N{\mathbb N}$Your counterexample is correct. Indeed, if
$$X_t=\sum_{j=0}^\infty a_j \ep_{t-j} \tag{1}\label{1}$$
for $t\in\N$, $\sum_{j=0}^\infty a_j^2=1$, and the $\ep_t$'s are iid zero-mean random variables with variance  $\si^2\in(0,\infty)$, then for $u\in\N$
$$EX_tX_{t+u}=\si^2\sum_{j=0}^\infty a_ja_{j+u}\to0$$
as $u\to\infty$, by the Cauchy--Schwarz inequality.
This contradicts the condition that $X_t=Y$ for all $t$, because then we would have $EX_tX_{t+u}=EX_1^2=\si^2>0$ for all $u\in\N$.

Generally (see e.g. Theorem 2, p. 263), a  Gaussian process $(X_t)$ is of the form \eqref{1} if and only if $(X_t)$ is a zero-mean stationary process with $EX_t^2<\infty$ which has an absolutely continuous spectral measure and its spectral
density $f$ satisfies the condition
$$\int_{-\pi}^\pi\ln f(t)\,dt>-\infty$$
(so that $f$ does not get too close to $0$ on any set of significant enough Lebesgue measure).
(In your counterexample, the spectral measure is degenerate, a multiple of the Dirac measure supported on the singleton set $\{0\}$, and hence not all absolutely continuous.)
