[proofplan] We use the defining outer-product representation of the central Wishart distribution. First represent each Wishart law as the law of a sum of independent products $XX^\top$ with $X \sim \mathcal{N}_p(0,\Sigma)$. Independence of $W_1$ and $W_2$ lets us replace their joint law by the product law of two independent Gaussian samples, and the sum then becomes the Wishart construction with sample size $n_1+n_2$. [/proofplan]

[step:Represent each Wishart law on an auxiliary Gaussian product space]Because the laws $W_p(n_1,\Sigma)$ and $W_p(n_2,\Sigma)$ are assumed to be defined, the covariance matrix $\Sigma$ is symmetric and positive semidefinite, so the Gaussian law $\mathcal{N}_p(0,\Sigma)$ on $(\mathbb{R}^p,\mathcal{B}(\mathbb{R}^p))$ is well-defined. Choose an auxiliary probability space $(\Omega,\mathcal{F},\mathbb{P})$ carrying random vectors \begin{align*} X_i : (\Omega,\mathcal{F},\mathbb{P}) &\to (\mathbb{R}^p,\mathcal{B}(\mathbb{R}^p)), \qquad 1 \le i \le n_1, \end{align*} and \begin{align*} Y_j : (\Omega,\mathcal{F},\mathbb{P}) &\to (\mathbb{R}^p,\mathcal{B}(\mathbb{R}^p)), \qquad 1 \le j \le n_2, \end{align*} such that the combined family $(X_1,\dots,X_{n_1},Y_1,\dots,Y_{n_2})$ is independent and every member has distribution $\mathcal{N}_p(0,\Sigma)$. Define the random matrices \begin{align*} S_X : (\Omega,\mathcal{F},\mathbb{P}) &\to (\mathbb{R}^{p \times p},\mathcal{B}(\mathbb{R}^{p \times p})) \\ \omega &\mapsto \sum_{i=1}^{n_1} X_i(\omega)X_i(\omega)^\top \end{align*} and \begin{align*} S_Y : (\Omega,\mathcal{F},\mathbb{P}) &\to (\mathbb{R}^{p \times p},\mathcal{B}(\mathbb{R}^{p \times p})) \\ \omega &\mapsto \sum_{j=1}^{n_2} Y_j(\omega)Y_j(\omega)^\top. \end{align*} The maps $S_X$ and $S_Y$ are independent because $S_X$ is a measurable function of $(X_i)_{i=1}^{n_1}$, $S_Y$ is a measurable function of $(Y_j)_{j=1}^{n_2}$, and the two Gaussian families are independent. By the definition of the central Wishart distribution, \begin{align*} S_X \sim W_p(n_1,\Sigma), \qquad S_Y \sim W_p(n_2,\Sigma). \end{align*} Therefore $W_1$ and $S_X$ have the same distribution, and $W_2$ and $S_Y$ have the same distribution.[/step]

[guided]The central Wishart distribution $W_p(n,\Sigma)$ is defined as the law of a sum of $n$ independent Gaussian outer products. The original random matrices $W_1$ and $W_2$ may live on an arbitrary probability space, and that space need not contain any Gaussian sample. Therefore we do not add Gaussian variables to the original space. Instead, we build an auxiliary model space whose random matrices have the required laws. Because the laws $W_p(n_1,\Sigma)$ and $W_p(n_2,\Sigma)$ are assumed to be defined, the covariance matrix $\Sigma$ is symmetric and positive semidefinite. Hence the Gaussian probability law $\mathcal{N}_p(0,\Sigma)$ on $(\mathbb{R}^p,\mathcal{B}(\mathbb{R}^p))$ is well-defined. For example, one may take the product probability space with $n_1+n_2$ coordinate maps, each distributed as $\mathcal{N}_p(0,\Sigma)$. Thus we may choose an auxiliary probability space $(\Omega,\mathcal{F},\mathbb{P})$ carrying random vectors $X_1,\dots,X_{n_1},Y_1,\dots,Y_{n_2}$ such that the combined family is independent and each random vector has distribution $\mathcal{N}_p(0,\Sigma)$. Define \begin{align*} S_X : (\Omega,\mathcal{F},\mathbb{P}) &\to (\mathbb{R}^{p \times p},\mathcal{B}(\mathbb{R}^{p \times p})) \\ \omega &\mapsto \sum_{i=1}^{n_1} X_i(\omega)X_i(\omega)^\top. \end{align*} By the defining representation of $W_p(n_1,\Sigma)$, this gives $S_X \sim W_p(n_1,\Sigma)$. Define also \begin{align*} S_Y : (\Omega,\mathcal{F},\mathbb{P}) &\to (\mathbb{R}^{p \times p},\mathcal{B}(\mathbb{R}^{p \times p})) \\ \omega &\mapsto \sum_{j=1}^{n_2} Y_j(\omega)Y_j(\omega)^\top. \end{align*} Again by the defining representation of the Wishart distribution, \begin{align*} S_Y \sim W_p(n_2,\Sigma). \end{align*} The random matrix $S_X$ is a measurable function of the Gaussian family $(X_i)_{i=1}^{n_1}$, and $S_Y$ is a measurable function of the Gaussian family $(Y_j)_{j=1}^{n_2}$. Since these two families are independent, $S_X$ and $S_Y$ are independent. Hence $S_X$ is a distributional copy of $W_1$, $S_Y$ is a distributional copy of $W_2$, and the pair $(S_X,S_Y)$ has the same product-law structure as the independent pair $(W_1,W_2)$.[/guided]

