the probability density function (PDF) of concatenation of two Gaussian variables Gaussian variable $x$ follows from $N(u_x,\sigma_x^2)$ and $y$ follows from $N(u_y,\sigma_y^2)$. Assume we have the vector $\bf{z}=[x,y]^T\in R^2$, then it seems that no matter whether $x$ and $y$ are independent or not, we always have that $\bf{z}$ also follows from the Gaussian distribution $N([u_x,u_y]^T, Cov([x,y]^T))$, where $Cov$ means the covariance. 
The above claim is reformulated from the last eight lines on the left column in page 4 of http://www.cs.bham.ac.uk/~axk/ICML_Flip_2013.pdf .
Could anyone show why the distribution of $\bf{z}$ is Gaussian? and how to get the related parameters?
 A: $Z=(X,Y)'$ is not necessary bivariate normal (it works if $X$ and $Y$ are independent). There are lots of counterexamples (e.g. link). Another example: $\psi_1$ and $\psi_2$ are independent standard normal r.v.s and
$$(X, Y)=(\psi_1,|\psi_2|)1\{\psi_1\ge 0\}+(\psi_1,-|\psi_2|)1\{\psi_1< 0\}$$
which is not bivariate normal. However, marginals $X$ and $Y$ are normal. 
On the other hand, $Z$ is a bivariate normal (multivariate, in general) if $Z=A\cdot X+\mu$ for some matrix $A: 2\times k$, a vector of independent standard normal r.v.s. $X=(X_1,\dots,X_k)'$, and a vector $\mu=(\mu_1,\mu_2)'$.
Edit: In your case $r_{ij}/\sigma\sim_{iid}\mathcal{N}(0, 1)$ and
$$\begin{bmatrix}u_i\\v_i\end{bmatrix}=
\begin{bmatrix}
h_{1} && h_{2} && \dots && h_{d} \\
x_{1} && x_{2} && \dots && x_{d}
\end{bmatrix}\cdot
\begin{bmatrix}r_{i1} \\ r_{i2} \\ \vdots \\ r_{id}
\end{bmatrix}
=\sigma^2A\cdot \frac{R_i}{\sigma^2}$$ 
s.t.
$$Cov(u_i,v_i)=\sigma^2\times AA'= \sigma^2\times
\begin{bmatrix} 
h'h && h'x \\
x'h && x'x
\end{bmatrix}
$$
