@User1865345 Also sounds good

@Firebug Theorem 2.2 from the paper is exactly what I was looking for. Thank you! Feel free to post it as an answer, I will accept it.

@User1865345 The key difference here is that normality has to be preserved for any normal random variable the transformation is applied to, not just one specific one. I'm not sure what you mean by weakening the requirements but I would also be interested in partial results, if you have any.

@whuber Exactly!

Unfortunately, the document is not accesible anymore. Could you give the title of the publication or provide a new link?

@whuber Could you elaborate on how to partition the matrix? All the structure it has is that it is positive semidefinite. I feel like some reparametrization would be necessary. How would that work?

@whuber While I think this is a good direcction, it's not very clear how to apply it because the standard formulas for Gaussian conditional distributions only apply if we marginalize out a subset of variables which is not the case here.

@user0 No, $\Sigma$ is an arbitrary positive semidefinite matrix

@whuber But how do I do that? $\mathbb{E}[x^2] \neq \mu^2$

@whuber Well x^2 follows a scaled non-central chi-square distribution so the expectation of $x^2$ should be $\sigma^2*(1+\mu^2)$ I think. So for $x^2 - s^2$ it is $(\sigma^2-1)\mu^2$.

@whuber Sure but how do I get an unbiased estimator of $s^2 - \mu^2$?

@whuber The problem is that I know neither $\mu$ nor $E[x]$ nor $E[X]$ and x is an unbiased estimator of $\mu$ with variance $\sigma^2$ but how can I get an unbiased estimate of $\mu^2$ or the variance?