[step:Transfer the joint distribution to the model sums]For a random matrix $R : (\Omega_R,\mathcal{F}_R,\mathbb{P}_R) \to (\mathbb{R}^{p \times p},\mathcal{B}(\mathbb{R}^{p \times p}))$, write $\mathcal{L}(R)$ for its law, meaning the pushforward probability measure $\mathbb{P}_R \circ R^{-1}$ on $\mathbb{R}^{p \times p}$. If $\mu$ and $\nu$ are probability measures on $\mathbb{R}^{p \times p}$, write $\mu \otimes \nu$ for their product probability measure on $\mathbb{R}^{p \times p} \times \mathbb{R}^{p \times p}$. Define the addition map \begin{align*} A : \mathbb{R}^{p \times p} \times \mathbb{R}^{p \times p} &\to \mathbb{R}^{p \times p} \\ (M,N) &\mapsto M+N. \end{align*} Since $W_1$ and $W_2$ are independent, the law of $(W_1,W_2)$ is the product of the laws of $W_1$ and $W_2$. Since $S_X$ and $S_Y$ are independent and have the same marginal laws as $W_1$ and $W_2$, respectively, the pairs $(W_1,W_2)$ and $(S_X,S_Y)$ have the same distribution. Applying the measurable map $A$ to both pairs gives \begin{align*} W_1+W_2 \sim S_X+S_Y. \end{align*}[/step]

[guided]We now pass from marginal distributional equality to equality of the sums. Marginal equality alone would not be enough: the distribution of a sum depends on the joint law. This is exactly where independence is used. For a random matrix $R : (\Omega_R,\mathcal{F}_R,\mathbb{P}_R) \to (\mathbb{R}^{p \times p},\mathcal{B}(\mathbb{R}^{p \times p}))$, write $\mathcal{L}(R)$ for its law, meaning the pushforward probability measure $\mathbb{P}_R \circ R^{-1}$. If $\mu$ and $\nu$ are probability measures on $\mathbb{R}^{p \times p}$, write $\mu \otimes \nu$ for their product probability measure on $\mathbb{R}^{p \times p} \times \mathbb{R}^{p \times p}$. The random matrices $W_1$ and $W_2$ are independent, so the law of the pair $(W_1,W_2)$ is \begin{align*} \mathcal{L}(W_1,W_2)=\mathcal{L}(W_1)\otimes \mathcal{L}(W_2). \end{align*} The constructed matrices $S_X$ and $S_Y$ are also independent, so \begin{align*} \mathcal{L}(S_X,S_Y)=\mathcal{L}(S_X)\otimes \mathcal{L}(S_Y). \end{align*} From the previous step, \begin{align*} \mathcal{L}(W_1)=\mathcal{L}(S_X), \qquad \mathcal{L}(W_2)=\mathcal{L}(S_Y). \end{align*} Therefore \begin{align*} \mathcal{L}(W_1,W_2)=\mathcal{L}(S_X,S_Y). \end{align*} Now define the measurable addition map \begin{align*} A : \mathbb{R}^{p \times p} \times \mathbb{R}^{p \times p} &\to \mathbb{R}^{p \times p} \\ (M,N) &\mapsto M+N. \end{align*} Applying the same measurable map to two random elements with the same distribution preserves equality in distribution. Hence \begin{align*} A(W_1,W_2) \sim A(S_X,S_Y), \end{align*} which is exactly \begin{align*} W_1+W_2 \sim S_X+S_Y. \end{align*}[/guided]

