We have a zero-mean Gaussian random field with length n and covariance matrix \bm{C}. Note that \bm{C} is the true noiseless covariance. Here we draw N=1404 realizations of the random field, and estimate a noisy covariance \bm{C}_N. As the Covariance is noisy, when estimating the inverse of the covariance matrix, we need to add a Hartlap correction, where the Hartlap factor is defined as h = \frac{N-n-2}{N-1}. We have \begin{equation} \bm{C}^{-1} = h \left\langle (\bm{C}_N)^{-1} \right\rangle. \end{equation} If we have a realization \bm{x}, the \chi^2 of this random realization is \begin{equation} \chi^2 = \bm{x}^{T} \cdot \bm{C}^{-1} \cdot \bm{x} = h ~ \bm{x}^{T} \cdot (\bm{C}_N)^{-1} \cdot \bm{x}\,. \end{equation} With many realizations of \bm{x}, we can sample the distribution of \chi^2.
However, the problem is that we do not know \bm{C} in real observations and we only have a noisy estimation \bm{C}_N. The question is — how to estimate the \chi^2 distribution for this Gaussian random field starting from the noisy covariance?
We fully reproduce the estimation process in the first section assuming \bm{C}_N is the true covariance. Note that the expectation of \bm{C_N} is the same as \bm{C}, but the expectation of its inverse is not the same as \bm{C}^{-1}: \begin{equation} \begin{split} \langle \bm{C}_N \rangle & = \bm{C} \\ \langle (\bm{C}_N)^{-1} \rangle & = \frac{1}{h} \bm{C}^{-1}\neq \bm{C}^{-1}\,. \end{split} \end{equation}
We draw N=1404 realizations of Gaussian random fields using \bm{C}_N and use them to estimate a noisy covariance \bm{C}'_N. We have \begin{equation} \bm{C_N}^{-1} = h \left\langle (\bm{C'}_N)^{-1} \right\rangle. \end{equation} For one random realization \bm{x'} drawn from the simulation with \bm{C}_N, the \chi^2 should be \begin{equation} \chi'^2 = h~ \bm{x'}^{T} \cdot (\bm{C}'_N)^{-1} \cdot \bm{x'} = \bm{x'}^{T} \cdot (\bm{C}_N)^{-1} \cdot \bm{x'}. \end{equation} By simulating different realizations of \bm{x}', we can sample the \chi'^2 distribution P(\chi^2) and use it to approximate P(\chi^2).
In addition to using \bm{C}_N to approximate the true covariance \bm{C}, we can also assume \frac{1}{h}\bm{C}_N is the true covariance. Note that the expectation of (\frac{1}{h}\bm{C_N}) is not the same as \bm{C}, but the expectation of its inverse is the same as \bm{C}^{-1}: \begin{equation} \begin{split} \langle \frac{1}{h} \bm{C}_N \rangle & \neq \bm{C} \\ \langle (\frac{1}{h} \bm{C}_N)^{-1} \rangle & = \bm{C}^{-1}\,. \end{split} \end{equation}
We draw N=1404 realizations of Gaussian random fields using \frac{1}{h}\bm{C}_N and use them to estimate a noisy covariance \bm{C}''_N. For one random realization \bm{x''} drawn from the simulation with \bm{C}_N, the \chi^2 should be \begin{equation} \chi''^2 = h ~ \bm{x''}^{T} \cdot (\bm{C}_N)^{-1} \cdot \bm{x''} = \bm{x''}^{T} \cdot (\bm{C})^{-1} \cdot \bm{x''} \end{equation} By simulating different realizations of \bm{x}'', we can sample the \chi''^2 distribution P(\chi^2) and use it to approximate P(\chi^2).