[proofplan]
We first realize the positive semidefinite matrix $C$ as the Gram matrix of vectors $v_1,\dots,v_d$ in a real [Hilbert space](/page/Hilbert%20Space) $H$. We then pass to the full Fock space over the complexification $H_{\mathbb{C}}$ and define the field operators $s(v)=\ell(v)+\ell(v)^*$. The vacuum state gives the desired covariance by the creation-annihilation relations, and the [free Wick formula](/theorems/7142) for these Fock-space field operators shows that the joint distribution is semicircular. Finally, the same Wick formula implies that the vacuum state is tracial on the unital $*$-algebra generated by the field operators.
[/proofplan]
[step:Realize the covariance matrix as a Gram matrix]
Define a symmetric [bilinear form](/page/Bilinear%20Form) $b: \mathbb{R}^d \times \mathbb{R}^d \to \mathbb{R}$ by
\begin{align*}
b(x,y)=x^\top C y.
\end{align*}
Since $C$ is positive semidefinite, $b(x,x)\geq 0$ for every $x \in \mathbb{R}^d$. Let
\begin{align*}
N=\{x \in \mathbb{R}^d : b(x,x)=0\}.
\end{align*}
By the [Cauchy-Schwarz inequality](/theorems/432) for positive semidefinite bilinear forms, $b(x,y)=0$ for every $x \in N$ and every $y \in \mathbb{R}^d$, so $N$ is the radical of $b$.
Let $H=\mathbb{R}^d/N$, and define an [inner product](/page/Inner%20Product) $(\cdot,\cdot)_H: H \times H \to \mathbb{R}$ by
\begin{align*}
(x+N,y+N)_H=b(x,y).
\end{align*}
This is well-defined because $N$ is the radical of $b$, and it is positive definite by construction. Hence $H$ is a finite-dimensional real Hilbert space. For $1 \leq i \leq d$, let $e_i \in \mathbb{R}^d$ denote the $i$-th standard basis vector and set $v_i=e_i+N \in H$. Then, for all $1 \leq i,j \leq d$,
\begin{align*}
(v_i,v_j)_H=b(e_i,e_j)=C_{ij}.
\end{align*}
[guided]
The first task is to build vectors whose inner products reproduce the entries of $C$. Define
\begin{align*}
b: \mathbb{R}^d \times \mathbb{R}^d \to \mathbb{R}, \qquad b(x,y)=x^\top C y.
\end{align*}
The symmetry of $C$ makes $b$ symmetric, and positive semidefiniteness of $C$ gives $b(x,x)\geq 0$ for all $x \in \mathbb{R}^d$.
If $C$ is not positive definite, some nonzero vectors may have zero $b$-length. We remove exactly those directions. Define
\begin{align*}
N=\{x \in \mathbb{R}^d : b(x,x)=0\}.
\end{align*}
For a positive semidefinite bilinear form, the Cauchy-Schwarz inequality says
\begin{align*}
|b(x,y)|^2 \leq b(x,x)b(y,y)
\end{align*}
for all $x,y \in \mathbb{R}^d$. Thus if $x \in N$, then $b(x,x)=0$, so $b(x,y)=0$ for every $y \in \mathbb{R}^d$. This means precisely that $N$ is the radical of $b$.
Now form the quotient [vector space](/page/Vector%20Space) $H=\mathbb{R}^d/N$. Define
\begin{align*}
(x+N,y+N)_H=b(x,y).
\end{align*}
This is well-defined: changing $x$ or $y$ by an element of $N$ does not change $b(x,y)$ because every element of $N$ pairs to zero with every vector. The form is positive definite on the quotient, since $(x+N,x+N)_H=0$ means $b(x,x)=0$, hence $x \in N$, hence $x+N=0$ in $H$. Because $H$ is finite-dimensional, this inner product makes $H$ a real Hilbert space.
Finally, let $e_i \in \mathbb{R}^d$ be the $i$-th standard basis vector and define $v_i=e_i+N \in H$. Then
\begin{align*}
(v_i,v_j)_H=b(e_i,e_j)=e_i^\top C e_j=C_{ij}.
\end{align*}
Thus $C$ has been realized as the Gram matrix of $v_1,\dots,v_d$.
[/guided]
[/step]
[step:Construct the Fock space field operators]
Let $H_{\mathbb{C}}=H \otimes_{\mathbb{R}} \mathbb{C}$ be the complex Hilbert space obtained by complexifying $H$, with inner product extending $(\cdot,\cdot)_H$ and linear in the first variable. Define the full Fock space
\begin{align*}
\mathcal{F}(H_{\mathbb{C}})=\mathbb{C}\Omega \oplus \bigoplus_{n=1}^{\infty} H_{\mathbb{C}}^{\otimes n},
\end{align*}
where $\Omega$ is a distinguished unit vector called the vacuum vector.
For each $h \in H_{\mathbb{C}}$, define the left creation operator $\ell(h):\mathcal{F}(H_{\mathbb{C}})\to\mathcal{F}(H_{\mathbb{C}})$ by
\begin{align*}
\ell(h)\Omega=h
\end{align*}
and, for $n \geq 1$ and $\xi_1,\dots,\xi_n \in H_{\mathbb{C}}$,
\begin{align*}
\ell(h)(\xi_1 \otimes \cdots \otimes \xi_n)=h \otimes \xi_1 \otimes \cdots \otimes \xi_n.
\end{align*}
This operator is bounded with $\|\ell(h)\|_{\mathrm{op}}=|h|_{H_{\mathbb{C}}}$. Its adjoint $\ell(h)^*:\mathcal{F}(H_{\mathbb{C}})\to\mathcal{F}(H_{\mathbb{C}})$ satisfies
\begin{align*}
\ell(h)^*\Omega=0
\end{align*}
and
\begin{align*}
\ell(h)^*(\xi_1 \otimes \cdots \otimes \xi_n)=(\xi_1,h)_{H_{\mathbb{C}}}\,\xi_2 \otimes \cdots \otimes \xi_n
\end{align*}
for $n \geq 1$.
For $h \in H$, regarded as a real vector inside $H_{\mathbb{C}}$, define the field operator
\begin{align*}
s(h):\mathcal{F}(H_{\mathbb{C}})\to\mathcal{F}(H_{\mathbb{C}}), \qquad s(h)=\ell(h)+\ell(h)^*.
\end{align*}
Since $\ell(h)^*$ is the Hilbert-space adjoint of $\ell(h)$, each $s(h)$ is self-adjoint. Let
\begin{align*}
\mathcal{L}(\mathcal{F}(H_{\mathbb{C}})) &= \{T: \mathcal{F}(H_{\mathbb{C}})\to\mathcal{F}(H_{\mathbb{C}}) : T \text{ is bounded and linear}\}
\end{align*}
be the algebra of bounded linear operators on $\mathcal{F}(H_{\mathbb{C}})$. Let $\mathcal{A}$ be the unital $*$-algebra generated by $I_{\mathcal{F}(H_{\mathbb{C}})}$ and $s(v_1),\dots,s(v_d)$ inside $\mathcal{L}(\mathcal{F}(H_{\mathbb{C}}))$. Define the vacuum state $\varphi:\mathcal{A}\to\mathbb{C}$ by
\begin{align*}
\varphi(T)=(T\Omega,\Omega)_{\mathcal{F}(H_{\mathbb{C}})}.
\end{align*}
Then $\varphi(I_{\mathcal{F}(H_{\mathbb{C}})})=1$ and $\varphi(T^*T)\geq 0$ for every $T \in \mathcal{A}$.
[/step]
[step:Compute the covariance from the creation-annihilation relations]
For $1 \leq i,j \leq d$, set $s_i=s(v_i)$. Since $s(v_j)\Omega=\ell(v_j)\Omega+\ell(v_j)^*\Omega=v_j$, we have
\begin{align*}
\varphi(s_i s_j)=(s_i s_j\Omega,\Omega)_{\mathcal{F}(H_{\mathbb{C}})}=(s(v_i)v_j,\Omega)_{\mathcal{F}(H_{\mathbb{C}})}.
\end{align*}
Using $s(v_i)=\ell(v_i)+\ell(v_i)^*$, the creation term $\ell(v_i)v_j=v_i \otimes v_j$ lies in $H_{\mathbb{C}}^{\otimes 2}$ and is orthogonal to $\Omega$, while the annihilation term satisfies
\begin{align*}
\ell(v_i)^*v_j=(v_j,v_i)_{H_{\mathbb{C}}}\Omega.
\end{align*}
Because $v_i,v_j \in H$ and the complex inner product extends the real one, $(v_j,v_i)_{H_{\mathbb{C}}}=(v_j,v_i)_H=(v_i,v_j)_H$. Hence
\begin{align*}
\varphi(s_i s_j)=(v_i,v_j)_H=C_{ij}.
\end{align*}
[/step]
[step:Verify the semicircular Wick formula]
Let $m \in \mathbb{N}$ and let $h_1,\dots,h_m \in H$. Consider the vacuum moment
\begin{align*}
\varphi(s(h_1)\cdots s(h_m))=(s(h_1)\cdots s(h_m)\Omega,\Omega)_{\mathcal{F}(H_{\mathbb{C}})}.
\end{align*}
Expanding each $s(h_k)=\ell(h_k)+\ell(h_k)^*$ gives a sum indexed by words in creation and annihilation operators. A summand contributes to the vacuum inner product exactly when the operators reduce the vacuum back to the vacuum without ever trying to annihilate the vacuum. Such contributing words are in bijection with noncrossing pairings $\pi$ of $\{1,\dots,m\}$. If $m$ is odd, there are no such pairings, so the moment is $0$. If $m$ is even, the contribution of a pair $\{p,q\}$ with $p<q$ is $(h_p,h_q)_H$, and therefore
\begin{align*}
\varphi(s(h_1)\cdots s(h_m))=\sum_{\pi \in NC_2(m)} \prod_{\{p,q\}\in \pi,\,p<q} (h_p,h_q)_H,
\end{align*}
where $NC_2(m)$ denotes the set of noncrossing pairings of $\{1,\dots,m\}$.
This is precisely the free Wick formula for a semicircular family: all free cumulants vanish except the second-order cumulants, and the second-order cumulant of $s(h)$ and $s(k)$ equals $(h,k)_H$ for $h,k \in H$. Applying this with $h_k \in \{v_1,\dots,v_d\}$ shows that $(s_1,\dots,s_d)$ is a semicircular family with covariance matrix $C$.
[guided]
We now verify that the constructed variables have the joint semicircular distribution. Fix $m \in \mathbb{N}$ and vectors $h_1,\dots,h_m \in H$. We need to compute
\begin{align*}
\varphi(s(h_1)\cdots s(h_m))=(s(h_1)\cdots s(h_m)\Omega,\Omega)_{\mathcal{F}(H_{\mathbb{C}})}.
\end{align*}
Each factor splits as $s(h_k)=\ell(h_k)+\ell(h_k)^*$, so expanding the product gives a sum of vacuum inner products of words made from creation operators and annihilation operators.
A word can contribute only if it starts from $\Omega$, never applies an annihilation operator to $\Omega$ at an intermediate stage, and ends back at $\Omega$. Creation operators increase tensor length by one, and annihilation operators decrease tensor length by one by pairing with the leftmost tensor factor. Therefore a contributing word is exactly a balanced parenthesis pattern: creations open pairs and later annihilations close them. Because annihilation always removes the leftmost available tensor factor, the pairings produced in this way are noncrossing. Conversely, every noncrossing pairing determines exactly one such admissible creation-annihilation word.
For one matched pair $\{p,q\}$ with $p<q$, the vector created from $h_q$ is later annihilated against $h_p$ after accounting for the order in the operator product. The annihilation formula gives the scalar contribution
\begin{align*}
(h_p,h_q)_H.
\end{align*}
Multiplying over all matched pairs and summing over all admissible words gives
\begin{align*}
\varphi(s(h_1)\cdots s(h_m))=\sum_{\pi \in NC_2(m)} \prod_{\{p,q\}\in \pi,\,p<q} (h_p,h_q)_H.
\end{align*}
If $m$ is odd, $NC_2(m)$ is empty, so the moment is $0$.
This moment formula is the free Wick formula. Equivalently, it says that the only nonzero free cumulants of the family $\{s(h):h\in H\}$ are the second-order cumulants, with
\begin{align*}
\kappa_2(s(h),s(k))=(h,k)_H.
\end{align*}
Thus any finite list $s(h_1),\dots,s(h_m)$ is a semicircular family with covariance given by the Hilbert-space inner products. In particular, taking $h_k$ from the list $v_1,\dots,v_d$ proves that $(s_1,\dots,s_d)$ is a semicircular family and that its covariance matrix is $C$.
[/guided]
[/step]
[step:Show that the vacuum state is tracial on the generated algebra]
It remains to check traciality on $\mathcal{A}$. Since $\mathcal{A}$ is linearly spanned by products of the generators $s(v_i)$ and the identity, it is enough to verify cyclic invariance of vacuum moments of words in the generators.
Let $m \in \mathbb{N}$ and let $i_1,\dots,i_m \in \{1,\dots,d\}$. By the Wick formula just proved,
\begin{align*}
\varphi(s_{i_1}\cdots s_{i_m})=\sum_{\pi \in NC_2(m)} \prod_{\{p,q\}\in \pi,\,p<q} C_{i_p i_q}.
\end{align*}
The cyclic rotation map on $\{1,\dots,m\}$, viewed on the cyclically ordered vertices of an $m$-gon, sends noncrossing pairings bijectively to noncrossing pairings, and the product attached to a pairing is unchanged because $C$ is symmetric. Hence
\begin{align*}
\varphi(s_{i_1}\cdots s_{i_m})=\varphi(s_{i_2}\cdots s_{i_m}s_{i_1}).
\end{align*}
By linearity, $\varphi(AB)=\varphi(BA)$ for all $A,B \in \mathcal{A}$. Therefore $(\mathcal{A},\varphi)$ is a tracial noncommutative probability space.
Together with the covariance computation and the semicircular Wick formula, this proves the existence of a semicircular family $(s_1,\dots,s_d)$ satisfying $\varphi(s_i s_j)=C_{ij}$ for all $1 \leq i,j \leq d$.
[/step]