[step:Identify the combined outer-product sum as a Wishart random matrix]By the definitions of $S_X$ and $S_Y$, \begin{align*} S_X+S_Y &= \sum_{i=1}^{n_1} X_iX_i^\top + \sum_{j=1}^{n_2} Y_jY_j^\top. \end{align*} Define a combined family of random vectors $(Z_k)_{k=1}^{n_1+n_2}$ by \begin{align*} Z_k : (\Omega,\mathcal{F},\mathbb{P}) &\to (\mathbb{R}^p,\mathcal{B}(\mathbb{R}^p)) \\ \omega &\mapsto \begin{cases} X_k(\omega), & 1 \le k \le n_1,\\ Y_{k-n_1}(\omega), & n_1+1 \le k \le n_1+n_2. \end{cases} \end{align*} The vectors $Z_1,\dots,Z_{n_1+n_2}$ are independent because the family $(X_i)_{i=1}^{n_1}$ is independent, the family $(Y_j)_{j=1}^{n_2}$ is independent, and the two families are independent of each other. Each $Z_k$ has distribution $\mathcal{N}_p(0,\Sigma)$. Therefore \begin{align*} S_X+S_Y = \sum_{k=1}^{n_1+n_2} Z_kZ_k^\top. \end{align*} By the definition of the central Wishart distribution, \begin{align*} S_X+S_Y \sim W_p(n_1+n_2,\Sigma). \end{align*} Combining this with $W_1+W_2 \sim S_X+S_Y$ gives \begin{align*} W_1+W_2 \sim W_p(n_1+n_2,\Sigma). \end{align*} This is the desired reproductive property.[/step]

[guided]We now identify the model sum as a single Wishart random matrix with the combined sample size. By the definitions of the maps $S_X : \Omega \to \mathbb{R}^{p \times p}$ and $S_Y : \Omega \to \mathbb{R}^{p \times p}$, for every $\omega \in \Omega$, \begin{align*} (S_X+S_Y)(\omega) &= \sum_{i=1}^{n_1} X_i(\omega)X_i(\omega)^\top + \sum_{j=1}^{n_2} Y_j(\omega)Y_j(\omega)^\top. \end{align*} Define random vectors $(Z_k)_{k=1}^{n_1+n_2}$ on the same probability space by \begin{align*} Z_k : (\Omega,\mathcal{F},\mathbb{P}) &\to (\mathbb{R}^p,\mathcal{B}(\mathbb{R}^p)) \\ \omega &\mapsto \begin{cases} X_k(\omega), & 1 \le k \le n_1,\\ Y_{k-n_1}(\omega), & n_1+1 \le k \le n_1+n_2. \end{cases} \end{align*} This definition is well-formed because all the $X_i$ and $Y_j$ have already been constructed on the common probability space $(\Omega,\mathcal{F},\mathbb{P})$. The family $(Z_k)_{k=1}^{n_1+n_2}$ is independent because it is just a relabelling of the independent combined family $(X_1,\dots,X_{n_1},Y_1,\dots,Y_{n_2})$. Each $Z_k$ has distribution $\mathcal{N}_p(0,\Sigma)$: for $1 \le k \le n_1$ this is the distribution of $X_k$, and for $n_1+1 \le k \le n_1+n_2$ this is the distribution of $Y_{k-n_1}$. With this relabelling, the two sums combine into one sum: \begin{align*} S_X+S_Y = \sum_{k=1}^{n_1+n_2} Z_kZ_k^\top. \end{align*} The defining representation of the central Wishart distribution now applies with sample size $n_1+n_2$, since the random vectors $Z_1,\dots,Z_{n_1+n_2}$ are independent and each has distribution $\mathcal{N}_p(0,\Sigma)$. Therefore \begin{align*} S_X+S_Y \sim W_p(n_1+n_2,\Sigma). \end{align*} The previous step proved \begin{align*} W_1+W_2 \sim S_X+S_Y. \end{align*} Combining these two distributional equalities gives \begin{align*} W_1+W_2 \sim W_p(n_1+n_2,\Sigma). \end{align*} This is precisely the reproductive property asserted in the theorem.[/guided]