[proofplan] We compute the conditional density by dividing the joint Gaussian density of $(X_1,X_2)$ by the marginal density of $X_2$. The main algebraic step is the Schur-complement factorization of $\Sigma$, which gives both the determinant identity and the completed-square form of the quadratic exponent. Once the joint density separates into a normal density in $x_1$ times a normal density in $x_2$, the quotient is exactly the stated multivariate normal density. [/proofplan]

[step:Show that the Schur complement is positive definite]Since $\Sigma$ is symmetric positive definite, the principal block $\Sigma_{22}$ is symmetric positive definite. Indeed, for every nonzero vector $v \in \mathbb{R}^{p_2}$, \begin{align*} v^\top \Sigma_{22} v = \begin{pmatrix} 0 \\ v \end{pmatrix}^{\top} \Sigma \begin{pmatrix} 0 \\ v \end{pmatrix} >0. \end{align*} Thus $\Sigma_{22}$ is invertible. Define the Schur complement matrix $S \in \mathbb{R}^{p_1\times p_1}$ by \begin{align*} S := \Sigma_{11}-\Sigma_{12}\Sigma_{22}^{-1}\Sigma_{21}. \end{align*} For any nonzero vector $u \in \mathbb{R}^{p_1}$, define \begin{align*} v := -\Sigma_{22}^{-1}\Sigma_{21}u \in \mathbb{R}^{p_2}. \end{align*} Then \begin{align*} \begin{pmatrix} u \\ v \end{pmatrix}^{\top} \Sigma \begin{pmatrix} u \\ v \end{pmatrix} &= u^\top \Sigma_{11}u+u^\top\Sigma_{12}v+v^\top\Sigma_{21}u+v^\top\Sigma_{22}v \\ &= u^\top \Sigma_{11}u-u^\top\Sigma_{12}\Sigma_{22}^{-1}\Sigma_{21}u -u^\top\Sigma_{12}\Sigma_{22}^{-1}\Sigma_{21}u +u^\top\Sigma_{12}\Sigma_{22}^{-1}\Sigma_{21}u \\ &= u^\top S u. \end{align*} Since $u \neq 0$, the vector $\begin{pmatrix}u\\v\end{pmatrix}$ is nonzero, so positive definiteness of $\Sigma$ gives $u^\top S u>0$. Therefore $S=\Sigma_{1\mid 2}$ is symmetric positive definite.[/step]

[guided]The conditional covariance matrix is the Schur complement of $\Sigma_{22}$ in $\Sigma$, so we first verify that it is a legitimate covariance matrix. Because $\Sigma$ is symmetric positive definite, every principal covariance block is positive definite. For $v \in \mathbb{R}^{p_2}$ with $v \neq 0$, insert $v$ into the second block and zero into the first block: \begin{align*} v^\top \Sigma_{22}v = \begin{pmatrix}0\\v\end{pmatrix}^{\top} \Sigma \begin{pmatrix}0\\v\end{pmatrix} >0. \end{align*} Hence $\Sigma_{22}$ is symmetric positive definite and therefore invertible. Now define \begin{align*} S := \Sigma_{11}-\Sigma_{12}\Sigma_{22}^{-1}\Sigma_{21}. \end{align*} To prove that $S$ is positive definite, fix a nonzero vector $u \in \mathbb{R}^{p_1}$. The correct second block to pair with $u$ is the value that cancels the mixed terms, namely \begin{align*} v := -\Sigma_{22}^{-1}\Sigma_{21}u. \end{align*} Substituting the vector $\begin{pmatrix}u\\v\end{pmatrix}$ into the quadratic form of $\Sigma$ gives \begin{align*} \begin{pmatrix} u \\ v \end{pmatrix}^{\top} \Sigma \begin{pmatrix} u \\ v \end{pmatrix} &= u^\top \Sigma_{11}u+u^\top\Sigma_{12}v+v^\top\Sigma_{21}u+v^\top\Sigma_{22}v \\ &= u^\top \Sigma_{11}u-u^\top\Sigma_{12}\Sigma_{22}^{-1}\Sigma_{21}u -u^\top\Sigma_{12}\Sigma_{22}^{-1}\Sigma_{21}u +u^\top\Sigma_{12}\Sigma_{22}^{-1}\Sigma_{21}u \\ &= u^\top S u. \end{align*} The vector $\begin{pmatrix}u\\v\end{pmatrix}$ is nonzero because its first block is $u \neq 0$. Since $\Sigma$ is positive definite, the left-hand side is strictly positive. Hence $u^\top S u>0$ for every nonzero $u \in \mathbb{R}^{p_1}$, which proves that $S=\Sigma_{1\mid 2}$ is positive definite.[/guided]

