I am trying to compute the distribution of 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. For the unidimensional case, I've seen that the Wiener integral distribute as $\int_0^t g(s)dBs \sim \mathcal{N}(0, \int_0^t g^2(s) ds)$ and I will like to use this to determine the distribution of the d-dimensional case. It is clear that the mean vector of $X_t^\epsilon(x)$ should be $e^{-Qt}x$, however I am struggling with finding the covariance of the d-dimensional Wiener integral associated.
Any references or ideas on how to proceed are very much welcomed.
Edit: In this case, I am assuming that $Q$ is a $d\times d$ deterministic matrix, which is exponentially stable (the real components of its eigenvalues are all positive) and $\epsilon$ is a scalar.
We can solve the PDE of the density $$\frac{\partial p(y,t)}{\partial t}=\sum_{k=1}^dQ_k\frac{\partial}{\partial y_k}(y_kp(y,t))+\sum_{k=1}^d\frac{\epsilon_k^2}{2}\frac{\partial^2 p(y,t)}{\partial y_k^2}, \ \ \ y \in \mathbb{R}^d$$ with initial condition $p(y,0)=\prod_{k=1}^d\delta(y_k-x_k)$. This yields a system of uncorrelated Gaussian random variables such that $$X_{k,t}(x_k)\sim \mathcal{N}\bigg(x_ke^{-Q_k t},\frac{\epsilon_k^2}{2Q_k}(1-e^{-2Q_k t})\bigg),k=\{1,2,...,d\}$$ or in matrix form $$X_t(x)\sim \mathcal{N}_d(\mu_t,\Sigma_t)$$ where $$\mu_t=\begin{bmatrix} x_1e^{-Q_1 t} \\ x_2e^{-Q_2 t} \\ ... \\ x_de^{-Q_d t} \\ \end{bmatrix} $$ $$\Sigma_t=\begin{bmatrix} \frac{\epsilon_1^2}{2Q_1}(1-e^{-2Q_1 t}) & 0 & ... & 0 \\ 0 & \frac{\epsilon_2^2}{2Q_2}(1-e^{-2Q_2 t}) & ... & 0 \\ 0 & 0 & ... & \frac{\epsilon_d^2}{2Q_d}(1-e^{-2Q_d t}) \\ \end{bmatrix} $$ To see that the covariance matrix is like this without solving the PDE, notice that the individual processes $X_{k,t}$ can be solved as in the monodimensional case using Ito and $$\mathbb{E}[(X_{k,t}-\mu_{k,t})(X_{j,t}-\mu_{j,t})|x]=(\epsilon_k\epsilon_j)\mathbb{E}\bigg[\bigg(\int_0^te^{-Q_k(t-s)}dW_{k,s}\bigg)\bigg(\int_0^te^{-Q_j(t-s)}dW_{j,s}\bigg)\bigg|x\bigg]$$ where, if $k \neq j$ then $W_{k,t}$ and $W_{j,t}$ are independent Brownian motions, which implies $\mathbb{E}[(X_{k,t}-\mu_{k,t})(X_{j,t}-\mu_{j,t})|x]=0$. However, if $k=j$ then $$\mathbb{E}[(X_{k,t}-\mu_{k,t})^2|x]=(\epsilon_k)^2\mathbb{E}\bigg[\bigg(\int_0^te^{-Q_k(t-s)}dW_{k,s}\bigg)^2\bigg|x\bigg]$$ Using Ito isometry: $$\mathbb{E}[(X_{k,t}-\mu_{k,t})^2|x]=(\epsilon_k)^2\int_0^te^{-2Q_k(t-s)}ds=\frac{\epsilon_k^2}{2Q_k}(1-e^{-2Q_k t})$$
If $Q$ is a $d \times d$ deterministic matrix, the SDE system becomes
$$dX_{k,t}(x)=\sum_{j=1}^dQ_{k,j}X_{j,t}dt+\epsilon_kdW_{k,t}, \ \ \ k = \{1,2,...,d\}$$
the PDE of the density becomes $$\frac{\partial p(y,t)}{\partial t}=-\sum_{k=1}^d\frac{\partial}{\partial y_k}\bigg(p(y,t)\sum_{j=1}^dQ_{k,j}y_{j}\bigg)+\sum_{k=1}^d\frac{\epsilon_k^2}{2}\frac{\partial^2 p(y,t)}{\partial y_k^2}, \ \ \ y \in \mathbb{R}^d$$
and the solution depends on the entries of $Q$. In the monodimensional case or 'vector' case $Q=[Q_1,Q_2,...,Q_d]$ $$dX_{k,t}(x)=-Q_kX_kdt+\epsilon_kdW_{k,t}$$ it is sufficient to set $Q_k>0$ to obtain an O-U system with a stable Gaussian distribution for $t \to \infty$.