[proofplan] We construct the Gaussian Hilbert space isometry explicitly from a countable orthonormal basis of $H$. The strategy is: (1) choose an orthonormal basis $(e_i)$ using separability, (2) build an i.i.d. standard Gaussian sequence $(X_i)$ on some probability space, (3) define $I$ on basis elements as $I(e_i) = X_i$ and extend by linearity to all of $H$ via $L^2$-convergent series, (4) verify the distribution, isometry, and closedness properties. [/proofplan] [step:Choose an orthonormal basis and build an i.i.d. Gaussian sequence] Since $H$ is a separable Hilbert space, there exists a countable orthonormal basis $(e_i)_{i=1}^\infty$ (or $(e_i)_{i=1}^d$ if $\dim H = d < \infty$; we treat the infinite-dimensional case, the finite case being identical with finite sums). Let $(\Omega, \mathcal{F}, \mathbb{P})$ be a probability space supporting a sequence $(X_i)_{i=1}^\infty$ of independent random variables, each distributed as $X_i \sim N(0,1)$. Such a probability space exists by the Kolmogorov Extension Theorem: take $\Omega = \mathbb{R}^\infty$, $\mathcal{F} = \mathcal{B}(\mathbb{R})^{\otimes \infty}$, and $\mathbb{P} = \bigotimes_{i=1}^\infty \gamma_1$, where $\gamma_1$ is the standard Gaussian measure on $(\mathbb{R}, \mathcal{B}(\mathbb{R}))$. [guided] We need a probability space carrying infinitely many independent standard Gaussians. Why does such a space exist? The Kolmogorov Extension Theorem guarantees that for any consistent family of finite-dimensional distributions, there exists a probability space and a stochastic process realising them. Here the finite-dimensional distributions are simply products of the standard Gaussian measure $\gamma_1$ on $\mathbb{R}$, which are trivially consistent under projection. Concretely, we take $\Omega = \mathbb{R}^\infty$ with the product $\sigma$-algebra $\mathcal{F} = \mathcal{B}(\mathbb{R})^{\otimes \infty}$ and the product measure $\mathbb{P} = \bigotimes_{i=1}^\infty \gamma_1$. The coordinate projections $X_i(\omega) = \omega_i$ are then i.i.d. $N(0,1)$. Since $H$ is separable, it admits a countable (possibly finite) orthonormal basis $(e_i)_{i=1}^\infty$. Every $f \in H$ has a unique expansion $f = \sum_{i=1}^\infty a_i e_i$ with $a_i = (f, e_i)_H$ and $\sum_{i=1}^\infty a_i^2 = \|f\|_H^2 < \infty$ by Parseval's identity. [/guided] [/step] [step:Define $I$ on basis elements and extend by $L^2(\mathbb{P})$-convergent series] Set $I(e_i) := X_i$ for each $i \geq 1$. For a general element $f = \sum_{i=1}^\infty a_i e_i \in H$ with $a_i = (f, e_i)_H$, define \begin{align*} I(f) := \sum_{i=1}^\infty a_i X_i. \end{align*} We verify $L^2(\mathbb{P})$-convergence of this series. Let $S_n := \sum_{i=1}^n a_i X_i$. Since the $X_i$ are independent with $\mathbb{E}[X_i] = 0$ and $\mathbb{E}[X_i^2] = 1$, for $m < n$: \begin{align*} \mathbb{E}\!\left[(S_n - S_m)^2\right] = \mathbb{E}\!\left[\left(\sum_{i=m+1}^n a_i X_i\right)^2\right] = \sum_{i=m+1}^n a_i^2, \end{align*} where the cross terms vanish because $\mathbb{E}[X_i X_j] = 0$ for $i \neq j$ (independence and zero mean). Since $\sum_{i=1}^\infty a_i^2 = \|f\|_H^2 < \infty$, the partial sums form a Cauchy sequence in $L^2(\mathbb{P})$. By completeness of $L^2(\mathbb{P})$, the series converges in $L^2(\mathbb{P})$ to a random variable $I(f)$. [guided] The key calculation is why the cross terms vanish. For $i \neq j$, independence of $X_i$ and $X_j$ gives $\mathbb{E}[X_i X_j] = \mathbb{E}[X_i]\mathbb{E}[X_j] = 0 \cdot 0 = 0$. So expanding the square of the partial sum: \begin{align*} \mathbb{E}\!\left[\left(\sum_{i=m+1}^n a_i X_i\right)^2\right] = \sum_{i=m+1}^n \sum_{j=m+1}^n a_i a_j \mathbb{E}[X_i X_j] = \sum_{i=m+1}^n a_i^2 \cdot 1 = \sum_{i=m+1}^n a_i^2. \end{align*} This is the tail of the convergent series $\sum a_i^2 = \|f\|_H^2$, so $\sum_{i=m+1}^n a_i^2 \to 0$ as $m, n \to \infty$. The partial sums $(S_n)$ are therefore Cauchy in the Hilbert space $L^2(\Omega, \mathcal{F}, \mathbb{P})$, which is complete, so $S_n \to I(f)$ in $L^2(\mathbb{P})$ for some limit $I(f)$. Note that linearity of $I$ is built into the definition: if $f = \sum a_i e_i$ and $g = \sum b_i e_i$, then $\alpha f + \beta g = \sum (\alpha a_i + \beta b_i) e_i$, so $I(\alpha f + \beta g) = \sum (\alpha a_i + \beta b_i) X_i = \alpha \sum a_i X_i + \beta \sum b_i X_i = \alpha I(f) + \beta I(g)$, where the last equality holds in $L^2(\mathbb{P})$ by linearity of the $L^2$ limit. [/guided] [/step] [step:Verify $I(f) \sim N(0, \|f\|_H^2)$ and $\mathbb{E}[I(f)I(g)] = (f,g)_H$] **Distribution.** Each partial sum $S_n = \sum_{i=1}^n a_i X_i$ is a finite linear combination of independent $N(0,1)$ variables, hence $S_n \sim N(0, \sum_{i=1}^n a_i^2)$. The characteristic function of $S_n$ is \begin{align*} \phi_{S_n}(u) = \mathbb{E}[e^{iuS_n}] = \exp\!\left(-\frac{u^2}{2}\sum_{i=1}^n a_i^2\right). \end{align*} Since $S_n \to I(f)$ in $L^2(\mathbb{P})$, convergence in $L^2$ implies convergence in probability, which implies convergence in distribution. By the continuity theorem for characteristic functions: \begin{align*} \phi_{I(f)}(u) = \lim_{n \to \infty} \phi_{S_n}(u) = \exp\!\left(-\frac{u^2}{2}\|f\|_H^2\right), \end{align*} which is the characteristic function of $N(0, \|f\|_H^2)$. Since the characteristic function determines the distribution uniquely, $I(f) \sim N(0, \|f\|_H^2)$. **Covariance.** For $f = \sum a_i e_i$ and $g = \sum b_i e_i$, the $L^2(\mathbb{P})$-convergence of $I(f) = \sum a_i X_i$ and $I(g) = \sum b_i X_i$ together with the Cauchy-Schwarz inequality gives: \begin{align*} \mathbb{E}[I(f)I(g)] = \lim_{n \to \infty} \mathbb{E}\!\left[\left(\sum_{i=1}^n a_i X_i\right)\left(\sum_{j=1}^n b_j X_j\right)\right] = \lim_{n \to \infty} \sum_{i=1}^n a_i b_i = \sum_{i=1}^\infty a_i b_i = (f,g)_H, \end{align*} where we used $\mathbb{E}[X_i X_j] = \delta_{ij}$ and Parseval's identity $(f,g)_H = \sum_{i=1}^\infty (f,e_i)_H(g,e_i)_H = \sum_{i=1}^\infty a_i b_i$. [guided] Why does convergence in $L^2(\mathbb{P})$ suffice to pass to the limit in the covariance computation? We need $\mathbb{E}[I(f)I(g)] = \lim_n \mathbb{E}[S_n^f S_n^g]$ where $S_n^f = \sum_{i=1}^n a_i X_i$ and $S_n^g = \sum_{i=1}^n b_i X_i$. Write: \begin{align*} |\mathbb{E}[I(f)I(g)] - \mathbb{E}[S_n^f S_n^g]| &= |\mathbb{E}[I(f)I(g) - S_n^f S_n^g]| \\ &= |\mathbb{E}[I(f)(I(g) - S_n^g)] + \mathbb{E}[(I(f) - S_n^f) S_n^g]| \\ &\leq \|I(f)\|_{L^2} \|I(g) - S_n^g\|_{L^2} + \|I(f) - S_n^f\|_{L^2} \|S_n^g\|_{L^2}, \end{align*} by the Cauchy-Schwarz inequality in $L^2(\mathbb{P})$. Since $\|I(f)\|_{L^2} = \|f\|_H < \infty$, $\|S_n^g\|_{L^2} \leq \|g\|_H$, and both $\|I(g) - S_n^g\|_{L^2} \to 0$ and $\|I(f) - S_n^f\|_{L^2} \to 0$ by $L^2$ convergence, the right-hand side tends to $0$. For the distribution claim, we use a standard fact: $L^2$ convergence implies convergence in probability (by Chebyshev's inequality: $\mathbb{P}(|S_n - I(f)| > \varepsilon) \leq \varepsilon^{-2}\mathbb{E}[(S_n - I(f))^2] \to 0$), which implies convergence in distribution. The characteristic function $\phi_{S_n}(u) \to \phi_{I(f)}(u)$ for each fixed $u$ by the definition of convergence in distribution. Since the limit $\exp(-u^2 \|f\|_H^2/2)$ is continuous at $u = 0$, the Levy Continuity Theorem identifies $I(f)$ as $N(0, \|f\|_H^2)$. [/guided] [/step] [step:Verify $I$ is an isometry and $S = I(H)$ is a closed Gaussian subspace] **Isometry.** For any $f \in H$: \begin{align*} \|I(f)\|_{L^2(\mathbb{P})}^2 = \mathbb{E}[I(f)^2] = (f,f)_H = \|f\|_H^2, \end{align*} using the covariance identity from the previous step with $g = f$. Hence $I: H \to L^2(\Omega, \mathcal{F}, \mathbb{P})$ is a linear isometry. **Closedness.** The image $S := I(H) \subseteq L^2(\Omega, \mathcal{F}, \mathbb{P})$ is closed. To see this, let $(I(f_n))_{n \geq 1}$ be a sequence in $S$ converging to some $Z \in L^2(\mathbb{P})$. Since $I$ is an isometry, $\|f_n - f_m\|_H = \|I(f_n) - I(f_m)\|_{L^2(\mathbb{P})} \to 0$, so $(f_n)$ is Cauchy in $H$. By completeness of $H$, there exists $f \in H$ with $f_n \to f$. By the isometry, $I(f_n) \to I(f)$ in $L^2(\mathbb{P})$. Uniqueness of limits gives $Z = I(f) \in S$. **Gaussian subspace.** Every element of $S$ is a centered Gaussian random variable. To verify this, let $Z = I(f)$ for some $f \in H$. We showed $I(f) \sim N(0, \|f\|_H^2)$, which is centered Gaussian. Moreover, any finite linear combination $\sum_{k=1}^m c_k I(f_k) = I(\sum_{k=1}^m c_k f_k)$ is also centered Gaussian (since $\sum c_k f_k \in H$), so $S$ is a Gaussian subspace of $L^2(\Omega, \mathcal{F}, \mathbb{P})$. [guided] Why is the closedness argument valid? The crucial input is that $I$ is an isometry, which means it preserves distances: $\|I(f) - I(g)\|_{L^2} = \|f - g\|_H$. This converts $L^2$-Cauchy sequences in the image back to Cauchy sequences in $H$. Since $H$ is a Hilbert space (hence complete), the preimage sequence converges, and the isometry pushes that convergence forward into $L^2$. A subtle point: why is $S$ a Gaussian subspace, not merely a set of Gaussian random variables? The definition requires that every finite linear combination of elements of $S$ is again Gaussian. This holds because $I$ is linear: $\sum c_k I(f_k) = I(\sum c_k f_k) \sim N(0, \|\sum c_k f_k\|_H^2)$. What about $L^2$ limits of elements of $S$? If $I(f_n) \to Z$ in $L^2$, we showed $Z = I(f) \in S$, so $Z$ is itself Gaussian. More generally, any element of the $L^2$-closure of a Gaussian subspace is Gaussian (since $L^2$ convergence implies convergence in distribution, and limits in distribution of centered Gaussians are centered Gaussian). Since $S$ is already closed, this is automatic. [/guided] [/step] [step:Verify the almost sure linearity identity] For $\alpha, \beta \in \mathbb{R}$ and $f, g \in H$, we must show $I(\alpha f + \beta g) = \alpha I(f) + \beta I(g)$ almost surely. Write $f = \sum a_i e_i$ and $g = \sum b_i e_i$. Then $\alpha f + \beta g = \sum (\alpha a_i + \beta b_i) e_i$, so: \begin{align*} I(\alpha f + \beta g) = \sum_{i=1}^\infty (\alpha a_i + \beta b_i) X_i = \alpha \sum_{i=1}^\infty a_i X_i + \beta \sum_{i=1}^\infty b_i X_i = \alpha I(f) + \beta I(g), \end{align*} where all three series converge in $L^2(\mathbb{P})$, and the rearrangement is justified because $L^2$-convergent series are linear. In particular, the equality $I(\alpha f + \beta g) = \alpha I(f) + \beta I(g)$ holds in $L^2(\mathbb{P})$, hence holds $\mathbb{P}$-almost surely (two $L^2$-equal random variables agree a.s. after choosing appropriate representatives). This completes the proof: we have constructed a probability space $(\Omega, \mathcal{F}, \mathbb{P})$, a Gaussian subspace $S = I(H) \subseteq L^2(\Omega, \mathcal{F}, \mathbb{P})$, and a linear isometry $I: H \to S$ satisfying $I(f) \sim N(0, \|f\|_H^2)$ and $\mathbb{E}[I(f)I(g)] = (f,g)_H$. [/step]