[proofplan]
We prove both implications. Independence first forces all mixed second moments to factor, so every entry of the cross-covariance block $\Sigma_{12}$ is zero. Conversely, if $\Sigma_{12}=0$, then the characteristic function of the joint vector $(X_1,X_2)$ factors into the product of the marginal characteristic functions of $X_1$ and $X_2$. The uniqueness theorem for characteristic functions then identifies the joint law with the product of the marginal laws, which is precisely independence.
[/proofplan]
[step:Use independence to annihilate the mixed covariance block]
Assume that $X_1$ and $X_2$ are independent. For $1 \le i \le p_1$ and $1 \le j \le p_2$, let $X_{1,i}: \Omega \to \mathbb{R}$ and $X_{2,j}: \Omega \to \mathbb{R}$ denote the $i$-th and $j$-th coordinate random variables of $X_1$ and $X_2$, respectively. Since $X$ is multivariate normal, all coordinate random variables have finite second moments.
The $(i,j)$-entry of $\Sigma_{12}$ is
\begin{align*}
(\Sigma_{12})_{ij}
&= \operatorname{Cov}(X_{1,i},X_{2,j}) \\
&= \mathbb{E}\big[(X_{1,i}-\mathbb{E}[X_{1,i}])(X_{2,j}-\mathbb{E}[X_{2,j}])\big].
\end{align*}
Independence of the random vectors $X_1$ and $X_2$ implies independence of the real-valued random variables $X_{1,i}$ and $X_{2,j}$. Therefore
\begin{align*}
\mathbb{E}[X_{1,i}X_{2,j}]
=
\mathbb{E}[X_{1,i}]\,\mathbb{E}[X_{2,j}],
\end{align*}
and hence
\begin{align*}
(\Sigma_{12})_{ij}
&= \mathbb{E}[X_{1,i}X_{2,j}]
- \mathbb{E}[X_{1,i}]\,\mathbb{E}[X_{2,j}] \\
&= 0.
\end{align*}
Since this holds for every $1 \le i \le p_1$ and $1 \le j \le p_2$, we have $\Sigma_{12}=0$.
[/step]
[step:Factor the joint characteristic function when the covariance block vanishes]
Assume now that $\Sigma_{12}=0$. Since $\Sigma$ is symmetric, $\Sigma_{21}=\Sigma_{12}^{\top}=0$. Define the joint characteristic function
\begin{align*}
\phi_X: \mathbb{R}^{p_1} \times \mathbb{R}^{p_2} &\to \mathbb{C} \\
(t_1,t_2) &\mapsto \mathbb{E}\big[\exp(i(t_1^{\top}X_1+t_2^{\top}X_2))\big].
\end{align*}
For a multivariate normal vector $X \sim \mathcal{N}_p(\mu,\Sigma)$, its characteristic function is
\begin{align*}
\phi_X(t)
=
\exp\left(i t^{\top}\mu-\frac{1}{2}t^{\top}\Sigma t\right),
\qquad t \in \mathbb{R}^p.
\end{align*}
Taking $t=(t_1^\top,t_2^\top)^\top \in \mathbb{R}^p$, we compute
\begin{align*}
t^{\top}\mu
&= t_1^{\top}\mu_1+t_2^{\top}\mu_2, \\
t^{\top}\Sigma t
&=
t_1^{\top}\Sigma_{11}t_1
+t_1^{\top}\Sigma_{12}t_2
+t_2^{\top}\Sigma_{21}t_1
+t_2^{\top}\Sigma_{22}t_2 \\
&=
t_1^{\top}\Sigma_{11}t_1
+t_2^{\top}\Sigma_{22}t_2.
\end{align*}
Thus
\begin{align*}
\phi_X(t_1,t_2)
&=
\exp\left(i t_1^{\top}\mu_1-\frac{1}{2}t_1^{\top}\Sigma_{11}t_1\right)
\exp\left(i t_2^{\top}\mu_2-\frac{1}{2}t_2^{\top}\Sigma_{22}t_2\right).
\end{align*}
[guided]
The goal is to show that the joint distribution separates into two independent pieces. For Gaussian vectors, the characteristic function encodes exactly the mean vector and covariance matrix, so the block diagonal form of $\Sigma$ should force the exponent to split.
Define
\begin{align*}
\phi_X: \mathbb{R}^{p_1} \times \mathbb{R}^{p_2} &\to \mathbb{C} \\
(t_1,t_2) &\mapsto \mathbb{E}\big[\exp(i(t_1^{\top}X_1+t_2^{\top}X_2))\big].
\end{align*}
The characteristic function formula for $X \sim \mathcal{N}_p(\mu,\Sigma)$ gives
\begin{align*}
\phi_X(t)
=
\exp\left(i t^{\top}\mu-\frac{1}{2}t^{\top}\Sigma t\right),
\qquad t \in \mathbb{R}^p.
\end{align*}
We now insert the partitioned vector $t=(t_1^\top,t_2^\top)^\top$. The linear part separates as
\begin{align*}
t^{\top}\mu
=
t_1^{\top}\mu_1+t_2^{\top}\mu_2.
\end{align*}
For the quadratic part, block multiplication gives
\begin{align*}
t^{\top}\Sigma t
&=
t_1^{\top}\Sigma_{11}t_1
+t_1^{\top}\Sigma_{12}t_2
+t_2^{\top}\Sigma_{21}t_1
+t_2^{\top}\Sigma_{22}t_2.
\end{align*}
The hypothesis $\Sigma_{12}=0$ removes the mixed term $t_1^{\top}\Sigma_{12}t_2$. Since $\Sigma$ is symmetric, $\Sigma_{21}=\Sigma_{12}^{\top}=0$, so the other mixed term also vanishes. Hence
\begin{align*}
t^{\top}\Sigma t
=
t_1^{\top}\Sigma_{11}t_1
+t_2^{\top}\Sigma_{22}t_2.
\end{align*}
Substituting these two identities into the Gaussian characteristic function yields
\begin{align*}
\phi_X(t_1,t_2)
&=
\exp\left(i t_1^{\top}\mu_1-\frac{1}{2}t_1^{\top}\Sigma_{11}t_1\right)
\exp\left(i t_2^{\top}\mu_2-\frac{1}{2}t_2^{\top}\Sigma_{22}t_2\right).
\end{align*}
This is the desired factorisation: one factor depends only on $t_1$, and the other depends only on $t_2$.
[/guided]
[/step]
[step:Identify the two factors as the marginal characteristic functions]
Define the marginal characteristic functions
\begin{align*}
\phi_{X_1}: \mathbb{R}^{p_1} &\to \mathbb{C} \\
t_1 &\mapsto \mathbb{E}\big[\exp(i t_1^{\top}X_1)\big],
\end{align*}
and
\begin{align*}
\phi_{X_2}: \mathbb{R}^{p_2} &\to \mathbb{C} \\
t_2 &\mapsto \mathbb{E}\big[\exp(i t_2^{\top}X_2)\big].
\end{align*}
The marginal distributions of a multivariate normal vector are multivariate normal, so
\begin{align*}
X_1 \sim \mathcal{N}_{p_1}(\mu_1,\Sigma_{11}),
\qquad
X_2 \sim \mathcal{N}_{p_2}(\mu_2,\Sigma_{22}).
\end{align*}
Therefore, for every $t_1 \in \mathbb{R}^{p_1}$ and $t_2 \in \mathbb{R}^{p_2}$,
\begin{align*}
\phi_{X_1}(t_1)
&=
\exp\left(i t_1^{\top}\mu_1-\frac{1}{2}t_1^{\top}\Sigma_{11}t_1\right), \\
\phi_{X_2}(t_2)
&=
\exp\left(i t_2^{\top}\mu_2-\frac{1}{2}t_2^{\top}\Sigma_{22}t_2\right).
\end{align*}
Combining this with the previous step gives
\begin{align*}
\phi_X(t_1,t_2)=\phi_{X_1}(t_1)\phi_{X_2}(t_2),
\qquad
(t_1,t_2)\in \mathbb{R}^{p_1}\times\mathbb{R}^{p_2}.
\end{align*}
[/step]
[step:Convert characteristic function factorisation into independence]
Let $\nu$ denote the joint law of $(X_1,X_2)$ on $\mathbb{R}^{p_1}\times\mathbb{R}^{p_2}$, let $\nu_1$ denote the law of $X_1$ on $\mathbb{R}^{p_1}$, and let $\nu_2$ denote the law of $X_2$ on $\mathbb{R}^{p_2}$. The characteristic function of the product measure $\nu_1\otimes\nu_2$ is
\begin{align*}
(t_1,t_2)
&\mapsto
\int_{\mathbb{R}^{p_1}\times\mathbb{R}^{p_2}}
\exp(i(t_1^{\top}x_1+t_2^{\top}x_2))\,d(\nu_1\otimes\nu_2)(x_1,x_2) \\
&=
\left(\int_{\mathbb{R}^{p_1}}\exp(i t_1^{\top}x_1)\,d\nu_1(x_1)\right)
\left(\int_{\mathbb{R}^{p_2}}\exp(i t_2^{\top}x_2)\,d\nu_2(x_2)\right) \\
&=
\phi_{X_1}(t_1)\phi_{X_2}(t_2).
\end{align*}
By the factorisation established above, this equals the characteristic function of $\nu$ for every $(t_1,t_2)\in\mathbb{R}^{p_1}\times\mathbb{R}^{p_2}$. By the uniqueness theorem for characteristic functions, $\nu=\nu_1\otimes\nu_2$. Therefore the joint law of $(X_1,X_2)$ is the product of its marginal laws, which is exactly the independence of $X_1$ and $X_2$.
Combining this implication with the first step proves that $X_1$ and $X_2$ are independent if and only if $\Sigma_{12}=0$.
[/step]