Let $ X\in\mathbb{R}^{(N+1)n}, U\in\mathbb{R}^{Nm}$ be two random variables, with $N,n,m >0$ positive integers. The covariances of $X$ and $U$ are $\Sigma_{X} = (I + BK)\mathcal{S}(I + BK)^\intercal$ and $\Sigma_U = K\mathcal{S} K^\intercal$, respectively, and the cross-covariance between $X$ and $U$ is $\Sigma_{XU} = (I + BK)\mathcal{S}K^\intercal$, where $\mathcal{S} \in\mathbb{R}^{(N+1)n\times (N+1)n}, K\in\mathbb{R}^{Nm\times (N+1)n}, B \in\mathbb{R}^{(N+1)n\times Nm}$. For reference, all covariance matrices are positive definite. I am interesting in finding the square root matrix of the combined random vector $Z := [X^\intercal, U^\intercal]^\intercal \in \mathbb{R}^{(N+1)n + Nm}$, that is, $\Sigma_{Z}^{1/2}$.
For $X$ and $U$ individually, it is very simple to the find the square root matrix, that is, matrices $L$ and $M$ such that $LL^\intercal = \Sigma_X$ and $MM^\intercal = \Sigma_U$. It is clear that
$$ L = (I + BK)\mathcal{S}^{1/2} \quad \textrm{and} \quad M = K\mathcal{S}^{1/2}, $$
where I am assuming that $\mathcal{S}^{1/2}$ can be easily computed (via Cholesky decomposition). I am struggling, though, to compute $\Sigma_{Z}^{1/2}$. The random vector $Z$ has the covariance matrix
$$ \Sigma_{Z} = \begin{bmatrix} (I + BK)\mathcal{S}(I + BK)^\intercal & (I + BK)\mathcal{S} K^\intercal \\ K\mathcal{S}(I+BK)^\intercal & K\mathcal{S}K^\intercal \end{bmatrix}. $$
However, I am not sure how to find a matrix $P$ such that $PP^\intercal = \Sigma_{Z}$, due to the presence of the cross-covariance terms.
Solved this problem:
$$ \Sigma_{Z} = \begin{bmatrix} I + BK \\ K \end{bmatrix} \mathcal{S} \begin{bmatrix} (I + BK)^\intercal & K^\intercal \end{bmatrix} $$ is the decomposition. From here, it's straightforward to compute the square root matrix, i.e.,
$$ \Sigma_{Z}^{1/2} = \begin{bmatrix} I + BK \\ K \end{bmatrix} \mathcal{S}^{1/2}, $$ which satisfies $\Sigma_{Z}^{1/2}(\Sigma_{Z}^{1/2})^\intercal = \Sigma_{Z}$.