For two multivariate normal vectors $\mathbf{x}\sim\mathcal N(\mathbf{\mu_x},\Sigma_{\mathbf{x,x}})$ and $\mathbf{Y}\sim\mathcal N(\mathbf{\mu_Y},\Sigma_{\mathbf{Y,Y}})$ we have that, given their covariance $\Sigma_{\mathbf{x,Y}}$:
$$E[\mathbf x|\mathbf Y]=\mathbf{\mu_x}+\Sigma_{\mathbf{x,Y}}\Sigma_{\mathbf{Y,Y}}^{-1}(\mathbf{Y-\mu_Y})$$
$$\text{Var}[\mathbf x|\mathbf Y]=\Sigma_{\mathbf{x,x}}-\Sigma_{\mathbf{x,Y}}\Sigma_{\mathbf{Y,Y}}^{-1}\Sigma_{\mathbf{Y,x}}$$
But the covariance is given in terms of the second moment:
$$\text{Var}[\mathbf x|\mathbf Y]=E[\mathbf {xx}^T|\mathbf Y]-E[\mathbf x|\mathbf Y]E[\mathbf {x}^T|\mathbf Y]$$
So
$$E[\mathbf {xx}^T|\mathbf Y] = \text{Var}[\mathbf x|\mathbf Y] + E[\mathbf x|\mathbf Y]E[\mathbf {x}^T|\mathbf Y]\\
=\Sigma_{\mathbf{x,x}}-\Sigma_{\mathbf{x,Y}}\Sigma_{\mathbf{Y,Y}}^{-1}\Sigma_{\mathbf{Y,x}}+
(\mathbf{\mu_x}+\Sigma_{\mathbf{x,Y}}\Sigma_{\mathbf{Y,Y}}^{-1}(\mathbf{Y-\mu_Y}))
(\mathbf{\mu_x}+\Sigma_{\mathbf{x,Y}}\Sigma_{\mathbf{Y,Y}}^{-1}(\mathbf{Y-\mu_Y}))^T$$
I haven't checked, but it seems this all stems from the simple fact that:
$$\mathbf x=\mathbf{\mu_x}+\Sigma_{\mathbf{x,Y}}\Sigma_{\mathbf{Y,Y}}^{-1}(\mathbf{Y-\mu_Y})+\mathbf e\\
E[\mathbf {xx}^T|\mathbf Y]=(\mathbf{\mu_x}+\Sigma_{\mathbf{x,Y}}\Sigma_{\mathbf{Y,Y}}^{-1}(\mathbf{Y-\mu_Y}))(\mathbf{\mu_x}+\Sigma_{\mathbf{x,Y}}\Sigma_{\mathbf{Y,Y}}^{-1}(\mathbf{Y-\mu_Y}))^T+E[\mathbf {ee}^T|\mathbf Y]
$$
So we could've started directly here, and then $E[\mathbf {ee}^T|\mathbf Y] = \text{Var}[\mathbf x|\mathbf Y]$.