Let $d_1,d_2 \in \mathbb{N}$. Let $(B_t)_{t \in [0, \infty)}$ be a standard $d_2$-dimensional Brownian Motion. Let $\mu$ be a constant ($d_1 \times d_1$) - matrix and let $ \sigma $ be a constant ($d_1 \times d_2$) - matrix.
I believe I have shown that the unique strong solution $(X_{t}^{x},B_{t})_{t \in (0,\infty]}$ to the equation
$$ dX_{t}^{x} = \mu X_{t}^{x} dt + \sigma dB_t \ \ \ \ \ \ in \ \ (0, \infty) $$ $$ X_{0}^{x} = x $$
Is given by
$$ X_{t}^{x} = exp(-t \mu) ( \int_{0}^{t} exp(s \mu ) \sigma dB_s + x) $$
The question then asks for the distribution of $X_{t}^{x}$ for each $x \in \mathbb{R}^d $ and $t \in [0, \infty)$. My thoughts are that, since this is a deterministic integrand with respect to Brownian motion, the integral must be a Gaussian vector and so its distribution will be determined entirely by its mean and covariance. But then I have trouble finding the covariance of a vector like this - any help would be much appreciated!
The covariance should be given by the multidimensional Ito-isometry
$$ C_t = \mathbb E[(X_t^x-e^{-\mu t}x)(X_t^x-e^{-\mu t}x)^T] = e^{-2\mu t}\mathbb E[\int_0^t e^{s\mu} \sigma dB_s (\int_0^t e^{s\mu}\sigma dB_s)^T] = e^{-2\mu t} \int_0^t e^{2\mu s} \sigma \sigma^T ds = \frac1{2\mu}\bigg(1-e^{-2\mu t} \bigg)\sigma \sigma^T $$