Background
Consider a multivariate Gaussian dataset $\mathbf{Y}$ with observations on $k$ individuals (rows) over $m$ variables (columns). The variables have covariance $\boldsymbol{\Sigma}$ (an $m\times m$ matrix) and the individuals have a known (expected) covariance structure described by the $k\times k$ matrix $\mathbf{V}$. Now assume $\mathbf{Y}$ is arranged into an $n$-length vector $\mathbf y$ such that its covariance is described by $\mathbf{W}=\mathbf{\Sigma}\otimes\mathbf{V}$. Assuming no missing data, $n=km$ and $\mathbf{W}$ is a block matrix of dimension $n\times n$, with block $i,j$ composed of $\Sigma_{ij}\mathbf{V}$. The log-likelihood for $\mathbf{y}$ given $\boldsymbol{\Sigma}$ and $\boldsymbol{\mu}$ is then \begin{equation} \operatorname{log}(L)=-\frac{1}{2}\left((n)\operatorname{log}(2\pi)+\operatorname{log}|\mathbf{W}|+(\mathbf{y}-\mathbf{X}\boldsymbol{\hat\mu})'\mathbf{W}^{-1}(\mathbf{y}-\mathbf{X}\boldsymbol{\hat\mu})\right) \end{equation} (for the restricted log-likelihood, replace $(n)$ with $(n-m)$ and subtract $\frac{1}{2}\operatorname{log}|\mathbf{X}'\mathbf{W}^{-1}\mathbf{X}|$).
Then $\mathbf{X}$ is the $n\times m$ matrix indicating which observations are from which variable ($\mathbf{X}_{ij}=1$ if the $i$th of variable $j$ corresponds to observation $y_u$ and $\mathbf{X}_{ij}=0$ otherwise) and $\boldsymbol{\hat\mu}$ is an $m$-length vector of estimated root states for each variable, given by \begin{equation} \boldsymbol{\hat\mu}=(\mathbf{X}'\mathbf{W}^{-1}\mathbf{X})^{-1}\mathbf{X}'\mathbf{W}^{-1}\mathbf{Y} \,.\end{equation}
Assuming complete (i.e. no missing) data, $\boldsymbol{\Sigma}$ has a familiar closed-form solution, but for my purposes I often need to numerically estimate the maximum likelihood (or restricted maximum likelihood) parameters to find $\boldsymbol{\hat\Sigma}$.
Question
Supposing $\mathbf{Z}$ is multivariate Gaussian and, like $\mathbf{Y}$, is arranged into a vector $\mathbf{z}$, but unlike $\mathbf{Y}$, $\mathbf{Z}$ is (potentially) composed of variables that may be entirely real, entirely complex, or both -- is there a way to implement a log-likelihood function to allow for maximum likelihood estimation of $\boldsymbol{\Sigma}$ for $\mathbf{Z}$ regardless of whether each variable is real or complex?
To motivate the question, and to clarify the nature of $\mathbf{Z}$, I am essentially attempting to model the roots of a set of nonlinear equations as random Gaussian variables. Sometimes the roots are all real, some roots are all complex, and I expect some might switch between the two (I'm not sure if the roots could also be purely imaginary).
I don't have a background in complex numbers, and my math/stats background is largely self-taught, so apologies if my notation is incorrect. I'm bringing my question here only after several weeks of unsuccessful attempts in finding the answer myself via the limited complex multivariate Gaussian literature that's out there. However, it does seem safe to conclude that I can't simply treat each complex variable as two real Gaussian variables and proceed with the real multivariate Gaussian distribution?
Also, I'm unsure if $\mathbf{Z}$ would meet the requirement for "circular symmetry" but it seems like the answer would be "no."