[step:Factor the covariance matrix and compute its determinant] Define the matrix $A \in \mathbb{R}^{p_1\times p_2}$ by \begin{align*} A := \Sigma_{12}\Sigma_{22}^{-1}. \end{align*} Since $\Sigma$ is symmetric, $\Sigma_{21}=\Sigma_{12}^\top$; since $\Sigma_{22}$ is symmetric positive definite, $\Sigma_{22}^{-1}$ is symmetric. Hence \begin{align*} A^\top = \Sigma_{22}^{-1}\Sigma_{21}. \end{align*} Then \begin{align*} \Sigma &= \begin{pmatrix} I_{p_1} & A \\ 0 & I_{p_2} \end{pmatrix} \begin{pmatrix} S & 0 \\ 0 & \Sigma_{22} \end{pmatrix} \begin{pmatrix} I_{p_1} & 0 \\ A^\top & I_{p_2} \end{pmatrix}. \end{align*} Taking determinants and using that both triangular factors have determinant $1$, we get \begin{align*} \det \Sigma = \det S \, \det \Sigma_{22}. \end{align*} [/step]

[step:Complete the square in the joint Gaussian density]For $m \in \{p,p_1,p_2\}$, let $\mathcal{L}^m$ denote $m$-dimensional Lebesgue measure on $\mathbb{R}^m$. Let $f_X: \mathbb{R}^{p} \to [0,\infty)$ be the density of $X$ with respect to $\mathcal{L}^{p}$. For $x_1 \in \mathbb{R}^{p_1}$ and $x_2 \in \mathbb{R}^{p_2}$, define \begin{align*} y_1 &:= x_1-\mu_1, \\ y_2 &:= x_2-\mu_2. \end{align*} Using the factorization of $\Sigma$ from the previous step, the inverse quadratic form is \begin{align*} \begin{pmatrix} y_1 \\ y_2 \end{pmatrix}^{\top} \Sigma^{-1} \begin{pmatrix} y_1 \\ y_2 \end{pmatrix} = (y_1-Ay_2)^\top S^{-1}(y_1-Ay_2)+y_2^\top\Sigma_{22}^{-1}y_2. \end{align*} Therefore the joint density is \begin{align*} f_X(x_1,x_2) &= \frac{1}{(2\pi)^{p/2}(\det\Sigma)^{1/2}} \exp\left( -\frac{1}{2} \begin{pmatrix} y_1 \\ y_2 \end{pmatrix}^{\top} \Sigma^{-1} \begin{pmatrix} y_1 \\ y_2 \end{pmatrix} \right) \\ &= \frac{1}{(2\pi)^{p_1/2}(\det S)^{1/2}} \exp\left( -\frac{1}{2}(x_1-\mu_1-A(x_2-\mu_2))^\top S^{-1} (x_1-\mu_1-A(x_2-\mu_2)) \right) \\ &\qquad\cdot \frac{1}{(2\pi)^{p_2/2}(\det\Sigma_{22})^{1/2}} \exp\left( -\frac{1}{2}(x_2-\mu_2)^\top\Sigma_{22}^{-1}(x_2-\mu_2) \right). \end{align*}[/step]

[guided]The goal is to isolate all dependence on $x_1$ into one Gaussian density. Define \begin{align*} y_1 &:= x_1-\mu_1, \\ y_2 &:= x_2-\mu_2, \end{align*} so that the exponent in the joint density is the quadratic form \begin{align*} \begin{pmatrix} y_1 \\ y_2 \end{pmatrix}^{\top} \Sigma^{-1} \begin{pmatrix} y_1 \\ y_2 \end{pmatrix}. \end{align*} The factorization \begin{align*} \Sigma &= \begin{pmatrix} I_{p_1} & A \\ 0 & I_{p_2} \end{pmatrix} \begin{pmatrix} S & 0 \\ 0 & \Sigma_{22} \end{pmatrix} \begin{pmatrix} I_{p_1} & 0 \\ A^\top & I_{p_2} \end{pmatrix} \end{align*} gives an explicit inverse by inverting the three factors in reverse order: \begin{align*} \Sigma^{-1} &= \begin{pmatrix} I_{p_1} & 0 \\ -A^\top & I_{p_2} \end{pmatrix} \begin{pmatrix} S^{-1} & 0 \\ 0 & \Sigma_{22}^{-1} \end{pmatrix} \begin{pmatrix} I_{p_1} & -A \\ 0 & I_{p_2} \end{pmatrix}. \end{align*} The rightmost factor sends $\begin{pmatrix}y_1\\y_2\end{pmatrix}$ to $\begin{pmatrix}y_1-Ay_2\\y_2\end{pmatrix}$. Therefore the quadratic form becomes \begin{align*} \begin{pmatrix} y_1 \\ y_2 \end{pmatrix}^{\top} \Sigma^{-1} \begin{pmatrix} y_1 \\ y_2 \end{pmatrix} = (y_1-Ay_2)^\top S^{-1}(y_1-Ay_2)+y_2^\top\Sigma_{22}^{-1}y_2. \end{align*} This is the completion of the square: the first term is a quadratic form in \begin{align*} x_1-\mu_1-A(x_2-\mu_2), \end{align*} and the second term depends only on $x_2$. Using the determinant identity $\det\Sigma=\det S\,\det\Sigma_{22}$ and $p=p_1+p_2$, the normalizing constant also splits: \begin{align*} \frac{1}{(2\pi)^{p/2}(\det\Sigma)^{1/2}} = \frac{1}{(2\pi)^{p_1/2}(\det S)^{1/2}} \frac{1}{(2\pi)^{p_2/2}(\det\Sigma_{22})^{1/2}}. \end{align*} Substituting the completed square and the split normalizing constant into the joint density gives \begin{align*} f_X(x_1,x_2) &= \frac{1}{(2\pi)^{p_1/2}(\det S)^{1/2}} \exp\left( -\frac{1}{2}(x_1-\mu_1-A(x_2-\mu_2))^\top S^{-1} (x_1-\mu_1-A(x_2-\mu_2)) \right) \\ &\qquad\cdot \frac{1}{(2\pi)^{p_2/2}(\det\Sigma_{22})^{1/2}} \exp\left( -\frac{1}{2}(x_2-\mu_2)^\top\Sigma_{22}^{-1}(x_2-\mu_2) \right). \end{align*}[/guided]

