Covariance matrix after projection on gaussian state I cannot understand the proof in Eisert article (https://arxiv.org/abs/quant-ph/0204052) about finding the covariance matrix of a state after projecting some of its modes onto a gaussian state.
We have the characteristic function of the remaining (unprojected) modes being
\begin{equation}
\chi(\xi_A)=\frac{1}{\pi^2}\int d\xi_5\dots d\xi_8 e^{-\xi^T\Gamma'\xi/2}e^{-\xi_B^T \gamma\xi_B}
\end{equation}
where $\xi=(\xi_A, \xi_B)$, the labels $A$ and $B$ referring respectively to the subsystem that ''remains'' and to the subsystem that is projected to a gaussian state with covariance matrix $\gamma=\mathrm{diag}(1/d,d,1/d,d )$ and where
\begin{equation}
\Gamma'=\begin{pmatrix}C_1 & C_3\\C_3^T &C_2 \end{pmatrix}
\end{equation}
He then states that the resulting covariance matrix, after carrying out the integration, is 
\begin{equation}
M_d=C_1-C_3(C_2+\gamma^2)^{-1}C_3^T
\end{equation}
but I have no idea of how he manages to carry out the integration, since the term $\xi\Gamma'\xi$ has mixed terms in it reading $\xi_B^T C^T \xi_A$ and $\xi_A^T C \xi_B$ and I don't know how to treat them.
Any hint?
 A: As suggested by @flippiefanus, it’s a case of completing the square. (Or, alternatively, it is a well known formula for Gaussian integration, but this sounds a bit like cheating)
If I’m not mistaken, you made a typo in your first equation, which is equation (8) in the Eisert, Scheel, Plenio article. You should replace $γ$ by $\frac{γ^2}2$. We have then
\begin{align}
χ(ξ_A)=\frac1{π^2}∫dξ_B\exp(-\tfrac12ξ_A^TC_1ξ_A 
   -\tfrac12ξ_B^TC_3^Tξ_A-\tfrac12ξ_A^TC_3ξ_B-\tfrac12ξ_B^T(C_2+γ^2)ξ_B)\\
=\frac{\exp(-\tfrac12ξ_A^TC_1ξ_A )}{π^2}∫dξ_B\exp( 
   -\tfrac12ξ_B^TC_3^Tξ_A-\tfrac12ξ_A^TC_3ξ_B-\tfrac12ξ_B^T(C_2+γ^2)ξ_B).
\end{align}
The idea is to express the argument of the integrand as $-\tfrac12(ξ_B+Δ)^TM(ξ_B+Δ)+\tfrac12Δ^TMΔ$, the “$+Δ$” translation of $\xi_B$ being irrelevant in the integration. This expression becomes
$$-\tfrac12ξ_B^T M Δ - \tfrac12 Δ^T M ξ_B  - \tfrac12 ξ_B^T M ξ_B $$
which, when identified with the argument of the exponential, leads to
\begin{align}
  M& =C_2+γ^2 & M\Delta&=C_3^Tξ_A & Δ^T&=ξ_A^TC_3M^{-1}\\
  && && &=ξ_A^TC_3(C_2+γ^2)^{-1}
\end{align} 
We have then
\begin{align}
χ(ξ_A)&=\frac{\exp(-\tfrac12ξ_A^TC_1ξ_A )}{π^2}
     ∫dξ_B\exp( -\tfrac12(ξ_B+Δ)^TM(ξ_B+Δ)+\tfrac12Δ^TMΔ)\\
 &=\frac{\exp(-\tfrac12ξ_A^TC_1ξ_A +\tfrac12Δ^TMΔ)}{π^2}
    ∫dξ_B\exp( -\tfrac12(ξ_B+Δ)^TM(ξ_B+Δ))\\
 &=\frac{\exp(-\tfrac12ξ_A^TC_1ξ_A +\tfrac12Δ^TMΔ)}{π^2}
    \underbrace{∫dξ_B\exp( -\tfrac12ξ_B^TMξ_B)}_{\text{scalar independent of $ξ_A$}}\\
 &∝\exp(-\tfrac12ξ_A^TC_1ξ_A +\tfrac12ξ_A^TC_3(C_2+γ^2)^{-1}C_3^Tξ_A)\\
 &=\exp(-\tfrac12ξ_A^T
       \underbrace{\left(C_1+C_3(C_2+γ^2)^{-1}C_3^T\right)}_{M_d}
       ξ_A)
\end{align} 
