[proofplan]
We first use positivity and traciality to show that the covariance entries are real and symmetric. Then we test the covariance matrix against an arbitrary real vector $a\in\mathbb R^d$. From $a$ we form the self-adjoint element $x=\sum_i a_i s_i$ of the noncommutative probability space. Positivity of the state gives $\varphi(x^*x)\geq 0$, and self-adjointness gives $x^*x=x^2$. Expanding $x^2$ by bilinearity identifies this nonnegative number with the quadratic form $a^\top C a$, which is exactly positive semidefiniteness.
[/proofplan]
[step:Show that the covariance entries are real and symmetric]
For $1\leq i,j\leq d$, traciality gives
\begin{align*}
C_{ij}=\varphi(s_i s_j)=\varphi(s_j s_i)=C_{ji}.
\end{align*}
Since $s_i$ and $s_j$ are self-adjoint, $(s_i s_j)^*=s_j s_i$. A positive linear functional on a $*$-algebra satisfies $\varphi(y^*)=\overline{\varphi(y)}$ for every $y\in\mathcal A$. Applying this to $y=s_i s_j$ gives
\begin{align*}
\overline{C_{ij}}=\overline{\varphi(s_i s_j)}=\varphi((s_i s_j)^*)=\varphi(s_j s_i)=C_{ji}.
\end{align*}
Together with $C_{ij}=C_{ji}$, this implies $\overline{C_{ij}}=C_{ij}$. Hence every entry $C_{ij}$ is real, and the matrix $C$ is symmetric.
[/step]
[step:Form the self-adjoint linear combination associated to a real vector]
Fix an arbitrary vector $a=(a_1,\dots,a_d)\in\mathbb R^d$. Define the element $x\in\mathcal A$ by
\begin{align*}
x=\sum_{i=1}^{d} a_i s_i.
\end{align*}
Since each $a_i\in\mathbb R$ and each $s_i$ is self-adjoint, the involution on $\mathcal A$ gives
\begin{align*}
x^*=\left(\sum_{i=1}^{d} a_i s_i\right)^*=\sum_{i=1}^{d} a_i s_i^*= \sum_{i=1}^{d} a_i s_i=x.
\end{align*}
Thus $x$ is self-adjoint.
[guided]
Fix an arbitrary vector $a=(a_1,\dots,a_d)\in\mathbb R^d$. To prove that $C$ is positive semidefinite, we must show that the quadratic form determined by $C$ is nonnegative on this arbitrary vector. The natural element of $\mathcal A$ corresponding to $a$ is the real linear combination of the semicircular variables. Define $x\in\mathcal A$ by
\begin{align*}
x=\sum_{i=1}^{d} a_i s_i.
\end{align*}
We need $x$ to be self-adjoint because positivity of the state applies to elements of the form $y^*y$. The coefficients are real, so the involution does not conjugate them into anything new. Since $s_i^*=s_i$ for every $i$, we compute
\begin{align*}
x^*=\left(\sum_{i=1}^{d} a_i s_i\right)^*=\sum_{i=1}^{d} a_i s_i^*=\sum_{i=1}^{d} a_i s_i=x.
\end{align*}
Therefore $x$ is self-adjoint.
[/guided]
[/step]
[step:Apply positivity of the state to obtain a nonnegative second moment]
Because $\varphi$ is positive, for every $y\in\mathcal A$ one has $\varphi(y^*y)\geq 0$. Applying this with $y=x$ gives
\begin{align*}
\varphi(x^*x)\geq 0.
\end{align*}
Since $x^*=x$, this becomes
\begin{align*}
\varphi(x^2)\geq 0.
\end{align*}
[/step]
[step:Expand the second moment as the covariance quadratic form]
Using multiplication in $\mathcal A$ and linearity of $\varphi$, we expand
\begin{align*}
\varphi(x^2)=\varphi\left(\left(\sum_{i=1}^{d} a_i s_i\right)\left(\sum_{j=1}^{d} a_j s_j\right)\right).
\end{align*}
Distributivity in $\mathcal A$ gives
\begin{align*}
\left(\sum_{i=1}^{d} a_i s_i\right)\left(\sum_{j=1}^{d} a_j s_j\right)=\sum_{i=1}^{d}\sum_{j=1}^{d} a_i a_j s_i s_j.
\end{align*}
Applying linearity of $\varphi$ and the defining relation $C_{ij}=\varphi(s_i s_j)$ yields
\begin{align*}
\varphi(x^2)=\sum_{i=1}^{d}\sum_{j=1}^{d} a_i a_j \varphi(s_i s_j)=\sum_{i=1}^{d}\sum_{j=1}^{d} a_i C_{ij} a_j.
\end{align*}
[/step]
[step:Conclude positive semidefiniteness of the covariance matrix]
Combining the nonnegativity from positivity of $\varphi$ with the expansion of $\varphi(x^2)$ gives
\begin{align*}
\sum_{i=1}^{d}\sum_{j=1}^{d} a_i C_{ij} a_j=\varphi(x^2)\geq 0.
\end{align*}
The vector $a\in\mathbb R^d$ was arbitrary. Hence the quadratic form associated to $C$ is nonnegative on every vector in $\mathbb R^d$, so $C$ is positive semidefinite.
[/step]