[step:Identify the marginal density of $X_2$] Define $f_{X_2}: \mathbb{R}^{p_2} \to [0,\infty)$ by \begin{align*} f_{X_2}(x_2) := \int_{\mathbb{R}^{p_1}} f_X(x_1,x_2)\,d\mathcal{L}^{p_1}(x_1). \end{align*} From the product form above and the normalization of the $\mathcal{N}_{p_1}(\mu_1+A(x_2-\mu_2),S)$ density, \begin{align*} f_{X_2}(x_2) = \frac{1}{(2\pi)^{p_2/2}(\det\Sigma_{22})^{1/2}} \exp\left( -\frac{1}{2}(x_2-\mu_2)^\top\Sigma_{22}^{-1}(x_2-\mu_2) \right). \end{align*} Thus $X_2 \sim \mathcal{N}_{p_2}(\mu_2,\Sigma_{22})$. [/step]

[step:Divide the joint density by the marginal density] For each fixed $x_2 \in \mathbb{R}^{p_2}$, the marginal density $f_{X_2}(x_2)$ is strictly positive. Define the map $q: \mathbb{R}^{p_1}\times\mathbb{R}^{p_2}\to[0,\infty)$ by \begin{align*} q(x_1,x_2) := \frac{f_X(x_1,x_2)}{f_{X_2}(x_2)}. \end{align*} This quotient defines a regular conditional density of $X_1$ given $X_2$. Indeed, for every Borel set $B \subset \mathbb{R}^{p_1}$ and every Borel set $C \subset \mathbb{R}^{p_2}$, Tonelli's theorem applies because the integrand is non-negative, and the definition of $f_{X_2}$ gives \begin{align*} \mathbb{P}(X_1\in B, X_2\in C) &= \int_C \int_B f_X(x_1,x_2)\,d\mathcal{L}^{p_1}(x_1)\,d\mathcal{L}^{p_2}(x_2) \\ &= \int_C \int_B q(x_1,x_2)\,d\mathcal{L}^{p_1}(x_1)\, f_{X_2}(x_2)\,d\mathcal{L}^{p_2}(x_2). \end{align*} Thus the probability measure $B\mapsto \int_B q(x_1,x_2)\,d\mathcal{L}^{p_1}(x_1)$ is a version of the conditional law of $X_1$ given $X_2=x_2$. Using the product decomposition of $f_X$ and the formula for $f_{X_2}$, the factors depending only on $x_2$ cancel, leaving \begin{align*} q(x_1,x_2) = \frac{1}{(2\pi)^{p_1/2}(\det S)^{1/2}} \exp\left( -\frac{1}{2}(x_1-\mu_1-A(x_2-\mu_2))^\top S^{-1} (x_1-\mu_1-A(x_2-\mu_2)) \right). \end{align*} Since $A=\Sigma_{12}\Sigma_{22}^{-1}$ and $S=\Sigma_{11}-\Sigma_{12}\Sigma_{22}^{-1}\Sigma_{21}$, this is the density of \begin{align*} \mathcal{N}_{p_1}\left( \mu_1+\Sigma_{12}\Sigma_{22}^{-1}(x_2-\mu_2), \Sigma_{11}-\Sigma_{12}\Sigma_{22}^{-1}\Sigma_{21} \right). \end{align*} Therefore, for every $x_2 \in \mathbb{R}^{p_2}$, \begin{align*} X_1 \mid X_2=x_2 \sim \mathcal{N}_{p_1}(\mu_{1\mid 2}(x_2),\Sigma_{1\mid 2}), \end{align*} as claimed. [/step]