I have been working with the d-dimensional Ornstein-Uhlenbeck process $dX_t^{\epsilon}=-QX_t^{\epsilon}dt + \epsilon dB_t$, $X_0^{\epsilon}=x \in \mathbb{R}$. I know that, by Ito's formula, you can obtain its strong solution as $$X_t^\epsilon(x)=e^{-Qt}x+\epsilon \int_0^te^{-Q(t-s)}dBs.$$
Where $\int_0^te^{-Q(t-s)}dBs$ is a Wiener integral and so $\int_0^te^{-Q(t-s)}dBs\sim \mathcal{N}(0,\int_0^te^{-2Q(t-s)}ds)$. Also, I assumed that $x, \epsilon$ and $Q$ are fixed, and $Q$ is invertible and with all of its eigenvalues with positive real component. With this, I computed $$\int_0^te^{-2Q(t-s)}ds=\frac{Q^{-1}}{2}(Id - e^{-2Qt})$$ and so, by using the properties of normal distributions $$X_t^\epsilon(x) \sim \mathcal{N}(e^{-Qt}x, \epsilon^2 \frac{Q^{-1}}{2}(Id - e^{-2Qt}))$$
My question here is, how can I determine that $\epsilon^2 \frac{Q^{-1}}{2}(Id - e^{-2Qt})$ is indeed a covariance matrix? (i.e that it is positive semi-definite and symmetric)
EDIT: I think that what makes no sense to me, is that supposedly one could obtain the invariant meassure $\mu$ by making $t\to \infty$. However, by the hypothesis on $Q's$ eigenvalues one gets that $$\epsilon^2 \frac{Q^{-1}}{2}(Id - e^{-2Qt}) \to \frac{\epsilon^2 Q^{-1}}{2}$$
And so, $\mu \sim \mathcal{N}(0, \frac{\epsilon^2 Q^{-1}}{2})$, however without any more restrictions for $Q$, $\frac{\epsilon^2 Q^{-1}}{2}$ would not be a covariance matrix.