Brownian motion - characteristic function Let me remind first the construction of Brownian motion.
Fix a vector $x \in \mathbb{R}^n$ and define
$p(t,x,y) := (2\pi t )^{-\frac{n}{2}} \cdot \exp{\left( - \frac{|x-y|^2}{2t} \right)},$ for $y \in \mathbb{R}^n$, $t >0$.
Then for $0 \leq t_1 \leq t_2 \leq \ldots \leq t_k$ we define a measure $\nu_{t_1, \ldots, t_k}$ on $\mathbb{R}^{nk}$ by
$$\nu_{t_1, \ldots, t_k}(F_1 \times \ldots \times F_k)= \int \limits_{F_1 \times \ldots \times F_k}^{} p(t_1, x, x_1)\prod_{j=1}^{k-1}p(t_{j+1}-t_j,x_{j}, x_{j+1}) dx_1\ldots dx_k,$$ where
the following conventions are used $dy=dy_1\ldots dy_k$ for Lebesgue measure and $p(0, x, y)dy=\delta_x(y)$ (the unit point mass at $x$).
By Kolmogorov's extension theorem applied to the probability measures $\nu_{t_1, \ldots, t_k}$ (which easily satisfies all the assumtions of that theorem) there exists a probability space $(\Omega, \mathcal{F}, P^x)$ and a stochastic process $\{B_t\}_{t \geq 0}$ on $\Omega$ such that the finite distributions of $B_t$ are given by
$$ (\star) \ \ P^x(B_{t_1} \in F_1, \ldots, B_{t_k} \in F_k)= \int \limits_{F_1 \times \ldots \times F_k}^{} p(t_1, x, x_1)\prod_{j=1}^{k-1}p(t_{j+1}-t_j,x_{j}, x_{j+1}) dx_1\ldots dx_k.$$
Problem 
I want to show that for the random variable $Z = (B_{t_1}, \ldots, B_{t_k}) \in \mathbb{R}^{nk}$ there exist a vector $M \in \mathbb{R}^{nk}$ and a non-negative definite matrix $C \in \text{M}_{nk \times nk}(\mathbb{R})$ such that
$$E^x\left[\exp\left(i\left<u, Z \right> \right)\right]= \exp\left( -\frac{1}{2} \left<Cu, u\right> + i \left<u, M\right> \right),$$
for all $u=(u_1, \ldots, u_{nk}) \in \mathbb{R}^{nk}$ (left hand side stands for the characteristic function of the random variable $Z$). Moreover, $M$ is the mean value of $Z$ and $C$ is a covariance matrix of $Z$.
I was trying to calculate it explicitly by writing the left hand side by its definition and applying the ($\star$) formula, but the integral which I've got is not so nice and I cannot conclude what I want to.
This exercise is nothing else, but showing that $B_t$ is a Gaussian process (so it is not hard to guess what should be $M$ and $C$ in this case). 
I hope you know some tricks how to compute such an integral in an easy way. Thanks in advance for any help.
 A: Denote $t_0 = 0$ and $x_0 = 0 $ for notational convenience. Then
$$
   \mathbb{E}( \exp(i \langle u, Z \rangle ) = \int_{\mathbb{R}^{n k}} \exp( i \langle u, x \rangle ) \prod_{j=1}^k p(t_{j}-t_{j-1}, x_{j} - x_{j-1}) \mathrm{d}x_1 \cdots \mathrm{d} x_k
$$
Now change variables $x_k = \sum_{i=1}^k s_k$. This has unit Jacobian, you get
$$
\mathbb{E}( \exp(i \langle u, Z \rangle ) =  \int_{\mathbb{R}^{n k}}  
 \prod_{k=1}^k \exp( i \, s_j \sum_{m=j}^k u_q ) ) \prod_{j=1}^k p(t_{j}-t_{j-1}, s_j) \mathrm{d}s_1 \cdots \mathrm{d} s_k
$$
Now integration with respect to each $s_j$ can be carried out independently and will produce $\exp( Q_j(u))$, where $Q_j$ is a quadratic multivariate polynomial in $u$. Hence the characteristic function of k-dimensional time-slice distribution of this Brownian motion process is going to be $\exp(Q(u))$. 
Now $Q(u) = \langle C u, u \rangle  + i \langle u, M \rangle $, where $C$ is covariance matrix and $M$ is mean vector. The free term will be missing, because c.f. is one for zero vector $u$.
