Joint probability distribution of the Ornstein-Uhlenbeck process

The best known stationary Markov process is the Ornstein-Uhlenbeck process fully characterized via:
$$
p_1 (x_1) = frac{1}{sqrt{2pi}} e^{-frac{1}{2}x_1^{2}}, quad p_{1|1}(x_2|x_1) = frac{1}{sqrt{2pi (1 – e^{-2tau})}} exp left[ -frac{(x_2 – x_1 e^{-tau})^2} {2(1-e^{-2tau})} right],
$$
$p_1$ is the probability to obtain the value $x_1$ at time $t_1$, while $p_{1|1}$ is the transtion probability, i.e., the conditioned probability of obataining $x_2$ at $t_2$ given $x_1$ at $t_1$. The transition probability depends only on the time difference $tau = t_2 – t_1 > 0$. This process was originally derived to determine the velocity of a Brownian particle and is essentially the only process which is Gaussian, Markovian and stationary.

How to calculate the joint probability distribution $p(x_1,x_2)$ of the Ornstein-Uhlenbeck process?

-> How to calculate the autocorrelation function of the Ornstein-Uhlenbeck process via
$$
langle langle x_1x_2ranglerangle = int dx_1 int dx_2 x_1 x_2 p(x_1x_2)
$$

Joint probability distribution $p(x_1,x_2)$ of the Ornstein-Uhlenbeck process. Using the definition,
$$
p(x_1, x_2) = p(x_2|x_1) p(x_1)
= frac{1}{sqrt{2pi (1 – e^{-2tau})}} exp left[ -frac{(x_2 – x_1 e^{-tau})^2} {2(1-e^{-2tau})} right] . frac{1}{sqrt{2pi}} e^{-frac{1}{2}x_1^{2}}
= frac{1}{2pi sqrt{(1 – e^{-2tau})}} exp left[ -frac{1}{2} left( frac{(x_2 – x_1 e^{-tau})^2} {2(1-e^{-2tau})} – x_1^2 right) right]
$$

How to proceed further or this is the final solution??